*Article* **The Symmetry of Pairing and the Electromagnetic Properties of a Superconductor with a Four-Fermion Attraction at Zero Temperature**

#### **Przemyslaw Tarasewicz**

Faculty of Pharmacy, Collegium Medicum in Bydgoszcz Nicolaus Copernicus University in Toru ´n, ul. Jagiello ´nska 13, 85-067 Bydgoszcz, Poland; tarasek1@cm.umk.pl; Tel.: +48-52-585-3428; Fax: +48-52-585-3308

Received: 20 September 2019; Accepted: 22 October 2019; Published: 2 November 2019

**Abstract:** Properties of a fermion system at zero temperature are investigated. The physical system is described by a Hamiltonian containing the BCS interaction and an attractive four-fermion interaction. The four-fermion potential is caused by attractions between Cooper pairs mediated by the phonon field. In this paper, the BCS interaction is assumed to be negligible and the four-fermion potential is the only one that acts in the system. The effect of the pairing symmetry used in the four-fermion potential on some zero-temperature properties is studied. This especially concerns the electromagnetic response of the system to an external magnetic field. It turns out that, in this instance, there are serious differences between the conventional BCS system and the one investigated in this paper.

**Keywords:** superconductivity; four-fermion attraction; Meissner effect; fermion quartets

**PACS:** 74.20.-z; 74.20.Fg

#### **1. Introduction**

Recent decades have been witnessing many new phenomena discovered in solid-state physics (e.g., high temperature superconductivity, heavy fermion superconductors and unusual properties of <sup>3</sup>*He*). These discoveries are intriguing and lead to many new questions on the nature of superconductivity and superfluidity. The problem of the internal symmetry of gap parameters stands for one such an example. Furthermore, we still do not know how many particles are engaged in constituting the fundamental clusters responsible for the occurrence of new phases. Is it possible that only two-particle entities are relevant to these systems (e.g., Cooper pairs) or may three or four-particle structures also be taken into account? The last question was addressed in [1]. This issue becomes more and more interesting due to some recent discoveries and suggestions. To give some examples, let us start with a work by Schneider and Keller [2] who measured the relationship between the critical temperature and zero temperature condensate density for some cuprates and Chevrel-phases superconductors. They noticed that the experimental data for *YBa*2*Cu*3*O*6.602, for example, ca be associated with the behavior of a dilute Bose gas. As a consequence, they put the suggestion that Bose condensation of weakly interacting fermion pairs could be a possible mechanism of transition from the normal to the superconducting state. Moreover, a subsequent discovery of Bunkov et al. [3] points to the presence of fermion quartets in <sup>3</sup>*He*. They investigated the problem of the influence of spatial disorder on the order parameter in superfluid <sup>3</sup>*He*. They followed Volovik [4] and suggested that unusual spectra of <sup>3</sup>*He* in an aerogel could be explained by a process in which impurities tend to destroy the anisotropic correlations of the order parameter. Especially interesting is that the correlation of higher symmetry can survive (e.g., four-particle correlations). Two papers [5,6] reported on a discovery of the half-*h*/2*e* magnetic flux quanta coexisting together with the usual ones in SQUIDs fabricated of bicrystalline *YBa*2*Cu*3*O*7−*<sup>δ</sup>* films. As is widely known, such a circumstance corresponds to the

presence of fermion quartets in a physical system. This in turn leads to taking the interplay between Cooper pairs and quartets into consideration. It should be mentioned that there were already some attempts for introducing the four-fermion interactions in the field of nuclear physics [7].

The problem of symmetry of pairing in superconductors has a long history. In classical papers [8,9], this topic was established and studied to a large extent. As is known in the classical BCS theory, the BCS potential binds electrons into Cooper pairs in which pairing has the *s*-wave symmetry, i.e., the energetic gap that opens at the Fermi level is isotropic in the momentum space. This kind of pairing symmetry is displayed mostly by metallic superconductors. However, there have appeared novel superconductors, e.g., copper oxides and heavy-fermion superconductors. In many of these materials, the order parameter turned out to be anisotropic in the momentum space, i.e., it vanished along some directions. In cuprates *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing occurs [10] whereas the mixed pairing symmetry *<sup>s</sup>* <sup>+</sup> *id* is frequently present in heavy-fermion systems [11]. In this paper, we deal with three cases: the pure *s*-wave pairing, the pure *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing and the mixed *<sup>s</sup>* <sup>+</sup> *idx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing symmetries. We also consider two kinds of electronic bands, i.e., the band with the rectangular density of states (DOS) with the parabolic dispersion relation (the two-dimensional system) and the band with the cosine dispersion relation (the one- and two-dimensional systems). We calculate the order parameter, the ground state energy per lattice site and the electromagnetic kernel as the quantity describing the response of the system to an external magnetic field.

#### **2. The Model**

In this paper, some ground state properties of a fermion system, especially in the case of the half-filled band case, are investigated. The Hamiltonian of the system reads

$$H = H\_{\rm BCS} + V\_{\rm 4} \tag{1}$$

where

$$H\_{\rm BCS} = \sum\_{\mathbf{k}\sigma} \xi\_{\mathbf{k}} a\_{\mathbf{k}\sigma}^{\*} a\_{\mathbf{k}\sigma} - N^{-1} \sum\_{\mathbf{k}, \mathbf{k}'} G\_{\mathbf{k} \mathbf{k}'} a\_{\mathbf{k} +}^{\*} a\_{-\mathbf{k} -}^{\*} a\_{-\mathbf{k}' -} a\_{\mathbf{k}' +} \cdots$$

The reduced four-fermion interaction *V*<sup>4</sup> has a form

$$V\_4 = -N^{-1} \sum\_{\mathbf{k}, \mathbf{k'}} g\_{\mathbf{k} \mathbf{k'}} b\_{\mathbf{k}}^\* b\_{-\mathbf{k}}^\* b\_{-\mathbf{k'}} b\_{\mathbf{k'}} \tag{2}$$

where

$$b\_{\mathbf{k}} = a\_{\mathbf{k} + a\_{\mathbf{k} - 1}}$$

and operators *a***k***σ*, *a*<sup>∗</sup> **<sup>k</sup>***<sup>σ</sup>* denote the annihilation and creation operators, respectively. *N* is the number of lattice sites and *ξ***<sup>k</sup>** denotes the one-electron energy counted with respect to the chemical potential *μ*. The coupling functions *G***kk** and *g***kk** are assumed to be nonzero in the whole band as it should be in many novel superconductors that are narrow-band systems. As we remember in classical superconductors they are restricted to a narrow shell around the Fermi level. Moreover, it is assumed that both of potentials are attractive so the system is a mixture of fermion pairs with opposite momenta and spins as well as quartets made of these pairs.

In [12], we derived the Hamiltonian (1) by making use of a Fröhlich type transformation of the electron-phonon Hamiltonian. Moreover, we were able to find the sign and the approximated form of coupling constant of the four-fermion interaction.

$$\mathcal{S}\mathbf{k}\_{\mathbb{F}}\mathbf{k}\_{\mathbb{F}} \approx \frac{D\_{\mathbf{k}\_{\mathbb{F}}}^{6}}{\hbar^{5}\omega\_{\mathbf{k}\_{\mathbb{F}}}^{5}},$$

where *D***k***<sup>F</sup>* and *h*¯ *ω***k***<sup>F</sup>* are the electron-phonon coupling constant and the phonon energy at Fermi level, respectively. The relationship between this coupling constant and BCS coupling constant was estimated to be

$$\mathcal{g}\_{\mathbf{k}\_{\mathbf{F}}\mathbf{k}\_{\mathbf{F}}} = \frac{G\_{\mathbf{k}\_{\mathbf{F}}\mathbf{k}\_{\mathbf{F}}}^3}{\hbar^2 \omega\_{\mathbf{k}\_{\mathbf{F}}}^2}$$

This relationship points to the weak character of the four-fermion interaction in comparison to familiar Cooper pairing. In classical superconductors, quartets cannot be observed due to weak BCS coupling and strong Coulomb repulsion between electrons. However, in materials with strong coupling between electrons and phonons, it seems to be quite reasonable to expect that this phenomenon could be visible and recognizable especially from surveys of magnetic flux quanta. It is widely accepted that if half-*h*/2*e* magnetic flux quanta appear among usual ones it points to the existence of quartets in the investigated system. As was established in Introduction these exotic half-*h*/2*e* magnetic flux quanta are present in some novel superconductors. These materials possess narrow conduction bands and the antiadiabatic limit ( *<sup>h</sup>*¯ *<sup>ω</sup> <sup>t</sup>* 1), where *<sup>h</sup>*¯ *<sup>ω</sup> <sup>t</sup>* is the energy of an optic phonon while *t* is the hopping parameter. In this case the starting Hamiltonian could be as follows

$$H = H\_0 + H\_{I\prime}$$

$$H\_0 = \sum\_{i \neq j \neq \sigma} t\_{ij} c\_{i\sigma}^\* c\_{j\sigma} + (E - \mu) \sum\_i n\_i + \sum\_{i,j} \mathcal{U}\_{ij} n\_i n\_j$$

$$H\_I = \sum\_{i,j} g\_{i\bar{j}}^s n\_i (b\_{\bar{j}}^\* + b\_{\bar{j}}) + \sum\_{i,j} g\_{i\bar{j}}^p n\_i n\_j (b\_i^\* + b\_{\bar{i}}),$$

where *c*∗ *<sup>i</sup><sup>σ</sup>* and *ci<sup>σ</sup>* represent creation and annihilation operators of a local electron from the narrow band. *<sup>σ</sup>* = ± denotes the spin of electrons and *<sup>i</sup>* and *<sup>j</sup>* refer to lattice sites. *ni* = *ni*<sup>+</sup> + *ni*<sup>−</sup> is the number operator. *b*∗ *<sup>i</sup>* and *bi* are creation and annihilation operators for a local phonon residing on the *i*-th site. *tij* denotes the hopping integral for electrons in the band. *E* refer to the site energy of electrons. *Uij* is the inter-site Coulomb interaction (*<sup>i</sup>* <sup>=</sup> *<sup>j</sup>*) and the local one (*<sup>i</sup>* <sup>=</sup> *<sup>j</sup>*). *<sup>g</sup><sup>s</sup> ij* and *<sup>g</sup><sup>p</sup> ij* are the coupling constants of the inter-site pnonon–electron potential and the inter-site phonon-electron-electron potential (*i* = *j*). For *i* = *j* one gets the local interactions of electrons with phonons. The second term in *HI* is based on the argument given by Hirsch et al. [13] that phonons can affect an electron pair sitting on the same lattice site. This is associated with the fact that such a pair affects the shape of an orbital that expands in the space. We extend this argument that two interacting electrons residing on two sites being very close to one another can be influenced by phonons. The term *HI* can be removed via making use of the following generalized Lang–Firsov type transformation *U<sup>S</sup>* with the generator

$$S = \sum\_{i,j} \left( \frac{\mathcal{S}\_{ij}^{\boldsymbol{\varsigma}}}{\hbar \omega} n\_i + \frac{\mathcal{S}\_{ij}^{\boldsymbol{P}}}{\hbar \omega} n\_i n\_j \right) \left( b\_j^{\boldsymbol{\varsigma}} - b\_j \right)^2$$

As a consequence there appear some new purely electronic nonlocal terms, namely, three- and four-electron interactions. We perceive one more possibility to derive a Hamiltonian containing the four-fermion interaction. If one considers the so-called Penson–Kolb–Hubbard model which is described by the following Hamiltonian

$$H\_{pkh} = (E - \mu) \sum\_{i\sigma} n\_{i\sigma} + lI \sum\_{i} n\_{i+} n\_{i-} + W \sum\_{i} n\_{i} n\_{i+1} + $$

$$-t^p \sum\_{i} (c\_{i+}^\* c\_{i-}^\* c\_{i+1-} c\_{i+1+} + c\_{i+1+}^\* c\_{i+1-}^\* c\_{i-} c\_{i+}) - t^s \sum\_{i\sigma} (c\_{i\sigma}^\* c\_{i+1\sigma} + c\_{i+1\sigma}^\* c\_{i\sigma}),$$

where *U* and *W* are coupling constants for intra-site and inter-site Coulomb interactions. −*t <sup>s</sup>* and <sup>−</sup>*t<sup>p</sup>* are single-electron hopping and pair-hopping integrals, respectively. We admit different signs of *U* and *W*. If the pair-hopping term is weak in comparison to the other potentials, then one can treat it as a perturbation and find the expansion in powers of *tp*. We suspect that the four-electron term would appear in the second order of this expansion.

One needs to mention that the problem of the four-fermion interaction in superconductors was investigated in a series of papers by Szcze¸sniak et al. [ ´ 14–16]. In some of them, the Eliashberg formalism was applied in order to deal with cuprates as systems with strong coupling of electrons to phonons.

In this paper, we neglect the BCS interaction. This could seem to be strange to do so, however, in [17], it was discovered that for the sufficiently strong four-fermion attraction, the BCS interaction gets turned off. The system in the superconducting state behaves like it would be a condensate of fermion quartets only. The Cooper pairs disappear because they are forced to combine and form fermion quartets. Of course, we are aware of the fact that the scenario based on the four-fermion interaction as the only potential in the system cannot be valid for high-temperature superconductors. If it were the case, then the half-*h*/2*e* magnetic flux quanta would be dominating over the conventional ones. So far, this situation has not occurred. Moreover, the four-fermion interaction leads to the first order phase transition whereas novel materials undergo rather the second one. The critical temperatures are lower than those resulting from the BCS theory with the same value of the coupling function [18,19]. This suggests the weaker nature of the four-fermion attraction in comparison to the BCS one. The best solution of this problem is to consider the system with both of potentials and more realistic dispersion relations as it is made in this paper. In [18,19], we make use of the parabolic dispersion relation with the density of states for free electrons that does not seem to be adequate for novel superconductors. Therefore this work describes the limiting case only. The difficulties significantly grow in the case of presence of two attractive potentials in the system.

#### **3. The Mean Field Approximation**

In this section, we would like to show the structure of the spectrum of the extended BCS model represented by the Hamiltonian (1). In that paper, we try to address mainly the problem of electrodynamics of the Hamiltonian *H* defined in the Introduction. The mean field method is used in order to obtain the ground state energy and two-order parameters for Cooper's pairs and quartets, respectively. Both gaps are assumed to be complex at this stage. Unfortunately, the final expressions turns out to be very complicated as we shall see later. Our mean field Hamiltonian reads

$$\begin{split} H\_{\mathrm{M}} &= \sum\_{\mathbf{k}>0} H\_{\mathrm{M}\mathbf{k}} \\ &= \sum\_{\mathbf{k}>0} \left( \mathfrak{J}\_{\mathbf{k}} \sum\_{\sigma} (n\_{\mathbf{k}\sigma} + n\_{-\mathbf{k}\sigma}) - \Delta\_{\mathrm{G}\mathbf{k}} (a\_{\mathbf{k}}^{\*} + a\_{-\mathbf{k}}^{\*}) - \Delta\_{\mathrm{G}\mathbf{k}}^{\*} (a\_{\mathbf{k}} + a\_{-\mathbf{k}}) - 2\Delta\_{\mathrm{g}\mathbf{k}}^{\*} \beta\_{\mathbf{k}} - 2\Delta\_{\mathrm{g}\mathbf{k}} \beta\_{\mathbf{k}}^{\*} + \mathbb{C}\_{\mathbf{k}} \right), \end{split} \tag{3}$$

where

$$a\_{\mathbf{k}} = a\_{-\mathbf{k}-} a\_{\mathbf{k}+} \qquad , \qquad \beta\_{\mathbf{k}} = b\_{-\mathbf{k}} b\_{\mathbf{k}}.$$

$$\mathbf{C\_{k}} = \Delta\_{\mathbf{G}\mathbf{k}} \sigma\_{\mathbf{k}}^{\*} + \Delta\_{\mathbf{G}\mathbf{k}}^{\*} \sigma\_{\mathbf{k}} + \Delta\_{\mathbf{g}\mathbf{k}} \tau\_{\mathbf{k}}^{\*} + \Delta\_{\mathbf{g}\mathbf{k}}^{\*} \tau\_{\mathbf{k}}$$

$$\Delta\_{\mathbf{G}\underline{p}} := N^{-1} \sum\_{\underline{k}'} \mathbf{G\_{kk'}} \sigma\_{\mathbf{k}'} \qquad , \qquad \Delta\_{\mathbf{g}\mathbf{k}} := N^{-1} \sum\_{\underline{k}'} \mathbf{g\_{kk'}} \tau\_{\mathbf{k}'} \tag{4}$$

with *<sup>σ</sup>G***<sup>k</sup>** <sup>=</sup> Tr *<sup>e</sup>*−*βHM***<sup>k</sup>** *<sup>α</sup>***<sup>k</sup>** Tr *<sup>e</sup>*−*βHM***<sup>k</sup>** and *<sup>τ</sup>g***<sup>k</sup>** <sup>=</sup> Tr *<sup>e</sup>*−*βHM***<sup>k</sup>** *<sup>β</sup>***<sup>k</sup>** Tr *<sup>e</sup>*−*βHM***<sup>k</sup>** in practice. We followed the standard procedure of Bogolyubov et al. The details of the procedure can be found in [20] as well. The sum in the Hamiltonian *HM* is over {**k** : **k** > 0} denoting the set of all 1-fermion momenta restricted to a definite half-space of *R*3.

Now it remains to diagonalize *HM***k**. *HM***<sup>k</sup>** acts in the 16-dimensional space *M***<sup>k</sup>** spanned by the vectors

$$\left| \left| n\_1 n\_2 n\_3 n\_4 \right\rangle := (a\_{\mathbf{k}+}^\*)^{n\_1} (a\_{\mathbf{k}-}^\*)^{n\_2} (a\_{-\mathbf{k}+}^\*)^{n\_3} (a\_{-\mathbf{k}-}^\*)^{n\_4} \left| 0 \right\rangle \tag{5}$$

where *ni* = 0, 1 for *i* = 1, 2, 3, 4. The diagonalization can be done by the method presented in [21,22]. It consists in splitting *M***<sup>k</sup>** into invariant subspaces with fixed eigenvalues of the spin projection operator 2*S***<sup>k</sup>** and two seniorities Λ**k**+, Λ**k**<sup>−</sup> defined as

$$\Delta S\_{\mathbf{k}} = \sum\_{a=\pm 1} \sum\_{\sigma=\pm 1} \sigma n\_{a\mathbf{k},\sigma} \quad , \quad \Lambda\_{\mathbf{k}\sigma} = n\_{\mathbf{k},\sigma} - n\_{-\mathbf{k},-\sigma} \quad , \quad \sigma = \pm \tag{6}$$

which commute with *HM***k**, namely,

$$\left[H\_{M\mathbf{k}\prime} \mathcal{Z} \mathbf{S}\_{\mathbf{k}}\right] = 0 \quad , \ \left[H\_{M\mathbf{k}\prime} \Lambda\_{\mathbf{k},\sigma}\right] = 0. \tag{7}$$

*M***<sup>k</sup>** splits into nine such invariant subspaces *M***k***<sup>i</sup>* (*i* = 1, 2, ...9) and due to the relations (7) diagonalization of *HM***<sup>k</sup>** can be performed separately in each of them.

	- 1. |1010 2*s* = 2 *λ*<sup>+</sup> = 1 *λ*<sup>−</sup> = −1 *E***<sup>k</sup>** = 2*ξ***<sup>k</sup>** 2. |0101 2*s* = −2 *λ*<sup>+</sup> = −1 *λ*<sup>−</sup> = 1 *E***<sup>k</sup>** = 2*ξ***<sup>k</sup>** 3. |1100 2*s* = 0 *λ*<sup>+</sup> = 1 *λ*<sup>−</sup> = 1 *E***<sup>k</sup>** = 2*ξ***<sup>k</sup>** 4. |0011 2*s* = 0 *λ*<sup>+</sup> = −1 *λ*<sup>−</sup> = −1 *E***<sup>k</sup>** = 2*ξ***<sup>k</sup>**
	- 5. |1000 |1110 <sup>2</sup>*<sup>s</sup>* = <sup>1</sup> *<sup>λ</sup>*<sup>+</sup> = <sup>1</sup> *<sup>λ</sup>*<sup>−</sup> = <sup>0</sup> *<sup>E</sup>***k**<sup>±</sup> 6. |0001 |0111 <sup>2</sup>*<sup>s</sup>* = −<sup>1</sup> *<sup>λ</sup>*<sup>+</sup> = −<sup>1</sup> *<sup>λ</sup>*<sup>−</sup> = <sup>0</sup> *<sup>E</sup>***k**<sup>±</sup> 7. |0010 |1011 <sup>2</sup>*<sup>s</sup>* = <sup>1</sup> *<sup>λ</sup>*<sup>+</sup> = <sup>0</sup> *<sup>λ</sup>*<sup>−</sup> = −<sup>1</sup> *<sup>E</sup>***k**<sup>±</sup> 8. |0100 |1101 <sup>2</sup>*<sup>s</sup>* = −<sup>1</sup> *<sup>λ</sup>*<sup>+</sup> = <sup>0</sup> *<sup>λ</sup>*<sup>−</sup> = <sup>1</sup> *<sup>E</sup>***k**<sup>±</sup>

where *<sup>E</sup>***k**<sup>±</sup> <sup>=</sup> <sup>2</sup>*ξ***<sup>k</sup>** <sup>±</sup> *EG***<sup>k</sup>** with *EG***<sup>k</sup>** = (*ξ*<sup>2</sup> **<sup>k</sup>** + |Δ*G***k**| <sup>2</sup>)1/2. In each of the subspaces *M***k***<sup>i</sup>* (*i* = 5, 6, 7, 8) the eigenproblem of *HM***<sup>k</sup>** reduces to that of the matrix

$$
\begin{pmatrix}
\mathbf{^XF}\_{\mathbf{k}} & \Delta\_{\mathbf{Gk}}^{\*} \\
\Delta\_{\mathbf{Gk}} & \mathbf{3}\_{\mathbf{g}}^{x}
\end{pmatrix}
\tag{8}
$$

The eigenvectors of *HM***<sup>k</sup>** in these subspaces have the form

$$\left| \left| E\_{\mathbf{k} \pm} \right> \right| = c\_{\mathbf{k} \pm} \left| m\_1 n\_2 n\_3 n\_4 \right> + d\_{\mathbf{k} \pm} \left| m\_1 m\_2 m\_3 m\_4 \right> \tag{9}$$

where <sup>4</sup> ∑ *i*=1 *ni* <sup>=</sup> 1, <sup>4</sup> ∑ *i*=1 *mi* = 3 and |*c***k**±| <sup>2</sup> <sup>=</sup> (*ξ***<sup>k</sup>** <sup>∓</sup> *EG***k**)<sup>2</sup> <sup>|</sup>Δ*G***k**|<sup>2</sup> + (*ξ***<sup>k</sup>** <sup>∓</sup> *EG***k**)<sup>2</sup> , <sup>|</sup>*d***k**±| <sup>2</sup> <sup>=</sup> <sup>|</sup>Δ*G***k**<sup>|</sup> 2 |Δ*G***k**|<sup>2</sup> + (*ξ***<sup>k</sup>** ∓ *EG***k**)<sup>2</sup>

(C) There is one 4-dimensional common subspace *M***k**<sup>9</sup> of 2*S***k**, Λ**k**<sup>±</sup> and *HM***<sup>k</sup>** spanned by the vectors |0000, |1001, |0110, |1111. The eigenvalue 2*S***k**, <sup>Λ</sup>**k**<sup>±</sup> in *<sup>M</sup>***k**<sup>9</sup> is 2*<sup>s</sup>* = *<sup>λ</sup>*<sup>+</sup> = *<sup>λ</sup>*<sup>−</sup> = 0. Let us use the following basis in *M***k**9:

$$\{|0000\rangle, -|1001\rangle, |0110\rangle, -|1111\rangle\}\tag{10}$$

and denote the projector on *M***k**<sup>9</sup> by *P***k**9. Therefore, in the basis (10) we have

$$P\_{\mathbf{k}\mathfrak{g}}H\_{\mathbf{k}}P\_{\mathbf{k}\mathfrak{g}} = \begin{pmatrix} 0 & \Delta\_{\mathbf{G}\mathbf{k}}^{\*} & \Delta\_{\mathbf{G}\mathbf{k}}^{\*} & 2\Delta\_{\mathbf{g}\mathbf{k}}^{\*} \\ \Delta\_{\mathbf{G}\mathbf{k}} & 2\mathfrak{g}\_{\mathbf{k}} & 0 & \Delta\_{\mathbf{G}\mathbf{k}}^{\*} \\ \Delta\_{\mathbf{G}\mathbf{k}} & 0 & 2\mathfrak{g}\_{\mathbf{k}} & \Delta\_{\mathbf{G}\mathbf{k}}^{\*} \\ 2\Delta\_{\mathbf{g}\mathbf{k}} & \Delta\_{\mathbf{G}\mathbf{k}} & \Delta\_{\mathbf{G}\mathbf{k}} & 4\mathfrak{g}\_{\mathbf{k}} \end{pmatrix} \tag{11}$$

Consecutively, the index **k** will be suppressed in this section. In terms of the unknown *x* = *E* − 2*ξ*, the secular equation for the eigenvalues *E* of *P***k**9*H***k***P***k**<sup>9</sup> takes the form

$$(\mathbf{x})(\mathbf{x}^3 + a\_2\mathbf{x}^2 + a\_1\mathbf{x} + a\_0) = \mathbf{0},\tag{12}$$

where

$$\begin{aligned} a\_0 &= -4(\Delta\_{\mathcal{S}} \Delta\_{G}^\* \Delta\_{G}^\* + \Delta\_{\mathcal{S}}^\* \Delta\_{G} \Delta\_{G}), \\ a\_1 &= -4(|\Delta\_{\mathcal{S}}|^2 + |\Delta\_{G}|^2 + \xi^2), \\ a\_2 &= 0. \end{aligned}$$

One root of Equation (12) is thus equal *x*(2) = 0, yielding *E*(2) = 2*ξ*. The remaining three can be determined by passing to the standard form of a 3rd-order equation:

$$x^3 + 3xp + 2q = 0\tag{13}$$

with

$$\begin{array}{l} 3p = a\_1 = -4(|\Delta\_{\mathcal{R}}|^2 + |\Delta\_G|^2 + \mathfrak{z}^2) \\ 2q = -4(\Delta\_{\mathcal{R}}\Delta\_G^\*\Delta\_G^\* + \Delta\_{\mathcal{R}}^\*\Delta\_G\Delta\_G) \end{array}$$

Equation (13) has the solutions

$$y\_k = 2(-1)^k \sqrt{r} \cos\left(\frac{\varphi}{3} + \frac{k\pi}{3}\right) \quad , \ k = 0, \pm 1 \tag{14}$$

where *<sup>r</sup>* <sup>=</sup> <sup>−</sup>*p*, cos *<sup>ϕ</sup>* <sup>=</sup> *tr*<sup>−</sup> <sup>3</sup> <sup>2</sup> with *t* = −*q*. These solutions yield the corresponding eigenvalues of *HM***<sup>k</sup>** in *M***k**9, viz.,

$$E^{(k)} = \mathfrak{Z}\_{\mathfrak{\*}}^{x} + y\_k \tag{15}$$

Let *u*(*j*), *v* (*j*) <sup>1</sup> , *v* (*j*) <sup>2</sup> , *<sup>s</sup>*(*j*) denote the components of the eigenvectors <sup>|</sup>*E*(*j*) of *<sup>P</sup>*9*HP*<sup>9</sup> in the basis (10):

$$P\_{\mathfrak{Y}}HP\_{\mathfrak{Y}}|E^{(j)}\rangle = E^{(j)}|E^{(j)}\rangle \tag{16}$$

along with

$$|E^{(j)}\rangle = u^{(j)}|0000\rangle - v\_1^{(j)}|1001\rangle + v\_2^{(j)}|0110\rangle - s^{(j)}|1111\rangle, \; j = 0, \pm 1, 2\tag{17}$$

One finds

$$|E^{(2)}\rangle = |2\xi\rangle = \frac{1}{\sqrt{2}}(|1001\rangle + |0110\rangle)$$

and the components of the remaining three eigenvectors <sup>|</sup>*E*(*j*) equal

$$|u^{(j)}|^2 = \frac{|a^{(j)}|^2}{D^{(j)}} \ , \ \ v\_1^{(j)} = v\_2^{(j)} \ , \ \ |v\_1^{(j)}|^2 = \frac{|b^{(j)}|^2}{D^{(j)}} , \ \ |s^{(j)}|^2 = \frac{|c^{(j)}|^2}{D^{(j)}} \ \ \ j = 0 \pm 1 \tag{18}$$

where

$$|a^{(i)}|^2 = (2\Lambda\_\mathcal{\mathbb{X}}^\*\Lambda\_\mathcal{G} - (4\mathfrak{z} - E^{(i)})\Lambda\_\mathcal{G}^\*)(2\Lambda\_\mathcal{\mathbb{X}}\Lambda\_\mathcal{G}^\* - (4\mathfrak{z} - E^{(i)})\Lambda\_\mathcal{G})(E^{(i)} - 2\mathfrak{z})^2,\tag{19}$$

$$|b^{(i)}|^2 = 4(\Delta\_{\mathcal{g}}^\* \Delta\_G \Delta\_G + \Delta\_{\mathcal{G}} \Delta\_G^\* \Delta\_G^\* + |\Delta\_G|^2 (E^{(i)} - 2\mathfrak{g}^\chi))^2,\tag{20}$$

$$|c^{(i)}|^2 = (E^{(i)}\Delta\_G + 2\Delta\_{\mathfrak{J}}\Delta\_G^\*)(E^{(i)}\Delta\_G^\* + 2\Delta\_{\mathfrak{J}}^\*\Delta\_G)(E^{(i)} - 2\mathfrak{J})^2,\tag{21}$$

$$D^{(i)} = |a^{(i)}|^2 + 2|b^{(i)}|^2 + |c^{(i)}|^2 \tag{22}$$

#### **4. The Case with the Zero BCS Potential**

As it has been stated in Introduction, we are interested in some zero-temperature properties of *H* in which the BCS interaction is discarded. Especially, this concerns the electromagnetic response to a weak external magnetic field. The objective is to show how the electromagnetic kernel depends on the symmetry of a quartet order parameter and on the type of one-electron density of states in the band. Before we do that, we will give some details regarding the spectrum and ground state properties of such a system. The mean field Hamiltonian is

$$\overline{H}\_{\mathcal{g}} = \sum\_{\mathbf{k} > 0} \left[ \xi\_{\mathbf{k}} \sum\_{\sigma} (n\_{\mathbf{k}\sigma} + n\_{-\mathbf{k}\sigma}) - 2\Delta\_{\mathbf{g}\mathbf{k}} (\beta\_{\mathbf{k}} + \beta\_{\mathbf{k}}^{\*}) + \mathbf{C}\_{\mathbf{k}} \right] = \sum\_{\mathbf{k} > 0} \overline{H}\_{\mathcal{g}\mathbf{k}\prime} \tag{23}$$

where Δ*g***<sup>k</sup>** is assumed to be real. Its eigenstructure has following form:


It is worth noting that the vector number 15 is the ground state vector |*G***<sup>k</sup>** of the system for momentum **k** and *E***<sup>k</sup>** = *ξ*2 **<sup>k</sup>** + <sup>Δ</sup><sup>2</sup> **<sup>k</sup>**, *<sup>u</sup>*<sup>2</sup> **<sup>k</sup>** <sup>=</sup> <sup>1</sup> 2 1 + *<sup>ξ</sup>***<sup>k</sup>** *E***k** and *v*<sup>2</sup> **<sup>k</sup>** <sup>=</sup> <sup>1</sup> 2 <sup>1</sup> <sup>−</sup> *<sup>ξ</sup>***<sup>k</sup>** *E***k** .

By making use of definition (4) we obtain the gap equation, namely

$$
\Delta\_{\mathbf{gk}} = N^{-1} \frac{1}{2} \sum\_{\mathbf{k'} > 0} g\_{\mathbf{k} \mathbf{k'}} \langle G | \beta\_{\mathbf{k'}} | G \rangle = N^{-1} \sum\_{\mathbf{k'}} g\_{\mathbf{k} \mathbf{k'}} \frac{\Delta\_{\mathbf{g} \mathbf{k'}}}{2E\_{\mathbf{k'}}}. \tag{24}$$

The chemical potential *μ* can be found from the expression for the average number of particles per lattice site

$$n = N^{-1} \frac{1}{2} \sum\_{\mathbf{k} > 0} \left\langle G \left| (n\_{\mathbf{k}+} + n\_{\mathbf{k}-} + n\_{-\mathbf{k}+} + n\_{-\mathbf{k}-}) \right| G \right\rangle = \frac{2}{N} \sum\_{\mathbf{k}} v\_{\mathbf{k}}^2 \tag{25}$$

where |*G* = ' **k**>0 |*G***<sup>k</sup>** is the total ground state vector. The ground state energy per lattice site reads

$$\frac{E\_G}{N} = \frac{1}{N} \sum\_{\mathbf{k}} (\xi\_\mathbf{k} - E\_\mathbf{k} + \frac{\Lambda\_\mathbf{k}^2}{2E\_\mathbf{k}}).\tag{26}$$

The above equations will be used throughout this paper.

#### **5. Interaction with External Electromagnetic Field at Zero Temperature**

This section is to large extent based on perturbation theory used in the book by Abrikosov [23]. Our objective is to demonstrate the existence of the Meissner–Ochsenfeld effect in the system described by the Hamiltonian *Hf* = *Hg* + *H* , where *H* represents the perturbation due to a weak static external electromagnetic field described by the vector potential **A**(**r**). Thus

$$H' = \frac{1}{2m} \int \mathbf{d}^3 r \psi^\*(\mathbf{r}) \left[ (-i\hbar \nabla - e\mathbf{A}/c)^2 - (-i\hbar \nabla)^2 \right] \psi(\mathbf{r}) .$$

After making an expansion of the field operators *ψ*∗(**r**), *ψ*(**r**) one obtains

$$H' = H\_1' + H\_2'$$

with

$$H\_1' = -\frac{e\hbar}{2mc} \sum\_{\mathbf{k}, \mathbf{q}, r} \mathbf{a}(\mathbf{q}) \cdot (\mathbf{q} + 2\mathbf{k}) a\_{\mathbf{k} + \mathbf{q}r}^\* a\_{\mathbf{k}r} \tag{27}$$

and

$$H\_2' = \frac{c^2}{2mc^2} \sum\_{\mathbf{k}, \mathbf{q}, \mathbf{q}', \sigma} \mathbf{a}(\mathbf{q}) \mathbf{a}(\mathbf{q}') a\_{\mathbf{k}\sigma}^\* a\_{\mathbf{k}-\mathbf{q}-\mathbf{q}'\sigma} \tag{28}$$

where

$$\mathbf{a}(\mathbf{k'} - \mathbf{k}) = 1/N \int d^3 r \mathbf{A}(\mathbf{r}) \exp\left[-i(\mathbf{k'} - \mathbf{k}) \cdot \mathbf{r}\right].$$

Now we would like to find corrections to the ground state energy of the system *EG* up to second order in **a**(**q**). *EG* is connected with the eigenenergy 2*ξ***<sup>k</sup>** − 2*E***k**, where *E***<sup>k</sup>** = *ξ*2 **<sup>k</sup>** + <sup>Δ</sup><sup>2</sup> **<sup>k</sup>** . To this end, let us replace **q** with **k** − **k** in *H* <sup>1</sup>. In terms of new creation and annihilation operators *a*<sup>∗</sup> **k***i* , *a***k***<sup>i</sup>* with *i* = 1, 2, 3, 4 (defined as *a*∗ **<sup>k</sup>**<sup>1</sup> := *a*<sup>∗</sup> **<sup>k</sup>**+, *a*<sup>∗</sup> **<sup>k</sup>**<sup>2</sup> := *a*<sup>∗</sup> **<sup>k</sup>**−, *<sup>a</sup>*<sup>∗</sup> **<sup>k</sup>**<sup>3</sup> := *a*<sup>∗</sup> −**k**+, *<sup>a</sup>*<sup>∗</sup> **<sup>k</sup>**<sup>4</sup> := *a*<sup>∗</sup> −**k**−) *<sup>H</sup>* <sup>1</sup> assumes the form

$$\begin{split} H\_{1}^{\prime} &= -\frac{c\hbar}{2mc} \sum\_{\begin{subarray}{c} \mathbf{k} > 0 \\ \mathbf{k}^{\prime} > 0 \end{subarray}} \Big\{ \mathbf{a} (\mathbf{k}^{\prime} - \mathbf{k}, t) \cdot (\mathbf{k}^{\prime} + \mathbf{k}) (a\_{\mathbf{k}^{\prime} \mathbf{k}}^{\*} a\_{\mathbf{k} 1} - a\_{\mathbf{k} 4}^{\*} a\_{\mathbf{k}^{\prime} 4}) + \cdots \\ &+ \mathbf{a} (-\mathbf{k}^{\prime} - \mathbf{k}, t) \cdot (-\mathbf{k}^{\prime} + \mathbf{k}) (a\_{\mathbf{k}^{\prime} \mathbf{} 3}^{\*} a\_{\mathbf{k} 1} - a\_{\mathbf{k} 4}^{\*} a\_{\mathbf{k}^{\prime} 2}) + \\ &+ \mathbf{a} (\mathbf{k}^{\prime} + \mathbf{k}, t) \cdot (\mathbf{k}^{\prime} - \mathbf{k}) (a\_{\mathbf{k}^{\prime} \mathbf{} 1}^{\*} a\_{\mathbf{k} 3} - a\_{\mathbf{k} 2}^{\*} a\_{\mathbf{k}^{\prime} 4}) + \\ &+ \mathbf{a} (-\mathbf{k}^{\prime} + \mathbf{k}, t) \cdot (-\mathbf{k}^{\prime} - \mathbf{k}) (a\_{\mathbf{k}^{\prime} 3}^{\*} a\_{\mathbf{k} 3} - a\_{\mathbf{k} 2}^{\*} a\_{\mathbf{k}^{\prime} 2}) \Big\}, \end{split}$$

where {**<sup>k</sup>** : **<sup>k</sup>** <sup>&</sup>gt; <sup>0</sup>} denotes the set of all 1-fermion momenta restricted to a definite half-space of *<sup>R</sup>*3. The expression above is a result of a rearrangement of terms entering (27) and switching **k** and **k** in some terms. Thus, we are able to calculate the first order corrections, viz.,

$$E\_1^{(1)} = \langle G | H\_1' | G \rangle = \begin{cases} 0, & \mathbf{k}' = \mathbf{k};\\ 0, & \mathbf{k}' \neq \mathbf{k}, \end{cases} \tag{29}$$

and

$$E\_2^{(1)} = \langle G \vert H\_2' \vert G \rangle = \frac{\varepsilon^2}{2mc^2} \sum\_{\mathbf{k} \neq \mathbf{q}' \vert \sigma} \mathbf{a}(\mathbf{q}) \mathbf{a}(\mathbf{q}') \langle G \vert a\_{\mathbf{k}\sigma}^\* a\_{\mathbf{k} - \mathbf{q} - \mathbf{q}' \vert \sigma} \vert G \rangle = \frac{\varepsilon^2}{2mc^2} \sum\_{\mathbf{k} \neq \mathbf{q}} 2v\_\mathbf{k}^2 \mathbf{a}(\mathbf{q}) \mathbf{a}(-\mathbf{q}). \tag{30}$$

In (30) operators are in the original form, it is, with momentum indices belonging to whole momentum space. Now the calculation of the second order term *E*(2) <sup>1</sup> can be performed, viz.,

$$E\_1^{(2)} = \sum\_{m \neq G} \frac{|\langle G | H\_1' | m \rangle|^2}{E\_G - E\_m}.$$

where |*m* and *Em* denote vectors representing excited states constructed from the vectors placed in Section 4 and their energies, respectively. *EG* is the ground state energy. To get this correction the eigenstructure from Section 4 for momenta **k** and **k** must be exploited. It is found that the only nonzero contributions are as follows:


where *u*<sup>2</sup> **<sup>k</sup>** <sup>=</sup> <sup>1</sup> <sup>2</sup> (<sup>1</sup> <sup>+</sup> *<sup>ξ</sup>***<sup>k</sup>** *<sup>E</sup>***<sup>k</sup>** ) and *<sup>v</sup>*<sup>2</sup> **<sup>k</sup>** <sup>=</sup> <sup>1</sup> <sup>2</sup> (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>***<sup>k</sup>** *<sup>E</sup>***<sup>k</sup>** ). After summing all nonzero contributions we obtain

$$\begin{split} \mathbf{E}\_{1}^{(2)} &= \begin{pmatrix} \frac{\varrho \mathbf{h}}{2mc} \end{pmatrix}^{2} \sum\_{\mathbf{k} \mathbf{q}} \begin{pmatrix} \frac{(u\_{\mathbf{k}+\mathbf{q}} v\_{\mathbf{k}})^{2}}{\xi\_{\mathbf{k}+\mathbf{q}} - \frac{\mathbf{r}}{\mathbf{q}} - 2(E\_{\mathbf{k}} + E\_{\mathbf{k}+\mathbf{q}})} + \\ \frac{(v\_{\mathbf{k}+\mathbf{q}} u\_{\mathbf{k}})^{2}}{\xi\_{\mathbf{k}} - \frac{\mathbf{r}}{\mathbf{q}} - 2(E\_{\mathbf{k}} + E\_{\mathbf{k}+\mathbf{q}})} \end{pmatrix} \mathbf{[a(q) \cdot (2\mathbf{k} + \mathbf{q})]} \mathbf{[a(-q) \cdot (2\mathbf{k} + \mathbf{q})]}. \end{split} \tag{31}$$

We exploited **a**∗(**q**) = **a**(−**q**) here. Furthermore the sum ∑ **k**>0 **k** >0 was replaced with the unrestricted .

sum ∑ **kq**

It is well known that the Hamiltonian of a system in an external magnetic field fulfils the following relation:

$$
\delta H = -\frac{1}{c} \int \mathbf{j} \delta \mathbf{A} d^3 r\_{\nu}
$$

where **j** is the current density operator. Averaging this relation over a given state, we obtain

$$
\delta E = -\frac{1}{c} \int \langle \mathbf{j} \rangle \delta \mathbf{A} d^3 r. \tag{32}
$$

Next we obtain expression (32) in Fourier representation, viz.,

$$
\delta E = -\frac{N}{c} \sum\_{\mathbf{q}} \mathbf{j}\_{\mathbf{q}} \delta \mathbf{a} (-\mathbf{q})\_{\mathbf{q}}
$$

with

$$\mathbf{j\_{q}} = \frac{1}{N} \int d^{3}r \mathbf{j} e^{i\mathbf{q}\cdot\mathbf{r}}.$$

In order to get the expression for **jq** the corrections to the ground state energy are differentiated with respect to the vector potential. It yields

$$\begin{split} \delta E &= -2\left[\frac{e^{2}}{mc^{2}} \sum\_{\mathbf{k}\mathbf{q}} v\_{\mathbf{k}}^{2} \mathbf{a}(\mathbf{q}) + \left(\frac{e\hbar}{2mc}\right)^{2} \sum\_{\mathbf{k}\mathbf{q}} \left(\frac{(u\_{\mathbf{k}+\mathbf{q}}v\_{\mathbf{k}})^{2}}{\xi\_{\mathbf{k}+\mathbf{q}}-\xi\_{\mathbf{k}}-2(E\_{\mathbf{k}}+E\_{\mathbf{k}+\mathbf{q}})} + \\ &+ \frac{(v\_{\mathbf{k}+\mathbf{q}}u\_{\mathbf{k}})^{2}}{\xi\_{\mathbf{k}}-\xi\_{\mathbf{k}+\mathbf{q}}-2(E\_{\mathbf{k}}+E\_{\mathbf{k}+\mathbf{q}})} \right) [\mathbf{a}(\mathbf{q})\cdot(2\mathbf{k}+\mathbf{q})](2\mathbf{k}+\mathbf{q}) \right] \delta \mathbf{a}(-\mathbf{q}). \end{split} \tag{33}$$

Now having the expression above for *δE* at disposal the current density can be obtained,

$$\begin{split} \mathbf{j\_{q}} &= \quad -2\frac{\varepsilon^{2}}{m c N} \sum\_{\mathbf{k}} v\_{\mathbf{k}}^{2} \mathbf{a}(\mathbf{q}) - \frac{2}{c N} \left(\frac{\varepsilon \hbar}{2m}\right)^{2} \sum\_{\mathbf{k}} \left(\frac{(u\_{\mathbf{k}+\mathbf{q}} v\_{\mathbf{k}})^{2}}{\xi\_{\mathbf{k}+\mathbf{q}} - \xi\_{\mathbf{k}} - 2(E\_{\mathbf{k}} + E\_{\mathbf{k}+\mathbf{q}})} + \\ & \quad + \frac{(v\_{\mathbf{k}+\mathbf{q}} u\_{\mathbf{k}})^{2}}{\xi\_{\mathbf{k}} - \xi\_{\mathbf{k}+\mathbf{q}} - 2(E\_{\mathbf{k}} + E\_{\mathbf{k}+\mathbf{q}})} \right) [\mathbf{a}(\mathbf{q}) \cdot (2\mathbf{k} + \mathbf{q})] (2\mathbf{k} + \mathbf{q}). \end{split} \tag{34}$$

For **q** tending to zero the relation for current density takes the following form

$$\mathbf{j}\_0 = -\mathbf{K}(0,0)\mathbf{a}(0),\tag{35}$$

with

$$K(0,0) = 2\frac{e^2}{mcN} \sum\_{\mathbf{k}} v\_{\mathbf{k}}^2 - \frac{8}{d} \frac{1}{cN} \left(\frac{e\hbar}{2m}\right)^2 \sum\_{\mathbf{k}} k^2 \frac{2(\mu\_{\mathbf{k}} v\_{\mathbf{k}})^2}{4E\_{\mathbf{k}}}.$$

where *d* is equal to the dimension of the system. The Meissner effect holds if *K*(0, 0) = *const* > 0. This is known as the Schafroth's criterion. Two zeros that stand as the arguments in *K*(0, 0) correspond to static external field (*ω* = 0) and *q* = 0. Finally, by making use of Equation (25) one obtains

$$K(0,0) = \frac{e^2}{mc}n - \frac{8}{d}\frac{1}{cN}\left(\frac{e\hbar}{2m}\right)^2 \sum\_{\mathbf{k}} k^2 \frac{2(\mu\_{\mathbf{k}}v\_{\mathbf{k}})^2}{4E\_{\mathbf{k}}}.\tag{36}$$

The Equation (36) is the basis for the investigation of the response of the system under the study to a weak external magnetic field. We shall make it for three different symmetries of the order parameter: the pure *<sup>s</sup>*-wave pairing, the *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing and the mixed *<sup>s</sup>* <sup>+</sup> *idx*<sup>2</sup>+*y*<sup>2</sup> -wave pairing cases. The *p*-wave case is not relevant to the problem of four-fermion potential in the form used in this paper because this type of pairing refers to the so-called equal-spin pairing but fermion quartets with the same |**k**| cannot manifest the pairing of this kind of symmetry. In a more general case, when the pairs and the quartets are present in the system, one can consider combinations of different types of pairing including the *p*-wave pairing as well.

#### **6. The Scenario with the Rectangular Density of States**

At first, let us consider our system with the rectangular DOS, defined as:

$$\rho(\xi) = \begin{cases} \frac{1}{\Box \xi'} & \xi \in (-\mu, D\varepsilon - \mu); \\ 0, & \text{otherwise}, \end{cases} \tag{37}$$

where *De* is the bandwidth. Let us begin with the order parameter with the pure *s*-wave pairing symmetry. In this case the coupling function *g***kk** = *g* > 0 and as a consequence Δ*g***<sup>k</sup>** = Δ*g*. In the thermodynamic limit the sums in Equations (24)–(26) and (36) turn into integrals, i.e.,

$$1 = \frac{\mathcal{g}}{2D\varepsilon} \int\_{-\mu}^{D\varepsilon-\mu} \frac{d\mathfrak{g}}{E} = \frac{\mathcal{g}}{2D\varepsilon} \left( \mathrm{arc}\,\sinh\frac{D\varepsilon-\mu}{\Delta\_{\mathfrak{g}}} + \mathrm{arc}\,\sinh\frac{\mu}{\Delta\_{\mathfrak{g}}} \right) \tag{38}$$

$$n = \frac{1}{D\varepsilon} \int\_{-\mu}^{D\varepsilon-\mu} \mathrm{d}\xi \left(1 - \frac{\overline{\xi}}{E}\right) = 1 - \frac{1}{D\varepsilon} \left(\sqrt{(D\varepsilon-\mu)^2 + \Delta\_{\overline{\xi}}^2} - \sqrt{\mu^2 + \Delta\_{\overline{\xi}}^2}\right),\tag{39}$$

$$\begin{split} \frac{^{E\_{\mathcal{C}}}}{^{N}} &= \frac{1}{D\varepsilon} \int\_{-\mu}^{D\varepsilon-\mu} \mathrm{d}\xi \left( \xi - E + \frac{\Delta\_{\mathcal{S}}^{2}}{2E} \right) \\ &= \frac{1}{2} (D\varepsilon - 2\mu) - \frac{1}{2D\varepsilon} \left( (D\varepsilon - \mu) \sqrt{(D\varepsilon - \mu)^{2} + \Delta\_{\mathcal{S}}^{2}} + \mu \sqrt{\mu^{2} + \Delta\_{\mathcal{S}}^{2}} \right) - \frac{\Delta\_{\mathcal{S}}^{2}}{\mathcal{S}} \end{split} \tag{40}$$

and finally

$$\begin{split} K(0,0) &= \frac{\varepsilon^{2}}{m\varepsilon} \Big[ n - \frac{1}{D\varepsilon} \int\_{-\mu}^{D\varepsilon-\mu} \mathrm{d}\_{\varepsilon}^{x} (\zeta + \mu) \frac{\Delta\_{\mathfrak{x}}^{2}}{4\mathcal{E}^{3}} \Big] \\ &= \frac{\varepsilon^{2}}{m\varepsilon} \Big[ n - \frac{\Lambda\_{\mathfrak{x}}^{2}}{4D\varepsilon} \Bigg[ \frac{1}{\sqrt{\mu^{2} + \Lambda\_{\mathfrak{x}}^{2}}} - \frac{1}{\sqrt{(D\varepsilon-\mu)^{2} + \Lambda\_{\mathfrak{x}}^{2}}} + \frac{\mu}{\Lambda\_{\mathfrak{x}}^{2}} \left( \frac{D\varepsilon-\mu}{\sqrt{(D\varepsilon-\mu)^{2} + \Lambda\_{\mathfrak{x}}^{2}}} + \frac{\mu}{\sqrt{\mu^{2} + \Lambda\_{\mathfrak{x}}^{2}}} \right) \Bigg] \end{split} \tag{41}$$

where *E* = *ξ*<sup>2</sup> + Δ<sup>2</sup> *<sup>g</sup>*. The one-electron dispersion relation counted with respect to *μ* is as follows *<sup>ξ</sup>***<sup>k</sup>** <sup>=</sup> *<sup>h</sup>*¯ <sup>2</sup>**k**<sup>2</sup> <sup>2</sup>*<sup>m</sup>* − *μ* while the system is two-dimensional. The Equations (38)–(40) are the same as in the BCS model with the rectangular DOS; However, the Equation (41) differs from the BCS counterpart that has the following form:

$$K(0,0) = \frac{e^2 n}{mc}.$$

As one can notice, at the zero temperature, there is only the so-called diamagnetic part in the electromagnetic kernel. At higher temperatures, the paramagnetic term appears and increases up to the transition temperature *Tc*. At this temperature, both terms cancel each other, which means the transition of a superconducting system goes back to the normal state. This happens so independently of the symmetry of pairing. The situation is significantly different if we replace the BCS interaction with the four-fermion attraction. By the appearance of the paramagnetic term, the Meissner effect gets weaker even at zero temperature. This was discovered in 2004 [24]. In that paper, the electromagnetic kernel at zero temperature reads

$$K(0,0) = \frac{3}{4} \frac{e^2 n}{mc}.\tag{42}$$

Now, for simplicity sake, let us consider the half-filled band case, i.e., *n* = 1. Then the Equations (38)–(41) yield *μ* = *De* <sup>2</sup> and

$$\frac{E\_G}{N} = -\frac{1}{2}\sqrt{\frac{D\varepsilon^2}{4} + \Lambda\_\mathcal{S}^2} - \frac{\Lambda\_\mathcal{S}^2}{\mathcal{S}}, \quad \Delta\_\mathcal{S} = \frac{D\varepsilon}{\sinh\frac{D\varepsilon}{\mathcal{S}}}, \quad K(0,0) = \frac{\varepsilon^2}{mc} \left[1 - \frac{1}{8}\frac{D\varepsilon}{\sqrt{\frac{D\varepsilon^2}{4} + \Lambda\_\mathcal{S}^2}}\right], \tag{43}$$

Note that in the limit of the shrinking band *De* → 0 (the atomic limit) one obtains the BCS result *K*(0, 0) = *<sup>e</sup>*<sup>2</sup> *mc* . In the opposite regime <sup>Δ</sup>*<sup>g</sup> De* <sup>2</sup> we obtain the result (42). This means that the external magnetic field penetrates a superconductor deeper in the case of the system with the four-fermion attraction and the wider band than in the case of the system with the BCS potential or than with the four-fermion potential but with a very narrow band.

Now, let us deal with the scenario where *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing symmetry. Let us take the coupling function and the order parameter in the form *g***kk** = *g*(*k*<sup>2</sup> *<sup>x</sup>* <sup>−</sup> *<sup>k</sup>*<sup>2</sup> *y*)(*k x* <sup>2</sup> <sup>−</sup> *<sup>k</sup> y* 2 ) and Δ*g***<sup>k</sup>** = Δ*g*(*k*<sup>2</sup> *x* − *k*2 *<sup>y</sup>*), respectively. In angle representation one gets *g*(*φ*, *φ* ) = *g* cos 2*φ* cos 2*φ* and Δ*g*(*φ*) = Δ*<sup>g</sup>* cos 2*φ*. The Equations (24)–(26) and (36) take the following form:

$$1 = \frac{\mathcal{S}}{4\pi D\varepsilon} \int\_0^{2\pi} \mathbf{d}\phi \cos^2 2\phi \left( \text{arc}\,\sinh\frac{D\varepsilon - \mu}{\Delta\_{\mathcal{S}} |\cos 2\phi|} + \text{arc}\,\sinh\frac{\mu}{\Delta\_{\mathcal{S}} |\cos 2\phi|} \right),\tag{44}$$

$$m = 1 - \frac{1}{2\pi D \varepsilon} \int\_0^{2\pi} \mathrm{d}\phi \left( \sqrt{(De - \mu)^2 + \Delta\_\mathcal{S}^2 \cos^2 2\phi} - \sqrt{\mu^2 + \Delta\_\mathcal{S}^2 \cos^2 2\phi} \right), \tag{45}$$

$$\frac{E\_{\mathcal{C}}}{N} = \frac{1}{2}(D\varepsilon - 2\mu) - \frac{1}{4\pi D\varepsilon} \int\_{0}^{2\pi} \mathrm{d}\phi \left( (D\varepsilon - \mu)\sqrt{(D\varepsilon - \mu)^{2} + \Lambda\_{\mathcal{S}}^{2}\cos^{2}2\phi} + \mu\sqrt{\mu^{2} + \Lambda\_{\mathcal{S}}^{2}\cos^{2}2\phi} \right) - \frac{\Lambda\_{\mathcal{S}}^{2}}{\mathcal{S}} \tag{46}$$

and finally

$$\begin{split} K(0,0) &= \frac{\varepsilon^2}{mc} \Big[ n - \frac{\Lambda\_{\mathcal{S}}^2}{8\pi D\varepsilon} \Big[ \int\_0^{2\pi} \mathrm{d}\phi \cos^2 2\phi \left( \frac{1}{\sqrt{\mu^2 + \Lambda\_{\mathcal{S}}^2 \cos^2 2\phi}} - \frac{1}{\sqrt{(D\varepsilon - \mu)^2 + \Lambda\_{\mathcal{S}}^2 \cos^2 2\phi}} \right) + \\ &+ \frac{\mu}{\Lambda\_{\mathcal{S}}^2} \Big[ \mathrm{d}\phi \Big( \frac{D\varepsilon - \mu}{\sqrt{(D\varepsilon - \mu)^2 + \Lambda\_{\mathcal{S}}^2 \cos^2 2\phi}} + \frac{\mu}{\sqrt{\mu^2 + \Lambda\_{\mathcal{S}}^2 \cos^2 2\phi}} \Big) \Big] \Big]. \end{split} \tag{47}$$

In the case of the half-filled band (*n* = 1) the Equations (44)–(47) take a simpler form, namely:

$$\mu = \frac{D\varepsilon}{2}, \qquad 1 = \frac{\mathcal{S}}{2\pi D\varepsilon} \int\_0^{2\pi} \mathbf{d}\phi \cos^2 2\phi \text{arc } \sinh \frac{D\varepsilon}{2\Lambda\_{\mathcal{S}} |\cos 2\phi|} \tag{48}$$

$$\frac{E\_G}{N} = -\frac{1}{4\pi} \int\_0^{2\pi} \mathrm{d}\phi \sqrt{\frac{D\epsilon^2}{4} + \Lambda\_\mathcal{S}^2 \cos^2 2\phi} - \frac{\Lambda\_\mathcal{S}^2}{\mathcal{S}} \tag{49}$$

and

$$K(0,0) = \frac{e^2}{mc} \left[ 1 - \frac{De}{16\pi} \int\_0^{2\pi} \frac{d\phi}{\sqrt{\frac{Dc^2}{4} + \Lambda\_\beta^2 \cos^2 2\phi}} \right]. \tag{50}$$

The integrals in the above equations can be expressed by the elliptic integrals which can be calculated numerically what is made later in the work. One can try to obtain some analytical results, namely, assume that <sup>Δ</sup>*<sup>g</sup> De* <sup>2</sup> then one can apply the following expansions:

$$\sqrt{\frac{D\varepsilon^2}{4} + \Lambda\_\mathcal{\mathcal{S}}^2 \cos^2 2\phi} = \frac{D\varepsilon}{2} + \frac{\Lambda\_\mathcal{\mathcal{S}}^2 \cos^2 2\phi}{D\varepsilon} + \dots$$

$$\frac{1}{\sqrt{\frac{D\varepsilon^2}{4} + \Lambda\_\mathcal{\mathcal{S}}^2 \cos^2 2\phi}} = \frac{2}{D\varepsilon} - \frac{4}{D\varepsilon^3} \Lambda\_\mathcal{\mathcal{S}}^2 \cos^2 2\phi + \dots$$

If we limit ourselves to the first order terms and substitute them into the Equations (49) and (50) we will obtain after integration:

$$\frac{E\_G}{N} \approx -\frac{1}{4} D \varepsilon - \frac{1}{4} \frac{\Delta\_\mathcal{S}^2}{Dc} - \frac{\Delta\_\mathcal{S}^2}{g} \tag{51}$$

and

$$K(0,0) \approx \frac{e^2}{mc} \left[\frac{3}{4} + \frac{1}{4} \frac{\Lambda\_\chi^2}{Dc^2}\right].\tag{52}$$

Notice that the neglecting of the second term in (52) gives the same result as (42). At last let us deal with the mixed *<sup>s</sup>* <sup>+</sup> *idx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing symmetry. We take the coupling function in the form *g***kk** = *g*<sup>0</sup> + *g*1(*k*<sup>2</sup> *<sup>x</sup>* <sup>−</sup> *<sup>k</sup>*<sup>2</sup> *y*)(*k x* <sup>2</sup> <sup>−</sup> *<sup>k</sup> y* 2 ), where *g*<sup>0</sup> and *g*<sup>1</sup> are positive real numbers. The gap parameter is complex and reads Δ*g***<sup>k</sup>** = Δ*<sup>s</sup>* + *i*Δ*d*(*kx* <sup>2</sup> <sup>−</sup> *ky* <sup>2</sup>) and the equation for it at zero temperature is as follows

$$
\Delta\_{\text{gk}} = \frac{1}{N} \sum\_{\mathbf{k'}} \frac{\Delta\_{\text{gk'}}}{2E\_{\mathbf{k'}}} \lambda
$$

where *E***<sup>k</sup>** = *ξ***<sup>k</sup>** + Δ<sup>2</sup> *<sup>s</sup>* + Δ<sup>2</sup> *<sup>d</sup>*(*kx* <sup>2</sup> <sup>−</sup> *ky* <sup>2</sup>)2. This equation splits up to two coupled equations for Δ*<sup>s</sup>* and Δ*d*, namely,

$$1 = \frac{\mathcal{G}0}{N} \sum\_{\mathbf{k}} \frac{1}{2E\_{\mathbf{k}}}, \quad 1 = \frac{\mathcal{G}1}{N} \sum\_{\mathbf{k}} \frac{(k\_{\mathbf{x}}^{\ \ 2} - k\_{\mathbf{y}}^{\ \ 2})^2}{2E\_{\mathbf{k}}}.$$

In the angle representation the factor *k*<sup>2</sup> *<sup>x</sup>* <sup>−</sup> *<sup>k</sup>*<sup>2</sup> *<sup>y</sup>* is replaced by cos 2*φ*. In this representation our problem can be expressed by the set of equations:

$$1 = \frac{\mathcal{G}^0}{4\pi D\varepsilon} \int\_0^{2\pi} \mathbf{d}\phi \left(\text{arc}\,\sinh\frac{D\varepsilon-\mu}{\sqrt{\Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi}} + \text{arc}\,\sinh\frac{\mu}{\sqrt{\Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi}}\right),\tag{53}$$

$$1 = \frac{\mathcal{G}\_1}{4\pi De} \int\_0^{2\pi} d\phi \cos^2 2\phi \left( \text{arc}\,\sinh\frac{D\varepsilon-\mu}{\sqrt{\Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi}} + \text{arc}\,\sinh\frac{\mu}{\sqrt{\Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi}} \right), \tag{54}$$

$$m = 1 - \frac{1}{2\pi D \varepsilon} \int\_0^{2\pi} \mathbf{d}\phi \Big( \sqrt{(D\varepsilon - \mu)^2 + \Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi} - \sqrt{\mu^2 + \Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi} \Big), \tag{55}$$

$$\frac{E\_G}{N} = \frac{1}{2}(De - 2\mu) - \frac{\Delta\_d^2}{\mathcal{g}\_1} - \frac{\Delta\_s^2}{\mathcal{g}\_0} -$$

$$- \frac{1}{4\pi D\varepsilon} \int\_0^{2\pi} \mathrm{d}\phi \left( (De - \mu)\sqrt{(De - \mu)^2 + \Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi} + \mu\sqrt{\mu^2 + \Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi} \right) \tag{56}$$

and

$$\begin{cases} K(0,0) = \quad \frac{\varepsilon^2}{mc} \Big[ n - \frac{1}{8\pi Dc} \Bigg[ \frac{2\pi}{\ell} \operatorname{d}\phi (\Lambda\_s^2 + \Lambda\_d^2 \cos^2 2\phi) \left( \frac{1}{\sqrt{\mu^2 + \Lambda\_s^2 + \Lambda\_d^2 \cos^2 2\phi}} - \frac{1}{\sqrt{(D\epsilon - \mu)^2 + \Lambda\_s^2 + \Lambda\_d^2 \cos^2 2\phi}} \right) + \operatorname{d}\phi (\Lambda\_s^2 - \mu) \frac{1}{\sqrt{(D\epsilon - \mu)^2 + \Lambda\_s^2 + \Lambda\_d^2 \cos^2 2\phi}} \Bigg] \tag{57} \\ \qquad + \operatorname{\mu} \int\_0^T \operatorname{d}\phi \left( \frac{D\epsilon - \mu}{\sqrt{(D\epsilon - \mu)^2 + \Lambda\_s^2 + \Lambda\_d^2 \cos^2 2\phi}} + \frac{\mu}{\sqrt{\mu^2 + \Lambda\_s^2 + \Lambda\_d^2 \cos^2 2\phi}} \right) \Bigg]. \end{cases} \tag{57}$$

Because we are interested in the half-filled band case the above equations take a much simpler form, i.e., *μ* = *De* <sup>2</sup> and

$$1 = \frac{\mathcal{G}\_0}{2\pi D\varepsilon} \int\_0^{2\pi} \text{arc}\,\text{sinh}\frac{D\varepsilon}{2\sqrt{\Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi}} \text{d}\phi,\tag{58}$$

$$\mathcal{L} = \frac{\mathcal{S}\_1}{2\pi D\varepsilon} \int\_0^{2\pi} \cos^2 2\phi \text{arc}\,\sinh\frac{D\varepsilon}{2\sqrt{\Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi}} \text{d}\phi,\tag{59}$$

$$\frac{E\_G}{L} = -\frac{1}{4\pi} \int\_0^{2\pi} \mathbf{d}\phi \sqrt{\frac{D\epsilon^2}{4} + \Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi} - \frac{\Delta\_s^2}{\mathcal{S}\_0} - \frac{\Delta\_d^2}{\mathcal{S}\_1} \tag{60}$$

and finally

$$K(0,0) = \frac{c^2}{mc} \left[ 1 - \frac{Dc}{16\pi} \int\_0^{2\pi} \mathrm{d}\phi \frac{1}{\sqrt{\frac{Dc^2}{4} + \Delta\_s^2 + \Delta\_d^2 \cos^2 2\phi}} \right].\tag{61}$$

Once again, let us try to find some analytical expressions for *EG <sup>N</sup>* and *K*(0, 0). Note that in the limit of *De* <sup>→</sup> 0 one obtains the BCS result *<sup>K</sup>*(0, 0) = *<sup>e</sup>*<sup>2</sup> *mc* . Therefore, the narrower the band is the more similar to the BCS superconductor the system is. Reversely, for the sufficiently wide band, i.e., *De* <sup>2</sup> , Δ*<sup>s</sup>* Δ*<sup>d</sup>* one has the following expansions:

$$
\sqrt{\frac{Dc^2}{4} + \Lambda\_s^2 + \Lambda\_d^2 \cos^2 2\phi} = \sqrt{\frac{Dc^2}{4} + \Lambda\_s^2} + \frac{1}{2} \frac{\Lambda\_d^2 \cos^2 2\phi}{\sqrt{\frac{Dc^2}{4} + \Lambda\_s^2}} + \dots
$$

$$
\frac{1}{\sqrt{\frac{Dc^2}{4} + \Lambda\_s^2 + \Lambda\_d^2 \cos^2 2\phi}} = \frac{1}{\sqrt{\frac{Dc^2}{4} + \Lambda\_s^2}} - \frac{1}{2} \frac{1}{(\frac{Dc^2}{4} + \Lambda\_s^2)^{\frac{1}{2}}} \Delta\_d^2 \cos^2 2\phi + \dots
$$

It suffices to limit ourselves to the zero and first order terms and substitute them into Equations (60) and (61). After integrating we obtain

$$\frac{E\_G}{N} \approx -\frac{1}{2}\sqrt{\frac{Dc^2}{4} + \Delta\_s^2} - \frac{1}{8}\frac{\Delta\_d^2}{\sqrt{\frac{Dc^2}{4} + \Delta\_s^2}} - \frac{\Delta\_d^2}{\mathcal{S}^1} - \frac{\Delta\_s^2}{\mathcal{S}^0} \tag{62}$$

and finally

$$K(0,0) \approx \frac{c^2}{mc} \left[ 1 - \frac{1}{8} \frac{De}{\sqrt{\frac{Dc^2}{4} + \Lambda\_s^2}} \left( 1 - \frac{1}{4} \frac{\Lambda\_d^2}{\frac{Dc^2}{4} + \Lambda\_s^2} \right) \right].\tag{63}$$

The last equation can be simplified if we neglect the term proportional to Δ<sup>2</sup> *<sup>d</sup>* because it is much lesser than the unity. Therefore, this yields

$$K(0,0) \approx \frac{e^2}{mc} \left[ 1 - \frac{1}{8} \frac{De}{\sqrt{\frac{Dc^2}{4} + \Delta\_s^2}} \right]. \tag{64}$$

The kernel (64) fulfills the double inequality

$$
\frac{3}{4} \frac{e^2}{mc} < K(0,0) < \frac{e^2}{mc}.\tag{65}
$$

We have analyzed the situation in which the component with the *s*-wave pairing symmetry and the bandwidth are much greater than the *dx*<sup>2</sup>−*y*<sup>2</sup> -wave component. If we admit that *De* <sup>2</sup> Δ*s*, Δ*<sup>d</sup>* then the ground state energy and the electromagnetic kernel read

$$\frac{E\_G}{N} \approx -\frac{1}{4} D \varepsilon - \frac{1}{2} \frac{\Delta\_s^2}{D\varepsilon} - \frac{1}{4} \frac{\Delta\_d^2}{D\varepsilon} - \frac{\Delta\_d^2}{g\_1} - \frac{\Delta\_s^2}{g\_0} \tag{66}$$

and

$$K(0,0) \approx \frac{e^2}{mc} \left[ \frac{3}{4} + \frac{1}{2} \frac{\Delta\_s^2}{De^2} + \frac{1}{4} \frac{\Delta\_d^2}{De^2} \right]. \tag{67}$$

Looking at the Equation (67) one can notice that the inequality (65) holds.

#### **7. The Scenario with the Tight-Binding Model—The One-Dimensional Case**

In this section we shall do the same kind of calculations but instead of the less realistic model with rectangular density of states we shall consider the one-dimensional model in the frame of the tight-binding approximation. In the one-dimensional system the dispersion relation has the form *ξ***<sup>k</sup>** = −2*t* cos *k* − *μ*. As is known in one-dimensional systems due to strong quantum fluctuations there is no long-range order even at *T* = 0 but one can assume that the investigated chain interacts with environment that stabilizes the superconductivity in the chain. One can say that we have to deal with the effective model. In the one-dimensional case one can take under considerations the *s*-wave pairing

only. It is no use to investigate the *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing case because the momentum has only one component. The gap equation reads

$$1 = \frac{g}{2\pi} \int\_0^{\pi} \frac{\mathrm{d}k}{\sqrt{(2t\cos k + \mu)^2 + \Delta\_\S^2}} = \frac{g}{2\pi} \int\_{-2t - \mu}^{2t - \mu} \frac{\mathrm{d}\xi}{\sqrt{\xi^2 + \Delta\_\S^2}} \frac{1}{\sqrt{4t^2 - (\xi + \mu)^2}},\tag{68}$$

where after the introducing of the density of states this equation has been easily transformed to the second form. Note that at half filling (*μ* = 0) DOS has two the so-called Van Hove singularities at *ξ* = ±2*t*. Next, the equation for the average number of particles per lattice site is

$$m = 1 + \frac{1}{\pi} \int\_0^\pi d\mathbf{k} \frac{2t \cos k + \mu}{\sqrt{(2t \cos k + \mu)^2 + \Delta\_\mathcal{g}^2}} = 1 - \frac{1}{\pi} \int\_{-2t - \mu}^{2t - \mu} \frac{\mathbf{d}\_\mathcal{g}^\mathbf{x}}{\sqrt{4t^2 - (\mathfrak{f} + \mu)^2}} \frac{\mathfrak{f}}{\sqrt{\mathfrak{f}^2 + \Delta\_\mathcal{g}^2}}.\tag{69}$$

The ground state energy per lattice site is expressed via

$$\frac{E\_G}{N} = -\mu - \frac{1}{\pi} \int\_0^\pi \mathrm{d}k \sqrt{(2t\cos k + \mu)^2 + \Delta\_\mathcal{g}^2} + \frac{\Delta\_\mathcal{g}^2}{\mathcal{g}} = -\mu - \frac{1}{\pi} \int\_{-2t - \mu}^{2t - \mu} \frac{\sqrt{\xi^2 + \Delta\_\mathcal{g}^2}}{\sqrt{4t^2 - (\xi + \mu)^2}} \mathrm{d}\xi + \frac{\Delta\_\mathcal{g}^2}{\mathcal{g}}.\tag{70}$$

Finally, the electromagnetic kernel can be found from

$$\begin{split} K(0,0) &= \frac{e^2}{mc} \Big[ n - \frac{t\Lambda\_{\xi}^2}{2\pi} \int\_0^{\pi} \frac{k^2 \mathrm{d}k}{\left( (2t\cos k + \mu)^2 + \Lambda\_{\xi}^2 \right)^{3/2}} \Bigg] \\ &= \frac{e^2}{mc} \Big[ n - \frac{t\Lambda\_{\xi}^2}{2\pi} \int\_{-2t-\mu}^{2t-\mu} \mathrm{d}\xi \frac{\left( \arccos\left( -\frac{\frac{\pi}{\xi} + \mu}{2t} \right) \right)^2}{\sqrt{4t^2 - (\frac{\pi}{\xi} + \mu)^2}} \frac{1}{(\frac{\pi^2}{\xi^2 + \Lambda\_{\xi}^2)^{3/2}}} \Bigg], \end{split} \tag{71}$$

where *<sup>t</sup>* = *<sup>h</sup>*¯ <sup>2</sup> <sup>2</sup>*ma*<sup>2</sup> is the hopping parameter while *a* is the lattice constant that is put here equal the unity. In this section we are interested in the half-filled case (*n* = 1) as well. The above equations get

a simpler form. Let us start from the chemical potential *μ* = 0 and the gap equation, namely,

$$1 = \frac{\mathcal{S}}{\pi} \int\_0^{2t} \mathrm{d}\xi \frac{1}{\sqrt{4t^2 - \xi^2}} \frac{1}{\sqrt{\xi^2 + \Delta\_{\mathcal{S}}^2}}. \tag{72}$$

Note that in the one-dimensional case the bandwidth *De* = 4*t*. The ground state energy per lattice site is as follows

$$\frac{E\_G}{N} = -\frac{2}{\pi} \int\_0^{2t} \mathrm{d}\xi \frac{\sqrt{\xi^2 + \Delta\_{\xi}^2}}{\sqrt{4t^2 - \xi^2}} + \frac{\Delta\_{\xi}^2}{\mathcal{S}}.\tag{73}$$

Finally, the electromagnetic kernel reads

$$K(0,0) = \frac{e^2}{mc} \left[ 1 - \frac{t\Lambda\_{\frac{\pi}{\mathcal{S}}}^2}{2\pi} \int\_{-2t}^{2t} \mathrm{d}\_{\sharp}^2 \frac{\left(\arccos\left(-\frac{\xi}{2t}\right)\right)^2}{\sqrt{4t^2 - \xi^2}} \frac{1}{(\xi^2 + \Lambda\_{\frac{\pi}{\mathcal{S}}}^2)^{3/2}} \right]. \tag{74}$$

We can try to find an analytical solution of the gap Equation (72) in the *De* <sup>2</sup> Δ*<sup>g</sup>* regime (the strong coupling regime). In this case the expansion of <sup>1</sup> 4*t*<sup>2</sup> cos2 *k*+Δ<sup>2</sup> *g* <sup>≈</sup> <sup>1</sup> <sup>Δ</sup>*<sup>g</sup>* <sup>−</sup> <sup>1</sup> 2 4*t* <sup>2</sup> cos2 *k* Δ3 *g* to the first order is useful, then the Equation (72) reduces to a cubic equation that can be solved by the method from Section 3. This equation takes the form

$$
\Delta\_{\mathcal{S}}^3 - \frac{\mathcal{S}}{2} \Delta\_{\mathcal{S}}^2 + \frac{\mathcal{S}t^2}{2} = 0. \tag{75}
$$

After finding the proper solution of the Equation (75) one can make use of the similar expansions of 4*t*<sup>2</sup> cos<sup>2</sup> *k* + Δ<sup>2</sup> *<sup>g</sup>* <sup>≈</sup> <sup>Δ</sup>*<sup>g</sup>* <sup>+</sup> <sup>1</sup> 2 4*t* <sup>2</sup> cos2 *<sup>k</sup>* <sup>Δ</sup>*<sup>g</sup>* and <sup>1</sup> (4*t*<sup>2</sup> cos2 *k*+Δ<sup>2</sup> *<sup>g</sup>*)3/2 <sup>≈</sup> <sup>1</sup> Δ3 *g* <sup>−</sup> <sup>3</sup> 2 4*t* <sup>2</sup> cos<sup>2</sup> *k* Δ5 *g* and use them for the calculation of the ground state energy and the electromagnetic kernel, i.e.,

$$\frac{E\_G}{N} \approx -\Delta\_\% - \frac{t^2}{\Delta\_\%} + \frac{\Delta\_\%^2}{g} \tag{76}$$

and

$$K(0,0) \approx \frac{e^2}{mc} \left[ 1 - \frac{\pi^2}{6} \frac{t}{\Delta\_{\mathcal{g}}} + \left( 2\pi + \frac{3}{4} \right) \frac{t^3}{\Delta\_{\mathcal{g}}^3} \right]. \tag{77}$$

The last equation can be approximated by the neglecting of the third term in the square brackets, i.e.,

$$K(0,0) \approx \frac{e^2}{mc} \left[1 - \frac{\pi^2}{6} \frac{t}{\Delta\_{\%}}\right]. \tag{78}$$

It is visible, that the paramagnetic term is proportional to the hopping parameter *t*. The lesser it is the greater the kernel becomes. In the limit *t* → 0 the kernel tends to the BCS one. As one can convince oneself that analytical result (78) is in agreement with numerical calculation outcomes which can be found in Section 9.

#### **8. The Scenario with the Tight-Binding Model—The Two-Dimensional Case**

In this section we shall do the same kind of calculations but we shall consider the two-dimensional model in the frame of the tight-binding approximation. In the case of two-dimensional one we shall pay an attention on the simple square lattice with the dispersion relation counted with respect to *μ*, i.e., *ξ***<sup>k</sup>** = −2*t*(cos *kx* + cos *ky*) − *μ*, where *t* is the hopping parameter while the lattice constant is put equal to the unity. Again, we are interested in the study of the following symmetries of the gaps: the *s*-wave pairing and the mixed *<sup>s</sup>* <sup>+</sup> *idx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing ones. The pure *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing symmetry case is treated in a numerical manner. The expressions for this case can be obtained from the equations for the *<sup>s</sup>* <sup>+</sup> *idx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing symmetry case by putting <sup>Δ</sup>*<sup>s</sup>* <sup>=</sup> 0 in them.

Let us start with the gap equation for Δ*<sup>g</sup>* with the *s*-wave pairing symmetry. This reads

$$\mathbf{1} = \frac{\mathcal{g}}{2\pi^2} \int\_0^\pi \mathrm{d}k\_x \int\_0^\pi \mathrm{d}k\_y \frac{1}{\sqrt{(2t(\cos k\_x + \cos k\_y) + \mu)^2 + \Delta\_g^2}}\,\tag{79}$$

next the equation for the average number of particles per lattice site

$$m = 1 + \frac{1}{\pi^2} \int\_0^\pi \mathrm{d}k\_x \int\_0^\pi \mathrm{d}k\_y \frac{2t(\cos k\_x + \cos k\_y) + \mu}{\sqrt{(2t(\cos k\_x + \cos k\_y) + \mu)^2 + \Delta\_\mathcal{g}^2}} \tag{80}$$

and those for the ground state energy per lattice site and the electromagnetic kernel:

$$\frac{E\_G}{N} = -\mu - \frac{1}{\pi^2} \int\_0^\pi \mathrm{d}k\_x \int\_0^\pi \mathrm{d}k\_y \sqrt{(2t(\cos k\_x + \cos k\_y) + \mu)^2 + \Lambda\_\mathcal{S}^2} + \frac{\Lambda\_\mathcal{S}^2}{\mathcal{S}},\tag{81}$$

$$K(0,0) = \frac{e^2}{mc} \left[ n - \frac{t\Delta\_\text{g}^2}{4\pi^2} \int\_0^\pi \text{dk}\_x \int\_0^\pi \text{dk}\_y (k\_x^2 + k\_y^2) \frac{1}{((2t(\cos k\_x + \cos k\_y) + \mu)^2 + \Delta\_\text{g}^2)^{\frac{3}{2}}} \right].\tag{82}$$

These equations for the half-filled band case can be produced by putting *n* = 1 and *μ* = 0 and substituting these values into them. Next, in order to find some analytical results we follow the procedure used in the previous section, i.e., assume that *De* <sup>2</sup> Δ*g*, where *De* = 8*t* is the bandwidth in the two-dimensional system. We will exploit the following expansions to the first order terms:

$$\begin{aligned} \sqrt{4t^2(\cos k\_x + \cos k\_y)^2 + \Delta\_{\mathcal{S}}^2} &\approx \Delta\_{\mathcal{S}} + \frac{1}{2} \frac{4t^2(\cos k\_x + \cos k\_y)^2}{\Delta\_{\mathcal{S}}}, \\\\ \frac{1}{\sqrt{4t^2(\cos k\_x + \cos k\_y)^2 + \Delta\_{\mathcal{S}}^2}} &\approx \frac{1}{\Delta\_{\mathcal{S}}} - \frac{1}{2} \frac{4t^2(\cos k\_x + \cos k\_y)^2}{\Delta\_{\mathcal{S}}^3} \end{aligned}$$

and

$$\frac{1}{(4t^2(\cos k\_x + \cos k\_y)^2 + \Delta\_\mathcal{g}^2)^{\frac{3}{2}}} \approx \frac{1}{\Delta\_\mathcal{g}^3} - \frac{3}{2} \frac{4t^2(\cos k\_x + \cos k\_y)^2}{\Delta\_\mathcal{g}^5}$$

These expansions are now substituted into the Equations (79), (81) and (82) and after making the integration one obtains:

$$1 = \frac{\text{g}}{2} \left[ \frac{1}{\Lambda\_{\text{g}}} - 2 \frac{t^2}{\Lambda\_{\text{g}}^3} \right] \Leftrightarrow \Lambda\_{\text{g}}^3 - \frac{\text{g}}{2} \Lambda\_{\text{g}}^2 + \text{g}t^2 = 0,\tag{83}$$

$$\frac{E\_G}{N} \approx -\Delta\_{\mathcal{S}} - 2\frac{t^2}{\Delta\_{\mathcal{S}}} + \frac{\Delta\_{\mathcal{S}}^2}{\mathcal{g}} \tag{84}$$

.

and

$$K(0,0) \approx \frac{\epsilon^2}{mc} \left[ 1 - \frac{\pi^2}{6} \frac{t}{\Delta\_{\mathcal{g}}} + \frac{3}{2} \frac{t^3}{\Delta\_{\mathcal{g}}^3} \left( \frac{1}{2} + \frac{4}{3} \pi + \frac{\pi^2}{3} \right) \right].$$

If *<sup>t</sup>* <sup>Δ</sup>*<sup>g</sup>* 1 then we can approximate *K*(0, 0) by

$$K(0,0) \approx \frac{e^2}{mc} \left[1 - \frac{\pi^2}{6} \frac{t}{\Delta\_{\mathcal{S}}}\right]. \tag{85}$$

Note that (85) is the same as (78).

The pure *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing case will be investigated in another paper. We will look now at the case with mixed symmetry of pairing, namely, the *<sup>s</sup>* <sup>+</sup> *idx*<sup>2</sup>−*y*<sup>2</sup> . The order parameter has the form Δ*g***<sup>k</sup>** = Δ*<sup>s</sup>* + *i*Δ*d*(cos *kx* − cos *ky*) while its module |Δ*g***k**| = Δ2 *<sup>s</sup>* + Δ<sup>2</sup> *<sup>d</sup>*(cos *kx* − cos *ky*)2. The coupling function is *g***kk** = *g*<sup>0</sup> + *g*1(cos *kx* − cos *ky*)(cos *k <sup>x</sup>* − cos *k <sup>y</sup>*), where *g*0, *g*<sup>1</sup> are real positive numbers. Similarly, as it was in Section 6, we shall give the fundamental equations for our problem: the gap equations, the equation for the chemical potential, the ground state energy and the electromagnetic kernel. They read

$$1 = \frac{\mathcal{S}0}{2\pi^2} \int\_0^\pi \mathrm{d}k\_x \int\_0^\pi \mathrm{d}k\_y \frac{1}{\sqrt{(2t(\cos k\_x + \cos k\_y) + \mu)^2 + \Lambda\_s^2 + \Lambda\_d^2(\cos k\_x - \cos k\_y)^2}},\tag{86}$$

$$1 = \frac{\mathcal{S}1}{2\pi^2} \int\_0^\pi \mathrm{d}k\_x \int\_0^\pi \mathrm{d}k\_y \frac{(\cos k\_x - \cos k\_y)^2}{\sqrt{(2t(\cos k\_x + \cos k\_y) + \mu)^2 + \Delta\_s^2 + \Delta\_d^2 (\cos k\_x - \cos k\_y)^2}},\tag{87}$$

$$n = 1 + \frac{1}{\pi^2} \int\_0^\pi \mathrm{d}k\_x \int\_0^\pi \mathrm{d}k\_y \frac{2t(\cos k\_x + \cos k\_y) + \mu}{\sqrt{(2t(\cos k\_x + \cos k\_y) + \mu)^2 + \Delta\_s^2 + \Delta\_d^2(\cos k\_x - \cos k\_y)^2}},\tag{88}$$

$$\frac{E\_C}{N} = -\mu - \frac{1}{\pi^2} \int\_0^\pi \mathrm{d}k\_x \int\_0^\pi \mathrm{d}k\_y \sqrt{(2t(\cos k\_x + \cos k\_y) + \mu)^2 + \Lambda\_s^2 + \Lambda\_d^2 (\cos k\_x - \cos k\_y)^2} + \frac{\Lambda\_s^2}{g\_0} + \frac{\Lambda\_d^2}{g\_1} \tag{89}$$

and

$$K(0,0) = \frac{c^2}{mc} \left[ n - \frac{l}{4\pi^2} \int\_0^\pi d\mathbf{k}\_x \int\_0^\pi d\mathbf{k}\_y (k\_x^2 + k\_y^2) \frac{\Lambda\_s^2 + \Lambda\_d^2(\cos k\_x - \cos k\_y)^2}{((2t(\cos k\_x + \cos k\_y) + \mu)^2 + \Lambda\_s^2 + \Lambda\_d^2(\cos k\_x - \cos k\_y)^2)^{\frac{2}{2}}} \right]. \tag{90}$$

Because we are interested in the investigation of the half-filled band case it suffices to put *n* = 1. For this number of electrons one obtains the chemical potential *μ* = 0. In order to obtain those equations for the half-filled band case, one needs to substitute these two values into them. As has been already done for the pure *s*-wave pairing case we can make some analytical calculations. Let us assume that *De* <sup>2</sup> Δ*<sup>s</sup>* and Δ*<sup>d</sup>* Δ*s*. Again let us make use of the following expansions to the first order terms:

$$\sqrt{4t^2(\cos k\_x + \cos k\_y)^2 + \Delta\_s^2 + \Delta\_d^2(\cos k\_x - \cos k\_y)^2} \approx 0$$

$$\Delta\_s + \frac{1}{2} \frac{4t^2(\cos k\_x + \cos k\_y)^2 + \Delta\_d^2(\cos k\_x - \cos k\_y)^2}{\Delta\_s}$$

and

$$\frac{1}{(4t^2(\cos k\_x + \cos k\_y)^2 + \Delta\_s^2 + \Delta\_d^2(\cos k\_x - \cos k\_y)^2)^{\frac{3}{2}}} \approx 0$$

$$\frac{1}{\Delta\_s^3} - \frac{3}{2} \frac{4t^2(\cos k\_x + \cos k\_y)^2 + \Delta\_d^2(\cos k\_x - \cos k\_y)^2}{\Delta\_s^5}.$$

Those expansions are next substituted into the integrals in (89) and (90). After the integration, one obtains the expressions for the ground state energy per lattice site together with the expression for the electromagnetic kernel, namely,

$$\frac{E\_G}{N} \approx -\Delta\_\text{s} - \frac{4t^2 + \Delta\_d^2}{2\Delta\_\text{s}} + \frac{\Delta\_\text{s}^2}{\mathcal{g}\_0} + \frac{\Delta\_d^2}{\mathcal{g}\_1} \tag{91}$$

and

$$\begin{split} K(0,0) &= \quad \frac{\varepsilon^{2}}{m\varepsilon} \Big[ 1 - \frac{1}{4\pi^{2}} \frac{1}{\Lambda\_{e}} \Big[ \frac{2}{3}\pi^{4} - \frac{3}{2} \frac{4^{2} + \Lambda\_{2}^{2}}{\Lambda\_{e}^{2}} \Big( \frac{4}{3}\pi^{3} + \frac{1}{3}\pi^{4} + \frac{\pi^{2}}{2} \Big) \Big] - \\ & \quad - \frac{1}{4\pi^{2}} \frac{1}{\Lambda\_{e}} \frac{\Lambda\_{2}^{2}}{\Lambda\_{e}^{2}} \Big[ \frac{4}{3}\pi^{3} + \frac{\pi^{2}}{2} + \frac{\pi^{4}}{3} - \frac{3}{2} \frac{4^{2} + \Lambda\_{2}^{2}}{\Lambda\_{e}^{2}} \Big( \frac{4}{3}\pi^{2} + \frac{\pi^{2}}{2} + \frac{81}{32}\pi^{2} \Big) + 3 \frac{4^{2} - \Lambda\_{2}^{2}}{\Lambda\_{e}^{2}} \Big( \frac{4}{3}\pi^{3} + \frac{\pi^{2}}{2} \Big) \Big] \Big]. \end{split} \tag{92}$$

If *<sup>t</sup>* <sup>Δ</sup>*<sup>s</sup>* 1 and <sup>Δ</sup>*<sup>d</sup>* <sup>Δ</sup>*<sup>s</sup>* 1 then

$$K(0,0) = \frac{e^2}{mc} \left[1 - \frac{\pi^2}{6} \frac{t}{\Delta\_\mathfrak{s}}\right],\tag{93}$$

thus this is the same expression as the formulas (78) and (85) but with proviso that the hopping parameter is much lesser than the *s*-wave component of the gap for quartets.

#### **9. Numerical Results**

In this section, we would like to present some numerical results regarding some zero-temperature properties of the system represented by the Hamiltonian (1) with *VBCS* = 0. In the beginning, let us consider the scenario with the rectangular DOS. We will show the results for the *s*-wave and the *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairings. For the sake of clarity, let us denote the gap for the *<sup>s</sup>*-wave pairing case by <sup>Δ</sup>*<sup>s</sup>* whereas the gap for the *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing by <sup>Δ</sup>*d*. We shall be using the coupling constant *<sup>g</sup>* <sup>=</sup> 0.5 eV whereas the bandwidth *De* will be equal to three values: 0.5 eV, 1 eV, 2 eV. The results are shown in Table 1.

**Table 1.** The numerical results for the *s*-wave pairing case with the rectangular DOS for the half-filled band with *<sup>n</sup>* <sup>=</sup> 1 and *<sup>μ</sup>* <sup>=</sup> *De* <sup>2</sup> . We have used the Equation (43) for numerical calculations.


From Table 1, one can notice that for the increasing bandwidth, the gap Δ*<sup>s</sup>* clearly decreases. The ground state energy per lattice site is negative and its module increases. One needs to mention that the same results would hold in the BCS case. However, the electromagnetic kernel multiplied by *mc <sup>e</sup>*<sup>2</sup> behaves in another way, namely, it differs from the BCS result that is equal to the unity and its value diminishes with the increasing bandwidth. This means that for a superconductor with the wider band, the Meissner effect is weaker than for one with the narrower band. Of course, the penetration depth is greater in the former case than in the latter one.

The results for the investigated system with the pure *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing symmetry are placed in Table 2.

**Table 2.** The numerical results for the *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing case with the rectangular DOS for the half-filled band with *n* = 1 and *μ* = *De* <sup>2</sup> . We have used the Equations (48)–(50) for numerical calculations.


As one can notice, the behavior of the system with this symmetry of pairing is very similar to that of the former one. The only difference is that the *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing values are lesser than those for the *s*-wave pairing. Especially, it is visible at looking at the gaps.

Now let us discuss the results for the one-dimensional system with the cosine dispersion relation. These results are shown in Table 3.

**Table 3.** The numerical results for the one-dimensional system with the *s*-wave pairing symmetry and the cosine dispersion relation for the half-filled band. Note that *De* = 4*t* and *g* = 0.5 eV with *n* = 1 and *μ* = 0. We have used the Equations (72)–(74) for numerical calculations.


One can see the same tendencies as in the previous tables. The gap Δ decreases for the increasing bandwidth *De*. The ground state energy per lattice site is also negative and its module increases with the increasing bandwidth. The electromagnetic kernel decreases as well though for the bandwidths *De* = 1 eV and *De* = 2 eV this quantity is lesser than <sup>3</sup> 4 *e*2 *mc* . Let us add that the approximated Equation (75) has three real solutions for *g* = 0.5 eV and two hopping parameters *t* = 0.01 eV or *t* = 0.025 eV. For the former value of *t* one obtains: Δ*s*<sup>1</sup> = −0.00981 eV, Δ*s*<sup>2</sup> = 0.010211 eV and Δ*s*<sup>3</sup> = 0.2496 eV whereas for the latter one Δ*s*<sup>1</sup> = −0.023885 eV, Δ*s*<sup>2</sup> = 0.026437 eV and Δ*s*<sup>3</sup> = 0.247448 eV. Note that only Δ*s*<sup>3</sup> in both cases are proper because they fulfil the inequality 2*t* <sup>Δ</sup> 1. Additionally, for *t* = 0.025 eV the value of Δ*s*<sup>3</sup> is very close to that in Table 3 for *De* = 0.1 eV corresponding to *t* = 0.025 eV.

Concerning the two-dimensional system, the situation is very similar to that above. Let us start from the description of the scenario with the *s*-wave pairing symmetry. Table 4 contains the numerical outcomes for four values of the bandwidth *De* = 8*t*. As has already been stated, there are the same tendencies as in the previous cases. It is worth noting that for the bandwidth *De* = 1 eV the paramagnetic term in the electromagnetic kernel is greater than its values for the previously discussed variants of the investigated model. The penetration depth is obviously greater than corresponding values of this quantity. Moreover, the Equation (83) has three real solutions for two values of the hopping parameter *t* = 0.01 eV and *t* = 0.025 eV. They are as follows: for the former one one gets Δ*s*<sup>1</sup> = −0.013768 eV, Δ*s*<sup>2</sup> = 0.014573 eV and Δ*s*<sup>3</sup> = 0.249195 eV while for the latter Δ*s*<sup>1</sup> = −0.033217 eV, Δ*s*<sup>2</sup> = 0.038433 eV and Δ*s*<sup>3</sup> = 0.244785 eV. In this instance one obtains a very good agreement between Δ*s*<sup>3</sup> and those values for *De* = 0.08 eV and *De* = 0.2 eV from Table 4.

**Table 4.** The numerical results for the two-dimensional system with the *s*-wave pairing symmetry and the cosine dispersion relation for the half-filled band. Note that *De* = 8*t*, *g* = 0.5 eV, *μ* = 0. We have used the Equations (79), (81) and (82) with *μ* = 0 for numerical calculations.


In the end, we shall describe the pure *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing case. Here, the investigated quantities behave in a similar way to the corresponding *s*-wave pairing ones. We adopt the same dispersion relation as in the *s*-wave case. The results are shown in Table 5.

**Table 5.** The numerical results for the two-dimensional system with the *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing symmetry and the cosine dispersion relation for the half-filled band. Note that *De* = 8*t*, *g* = 0.5 eV, *μ* = 0. We have used the Equations (87), (89) and (90) with *μ* = 0 and Δ*<sup>s</sup>* = 0 for numerical calculations.


In this instance one can notice that the order parameter for the *s*-wave pairing case is visibly greater than in the *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing one. The same assertion concerns the electromagnetic kernel. However, the ground state energy per lattice site for the latter kind of pairing and *De* = 1 eV is slightly greater than that in the former kind of pairing.

#### **10. Conclusions**

The system with the four-fermion attraction has been investigated in this paper. The original model represented by the Hamiltonian (1) comprises the BCS interaction, however, it has been neglected here for simplicity. The study has been devoted to some zero-temperature values of some quantities in this system such as the order parameter, the ground state energy per lattice site and the electromagnetic response to an external magnetic field. An essential part of the study is to get an answer to the question of how the symmetry of the pairing affects those quantities. It turns out that the influence of this factor cannot be neglected. The difference between both kinds of pairings is significant independent of what type of dispersion relation we will use. Basically, all quantities with the *dx*<sup>2</sup>−*y*<sup>2</sup> -wave pairing are lesser than those with the *s*-wave pairing. Moreover, the electromagnetic kernel for the system in the frame of the tight-binding approximation (the cosine dispersion relation) is lesser than this kernel for the system with the rectangular DOS. In both of cases, the same bandwidth is used. The example with *De* = 1 eV is visible in Tables from the previous section. The interesting result is found for the system with the cosine dispersion relation. For a very small ratio *<sup>t</sup>* <sup>Δ</sup>*<sup>g</sup>* the electromagnetic kernel has approximately the same value independent of the dimension and type of pairing.

There are some open problems not undertaken in this paper. For example, properties of the system at finite temperatures are not studied here. From [18,19], it is known that there should be some differences between the conventional BCS system and the system studied in this paper. One can mention the character of the phase transition and the value of the critical temperature. Moreover, the incorporation of the BCS interaction can, to large extent, affect different properties of the system. Additionally, different types of pairings in such interactions can lead to many interesting phenomena.

**Funding:** This research received no external funding.

**Acknowledgments:** The author would like to thank R. Szcze¸sniak and W. Leo ´ ´ nski for their support during the writing of this text.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


c 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Symmetry* Editorial Office E-mail: symmetry@mdpi.com www.mdpi.com/journal/symmetry

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18