Machine learning prediction of superconducting pairing potential for quasi-one-dimensional disordered s-wave superconductors

Superconductivity in quasi-one-dimensional (Q1D) systems represents a fascinating frontier in condensed matter physics, where the interplay of quantum confinement, strong electron correlations, and disorder leads to a rich variety of physical phenomena markedly different from bulk superconductors[1-4]. The reduced dimensionality profoundly affects both single-particle and collective excitations, enhancing quantum fluctuations and making these systems particularly sensitive to perturbations. In clean 1D systems, the superconducting state is already precarious due to enhanced phase fluctuations, while the introduction of disorder creates an even more complex landscape where superconductivity can be locally suppressed or, counterintuitively, enhanced[5, 6]. When coupled with a system of higher dimensionality, Q1D systems exhibit a notable reduction in fluctuations[7, 8]. This phenomenon may facilitate the attainment of elevated critical temperatures, typical of Q1D systems.

The theoretical description of such systems typically relies on the Bogoliubov-de Gennes (BdG) formalism[9], which provides a self-consistent mean-field framework for treating inhomogeneous superconductivity. The BdG equations couple the superconducting pairing potential to the underlying electronic structure, requiring iterative numerical solutions that become increasingly demanding as system size grows. For disordered systems, this challenge is compounded by the need to examine many realizations of disorder to obtain statistically meaningful results. Traditional approaches often face severe computational limitations when addressing key questions about the spatial structure of the pairing potential, its statistical distribution, or the nature of the superconductor-insulator transition. These limitations become particularly acute when studying mesoscopic systems or when attempting to map out phase diagrams as functions of multiple parameters.Recent years have witnessed remarkable progress in applying machine learning (ML) techniques to problems in quantum many-body physics[10-13]. Neural networks and other ML architectures have demonstrated exceptional capability in learning complex mappings between system parameters and physical observables, often achieving accuracy comparable to conventional numerical methods while offering superior computational efficiency. In the context of disordered superconductors, ML approaches present several unique advantages: they can learn the underlying physical principles from training data and generalize them to new configurations without solving the full quantum problem; they enable rapid exploration of large parameter spaces that would be prohibitive for traditional methods; and they can identify and characterize rare but physically important configurations that might be missed in conventional sampling. As one of the most successful ML approaches, neural networks (NNs) are capable of addressing the complex non-linear correlation between disorder potentials and the emerging superconducting pairing potential while maintaining computational tractability for statistical studies of disorder effects. This combination of physical accuracy and scalability makes NNs the optimal choice for this class of quantum many-body problems.

This work presents a neural network approach for predicting solutions to the BdG equations for Q1D disordered superconductors. By learning from self-consistent solutions across various disorder configurations, the network successfully captures the complex relationship between disorder potential and superconducting pairing potential and reproduces both the average suppression of superconductivity and the formation of localized superconducting regions that characterize these systems.

A disordered s-wave superconductor is described by the following Hamiltonian[14-19]

Here \( \hat{c}_{i\sigma }^{\dagger } (\hat{c}^{\,}_{i\sigma }) \) creates (annihilates) an electron with spin \( \sigma \) at site \( i \) of the lattice; \( \hat{n}^{\,}_{i\sigma }=\hat{c}_{i\sigma }^{\dagger }\hat{c}^{\,}_{i\sigma } \); \( t^{\,}_{ij}=t \) is the tunneling amplitude between neighboring sites; \( \Delta^{\,}_i = - g \langle \hat{c}^{\,}_{i\uparrow } \hat{c}^{\,}_{i\downarrow} \rangle \) is the pairing potential with \( g > 0 \) being the on-site coupling constant; \( V_i \) is the disorder potential introduced at each lattice site and \( \mu^{\,}_i = \mu^{\,}_0 - g \langle n^{\,}_i\rangle /2 \) is the chemical potential with the site-dependent Hartree shift. We consider a quasi-one-dimensional lattice of size \( L \) with periodic boundary conditions, where quantum fluctuations are strongly suppressed[7, 8]. In what follows, all energy values are represented in units of the hopping amplitude \( t \), while distances are measured in terms of the lattice spacing \( a \).

The Hamiltonian (1) may be also considered within the BdG mean-field approach[9, 20] and the corresponding BdG equations are written as

where the matrices of operators \( \hat{H}^{\,}_0 \) and \( \hat{\Delta} \) have the size of \( L\times L \) and the following matrix elements,

The solution of the BdG equations (2) needs to be found in a self-consistent manner with the conditions for the pairing potential \( \Delta_i \) and the Hartree potential \( U_{i} \) at each lattice site \( i \)[14, 19],

At every lattice site, the disorder potential \( V_i \) is randomly sampled from a uniform distribution spanning \( [-V_0, V_0] \). Here, \( V_0 \) denotes the disorder strength, which is systematically varied over the range from zero to 0.75 in increments of 0.05.

In numerical calculations, we choose a coupling strength of \( g=1.5 \), which places the system in the strong-coupling regime with \( \Delta/\varepsilon^{\,}_F \approx 0.1 \), with \( \Delta \) being the pairing potential averaged over the system and \( \varepsilon^{\,}_F \) the Fermi energy. The chemical potential \( \mu \) is tuned to achieve an electron filling factor of \( n=0.85 \), avoiding the half-filling and related symmetry effects. The calculations are made for \( L=256 \) to prevent any finite-size effects for the chosen parameters.

The neural network architecture, illustrated in Figure 1c, consists of four fully connected hidden layers with 100 neurons each, utilizing ReLU activation functions. This configuration results in a model containing over 10⁴ trainable parameters, providing sufficient complexity to capture the complex relationship between disorder configurations and the superconducting pairing potential.

Figure 1

The histograms of disorder (panel a) and pairing potential (panel b) distributions used in the training dataset. The sketch of the neural network architecture is shown in panel c.

For training data generation, we performed self-con-sistent solutions of the BdG equations (2) for more than 200 distinct disorder realizations at each disorder strength \( V_0 \). The network's input for predicting the pairing potential \( \Delta_i \) at site \( i \) includes the disorder potential \( V_j \) for all sites \( j \) within the region \( |j-i| \leq R \), where \( R \) is a cut-off value to be chosen in accordance with the system's coherence length \( \xi \approx 6 \). This cutoff ensures the inclusion of all relevant spatial correlations while optimizing computational efficiency, as disorder effects beyond this range become negligible for the local pairing potential. We emphasize that the method's applicability is limited to regimes where the superconducting coherence length \( \xi \) remains finite. In systems with intrinsically large \( \xi \) or near the critical temperature \( T^{\,}_c \), the computational cost of generating the training dataset through exact BdG calculations becomes prohibitive. A promising alternative pathway for such regimes could involve using spatially-resolved experimental data[21], as the ground truth for training, thereby circumventing the need for full-scale theoretical simulations.

The training dataset statistics, shown in figure Figure 1(panels a and b), reveal important physical characteristics. The disorder potential histogram Figure 1a) follows the expected triangular distribution, while the pairing potential distribution (Figure 1b) exhibits a bimodal structure. The left broad peak corresponds to regions of suppressed superconductivity near the superconductor-insulator transition. The sharp peak on the right represents the pairing potential value for a clean system, and the high-energy tail indicates rare disorder-induced enhancements of superconductivity. The relative scarcity of these enhancement events suggests they may contribute disproportionately to the prediction error.

We implemented the Adam optimizer[22] with a scaling learning rate policy, varying the rate from 10^-3 to 10^-6 with early stopping during training. The optimization process continued until the mean-squared error (MSE) loss function value fell below 10^-3, ensuring robust convergence. This architecture maintains full physical interpretability while offering exceptional scalability to large systems with identical coupling and filling parameters, making it particularly suitable for studies of extended disordered superconductors. While the present study demonstrates the feasibility of neural network-based prediction using conventional digital hardware, significant further acceleration could be achieved through analog neuromorphic computing[23].

We evaluate the performance of the trained NN on a test dataset consisting of over 20 distinct disorder configurations for each disorder strength \( V^{\,}_0 \), totaling more than 50000 target values of the pairing potential \( \Delta^{\,}_i \). Crucially, these configurations were not included in the training process, ensuring an unbiased assessment of the model's predictive capabilities.

Figure 2

The MSE loss function for the test dataset as a function of the cut-off value \( R \).

To achieve precise modeling of how disorder affects the pair potential without sacrificing computational efficiency, it was essential to determine the selected cut-off parameter \( R \), which specifies the number of neighboring sites to consider for predictions at each site. Error analysis of the test dataset presented in Figure 2 shows that the prediction accuracy improves with increasing \( R \) until saturating below the target threshold 10^-3 for \( R \geq 8 \), which is consistent with the assumption that the effect of disorder at distances greater than \( \xi \) is insignificant. The results demonstrate a clear trade-off: while larger \( R \)-values should reduce the prediction error by incorporating more distant disorder potentials, they also exponentially increase the required training data size to properly sample the expanded configuration space. Through this analysis, we identified \( R = 8 \) as the optimal truncation value to reliably reproduce effects of disorder on superconductivity while remaining computationally tractable for large-scale statistical studies.

Figure 3

Average pairing potential prediction by NN for various disorder strengths. Red dots show values of numerical self-consistent solution and blue line corresponds to the NN predictions.

The primary focus of this analysis is the model's ability to reproduce the disorder-induced suppression of superconductivity — a phenomenon that lies beyond the scope of the Anderson theorem. Figure 3 demonstrates the evolution of the spatially averaged pairing potential \( \bar{\Delta} \) as a function of disorder strength \( V^{\,}_0 \). The self-consistent BdG calculations reveal a systematic suppression of \( \bar{\Delta} \) with increasing \( V_0 \), reflecting the gradual destruction of superconducting order. Remarkably, the neural network predictions exhibit perfect agreement with the BdG results, quantitatively capturing both the initial weakening of superconductivity at moderate disorder strengths and its slow progression toward a superconductor-insulator transition with the increasing \( V^{\,}_0 \).

We now examine in detail the NN's capacity to reproduce the spatial variations of the pairing potential \( \Delta^{\,}_i \) in disordered systems. Figure 4 presents characteristic spatial profiles of \( \Delta^{\,}_i \) and their statistical distributions for various disorder strengths \( V^{\,}_0 \). The physics of disordered superconductors predicts that the pairing potential distribution should evolve from a narrow peak in a clean system to a broad, asymmetric form in the presence of disorder. Our results confirm this expectation: the distribution becomes lognormal as disorder grows, with its mode shifting toward zero while developing an extended tail toward higher \( \Delta^{\,}_i \) values. This behavior reflects the well-known phenomenon of the formation of localized superconducting 'islands' embedded in an otherwise suppressed superconducting background~\cite{Sadovskii1997, Ghosal2001, Gant2010, neverov2022}.

Figure 4

Pairing potential predictions by NN for various disorder strengths at every site of the system (left) and the corresponding pairing potential distributions (right). Red dots show values of numerical self-consistent solution and blue lines correspond to the NN predictions.

Notably, while the spatial average \( \bar{\Delta} \) decreases monotonically with disorder (as shown in Figure 3), the NN successfully identifies local regions where the pairing potential exceeds its clean-system value. These areas of enhanced superconductivity correspond to rare but physically significant configurations where disorder potential fluctuations create favorable conditions for pair formation. Such localized superconducting regions are particularly interesting as they may host elevated local critical temperatures while remaining unable to sustain global phase coherence or supercurrent.

The NN's predictions show excellent agreement with BdG calculations for moderate disorder strengths \( V^{\,}_0< 0.5 \), accurately reproducing the main lognormal distribution. For stronger disorder, minor discrepancies are observed near zero value. These differences likely stem from insufficient sampling of rare disorder configurations in the training set. We anticipate that expanding the training dataset to include more realizations of strongly disordered systems would improve the network's performance in this regime by better capturing the statistics of pairing potential fluctuations.

Our investigation further examines the system's capability to predict pair potential correlations in the presence of disorder (Figure 5). At weak disorder strengths, the \( \Delta^{\,}_i\Delta^{\,}_j \) correlation function remains nearly constant, demonstrating minimal spatial fluctuations, a characteristic signature of disorder having negligible impact on superconducting coherence. As disorder increases, the correlation function develops distinctive features: a pronounced dip emerges at intermediate distances, reflecting the formation of localized superconducting clusters embedded within the suppressed background.

Figure 5

The normalized pair potential correlations \( \Delta^{\,}_i \Delta^{\,}_j \) as a function of distance between sites \( i \) and \( j \) for different values of disorder strength \( V^{\,}_0 \). Dots represent self-consistent numerical solutions, while solid lines correspond to the NN predictions.

The depth of this correlation dip quantitatively measures the relative strength of these localized clusters compared to the average pairing amplitude. The model successfully captures the long-range behavior of correlations extending to distances up to \( 2R-1 \), beyond which the NN does not take correlations into account. Notably, the difference between the exact self-consistent solution and the predictions of NN at distances \( |i-j|\ge 2R \) is negligible.

The demonstrated accuracy in reproducing the complex spatial and statistical properties of pairing potential confirms that the NN has effectively learned the essential physics of disorder effects in superconductors. This includes not just the average suppression of superconductivity, but also the more subtle phenomenon of its local enhancement, a capability that makes this approach particularly valuable for studying inhomogeneous superconducting systems.

A key advantage of this approach lies in its scalability. The computational cost of a self-consistent BdG solution for a single disorder realization scales as \( \mathcal{O}(K L^2) \) for sparse methods (or \( \mathcal{O}(K L^3) \) if exact diagonalization is used), where \( K \sim 10^2-10^3 \) is the number of self-consistency iterations. In contrast, the neural network's prediction for the full spatial profile of the pairing potential scales only as \( \mathcal{O}(L) \), since the trained model is applied independently at each lattice site. The neural network architecture may be, therefore, effectively generalized to larger systems without requiring retraining. This feature is particularly valuable for studying extended disordered superconductors where traditional BdG calculations would be computationally prohibitive.

Our study demonstrates that a neural network can accurately predict the spatial distribution of the superconducting pairing potential in disordered quasi-one-dimen-sional systems, bypassing the need for iterative Bogoliubov-de Gennes calculations. The model captures both the average suppression of superconductivity with increasing disorder strength and the formation of localized regions with enhanced pairing potential, reproducing key features observed in direct numerical solutions.

While the neural network performs well for moderate disorder strengths, small discrepancies emerge in strongly disordered regimes, likely due to limited training data for rare configurations. The method's computational efficiency and scalability make it a practical tool for studying large disordered systems, though further improvements could be achieved by expanding the training dataset to better capture extreme disorder effects.

This work is supported by a grant from the Ministry of Science and Higher Education of the Russian Federation (project No. 075-15-2025-010).

Corresponding author: Vyacheslav D. Neverov,

orcid.org/ 0000-0002-8999-8297, e-mail slavesta10@gmail.com.

The authors declare no competing financial or non-financial interests.