Skip to main content

Thank you for visiting You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

Pairing symmetries in the Zeeman-coupled extended attractive Hubbard model


By introducing the possibility of equal- and opposite-spin pairings concurrently, we show that the ground state of the extended attractive Hubbard model (EAHM) exhibits rich phase diagrams with a variety of singlet, triplet, and mixed parity superconducting orders. We study the competition between these superconducting pairing symmetries invoking an unrestricted Hartree–Fock–Bogoliubov–de Gennes (HFBdG) mean-field approach, and we use the d-vector formalism to characterize the nature of the stabilized superconducting orders. We discover that, while all other types of orders are suppressed, a non-unitary triplet order dominates the phase space in the presence of an in-plane external magnetic field. We also find a transition between a non-unitary to unitary superconducting phase driven by the change in average electron density. Our results serve as a reference for identifying and understanding the nature of superconductivity based on the symmetries of the pairing correlations. The results further highlight that EAHM is a suitable effective model for describing most of the pairing symmetries discovered in different materials.


In a conventional s-wave superconductor, described by BCS–Migdal–Eliashberg theory1,2,3,4,5,6, a phonon-mediated effective attraction causes the electrons to form spin-singlet (with total spin \(S=0\)) Cooper pairs with an isotropic s-wave orbital order parameter (OP) symmetry7. Most of the superconductors discovered in the early phase of the past century have this type of conventional OP symmetry8. On the contrary, the general consensus about the novel high \(T_{c}\) cuprate superconductors is that they are identified as strong candidates for unconventional d-wave superconductors9,10,11, which support the formation of spin-singlet Cooper pairs with an anisotropic d-wave orbital OP symmetry. In general, Cooper pairs can also be formed in the spin-triplet state, with total spin \(S=1\) and anisotropic orbital OP12,13. Recent studies revealed that a two dimensional chiral p-wave spin-triplet superconductor is a candidate to host Majorana fermions14,15,16. A pair of Majorana fermions, bound to topological defects, together known as Ising anyons, exhibit non-Abelian exchange statistics and such an object can be considered to be the potential building block for decoherence free quantum computation17. Regardless of their technological relevance, superconductors with triplet-superconductivity emerging due to the intrinsic properties of the material itself are quite rare in nature. While \(^{3}\text {He}\) was the first charge-less many body system where triplet-pairing state was realised18,19,20,21, it was \(\text {UPt}_{3}\) which was identified as the first spin-triplet superconductor in charged many body systems22,23. Strong candidate materials for triplet-superconductivity, so far, include Uranium based heavy-fermion superconductors \(\text {UPt}_{3}\), \(\text {UGe}_{2}\), URhGe, UCoGe24, organic superconductor \(\text {(TMTSF)}_{2}\text {PF}_{6}\)25 and probably recently discovered heavy-fermion superconductor \(\text {UTe}_{2}\)26,27,28. We would also like to add that Strontium Ruthenate (Sr2RuO4), which for its most part was understood as a strong candidate for odd-parity superconductor, has undergone immense scrutiny. The exact nature of the SC order in Sr2RuO4 has remained an open question for many years now, with a recent study ruling out the existence of a pure odd-parity superconducting order, and explaining that the Knight shifts originate from quasiparticles rather than the condensate itself29. Other studies impose constraints on the possible order of Sr2RuO4 superconducting state—stress-induced splitting of onset temperatures for SC and time-reversal symmetry (TRS) breaking has also been reported30, and ultrasound measurements suggests the existence of a multi-component order parameter31.

In conventional superconductors TRS allows energetically degenerate electron states \(|\mathbf {k}\uparrow \rangle\) and \(|-\mathbf {k}\downarrow \rangle\) to form spin-singlet Cooper pairs. On the other hand, parity symmetry (inversion symmetry in this scenario) along with TRS is needed to form odd-parity spin-triplet Cooper pairs32,33,34,35. In case of broken inversion symmetry, there is also the possibility of forming Cooper pairs with mixed-parity superconducting (SC) states32,35,36,37. Possible occurrence of mixed parity SC states have been reported in several theoretical studies on different contexts, including quantized vortices of superfluid \(^{3}\text {He}\)38, helical mixed parity SC phase relevant to the unconventional SC phases of \(\text {UPt}_{3}\)39, Larkin–Ovchinnikov–Fulde–Ferrell phase with inhomogeneous order parameter40,41, two-dimensional superconductors with broken inversion symmetry36, noncentrosymmetric superconductors42, vortex phase of a d-wave superconductor in presence of paramagnetic effects43, SC topological insulators in presence of surface Dirac fermions44, and disordered monolayer transition metal dichalcogenides45. While the search for unconventional pairing mechanism has not been an easy path, understanding unconventional pairing symmetries of the OP has been of utmost importance since it certainly gives us clues about the possible pairing mechanisms. Therefore, understanding and characterizing the unconventional OP symmetries of different spin-triplet and mixed-parity SC states remains a problem of crucial importance.

Most interestingly, the response of a superconducting state to an external magnetic field can reveal significant details about the pairing state of the Cooper pairs. For example, unlike the fully gapped s-wave superconductors, \(d_{x^2-y^2}\)-wave superconductors exhibit gapless quasiparticle spectrum, which allows even a weak amount of Zeeman field to spin-polarize the quasiparticle excitations, resulting in a stronger response to an external magnetic field46,47,48. For a two dimensional system it is quite interesting to study the effect of an in-plane (parallel to the plane) magnetic field, as the orbital motions of the electrons do not couple with the external magnetic field at all, and the Zeeman coupling of the electronic spins to the magnetic field is enough to study the response of the superconducting state. The effect of an in-plane magnetic field has been studied in various contexts, including quasi-2D systems like cuprates, emergence of a dissipative state in a superconducting \(\text {Mo}_{0.79}\text {Ge}_{0.21}\) nanostrip, anisotropy of the upper critical field in \(\text {Sr}_{2}\text {RuO}_{4}\), and magnetic field driven nodal topological superconductivity in monolayer transition metal dichalcogenides49,50.

In this article we provide a detailed study of the competition between different superconducting pairings and their response to non-zero temperature and Zeeman field. We consider an extended attractive Hubbard model, with on-site and nearest neighbour attractive interaction on a square lattice, and treat the many-body interaction terms using mean-field analysis. We decouple the interaction terms into pairing channels without imposing any symmetry constraint on the superconducting pairing correlations. By doing so, it allows one to compare different superconducting ground states with different symmetries and spin pairings, i.e. a singlet and three components of triplet states. We employ the d-vector definition of the triplet order parameter to keep track of the triplet pairings throughout our analysis. We demonstrate that in this model numerous superconducting states with various symmetries can be stabilized and studied. Therefore, the extended attractive Hubbard model can be a universal effective model for studying unconventional, particularly triplet, superconductivity. The rest of the article is organized as follows. We start by presenting the model and defining all the parameters involved, and subsequently discussing the results: In “Determination of relative phase angle between different pairing correlations” section, using energy minimization arguments we self-consistently reduce the exploration parameter space and fix certain relative phase angles between different pairing correlations. Followed by “Competition among the different triplet states” section, in which we explore the competition between the three kinds of triplet pairings at non-zero magnetic field to lift the spin degeneracy. In “Ground State phase diagrams” section, we present the ground state phase diagram of nearest-neighbour attractive strength versus particle density. Next in “Finite temperature study” and “Effect of Zeeman coupling” sections, we study the behaviour of different superconducting states at non-zero temperatures and Zeeman field respectively. Here we will encounter how different triplet and singlet orders react to the Zeeman field. Then we do a further characterization of superconducting phases that are stabilized in our study using symmetry and d-vector approach in “Characterization of phases” section. Finally, we conclude by summarising our findings. All the methods we employ in this work, to analyse the model, are presented in the “Methods” section at the end of the article. In the same section, we also discuss the d-vector formalism and its usefulness in characterizing triplet order parameters.

Results and discussion

For superconductors like high-temperature Cuprates, the electronic band is quite narrow i.e., the orbitals have a small overlap between adjacent atoms. We therefore make use of the tight-binding Hamiltonian, which is either constructed from atomic orbitals or from Wannier orbitals, to study narrow band behaviours arising from electronic correlation effects. We, therefore, consider the Extended Attractive Hubbard Model, defined on a square lattice to be the doorway to address our problem. In real space, the Hamiltonian corresponding to the Extended Attractive Hubbard Model on a two-dimensional square lattice can be written as,

$$\begin{aligned} \mathcal {H}= \mathcal {H}_{\text {kin}} + \mathcal {H}_{\mu } + \mathcal {H}_{\text {int}}^{\text {onsite}} + \mathcal {H}_{\text {int}}^{\text {nn}} + \mathcal {H}_{\text {B}}, \end{aligned}$$

where, \(\mathcal {H}_{\text {kin}}\) represents the kinetic energy of the electrons, \(\mathcal {H}_{\mu }\) is the chemical potential term, \(\mathcal {H}_{\text {int}}^{\text {onsite}}\) is the on-site attractive interaction between two different spin projections, \(\mathcal {H}_{\text {int}}^{\text {nn}}\) is the nearest-neighbor (nn) attractive interaction which would be responsible for all the unconventional superconducting phases, and finally \(\mathcal {H}_{\text {B}}\) contains the Zeeman coupling of electron spin to an external magnetic field B. The corresponding operator terms, written in the second quantization notation are given as,

$$\begin{aligned} \mathcal {H}_{\text {kin}}&= -t\sum \limits _{\langle ij\rangle ,\sigma }c_{i\sigma }^{\dagger }c^{}_{j\sigma } + \text {H.c.}, \end{aligned}$$
$$\begin{aligned} \mathcal {H}_{\mu }&= -\mu \sum \limits _{i\sigma }n_{i\sigma } = -\mu \sum \limits _{i\sigma }c_{i\sigma }^{\dagger }c^{}_{i\sigma }, \end{aligned}$$
$$\begin{aligned} \mathcal {H}_{\text {int}}^{\text {onsite}}&= -U\sum \limits _{i}n_{i\uparrow }n_{i\downarrow } = -U\sum \limits _{i}c_{i\uparrow }^{\dagger }c^{}_{i\uparrow } c_{i\downarrow }^{\dagger } c^{}_{i\downarrow }, \end{aligned}$$
$$\begin{aligned} \mathcal {H}_{\text {int}}^{\text {nn}}&= -V\sum \limits _{\langle ij\rangle }n_{i}n_{j} = -V\sum \limits _{{\mathop {\sigma ,\sigma ^{\prime }}\limits ^{\langle ij\rangle }}} c_{i\sigma }^{\dagger }c^{}_{i\sigma } c_{j\sigma ^{\prime }}^{\dagger }c^{}_{j\sigma ^{\prime }}, \end{aligned}$$
$$\begin{aligned} \mathcal {H}_{\text {B}}&= -B\sum \limits _{i}(n_{i\uparrow }-n_{i\downarrow }) = -B\sum \limits _{i}(c_{i\uparrow }^{\dagger }c^{}_{i\uparrow } - c_{i\downarrow }^{\dagger }c^{}_{i\downarrow }). \end{aligned}$$

In the above, \(c_{i\sigma }^{\dagger }\)(\(c_{i\sigma }\)) creates (annihilates) an electron with spin \(\sigma\) at \(i{\text {th}}\) site of the lattice. The sum over \(\langle ij\rangle\) represents the sum over the nearest-neighbor sites of the square lattice. The operator \(n_{i\sigma }\) is the electron occupation number operator at \(i{\text {th}}\) site for \(\sigma\) spin-projection. The total electron occupation number at \(i{\text {th}}\) site is, therefore, defined as \(n_{i} = n_{i\uparrow } + n_{i\downarrow }\). The parameter U is the on-site attractive interaction strength, whereas V is the inter-site attractive interaction strength. The z-component of the magnetic field is given by B, which lies in the lattice plane. t is the usual nearest-neighbor hopping amplitude and finally, \(\mu\) is the chemical potential as we work in the grand canonical ensemble.

Looking closely at the interaction terms individually and expanding them in second quantization notation we have:

$$\begin{aligned} \mathcal {H}_{\text {int}}^{\text {onsite}}&= -U\sum \limits _{i}n_{i\uparrow }n_{i\downarrow } = -U\sum \limits _{i}c_{i\uparrow }^{\dagger }c_{i\uparrow } c_{i\downarrow }^{\dagger }c_{i\downarrow } = -U\sum \limits _{i}c_{i\uparrow }^{\dagger }c_{i\downarrow }^{\dagger } c_{i\downarrow }c_{i\uparrow }, \end{aligned}$$

where we rearranged the Fermionic operators in particle–particle channel or so called pairing channel, using anti-commutation rules. Similarly for the inter-site term:

$$\begin{aligned} \mathcal {H}_{\text {int}}^{\text {nn}}&= -V\sum \limits _{\langle ij\rangle }n_{i}n_{j} = -V\sum \limits _{i,\delta }n_{i}n_{i+\delta } = -V\sum \limits _{i,\delta }(n_{i\uparrow }+n_{i\downarrow }) (n_{i+\delta \uparrow }+n_{i+\delta \downarrow }) \nonumber \\&=-V\sum \limits _{i,\delta } (n_{i\uparrow }n_{i+\delta \uparrow } +n_{i\downarrow }n_{i+\delta \downarrow } +n_{i\uparrow }n_{i+\delta \downarrow } +n_{i\downarrow }n_{i+\delta \uparrow } ) \nonumber \\&= \underbrace{-V\sum \limits _{i,\delta }c_{i\uparrow }^{\dagger }c_{i+\delta \uparrow }^{\dagger }c_{i+\delta \uparrow }c_{i\uparrow }}_A~~~ \underbrace{-V\sum \limits _{i,\delta }c_{i\downarrow }^{\dagger }c_{i+\delta \downarrow }^{\dagger } c_{i+\delta \downarrow }c_{i\downarrow }}_B ~~~ \underbrace{-V\sum \limits _{i,\delta }c_{i\uparrow }^{\dagger }c_{i+\delta \downarrow }^{\dagger } c_{i+\delta \downarrow }c_{i\uparrow }}_C \nonumber \\&\underbrace{\quad -V\sum \limits _{i,\delta }c_{i+\delta \uparrow }^{\dagger }c_{i\downarrow }^{\dagger } c_{i\downarrow }c_{i+\delta \uparrow }}_D, \end{aligned}$$

where \(\delta\) is an index which denotes unit vectors along the direction of two nearest neighbors, namely \(\delta = +{\hat{x}}\) and \(+{\hat{y}}\). Note that we have grouped terms as A, B, C, and D for the ease of reference. Now terms A and B lead to the possibility of triplet solution with \(S_{z}=\pm 1\) (ESP states), while terms C and D lead to singlet and triplet solution with \(S_{z}=0\) (OSP states).

Notice that this model describes the many-particle interaction which in itself is very hard to analyse. We therefore make the mean-field approximation and via the Bogoliubov-de Gennes framework, we solve the Hamiltonian computationally, and present the results in the following sections. The details to the BdG framework are given in “Bogoliubov-de Gennes method” section.

Determination of relative phase angle between different pairing correlations

Before we discuss the competition among different superconducting states, we consider the possibility of obtaining different triplet states from different combinations of pairing mean field parameters.

Triplet combinations from \(\Delta ^{\uparrow \downarrow }(\mathbf {k})\)

  1. 1.

    \(p_{x}\) state (and correspondingly \(p_y\) state): \(\Delta _{x}^{+}=-\Delta _{x}^{-} \in {\mathbb {C}};~~~~~~\Delta _{y}^{+}= \Delta _{y}^{-}=0\).

  2. 2.

    \(p_{x}\pm ip_{y}\) state : \(\Delta _{x}^{+}=-\Delta _{x}^{-} = C \in {\mathbb {C}}; ~~~~~ \Delta _{y}^{+}=-\Delta _{y}^{-}=iC\).

In general, \(p_{x}\) and \(p_{y}\) can be of the form: \(p_{x}+e^{i\theta }p_{y}\). With the help of energetics and self-consistency we have checked that \(\theta = \frac{\pi }{2}\left( -\frac{\pi }{2}\right)\) corresponds to the most stable solution, i.e. it takes the form \(p_{x}\pm ip_{y}\), independent of the external parameters. The energetics suggests that for all external model parameters, the most stable relative phase angle between \(p_x\) and \(p_y\) states is \(\pi /2\). This allows us to reduce our exploration space \(\{ \Delta \}\). Similarly, we fix the phase relation among, \(\Delta _{x}^+,\Delta _x^-,\Delta _y^+\) and \(\Delta _y^-\).

Triplet combinations from \(\Delta ^{\uparrow \uparrow }(\mathbf {k})\) and \(\Delta ^{\downarrow \downarrow }(\mathbf {k})\)

  1. 1.

    \(\Delta ^{\uparrow \uparrow }(k)=2iV\left[ \Delta _{x}^{\uparrow }\sin (kx)+ \Delta _{y}^{\uparrow }\sin (ky)\right]\).

  2. 2.

    \(\Delta ^{\downarrow \downarrow }(k)=2iV\left[ \Delta _{x}^{\downarrow }\sin (kx)+ \Delta _{y}^{\downarrow }\sin (ky)\right]\).

Again by comparing energies we find the phase relations between \(\Delta _{x}^{\uparrow }\) and \(\Delta _{y}^{\uparrow }\) (and correspondingly, between \(\Delta _{x}^{\downarrow }\) and \(\Delta _{y}^{\downarrow }\)) i.e., \(\Delta _{y}^{\uparrow } = e^{i\Phi _{1}}\Delta _{x}^{\uparrow }\) (and \(\Delta _{y}^{\downarrow } = e^{i\Phi _{1}}\Delta _{x}^{\downarrow }\)). We find that \(\Phi _{1} = \frac{\pi }{2}\) minimizes the average energy \(\langle E\rangle\) (Fig. 1e) independent of the external parameters. Note that we also allowed \(\uparrow\) and \(\downarrow\) parring correlation functions to have a relative phase of \(\Phi _2\), and our results are independent of this relative phase. These phase relations reduce our parameter space considerably. Though these plots (Fig. 1e) are generated on the basis of mere energetics, we also checked the validity of these relations in self-consistent solutions. Notice that these relations only hold for the correlation functions \(\{\Delta _{x}^{\uparrow }, \Delta _{y}^{\uparrow }\}\) (and \(\{\Delta _{x}^{\downarrow }, \Delta _{y}^{\downarrow }\}\)) and we cannot, at this stage, comment on the phase difference between \(\Delta ^{\uparrow \uparrow }(\mathbf {k})\) and \(\Delta ^{\downarrow \downarrow }(\mathbf {k})\) because these results are independent of \(\Phi _2\).

Figure 1
figure 1

(a–d) The plots show the variation of average energy \(\langle E\rangle\) corresponding to three triplet states with hopping amplitude t at zero temperature for different magnitudes of the magnetic field, (a) \(B=0\), (b) \(B=(0.2/1.8) V\), (c) \(B=(0.5/1.8) V\) and combinations of pairing mean \(B=(2.0/1.8) V\). (e) Plot shows the variation of average energy \(\langle E\rangle\) with relative phase angle \(\Phi _{1}\), where \(\Delta _{x}^{\uparrow } = 0.5,\Delta _{y}^{\uparrow } = 0.5e^{i\Phi _{1}}, \Delta _{x}^{\downarrow } = 0.5e^{i\Phi _{2}}\) and \(\Delta _{y}^{\downarrow } = 0.5e^{i(\Phi _1+\Phi _2)}\). The variation shown here is independent of the choice of external parameters like the chemical potential \(\mu\), hopping parameter t, magnetic field \(\mathbf {B}\), and the relative phase between \(\Delta _{\delta }^{\uparrow }\) and \(\Delta _{\delta }^{\downarrow }\) i.e., \(\Phi _2\). This suggests that, independent of the external parameters, the most stable relative phase between \(\Delta _{x}^{\sigma }\) and \(\Delta _{y}^{\sigma }\) is \(\Phi _1 = \pi /2\). (f,g) The plots shows the variation of average energy \(\langle E\rangle\) for the three triplet states with the magnetic field B, for (f) \(t=0\) and (g) \(t=(1/1.8) V\).

Competition among the different triplet states

In this section we compare the energies of the three triplet states: \(|S_{z}=0\rangle\) \(\equiv \frac{1}{\sqrt{2}}\left( |\uparrow \downarrow \rangle + |\downarrow \uparrow \rangle \right)\), \(|S_{z}=1\rangle\) \(\equiv |\uparrow \uparrow \rangle\), and \(|S_{z}=-1\rangle\) \(\equiv |\downarrow \downarrow \rangle\). We expect all three states \(|S_{z}=0\rangle\), \(|S_{z}=1\rangle\), and \(|S_{z}=-1\rangle\) to be degenerate in the absence of magnetic field. Furthermore, we expect this degeneracy to be lifted in the presence of magnetic field, and one of the equal spin pairing (ESP) states \(|S_{z}=1\rangle\) or \(|S_{z}=-1\rangle\) will have lower energy in the presence of an in-plane magnetic field depending on its orientation. In our formulation, singlet and \(|S_{z}=0\rangle\) triplet component enter into the Hamiltonian through the superconducting pairing correlation \(\Delta ^{\uparrow \downarrow }(\mathbf {k})\); and \(|S_{z}=1\rangle\), \(|S_{z}=-1\rangle\) triplet contributions enter through the triplet superconducting correlations \(\Delta ^{\uparrow \uparrow }(\mathbf {k})\) and \(\Delta ^{\downarrow \downarrow }(\mathbf {k})\). Therefore, we expect this behaviour to translate to the superconducting pairing correlations as well.

The analysis is summarised in Fig. 1a–d,f–g. In Fig. 1a–d, we plot the variation of average system energy corresponding to all possible triplet states (both OSP and ESP) with increasing amplitude of nearest-neighbor hopping parameter t. It is to be noted that, in Fig. 1, all energies are measured in units of V/1.8. This slightly unusual choice of basic energy scale is unavoidable as in Fig. 1 we are trying to connect the discussion to the atomic limit by presenting the effect of variations in the hopping parameter. On the other hand, for the rest of the article, we persistently use \(t=1\) to be the unit for the remaining parameters. As we focus on the energies of the triplet states, we pick a specific point in the parameter space (\(U=0\), \(\mu = -(1.5/1.8)V\)) where we found pure-triplet (\(p_{x}\pm ip_{y}\) type) state being stabilized in our previous study51. As the interaction strength V is kept fixed while increasing nearest-neighbor hopping amplitude t, effectively we are moving from a strong-coupling limit to a weak-coupling limit.

In the strong-coupling limit, where the electron-pairs are tightly bound, we can appropriately apply the well-known physics of spin-singlet and spin-triplet pairing to the spin degrees of freedom associated with the electron-pairs. In this limit, absence of magnetic field leads to degenerate states. In fact, absence of magnetic field provides the freedom to transform ESP states into specific OSP states, via a suitable choice of quantization axis due to rotational symmetry. In the presence of a magnetic field this symmetry is broken and therefore the degeneracy is lifted. In the weak-coupling limit, the system gains energy via de-localization of electrons, and characterization of different triplet states in terms of local spin operators is not valid. Note that we use the phase-lock between different pairing correlations obtained in the previous section—this reduces the exploration parameter space significantly and gives us the freedom to use only magnitudes of different correlation functions. In Fig. 1a, at \(t=0\), we choose these magnitudes in a manner so that all the triplet states, both OSP and ESP, remain degenerate. With the increase of nearest-neighbour hopping amplitude t we note that the OSP state (\(\uparrow \downarrow\) in figure) becomes energetically favourable. Furthermore, at any finite value of Zeeman-coupling, \(\uparrow \uparrow\)-type ESP state remain energetically favourable as we expected, while the energy of \(\downarrow \downarrow\)-type ESP state in comparison to OSP state behave differently in different coupling regimes. At moderate values of Zeeman-coupling, the system prefers OSP state to \(\downarrow \downarrow\)-type ESP state in both weak and strong coupling limits (Fig.1b,c). However, at a larger value of Zeeman-coupling, OSP state is completely disfavoured in the weak-coupling limit (Fig.1d). One should take a note of the fact that this particular observation is validated by self-consistent solutions presented in the forthcoming sections.

Figure 1f,g completes this analysis, as it captures the variation of energy of the triplet states with increasing values of Zeeman-coupling in the two limiting cases (\(t = 0\) and \(t=(1/1.8) V\)). At strong-coupling (\(t=0\)) limit, OSP state is of lower energy compared to \(\downarrow \downarrow\)-type ESP state below a critical Zeeman-coupling (\(B\approx (2/1.8) V\)). Beyond such critical value of Zeeman-coupling OSP state becomes completely disfavoured energetically with respect to the ESP states. On the other hand, at weak-coupling limit (\(t=(1/1.8) V\)), energy of OSP state remains comparable with \(\downarrow \downarrow\)-type ESP state at moderate values of Zeeman-coupling, and eventually becomes energetically disfavoured at higher values of Zeeman-coupling. On one hand this section provides us with the idea of how OSP and ESP states behave at two different coupling regimes. On the other hand, it provides us an insight into the nature of triplet states and which type is likely to be stabilized at different values of Zeeman-coupling.

Ground state phase diagrams

In the previous work51 we explored the possibility of stabilizing OSP superconducting solutions in the framework of extended attractive Hubbard model defined on a square lattice in the absence of a magnetic field. The results were based on the set of order parameters {\(\Delta\)} = { \(\Delta _{s},\Delta _{s^{*}},\Delta _{d_{x^2-y^2}} \Delta _{p_{x}},\Delta _{p_{y}}\) }, which were obtained using the standard definition of superconducting OPs, defined in literature51, from the original set of pairing correlations { \(\Delta , \Delta _{\delta }^{+}, \Delta _{\delta }^{-}\) }. Explicitly writing {\(\Delta\)} we have,

$$\begin{aligned} \Delta _{s}&= \Delta ,\nonumber \\ \Delta _{d_{x^{2}-y^{2}}/s^{*}}&= \left[ (\Delta _{x}^{+}+\Delta _{x}^{-}) \mp (\Delta _{y}^{+}+\Delta _{y}^{-}) \right] /4 ,\nonumber \\ \Delta _{p_{x}/p_{y}}&= \left[ \Delta _{x/y}^{+} - \Delta _{x/y}^{-} \right] /2 . \end{aligned}$$

Note that \(\delta\) is an index for the unit vectors \(\delta = {\hat{x}},{\hat{y}}\) along two independent directions on a square lattice. So { \(\Delta , \Delta _{\delta }^{+}, \Delta _{\delta }^{-}\) } is actually a shorthand notation for { \(\Delta , \Delta _{x}^{+}, \Delta _{y}^{+}, \Delta _{x}^{-}, \Delta _{y}^{-}\) }.

In this article, we also include the possibility of the ESP-type superconducting solutions along with OSP-type solutions. Thus our new set of pairing correlations include four more quantities, namely \(\Delta _{x}^{\uparrow }\), \(\Delta _{y}^{\uparrow }\), \(\Delta _{x}^{\downarrow }\), and \(\Delta _{y}^{\downarrow }\). So, our new set of order parameters becomes

$$\begin{aligned} \{ \Delta \} = \overbrace{ \{ \Delta ,\quad \Delta _{\delta }^{+},\quad \Delta _{\delta }^{-} }^{\text {OSP}},\quad \overbrace{ \Delta _{\delta }^{\uparrow },\quad \Delta _{\delta }^{\downarrow } }^{\text {ESP}} \}, \end{aligned}$$

where, again we have used a shorthand notation as explained above. Note that, in Eq. (10), we have grouped the pairing correlations in two sets, namely OSP and ESP. The OSP group contains pairing correlations (\(\Delta\),\(\Delta _{\delta }^{+}\),\(\Delta _{\delta }^{-}\)) that allow the OSP-type triplet states to exist, while the ESP group contains pairing correlations (\(\Delta _{\delta }^{\uparrow }\),\(\Delta _{\delta }^{\downarrow }\)) that give rise to ESP-type triplet states. We follow an unrestricted approach, as adopted in Ref.51, where OSP group of pairing correlations are allowed to take such forms where they can either make the singlet or the OSP-type triplet, or a mixture of them them to be stabilized in different parameter regimes. On the other hand, ESP group of pairing correlations allow only ESP-type triplet states to be stabilized because of their symmetrized-spin wavefunction (and therefore the orbital part of the wavefunction which is the superconducting gap function, should be anti-symmetric) by definition.

Having defined our order parameters in the above manner, we set out to study the variation of magnitudes of these order parameters with varying chemical potential \(\mu\). Average electron density per site \(\langle n\rangle\) is calculated for each value of \(\mu\). This allows one to plot magnitudes of the order parameters, defined in Eq. (10), with varying \(\langle n\rangle\). In the self-consistent approach, if we start with a random initial configuration of \(\{ \Delta \}\) in the absence of magnetic field, we find individual magnitudes of the converged solutions to be arbitrarily different for different initial runs due to degeneracy. This prevents us to plot a smooth variation of the OPs with \(\langle n\rangle\). In the absence of Zeeman coupling this is quite expected. As in the absence of Zeeman coupling, there is no fixed quantization axis present in the system, possible OSP and ESP-type triplet states can form different linear combinations of corresponding pairing correlations, leaving us with a degenerate set of solutions for each value of \(\mu\). In this case the d-vector formalism, defined in the “d-vector formalism” section, comes to our rescue. We find magnitudes of \(d_{0}(\mathbf {k})\) and \(\mathbf {d}(\mathbf {k})\), averaged over the whole Brillouin zone, to be the best suited order parameters for the singlet and the triplet states respectively, in this scenario. Thus we define our singlet and triplet order parameters in the following way:

$$\begin{aligned} \Delta _{\text {SP}} = \frac{1}{N_{s}}\sum \limits _{\mathbf {k}} |d_{0}(\mathbf {k})|\quad \text {and}\quad \Delta _{\text {TP}} = \frac{1}{N_{s}}\sum \limits _{\mathbf {k}} |\mathbf {d}(\mathbf {k})|. \end{aligned}$$
Figure 2
figure 2

The color in the background indicates the nature of the SC order as, red: pure-singlet (SP), green: pure-triplet (TP), and blue: singlet–triplet mixture (MP). The gray background represents non-superconducting region (NS). (a–d) Variation of singlet \(\Delta _{SP}\) and triplet \(\Delta _{TP}\) order parameter magnitudes with electron density \(\langle n\rangle\) at zero temperature for onsite interaction \(U=t\) and inter-site interaction strength: (a) \(V = 0.4t\), (b) \(V = 1.6t\), (c) \(V = 3.0t\) and (d) \(V = 4.0t\). (e) Plot (e) shows the V-\(\langle n\rangle\) ground state phase diagram at \(B=0\), \(U=t\). (f,g) Variation of singlet \(\Delta _{SP}\) and triplet \(\Delta _{TP}\) order parameter magnitudes with temperature T in the absence of Zeeman-coupling for on-site interaction \(U=t\), inter-site interaction \(V=3t\), and average electron density: (f) \(\langle n\rangle =1.0\) and (g) \(\langle n\rangle =0.44\).

Using the newly defined singlet and triplet order parameters in Eq. (11), we study the variation of these order parameters with average electron density per site \(\langle n\rangle\). In Fig. 2a–d, we plot the variation of singlet order parameter \(\Delta _{SP}\) and triplet order parameter \(\Delta _{TP}\) with \(\langle n\rangle\) for different values of inter-site attraction V. We keep on-site attraction at a fixed value (\(U=t\)), as it is the inter-site attraction V that plays a role in stabilizing different unconventional superconducting orders as deduced from our previous works51. When inter-site attraction V is small (\(V=0.4t\)) compared to on-site attraction U, it is the singlet phase that occupies the whole density region (Fig. 2a). At \(V=1.6t\), singlet phase occurs both at lower density region and near half-filling (Fig. 2b). Though orbital information of the order parameter is averaged out in defining the singlet and triplet order parameters, it’s easy to suggest that the singlet phase near half-filling is \(d_{x^2-y^2}\) type, while the singlet phase near low-density region is of \(s+s^{*}\) type based off our previous work51. Pure triplet phase occurs near quarter-filling, while singlet–triplet mixed parity phase occurs between pure triplet phase and singlet phase of \(d_{x^2-y^2}\) type. The pure triplet phase consists of all three types of triplet pairings due to absence of magnetic field and therefore it is meaningless to give it a specific name. It is important to note that pure-\(d_{x^2-y^2}\) type singlet phase occurs at low-V region, when V is still larger than on-site interaction U. At higher value of inter-site attraction (\(V=3t\)), pure-\(d_{x^2-y^2}\) type singlet phase vanishes (Fig. 2c), while singlet–triplet mixed phase occupies the corresponding density region (near half-filling). At \(V=3t\) a small window of singlet–triplet mixed phase occurs between \(s+s^{*}\) type singlet and pure triplet phase. If we increase V further (\(V=4t\)), most of the density region is occupied by dominant singlet–triplet mixed phase, while \(s+s^{*}\) type singlet phase still occurs at low density region (Fig. 2d).

This variation of order parameters with average electron density motivates us to draw a V-\(\langle n\rangle\) ground state phase diagram in the absence of a magnetic field (Fig. 2e). We notice that when the inter-site interaction V is less than the on-site interaction, i.e. \(V < t\), the whole density region is occupied by the pure singlet phase (red in color). This is expected as a dominant on-site attraction is known to stabilize singlet s-wave SC order. While U is stabilizing singlet s-wave SC order, V is also playing an important role in stabilizing the singlet \(s^{*}\) and the unconventional \(d_{x^2-y^2}\) SC orders in this region along with s-wave. When the lower density region prefers singlet s and \(s^{*}\) SC orders to other phases, near half-filling it is the singlet \(d_{x^2-y^2}\)-wave that prevails. Notice that these inferences come from a combination of the present work and the previous work51. As inter-site interaction V becomes larger compared to on-site interaction U, it not only stabilizes singlet \(s^{*}\) and \(d_{x^2-y^2}\) SC orders, it also helps pure triplet and mixed-parity phases to get stabilized at different density regimes. The pure triplet phase (green in color) occurs in a delta-shaped region near quarter filling (\(0.4 \le \langle n\rangle \le 0.6\)), when the inter-site attraction V is higher than the on-site attraction U but not as large as \(V=4t\). With V getting larger, the mixed-parity (blue in color) phase becomes more stable, specially in the higher density region. It’s interesting to note that \(d_{x^2-y^2}\)-type singlet phase is prone to get stabilized near half-filling region and pure triplet phase has a tendency to get stabilized near quarter-filling. But depending on the strength of inter-site interaction the system gains energy by stabilizing a new mixed-parity (possibly \(d+p\)-type extrapolating from our previous study) state, near the higher density region. It’s also important to note that, in the absence of a magnetic field, the phase diagram is very similar to our previous study51. The only difference is that in the triplet region, we expect the phase to be in a superposition of all possible triplet orders, which was not possible earlier because of the limitations of our exploration space. In the section where we study the effects of magnetic field, “Effect of Zeeman coupling” section, the phase diagram will change drastically since now we have the possibility of a competition between different triplet phases. Given that the values of V/t and B/t used in the study are unusually large, we do not claim that the model can be directly applicable to a known superconductor. Nevertheless, our results become relevant to systems with nearly flat bands that lead to a very small effective t. The model may also be realized in optical lattices of ultracold atomic gases. Furthermore, while it appears that the interesting competition between different SC phases begins only for \(V>t\) (see Fig. 2c), it is important to note that the interesting phases in fact show up for \(V > U\)51. The choice \(U=1\) is merely for computational convenience as the lattice sizes required to provide stable results are smaller for larger values of U and V (see “Methods” section). In fact, one of the possible realizations of inter-site attractive interactions occurs via \(t-J\) model where the onsite term is repulsive. This leads to a already well studied problem of competition between magnetism and superconductivity. Given that our aim here is to study the competition between different SC states, we avoided the scenario where magnetism leads to further complexity. We again remind the reader that all the results discussed here are obtained within a mean-field approach.

Figure 3
figure 3

The color in the background indicates the nature of the SC order as, red: pure-singlet (SP), green: pure-triplet (TP), blue: singlet–triplet mixture (MP), and gray: non-superconducting region (NS). White region represents phase separation (PS). (a) T-\(\langle n\rangle\) ground state phase diagram, in the absence of a magnetic field, at \(U=t\) and \(V=3t\). (b) V-\(\langle n\rangle\) ground state ground state phase diagram at finite magnetic field, \(B=0.1t\) and \(U=t\). (c) B-\(\langle n\rangle\) ground state phase diagram at zero temperature, \(V=2t\), \(U=t\). (d) B-T phase diagram for \(U=t\) and \(V=3t\) at half-filling, \(\langle n\rangle = 1\). Note that the phases are expected to be symmetric about half-filling, \(\langle n\rangle = 1\).

Finite temperature study

After exploring the stability of singlet, triplet, and mixed-parity states in the ground state of the extended attractive Hubbard model on a square lattice, we intend to study the behavior of the corresponding SC states at finite temperatures. General consensus is that thermal excitations destroy superconductivity. So the expectation is that the magnitudes of the superconducting order parameters, defined in our system, will decrease and eventually vanish when temperature is increased beyond a critical value. It is well-known that beyond this critical temperature \(T_{c}\) superconductors make a transition from superconducting state to the normal state. As we limit ourselves to the framework of mean-field theory, it is not our aim to infer about the specific values of \(T_{c}\) we obtain in our calculations. We rather intend to study how different SC orders react to the onset of finite temperature.

In Fig. 2f,g, we plot the variation of singlet \(\Delta _{SP}\) and triplet \(\Delta _{TP}\) order parameter with temperature for two different points in parameter space taken from the ground state V-\(\langle n\rangle\) phase diagram. In Fig. 2f, we start with a mixed-parity state exactly at half-filling \(\langle n\rangle =1\). While the red dashed line represents the singlet order parameter \(\Delta _{SP}\), the green dotted line represents triplet order parameter \(\Delta _{TP}\). As temperature is increased, the triplet order parameter starts to fall off and around \(T\approx 0.38\) it vanishes completely. Interestingly, the singlet order parameter \(\Delta _{SP}\) seems to remain constant, and begins to reduce only when \(\Delta _{TP}\) has vanished. On a careful analysis, this turns out to be a simple quantitative effect. In the MP state, the total gap has contributions from both singlet and triplet order parameters. At low but finite T, excitations across the SC gap suppress the smaller OP relatively strongly compared to the larger OP. Consequently, the larger OP appear to remain unaffected as the smaller OP completely vanishes. We have explicitly checked that if the two OPs are comparable in magnitude than both are simultaneously affected, and in that case the transition is directly from a MP SC state to a non-SC state. At \(T\approx 0.38\) the system makes a transition from a mixed-parity SC state to a pure \(d_{x^2-y^2}\) type singlet superconducting state. If temperature is increased further \(\Delta _{SP}\) starts to fall off and it completely vanishes at \(T\approx 0.68\). Thus Fig. 2f depicts how the system makes a transition from a mixed-parity superconducting state to a pure singlet superconducting state before making a transition to a non-superconducting state. In Fig. 2g, we start with a different point of the parameter space, again taken from the ground state V-\(\langle n\rangle\) phase diagram. In this case our zero temperature initial choice is at density \(\langle n\rangle = 0.44\), where we find pure triplet state being stable. Notably with the increase of temperature the triplet order parameter \(\Delta _{TP}\) falls off exactly the way it did for the previous case and the temperature where \(\Delta _{TP}\) completely vanishes is exactly same as the previous case, i.e., \(T\approx 0.38\). But in this case the triplet SC state makes a direct transition to non-superconducting state.

Finally we summarize these results in Fig. 3a, where we present the T-\(\langle n\rangle\) phase diagram to complete the analysis of the finite temperature effects. In the T-\(\langle n\rangle\) phase diagram we find that every superconducting phases at zero temperature makes a transition from superconducting to non-superconducting phase when temperature is increased. While the transition temperature for singlet phase takes a dome-shaped feature as a function of density, for triplet and mixed-parity states the transition temperature does not have such sharp feature with varying density. As discussed earlier, we find four different sectors of such transition. In three sectors the system makes a direct transition from superconducting state to non-superconducting state retaining the specific form of pairing symmetry, while the mixed-parity state in the density region \(0.5\le \langle n\rangle \le 1.0\) makes a transition to singlet superconducting state before making a transition to non-superconducting state.

Effect of Zeeman coupling

After analysing the behaviour of unconventional superconducting phases at zero and finite temperature in the absence of a magnetic field, we study the effects of the magnetic field on these phases in different regions of the parameter space. We limit our study to the effect of Zeeman-coupling to an external magnetic field. Thus, we ignore the effects of an external magnetic field on orbital degrees of freedom and confine ourselves to investigate the effects of the magnetic field on the spin degrees of freedom only. An in-plane magnetic field serves our purpose. Motivated by the zero temperature V-\(\langle n\rangle\) phase diagram presented in Fig. 2e, we take interest in knowing how this phase diagram changes with the onset of a Zeeman field. Figure 3b shows the V-\(\langle n\rangle\) phase diagram in the the presence of magnetic field with \(B = 0.1t\). The first thing we notice is the appearance of a non-superconducting region (gray in color), which is quite expected as the magnetic field destroys superconductivity. It is interesting to note that the non-superconducting region appears at lower values of V suggesting that for the superconductivity to survive it really is a competition between the inter-site interaction strength and the magnetic field strength.

Figure 4
figure 4

Variation of average electron density \(\langle n\rangle\) with chemical potential \(\mu\) for, (a) \(B=0\) and \(B=0.1t\). The discontinuity in (b) corresponds to a region of phase separation.

From the ground state V-\(\langle n\rangle\) phase diagram in the absence of a magnetic field Fig. 2e, we notice that when inter-site interaction V is less than the on-site interaction U, the whole density region is occupied by singlet superconducting phase. When the Zeeman field is switched on and the strength of the Zeeman field is maintained at \(B = 0.1t\), it turns out to be strong enough to destroy superconductivity when V is approximately less than 0.8t. So any superconducting phase occurring below \(V\approx 0.8t\) is destroyed by a Zeeman field of strength \(B = 0.1t\). In the absence of a magnetic field, there is a pure \(d_{x^2-y^2}\) type singlet phase, which is stable near half-filling. It appears that presence of an external magnetic field of strength \(B = 0.1t\), destroys such pure \(d_{x^2-y^2}\) type superconducting phase, as for any strength of V (within the scale provided in Fig. 3b), only the singlet–triplet mixed state gets stabilized. Most interestingly, the region, where pure triplet superconducting phase is stabilized in V-\(\langle n\rangle\) phase diagram, gets larger in presence of the Zeeman-coupling suggesting that the pure triplet superconducting phase gets enhanced by the onset of a magnetic field, this is due to the allowed ESP superconducting correlations.

Another important aspect we notice is that there is a region of phase separation appearing between the pure singlet phase and the pure triplet phase when V is approximately greater than 2t. To look deeper and confirm this aspect of phase separation, we study the variation of average electron density per site \(\langle n\rangle\) with the chemical potential \(\mu\) at two values of Zeeman-coupling, (a) \(B=0\), and (b) \(B=0.1t\) presented in Fig. 4. For the purpose of the above-mentioned illustration we keep the on-site interaction at \(U=t\), while the inter-site interaction is kept fixed at \(V=4t\). While in the absence of a Zeeman-coupling we see a continuous variation of \(\langle n\rangle\) with \(\mu\), in the presence of a Zeeman-coupling, we clearly notice a discontinuity in the electron density at \(\mu \approx - 2.5t\). The discontinuity in the density is of the order \(\delta _{\langle n\rangle }\approx 0.1\), at \(V=4t\) (Fig. 3b).

To completely understand the effect of Zeeman-coupling on the superconducting phases, we present, in Fig. 3c, a zero temperature B-\(\langle n\rangle\) phase diagram at \(U=t\) and \(V=2t\). At \(B=0\), \(U=t\) and \(V=2t\), singlet, triplet, and mixed-parity superconducting phases occupy different sectors of density regions, as observed in the ground state V-\(\langle n\rangle\) phase diagram in the absence of the Zeeman-coupling to an external magnetic field. As the strength of the magnetic field is increased it appears that the pure triplet superconducting state becomes favourable. In fact, the pure singlet superconducting phase disappears from the whole density region when the strength of the magnetic field becomes approximately greater than 0.23t, while the mixed-parity state completely disappears from the whole density region when B becomes larger than \(\sim 0.4t\). Eventually at \(B>0.4t\) only pure triplet superconducting phase survives throughout the whole density region. Though we must take a note of the fact that a non-superconducting region (gray in color) prevails over the other phases at lower electron densities. The information we gather from Fig. 3c is of twofold nature. Firstly, it is clear that Zeeman-coupling to an external magnetic field at low strength favours pure triplet superconductivity. Secondly, the disappearance of pure singlet and mixed-parity state indicates the possibility of disappearance of OSP type superconducting states and appearance of ESP type superconducting states beyond certain strength of the external magnetic field. As OSP type states give rise to singlet states like s,\(s^{*}\) and \(d_{x^2-y^2}\) etc., and ESP type states give rise to pure triplet superconducting phases, the above-mentioned behavior is well-explained.

We conclude this section by looking at the effect of Zeeman-coupling at finite temperatures. Figure 3d shows the B-T phase diagram to describe these effects at a specific point of the parameter space, namely \(U=t\), \(V=3t\) and \(\mu =0\). As expected, we see that for high enough temperatures, the normal state takes over for all kinds of superconducting phases. Also notice that at higher magnetic field strength, the superconductivity is destroyed at lower critical temperatures, in contrast to the critical temperatures at lower magnetic field strengths.

Characterization of phases

Our main aim in this study is to provide an identification, based on explicit energy minimization, of distinct-symmetry superconducting order parameters in the presence of a Zeeman coupling. So far we have classified these superconducting states in terms of opposite spin pairing (OSP) and equal spin pairing (ESP) states, leading to a description of minimum energy solutions in terms of pure singlet, pure triplet, and mixed-parity superconducting states. While, the main classification is done based on the nature of the spin state of the system we do mention different pairing symmetries, such as s-wave, \(s^{*}\)-wave, \(p_{x}+ip_{y}\)-wave and \(d_{x^2-y^2}\)-wave to define the nature of the orbital part of the gap function51,52. In the case of pure singlet, or pure triplet phase the nature of the orbital part can be easily understood, while for mixed-parity states such a separation does not exist. In Fig. 5 we look at various quantities related to the superconducting order parameters at zero temperature with varying average electron density per site \(\langle n\rangle\). Results presented in the first column of panels in Fig. 5 correspond to the case when the Zeeman field \(B = 0.12t\). On the other hand, the second column represents data for a relatively higher value of the Zeeman field, \(B = 0.6t\). The on-site attraction potential U and the inter-site attraction potential V are kept fixed at \(U = t\) and \(V = 2t\), for all the plots in Fig. 5.

Figure 5
figure 5

The color in the background indicates the nature of the SC order as, cyan: time reversal symmetry preserving (TRP) SC order, magenta: time reversal symmetry breaking (TRB) SC order, and gray: non-superconducting (NS) region. (a,b) Variation of \(\Delta _{SP}\) and \(\Delta _{TP}\) as a function of average electron density \(\langle n\rangle\) for (a) \(B = 0.12t\), and (b) \(B = 0.6t\). Inset in (b) displays orbital part of the singlet component of SC order plotted in (a). Note that the presence of s-wave and \(s^{*}\)-wave in pure-singlet phase, and the presence of \(d_{x^{2}-y^{2}}\)-wave character in the mixed-parity phase is confirmed from this plot. (c,d) SC pair correlation averages \(\Delta _{x}^{\uparrow }\), \(\Delta _{x}^{\downarrow }\), \(\Delta _{y}^{\uparrow }\), and \(\Delta _{y}^{\downarrow }\) as a function of \(\langle n\rangle\) for (c) \(B = 0.12t\), and (d) \(B = 0.6t\). (e,f) Variation of spin-resolved average electron densities (\(\langle n_{\uparrow }\rangle\) and \(\langle n_{\downarrow }\rangle\)) and \(\Delta _{\mathbf {q}}\), representing the non-unitarity of the SC order for (e) \(B = 0.12t\), and (f) \(B = 0.6t\). For clarity \(\Delta _{\mathbf {q}}\) is scaled by appropriate factors in (e,f). We have used \(U = t\), \(V = 2t\), and \(T = 2\times 10^{-4}t\) for all the data shown in this figure. Since there are multiple quantities displayed in each panel, a common y-axis label is not possible. The quantities plotted are directly mentioned with the help of legends.

Figure 6
figure 6

Angular variation of the superconducting energy gaps corresponding to the two pseudo-spin degrees of freedom, for different values of Zeeman field strength \(B = 0\), \(B = 0.12t\), and \(B = 0.6t\). The values of other parameters are fixed as \(U = t\), \(V = 2t\), \(\langle n\rangle \approx 0.6\) for all the panels in this figure. (a,b) Superconducting energy gaps, corresponding to two pseudo-spins, at the \(\uparrow\)-spin Fermi surface in the presence of magnetic field. (c) \(\uparrow\)-spin Fermi surface for different values of Zeeman field strength \(B = 0\), \(B = 0.12t\), and \(B = 0.6t\). (d,e) Superconducting energy gaps, corresponding to two pseudo-spins, at the \(\downarrow\)-spin Fermi surface in the presence of magnetic field. (f) \(\downarrow\)-spin Fermi surface for different values of Zeeman field strength \(B = 0\), \(B = 0.12t\), and \(B = 0.6t\).

In Fig. 5a, we plot the singlet and triplet order parameters as defined previously. We observe three different sectors in the whole density profile: pure singlet at lower density region, pure triplet near quarter filling, and mixed-parity solution around half filling. Further characterization of the singlet phases, based on the definitions of s-wave, \(s*\)-wave, and d-wave are shown in the inset of Fig. 5b. Note that near half filling, where the singlet phase shows d-type pairing (inset of Fig. 5b), it is a mixed parity state with \(d_{x^2-y^2}\) as a singlet component. The presence of s-wave and \(s^{*}\)-wave character in pure singlet phase and \(d_{x^{2}-y^{2}}\)-wave character in the mixed-parity phase is consistent with our previous findings when the Zeeman field is absent51. These pairing symmetries do not change with the addition of magnetic field, however, their relative magnitudes do change as discussed in the previous sections. Finite values of singlet order parameter \(\Delta _{SP}\) and triplet order parameter \(\Delta _{TP}\) at different sectors of the doping regime confirm the existence of pure-singlet, pure-triplet, and mixed-parity phases in the system at a lower strength of the magnetic field, \(B = 0.12t\). The background color indicates the time reversal symmetry of the superconducting phase (d-vector formalism also allows us to probe the time reversal symmetry in a superconducting phase, refer “d-vector formalism” section). Notice that even though Fig. 5a is generated for non-zero magnetic field \(B=0.12t\), we see the existence of time reversal invariant (TRP-SC) superconducting phase appearing for lower densities. This is the enduring s-wave superconducting phase which is present in the absence of magnetic field at lower densities (Fig. 3). We can see from Fig. 5b that as we increase the magnetic field \(B=0.6t\), the remnant time reversal invariant phase is destroyed. Of course, this transition is smooth with the increase of Zeeman field. We also notice that the mixed parity phase breaks time reversal, in order to take advantage of the external magnetic field. At a very low electron density, superconductivity completely vanishes, giving rise to the non-superconducting (NS) region. Fig. 5a,b shows that, in the presence of a magnetic field, it is the ESP-type correlation that becomes more favourable, although at low enough Zeeman fields, other correlations can also survive.

In Fig. 5c,d we plot the magnitudes of the ESP-type correlations, \(\{\Delta _{x}^{\uparrow }, \Delta _{x}^{\downarrow }, \Delta _{y}^{\uparrow }, \Delta _{y}^{\downarrow }\}\), for the cases presented in Fig. 5a,b respectively. we observe that in the pure triplet state, near quarter filling, \(\Delta _{x}^{\uparrow }\) and \(\Delta _{y}^{\uparrow }\) become equal in magnitude. Similarly magnitudes of \(\Delta _{x}^{\downarrow }\) and \(\Delta _{y}^{\downarrow }\) also become equal, while \(\Delta _{\delta }^{\uparrow }\) and \(\Delta _{\delta }^{\downarrow }\) acquire different values due to the Zeeman coupling. This observation was also inferred from “Determination of relative phase angle between different pairing correlations” section. Interestingly, even in the presence of Zeeman field, this relation among the triplet ESP correlations appears to change as the \(d_{x^2-y^2}\) singlet component becomes dominant near half filling. The character of \(d_{x^2-y^2}\) breaks the symmetry between the x and y components of the correlation functions, and as \(d_{x^2-y^2}\) wave state becomes more dominant near half filling, this character also shows up in the ESP pairing correlations, forcing \(\Delta _{x}^{\uparrow }\) and \(\Delta _{x}^{\downarrow }\) to become equal in magnitude and similarly for the y component as well. Now as we increase the magnetic field, the pure triplet phase takes over most of the density profile and naturally, only time reversal breaking character of the ESP pairing correlations remains. There is however, an interesting observation that exactly at half filling the \(\uparrow\) and \(\downarrow\) symmetry is forced upon the triplet phase, even in the presence of strong Zeeman field. Since this formalism is intrinsically particle-hole symmetric, at half-filling we see this symmetry showing up as the symmetry between \(\uparrow\) and \(\downarrow\) spins which is why all the ESP correlations functions, even in the presence of Zeeman field are forced to take the same values.

In Fig. 5e spin-resolved average electron densities are plotted along with the magnitude of the \(\mathbf{q}\) vector averaged over the Brillouin zone, \(\Delta _\mathbf{q} \equiv \sum _\mathbf{k} |\mathbf{q}|\). Here \(\mathbf{q}\)-vector is defined in terms of the \(\mathbf{d}\)-vector as \(\mathbf{q}=i~\mathbf{d}\times \mathbf{d}^*\). The \(\mathbf{q}(\mathbf{k})\)-vector13 physically represents the net spin average present in the pairing state with momentum \(\mathbf {k}\). This does not always ensure a net total spin moment, averaged over the Fermi surface. However, it entails the fact that the pair correlation for \(\uparrow\)-spin electrons is different than that of the \(\downarrow\)-spin electrons, thus non-unitarity of a SC order is usually associated with the breaking of time reversal symmetry (TRS)53. Non-unitarity also implies the opening of two distinct SC energy gaps, driven by the breaking of TRS. In the pure triplet region we observe the possibility of the formation of an effective magnetic moment, however, interestingly enough, even in the presence of (small) magnetic field, in the regions where singlet superconducting phases are present, the spin-resolved electron densities are always equal and thus the magnetic moment is always zero. The magnitude of the \(\mathbf{q}\) vector, \(\Delta _\mathbf{q}\) tells us that the gap function is non-unitary (unitary when \(\Delta _\mathbf{q}=0\)) in nature when the triplet order parameter is finite. Note that non-unitarity of triplet gap function is defined in “d-vector formalism” section. From Fig. 5e,f, we see that non-unitarity is enhanced due to Zeeman field which is in correspondence to the enhancement of the pure triplet phase. It is again interesting to note that quarter filling seems to be the most favourable region for a non-unitary phase to exist and as we go close to half-filling, the non-unitarity is destroyed, consequently, the triplet superconducting phase becomes unitary exactly at half-filling due to the symmetry of the Hamiltonian.

We conclude the section by presenting, in Fig. 6, the angular variation of the superconducting energy gaps corresponding to the two pseudo-spin degrees of freedom, in the presence of the Zeeman field. We choose a representative point at \(U = t\), \(V = 2t\), and \(\langle n\rangle \approx 0.6\). In the presence of the Zeeman field, the non-interacting Fermi surface splits into two, resulting in separate Fermi surfaces for \(\uparrow\)-spin (see Fig. 6c) and \(\downarrow\)-spin (see Fig. 6f) electrons. Interactions lead to a gap opening at these Fermi surfaces, resulting in two gaps as the solution of a \(4 \times 4\) matrix at a given \(\mathbf{k}_F\)-point leads to two pairs of \(\pm E(\mathbf{k_F})\) eigenvalues. Although it may not be straightforward to measure these two gaps in experiments, we nevertheless provide the description in terms of two distinct gaps. In Fig. 6a,b we plot the magnitude of these two distinct superconducting energy gaps at \(\uparrow\)-spin Fermi surface. Similarly, in Fig. 6d,e we plot the magnitude of gaps at \(\downarrow\)-spin Fermi surface. This angular variation of the superconducting energy gaps encodes the symmetry of the underlying superconducting state: while at \(B = 0\), a \(d_{x^2-y^2}\)-wave character is prominent from the symmetry of the gap, at \(B = 0.12t\) the chiral \(p_{x}\pm i p_{y}\) seems to be the dominant symmetry present in the gap structure. Interestingly, at a higher value of the magnetic field, \(B = 0.6t\), one of the SC gaps approaches an isotropic form (see Fig. 6a,d) while the other one retains the \(p_{x}\pm i p_{y}\) character see Fig. 6b,e. Note that each of the two gaps at the up and down spin Fermi surfaces contains a linear combination of various OPs. Therefore, one can conclude that for larger values of the magnetic field, one of these linear combinations is such that all angular dependence gets canceled out. A slight difference in anisotropy between gap structure in panels (b,e) can also be noticed for larger B. This is related to the B-induced difference in the shapes of the up and the down spin Fermi surfaces.


Using the unrestricted BdG mean-field approach to the EAHM, we study the competition among different types of superconducting orders. The unique aspect of our study is that we rely on energetics for obtaining the superconducting solutions with different order parameter symmetries. We find rich ground state phase diagrams with exotic phases such as unitary, non-unitary, pure triplet and mixed parity states. Based on energetics, we provide an understanding of specific phase relations between different superconducting order parameters. We find an interesting competition between the anti-aligned, w.r.t. the Zeeman field, \(\downarrow \downarrow\) ESP state and the OSP state. At moderate values of the Zeeman field, the system prefers the OSP state however, at larger values of Zeeman-coupling, OSP state is completely disfavoured in the weak-coupling limit. The \(\uparrow \uparrow\) ESP state is favoured in the presence of Zeeman field. We use the d-vector formalism to further distinguish and characterize the different triplet phases. Pure triplet phase is stable near quarter filling and singlet–triplet mixed phase is favoured near half-filling. The definition of the triplet order parameters being rotationally invariant allows us to study their behaviour even in the presence of magnetic field consistently. Finite temperature study shows transitions from a mixed parity ground state to a pure singlet phase of \(d_{x^2-y^2}\)-type, and finally to a non-superconducting phase. In the presence of Zeeman field, the pure triplet phase dominates the phase diagram, suppressing the pure singlet and mixed parity phases. It is interesting to note that the mixed parity phase is relatively more robust against the Zeeman field as compared to the pure singlet phase. Finally, we discuss the key distinction between different phases in terms of the angular shapes of the gap functions that are experimentally measurable. Our results will be particularly useful in understanding the effect of planar magnetic fields on the nature of superconducting states stabilized in quasi two-dimensional systems, such as atomically thin conducting interfaces between insulating oxides.


Mean-field decoupling in pairing channel

Both the interaction terms, the on-site interaction and the inter-site interaction term, are represented by two-body operator terms in the second quantization formalism. To reduce the complexity of these terms, we consider the mean-field treatment, where the many-body interaction effects are mimicked via the interaction of a single electron system with a mean-field created by the aggregated effect of rest of the electrons. Mathematically speaking, such a two-body operator term \({\hat{X}}{\hat{Y}}\) can be approximated as

$$\begin{aligned} {\hat{X}}{\hat{Y}}&= \big \{\langle {\hat{X}}\rangle + ({\hat{X}} - \langle {\hat{X}}\rangle ) \big \}\big \{\langle {\hat{Y}}\rangle + ({\hat{Y}} - \langle {\hat{Y}}\rangle )\big \} \nonumber \\&= \langle {\hat{X}}\rangle \langle {\hat{Y}}\rangle + \langle {\hat{X}}\rangle ({\hat{Y}}-\langle {\hat{Y}}\rangle ) + \langle {\hat{Y}}\rangle ({\hat{X}}-\langle {\hat{X}}\rangle ) + ({\hat{X}}-\langle {\hat{X}}\rangle )({\hat{Y}}-\langle {\hat{Y}}\rangle ) \nonumber \\&\approx {\hat{X}}\langle {\hat{Y}}\rangle + {\hat{Y}}\langle {\hat{X}}\rangle - \langle {\hat{X}}\rangle \langle {\hat{Y}}\rangle , \end{aligned}$$

where we have neglected the second order correction contribution. In practice, we can apply this idea of mean-field decoupling to the interaction terms of the Hamiltonian using different channels, such as density channel or pairing channel. Where density channel decoupling leads to the mean-fields of the form \(\langle c^{\dagger }c\rangle\), and pairing channel decoupling leads to mean-fields of the form \(\langle c^{\dagger }c^{\dagger }\rangle\). As we are specifically interested in superconducting solutions, we carefully choose pairing channel decoupling for our system and allow the superconducting correlation functions to take non-zero values to minimize energy.

Applying the mean field approximation described above on the interaction Hamiltonian (Eq. 7), and treating \({\hat{X}} = c_{i\uparrow }^{\dagger }c_{i\downarrow }^{\dagger }\) and \({\hat{Y}} = c_{i\downarrow }c_{i\uparrow }\), we get:

$$\begin{aligned} \mathcal {H}_{\text {MF}}^{\text {onsite}}&= -U\sum \limits _{i} \left[ c_{i\uparrow }^{\dagger }c_{i\downarrow }^{\dagger }\langle c_{i\downarrow }c_{i\uparrow }\rangle + c_{i\downarrow }c_{i\uparrow }\langle c_{i\uparrow }^{\dagger }c_{i\downarrow }^{\dagger }\rangle \right] \nonumber \\&\quad +U\sum \limits _{i} \left[ \langle c_{i\uparrow }^{\dagger }c_{i\downarrow }^{\dagger }\rangle \langle c_{i\downarrow }c_{i\uparrow }\rangle \right] = -U\sum \limits _{i} \left[ \Delta _{i}c_{i\uparrow }^{\dagger }c_{i\downarrow }^{\dagger } + \text {H.c.} - |\Delta _{i}|^2 \right] , \end{aligned}$$

where, \(\Delta _{i} = \langle c_{i\downarrow }c_{i\uparrow }\rangle\) is defined as the on-site pairing correlation at \(i{\text {th}}\) site. Similarly, applying the mean field approximation to the terms A, B, C and D of Eq. (8), and decoupling into all possible pairing channels, we get:

$$\begin{aligned} \mathcal {H}_{\text {MF}}^{\text {nn}}&= -V\sum \limits _{i,\delta } \left[ \Delta _{i\delta }^{\uparrow } c_{i\uparrow }^{\dagger }c_{i+\delta \uparrow }^{\dagger } + \text {H.c.} - |\Delta _{i\delta }^{\uparrow }|^2 \right] -V\sum \limits _{i,\delta } \left[ \Delta _{i\delta }^{\downarrow } c_{i\downarrow }^{\dagger }c_{i+\delta \downarrow }^{\dagger } + \text {H.c.} - |\Delta _{i\delta }^{\downarrow }|^2 \right] \nonumber \\&\quad -V\sum \limits _{i,\delta } \left[ \Delta _{i\delta }^{+} c_{i\uparrow }^{\dagger }c_{i+\delta \downarrow }^{\dagger } + \text {H.c.} - |\Delta _{i\delta }^{+}|^2 \right] \nonumber \\&\quad -V\sum \limits _{i,\delta } \left[ \Delta _{i\delta }^{-} c_{i+\delta \uparrow }^{\dagger }c_{i\downarrow }^{\dagger } + \text {H.c.} - |\Delta _{i\delta }^{-}|^2 \right] , \end{aligned}$$

where, we have defined the inter-site pairing correlations as:

$$\begin{aligned} \Delta _{i\delta }^{\uparrow } = \langle c_{i+\delta \uparrow }c_{i\uparrow }\rangle ;&\quad \Delta _{i\delta }^{\downarrow } = \langle c_{i+\delta \downarrow }c_{i\downarrow }\rangle \nonumber \\ \Delta _{i\delta }^{+} = \langle c_{i+\delta \downarrow }c_{i\uparrow }\rangle ;&\quad \Delta _{i\delta }^{-} = \langle c_{i\downarrow }c_{i+\delta \uparrow }\rangle . \end{aligned}$$

After combining all the terms together, our effective mean-field Hamiltonian in real space looks like:

$$\begin{aligned} \mathcal {H}_{\text {MF}}&= -t\sum \limits _{\langle ij\rangle ,\sigma }\left[ c_{i\sigma }^{\dagger }c_{j\sigma } + \text {H.c.}\right] -\mu \sum \limits _{i\sigma }c_{i\sigma }^{\dagger }c_{i\sigma } -B\sum \limits _{i}(c_{i\uparrow }^{\dagger }c_{i\uparrow }-c_{i\downarrow }^{\dagger }c_{i\downarrow }) -U\sum \limits _{i} \left[ \Delta _{i}c_{i\uparrow }^{\dagger }c_{i\downarrow }^{\dagger } + \text {H.c.} - |\Delta _{i}|^2 \right] \nonumber \\&\quad -V\sum \limits _{i,\delta } \left[ \Delta _{i\delta }^{\uparrow } c_{i\uparrow }^{\dagger }c_{i+\delta \uparrow } ^{\dagger } + \text {H.c.} - |\Delta _{i\delta }^{\uparrow }|^2 \right] -V\sum \limits _{i,\delta } \left[ \Delta _{i\delta }^{\downarrow } c_{i\downarrow }^{\dagger }c_{i+\delta \downarrow }^{\dagger } + \text {H.c.} - |\Delta _{i\delta }^{\downarrow }|^2 \right] \nonumber \\&\quad -V\sum \limits _{i,\delta } \left[ \Delta _{i\delta }^{+} c_{i\uparrow }^{\dagger }c_{i+\delta \downarrow }^{\dagger } + \text {H.c.} - |\Delta _{i\delta }^{+}|^2 \right] -V\sum \limits _{i,\delta } \left[ \Delta _{i\delta }^{-} c_{i+\delta \uparrow }^{\dagger }c_{i\downarrow }^{\dagger } + \text {H.c.} - |\Delta _{i\delta }^{-}|^2 \right] . \end{aligned}$$

This Hamiltonian allows for spatially varying superconducting solutions since all the correlations are site dependent. However, working with the periodic boundaries and a uniform structure of superconducting correlation functions so that the translational invariance in preserved has a nice advantage. It allows us to block diagonalize the Hamiltonian by working in the Fourier space which drastically simplifies the problem. This is what we do in the next section.

Effective Hamiltonian in momentum space

Observing the effective mean-field Hamiltonian in real space Eq. (16) it is apparent that the effective Hamiltonian is a function of

  1. 1.

    a set of external parameters: t, \(\mu\), U, V, B, and

  2. 2.

    a set of complex superconducting configurations (pairing correlations): \(\Delta _{i}\), \(\Delta _{i\delta }^{+}\), \(\Delta _{i\delta }^{-}\), \(\Delta _{i\delta }^{\uparrow }\), and \(\Delta _{i\delta }^{\downarrow }\).

In this article, we treat a clean system on a square lattice. Thus we exclude the possibility of inhomogeneity in our solutions, i.e., we expect our superconducting solutions to respect translational symmetry of the underlying lattice. This particular choice of translational symmetry makes Bloch basis to be more suitable to describe the effective Hamiltonian. Thus we apply the following transformation from Wannier basis to Bloch basis:

$$\begin{aligned} c_{i\sigma }&= \frac{1}{\sqrt{N_{s}}}\sum \limits _{\mathbf {k}} e^{-i\mathbf {k}\cdot \mathbf {r}_{i}} c_{\mathbf {k}\sigma }\text { and }\nonumber \\ c_{i\sigma }^{\dagger }&= \frac{1}{\sqrt{N_{s}}}\sum \limits _{\mathbf {k}} e^{i\mathbf {k}\cdot \mathbf {r}_{i}} c_{\mathbf {k}\sigma }^{\dagger }, \end{aligned}$$

where \(N_s\) is the total number of lattice points. This essentially block diagonalize the Hamiltonian and simplifies the problem. In Bloch basis, the effective mean-field Hamiltonian looks like:

$$\begin{aligned} \mathcal {H}_{\text {MF}}&= \sum \limits _{\mathbf {k}}\left[ \epsilon ^{\uparrow }(\mathbf {k})c_{\mathbf {k}\uparrow }^{\dagger }c_{\mathbf {k}\uparrow } + \epsilon ^{\downarrow }(\mathbf {k})c_{\mathbf {k}\downarrow }^{\dagger }c_{\mathbf {k}\downarrow }\right] +\sum \limits _{\mathbf {k}}\left[ \Delta ^{\uparrow \downarrow }(\mathbf {k}) c_{\mathbf {k}\uparrow }^{\dagger }c_{-\mathbf {k}\downarrow }^{\dagger } + \text {H.c.} \right] +\sum \limits _{\mathbf {k}}\left[ \Delta ^{\uparrow \uparrow }(\mathbf {k}) c_{\mathbf {k}\uparrow }^{\dagger }c_{-\mathbf {k}\uparrow }^{\dagger } + \text {H.c.} \right] \nonumber \\&\quad +\sum \limits _{\mathbf {k}}\left[ \Delta ^{\downarrow \downarrow }(\mathbf {k}) c_{\mathbf {k}\downarrow }^{\dagger }c_{-\mathbf {k}\downarrow }^{\dagger } + \text {H.c.} \right] + N_{s}U|\Delta |^2 + N_{s}V\sum \limits _{\delta }\left[ |\Delta _{\delta }^{+}|^2 + |\Delta _{\delta }^{-}|^2 + |\Delta _{\delta }^{\uparrow }|^2 + |\Delta _{\delta }^{\downarrow }|^2 \right] , \end{aligned}$$


$$\begin{aligned} \epsilon ^{\sigma }(\mathbf {k})&= -2t[\cos (k_x)+\cos (k_y)] -\mu - (-1)^\sigma B \nonumber \\ \Delta ^{\uparrow \downarrow }(\mathbf {k})&= -U\Delta -V\sum \limits _{\delta } \left[ \Delta _{\delta }^{+}e^{-i(\mathbf {k}\cdot \varvec{\delta })} + \Delta _{\delta }^{-}e^{i(\mathbf {k}\cdot \varvec{\delta })} \right] \nonumber \\ \Delta ^{\sigma \sigma }(\mathbf {k})&= -V\sum \limits _{\delta }\Delta _{\delta }^{\sigma }e^{-i(\mathbf {k}\cdot \varvec{\delta })}, \end{aligned}$$

here \(\sigma ={0,1}\) represent the spins \(\{\uparrow ,\downarrow \}\). To make use of translational invariance, we assumed spatial homogeneity of superconducting pairing correlations. So the set of variables: {\(\Delta _{i},\Delta _{i\delta }^{\uparrow },\Delta _{i\delta }^{\downarrow },\Delta _{i\delta }^{+}\) and \(\Delta _{i\delta }^{-}\)} becomes independent of site index i and our new set of variables becomes: {\(\Delta ,\Delta _{\delta }^{\uparrow },\Delta _{\delta }^{\downarrow },\Delta _{\delta }^{+}\) and \(\Delta _{\delta }^{-}\)}. In this effective Hamiltonian, \(\epsilon ^{\uparrow }(\mathbf {k})\) describes the kinetic energy of \(\uparrow\)-spin electrons with respect to an effective chemical potential \(\mu _{\text {eff}}^{\uparrow } = \mu +B\). Similarly \(\epsilon ^{\downarrow }(\mathbf {k})\) describes the kinetic energy of \(\downarrow\)-spin electrons with respect to an effective chemical potential \(\mu _{\text {eff}}^{\downarrow } = \mu -B\). It is evident that the presence of a magnetic field breaks the spin-degeneracy in the system via Zeeman coupling. It’s worth mentioning that we work in the regime where magnetic field couples only with spin-degree of freedom. Orbital degree of freedom remains completely unperturbed by the presence of magnetic field, which leaves us with the freedom to ignore the Peierls substitution. In the above mentioned Hamiltonian Eq. (18) \(\Delta ^{\uparrow \uparrow }(\mathbf {k})\) and \(\Delta ^{\downarrow \downarrow }(\mathbf {k})\) represent superconducting correlations corresponding to equal spin pairing (ESP) states, while \(\Delta ^{\uparrow \downarrow }(\mathbf {k})\) represents superconducting correlations corresponding to opposite spin pairing (OSP) states. While the on-site attraction strength U controls OSP states only, which is evident from the definitions of these pairing correlations Eq. (19), inter-site attraction V controls both OSP and ESP states. This very fact suggests that V is meant to play an important role in stabilizing non-trivial superconducting orders with variety of pairing symmetries.

Bogoliubov-de Gennes method

The mean field Hamiltonian \(\mathcal {H}_{\text {MF}}\) Eq. (18) can be written in terms of Nambu spinors \({\Psi }_{\mathbf {k}}\) as:

$$\begin{aligned} \mathcal {H}_{\text {MF}} = \sum \limits _{\mathbf {k}}{\Psi }_{\mathbf {k}}^{\dagger }\mathcal {H}_{\text {BdG}}(\mathbf {k}){\Psi }_{\mathbf {k}}, \end{aligned}$$

where, Nambu spinor \({\Psi }_{\mathbf {k}}\) is a column vector of the form, \({\Psi }_{\mathbf {k}} = (c^{}_{\mathbf {k}\uparrow }\quad c_{-\mathbf {k}\downarrow }^{\dagger }\quad c^{}_{\mathbf {k}\downarrow }\quad c_{-\mathbf {k}\uparrow }^{\dagger } )\) and \(\mathcal {H}_{\text {BdG}}(\mathbf {k})\) is a \(4\times 4\) Hamiltonian matrix of the form:

$$\begin{aligned} \mathcal {H}_{\text {BdG}}(\mathbf {k})= \begin{pmatrix} \epsilon ^{\uparrow }(\mathbf {k}) &{} \Delta ^{\uparrow \downarrow }(\mathbf {k}) &{} 0 &{} \Delta _{s}^{\uparrow \uparrow }(\mathbf {k}) \\ \quad \\ (\Delta ^{\uparrow \downarrow }(\mathbf {k}))^{*} &{} -\epsilon ^{\downarrow }(-\mathbf {k}) &{} (\Delta _{s}^{\downarrow \downarrow }(\mathbf {k}))^{*} &{} 0 \\ \quad \\ 0 &{} \Delta _{s}^{\downarrow \downarrow }(\mathbf {k}) &{} \epsilon ^{\downarrow }(\mathbf {k}) &{} -(\Delta ^{\uparrow \downarrow }(-\mathbf {k})) \\ \quad \\ (\Delta _{s}^{\uparrow \uparrow }(\mathbf {k}))^{*} &{} 0 &{} -(\Delta ^{\uparrow \downarrow }(-\mathbf {k}))^{*} &{} -\epsilon ^{\uparrow }(-\mathbf {k}) \end{pmatrix}. \end{aligned}$$

\(\mathcal {H}_{\text {BdG}}(\mathbf {k})\) is known as the mean field Bogoliubov-de Gennes Hamiltonian. We can diagonalize the BdG Hamiltonian by defining new Fermionic quasi-particle operators that mix the electronic operators as,

$$\begin{aligned} c_{\mathbf {k}\uparrow }&= \sum ^{'}_{\mathbf {k},\alpha } u_{\mathbf {k}\uparrow }^{\alpha } \gamma _{\mathbf {k}\alpha } + (v_{\mathbf {k}\uparrow }^{\alpha })^{*} \gamma _{-\mathbf {k}\alpha }^{\dagger } \nonumber \\ c_{\mathbf {k}\downarrow }&= \sum ^{'}_{\mathbf {k},\alpha } u_{\mathbf {k}\downarrow }^{\alpha } \gamma _{\mathbf {k}\alpha } + (v_{\mathbf {k}\downarrow }^{\alpha })^{*} \gamma _{-\mathbf {k}\alpha }^{\dagger }. \end{aligned}$$

The prime on the summation is a restriction to include only those states in the summation that will lead to a positive energy excitation spectrum of the diagonal Hamiltonian. Using this transformation, the Hamiltonian is diagonalized in the following form,

$$\begin{aligned} \mathcal {H}_{\text {MF}} = E_{g} + \sum \limits _{\mathbf {k}\alpha }E_{\mathbf {k}\alpha } \gamma _{\mathbf {k}\alpha }^{\dagger }\gamma _{\mathbf {k}\alpha }, \end{aligned}$$

here \(\{E_{\mathbf {k}\alpha }\}\) is the excitation spectrum for the Bogoliubov quasiparticles. The physical constraint of non-negativity on the excitation energies is implemented by discarding the negative energy states from the definition of the Bogoliubov transformation. This results in \(\{E_{\mathbf {k}\alpha }\}\) to be a set of positive eigenvalues, and \(\gamma _{\mathbf {k}\alpha }^{\dagger }(\gamma _{\mathbf {k}\alpha })\) creates (annihilates) a Bogoliubov quasiparticle with momentum \(\mathbf {k}\) and pseudo-spin \(\alpha\). The BdG quasiparticles describe elementary excitations of the condensate, \(E_{g}\) being the ground state energy of the condensate.

d-vector formalism

In the absence of magnetic field, all the triplet superconducting solutions should be degenerate. However, while performing self-consistency, a randomized initial configuration should converge to a configuration with minimum energy. Due to degeneracy in the absence of magnetic field, whichever initial triplet configuration we provide, the self-consistency approach tends to throw out an arbitrary configuration. To bypass this issue, we divide all the superconducting configurations into different classes based on the characteristics we are most interested in, namely, singlet and triplet pairing. We do so by employing the d-vector formalism.

The most general superconducting gap function can be represented in a matrix form as,

$$\begin{aligned} \Delta (\mathbf{k}) = \begin{pmatrix} \Delta ^{\uparrow \uparrow }(\mathbf{k})&{} \Delta ^{\uparrow \downarrow }(\mathbf{k})\\ \Delta ^{\downarrow \uparrow }(\mathbf{k}) &{} \Delta ^{\downarrow \downarrow }(\mathbf{k}) \end{pmatrix}. \end{aligned}$$

In this, the elements correspond to the spin state of the electrons that constitute the Cooper pair. Furthermore, we may write the Cooper pair wavefunction as,

$$\begin{aligned} \psi _{\sigma \sigma '}(\mathbf{k}) = [ \Delta _0(\mathbf{k}) + \mathbf{d}(\mathbf{k}) \cdot {\varvec{\sigma }} ]i(\sigma _y)_{\sigma \sigma '}, \end{aligned}$$

where, \(\psi _{\sigma \sigma '}(k)\) is the Cooper pair wavefunction with \(\sigma \sigma '\) pairing. The first term, \(\Delta _0(\mathbf{k})\), is the superconducting gap function for the singlet type of superconductor, whereas, the three-components of the triplet superconducting gap function is related to the complex vector \(\mathbf{d}(\mathbf{k})\). The advantage of writing it in this form is that a rotation of spin quantization axis in spin space would be equivalent to a 3D rotation of \(\mathbf{d}(\mathbf{k})\) vector. A rotation of \(\mathbf{d}(\mathbf{k})\) vector, would in turn adjust the superconducting gap functions corresponding to the three spin triplet components accordingly, and therefore, it makes it easier to track them. Now, since the length of \(\mathbf{d}(\mathbf{k})\) vector (and averaged over \(\mathbf{k}\) space) will remain invariant of spin-quantization axis, we may use this quantity to track down the overall triplet component of superconductivity. We therefore use,

$$\begin{aligned} \Delta _{\text {SP}}&= \frac{1}{N_s} \Sigma _\mathbf{k} |\Delta _0(\mathbf{k})|, \end{aligned}$$
$$\begin{aligned} \Delta _{\text {TP}}&= \frac{1}{N_s} \Sigma _\mathbf{k} |\mathbf{d}(\mathbf{k})|, \end{aligned}$$

as our definitions for the singlet superconducting gap function and triplet superconducting gap functions, respectively.

This formalism also gives us an handle on several properties of the superconducting state for example unitarity of the superconducting wavefunction and time reversal symmetry. A superconducting gap function, in matrix form, \(\Delta (\mathbf{k})\) is called unitary if the product \(\Delta (\mathbf{k})\Delta ^{\dagger }(\mathbf{k})\) is proportional to the unit matrix, otherwise the superconducting gap function, in matrix form, is known as non-unitary. With this definition, it is clear that only triplet pairing matrices can be non-unitary since,

$$\begin{aligned} \Delta ^{\dagger }(\mathbf{k})\Delta (\mathbf{k}) = |\mathbf{d(\mathbf{k})}|^2 {\mathbb {I}} + \mathbf {q} \cdot \mathbf {\sigma }, \end{aligned}$$

where \(\mathbf{{q}} = i (\mathbf {d}\times \mathbf {d}^*)\). The measure of \(\Delta _\mathbf{q}\equiv \sum _\mathbf{k}|\mathbf{q}|\) tells us if the superconducting phase is unitary. Also, the time reversal symmetry can easily be checked by,

$$\begin{aligned} \mathbf{d}(\mathbf{k}) {\mathop {=}\limits ^{?}} -\mathbf{d}(\mathbf{k})^* ~~~~~~~~~~~~~~\text {(Time reversal invariant if equal)}. \end{aligned}$$

Therefore, using the d-vector approach we can infer many properties of the superconducting wavefunction.

Self-consistency and minimization

In this section we lay out how the method of self-consistency is implemented computationally. We set electronic hopping parameter \(t = 1\) as the basic energy scale, then we are left with four independent external parameters in the Hamiltonian, viz., \(\mu , U, V, B\). Corresponding to these, we will obtain a set of self-consistent SC pairing correlations, \(\{\Delta \}\) that defines the SC OP. To solve the BdG equations numerically, an initial guess of \(\{\Delta \}\) is fed into the Hamiltonian at some external parameter value. This Hamiltonian is diagonalized and eigenvalues and eigenvectors are calculated. The obtained eigenspectrum is used to further redefine the Hamiltonian and is re-diagonalized. This set of steps is labelled as an iteration and the cycle of iteration is repeated until the averages converge within a specified error, which was set to \(10^{-5}\) in our calculations. The ground state (at \(T=0\)) energy is calculated by summing over all the positive energy states. Therefore, the problem now reduces to minimizing the total energy w.r.t. the set \(\{\Delta \}\) of pairing correlations. Since we do not fix any constraint on the nature of the pairing correlation, energy minimization shall result in the most energetically favourable pairing correlation symmetry that may either by singlet like or triplet like or even some strange combination of these two. This method of calculating averages self-consistently is actually equivalent to the energy minimization. We use \(64 \times 64\) k-grid for \(V \ge 2t\), and a \(512 \times 512\) k-grid for the cases \(V < 2t\).

Code availability

All the codes are written from scratch, in Python or Fortran. Plots are generated using Python scripts.


  1. 1.

    Cooper, L. N. Bound electron pairs in a degenerate fermi gas. Phys. Rev. 104, 1189 (1956).

    ADS  CAS  MATH  Google Scholar 

  2. 2.

    Bardeen, J., Cooper, L. N. & Schrieffer, J. R. Theory of superconductivity. Phys. Rev. 108, 1175 (1957).

    ADS  MathSciNet  CAS  MATH  Google Scholar 

  3. 3.

    Bardeen, J., Cooper, L. N. & Schrieffer, J. R. Microscopic theory of superconductivity. Phys. Rev. 106, 162 (1957).

    ADS  MathSciNet  CAS  MATH  Google Scholar 

  4. 4.

    Migdal, A. Interaction between electrons and lattice vibrations in a normal metal. Sov. Phys. JETP 7, 996–1001 (1958).

    MathSciNet  Google Scholar 

  5. 5.

    Eliashberg, G. Interactions between electrons and lattice vibrations in a superconductor. Sov. Phys. JETP 11, 696–702 (1960).

    MathSciNet  MATH  Google Scholar 

  6. 6.

    Eliashberg, G. Temperature green function for electrons in a superconductor. Sov. Phys. JETP 12, 1000–1002 (1961).

    Google Scholar 

  7. 7.

    Tinkham, M. Introduction to Superconductivity (Courier Corporation, 2004).

    Google Scholar 

  8. 8.

    Bennemann, K. & Ketterson, J. B. History of superconductivity: Conventional, high-transition temperature and novel superconductors. In Superconductivity (eds Bennemann, K. & Ketterson, J. B.) 3–26 (Springer, 2008).

    Google Scholar 

  9. 9.

    Coffey, D. & Coffey, L. Quasiparticle decay effects in the superconducting density of states: Evidence for d-wave pairing in the cuprates. Phys. Rev. Lett. 70, 1529 (1993).

    ADS  CAS  PubMed  Google Scholar 

  10. 10.

    Tsuei, C. & Kirtley, J. Phase-sensitive evidence for d-wave pairing symmetry in electron-doped cuprate superconductors. Phys. Rev. Lett. 85, 182 (2000).

    ADS  CAS  PubMed  Google Scholar 

  11. 11.

    Oates, D., Park, S.-H. & Koren, G. Observation of the nonlinear Meissner effect in ybco thin films: Evidence for a d-wave order parameter in the bulk of the cuprate superconductors. Phys. Rev. Lett. 93, 197001 (2004).

    ADS  CAS  PubMed  Google Scholar 

  12. 12.

    Maki, K., Kee, H. Y. & Morita, Y. Triplet superconductivity in a nutshell. J. Supercond. Novel Magn. 22, 71–74 (2009).

    CAS  Google Scholar 

  13. 13.

    Mackenzie, A. P. & Maeno, Y. The superconductivity of sr 2 ruo 4 and the physics of spin-triplet pairing. Rev. Mod. Phys. 75, 657 (2003).

    ADS  CAS  Google Scholar 

  14. 14.

    Kopnin, N. & Salomaa, M. Mutual friction in superfluid he 3: Effects of bound states in the vortex core. Phys. Rev. B 44, 9667 (1991).

    ADS  CAS  Google Scholar 

  15. 15.

    Volovik, G. Fermion zero modes on vortices in chiral superconductors. J. Exp. Theor. Phys. Lett. 70, 609–614 (1999).

    CAS  Google Scholar 

  16. 16.

    Sarma, S. D., Nayak, C. & Tewari, S. Proposal to stabilize and detect half-quantum vortices in strontium ruthenate thin films: Non-abelian braiding statistics of vortices in a p x+ i p y superconductor. Phys. Rev. B 73, 220502 (2006).

    ADS  Google Scholar 

  17. 17.

    Beenakker, C. Search for majorana fermions in superconductors. Annu. Rev. Condens. Matter Phys. 4, 113–136 (2013).

    ADS  CAS  Google Scholar 

  18. 18.

    Thouless, D. J. Perturbation theory in statistical mechanics and the theory of superconductivity. Ann. Phys. 10, 553–588 (1960).

    ADS  MathSciNet  MATH  Google Scholar 

  19. 19.

    Brueckner, K., Soda, T., Anderson, P. W. & Morel, P. Level structure of nuclear matter and liquid he 3. Phys. Rev. 118, 1442 (1960).

    ADS  MathSciNet  CAS  Google Scholar 

  20. 20.

    Anderson, P. W. & Morel, P. Generalized Bardeen–Cooper–Schrieffer states and the proposed low-temperature phase of liquid he3. Phys. Rev. 123, 1911 (1961).

    ADS  MathSciNet  CAS  Google Scholar 

  21. 21.

    Leggett, A. J. A theoretical description of the new phases of liquid he3. Rev. Mod. Phys. 47, 331 (1975).

    ADS  CAS  Google Scholar 

  22. 22.

    Tou, H. et al. Odd-parity superconductivity with parallel spin pairing in upt 3: Evidence from 195 pt knight shift study. Phys. Rev. Lett. 77, 1374 (1996).

    ADS  CAS  PubMed  Google Scholar 

  23. 23.

    Tou, H. et al. Nonunitary spin-triplet superconductivity in upt 3: Evidence from p 195 t knight shift study. Phys. Rev. Lett. 80, 3129 (1998).

    ADS  CAS  Google Scholar 

  24. 24.

    Aoki, D., Ishida, K. & Flouquet, J. Review of u-based ferromagnetic superconductors: Comparison between uge2, urhge, and ucoge. J. Phys. Soc. Jpn. 88, 022001 (2019).

    ADS  Google Scholar 

  25. 25.

    Lee, I. et al. Triplet superconductivity in an organic superconductor probed by nmr knight shift. Phys. Rev. Lett. 88, 017004 (2001).

    ADS  PubMed  Google Scholar 

  26. 26.

    Ran, S. et al. Nearly ferromagnetic spin-triplet superconductivity. Science 365, 684–687 (2019).

    ADS  CAS  PubMed  Google Scholar 

  27. 27.

    Metz, T. et al. Point-node gap structure of the spin-triplet superconductor ute 2. Phys. Rev. B 100, 220504 (2019).

    ADS  CAS  Google Scholar 

  28. 28.

    Sundar, S. et al. Coexistence of ferromagnetic fluctuations and superconductivity in the actinide superconductor ute 2. Phys. Rev. B 100, 140502 (2019).

    ADS  CAS  Google Scholar 

  29. 29.

    Chronister, A. et al. Evidence for even parity unconventional superconductivity in sr2ruo4. Proc. Natl. Acad. Sci. 118, e2025313118 (2021).

    CAS  PubMed  Google Scholar 

  30. 30.

    Grinenko, V. et al. Split superconducting and time-reversal symmetry-breaking transitions in sr2ruo4 under stress. Nat. Phys. 17, 748–754 (2021).

    CAS  Google Scholar 

  31. 31.

    Benhabib, S. et al. Ultrasound evidence for a two-component superconducting order parameter in sr 2 ruo 4. Nat. Phys. 17, 194–198 (2021).

    CAS  Google Scholar 

  32. 32.

    Sigrist, M. & Rice, T. Symmetry classification of states in high temperature superconductors. Z. Phys. B Condens. Matter 68, 9–14 (1987).

    ADS  Google Scholar 

  33. 33.

    Tsuei, C. & Kirtley, J. Pairing symmetry in cuprate superconductors. Rev. Mod. Phys. 72, 969 (2000).

    ADS  CAS  Google Scholar 

  34. 34.

    Annett, J. F. Unconventional superconductivity. Contemp. Phys. 36, 423–437 (1995).

    ADS  CAS  Google Scholar 

  35. 35.

    Annett, J. F. Symmetry of the order parameter for high-temperature superconductivity. Adv. Phys. 39, 83–126 (1990).

    ADS  CAS  Google Scholar 

  36. 36.

    Gorkov, L. P. & Rashba, E. I. Superconducting 2d system with lifted spin degeneracy: Mixed singlet-triplet state. Phys. Rev. Lett. 87, 037004 (2001).

    ADS  CAS  Google Scholar 

  37. 37.

    Lepori, L., Giuliano, D., Nava, A. & Perroni, C. Interplay between singlet and triplet pairings in multi-band two-dimensional oxide superconductors. Preprint at (2021).

  38. 38.

    Salomaa, M. & Volovik, G. Quantized vortices in superfluid he3. Rev. Mod. Phys. 59, 533 (1987).

    ADS  CAS  Google Scholar 

  39. 39.

    Mineev, K. & Samokhin, K. Helical phases in superconductors. Zh. Eksp. Teor. Fiz 105, 747–763 (1994).

    CAS  Google Scholar 

  40. 40.

    Romano, A., Cuoco, M., Noce, C., Gentile, P. & Annunziata, G. Field-induced transition from chiral spin-triplet to mixed-parity Fulde–Ferrell–Larkin–Ovchinnikov superconductivity. Phys. Rev. B 81, 064513 (2010).

    ADS  Google Scholar 

  41. 41.

    Matsuo, S., Shimahara, H. & Nagai, K. Order parameter mixing effecting the Fulde–Ferrell state. J. Phys. Soc. Jpn. 63, 2499–2502 (1994).

    ADS  CAS  Google Scholar 

  42. 42.

    Fujimoto, S. Unambiguous probe of parity mixing of cooper pairs in noncentrosymmetric superconductors. Phys. Rev. B 79, 220506 (2009).

    ADS  Google Scholar 

  43. 43.

    Lebed, A. G. Cooper pairs with broken parity and spin-rotational symmetries in d-wave superconductors. Phys. Rev. Lett. 96, 037002 (2006).

    ADS  CAS  PubMed  Google Scholar 

  44. 44.

    Mizushima, T., Yamakage, A., Sato, M. & Tanaka, Y. Dirac-fermion-induced parity mixing in superconducting topological insulators. Phys. Rev. B 90, 184516 (2014).

    ADS  Google Scholar 

  45. 45.

    Möckli, D. & Khodas, M. Robust parity-mixed superconductivity in disordered monolayer transition metal dichalcogenides. Phys. Rev. B 98, 144518 (2018).

    ADS  Google Scholar 

  46. 46.

    Ekin, J. et al. Correlation between d-wave pairing behavior and magnetic-field-dependent zero-bias conductance peak. Phys. Rev. B 56, 13746 (1997).

    ADS  CAS  Google Scholar 

  47. 47.

    Kashiwaya, H. et al. Zeeman magnetic field responses of d-wave superconductors observed by tunneling spectroscopy. J. Phys. Chem. Solids 67, 350–352 (2006).

    ADS  CAS  Google Scholar 

  48. 48.

    Prohammer, M. & Carbotte, J. Upper critical field of s-and d-wave superconductors with anisotropic effective mass. Phys. Rev. B 42, 2032 (1990).

    ADS  CAS  Google Scholar 

  49. 49.

    Wang, Y.-L. et al. Parallel magnetic field suppresses dissipation in superconducting nanostrips. Proc. Natl. Acad. Sci. 114, E10274–E10280 (2017).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  50. 50.

    He, W.-Y. et al. Magnetic field driven nodal topological superconductivity in monolayer transition metal dichalcogenides. Commun. Phys. 1, 1–7 (2018).

    ADS  CAS  Google Scholar 

  51. 51.

    Nayak, S. & Kumar, S. Exotic superconducting states in the extended attractive Hubbard model. J. Phys. Condens. Matter 30, 135601 (2018).

    ADS  PubMed  Google Scholar 

  52. 52.

    Batra, N., Nayak, S. & Kumar, S. Topological transitions in a model for proximity-induced superconductivity. Phys. Rev. B 100, 214517 (2019).

    ADS  CAS  Google Scholar 

  53. 53.

    Sigrist, M. & Ueda, K. Phenomenological theory of unconventional superconductivity. Rev. Mod. Phys. 63, 239 (1991).

    ADS  CAS  Google Scholar 

Download references


We acknowledge the use of Computing Facility at IISER Mohali.

Author information




S.K. and S.N. conceived the project. S.N. and N.B. formulated the methodology and ran preliminary simulations for consistency. S.N. performed the simulations and arranged the results, as presented in the final manuscript. All authors contributed to analyzing the results and writing the final manuscript. S.K. supervised the project.

Corresponding author

Correspondence to Sanjeev Kumar.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Nayak, S., Batra, N. & Kumar, S. Pairing symmetries in the Zeeman-coupled extended attractive Hubbard model. Sci Rep 11, 22724 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.


Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing