1 Introduction

Since its discovery by Kamerling-Onnes [1] at the beginning of the twentieth century superconductivity became a particularly important field of physics with wide technical applications. The identification of superconductors as perfect diamagnets, i.e. the Meissner effect [2] was discovered in the thirties and was soon followed by the beautiful phenomenological electromagnetic theory of London [3, 4]. The fundamental theoretical breakthrough however is due to Bardeen, Cooper and Schrieffer [5]. They have shown that the origin of the superconductive phase transition lies in the correlation between electrons of opposite momenta and spin resulted from phonon exchange. While the most important and striking feature is the absence of resistance below a critical temperature, there is no deep understanding of it. Within the present unsatisfactory status of non-equilibrium statistical mechanics irreversibility (dissipation) is introduced “by hand”, although there are significant recent progresses in understanding the treatment of open systems [6, 7]. Therefore, the understanding at least of the equilibrium properties is a central point. The most important in this respect is the theory of the Meissner effect to which we devote our discussion.

The standard modeling of the BCS idea within self-consistent electron theories [8,9,10] offered no convincing results for the Meissner effect. The failure is due to the incorrect implementation of electromagnetism. Standard condensed matter theory includes just Coulomb interactions of charged particles. This offers no basis for the treatment of magnetic phenomena. One had to introduce magnetic interaction between localized spins to understand ferro-magnetism. The (dia)-magnetism resulted form the electronic currents remained still ignored. However, at least on the mean-field level this could be corrected by the introduction of the average internal magnetic field based on physical intuition (see [11]).

Starting from the non-relativistic quantum electrodynamics a current–current magnetic interaction may be derived [12, 13]. We show in this paper that its inclusion improves the standard bulk mean field theory explaining within this frame the perfect diamagnetism. In this sense we contradict our recent skepticism [13].

2 The standard theory of the Meissner effect in bulk

The revelation of BCS that the source of superconductivity lies in the correlation of electrons of opposite momenta and spin due to interaction with phonons was decisive for the theory. However, to pursue this idea within the many-body theory including phonons seems difficult. Therefore, one tried to construct pure electron many-body theories with a built-in “potential” giving rise to such correlations and implicitly to a superconducting phase transition.

Such a model Hamiltonian is due to Rickayzen [8, 9] that we describe here. Its advantage is the explicit coordinate-space formulation. The specific version of Bogolyubov-de Gennes [10] is included in this frame.

Since one is concentrated on equilibrium properties in a \(\mu \)-system one considers \(\mathcal{H}\equiv H-\mu N\) instead of the Hamiltonian H. This is compulsory if one wants to treat phase transitions like Bose condensation or superconductivity breaking particle fermion number conservation. Rickayzen introduces besides the kinetic energy in the presence of classical time-independent given vector \({{\mathbf \mathcal{A}}}(x)\) and scalar \(\phi (x)\) potentials, an electron–electron interaction by a “correlating” potential \(W(x)\), that produces no bound states, but gives rise to correlations of BCS type. The second quantized wave functions of the electrons \(\psi _{\sigma }(x)\) and \(\psi _{\sigma }^+ (x)\) obey the usual fermionic anti-commuting rules, while the field \(\mathcal{A}(x)\) is a classical one and the operator \(\mathcal{H}\) (in the absence of an external scalar potential is defined as

$$\begin{aligned}&\mathcal{H} \equiv H-\mu N \nonumber \\&\quad = \sum _{\sigma =\pm \frac{1}{2}} \int {\text {d}}x\psi _{\sigma }^{+}(x) \bigg \{ \frac{1}{2m}\left( -\imath \hbar \nabla -\frac{e}{c}\mathbf \mathcal{{A}}(x)\right) ^{2} \nonumber \\&\quad \quad -\mu \bigg \} \psi _{\sigma }(x) +\frac{1}{2}\int {\text {d}}x\int {\text {d}}{x}'W(x-{x}')\nonumber \\&\quad \quad \times \left[ \psi _{\frac{1}{2}}^{+}(x)\psi _{-\frac{1}{2}}^{+}({x}')\psi _{\frac{1}{2}}({x}')\psi _{-\frac{1}{2}}(x)+h.c.\right] . \end{aligned}$$
(1)

This Hamiltonian is invariant against time-independent gauge transformations of the vector potential \({\mathbf \mathcal{A}}({\mathbf {x}})\). We shall keep the discussion in a rather general frame without choosing a definite potential \(W(x)\) and ignore the Coulomb interactions (including a positive background) since they play no role in the next steps to follow. Further, one resorts to a self-consistent Hartree–Fock approximation including anomalous averages like \(\langle \psi _{\frac{1}{2}}^{+}(\mathbf {x})\psi _{-\frac{1}{2}}^{+}(\mathbf {x'})\rangle \).

Then, the self-consistent Hamiltonian is

$$\begin{aligned} \mathcal{H}_{s.c.} =&\sum _{\sigma =\pm \frac{1}{2}} \displaystyle \int {\text {d}}x\psi _{\sigma }^{+}(x) \bigg \{ \frac{1}{2m}\left( -\imath \hbar \nabla -\frac{e}{c}\mathbf \mathcal{{A}}(x)\right) ^{2}\nonumber \\&\quad -\mu \bigg \} \psi _{\sigma }(x) +\frac{1}{2} \displaystyle \int {\text {d}}x\displaystyle \int {\text {d}}{x}'W(x-{x}') \nonumber \\&\quad \times \left[ \langle \psi _{\frac{1}{2}}^{+}(x)\psi _{-\frac{1}{2}}^{+}({x}')\rangle \psi _{-\frac{1}{2}}({x}')\psi _{\frac{1}{2}}(x)\right. \nonumber \\&\quad + \left. \langle \psi _{-\frac{1}{2}}({x}')\psi _{\frac{1}{2}}(x)\rangle \psi _{\frac{1}{2}}^{+}(x)\psi _{-\frac{1}{2}}^{+}({x}')\right. \nonumber \\&\quad \left. -\langle \psi _{\frac{1}{2}}^{+}(x)\psi _{-\frac{1}{2}}^{+}({x}')\rangle \langle \psi _{-\frac{1}{2}}({x}')\psi _{\frac{1}{2}}(x)\rangle \right] . \end{aligned}$$
(2)

One can show within this frame, that in the absence of the field \({\mathbf \mathcal{A}}({\mathbf {x}})\) a phase transition may occur below a critical temperature, provided the potential \(W({\mathbf {x}})\) ensures a non-vanishing solution for the symmetry breaking gap parameter (\(\varDelta ({\mathbf {k}})\ne 0\)) of the largely described “gap equation” we do not give here. This condition is equivalent to the vanishing of the first derivative of the free energy with respect to \(\varDelta ({\mathbf {k}})\).

The next step is to use equilibrium linear response theory to get the relationship between the average \(\langle {\mathbf {j}}(x)\rangle \) of the current density operator

$$\begin{aligned} \mathbf {j}(x)\equiv \frac{e}{2m}\psi ^{+}(x)\left( -\imath \hbar \nabla +\frac{e}{c}{\mathbf \mathcal{A}}(x)\right) \psi ({x}')+h.c. \end{aligned}$$
(3)

and the weak vector potential \({\mathbf \mathcal{A}}\) in the Coulomb gauge \(\nabla {\mathbf \mathcal{A}}(x)=0 \, {\text {and}} \, \phi (x)=0\).

A peculiarity of the linear response within self-consistent theories is that the deviation of the averages from their equilibrium values constitute so called induced perturbations. This means, that the true perturbation (to first order in \({\mathbf \mathcal{A}}\)) is

$$\begin{aligned} \mathcal{H}'_{s.c.} =&-\frac{1}{c} \displaystyle \int {\text {d}}x\mathbf {j}(x)\mathbf {\mathcal{A}}(x) + \displaystyle \int {\text {d}}x\int {\text {d}}{x}'W(x-{x}') \nonumber \\&\quad \times \left[ \eta (x,{x}')\psi _{-\frac{1}{2}}({x}')\psi _{\frac{1}{2}}(x)+h.c.\right] \end{aligned}$$
(4)

with

$$\begin{aligned} \eta (\mathbf {x},\mathbf {x}')\equiv \langle \psi _{\frac{1}{2}}^{+}(\mathbf {x})\psi _{-\frac{1}{2}}^{+}(\mathbf {x'})\rangle -\langle \psi _{\frac{1}{2}}^{+}(\mathbf {x})\psi _{-\frac{1}{2}}^{+}(\mathbf {x'})\rangle _{0} \end{aligned}$$
(5)

being the deviations of the anomalous averages from their values in the absence of the field \({\mathbf \mathcal{A}}\). The resulting linear relationship between the Fourier transforms of the two transverse vectors reads as

$$\begin{aligned} \langle \tilde{j}_{\mu }({\mathbf {k}})\rangle =\kappa (k)\tilde{\mathcal{A}}_{\mu }({\mathbf {k}}); \qquad (\mu =1,2,3), \end{aligned}$$
(6)

with the scalar coefficient \(\kappa (k)\) in an infinite homogeneous, isotropic system being a function of \(k=\sqrt{{{\mathbf {k}}}^2 }\). Its explicit expression for any k was calculated explicitly [8] by neglecting the induced perturbations.

One is interested in this relationship only for small wave vectors (slowly varying behavior in the coordinate space!). If in the absence of the perturbation one had an anomalous superconducting phase with a non-vanishing gap \(\varDelta ({\mathbf {k}})\) at \({\mathbf {k}}=0\) one gets (without induced perturbations [8]) that \(\kappa (0)\) is finite and strictly negative

$$\begin{aligned} \kappa (0)=-\frac{1}{c\varLambda }< 0 . \end{aligned}$$
(7)

The same result is described in the frame of the anomalous Green functions, under similar approximations by Schrieffer [14].

It may be shown however [15], that under the condition of stability of the superconducting phase (which amounts to a non-negative second derivative of the free energy with respect to the gap parameter \(\varDelta ({\mathbf {k}})\) at \({\mathbf {k}}=0\)) \(\varLambda \) is indeed positive and the contribution produced by \(\eta (x,{x}')\) does not change Rickayzen’s result. For \({\mathbf {k}}\ne 0\) no sound results are available, since it requires the full calculation considering the induced perturbation terms.

So far, the calculations are OK. The problem resides in the interpretation of these results. According to the old interpretation (see [16]) Eqs. 6, 7 have to be compared to the second London equation that reads as

$$\begin{aligned} \nabla \times {\mathbf {i}}=-\frac{1}{c\varLambda }{\mathbf {B}}\qquad {\text {or}} \qquad \mathbf {i}=-\frac{1}{c\varLambda }{\mathbf { A}} , \end{aligned}$$

with \(\mathbf {i}(x)\) being the macroscopic super-current density and \({\mathbf {A} }(x)\) the total macroscopic field in the superconductor. Thus the Meissner effect would be explained within this model Hamiltonian. Moreover, the parameter \(\varLambda \) might be interpreted as the London penetration length.

Unfortunately, such a reasoning is misleading for several reasons. If Meissner effect occurs, the total magnetic field \({\mathbf {B}}\) as well as the super-current are identically null in the bulk whenever the source of the external field lies outside the system ! The reference to the London equation is based on the assumed behavior of the kernel \(\kappa (k)\) for \(k\ne 0\). (See more in the discussion at the end of the paper.) We intend to proof something less ambitious, but implying no further assumptions about the kernel.

On the other hand, the classical field \({\mathbf \mathcal{A}}\) in this model has undefined sources and no magnetic field produced by the electrons is present at all in this theory. Therefore it cannot be identified with the total macroscopic magnetic vector potential \({\mathbf { A}}\). Besides, this Hamiltonian is invariant with respect to time independent gauge transformations and therefore \({\mathbf {\mathcal{A}}}\) has to be identified with the external field. In the frame of the non-relativistic QED the Hamiltonian is defined in a fixed gauge, namely in the Coulomb gauge for the quantized (“radiation”) electromagnetic field.

The primary task of the microscopical theory in the bulk is just to show, that in thermal equilibrium a homogeneous constant external magnetic field is compensated completely by the internal one produced by the electrons. Obviously, some ingredients are still missing. One is at the range of validity of the ordinary quantum mechanics, that takes no magnetic interactions between the electrons into account. Steps towards the non-relativistic quantum electrodynamics of charged particles are compulsory. This criticism of the electromagnetic aspects is pertinent also to the Ginzburg–Landau non-linear theory of superconductivity [17].

3 The current–current interaction and the theory of superconductivity

Already 100 years ago Darwin [18] has argued in the frame of the classical electrodynamics of point-like electrons, that up to order \(1/c^2\) one might separate the motion of the particles from that of the electromagnetic field. From this separation emerges a magnetic electronic current–current interaction. Since the classical electrodynamics of point-like charged particles is a vicious theory having neither Lagrangian, nor Hamiltonian formulation, his derivation lacked any rigor, nevertheless it contains a grain of truth and was considered later also by Landau and Lifshitz [19]. They have shown that Darwin implicitly has chosen a certain very strange non-linear choice of gauge. Such a gauge however imposes constraints on the velocities of the particles and therefore no ordinary canonical formalism is allowed. See in this context Dirac’s theory of canonical formalism with constraints [20, 21].

Actually one needs a new analysis of the problem in the frame of the non-relativistic QED (see [22, 23] for a modern presentation). A natural choice of the gauge is here the Coulomb one. Indeed, in this gauge one has to do only with the two physical degrees of freedom of the photons and one is free from artificial constraints.

The non-relativistic QED is defined by the Hamiltonian (here just for electrons and photons) in the presence of time-dependent external electric and magnetic potentials

$$\begin{aligned} H^{QED}(t)= & {} \sum _{\mathbf {q},\lambda }\hbar \omega _{q}b_{\mathbf {q},\lambda }^{+}b_{\mathbf {q},\lambda } \nonumber \\&+\int d{\mathbf {x}} \mathcal{N} \left[ \frac{1}{2m} \left( \imath \hbar \nabla \psi ^+({\mathbf {x}}) \right. \right. \nonumber \\&\left. + \frac{e}{c}({\mathbf {A}}_\bot ({\mathbf {x}}) + {\mathbf {A}}_{ext}({\mathbf {x}},t))\psi ^+({\mathbf {x}})\right) \nonumber \\&\times \left. \left( \frac{\hbar }{\imath }\nabla \psi ({\mathbf {x}}, t) +\frac{e}{c}({\mathbf {A}}_\bot ({\mathbf {x}})+{\mathbf {A}}_{ext}({\mathbf {x}},t)))\psi ({\mathbf {x}})\right) \right] \nonumber \\&+ \frac{1}{2}\int {\text {d}}\mathbf {x}\int d{\mathbf {x}}'\psi ^{+}(\mathbf {x})\psi ^{+}(\mathbf {x}')\frac{e^{2}}{|\mathbf {x}-\mathbf {x}'|}\psi (\mathbf {x}')\psi (\mathbf {x})\nonumber \\&+ e\int {\text {d}}\mathbf {x} V_{ext}(\mathbf {x},t)\psi ^+(\mathbf {x})\psi (\mathbf {x}) . \end{aligned}$$
(8)

with the transverse radiation field being

$$\begin{aligned} \mathbf {A}_{\bot }(\mathbf {x})=\sum _{\lambda =1,2}\sqrt{\frac{\hbar {c}}{\varOmega }}\sum _{\mathbf {q}}\frac{1}{\sqrt{|\mathbf {q}|}}\mathbf {e}_{\mathbf {q}}^{(\lambda )} e^{-\imath \mathbf {q}\mathbf {x}}\left( b_{\mathbf {q},\lambda }+ b_{-\mathbf {q},\lambda }^{+}\right) , \nonumber \\ \end{aligned}$$
(9)

while \({\mathbf {A}}_{ext}(x,t)\) is a classical external vector potential. The unit vectors \(\mathbf {e}_{\mathbf {q}}^{(\lambda )}\) are orthogonal to the wave vector \(\mathbf {q}\) as well as to each other and according to the general recipe of second quantization a normal ordering \(\mathcal{N}[... ]\) had to be introduced also with respect to the photon creation and annihilation operators \(b_{\mathbf {q},\lambda }^{+} \), \(b_{\mathbf {q}}\) in the Hamiltonian and the photon frequency is \(\omega _q= c|q|\).

This non-relativistic QED has been widely used in the quantum optics of atoms and solids (see a recent book [24] about).

Although in their basic paper “de Haas-van Alphen Effect and the Specific Heat of an Electron Gas”, Holstein, Norton and Pincus [25] already have shown, that the non-relativistic QED current–current interaction is the basic piece for understanding diamagnetism, their ideas and proofs are ignored in the treatment of superconductivity.

In order to obtain a pure electronic Hamiltonian one must restrict the discussion to the subspace of states without photons. A good object to perform such a discussion is either the S-matrix, or the theory of Green functions in terms of the Feynman diagrams. As it was shown recently [12, 13] the necessary steps are to neglect time-retardation in the photon propagator and ignore the “sea-gull” vertices . This restricts the validity of the resulting theory to order \(1/c^2\).

The resulting electronic Hamiltonian (here with time-independent external potentials) is

$$\begin{aligned} \mathbf{H}= & {} \int d{x} \psi ^{+}(x)\left[ \left( \frac{\hbar }{\imath }\nabla -\frac{e}{c}\mathbf {A}_{ext}(x)\right) ^2 + V_{ext}(x)-\mu \right] \psi (x)) \nonumber \\&+ \frac{1}{2} \int {\text {d}}x\int {\text {d}}{x}'\frac{\mathcal{N}\left[ \rho (x)\rho ({x}')\right] }{|x-{x}'|} \nonumber \\&- \frac{1}{2} \int {\text {d}}x\int {\text {d}}{x}'\frac{\mathcal{N}\left[ \mathbf {j}_\bot (x)\mathbf {j}_\bot ({x}')\right] }{c^2|x-{x}'|} . \end{aligned}$$
(10)

Here \(\rho (x)\) denotes the charge density operator

$$\begin{aligned} \rho (x)=e\psi ^{+} (x)\psi (x), \end{aligned}$$
(11)

while \(\mathbf {j}_\bot (\mathbf {x})\) denotes the transverse part of the current density operator

$$\begin{aligned} \mathbf {j}(x) = \frac{e}{2m}\left( \psi ^{+}(x) \left( \frac{\hbar }{\imath }\nabla \!-\!\frac{e}{c}\mathbf {A}_{ext}(x)\!\!\right) \psi (x) +\!h.c.\right) . \end{aligned}$$
(12)

It is essential that the kinetic energy contains only the external magnetic vector potential \(\mathbf {A}_{ext}(x)\). As we shall see later, the internal magnetic field is directly generated by the current densities. The presence of the covariant derivatives with respect to the external vector potential ensures the invariance with respect to time-independent gauge transformations of the external fields as it holds in the QED Hamiltonian Eq. 8 defined in the Coulomb gauge for the radiation field.

The current–current term is nothing else but the well-known Biot–Savart law resulting from the exchange of transverse photons responsible for magnetic forces. One may argue that due to the smallness of the velocities in the condensed matter such an \(1/c^2\) term may be neglected. This is obviously false. Our everyday experience teaches us, that a macroscopic number of slow electrons may create enormous magnetic fields. We shall show that exactly this term is necessary to complete the theory of the preceding Section leading to the exact cancellation of a total constant magnetic field in the bulk.

Modifying accordingly the Hamiltonian of Eq. 1 and taking into account the presence of just a time independent external field \({\mathbf {A}}_{ext}\) the complete BCS-Hamiltonian looks now as

$$\begin{aligned} \mathcal{H}^{new}= & {} \sum \limits _{\sigma =\pm \frac{1}{2}} \displaystyle \int {\text {d}}x\psi _{\sigma }^{+}(x) \bigg \{ \frac{1}{2m}\left( -\imath \hbar \nabla -\frac{e}{c}\mathbf {A}_{ext}(x)\right) ^{2} \nonumber \\&\quad \left. -\mu \bigg \} \psi _{\sigma }(x) +\frac{1}{2}\displaystyle \int {\text {d}}x\displaystyle \int {\text {d}}{x}'W(x-{x}') \right. \nonumber \\&\quad \times \left[ \psi _{\frac{1}{2}}^+(x)\psi _{-\frac{1}{2}}^+({x}')\psi _{\frac{1}{2}}({x}')\psi _{-\frac{1}{2}}(x)+h.c.\right] \nonumber \\&\quad -\frac{1}{2} \displaystyle \int {\text {d}}x\displaystyle \int {\text {d}}{x}'\frac{\mathcal{N}\left[ \mathbf {j}_\bot (x)\mathbf {j}_\bot ({x}')\right] }{c^2|x-{x}'|}. \end{aligned}$$
(13)

Then, within the mean-field approximation of the current–current term in Eq.13 the modification of Eq.2 is

$$\begin{aligned} \mathcal{H}_{s.c.}^{new}= & {} \sum \limits _{\sigma =\pm \frac{1}{2}} \displaystyle \int {\text {d}}x\psi _{\sigma }^{+}(x) \bigg \{ \frac{1}{2m}\left( -\imath \hbar \nabla -\frac{e}{c}\mathbf {A}_{ext}(x)\right) ^{2} \nonumber \\&- \left. \mu \bigg \} \psi _{\sigma }(x) +\frac{1}{2}\displaystyle \int {\text {d}}x\displaystyle \int {\text {d}}{x}'W(x-{x}') \right. \nonumber \\&\times \left[ \langle \psi _{\frac{1}{2}}^{+}(x)\psi _{-\frac{1}{2}}^{+}({x}')\rangle \psi _{-\frac{1}{2}}({x}')\psi _{\frac{1}{2}}(x)\right. \nonumber \\&+\left. \langle \psi _{-\frac{1}{2}}({x}')\psi _{\frac{1}{2}}(x)\rangle \psi _{\frac{1}{2}}^{+}(x)\psi _{-\frac{1}{2}}^{+}({x}')\right. \nonumber \\&\left. -\langle \psi _{\frac{1}{2}}^{+}(x)\psi _{-\frac{1}{2}}^{+}({x}')\rangle \langle \psi _{-\frac{1}{2}}({x}')\psi _{\frac{1}{2}}(x)\rangle \right] \nonumber \\&- \displaystyle \int {\text {d}}x\displaystyle \int {\text {d}}{x}'\frac{{\mathbf {j}}_\bot (x)\langle {\mathbf {j}}_\bot ({x}')\rangle }{c^2 |x-{x}'|} . \end{aligned}$$
(14)

In the above expression one may already identify the internal transverse vector field

$$\begin{aligned} {\mathbf {A}}_{int}({x},t)= \int {\text {d}}{x}'\frac{\langle {\mathbf {j}}_\bot ({x}',t)\rangle }{c|x-{x}'|} \end{aligned}$$
(15)

or in Fourier transform

$$\begin{aligned} \tilde{A}_{int}^{\mu }({\mathbf {k}})=\frac{4\pi }{ck^2 }\langle \tilde{j}_{\mu }({\mathbf {k}})\rangle . \end{aligned}$$
(16)

As before we choose the external vector potential \({\mathbf {A}}_{ext}\) again to be transverse (\(\nabla {\mathbf {A}}_{ext}=0)\). Its longitudinal part would have been anyway irrelevant for the magnetic field \( {\mathbf {B}}_{ext}=\nabla \times {\mathbf {A}}_{ext} \). Thus all the vectors are implicitly transverse an no need to mention it in the notations.

The terms of first order in \({\mathbf {A}}_{ext}\) differ from those of the previous approach only by the replacement

$$\begin{aligned} {\mathbf {j}{\mathbf \mathcal{A}} \quad \Rightarrow \quad \mathbf {j}}({\mathbf {A}}_{ext}+ {\mathbf {A}}_{int}) \end{aligned}$$

in the perturbation Eq. 4. This could have been obtained also on intuitive physical grounds [11]. However, it is important to stress, that this replacement would not be correct overall in the Hamiltonians of Eqs. 1, 2 , but just in the first order terms.

Therefore, due to the new induced term now we get instead of Eq. 6 with the same previously defined \(\kappa (k)\) the relationship

$$\begin{aligned} \langle \tilde{j}^{\mu }({\mathbf {k}})\rangle =\kappa (k)\left( \tilde{A}_{ext}^{\mu }({\mathbf {k}})+ \tilde{A}_{int}^{\mu }({\mathbf {k}})\right) ; \quad (\mu =1,2,3) . \end{aligned}$$
(17)

Multiplying this equation by \(\frac{4\pi }{ck^2}\) one recovers on the left-hand side again the internal field \(\tilde{A}^{int}_{\mu }({\mathbf {k}})\), therefore

$$\begin{aligned} {\displaystyle \tilde{A}_{int}^{\mu }({\mathbf {k}})=\frac{\frac{4\pi }{ck^2}\kappa (k)}{1-\frac{4\pi }{ck^2 }\kappa (k)}\tilde{A}_{ext}^{\mu }({\mathbf {k}}); \qquad (\mu =1,2,3) .} \end{aligned}$$
(18)

For the internal and external magnetic fields holds similarly

$$\begin{aligned} {\tilde{B}}^{\mu }_{int}({\mathbf {k}})=\frac{\frac{4\pi }{ck^2 }\kappa (k)}{1-\frac{4\pi }{ck^2 } \kappa (k)}{\tilde{B}}^{\mu }_{ext}({\mathbf {k}}); \quad (\mu =1,2,3). \end{aligned}$$
(19)

This last equation has been explicitly obtained by Tinkham [11] as already mentioned and may be also obtained through Zubarev’s early reasoning [26] based just on the macroscopic Maxwell equations, without any reference to a Hamiltonian. He obtains the linear response to the field \({\mathbf {B}}\) from the known linear response to the field \({\mathbf {H}}\).

Since it was proven already [8, 15] that under the condition of a stable superconductive phase \(\kappa (0)\) is finite, from Eq. 19 follows

$$\begin{aligned} {\tilde{B}}^{\mu }_{int}(0)=-{\tilde{B}}^{\mu }_{ext}(0) . \end{aligned}$$
(20)

This proves that in the frame of this modified s.c. BCS model in the stable superconductive phase no constant magnetic field in the bulk may survive! (Actually the vanishing of the constant magnetic field in real space still needs the easily accepted assumption that a homogeneous external field induces also a homogeneous magnetic field in the bulk.) The internal magnetic field of the electrons compensates fully the applied external magnetic field. The finite value of the coefficient \(\kappa (0)\) of Eq. 7 remains the criterion for the Meissner effect [16], however, without supplementary assumptions and approximations it cannot be related to the London penetration length and its sign is irrelevant for the proof.

4 Conclusions

By improving the standard bulk model of BCS superconductivity in its electromagnetic aspects, we have proven that in its mean-field approximation, if the superconductive phase transition is stable, the Meissner effect results in the sense of the perfect compensation of a constant external magnetic field by the internal magnetic field i.e. we have an ideal diamagnet. This confirms the intuitive argumentation of [11] Over the whole derivation no definite “correlating potential” \(W(x)\) was considered, but just the general requirement of the minimum of the free energy with a non-vanishing gap parameter. Our result implies no supplementary assumptions nor new approximations, but just the implementation of the ordinary non-relativistic QED up to order \(1/c^2\).

The true relevance of the Hamiltonians Eqs. 13,14 manifests itself in the non-linear “sea-gull” term \(\frac{e^2}{2mc^2}\psi ^+\)\( \psi {\mathbf \mathcal{A}}^2\) , where the replacement of the vector potential \({\mathbf \mathcal{A}}\) by the total average vector potential \({\mathbf {A}}_{ext}+{\mathbf {A}}_{int} \), defying simple intuition, would be wrong. Of course, beyond the mean field approximation and linear response this aspect must be even more important.

Of course, a microscopic derivation of the London equations remains highly desirable, but within this bulk approach it is an almost impossible task. One needs anyway to consider a specific potential \(W(x)\) and calculate the correct linear response for any k including also the induced terms. One must consider then a local external magnetic field embedded in the bulk (a solenoid). On the other hand, a local relation in the \({\mathbf {k}}\)-space does not imply a local (or at least an almost local) relation in the \(x\)-space. One might consider alternatively the boundary problem (in a half space) with a homogeneous magnetic external field. Again linear response in coordinate space without strongly localized correlations leads nowhere. It seems plausible that one needs one more ingredient, destroying spatial long-range correlations typical for phase transitions, without destroying the phase transition itself.

In this context it is useful to remind the situation with the Debye electric field penetration in a semiconductor. It may be obtained in a model of free classical electrons with a s.c. potential in equilibrium (Maxwell–Boltzmann distribution) within the linear approximation. However, a quantum mechanical derivation within linear response needs short range correlations in real space, that the free electron model cannot deliver. Here also some ingredient is necessary, although it is easier to imagine one since no phase transition is involved. Nevertheless, an analogous relation to Eq. 19 may be derived for the electric field with a non-vanishing kernel at \(k=0\) showing that in equilibrium no constant electric field survives in bulk.

To conclude, with the implementation of the current–current interaction, the standard mean field theory [9] can explain ideal diamagnetism in bulk, but not yet the London equation.