Functional Form of the Superconducting Critical Temperature from Machine Learning

TL;DR

This paper uses machine learning to improve the analytical prediction of superconducting critical temperatures, revealing limitations of existing models and suggesting new descriptors for better accuracy in superconductor discovery.

Contribution

It introduces a machine learning approach that refines analytical formulas for $T_c$, highlighting the need for new descriptors beyond traditional electron-phonon parameters.

Findings

Improved analytical expression for $T_c$ over Allen-Dynes fit.

Reasonable prediction for high-pressure H$_3$S superconductor.

Identification of the need for Fermi surface descriptors.

Abstract

Predicting the critical temperature $T_c$ of new superconductors is a notoriously difficult task, even for electron-phonon paired superconductors for which the theory is relatively well understood. Early attempts by McMillan and Allen and Dynes to improve on the weak-coupling BCS formula led to closed-form approximate relations between $T_c$ and various measures of the phonon spectrum and the electron-phonon interaction appearing in Eliashberg theory. Here we propose that these approaches can be improved with the use of machine learning algorithms. As an initial test, we train a model for identifying low-dimensional descriptors using the $T_c < 10$ K data tested by Allen and Dynes, and show that a simple analytical expression thus obtained improves upon the Allen-Dynes fit. Furthermore, the prediction for the recently discovered high $T_c$ material H$_3$S at high pressure is quite reasonable. Interestingly, $T_c$'s for more recently discovered superconducting systems with a more two-dimensional electron-phonon coupling, which do not follow Allen and Dynes' expression, also do not follow our analytic expression. Thus, this machine learning approach appears to be a powerful method for highlighting the need for a new descriptor beyond those used by Allen and Dynes to describe their set of isotropic electron-phonon coupled superconductors. We argue that this machine learning method, and its implied need for a descriptor characterizing Fermi surface properties, represents a promising new approach to superconductor materials discovery which may eventually replace the serendipitous discovery paradigm begun by Kamerlingh Onnes.

Figures

Equations

$$T_c\simeq 1.14{\omega }_D{e}^{-\frac{1}{\lambda }},$$

$$T_c\simeq \frac{{\omega }_D}{1.45}\exp \left(-\frac{1.04(1+\lambda )}{\lambda -{\mu }^{\ast }(1+0.62\lambda )}\right),$$

$$T_c=\frac{f_1{f}_2{\omega }_{\log }}{1.20}\exp \left(-\frac{1.04(1+\lambda )}{\lambda -{\mu }^{\ast }(1+0.62\lambda )}\right),$$

$${T_c}^{\mathrm{SISSO}}=0.09525\frac{{\lambda }^4{\omega }_{\log }}{{\lambda }^3+\sqrt{{\mu }^{\ast }}}.$$

Full text

I Introduction

Discovery of new superconductors has historically proceeded largely serendipitously, with guidance from rules of thumb (such as Matthias’ e/a ratio) rather than many-body and ab-initio theory. The space of possible materials to search for new superconductors is vast, considering that many discoveries in the last thirty years are multinary compounds. Thus, it is desirable to appeal to recent computational developments, aided by theory, to assist this process. The history of ab-initio and materials-genome type approaches to superconducting materials discovery has recently been reviewed by Norman Norman (2016), Pickett Pickett (2017), and Duan et alDuan et al. (2019).

While initially, success in prediction (as opposed to analysis after discovery, i.e., postdiction) was rare to nonexistent, more recently the potential for theory to aid in the discovery of new high-temperature superconductors was dramatically demonstrated by the prediction and subsequent discovery, in 2015, of superconductivity at $T_{c}=200$ K in H3S at one $\sim$150\text{,}\mathrm{G}\mathrm{P}\mathrm{a}$$ pressure Drozdov et al. (2015). This experiment shattered the assumed ceiling for $T_{c}$ in electron-phonon superconductors and was followed by the recent discovery of superconductivity in compressed lanthanum hydride at 215-250 K Drozdov et al. (2018); Somayazulu et al. (2019), also preceded by a theoretical prediction Liu et al. (2017); Peng et al. (2017). Recent computational approaches to hydride superconductivity have been reviewed in Refs. Boeri, 2019; Pickett and Eremets, 2019.

Despite these undeniable successes and the demonstration that the old assumed limit of 35-40 K for $T_{c}$ due to the exchange of phonons is incorrect, these experiments do not provide a clear strategy to optimize $T_{c}$ in the vast phase space of materials. This is at least partially due to an inability to identify the correct materials descriptors, parameters directly reflecting the underlying mechanism of superconductivity. For some classes of materials, e.g., thermoelectrics, considerable progress has been made in high-throughput approaches identifying simple observables recorded in databases that contribute to a material’s figure of merit Curtarolo et al. (2013). For superconductivity, however, such approachesStanev et al. (2018) are considerably more difficult, both because the theory is more complex, and the figure of merit, $T_{c}$, depends extremely sensitively on the underlying interactions.

This last difficulty is clear already from the Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity Bardeen et al. (1957), among whose great successes was the proof that for weak attractive interactions, fermions pair with an instability that corresponds to an essential singularity in the dimensionless coupling constant $\lambda$, leading to the well-known expression,

[TABLE]

where $\omega_{D}$ is the Debye frequency. BCS theory is successful because it predicts superconducting properties accurately in terms of measured $T_{c}$’s, but the essential singularity alone suggests that accurate calculations will be difficult. Besides, Eq. (1) is strictly valid only in the weak coupling limit $\lambda\ll 1$ and if the Coulomb interaction is neglected.

The inadequacy of the BCS expression for $T_{c}$ was already clear by the late 1960’s, when McMillan McMillan (1968) introduced an improved formula based on Eliashberg theory Eliashberg (1960), relating $T_{c}$ to a small number of physical quantities calculated from the effective electron-phonon interaction $\alpha^{2}F(\omega)$ that could in principle be extracted from tunneling data McMillan and Rowell (1965),

[TABLE]

where $\mu^{\ast}$ is the Coulomb pseudopotential. This approximate formula led McMillan to predict a maximum $T_{c}$ in a given class of materials. Dynes Dynes (1972) later replaced the prefactor $\omega_{D}$/1.45 of the McMillan equation (2) with $\langle\omega\rangle/1.20$, where $\langle\omega\rangle$ is the first moment of the distribution $g(\omega)=2/(\lambda\omega)\alpha^{2}F(\omega)$.

Based on a reanalysis of Eliashberg theory and newly available computational checks in special cases, Allen and Dynes Allen and Dynes (1975) proposed an alternate approximate formula,

[TABLE]

where $f_{1}$ and $f_{2}$ are correction factors that depend on $\lambda,\mu^{\ast},\omega_{\log}$, and $\bar{\omega}_{2}$, where $\bar{\omega}_{2}$ is the second moment of the phonon density of states $F(\omega)$. They showed that the expression (3) fit the $T_{c}$ of a variety of superconductors known at the time, using data derived from tunneling, and that it implied the absence of any maximum $T_{c}$, except that caused by the competition between $\lambda$ and $\omega_{\log}\equiv\exp\langle\ln\omega\rangle$, where the average is taken over $g(\omega)$. Unlike the McMillan expression, which saturates to a constant value as $\lambda\rightarrow\infty$, the Allen-Dynes equation obeys an asymptotic result of Eliashberg theory, that $T_{c}\sim\sqrt{\lambda}$ as $\lambda\rightarrow\infty$ with other parameters fixed.

The Allen-Dynes equation has played a crucial role in the discussion of high-temperature superconductivity and indeed is often used to extract quoted values of $\lambda$ in the literature for materials where tunneling data is not available. Nevertheless, it is important to recall that it has been derived from Eliashberg theory, which itself is implemented with various approximations, e.g., the momentum dependence of the electron-phonon interaction was often neglected in early studies. The full evaluation of the Eliashberg equation is computationally expensive and not currently suitable for high-throughput superconductor discovery. It would be highly desirable to develop an expression for $T_{c}$ that generalizes the Allen-Dynes equation and is applicable over an extensive range of parameters to guide such searches.

In this letter, we use modern machine learning techniques to critically examine the Allen-Dynes equation in the context of recently discovered materials. These techniques are analytical in nature, meaning they search for analytical relations between a minimal set of features, i.e., physical parameters, and the desired properties. Specifically, we apply the Sure-Independence Screening and Sparsifying Operator (SISSO) method Ouyang et al. (2018) to estimate $T_{c}$ from $\lambda$, $\mu^{\ast}$, and $\omega_{\log}$ with the goal to obtain an equation of similar or enhanced performance to the one proposed by Allen and Dynes Allen and Dynes (1975). We find that we can improve on the Allen-Dynes fit to strong-coupling superconductors, with a smaller set of descriptors. More interestingly, the approach identifies outliers like MgB2, $T_{c}$=39 K, which suggests the importance of new physics essential to high $T_{c}$ that needs to be incorporated in an improved formula to guide the search for new electron-phonon superconductors in materials space.

II Methods

To generate models for predicting $T_{c}$, we apply recently developed methods of equation-based machine learning, subject to physical constraints. In the SISSO approach, the predictive models are expressed as analytical formulas relating physical quantities with algebraic operations such as addition and exponentiation. Given a tabulated set of scalar-valued physical quantities, or features, the SISSO method constructs additional features by iteratively applying operations from a specified set, e.g., +, $\times$, exp, $\sqrt{~{}}$, 2.

To pinpoint the best equations, the SISSO method employs the sure-independence screening (SIS) method and the sparse-solution algorithm using sparsifying operators (SO) in tandem. After constructing the feature space, the SIS method selects a subspace of features with the largest linear correlation with the target property ($T_{c}$), i.e., the largest absolute value of their dot product. The SO step then evaluates all possible combinations of features from the SIS subspace, yielding the optimal least-squares solution and residual. With such a vast feature space, the combinatorial optimization in each SO step relies on $L_{0}$ regularization, which penalizes the number of non-zero coefficients. Combined with one numerical prefactor, fit from available data, each feature is used to generate one predictive model.

We benchmark the performance of different models identified by SISSO using leave-one-out cross-validation. Given $N$ available data points, each model is repeatedly fit using $N-1$ points and evaluated with the excluded point. The average evaluation error across $N$ iterations, where each point is tested once, is the leave-one-out cross-validation error. This method can help to maximize the transferability of a model by reducing “overfitting”, i.e., models that exhibit low root-mean-square error in predictions on the training data but very high root-mean-square error in the testing data.

We apply the SISSO method to estimate $T_{c}$ from $\lambda$, $\mu^{\ast}$, and $\omega_{\log}$ to obtain an equation of similar performance to the one proposed by Allen and Dynes Allen and Dynes (1975). We use the values of $\lambda$, $\mu^{\ast}$, and $\omega_{\log}$, and the target property, $T_{c}$, from the data for 29 superconducting materials provided by Allen and Dynes (Table I in Ref. Allen and Dynes (1975)). Next, we apply the SISSO method with 4 binary operators (+, $-$, $\times$, /) and 7 unary operators (exp, log, $\sqrt{}$, $\sqrt[3]{}$, -1, 2, 3) three times to generate 3,414,094 features. Fig. 1 shows the rapid growth of the feature space with the number of iterations. Of the initial feature space, we select the equations with the highest linear correlation to $T_{c}$ using sure-independence screening with a minimum correlation magnitude (inner product) of 0.5. To further reduce the number of features and eliminate unphysical equations, we apply constraints. We select equations that are linearly proportional to $\omega_{\log}$ and obey the proper $\lambda\rightarrow 0$ limiting behavior. Additionally, we filter for equations that are strictly positive, real, finite, continuous, and monotonic across the relevant training and testing feature spaces. To evaluate the generalizability and performance of these equations, we compute the error against a testing set of 9 superconductors Arnold et al. (1980); Mitrović et al. (1984); Kim et al. (2006); Kwo and Geballe (1981); Bending et al. (1987); Kihlstrom et al. (1984); Nishio et al. (1993); Wolf and Noer (1979); Kihlstrom et al. (1985), shown in green in Fig. 2, using leave-one-out cross validation.

III Results

Fig. 2 illustrates the main proof-of-principle result that machine learning can provide an analytic equation of similar performance to the Allen-Dynes equation. The equation-based machine learning uses the values of $\lambda,\omega_{\log}$, and $\mu^{\ast}$ of the 29 materials in Table I of Allen and Dynes Allen and Dynes (1975), and neglects the average frequency $\bar{\omega}_{2}\equiv\langle\omega^{2}\rangle^{1/2}$ taken over the $g(\omega)$ distribution that is also used in the Allen-Dynes equation. The SISSO method and subsequent physical constraints lead to the optimal equation,

[TABLE]

Importantly, Eq. (4) emerged from our approach with the smallest RMSE even before any of the physical constraints summarized in Fig. 1 were applied. Fig. 2(a) compares the performance of this equation with the modified McMillan and Allen-Dynes equations for the measured $T_{c}$’s of the 29 materials that train the model. The root-mean-square error (RMSE) of this equation evaluated on the training data is 0.25 K, significantly smaller than the RMSE of 0.92 K for the modified McMillan equation, and also slightly lower than the RMSE of 0.30 K for the Allen-Dynes equation. This result is impressive given the use of only 3 parameters and a single numerical coefficient compared to 3 parameters and 4 coefficients for the modified McMillan and 4 parameters and 7 coefficients for the Allen-Dynes equation.

Figure 2(b) shows the testing of Eq. (4) for a variety of other superconductors, mostly of higher $T_{c}$. Because $\mu^{\ast}$ data were not available for these materials, we adopt a constant value of $\mu^{\ast}=0.1$. This procedure introduces some unknown error into the analysis, but despite this, the fit to the new materials is rather good, with an RMSE of only 3.2 K (9.1%).

It is important to note that Eq. (4) is not derived from any physical theory and therefore may contain some terms that may make no physical sense, e.g., the appearance of the $\sqrt{\mu^{\ast}}$ term, which may be a proxy for a constant term due to the small range of data and the paucity of features at this level of learning. The limit $T_{c}\rightarrow 0$ as $\lambda\rightarrow 0$ in Eq. (4) even at nonzero $\mu^{*}$ may reflect the lack of data at small coupling. Also, Eq. (4) increases monotonically with $\lambda$, with linear behavior at very high couplings. This behavior violates the asymptotic limit of Eliashberg theory, $T_{c}\sim\sqrt{\lambda}$, built into the Allen-Dynes equation Allen and Dynes (1975). Again, this disagreement with physics is due to the absence of data points, either in the training or the testing set, which deviate significantly from the linear behavior predicted by Eq. (4).

Fig. 3 shows the functional behavior $T_{c}(\lambda)$ of the 100 highest-scored equations discovered by SISSO; it is clear that almost all of these equations are equally valid over the range of $\lambda$ values where data exist. This highlights the need for measurements to determine the materials parameters $\lambda$, $\omega_{\log}$, and $\mu^{\ast}$ reliably for both very low $T_{c}$ materials, as well as for some of the recently discovered higher-$T_{c}$ systems.

Fig. 2(b) also shows some dramatic failures of the learned equation, namely for MgB2 and NbS2. The probable reasons for these failures are both revealing and reassuring. The point labeled MgB2(1) with a predicted $T_{c}$ of 10 K is one where $\omega_{\log}$, a logarithmic average of the electron-phonon interaction function $\alpha^{2}F/\omega$, was determined from a specific heat measurement of the Debye frequency $\omega_{D}$, which depends only on the phonon density of states $F(\omega)$. Relating the Debye frequency with $\omega_{\log}$ neglects the difference between the two distributions Allen and Dynes (1975). This assumption is particularly poor in MgB2, where high-frequency phonons couple anomalously strongly. In addition, $\lambda$ was determined from standard expressions for the high-temperature resistivity of a 3D metal. It is well known that MgB2 has strong 2D character, and that the full momentum and band dependence of the Eliashberg function $\lambda_{n{\bf k},n{\bf k^{\prime}}}$ must be accounted for to obtain reasonable values for $T_{c}$ from first principles Margine and Giustino (2013). It is interesting to note that if one uses the higher value of $\lambda$ obtained from Ref. Margine and Giustino, 2013 in Eq. (4), one obtains data point MgB2(2), with the significantly enhanced predicted $T_{c}$ of 20 K, but still far from the measured value of 40 K and even further from the full Eliashberg calculation of 50 K Margine and Giustino (2013).

These discrepancies indicate, not surprisingly, that a machine trained on a database of nearly isotropic low-$T_{c}$ superconductors cannot capture the physics of highly anisotropic higher-$T_{c}$ materials using the simple averaged descriptors chosen by Allen and Dynes. The same principle apparently applies to NbS2, which while having a low-$T_{c}$ is quite 2-dimensional. However, Eq. (4) may have significant predictive power extrapolated to higher-$T_{c}$ 3D systems. For example, if we take values of $\lambda$ and $\omega_{\log}$ calculated from first principles for H3S at 140 GPa pressure Duan et al. (2014); Komelj and Krakauer (2015), the predicted $T_{c}$ from Eq. (4) is 262 K, compared to the measured value of 203 K. This result is similar to the result obtained from Allen-Dynes using calculated parameters, but substantially higher than the modified McMillan equation used in Refs. Duan et al., 2014; Komelj and Krakauer, 2015.

IV Conclusion

We have demonstrated that machine learning can discover equations that describe the dependence of superconducting $T_{c}$’s on moments of distributions of phonon frequencies and electron-phonon couplings, as used originally by Allen and Dynes in their attempt to understand the systematics of $T_{c}$ in the framework of Eliashberg theory. While the method is quite successful in predicting known superconductors of the same general type as the original Allen-Dynes dataset, with fewer parameters, the existence of a few anomalous outliers suggests that the use of such methods for high-throughput materials discovery will require new descriptors that capture anomalous features, e.g., the anisotropy of the electron-phonon interactions and unusual electronic states that take advantage of them. A natural modern extension of the philosophy of Allen and Dynes is then to calculate from first principles a few key measures of electronic structure crucial for superconductivity, together with the moments discussed above, and apply machine-learning methods as described here. We anticipate that this approach will allow a much more efficient and thorough investigation of materials space than current approaches that rely on fully anisotropic Eliashberg calculations for each material.

V Acknowledgements

We are grateful to L. Boeri, P. Allen, and W. Pickett for valuable discussions. The work by SRX and RGH was supported by the National Science Foundation under grants Nos. PHY-1549132 and DMR-1609306. The work by GRS was performed under the auspices of the Office of Sciences, United States Department of Energy, contract numbers DE-FG02-86ER45268.