Structure-property relations of silicon oxycarbides studied using a machine learning interatomic potential

Silicon oxycarbides show outstanding versatility due to their highly tunable composition and microstructure. Consequently, a key challenge is a thorough knowledge of structure-property relations in the system. In this work, we fit an atomic cluster expansion potential to a set of actively learned DFT training data spanning a wide configurational space. We demonstrate the ability of the potential to produce realistic amorphous structures and rationalize the formation of different morphologies of the turbostratic free carbon phase. Finally, we relate the materials stiffness to its composition and microstructure, finding a delicate dependence on Si-C bonds that contradicts commonly assumed relations to the free carbon phase.


Introduction
Silicon oxycarbides (Si-O-Cs) are highly versatile materials that combine remarkable structural and functional properties.Among them are high temperature resistance, good mechanical strength [1,2], great creep-and corrosionresistance [3,4], as well as piezoresistivity and the ability to reversibly store Li + , Na + and K + [5,6,7,8].These properties make them interesting for applications in very different fields, such as protective coatings [9], energy storage [10] and biomedicine [11].Despite the plethora of desirable properties and intensive research, open questions about structural features and their relation even to basic characteristics such as the Young's modulus (E) remain.From NMR measurements it is known that Si-O-Cs consist of corner-shared SiO x C 4−x tetrahedra with car-bidic sp 3 -hybridized carbon [12,13] and a segregated secondary phase with sp 2 -hybridized turbostratic carbon [14,15].The detailed nanostructure, however, remains elusive and sensitively depends on the composition, precursor structure and processing conditions.For example, it is still unclear in which form the turbostratic carbon is present in Si-O-C.Scarmi et al. [16] and Saha et al. [17] argued that interpenetrating networks of graphenelike carbon and silica-rich mixed tetrahedra domains are formed, based on the high creep resistance of the material.In contrast, Widgeon et al. [12] found that a model with graphitic inclusions embedded in a silica rich matrix result in a better match of the mass fractal dimensions of mixed tetrahedra.Here, atomistic simulations may help in understanding structure formation and structureproperty relations of Si-O-C compounds at the nanoscale.As shown for example in a series of studies by Kroll [18,19,20], ab-initio molecular dynamics (MD) simulations can be employed to investigate structural details, energetics and elastic properties of Si-O-C, but are limited to structures consisting of a few hundred atoms and simulation times on the order of tens of picoseconds.Consequently, the heterogeneity on the nanoscale intrinsic to Si-O-C can not be reproduced in full extend.Large-scale MD simulations, on the other hand, require suitable interatomic potentials and the complex nature of the strongly directional covalent bonds is hard to capture in empirical formulas like bond order potentials.Recent studies showed that the ReaxFF framework [21] allows to study specific aspects of the Si-O-C system.For example, Newsome et al. simulated the oxidation of silicon car-bide [22].Soria et al. investigated organic molecules on silicon surfaces [23], Gao et al. and Ponomarev et al. simulated the pyrolysis of specific polymers to amorphous Si-O-C [24,25].These potentials, however, have a limited scope and require a reparametrization for each application.Furthermore, only two out of the four parameter sets are publicly available [22,23].Here, modern machine learning interatomic potentials (MLIPs) offer an alternative approach to describe complex systems over a wide compositional and structural range at similar computational cost, but at the expanse of requiring more training data.Recent studies have shown the successful application of MLIPs to carbon [26,27,28,29,30,31,32,33], silicon [34,35,36,37,38,39], SiO 2 [40,41,42,43] and SiC [44,45].In this work, we present an Atomic Cluster Expansion (ACE) potential for the Si-O-C system fitted to an extensive database generated using the active learning (AL) capabilities of Moment Tensor potentials (MTPs) [46] and ACE [47].We show that the potential achieves a high accuracy for a wide compositional range and that it can be used to produce realistic amorphous Si-O-C structures.We investigate structural features, formation energies and Young's moduli for samples with varying compositions and precursor configurations.Thereby, we find that the structure model containing graphene like sheets and the graphitic inclusions in a silica matrix are both likely to describe the structure of Si-O-C, but occur at different stages of processing.Furthermore, we establish relations of Young's modulus to the fraction of Si-C bonds in mixed SiO x C 4−x tetrahedra and the silica volume.

Training and testing data
The training data for the potential was produced as follows.
Initial Si, O and Si-O structures were taken from an ACE potential previously fitted for the Si-O system [48].Pure C structures were generated using the ASE package [49] and SiC structures were taken from the materials project database [50].Si-O-C structures with varying compositions were generated with two different procedures.Firstly, structures were produced to match the expected coordination, i.e. 4 1: Compositions of hydrogen stripped polymers, sorted from low to high carbon content.These compositions were used to generate polymer derived (PD), bulk fragment (BF), isolated atoms (Ats) and graphite flakes and isolated atoms (GrAts) samples (See Fig. 1a).
for Si and C and 2-fold for O [51,52].Additional structures were then generated using an AL procedure, schematically shown in Fig. 1.
We employed PACKMOL [53], which implements an algorithm to densely pack structural fragments while keeping adjustable minimal distances between atoms [54].With this approach structures based on polymers inspired from Polymethylsilsesquioxane (PMSQ), the polyorganosiloxanes RD-212, RD-684 and SILRES-604 and a fully artificial polymer made from Si 2 O 2 C monomers were generated with removed H atoms.In the following we will use the compositions, instead of polymer names as shown in Tab. 1. Stripping H from the polymers avoids a massive extension of the necessary training data, as the configuration space for a quarternary system is much larger than for a ternary, but keeps the polymeric backbone.We expect that this has a negligible influence for final structures produced with pyrolysis temperatures of 1273 K to 1523 K and higher, because H evaporates during synthesis and the remaining amount of H atoms in them is very small [55,56,57].In order to cover a wider variety of compositions and structural features, we also produced configurations made from isolated atoms (Ats) and bulk fragments (BFs) SiO 2 , SiC 4 , CSi 4 and C 6 as shown in Fig. 1a.For the structures generated directly from Ats and BFs the composition was varied, with Si and C concentrations ranging from 0 % to 100 % and the oxygen concentration from 0 to 2.2 times the Si concentration.The AL capabilities of MTPs [58] and ACE potentials [47] were employed to iteratively generate new Si-O-C configurations Structures for AL and samples used in the analysis were produced using a cook and quench simulations depicted in (c).The compressions step shown as dotted line was only applied for PD and GrAts structures, where it was necessary to obtain nonporous bulk samples.For structures made from BFs or Ats it was not necessary, because the initial packing already leads to reasonable densities.
as shown in Fig. 1b.Here, we started with MTPs because the AL capabilities for ACE were implemented only recently [47].Furthermore, we covered a wider variety of structures by using both codes.The training data set was considered complete when no new structures, that had a maximum density-functional theory (DFT) force of less than 150 eV/ Å, were discovered in cook and quench simulations with temperatures up to 3000 K and pressures up to 200 GPa.For the test data set the test data of the Si-O potential [48] was supplemented with a separately created set of Si-O-C structures.These Si-O-C test structures were created with PACKMOL at varying densities and compositions.Consequently, the atoms were randomly displaced.Finally, the training and test data sets were filtered to not contain structures with a maxi-mum force of more than 150 eV/ Å, a minimal distance between atoms smaller than 0.6 Å or greater than 4 Å or an energy of more than 20 eV/atom above the convex hull.

Details of potentials
For the potentials we used a cutoff of 5 Å.The ACE and MTPs were fitted using the pacemaker [39,59] and MLIP packages [46].For the AL process with MTPs a level 26 potential was employed.The intermediate ACE potentials for AL were fitted with the triple embedding ρ 0.5 1 +ρ 1 2 +ρ 2 3 and a total of 2325 basis functions.In principle, the most straightforward way of increasing the accuracy of ACE potentials is to increase the number of basis functions.However, this considerably increases the computational cost [59].Instead, we tested different embedding terms for the fi-nal potential, which is computationally very cheap.Here, we found a highly nonlinear sum of expansions n i ϕ α i i with n = 10 and exponents α i 0.125, 0.25, 0.375, 0.5, 0.75, 0.875, 1, 1.25, 1.5 and 2 to result in the best testing errors.

SiOC sample structures
The structure samples used to analyze formation energies, structural features and elastic properties were produced in a cook and quench process, as schematically shown in Fig. 1c.Here we tested two different degrees of freedom, the influence of the composition and the effect of the precursor on the structure and properties of the final sample.The compositions of the structures correspond to the 5 polymeric compositions shown in Tab. 1.As precursors, we employed 4 different types of building blocks shown in Fig. 1a.The graphite sheets in GrAts structures consist of 160 atom 2-layer graphite.The resulting 20 structures were used in cook and quench simulations with annealing temperatures of 1000, 1500 and 2000 K to obtain a total of 60 sample structures.Here, the annealing time was 1 ns.As shown in the supplemental material (section S1) this time is sufficient to reach a steady state regarding different structural features.The employed quench rate was 1 × 10 12 K/s.PD and GrAts based structures required an additional compression step before the cook and quench process to obtain nonporous initial structures.For this purpose, they were equilibrated at 500 K with an applied isotropic pressure of 10 GPa for 10 ps before heating them up to the annealing temperature.

Simulations
DFT calculations were carried out with the same settings as used for a silica potential previously fitted [43] to keep the training data consistent.The plane wave code VASP [60,61,62] with projector-augmented wave [63] pseudopotentials and the SCAN [64] meta-GGA exchange-correlation were employed with a plane-wave cutoff of 900 eV and a k-spacing of 0.23 Å−1 .Classical MD simulations were carried out with LAMMPS [65], applying GPU accelerated KOKKOS versions where possible.If not otherwise noted an NPT ensemble with isotropic 0 Pa pressure, Nosé-Hoover thermo-and barostats and a timestep of 1 fs were employed in the simulations.

Results
In the following we will start by shortly presenting the results of the applied AL procedure and evaluate the newly developed MLIP.Then we will discuss the structure of the produced Si-O-C samples based on their composition and precursors going from low to high C contents.Finally, we will relate structural features to the Young's modulus of the samples.

Training and performance of the potential
The flexibility of the ACE formalism allows for an accurate description of highly complex materials, under the condition that similar atomic configurations have been part of the data used in their training procedure.In this work, an AL procedure was employed to ensure that the training data covers a large phase space volume.Here, we started with MTPs as implemented in the MLIP package [46], and continued with ACE potentials for which AL capabilities were implemented only recently [47].
We included structures with varying compositions and at high temperatures and pressures, as well as different defective structures that can form during MD simulations.This makes the potential applicable to a wide area of problems.Fig. 2a and 2b show the resulting distribution of the training data in terms of atomic energies and volumes.
The mostly actively learned Si-O-C structures are scattered widely in this 2D representation of phase space.The pure elements and SiC structures were mostly made by hand.In combination with less compositional degrees of freedom this leads to a comparatively narrow distribution for them.To prevent the occurrence of unphysically large forces on single atoms during MD simulations at high temperatures or pressures we found it helpful to add SiC structures with a high density, leading to a second area with aggregated SiC structures around 4 eV/atom above the convex hull.Energies and forces predicted using the potential agree well with those calculated using DFT over the whole range of structures, as shown in the scatter plots in Fig. 2c and 2d.For energies an RMSE of 24 and 36 meV/atom for the training and test sets were obtained.For forces the RMSEs were 479 and 650 meV/ Å respectively.We want to note that the higher errors for the testing set are a result from the different distribution of the structures regarding their composition and energy and not from overfitting as shown by the continuous decrease of testing errors in Fig. 2e and 2f.

Si
The microstructure and properties of Si-O-C compounds depend strongly on the precursor material and processing conditions like pyrolysis temperature [7].Time and lengths scales of the experimental processes can not be reproduced directly in MD simulations.Instead, we used different building blocks to produce a variety of microstructures in cook and quench simulations as previously described.Examples of these initial structures are shown in the upper row of Fig. 3 Exemplary samples produced via the cook and quench protocol described in section 2.3 with an annealing temperature of 1500 K are shown in the row below.Qualitatively, the PD, BF and Ats structures are very similar.Major differences are only found for GrAts.Here, large graphite areas are still present in the final structure.From a purely thermodynamic viewpoint this is surprising, as one would expect that Si 0.4 O 0.4 C 0.2 decomposes into SiO 2 and SiC, but at the tested temperatures and timescales the kinetics do not allow for such a phase separation.To quantify differences between the structures we calculated the Voronoi volume fractions occupied by 3-fold coordinated C atoms, i.e. the free carbon phase and SiO 4 tetrahedra, i.e. the silica phase.The results are shown in Fig. 4.Only the structures based on GrAts contain a meaningful amount of free carbon.As discussed previously there is no thermodynamic driving force for the formation of graphene or graphite like carbon, so its volume fraction is determined by the kinetics of the system and this behavior can be expected.In the GrAts structure the volume fraction stays constant between annealing temperatures 1000 and 1500 K.At 2000 K a slight decrease can be observed, showing that the graphite phase is kinetically stabilized up to Temperatures greater than 1500 K. Regarding the silica volume fraction a continuous increase at increasing temperatures can be observed for all structures, indicating a higher degree of phase separation.However, the data shows large differences with regard to the absolute values.The largest amount is observed for the GrAts structure, followed by BFs, Ats and finally the PD configuration.Again this observation can be rationalized by the thermodynamics and kinetics of the system.In the GrAts structure the unpaired Si and O atoms can quickly form silica without interference of the carbon, because it is already bound in the graphite phase.Si-C can only form in the relatively small interface regions or requires long range diffusion.In the BF based structure SiO 2 , SiC 4 and Si 4 C units are already present, only requiring rotations and minor rearrangement to form the thermo-dynamically favored products, while the diffusion paths necessary to achieve phase separation in the atom based structure are long and intermediate bonds can form and need to brake again.Similarly, the PD structure contains Si bonded to C and O, which needs to break before forming pure SiO 4 or SiC 4 tetrahedra.The right column of Fig. 4 shows the formation energy of samples with respect to αquartz, β-SiC and graphite.In the case of Si 0.4 O 0.4 C 0.2 they are very close in energy, despite the significant structural differences.Especially for the structures annealed at 2000 K the difference of nearly 0 for may irritate, when considering that the system still contains a considerable amount of free carbon in the case GrAts, but this energetically unfavorable state is apparently compensated by the high fraction of very favorable silica in the system.should form.Consequently, the formation energies of structures based on GrAts are considerably lower than those of the other precursors for Si 0.25 O 0.5 C 0.25 .The difference between the Ats and GrAts structures with an annealing temperature of 1000 K is about 300 meV/atom, i.e. the driving forces for phase separation are very high.Since the structures are in a steady state, this also indicates a very low mobility.For an annealing temperature of 2000 K the difference between precursors becomes much smaller, reaching around 100 meV/atom.It is expected that the differences shrink, because the increased mobility at higher temperatures allows coming closer to the thermodynamic equilibrium.The free carbon and silica volume in Si 0.25 O 0.5 C 0.25 samples increases significantly with increasing annealing temperatures.The least changes are observed for the GrAts structure, which is already close to the phase separation.Generally, the arguments regarding the kinetics of the different precursors discussed for Si 0.4 O 0.4 C 0.2 also apply for Si 0.25 O 0.5 C 0.25 , so similar trends can be observed, with the main differences determined by the decomposition products.1a), which could ease the formation of graphite.However, this can not be observed in our data.Generally the three compounds show high free carbon contents and qualitatively similar temperature dependence for both, the silica and free carbon fractions.The energetically most favorable sample is the one made from GrAts.Similar to Si 0.4 O 0.4 C 0.2 this can be explained by the high degree of phase separation, which corresponds to the thermodynamic equilibrium.Si 0.125 O 0.125 C 0.75 and Si 0.121 O 0.121 C 0.758 behave similar.

Relation to model structures
As discussed previously, two models for the nanostructure of Si-O-C exist.One suggests silica rich nanodomains, which are separated by an interconnected graphene like carbon network [16,17].The other describes the structure as graphitic inclusions in a silica rich matrix [12].The results from our simulations are shown in Fig. 5.For structures that should contain graphite upon decomposition we found the GrAts models to be energetically favorable, favoring the latter model.However, in the other samples no graphitic carbon could be found and instead graphene like layers spread through the system, pointing to the former.This allows us to conclude that both models are representative for two distinct stages.In early stages of structure formation the structure is likely described by the first model, because the slow kinetics in the system prevent the formation of graphitic inclusions.With increasing temperatures and pyrolysis times the structure evolves towards the latter model argued for by Widgeon et al., as it is lower in energy.

Elastic properties
The structure of Si-O-C is highly tunable, depending on the pyrolysis conditions and precursors.An understanding of structureproperty relations can therefore guide the search for processing routes that best match specific needs.Here, we investigate the structural features influencing the elastic properties of Si-O-C.For this purpose we calculated the elastic tensors of our samples using pymatgen [72] and derived the Young's modulus with the Lamè constants µ = C 44 and λ = C 12 , as valid for isotropic materials.Here, the assumption of isotropy is very accurate with differences between C ij that should be equivalent due to symmetry in the range of ±3 %.Only the GrAts structures show larger anisotropies, due to the limited amount of randomly orientated graphite flakes, leading to differences in the range of ±15 %.To reduce resulting errors we averaged over the supposedly equivalent directions.The deformations applied to obtain elastic tensors are not fully reversible because the relaxation of atomic positions leads to small energy barriers on the path back to the initial state that cannot be overcome in static relaxations.Corresponding energy-strain relations can be found in the supplemental material.Fig. 6 shows the dependence of E on the SiO 4 volume, the average amount of carbon in SiO 4−x C x tetrahedra and the free carbon volume in the samples.Generally the modulus ranges from 70 GPa, which corresponds to E of silica glass [73] [66,67,68,69,3,70] as collected in [70] and from [71].
Therefore, we conclude that the amount of free carbon only weakly influences E, but at similar compositions more free carbon is equivalent to lower amounts of the strong Si-C bonds, indirectly leading to a lower stiffness.

Conclusion
Using multiple iterations of AL techniques implemented for MTPs and ACE potentials we produced a highly diverse set of amorphous Si-O-C structures, spanning a large area of phase space.With the converged dataset, we fitted a nonlinear ACE potential to energies and forces calculated with DFT and showed that the potential can accurately reproduce them.Due to the diversity of the training data, the potential has a large applicability range and can describe the formation of Si-O-C compounds and their properties in a large temperature and pressure range.
Applying the potential, we produced amorphous Si-O-C samples with various compositions in a cook and quench procedure.Here, we tested the influence of different initial structures on the final configuration.We found that manually added graphite agglomerates are kinetically stabilized in MD simulation times and lead to energetically favorable states compared to graphene like sheets, if the structure contains excess C.However, their formation is kinetically hindered, and we assume that an interpenetrating network of silica rich domains and graphene like sheets is formed as an intermediate step during synthesis of amorphous Si-O-C.
To test and establish structure-property rela-tions in Si-O-C we calculated Young's moduli of the samples.We found that the silica volume fraction and the average amount of carbon in mixed SiO 4−x C x tetrahedra correlate well with the stiffness across all samples.Silica reduces the stiffness and high amounts of Si-C bonds increase it.Furthermore, the free carbon volume fraction is a good indicator of E for samples with similar compositions, but not in general.
We believe that the presented potential will allow detailed investigations of structure formation and properties in Si-O-C.In a broader context, we have shown that modern MLIPs can be employed to study multi element material systems that form highly complex structures.Therefore, they can greatly benefit the understanding of glassy and ceramic materials on the atomistic scale.

Figure 1 :
Figure 1: Si-O-C structure generation and active-learning strategy.Most initial structures for the AL process were generated by densely packing smaller building blocks using the PACKMOL program.These building blocks and the nomenclature used to describe them throughout this work are shown in (a).Si is shown in beige, O in red and C in gray.PD structures are stripped of H atoms. GrAts building blocks were used to produce large sample structures, but not in the AL process.The AL procedure employed to iteratively improve the training data is schematically shown in (b).Structures for AL and samples used in the analysis were produced using a cook and quench simulations depicted in (c).The compressions step shown as dotted line was only applied for PD and GrAts structures, where it was necessary to obtain nonporous bulk samples.For structures made from BFs or Ats it was not necessary, because the initial packing already leads to reasonable densities.

Figure 2 :
Figure 2: Training data and performance of the potential.Training data resulting from the AL procedure (a) with the black box indicating the magnified region shown in (b).Error of energies as function of the formation energy distance from the convex hull (c) and correlation of forces (d) as calculated with DFT and the ACE potential for the training and test datasets.Training and test RMSEs of energies (e) and forces (e) during optimization of the potential.

Figure 3 :
Figure 3: Si 0.4 O 0.4 C 0.2 assembled precursors and final structures.Shown are samples as packed with PACKMOL (upper row) and after the cook and quench process (lower row) with an annealing temperature of 1500 K for all structure prototypes.The PD, BF and Ats structures appear to be very similar upon visual inspection.The graphite like flakes in the GrAts based structure are still present after processing, showing that they are kinetically stabilized, as one would expect a decomposition into SiO 2 and SiC from an energetic viewpoint.
0.25 O 0.5 C 0.25 The carbon content of Si 0.25 O 0.5 C 0.25 is very close to the previously discussed Si 0.4 O 0.4 C 0.2 , but the different Si:O ratio leads to a different thermodynamic situation.In equilibrium Si 0.25 O 0.5 C 0.25 should split into a pure silica and a graphite phase, while no SiC

Figure 4 :
Figure 4: Structure analysis and formation energies of Si-O-C samples.The left column shows the free carbon volume fraction, defined as Voronoi volume of all 3-fold coordinated C atoms, the center column shows the silica volume fraction (Voronoi volume of SiO 4 tetrahedra.The right column contains the formation energies of samples with respect to α-quartz, β-SiC and graphite.Rows correspond to different sample compositions.Si 0.125 O 0.125 C 0.75 and Si 0.121 O 0.121 C 0.758 are shown in the same row, as a direct comparison of them is interesting especially in the case of PD structures due to their compositional similarity.
0.25 O 0.25 C 0.5 , Si 0.125 O 0.125 C 0.75 and Si 0.121 O 0.121 C 0.758 Si 0.25 O 0.25 C 0.5 , Si 0.125 O 0.125 C 0.75 and Si 0.121 O 0.121 C 0.758 contain Si and O in a 1:1 ratio, but, in contrast to Si 0.4 O 0.4 C 0.2 , also excess carbon.Thus we expect the formation of a varying amount of free carbon phase in the system and not just SiO 2 and SiC.The compositional similarity of Si 0.125 O 0.125 C 0.75 and Si 0.121 O 0.121 C 0.758 allows us to investigate the influence of different polymer-like precursor structures.In the Si 0.121 O 0.121 C 0.758 polymer C 6 rings are stacked very closely (cf.Fig.

Figure 5 :
Figure 5: Structure formation of free carbon phase in Si-O-C.For low pyrolysis temperatures and times graphene like carbon forms as argued by Scarmi et al.[16] and Saha et al.[17] based on high creep resistance.At higher temperatures energetically more favorable graphite like agglomerates form, which are consistent with the mass fractal dimension as found by Widgeon et al.[12].The figures show the free carbon phase as found for Si 0.25 O 0.25 C 0.5 PD sample annealed at 2000 K (a), GrAts sample annealed at 2000 K (b) and an additional PD sample annealed at 3500 K (c).Their formation energies suggests that the Si-O-C system would evolve from PD to GrAts structures, so we consider it as an intermediate step even though timescales accessible in MD do not allow a direct observation.

Figure 6 :
Figure 6: Dependence of E on structural features.Calculated E values for samples in dependence of silica volume fraction (a), average amount of carbon in mixed tetrahedra (b) and free carbon volume fraction (c).Experimental data in (c) is taken from[66,67,68,69,3,70] as collected in[70] and from[71].