BCS-BEC crossover and superconductor-insulator transition in Hopf-linked Graphene layers: Hopfene

We employ the attractive Hubbard Hamiltonian [61]

$$\begin{eqnarray}&&\hat{H}=-\displaystyle \sum _{{ij}\sigma }{t}_{{ij}}{\hat{\psi }}_{i\sigma }^{\dagger }{\hat{\psi }}_{j\sigma }-\mu \displaystyle \sum _{i\sigma }{\hat{n}}_{i\sigma }-U\displaystyle \sum _{i}{\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow },\end{eqnarray} \tag{ 1 }$$

where ${\hat{\psi }}_{i\sigma }^{\dagger }$ and ${\hat{\psi }}_{i\sigma }$ are creation and annihilation operators at the atomic lattice site i for spin σ, respectively, t_ij is the transfer energy between the lattice sites i and j, μ is chemical potential, ${\hat{n}}_{i\sigma }$ is the number operator, and U is the attractive interaction strength. The real crystal of Hopfene would contain defects and disorder of impurities. In particular, we think it is important to consider new type of defects, named topological defects, which break translational long-range order, while preserving the strong σ-bonds (figure 1 (e)), because Hopfene is made of soft Graphene layers against deformations. The structure will be maintained as far as the topological long-range order is preserved, even if the perfect translational symmetry is broken. These topological defects would not affect s-wave superconductivity, as far as the radius of the Cooper pair is much larger than the size of the topological defects. In this paper, we focus on (2,2)-Hopfene (figure 1) and restrict the electron transfers within the nearest neighbour sites [44, 54]. In (2,2)-Hopfene, the numbers of sub-lattice atoms for all sublattices are the same ${ \mathcal N }={{ \mathcal N }}_{{\rm{A}}}={{ \mathcal N }}_{{\rm{B}}}={{ \mathcal N }}_{{\rm{C}}}={{ \mathcal N }}_{{\rm{D}}}$, so that it is convenient to split the sum over all lattice sites into the sum within sub-lattices as

$$\begin{eqnarray}&&\displaystyle \sum _{i}=\displaystyle \sum _{{\rm{S}}={\rm{A}},{\rm{B}},{\rm{C}},{\rm{D}}}\ \displaystyle \sum _{i\in {\rm{S}}}^{{ \mathcal N }}.\end{eqnarray} \tag{ 2 }$$

This means that we can assign a number from 1 to ${ \mathcal N }$ for each atom with a sublattice index from A to D to describe the lattice coordinate of an atom. The total number of lattice sites is ${{ \mathcal N }}_{\mathrm{tot}}=4{ \mathcal N }=8{{ \mathcal N }}_{{\rm{U}}}$ with the total number of the unit cell, ${{ \mathcal N }}_{{\rm{U}}}$, and the number of sub-lattice atoms for each sublattice in the unit cell is 2 and the total number of atoms in the unit cell is 8. Similarly, it is convenient to describe ${\hat{\psi }}_{i\sigma }$ with a sublattice index as ${\hat{\psi }}_{i\sigma }={\hat{\psi }}_{i\sigma }^{{\rm{S}}}$, and consequently we can define Fourier transformation and its inverse as

$$\begin{eqnarray}&&{\hat{\psi }}_{{\bf{k}}\sigma }^{{\rm{S}}}=\displaystyle \frac{1}{\sqrt{{ \mathcal N }}}\displaystyle \sum _{i\in {\rm{S}}}^{{ \mathcal N }}{\hat{\psi }}_{i\sigma }^{{\rm{S}}}{{\rm{e}}}^{-i{\bf{k}}\cdot {{\bf{x}}}_{i}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\hat{\psi }}_{i\sigma }^{{\rm{S}}}=\displaystyle \frac{1}{\sqrt{{ \mathcal N }}}\displaystyle \sum _{{\bf{k}}}^{{ \mathcal N }}{\hat{\psi }}_{{\bf{k}}\sigma }^{{\rm{S}}}{{\rm{e}}}^{i{\bf{k}}\cdot {{\bf{x}}}_{i}},\end{eqnarray} \tag{ 4 }$$

respectively. We will describe the wave number vectors, k, equivalent to the momentum p = ℏk, in the unit of $2\pi ({a}^{-1},{a}^{-1},{c}^{-1})$ and ℏ is the Dirac constant, which is the Plank constant divided by 2π.

Alternatively, we can also define the field-vector operator, ${\hat{{\boldsymbol{\psi }}}}_{i\sigma }$, with a sub-lattice index similar to a spinor description as ${\hat{{\boldsymbol{\psi }}}}_{i\sigma }=\left({\hat{\psi }}_{i\sigma }^{{\rm{A}}},{\hat{\psi }}_{i\sigma }^{{\rm{B}}},{\hat{\psi }}_{i\sigma }^{{\rm{C}}},{\hat{\psi }}_{i\sigma }^{{\rm{D}}}\right)=\left({a}_{i\sigma },{b}_{i\sigma },{c}_{i\sigma },{d}_{i\sigma }\right)$, and its Fourier transform is ${\hat{{\boldsymbol{\psi }}}}_{{\bf{k}}\sigma }\,=\left({a}_{{\bf{k}}\sigma },{b}_{{\bf{k}}\sigma },{c}_{{\bf{k}}\sigma },{d}_{{\bf{k}}\sigma }\right)$. This simplification was possible because of the same number of lattice sites among sub-lattices. If this was not the case, it is required to define two coordination vectors: one for the unit cell and the other for the relative vectors inside the cell.

Below, we will examine the 1st term of $\hat{H}$ as the kinetic energy Hamiltonian ${\hat{H}}_{\mathrm{kin}}$, the 2nd term as the chemical potential Hamiltonian ${\hat{H}}_{\mathrm{chem}}$, and the last term as the interaction Hamiltonian ${\hat{H}}_{\mathrm{int}}$ in more details.

${\hat{H}}_{\mathrm{kin}}$ represents delocalised π-electrons in the carbon network, and we expect 1 electron per atom without doping. The other 3 electrons among 4 electrons in the outer shell of a carbon atom are forming strong σ bonds. We assume the nearest neighbour transfer energy of −t within Graphene layers and that of −t_H for Hopf-links. Then, we obtain

$$\begin{eqnarray}{\hat{H}}_{\mathrm{kin}}=\displaystyle \sum _{{\bf{k}}\sigma }{\hat{{\boldsymbol{\psi }}}}_{{\bf{k}}\sigma }^{\dagger }\left(\begin{array}{cccc}0 & {h}_{\mathrm{AB}} & {h}_{\mathrm{AC}} & {h}_{\mathrm{AD}}\\ {h}_{\mathrm{AB}}^{* } & 0 & {h}_{\mathrm{BC}} & {h}_{\mathrm{BD}}\\ {h}_{\mathrm{AC}}^{* } & {h}_{\mathrm{BC}}^{* } & 0 & {h}_{\mathrm{CD}}\\ {h}_{\mathrm{AD}}^{* } & {h}_{\mathrm{BD}}^{* } & {h}_{\mathrm{CD}}^{* } & 0\end{array}\right){\hat{{\boldsymbol{\psi }}}}_{{\bf{k}}\sigma },\end{eqnarray} \tag{ 5 }$$

where we have defined

$$\begin{eqnarray*}\begin{array}{rcl}{h}_{\mathrm{AB}} & = & -t\left({{\rm{e}}}^{{{ik}}_{x}a/3}+2{{\rm{e}}}^{-{{ik}}_{x}a/6}\cos ({k}_{z}c/2)\right),\\ {h}_{\mathrm{CD}} & = & -t\left({{\rm{e}}}^{{{ik}}_{y}a/3}+2{{\rm{e}}}^{-{{ik}}_{y}a/6}\cos ({k}_{z}c/2)\right),\\ {h}_{\mathrm{AC}} & = & -{t}_{{\rm{H}}}\left({{\rm{e}}}^{{{ik}}_{x}a/6}{{\rm{e}}}^{{{ik}}_{y}a/3}+{{\rm{e}}}^{-{{ik}}_{x}a/3}{{\rm{e}}}^{-{{ik}}_{y}a/6}+2{{\rm{e}}}^{{{ik}}_{x}a/6}{{\rm{e}}}^{-{{ik}}_{y}a/6}\cos ({k}_{z}c/2)\right),\\ {h}_{\mathrm{AD}} & = & -{t}_{{\rm{H}}}\left({{\rm{e}}}^{{{ik}}_{x}a/6}{{\rm{e}}}^{-{{ik}}_{y}a/3}+{{\rm{e}}}^{-{{ik}}_{x}a/3}{{\rm{e}}}^{{{ik}}_{y}a/6}+2{{\rm{e}}}^{{{ik}}_{x}a/6}{{\rm{e}}}^{{{ik}}_{y}a/6}\cos ({k}_{z}c/2)\right),\\ {h}_{\mathrm{BC}} & = & -{t}_{{\rm{H}}}\left({{\rm{e}}}^{-{{ik}}_{x}a/6}{{\rm{e}}}^{{{ik}}_{y}a/3}+{{\rm{e}}}^{{{ik}}_{x}a/3}{{\rm{e}}}^{-{{ik}}_{y}a/6}+2{{\rm{e}}}^{-{{ik}}_{x}a/6}{{\rm{e}}}^{-{{ik}}_{y}a/6}\cos ({k}_{z}c/2)\right),\\ {h}_{\mathrm{BD}} & = & -{t}_{{\rm{H}}}\left({{\rm{e}}}^{-{{ik}}_{x}a/6}{{\rm{e}}}^{-{{ik}}_{y}a/3}+{{\rm{e}}}^{{{ik}}_{x}a/3}{{\rm{e}}}^{{{ik}}_{y}a/6}+2{{\rm{e}}}^{-{{ik}}_{x}a/6}{{\rm{e}}}^{{{ik}}_{y}a/6}\cos ({k}_{z}c/2)\right).\end{array}\end{eqnarray*}$$

If we define the Hamiltonian matrix for Graphene as

$$\begin{eqnarray*}{{ \mathcal H }}_{{\rm{G}}}(h)=\left(\begin{array}{cc}0 & h\\ {h}^{* } & 0\end{array}\right),\end{eqnarray*}$$

and that of Hopf-links as

$$\begin{eqnarray*}{h}_{{\rm{H}}}=\left(\begin{array}{cc}{h}_{\mathrm{AC}} & {h}_{\mathrm{AD}}\\ {h}_{\mathrm{BC}} & {h}_{\mathrm{BD}}\end{array}\right),\end{eqnarray*}$$

then, it can be summarised as

$$\begin{eqnarray}{\hat{H}}_{\mathrm{kin}}=\displaystyle \sum _{{\bf{k}}\sigma }{\hat{{\boldsymbol{\psi }}}}_{{\bf{k}}\sigma }^{\dagger }\left(\begin{array}{cc}{{ \mathcal H }}_{{\rm{G}}}({h}^{\mathrm{AB}}) & {{ \mathcal H }}_{{\rm{H}}}\\ {{ \mathcal H }}_{{\rm{H}}}^{* } & {{ \mathcal H }}_{{\rm{G}}}({h}^{\mathrm{CD}})\end{array}\right){\hat{{\boldsymbol{\psi }}}}_{{\bf{k}}\sigma }.\end{eqnarray} \tag{ 6 }$$

In the limit of the absence of the Hopf-link transfer, ${t}_{{\rm{H}}}\to 0$, the energy dispersion reduces to that of Graphene for Dirac Fermions

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{AB}}({\bf{k}})=\pm \sqrt{| {h}_{\mathrm{AB}}{| }^{2}}=\pm \sqrt{1+4{\cos }^{2}({k}_{z}c/2)+4\cos ({k}_{x}a/2)\cos ({k}_{z}c/2)}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{CD}}({\bf{k}})=\pm \sqrt{| {h}_{\mathrm{CD}}{| }^{2}}=\pm \sqrt{1+4{\cos }^{2}({k}_{z}c/2)+4\cos ({k}_{y}a/2)\cos ({k}_{z}c/2)},\end{eqnarray} \tag{ 8 }$$

with t = 2.8 eV [29, 33, 34]. In general, ${t}_{{\rm{H}}}\ne 0$, we need to diagonalise ${\hat{H}}_{\mathrm{kin}}$ numerically at each k, and we can obtain the band diagram with 4 energy dispersions, ₁ (k) − ₄ (k), labelled from the bottom. After the diagonalisation, the state will carry the band index, α, instead of S, and we obtain

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{kin}}=\displaystyle \sum _{{\bf{k}}\sigma \alpha }{\epsilon }_{\alpha }{\hat{\psi }}_{{\bf{k}}\sigma \alpha }^{\dagger }{\hat{\psi }}_{{\bf{k}}\sigma \alpha }.\end{eqnarray} \tag{ 9 }$$

Figure 2 shows examples of calculated tight-binding band structures. At t_H = 0, electrons in Graphene layers cannot transfer from each other through Hopf-links, so that electrons are confined within Garphene layers. In this case, the band dispersions are flat against the directions perpendicular to the Graphene layers, so that we have flat-bands along k_x, from (0,0.5,1) to (0.5,0.5,1), and k_y, from (0,0,1) to (0,0.5,1) (figure 2 (a)). There are 10 Dirac points at (±1, 0, ±1/3), (0, ±1, ±1/3), and (0,0, ±2/3). Among these points, at (0,0, ±2/3), Dirac points coming from horizontal and vertical Graphene layers are degenerate (figure 2(a)). Therefore, we can recognise 2 dispersion curves are completely overlapping around the Dirac point (figure 2(a)). At ${t}_{{\rm{H}}}\ne 0$, Dirac point at (0, 0, 2/3) split into 2 Dirac points (figures 2(b) and (c)) due to the formation of bonding and anti-bonding states through Hopf-links. We have compared these tight-binding results with the first-principle calculations [55] and confirmed reasonable agreements and extracted t_H = 1.5 eV. All calculations below are based on this parameter with t = 2.8 eV [29, 33, 34] and t_H = 1.5 eV.

Zoom In Zoom Out Reset image size

Due to its analytic nature of the tight-binding model, it is easy to calculate energy band dispersions, _α = _α(k_x, k_y, k_z). However, we must fix 1 component in order to describe in a 3D form (figure 3). We can recognise Dirac points in the (k_x, k_z) plot at fixed k_y = 0 (figure 3(a)), and the dispersion is symmetric under the exchange of k_x and k_y, so that exactly the same energy dispersion is obtained in the (k_y, k_z) plot at k_x = 0 (not shown). At k_z = 0, we do not see Dirac points, and a simple band interception was found at the zone edge of (1, 1, 0). We also found the broken particle-hole symmetry, which means _α (k_x, k_y, k_z) is symmetric against  = 0 by the change of the sign. The particle-hole symmetries were found in Graphene layers, as seen in _AB and _CD of equations (7) and (8) while it was broken at ${t}_{{\rm{H}}}\ne 0$. This could be understood simply by diagonalising 2 and 3 sites models as shown in table 1. The simplest 2 sites model will give two energy states, corresponding to particle and hole state with the energy of t and −t, respectively. If we introduce 1 more site allowing equivalent transfers to these 2 sites, mimicking Hopf-links, we obtain the asymmetric eigenstates, $-t-2{t}_{{\rm{H}}}^{2}/t$, $2{t}_{{\rm{H}}}^{2}/t$, and t, in the limit of $t\gg {t}_{{\rm{H}}}$. Therefore, the energy levels depend on the sign of t and the particle-hole symmetry is broken. On the other hand, the sign change of t_H does not break the particle-hole symmetry. This comes from the fact that 2 transfers through Hopf-link are required to come back to the original Graphene sites, so that the sign change of t_H cannot contribute. More generally in Hopfene, even number of transfers thorough Hopf-links are required to come back, so that the obtained energy bands _α (k_x, k_y, k_z) did not change by ${t}_{{\rm{H}}}\to -{t}_{{\rm{H}}}$. Please also note that the spacial reflection symmetry ${\bf{k}}\to -{\bf{k}}$ (_α (k) = _α (−k)) and the time reversal symmetry $t\to -t$ are still maintained without applying magnetic fields. It is also symmetric with the exchange of k_x and k_y in (2,2)-Hopfene, due to the symmetry between horizontal and vertical directions, but it would not be the case in (n, m)-Hopfene with $n\ne m$.

Zoom In Zoom Out Reset image size

Table 1.  Two and three sites models for the impacts of signs in transfer energies.

Model	Hamiltonian	Eigenvalues	Eigenvalues ($t\gg {t}_{{\rm{H}}}$)
2 sites	$\left(\begin{array}{cc}0 & -t\\ -t & 0\end{array}\right)$	±t	±t
2 sites + Hopf-link	$\left(\begin{array}{ccc}0 & -t & -{t}_{{\rm{H}}}\\ -t & 0 & -{t}_{{\rm{H}}}\\ -{t}_{{\rm{H}}} & -{t}_{{\rm{H}}} & 0\end{array}\right)$	t, $-\tfrac{1}{2}\left(t\pm \sqrt{{t}^{2}+8{t}_{{\rm{H}}}^{2}}\right)$	$-t-\tfrac{2{t}_{{\rm{H}}}^{2}}{t}$, $\tfrac{2{t}_{{\rm{H}}}^{2}}{t}$, t
Negative transfer	$\left(\begin{array}{ccc}0 & -t & -{t}_{{\rm{H}}}\\ t & 0 & -{t}_{{\rm{H}}}\\ {t}_{{\rm{H}}} & -{t}_{{\rm{H}}} & 0\end{array}\right)$	−t, $\tfrac{1}{2}\left(t\pm \sqrt{{t}^{2}+8{t}_{{\rm{H}}}^{2}}\right)$	−t, $-\tfrac{2{t}_{{\rm{H}}}^{2}}{t}$, $t+\tfrac{2{t}_{{\rm{H}}}^{2}}{t}$
Negative Hopf-link	$\left(\begin{array}{ccc}0 & -t & {t}_{{\rm{H}}}\\ -t & 0 & {t}_{{\rm{H}}}\\ {t}_{{\rm{H}}} & {t}_{{\rm{H}}} & 0\end{array}\right)$	t, $-\tfrac{1}{2}\left(t\pm \sqrt{{t}^{2}+8{t}_{{\rm{H}}}^{2}}\right)$	$-t-\tfrac{2{t}_{{\rm{H}}}^{2}}{t}$, $\tfrac{2{t}_{{\rm{H}}}^{2}}{t}$, t

Another interesting feature of the band diagram of Hopfene is the existence of nodal lines [39] (figure 4). This is coming from the original flat dispersion relationships of _AB and _CD against k_y and k_x, which are directions perpendicular to horizontal AB-layers and vertical CD-layers, respectively. Therefore, the bands close to the Dirac points are also dispersionless along these directions, so that Dirac points are connected as nodal lines. Once nodal lines were formed, they cannot be eliminated upon the adiabatic introduction of t_H, since nodal lines are topologically protected. In the example of figure 3, nodal lines have dispersion along k_x = k_y − 1 at k_z = 1/2 but the lines are in perfect agreement with the adjacent bands (figure 3(b)), confirming the degeneracies and topological features of band diagrams.

Zoom In Zoom Out Reset image size

In the presence of t_H, electronic states between horizontal and vertical Graphene layers are hybridised to form topologically unique Fermi surfaces [54]. The band structures can be highlighted by the existence of 3D, 2D, and 1D Weyl and Dirac Fermions [54] with flat-band features (figure 2). The Fermi energy of the non-doped (2,2)-Hopfene was E_F = 1.1 eV, but this could be changed by doping alkali-atoms [15, 56, 57] or NH₃ [71] at interstitial sites due to the large spacing in the unit cell [44, 54]. Here, we describe typical structures of Fermi surfaces in (2,2)-Hopfene.

3D Weyl Fermions exist around at the Dirac point with the energy of E_d = −1.0 eV, where ₃ and ₂ have linear dispersions along all 3 directions, k_x, k_y, and k_z, and therefore, the spatial quantum confinements are strongest [54]. Figure 5 shows the Fermi surface of 3D Weyl Fermions below (E_F = −0.8 eV, figures 5(a) and (c)) and above (E_F = −1.2 eV, figures 5(b) and (d)) this point at E_d = −1.0 eV. Fermi surfaces of 3D Weyl Fermions are small, because of the limited DOS due to strong confinements and linear dispersions. The reduced DOS is not favourable for superconductivity [67–70], as we shall see in the following sections for Hopfene. The Fermi surface of ₁ in figure 5(d) is particularly unique, because the 2 edges of the Fermi surface are connected to the Fermi surfaces of ₂ by point-contacts. This is originated from the nodal lines, shown in figure 4. Due to the close proximity of the Dirac point to the nodal lines, the Fermi surfaces of ₁ were terminated on the nodal lines, so that Fermi surfaces of ₁ and ₂ were connected by points without overlapping each other. This topological feature is maintained to a certain extent, and the filling of these bands were competing: as we increase E_F, ₁ is eventually occupied completely, leading to the reduction of the size of the Fermi surface of lower valley (Figure 5(d)) for 3D Weyl Fermions, and the flat-bands of ₂ is taking over the area, while further increase of E_F induces the reduction of the volume of Fermi surfaces of ₂, and the Fermi surface of the upper valley (figure 5(c)) starts to take over.

Zoom In Zoom Out Reset image size

We also found the weaker confinements at the 1D Dirac point of E_d = 2.7 eV, where the confinement is strong only along k_x direction (figure 6). The Fermi surfaces located these valleys are more spreading along vertical and diagonal directions, since band dispersions along these directions are flat (figures 6(b) and (c)). We can also recognise that the close proximity of nodal lines near E_d due to the point contacts of the Fermi surfaces of ₄ with ₃, at (0, 0, 1), which is equivalent to (0, 0, −1).

Zoom In Zoom Out Reset image size

A similar point-contacts of the Fermi surfaces were also found at E_d = −4.3 eV (figure 7), where 4 Dirac points are found and each Fermi surface of the upper valleys (₂) have 4 point contacts to the Fermi surfaces of the lower band (₁). As we increase E_F, the flat-bands of ₁ (figure 7(a)) are gradually filled, leading to the reduced Fermi surfaces, and eventually the Fermi surfaces of another flat-bands of ₂ (figure 7(b)) start to grow. These changes of the Fermi surfaces upon increasing E_F will give important consequences on the evolutions of superconductivity using these valleys, and the larger DOS is favourable for larger superconducting order parameters, as expected from the BCS theory [3, 4, 6, 61, 62].

Zoom In Zoom Out Reset image size

The DOS per unit volume per spin for the band α is given by

$$\begin{eqnarray}&&{D}_{V\sigma }^{\alpha }(\epsilon )=\displaystyle \frac{1}{V}\displaystyle \sum _{{\bf{k}}}\delta \left(\epsilon -{\epsilon }^{\alpha }({\bf{k}})\right)\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{2}{{v}_{{\rm{U}}}}{\int }_{-1/2}^{1/2}{{dk}}_{x}{\int }_{-1/2}^{1/2}{{dk}}_{y}{\int }_{-1/2}^{1/2}{{dk}}_{z}\delta \left(\epsilon -{\epsilon }^{\alpha }({\bf{k}})\right)\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{2}{{v}_{{\rm{U}}}}{\int }_{-1}^{1}\displaystyle \frac{{{dk}}_{x}}{2}{\int }_{-1}^{1}\displaystyle \frac{{{dk}}_{y}}{2}{\int }_{-1}^{1}\displaystyle \frac{{{dk}}_{z}}{2}\delta \left(\epsilon -{\epsilon }^{\alpha }({\bf{k}})\right),\end{eqnarray} \tag{ 12 }$$

where V is the volume of the system, v_U = abc is the volume of the unit cell, a factor of 2 in equation (11) is coming from the number of each sub-lattice atoms in the unit cell, and equation (12) is described by the integration in the extended Brillouin zone of k_x ∈ (−2π/a, 2π/a), k_y ∈ (−2π/a, 2π/a), and k_z ∈ (−2π/c, 2π/c). The DOS per unit cell per spin for the band α is

$$\begin{eqnarray}&&{D}_{{\rm{U}}\sigma }^{\alpha }(\epsilon )={v}_{{\rm{U}}}{D}_{V\sigma }^{\alpha }(\epsilon ).\end{eqnarray} \tag{ 13 }$$

The total number of atoms per unit cell is 4 × 2 = 8, so that the DOS per atom per spin for the band α is given by

$$\begin{eqnarray}&&{D}_{0\sigma }^{\alpha }(\epsilon )=\displaystyle \frac{1}{4}{\int }_{-1}^{1}\displaystyle \frac{{{dk}}_{x}}{2}{\int }_{-1}^{1}\displaystyle \frac{{{dk}}_{y}}{2}{\int }_{-1}^{1}\displaystyle \frac{{{dk}}_{z}}{2}\delta \left(\epsilon -{\epsilon }^{\alpha }({\bf{k}})\right),\end{eqnarray} \tag{ 14 }$$

and the total DOS per atom per spin is

$$\begin{eqnarray}&&{D}_{0\sigma }(\epsilon )=\displaystyle \sum _{\alpha =1}^{4}{D}_{0\sigma }^{\alpha }(\epsilon ),\end{eqnarray} \tag{ 15 }$$

which satisfies the sum rule

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }d\epsilon {D}_{0\sigma }(\epsilon )=\displaystyle \sum _{\alpha =1}^{4}\displaystyle \frac{1}{4}=1,\end{eqnarray} \tag{ 16 }$$

meaning that each atom can accommodate 1 electron per spin, so that the maximum number of 2 electrons can occupy an atom including spin.

Figure 8 shows the calculated DOS of non-interacting electrons at t_H = 1.5 eV. The band structure predicts semi-metal characteristics with peaks, coming from flat-bands. The point contact of ₁ and ₃ bands at  = −1.0 eV is showing 3D Dirac point, where the peak of DOS of the band ₂ is coming from the flat-band feature (figure 5).

Zoom In Zoom Out Reset image size

Finally, we will focus on the impact of the attractive interaction Hamiltonian,

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{int}}=-U\displaystyle \sum _{{\rm{S}}={\rm{A}},{\rm{B}},{\rm{C}},{\rm{D}}}\ \displaystyle \sum _{i\in {\rm{S}}}{\hat{n}}_{i\uparrow }^{{\rm{S}}}{\hat{n}}_{i\downarrow }^{{\rm{S}}},\end{eqnarray} \tag{ 17 }$$

where ${\hat{n}}_{i\sigma }^{{\rm{S}}}$ is the number operator at the site i in the sub-lattice S. We introduce Fourier transformation and its inverse of ${\hat{n}}_{i\sigma }^{{\rm{S}}}$ as

$$\begin{eqnarray}&&{\hat{\rho }}_{{\bf{q}}\sigma }^{{\rm{S}}}=\displaystyle \sum _{i\in {\rm{S}}}^{{ \mathcal N }}{\hat{n}}_{i\sigma }^{{\rm{S}}}{{\rm{e}}}^{-i{\bf{q}}\cdot {{\bf{x}}}_{i}}=\displaystyle \sum _{{\bf{k}}}{\hat{\psi }}_{{\bf{k}}\sigma }^{{\rm{S}}\dagger }{\hat{\psi }}_{{\bf{k}}+{\bf{q}}\sigma }^{{\rm{S}}}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\hat{n}}_{i\sigma }^{{\rm{S}}}=\displaystyle \frac{1}{{ \mathcal N }}\displaystyle \sum _{{\bf{q}}}{\hat{\rho }}_{{\bf{q}}\sigma }^{{\rm{S}}}{{\rm{e}}}^{i{\bf{q}}\cdot {{\bf{x}}}_{i}}.\end{eqnarray} \tag{ 19 }$$

Inserting these into ${\hat{H}}_{\mathrm{int}}$, we obtain

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{int}}=-\displaystyle \frac{U}{{ \mathcal N }}\displaystyle \sum _{{\rm{S}}={\rm{A}},{\rm{B}},{\rm{C}},{\rm{D}}}\ \displaystyle \sum _{i\in {\rm{S}}}{\hat{\rho }}_{{\bf{q}}\uparrow }^{{\rm{S}}}{\hat{\rho }}_{-{\bf{q}}\downarrow }^{{\rm{S}}},\end{eqnarray} \tag{ 20 }$$

which shows the interaction works within the sub-lattice, because of the local nature of the Hubbard interaction. We examine how the interaction works among sub-bands after the diagonalisation of the kinetic Hamiltonian as

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{kin}}=\displaystyle \sum _{{\rm{S}},{{\rm{S}}}^{{\prime} }}\ \displaystyle \sum _{{\bf{k}},\sigma }{\hat{\psi }}_{{\bf{k}}\sigma }^{{\rm{S}}}{{ \mathcal H }}_{\mathrm{kin}}^{{\mathrm{SS}}^{{\prime} }}({\bf{k}}){\hat{\psi }}_{{\bf{k}}\sigma }^{{{\rm{S}}}^{{\prime} }}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&=\ \displaystyle \sum _{\alpha }\ \displaystyle \sum _{{\bf{k}},\sigma }{\hat{\psi }}_{{\bf{k}}\sigma }^{\alpha }{{ \mathcal H }}_{\mathrm{kin}}^{\alpha \alpha }({\bf{k}}){\hat{\psi }}_{{\bf{k}}\sigma }^{\alpha },\end{eqnarray} \tag{ 22 }$$

where the Hamiltonian matrix

$$\begin{eqnarray}&&{{ \mathcal H }}_{\mathrm{kin}}^{\alpha {\alpha }^{{\prime} }}({\bf{k}})={U}_{\alpha {\rm{S}}}^{\dagger }({\bf{k}}){{ \mathcal H }}_{\mathrm{kin}}^{{\mathrm{SS}}^{{\prime} }}({\bf{k}}){U}_{{{\rm{S}}}^{{\prime} }{\alpha }^{{\prime} }}({\bf{k}})\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&=\ {\epsilon }_{\alpha }({\bf{k}}){\delta }_{\alpha {\alpha }^{{\prime} }},\end{eqnarray} \tag{ 24 }$$

was diagonalised by the unitary matrix U_{S α} . It is not trivial whether the interaction is still within the intra-bands, because U_{S α} is included in the Fourier transform of the density operator as

$$\begin{eqnarray}&&{\hat{\rho }}_{{\bf{q}}\sigma }^{{\rm{S}}}=\displaystyle \sum _{{\bf{k}}}\displaystyle \sum _{\alpha {\alpha }^{{\prime} }}{U}_{{\rm{S}}\alpha }^{* }({\bf{k}}){U}_{{\rm{S}}{\alpha }^{{\prime} }}({\bf{k}}+{\bf{q}}){\hat{\psi }}_{{\bf{k}}\sigma }^{\alpha \dagger }{\hat{\psi }}_{{\bf{k}}+{\bf{q}}\sigma }^{{\alpha }^{{\prime} }},\end{eqnarray} \tag{ 25 }$$

and we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\hat{H}}_{\mathrm{int}} & = & -\displaystyle \frac{U}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}{k}^{{\prime} }q}\ \displaystyle \sum _{{\rm{S}}}\ \displaystyle \sum _{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}{\alpha }_{4}}\\ & & {U}_{{\rm{S}}{\alpha }_{1}}^{* }({\bf{k}}){U}_{{\rm{S}}{\alpha }_{3}}^{* }\left(-{k}^{{\prime} }\right){U}_{{\rm{S}}{\alpha }_{4}}(-({k}^{{\prime} }+q)){U}_{{\rm{S}}{\alpha }_{2}}({\bf{k}}+{\bf{q}})\\ & & {\hat{\psi }}_{{\bf{k}}\uparrow }^{{\alpha }_{1}\dagger }{\hat{\psi }}_{-{k}^{{\prime} }\downarrow }^{{\alpha }_{3}\dagger }{\hat{\psi }}_{-({k}^{{\prime} }+q)\downarrow }^{{\alpha }_{4}}{\hat{\psi }}_{{\bf{k}}+{\bf{q}}\uparrow }^{{\alpha }_{2}}.\end{array}\end{eqnarray} \tag{ 26 }$$

We show the products of unitary matrices do not contribute in the end. We utilise the fundamental character of the Unitary matrix, ${U}^{-1}={U}^{\dagger }$, namely

$$\begin{eqnarray}&&{\left({U}^{-1}({\bf{k}})\right)}_{{{\rm{S}}}_{1}{{\rm{S}}}_{2}}={U}_{{{\rm{S}}}_{2}{{\rm{S}}}_{1}}^{* }({\bf{k}}),\end{eqnarray} \tag{ 27 }$$

and the reflection symmetry of the system

$$\begin{eqnarray}&&{U}_{{{\rm{S}}}_{1}{{\rm{S}}}_{2}}(-{\bf{k}}),={U}_{{{\rm{S}}}_{1}{{\rm{S}}}_{2}}^{* }({\bf{k}}).\end{eqnarray} \tag{ 28 }$$

In the Cooper channel, the contribution from the zero center-of-mass momentum of the Cooper pair dominates due to BEC nature of superconductivity, so that the pairing of k' = k within the sub-babnds in ${\hat{H}}_{\mathrm{int}}$ dominates. Then, we consider

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{{\alpha }_{1}={\alpha }_{3}}{U}_{{\rm{S}}{\alpha }_{1}}^{* }({\bf{k}}){U}_{{\rm{S}}{\alpha }_{3}}^{* }(-{\bf{k}}) & = & \displaystyle \sum _{{\alpha }_{1}}{\left(U{\left({\bf{k}}\right)}^{-1}\right)}_{{\alpha }_{1}{\rm{S}}}{\left(U({\bf{k}})\right)}_{{\rm{S}}{\alpha }_{1}}\\ & = & {\left(U({\bf{k}})U{\left({\bf{k}}\right)}^{-1}\right)}_{\mathrm{SS}}\\ & = & 1.\end{array}\end{eqnarray} \tag{ 29 }$$

Similarly, we confirm

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{{\alpha }_{4}={\alpha }_{2}}{U}_{{\rm{S}}{\alpha }_{4}}^{* }(-({\bf{k}}+{\bf{q}})){U}_{{\rm{S}}{\alpha }_{2}}^{* }({\bf{k}}+{\bf{q}}) & = & \displaystyle \sum _{{\alpha }_{2}}{\left(U({\bf{k}}+{\bf{q}})\right)}_{{\rm{S}}{\alpha }_{2}}{\left(U{\left({\bf{k}}+{\bf{q}}\right)}^{-1}\right)}_{{\alpha }_{2}{\rm{S}}}\\ & = & {\left(U({\bf{k}}+{\bf{q}})U{\left({\bf{k}}+{\bf{q}}\right)}^{-1}\right)}_{\mathrm{SS}}\\ & = & 1.\end{array}\end{eqnarray} \tag{ 30 }$$

Therefore, the dominant interaction would be the paring within sub-bands

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{Cooper}}=-\displaystyle \frac{U}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}{k}^{{\prime} }}\ \displaystyle \sum _{\alpha }{\hat{\psi }}_{{\bf{k}}\uparrow }^{\alpha \dagger }{\hat{\psi }}_{-{\bf{k}}\downarrow }^{\alpha \dagger }{\hat{\psi }}_{-{k}^{{\prime} }\downarrow }^{\alpha }{\hat{\psi }}_{{{\bf{k}}}^{{\prime} }\uparrow }^{\alpha },\end{eqnarray} \tag{ 31 }$$

and the factor of unitary matrices become unity.

Another dominant source of interaction is the Hartree–Fock channel at q = 0

$$\begin{eqnarray}\begin{array}{rcl}{\hat{H}}_{\mathrm{HF}} & = & -\displaystyle \frac{U}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}{k}^{{\prime} }}\ \displaystyle \sum _{{\rm{S}}}\ \displaystyle \sum _{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}{\alpha }_{4}}\\ & & {U}_{{\rm{S}}{\alpha }_{1}}^{* }({\bf{k}}){U}_{{\rm{S}}{\alpha }_{2}}({\bf{k}}){U}_{{\rm{S}}{\alpha }_{3}}^{* }\left(-{k}^{{\prime} }\right){U}_{{\rm{S}}{\alpha }_{4}}(-({k}^{{\prime} }))\\ & & {\hat{\psi }}_{{\bf{k}}\uparrow }^{{\alpha }_{1}\dagger }{\hat{\psi }}_{{\bf{k}}\uparrow }^{{\alpha }_{2}}{\hat{\psi }}_{-{k}^{{\prime} }\downarrow }^{{\alpha }_{3}\dagger }{\hat{\psi }}_{-{k}^{{\prime} }\downarrow }^{{\alpha }_{4}},\end{array}\end{eqnarray} \tag{ 32 }$$

for which the main contribution is coming from the intra-band interaction to conserve the number of electrons in each sub-band. Then, we consider

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{{\alpha }_{1}={\alpha }_{2}}{U}_{{\rm{S}}{\alpha }_{1}}^{* }({\bf{k}}){U}_{{\rm{S}}{\alpha }_{2}}({\bf{k}}) & = & \displaystyle \sum _{{\alpha }_{1}}{\left(U{\left({\bf{k}}\right)}^{-1}\right)}_{{\alpha }_{1}{\rm{S}}}{\left(U({\bf{k}})\right)}_{{\rm{S}}{\alpha }_{1}}\\ & = & {\left(U({\bf{k}})U{\left({\bf{k}}\right)}^{-1}\right)}_{\mathrm{SS}}\\ & = & 1.\end{array}\end{eqnarray} \tag{ 33 }$$

Similarly, we confirm

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{{\alpha }_{3}={\alpha }_{4}}{U}_{{\rm{S}}{\alpha }_{3}}^{* }\left(-{k}^{{\prime} }\right){U}_{{\rm{S}}{\alpha }_{4}}\left(-{k}^{{\prime} }\right) & = & \displaystyle \sum _{{\alpha }_{3}}{\left(U{\left(-{k}^{{\prime} }\right)}^{-1}\right)}_{{\alpha }_{3}{\rm{S}}}{\left(U\left(-{k}^{{\prime} }\right)\right)}_{{\rm{S}}{\alpha }_{3}}\\ & = & {\left(U\left(-{k}^{{\prime} }\right)U{\left(-{k}^{{\prime} }\right)}^{-1}\right)}_{\mathrm{SS}}\\ & = & 1.\end{array}\end{eqnarray} \tag{ 34 }$$

Therefore, we obtain

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{HF}}\approx -\displaystyle \frac{U}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}{k}^{{\prime} }}\ \displaystyle \sum _{\alpha }{\hat{\psi }}_{{\bf{k}}\uparrow }^{\alpha \dagger }{\hat{\psi }}_{{\bf{k}}\uparrow }^{\alpha }{\hat{\psi }}_{-{k}^{{\prime} }\downarrow }^{\alpha \dagger }{\hat{\psi }}_{-{k}^{{\prime} }\downarrow }^{\alpha },\end{eqnarray} \tag{ 35 }$$

and the unitary matrix U_{S α}(k) does not contribute to expectation values.

The unitary transformation does not change the total number of electrons

$$\begin{eqnarray}&&\hat{N}=\displaystyle \sum _{{\bf{k}}\sigma {\rm{S}}}\ {\hat{\psi }}_{{\bf{k}}\sigma }^{{\rm{S}}\dagger }{\hat{\psi }}_{{\bf{k}}\sigma }^{{\rm{S}}}=\displaystyle \sum _{{\bf{k}}\sigma \alpha }\ {\hat{\psi }}_{{\bf{k}}\sigma }^{\alpha \dagger }{\hat{\psi }}_{{\bf{k}}\sigma }^{\alpha },\end{eqnarray} \tag{ 36 }$$

either. We can also define the average number of electrons per atom for each sub-lattice

$$\begin{eqnarray}&&{n}_{{\rm{S}}}=\displaystyle \frac{1}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}\sigma }\ \langle {\hat{\psi }}_{{\bf{k}}\sigma }^{{\rm{S}}\dagger }{\hat{\psi }}_{{\bf{k}}\sigma }^{{\rm{S}}}\rangle ,\end{eqnarray} \tag{ 37 }$$

which is n_S ∈ (0,2) and the filling fraction is defined by n_S/2∈ (0,1). The average total number of electrons per site among all 4 sub-lattices is

$$\begin{eqnarray}&&n=\displaystyle \frac{1}{4}\displaystyle \sum _{{\rm{S}}={\rm{A}},{\rm{B}},{\rm{C}},{\rm{D}}}{n}_{{\rm{S}}}.\end{eqnarray} \tag{ 38 }$$

At zero temperature, T = 0, we use the BCS wavefunction [3, 4, 6, 7, 61]

$$\begin{eqnarray}&&| \mathrm{BCS}\rangle =\displaystyle \prod _{{\bf{k}}\alpha }\left[{u}_{{\bf{k}}}^{\alpha }+{v}_{{\bf{k}}}^{\alpha }{\hat{\psi }}_{{\bf{k}}\uparrow }^{\alpha \dagger }{\hat{\psi }}_{-{\bf{k}}\downarrow }^{\alpha \dagger }\right]| 0\rangle ,\end{eqnarray} \tag{ 39 }$$

where ${u}_{{\bf{k}}}^{\alpha }$ and ${v}_{{\bf{k}}}^{\alpha }$ are described as

$$\begin{eqnarray}&&{u}_{{\bf{k}}}^{\alpha }=\cos \left(\displaystyle \frac{{\theta }_{{\bf{k}}}^{\alpha }}{2}\right)\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{v}_{{\bf{k}}}^{\alpha }={{\rm{e}}}^{i\phi }\sin \left(\displaystyle \frac{{\theta }_{{\bf{k}}}^{\alpha }}{2}\right)\end{eqnarray} \tag{ 41 }$$

with the superconducting U(1) phase, ϕ, to satisfy the normalisation condition

$$\begin{eqnarray}&&| {u}_{{\bf{k}}}^{\alpha }{| }^{2}+| {v}_{{\bf{k}}}^{\alpha }{| }^{2}=1\end{eqnarray} \tag{ 42 }$$

for $\langle \mathrm{BCS}| \mathrm{BCS}\rangle =1$. Then, we can calculate the thermodynamic potential energy [3, 4, 6, 7, 61, 72]

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{0}={{\rm{\Omega }}}_{0}(T=0,V,\mu )\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&=\langle {\hat{H}}_{\mathrm{kin}}+{\hat{H}}_{\mathrm{chem}}+{\hat{H}}_{\mathrm{int}}\rangle \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}\alpha }{\xi }_{{\bf{k}}\alpha }^{0}\left(1-{\theta }_{{\bf{k}}}^{\alpha }\right)-U\displaystyle \sum _{\alpha }{\left(\displaystyle \frac{{n}_{\alpha }}{2}\right)}^{2}-\displaystyle \frac{1}{U}\displaystyle \sum _{\alpha }| {{\rm{\Delta }}}_{\alpha }{| }^{2},\end{eqnarray} \tag{ 45 }$$

where ${\xi }_{{\bf{k}}\alpha }^{0}={\epsilon }_{{\bf{k}}\alpha }-\mu$. The electron density, n_α, and the superconducting energy gap, Δ_α, for the band α are self-consistently determined by the density equation

$$\begin{eqnarray}&&{n}_{\alpha }=\displaystyle \frac{1}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}}\left(1-\cos {\theta }_{{\bf{k}}}^{\alpha }\right),\end{eqnarray} \tag{ 46 }$$

and the gap equation

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\alpha }={{\rm{e}}}^{i\phi }\displaystyle \frac{U}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}}\displaystyle \frac{1}{2}\sin {\theta }_{{\bf{k}}}^{\alpha }.\end{eqnarray} \tag{ 47 }$$

By calculating the functional derivative of Ω₀, the saddle-point solution is given by

$$\begin{eqnarray}&&\displaystyle \frac{\delta {{\rm{\Omega }}}_{0}}{\delta {\theta }_{{\bf{k}}}^{\alpha }}=0,\end{eqnarray} \tag{ 48 }$$

which yields

$$\begin{eqnarray}&&{\theta }_{{\bf{k}}}^{\alpha }={\tan }^{-1}\left(\displaystyle \frac{| {{\rm{\Delta }}}_{\alpha }| }{{\xi }_{{\bf{k}}\alpha }}\right),\end{eqnarray} \tag{ 49 }$$

where ${\xi }_{{\bf{k}}\alpha }={\xi }_{{\bf{k}}\alpha }^{0}-{{Un}}_{\alpha }/2$. By using the gapped quasi-particle energy spectrum

$$\begin{eqnarray}&&{E}_{{\bf{k}}\alpha }=\sqrt{{\xi }_{{\bf{k}}\alpha }^{2}+| {{\rm{\Delta }}}_{\alpha }{| }^{2}},\end{eqnarray} \tag{ 50 }$$

the density and gap equations become

$$\begin{eqnarray}&&{n}_{\alpha }=\displaystyle \frac{1}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}}\left(1-\displaystyle \frac{{\xi }_{{\bf{k}}\alpha }}{{E}_{{\bf{k}}\alpha }}\right)\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}&&| {{\rm{\Delta }}}_{\alpha }| =\displaystyle \frac{U}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}}\displaystyle \frac{| {{\rm{\Delta }}}_{\alpha }| }{2{E}_{{\bf{k}}\alpha }},\end{eqnarray} \tag{ 52 }$$

respectively, and we obtain the saddle-point free energy

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{0}=\displaystyle \frac{1}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}}{\xi }_{{\bf{k}}\alpha }^{0}\left(1-\displaystyle \frac{{\xi }_{{\bf{k}}\alpha }}{{E}_{{\bf{k}}\alpha }}\right)-U\displaystyle \sum _{\alpha }{\left(\displaystyle \frac{{n}_{\alpha }}{2}\right)}^{2}-\displaystyle \frac{1}{U}\displaystyle \sum _{\alpha }| {{\rm{\Delta }}}_{\alpha }{| }^{2}.\end{eqnarray} \tag{ 53 }$$

The gap equation always has a normal state solution of $| {{\rm{\Delta }}}_{\alpha }| =0$, while a superconducting solution has none-zero gap energy, $| {{\rm{\Delta }}}_{\alpha }| \ne 0$. It is not always obvious whether a superconducting solutions has a lower energy, as discovered in Graphene [67–69], due to vanishing DOS at Dirac points. In the grand canonical ensemble with variable number of electrons, the state with the lowest thermodynamic potential energy Ω₀(T = 0, V, μ) will be realised at the given chemical potential, μ. For the fixed number of electrons, N_e, or equivalently, at the fixed density, $n={N}_{{\rm{e}}}/{{ \mathcal N }}_{\mathrm{tot}}$, on the other hand, we need to find out a state with the lowest Helmholtz free energy, whose zero-temperature form is given by the internal energy

$$\begin{eqnarray}&&{F}_{0}={F}_{0}(T=0,V,n)\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&=\ \langle {\hat{H}}_{\mathrm{kin}}+{\hat{H}}_{\mathrm{int}}\rangle .\end{eqnarray} \tag{ 55 }$$

The difference between Ω₀(T = 0, V, μ) and F₀(T = 0, V, n) was crucial to identify the correct ground state.

At finite temperature, $T\ne 0$, we can use a standard Bogoliubov transformations [3, 4, 6, 7, 61, 72] to obtain the self-consistency conditions

$$\begin{eqnarray}&&{n}_{\alpha }=\displaystyle \frac{1}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}}\left(1-\displaystyle \frac{{\xi }_{{\bf{k}}\alpha }}{{E}_{{\bf{k}}\alpha }}\tanh \left(\displaystyle \frac{\beta {E}_{{\bf{k}}\alpha }}{2}\right)\right)\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&| {{\rm{\Delta }}}_{\alpha }| =\displaystyle \frac{U}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}}\displaystyle \frac{| {{\rm{\Delta }}}_{\alpha }| }{2{E}_{{\bf{k}}\alpha }}\tanh \left(\displaystyle \frac{\beta {E}_{{\bf{k}}\alpha }}{2}\right),\end{eqnarray} \tag{ 57 }$$

where β = 1/(k_BT) with the Boltzmann constant of k_B, and the thermodynamic potential becomes

$$\begin{eqnarray}&&{\rm{\Omega }}={\rm{\Omega }}(T,V,\mu )\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{2}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}}{E}_{{\bf{k}}\alpha }\displaystyle \frac{1}{{{\rm{e}}}^{\beta {E}_{{\bf{k}}\alpha }}+1}-\displaystyle \frac{1}{{ \mathcal N }}\displaystyle \sum _{{\bf{k}}}\left({E}_{{\bf{k}}\alpha }-{\xi }_{{\bf{k}}\alpha }\right)+U\displaystyle \sum _{\alpha }{\left(\displaystyle \frac{{n}_{\alpha }}{2}\right)}^{2}+\displaystyle \frac{1}{U}\displaystyle \sum _{\alpha }| {{\rm{\Delta }}}_{\alpha }{| }^{2},\end{eqnarray} \tag{ 59 }$$

whose first term vanishes towards T = 0, corresponding to the disappearance of a quasi-particle, which is the main reason for zero resistivity and the Meissner effect, because the absence of single particle excitation leads to the macroscopic coherence, characterised by the current proportional to the vector potential [3–6].

Numerical solutions of self-consistency equations at T = 0 are shown in figure 9. The non-doped (2, 2)-Hopfene is half-filled (n = 1), and it is semi-metal at U = 0 (figures 8 and 9(a)). With the introduction of $U\ne 0$, Δ₃ developed with the peak at U ∼ 3 eV (figure 9(b)) and it decreased upon increasing U. This was attributed to the reduction of the filling fraction of n₃ due to the decrease of the effective lower band energy of ξ_k2 = _k2 − Un₂/2 by the Hartree–Fock interaction. Consequently, the _k3 band becomes empty at U > U_MI ∼ 4 eV, and Δ₃ vanishes, and we expect a Metal-Insulator transition in the normal state at U_MI due to the opening of the semiconductor gap. With sufficient doping, such significant changes in the filling fractions were not found in a system (figure 9(b)). On the other hand, Δ₂ remained to be zero at n = 1 until U_c ∼ 1 eV, showing that the system is very close to the quantum critical point [67–69]. At U < U_c, the free energy is lower by reducing ξ_k2 without opening Δ₂, which is similar to the finding of the quantum critical behaviour in 2D Graphene. However, in Hopfene, the DOS is finite at μ ∼ 1 eV (figures 8 and 9(e)) and the reduction in the mean-field Hartree–Fock energy overcame the emergence of superconducting order. Surprisingly, further increase of U > U_MI allowed the increase of Δ₂, although the increase is significantly smaller than that in the doped case at n = 1.25 (figure 9(d)). This increase of Δ₂ at n = 1 is similar to the pioneering work found in [65], where the superconducting order is expected in a semiconductor with a small band gap, if the attractive interaction between electrons in conduction band and holes in valence band is strong enough to overcome the semiconductor gap. An excitonic strong pairing and condensation were experimentally found in Graphite under extremely high magnetic fields in the quantum limit [66]. In our case, the attractive interaction is working within the sub-bands, and the opening of the superconducting gap Δ₂ is smaller than the semiconductor gap ∼U. Consequently, the energy gain by the superconducting order at n = 1 was significantly limited.

Zoom In Zoom Out Reset image size

The condensation energies were calculated by the difference of the Helmholtz free energy between superconducting and normal state, Δ F = F_super − F_normal (figure 10). At n = 1, we estimated solely from the BCS theory (figure 10(a)), because the band is almost completely filled and the impact of strong correlation is limited. In fact, Δ F was order of magnitude smaller than that in the doped case (figure 10(b)), so that the system at n = 1 is closer to a band insulator. In particular, Δ F ∼ 0 at U ∼ U_MI, where Δ₃ is also close to zero. On the other hand, the BCS theory overestimates Δ F for n = 1.25 at the strong coupling limit. It is beyond the scope of this paper to include the impact of strong correlation precisely, because it is quite challenging to include small energy scale of the order t²/U, which is similar to the problem of Mott transition [61]. Here, we roughly estimate the impact of strong correlation by considering a Gutzwiller-BCS variational wavefunction [64, 73, 74]

$$\begin{eqnarray}&&| \mathrm{GBCS}\rangle =\displaystyle \prod _{i}\left[1+(g-1){n}_{i\uparrow }{n}_{i\downarrow }\right]| \mathrm{BCS}\rangle ,\end{eqnarray} \tag{ 60 }$$

by which local spatial correlation is included in the limit of $g\to \infty$ [64, 74]. In this limit, the electrons are spreading completely in the momentum space to form local singlet of preformed electrons. The impact of the correlation on the free energy in a superconducting state is actually not significant, while the change in the normal state is significant [64]. Therefore, we estimated the normal state free energy simply by assuming all electrons to form locally singlet pairs with the interaction energy of −Un/2, while the kinetic energy was calculated by assuming local paring corresponds to the spreading in the entire sub-band in momentum space [64, 74]. The calculated condensation energy is considered to be a simplified BEC model, which is shown in figure 10(b). We could somehow include the energy scale of t²/U, and the true condensation energy would show a crossover from BCS to BEC.

Zoom In Zoom Out Reset image size

We have also examined the filling dependence of superconductivity in (2,2)-Hopfene (figure 11). We found significant decreases of superconducting order parameters at every quarterly fillings (figures 11(a), (c), and (e)), because of the complete filling of the lower sub-bands. These observations are consistent with the filling fractions found in figures 11(b), (d), and (f), which show the subsequent completion of filling at n = 0.5, 1, 1.5 and 2, so that the system is close to quantum critical points towards a semiconductor or an insulator. At n = 1, it was quite marginal whether the superconducting order remains, as we have seen in figures 9 and 10, and according to our numerical calculations, the system remains to be superconducting at T = 0. The robust superconductivity at n = 1 was coming from the sharp peak of the sub-band ₃ at the bottom edge of the band, because of the flat-band (figure 5(c)). At n = 0.5, the superconducting order parameters are reduced even further, which was coming from the stronger confinement of 3D Weyl Fermions at the top edge of the sub-band ₁ (figure 5(d)). Therefore, the system at n = 0.5 in the strong coupling is very close to the situation of the reduced superconductivity found in 2D Dirac Fermions [67–69]. In our case, the system is semi-metal with the finite DOS from sub-bands ₁ and ₂ at the weak coupling, so that the superconducting order is still finite and there is no finite critical attractive interaction strength to induce the superconductivity. In that sense, from fundamental point of view, superconductivity in 3D Hopfene is less significant compared with the case in 2D Graphene. This is consistent with the common understanding, that order parameters are more robust against quantum and thermal fluctuations in higher dimensions, where mean-field pictures are validated [11, 13, 35, 37]. Nevertheless, at strong coupling, the Hartree–Fock interaction reduces the sub-band energy, effectively, leading to an isolation of ₁ and opening of a semiconductor gap.

Zoom In Zoom Out Reset image size

The superconducting energy gaps were significantly affected by DOS (figures 11(e) and 8) at weak coupling, as expected from the BCS theory. On the other hand, in the strong coupling, the gaps are smooth functions of fillings (figures 11(a) and (c)), showing the states are spreading to the entire Brillouin zone of the momentum space in order to form local singlet pairs, so that the system is close to BEC.

We have also examined the temperature dependence of superconducting order parameters by solving the self-consistency conditions, equations (56) and (57) simultaneously (Figure 12). We found continuous reductions of order parameters as we increase the temperature, so that the phase transitions are of the second order. This is expected in a 3D material, for which the mean field theory should be appropriate. Superconductivity is considerably reduced at n = 1, while it is order of magnitude larger at n = 1.25, consistent with previous results. Obviously, the superconducting transition temperatures are too large if we assume U > 2 eV. In practice, superconductivity in carbon allotropes [15, 16, 56, 57, 60] were found only at temperatures much lower than room temperature. Therefore, a weak coupling condition U < 1 eV is appropriate, when we consider realistic experiments. Nevertheless, it would be valuable to consider a crossover [61] or a possible quantum phase transition [19, 64] at least theoretically. Stronger coupling would be expected to excitonic systems [66], for which the model needs to be modified to allow the inter-band interactions.

Zoom In Zoom Out Reset image size

Finally, we have made phase diagrams of (2, 2)-Hopfene based on an attractive Hubbard model (figure 13). At n = 1, which would be expected for non-doped Hopfene, we expect BCS superconductivity with reduced transition temperatures. In the weak coupling limit, only Δ₃ opens up. At moderate interaction strengths, 1 eV < U < 4 eV, both Δ₂ and Δ₃ open at low temperature. At strong coupling limit U > 4 eV, we expect insulator to superconductor transition at high temperature, because of the opening of a semiconductor gap, solely from the attractive Hartree–Fock interaction, which is different from Mott transition coming from strong correlation. Obviously, the transition temperature of the order 1000 K is not realistic at all, since Hopfene would be burned out. Nevertheless, we hope there is some theoretical implications of this result, since our results suggest competitions of a semiconductor band gap and a superconducting energy gap with an increase of attractive interaction strength, which is similar to a problem suggested by Nozières and Pistolesi [65]. At n = 1.25 with sufficient doping, we confirmed a smooth crossover from BCS to BEC (figure 13(b)). The transition temperature of BEC was simply estimated from the condensation energy, so that it is not accurate quantitatively. Nevertheless, the result suggests the existence of the optimum interaction strength at around U ∼ 2 eV, so that the infinitely strong attraction would not help to increase transition temperature, which is consistent with previous publications [61]. The gap formation is important at weak coupling to expect superconductivity, while the phase coherence is required to expect BEC of local singlet pairs. The sharp increase in DOS by a flat-band will also help to increase the transition temperature at practical weak coupling, while it is less effective due to the spreading in momentum space at strong coupling. The impact of doping on superconductivity in Hopfene is significant and the system is located to the close proximity of the quantum critical point, where the transition temperature is reduced to zero, due to the competition between the opening of semiconductor and superconducting gaps.

Zoom In Zoom Out Reset image size