Modeling phase separation in polymer-derived carbonitride ceramics through extended machine learning molecular dynamics

Abstract: Polymer-derived ceramics combine the thermal stability of ceramics with the versatile properties of carbon domains, but modeling their atomic-scale evolution during processing remains elusive due to the limitations of traditional computational methods. To address this issue, here we develop and apply a machine learning interatomic potential for silicon carbonitride-based (Si-C-N-H) systems, trained on a diversified database of over 9000 configurations -including amorphous models, high-temperature states, surfaces, and crystal structure predictions - to capture the full complexity of these materials. This potential enables large-scale molecular dynamics simulations of 8000-atom systems revealing the atomic-scale evolution of the polymer-derived ceramic during thermal treatment. A key result of this work is the occurrence of a phase separation where carbon domains progressively nucleate from the amorphous SiCN matrix during thermal processing, forming distinct graphene-like sheets while preserving the integrity of the ceramic network. The resulting models reproduce the experimental atomic pair distribution functions with exceptional fidelity, validating our approach and providing microscopic explanations for the material unique combination of ceramic and graphitic properties. In this process, defective 5- and/or 7-member carbon rings, mediate the transformation to stable 6-member aromatic structures. These findings offer new atomic-scale insights into the thermal stability and structural transformation pathways of polymer-derived ceramics, while our methodology opens avenues for studying complex amorphous systems with first-principles accuracy at experimentally relevant scales.

• The paper develops an E(3)-equivariant MACE potential trained on 9,443 diverse DFT configurations, achieving 12.4 meV/atom energy and 149 meV/Å force RMSE while maintaining consistent transferability across amorphous, high-temperature, surface, and crystalline environments.
• The paper uses 8,000-atom, nanosecond-scale molecular dynamics to show that roughly 38–40% of carbon forms graphene-like sheets chemically bonded to a mixed Si-C-N tetrahedral matrix, with the best model reaching an experimental X-ray PDF RMSE of 0.260.
• The paper identifies defect-mediated carbon growth in which transient 3-, 5-, and 7-member rings reorganize into aromatic six-member rings, while extended annealing promotes sheet growth and improves structural agreement beyond prior 400-atom first-principles simulations.

Motivation and context

Polymer-derived ceramics (PDCs) obtained by pyrolysis of preceramic polymers such as polysilazanes exhibit a nanoscale biphasic microstructure: an amorphous SiC$_x$N$_{4-x}$ tetrahedral network coexisting with a "free carbon" phase whose degree of organization increases with pyrolysis temperature. This free carbon phase governs key functional properties—electrical conductivity, diffusion-barrier behavior against crystallization of Si$_3$N$_4$/SiC, and catalytic charge transfer in metal-containing nanocomposites—yet its nucleation and growth mechanisms at the atomic scale remain poorly characterized. Experimental probes each have blind spots: TEM resolves graphitic domains but not their formation pathways, Raman cannot distinguish small sp$^2$ clusters from extended sheets, and solid-state NMR lacks long-range information.

Prior computational work was similarly constrained. Reverse Monte Carlo and Tersoff-potential studies confirmed phase separation but lacked chemical reactivity; first-principles molecular dynamics (FPMD) by Kroll on 120-atom Si$_{40}$C$_{40}$N$_{40}$ cells captured clustering but over only tens of picoseconds; DFTB studies proposed cluster-versus-layered free carbon morphologies without dynamics. The authors' own recent FPMD study of a polyvinylsilazane-derived Si$_{32}$C$_{25}$N$_{4-x}$0H$_{4-x}$1 ceramic validated 400-atom models against experimental pair distribution functions (PDFs), showing that poorly spatially correlated carbon layers best match experiment—but FPMD system sizes and timescales precluded observing domain growth. Existing machine-learning interatomic potentials (MLIPs) for Si-C-N-H based on moment tensor potentials (MTP) enabled pyrolysis simulations but were not designed to systematically investigate free-carbon nucleation, growth, or quantitative comparison of the final phase-separated structure with experiment.

Training database construction

The central methodological contribution is a training database of 9443 configurations spanning seven classes, all recomputed at a single consistent DFT level (CP2K, PBE, GTH pseudopotentials, DZVP basis, 1000 Ry cutoff, Grimme D3 dispersion):

Class	Composition	Atoms	Configurations
SiCNH-periodic	SiCNH	200–400	4033
SiCNH-periodic-4000K	SiCNH	400	1497
SiCNH-surface	SiCNH	200–400	2000
Amorphous-Carbon-200/250	C	216	249 each
Amorphous-Silicon	Si	512	250
Crystal-Structure-Prediction	SiCNH variants	46–100	1165

Three design choices are notable. First, the high-temperature class (up to 4000 K) samples far-from-equilibrium states expected during phase separation. Second, amorphous carbon models at two densities (2.00 and 2.50 g/cm$_{4-x}$2, containing sp/sp$_{4-x}$3/sp$_{4-x}$4 mixtures) and amorphous silicon ensure transferability to compositionally extreme local environments encountered in nascent carbon domains. Third, rather than using crystal structure prediction (USPEX) to find ground states, the authors exploit the structural diversity of evolutionary-search trajectories across eight compositions (reference, C-rich/poor, N-rich/poor, H-rich, C-H-rich), including fixed low-density runs at 1.25× optimized volume, to populate rare bonding motifs absent from FPMD sampling. A multidimensional-scaling projection using Valle-Oganov fingerprints confirms that the CSP configurations fill regions of configuration space disjoint from the amorphous classes.

Potential architecture, training protocol, and accuracy

The potential is an E(3)-equivariant MACE message-passing network with two layers, 128 scalar and 128 vector features, correlation order $_{4-x}$5 (four-body interactions), spherical harmonics up to $_{4-x}$6, a 5.5 Å cutoff (11 Å effective receptive field), trained in five progressive stages with a two-stage loss schedule per stage (force-weighted optimization followed by energy refinement). Twenty percent of every class is held out as a test set throughout.

On the unseen test set, the final model achieves RMSEs of 12.4 meV/atom for energies and 149 meV/Å for forces (relative force RMSE 10.2%). Per-class errors are below 10 meV/atom for energies everywhere; force errors range from 51 meV/Å (amorphous silicon) to 218 meV/Å (4000 K class) and 373 meV/Å (CSP class), but relative force errors remain within 8.5–12.6% for all classes, indicating that elevated absolute errors simply reflect higher-energy configurations rather than degraded transferability. This level of accuracy, particularly on out-of-distribution crystalline environments, supports the claim that the database diversification strategy yields a robust potential suitable for large-scale reactive dynamics.

Large-scale simulations and structural validation

Four 8000-atom models (~4.5 nm cubic cells, experimental density 2.17 g/cm$_{4-x}$7) were generated from distinct initializations: a Packmol random distribution (SiCNH8000-Random), a ceramicNetworkBuilder network favoring small carbon motifs (SiCNH8000-CNB), a network with 46 pre-inserted carbon rings representing 32.5% of carbon atoms (SiCNH8000-CL), and an extended-cycle variant annealed for 2 ns at 2200 K (SiCNH8000-CNB_2200). Simulations ran in LAMMPS with Kokkos+CUDA on NVIDIA H100 hardware, NVT ensemble, 0.5 fs timestep, Nosé-Hoover thermostat, with thermal cycles reaching nanosecond total duration—a scale inaccessible to FPMD by roughly three orders of magnitude in both size and time.

Validation against experimental X-ray PDFs shows all four MLIP-MD models outperform the previous best 400-atom FPMD model. After optimizing amplitude and $_{4-x}$8-scaling factors (the latter ~0.96, reflecting the known slight bond-length overestimation of the DFT setup and probable closed porosity in the measured sample), the PDF RMSE improves monotonically with increasing carbon organization:

Model	$_{4-x}$9 scaling	PDF RMSE
FPMD SiCNH400-CL-3	0.977	0.521
SiCNH8000-Random	0.958	0.311
SiCNH8000-CNB	0.961	0.315
SiCNH8000-CL	0.962	0.289
SiCNH8000-CNB_2200	0.961	0.260

The improvement is most pronounced at the ~3.00 Å peak (N-(Si)-N, Si-(N)-Si, C-(C-C)-C correlations), previously underestimated in FPMD models, and at medium-range correlations beyond 4 Å. The authors attribute this to ns-scale relaxation eliminating irregular environments inherited from initialization—an interpretation supported by the systematic narrowing of partial PDF peaks relative to FPMD counterparts. A caveat applies: the model density exceeds the experimental value by ~13%, tied to the same assumed underestimation of measured density due to closed porosity; this assumption is carried over from the FPMD work rather than independently verified here.

Local structure and phase separation

Across all models, 87–90% of Si atoms are fourfold coordinated, overwhelmingly in mixed SiC$_3$0N$_3$1 tetrahedra (unmixed Si-C$_3$2 ≈ 3%, Si-N$_3$3 < 10%), with 3–4% of Si involving Si-Si bonds. Nitrogen is predominantly threefold and bonded to silicon. Carbon splits between fourfold coordination (~52–58%, matrix-resident, frequently H-stabilized as C-Si$_3$4H at ~20%) and threefold coordination (~37–47%, free-carbon domains). In the best models, the C-C$_3$5 environment exceeds 15%, marking extended sheets whose edge atoms remain sp$_3$6-hybridized in C-SiC$_3$7 environments—meaning the carbon domains are chemically bonded to the matrix, not merely held by van der Waals contact. Notably, the C-Si$_3$8C environment (~10% in SiCNH8000-CNB) nearly vanishes after extended annealing, indicating limited stability of this motif. Hydrogen resides mostly on carbon (~60–65%), stabilizing matrix carbons unable to complete four Si bonds and capping network continuity.

The phase separation reaches a chemically consistent stoichiometry in SiCNH8000-CNB_2200: ~38–40% of carbon atoms reside in graphitic sheets and ~60% in the matrix, giving a matrix composition close to Si$_3$9C$_4$0N$_4$1, i.e., approximately 10 C + 20 (SiC)$_4$2(Si$_4$3N$_4$4)$_4$5. That the simulated decomposition reproduces a physically sensible SiC/Si$_4$6N$_4$7 mixture ratio constitutes an independent consistency check beyond PDF agreement.

Nucleation and growth mechanisms

Tracking carbon motifs (single atoms, dimers, linear chains, branched chains, sheets) through the thermal cycles reveals a clear trajectory: dimers and linear chains are consumed during high-temperature plateaus while single atoms (~50–60%) and sheets grow. In SiCNH8000-CNB_2200, sheet fraction rises from 6% to ~17% at 1800 K and continues growing at 2200 K until slowing after ~300 ps; the number of sheets stays constant while their size grows, implying layer growth rather than new-nucleus proliferation. Ring statistics show 3- and 4-member rings disappearing within ~50 ps, sustained 5-member rings, and rapid 6-member ring accumulation driving sheet formation.

The atomistic mechanisms identified are defect-mediated throughout:

• De novo ring formation: dispersed single atoms and dimers coalesce into linear hexamer chains via Si-C/C-N bond breaking compensated by new Si-C and Si-N bonds; chain closure forms a 6-member aromatic ring that then captures peripheral atoms.
• Ring-size correction: 7-member rings shed one carbon through a transient 3-member ring intermediate (the ejected carbon forming a C-C$_4$8Si$_4$9 tetrahedron before rejoining the matrix), while 5-member rings are expanded to 6 members by intercalation of a carbon delivered through an analogous 3-member-ring pathway.
• Sheet coalescence: adjacent sheets linked by lateral bonds merge through shared-atom reorganization, with initially out-of-plane atoms migrating to peripheries before closure into a larger combined sheet.

These pathways parallel metadynamics results for graphene nucleation on SiC(0001) surfaces, where under-coordinated carbon reorganizes through chains, defective rings, and aromatic seeds, suggesting a generic defect-mediated mechanism for sp$^2$0 carbon ordering in silicon-containing environments.

Limitations and open questions

Several caveats bound the conclusions. The phase-separated structure emerges only after prolonged annealing at 2200 K, well above the experimental pyrolysis temperature (<1400 °C), so the simulation accesses an accelerated, thermodynamically driven pathway whose kinetic correspondence to real processing is asserted by PDF agreement rather than demonstrated directly. The $^2$1-scaling factor of ~0.96 and the associated 13% density discrepancy rest on the assumption of helium-pycnometry underestimation due to closed porosity, which remains unverified. Initial-condition dependence persists: the Random and CNB models at 1800 K do not reach full phase separation, so the claimed convergence toward the experimental structure holds specifically for sufficiently long, sufficiently hot annealing from appropriately tuned initializations. Finally, the free-carbon morphology achieved corresponds to poorly correlated, non-stacked sheets; extension to turbostratic stacking in carbon-richer PDCs is identified by the authors as the natural next application but is not established here.

Conclusion

This work delivers a MACE-based interatomic potential for Si-C-N-H systems trained on a deliberately diversified 9443-configuration database combining FPMD amorphous trajectories, high-temperature states, surfaces, elemental amorphous references, and USPEX-derived crystalline diversity, achieving 12.4 meV/atom and 149 meV/Å test accuracy with uniform relative force errors across all classes. Applied to 8000-atom, nanosecond-scale MD, it produces the first atomic-scale dynamical picture of free-carbon nucleation and growth in polymer-derived SiCN, yielding phase-separated models—~40% of carbon in graphene-like sheets chemically anchored to a mixed-tetrahedra SiCN matrix—that reproduce experimental PDFs substantially better than prior FPMD models, with the best agreement (RMSE 0.260) obtained after extended 2200 K annealing. The identification of transient 3-, 5-, and 7-member rings as obligatory intermediates in aromatic 6-member ring formation provides a concrete mechanistic account of free-carbon structuring, and the methodology establishes a transferable framework for first-principles-accuracy modeling of complex amorphous ceramics at experimentally relevant scales.