Magnetic flux noise in superconducting qubits and the gap states continuum

In the present study we investigate the selected local aspects of the metal-induced gap states (MIGSs) at the disordered metal–insulator interface, that were previously proposed to produce magnetic moments responsible for the magnetic flux noise in some of the superconducting qubit modalities. Our analysis attempts to supplement the available studies and provide new theoretical contribution toward their validation. In particular, we explicitly discuss the behavior of the MIGSs in the momentum space as a function of the onsite energy deviation, that mimics random potential disorder at the interface in the local approximation. It is found, that when the difference between the characteristic electronic potentials in the insulator increases, the corresponding MIGSs become more localized. This effect is associated with the increasing degree of the potential disorder that was earlier observed to produce highly localized MIGSs in the superconducting qubits. At the same time, the presented findings show that the disorder-induced localization of the MIGSs can be related directly to the decay characteristics of these states as well as to the bulk electronic properties of the insulator. As a result, our study reinforces plausibility of the previous corresponding investigations on the origin of the flux noise, but also allows to draw future directions toward their better verification.

by the Ruderman-Kittel-Kasuya-Yosida (RKKY) interactions at the superconductor-insulator interface 21,22 . In what follows, the RKKY interactions were also adopted in the spin-cluster model by De 23 . Yet another proposal by Wu and Yu suggested that the noise emerges from the hyperfine interactions of the relaxing surface spins 24 . To this end, Choi et al. claimed increased role of the metal-induced gap states that become localized at the superconductor-insulator interface by the potential disorder 25 .
In general, the described above theories follow the most important experimental findings on the magnetic flux noise. Specifically, they render the universal and weak scalability of the noise with the overall size of the SQUID 8,25 . However, they also independently tackle other aspects of the discussed effect. In details, the model of Koch et al. 20 is supported by the recent experimental studies that consider qubit coherence optimization based on SQUID geometry 26 . Some other models emphasize the interactions between surface spins [21][22][23][24] , in agreement with the experimental observations provided in 27 . The model of Choi et al. 25 , concentrates on describing correctly the areal spin density from the susceptibility measurements 28 . Finally, the latest theoretical approaches attempt to include the role of an extrinsic effects coming from the surface adsorbents 29,30 . As a result, they altogether suggest complex nature of the noise that may arise from the interplay of a various intrinsic and extrinsic effects.
In this context, there is a strong motivation to verify each of the theoretically postulated contributions to the origin of the flux noise. This will allow to account for the most important physical parameters that govern the noise, toward its further reduction. In the present study, of special interest is the role of the metal-induced gap states (MIGSs), as proposed by Choi et al. in 25 . The importance of the MIGSs steams from the fact that they are considered to inevitably appear whenever metal-insulator junction (MIJ) is created 31,32 . Thus, they should be also present at the Josephson junctions that build SQUIDs. Moreover, the MIGSs are likely to localize at the discussed interface, as they constitute direct continuum of the metal states that decay into the first few layers of the insulator 33 . According to Choi et al. 25 , the mentioned localization of the MIGSs notably contributes to the flux noise (giving rise to the paramagnetic local moments with the observed areal density 28 ) and should occur due to the presence of the potential disorder at the metal-insulator interface. Such interplay between the interfacial disorder and the MIGSs localization is the central point of the approach presented by Choi et al. 25 . However, we argue here that the universal and inherent nature of this relation has not been yet explicitly recognized, in order to further reinforce importance of the MIGSs concept in explaining the magnetic flux noise origin.
With respect to the above, in this study we attempt to present new insight into the localization of the MIGSs, as compared to the study already presented in Choi et al. 25 . In particular, contrary to the supercell approach by Choi et al. 25 , our theoretical analysis is intended to directly examine relation between the potential disorder at the unit cell level and the MIGSs localization, near the metal-insulator interface. To do so, we use the so-called complex band structure method, instead of following the isotropically-averaged (in the momentum space) density of states scenario similar to the Cardona-Christiansen approximation 34 and adopted by Choi et al. 25 . In this manner, we are able to directly investigate the behavior of the MIGSs in the momentum space, when the onsite potentials varies. Hence, the behavior of interest is determined in terms of the insulator band structure (not the entire interface like in case of Choi et al. 25 ), so that the localization effect is explicitly related to the intrinsic properties of the interface. Note that this approach follows in spirit considerations presented by Tersoff for the Fermi level pinning at the metal-semiconductor junctions 33 . Moreover, since our calculations are employed at the local level they demonstrate importance of the relative energy difference and its fluctuations within unit cells, in correspondence to the fundamental features of the insulating layer. Therefore, the presented here study can be viewed not only as a supplementary analysis to the investigations conducted by Choi et al. 25 , but also first step in developing fully-anisotropic theoretical approach that explains disorder-induced MIGSs localization at the MIJs, in direct relation to the main characteristics of these states (e.g. their decay characteristics not considered in the study of Choi et al. 25 ).

Theoretical model
To investigate the localization effect of the MIGSs near the MIJ, we concentrate our attention on the insulator region in a benchmark CsCl structure. The electronic properties of such structure are described in the framework of the mean-field Hubbard Hamiltonian given as: where ε i is the electronic potential of the s-type orbital at the ith site, equal to -4 and 2 eV for the undisturbed Cs and Cl atoms, respectively. In what follows, the c i,σ ( c † i,σ ) operator creates (anihilates) electron at the ith site with the spin σ =↑, ↓ , whereas j denotes nearest-neighbors (NN) or next-nearest-neighbors (NNN) of i. Hence, the t i,j is the hopping energy set to − 0.5 eV for both the NN and NNN cases. Note that the tight-binding parameters are adopted here directly from the study of Choi et al. 25 , where above values were chosen in order to set the van Hove singularity away from the conduction band edge at the MIJ, therefore keeping the density of states relatively realistic in the vicinity of the insulators band gap. To this end, the U parameter describes the on-site Coulomb repulsion, with n i,σ = c † i,σ c i,σ being the number operator. We take U = 3.25 eV, to account for the well localized electrons, as suggested in 25 . Note, that above assumptions allow us to conduct our analysis on the same footing with the discussion provided by Choi et al. (see 25 for more details). Moreover, the assumed U value is much greater than the superconducting pairing gap, hence our approach is expect to be valid also in the case when the attached metal is superconducting. Nonetheless, it is instructive to note here that we do not explicitly consider attached metal, by arguing the fact that MIGS should constitute intrinsic property of the insulator, according to  32,33,[35][36][37] . This is to say, the interfacial behavior of MIGS is described by the insulator complex band structure (CBS), that can be generalized to the interface 32 . Such approach is possible, since MIGS constitute direct analytical continuation of the propagating states in a metal 37,38 . It also allows to trace, hitherto not considered, canonical aspects of the MIGS localization at the MIJ. In particular, the MIGSs are calculated here by adopting the CBS method, that requires us to solve the following generalized eigenvalue problem 32,39 : In Eq. (2), the H n and H n+1 matrices are the component Hamiltonian terms for the reference unit cell and its coupling to the neighboring cells, respectively. In this manner, the ψ n and ψ n+1 describes the wave function coefficients accordingly associated with the n and n + 1 cell. These coefficients are chosen so that they satisfy the following phase relation ψ n+1 = ϑψ n , where ϑ denotes the generalized Bloch phase factor. The problem of Eq. (2) is solved in the self-consistent iterative manner, to match the paramagnetic solutions of Eq. (1). Note, that the paramagnetic behavior of the insulator is dictated by the experimental results of Sendelbach et al. 28 . For more technical details on the CBS method we refer to the studies already mentioned above i.e. 32,39 .
In general, the Eq. (2) produces the pairs of ϑ and 1/ϑ eigenvalues that are related to each other by the timereversal symmetry. Herein, we consider the behavior of such solutions in the momentum space ( k ). When |ϑ| = 1 the solutions correspond to the typical propagating states. On the other hand, when |ϑ| < 1 some of these eigenvalues can be interpreted as a MIGSs. Specifically, the MIGSs are the states with |ϑ| < 1 that appear in the energy gap and create characteristic complex band loops, joining the maximum of the valence band with the conduction band minimum. Of particular importance to the present discussion, is the point in energy space where the MIGSs exhibit the largest density. According to Choi et al., these are the states most susceptible to the localization that should provide biggest contribution to the flux noise. By recalling the sum rule on the density of states 40,41 , the density of MIGSs must be derived from the contributions of the conduction and valence bands of the insulator. In this context, the highest density of MIGSs is expected at the point where the conduction and valence bands contributes equally to the MIGSs (the charge neutrality point, also called the branch point). By following Allen 42,43 and Tersoff 33 , this point can be located within the energy gap by using the following cellaveraged Green's function: where the parameter N gives the number of the unit cells in the considered system, n counts all the accessible eigenvalues in the momentum space ( E n,k ), and η takes infinitesimal positive value to differentiate between the advanced and retarded Green's function. Moreover, the R is the lattice vector of the CsCl crystal structure. Since the function of Eq. (3) varies its sign along with the bonding character of the considered eigenvalues, it yields the location of the charge neutrality point at the value of 0. The Im[k] value at this energy level is then interpreted as a characteristic decay rate ( κ ) of the MIGSs within the insulator energy gap; according to the wave function decay defined as e −κa , with a being the lattice constant. Hence, it provides the effective measure of how well given decaying states are localized near the MIJ.
In the context of the above, the present study attempts to analyze how the characteristic κ value changes as a function of the potential disorder at the interface. However, to provide the most in-depth insight into the analyzed processes, their behavior is related here to the fundamental aspects of the potential disorder at the local level. In particular, when the random potential disorder is induced it varies the onsite energies ( ε i ) at the two inequivalent sites (the Cs and Cl atoms) within each of the unit cells that build the insulator. To characterize such potential deviations at the local level, we introduce energy deviation parameter: is the unperturbed onsite energy at the Cs (Cl) atom, as given above, and ε ′′ Cs ( ε ′′ Cl ) denotes its perturbed counterpart. In this manner, the above parameter describes respective differences between perturbed and unperturbed potentials at the Cs and Cl atoms within the unit cell. In particular, when = 0 eV, it corresponds to the unperturbed CsCl structure. On the other hand, when is positive (negative) and takes on increasing values, the potential distribution is perturbed and the asymmetry between Cs and Cl sites decreases (increases).
We note that parameter, in our opinion, is of particular importance to the behavior of the MIGSs, since mentioned energy difference models region in the vicinity of the energy gap, as well as the gap itself. Therefore, it can be expected that this difference has particular influence on the behavior of the MIGSs. This fact is additionally reinforced by already mentioned observations made by Choi et al., that suggest states located within the band gap or at its edges to be most susceptible to the localization. Note, however, that the presented approach is not limited to the local picture and can be extended further to analyze processes of interest in the framework of a large scale calculations e.g. within the approach derived from the Anderson localization model, as presented in 25 . (2)

Numerical results
In Fig. 1A we depict the CBS solutions in the complex plane of ϑ , that has been calculated for the unperturbed CsCl structure by using Eq. (2). For convenience, the projection of these results on the Re[k] (left panel) and Im[k] axes is presented in the right and left panel of Fig. 1B, respectively. Altogether these are the reference results where deviations of the onsite energies within the insulator unit cell ( ) are set to zero (the CsCl structure is assumed to be ideal). This is to say, ε ′ Cs = ε ′′ Cs and ε ′ Cl = ε ′′ Cl in case of the results depicted in Fig. 1A,B. The purpose of Fig. 1A,B is to conveniently introduce reader to the employed CBS representation and briefly discuss its characteristics, prior to the main discussion. In general, the CBS, like the one in Fig. 1A,B, presents set of three different types of states, that are characterized by the complicated functional behavior. The propagating states correspond to the solutions for |ϑ| = 1 and real k, composing characteristic unit circles in the complex plane of ϑ (see Fig. 1A). For clarity, in Fig. 1A we depict only states for Re[ϑ] ∈ �−1, 1� ∩ Im[ϑ] ∈ �0, 1� , due to the symmetry of the CBS solutions about the Im[ϑ] axis. Note that these states are valid only deep into the bulk of the ideal periodic system and correspond to the typical electronic band structure with its characteristic features, e.g. the indirect band gap of the CsCl structure as visible in Fig. 1B. For the purpose of this study, however, we are only interested in the evanescent states that appear when the perfect periodicity of the system is broken e.g. at the surfaces, scattering incidents or the discussed here interfaces 44 . In what follows, we distinguish two major types of the evanescent states that can be observed within our theoretical model. First, we note the evanescent states that appear at the boundaries of the adopted Brillouin zone and are depicted by the green color ( |ϑ| < 1 ∩ Im[ϑ] = 0 and k is real). These states exhibit exponentially decaying character inside the unit circles defined by the propagating states. Next, we observe the gap states of interest, that are given in the red color ( |ϑ| < 1 ∩ Im[ϑ] � = 0 and k is complex). Contrary to the boundary states they decay in the oscillatory manner and compose already mentioned complex band loops, that join maximum of the valence band with the minimum of the conduction band and can be interpreted as the MIGSs. To elucidate the presented results even further we remind that each of the above solutions appears in pairs i.e. the pairs of ϑ and 1/ϑ eigenvalues. In details, the gap states are symmetric about the Im[ϑ] axis (like the propagating states), whereas the edge states exhibit symmetry about the Re[ϑ] axis (see Fig. 1A). The component solutions of a given pair are interpreted as the right and left propagating As mentioned earlier, we argue here that the deviation is the local manifestation of the potential disorder at the interface.
By inspecting the results given in Figs. 2 and 3, we can observe that the functional behavior of all states in the energy-, ϑ -and k-space changes notably as the parameter varies. Of particular importance to our discussion is the behavior of the gap states. One can notice that, as the parameter decreases toward lower negative values, the total magnitude of the complex band loops increases. In details, the size of the complex loops increases in the energy-, ϑ -and k-space. This behavior occurs due to the fact that the onsite energies at the Cs and Cl atoms determine the upper and lower energy limits of the complex band loops. Therefore, when their relative energy difference changes the size of the complex band loops changes accordingly. It can be then qualitatively argued that strongly influences the MIGSs at the MIJ, including the pivotal position of the branch points associated with the complex band loops (see Eq. 3). However, to investigate how the parameter is related to the localization strength of the MIGSs, their behavior must be analyzed in the momentum space.
In the context of the above, we adopt the cell-averaged Green's function method of Eq. (3) to elucidate the relation between the MIGSs localization strength and the parameter. Herein, we use the numerical techniques previously adopted in 37 . Specifically, by solving the Eq. (3) for each of the presented cases, we search for the energy level ( E B ) that corresponds to the location of the mentioned branch point at the given complex band loop. Next, we extract the corresponding value of the characteristic decay rate parameter, in each of the considered cases, by using the following relation κ = − Re[log(ϑ(E B ))/a] . The values of κ , determined in the described way, are presented in Fig. 4A as a function of ∈ [− 5 eV, 5 eV] . Note that the decay rate parameter is expressed there in the units of the inverse unit cell size of the CsCl structure (1/d). In correspondence to the behavior of the complex band loops, the decay rate κ decreases, almost linearly, along with the increase of the parameter. In this manner, the highest values of κ are obtained for the gap states with the biggest magnitude in the ϑ -and k-space (see Figs. 2 and 3, respectively). However, to derive the localization characteristics of the MIGSs from the results presented in Fig. 4A, it is instructive to recall the fact that the wave function decay is given by e −κa . With respect to this relation, the most localized states are the ones described by the the highest values of κ . This statement is additionally reinforced by the fact that decay rate is simply the inverse of the decay length ( ), namely = 1/κ as depicted in Fig. 4B. In other words, the MIGSs characterized by the largest κ values, exhibit the shortest decay lengths (see Fig. 4).
The results presented in Fig. 4A,B allow us to draw several important observations regarding the localization of the MIGSs at the MIJ. Of particular importance is the fact that the provided results explicitly show direct relation www.nature.com/scientificreports/ between the MIGSs localization strength and the onsite energy deviation parameter ( ). By arguing the fact that parameter is the local manifestation of the potential disorder at the MIJ, one can observe that certain interfacial variations of the potential may cause MIGSs to become more localized. In reference to the study of Choi et al. 25 , the energy distribution at the perturbed MIJ can be given as P( where E 0 denotes the unperturbed onsite energy and δ is the standard deviation. In his work Choi et al. defines additionally the disorder degree as R = 2δ/W , with W being the metal bandwidth. In the context of these relations, Choi et al. show that the localization of the MIGSs increases together with the increase of the R parameter, that is proportional to the standard deviation. This observations can be related to our study, by saying that represents selected changes of the potential distribution within given unit cell, when the δ parameter increases.

Summary and conclusions
In summary, by using the complex band structure method we have determined the localization behavior of the MIGSs in the momentum space, with respect to the local onsite energy deviation ( ) at the MIJs and bulk properties of the benchmark CsCl insulator. The analysis was performed to supplement findings by Choi et al. 25 on the origin of flux noise in superconducting qubits. Specifically, we have found that the MIGSs localization increases along with the increase of the energy deviation. This behavior is related in our analysis to the fact that difference between the Cs and Cl onsite energies increases with respect to the unperturbed case. In terms of the MIJs, the described localization increase corresponds directly to the potential barrier growth at the interface. This is to say, higher difference of the onsite energies causes bigger energy separation between the band edges of the propagating states that give rise to a given MIGS. As a result, the discussed parameter has been related here to the standard deviation within the disorder degree defined by Choi et al. 25 . This allows to suggest that the calculated MIGSs correspond to the highly localized states observed by Choi et al. 25 under substantial degree of the interfacial disorder.
In what follows, in our analysis we confirm that the variations of the potential at the interface can cause stronger localization of the MIGSs, that appears as an effect directly related to the inherent electronic properties of the insulator material. At the same time, the obtained results validate, to some extent, the proposal of Choi et al., saying that the MIGSs can play an important role in describing the magnetic flux noise in superconducting qubits. However, it should be noticed that further investigations are desirable, in order to combine our approach with the large scale modeling techniques, toward more comprehensive analysis of the discussed processes. In particular, it should be possible to adapt our theoretical techniques within the calculations that tackle realistic metal-insulator junctions with the random potential disorder. As least two directions in this respect can be listed: (1) the Anderson-derived scenario as considered by the Choi et al. 25 or (2) the technique based on the renormalization group approach 45 . Such investigations are expected to allow to relate the decay characteristics of the MIGSs, not only to the potential fluctuations, but also the experimentally observed areal density of the paramagnetic spins in SQUIDs. As a result, further verification of the flux noise origin in terms of MIGSs should be feasible. Simultaneously, we note that the described improvement of the theoretical techniques may allow to consider still open problems such as the role of interactions between the paramagnetic spins in producing the flux noise. To this end, the presented results indicate potential experimental routes toward suppression of the noise. In particular, they suggest avenues aimed at the reduction of the potential disorder at the MIGSs (e.g. by producing interfaces via epitaxial growth) or engineering Fermi level position within insulator toward regions less sensitive against disorder (the band gap edges).