Assessing Thermodynamic Selectivity of Solid-State Reactions for the Predictive Synthesis of Inorganic Materials

Synthesis is a major challenge in the discovery of new inorganic materials. Currently, there is limited theoretical guidance for identifying optimal solid-state synthesis procedures. We introduce two selectivity metrics, primary and secondary competition, to assess the favorability of target/impurity phase formation in solid-state reactions. We used these metrics to analyze 3520 solid-state reactions in the literature, ranking existing approaches to popular target materials. Additionally, we implemented these metrics in a data-driven synthesis planning workflow and demonstrated its application in the synthesis of barium titanate (BaTiO3). Using an 18-element chemical reaction network with first-principles thermodynamic data from the Materials Project, we identified 82985 possible BaTiO3 synthesis reactions and selected 9 for experimental testing. Characterization of reaction pathways via synchrotron powder X-ray diffraction reveals that our selectivity metrics correlate with observed target/impurity formation. We discovered two efficient reactions using unconventional precursors (BaS/BaCl2 and Na2TiO3) that produce BaTiO3 faster and with fewer impurities than conventional methods, highlighting the importance of considering complex chemistries with additional elements during precursor selection. Our framework provides a foundation for predictive inorganic synthesis, facilitating the optimization of existing recipes and the discovery of new materials, including those not easily attainable with conventional precursors.


Introduction
The predictive synthesis of inorganic materials remains a grand challenge in modern chemistry and materials science. 1 Unlike organic synthesis, which is often described via discrete reaction steps or mechanisms, inorganic materials synthesis reactions cannot be deconstructed into elementary steps, 2,3 hindering the analogous development of retrosynthetic analysis techniques 4 and computer-aided synthesis planning. 5 This lack of successful mechanistic models has made the synthesis of predicted new materials a critical bottleneck in high-throughput computational materials design efforts, 6 with many proposed materials having yet to be successfully synthesized. [7][8][9] While there are numerous inorganic synthesis methods (e.g., hydrothermal, mechanochemical, sol-gel, intercalation, etc.), 10 we limit the scope of this work to bulk solid-state synthesis via powder reactions. This choice has been motivated by the straightforward and scalable nature of working with bulk powders, which makes solid-state synthesis suitable for applications ranging from one-off laboratory synthesis to industrial mass manufacturing. In powder reactions, product formation proceeds via nucleation and growth at interfacial contact areas in the powder mixture (Figure 1a). 11 The equilibrium phases of the reacting system can be predicted by constructing a convex hull in free energy and composition space, where the composition axis is a mixing ratio between the two precursor compositions (Figure 1b). 12 Here, we calculate the convex hull exclusively using normalized compositions and energies (i.e., on a per-atom basis). This construction, which we refer to herein as the "interface reaction hull", is a subsection of the compositional phase diagram for binary systems and a "quasibinary" two-dimensional slice of the full phase diagram for chemical systems with three or more elements. The exact product species, and the sequence in which they appear, cannot be predicted with thermodynamics alone; to do so requires intimate knowledge of the kinetic rates of all physically feasible reactions. However, a commonly adopted theoretical simplification assumes that the reaction product(s) with the most negative pairwise reaction energy will be the first to nucleate and grow as a powder mixture is heated. 11,13 This hypothesis is based on two principles: 1) the random packing of solid crystallites results in very few locations where three or more particles are simultaneously in contact, and 2) the activation energy barrier to nucleation scales inversely with free energy as ∆G † ∝ γ 3 ∆G 2 rxn . The surface energy, γ, is particularly important when comparing the feasibility of reactions with similar energies. 14 Despite this, many solid-state reactions are likely not nucleation-limited but rather transport-limited due to the relative sluggishness of solid-state diffusion and the often large driving forces of these reactions (10-100 kJ/mol). 15 Figure 1: Modeling chemical reactions at heterogeneous solid interfaces in a binary/quasibinary chemical system. (a) Cartoon model of a powder reaction between the hypothetical precursors: α (orange) and β (blue). The nucleation of a new phase, γ (green), occurs at the α|β interface according to the reaction α + β −−→ γ. (b) A possible reaction pathway for the powder system in which a secondary reaction of the equilibrium phase (γ) yields an impurity phase, δ (red). The interface reaction hull (top) shows available interfacial reactions and their corresponding Gibbs free energies, G, and mixing ratios, x. The one-dimensional spatial model (bottom) shows reaction steps beginning from an equal mixture of α and β. The impurity phase, δ, may be kinetically retained in a local equilibrium state; however, with infinite time, the system should approach the global equilibrium state composed entirely of γ.
The complex interplay between thermodynamics and kinetics makes solid-state synthesis prone to the unpredictable and undesirable formation of impurity phases. 13 A classic example of a nonselective synthesis is that of the prototypical multiferroic bismuth ferrite, BiFeO 3 , via the standard reaction from binary oxides: Bi 2 O 3 + Fe 2 O 3 −−→ 2 BiFeO 3 . This reaction typically yields impurity phases Bi 2 Fe 4 O 9 and Bi 25 FeO 39 , which are challenging to isolate and remove. 16,17 Unfortunately, the presence of an impurity phase is difficult to predict a priori, and is typically attributed to "kinetic" factors or changes in phase equilibria related to precursor purity, morphology, volatility, or processing conditions. To optimize the performance of solid-state reactions and maximize conversion to the desired target, the experimentalist frequently relies on intuition and heuristic rules to choose the 1) precursor compositions (typically off-the-shelf binary phases such as carbonates, oxides, etc.), 2) grinding/milling protocol, 3) synthesis annealing temperature, 4) synthesis atmosphere (e.g., vacuum, flowing O 2 ), 5) synthesis time, and 6) cooling protocol. Heuristics include well-known rules such as Tamman's Rule for estimating reaction onset temperature (i.e., two-thirds the melting temperature of the precursor with the lowest melting point), 18 as well as "chemical intuition" or human-biased experimental protocols (e.g., selecting synthesis times based on common increments, such as four hours, eight hours, etc.). 19 Unfortunately, these heuristics may be insufficient to achieve successful synthesis on the first attempt(s), necessitating follow-up experiments that can be time-intensive and costly. In the worst cases, human-biased heuristics lead to lower success rates than randomly-generated experimental protocols. 20 The a priori calculation of reaction selectivity in solid-state synthesis permits the ranking of synthesis approaches based on their thermodynamic likelihood of success, thereby circumventing the current time-consuming trial-and-error (Edisonian) approach. The assessment of reaction selectivity is particularly relevant in the proposal of optimal synthesis precursors; 21 in several cases, improved navigation of the phase diagram was shown to lead to a more practical synthesis. 11,[22][23][24][25] However, no solid-state reaction selectivity metric has been formally established. In recent work, 14 Aykol et al. demonstrated a computational workflow for ranking solid-state synthesis reactions by two performance metrics: 1) a catalytic nucleation barrier factor incorporating structural similarity and epitaxial matching, and 2) the number of known competing phases. These metrics perform well in rationalizing successful syntheses in the literature but lack generality; for example, the nucleation metric is derived assuming all reactions are nucleation-limited, which, as discussed, is not true for many solid-state reactions. Additionally, while a metric based on the total number of competing phases is significant as it hints at a measure of reaction selectivity, such a scheme does not account for the relative stability of these competing phases. A count-based selectivity metric is also biased by how many phases are known to exist at the present time and the extent to which various structural configurations (e.g., disordered or defective phases) have been enumerated within the data.
In this work, we address the longstanding issue of assessing the selectivity of solid-state reactions by deriving two complementary thermodynamic metrics measuring the degree of phase competition from the interface reaction model. We incorporate these competition metrics into a computational synthesis planning workflow for identifying and ranking synthesis reactions, which builds upon the high-throughput reaction enumeration tools we previously developed for constructing solid-state chemical reaction networks 26 from large materials databases such as the Materials Project. 27 Our selectivity metrics, computational workflow, literature analysis, and experimental findings yield a framework for generating more optimal and efficient solid-state synthesis routes, providing a foundation for the predictive synthesis of inorganic materials. The suggestion of nonstandard precursors, particularly those involving additional elements beyond those in the target composition, expands the synthetic capabilities of the solid-state approach.

Derivation of Selectivity Metrics
The Interface Reaction Hull To construct the interface reaction hull (Figure 1b), one begins with thermodynamic data for the reacting system (i.e., a set of relevant phases and their compositions and energies). In this work, we acquire formation enthalpies, ∆H f , from the Materials Project database 27 and extend them to Gibbs free energies of formation, ∆G f (T ), through the use of a prior machine learning model 28 and supplemental experimental thermochemistry data 29 (see Methods). For systems with two elements, the interface reaction hull is equivalent to the binary compositional phase diagram, where each vertex represents a single phase. However, for systems with three or more elements, the non-precursor vertices include both single phases and mixtures of phases. These mixtures are stoichiometric combinations of phases representing the products of balanced reactions of the precursors. The balanced reactions can be determined via 1) computing slices of the full compositional phase diagram along the tie-line connecting the precursors 12 or 2) combinatorial reaction enumeration. 26 More specifically, the maximum number of products for a particular vertex is one less than the number of elements in the system. The interface reaction hull is thus generalized such that all non-precursor vertices correspond to reactions with coordinates given by the atomic mixing ratio of precursors, x, and the Gibbs free energy of the reaction, ∆G rxn . This model can be further generalized to environmental conditions other than fixed temperature and pressure by constructing the hull with the appropriate thermodynamic potential. For example, in open systems, one would use the grand potential energy, Φ. Note that in these systems, the hull vertices may include additional open elemental reactants/products (e.g., O 2 ) that do not factor into the determination of x.
As a model for solid-state reactions, the interface reaction hull construction also rationalizes the formation of impurity phases. To demonstrate this, we revisit the binary system in Figure 1. In this system, the target phase (γ) is predicted to form first because it is the phase with the highest driving force of formation (most negative ∆G rxn ), irrespective of the average composition of the total system. 12 When the γ phase forms, however, it also introduces two additional interfaces: α|γ and γ|β. These secondary interfaces can be modeled via the construction of new interface reaction hulls or, more simply, by splitting the original hull into two subsections (i.e., to the left and the right of the target). Figure 1b suggests that the γ|β interface should produce an impurity phase via the exergonic reaction, γ + β −→ δ. Hence, the full conversion of reactants to the target phase is impeded while δ persists. In a "one-dimensional" solid-state reaction (Figure 1b, bottom), local thermodynamic equilibrium may be achieved when the system reaches a state in which all stable product phases on the interface reaction hull have formed and the growth of the product layer(s) slows down until it ceases entirely. This situation has been observed in previous experimental studies on diffusion couples. 30,31 This observed mixture of products may be kinetically "stable" (i.e., with locally stable interfaces), but it is not the global equilibrium state of the system. Rather, the global equilibrium state is the combination of phases that minimizes the free energy given the composition of the entire powder mixture; in Figure 1b, this corresponds to entirely γ. In powder reactions, access to the global equilibrium state is often provided by re-grinding steps in which new interfaces are exposed and mixed to facilitate complete conversion to the equilibrium products. However, in reacting systems with significant phase competition and slow transport, high temperatures and long heating times may be necessary but impractical; pure target synthesis may not be achievable if the desired products are unstable at high temperatures. These situations can be avoided entirely by proposing alternative precursors that are more selective (i.e., those with interface reaction hulls containing few to no competing phases).

Measuring Phase Competition
To predict the thermodynamic selectivity of a solid-state reaction, we propose two complementary metrics for assessing phase competition using the interface reaction hull: primary (C 1 ) and secondary (C 2 ) competition. Although both metrics measure the relative energetic favorability of competing reactions, they model different mechanisms for impurity formation.
The primary competition measures the favorability of competing reactions of the original precursors, while the secondary competition measures the favorability of subsequent competing reactions between the precursors and target phase(s). The origin of the two competition mechanisms is illustrated in Figure 2. In this simple interface reaction hull for a binary (two-element) system, two interface reactions can form a competing phase (magenta square). The primary reaction (yellow arrow) occurs at the interface between the two precursors α and β, whereas the secondary reaction (green arrow) occurs between the target phase (pink star) and the remaining β precursor, leading to a smaller driving force (arrow length). The coordinates of the competing phase, which must lie within the illustrated bounds (gray triangle) if the target phase is to be stable, determine the relative values of the primary and secondary reaction energies.
Primary competition, C 1 , is measured via calculation of the relative thermodynamic advantage of the most exergonic competing reaction from the original precursors, as assessed through an energy difference: (1) Here, ∆G rxn is the energy of the target synthesis reaction, and ∆G c i are the energies of possible competing reactions from the precursors. Lower C 1 values are favorable and result in more selective target formation. When C 1 is positive, the target reaction is less energetically favorable than the competing reaction with the greatest driving force (most negative energy), suggesting that a competing phase is likely to form. On the other hand, when C 1 is negative, the target reaction is predicted to have the greatest driving force of any reaction on the hull.
By considering only the single most competitive reaction, this functional form avoids the aforementioned bias related to using the total number of competing reactions. When no exergonic competing reactions are predicted for an interface, the competing reaction energy term is assigned a value of zero, representing the scenario in which the precursors do not react (e.g., α −−→ α). This results in the limiting condition: C 1 ≥ ∆G rxn .
Secondary competition, C 2 , assesses the favorability of impurity phase formation via secondary reactions between the target and precursor(s). This metric is important and distinct from primary competition because it measures the relative stability of the products of the target synthesis reaction with respect to decomposition into the competing phase(s).
Their relative stability can be measured by computing the "inverse distance to the hull" ( Figure 2), which for systems with one competing phase, is equivalent to the secondary reaction energy, ∆G d , at the precursor-target interface.
In an interface reaction hull with several competing phases, a sequence of multiple secondary reactions may occur ( Figure 3). When the target phase is formed, it introduces two new precursor-target interfaces that divide the hull into subsections to the left and right of the target. A secondary reaction may occur in either subsection, exposing another two interfaces. If, at either interface, there is a remaining driving force to form additional competing phases, then this process may continue in a recursive fashion until all possible secondary reactions have occurred. There are multiple ways to draw a feasible secondary reaction sequence (Figures 3b,c). Consider a particular secondary reaction sequence indexed j. This sequence has a total energy given by the sum of the energies of its n steps, where the number of reaction steps (n) in the sequence also equals the number of non-reactant (interior) vertices in the hull subsection.
If every secondary reaction step is required to be the one with the minimum energy (i.e., largest driving force), then only one unique reaction sequence exists in the left and right hull subsections. However, one must consider the situation where the minimum-energy principle does not hold due to kinetic limitations; this applies in particular to hulls where all secondary reactions have similar magnitude driving forces or a particular phase is kinetically limited from forming, perhaps due to an overall small driving force. Therefore, we must consider the alternative secondary reaction sequences shown in Figures 3b,c. These alternative sequences are not necessarily less favorable; while each alternative sequence may feature a first reaction step with a smaller driving force (i.e., small ∆G d j,1 ), the latter steps may have larger magnitude energies, resulting in comparable total energy (∆G 2,j ) for the particular sequence.
To encompass all combinatorial possibilities in our estimation of secondary competition, we choose to compute the mean total energy of all feasible secondary reaction sequences: Determining the total number of unique secondary reaction sequences, N , is mathematically equivalent to calculating the total number of full binary trees with n interior nodes, which yields the Catalan number sequence, u n = 1, 1, 2, 5, 14, 42, 132, 429, ...(n = 0, 1, 2, ...). 32 Using this connection to the Catalan numbers, we developed a non-recursive algorithm for calculating ∆G 2 that is significantly faster than the equivalent recursive solution (see Methods).
Finally, we formulate the C 2 metric such that it accounts for all possible reaction sequences in either hull subsection to the left (L) and right (R) of the target phase: The negative factor is included so that a lower C 2 value corresponds to a more favorable selectivity. Because our definition of a secondary reaction assumes that ∆G d ≤ 0, the secondary competition metric obeys the limiting behavior: C 2 ≥ 0.
We note that both C 1 and C 2 implicitly assume that the target phase is thermodynamically stable ("on the hull") under the conditions for which the equilibrium phase diagram is derived. However, the competition metrics are still calculable for a metastable phase by manually decreasing its energy until it becomes stable.

Application to Experimental Literature
Using the competition metrics C 1 and C 2 , we now assess the selectivities of solid-state re-   source of error in our reaction energy calculations is likely the disagreement between the assumed and actual gas partial pressures; it is often challenging to model the actual environmental conditions of synthesis, as it requires that they be both 1) accurately reported and 2) correctly extracted from the text. Another known source of error is systematic challenges in estimating GGA-calculated formation energies of carbonate compounds. 14 However, this has been addressed and partially mitigated with a fitted energy correction (see Methods).
To quantitatively assess thermodynamic optimality, we follow a similar approach to prior works 24,26 and define a cost function, Γ, that combines the driving force and reaction selectivity evaluated at a particular set of conditions (i.e., temperature and atmosphere). In this work, we opt to use a simple linear weighted summation of the reaction's energy and competition scores, where ∆E rxn is the reaction energy (either ∆G rxn or ∆Φ rxn ), and x 0 , x 1 , and x 2 are userdefined weights for each parameter. Due to the different scaling of each parameter, we find that a non-equal weighting of x 0 = 0.10, x 1 = 0.45, and x 2 = 0.45 produces reasonably diverse results that do not significantly favor one parameter over another. We note that this selection is arbitrary and subject to further optimization.
Unfortunately, closed and open reactions cannot be rigorously compared due to their different energy scales (i.e., Gibbs free energies vs. grand potential energies). This is further evidenced by differences in the reaction metric distributions (Figures 4a,b).  from the text). Given the extent of coverage of our thermodynamic data, we limit our analysis to only the targets for which we have at least five recipes successfully calculated (see Methods). In the following sections, we select several targets with costs below and above the median Γ t value, analyzing the factors leading to their optimal and suboptimal synthesis recipes, respectively.

Optimal Literature Recipes
Of the 40 most popular targets in the literature dataset ( Figure 4c), twenty have synthesis recipes with an average cost value (Γ t ) below the dataset's median (-0.089), indicating generally favorable thermodynamic optimality. A selection of these targets and their highest/lowest-ranked recipes are provided in Table 1, along with other selected reactions of interest. DOIs and raw (untransformed) costs for each reaction are provided in the Supporting Information.
For many targets in Table 1, the conventional reaction involving off-the-shelf binary precursors is predicted to be thermodynamically optimal. For example, in the synthesis of the spinel ZnGa 2 O 4 , the reaction between the binary oxides, ZnO is already perfectly selective (C 1 = ∆E rxn and C 2 = 0) over all temperatures in the dataset.
Furthermore, these precursors appear in all 17 calculated literature recipes for this target. in an open oxygen environment favor the production of BiVO 4 and BiFeO 3 . Finally, for Y 3 Al 5 O 12 , there appears to be a tradeoff between selectivity (which is optimal at lower temperatures) and driving force (which is optimal at high temperatures); it appears that intermediate temperatures (600 ℃) in a closed environment result in the most suitable compromise. We note, however, that there are many reasons to use specific processing conditions outside of the pursuit of target phase purity (e.g., improved density, annihilation of defects, optimization of crystallite sizes, etc.). These alternative reasons may explain some of the variability of conditions reported in the literature.
Interestingly, several targets in Table 1 feature synthesis recipes that appear extremely unfavorable. These often involve elemental precursors, such as Bi + 0.
Referencing the original articles from which these reactions were sourced 34,35 suggests that these precursors were not actually used, and LiOH (particularly the monohydrate) also appears to offer some advantage in the synthesis of LiNiO 2 , especially in oxygen at low temperatures (480-600 ℃).

Suboptimal Literature Recipes
Reactions can still be successful despite high C 1 and C 2 values. However, our thermodynamic assessment suggests that these reactions are suboptimal and likely require some combination Table 1: Thermodynamic analysis of optimal experimental synthesis recipes for selected popular targets in the literature. The ranking of reactions is determined by the transformed cost, Γ t . The highest-and lowest-ranked reactions are shown along with selected reactions of interest. Raw costs and DOIs are available in the Supporting Information.

Reaction
Temp. Energy Open 412  of 1) long heating times to promote thermodynamic equilibrium, 2) follow-up regrinding steps, or 3) fine-tuning of temperature and reaction atmosphere. In Table 2, we highlight several popular target materials that are associated with higher-than-average costs for their synthesis recipes.
Our findings suggest that the targets in Table 2 require a more judicious synthesis due to inherently greater competition in their phase space. Many targets are ranked poorly because the conventional reaction is predicted to be suboptimal. This appears to be true for   36 The use of PbCO 3 as an alternative appears to perform similarly, albeit with a slightly more favorable driving force and lower C 1 .
For Ca 3 (PO 4 ) 2 , Li 2 MnO 3 , and LiFePO 4 , the highest-ranked synthesis recipes appear to be nearly thermodynamically optimal already (i.e., low C 1 and C 2 ). The optimal recipe for and open oxygen conditions appears to be the most favorable; however, with this route, the driving force is still somewhat low and features a moderate amount of phase competition. The BaTiO 3 system will be examined extensively in the following sections as an experimental case study for assessing reaction selectivity.
The literature reaction data reveal an inherent tradeoff between driving force and selectivity. Reactions with only elemental precursors (e.g., Li metal, O 2 gas) typically have significantly greater driving forces but also much higher competition scores (i.e., lower se- Determining an optimal synthesis reaction can be formulated as an optimization problem of simultaneous minimization of the three reaction metrics: ∆G rxn , C 1 , and C 2 . A common approach for multi-objective optimization in synthesis planning is identifying the Pareto front. 14 Here, we calculate a three-dimensional Pareto front for the BaTiO 3 synthesis reactions ( Table 3)    byproducts (e.g., Ba 2 TiO 4 ), or refractory precursors (e.g., TiC). Some of these suggested reactions can be removed easily with user-applied filters that account for specific experi-mental restrictions. For example, one can supply a list of available precursor compositions (i.e., "off-the-shelf" phases) or composition types to be avoided (e.g., sulfides, acids, etc.).
In our provided code (see Methods), we support functionality for the former by including a list of hundreds of common precursors compiled from the catalogs of chemical suppliers.
Filtering by these commonly available precursors (e.g., be screened for reactivity, volatility, safety, and material costs. These challenges can be mitigated through the use of additional data or models; for example, reactivity can be approximated through surrogate data, such as defect formation energies or physical properties (melting points, hardness, etc.). 42 We selected nine reactions from Figure 5 to test experimentally. The calculated thermodynamic metrics for these reactions are provided in Table 4. We intentionally selected reactions spanning various precursor chemistries, free energies, and competition scores. To prioritize the study of impurity-forming reactions and avoid the aforementioned practicality challenges, we did not explicitly include any reactions on the Pareto front. The conventional synthesis route, BaCO 3 +TiO 2 −−→ BaTiO 3 +CO 2 , was chosen as a baseline reference (Expt. 1). A positive energy is calculated for this reaction (∆G rxn =+0.042 eV/atom at T = 600 ℃); however, this is likely due to residual uncorrected error in the calculated energy of BaCO 3 , which is not included in the NIST-JANAF dataset. We did not test the analogous reaction with BaO precursor due to its hygroscopic nature, which makes it difficult to handle.
However, we did test the alternative reaction from barium peroxide, BaO 2 (Expt. 2). We included two reactions that form BaTiO 3 directly from at least one other ternary phase (Ex-pts. 3,4). A reaction with a Ti metal precursor was selected due to its extremely high C 2 score (Expt. 5). Two metathesis reactions (Expts. 6, 7) were selected due to their predicted high performance, including the unconventional use of a sulfide precursor (BaS). The final two reactions (Expts. 8,9) were selected for being endergonic (∆G rxn > 0) to validate the accuracy of our free energy predictions. We note that several experiments feature precursors that are not easily purchasable from a chemical supplier (e.g., Ba 2 TiO 4 , Na 2 TiO 3 ); these phases were synthesized following recipes reported in the literature (see Methods). The nine synthesis experiments were completed using a gradient furnace (see Methods), 43 allowing the observation of reaction products over a wide range of temperatures (∼200-1000 ℃) with ex post facto SPXRD. The experimental results are summarized in Figure 6.
Mole fractions of each phase were determined via Rietveld refinement of SPXRD patterns acquired at various positions (temperatures) along the length of the sample after heating and subsequent cooling to ambient temperature. The ex post facto phase fraction plots can be interpreted similarly to those constructed with in situ data, as each effectively illustrates the reaction pathway during heating. More rigorously, however, each point in Figure 6 is pseudoindependent and best interpreted as the result of an isothermal (ex situ) reaction at its associated temperature. Shorter reaction times were selected to ensure capture of the onset of short-lived intermediate phases critical to assessing reaction pathway selectivities (Table  S2). Selected Rietveld analysis results for each experiment are provided in the Supporting Information ( Figures S4-S12). (ai) Mole fractions of observed phases for Experiments 1-9, as determined through Rietveld refinements of ex post facto SPXRD data. Phase types are distinguished by shape: precursors (circles), targets (diamonds), byproducts (triangles), and impurities (exes). Background shading denotes the total mole fractions of precursor (gray), impurity (pink), and target/byproduct (green). The median total reaction time was ∼67 min; exact times for each experiment are provided in Table S2.
The observed reaction pathways demonstrate significant variation in target and impurity formation. Visualizing the interface reaction hulls of the selected experiments helps to rationalize their predicted and observed performance ( Figure S13). The most complex of these pathways is that of Expt. 5, which features the formation of BaTiO 3 at an intermediate temperature range (400-500 ℃) before impurities Ba 2 TiO 4 , TiH 2 , and TiO 2 begin to dominate. Indeed, Expt. 5 exhibits both the largest driving force (-0.530 eV/atom) and C 2 value (0.534 eV/atom) of any reaction, supporting the observation of a complex reaction pathway containing many impurity phases.
The conventional synthesis reaction between BaCO 3 and TiO 2 (Expt. 1) was largely incomplete after 60 min at 1000 ℃ and exhibited significant formation of BaTi 2 O 5 . The interface reaction hull for this system ( Figure S13a) suggests that the formation of BaTi 2 O 5 is the most favorable reaction outcome, supporting our observations. In Expt. 2, the reaction of TiO 2 with BaO 2 appears to significantly decrease the reaction onset temperature but also features substantial impurity formation. In three of the experiments (Expts. 4, 6, and 7), the BaTiO 3 synthesis reaction is the most favorable reaction on the hull, resulting in C 1 < 0.
Notably, the use of ternary precursor(s) in Expts. 3 and 4 (Ba 2 TiO 4 , BaTi 2 O 5 ) results in low (or zero) C 2 , but also very little driving force. As a result, we observe near-perfect selectivity  the calculated reaction metrics (∆G rxn , C 1 , and C 2 ). The full set of all (3x3) pairwise correlation plots is available in Figure S14. We observe that the calculated reaction energy (∆G rxn ) correlates most strongly with the minimum amount of precursor remaining (P ) at the conclusion of the experiment (Figure 7b). With infinite reaction times, we would theoretically expect this distribution to resemble a step function: P = 0 for reactions with ∆G rxn < 0, and P = 1 for those with ∆G rxn > 0. In our work, the distribution is less defined, given the shorter reaction times and different chemistries explored. The coordinates of Expts. 1 and 2 appear to deviate the most from a step-like distribution. Difficulty in modeling the energetics of carbonate reactions was previously discussed and likely explains the deviation of Expt. 1. In contrast, the deviation of Expt. 2 is likely kinetic in nature, as the reaction appears to stall ( Figure 6b); this may suggest that the maximum temperature (750-800 ℃) is too low to achieve sufficient reaction completion.
The primary (C 1 ) and secondary (C 2 ) competition metrics display negative and positive correlations with the maximum amounts of target (T ) and impurity (I) phases formed, respectively (Figures 7c,d). We note that these trends are not strictly obeyed in a monotonic fashion. Still, the correlations are significant and can be theoretically rationalized by their derivation from the interface reaction hull. It is reasonable that C 1 correlates most strongly with T , as this competition metric effectively measures the relative favorability of the target reaction over competing reactions. Similarly, it is sensible that C 2 should correlate most strongly with I given its derivation as a measure of the relative stability of impurity phases with respect to the target phase. To be precise, the experimental correlation between C 2 and I suggests that one should consider not only the sum of the inverse hull distance energies for all competing phases but the total energy of the entire secondary reaction sequence containing them (Equation 2). By definition, this quantity includes, and thus will always be larger than or equal to, the sum of all competing inverse hull distance values for a particular interface reaction hull. Regarding the functional form of C 2 , the question arises as to whether a more simple summation of the maximum energy secondary reactions in the left and right hull subsections is a sufficient measure for secondary competition. While the maximum energy secondary reactions tend to account for much of the value of C 2 , on average, the full C 2 metric is 0.077 eV/atom greater than considering the maximum energy secondary reactions alone ( Figure S15), suggesting that C 2 is a more conservative metric. For completeness, we also tested an alternative formulation of the secondary competition metric using the enclosed "area" to the hull. While this metric correlates with C 2 , its calculation is more numerically unstable, and its units are less interpretable. Hence, we generally recommend the approach of modeling full secondary reaction sequences, which is straightforward to implement using our secondary competition algorithm (see Methods).
In a previous study, 24 we suggested a solid-state reaction selectivity metric based on the difference in elemental chemical potentials between precursors and targets, measured by a distance along the chemical potential diagram (i.e., the "total chemical potential distance").
This metric was used to rationalize the unique selectivity of the Na-based precursor in synthesizing pyrochlore Y 2 Mn 2 O 7 from YOCl and AMnO 2 (A=Li, Na, K). While straightforward to compute using just the chemical potential diagram, the distance metric operates in the space of chemical potentials rather than reaction energies, rendering it less intuitive and more difficult to precisely discern the specific competing reactions. From the results here, we generally recommend that C 2 be used instead of total chemical potential distance where possible. Both selectivity metrics capture similar characteristics of the competing phase space: each effectively involves the summation of competing phase stabilities through inverse hull distances or the corresponding chemical potential stability ranges. More precisely, both quantities are correlated ( Figure S16) because chemical potentials are mathematical derivatives of the convex hull in energy-composition space. However, the total chemical potential distance is biased, particularly by competing phases with defective elemental-like compositions (e.g., Mg 149 Cl). Due to the increased weight of the entropic (−T S) term in the definition of Gibbs free energy, chemical potential diagrams featuring these compositions as competing phases may yield very high (unfavorable) total chemical potential distance values for synthesis reactions.
We acknowledge that while C 1 and C 2 are meant to capture different, independent mechanisms by which competing phases form, these metrics are at least partially correlated due to the geometric constraints of the convex hull. In particular, one situation is geometrically limited from occurring: high C 2 and low C 1 ( Figure S17). Stated explicitly: if a competing phase lies significantly below the tie-line formed by the target and a precursor (i.e., high C 2 ), then both the target and that competing phase necessarily have similar reaction energies to form from the precursors, leading to high C 1 . In general, however, this restriction does not make the metrics redundant; while there is some correlation between the two, the correlation is not particularly strong ( Figure S16). Therefore, we generally recommend the tandem use of both selectivity metrics.
The major limitation of our current synthesis planning workflow is the assumption that optimal synthesis reactions can be predicted with the thermodynamic energy landscape alone. While this is not the case for all chemistries, we show that, at least for chemical systems that exhibit practical solid-state reaction kinetics, the energy landscape alone can provide much of the rationalization for the observation of impurity phases. To say that impurity and secondary phases are inherently "kinetic" products is a misnomer. Rather, these phases may be the thermodynamic minima of smaller "local" interface systems, distinct from the thermodynamic products of the entire reaction mixture (the global thermodynamic solution). Furthermore, these impurities are often not easily convertible to final products without long-range mass transport or intervention (e.g., via re-grinding or subsequent heating). This explains why impurities are often pervasive and challenging to remove in chemical systems with lower driving forces and/or slower kinetics (e.g., BiFeO 3 ).
Although our current study focuses on the synthesis of oxides, we expect our synthesis planning approach to be suitable to other chemistries where solid-state synthesis can be employed. This includes the chemistries of most ionic compounds: halides, chalcogenides, pnictides, some silicides/carbides, etc. Still, one must ensure that there is enough thermodynamic data available to accurately model phase competition in the chemical system of interest. This is generally true for oxide compounds due to their high prevalence in literature and thermodynamic data; for example, currently, ∼53% of the nearly one-hundred and fifty-thousand compounds in the Materials Project contain oxygen. While the predictive accuracy is currently greatest for oxides, we expect our approach to grow in accuracy and general applicability as computed materials databases grow in size and chemical complexity.
Currently, our workflow focuses exclusively on optimizing product purity; however, there are many issues one must consider when designing a synthesis recipe for a target compound: material cost, safety concerns, stability in air, handling challenges, availability of precursors, etc. Many routes suggested involve the formation of byproducts that are not easily removable from the product mixture (e.g., the formation of BaTiO 3 with byproduct Ba 2 TiO 4 ). To this point, one should be thoughtful in designing criteria by which to filter recommended synthetic routes. For example, one can prioritize the formation of only gaseous byproducts (e.g., O 2 or CO 2 ) or those easily removable by a solvent (e.g., NaCl). The cost function used to rank reactions can be modified to include other reaction metrics of interest, such as the estimated economic cost of the precursor materials. While not explicitly demonstrated here, the synthesis planning workflow can also be extended for application in multi-step syntheses, allowing one to retrosynthetically sequence reactions to a target material beginning with purchasable, "off-the-shelf" precursors.

Conclusions
Using the interface reaction model for powder reactions, we proposed two thermodynamic selectivity metrics for solid-state reactions: primary (C 1 ) and secondary (C 2 ) competition.
To systematically and critically examine the effectiveness of our metrics, we analyzed ex- X-ray diffraction reveals that C 1 and C 2 correlate with the maximum amounts of target and impurity formed, respectively.
The main advantage of our approach compared to recent, existing approaches 14,25 is the ability to simultaneously consider a wide range of chemistries, including those with unconventional additional elements. These so-called hyperdimensional chemistries 40 allow one to bypass commonly encountered intermediates in target systems with many competing phases. This was demonstrated particularly for the BaTiO 3 system studied in this work and is relevant for many other materials in the literature that are conventionally synthesized with theoretically suboptimal precursors (e.g., Na 2 Ti 3 O 7 , NaTaO 3 , LiMn 2 O 4 , etc.).
We anticipate that the selectivity metrics presented here and our computational synthesis planning workflow will significantly reduce the synthesis bottleneck, providing more rapid development of synthesis approaches for new, predicted materials. Our workflow provides a theoretical rationale for using certain precursors and synthesis conditions over other options, which promises to optimize existing synthesis procedures for current technologically important materials.
We envision our approach to be particularly useful in aiding high-throughput automated laboratory exploration efforts. 45 Predictions can be used to design and downselect the synthesis reactions tested, reducing the cost and current trial-and-error approach to inorganic materials synthesis. The future inclusion of models for the kinetic behavior of reactions, such as estimates of the reactivity of precursors based on solid-state diffusivities, will further enhance predictions.

Thermodynamic Data
Gibbs free energies of formation, ∆G f (T ), were acquired or approximated in a similar approach to those of previous works. 24,26 We acquired experimental ∆G f (T ) values from the NIST-JANAF thermochemical tables 29 where available. Experimental values were limited to compounds with low melting points (i.e., T m ≤ 1500 ℃), as these systems demonstrate more complex phase change behavior over the temperature range studied here. For predominantly solid compounds (i.e., those with melting points above this threshold), as well as for all other phases not available in the NIST-JANAF thermochemical tables, we estimated ∆G f (T ) using the machine-learned Gibbs free energy descriptor identified by Bartel et al. 28 This descriptor was applied using formation enthalpies, ∆H f (T = 298 K), acquired from the  Figure S18).

Synthesis Planning Workflow
The synthesis reaction calculation and ranking procedure was implemented as a Pythonbased workflow in the existing reaction-network package. 26 The code is available on GitHub at https://github.com/materialsproject/reaction-network. The workflow was constructed and launched on computing resources using the jobflow 46 and fireworks 47 workflow packages.
The synthesis planning workflow consists of three sequential steps. First, phases and their formation energies for the chemical system of interest are acquired as previously described.
The total number of phases can be optionally reduced by setting a threshold for the maximum energy above hull (∆G hull ). In this work, we used a moderately large threshold of ∆G hull ≤ 50 meV/atom, evaluated at ambient temperature (T = 300 K). Second, reaction enumeration is performed for the acquired phases using the combinatorial and free energy minimization approaches described in our previous work on solid-state reaction networks. 26 Note that the combinatorial approach allows one to identify reaction product combinations above the hull (i.e., "metastable" products), which makes the analysis more robust to numerical error in the thermodynamic data.

Secondary Competition Algorithm
The secondary competition score, C 2 , is defined as the negative sum of the mean secondary reaction sequence energies to the left and right of the target on the interface reaction hull (Equation 4). One approach for acquiring these quantities involves using a recursive algorithm to identify all possible sequences and their energies. However, this strategy is too slow for the high-throughput calculation of C 2 in systems with many competing reactions.
Instead, we have identified a non-recursive algorithm that takes advantage of the connection between this problem and the recursive construction of binary trees via the use of the Catalan number sequence. Our algorithm reformulates the sum of all secondary reaction sequence energies as a sum of individual secondary reaction energies weighted by their multiplicities, i.e., the total number of appearances of a particular reaction within the set of all possible secondary reaction sequences. The energy of any reaction indexed k can be calculated geometrically as the altitude, h k , of the triangle formed by its product vertex and two reactant vertices on the interface reaction hull. We find that the altitude multiplicity, m h k , is determined to be the product of three Catalan numbers, u n , such that where n l and n r refer to the number of interior vertices (i.e., within the triangle) to the left and right of the vertex of interest, respectively, and n is the total number of interior vertices for the entire hull subsection. For example, for secondary reactions between nearest neighbors, n l = 0 and n r = 0, resulting in an altitude multiplicity of m h = u n−1 .
The mean secondary reaction sequence energy for the hull subsection can then be calculated as where the sum occurs over all of the V unique reaction energies (altitudes), which is the number of unique triangles that can be constructed for the hull subsection, including the two exterior vertices: V = n+2 3 . The total number of unique secondary reaction sequences equals the corresponding Catalan number, N = u n . Finally, once this process has been performed for both the left and right hull subsections, the secondary competition (C 2 ) can be calculated via Equation 4.

Literature Reactions
Solid-state literature reactions studied in this work were acquired from the text-mined dataset of 31,782 inorganic materials synthesis recipes originally extracted from the literature by Kononova et al. 33 and available at https://github.com/CederGroupHub/text-mined-s ynthesis_public (version 2020-07-13). The original dataset was filtered down to 8,530 reactions that contain: 1) precursors composed of ≤ 2 solids and ≤ 1 elemental gases (i.e., O 2 , H 2 , and N 2 ), 2) no elements with an atomic number greater than 94 or for which the Gibbs free energy descriptor does not apply (e.g., Ne, Ar, Pm, Ra), 3) ten or fewer total elements due to limitations in the convex hull algorithm, and 4) no ions. Finally, these reactions were required to be stoichiometrically balanceable after adjusting compositions for hydrates and fractional formulas. For reactions containing variable compositions with one open variable (e.g., Nd 1-x Sr x CoO 3 ), we attempted to substitute all extracted values of x and retained the reactions that could be successfully balanced.
Competition metrics and free energies were assessed for each of the remaining reactions.
For the enumerated competing reactions, metastable phases were considered up to a maximum threshold of ∆G hull = 50 meV/atom, evaluated at ambient temperature (T = 300 K). Interface reaction hulls were constructed at the maximum temperature reported during synthesis, T syn . If this was not provided, a temperature of 800 ℃ was assumed. Formation energies, ∆G f (T syn ), were assigned based on the ground-state energy for a given composition; i.e., we selected the lowest available formation energy of all polymorphs with the composition of interest. For increased accuracy, we did not include a reaction if any of its entries were missing from our thermodynamic data. For reactions with an open gas (i.e., O 2 , H 2 , N 2 ), we assigned a chemical potential of µ gas = 0 eV (i.e., standard state at T syn ) for that element.
For reactions completed in air, we assumed an O 2 partial pressure of 0.21 atm and thus assigned a chemical potential of µ O = 1 2 k b T syn ln(0.21) eV. Finally, we removed duplicates with the same reaction equation and temperature/environment, as well as identity reactions (e.g., A → A). These filtering steps yielded a total of 3,520 unique literature reactions.

Precursor Materials
Precursors for all experiments were purchased from chemical providers or prepared via known  BaS was prepared using BaSO 4 and activated carbon (C, J.T.Baker 99.9%). 50 The chemicals were mixed, ground using a mortar and pestle, pressed into a 0.5" diameter pellet with two tons of force, and heated in an alumina boat at 1100 ℃ for 7-10 min in air, with a heating rate of 10 ℃/min and natural cooling in the furnace. The product was phase pure with no detectable impurities. Na 2 TiO 3 was prepared using stoichiometric amounts of sodium hydroxide (NaOH, Fisher Scientific 99.9%) and anatase TiO 2 , with a slight excess of NaOH. 51 The chemicals were mixed, ground using a mortar and pestle, and heated in an alumina boat at 500 ℃ for 2 hrs with a heating rate of 10 ℃/min. The product was mostly phase pure with minor impurities. The sodium titanate peaks are best fit by a cubic α-Na 2 TiO 3 structure with a small crystallite size. A minor amount of unreacted anatase TiO 2 was present in the product (∼1 mol %). Na 2 CO 3 also appears to be present as an impurity (∼11 mol %); we suspect this is due to contamination of the NaOH precursor via reaction with CO 2 in the air.

Ex Post Facto SPXRD Reactions
Synchrotron powder X-ray diffraction (SPXRD) data were collected in transmission (i.e., Debye-Scherrer) geometry on beamline 28-ID-2 (XPD) at the National Synchrotron Light Source-II (NSLS-II). Data were collected on a 2D area detector (PerkinElmer XRD 1621, 2048x2048 pixel array, 200x200 µm pixel size) at a sample-to-detector distance of 1407.1 mm using an incident X-ray energy of 68.12 keV (λ = 0.182Å) with a 0.60 x 0.20 mm beam size.
A total acquisition time of 1 s was used, summing five subframes collected for 0.2 s each.
Samples were packed into 1.1 mm OD/0.9 mm ID quartz capillaries. The capillary ends were filled with a 3 mm plug of powder silicon (Si, Strem 99.0%) followed by a cap of recycled silicon dioxide (SiO 2 ). To account for possible gas production, the capillaries for Expts. 1, 2, 5, 7, and 8 were left unsealed, and a moderate vacuum (P gage = −20 in Hg) was pulled on the samples during heating. All other sample capillaries (Expts. 3, 4, 6, and 9) were flame-sealed under argon.
Experiments were carried out in the gradient furnace described in Ref. 43, which heats samples to different temperatures across a range of spatial positions on the capillary ( Figure   S19). The furnace was operated with a Eurotherm 2408 temperature controller and a TDK Lambda 900W (30V/30A) power supply. Furnace heating elements were wound from resistive wire (Kanthal A-1, #24 awg). A K-type thermocouple (stainless steel, 0.01" OD) placed at an intermediate position along the sample was used as input for PID control of the furnace. We performed experiments in three temperature ranges with setpoints of T H = 550 ℃ (Expts. 1,3,4,8), T L1 = 450 ℃ (Expts. 2,5,6,9), and T L2 = 400 ℃ (Expt. 7). This choice was motivated by differences in reactivity among the samples. Position-dependent temperatures were determined using a fit of measured in situ lattice expansion from NaCl/Si and Al 2 O 3 /MgO standards ( Figure S20). Using the root-mean-square error of the curve fit, the estimated uncertainty for each temperature point is 11.4 ℃ (T H ), 7.9 ℃ (T L1 ), and 10.9 ℃ (T L2 ). The experiments spanned a total temperature range of 189-1064 ℃.
The median total time of each experiment was ∼67 min. Heating, holding, and cooling times varied among experiments due to differences in sample heat capacities, reactivities, and the sizes of investigated temperature windows; specifically, unreactive samples (Expts. 8,9) and samples with smaller studied temperature windows (Expts. 3,4) were held at elevated temperatures for shorter times. These differences are accounted for in our analysis via the use of relative metrics (i.e., mole fraction) and normalization by reaction progress.
The exact heating, holding, and cooling times for each sample are shown in Table S2.

Quantitative Phase Analysis of Powder Diffraction Data
Quantitative analysis of synchrotron powder X-ray powder diffraction data was carried out with the Rietveld method using either the TOPAS v6 (Expts. 1-3, 5, 7-9) or GSAS-II (Expts. 4, 6) software packages. 52,53 Atomic displacement parameters were fixed to B = 1Å 2 , and peak broadening was primarily modeled via crystal size broadening using a Lorentzian function. Site occupancies were fixed at one.

TOC Graphic
Synopsis Selectivity analysis of massive inorganic chemical reaction networks identifies alternative synthesis recipes for solid-state materials (e.g., BaTiO 3 ) that outperform conventional approaches.
Power transform of costs in literature reactions Figure S1: Power transformation applied to literature reaction costs, Γ. The power transform monotonically transforms the reaction costs for closed and open reactions such that they resemble standard normal distributions, allowing for better comparison between the two datasets.

S-2
BaO|TiO 2 interface reaction hull Figure S2: Calculated interface reaction hull for BaO|TiO 2 at T = 600 ℃. Blue squares denote experimentally known ternary phases in the system. Ba 2 TiO 4 is the most commonly observed intermediate/impurity in standard syntheses of BaTiO 3 from BaO (or BaCO 3 ) and TiO 2 ; it is also predicted to be the phase with the greatest driving force to form. Table S1: Selected experimental reactions to BaTiO 3 and their associated grand potential energies (µ O = −0.1001 eV), ∆Φ rxn (T = 600 ℃), primary competition scores, C 1 , secondary competition scores, C 2 , and costs, Γ.
(eV/at) (eV/at) (eV/at) (eV/at)   Figure 5, reactions are plotted on a shared axis of reaction energy, ∆Φ rxn , and on independent axes of (a) primary competition, C 1 , and (b) secondary competition, C 2 . Orange squares represent reactions on the three-dimensional Pareto frontier of ∆Φ rxn , C 1 , and C 2 . Blue diamonds indicate selected reactions experimentally tested in this work.

S-5
Selected Rietveld refinements for BaTiO 3 experiments Figure S4: Selected Rietveld refinement for Experiment 1. The observed pattern represents ex post facto synchrotron powder X-ray diffraction data (SPXRD) captured following reaction at T = 1040 ℃, which corresponds to the temperature with the highest BaTiO 3 yield. Phase fractions are shown in units of mole percent.

S-11
Figure S10: Selected Rietveld refinement for Experiment 7. The observed pattern represents ex post facto SPXRD data captured following reaction at T = 727 ℃, which corresponds to the temperature with the highest BaTiO 3 yield. Phase fractions are shown in units of mole percent.

S-12
Figure S11: Selected Rietveld refinement for Experiment 8. The observed pattern represents ex post facto SPXRD data captured following reaction at T = 603 ℃. This value corresponds to the median temperature studied since BaTiO 3 was not observed at any temperature. Phase fractions are shown in units of mole percent.

S-13
Figure S12: Selected Rietveld refinement for Experiment 9. The observed pattern represents ex post facto SPXRD data captured following reaction at T = 594 ℃. This value corresponds to the median temperature studied since BaTiO 3 was not observed at any temperature. Phase fractions are shown in units of mole percent.

S-14
Interface reaction hulls for selected BaTiO 3 experiments Figure S13: Interface reaction hulls for selected BaTiO 3 experiments. (a-i) Hulls for Experiments 1-9, as extracted from the full reaction network calculated during the synthesis planning workflow. Red diamonds mark the selected reactions of interest. All hulls are plotted on a uniform energy scale to facilitate comparison.

S-15
All pairwise correlation plots between reaction metrics and experimental outcomes Figure S14: Pairwise correlation plots between reaction metrics and experimental outcomes. Plots of minimum remaining precursor (P ), maximum target (T ), and maximum impurity (I) formed as a function of reaction energy, primary competition, and secondary competition. The target and impurity plots have been normalized by the amount of precursor consumed (1 − P ).

Correlations between selectivity metrics
Figure S15: Pairwise correlations of alternative secondary competition metrics. Secondary competition (max) is defined as the sum of the energies of only the secondary reactions with the highest driving forces on either side of the target. Secondary competition (area) is the enclosed area of the interface reaction hull. Figure S16: Pairwise correlations of selectivity metrics. The total chemical potential distance is calculated using the methodology outlined in Ref. S1.

S-18
Figure S17: Examples of varying C 1 and C 2 indicating their partial correlation.
(a, b) When secondary competition (C 2 ) is low, primary competition (C 1 ) can either be low or high. c) However, when C 1 is low, C 2 can not be high due to the geometric constraints of the hull. d) There are no restrictions for both C 1 and C 2 to be high. The pictured device has the same specifications as described in Ref. S2. The heating wires are wound with variable pitch, resulting in a wide temperature profile over the powder sample (capillary tube). A thermocouple is used as a setpoint for determining the power supplied to the heating elements.