Anomalous isotope effect in BCS superconductors with two boson modes

The isotope effect in the superconducting transition temperature is anomalous if the isotope coefficient α < 0 or α > 1/2. In this work, we show that such anomalous behaviors can naturally arise within the Bardeen–Cooper–Schrieffer framework if both phonon and non-phonon modes coexist. Different from the case of the standard Eliashberg theory (with only phonon) in which α ⩽ 1/2, the isotope coefficient can now take arbitrary values in the simultaneous presence of phonon and the other non-phonon mode. In particular, most strikingly, a pair-breaking phonon can give rise to large isotope coefficient α > 1/2 if the unconventional superconductivity is mediated by the lower frequency non-phonon boson mode. Based on our studies, implications on several families of superconductors are discussed.


Introduction
Isotope effect [1,2] is a cornerstone of the Bardeen-Cooper-Schrieffer (BCS) theory [3] for phonon mediated superconductors. It describes the change of transition temperature T c caused by isotope substitution. The standard BCS theory gives T c = 1.13Ωe −1/λ and thus predicts the isotope coefficient where M is the ion mass, Ω is the Debye frequency and λ is the dimensionless isotope independent electron-phonon coupling constant. Including the Coulomb pseudopotential [4] μ * = μ/[1 + μ ln(E F /Ω D )] (μ is defined at the Fermi energy E F ) will reduce the isotope coefficient even to negative values, as clarified in the more elaborated Eliashberg theory [5][6][7][8].
In real materials, the anomalous isotope effect (defined as α < 0 or α > 1/2) has been observed in many experiments. (Although the negative isotope coefficient α < 0 is a standard behavior predicted by the Eliashberg theory, it can only occur at very low T c and thus is inconsistent with the experiments. Therefore, people still prefer to call α < 0 anomalous.) Among all these materials, cuprates may be the most systematically studied in the past thirty years [9]. In cuprates, the isotope coefficient of the O-atoms in the CuO 2 -plane [10] is found to nearly vanish at optimal doping [11] and to increase with decreasing T c either upon under-or over-doping, to values even larger than 1/2 [12,13]. Such an interesting observation has stimulated many theoretical works on the role of phonons in the superconductivity mechanism of cuprates [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27], although spin and other types of fluctuations are also widely observed to be closely related to superconductivity [28][29][30][31]. In fact, cuprates are not the only material having α > 1/2. Such an anomalous behavior has also been observed in some iron- [32] and C 60 -based superconductors [33][34][35]. In Sr 2 RuO 4 , even a similar scaling behavior of α versus T c was reported and drops to negative values near the maximum T c . [36] Such a negative isotope coefficient cannot be explained by the standard Eliashberg theory where α can drop to negative values but only as T c diminishes [8]. Similarly, the inverse isotope effect has also been observed in iron-based superconductors [37] and PdH [38]. In the latter case, the anomalous isotope effect has been attributed to anharmonic phonon effect [39][40][41].
How to understand these anomalous isotope effect in a unified way is a longstanding problem. This issue has been explored rather extensively in the literature. Among these studies, for cuprates in particular, some material dependent properties such as Van Hove singularity [14,18], pseudogap [19], anharmonic phonon effect [16,22], or bipolaron [13,26,27] are proposed to be responsible for the anomalous isotope effect. A somewhat more 'universal' approach was to consider a pair-breaking non-phonon mode which is found to cause large α > 1/2 [15,21,24]. However, most of these studies assumed phonon mediated superconductivity with λ ph > 0 (henceforth the subscript ph/nph stands for phonon/non-phonon in this work). But for unconventional superconductivity (other than uniform s-wave pairing), λ ph can be negative for a given phonon mode [25,42]. See the appendix for several typical phonon modes in cuprates for example. Another aspect for some of these works may be the heavy dependence on the material details, e.g. complex electron band structures or electron-boson coupling functions α 2 F, to solve the Eliashberg equations. But a universal understanding of the anomalous isotope effect is still lacking.
Motivated by these experimental and theoretical progresses, we ask a somewhat very simple question: what happens in the BCS framework (not only with s-wave pairing) including two kinds of single-frequency bosons Ω 1,2 with the electron-boson couplings λ 1,2 ? Thanks to the small number of model parameters, a thorough study becomes possible. Interestingly, all types of values of α can appear within this simple approach. In particular, we find a pair-breaking phonon can give rise to α > 1/2 as long as the other non-phonon mode has a lower frequency. Implications on different superconducting materials will be discussed.

BCS theory with two boson modes
Boson mediated interactions are retarded. One key point of the standard BCS theory [3] is to simplify the frequency dependency of the boson mediated pairing interactions into the momentum space within an energy shell: with Θ the step function and ε k the normal state band dispersion. Following this idea directly, including two boson modes, we have the interaction Hamiltonian H I = − kk V kk Δ † k Δ k with Δ k the pairing operator and V kk the pairing interaction given by where the pairing interactions have been decoupled into different symmetry channels (e.g. s-, p-or d-wave etc.) labeled by the subscript m with f m (k) defined as the form factor (the mth eigenvector of V kk ). V m1 and V m2 are effective interactions (positive for attractive ones) mediated by the Ω 1 and Ω 2 modes, respectively. U m is the instantaneous interaction. Without lost of generality, Ω 1 Ω 2 is always assumed in the following. At T = T c , we have the linearized gap equation where N k is the number of momentum k. As a result of the three-piecewise behavior of V kk in equation (2), Δ(k) can be approximated as a three-piecewise function. Then, the momentum summation (replaced by energy integral with constant density of states) can be performed in three regimes: (0, Ω 1 ), (Ω 1 , Ω 2 ) and (Ω 2 , E F ), respectively. For the m-wave pairing (it means the pairing function is f m ), the kernel becomes (assuming T c Ω 1 ) where κ = 2e γ /π ≈ 1.13, λ 1,2 = N V m1,2 f 2 m , and μ = N U m f 2 m with N the electron density of states and the average is performed on the Fermi surface. T c is then determined by letting the largest eigenvalue of K Figure 1. T c and α 1,2 as functions of λ 1 and λ 2 by setting Ω 2 = 2Ω 1 and μ * = 0, 0.15, −0.15, respectively. For clarity, we divide the phase diagrams into different regimes: α < 0, 0 < α < 1/2 and α > 1/2 which are denoted as A, B and C, respectively. Since α 1,2 ranges from −∞ to ∞, we plot tanh(α 1,2 ) instead.
to be 1, i.e. det(K − I) = 0, giving rise to the T c -formula: where μ * is the Coulomb pseudopotential defined at Ω 2 as usual, i.e. [4]. It can be easily checked by setting λ 1 = 0 or λ 2 = 0, equation (5) does reduce to the standard BCS result. In fact, equation (5) can be rewritten in a more familiar way: is the 'pseudopotential' defined at Ω 1 and contributed by (μ * − λ 2 ) defined at Ω 2 . (Here, we slightly generalize the concept of the pseudopotential to include both instantaneous and retarded interactions above a given frequency, which can also be understood as the ladder approximation in the pairing channel.) Next, we turn to the isotope effect. From equation (5), the isotope coefficient α 1 or α 2 can be obtained exactly, corresponding to Ω 1 or Ω 2 as the phonon mode, respectively: and Clearly, if both Ω 1 and Ω 2 are phonons, the total isotope coefficient α 1 + α 2 is always less than 1/2 in agreement with the Eliashberg theory [8]. Now, let us suppose only one of them is phonon. In figure 1, we plot T c and α 1,2 as functions of λ 1 and λ 2 by fixing Ω 2 = 2Ω 1 and μ * = 0, ±0.15, respectively. Here, μ * < 0 is also considered in order to describe instantaneous interaction induced superconductivity (e.g. negative-U Hubbard model). For clarity, we divide the phase diagrams into several regimes: Let us firstly examine α 1 if Ω 1 is the phonon, as shown in figures 1(b), (e) and (h). Only regimes A and B appear, corresponding to α 1 1/2, which can be seen directly in its expression equation (7). An interesting feature is that as T c decreases α 1 drops to negative values and finally diverges logarithmically: α 1 ∼ − ln 2 (T c ), unless in the critical point (λ 1 = 0 and λ 2 = μ * ). This feature is already captured by the standard Eliashberg theory [8] when λ 1 > 0, as a result of the poorer screening of μ * 1 caused by increasing Ω 1 . On the other hand, if λ 1 < 0, the negative α 1 is expected as a result of its pair-breaking effect directly.
Our new result is in α 2 if the higher frequency boson mode Ω 2 is the phonon, as shown in figures 1(c), (f) and (i). All three regimes can be found for nonzero μ * . Large α 2 > 1/2 (regime-C) is a ubiquitous feature for λ 1 λ 2 < 0. Sandwiched between them is the regime of normal isotope coefficient 0 < α 2 < 1/2 (regime-B). For μ * = 0, there is also a regime of α 2 < 0 (regime-A). In fact, as seen from equation (8), α 2 ∼ λ 2 (λ 2 − 2 μ * )ln 2 (T c ) as T c → 0. When λ 2 > 0, the large α 2 originates from the pair-breaking effect of the Ω 1 mode, as being discussed in references [15,21,24]. Astonishingly, we have found another regime having large α 2 for λ 2 < 0. At first glance, this seems to be impossible since in this case the phonon is harmful to superconductivity. How can a 'repulsive' or pair-breaking phonon cause a large positive isotope effect? The answer is: the 'pseudopotential' μ * 1 (equation (6)) contributed by (μ * − λ 2 ) can be reduced by increasing Ω 2 , hence, leading to higher T c . In figure 2, we plot several typical scalings of α 2 versus T c along different paths as shown in the inset. Line cut L1 is described by the standard Eliashberg theory by setting λ 1 = 0. Both L2 and L3 give large α 2 as T c decreases, while L2 also have negative α 2 < 0 in the high T c regime. Along L4, α 2 is always negative and greatly enlarges the negative isotope coefficient regime of the standard Eliashberg theory (L1). Discussions of these theoretical results in light of real materials are left to the next section.
In the above, we have found a universal scaling as T c → 0 by tuning λ 1,2 to suppress the superconductivity. In practice, there is another theoretical possibility to get T c → 0 by tuning Ω 1 → 0 when λ 1 > 0 and λ 2 − μ * < 0. In this case, we get the standard BCS result α 1 = 1/2 and α 2 ∼ ln −2 (Ω 1 ) → 0. From the above discussions, we have seen that equations (7) and (8) give a full description of all possible values of the isotope coefficient in a BCS superconductor with two boson modes. There are mainly three approximations in the above theory: (1) We treat the boson mediated retarded interactions in momentum space directly. The frequency dependence of the gap functions is changed into momentum dependence effectively. This is just the standard BCS approximation [3]. (2) We have ignored band renormalization effect caused by the boson modes. (3) We have assumed the pairing interactions and gap functions to be approximated by the three-piecewise functions. In order to justify our approximations, we have performed numerical studies of the Eliashberg theory and obtained similar results, indicating the above BCS picture indeed works and captures the main physics qualitatively. See the appendix for more details. In particular, the band renormalization effect can be approximately included in our BCS treatment. The only difference is the kernel now becomes where Z 1 ≈ 1 + λ 1z + λ 2z and Z 2 ≈ 1 + λ 2z . It should be emphasized that for phonons λ z = λ s and for magnetic modes λ z = −λ s where λ s denotes the electron-boson coupling constant in the uniform s-wave paring channel. (See the appendix for more details.) We have checked that including the band renormalization effect does not change the above BCS results qualitatively. The derivations of T c and α 1,2 are straightforward and thus left out. Their analytical expressions can be found in the appendix. Another approximation in our study is that both boson modes are taken as Einstein modes, i.e. single-energy modes. For continuous boson modes, Ω 1 and Ω 2 can be understood as their representative energies. In most materials, the phonon and non-phonon modes range in different energy regions, and hence, we expect the above Einstein mode approximations are qualitatively correct.

Summary and discussions
In summary, we have found the anomalous isotope effect (α < 0 or α > 1/2) can be explained in the BCS theory when both phonon and non-phonon modes coexist. If the phonon frequency is lower, α 1/2. But if the phonon has a higher frequency, any values of α can be obtained. Interestingly, we have obtained a scaling behavior α → ± ln 2 (T c ) as T c → 0 by tuning λ 1,2 to suppress the superconductivity. Most strikingly, when the phonon mode has a higher frequency, α can be larger than 1/2 even if it is pair-breaking. Finally, as possible applications, we provide some remarks about several superconducting materials from the viewpoint of the BCS framework. (1) Cuprates. Low energy boson modes have been widely observed in many different experiments in cuprates, including mainly two candidates: phonon and magnetic modes [29,30,43]. We have listed several candidate boson modes in the appendix. Taking different experiments together, roughly speaking, the 70 meV kink observed in angle-resolved photoemission spectral (ARPES) can be assigned to the breathing phonon [44][45][46] and the low energy (10-60 meV including the famous 41 meV resonance [47]) [29][30][31] are dominated by antiferromagnetic (AF) excitations with the hour glass dispersion [48], although phonon may also have some contributions [20,[49][50][51][52][53]. Since the breathing phonon is against d-wave superconductivity [17] but the AF fluctuation can mediate d-wave superconductivity [42,[54][55][56], within our BCS picture, they correspond to 10 meV < Ω 1 < 60 meV, Ω 2 ∼ 70 meV, λ 1 > 0 and λ 2 < 0. For a rough but specific estimation, if we choose Ω 1 ∼ 35 meV, Ω 2 ∼ 70 meV, λ 2 ∼ −0.1, then the maximal T c ∼ 80 K at optimal doping corresponds to λ 1 ∼ 0.8. Upon doping away from the optimal doping, λ 1 drops to reduce T c and enhance α 2 and finally leads to the scaling behavior α 2 ∼ ln 2 (T c ) as T c → 0. This behavior is similar to the experimental observations [12,13]. In our theory, adding μ * does not change the qualitative behavior and thus is not in contradictory with the resonating valence bond (RVB) theory [57] which corresponds to an attractive pseudopotential μ * < 0. But based on our picture, the lower energy boson mode is necessary to obtain α 2 > 1/2. In particular, for La 2−x Sr x CuO 4 near 1/8 doping [58] compared with YBa 2 Cu 3 O y [59], stronger charge fluctuation results in weaker λ 1 and leads to smaller T c and larger α 2 , also in agreement with the experiment. (2) Sr 2 RuO 4 . Ferromagnetic (FM) fluctuations are widely believed to mediate the superconductivity in Sr 2 RuO 4 [60]. The magnetic mode energy is found to be less than 15 meV [61], much less than the O-phonon frequency around 50 meV [36]. Therefore, it is in the same parameter regime as cuprates. As a result, its α versus T c shows similar behavior as cuprates except α 2 < 0 for higher T c samples [36]. (3) Iron-based superconductors are found to be similar to cuprates in the sense that AF fluctuations are closely related to superconductivity (either s ± -or d-wave). [62,63] If we take the AF fluctuation as Ω 1 ∼ 15 meV [64,65] and phonon as Ω 2 ∼ 40 meV [66], we can obtain both α 2 < 0 [37] or α 2 > 1/2 [32]. However, further systematic experiments of the isotope effect upon doping are needed to pin down the role of phonons.
(4) C 60 -based superconductors. In fullerides superconductors A 3 C 60 (A stands for K, Rb, Cs), phonon mediated s-wave pairing is widely accepted [33]. α > 1/2 has been reported [34,35] and explained by the breakdown of Migdal theorem [67,68]. Nevertheless, our theory provides another possibility: existence of a lower frequency non-phonon mode can also lead to large α. Interestingly, in A15-Cs 3 C 60 superconductivity is found to be near the AF parent [69] such that the spin fluctuation may also play some role in it [70].
Although the above discussions are based on the BCS framework, we hope the qualitative conclusions should be the leading effect even beyond the weak-coupling limit. Therefore, we emphasize that the isotope effect can be strongly affected by the non-phonon boson mode.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.

Appendix A. Eliashberg theory
In this section, we make a benchmark for our BCS approach by numerically solving the Eliashberg equations for general pairing symmetries. At first, for simplicity, we consider only one phonon mode and give a self-contained derivation. The Eliashberg theory is based on the self energy in Nambu space [8,71,72], where G/D are the electron/phonon propagators, g is the electron-phonon vertex, and p(p ) stands for both momentum and frequency. σ 3 is the third Pauli matrix. By choosing the ansatz: and comparing two sides of equation (A.1), we get (particle-hole symmetry is assumed here and can be generalized straightforwardly) Then, singular mode decomposition is performed for |g| 2 D such that |g(p − p )| 2 = m g 2 m f m (p)f m (p ) where f m (p) are form factors in different symmetry channels. Next, we take two other ansatzs [72]:   where Ω is the solid angle (not confused with Debye frequency Ω D or Ω 1,2 ) and In the following, we neglect the tilde symbol inλ m for simplicity. Before going on, a short discussion on the Coulomb pseudopotential is given. Coulomb pseudopotential μ is not others but a boson mode with infinite frequency such that its λ m is frequency independent. Therefore, it has no contribution to Z but has to be included in the self-consistent equation of Δ m . In practice, the Coulomb pseudopotential may also be defined at a middle frequency ω c satisfying [4,8,73]. In practice, equation (A.10) is firstly solved to obtain Z(ω n ) numerically. Then, T c is obtained by finding the largest eigenvalue of the kernel to be 1 and Δ m is given by the eigenvector. After obtaining imaginary frequency data Z(iω n ) and Δ(iω n ), we perform the analytical continuation using the Padé approximation [74]. The results are shown in figure A1 by setting λ iz = |λ i | in (b) and λ iz = 0 in (c), respectively. Z(ω) and Δ(ω) are found to show drastic change near Ω 1 and Ω 2 , supporting our three-piecewise approximation in the main text. As a further benchmark, we also present the phase diagrams on the λ 1 − λ 2 plane in figure A2, which are in quite good agreement with the BCS theory.

Appendix B. Modified BCS theory
In this section, we show the results of the modified BCS theory by considering the band renormalization effect. In this case, the kernel is given by equation (10) in the main text. Following the method in the main text, T c and α 1,2 can be obtained as follows and and Notice that the scaling behavior of α 1,2 ∼ ± ln 2 (T c ) as T c → 0 keeps unchanged.
In this section, we list several typical phonon and magnetic modes in cuprates in table C1 together with three instantaneous interactions which can be taken as the pseudopotentials. All phonon modes have positive λ s and all magnetic modes have negative λ s . Therefore, for conventional uniform s-wave superconductors, phonon can mediate superconductivity but the magnetic modes only cause pair-breaking. Quite differently, for unconventional superconductors, both phonon and magnetic modes can be either positive or negative depending on different pairing symmetries. For the d x 2 −y 2 -wave pairing in cuprates, B 1g -buckling phonon mode has a positive λ d due to its form factor cos 2 (q x /2) + cos 2 (q y /2) and has been used as one candidate of the pairing mechanism [17,20,23,25]. However, the buckling mode requires the mirror symmetry breaking [17] and does not exist in single layer cuprates. Differently, the breathing phonon mode always exists and has been evidenced in ARPES experiments as the 70 meV kink [44,46], which is in fact against d-wave SC since its λ d < 0 as a result of its form factor sin 2 (q x /2) + sin 2 (q y /2) . [17,25] Besides, the Holstein phonon has no direct d-wave component and thus can be neglected in the leading order approximation without considering its coupling to other interaction channels.
On the other hand, the magnetic fluctuations are widely observed [29][30][31] and believed to be closely related to the d-wave superconductivity, including the spin fluctuation mechanism [42,[54][55][56] and the emergent effective SO(5) symmetry of the t-J model [75,76]. Combining most experiments, especially neutron and ARPES, it is reasonable to assume its energy ranging from 10 meV to 60 meV [29][30][31]43]. For simplicity, we have considered two extreme cases: AF and FM fluctuations with opposite λ d .
Finally, in order to include the instantaneous interactions, we also consider three interaction terms: Hubbard, Heisenberg exchange and Coulomb interaction between nearest neighboring sites. Although the Hubbard term has no d-wave pairing interaction directly, the Heisenberg term does have attractive component in the d-wave pairing channel, which in fact plays the essential role of pairing in the RVB theory [57,[77][78][79]. In addition, the nearest neighbor Coulomb interaction (which already exists in the t-J model) leads to a repulsive pairing interaction and thus should contribute to a positive Coulomb pseudopotential.