Machine learning modeling of superconducting critical temperature

Superconductivity has been the focus of enormous research effort since its discovery more than a century ago. Yet, some features of this unique phenomenon remain poorly understood; prime among these is the connection between superconductivity and chemical/structural properties of materials. To bridge the gap, several machine learning schemes are developed herein to model the critical temperatures ($T_{\mathrm{c}}$) of the 12,000+ known superconductors available via the SuperCon database. Materials are first divided into two classes based on their $T_{\mathrm{c}}$ values, above and below 10 K, and a classification model predicting this label is trained. The model uses coarse-grained features based only on the chemical compositions. It shows strong predictive power, with out-of-sample accuracy of about 92%. Separate regression models are developed to predict the values of $T_{\mathrm{c}}$ for cuprate, iron-based, and"low-$T_{\mathrm{c}}$"compounds. These models also demonstrate good performance, with learned predictors offering potential insights into the mechanisms behind superconductivity in different families of materials. To improve the accuracy and interpretability of these models, new features are incorporated using materials data from the AFLOW Online Repositories. Finally, the classification and regression models are combined into a single integrated pipeline and employed to search the entire Inorganic Crystallographic Structure Database (ICSD) for potential new superconductors. We identify more than 30 non-cuprate and non-iron-based oxides as candidate materials.


INTRODUCTION
Superconductivity, despite being the subject of intense physics, chemistry and materials science research for more than a century, remains among one of the most puzzling scientific topics [1]. It is an intrinsically quantum phenomenon caused by a finite attraction between paired electrons, with unique properties including zero DC resistivity, Meissner and Josephson effects, and with an ever-growing list of current and potential applications. There is even a profound connection between phenomena in the superconducting state and the Higgs mechanism in particle physics [2]. However, understanding the relationship between superconductivity and materials' chemistry and structure presents significant theoretical and experimental challenges. In particular, despite focused research efforts in the last 30 years, the mechanisms responsible for high-temperature superconductivity in cuprate and iron-based families remain elusive [3,4].
Recent developments, however, allow a different approach to investigate what ultimately determines the superconducting critical temperatures (T c ) of materials. Extensive databases covering various measured and calculated materials properties have been created over the years [5][6][7][8][9]. The shear quantity of accessible information also makes possible, and even necessary, the use of data-driven approaches, e.g., statistical and machine learning (ML) methods [10][11][12][13]. Such algorithms can be developed/trained on the variables collected in these databases, and employed to predict macroscopic properties such as the melting temperatures of binary compounds [14], the likely crystal structure at a given composition [15], band gap energies [16,17] and density of states [16] of certain classes of materials.
Taking advantage of this immense increase of readily accessible and potentially relevant information, we develop several ML methods modeling T c from the complete list of reported (inorganic) superconductors [18]. In their simplest form, these methods take as input a number of predictors generated from the elemental composition of each material. Models developed with these basic features are surprisingly accurate, despite lacking information of relevant properties, such as space group, electronic structure, and phonon energies. To further improve the predictive power of the models, as well as the ability to extract useful information out of them, another set of features are constructed based on crystallographic and electronic information taken from the AFLOW Online Repositories [19][20][21][22].
Application of statistical methods in the context of superconductivity began in the early eighties with simple clustering methods [23,24]. In particular, three "golden" descriptors confine the sixty known (at the time) superconductors with T c > 10 K to three small islands in space: the averaged valence-electron numbers, orbital radii differences, and metallic electronegativity differences. Conversely, about 600 other superconductors with T c < 10 K appear randomly dispersed in the same space. These descriptors were selected heuristically due to their success in classifying binary/ternary structures and predicting stable/metastable ternary quasicrystals. Recently, an investigation stumbled on this clustering problem again by observing a threshold T c closer to log T thres c ≈ 1.3 T thres c = 20 K [25]. Instead of a heuristic approach, random forests and simplex fragments were leveraged on the structural/electronic properties data from the AFLOW Online Repositories to find the optimum clustering descriptors. A classification model was developed showing good performance. Separately, a sequential learning framework was evaluated on superconducting materials, exposing the limitations of relying on randomguess (trial-and-error) approaches for breakthrough discoveries [26]. Subsequently, this study also highlights the impact machine learning can have on this particular field. In another early work, statistical methods were used to find correlations between normal state properties and T c of the metallic elements in the first six rows of the periodic table [27]. Other contemporary work hones in on specific materials [28,29] and families of superconductors [30,31] (see also Ref. [32]).
Whereas previous investigations explored several hundred compounds at most, this work considers more than 16, 000 different compositions. These are extracted from the SuperCon database, which contains an exhaustive list of superconductors, including many closely-related materials varying only by small changes in stoichiometry (doping plays a significant role in optimizing T c ). The order-of-magnitude increase in training data (i ) presents crucial subtleties in chemical composition among related compounds, (ii ) affords family-specific modeling exposing different superconducting mechanisms, and (iii ) enhances model performance overall. It also enables the optimization of several model construction procedures. Large sets of independent variables can be constructed and rigorously filtered by predictive power (rather than selecting them by intuition alone). These advances are crucial to uncovering insights into the emergence/suppression of superconductivity with composition.
As a demonstration of the potential of ML methods in looking for novel superconductors, we combined and applied several models to search for candidates among the roughly 110, 000 different compositions contained in the Inorganic Crystallographic Structure Database (ICSD). The framework highlights 35 compounds with predicted T c 's above 20 K for experimental validation. Of these, some exhibit interesting chemical and structural similarities to cuprate superconductors, demonstrating the ability of the ML models to identify meaningful patterns in the data. In addition, most materials from the list share a peculiar feature in their electronic band structure: one (or more) flat/nearly-flat bands just below the energy of the highest occupied electronic state. The associated large peak in the density of states (infinitely large in the limit of truly flat bands) can lead to strong electronic in-stability, and has been discussed recently as one possible way to high-temperature superconductivity [33,34].

RESULTS
Data and predictors. The success of any ML method ultimately depends on access to reliable and plentiful data. Superconductivity data used in this work is extracted from the SuperCon database [18], created and maintained by the Japanese National Institute for Materials Science. It houses information such as the T c and reporting journal publication for superconducting materials known from experiment. Assembled within it is a uniquely exhaustive list of all reported superconductors, as well as related non-superconducting compounds.
From SuperCon, we have extracted a list of approximately 16, 400 compounds, of which 4, 000 have no T c reported (see Methods for details). Of these, roughly 5, 700 compounds are cuprates and 1, 500 are iron-based (about 35% and 9%, respectively), reflecting the significant research efforts invested in these two families. The remaining set of about 8, 000 is a mix of various materials, including conventional phonon-driven superconductors (e.g., elemental superconductor, A15 compounds), known unconventional superconductors like the layered nitrides and heavy fermions, and many materials for which the mechanism of superconductivity is still under debate (such as bismuthates and borocarbides). The distribution of materials by T c for the three groups is shown in Figure 2a.
Use of this data for the purpose of creating ML models can be problematic. Training a model only on superconductors can lead to significant selection bias that may render it ineffective when applied to new materials [35]. Even if the model learns to correctly recognize factors promoting superconductivity, it may miss effects that strongly inhibit it. To mitigate the effect, we incorporate about 300 materials found by H. Hosono's group not to display superconductivity [36]. However, the presence of non-superconducting materials, along with those without T c reported in SuperCon, leads to a conceptual problem. Surely, some of these compounds emerge as nonsuperconducting "end-members" from doping/pressure studies, indicating no superconducting transition was observed despite some efforts to find one. However, a transition may still exist, albeit at experimentally difficult to reach or altogether inaccessible temperatures (for most practical purposes below 10 mK) [37]. This presents a conundrum: ignoring compounds with no reported T c disregards a potentially important part of the dataset, while assuming T c = 0 K prescribes an inadequate description for (at least some of) these compounds. To circumvent the problem, materials are first partitioned in two groups by their T c , above and below a threshold temperature (T sep ), for the creation of a classification model. Compounds with no reported critical temperature can be classified in the "below-T sep " group without the need to specify a T c value (or assume it is zero).
For most materials, the SuperCon database provides only the chemical composition and T c . To convert this information into meaningful features/predictors (used interchangeably), we employ the Materials Agnostic Platform for Informatics and Exploration (Magpie) [38]. Magpie computes a set of attributes for each material, including elemental property statistics like the mean and the standard deviation of 22 different elemental properties (e.g., period/group on the periodic table, atomic number, atomic radii, melting temperature), as well as electronic structure attributes, such as the average fraction of electrons from the s, p, d and f valence shells among all elements present.
The application of Magpie predictors, though appearing to lack a priori justification, expands upon past clustering approaches by Villars and Rabe [23,24]. They show that, in the space of a few judiciously chosen heuristic predictors, materials separate and cluster according to their crystal structure and even complex properties such as high-temperature ferroelectricity and superconductivity. Similar to these features, Magpie predictors capture significant chemical information, which plays a decisive role in determining structural and physical properties of materials.
Despite the success of Magpie predictors in modeling materials properties [38], interpreting their connection to superconductivity presents a serious challenge. They do not encode (at least directly) many important properties, particularly those pertinent to superconductivity. Incorporating features like lattice type and density of states would undoubtedly lead to significantly more powerful and interpretable models. Since such information is not generally available in SuperCon, we employ data from the AFLOW Online Repositories [19][20][21][22]. The materials database houses nearly 170 million properties calculated with the software package AFLOW [6,[39][40][41][42][43][44][45][46][47]. It contains information for the vast majority of compounds in the ICSD [5]. Although the AFLOW Online Repositories contain calculated properties, the DFT results have been extensively validated with ICSD records [17,25,[48][49][50][51].
Unfortunately, only a small subset of materials in SuperCon overlaps with those in the ICSD: about 800 with finite T c and less than 600 are contained within AFLOW. For these, a set of 26 predictors are incorporated from the AFLOW Online Repositories, including structural/chemical information like the lattice type, space group, volume of the unit cell, density, ratios of the lattice parameters, Bader charges and volumes, and formation energy (see Supplementary Materials). In addition, electronic properties are considered, including the density of states near the Fermi level as calculated by AFLOW. Previous investigations exposed limitations in applying ML methods to a similar dataset in isolation [25]. Instead, a framework is presented here for combining models built on Magpie descriptors (large sampling, but features limited to compositional data) and AFLOW features (small sampling, but diverse and pertinent features).
Once we have a list of relevant predictors, various ML models can be applied to the data [52,53]. All ML algorithms in this work are variants of the random forest method [54]. Fundamentally, this approach combines many individual decision trees, where each tree is a nonparametric supervised learning method used for modeling either categorical or numerical variables (i.e., classification or regression modeling). A tree predicts the value of a target variable by learning simple decision rules inferred from the available features (see Figure 1 for an example).
Random forest is one of the most powerful, versatile, and widely-used ML methods [55]. There are several advantages that make it especially suitable for this problem. First, it can learn complicated non-linear dependencies from the data. Unlike many other methods (e.g., linear regression), it does not make any assumptions about the relationship between the predictors and the target variable. Second, random forests are quite tolerant to heterogeneity in the training data. It can handle both numerical and categorical data which, furthermore, does not need extensive and potentially dangerous preprocessing, such as scaling or normalization. Even the presence of strongly correlated predictors is not a problem for model construction (unlike many other ML algorithms). Another significant advantage of this method is that, by combining information from individual trees, it can estimate the importance of each predictor, thus making the model more interpretable. However, unlike model construction, determination of predictor importance is complicated by the presence of correlated features. To avoid this, standard feature selection procedures are employed along with a rigorous predictor elimination scheme (based on their strength and correlation with others). Overall, these methods reduce the complexity of the models and improve our ability to interpret them. Classification models. As a first step in applying ML methods to the dataset, a sequence of classification models are created, each designed to separate materials into two distinct groups depending on whether T c is above or below some predetermined value. The temperature that separates the two groups (T sep ) is treated as an adjustable parameter of the model, though some physical considerations should guide its choice as well. Classification ultimately allows compounds with no reported T c to be used in the training set by including them in the below-T sep bin. Although discretizing continuous variables is not generally recommended, in this case the benefits of including compounds without T c outweigh the potential information loss.
In order to choose the optimal value of T sep , a series of random forest models are trained with different threshold temperatures separating the two classes. Since setting T sep too low or too high creates strongly imbalanced classes (with many more instances in one group), it is important to compare the models using several different metrics. Focusing only on the accuracy (count of correctly-classified instances) can lead to deceptive re- Example of a single decision tree used to classify materials depending on whether Tc is above or below 10 K. A tree can have many levels, but only the three top are shown. The decision rules leading to each subset are written inside individual rectangles. The subset population percentage is given by "samples", and the node color/shade represents the degree of separation, i.e., dark blue/orange illustrates a high proportion of Tc > 10 K/Tc < 10 K materials (the exact value is given by "proportion"). A random forest consists of a large numbercould be hundreds or thousands -of such individual trees.
sults. Hypothetically, if 95% of the observations in the dataset are in the below-T sep group, simply classifying all materials as such would yield a high accuracy (95%), while being trivial in any other sense. To avoid this potential pitfall, three other standard metrics for classification are considered: precision, recall, and F 1 score. They are defined using the values tp, tn, f p, and f n for the count of true/false positive/negative predictions of the model: where positive/negative refers to above-T sep /below-T sep . The accuracy of a classifier is the total proportion of correctly-classified materials, while precision measures the proportion of correctly-classified above-T sep superconductors out of all predicted above-T sep . The recall is the proportion of correctly-classified above-T sep materials out of all truly above-T sep compounds. While the precision measures the probability that a material selected by the model actually has T c > T sep , the recall reports how sensitive the model is to above-T sep materials. Maximizing the precision or recall would require some compromise with the other, i.e., a model that labels all materials as above-T sep would have perfect recall but dismal precision. To quantify the trade-off between recall and precision, their harmonic mean (F 1 score) is widely used to measure the performance of a classification model. With the exception of accuracy, these metrics are not symmetric with respect to the exchange of positive and negative labels. For a realistic estimate of the performance of each model, the dataset is randomly split (85%/15%) into training and test subsets. The training set is employed to fit the model, which is then applied to the test set for subsequent benchmarking. The aforementioned metrics (Equations 1-4) calculated on the test set provide an unbiased estimate of how well the model is expected to generalize to a new (but similar) dataset. With the random forest method, similar estimates can be obtained intrinsically at the training stage. Since each tree is trained only on a bootstrapped subset of the data, the remaining subset can be used as an internal test set. These two methods for quantifying model performance usually yield very similar results.
With the procedure in place, the models' metrics are a b c d

FIG. 2.
SuperCon dataset and classification model performance. (a) Histogram of materials categorized by Tc (bin size is 2 K, only those with finite Tc are counted). Blue, green, and red denote "low-Tc", iron-based, and cuprate superconductors, respectively. In the inset: histogram of materials categorized by ln (Tc) restricted to those with Tc > 10 K. (b) Performance of different classification models as a function of the threshold temperature (Tsep) that separates materials in two classes by Tc. Performance is measured by accuracy (gray), precision (red), recall (blue), and F1 score (purple). The scores are calculated from predictions on an independent test set, i.e., one separate from the dataset used to train the model. In the inset: the dashed red curve gives the proportion of materials in the above-Tsep set. (c) Accuracy, precision, recall, and F1 score as a function of the size of the training set with a fixed test set. (d) Accuracy, precision, recall, and F1 as a function of the number of predictors.
evaluated for a range of T sep and illustrated in Figure 2b.
The accuracy increases as T sep goes from 1 K to 40 K, and the proportion of above-T sep compounds drops from above 70% to about 15%, while the recall and F 1 score generally decrease. The region between 5 − 15 K is especially appealing in (nearly) maximizing all benchmarking metrics while balancing the sizes of the bins. In fact, setting T sep = 10 K is a particularly convenient choice. It is also the temperature used in Refs [23,24] to separate the two classes, as it is just above the highest T c of all elements and pseudoelemental materials (solid solution whose range of composition includes a pure element). Here, the proportion of above-T sep materials is approximately 38% and the accuracy is about 92%, i.e., the model can correctly classify nine out of ten materialsmuch better than random guessing. The recall -quantifying how well all above-T sep compounds are labeled and, thus, the most important metric when searching for new superconducting materials -is even higher. (Note that the models' metrics also depend on random factors such as the composition of the training and test sets, and their exact values can vary.) The most important factors that determine the model's performance are the size of the available dataset and the number of meaningful predictors. As can be seen in Figure 2c, all metrics improve significantly with the increase of the training set size. The effect is most dramatic for sizes between several hundred and few thousands instances, but there is no obvious saturation even for the largest available datasets. This validates efforts herein to incorporate as much relevant data as possible into model training. The number of predictors is another very important model parameter. In Figure 2d, the accuracy is calculated at each step of the backward feature elimination process. It quickly saturates when the number of predictors reaches 10. In fact, a model with only 5 predictors achieves almost 90% accuracy. (Note that these are the five most informative predictors, selected by the model out of the full list of 145 ones.) For an understanding of what the model has learned, an analysis of the chosen predictors is needed. In the random forest method, features can be ordered by their importance quantified via the so-called Gini importance or "mean decrease in impurity" [52,53]. For a given feature, it is the sum of the Gini impurity [56] over the number of splits that include the feature, weighted by the number of samples it splits, and averaged over the entire forest. Due to the nature of the algorithm, the closer to the top of the tree a predictor is used, the greater number of predictions it impacts.
Although correlations between predictors do not affect the model's ability to learn, it can distort importance estimates. For example, a material property with a strong effect on T c can be shared among several correlated predictors. Since the model can access the same information through any of these variables, their relative importances are diluted across the group. To reduce the effect and limit the list of predictors to a manageable size, the backward feature elimination method is employed. The process begins with a model constructed with the full list of predictors, and iteratively removes the least significant one, rebuilding the model and recalculating importances with every iteration. (This iterative procedure is necessary since the ordering of the predictors by importance can change at each step.) Predictors are removed until the accuracy drops by no more than 2%, reducing the full list of 145 down to 5. Furthermore, two of these predictors are strongly correlated with each other, and we remove the less important one. This has a negligible impact on the model performance, yielding four predictors total (see Table 1) with an above 90% accuracy score -only slightly worse than the full model. Scat-TABLE 1. The most relevant predictors and their importances for the classification and general regression models. avg(x) and std(x) denote the composition-weighted average and standard deviation, respectively, calculated over the vector of elemental values for each compound [38]. For the classification model, all predictor importances are quite close. ter plots of the pairs of the most important predictors are shown in Figure 3, where blue/red denotes whether the material is in the below-T sep /above-T sep class. For comparison, we create another classifier based on the average number of valence electrons, metallic electronegativity differences, and orbital radii differences, i.e., the predictors used in Refs. [23,24] to cluster materials with T c > 10 K. A classifier built only with these three predictors is less accurate than both the full and the truncated models presented herein, but comes quite close: the full model has about 3% higher accuracy and F 1 score, while the truncated model with four predictors is less that 2% more accurate. The rather small (albeit not insignificant) differences demonstrates that even on the scale of the entire SuperCon dataset, the predictors used by Villars and Rabe [23,24] capture much of the relevant chemical information for superconductivity. Regression models. After constructing a successful classification model, we now move to the more difficult challenge of predicting T c . Creating a regression model may enable better understanding of the factors controlling T c of known superconductors, while also serving as an organic part of a system for identifying potential new ones. Leveraging the same set of elemental predictors as the classification model, several regression models are presented focusing on materials with T c > 10 K. It avoids the problem of materials with no reported T c with the assumption that, if they were to exhibit superconductivity at all, their critical temperature would be below 10 K. Another problem is that the T c 's are unevenly distributed over the T c axis (see Figure 2a). To avoid this, ln (T c ) is used as the target variable instead of T c (Figure 2a inset), which creates a more uniform distribution and is also considered a best practice when the range of a target variable covers more than one order of magnitude (as in the case of T c ). Following this transformation, the dataset is parsed randomly (85%/15%) into training and test subsets (similarly performed for the classification model).
Present within the dataset are distinct families of superconductors with different driving mechanisms for superconductivity, including cuprate and iron-based hightemperature superconductors, with all others denoted "low-T c " for brevity (no specific mechanism in this group). Surprisingly, a single regression model does reasonably well among the different families -benchmarked on the test set, the model achieves R 2 ≈ 0.88 (Figure 4a). It suggests that the random forest algorithm is flexible and powerful enough to automatically separate the compounds into groups and create group-specific branches with distinct predictors (no explicit group labels were used during training and testing). As validation, three separate models are constructed trained only on a specific family, namely the low-T c , cuprate, and iron-based superconductors, respectively. Benchmarking on mixed-family test sets, the models performed well on compounds belonging to their training set family while demonstrating no predictive power on the others. Figures 4b-d illustrate a cross-section of this comparison. Specifically, the model trained on low-T c compounds dramatically underestimates the T c of both high-temperature superconducting families (Figures 4b and c), even though this test set only contains compounds with T c < 40 K. Conversely, the model trained on the cuprates tends to overestimate the T c of low-T c (Figure 4d) and iron-based (Figure 4e) superconductors. This is a clear indication that superconductors from these groups have different factors determining their T c . Interestingly, the family-specific models do not perform better than the general regression containing all the data points: R 2 for the low-T c materials is about 0.85, for cuprates is just below 0.8, and for iron-based compounds is about 0.74. In fact, it is a purely geometric effect that the combined model has the highest R 2 . Each group of superconductors contributes mostly to a distinct T c range, and, as a result, the combined regression is better determined over longer temperature interval.
In order to reduce the number of predictors and increase the interpretability of these models without significant detriment to their performance, a backward feature elimination process is again employed. The procedure is very similar to the one described previously for the classification model, with the only difference being that the reduction is guided by R 2 of the model, rather than the accuracy (the procedure stops when R 2 drops by 3%).
The most important predictors for the four models (one general and three family-specific) together with their importances are shown in Tables 1 and 2. Differences in important predictors across the family-specific models reflect the fact that distinct mechanisms are responsible for driving superconductivity among these groups. The list is longest for the low-T c superconductors, reflecting the eclectic nature of this group. Similar to the general regression model, different branches are likely created for distinct sub-groups. Nevertheless, some important predictors have straightforward interpretation. As illustrated in Figure 5a, low average atomic weight is a necessary (albeit not sufficient) condition for achieving high T c among the low-T c group. In fact, the maximum T c for a given weight roughly follows 1/ √ m A .
Mass plays a significant role in conventional superconductors through the Debye frequency of phonons, leading to the well-known formula T c ∼ 1/ √ m, where m is the ionic mass. Other factors like density of states are also important, which explains the spread in T c for a given m A . Outlier materials clearly lying above the ∼ 1/ √ m A line include bismuthates and chloronitrates, suggesting the conventional electron-phonon mechanism is not driving superconductivity in these materials. Indeed, chloronitrates exhibit a very weak isotope effect [57], though some unconventional electron-phonon coupling could still be relevant for superconductivity [58]. Another important feature for low-T c materials is the average number of valence electrons. This recovers the empirical relation first discovered by Matthias more than sixty years ago [59]. Such findings validate the ability of ML approaches to discover meaningful patterns that encode true physical phenomena. Similar T c -vs.-predictor plots reveal more interesting and subtle features. A narrow cluster of materials with T c > 20 K emerges in the context of the mean covalent radii of compounds -another important predictor for low-T c superconductors. The cluster includes (leftto-right) alkali-doped C 60 , MgB 2 -related compounds, and bismuthates. The sector likely characterizes a region of strong covalent bonding and corresponding highfrequency phonon modes that enhance T c (however, frequencies that are too high become irrelevant for superconductivity). Another interesting relation appears in the context of the average number of d valence electrons. Figure 5c illustrates a fundamental bound on T c of all non-cuprate and non-iron-based superconductors.
A similar limit exists for cuprates based on the average number of unfilled orbitals (Figure 5d). It appears to be quite rigid -several data points found above it on inspection are actually incorrectly recorded entries in the database and were subsequently removed. The connection between T c and the average number of unfilled orbitals [60] may offer new insight into the mechanism for superconductivity in this family. Known trends include higher T c 's for structures that (i ) stabilize more than one superconducting Cu-O plane per unit cell and (ii ) add more polarizable cations such as Tl 3+ and Hg 2+ between these planes. The connection reflects these observations, since more copper and oxygen per formula unit leads to lower average number of unfilled orbitals (one for copper, two for oxygen). Further, the lower-T c cuprates typically consist of Cu 2− /Cu 3− -containing layers stabilized by the addition/substition of hard cations, such as Ba 2+ and La 3+ , respectively. These cations have a large number of unfilled orbitals, thus increasing the compound's average. Therefore, the ability of betweensheet cations to contribute charge to the Cu-O planes may be indeed quite important. The more polarizable the A cation, the more electron density it can contribute to the already strongly covalent Cu 2+ -O bond. Including AFLOW. The models described previously demonstrate surprising accuracy and predictive power, especially considering the difference between the relevant energy scales of most Magpie predictors (typically in the range of eV) and superconductivity (meV scale). This disparity, however, hinders the interpretability of the models, i.e., the ability to extract meaningful physical correlations. Thus, it is highly desirable to create accurate ML models with features based on measurable macroscopic properties of the actual compounds (e.g., crystallographic and electronic properties) rather than composite elemental predictors. Unfortunately, only a small subset of materials in SuperCon is also included in the ICSD: about 1, 500 compounds in total, only about 800 with finite T c , and even fewer are characterized with ab initio calculations. In fact, a good portion of known superconductors are disordered (off-stoichiometric) materials and notoriously challenging to address with DFT calculations. Currently, much faster and efficient methods are becoming available [40] for future applications.
To extract suitable features, data is incorporated from the AFLOW Online Repositories -a database of DFT calculations managed by the software package AFLOW. It contains information for the vast majority of compounds in the ICSD and about 550 superconducting materials. In Ref. 25, several ML models using a similar set of materials are presented. Though a classifier shows good accuracy, attempts to create a regression model for T c led to disappointing results. We verify that using Magpie predictors for the superconducting compounds in the ICSD also yields an unsatisfactory regression model. The 2. The most significant predictors and their importances for the three material-specific regression models. avg(x), std(x), max(x) and frac(x) denote the composition-weighted average, standard deviation, maximum, and fraction, respectively, taken over the elemental values for each compound. l 2 -norm of a composition is calculated by ||x||2 = issue is not the lack of compounds per se, as models created with randomly drawn subsets from SuperCon with similar counts of compounds perform much better. In fact, the problem is the chemical sparsity of superconductors in the ICSD, i.e., the dearth of closely-related compounds (usually created by chemical substitution). This translates to compound scatter in predictor space -a challenging learning environment for the model.
The chemical sparsity in ICSD superconductors is a significant hurdle, even when both sets of predictors (i.e., Magpie and AFLOW features) are combined via feature fusion. Additionally, this approach alone neglects the majority of the 16, 000 compounds available via Super-Con. Instead, we constructed separate models employing Magpie and AFLOW features, and then judiciously combined the results to improve model metrics -known as late or decision-level fusion. Specifically, two independent classification models are developed, one using the full SuperCon dataset and Magpie predictors, and another based on superconductors in the ICSD and AFLOW predictors. Such an approach can improve the recall, for example, in the case where we classify "high-T c " superconductors as those predicted by either model to be above-T sep . Indeed, this is the case here where, separately, the models obtain a recall of 40% and 66%, respectively, and together achieve a recall of about 76% (subject to small fluctuations due to variations in the test sets). In this way, the models' predictions complement each other in a constructive way such that above-T sep materials missed by one model (but not the other) are now accurately classified. Searching for new superconductors in the ICSD. As a final proof of concept demonstration, the classification and regression models described previously are integrated in one pipeline and employed to screen the en-tire ICSD database for candidate "high-T c " superconductors. (Note that "high-T c " is a simple label, the precise meaning of which can be adjusted.) Similar tools power high-throughput screening workflows for materials with desired thermal conductivity and magnetocaloric properties [51,61]. As a first step, the full set of Magpie predictors are generated for all compounds in SuperCon. A classification model similar to the one presented above is constructed, but trained only on materials in Super-Con and not in the ICSD (used as an independent test set). The model is then applied on the ICSD set to create a list of materials with predicted T c above 10 K. Opportunities for model benchmarking are limited to those materials both in the SuperCon and ICSD datasets, though this test set is shown to be problematic. The set includes  6. DOS of four compounds identified by the ML algorithm as potential materials with Tc > 20 K. The partial DOS contributions from s, p and d electrons and total DOS are shown in blue, green, red, and black, respectively. The large peak just below EF is a direct consequence of the flat band(s) present in all these materials. These images were generated automatically via AFLOW [43]. In the case of substantial overlap among k-point labels, the right-most label is offset below.
about 1, 500 compounds, though T c is reported for only about half of them. The model achieves an impressive accuracy of 0.98, which is overshadowed by the fact that 96.6% of these compounds belong to the T c < 10 K class. The precision, recall, and F 1 scores are about 0.74, 0.66, and 0.70, respectively. These metrics are lower than the estimates calculated for the general classification model, which is not unexpected given that this set cannot be considered randomly selected. Nevertheless, the performance suggests a good opportunity to identify new candidate superconductors. Next in the pipeline, the list is fed into a random forest regression model (trained on the entire SuperCon database) to predict T c . Filtering on the materials with T c > 20 K, the list is further reduced to about 2, 000 compounds. This count may appear daunting, but should be compared with the total number of compounds in the database -about 110, 000. Thus, the method selects less than two percent of all materials, which in the context of the training set (containing more than 20% with "high-T c "), suggests that the model is not overly biased toward predicting high critical temperatures.
The vast majority of the compounds identified as candidate superconductors are cuprates, or at least compounds that contain copper and oxygen. There are also some materials clearly related to the iron-based superconductors. The remaining set has 35 members, and is composed of materials that are not obviously connected to any high-temperature superconducting families (see Table 3) [62]. None of them is predicted to have T c in excess of 40 K, which is not surprising, given that no such instances exist in the training dataset. All contain oxygen -also not a surprising result, since the group of known superconductors with T c > 20 K is dominated by oxides.
The list comprises several distinct groups. Especially interesting are the compounds containing heavy metals (such as Au, Ir, Ru), metalloids (Se, Te), and heavier post-transition metals (Bi, Tl), which are or could be pushed into interesting/unstable oxidation states. Charge doping and/or pressure may be needed to drive these materials into a superconducting state. The most surprising and non-intuitive of the compounds in the list are the silicates and the germanates. These materials form corner-sharing SiO 4 or GeO 4 polyhedra, not unlike quartz glass, and also have counter cations with full or empty shells such as Cd 2 + or K + . Converting these insulators to metals (and possibly superconductors) likely requires significant charge doping. However, the similarity between these compounds and cuprates is meaningful. In compounds like K 2 CdSiO 4 or K 2 ZnSiO 4 , K 2 Cd (or K 2 Zn) unit carries a 4+ charge that offsets the (SiO 4 ) 4− (or (GeO 4 ) 4− ) charges. This is reminiscent of the way Sr 2 balances the (CuO 4 ) 4− unit in Sr 2 CuO 4 . Such chemical similarities based on charge balancing and stoichiometry were likely identified and exploited by the ML algorithms. The electronic properties calculated by AFLOW offer additional insight into the results of the search, and suggest a possible connection among these candidate. Plotting the electronic structure of the potential superconductors exposes an extremely peculiar feature shared by almost all -one or several (nearly) flat bands just below the energy of the highest occupied electronic state. Such bands lead to a large peak in the DOS (see Figure 6) and can cause a significant enhancement in T c . Peaks in the DOS elicited by van Hove singularities can enhance T c if sufficiently close to E F [63][64][65]. However, note that unlike typical van Hove points, a true flat band creates divergence in the DOS (as opposed to its derivatives), which in turn leads to a critical temperature dependence linear in the pairing interaction strength, rather than the usual exponential relationship yielding lower T c [33]. Additionally, there is significant similarity with the band structure and DOS of layered BiS 2 -based superconductors [66].
This band structure feature came as the surprising result of applying the ML model. It was not sought for, and, moreover, no explicit information about the electronic band structure has been included in these predictors. This is in contrast to the algorithm presented in Ref. 30, which was specifically designed to filter ICSD compounds based on several preselected electronic structure features.
While at the moment it is not clear if some (or indeed any) of these compounds are really superconducting, let alone with T c 's above 20 K, the presence of this highly unusual electronic structure feature is encouraging. Attempts to synthesize several of these compounds are already underway.

DISCUSSION
Herein, several machine learning tools are developed to study the critical temperature of superconductors. Based on information from the SuperCon database, initial coarse-grained chemical features are generated using the Magpie software. As a first application of ML methods, materials are divided into two classes depending on whether T c is above or below 10 K. A non-parametric random forest classification model is constructed to predict the class of superconductors. The classifier shows excellent performance, with out-of-sample accuracy and F 1 score of about 92%. Next, several successful random forest regression models are created to predict the value of T c , including separate models for three material sub-groups, i.e., cuprate, iron-based, and "low-T c " compounds. By studying the importance of predictors for each family of superconductors, insights are obtained about the physical mechanisms driving superconductivity among the different groups. With the incorporation of crystallographic-/electronic-based features from the AFLOW Online Repositories, the ML models are further improved. Finally, we combined these models into one integrated pipeline, which is employed to search the entire ICSD database for new inorganic superconductors. The model identified about 30 oxides as candidate materials. Some of these are chemically and structurally similar to cuprates (even though no explicit structural information was provided during training of the model). Another feature that unites almost all of these materials is the presence of flat or nearly-flat bands just below the energy of the highest occupied electronic state.
In conclusion, this work demonstrates the important role ML models can play in superconductivity research. Records collected over several decades in SuperCon and other relevant databases can be consumed by ML models, generating insights and promoting better understanding of the connection between materials' chemistry/structure and superconductivity. Application of sophisticated ML algorithms has the potential to dramatically accelerate the search for candidate high-temperature superconductors.

METHODS
Superconductivity data. The SuperCon database consists of two separate subsets: "Oxide & Metallic" (inorganic materials containing metals, alloys, cuprate high-temperature superconductors, etc.) and "Organic" (organic superconductors). Downloading the entire inorganic materials dataset and removing compounds with incompletely-specified chemical compositions leaves about 22, 000 entries. In the case of multiple records for the same material, the reported material's Tc's are averaged, but only if their standard deviation is less than 5 K, and discarded otherwise. This brings the total down to about 16, 400 compounds, of which around 4, 000 have no critical temperature reported. Each entry in the set contains fields for the chemical composition, Tc, structure, and a journal reference to the information source. Here, structural information is ignored as it is not always available.
There are occasional problems with the validity and consistency of some of the data. For example, the database includes some reports based on tenuous experimental evidence and only indirect signatures of superconductivity, as well as reports of inhomogeneous (surface, interfacial) and nonequilibrium phases. Even in cases of bona fide bulk superconducting phases, important relevant variables like pressure are not recorded. Though some of the obviously erroneous records were removed from the data, these issues were largely ignored assuming their effect on the entire dataset to be relatively modest. The data cleaning and processing is carried out using the Python Pandas package for data analysis [67]. Chemical and structural features. The predictors are calculated using the Magpie software [68]. It computes a set of 145 attributes for each material, including: (i) stoichiometric features (depends only on the ratio of elements and not the specific species); (ii) elemental property statistics: the mean, mean absolute deviation, range, minimum, maximum, and mode of 22 different elemental properties (e.g., period/group on the periodic table, atomic number, atomic radii, melting temperature); (iii) electronic structure attributes: the average fraction of electrons from the s, p, d and f valence shells among all elements present; and (iv ) ionic compound features that include whether it is possible to form an ionic compound assuming all elements exhibit a single oxidation state.
ML models are also constructed with the superconducting materials in the AFLOW Online Repositories. AFLOW is a high-throughput ab initio framework that manages density functional theory (DFT) calculations in accordance with the AFLOW Standard [21]. The Standard ensures that the calculations and derived properties are empirical (reproducible), reasonably well-converged, and above all, consistent (fixed set of parameters), a particularly attractive feature for ML modeling. Many materials properties important for superconductivity have been calculated within the AFLOW framework, and are easily accessible through the AFLOW Online Repositories. The features are built with the following properties: number of atoms, space group, density, volume, energy per atom, electronic entropy per atom, valence of the cell, scintillation attenuation length, the ratios of the unit cell's dimensions, and Bader charges and volumes. For the Bader charges and volumes (vectors), the following statistics are calculated and incorporated: the maximum, minimum, average, standard deviation, and range.
Machine learning algorithms. Once we have a list of relevant predictors, various ML models can be applied to the data [52,53]. All ML algorithms in this work are variants of the random forest method [54]. It is based on creating a set of individual decision trees (hence the "forest"), each built to solve the same classification/regression problem. The model then combines their results, either by voting or averaging depending on the problem. The deeper individual tree are, the more complex the relationships the model can learn, but also the greater the danger of overfitting, i.e., learning some irrelevant information or just "noise". To make the forest more robust to overfitting, individual trees in the ensemble are built from samples drawn with replacement (a bootstrap sample) from the training set. In addition, when splitting a node during the construction of a tree, the model chooses the best split of the data only considering a random subset of the features.
The random forest models above are developed using scikitlearn -a powerful and efficient machine learning Python library [69]. Hyperparameters of these models include the number of trees in the forest, the maximum depth of each tree, the minimum number of samples required to split an internal node, and the number of features to consider when looking for the best split. To optimize the classifier and the combined/family-specific regressors, the GridSearch function in scikit-learn is employed, which generates and compares candidate models from a grid of parameter values. To reduce computational expense, models are not optimized at each step of the backward feature selection process.
Prediction errors of the regression models. Previously, several regression models were described, each one designed to predict the critical temperatures of materials from different superconducting groups. These models achieved an impressive R 2 score, demonstrating good predictive power for each group. However, it is also important to consider the accuracy of the predictions for individual compounds (rather than on the aggregate set), especially in the context of searching for new materials. To do this, we calculate the prediction errors for about 300 materials from a test set. Specifically, we consider the difference between the logarithm of the predicted and measured critical temperature [ln(T meas c )−ln(T pred c )] normalized by the value of ln(T meas c ) (normalization compensates the different Tc ranges of different groups). The models show comparable spread of errors. The histograms of errors for the four models (combined and three group-specific) are shown in Fig. 7. The errors approximately follow a normal distribution, centered not at zero but at a small negative value. This suggests the models are marginally biased, and on average tend to slightly underestimate Tc. The variance is comparable for all models, but largest for the model trained and tested on iron-based materials, which also shows the smallest R 2 . Performance of this model is expected to benefit from a larger training set.