Ion correlations explain kinetic selectivity in diffusion-limited solid-state synthesis reactions

MainThe demand for new inorganic materials to improve energy technologies has led to the emergence of powerful, data-driven materials design. However, realizing in silico-designed materials through inorganic synthesis remains a challenge, lagging behind organic synthesis due to the absence of a general mechanistic model for solid-state reactions1. Current state-of-the-art atomistic modelling of solid-state synthesis describes reaction behaviour in terms of bulk thermodynamic properties from high-throughput databases of density functional theory (DFT) calculations such as the Materials Project (MP)2. Prominent examples include reaction networks3 (which identify thermodynamically favourable pathways) and active learning algorithms4 (which iteratively refine synthesis recipes based on experimental feedback). However, predictions based solely on thermodynamics can be inaccurate, particularly in systems with competing phases of similar formation energies5. In such cases, the limited transport of essential constituents may prevent the formation of globally stable products, hindering the attainment of thermodynamic equilibrium. Prior attempts to understand such effects have led to the use of empirical rate expressions6,7, which fit ‘effective rate constants’ from the degree of conversion of the reactants, making them post hoc and not predictive of solid-state reaction products. Computational studies of these reactions largely rely on phase-field methods, which require solving coupled, multiphysics partial differential equations. Prior work has focused on simplified scenarios, including single-step reactions without intermediates8, sintering-driven morphology evolution9 and the growth of predefined nucleated phases in melts10. Although informative, these studies depend on experimental or modelled inputs and, thus, require well-characterized systems to resolve diffusion and nucleation kinetics11.In general, solid-state reaction kinetics can be broken down into nucleation-limited and diffusion-limited regimes1,12. Although the former provides an insight into the first phase(s) that forms, which is useful in, for example, thin-film synthesis, such information does not conclusively predict the bulk distribution of products in powder reactions, which proceeds by the diffusion-controlled transfer of the precursor constituents to the reaction zone (Fig. 1). We hypothesize that the synthesis evolution of such systems can be described as an optimization of the local energy under the time-dependent constraint of available ionic fluxes through a defective, liquid-like interphase with the same stoichiometry as candidate, nucleating phases. Here local energy denotes the energy at a prescribed reaction interface—that is, the energy that can be reduced via the reaction of solid reactants without explicitly enforcing global mass constraints. Simply put, we propose that the growth of the crystalline nucleus in a pairwise powder reaction is governed by the ionic transport of its constituents through a liquid-like interphase of the same composition as the nucleus.Fig. 1: Reactive interface between two precursors, one of which is cation A rich (α) and the other is cation B rich (β), which react to form an interphase γ.The alternative text for this image may have been generated using AI.Full size imageThe diffusion of A towards the B-rich precursor through the disordered/liquid-like region of the interface will determine the local availability of A and B ions, which, in turn, controls which phases can form and at what rate.To showcase our model, we choose the Ba–Ti–O chemical space, which exhibits particularly competitive polymorphism (Fig. 2a). Technologically relevant Ba–Ti–O ternaries (for example, BaTiO3, Ba2Ti9O20, BaTi5O11 and BaTi2O5) are typically synthesized from BaCO3 and TiO2 precursors at varying ratios. The formation of ferroelectric BaTiO3 is well studied13,14, commonly via heating mixed powders at 1,000–1,300 °C. Although BaTiO3 is the primary product, Ba2TiO4 is energetically favoured at the TiO2/BaO interface (Fig. 2a). Impurities depend on synthesis conditions, with Ba2TiO4 dominant below 1,050 K and higher temperatures (~1,200 K) increasing the BaTiO3 yield but promoting secondary phases such as BaTi2O5.Fig. 2: Thermodynamic and kinetic landscape in the Ba–Ti–O chemical space.The alternative text for this image may have been generated using AI.Full size imagea, Thermodynamic hull of stability for BaO and TiO2 at 600 °C, computed using entries from the MP. Finite-temperature effects are accounted for using a machine learning estimator for the vibrational contribution to entropy37. The blue squares (except Ba3TiO5) are phases that have been experimentally observed in solid-state reactions, and the red ones are additional phases in the MP. The shape of the hull (that is, the phases and their relative formation energies) remains unchanged for all the temperatures considered in this work. b, Calculated effective diffusion rate constant (KD), which is a measure for the average flux of all ions through the product phase of a reaction, across the BaO–TiO2 reaction interface at a series of temperatures in the range of 1,000–1,750 K, ordered by increasing Ba:Ti ratio. KD peaks for phases with intermediate Ba:Ti ratio, and drops substantially in phases with Ba:Ti > 1 at typical synthesis temperatures. Error bars, where visible, in KD are computed using the standard deviation across five MD trajectories.Source dataPrevious studies14,15, particularly recent works5,16, have analysed selectivity in solid-state synthesis reactions based on first-principles thermodynamic data. Reference 5 suggested a thermodynamic threshold of 60 meV per atom above which the initial product formation can be reliably predicted. By contrast, for systems with multiple competing phases of comparable driving forces less than the threshold, the initial product was often determined by kinetics. For the Ba–Ti–O system, the formation energy difference between the product with the highest formation driving force (Ba2TiO4) and BaTiO3 and BaTi2O5 is ~51 meV per atom and 46 meV per atom, respectively. Hence, we expect synthesis outcomes in this system to be impacted by a temperature-dependent interplay between diffusive fluxes and reaction energies. Indeed, experimentally, the first observed product is generally Ba2TiO4 (ref. 16); however, products with a lower driving force are subsequently found, such as BaTi2O5 and BaTiO3. Using Ba–Ti–O as an exacting test case, we here present a general, predictive framework for solid-state synthesis outcomes that integrates rigorously computed ionic transport and first-principles thermodynamics with a cellular automata (ReactCA), a discrete computational model in which grid cells evolve based on local neighbour interactions (Fig. 3). As shown, predictions closely match four carefully characterized experiments across time and temperature.Fig. 3: Overview of ReactCA simulation framework for simulating solid-state powder synthesis reactions.The alternative text for this image may have been generated using AI.Full size imagea, Progression of one step of the reaction via pairwise interfacial reactions. b, Schematic illustrating the main stages of the framework: (i) the user specifies synthesis recipe to study (consisting of the chemical system to be considered, precursors and their amounts, reaction atmosphere and heating profile), based on which formation energies are obtained from the MP and machine learning estimators are used to compute the vibrational free energies, melting points and diffusion rate constants, and reactions are then enumerated and scored; (ii) a random initial distribution of particles (voxels) is generated, and the evolution rule is repeatedly applied using the scores assigned to each reaction to simulate the reaction; (iii) simulation steps are concatenated into a trajectory, which is analysed to determine phase fractions and the pathway followed over the course of the reaction19.Kinetics in the Ba–Ti–O phase spaceWe consider two reactants (here BaO and TiO2) and nine possible solid-state reaction products in the Ba–Ti–O space (Fig. 2a, blue squares). We use BaO instead of BaCO3, as it is recognized that the experimental precursor BaCO3 decomposes to BaO at ~1,100 K (ref. 17), before any ternary reaction occurs13. We calculate the flux of constituent ions from the chemical potential difference across the interface and transport coefficients of Ba2+, Ti4+ and O2− using Onsager analyses. These analyses are based on 5-ns molecular dynamics (MD) trajectories generated by machine learning interatomic potentials trained on ab initio molecular dynamics (AIMD) data for each liquid-like, non-crystalline analogue of nine possible products in the Ba–Ti–O system.Using ionic fluxes of both Ba2+ and Ti4+, effective diffusion rate constants (KD) are estimated for each considered liquid-like interphase product (Methods). Figure 2b shows KD for nine possible product compositions across the BaO/TiO2 interface for temperatures of 1,000–1,750 K. Interestingly, above 1,000 K, KD in Ti-rich phases is more than an order of magnitude higher than that in the Ba-rich phases. In the Ti-rich phases, KD increases by around an order of magnitude with every 250 K rise in temperature, and plateaus to a similar value for most phases at 1,750 K. On the other hand, KD only increases by an order of magnitude in Ba-rich phases (Ba:Ti ratio > 1) on raising temperature by 750 K.Kinetics-informed cellular automaton simulations of Ba–Ti–O solid-state synthesis reactionsTo compare temperature and time-dependent evolution of products with synthesis experiments reported elsewhere13,15,16,18, we use the recently developed cellular automaton simulation framework ReactCA19. ReactCA was designed to simulate solid-state reactions by modelling a three-dimensional grid of cells that evolve based on customizable local rules, incorporating both thermodynamic and kinetic inputs through a proxy for the reaction rate (a ‘scoring’ function). Here we extend the ReactCA framework by allowing the scoring function to depend on the instantaneous growth rate, which is a function of a calculated effective ionic diffusion constant of the liquid-like product phase (KD) at temperature T, a modified thermodynamic driving force (\(\frac{\Delta {G}^{* }}{{k}_{{\rm{B}}}T}\)) and a heuristic for Tammann’s rule20. Below the Tammann temperature, reaction rates are low but possible. Above it, rates increase with temperature due to both diffusion and thermodynamic contributions; however, at high temperatures, the saturation of diffusion rates shifts the balance in favour of thermodynamics-controlled outcomes. We simulate precursor Ba:Ti stoichiometries from 1:5 to 1:1 using the same heating profiles as the experiments (Figs. 4 and 5)15,16,18,21. Simulations using a scoring function that excludes the effective diffusion rates, that is, reaction rates based only on thermodynamics and Tamman’s rule are shown in the same figure for comparison.Fig. 4: Reaction simulations using BaO and TiO2 precursors with different ratios using ReactCA, informed by both diffusive fluxes and reaction thermodynamics versus only thermodynamics.The alternative text for this image may have been generated using AI.Full size imagea,b, Reaction 1 (BaO + TiO2) (a) and reaction 2 (BaO + 5TiO2) (b). The shade of the grey background distinguishes the different regimes of reaction selectivity: the light grey region signifies the activation-controlled regime, the darker grey region signifies the kinetics-controlled regime and the dark region signifies the thermodynamics-controlled regime.Source dataReaction 1 (Fig. 4a) represents the prototypical 1:1 BaO to TiO2 synthesis reaction, involving heating to 1,500 K followed by a brief annealing period16,21, resulting in BaTiO3 as the predominant product. In our simulation, the first products to form at low temperatures are Ba2TiO4 and BaTiO3. This agrees excellently with experimental results showing Ba2TiO4 as the primary impurity (up to 27 mol%)13,16 when the reaction was performed between 1,000 and 1,050 K. At higher temperatures, we find that Ba2TiO4 is consumed in favour of BaTi2O5, which emerges as the major impurity past 1,200 K, matching the data in refs. 15,16. With continued heating, BaTi2O5 is mostly consumed but persists as the primary impurity (100 min at 1,300 K, giving almost phase-pure BaTiO3. Performing the same simulation without incorporating diffusion rates yields a qualitatively different outcome. Results from simulations that incorporate only the thermodynamics-based data as well as Tamman’s rule show a mixture of BaTiO3 (~70 mol%), Ba2TiO4 (~15 mol%) and BaTi2O5 (~20 mol%) after annealing at 1,500 K, with Ba2TiO4 and BaTi2O5 amounts peaking at similar temperatures. In particular, BaTiO3 accounts for only 70% of the product, in stark contrast to the nearly phase-pure BaTiO3 observed experimentally.Reaction 2 (Fig. 4b) follows the experiments conducted earlier22 in which a 1:5 ratio of BaO and TiO2 was reacted at a series of temperatures (1,250–1,500 K). At low temperatures ( BaTi4O9 > Ba2Ti9O20. In our simulation (Fig. 5a), small amounts of BaTiO3 and BaTi2O5 initially form, and quickly convert to Ba2Ti9O20, BaTi5O11 and BaTi4O9. In particular, our results successfully reproduce the experimentally observed trend in the initial formation rates, following the sequence BaTi5O11 > BaTi4O9 > Ba2Ti9O20. Furthermore, and in qualitative agreement with experiments, the growth of BaTi5O11 and BaTi4O9 stagnates during the 1,400 K anneal and finally declines above 1,400 K, whereas Ba2Ti9O20 is continually produced at a positive rate. In ref. 18, trace quantities of Ba4Ti13O30 and some unreacted rutile (TiO2) were also detected when the reaction was performed below 1,400 K, which are also observed in our simulations.Fig. 5: Reaction simulations using BaO and TiO2 precursors with different ratios using ReactCA, informed by both diffusive fluxes and reaction thermodynamics versus only thermodynamics.The alternative text for this image may have been generated using AI.Full size imagea,b, Reaction 3 (2BaO + 9TiO2) (a) and reaction 4 (BaO + 2TiO2) (b). The shade of grey background distinguishes the different regimes of reaction selectivity: the light grey region signifies the activation-controlled regime, where kinetics is slow for all compositions; the darker grey region signifies the kinetics-controlled regime, where kinetics dominates the phase selectivity among thermodynamically competitive phases; and the dark region signifies the thermodynamics-controlled regime, where kinetics is fast for all compositions.Source dataReference 24 provided a recipe to synthesize phase-pure BaTi2O5 by annealing a 1:2 BaO:TiO2 precursor ratio at three successive temperature steps: ~1,200 K, 1,300 K and 1,400 K, which we simulate through reaction 4 (Fig. 5b). After the first step, their sample primarily contained BaTiO3 and some amount of BaTi2O5. After the second annealing step, the BaTiO3 content was reduced and BaTi2O5 increased proportionately. After the third annealing step, close to phase-pure BaTi2O5 was observed, with BaTiO3 and Ba6Ti17O40 appearing as impurity phases. On longer heating time above 1,400 K, BaTiO3 and Ba6Ti17O40 (the equilibrium phases above 1,473 K) increase in phase fraction24. Our simulations show excellent agreement; we show a rapid formation of BaTiO3 as the majority phase in the first annealing step, which is overtaken by BaTi2O5 in the second annealing step and we recover phase-pure BaTi2O5 with trace ( 1,700 K), close to the melting point of the system, kinetics is fast for all ionic species, and the mixture of phases closest to the globally set composition on the thermodynamic hull form to establish thermodynamic equilibrium. This insight agrees with existing theories of solid-state reactivity1,32,33,34. For example, BaTiO3 and BaTi2O5, which exhibit large negative formation energies on the Ti-rich side of the hull and relatively fast kinetics (Fig. 2), tend to form at lower temperatures (reactions 2 and 3). Raising the temperature above 1,700 K, corresponding to the thermodynamic regime, promotes the formation of BaTi5O11 and Ba2Ti9O20, which is thermodynamically the most favourable product mixture for the given precursor ratio.Finally, we assess the impact of structural similarity in solid-state synthesis. Several prior works5,18,35 argue that structural similarity or ‘templating’ between the precursors, intermediates and products promote the formation of specific phases. Indeed, in ref. 18, it was hypothesized that Ba2Ti9O20 forms above 1,373 K only when the intermediate BaTi5O11 provides a lowering of interfacial and strain energies by the templating of the Ba2Ti9O20 phase onto BaTi5O11. Here we reproduce the product distribution of ref. 18 without incorporating any structural similarity effects. Specifically, we find similar relative rates of formation for the three ternary phases in the final product distribution (BaTi5O11 > BaTi4O9 > Ba2Ti9O20) during the early stages of reactions 2 and 3 (Figs. 4 and 5) and that, importantly, Ba2Ti9O20 is only thermodynamically stable above 1,373 K, as indicated by the BaO–TiO2 phase diagram36. The modest driving force for its formation (Fig. 2) accounts for its appearance only at temperatures above 1,373 K in both experiments in ref. 18 and our simulations, as well as its absence in the shorter reaction times used by others16. To generally analyse the structural effects in the Ba–Ti–O system, we included the nucleation rate computed using the approach in ref. 35 Their metric shows that BaTiO3 and BaTi2O5 exhibit the lowest nucleation barriers and should dominate at low to moderate temperatures (Supplementary Table 1). However, Ba2TiO4 often forms as a low-temperature intermediate in both experiments and simulations, which cannot be explained by this approach. Hence, although structural similarity serves as a valuable conceptual framework for synthesis reactions that proceed topochemically, we find that the proposed approach better explains product outcomes when the reactants undergo substantial structural transformation and/or decompose entirely. We expect our method to be well suited to study this thermodynamics–kinetics interplay in more diffusion-limited solid-state synthesis reactions, with the ability to perform such simulations in high throughput to search for better recipes of both existing and unrealized materials that are possibly kinetically prevented from forming.MethodsTo capture ion dynamics effects in dense liquid-like or non-crystalline phases, we adapt the framework of ref. 38 to obtain both self- and correlated diffusive transport coefficients in the form of the Onsager transport matrix. Robust estimates of transport coefficients are obtained from nanosecond-long MD trajectories, using an atomic cluster expansion (ACE)-based machine learning interatomic potentials39, trained on 75 AIMD trajectories from non-crystalline systems. The resulting kinetic data, as well as thermodynamic data from the MP, are implemented in the ReactCA simulation framework19, which allows for describing both spatial and temporal phase evolution over the course of a prescribed solid-state reaction. Below, we describe the components of the simulation framework in more detail.Accelerating AIMD using ACEWe use the MPMorph40,41 workflow, as implemented in the atomate2 Python package, to generate amorphous configurations for training ACE. Sample random amorphous structures of the desired composition are first generated using packmol. The MPMorph workflow begins by adjusting the volume of this structure to 0.8 and 1.2 times the initial volume, followed by a 4-ps NVT AIMD run to determine the equation of state at the given temperature. A trial volume is then selected based on the equation of state for another 4-ps NVT AIMD run to check for energy convergence. If the energy converges, a 20-ps AIMD ‘production’ run is performed to equilibrate the structure at this volume. If not, the workflow continues to iteratively rescale the volume until energy convergence is achieved. Once a converged volume is found, it is used for the 20-ps production run. The workflow is run on high-performance computing resources using the jobflow and fireworks Python packages.Compositions according to BaTiO3, Ba2TiO4, BaTi2O5 and BaTi5O11 were run at 1,000, 1,250 and 1,500 K each, and frames were sampled every 100 fs to ensure sufficient variety in the sampled configurations. The sampled frames were used in addition to MP data in the Ba–Ti–O chemical space to train ACE potentials using the pacemaker42,43 Python package. The ACE potential is parameterized with 500 basis functions per element, considering neighbourhood lists within 5 Å. Following ref. 44, a higher-order Finnis–Sinclair embedding function is used. Training is performed by hierarchically growing the potential using a power-order based ‘ladder-fitting’ approach43 with a higher weight on the loss contribution from forces (99%) for 2,000 iterations. Once convergence of the training was obtained, the potential was fine-tuned by reducing the force loss weight to 90% for another 2,000 iterations and keeping the shape and complexity of the potential fixed.Active learning was performed using the extrapolation grade and high-temperature MD approach described in ref. 45. We generate randomly packed structures using packmol for compositions corresponding to BaTi5O11, Ba2Ti9O20, BaTi4O9, Ba4Ti13O30,Ba6Ti27O40, BaTi2O2, BaTiO3, Ba2TiO4 and Ba3TiO5 as input for ACE-based NVT MD simulations for 100 ps (100,000 steps) at 2,000 K. We use the extrapolation grade, which quantifies the deviation of the test configurations sampled through ACE-MD from those encountered in the training data to filter out structures that are sufficiently different from the training data. During the high-temperature MD simulations, all configurations encountered by the model with an extrapolation grade of more than 5 are considered for active learning, and the MD simulation is terminated if the model encounters a configuration with extrapolation grade exceeding 100. A subset of the collected structures is determined using the D-optimality criterion45 to ensure diversity in the active learning data, and DFT single-point calculations are performed to compute the DFT energy and forces. The trained potential is then retrained with the active learning data for an additional 1,000 iterations, and this active learning loop is repeated four times. In total, including the data acquired during active learning, the training dataset consisted of 4,628 structures containing 492,216 atoms.It is important to highlight that the objective here is not to develop a potential capable of addressing a broad variety of tasks (such as crystal relaxation, property calculations and so on) within the Ba–Ti–O chemical space, as has been successfully achieved with simpler systems. Specifically, the goal of the trained potential is to obtain longer-time NVT MD trajectories with similar-to-DFT accuracy on the amorphous liquid-like configurations. With this purpose in mind, benchmarks for the quality of the trained model are discussed in Supplementary Section 3.Onsager transport framework for a solid-state reaction interfaceIn the framework of linear irreversible thermodynamics, the diffusive flux induced in an ion i due to a driving force can be written as27,38$${J}_{i}=-\mathop{\sum }\limits_{j}{L}_{ij}\nabla {\widetilde{\mu }}_{j}.$$ (1) This expression considers the effect of driving forces (\(\{\nabla {\widetilde{\mu }}_{j}\}\)) on all ions present in the system to the flux experienced by ion i through the Onsager transport matrix (L = [Lij]). The transport coefficient Lij is then computed from an MD trajectory using the differential form of the Green–Kubo relations27,38:$${L}_{ij}=\frac{1}{6{k}_{{\rm{B}}}TV}\mathop{{\lim}}\limits_{t\to \infty }\frac{\mathrm{d}}{{\mathrm{d}}t} .$$ (2) Here V and T are the volume and temperature of the system, respectively; \({{\bf{r}}}_{\alpha }^{i}(t)-{{\bf{r}}}_{\alpha }^{i}(0)\) is the displacement of the αth particle of species i at time t. The self-transport coefficient, \({L}_{ii}^{{\rm{self}}}\), can be computed using a similar relation:$${L}_{ii}^{\mathrm{self}}=\frac{1}{6{k}_{{\rm{B}}}TV}\mathop{{\text{lim}}}\limits_{t\to \infty }\frac{{\rm{d}}}{{\rm{d}}t} .$$ (3) The self-transport coefficient of an ion i is related to the Nernst–Einstein estimate of diffusivity, also called the self-diffusion coefficient (Di) through the following relation:$${L}_{ii}^{{\rm{self}}}=\frac{{{D}}_{i}{c}_{i}}{{k}_{{\rm{B}}}T},$$ (4) where ci is the concentration of species i. This self-term exactly describes ionic mobility for dilute systems, where diffusion is ideal. Practically, Lij is the slope of a linear fit (the diffusive regime) of the time correlation function, \(\), with time t. Hence, a linear regime of at least 200 ps (200,000 steps) is used in the time correlation function versus time plot to compute all the transport coefficients. More details on how the transport coefficients are fit, as well as other caveats when using this theory on solid-state reactions, are provided in Supplementary Section 4.Deviations in diffusion due to correlated ion movementThe deviation in effective diffusion coefficients can be studied through the cross-ion correlations, which are quantified by the off-diagonal terms of L, and the distinct ion correlations, which are computed as$${L}_{ii}^{{\rm{distinct}}}={L}_{ii} - {L}_{ii}^{{\rm{self}}}.$$ (5) In this work, we focused on the effect of correlations amongst the cation pairs Ba–Ba and Ti–Ti through the corresponding distinct ion correlation terms. To compare the correlation terms between phases, we calculate the distinct correlation normalized by the self-part of the correlation function, which we call the correlation factor (f):$${f}_{i}=\frac{{L}_{ii}^{{\rm{distinct}}}}{{L}_{ii}^{{\rm{self}}}}=\frac{{L}_{ii}}{{L}_{ii}^{{\rm{self}}}}-1.$$The correlation factor typically exhibits values between –1 and +1, with values close to 0 implying uncorrelated motion between ions of type i. Additionally, we analyse the local coordination environment within the liquid-like phases by counting all geometrically feasible bonds for each ion within a cut-off radius of 5 Å. This approach captures spatial correlations or bonding interactions within both first and second coordination shells. Both shells are included because cations typically interact through bridging O ions. This enumeration is performed using the CrystalNN algorithm implemented in pymatgen46,47.Estimating reaction kinetics from atomistic transportThe liquid-like structures generated when training the ACE potential serve as inputs for NVT MD simulations using the ACE potential. These simulations are conducted for 5 ns with a 1-fs time step at temperatures of 750, 1,000, 1,250, 1,500 and 1,750 K, using the Langevin thermostat with a friction factor of 0.01 fs−1. To obtain better statistics on the computed transport coefficients, five frames from the production run of the ACE-based MPMorph flow are used as starting configurations for the NVT MD runs, leading to five replicate MD runs for each composition.Following ref. 48, as the first order, we approximate the driving forces \(\nabla {\widetilde{\mu }}_{j}\) in equation (1), which are electrochemical potential gradients, by the bulk chemical potential gradient for species j, ∇μj. These gradients are obtained through the relevant chemical potential diagram for the chemical system, which is built using pymatgen47 with data for all crystalline phases from the MP (Supplementary Fig. 5). We approximate the gradient in chemical potential for species j between the precursors (or in general, the reactants) of a solid-state reaction by the shortest distance on chemical potential diagram between domains corresponding to the precursors (reactants), say α and β, divided by the thickness of the product layer (h in Supplementary Fig. 8):$$\nabla {\widetilde{\mu }}_{j}\approx \nabla {\mu }_{j}=\frac{\min \left(\,{\mu }_{j}^{\alpha }-{\mu }_{j}^{\beta }\right)}{h}.$$ (6) This approach ensures that when the phases are in thermodynamic equilibrium, the reaction driving force and the resulting flux are zero30. The effective ‘diffusion rate constant’ (KD) is defined as follows:$${K}_{{\rm{D}}}=\mathop{\sum }\limits_{i}\mathop{\sum }\limits_{j}\frac{\left|{L}_{ij}^{\gamma }{V}^{\,\gamma }\times \min \left(\,{\mu }_{j}^{\alpha }-{\mu }_{j}^{\,\beta }\right)\right|}{{n}_{i}^{\gamma }{N}_{{\rm{A}}}}.$$ (7) where Vγ is the molar volume of the product γ, \({n}_{i}^{\gamma }\) is the fraction of atoms of species i in γ and NA is the Avogadro number. We assume a core–shell model for the formation of the product (the derivation is shown in Supplementary Section 6), with the growth of γ supported by the transport of Ba2+ and Ti4+ across the interface. In particular, the growth rate varies with the thickness of the product phase, and decreases as the reaction proceeds. For a purely diffusion-limited case for the geometry shown in Supplementary Fig. 8, the rate equation can be solved to give49$$1-{(1-y)}^{1/3}=\frac{\sqrt{2{K}_{{\rm{D}}}t}}{{r}_{0}},$$ (8) where y (1 ≥ y ≥ 0) is the degree of completion of the reaction and r0 is the radius of the precursor powder particle, reminiscent of the Jander equation, which is commonly used to fit kinetic models to the solid-state reaction data6,7. In this study, we assume that powder particles are well mixed and uniform in size, allowing us to analyse the kinetic feasibility of a solid-state reaction using the term \(\frac{\sqrt{{K}_{{\rm{D}}}}}{{r}_{0}}\), where \({r}_{{\rm{B}}{\rm{a}}{\rm{O}}}={r}_{{\rm{T}}{\rm{i}}{{\rm{O}}}_{2}}={r}_{0}\). The parameter r0 is set uniformly for all reactions in our simulations. Specifically, we assign r0 a value of 10−1 µm, which is representative of particle sizes commonly used in solid-state synthesis.We emphasize that this formulation is not applicable for all solid-state reactions, as the \({\widetilde{\mu }}_{j}\approx {\mu }_{j}\) approximation does not hold when species exhibit variable oxidation states, or there exists a substantial amount of charge transfer between species, both of which are known to occur in many solid-state reactions involving transition metals. In the present study, all phases considered have species in only a single oxidation state: +2 for Ba, +4 for Ti and –2 for O. Furthermore, the phenomenological expression for the flux (equation (1)) is defined under the centre-of-mass frame of reference, which imposes constraints on the total number of independent material fluxes in the system for incompressible systems. These constraints are described in Supplementary Section 4.Simulation of solid-state synthesis reactions with ReactCAA cellular automaton framework was recently developed (ReactCA)19 to simulate the evolution of mixtures of powders over the course of solid-state reactions. The reaction vessel is discretized into voxels (cells) and each voxel represents a powder particle of a given phase. At a given temperature, all possible reactions that can occur between any other powders (phases) in the system are enumerated combinatorially with an approach developed as part of a previous work3, and are scored using a scoring function. The system is initialized with a random distribution of two phases (in this case, representing two well-mixed precursor powders of uniform size). At every step of the simulation, two reactants are selected, and a reaction is probabilistically selected from the enumerated reactions based on their scores. Then, the reactants are probabilistically replaced with products based on the stoichiometry of the reaction. This process is repeated millions of times throughout the simulation to accurately model the evolution of phases within the reaction vessel. In the original publication19, reaction scores (that is, rates) between neighbouring particles were estimated based on the reaction thermodynamics (as estimated by zero-temperature formation energies from the MP in conjunction with a machine learning estimate of the vibrational Gibbs free energy37) and machine learning estimates of the melting point of the precursor materials50 (a preliminary proxy for the kinetic facility of the reaction). In this work, given a reaction α + β → γ, we incorporate the diffusive fluxes using a new scoring function S for ReactCA, as follows.$$\begin{array}{l}\begin{array}{l}S={\sigma }_{1}\left(\frac{{K}_{\mathrm{D}}}{{r}_{0}^{2}s}\times \frac{\Delta {G}^{\ast }}{{k}_{\mathrm{B}}T}\right)\times {\sigma }_{2}\left(\frac{T}{{T}_{{\rm{m}},\mathrm{reactant}}}\right)\\ {\sigma }_{1}(x)=\frac{1}{3}\log [1+\exp (ax)]\\ {\sigma }_{2}(x)=\frac{1}{2}\log [1+\exp (bx-c)]\\ \Delta {G}^{\ast }=1+\mathrm{erf}(-d(\Delta {G}_{{\mathrm{rxn}}}+e))\\ {K}_{\mathrm{D}}=\mathop{\sum }\limits_{i}\mathop{\sum }\limits_{j}\frac{\left|{L}_{ij}^{\gamma }{V}^{\,\gamma }\times {\mathrm{min}}\left(\,{\mu }_{j}^{\alpha }-{\mu }_{j}^{\,\beta }\right)\right|}{{n}_{i}^{\gamma }{N}_{\mathrm{A}}}\end{array}\end{array}$$ (9) The form of this scoring function is an extension to the scorer present in ReactCA19, using the same constants a, b, c, d, e. The factor s is set equal for all reactions considered, to appropriately scale the scores to be stable in the CA simulation and to ensure the reaction between the precursors initiate at temperatures above the decomposition temperature of BaCO3 (s = 5 × 1013). This scoring function has a ‘soft’ activation towards Tammann’s rule20, implemented through the σ2 part of the function, and a ‘switching-on’ behaviour towards reactions with negative reaction free energies (ΔGrxn) through the error function. For more details on these aspects of the scoring function, we refer to ref. 19. In this work, we implement explicit, quantitative transport diffusion constants in the liquid-like product phase (KD) through a soft activation with the σ1 function, which is also a soft-plus function. This function also takes the modified ΔGrxn and the temperature into account. Hence, below the Tammann temperature (\(\frac{T}{{T}_{{\rm{m}},{\rm{r}}{\rm{e}}{\rm{a}}{\rm{c}}{\rm{t}}{\rm{a}}{\rm{n}}{\rm{t}}}}\le 0.67\)), the scores are non-zero, but only slightly positive. Above the Tammann temperature, reactions with a slightly positive ΔGrxn or with low ionic fluxes have a non-zero but small score assigned to them. With increasing reaction temperature, the σ2 part of the function increases, thereby elevating the reaction rate. Similarly, raising the temperature also increases the diffusion rate constant KD (often in a superlinear fashion), which generally boosts the output of the σ1 part of the scoring function. However, as demonstrated in this work, KD approaches a saturation point with increasing temperature, which, in turn, causes the output of σ1 to also level off or even decrease at very high temperatures. This effect leads to the thermodynamic regime of phase selectivity when the reaction temperatures are close to or above the melting point. Finally, the scoring function maintains backward compatibility with the old scoring function, which is particularly important when kinetic transport parameters are not available for some or all phases in the system.DFT calculationsDFT single-point calculations and AIMD calculations were performed using the Vienna ab initio simulation package51,52 and the Perdew–Burke–Ernzerhof53 formulation of generalized gradient approximation with projector augmented-wave potentials54,55. Furthermore, to minimize computation requirements, all AIMD calculations were performed using the Γ point only and were non-spin polarized. For more details on the parameters used, we refer to the MPMorphMDSet class in atomate2. To construct the high-temperature phase diagrams, a machine learning estimator for the vibrational contribution to entropy37 was used to estimate the finite-temperature vibrational Gibbs free energies for all phases.