key: cord-0981150-jcegv8on authors: Wang, Da; Xu, Jun-Qi; Zhang, Hai-Jun; Wang, Qiang-Hua title: Anisotropic scattering caused by apical oxygen vacancies in thin films of overdoped high-temperature cuprate superconductors date: 2021-06-17 journal: Physical review letters DOI: 10.1103/physrevlett.128.137001 sha: 21ffd1b2919c9db6d692c2006f91e28c9b8833a3 doc_id: 981150 cord_uid: jcegv8on There is a hot debate on the anomalous behavior of superfluid density $rho_s$ in overdoped La$_{2-x}$Sr$_x$CuO$_4$ films in recent years. The linear drop of $rho_s$ at low temperatures implies the superconductors are clean, but the linear scaling between $rho_s$ (in the zero temperature limit) and the transition temperature $T_c$ is a hallmark of the dirty limit in the Bardeen-Cooper-Schrieffer (BCS) framework [I. Bozovic et al., Nature 536, 309 (2016)]. This dichotomy motivated exotic theories beyond the standard BCS theory. We show, however, that such a dichotomy can be reconciled naturally by the role of increasing anisotropic scattering caused by the apical oxygen vacancies. Furthermore, the anisotropic scattering also explains the"missing"Drude weight upon doping in the optical conductivity, as reported in the THz experiment [F. Mahmood et al., Phys. Rev. Lett. 122, 027003 (2019)]. Therefore, the overdoped cuprates can actually be described consistently by the $d$-wave BCS theory with the unique anisotropic scattering. Introduction.-For a long time, the overdoped cuprate superconductors were believed to be described quite well by the Bardeen-Cooper-Schrieffer (BCS) theory [1, 2] . However, in 2016, such a belief was challenged by the measurement of the superfluid density ρ s using mutual inductance technique on a large amount of high quality overdoped La 2−x Sr x CuO 4 (LSCO) films [3] . Two seemingly contradicting results were reported: ρ s (0)−ρ s (T ) ∝ T and ρ s (0) ∝ T c where ρ s (0) is the zero temperature value of ρ s (T ) and T c is the transition temperature. Within the conventional BCS theory, the former scaling law indicates the d-wave superconductors are very clean, but the latter is a result of dirty BCS superconductors [4] . This dichotomy regarding the dirtiness/cleanness of the overdoped cuprates motivated exotic theoretical investigations [5] [6] [7] [8] [9] [10] [11] [12] . The superfluid density has also been measured by THz optical conductivity experiment [13] , and is quantitatively consistent with the mutual inductance measurement [3] . In the meantime, the Drude fitting of the optical conductivity σ 1 (ν) indicates the dirty limit [14, 15] . So the same dichotomy also exists in optical measurements. Moreover, as yet another puzzle, σ 1 (ν → 0) should be identical to, but is fitted to be significantly smaller than the dc conductivity σ dc [13] . This superficial loss of Drude weight seems to increase with overdoping. In this Letter, we propose a scenario to resolve all of the above mysteries. We realize that in addition to the conventional isotropic scattering rate Γ s [14, 15] , the apical oxygen vacancies cause an anisotropic scattering rate Γ d cos 2 (2θ), with θ the azimuthal angle of the quasiparticle momentum relative to the antinodal direction. Since the low-energy nodal quasiparticles are barely affected by Γ d , they reduce the superfluid density linearly in temperature if in addition Γ s → 0. But the total scattering rate Γ θ = Γ s + Γ d cos 2 (2θ) determines the typical behavior ρ s (0) ∝ T c in the dirty limit. Meanwhile, the strong anisotropy in Γ θ causes a continuous distribution of Lorentzian components in σ(ν), so that σ 1 (ν → 0) would be underestimated by extrapolation from finite frequencies if a single isotropic scattering rate were assumed instead [13] . Oxygen vacancies and "cold spot" model.-From both early [16, 17] and recent [18] studies, apical oxygen vacancies were found to be an important factor to prevent obtaining (high-quality) overdoped LSCO samples. Therefore, the effect of apical oxygen vacancies deserves study in overdoped cuprates carefully. This is another motivation of the present work. In Fig. 1 , we plot a configuration of one apical oxygen atom above the CuO 2 -plane (only Cu-sites plotted). Notice that since the carriers have d x 2 −y 2 -like orbital content [19] , there is no coupling between the oxygen and the d x 2 −y 2 -orbital just below it. Therefore, the leading couplings are with the next nearest neighboring sites, giving the hopping integrals switching signs alternatively, as shown in Fig. 1 . First, if there is no oxygen vacancy, the second order process of these outof-plane hoppings gives an additional band dispersion −4κ(cos k x − cos k y ) 2 where κ = t 2 o /ε o with t o and ε o the hopping amplitude and energy level distance between the apical oxygen and d x 2 −y 2 orbitals [20] . Such a disper-sion has already enters the carrier band. Next, we consider one oxygen vacancy at the origin. Relative to the uniform background, the vacancy now leads to an additional term H i = δδ σ κf δ f δ c † δσ c δ σ where f δ = 1(−1) for δ = ±x(±ŷ) and σ denotes spin, which gives rise to the following scattering Hamiltonian where k and k are momentums, and N is the number of copper atoms. Since the oxygen vacancies are out-of-plane and only couple to next nearest neighboring d x 2 −y 2 -orbitals, the scattering strength is anticipated to be small. Under the Born approximation [4] , H i gives a scattering rate Γ(k x , k y ) ∝ Γ d (cos k x − cos k y ) 2 /4. For analytical convenience but without loss of qualitative physics, throughout this work, we further assume circular fermi surface and use wide band approximation, so that the scattering rate is only angle-dependent, i.e. where we have added the isotropic scattering rate Γ s arising from, e.g., in-plane impurities. In the following calculations, Γ s and Γ d are our model parameters. In fact, this kind of scattering rate has been proposed to explain the transport phenomena in underdoped cuprates, called "cold spot" model [21] , where the Γ d -term is attributed to the pair fluctuations or interaction effects, which hence is expected to decrease upon doing. But here, in our case, the Γ d -term is caused by oxygen vacancies and should increase upon doping, according to both early [16, 17] and recent [18] experiments. We now look at the effect of the anisotropic scattering Γ d on single particle excitations. By choosing the d-wave pairing function ∆(θ) = ∆ 0 cos(2θ) and fixing Γ s,d = 0.2∆ 0 , we calculate the density of states (DOS) ρ(ω) = dθ 2π ρ(ω, θ), where ρ(ω, θ) is the partial DOS (PDOS) given by where N is the normal state DOS and E θ = ε 2 + ∆(θ) 2 . ρ(ω) and ρ(ω, θ) are plotted in Fig. 2 (a)-(d). In the calculations, we perform the integrals over ε and θ using standard numerical integral technique. Fig. 2(a) is the typical behavior of the conventional dwave superconductor with isotropic scattering rate: zero energy cusp and coherence peak are both smeared out by isotropic scattering rate. Differently, the Γ d -scattering shown in Fig. 2 (b) only suppresses the coherence peak but does not change the zero energy cusp, as a result of the vanishing of scattering rate for nodal quasiparticles. Accordingly, we can divide the energy space into two regions: lower energy "clean" one and higher energy "dirty" one. The above picture is also seen clearly from the PDOS in Fig. 2(d) , where the spectral becomes smeared out for high energy quasiparticles (near the antinodes) but remains sharp for low energy ones (near the nodes). The dirty BCS theory for d-wave superconductors [22] is a direct generalization of the Abrikosov-Gor'kov (AG) theory [23] , since the potential scattering here causes pair-breaking which is similar to the magnetic impurities in s-wave superconductors. The gap ∆ 0 is determined by the self-consistent equation where λ is the BCS coupling constant, ω n = (2n − 1)πT the Matsubara frequency, and φ θ = cos 2θ the d-wave form factor. [24] In numerical calculations, we use λ = 0.3 and Ω = 800, which gives T c0 = 1.134Ωe −2/λ = 1.15 for clean superconductors (with this setup, 0.87T c0 is taken as the energy unit). By letting ∆ 0 → 0, we obtain the generalized T c -formula where ψ(z) is the digamma function. The results of T c vs Γ s and Γ d are shown in Fig. 2 (e). A stronger Γ d is needed to kill superconductivity as the low energy excitations are less affected. In Fig. 2 (f), we show the results of 2∆ 0 /T c . It is interesting to find that the anisotropic scattering drives the ratio farther away from the BCS prediction 4.28 in the clean limit to values as large as ∼ 10. In the literature, a large gap-T c ratio is often taken as an indication of strong coupling superconductivity. Here, the anisotropic scattering gives another interpretation. Superfluid density.-After recognizing the clean-dirty dichotomy caused by oxygen vacancies at the single particle level, we turn to discuss the superfluid density, which can be obtained as [25] where v F is the Fermi velocity. Here, the current vertex correction can be proven to vanish in the non-crossing approximation as a result of the factorized scattering potential (Eq. 1), as shown in the Supplementary Material [26] . In Fig. 6 (a) and (b), ρ s vs T is plotted for pure Γ s and Γ d scatterings, respectively. The most obvious difference is that any nonzero Γ s causes power law temperature dependence of ρ s (T ), but Γ d preserves the linear dependence as in the clean limit. This is already anticipated from the DOS feature in Fig. 2 (b) because the normal fluid density ρ n is roughly determined by quasiparticles and thus satisfies ρ n (T ) − ρ n (0) ∝ T , leaving the superfluid density satisfying ρ In real materials, both Γ s and Γ d are expected to coexist. In order to make quantitative comparisons with the experiment, we renormalize T and ρ s by T c and ρ s (0), respectively, as shown in Fig. 6 (c). We have superposed the experimental data of all 100 samples with T c ranging from 41.6 to 5.1K [3] , shown as the red dots. Four typical theoretical curves are plotted: (Γ s , Γ d ) = (0, 0), (0.1, 0), (1, 0), and (0, 1.2). As can be seen, almost all the experimental data reside within the region enclosed by the curves of (0.1, 0) (green line) and (0, 1.2) (blue line), while if Γ s is large and dominant, as the (1, 0)line shows, the renormalized ρ s -T curve bends more significantly, which is at odds with the experimental data. These observations indicate the more important role of Γ d . As an example, the experimental data (open circles) with T c = 5.1K are shown in the inset, which is clearly very linear. We can fit the data with pure Γ d = 0.75 (line) quite well. Another key observation of the experiment [3] is the linear relationship between ρ s (0) vs T c . To see that, we plot our results in Fig. 6(d) . Pure Γ s and Γ d scatterings correspond to blue and red curves, which enclose the physical region (shaded region in the inset), in which ρ s (0) is always roughly proportional to T c although not exactly. Therefore, this is a hallmark of the dirty BCS superconductors regardless of isotropic (Γ s ) or anisotropic (Γ d ) scatterings. Furthermore, we have also fix Γ s = 0.05 and 0.1 and tune Γ d to suppress T c . Both of them show the bending behavior near T c → 0, showing much better behavior than the separate scattering cases in view of the experimental data . The same as (a) but for the superconducting state conductivity σ1. The inset shows the result for isotropic scattering, revealing the universal conductivity in the zero frequency limit. (c) σ1/σ1n versus ν, with σ1n and σ1 from (a) and (b), respectively. The inset is for isotropic scattering. In (a)-(c), the line color is associated with the scattering parameter setting identically, and in (b) and (c), the insets share the same axis labels as for the main panels. (d) Best fits (lines) to the experimental data (symbols) with Tc = 17.5K measured at 22K (blue) and 1.6K (red) [13] . The dashed lines are best fits with Γs alone, and the solid lines include both Γs and Γ d . optical conductivity σ 1 (ν) can be obtained as [27] where f (ω) is the Fermi distribution function, G(ω, ε, θ) is the Green's function in the Nambu space: with σ 0,1,3 the identity and Pauli matrices. By numerical calculations, we obtain the frequency dependence of σ 1 (ν) in both normal and superconducting states (with ∆ 0 = 1 as the energy unit), as shown in Fig. 7(a) and (b) , respectively. In the normal state, as Γ d increases, σ 1 (ν) becomes more and more broadened to show long tail behavior. This is an important feature obtained in the "cold spot" model [21] . In the superconducting state, a nonzero Γ d breaks the universal conductivity [28] , which applies only in the case of pure Γ s scattering as shown in the inset of Fig. 7(b) . Furthermore, we renormalize the superconducting state σ 1 (ν) by the normal state value σ 1n (ν) to obtain the results shown in Fig. 7(c) . The resulted curves show similar behavior to the DOS at low frequencies: nearly linear νdependence if Γ d dominates which is quite different from the case of pure Γ s scattering (inset) with quadratic νdependence. This can be used to examine the types of scattering in future experiments. The most obvious feature when Γ d dominates is the behavior of σ 1 (ν) is no longer a simple Lorentian-like Drude peak (as in the case of isotropic scattering), since the quasiparticles experience angle dependent scattering rates. In fact, the conductivity should be an integration over a continuous distribution of Lorentzian components. Since the Lorentzian is sharper and higher for smaller scattering rates, the overall line shape of σ 1 (ν) should be steeper as the frequency approaches zero. Therefore, if the finite frequency data are used to perform the fit to a single Lorentzian-like Drude peak, as in the THz experiment, the extrapolated value of σ 1 (ν → 0) will be lower than the actual one, leading to superficial "missing" of the Drude weight. As an example, in Fig. 7 (d) we fit the experimental data with T c = 17.5K (symbols) in both normal (blue) and superconducting (red) states [13] . The best fits with only Γ s (dashed lines) are found to underestimate σ 1 (ν → 0) and the spectral weight, as compared to the result with both Γ s and Γ d (solid lines). Since the pairing gap ∆ 0 is not available in the experiment, we take it also as a parameter, and the fitted value is 1.11THz. More systematic fitting for a series of overdoped samples can be found in the Supplemental Material [26] , which supports our main point that apical oxygen vacancy increases with overdoping, leading to larger Γ d scattering. Summary and discussions.-In summary, we have found that the apical oxygen vacancies give rise to an anisotropic scattering rate Γ d cos 2 (2θ) which causes a clean-dirty dichotomy for low-high energy excitations. This provides a natural explanation of the anomalous behavior of the superfluid density and THz optical conductivity in the experiments. Therefore, we conclude that the superconducting states of overdoped cuprates are still captured by the BCS theory, as long as the anisotropic scattering rate is adequately considered. Finally, we make some remarks. (1) In this work, we have omitted many material-dependent details (such as specific tight-binding models, pairing interactions and interaction effects) in order to obtain universal results. The doping dependence enters the problem via the scattering rate Γ d , which increases with doping. (2) The anisotropic scattering rate has been reported in angledependent magnetoresistivity experiments in Nd-LSCO [29] and Tl 2 Ba 2 CuO 6+δ [30] . But to the best of our knowledge, there has not been such report in overdoped LSCO films. A systematic study of the anisotropic scattering rate can check our model and finally answer the question on whether overdoped cuprates are dirty BCS superconductors or not. (3) Our prediction suggests that it is necessary to extend the optical measurements down to even lower frequency (e.g. GHz) in order to obtain accurate behavior of the optical conductivity. (4) About superfluid density, there are two other kinds of scalings reported in literature: Uemura's law [31] ρ s (0) ∝ T c and Homes' law ρ s (0) ∝ σ dc T c [32] . The former is obtained in underdoped cuprates, where phase fluctuations are very strong and cannot be neglected. [33] . On the other hand, Homes' law can be explained by dirty BCS superconductors with pair-preserving impurities only [4, 34] . Here, in the overdoped LSCO films, the oxygen vacancies are pairbreaking, hence Homes' law does not apply. (5) Recently, the apical oxygen vacancies have also been observed in optimally doped YBa 2 Cu 3 O 7−x [35] . Then, according to our study, σ 1 (ν) at low frequency should be non-standard Drude-like, which may be consistent with the early microwave experiment [36] . Furthermore, if the apical oxygen vacancies are widely present in different families of cuprates, then the effect of Γ d scattering should be considered carefully in future studies. Anisotropic scattering rate: Born approximation As explained in the main text, since the impurity is out-of-plane, its intensity is expected to be quite small such that the Born approximation is acceptable. The self energy under the Born approximation [4] as shown in Fig. 5 is where G(k, iω n ) = (iω n − ε k ) −1 . As usual, for simplicity, we adopt the wide band approximation with circular Fermi surface and keep the form factor (cos k x − cos k y ) with only the angle dependence 2 cos(2θ), we obtain Σ(θ, iω n ) = 256κ 2 n imp N cos 2 (2θ) dε dθ 2π by analytic continuation, it gives the desired scattering rate Γ θ = Γ d cos 2 (2θ) with Γ d = 128πκ 2 n imp N . Together with the usual isotropic scattering rate Γ s , we get Eq. 2 in the main text. The retarded Green's function in the superconducting state can be written in the Nambu space as G −1 R (ω, k) = (ω + iΓ k )σ 0 − ε k σ 3 − ∆ k σ 1 , where σ 0,1,3 are identity and Pauli matrices. Keeping only the angle dependence and after integration over ε, we obtain The angle-dependent partial density of states (PDOS) ρ(ω, θ) is defined as the imaginary part of G R (ω, θ), hence, given by which is Eq. 3 in the main text. The angle integral of ρ(ω, θ) then gives the total density of states (DOS) ρ(ω) = dθ 2π ρ(ω, θ). It is interesting to notice that near the nodal region θ = π/4 − δθ, the d-wave pairing ∆ θ = ∆ 0 cos(2θ) ≈ ∆ 0 δθ vanishes linearly with δθ, while the Γ d -scattering vanishes quadratically Γ θ ≈ Γ d δθ 2 . Therefore, the Γ d -scattering has little effect on the low energy quasiparticles, which is quite different from the isotropic scattering and causes the clean-dirty dichotomy as we discussed in the main text. Within the BCS framework, we assume a pairing interaction where φ k = (cos k x − cos k y )/2, correspondingly φ θ = cos(2θ). Then, the gap function ∆ k = ∆ 0 φ k is determined by the self-consistent condition, which gives Eq. 4 in the main text by defining λ = V N . In the last line, we have also adopted the wide band approximation with circular Fermi surface. Next, by letting ∆ 0 → 0, we can determine T c . The calculation is in parallel to Ref. 23 . The above self-consistent equation becomes where we have added and subjected the same term (with Γ θ = 0) such that the summation of the first two terms converges and upper limit can now be pushed to infinity. The last term is the same as the clean case without disorder except T c0 replaced by T c . Using the clean case result T c0 = αΩe −2/λ for d-wave pairing, we obtain substituting it into Eq. 19, we obtain which is Eq. 5 in the main text. Next, let us consider the superfluid density. The calculation mainly follows Ref. 25 . The superfluid density is given by two diagrams as shown in Fig. 6 . The first diagram is from the paramagnetic current J P k = e∇ε k d † k d k and the second one is from the diamagnetic current J d k = e 2 ∇ 2 ε k d † k d k . Put these two contributions together, we obtain the superfluid density where we keep only the xx-component. Similar to previous sections, the k-summation is replaced by integrals over ε and θ under wide band approximation with circular Fermi surface, giving rise to which is Eq. 6 in the main text. The optical conductivity is contributed only by the paramagnetic current as shown in Fig. 7 . The current-current correlation function K(iν n ) is given by where the spectral representation has been used. Now the summation of iω n can be completed and gives After analytic continuation iν n → ν + i0 + , the optical conductivity σ 1 (ν) can be obtained as Finally, we replace the k-summation by energy and angle integrals and arrive at the final expression which is Eq. 7 in the main text. [13] in the normal state (a) and superconducting state (b), respectively, with (solid lines) and without (dashed lines) the anisotropic scattering rate Γ d . In the superconducting state, we fix ∆0 = 2Tc as a rough estimate of the pairing gap. For comparison, we also present the fitted results (dotted lines) with ∆0 = 0 as performed in Ref. 13 . In the inset of (a), we plot the best fitted parameters (Γs, Γ d ) versus Tc, respectively. On the basis of Eq. 27, we fit the experimental data of σ 1 from Ref. 13 . We read the data points from ν = 0.4 to 1.6THz by every 0.1THz. The best fits are shown in Fig. 8(a) for normal state and in Fig. 8(b) for superconducting state, respectively, with (solid lines) and without (dashed or dotted lines) Γ d . For the superconducting state, we have fixed ∆ 0 = 2T c as a rough estimate of the pairing gap. Compared with pure Γ s -fitting, the results using both Gs and Gd agree better to most of the data. As anticipated, the pure-Γ s fitting underestimates both σ 1 (ν → 0) and the spectral weight (enclosed area). We further show the best fitted parameters (Γ s , Γ d ) obtained from the normal state data in the inset of Fig. 8(a) , which shows clearly that Γ d becomes stronger upon overdoping (reducing T c ) while Γ s barely changes, supporting our main point that the apical oxygen vacancies become more and more important in overdoped LSCO upon doping. In order to determine the position of the oxygen vacancies theoretically, we performed first principle calculations on La 2−x Sr x CuO 4−y within the local density approximation (LDA), which has been found to work very well in the overdoped cuprates [37] . By constructing two kinds of super unit cells, 2 × 2 × 1 and 2 × 2 × 2, with one oxygen vacancy, and optimizing the positions of the doped Sr atoms, we obtain the total energy difference between two kinds of oxygen vacancies ∆E = E apical − E planar as shown in Fig. 9(a) . Without Sr doping, the planar oxygen vacancy E planar is lower than the apical one E apical , but upon Sr doping, their difference reduces and finally reverses at high doping levels. Such a tendency becomes stronger when we enlarge the super unit cell from 2 × 2 × 1 to 2 × 2 × 2, indicating the oxygen vacancies in overdoped LSCO are very likely to be apical rather than planar ones, in agreement with the Raman experiment [18] . This behavior can be roughly understood from the density of states (DOS), as shown in Fig. 9(b) . As can be seen, the planner oxygen vacancy reduces the DOS near the Fermi level at x = 0 and thus saves more energy. Instead, the apical vacancy enhances the DOS near the Fermi level. Upon hole doping caused by Sr, the DOS of the apical case is reduced and thus finally becomes preferred than the planar case. From both the experiment and our LDA calculations, we have found the oxygen vacancies are dominated by apical ones in overdoped LSCO. However, at small doping, the planar ones are more stable. Therefore, we may also ask what happens for the planar oxygen vacancies? At first, one planar oxygen vacancy breaks the Zhang-Rice singlets on its two neighboring sites. To describe such an effect, we consider a model by adding two onsite potentials given by where the subscript x, y labels the position of the vacancy is on the x-bond or y-bond, hence, d † 0 , d † x and d † y generate electrons on (0, 0), (x, 0) and (ŷ, 0), respectively. We call this scattering process as type-I. Its amplitude |V kk | 2 is given by Such a disorder scattering is difficult to treat exactly as it depends on the transferred momentum (k x,y − k x,y ) (but not k and k independently like the apical case) and breaks the rotational symmetry. Nevertheless, let us qualitatively consider its role on the d-wave superconducting properties. If we keep only the dominant scattering processes, they occur at k x ≈ k x (for x-bond vacancies)and k y ≈ k y (for y-bond vacancies). These scatterings preserve the pairing sign and are expected to be not pair-breaking, in some sense like the point disorder in s-wave superconductors. Therefore, the standard nonmagnetic disorder effects on the conventional dirty s-wave superconductors are expected: such disorders should have little effect on T c but reduce the superfluid density. Furthermore, these scatterings do not vanish for nodal quasiparticles, hence, the scattering rate Γ is finite in this region (like Γ s ) which cannot explain the linear T -dependence of the superfluid density. On the other hand, one planar oxygen vacancy can also change the hopping on this bond (called type-II process), which can be described by the following Hamiltonian The symbols are similar to the above. The amplitude of this scattering potential is |V II kk | 2 = 4V 2 cos 2 k x,y + k x,y 2 (31) which is largest at k x ≈ −k x (for x-bond vacancies) and k y ≈ −k y (for y-bond vacancies). These scatterings are also not pair-breaking and thus are expected to be similar to the above type-I scatterings. Put these two processes together, we conclude that the planar oxygen vacancies behave more like pair-preserving disorders in d-wave superconductors and cannot explain the superfluid density experiment. Now let us examine the effect of the current vertex correction caused by the anisotropic disorder scattering. We label the bare current vertex by v k and the dressed vertex by Λ k . They are connected in a self-consistent way as shown in Fig. 10 , where we have neglected all crossing diagrams which are at least of order n 2 imp κ 4 and are expected to be more relevant to localization [38] . In fact, the resistivity in the overdoped LSCO becomes smaller and smaller upon doping, indicating the system are more and more itinerant than localization [3] . From Fig. 10 , we can write down the self-consistent equation For a general scattering V kk , the above equation can be solved to give rise to an additional factor (1 − cos θ kk ) in the so-called transport scattering rate τ −1 tr which enters σ 1 (ν) [4] . But here, in our model, it becomes quite simple because the disorder scattering potential V kk is already factorized V kk = 16κ cos 2θ k cos 2θ k such that the summation over k can be completed exactly. Moreover, notice that Λ k is time reversal (k → −k ) odd, while |V kk | 2 G 2 (iω n , k ) is time reversal even, the integral over k must be zero, which means Λ k = v k . Therefore, in our disorder model, the current vertex correction totally vanishes, if the crossing diagrams can be safely neglected. Methods of quantum field theory in statistical physics (NJ Physica C: Superconductivity and its Applications npj Quantum Mater Exactly speaking, the frequency summation should be bounded with a soft-cutoff (with the form of a boson propagator) instead of the hard-cutoff used here. However, only the latter can be mapped to the momentum cutoff strategy as in the BCS treatment Introduction to Many-Body Physics See Supplemental Material [url] for calculation details and further discussions, which includes Refs In this supplementary material, we first derive the scattering potential caused by the apical oxygen vacancies and then the anisotropic scattering rate Γ d cos 2 (2θ) within the Born approximation. The effects of this scattering rate on some superconducting properties, including density of states, gap, T c , superfluid density, and optical conductivity, are derived in a self-contained way. Finally, some further discussions are given, including first principle calculations to determine the position of the oxygen vacancies (apical or planar), the effect of planar oxygen vacancies (if they exist), and the current vertex corrections. At first, suppose there is no oxygen vacancy, we have the following translational invariant Hamiltonianwhere d † i generates the in-plane electrons effectively at the Cu-site and p † i generates the ones at the apical O-site. The sign change of the out-of-plane hopping (with amplitude t o ) between ±x and ±ŷ are caused by the (x 2 − y 2 ) orbital character of the in-plane carriers (either d x 2 −y 2 orbital electrons or Zhang-Rice singlets), see Fig. 1 in the main text for clarity. We have omitted the spin indices here for simplicity because they are irrelevant to obtain the scattering potential. After integrating out the p-electrons (or by a second order perturbation), we obtainwhere κ = where the first term ε kd is from the pure d-electrons, and the second term comes from the coupling to apical oxygen orbitals. In together, ε k is taken as the free band without any disorder. Next, let us consider one apical oxygen vacancy at the origin. Relative to the translational invariant background, the single vacancy is expressed by the impurity Hamiltonianwhose Fourier transformation then gives a scattering potentialwhich is Eq. 1 in the main text.