The unbalanced phonon-induced superconducting state on a square lattice beyond the static boundary

TL;DR

This study investigates the induction of superconductivity on a square lattice via electron-phonon interactions, revealing a lower unbalance parameter threshold than the static boundary and significant deviations from BCS theory predictions.

Contribution

The paper demonstrates superconductivity induction at unbalance parameters below the static boundary using Eliashberg equations with momentum-dependent interactions, extending understanding beyond BCS theory.

Findings

Superconductivity occurs for unbalance parameter $\, extless 0.42$

Thermodynamic functions differ markedly from BCS predictions

The static boundary $\, extgreater 0.93$ is not reached in the dynamic case

Abstract

The paper presents our verification of induction of the superconducting state on a square lattice by the linear electron-phonon interaction for values of the unbalance parameter ($\gamma=\lambda_{D}/\lambda_{ND}$) less than $\gamma_{C}=0.42$. Symbols $\lambda_{D}$ and $\lambda_{ND}$ denote the values of the coupling constant in the diagonal and the non-diagonal channel of the self-energy. Calculations were carried out using the Eliashberg equations, in which the order parameter ($\Delta_{\bf k}\left(i\omega_{n}\right)$) and the wave function renormalising factor ($Z_{\bf k}\left(i\omega_{n}\right)$) depend explicitly on the Matsubara frequency ($\omega_{n}$) and the wave vector (${\bf k}$). The value of $\gamma_{C}$ in the static boundary ($\Delta_{\bf k}\left(i\omega_{n}\right)\rightarrow \Delta_{\bf k}\left(i\omega_{n=1}\right)$), equal to ($0.93$), is significantly greater than the obtained limit value of $0.42$. Values of the thermodynamic functions of the superconducting state determined for our assumptions are significantly different from the values calculated in accordance with the BCS theory. The results were obtained for the electron-phonon interaction function explicitly dependent on the momentum transfer between electron states.

Tables

Table 1. Table 1: Parameters λ 𝜆 \lambda , ω ln subscript 𝜔 ln \omega_{{\rm ln}} , and ω 2 subscript 𝜔 2 \omega_{2} represent the electron–phonon coupling constant, the logarithmic phonon frequency, and the second moment of the normalized weight function, respectively. Functions f 1 subscript 𝑓 1 f_{1} and f 2 subscript 𝑓 2 f_{2} are the strong–coupling correction function and the shape correction function, respectively.

Quantity	Value
$\lambda =2{\int }_0^{+\infty }𝑑\Omega \frac{{\alpha }^2\left(\Omega \right)F\left(\Omega \right)}{\Omega }$	$2.06$
${\omega }_{\ln }=\exp \left[\frac{2}{\lambda }{\int }_0^{+\infty }𝑑\Omega \frac{{\alpha }^{2}F\left(\Omega \right)}{\Omega }\ln \left(\Omega \right)\right]$	$0.109$ $t$
$\sqrt{{\omega }_2}={\left[\frac{2}{\lambda }{\int }_0^{+\infty }𝑑\Omega {\alpha }^{2}F\left(\Omega \right)\Omega \right]}^{1/2}$	$0.208$ $t$
$f_1={\left[1+{\left(\frac{\lambda }{{\Lambda }_1}\right)}^{\frac{3}{2}}\right]}^{\frac{1}{3}}$	-
$f_2=1+\frac{\left(\frac{\sqrt{{\omega }_2}}{{\omega }_{\ln }}-1\right){\lambda }^2}{{\lambda }^2+{\Lambda }_2^2}$	-
${\Lambda }_1=2.46$	-
${\Lambda }_2=1.82\left(\sqrt{{\omega }_2}/{\omega }_{\ln }\right)$	-

Equations

$${\phi }_{\mathbf{k}}\left(i{\omega }_n\right)=\frac{1}{\beta N}\sum _{m\mathbf{q}}{K}_{\mathbf{q}}\left({\omega }_n-{\omega }_m\right)\frac{{\phi }_{\mathbf{k}-\mathbf{q}}\left(i{\omega }_m\right)}{D_{\mathbf{k}-\mathbf{q}}\left(i{\omega }_m\right)},$$

$$Z_{\mathbf{k}}\left(i{\omega }_n\right)=1+\frac{\gamma }{\beta N}\sum _{m\mathbf{q}}\frac{{\omega }_m}{{\omega }_n}{K}_{\mathbf{q}}\left({\omega }_n-{\omega }_m\right)\frac{Z_{\mathbf{k}-\mathbf{q}}\left(i{\omega }_m\right)}{D_{\mathbf{k}-\mathbf{q}}\left(i{\omega }_m\right)},$$

$$Z_{m_0}\left({\mathbf{q}}_0\right)=1+\sum _{L=1}^{+\infty }\frac{1}{{\left(\beta N\right)}^L}\sum _{{\mathbf{q}}_1\sim {\mathbf{q}}_L=0^{+}}^{\pi }\sum _{m_1\sim m_L=1}^M\frac{{\omega }_{m_L}}{{\omega }_{m_0}}\prod _{j=1}^L{J}_{m_{j-1}{m}_j}\left({\mathbf{q}}_{j-1},{\mathbf{q}}_j\right)d_{m_j}^{-1}\left({\mathbf{q}}_j\right),$$

$$\rho \left(\epsilon \right)=\frac{1}{N}\sum _{k_x,k_y}\delta \left(\epsilon -{\epsilon }_{k_x,k_y}\right)={\rho }_0{\int }_{-\pi }^{\pi }d{k}_x{\int }_{-\pi }^{\pi }d{k}_y\delta \left(\epsilon -{\epsilon }_{k_x,k_y}\right).$$

$$F\left(\omega \right)=\frac{1}{N}\sum _{q_x,q_y}\delta \left(\omega -{\omega }_{q_x,q_y}\right)=F_0{\int }_{-\pi }^{\pi }d{q}_x{\int }_{-\pi }^{\pi }d{q}_y\delta \left(\omega -{\omega }_{q_x,q_y}\right).$$

$${\alpha }^{2}F\left(\mathbf{k},{\mathbf{k}}^{\prime },\omega \right)=\rho \left(0\right)\sum _{\text{mbox}\mathbf{q}}|g\left(\mathbf{q},\mathbf{k},{\mathbf{k}}^{\prime }\right){|}^2\delta \left(\omega -{\omega }_{\text{mbox}\mathbf{q}}\right).$$

$${\alpha }^{2}F\left(\omega \right)=\frac{1}{N^2}\sum _{\text{mbox}\mathbf{k},\text{mbox}{\mathbf{k}}^{\prime }}{w}_{\text{mbox}\mathbf{k}}{w}_{\text{mbox}{\mathbf{k}}^{\prime }}{\alpha }^{2}F\left(\mathbf{k},{\mathbf{k}}^{\prime },\omega \right),$$

$${\alpha }^{2}F\left(\omega \right)=\rho \left(0\right)\sum _{\text{mbox}\mathbf{q}}{g}_{\text{mbox}\mathbf{q}}^2\delta \left(\omega -{\omega }_{\text{mbox}\mathbf{q}}\right).$$

$$k_B{T}_C=f_1{f}_2\frac{{\omega }_{\ln }}{1.2}\exp \left[\frac{-1.04\left(1+\lambda \right)}{\lambda }\right],$$

$$R_{\Delta }$$

$$R_C$$

$$R_H$$

$$R_{\Delta }$$

$$R_C$$

$$R_H$$

$${\phi }_n$$

$$Z_n$$

Full text

Physical properties of the phonon–induced superconducting state on a square lattice were determined in the present work within the Eliashberg formalism. The initial results, obtained in the limit of momentum approximation (which neglects the self–consistency with respect to Matsubara frequency), were already presented by us elsewhere Szewczyk et al. (2018). There we showed that the balanced superconducting state on a square lattice cannot arise in the case of the constant value of the electron–phonon coupling function ($g_{\mbox{{\scriptsize{\bf q}}}}\sim g$). By the balanced superconducting state we mean such a state, for which the coupling constant $\lambda_{D}$ in the diagonal channel of the self–energy matrix ($M_{\mbox{{\scriptsize{\bf k}}}}\left(i\omega_{n}\right)$) is equal to the coupling constant in its non–diagonal channel $\lambda_{ND}$. Please notice that the diagonal elements of the self–energy determine, inter alia, the influence of the electron–phonon interaction on the value of the effective electron mass. The off–diagonal elements of the $M_{\mbox{{\scriptsize{\bf k}}}}\left(i\omega_{n}\right)$ are directly related to the order parameter of the phase transition between the superconducting and the normal state.

At present, we know many physical mechanisms which can lead to the unbalanced phonon–induced superconducting state. Particular attention should be paid to the spin–fluctuation (paramagnon) scattering, which contributes to the effective coupling constant in the diagonal channel in the opposite way than it does in the antidiagonal channel ($s$–wave symmetry) Rietschel and Winter (1979); Rietschel et al. (1980); Fay and Appel (1980); Daams et al. (1981); Fay and Appel (1982); Leavens and MacDonald (1983). The unbalance of the superconducting state can also be generated by the non–conventional terms of the electron–phonon interaction, e.g. the electron–electron–phonon Szczȩśniak (2012) or, similarly, the hole–hole–phonon interaction Durajski and Szczȩśniak (2018). Additionally, they induce the asymmetry of the electron density of states, the pseudo–gap in the electron density of states, or the anomalous dependence of the energy gap on doping Szczȩśniak et al. (2017). These effects can be of particular importance for the correct description of the superconducting state in cuprates Bednorz and Müller (1986, 1988); Dagotto (1994), especially since they are observed experimentally Renner et al. (1998a, b); Fischer et al. (2007). At the level of the conventional Eliashberg formalism, the unbalanced superconducting state arises also when the non–phonon contribution from the $d$–wave type symmetry is taken into account in the $\alpha^{2}_{\mbox{{\scriptsize{\bf k}}}}F\left(\omega\right)$ Eliashberg function modeling the electron–phonon interaction Rieck et al. (1990); Jiang et al. (1993); Zasadzinski et al. (2003). Nevertheless the considered model is too simple for the satisfactory explanation of even some of the anomalous properties of the superconducting state in cuprates.

The purpose of our investigation, the results of which are gathered in the present work, was to find information, as precise as possible, on the physical properties of the phonon–induced superconducting state on a square lattice. Therefore we solved the Eliasberg equations Eliashberg (1960) with respect both to the electron wave vector and to the Matsubara frequency in the fully self–consistent way. We took into consideration the electron–phonon coupling function explicitly dependent on the phonon transfer of momentum between electron states Allan et al. (2017). It is worth noticing that the Eliashberg equations on a square lattice are for the first time analysed in such a precise manner.

The obtained results allowed us also to comment on the problem of the significance of the electron–phonon interaction in the context of determination of the pairing mechanism in cuprates Dagotto (1994); Szczȩśniak (2012). This issue is discussed at greater length in the concluding part of the work.

In its most advanced form, the analysis of the phonon–induced superconducting state is performed within the Eliashberg formalism, which is based on the self–consistent calculation of the matrix elements of the thermodynamic Green’s function: $G_{\mbox{{\scriptsize{\bf k}}}}\left(i\omega_{n}\right)=\left<\left<\Psi_{\mbox{{\scriptsize{\bf k}}}}|\Psi^{\dagger}_{\mbox{{\scriptsize{\bf k}}}}\right>\right>_{i\omega_{n}}$. The symbol $\Psi^{\dagger}_{\mbox{{\scriptsize{\bf k}}}}$ denotes the Nambu spinor with Hermitian coupling: $\Psi^{\dagger}_{\mbox{{\scriptsize{\bf k}}}}=\left(c^{\dagger}_{\mbox{{\scriptsize{\bf k}}}\uparrow},c_{-\mbox{{\scriptsize{\bf k}}}\downarrow}\right)$, and $\omega_{n}=\frac{\pi}{\beta}\left(2n-1\right)$ stands for the fermionic Matsubara frequency ($\beta=1/k_{B}T$). The quantity $c^{\dagger}_{\mbox{{\scriptsize{\bf k}}}\sigma}$ is the creation operator of the electron state with momentum ${\bf k}$ and spin $\sigma\in\{\uparrow,\downarrow\}$. The explicit form of the self-energy $M_{\mbox{{\scriptsize{\bf k}}}}\left(i\omega_{n}\right)$ is determined as a rule by application of the Hamiltonian which models the linear coupling between the electron and the phonon subsystems Fröhlich (1952): $H=\sum_{\mbox{{\scriptsize{\bf k}}}}\varepsilon_{\mbox{{\scriptsize{\bf k}}}}\Psi^{\dagger}_{\mbox{{\scriptsize{\bf k}}}}\tau_{3}\Psi_{\mbox{{\scriptsize{\bf k}}}}+\sum_{\mbox{{\scriptsize{\bf q}}}}\omega_{\mbox{{\scriptsize{\bf q}}}}b^{\dagger}_{\mbox{{\scriptsize{\bf q}}}}b_{\mbox{{\scriptsize{\bf q}}}}+\frac{1}{\sqrt{N}}\sum_{\mbox{{\scriptsize{\bf k}}}\mbox{{\scriptsize{\bf q}}}}g_{\mbox{{\scriptsize{\bf q}}}}\Psi^{\dagger}_{\mbox{{\scriptsize{\bf k}}}+\mbox{{\scriptsize{\bf q}}}}\tau_{3}\Psi_{\mbox{{\scriptsize{\bf k}}}}\phi_{\mbox{{\scriptsize{\bf q}}}}$, where $\varepsilon_{\mbox{{\scriptsize{\bf k}}}}$ is the electronic dispersion relation and $\phi_{\mbox{{\scriptsize{\bf q}}}}=b_{\mbox{{\scriptsize{\bf q}}}}+b^{\dagger}_{-\mbox{{\scriptsize{\bf q}}}}$. The quantity $b^{\dagger}_{\mbox{{\scriptsize{\bf q}}}}$ denotes the creation operator of the phonon state with momentum ${\bf q}$. The symbol $\omega_{\mbox{{\scriptsize{\bf q}}}}$ represents phonon dispersion relation ($\omega_{\mbox{{\scriptsize{\bf q}}}}=\omega_{-\mbox{{\scriptsize{\bf q}}}}$). The relatively easy quantum operator calculations bring about the exact result: $M_{\mbox{{\scriptsize{\bf k}}}}\left(i\omega_{n}\right)=\frac{1}{N}\sum_{\mbox{{\scriptsize{\bf q}}}\mbox{{\scriptsize{\bf q}}}^{\prime}}g_{\mbox{{\scriptsize{\bf q}}}}g_{\mbox{{\scriptsize{\bf q}}}^{\prime}}\tau_{3}\left<\left<\Psi_{\mbox{{\scriptsize{\bf k}}}-\mbox{{\scriptsize{\bf q}}}}\phi_{\mbox{{\scriptsize{\bf q}}}}|\Psi^{\dagger}_{\mbox{{\scriptsize{\bf k}}}+\mbox{{\scriptsize{\bf q}}}^{\prime}}\phi_{\mbox{{\scriptsize{\bf q}}}^{\prime}}\right>\right>_{i\omega_{n}}\tau_{3}$. Further analysis of the problem is based on the Migdal approximation and on the Wick’s theorem Migdal (1958); Fetter and Walecka (1971), by virtue of which the thermodynamic Green’s function $\left<\left<\Psi_{\mbox{{\scriptsize{\bf k}}}-\mbox{{\scriptsize{\bf q}}}}\phi_{\mbox{{\scriptsize{\bf q}}}}|\Psi^{\dagger}_{\mbox{{\scriptsize{\bf k}}}+\mbox{{\scriptsize{\bf q}}}^{\prime}}\phi_{\mbox{{\scriptsize{\bf q}}}^{\prime}}\right>\right>_{i\omega_{n}}$ is approximated by the product of the $G_{\mbox{{\scriptsize{\bf k}}}}\left(i\omega_{n}\right)$ function and the phonon propagator for the non–interacting phonons. The procedure of reaching self–consistency allows to obtain the non–linear Eliashberg equations, which determine the thermodynamics of the superconducting phase on the quantitative level. For the half-filled electon band, at the occurrence of the slight particle–hole symmetry violation, the Eliashberg equations take the form:

[TABLE]

[TABLE]

The symbol $\varphi_{\bf k}\left(i\omega_{n}\right)$ denotes the order parameter function, the order parameter being defined by the $\Delta_{\bf k}\left(i\omega_{n}\right)=\varphi_{\bf k}\left(i\omega_{n}\right)/Z_{\bf k}\left(i\omega_{n}\right)$ ratio, where $Z_{\bf k}\left(i\omega_{n}\right)$ represents the wave function renormalisation factor. The pairing kernel of the electron–phonon interaction is given by the formula: $K_{\bf q}\left(\omega_{n}-\omega_{m}\right)=2g^{2}_{\bf q}\frac{\omega_{\mbox{{\scriptsize{\bf q}}}}}{\left(\omega_{n}-\omega_{m}\right)^{2}+\omega^{2}_{\mbox{{\scriptsize{\bf q}}}}}$. Additionally, $D_{\bf k}\left(i\omega_{n}\right)=\left(\omega_{n}Z_{\bf k}\left(i\omega_{n}\right)\right)^{2}+\varepsilon^{2}_{{\bf k}}+\varphi^{2}_{\bf k}\left(i\omega_{n}\right)$. Usually the Eliashberg equations are solved in the isotropic approximation: $\varphi_{\bf k}\left(i\omega_{n}\right)\sim\varphi\left(i\omega_{n}\right)$ and $Z_{\bf k}\left(i\omega_{n}\right)\sim Z\left(i\omega_{n}\right)$, while the pairing kernel is transformed in the following way: $K_{\bf q}\left(\omega_{n}-\omega_{m}\right)\sim 2\int^{\omega_{D}}_{0}d\omega^{\prime}\frac{\alpha^{2}F\left(\omega^{\prime}\right)}{\left(\omega_{n}-\omega_{m}\right)^{2}+\omega^{\prime 2}}$, where $\omega_{D}$ is the Debye phonon frequency. The spectral function $\alpha^{2}F\left(\omega\right)$ for a specific physical system can be determined by means of the DFT method. This approach works perfectly e.g. for the recently discovered electron–phonon hydrogen–containing superconductors Li et al. (2014); Duan et al. (2014); Durajski et al. (2015, 2016); Durajski and Szczȩśniak (2017); Szczȩśniak and Durajski (2018); Kostrzewa et al. (2018); Kruglov et al. (2018), characterised by the extremely high values of the critical temperature: ${\rm H_{2}S}$ ($T_{C}=150$ K for $p=150$ GPa), ${\rm H_{3}S}$ ($T_{C}=203$ K for $p=150$ GPa Drozdov et al. (2014, 2015), and ${\rm LaH_{10}}$ ($T_{C}=215$ K for $p=150$ GPa) Drozdov et al. (2019), ($T_{C}=260$ K for $p\in\left(180-200\right)$ GPa) Somayazulu et al. (2019).

Taking into account the results reported in Szewczyk et al. (2018), we supplemented the Eq. (2) with the $\gamma$ parameter, which determines the degree of unbalance of the electron–phonon coupling constants in the diagonal and the non–diagonal channel of self–energy.

We performed numerical calculations on the $N\times N$ momentum lattice, where $N=200$, and we took into account $200$ Matsubara frequencies. Therefore eight million equations had to be solved. To minimize the hardware requirements, the Eliashberg equations were written in the form which explicitly expressed all symmetries of solutions with respect to the wave vector and to the Matsubara frequency. Assuming that the auxiliary symbol $f_{\mbox{{\scriptsize{\bf k}}}}\left(i\omega_{n}\right)=f_{n}\left(k_{x},k_{y}\right)$ should be understood as $\varphi$ or $Z$, one can write that: $f_{n}\left(k_{x},k_{y}\right)=f_{n}\left(-k_{x},k_{y}\right)=f_{n}\left(k_{x},-k_{y}\right)=f_{n}\left(-k_{x},-k_{y}\right)$, and $f_{n}\left(k_{x},k_{y}\right)=f_{-n+1}\left(k_{x},k_{y}\right)$. The Eliashberg equations were solved by means of the AMCIA program, which generalises the numerical procedures used in another our work Szczȩśniak (2006a). The AMCIA numerical environment requires introduction of the initial values of the $\varphi_{\bf k}(i\omega_{n})$ and $Z_{\bf k}(i\omega_{n})$ functions. As far as the wave function renormalising factor is concerned, we took into account the first $25$ terms of the series:

[TABLE]

where: ${\bf q}_{j}=\left(q^{j}_{x},q^{j}_{y}\right)$ and $J_{m_{j-1}m_{j}}\left({\bf q}_{j-1},{\bf q}_{j}\right)=J_{m_{j-1}m_{j}}\left(q^{j-1}_{x},q^{j-1}_{y},q^{j}_{x},q^{j}_{y}\right)=K^{\prime}_{m_{j-1}m_{j}}\left(q^{j-1}_{x}+q^{j}_{x},q^{j-1}_{y}+q^{j}_{y}\right)+K^{\prime}_{m_{j-1}m_{j}}\left(q^{j-1}_{x}+q^{j}_{x},q^{j-1}_{y}-q^{j}_{y}\right)+K^{\prime}_{m_{j-1}m_{j}}\left(q^{j-1}_{x}-q^{j}_{x},q^{j-1}_{y}+q^{j}_{y}\right)+K^{\prime}_{m_{j-1}m_{j}}\left(q^{j-1}_{x}-q^{j}_{x},q^{j-1}_{y}-q^{j}_{y}\right)$, while $K^{\prime}_{m_{j-1}m_{j}}\left(q_{x},q_{y}\right)=K_{\bf q}\left(\omega_{m_{j-1}}-\omega_{m_{j}}\right)-K_{\bf q}\left(\omega_{m_{j-1}}+\omega_{m_{j}}\right)$. Additionally: $d_{m_{j}}\left({\bf q}_{j}\right)=\omega^{2}_{m_{j}}+\varepsilon^{2}_{{\bf q}_{j}}$. We choose positive non-zero numbers as initial values for the order parameter function. It should be mentioned that the number of the momentum lattice points and the number of Matsubara frequencies taken into account in our calculations provide for the convergence of the obtained results (see Appendix A).

The energy of electron states was modeled by using the tight–binding description on the ion equilibrium positions assuming the nearest–neighbor ($t$) and next–nearest neighbor ($t^{\prime}$) hopping Dagotto (1994): $\varepsilon_{\mbox{{\scriptsize{\bf k}}}}=-2t\left[\cos\left(k_{x}\right)+\cos\left(k_{y}\right)\right]+4t^{\prime}\cos\left(k_{x}\right)\cos\left(k_{y}\right)$, where the lattice constant is taken as unity. In the present work we assumed that energy is expressed in units of hopping integral $t$, and $t^{\prime}=0.1$ $t$. From the physical point of view, the low value of the $t^{\prime}$ integral indicates the relatively small violation of the particle–hole symmetry. Therefore it can be assumed for the considered case that the values of the chemical potential and the energy shift functions are equal to zero. For higher values of $t^{\prime}$ the full system of Eliashberg equations should be analysed according to the scheme presented in References Carbotte (1990); Radtke and Norman (1994). We took into account the acoustic phonons stemming from the nearest–neighbor spring: $\omega_{\mbox{{\scriptsize{\bf q}}}}=\omega_{0}\sqrt{2-\cos{q_{x}}-\cos{q_{y}}}$. In model calculations it was assumed that $\omega_{0}=0.15$ $t$, so that the maximum phonon frequency $\omega_{D}=2\omega_{0}=0.3$ $t$. In the case of a square lattice, the electron–phonon matrix elements are given by Allan et al. (2017): $g_{\mbox{{\scriptsize{\bf q}}}}=g_{0}|{\bf q}|\sqrt{1/\omega_{\mbox{{\scriptsize{\bf q}}}}}$, and $g_{0}=0.031$ $t^{3/2}$. The value of the $g_{\mbox{{\scriptsize{\bf q}}}}$ function does not exceed $0.3$ $t$ in the considered case. Let us notice that the function $g_{\mbox{{\scriptsize{\bf q}}}}$ reproduces the results reported in Madelung (1978); Coleman (2015); Giustino (2017), and has the structure predicted by Bloch Bloch (1928). It should be noticed, however, that in a general case the electron–phonon coupling function depends explicitly also on the ${\bf k}$ vector: $g_{\mbox{{\scriptsize{\bf k}}},\mbox{{\scriptsize{\bf k}}}+\mbox{{\scriptsize{\bf q}}}}$.

As the first step, let us discuss the properties of the superconducting state in the isotropic approximation and find the thermodynamic properties of the superconducting phase by means of analitical formulae. The obtained results shall also make us sure that the initial values of parameters ($\omega_{0}$ and $g_{0}$) were chosen in the physically acceptable way.

The electron density of states function is described by the formula:

[TABLE]

The symbol $\rho_{0}$ denotes the normalisation constant, which can be estimated from the condition: $\int^{W_{u}}_{W_{d}}d\varepsilon\rho\left(\varepsilon\right)=1$, where $W_{d}$ and $W_{u}$ are the upper and the lower limits of the electron band, equal to $W_{u}=4.6$ $t$ and $W_{d}=-4.1$ $t$, respectively. Additionally, $\rho_{0}=0.0254144$. The profile of the electron density of states function is shown in Fig. 1 (a). The value of the considered function at the Fermi level is relatively high ($\rho\left(0\right)=0.18724$ $t^{-1}$). This results from the presence of the van Hove singularity in the electron density of states van Hove (1953); Goicochea (1994); Szczȩśniak and Grabiński (2002). It is relatively close to the Fermi energy level due to the fact that the hopping integral $t$ is significantly greater than $t^{\prime}$ (the slight violation of the particle–hole interaction). The phonon density of states is calculated in the similar manner:

[TABLE]

The $F_{0}$ value was determined from the normalisation condition: $\int^{\omega_{D}}_{0}d\omega F\left(\omega\right)=1$ and found to be $F_{0}=0.026079$. The phonon density of states function is plotted in Fig. 1 (b). The maximum value of the function $F\left(\omega\right)$ corresponding to the energy $\omega\sim 0.21$ $t$ is associated with the flattening of the phonon dispersion relation. For $\omega\rightarrow 0$ we obtained $F\left(\omega\right)\rightarrow 0$. As expected, the phonon density of states function vanishes for the frequency $\omega_{D}$. In the case of the Eliashberg function, we started from the general formula:

[TABLE]

Then the spectral function $\alpha^{2}F\left({\bf k},{\bf k}^{\prime},\omega\right)$ was averaged over the Fermi surface in order to achieve the isotropic Eliashberg function:

[TABLE]

where: $w_{\mbox{{\scriptsize{\bf k}}}}=\delta\left(\varepsilon_{\mbox{{\scriptsize{\bf k}}}}\right)/\rho\left(0\right)$. Simple calculations led to the result:

[TABLE]

We presented the isotropic Eliashberg function in Fig. 1 (c). While comparing its profile with the profile of the phonon density of states function, one can notice a difference in their shapes. The difference characterises the contribution from the squares of the electron–phonon coupling and the electron density of states functions.

We assessed the walue of the critical temperature from the Allen–Dynes formula Allen and Dynes (1975):

[TABLE]

however we did not take into account the depairing electron correlations modeled by Coulomb pseudopotential $\mu^{\star}$ Morel and Anderson (1962) in our analysis. This correlations, if included, shift the critical temperature value down in the considered model, but they usually do not extinguish totally the superconducting phase, because the pseudopotential does not take very high values ($\mu^{\star}\sim 0.1$). The symbols used in Eq. (9) are defined in Tab. 1.

The dimensionless coefficients corresponding to those used in the BCS theory can be also calculated in a relatively easy way Bardeen et al. (1957a, b):

[TABLE]

Values of these coefficients calculated within the BCS theory are equal to $3.53$, $1.43$, and $0.168$, respectively. In the considered case, however, when we can use the explicit form of the Eliashberg function, the coefficients $R_{\Delta}$, $R_{C}$ and $R_{H}$ should be assessed on the basis of the formulae Carbotte (1990):

[TABLE]

For the values of initial parameters assumed in the work we obtained: $k_{B}T_{C}=0.029$ $t$. Additionally: $R_{\Delta}=5.5$, $R_{C}=2.64$ oraz $R_{H}=0.135$. It means that the physical system is characterised by a strong electron–phonon coupling, what can be also concluded from the value of $\lambda$ (see Tab. 1). If $t=250$ meV, then the critical temperature is equal to $84.19$ K.

Let us proceed to the analysis of the thermodynamic properties of the superconducting state using the full Eliashberg equations. It occurs that the achieved results are significantly different from the results obtained within the isotropic approximation.

The dependence between the averaged value of the order parameter $\left<\Delta(i\omega_{n=1})\right>=\frac{1}{N}\sum_{\mbox{{\scriptsize{\bf k}}}}\Delta_{\mbox{{\scriptsize{\bf k}}}}(i\omega_{n=1})$ and the value of the unbalance parameter is plotted in Fig. 2 (a). It can be easily noticed that the superconducting phase ceases to exist for $\gamma\geq\gamma_{C}=0.42$. Let us remember that within the momentum approximation with a constant value of the electron–phonon coupling function ($g_{\mbox{{\scriptsize{\bf q}}}}\sim g$) we got $\gamma_{C}=0.94$ Szewczyk et al. (2018). A similar value, $\gamma_{C}=0.93$, is obtained from calculations in which the momentum transfer is explicitly taken into account in the electron–phonon coupling function (see Fig. 2 (b)). It should be strongly emphasized that our result does not mean that the balanced phonon–induced superconducting state cannot exists on a square lattice in general. Attention should be paid to the fact that, despite the relatively sophisticated analysis presented in the work, we applied also a series of approximations. They can significantly influence the final results. We would like to mention here, in particular, the lack of explicit dependence of the electron–phonon coupling function on the $\bf k$ vector, omitting the vertex corrections to the electron–phonon interaction in the analysis (and yet $\omega_{D}/W\sim 0.1$, where $W\sim 4$ $t$), or disregarding the anharmonic corrections, etc. Therefore our result should be considered as physically interesting, but acquired under some specific assumptions. Formulation of general conclusions with respect to the phonon–induced superconducting state on a square lattice demands for still deeper analysis in future.

The dependence of the electron effective mass ($m^{\star}_{e}$) on the value of the unbalance parameter $\gamma$ is presented in Fig. 2 (c), wherein $m^{\star}_{e}=\left<Z(i\omega_{n=1})\right>m_{e}$, and the symbol $m_{e}$ denotes the band mass of an electron. The obtained results prove that the superconducting phase on a square lattice ceases due to the fact that the electron effective mass takes high values ($\left[m^{\star}_{e}\right]_{\gamma=\gamma_{C}}=5.18$ $m_{e}$ and $\left[m^{\star}_{e}\right]_{\gamma=1}=9.48$ $m_{e}$), what results in the decrease in the critical temperature: $T_{C}\sim\exp\left(-m^{\star}_{e}\right)$ McMillan (1968); Allen and Dynes (1975). The values of the electron effective mass for the momentum approximation and the explicitly given $g_{\bf q}$ are significantly less (Fig. 2 (d)). Quite similar values of the electron effective mass we can observe for the momentum approximation ($\left[m^{\star}_{e}\right]_{\gamma=\gamma_{C}}\sim 4.73$ $m_{e}$ and $\left[m^{\star}_{e}\right]_{\gamma=1}=5.46$ $m_{e}$) under the additional assumption $g_{\mbox{{\scriptsize{\bf q}}}}\sim g$ Szewczyk et al. (2018). Let us notice that the above mentioned results agree with the fundamental result achieved by Cooper, who proved that the Fermi system is unstable towards the formation of bound Cooper pairs for any infinitesimally small coupling Cooper (1956). Cooper’s considerations dealt with the case of $\gamma=0$ ($m^{\star}_{e}=m_{e}$), for which we also observed an induction of the superconducting condensate.

As far as the results achieved within the isotropic approximation ($\Delta_{\bf k}\left(i\omega_{n}\right)\sim\Delta\left(i\omega_{n}\right)$ and $Z_{\bf k}\left(i\omega_{n}\right)\sim Z\left(i\omega_{n}\right)$) are concerned, the unbalanced Eliashberg equations were solved by Cappelluti and Ummarino in Cappelluti and Ummarino (2007). For the case of 3D system, the existence of the balanced superconducting state was confirmed, and even the unbalanced superconducting state of the $\gamma$ parameter value significantly exceeding unity was stated. The calculations which we performed for the isotropic approximation on a square lattice also predict induction of either the balanced or the unbalanced superconducting state for $\gamma>1$ (see Appendix A). The above results mean that the isotropic approximation, at least with respect to the square lattice, should be applied with considerable caution.

Having acquired the explicit solutions of the Eliashberg equations, one can determine the regions of the first Brillouin zone, for which the function $m^{\star}_{e}\left({\bf k}\right)$ contributes particularly greatly to the averaged value of the effective electron mass. We present diagrams showing the values of the wave function renormalisation factor versus the wave vector in Fig. 3 (left column). It can be seen that an increase in the effective electron mass is generated mainly by increase in the value of $m^{\star}_{e}\left({\bf k}\right)$ function in the whole Brillouin zone. Although at the boundaries of the Brillouin zone it is particularly high. The order parameter function $\varphi_{\bf k}\left(i\omega_{n=1}\right)$, which is plotted in Fig. 3 (right column), has a similar form.

Our results prove that only the unbalanced superconductiong state of $\gamma<\gamma_{C}=0.42$ can be induced on a square lattice. We selected four values of the unbalance parameter ([math], $0.12\gamma_{C}$, $0.24\gamma_{C}$, and $0.31\gamma_{C}$) in order to analyse thermodynamical properties of this state. The achieved results are plotted in Fig. 4 for both the order parameter and the wave function renormalising factor. It can be seen that in every case the dependence of the order parameter on temperature deviates considerably from the one predicted by the BCS theory Bardeen et al. (1957a, b). Taking into account the ratio $R_{\Delta}=2\Delta\left(0\right)/k_{B}T_{C}$ we get: $R^{0}_{\Delta}=3.37$, $R^{0.12\gamma_{C}}_{\Delta}=3.13$, $R^{0.24\gamma_{C}}_{\Delta}=3$, and $R^{0.31\gamma_{C}}_{\Delta}=2.92$. Please notice that the BCS theory predicts $R_{\Delta}=3.53$.

The influence of temperature on the value of the effective electron mass is presented in Fig. 4 (b). The dependence of $m^{\star}_{e}$ on $T$ is weak, as it is both in the isotropic case and in the momentum approximation Szczȩśniak (2006b); Szewczyk et al. (2018).

The thermodynamics of the superconducting state is fully determined by the values of both the order parameter and the wave function renormalising factor. Considering the averaged values of the discussed functions, we calculated the difference in free energy between the superconducting and the normal state Bardeen and Stephen (1964): $\Delta F/\rho\left(0\right)=-\left(2\pi/\beta\right)\sum^{M}_{n=1}[\sqrt{\omega^{2}_{n}+\left<\Delta_{n}\right>^{2}}-|\omega_{n}|][\left<Z^{\left(S\right)}_{n}\right>-\left<Z^{\left(N\right)}_{n}\right>\frac{|\omega_{n}|}{\sqrt{\omega^{2}_{n}+\left<\Delta_{n}\right>^{2}}}]$. Symbols $Z^{\left(S\right)}_{n}$ and $Z^{\left(N\right)}_{n}$ denote the values of the wave function renormalising factor in the superconducting and the normal state, respectively. The thermodynamic critical field should be determined from the formula: $H_{C}=\sqrt{-8\pi\Delta F}$, whereas the difference in the specific heat between the superconducting and the normal state can be found using the formula: $\Delta C=C^{S}-C^{N}=-T\frac{d^{2}\Delta F}{dT^{2}}$, where: $C^{N}=\gamma{T}$. The Sommerfeld constant is equal to: $\gamma=\frac{2}{3}\pi^{2}k_{B}^{2}\rho(0)\left(1+\lambda\right)$. The achieved results we showed in Fig. 5. Then we calculated the values of the non-dimensional thermodynamic ratios: $R_{C}=\left[C^{S}\left(T_{C}\right)-C^{N}\left(T_{C}\right)\right]/C^{N}\left(T_{C}\right)$ and $R_{H}=T_{C}C^{N}\left(T_{C}\right)/H_{C}^{2}\left(0\right)$. The considered results correspond to the universal constants in the BCS theory, the latter taking the values equal to $1.43$ and $0.168$, respectively Bardeen et al. (1957a, b). For the superconducting state on a square lattice, we arrived at the results which differ significantly from the ones occurring in the BCS theory: $R^{0}_{C}=2.4$, $R^{0}_{H}=0.227$; $R^{0.12\gamma_{C}}_{C}=1.89$, $R^{0.12\gamma_{C}}_{H}=0.283$; $R^{0.24\gamma_{C}}_{C}=1.28$, $R^{0.24\gamma_{C}}_{H}=0.413$; $R^{0.31\gamma_{C}}_{C}=0.89$, $R^{0.31\gamma_{C}}_{H}=0.595$.

To summarize, we proved that – under the assumptions made with respect to the electron–phonon coupling function, vertex corrections, and anharmonic effects – the balanced phonon–induced superconducting state cannot be generated on a square lattice. On the other hand, the linear electron–phonon interaction can induce the unbalanced superconducting state in cases for which the unbalance parameter takes a value less than $0.42$. This value is much lower than $\gamma_{C}$ estimated in the static boundary ($0.93$). It means that the dynamic effects modeled by the explicit dependence of the order parameter and the wave function renormalising factor on the Matsubara frequency are very important for the superconducting state on a square lattice. This result matched our expectations, given that the interaction between electrons is transmitted via slow phonons. The reason for a decrease in the critical temperature value with an increase in $\gamma$ parameter value is the anomalously high increase in the effective electron mass. It should be rightly stressed that the possibility of induction of the superconducting state on a square lattice most probably depends strongly on the form of the electron (phonon) dispersion relation, the electron–phonon coupling function, or the assumed approximations. Therefore one cannot claim that the electron–phonon interaction in low–dimensional systems (2D) is unimportant with respect to the superconducting state. Our result evidence that the analysis of the thermodynamical properties of the superconducting phase, if too cursory, e.g. within the isotropic approximation, can generate false results. This conclusion is of particular importance in the case of the most recent low–dimensional systems (decorated graphene Profeta et al. (2012); Pešić et al. (2014); Guzman et al. (2014); Kaloni et al. (2014), silicene Wan et al. (2013) or phosphorene Ge et al. (2015); Shao et al. (2014)), in which the superconducting state can be induced.

The performed numerical analysis proved that the unbalanced superconducting state is described by the thermodynamic parameters, which values differ significantly from the values predicted by the BCS theory.

Please notice that our results can make the understanding of the pairing mechanism in cuprates considerably easier Dagotto (1994); Szczȩśniak (2012). It is generally accepted fact that electrons in cuprates form a strongly correlated system Emery (1987); Littlewood et al. (1989); Hybertsen et al. (1990). Numerical calculations carried out for the Hubbard model demonstrate that the value of the on-site $U$ integral is equal to about $5$ eV Hybertsen et al. (1990). It should be strongly emphasised that the balanced linear electron–phonon interaction is too weak to induce the experimentally observed superconducting state (the full discussion of the considered problem at the ab initio level can be found in Bohnen et al. (2003)). On the other hand, there exists a good deal of experimental data which indicate that the interaction between electrons and phonons in cuprates is really significant. Let us recall the results obtained by means of the ARPES method, which demonstrated the existence of a break in the energy spectrum near the phonon energy Damascelli et al. (2003); Cuk et al. (2005). Additionally, the ARPES method made possible to discover the isotope effect of the real part of the self–energy Gweon et al. (2004). The isotope effect is also related to the critical temperature, what is particularly distinct in the strongly underdoped regions Franck (1994). Moreover, the vibrations of the crystal lattice modify the penetration depth and the results of Raman measurements Hofer et al. (2000); Schneider (2005).

Therefore an obvious suggestion comes into mind that the high–temperature superconducting state in cuprates may be induced by the strongly unbalanced electron–phonon interaction. This suggestion is quite probable for the reason that it is very difficult to prove the existence of the superconducting state of sufficiently high critical temperature value in pure electron models with $U>0$. Such a problem, however, does not occur in effective models with $U<0$. Our hypothesis is also confirmed by the results achieved by Kim and Tesanovic Kim and Tesanovic (1993), who showed, on the example of ${\rm La_{2-x}Sr_{x}CuO_{4}}$, ${\rm(Y_{1-x}Pr_{x})Ba_{2}Cu_{3}O_{7-y}}$ and ${\rm YBa_{2-x}La_{x}Cu_{3}O_{7}}$ compounds, that the strong Coulomb correlations do not suppress the phonon pairing mechanism over a wide doping range.

Appendix A Dependence of the order parameter value on the number of momentum lattice points and the number of Matsubara frequencies

We checked the number of momentum lattice points and the number of Matsubara frequencies for which the value of the order parameter is saturated in the temperature $k_{B}T_{0}=0.0001$ $t$. We plotted the results in Fig. 6. It can be assumed with good accuracy that the averaged value of the order parameter is saturated for $N\sim 79$ and $n\sim 67$, under the assumption that $n=200$ and $N=200$, respectively. It is worth remembering that the value of the order parameter will be saturated much faster for higher temperature values.

Appendix B The unbalanced superconducting state on a square lattice – the isotropic approximation

Let us discuss the issue of the unbalanced superconducting state within the isotropic approximation. The Eliashberg equations take the form:

[TABLE]

where: $p_{\mbox{{\scriptsize{\bf k}}}}\left(m\right)=\left(Z_{m}\omega_{m}\right)^{2}+\varepsilon^{2}_{\mbox{{\scriptsize{\bf k}}}}+\varphi_{m}^{2}$. The pairing kernel of the electron–phonon interaction is given by: $K^{\prime}\left(n-m\right)=\lambda\frac{\nu^{2}}{\left(n-m\right)^{2}+\nu^{2}}$, where $\lambda$ is the electron–phonon coupling constant. We used Kresin’s method of introducing the average phonon frequency: $\left<\Omega\right>\sim\omega_{D}$. In this case: $\nu=\beta\omega_{D}/2\pi$ Kresin (1984, 1987). The density of states per spin direction $\rho\left(\varepsilon\right)$ is reproduced by: $\rho\left(\varepsilon\right)=b_{1}\ln\left|\varepsilon/b_{2}\right|$, where $b_{1}=-0.04687t^{-1}$ and $b_{2}=21.17796t$ Szczȩśniak (2006b). Hence we obtain: $\sum_{\mbox{{\scriptsize{\bf k}}}}p_{\mbox{{\scriptsize{\bf k}}}}^{-1}\left(m\right)\simeq\int^{W}_{-W}d\varepsilon\rho\left(\varepsilon\right)p_{\varepsilon}^{-1}\left(m\right)$, where $W=4t$ is half of the band width.

The dependence of the order parameter on $\gamma$ is plotted in Fig. 7 (a). It can be seen that, in the case of the isotropic approximation, the superconducting state on the square lattice is induced even for high values of the unbalance parameter. The results concerning the effective electron mass are presented in Fig. 7 (b). The values of $m^{\star}_{e}$ in the isotropic approximation are not very high, as compared with the accurate results (see Fig. 2 (c)).