Autonomous materials synthesis via hierarchical active learning of nonequilibrium phase diagrams.

Autonomous experimentation enabled by artificial intelligence offers a new paradigm for accelerating scientific discovery. Nonequilibrium materials synthesis is emblematic of complex, resource-intensive experimentation whose acceleration would be a watershed for materials discovery. We demonstrate accelerated exploration of metastable materials through hierarchical autonomous experimentation governed by the Scientific Autonomous Reasoning Agent (SARA). SARA integrates robotic materials synthesis using lateral gradient laser spike annealing and optical characterization along with a hierarchy of AI methods to map out processing phase diagrams. Efficient exploration of the multidimensional parameter space is achieved with nested active learning cycles built upon advanced machine learning models that incorporate the underlying physics of the experiments and end-to-end uncertainty quantification. We demonstrate SARA’s performance by autonomously mapping synthesis phase boundaries for the Bi2O3 system, leading to orders-of-magnitude acceleration in the establishment of a synthesis phase diagram that includes conditions for stabilizing δ-Bi2O3 at room temperature, a critical development for electrochemical technologies.

INTRODUCTION

Artificial intelligence (AI) holds great promise for revolutionizing scientific fields as varied as biology (), chemistry (), physics (), and economics (). Much of AI’s impressive recent successes have been in data-rich applications: AlphaFold (), for example, uses a library of tens of thousands of existing protein data to build a highly successful model for protein folding. In fields that lack comparably vast data, however, a great opportunity lies in guiding the exploratory process itself to minimize the number of experiments that are required to achieve insights, i.e., active learning (AL) (), and thereby accelerate the pace of scientific discovery ().

These AI-guided efforts have shown great promise in the design of quantum experiments (), drug development (), and wind turbine control () and are of particular importance in materials research aimed at designing and optimizing functional materials that lie at the core of technological advances. High-throughput (HT) experimental synthesis and characterization of materials systems through thin-film deposition of inorganic composition spreads (), so-called libraries, present a promising avenue to rapidly explore a vast chemical, structural, and property space (). These methods have been well established for comprehensive synthesis of composition spaces with two to four components, where the resulting tens to thousands of materials can be evaluated via automated characterization. While this approach has been quite effective for identification of materials with desired properties, the opportunity for broader materials exploration to enable new technologies is highlighted by the limited exploration of synthesis conditions to date. The portion of the materials search space that has been explored is vanishingly small when considering the dynamic range of thermal processing conditions, which are inherent to processing-composition-structure-property (PCSP) relationships. The breadth of relevant thermal processing conditions makes exhaustive sampling untenable, and as HT experiments lower the time for an experiment cycle, human intervention in decision-making becomes a worsening bottleneck. Therefore, AI and AL are critical in both reducing the number of experiments to a more tractable scale and accelerating the process of making decisions to match the rate of incoming data (–).

Recently, HT experimentation and AL techniques have been combined in a closed-loop fashion, where an AI instance iteratively proposes a sequence of experiments to explore and discover new materials. These efforts include identifying phase-change materials via Bayesian AL (), the discovery of NiTi-based shape memory alloys with low thermal hysteresis (), the synthesis of BaTiO3-based piezoelectrics with large electrostrain (), the selective growth of carbon nanotubes (,), the search for perovskite-type materials for photovoltaic applications () and inorganic quantum dots (), maximizing hole mobility of organic solar cells (), and accelerating toughness optimization in additive manufacturing ().

Despite this progress, the current state of the art exhibits considerable limitations. Chiefly, most attempts at closed-loop cycles still rely heavily on human intervention, preventing them from reaching true autonomy in materials discovery. Furthermore, although AL guidance has recently been deployed to great effect for discovery of optical phase-change thin-film materials (), the search space was limited to presynthesized compositions using a single processing condition. The time and temperature scales relevant to nonequilibrium thermal processing of solid-state inorganic materials pose substantial problems for incorporating synthesis in an autonomous loop, although the utility of spanning synthesis and characterization in an AL framework has been demonstrated by several recent autonomous workflows for chemical synthesis (, –). The complexity and the degrees of freedom of the PCSP space are particularly challenging to incorporate in autonomous experimentation when considering metastable materials that form far from equilibrium at different, often unpredictable processing conditions (,). Commonly deployed off-the-shelf AL models are often not sufficient for achieving highly efficient learning and are frequently outperformed by random search with twice the number of samples (), a problem that is exacerbated with increased dimensionality of the search space ().

Expert human scientists navigate complex search spaces by incorporating their prior knowledge, such as physics-based models that underlie the acquired data. Incorporating such knowledge in AL often requires the development of new AI methods. Lastly, exploration via AL critically relies on uncertainty quantification in the not-yet-sampled regions of parameter space, which, for complex experimental workflows, requires error propagation. Arguably, the most immediate obstacle to accelerating experimental exploration via AL lies in the dual challenges of developing noise models for each type of experiment and integrating them into a computational framework for end-to-end uncertainty quantification. In aggregate, these challenges motivate the establishment of a framework that integrates AI methods at multiple scales to perform scientifically meaningful interpretation, modeling, and uncertainty quantification of multiple streams of incoming data.

Our vision of the Scientific Autonomous Reasoning Agent (SARA) () is to develop a fully autonomous HT materials discovery and exploration framework by integrating robotic HT materials synthesis () with AI instances to accelerate both materials synthesis and analysis. In particular, SARA aims to automate the representation, characterization, planning, optimization, and learning of materials knowledge in a fully integrated manner. To achieve this goal, we envision the deployment of agents, which individually specialize on specific subtasks but closely interact with each other to accelerate the discovery efforts. These agents include, but are not limited to, synthesis and probing robotics to conduct experiments and highly optimized, physics-based AI models that evaluate currently available data with their associated uncertainties and that drive AI-guided discovery.

In this work, we take strides toward realizing this vision and present a fully integrated, autonomous framework that iteratively maps out the synthesis phase boundaries of metastable compounds in a closed-loop fashion. To this end, we incorporate a system of nested () cycles harnessed by SARA’s specialized AI agents to synthesize and explore thin-film libraries with lateral gradient laser spike annealing (lg-LSA) (): An internal (highest frequency) autonomous loop iteratively proposes optimized property measurements of a given lg-LSA stripe using a hierarchy of optical characterization techniques, while an external autonomous loop proposes and executes the next lg-LSA synthesis via a model that aggregates knowledge obtained by inner loop iterations. This architecture can be readily expanded and nested into higher-level loops that, e.g., optimize thin-film deposition, materials systems, and quantum materials computation.

SARA’s nested synthesis, microscopy imaging, and reflectance spectroscopy loops driven by the specialized AIs with AL reflect the hierarchical nature of scientific discovery. A primary goal of studying PCSP relationships is the enumeration of all possible syntheses that yield unique materials, a knowledge base that must be built from synthesis phase diagrams over a broad range of synthesis techniques, multiple parameter spaces defined within each technique, and many experimental campaigns to map synthesis phase diagrams in those spaces. Coordination among the levels of hierarchy is critical for maximizing high-level knowledge generation from low-level experiments, which guides our development of nested AL algorithms that seamlessly incorporate task coordination and uncertainty propagation. This framework is extensible with respect to incorporation of additional levels of hierarchy and/or expansion of techniques, such as additional property measurements and on-the-fly quantum mechanical calculations (, ), that enrich the knowledge within a given level of hierarchy. Networking of capabilities and knowledge sources elevates the use of AI and AL from process optimization to accelerated scientific discovery, a grand vision of AI-assisted science.

RESULTS

Our goal is to explore synthesis phase diagrams, especially the relatively unexplored ultrafast-annealing region where metastable polymorphs of metal oxides are more likely to form. These metastable oxide materials often exhibit improved properties over thermodynamic ambient ground states and are relevant for countless industrial applications. The cubic high-temperature polymorph of ZrO2, for example, is frequently used as a thermal coating material (–) because of its low thermal conductivity, while the anatase phase of TiO2 has attracted interest as a photocatalytic material (–). These are only two of the most prominent examples of materials systems where metastable phases outperform their respective equilibrium counterparts.

Here, we study the Bi-O system, which exhibits a rich phase diagram with dozens of experimentally observed polymorphs. In particular, we focus on the Bi2O3 composition, for which five distinct crystalline phases are known (, ). The monoclinic α-Bi2O3 is the thermodynamic ground state at room temperature, while four high-temperature phases have been reported: tetragonal β-Bi2O3, body-centered cubic γ-Bi2O3, cubic δ-Bi2O3, and orthorhombic ϵ-Bi2O3. The metastable δ phase has attracted interest as a solid oxide electrolyte in fuel cells (): Because of its defective fluorite-type crystal structure with a high concentration of oxygen vacancies, δ-Bi2O3 has the highest oxygen ion conductivity of any solid oxide known to date. Unfortunately, it exhibits only a narrow thermodynamic stability window between 727° and 824°C, which has so far precluded its use on an industrial scale. Substitution of yttrium or rare earth oxides can stabilize δ-Bi2O3 to room temperature but leads to a degraded ion conductivity. Hence, efforts have been aimed at finding routes to retain phase-pure δ-Bi2O3 to ambient conditions ().

Our samples are deposited as amorphous thin films by reactive sputtering on a silicon substrate. For other materials systems, composition spreads can be similarly deposited, allowing the mapping of a composition gradient c(x) to the location x on the substrate. We process the thin-film libraries using lg-LSA to form and kinetically trap metastable phases during the quench to ambient conditions. In contrast to conventional methods for annealing thin-film samples, such as hot plate, furnace, and rapid thermal annealing (), lg-LSA allows a controlled and rapid thermal processing over a wide range of conditions in a spatially confined region of less than 1 mm, with quench rates of 104 to 107 K/s and peak temperatures Tp up to 1400°C (limited by melt of the silicon substrates). Scanning a laser beam with a bi-Gaussian–like power profile (see the backdrop in the left panel of Fig. 1) over the film allows a single lg-LSA stripe to produce a spatially inhomogeneous thermal profile T(x) (where x runs across the stripe). The duration of heating is characterized by a dwell time τ defined by the ratio of the laser full width at half maximum (FWHM) divided by the scan velocity of the laser (typical dwells range from 100 to 10,000 μs). Hence, at a given dwell time τ, a single lg-LSA experiment produces a continuous range of temperature conditions wherein phase transitions, including formation of the sought metastable phases, need to be detected with a speed and level of automation commensurate with this robotic synthesis procedure to fully capitalize and elevate high-throughput synthesis to high-throughput discovery of phase boundaries.

Fig. 1.

SARA’s closed-loop autonomous materials synthesis and discovery cycle.

Starting from a set of initially selected processing conditions, SARA synthesizes an lg-LSA stripe on a thin-film library and subsequently sends it to its characterization AI agent, XAI. (Left) Schematic illustration of the lg-LSA/camera setup with the bi-Gaussian power profile in the backdrop, laser to the left, camera to the right, and thin-film sample mounted on a stage. Using a hierarchy of characterization techniques, XAI analyzes the stripe to determine intricate changes in its optical properties. In particular, XAI first acquires a microscope image to determine the positions of likely phase boundaries, which informs the reflectance spectroscopy measurements. XAI’s physics-informed AL model accelerates the spectroscopy acquisition, resulting in an accurate gradient model of the lg-LSA stripe. (Middle) Microscope image, reflectance spectroscopy heatmap, and first four Legendre coefficients from the XAI representation for a representative lg-LSA stripe. Lastly, the gradients are fed into SARA’s synthesis AI agent, ΣAI, which generates a gradient phase boundary map and also proposes the next experimental processing conditions to improve the phase boundary with as few experiments as possible. (Right) Model gradient phase map showing high-gradient regions in yellow.

To reduce both computational and experimental cost, we need to autonomously map out the processing phase space {x, τ, T} with as few synthesis experiments and property measurements as possible. Because the lg-LSA is an irreversible method, a specific position x [and potentially its associated composition c(x) in the presence of a composition gradient] can only be annealed once, further emphasizing the need for optimizing the selection of the processing conditions. Once an lg-LSA stripe is processed, a conclusive structural characterization across the thermal gradient is possible with grazing-incidence high-intensity x-ray diffraction (XRD) to resolve the crystal structure (). However, access to synchrotron facilities capable of producing x-rays with appropriate wavelength, intensity, and micrometer-scale spatial resolution comprises an inherently limited resource that motivates development of alternative phase boundary detection methods. To address this issue, we developed a complementary technique based on microscopy imaging and optical spectroscopy to rapidly assess phase boundaries. We recently demonstrated that structural phase changes are directly associated with changes in the optical thin-film properties of transparent films, in particular the optical thickness nd (), where n is the refractive index and d is the film thickness. Essentially, the gradients of the optical measurements across an lg-LSA stripe provide a means to map out phase boundaries without explicit crystallographic phase identification, thereby producing an unlabeled processing phase diagram without costly XRD experiments.

Here, we put forth how SARA integrates lg-LSA synthesis and optical phase boundary detection in a hierarchical autonomous workflow by using characterization and synthesis agents, XAI (pronounced chi AI) and ΣAI (pronounced sigma AI), respectively, as illustrated in Fig. 1. Starting with an initial processing condition, SARA synthesizes an lg-LSA stripe on a thin-film library. Then, SARA uses its internal characterization agent XAI to probe the stripe using a set of optical techniques: (i) microscopy imaging to rapidly inspect the anneal stripe (see the top panel in “Optical characterization” in Fig. 1) and (ii) more elaborate, but costly, reflectance measurements (see “Reflectance spectroscopy” in Fig. 1). In particular, XAI uses the observed features from the micrograph as prior knowledge to guide and acquire an accurate reflectance map with as few measurements as possible. The gradients of the reflectance map are then fed into SARA’s synthesis AI agent ΣAI, which incorporates the reflectance gradient information of each lg-LSA stripe into a phase boundary map as a function of the parameters {x, τ, T}. The high-gradient regions of this map determine the boundaries between phase fields and produce an unlabeled processing phase diagram (see “Phase boundary mapping” in Fig. 1). ΣAI is also responsible for proposing the next most promising synthesis conditions to effectively explore the search space. We discuss XAI and ΣAI in detail below.

XAI: Accelerating data acquisition and characterization

XAI’s primary task is to construct an accurate reflectance spectroscopy map r(x, λ) of an lg-LSA sample annealed at Tp and τ while measuring it at as few positions x across the stripe as possible. Because the acquisition time for a single such measurement r(x,·) is around 4.5 s, an exhaustive scan across a stripe of 1.5 mm in 10−μm intervals requires more than 11 min, forming one of the main bottlenecks of our HT experimental setup. To accelerate the reflectance data acquisition, we propose an AL scheme that takes advantage of multimodal measurements and incorporates physical structure into a Gaussian process (GP) regression model to yield highly optimized data acquisition and analysis.

The overall workflow of the XAI cycle is outlined in Fig. 2A. In a first step, SARA captures a microscope image of an lg-LSA stripe to analyze the overall condition of the anneal and to extract key features. This single RGB (red-green-blue) image of a stripe is inherently throughput-matched to the lg-LSA synthesis, producing prior knowledge for the XAI’s AL cycle to accelerate reflectance measurements. A representative microscope image is shown in Fig. 2B. These micrographs can be used to rapidly assess the conditions and the integrity of the anneal. Obvious damage of the thin film such as delamination and scratches or contamination such as dust particles, residual lithography artifacts, and dirt can be easily detected, which invalidates the lg-LSA stripe and can trigger resynthesis. The incorporation of such automated quality control in the autonomous loop alleviates responsibility for the XAI loop to effectively respond to invalid data, a critical aspect of autonomous workflows for robust operation ().

Fig. 2.

The characterization AL loop to accelerate acquisition of the reflectance spectroscopy necessary for phase boundary detection.

(A) Overall workflow. A microscope image (B) is captured to extract the stripe features, which are fed as a scaling function to the XAI kernel. The core features are the LSA prior and the RGB transition prior, which are sums of generalized Gaussian functions, as shown in (C). The corresponding gradient peak positions are denoted as yellow vertical lines in (B). XAI takes these functions as prior knowledge to set up a stripe-specific kernel that facilitates rapid model convergence. The AL loop is performed iteratively on reflectance measurements r(x, λ) over positions x, which are expanded into Legendre polynomials to reduce the dimensionality (see the middle panel in Fig. 1). (D) Performance of different kernel designs, illustrating that our XAI kernel with both LSA (XAI + LSA) and LSA + RGB (XAI + RGB) priors outperforms other conventional kernels. (E) Performances of the different acquisition function (R, random; U, uncertainty; IU, integrated uncertainty; IGU, integrated gradient uncertainty sampling). The solid lines represent the XAI + RGB kernel, while the dashed lines correspond to the RBF kernel.

SARA proceeds by constructing a stripe-specific GP kernel that incorporates the underlying physics of both the lg-LSA and optical spectroscopy processes. Notably, the bi-Gaussian power profile produces stripes of nearly perfect lateral symmetry at steady state, with their centers reaching the corresponding peak temperatures Tp and the continuous variation in lateral thermal gradient mirrored on each side of the stripe. We incorporate this structure into the kernel of XAI by forcing its main component to be symmetric around the center of a stripe (see the “XAI” section in the “materials and methods” part). In addition, SARA extracts key features of the stripe texture from the micrograph to further improve the kernel design, i.e., by identifying the stripe center and by detecting systematic optical changes across the stripe that we associate with structural transitions (, ). These optical transitions are identified by peaks in the gradient signal across a stripe, the locations of which are shown as vertical yellow lines in Fig. 2B. Furthermore, the two outermost detected peaks in the gradient signal give an estimate of how wide the lg-LSA stripe is, i.e., where the unannealed, amorphous film ends and the crystallization begins. We use slightly broadened peaks in the RGB gradient signal (purple line in Fig. 2C) and the overall width of the lg-LSA stripe (red line in Fig. 2C) as the RGB and LSA prior, respectively. These two functions are then used to rescale the kernel of the GP in the XAI cycle. Lastly, we account for small thickness variations of the film across the stripe by adding a linear component to the kernel.

To improve the efficiency of XAI, the reflectance r(x, λ) at any position x is expanded in Legendre polynomials as a function of wavelength λ before it is fed into the GP. Because the reflectance varies smoothly with λ, the Legendre expansion can be truncated between the 10th and 20th order at essentially no loss in accuracy () (see fig. S1), which reduces the dimensionality from the 2046 measured photon wavelength to a compact space of 10 to 20 Legendre coefficients. For our system, we use 16 coefficients throughout. Figure 1 (bottom middle) shows the first four Legendre coefficients of the reflectance data and our GP model’s posterior predictive mean and uncertainty for those coefficients.

To demonstrate the advantage of our specialized XAI kernel with respect to a set of conventional kernels, we perform statistical benchmarks on 617 lg-LSA experiments at distinct conditions, (T, τ). For each of the stripes, we measure the reflectance at n randomly selected positions on a grid spaced 10 μm apart and use these measurements as inputs to a GP model with different kernels. The ground truth is exhaustively measured across the whole stripe ranging over 1.5 mm, corresponding to a total of 151 measurements. For a range of n, we repeat this test 32 times for every stripe with independent random locations and average the coefficient of determination R2 for each kernel on the exhaustive data. This reduces the statistical noise in the results to a negligible value. Furthermore, we benchmark every kernel with a range of length scales and select the best in terms of R2 score (see the “Error metrics” section). By construction, our benchmark disentangles the effects of AL and kernel design, and the kernel with the right inductive bias will express the data best, even if all measurements are random. The results of this benchmark are illustrated in Fig. 2D, showing the performance of the various kernels with respect to the number of random measurements n. The radial basis function (RBF) kernel performs poorly, barely reaching an R2 score of 0.8 within 37 measurements. The Matérn kernel performs better, requiring n = 25 to reach the same score. The XAI kernels perform best: Depending on whether prior knowledge from the microscope image is included (“XAI + LSA” with LSA prior only and “XAI + RGB” with LSA and gradient peak prior), we obtain an R2 value of 0.8 with as few as 16 sampling points. The precise modeling of the optical measurements, its incorporation into the AL model, and the model’s initialization with the RGB image prior knowledge all contribute to the fast learning rate at the onset of autonomous experimentation, as required for efficient AL.

Having designed the kernel for the XAI cycle, we turn our attention to the acquisition function, that is, the function that chooses the next measurement based on the available information. An important component of many performant acquisition functions is the reduction of uncertainty in a target variable. Here, we benchmark three different acquisition functions, two of which are nonstandard. In particular, we study uncertainty (U) sampling, which chooses the next measurement at the point of maximum uncertainty in the Legendre coefficients; integrated uncertainty (IU) sampling, which selects the point that minimizes the integrated uncertainty over the whole sampling domain; and integrated gradient uncertainty (IGU) sampling, which is similar to IU but reduces the overall uncertainty in the gradients of the model. The last strategy targets our quantity of interest, because the reflectance gradients are indicative of the phase boundaries in the processing phase diagram. For this reason, we quantify the error of the model in the gradients, rather than the error to the observed data. Because we cannot directly observe the gradients, we generate ground truth data by training our GP model on the exhaustive measurements and taking the derivative of the fitted model. We then record the R2 score of the derivatives of the model for each of the acquisition functions at every iteration.

In Fig. 2E, we show the performance of the various sampling strategies as a function of AL iteration i, using either the XAI + RGB kernel (solid lines) or the RBF kernel (dashed lines). The best performance is achieved with the stripe-specific, highly optimized XAI + RGB kernel in conjunction with IGU sampling, reaching an R2 score of 0.8 and 0.9 within 9 and 15 iterations, respectively. Note that random sampling with the best kernel design still outperforms the best sampling strategies with the worst kernel. Furthermore, the acquisition functions do not differ markedly with the XAI + RGB kernel, highlighting the importance of incorporating the problem structure into our AI model and AL cycle. Compared to random sampling with an RBF kernel, the best strategy accelerates the acquisition and characterization by a factor of 9.7 for an R2 of 0.8, approximately one order of magnitude.

ΣAI: Accelerating phase exploration and processing conditions

Once an lg-LSA stripe has been processed by XAI, its output reflectance gradient information is fed into the external synthesis AI agent, ΣAI. Its main task is threefold: (i) assemble the incoming data, (ii) propagate uncertainty from every lg-LSA experiment to predict the gradient signal and its uncertainty throughout the search space, and (iii) ultimately propose new conditions for the synthesis experiments. The overall workflow of this process, which integrates the techniques described below, is shown in Fig. 3A.

Fig. 3.

The synthesis AL loop to accelerate materials exploration.

(A) Overall external workflow. Starting from an initial set of conditions, an lg-LSA stripe is annealed and processed by XAI. The gradients are then fed into the ΣAI agent, which constructs a (preliminary) gradient phase map and proposes the next experimental conditions. (B) The transformation of the XAI reflectance gradients requires a rigorous error assessment and propagation. Because of the symmetric XAI kernel, only one side from the stripe center is sampled on a uniform temperature mesh. The errors propagated to the ΣAI stem from the variation in the peak temperature (peak error) and the gradient of the temperature profile (profile error). (C) Performance of different ΣAI acquisition strategies: random (R), uncertainty (U), stripe uncertainty (SU), and upper confidence bound (UCB) sampling. The solid and dashed lines correspond to GP regression with and without input uncertainty, respectively. (D) Gradient phase map of Bi2O3, where the peak ridges are highlighted with light lines. The phase regions are labeled a posteriori with selected XRD measurements, from low to high temperatures: (i) amorphous as-deposited; (ii) rearranged, densified amorphous; (iii) δ-Bi2O3; (iv) mixed-phase region of δ-Bi2O3 and β-Bi2O3; (v) pure β-Bi2O3; and, lastly, (vi) melt-quenched amorphous.

The optical data of an lg-LSA anneal are processed through the nested XAI loop, the output of which is the gradients of the reflectance spectroscopy across a stripe, g(x) = ∥∂(x,·)∥2. This spatial gradient information is then transformed onto a temperature scale based on the Gaussian-type temperature profile (x) shown in Fig. 3B (blue line). Because the XAI kernel is symmetric up to the linear term, the gradient information is symmetric about the peak temperature T so that we only need to sample g(x) along one side from the stripe center (orange crosses in Fig. 3B).

In principle, one single lg-LSA stripe would produce the complete temperature conditions between room temperature T and T at a given dwell time τ. Hence, the set of metastable materials and their transition conditions would be available from a single stripe if one selected a high Tp (e.g., 1400°C) and Tmin = Tr. In practice, the concomitant increase in temperature gradient with Tp would require progressively higher spatial resolution to characterize the full range of transitions and results in undesirably high uncertainty in the modeled temperature. With our experimental characterization technique, the spatial resolution is limited to approximately 10 μm, and thus, the design of lg-LSA synthesis conditions must be done under consideration of the position-dependent temperature variation within a single spectroscopy measurement, which makes the selection of Tp at a given τ a nontrivial decision based on the aggregate information that can be gained from the entire lg-LSA stripe.

Properly propagating the multiple sources of uncertainty from synthesis and characterization through the model of the phase boundary map is extremely important: In standard GP regression, the inputs are assumed to be free of noise, but accounting for such errors is crucial when dealing with experimental measurements. Here, we include and propagate the uncertainties of the inputs due to two sources. First, the peak temperature reached in an lg-LSA anneal can vary within an error range of up to = 25°C at 1400°C due to fluctuation in the laser power, even after reaching steady state (Peak error in Fig. 3B). Second, the temperature profile itself gives rise to an error proportional to the spatial rate of change σ(x) ∝∣∂(x)∣, as shown by the green area in Fig. 3B. Note that the error bars in Fig. 3B are not to scale and intended solely for a schematic illustration.

As opposed to the XAI model of optical spectroscopy, the gradient map in synthesis space has no analogous physics-based model, in part because too few synthesis phase diagrams are known, and their underlying features remain an open question. The large dynamic range of dwell time motivates its logarithmic sampling, and the distinct influence of temperature and dwell time on synthesis motivates independent parameterization of these two dimensions. While we aim to learn more structured representations of synthesis phase diagrams in future refinements of the SARA framework, for the purposes of the present work, we find a Matérn kernel, with separate length scales for the temperature and dwell time dimensions, to enable rapid model convergence in the ΣAI loop while remaining flexible with respect to the gradient map in synthesis phase space.

In contrast to the XAI cycle, there is more opportunity to incorporate structure based on prior knowledge into the acquisition function, rather than the kernel, of ΣAI. As shown in Fig. 3C, random sampling performs only slightly worse than more sophisticated acquisition methods like uncertainty sampling or upper confidence bound (UCB) sampling in terms of , a generalization of R2 that takes into account the heteroscedasticity of the data due to the propagation of uncertainty (see the “Error metrics” section). This behavior can be understood by considering the following: Every experiment at {Tp, τ} produces a range of temperatures T < T at which new information is obtained, thereby reducing the uncertainty not only at Tp but also in a wide range of temperatures below it. Hence, uncertainty sampling at Tp and τ alone is a poor strategy. To address this issue, we introduce stripe uncertainty (SU) sampling, which takes into account the uncertainty in the whole temperature range between Tmin and Tp. This strategy greatly improves performance, reaching within 15 iterations.

The plot in Fig. 3D shows the gradient heatmap of Bi2O3 from an exhaustive sampling of all 617 lg-LSA stripes, with gradient peaks at every value of τ highlighted in white. Note that these peaks are connected and form ridges that can be well interpreted as phase field boundaries. To label these phase fields, we selectively collect and analyze XRD data of lg-LSA stripes annealed at conditions close to the field centers (see the Supplementary Materials for details). Notably, only few XRD measurements suffice to label the phase map, once the phase boundaries have been determined via the reflectance data. The phase field (i) below approximately 350°C corresponds to the as-deposited amorphous film, while a slight gradient ridge separates it from (ii), a densified, amorphous regime. At approximately 500°C, there is a ridge that extends across the complete dwell range, corresponding to the crystallization onset of the δ phase in domain (iii). The boundary separating (iii) from (iv) at approximately 550°C corresponds to the onset of a two-phase region, where both the δ and β phases of Bi2O3 coexist, and above approximately 650°C, we observe phase-pure β-Bi2O3 in (v). The gradient ridge between (iv) and (v) is particularly weak and wide in T, and we consider the bump arising at around 103.5 μs to be an artifact within the measurement uncertainty. The phase field above approximately 810°C corresponds to amorphous Bi2O3 that reforms after quenching from melt [the bulk melting temperature of Bi2O3 is 817°C ()] and stretches out across all values of τ.

A representative evolution of the actively learned gradient phase map is shown in Fig. 4, with six snapshots from (A) to (F). Figure 4A shows the preliminary gradient map at iteration n = 3: The two gradient ridges spanning all dwells qualitatively correspond to the crystallization boundaries from either melt (v and vi in Fig. 3D) or the deposited, densified thin film (ii and iii). At n = 8 in Fig. 4B, we detect the onset of the two-phase region (iii and iv), and at n = 15 (Fig. 4C), we detect the phase-pure β-Bi2O3 boundary (iv and v). With only n = 25 iterations, we identify the last boundary, namely, the amorphous densification onset (i and ii; see Fig. 4D). At this point, the overall features of the gradient phase map are already qualitatively captured completely, and subsequent iterations merely refine the boundary locations (Fig. 4E), getting closer to the exhaustive phase map in Fig. 4F.

Fig. 4.

The evolution of the actively learned gradient phase map of the Bi2O3 system at selected number of iterations n.

(A to E) We use the stripe uncertainty (SU) acquisition strategy, starting from a randomly selected condition (Tp, τ)1. The gradient ridges are shown as white lines, and the conditions (Tp, τ) at which the experiments have been performed are shown as white crosses (note that not all crosses are shown, because the plots have been cropped to a smaller range than the range of experimental conditions). For each panel, the number of sampled conditions n is indicated at the top, together with the corresponding score. (F) Final exhaustively sampled phase map.

Two factors are crucial for ΣAI to achieve a factor of approximately 14 acceleration to reach compared to random sampling without propagation of input uncertainties. First, incorporating materials synthesis into our SARA discovery framework allows us to check for convergence of the phase diagram on the fly. Even with random sampling, the possibility of quantifying the progress and monitoring convergence in the gradient mapping informs us how well the phase space has been sampled, thereby substantially reducing the resource cost. Second, the comprehensive uncertainty propagation in conjunction with the stripe uncertainty acquisition function realizes the full potential of AI and AL and decreases the required samples to a fraction of the exhaustive measurements. SARA’s overall AL acceleration is the product of the acceleration factors of XAI and ΣAI due to the cycles’ nested design.

DISCUSSION

In conclusion, we have developed SARA, an AI-driven autonomous closed-loop materials discovery framework that integrates robotic materials synthesis with automated microscopy imaging and reflectance spectroscopy characterization. SARA incorporates a set of nested AL loops based on specialized physics-inspired GP regression models to synthesize, characterize, and iteratively explore nonequilibrium synthesis phase maps using HT lg-LSA thin-film processing. In particular, SARA tightly integrates the physics of the experiments and quantifies experimental uncertainties in both the inputs and the outputs of the model. We highlight SARA’s capabilities on the Bi2O3 system by showing that SARA reduces the time to map the system’s phase boundaries by more than two orders of magnitude, in contrast to random or exhaustive searches. In particular, SARA identifies the synthesis conditions that trap metastable δ-Bi2O3 at room temperature, a promising solid oxide electrolyte. While scaling-up synthesis is a challenge for future work, a flat-top laser profile could be applied to anneal Bi2O3 at these processing conditions and integrated into fabrication of thin-film solid oxide fuel cells in portable power applications, or other integrated solid-state microelectromechanical system (MEMS) devices (, ).

The speedup in synthesis and data acquisition achieved by SARA is a fundamental prerequisite for paving the path toward exploratory HT efforts with additional chemical degrees of freedom and extended processing parameters and when targeting property optimization. The gradient phase map construction can be extended to additional degrees of freedom, e.g., on composition spreads over a continuous range of chemistries. While our current lg-LSA synthesis is limited to inorganic thin films that are transparent to infrared (IR) laser radiation (), e.g., complex oxides, future efforts are aimed at incorporating IR-transparent substrates that will allow the autonomous synthesis and processing of broader materials classes, e.g., metals and alloys. Furthermore, techniques that enable the characterization of in situ lg-LSA processing would provide us the means to better understand the transformation kinetics of metastable phases and improve our design of physics-inspired AI models. SARA’s nested AI architecture also allows the incorporation of additional agents for multi-objective optimization efforts by including robotic measurements of target properties. In addition to phase boundary mapping, research objectives for which SARA would enable new modalities of materials design include the following: discovery of a synthesis condition for a not-yet-synthesized phase, extension of the optical spectroscopy to characterize visible absorption to identify syntheses of materials for solar energy applications, and incorporation of new performance characterization such as electrical conductivity measurements. These latter examples involve mapping of synthesis phase diagrams in the context of performance metrics for a target application, the central goal of studying PCSP relationships. SARA’s autonomous execution of these studies constitutes a grand vision of AI-assisted materials science.

MATERIALS AND METHODS

Experiments and measurements

Thin-film deposition

We used thermally oxidized (200 nm oxide), highly doped (p-type, 0.01 to 0.02 Ohm·cm) Si wafers with lithographically patterned gold alignment marks as substrates for our thin-film deposition. Radio frequency (RF) reactive sputtering from a Bi target in an atmosphere of 8-mTorr Ar and 2-mTorr O2 was used to deposit the Bi2O3 sample in a custom-built sputter system. The substrate was rotated while operating the target at an RF power of 20 W to create a 170-nm-thick film with <10% thickness variation.

Lateral gradient laser spike annealing

The lg-LSA anneals were conducted using a continuous-wave CO2 laser operating at λ = 10.6 μm and maximum power of 125 W, which was configured to produce a power profile with a bi-Gaussian shape (320-μm-wide lateral FWHM and 80-μm-long longitudinal FWHM). To reach steady state, each anneal was conducted on a 5-mm-long stripe, with peak temperatures ranging from 400° to 1300°C and processing dwell times between 250 μs and 10 ms. The stripes were located 2 mm apart to avoid thermal overlap between anneals. With this configuration, a 100-mm-diameter wafer offers space for a total of up to 625 stripes with distinct anneal conditions. Note that the dwell τ is related to the scan velocity v via the FWHM of the laser in the scan direction (longitudinal) through . τ is approximately the time scale during which the temperature is within 5% of the peak temperature (). To avoid potential location bias on the wafer arising from variations in film thickness, the anneal locations were randomized across the thin film with respect to Tp and τ. In total, we annealed 617 lg-LSA stripes on our Bi2O3 sample with 400 ≤ Tp ≤ 1300∘C and 250 ≤ Tp ≤ 10,000 μs.

Microscopy imaging

We used a Thorlabs complementary metal-oxide semiconductor camera (RGB channels with 1024 × 1280 pixels), which was aligned normal to the sample, together with a coaxial illumination using white light over a spot size of approximately 1 mm in diameter. The camera magnification was set to produce a field of view of approximately 1 mm horizontally, resulting in a spacing of 0.92 μm between pixels. The raw microscope images and the software to process them are available online (, ).

Reflectance spectroscopy

A white light source (400 < λ < 900 nm) was focused down to a single 10-μm-diameter spot using optical fibers to locally illuminate the sample, allowing spatially resolved reflectance measurements. We used a flame spectrometer from Ocean Optics to collect the reflectance spectroscopy with an optimized integration time of ≈4500 ms. The reflected light was measured from λ = 340 to λ < 1026 nm at 2046 discrete values. The reflectance data were calibrated and normalized with respect to a dark reference spectrum, and a spectrum from an Ag-coated mirror. For the exhaustive reflectance measurements, the optical fiber was scanned across an lg-LSA stripe over a range of 1.5 mm in 10-μm increments, leading to 151 samples per stripe. The raw reflectance data and the software to process are available online (, ).

X-ray diffraction

The XRD data were collected using the ID3B beamline at the Cornell High Energy Synchrotron Source (CHESS) with a 9.7-keV beam, which was focused to a spot size on the sample of 20 μm by 40 μm at a 2∘ angle of incidence. A Pilatus 300K detector was used to capture the diffracted signal. The XRD data were collected every 10 μm across a stripe with a 50-ms integration time for each frame. The 2D detector data were integrated along the χ direction using pyFAI ().

Computational methods

In the following, bold lowercase letters refer to vectors and bold uppercase letters refer to matrices. Given a collection of inputs X = [x1, …, x] of a function f, we let f(X) be the result of the application of f to each column of X, f(X) ≔ [f(x1), …, f(x)].

Gaussian processes

A GP is a distribution over functions whose finite-dimensional marginal distributions are multivariate normal. That is, for any sample f of a GP and any finite selection of inputs X, we have f(X) ∼ N(μ, Σ), for some mean vector μ and covariance matrix Σ. Analogous to the multivariate case, a GP is completely defined by its first and second moments: a mean function μ(·) and a covariance function κ(·,·), also known as a kernel. In particular, if f ∼ GP(μ, κ), then for any finite collection of inputs Xwhere κ(X, X) is the matrix whose (i, j)nth entry is κ(x, x). Fortunately, the posterior mean μp and posterior covariance κp of a GP conditioned on observations with normally distributed noise have closed forms and only require linear algebraic operationswhere for homoscedastic regression and σ is the SE of the target y. Because a GP’s behavior is chiefly determined by the kernel, its performance can be improved markedly by incorporating important problem structure into the kernel. For more background on GPs, see (). For the present work, we developed a GP framework in Julia () with which we implemented SARA’s AL technology ().

Active learning

The field of AL considers the problem of selecting data in an optimal way to reduce the total amount of data that is required to effectively train a model (, ). To this end, the notion of an acquisition function is important. An acquisition function a(X, y) depends on currently available data and outputs a suggested observation x*. For example, if f∣X, y denotes the posterior of f after having seen the data, and var(f∣X, y) is the posterior variance (itself a function), thendefines an acquisition strategy known as uncertainty sampling. Other acquisition functions are based on upper-confidence bounds, expected improvement, and probability of improvement. Overall, an important ingredient for AL is the quantification of uncertainty, which is a strength of Bayesian models. In the realm of Bayesian models, GPs are of particular importance because of their unique combination of flexibility, closed-form inference formulas, and uncertainty quantification. For these reasons, we chose to build SARA’s computational backbone on GPs.

Input noise

Because of the importance of uncertainty quantification for AL, it is critical to take all sources of uncertainty into account. In the case of SARA, it is crucial to account for errors in not only the measurements (i.e., model outputs) but also the experimental conditions (i.e., model inputs) due to intrinsic experimental uncertainties in the temperature profile. However, the general problem of posterior inference with input noise is intractable. For this reason, one needs to use approximate methods like variational approximations (, ) and Markov chain Monte Carlo () or methods that transform the problem of homoscedastic regression with input noise to one of heteroscedastic regression () without input noise (–). A particularly efficient technique is that of McHutchon and Rasmussen (), which is based on propagating the input uncertainty using a linear approximation of the standard posterior mean. According to this model, given the regular posterior mean μp(x), the input noise–corrected version can be computed by updatingin Eq. 2 for the GP posterior. Notably, we generalize the original work in making the input uncertainty σ(X) dependent on the input. This is possible because the non-constant uncertainties in SARA’s experimental process can be estimated well by physical considerations (see the “ΣAI” section for details). Lastly, note that Eq. 4 makes the approximate posterior uncertainty dependent on the values of the data via the posterior mean, not just the locations of the measurements.

XAI

The goal of XAI is to infer the reflectance r(x, λ) using the least number of measurement locations x as possible. Each measurement of the inner loop acquires the wavelength-dependent spectroscopic reflectance of the underlying thin film, that is, a vector whose entries correspond to reflectance intensities at a given wavelength.

To aid the efficiency of our model, we first reduce the dimension of the output by projecting it onto the basis of a small number (10 to 20) of Legendre polynomials. Because the signal is smooth as a function of wavelength, it admits a sparse approximation in this basis, allowing the compression of the signal with virtually no loss of information () [see also the Supplementary Materials]. The AL cycle then works on the dimensionality-reduced form of the reflectance data.

In the following, we describe the construction of the XAI kernel function, which integrates special structure of the data and is a critical part of XAI. In particular, the kernel incorporates (i) lateral symmetry, (ii) variance scaling based on RGB data, and (iii) asymptotically linear behavior. Starting with a Matérn 5/2 kernel k with a length scale l, we symmetrize it via ksym(x, y) = k(x, y) + k(x − c, y − c) around the stripe center c, which we estimate from the RGB images. We incorporate further information from the RGB images by scaling the kernel with the LSA or RGB prior function frgb shown in Fig. 2C. In particular, we use the peaks in the RGB gradient signal, slightly broaden them by a Gaussian with σ = 20 μm, and sum them to our RGB prior function (purple line in Fig. 2C). In addition, the overall width of the lg-LSA stripe gives rise to the LSA prior, which is a generalized Gaussian with a wide shape parameter of β = 4 and a scale parameter σ defined by the stripe width (red line in Fig. 2C). frgb is then given by a weighted sum of these two prior functions. This scaling constrains the search space, because we do not expect a lot of change in the underlying material if the experimental conditions (e.g., temperature) stay similar, and gives rise to the kernel frgb(x)ksym(x, y)frgb(y). Lastly, we incorporate an asymptotically linear behavior, due to thickness variations in the wafer, with the linear kernel kline(x, y) = x · y + b, where b is a constant that controls the variance of the bias term of the line. As a result, the XAI kernel for one Legendre coefficient is proportional to

For all the Legendre coefficients, we then use a GP with the kernel aX(x, y), where {a} are scaling coefficients that incorporate the different variances of the Legendre coefficients, to learn the reflectance map. This can also be interpreted as computing a vector-valued GP with the matrix-valued kernel KX(x, y) = diag (a) kX(x, y), where a is the vector of scaling coefficients. For a comprehensive review on matrix-valued kernels, see (). The length scale l of the Matérn kernel can be optimized via maximization of the marginal likelihood (). However, to make the reported results in Fig. 2D independent of this nonconvex optimization procedure, we ran the benchmarks using a range of fixed length scales and reported the best performing combination for each kernel.

Regarding the acquisition function, in addition to uncertainty sampling, we benchmark XAI using IU sampling, a policy that reduces the total variance over a set of potential measurement locations Z. In particular, IU is defined bywhere X is the set of inputs and y is the set of outputs of the model. Note that we can calculate the quantity var(f ∣ X, y, x*) because the standard posterior GP variance only depends on the measurement location, not the value y*. Lastly, we note that the derivative of a GP is also a GP (). Plugging the derivative GP into Eq. 6 yields IGU sampling, which achieves the best performance in the XAI acquisition benchmark (see Fig. 2).

ΣAI

ΣAI works to identify phase regions and their boundaries in the temperature–dwell time space and, more generally, the processing-composition space. The raw reflectance data cannot be used directly for this task because of two reasons. First, the data are measured as a function of position, not temperature. Therefore, we convert the stripe-specific reflectance function (x, λ) to the temperature domain using the temperature profile , yielding (T, λ). Second, the reflectance varies not only with the phase behavior but also with the film thickness across the wafer. For this reason, we calculate the L2-norm of the rate of change of the spectroscopic reflectance, which is invariant to linear thickness variations of the film. In particular, for all T < Tp, we want to inferd quantifies how much the spectroscopic reflectance changes as a function of temperature and dwell time and is a strong indicator of phase changes (). Estimating the phase boundaries then reduces to getting an accurate estimate of d over all (T, τ) (and potentially composition c). This is the goal of the ΣAI loop.

Crucially, experimental errors can occur in x and, therefore, in T, making it imperative to quantify the uncertainty due to these input errors and propagate them to ΣAI. Our benchmarks show that ignoring these uncertainties leads to a substantial deterioration in AL performance (see Fig. 3C).

To this end, we now discuss the intrinsic experimental uncertainties due to the temperature profile (x) of the laser. In particular, we compute the variance of the true temperature around the value predicted by the temperature profile as a function of position bywhere is the SE in the peak temperature and σ is the SE in the position. The first term quantifies the error at the peak temperature, which is largest at high temperatures (1400∘C) and falls off linearly with T. The second term quantifies uncertainties of the temperature profile, which not only are due to limited spatial resolution but also encompass random asymmetries in the profile of the laser. The form of term is derived using SE propagation techniques (). For the results reported here, = 25∘C and σ = 50 μm.

The expression for the temperature uncertainty in Eq. 8 is then used in conjunction with Eq. 4 to compute a GP that comprehensively quantifies the uncertainties in the Legendre coefficients of the optical reflectance as a function of temperature. To compute d in Eq. 7, one simply sums the squared derivatives of the GPs of the Legendre coefficients of the reflectance

Because we have access to the uncertainties in from the GP, we can use uncertainty propagation techniques on Eq. 9 to calculate a first-order uncertainty estimate of d(T, τ). For the outer loop, we used a two-dimensional Matérn 5/2 kernel with different length scales across each dimension. This allows the GP to learn independent sensitivity parameters of the experiment for the input dimensions. Note that the ΣAI benchmarks in Fig. 3C were carried out with fixed length scales to disentangle the effects of different acquisition functions and hyperparameter learning.

For ΣAI, we designed an acquisition strategy that incorporates the property that a single stripe generates data throughout a range of temperatures. In particular, given experimental conditions x that give rise to a stripe (Tp, τ, etc.), we sum the uncertainties of all relevant observations x that are in the set Stripe(x) of conditions on the stripe x. In particular, we propose stripe uncertainty sampling:

Notably, one can use the same principle to generalize other acquisition functions. We investigated a stripe upper-confidence bound sampling policy. However, it performed worse or equal to the simpler stripe uncertainty sampling policy above. The synergy of the comprehensive uncertainty quantification and the stripe sampling function yields considerable benefits, as displayed in Fig. 3C.

Error metrics

In our benchmarks of the kernels and acquisition functions for the inner loop, we used the coefficient of determination R2 to measure performance, defined bywhere μ(y) is the mean of the data y. The advantage of using R2 over other canonical measures like the mean-squared error is that it is dimensionless and easily interpretable as the proportion of the variance of the data that is explained by the model f.

As R2 weighs the deviation at every data point equally, it is not an ideal measure for heteroscedastic data, like the optical gradient data of ΣAI. For this reason, we use a generalization of R2, based on the log-likelihood of the heteroscedastic normal errors, to measure performance in the ΣAI benchmarks. In particular, the measure is given bywhere σ is the SD of the ith error. For SARA, the σ are the product of the comprehensive uncertainty quantification of the experimental process. Clearly, reduces to R2 if the noise variances are all equal. If they are not, is a better measure of misfit, as it weighs the residuals of more certain data points stronger than those with greater uncertainty. Notably, similar pseudo-R2 scores based on log-likelihoods are used throughout statistics and applied fields (–).