On the use of Monkhorst–Pack scheme to evaluate superconductivity and the issue of umklapp electron–phonon interactions

Abstract

The Monkhorst–Pack scheme is a method to save time in the days of slow computers. It excludes umklapp phonons with significant consequences. Its widespread application to evaluate superconductivity arises from the desire to reduce phonon contributions to solve a historical difficulty of the BCS theory. An alternative method turns out to be more accurate in Pb and Pd.

---

I. Introduction

In 1957, Bardeen, Cooper and Schrieffer published the famed BCS theory. It leads to a formula to relate T_D (Debye temperature) and T_c (transition temperature of superconductors).¹ In 1968, Ashcroft predicted metallic hydrogen (T_D ∼ 3500 K by estimation) might be a high-temperature superconductor.² In 2004, Ashcroft recommended hydrogen-dominated metallic alloys as a better alternative because the attainment of metallic states should be well within the current capabilities of diamond anvil cells but at pressures considerably lower than may be necessary for hydrogen.³

In 1970, Satterthwaite and Toepke found 8.05 ≤ T_c ≤ 8.35 K in thorium–hydrogen and thorium–deuterium, surprisingly high in the 1970s. It was in a Pyrex enclosure in ambient pressure, with no detectable isotope effect.⁴ In 1972, Stritzker and Buckel found T_c up to 9 and 11 K in palladium-hydrogen and palladium–deuterium, isotope effect reversed.⁵ In 2015, Drozdov and colleagues found T_c = 203 K in the sulfur hydride system in a diamond anvil cell under pressure.⁶ In 2020, Dias and colleagues found T_c ≃ 287 K in a carbonaceous sulfur hydride under pressure,⁷ retracted by the journal investigating if its style affected the result.

In 1972 and 1978, Ganguly and Papaconstanstopoulos investigated the reversed isotope effect in palladium hydrides theoretically.^8,9 In 2020, Semenoka and colleagues reported a rule to predict the maximum T_c for metal hydrides.¹⁰ Meanwhile, Flores–Livas and colleagues offered an up-to-date perspective on conventional superconductivity in hydrides at high pressures.¹¹ Very recently, Peng,¹² Chen,¹³ Dou,¹⁴ Duan,¹⁵ Li¹⁶ and colleagues published further results of theoretical studies of superconductivity in compressed hydrides.

In light of the fascinating progress outlined above, we need to consider if the Monkhorst–Pack (MP) scheme is applicable to evaluate the BCS theory. It is a method to save time in the days of slow computers^17–19 readily available from popular computer packages such as QUANTUM ESPRESSO.²⁰ It resembles Fourier expansion, but in terms of the so-called “stars”, representing groups of plane waves propagating in the symmetric directions of the lattice. It serves the purpose of “integrating periodic functions of a Bloch wave vector over either the entire Brillouin zone (BZ) or over specified portions”.¹⁹ Therefore, by design, it excludes umklapp phonons.

We do not know if anyone has ever extended the MP scheme to include umklapp phonons. On the other hand, in 1994, Savrasov, Savrasov and Andersen stated without ambiguity that, when they calculated electron–phonon coupling strengths and transport properties in Al, Nb and Mo, the phonon integration is over BZ.²¹ When applied to a semiconductor, their “tetrahedron” method is identical to the MP method.²² They treated Ta, Pb, Mo, V, Cu and Pd similarly.²³ Camargo–Martínez and colleagues used the MP scheme to study superconductivity in H₃S,²⁴ stating they had knowingly excluded umklapp contributions.^25,26 Du and colleagues used the MP scheme to study transition metallic hydrides²⁷ and made a similar statement.²⁸

We believe the practice in^21–28 amounts to an answer to the call for reducing phonon contributions to superconductivity to solve a historical difficulty. In the BCS theory, superconductivity arises from the interactions between electrons and phonons.¹ In the Drude theory, electrical resistivity arises from the same interactions.²⁹ In both cases, the e–p interactions do not depend on temperature, leaving other considerations to induce temperature dependence of the phenomena. Many tried but failed to find such e–p interactions: whenever the interactions were reasonable in the normal state, they always were too strong in the superconducting state, and vice versa.^30–42 In 1977, Peter, Ashkenazi and Dacorogna suggested multiplying the electron–phonon coupling constants by 0.7 in superconductors,³⁶ calling for a reduction of phonon contributions.

But why do we exclude the umklapp phonons? On account of the availability of the MP scheme? On account of the popularity of the practice? We need an adequate reason. In a superconductor, two electron states in a pair have identical occupation probability by symmetry. But they may also have different occupation probabilities when they accommodate electrons driven by normal and umklapp phonons. We devise an alternative scheme to pair superconducting electrons to avoid conflict.^43,44 We show it can solve the historical difficulty with high accuracy in Pb and Pd.

We arrange the article as follows. In Sections II to IV, we detail the issue of umklapp phonons and the MP scheme. We explain why the scheme excludes umklapp phonons. We describe the method of Carbotte and Dynes to include umklapp contributions. We also describe the pioneering procedure by Savrasov and Savrasov (S & S) to implement the MP scheme to evaluate superconductivity, together with its recent advances. In Sections VI to VIII, we discuss the historical and theoretical difficulties that oblige us to reduce phonon contributions to superconductivity. In Sections IX to XIII, we discuss the alternative pairing scheme to reduce phonon contributions. We demonstrate its high accuracy in Pb and Pd via detailed calculations. We place brief discussions and conclusions in Sections XIV to XVII.

II. Umklapp phonons: MP scheme

The MP scheme is an alternative to Fourier expansion to present the states of electrons or phonons in BZ in terms of the ‘stars’. In an article in 1973, titled “Mean value point in the Brillouin zone”,¹⁷ Baldereschi proposed the formulationwhich imitates the usual Fourier series for an arbitrary function in reciprocal space, f(q), with the harmonic waves in the Cartesian coordinates grouped into waves propagating in the symmetric directions of the lattice of a crystal. In eqn (1), we havewhere R₀ and R_m are coordinates of the atomic site and an atom on the m-th shell, m = 1, 2, 3, …, R_m runs over the m-th shell, and N_m number of atomic sites on that shell. In eqn (2)R_m − R₀ represents a set of vectors resembling rays from a point source. It customarily is referred to as a star. Baldereschi found the mean value point, q₀, that ensures A_m(q₀) = 0 for m ≤ 2 or 3. He argued eqn (1) could be reduced to f(q₀) ≃ f₀ because other terms were negligible in a rapidly converging series. The treatment simplifies phonon integration significantly to a surprisingly good effect.¹⁷

Shortly afterwards, also in 1973, Chadi and Cohen¹⁸ proved the following relation:

with δ_mn = 1 or 0 depending on if m = n, closely resembling the orthogonal relation in Fourier analysis. It is rather regretful that Chadi and Cohen did not carry on to find the expression of f_m in eqn (1). Instead, they replaced the mean value point of Baldereschi with a grid of “special points” to integrate f(q), see ref. 18 for further details of the algorithm.

In 1976 Monkhorst and Pack¹⁹ found the expression

and used it to supplement eqn (2) and (3) to establish a well-defined self-consistent system comparable with Fourier expansion. Monkhorst and Pack also proved:which leads through eqn (1), (2) and (4) to a discrete system comparable with the system of fast Fourier transform (FFT). In the Mokhorst–Pack scheme, f(q) in eqn (1) is evaluated over an uniform “coarse mesh” of q, or sample points, to reduce the load of computation. Afterwards, f(q) is evaluated over a “fine mesh” by interpolation or extrapolation.

The MP scheme is for “integrating periodic functions of a Bloch wave vector over either the entire Brillouin zone (BZ) or over specified portions”.¹⁹ Extending the MP scheme to include umklapp phonons will be challenging because we do not know the lattice vector involved in the umklapp process. In Fig. 1, we use the filled squares to mark the MP sample points and the dashed lines to mark the borders of BZ. We show the phonons are periodic only when their momentum is in the direction of the lattice vector. It will not be easy to determine the umklapp phonons under other circumstances.

Fig. 1 Normal and umklapp lead phonons, shaded triangle envelops Debye frequencies and momenta, vertical dashed lines mark borders of BZ. Top: In the (0.77, 0.25, 0.59) direction. Bottom: In the (0, 1, 1) direction, filled squares represent Monkhorst–Pack sample points.

III. Umklapp phonons: CD sphere

To illustrate the difficulty to extend the MP scheme, we present electron–phonon interactions schematically in Fig. 2, where the shaded circles are for spherical Fermi surfaces in the extended zone scheme, Γ the initial state of an electron. We also use A, B, C, D, E and F to present the end states of the electron when scattered by a phonon of momentum q ≤ q_D, q_D being the Debye momentum. Since Γ, A, and B are in the first BZ, we refer to the events as normal scatterings, others umklapp scatterings.

Fig. 2 Top: Fermi spheres (shaded circles, radius = k_F, first BZ in the S–W corner) and phonons on a spherical surface. An electron is scattered from Γ to A, B (normal scattering) and C, D, E, F (umklapp scattering, folded back to CC, DD, EE, FF in the first BZ). Bottom: Carbotte–Dynes phonon sphere (centre = Γ, radius = 2k_F).

We use CC, DD, EE, FF to mark images (replicas) of C, D, E, F, respectively, in the first BZ. A state and its image, say C and CC, have identical coordinates against the centres of the Fermi surfaces they are on. In the previous paragraph, we say a normal phonon drives an electron from Γ to C. Alternatively, we can say an umklapp phonon drives the electron from Γ to CC. Either way, the result is the same in terms of the strength and enumeration of the electron–phonon interactions.

It is apparent that, as long as we have normal phonons to scatter an electron from Γ to C, D, E and F, we can think we have umklapp phonons to scatter the electron from Γ to CC, DD, EE and FF. It is wrong to believe the destination states of the electron are restricted to A and B and the arc in between and conclude that just a portion of the Fermi surface is subject to bombardment by the scattered electrons.

To evaluate the strength of the electron–phonon interactions, we have to count the number of phonons. We have two options. We may choose to follow the open circle in the top part of Fig. 2 to trace the electron states, exemplified by A, B, C, D, E and F, and integrate the phonons. The radius of the circle must vary between 0 and q_D. We must keep identifying the electron states from several BZs. It will be a cumbersome and impractical task. Alternatively, we can choose to follow Carbotte and Dynes to integrate A, B, CC, DD, EE and FF, …, in a large sphere enveloping the entire Fermi surface,⁴⁵ shown in the bottom part of Fig. 2, referred to as the Carbotte–Dynes sphere, or CD sphere for short.

IV. Umklapp phonons: S & S procedure

The procedure of Savrasov and Savrasov (S & S) is an early attempt to implement the MP scheme to evaluate superconductivity and electrical resistivity. It starts from the well-established framework of density functional theory (DFT) to find phonons by solving the following matrix eigenvalue equation:for lattice vibrations. Eqn (6) is identical to eqn (9) in ref. 22 (ψ_m, ψ_n, A_n and ε marked somewhat differently). It assumes a Fourier expansion
where |k〉 is a wave function of the system (marked as |k_j〉 originally). It evaluates ε and V_eff numerically.²²

In 1996, S & S investigated eight metals in ambient.²³ Nowadays, DFT calculations can predict a significant portion of the properties of a given material. Given a modest set of elements, it is already possible to predict the thermodynamically stable crystal structures, not only in ambient but also spanning a vast pressure range to terapascals and beyond.⁴⁶

With the next step, S & S enters untested waters because it is about a “new, generally applicable method for ab initio calculation of the wave-vector dependent electron–phonon coupling”.²¹ It evaluates

g_{[small script l]}(q) = 〈k + q|δV|k〉 + … (7)

for the matrix element of electron–phonon coupling (interaction). Eqn (7) is identical to eqn (3) in ref. 21, save some presentational changes for clarity, [small script l] identifies phonon polarization, k and q stand for electron and phonon momenta. To the best of our knowledge, the authors of ref. 21 were the first to let V in eqn (7) be identical to V_eff in eqn (6). It is a significant advance because a pseudopotential is a phenomenological and task-specific entity, not universally applicable. For example, the empty-core model of Slater⁴⁷ and the muffin-tin model of Ashcroft⁴⁸ are entirely different and not supposed to be interchangeable.

The final step of the S & S procedure is to integrate the matrix element in eqn (7). They found

which is identical to eqn (1) in ref. 21 but in somewhat different notations. In eqn (8)ν is the phonon energy, N number of atoms in unit volume, ε and ε_F electron and Fermi energies, respectively, ω_{[small script l]} phonon frequency. If we let S(q) = 1 or 0 depending on if q is in the first phonon BZ, then eqn (8) becomes eqn (1) in ref. 21 in every detail. It amounts to a numerical scheme to exclude all the umklapp phonons. We will discuss the physics of α²F(ν) and phonon contributions to it in some detail in Section IX, with S(q) from a different scheme.

S. Y. Savrasov declared that, when applied to a semiconductor, his “tetrahedron method” is identical to the special-point scheme of Monkhorst and Pack.²² It is not clear whether or not S & S had adopted the MP scheme in ref. 23 but Fig. 1 there tells us they probably had employed a 4 × 4 × 4 MP grid to sample the phonons from within BZ for integration.

In Fig. 3 we use the curves to present the S & S α²F(ν) in Pd and Pb sampled from Fig. 1 in ref. 23. We also use the shaded silhouettes to present more accurate examples for comparison. The S & S evaluations are overly strong in both metals. We will further the comparison in the following sections.

Fig. 3 Electron–phonon spectral densities, α²F(ν), in Pd and Pb. Top: From S & S calculation (curve, λ = 0.340) and the alternative scheme (silhouette, λ = 0.199). Bottom: From S & S calculation (curve, λ = 1.55) and McMillan experiment (silhouette, λ = 1.33).⁴⁹ The S & S evaluations are overly high in both metals.

V. Umklapp phonons: recent advances

In 2019, Camargo-Martínez, González-Pedreros and Baquero studied the pressure effect on high-T_c superconductivity in H₃S.²⁴ In response to our enquiry in ref. 50, they confirmed that “the electron–phonon umklapp process has a vital relevance in explaining the behaviour of electrical resistance” but claimed “in the total absence of electrical resistance, …, perhaps an electron that undergoes the umklapp-process will contribute to electrical resistance rather than the superconducting phase”.²⁵ In short, they believe the umklapp process does not contribute much to superconductivity.

In a following publication,²⁶ González-Pedreros, Camargo-Martínez and Mesa continued to apply the MP scheme to evaluate superconductivity in D₃H. They cited,^25,50 claiming that “the presence of umklapp process is mainly associated with the electrical resistivity, which is measured in the normal state and not in the superconductor one”. They announced “Our calculations were considered in the superconducting state, and Umklapp processes have not been included”.

In 2021, Du and colleagues studied high-temperature superconductivity in transition metallic hydrides under pressure.²⁷ They stated that the umklapp phonons “play a key role in explaining the DC conductivity of normal metals”.²⁸ In response to our enquiry,⁵¹ they claimed: “the pairing of electrons with opposite momenta, caused by an exchange by virtual optical phonons depends on umklapp phonons if k_F ≪ q_D.” They also claimed: “if these quantities are of the same order of magnitude, or k_F ≫ q_D, then pairing occurs without a noticeable contribution of the umklapp phonons”.²⁸ Interestingly, we have k_F ≥ 1.145q_D when valence ≥3, true in virtually all metallic superconductors, proving the MP scheme is inapplicable.

VI. Pairing difficulty: history

In theory, there is no difference in the interactions between electrons and phonons in the superconducting and non-superconducting states. In history, however, whenever reasonably accurate superconductivity arisen from the first-principles calculation, electrical resistivity always turned out to be too weak and vice versa.

In an early example, Tomlinson and Carbotte in 1976 evaluated α²F(ν) for lead with the pseudopotential of Appappillai and Williams, which is a specification of the Heine–Abarenkov potential.³⁰ The agreement with tunnelling measurements was considered to be “very good”⁵² although the longitudinal numerical peak of α²F(ν) at 9 meV is more than 2 times stronger than the observed tunnelling peak. It is not clear if correct T_c would arise from the first-principles α²F(ν).

Soon after, in 1977, Tomlinson and Carbotte evaluated the electrical resistivity, ρ(T), for lead with the same potential.³¹ Between T = 4 and 295 K, their ρ(T) is about 75% of the observed values, as is shown clearly in the graphic portrayal of the result in 1981 by Eiling and Schiling.³² A similar discrepancy occurred in aluminium (Heine–Abarenkov potential tabulated by Harrison, discrepancy in electrical resistivity obscured by logarithmic scaling).^33–35

Also in 1977, Peter, Ashkenazi and Dacorogna studied the effects of anisotropy on the T_c of niobium.³⁶ Intriguingly they found that the electron–phonon coupling constants determined are probably too large and have to be multiplied by a factor of 0.7 (0.49 in coupling strength) to obtain the observed T_c.

Again in 1977, Harmon and Sinha evaluated α²F(ν) for niobium.³⁷ They borrowed the muffin-tin potential used in the calculation of the band structure. They found λ = 1.58, which is considerably larger than the commonly accepted values (∼1.0).

In 1979, Glötzel, Rainer and Schober evaluated α²F(ν) for vanadium, niobium, tantalum, molybdenum, tungsten, palladium, platinum and lead.³⁸ They carefully avoided any uncontrolled approximations. To find phonon dispersion, they used published Born-von Kármán fits to the force constants. To estimate the strength of the electron–phonon interaction, they adopted muffin-tin potentials developed for band structure calculations. The value of T_c turned out to be 2 to 3 times too high. They show T_c is sensitive to the Coulomb pseudopotential, μ*, but stick to a reasonable choice μ* = 0.13. They conclude their careful approach was incapable of reproducing the observed values of T_c.

In 1987, Al-Lehaibi, Swihart, Butler and Pinski evaluated both ρ(T) and α²F(ν) for tantalum with a muffin-tin potential from band calculation.³⁹ While ρ(T) was found to be slightly lower than experimentally observed, α²F(ν) exceeded the tunnelling values significantly, giving T_c = 7.01 K (4.5 K experimentally). This was considered a puzzle. A similar puzzle occurred when niobium was investigated.^40–42

VII. Pairing difficulty: inconsistent occupancy

There is a paradox disallowing contributions to superconductivity by some phonons. The topic was discussed before.^43,44 We revisit it here, with additional explanations, for the convenience of the reader. In Fig. 4 we show the Fermi surface in the first BZ in Fig. 2, and a pair of electron states on the surface, A = (k′, ↑) and FF = (−k′, ↓). Initially we have an electron in the state Γ = (k, ↑). Since an electron has two spins, we may also have Γ = (k, ↓).

Fig. 4 Schematic of a spherical Fermi surface and a pair of electrons with initial momenta k and −k, and destination momenta k′ and −k′. In normal and umklapp scattering the pair follows the paths indicated by the solid and dashed arrows respectively. The range of normal scattering is marked by the shaded wedges (half angle = 60°, valency = 2).

Letting Γ = (k, ↑), we have a normal phonon to drive the electron into A = (k′, ↑). Note the spin remains unchanged before and after scattering. By symmetry we also have a normal phonon to drive another electron from (−k, ↓) to (−k′, ↓). In Fig. 4 we show the paths of the two electrons as the two curved solid arrows, which are in association with a certain occupation probability of the destination pair state, say h(q) = 0.3, q = k − k′ being the momentum of the normal phonon.

Letting Γ = (k, ↓), we have an umklapp phonon to drive the electron into FF = (−k′, ↓). By symmetry we also have an umklapp phonon to drive another electron from (−k, ↑) to (k′, ↑). In Fig. 4 we show the paths of the two electrons as the two curved dashed arrows, which are in association with another occupation probability, say h(q) = 0.5, q = k + k′ being the momentum of the umklapp phonon, in contrast to the previous probability 0.3 in association with q = k − k′.

We should notice we have followed Mott and Jones to assume we always have a phonon to scatter an electron between any two states in a shell surrounding the Fermi surface.²⁹ For example, in Fig. 4, the states Γ, A and FF are all at the Fermi level. Another example is eqn (8), where we always have ε_F as the electron energy due to the delta function δ(ε − ε_F), assuming we have phonons to let δ(ν − ħω_{[small script l]}) = 1 at any values of q requested by the events.

VIII. Pairing difficulty: double occupancy

To supplement our discussion in the previous section, we follow BCS to construct the ground state made from paired valence electrons.¹ We haveto present states in the events shown as the solid and dashed arrows in Fig. 4, respectively, both start from k. Note the spins of the electrons remain unchanged before and after scattering. We also have:to present the states in the events shown as the solid and dashed curved arrows starting from −k. Note that the normal destination state in eqn (9) is identical to the umklapp destination state in eqn (10). The normal destination state in eqn (10) is identical to the umklapp destination state in eqn (9). Lettingbe pair generation operators to generate particles in the destination states in both eqn (9) and (10). The ground state in the BCS theory can be written as¹where h_N is the occupation probability of the destination state of normal scattering in eqn (9), h_U occupation probability of the destination state of umklapp scattering in eqn (10). In eqn (12), we generate just two pairs of quasi-particles explicitly, others written collectively as |Φ〉. Due to double Fermion occupations we haveleading through direct calculation towhich is not normalized unless h_N or h_U vanishes, assuming 〈Φ|Φ〉 = 1. Eqn (14) arises because the exclusion principle does not allow two fermions to occupy the same state. The values of h_N and h_U make little difference. It deems the scenario described in eqn (9) and (10) a forbidden event when we have say 〈Ψ|Ψ〉 = 0.9999 or 1.0001. This difficulty always will be there unless normal and umklapp phonons are arranged not to share the same range of destination states.

In eqn (13), generates two electrons scattered by normal and umklapp phonons respectively. We assume equal energy for the two electrons, assuring double occupancy of the electron state. We have no other choice because, by the algebra of the BCS theory, the electron states are always at the Fermi level due to the delta function δ(ε − ε_F) in eqn (8). The problem is usually latent but outstands in the pairing events. Many opted to exclude umklapp phonons to deal with it.^21–28 We show it enlightens us to find an alternative option in the next section.

IX. Alternative pairing: scheme

We letin eqn (8) to identify the phonons allowed to pair electrons. In eqn (15)θ is the angle of scattering (angle between k and k′ in Fig. 4) and θ_N maximum angle of normal scattering (half angle of the shaded wedges). For Fermi electrons and Debye phonons we have θ_N = 2 sin⁻¹(q_D/2k_F), where k_F is the Fermi momentum, q_D/2k_F = (4Z)^−1/3 and Z valency, giving θ_N = 78.1°, 69.0°, 51.8° for Z = 1, 2, 3 and so forth, assumed to equal the number of conducting per atom for simplicity. Experimentally, Z from measurement may be fractional, and θ_N depends on the direction of q because the phonon BZ is a polyhedron.

Eqn (8) and (15) exclude all the normal phonons, avoiding the paradox discussed in Sections VII and VIII from taking place. In addition, they exclude the umklapp phonons capable of driving the electron into the other hemisphere of the Fermi surface, avoiding a similar paradox. We understand eqn (8) and (15) as a reasonable proposition to be tested numerically.

In the BCS theory, most Cooper pairs are in a thin layer across the Fermi surface. It is out of the exclusion principle because the internal states are all occupied, leaving no place to accommodate scattered electrons for pairing. Now eqn (15) draws a map over that thin layer, telling us graphically where the scattered electrons can land without violating the exclusion principle.

X. Alternative pairing: phonons

In Fig. 5, we use the open circles to represent vibration frequencies of the Pd (top) and Pb (bottom) atoms excited by neutrons along the symmetric directions of the (fcc) lattice.^53,54 Theoretically, in simple metals, the Born-von Kármán theory manifests itself as a series of 3 × 3 matrices, with elementsevaluated over the atomic shells surrounding an atomic site (YX, ZX and ZY defined by symmetry). In eqn (16)r is the normalized radius of an atomic shell, r = (R_m − R₀)/a, R₀ and R_m are defined in eqn (2), r² = x² + y² + z², and a crystal constant. In eqn (16) and , all in joules or eV, arising from the Coulomb force in the radial and transverse directions respectively.⁵⁵

Fig. 5 Experimental and theoretical phonon dispersion for palladium (top) and lead (bottom), open circles from neutron scattering measurement, curves from Born–von Kármán theory, usual (fcc) reciprocal site conventions.

We find and with the help of the popular method of optimization from Hooke and Jeeves.⁵⁶ In the pattern-search stage of the method, we perturb the values of and to evaluate the penalty function measuring the r.m.s. difference between the theoretical and experimental phonon frequencies shown in Fig. 5. We refer to a perturbation as positive if it reduces that r.m.s difference. Otherwise, we refer to it as negative. In the pattern-move stage, we implement positive perturbations in the original directions, negative perturbations in opposite directions, magnitude proportional to their effects. We repeat the process until we can no longer improve the fitting between the circles and curves in Fig. 5.

We let V_eff be third-order polynomials between the shells and second-order polynomials within the first shell, with when r = 0. We list the values of V_eff and its derivatives in Pb in Table 1. We have , 0.020, 0.637, 0.894, 0.150, −1.137 and , −3.063, −2.038, 2.932, −0.094, −0.078 on the first six atomic shells in Pd, a = 0.389 nm.

Table 1 Pb lattice potential data^a

m	Site^b	N m	r × a	V eff
a V eff and its derivatives in eV, a = 0.495 nm. b One sample site per atomic shell, in 0.5a.
0	(0, 0, 0)	1	0.000	−2.868	—	—
1	(1, 1, 0)	12	0.350	−0.599	2.138	−17.062
2	(2, 0, 0)	6	0.495	0.034	0.002	−1.839
3	(2, 1, 1)	24	0.606	0.028	0.030	−0.175
4	(2, 2, 0)	12	0.700	0.001	−0.193	−2.417
5	(3, 1, 0)	24	0.783	−0.013	0.191	−0.173
6	(2, 2, 2)	8	0.857	0.029	0.000	−0.542
7	(3, 2, 1)	48	0.926	0.035	0.031	−0.050
8	(4, 0, 0)	6	0.990	0.040	−0.001	−0.101
9	(4, 1, 1)	36	1.050	0.039	0.000	0.420
10	(4, 2, 0)	24	1.107	0.009	−0.198	−0.355
11	(3, 3, 2)	24	1.161	−0.024	0.000	−0.002
12	(4, 2, 2)	24	1.213	−0.024	−0.005	−0.123
13	(4, 3, 1)	72	1.262	−0.007	0.124	0.020
14	(5, 2, 1)	48	1.356	0.027	0.000	0.001
15	(4, 4, 0)	12	1.400	0.027	0.000	0.000
16	(4, 3, 3)	48	1.443	0.029	0.017	−0.011
17	(4, 4, 2)	30	1.485	0.032	0.004	−0.193
18	(5, 3, 2)	72	1.526	0.017	−0.135	0.168
19	(6, 2, 0)	24	1.565	0.001	−0.008	−0.055
20	(5, 4, 1)	48	1.604	0.000	−0.012	0.117

---

In Fig. 6 we use the vertical lines to present the locations of the shells. We also use the shaded silhouette to present the pseudopotentials for the electron–phonon interactions. Although the curves and silhouettes are significantly different, especially within the first shell, we should note the values of V_eff between the shells do not affect the phonons. We can let, for example, V_eff ≡ 0 within the first shell (empty core model).⁴⁷

Fig. 6 Pseudopotentials for lattice vibrations, V_eff in eqn (16) (curves), and electron–phonon interactions, V in eqn (20) (silhouettes), in palladium (top, V₀ = 31.6 eV and r₀ = 0.324 nm) and lead (bottom, V₀ = 29.0 eV and r₀ = 0.397 nm for the silhouette),⁴³ vertical lines indicate radius of the atomic shells.

XI. Alternative pairing: resistivity

Next, we evaluate ρ(T) for electrical resistivity against temperatures. The formula involves a pseudopotential, V, for the interactions between phonons and electrons. We will use the experimental values of ρ(T) to determine V and compare it with say V_eff in eqn (16). We start with the Drude formula:where T represents temperature in Kelvin, N number of atoms in unit volume, m_e and e electron mass and charge respectively.⁵⁷ In eqn (17)τ represents the collision (relaxation) time,x = ν/k_BT, k_B being the Boltzmann constant. We also have the electron–phonon transport densitywhere ω_{[small script l]} is the circular frequency of the phonon, [small script l] identifies polarization, ε and ε_F electron energy and Fermi energy, respectively, q = k′−k, k and k′ stand for the initial and end states of an electron.³¹

In eqn (19) we have the matrix element:⁵²

which is the equivalence of eqn (7) in the many-particles formalism, M being the mass of the atom, V = V(r) atomic potential, and r distance from the atomic site in real space. The atomic potential is differentiated to measure via Taylor expansion the first-order effect of atomic vibration on the valence electrons. We assumewhen 0 ≤ r ≤ r₀, otherwise V(r) = 0 (muffin-tin model).⁴⁸ We will adjust V₀ and r₀ to fit the output of eqn (17) to experimental data.

In metals, the electron–phonon interactions take place just on the Fermi surface due to the exclusion principle. We convert the summation in eqn (19) into integration. We let dk′ = (dk′/dε)dε to let the integration accommodate the delta function, δ(ε − ε_F), to implement the physics. Sine k′² sin(θ)dθ = qdq = q(dq/dω_{[small script l]})dω_{[small script l]}, the summation in eqn (19) also accommodates δ(ν − ħω_{[small script l]}), where θ is the angle between k and k′ in Fig. 4. The output of eqn (19) depends on the direction of k. We average it over the Fermi surface and find the analytical expression for α_tr²F(ν) in ref. 43.

We perform a two-dimensional search to specify V(r) in eqn (21) to minimize the difference between the experimental and computational values of ρ(T) illustrated in Fig. 7. We find

in Pb, compared with V_eff in Table 1. We find V₀ = 31.6 eV, r₀ = 0.324 nm in Pd. We also illustrate the results in Fig. 6 for a comparison with V_eff in eqn (16).

Fig. 7 Electric resistivity from experimental measurement (circles), computation (solid curves) and S & S evaluation (filled squares). Top: In palladium, deviation = 0.26% on average against resistivity in lead at 300 K. Bottom: In lead, deviation = 0.14% on average.⁴³

XII. Alternative pairing: superconductivity

We test if our pairing scheme in eqn (15), which excludes an extensive range of phonons, ensures accurate superconductivity. We will perform a two-dimensional search to specify V(r) in eqn (21) to minimize the difference between the experimental and computational values of the superconducting tunnelling current shown in the bottom part of Fig. 8. We will compare the resulting values of V₀ and r₀ with the values in eqn (22).

Fig. 8 Superconductivity in Pb, from the alternative pairing scheme in eqn (15) to curtail phonon contributions, Δ₀ adjusted to ensure an accurate value of T_c. Top: α²F(ν) from experiment (silhouette) and computation (histogram). Bottom: Tunnelling conductance from experiment (open circles) and computation (small circles).

In the top part of Fig. 8, the shaded silhouette represents the experimental α²F(ν) from ref. 49 which arises from the measured values of the superconducting tunnelling current shown in the bottom part of the figure. Our numerical α²F(ν), shown as the histogram, arises from eqn (8) which in turn arises from S(q) in eqn (15) and g_{[small script l]}(q) in eqn (20)via V(r) in eqn (21). It enables us the solve the Eliashberg equation at T = 0 numerically provided we know the edge of the superconducting gap function, Δ₀.⁵² We vary the trial value of Δ₀ until we find the experimental value of T_c from our numerical solution to the Eliashberg equation at T > 0.

Eventually, we find the computational tunnelling conductance current shown as the solid curve in the bottom part of Fig. 8. We have λ = 1.464 for our numerical α²F(ν), compared with λ = 1.33 ± 0.02 in ref. 49. We find μ* = 0.0994 and Δ₀ = 1.442 meV, giving 2Δ₀/k_BT_c = 4.654, compared with μ* = 0.1 ± 0.02 and Δ₀ = 1.358 ± 0.004 in ref. 49. We find

which deviate little from the values in eqn (22) for electrical resistivity.

We find that, within the 0.3% error margin between the experimental and computational values of ρ(T), V₀ and r₀ in eqn (22) are not unique. On the other hand, we always find unique values of V₀ and r₀ in eqn (23) to fit values of the tunnelling currents, full of rich details about the phonons. Therefore, we find V₀ and r₀ in eqn (22) iteratively. We treat V₀ and r₀ from eqn (23) as trial values of their counter parts in eqn (22) and proceed to fit ρ(T) once again. Usually, one iteration is sufficient to secure proper values of V₀ and r₀ in eqn (22).

In the top part of Fig. 3, we use the curve to show the S & S α²F(ν) in Pd, λ = 0.35 (0.34 by our calculation).²³ It leads trough the Eliashberg equation at T = 0 to Δ₀ = 0.0286 meV when μ* = 0.13. Consequently, we find T_c = 0.19 K from the relation 2Δ₀/k_BT_c = 3.52, that is Pd becomes a superconductor. In comparison, with the umklapp pairing scheme in eqn (15), we find a different α²F(ν) shown as the shaded silhouette in the top part of Fig. 2, λ = 0.199. It leads through the Eliashberg equation to Δ₀ = 0, giving T_c = 0.

XIII. Alternative pairing: numerical scheme

The scheme in eqn (15) arises entirely from the theoretical point of view to avoid the pairing difficulty revealed and analysed in Sections VII and VIII. Now we test a pairing scheme arising from a numerical point of view. We will include or exclude phonons depending on if the output of the Eliashberg equation matches experimental data accurately.

In Eqn (8), we let S(q) = S_i when θ = θ_i, i = 0,1,2,…,20, θ_i being mesh points evenly distributed between θ = 0 and π/2. We always force S_i = 0 when θ > π/2. We follow the procedure of Hooke and Jeeves⁵⁶ to vary the values of S_i to minimize the difference between the experimental and computational values of the superconducting tunnelling conductance, shown as the small circles and the curve, respectively, in the bottom part of Fig. 9. We force S_i = 0 when θ < θ_N initially but let the procedure run free afterwards.

Fig. 9 Superconductivity in Pb from the alternative pairing scheme, with numerical S(q) in eqn (8) and same g_{[small script l]}(q) in eqn (8) and (19), conventions identical with Fig. 8.

We evaluate the matrix element in eqn (20) with the pseudopotential in eqn (21) and (22) so that the strength of the electron–phonon interactions is identical to that for electrical resistivity. We find the computational tunnelling conductance current shown as the solid curve in the bottom part of Fig. 9. We have λ = 1.726 for the numerical α²F(ν) shown as the histogtam in the top part of Fig. 9. We find μ* = 0.151, Δ₀ = 1.458 meV, giving 2Δ₀/k_BT_c = 4.706.

In Fig. 10, we use the solid curve to present S(q) for our numerical pairing scheme. It takes the value of S_i when θ = θ_i. Otherwise, we define its value by interpolation. It fits well the theoretical values of S(q) in eqn (15), shown as the shaded rectangle, apart from some minor disturbances between q = 0 and q_D. The slight slope of the curve around q_D is from the non-spherical geometry of the phonon BZ.

Fig. 10 Survival rate, S(q), in lead from theory (silhouette) and numerical calculation (curve).

XIV. Discussion: palladium

We find λ = 0.199 in Pd, α²F(ν) shown as the shaded silhouette in the top part of Fig. 3, giving T_c = 0 in accord with the experiment. We are not alone. In 1977, Papaconstantopoulos and co-workers found superconducting properties of 32 metals, from self-consistent band-structure calculations and measured values of the Debye temperature, in conjunction with the theory of Gaspari and Gyorffy, λ = 0.148 (0.281 without correction) in Pd.⁵⁸ Consequently, in 2010, Doubble and co-workers concluded that the observation of strong spin fluctuations in Pd does not lead to a large contribution to the linear specific heat.⁵⁹

S & S found λ = 0.340 in Pd, α²F(ν) shown as the curve in the top part of Fig. 3, giving T_c > 0 in contrast to the experiment. They, too, are not alone. In 1978 and 1979, Pinski, Allen and Butler found λ = 0.41 in Pd,^42,60 phonons integrated over BZ following a similar practice to investigate Nb.⁴⁰ In 2005, Takezawa, Nagara and Suzuki found λ = 0.377 from DFT, phonons integrated on an 8 × 8 × 8 mesh within BZ.⁶¹ In 2006, Sklyadneva and co-workers found λ = 0.40 also from DFT, phonons integrated over a dense mesh in the first BZ.⁶²

In 2020, Kawamura, Hizume and Ozaki found superconducting properties (or the lack of them) of almost all elements in the periodic table from DFT via first-principles calculations with necessary material informatics.⁶³ They integrated phonons in BZ with the optimized tetrahedron method described briefly by Savrasov.²² Interestingly, they found λ = 0.325 or 0.333 in Pd, which leads through their first-principles μ* and spin fluctuations to T_c = 0.

In 1979, Striktzer found Pd films, irradiated with He⁺ ions, became superconductors, T_c = 3.2 K.⁶⁴ He suggested that, in the case of λ = 0.34 or so, irradiation might have smeared the Fermi surface to defeat the effect of spin fluctuations and revive superconductivity. He was open to considering smaller λ because he mentioned superconductivity enhancement by irradiation.

In 2013, Hayashi and co-workers found λ = 0.35 ± 0.05 from an angle-resolved photo-emission study of electron structure and self-energy in non-superconducting Pd.⁶⁵ It is not clear if λ∼ 0.35 still stands in ref. 42,60–63 after the authors have extended their theory to include the umklapp phonons.

XV. Discussion: alkali metals

In Fig. 2, we present a spherical Fermi surface confined entirely within the first BZ without touching its boundary. The gaps between the Fermi surface and its replicas block umklapp scattering when thermal agitations are too weak to drive the electron to jump across the gaps. In the case of alkali metals, the plot in Fig. 2 is no longer schematic but realistic. In 1971, Ekin showed with careful experiment and analysis that, in potassium, the contribution to electrical resistivity from umklapp phonons drops rapidly below ∼2.5 K.⁶⁶

Consequently, the analysis in Sections VII and VIII may not apply to alkali metals. In the absence of the umklapp phonons, represented by the curved dashed arrows in Fig. 4, we will not encounter inconsistent or double occupancies of the paired electron states. It appears reasonable for Yan, Zhang, Wang and Yang to use the MP scheme to study superconductivity in dense Li.⁶⁷ The use of the MP scheme in ref. 24 and 27 might also be justifiable if the authors could prove similarities between hydrides and alkali metals.

Furthermore, in any metal, umklapp phonons might not drive scattered electrons to cover the entire Fermi surface. There might be some pockets of the electron states on the surface beyond the reach of umklapp scattering, leaving room for normal phonons to pair electrons, evidenced by the minor disturbances of S(q) below Debye momentum in Fig. 10.

Regrettably, in either alkali or any other metals, we do not have sufficient reason to follow S & S and Andersen^21–23 to assume V = V_eff in eqn (6) and (7). It is a great convenience, but we may not be in a position to enjoy it safely. It allows us to identify, for example, the roles of acoustic or optical phonons, but we cannot be sure about it. We still have to estimate T_c roughly from the values of T_D, λ and μ* (ref. 68) not very different from the practice of McMillan in 1968.⁶⁹

XVI. Discussion: BCS interaction

The interaction term of the reduced Hamiltonian, V_kk′ in eqn (2.14) in ref. 1, has no dependence on temperature in the BCS theory. It then leads through the entropy of the electron pair ensemble to T_c. An anonymous referee asks us to replace it with its counterpart in a replacement theory⁷⁰ where we must choose the interaction term from eqn (4)–(6) if T < T_c. Otherwise, we choose eqn (7)–(9). We do not know V_kk′ before we know T_c. We also do not know T_c before we know V_kk′. We have been puzzled by the dilemma. We believe it will be more appropriate for us to follow⁷⁰ or comment on it in future.

XVII. Conclusions

Peter, Ashkenazi and Dacorogna called for a reduction of phonon contributions to superconductivity to solve a historical difficulty of the BCS theory. It has led to widespread applications of the MP scheme, readily available from computer packages, to evaluate superconductivity clearly because the scheme excludes umklapp phonons. González-Pedreros, Camargo-Martínez and Baquero believe the umklapp process does not contribute to superconductivity. Du and colleagues believe the umklapp process is relevant if the Fermi momentum is much smaller than the Debye momentum. Otherwise, they believe pairing occurs without a noticeable contribution of the umklapp phonons.

Savrasov and Savrasov investigated the effect of the MP scheme before and after a metal becomes a superconductor. In their evaluation, superconductivity in Pb is overly strong in the absence of the umklapp phonons. The electrical resistivity becomes too weak. Furthermore, Pd becomes a superconductor in ambient erroneously. It also is with weak electrical resistivity. The lesson is clear: we must reduce phonon contributions when and only when the metal is superconducting to avoid unwanted weakening of electrical resistivity.

The BCS theory indeed obliges us to do so. It assigns definite occupation probabilities to the paired electron states, sometimes inconsistent with the occupation probabilities due to umklapp scattering. We find we must exclude all the normal phonons. We also have to exclude the umklapp phonons capable of driving an electron to the other hemisphere of the Fermi surface. It leads to accurate resistivity and superconductivity in Pb and Pd. In particular, it explains why Pd is not a superconductor in the bulk form.

Further verification of the approach against other classical metallic and (or) novel high-temperature superconductors, such as various hydrides under pressure, is necessary because it concerns the worth of many previous and future endeavours. It also is urgent in light of the rapid development in the field.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank the late Professor David George Walmsley for his inspiration, experience and contributions in this direction of study.

References

1. J. Bardeen, L. N. Cooper and J. R. Schrieffer, Phys. Rev., 1957, 108, 1175–1204.
2. N. W. Ashcroft, Phys. Rev. Lett., 1968, 21, 1748–1749.
3. N. W. Ashcroft, Phys. Rev. Lett., 2004, 92, 187002.
4. C. B. Satterthwaite and I. L. Toepke, Phys. Rev. Lett., 1970, 25, 742–743.
5. B. Stritzker and W. Buckel, Z. Phys. A, 1972, 257, 1–8.
6. A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov and S. I. Shylin, Nature, 2015, 525, 73–76.
7. E. Snider, N. Dasenbrock-Gammon, R. McBride, M. Debessai, H. Vindana, K. Vencatasamy, K. V. Lawler, A. Salamat and R. P. Dias, Nature, 2020, 586, 373–377.
8. B. N. Ganguly, Z. Phys. A, 1973, 265, 433–439.
9. D. A. Papaconstantopoulos, B. M. Klein, E. N. Economou, E. N. Economou and L. L. Boyer, Phys. Rev. B: Condens. Matter Mater. Phys., 1978, 7, 141–149.
10. D. V. Semenoka, I. A. Kruglovb, I. A. Savkind, A. G. Kvashnina and A. R. Oganova, Curr. Opin. Solid State Mater. Sci., 2020, 24, 100808.
11. J. A. Flores-Livas, L. Boeri, A. Sanna, G. Profeta, R. Arita and M. Eremets, Phys. Rep., 2020, 856, 1–78.
12. F. Peng, Y. Sun, C. J. Pickard, R. J. Needs, Q. Wu and Y. Ma, Phys. Rev. Lett., 2017, 119, 107001.
13. B. Chen, L. J. Conway, W. Sun, X. Kuang, C. Lu and A. Hermann, Phys. Rev. B, 2021, 103, 035131.
14. X. Dou, X. Kuang, W. Sun, G. Jiang, C. Lu and A. Hermann, Phys. Rev. B, 2021, 104, 224510.
15. Q. Duan, J. Shen, X. Zhong, H. Lu and C. Lu, Phys. Rev. B, 2022, 105, 214503.
16. S. Li, H. Wang, W. Sun, C. Lu and F. Peng, Phys. Rev. B, 2022, 105, 224107.
17. A. Baldereschi, Phys. Rev. B: Solid State, 1973, 7, 5212–5215.
18. D. J. Chadi and M. L. Cohen, Phys. Rev. B: Solid State, 1973, 8, 5747–5753.
19. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188–5192.
20. P. Giannozzi, S. Baroni and N. Bonini, et al. , J. Phys.: Condens. Matter, 2009, 21, 395502.
21. S. Y. Savrasov, D. Y. Savrasov and O. K. Andersen, Phys. Rev. Lett., 1994, 72, 372–375.
22. S. Y. Savrasov, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 16470–16486.
23. S. Y. Savrasov and D. Y. Savrasov, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 16487–16501.
24. J. A. Camargo-Martínez, G. I. González-Pedreros and R. Baquero, Supercond. Sci. Technol., 2019, 32, 125013.
25. G. I. González-Pedreros, J. A. Camargo-Martínez and R. Baquero, Supercond. Sci. Technol., 2021, 34, 098001.
26. G. I. González-Pedreros, J. A. Camargo-Martínez and F. Mesa, Sci. Rep., 2021, 11, 22618.
27. M. Du, Z. Zhang, H. Song, H. Yu, T. Cui, V. Z. Kresinc and D. Duan, Phys. Chem. Chem. Phys., 2021, 23, 6717–6724.
28. M. Du, Z. Zhang, H. Song, H. Yu, T. Cui, V. Z. Kresinc and D. Duan, Phys. Chem. Chem. Phys., 2022, 24, 1898–1899.
29. N. F. Mott and H. Jones, The Theory of the Properties of Metals and Alloys, 1936.
30. P. G. Tomlinson and J. P. Carbotte, Phys. Rev. B: Solid State, 1976, 13, 4738–4744.
31. P. G. Tomlinson and J. P. Carbotte, Can. J. Phys., 1977, 55, 751–760.
32. A. Eiling and J. S. Schiling, J. Phys. F: Met. Phys., 1981, 11, 623–639.
33. C. R. Leavens and J. P. Carbotte, Ann. Phys., 1972, 70, 338–377.
34. P. T. Truant and J. P. Carbotte, Can. J. Phys., 1974, 52, 618–623.
35. H. K. Leung, J. P. Carbotte and C. R. Leavens, J. Low Temp., 1976, 24, 25–34.
36. M. Peter, J. Ashkenazi and M. Dacorogna, Helv. Phys. Acta, 1977, 50, 267–278.
37. B. N. Harmon and S. K. Sinha, Phys. Rev. B: Solid State, 1977, 16, 3919–3924.
38. D. Glötzel, D. Rainer and H. R. Schober, Z. Phys. B, 1979, 35, 317–326.
39. A. Al-Lehaibi, J. C. Swihart, W. H. Butler and F. J. Pinski, Phys. Rev. B: Condens. Matter Mater. Phys., 1987, 36, 4103–4111.
40. W. H. Butler, H. G. Smith and N. Wakabayashi, Phys. Rev. Lett., 1977, 39, 1004–1007.
41. W. H. Butler, F. J. Pinski and P. B. Allen, Phys. Rev. B: Condens. Matter Mater. Phys., 1979, 19, 3708–3721.
42. F. J. Pinski, P. B. Allen and W. H. Butler, Phys. Rev. Lett., 1978, 41, 431–434.
43. X. H. Zheng and D. G. Walmsley, Phys. Scr., 2014, 89, 095803.
44. X. H. Zheng and D. G. Walmsley, Phys. Scr., 2019, 94, 085801.
45. J. P. Carbotte and R. C. Dynes, Phys. Rev., 1968, 172, 476–484.
46. P. J. Hasnip, K. Refson, M. I. J. Probert, J. R. Yates, S. J. Clark and C. J. Pickard, Phil. Trans. R. Soc. A, 2022, 372, 20130270.
47. J. C. Slater, Phys. Rev., 1937, 51, 846–851.
48. N. W. Ashcroft, Phys. Lett., 1966, 23, 48–50.
49. W. L. McMillan and J. M. Rowell, Phys. Rev. Lett., 1965, 14, 108–112.
50. X. H. Zheng and J. X. Zheng, Supercond. Sci. Technol., 2021, 34, 098002.
51. X. H. Zheng and J. X. Zheng, Phys. Chem. Chem. Phys., 2022, 24, 1896–1897.
52. J. P. Carbotte, Rev. Mod. Phys., 1990, 62, 1027–1157.
53. R. Stedman, L. Almqvist and G. Nilson, Phys. Rev., 1967, 162, 549–557.
54. A. P. Miiller and B. N. Brockhouse, Can. J. Phys., 1971, 49, 704–723.
55. X. H. Zheng and D. G. Walmsley, J. Low Temp. Phys., 2012, 166, 279–297.
56. R. Hooke and T. A. Jeeves, J. Assoc. Comput. Mach., 1961, 8, 212–229.
57. N. W. Ashcroft and N. D. Mermin, Solid State Physics, Brooks/Cole, 1976.
58. D. A. Papaconstantopoulos, B. M. Klein, A. R. Williams, V. L. Morruzzi and J. F. Janak, Phys. Rev. B: Solid State, 1977, 15, 4221–4226.
59. R. Doubble, S. M. Hayden, P. Dai, H. A. Mook, J. R. Thompson and C. D. Fros, Phys. Rev. Lett., 2010, 105, 027207.
60. F. J. Pinski and W. H. Butler, Phys. Rev. B: Condens. Matter Mater. Phys., 1979, 19, 6010–6015.
61. T. Takezawa, H. Nagara and N. Suzuki, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 012515.
62. I. Yu Sklyadneva, A. Leonardo, P. M. Echenique, S. V. Eremeev and E. V. Chulkov, J. Phys.: Condens. Matter, 2006, 18, 7923–7935.
63. M. Kawamura, Y. Hizume and T. Ozaki, Phys. Rev. B, 2020, 101, 134511.
64. B. Striktzker, Phys. Rev. Lett., 1979, 42, 1769–1793.
65. H. Hayashi, K. Shimada, J. Jiang, H. Iwasawa, Y. Aiura, T. Oguchi, H. Namatame and M. Taniguchi1, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 035140.
66. J. W. Ekin, Phys. Rev. Lett., 1971, 26, 1551–1553.
67. Y. Yan, Y. Zhang, Y. Wang and G. Yang, Solid State Commun., 2016, 225, 7–11.
68. X. H. Zheng, J. X. Zheng and D. G. Walmsley, Solid State Commun., 2019, 297, 17–20.
69. W. L. McMillan, Phys. Rev., 1968, 167, 331–344.
70. A. D. Arulsamy, Phys. Scr., 2029, 94, 055803.

---