Abstract
By introducing the possibility of equal and oppositespin 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) meanfield approach, and we use the dvector formalism to characterize the nature of the stabilized superconducting orders. We discover that, while all other types of orders are suppressed, a nonunitary triplet order dominates the phase space in the presence of an inplane external magnetic field. We also find a transition between a nonunitary 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.
Introduction
In a conventional swave superconductor, described by BCS–Migdal–Eliashberg theory^{1,2,3,4,5,6}, a phononmediated effective attraction causes the electrons to form spinsinglet (with total spin \(S=0\)) Cooper pairs with an isotropic swave orbital order parameter (OP) symmetry^{7}. Most of the superconductors discovered in the early phase of the past century have this type of conventional OP symmetry^{8}. On the contrary, the general consensus about the novel high \(T_{c}\) cuprate superconductors is that they are identified as strong candidates for unconventional dwave superconductors^{9,10,11}, which support the formation of spinsinglet Cooper pairs with an anisotropic dwave orbital OP symmetry. In general, Cooper pairs can also be formed in the spintriplet state, with total spin \(S=1\) and anisotropic orbital OP^{12,13}. Recent studies revealed that a two dimensional chiral pwave spintriplet superconductor is a candidate to host Majorana fermions^{14,15,16}. A pair of Majorana fermions, bound to topological defects, together known as Ising anyons, exhibit nonAbelian exchange statistics and such an object can be considered to be the potential building block for decoherence free quantum computation^{17}. Regardless of their technological relevance, superconductors with tripletsuperconductivity emerging due to the intrinsic properties of the material itself are quite rare in nature. While \(^{3}\text {He}\) was the first chargeless many body system where tripletpairing state was realised^{18,19,20,21}, it was \(\text {UPt}_{3}\) which was identified as the first spintriplet superconductor in charged many body systems^{22,23}. Strong candidate materials for tripletsuperconductivity, so far, include Uranium based heavyfermion superconductors \(\text {UPt}_{3}\), \(\text {UGe}_{2}\), URhGe, UCoGe^{24}, organic superconductor \(\text {(TMTSF)}_{2}\text {PF}_{6}\)^{25} and probably recently discovered heavyfermion superconductor \(\text {UTe}_{2}\)^{26,27,28}. We would also like to add that Strontium Ruthenate (Sr_{2}RuO_{4}), which for its most part was understood as a strong candidate for oddparity superconductor, has undergone immense scrutiny. The exact nature of the SC order in Sr_{2}RuO_{4} has remained an open question for many years now, with a recent study ruling out the existence of a pure oddparity superconducting order, and explaining that the Knight shifts originate from quasiparticles rather than the condensate itself^{29}. Other studies impose constraints on the possible order of Sr_{2}RuO_{4} superconducting state—stressinduced splitting of onset temperatures for SC and timereversal symmetry (TRS) breaking has also been reported^{30}, and ultrasound measurements suggests the existence of a multicomponent order parameter^{31}.
In conventional superconductors TRS allows energetically degenerate electron states \(\mathbf {k}\uparrow \rangle\) and \(\mathbf {k}\downarrow \rangle\) to form spinsinglet Cooper pairs. On the other hand, parity symmetry (inversion symmetry in this scenario) along with TRS is needed to form oddparity spintriplet Cooper pairs^{32,33,34,35}. In case of broken inversion symmetry, there is also the possibility of forming Cooper pairs with mixedparity superconducting (SC) states^{32,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 parameter^{40,41}, twodimensional superconductors with broken inversion symmetry^{36}, noncentrosymmetric superconductors^{42}, vortex phase of a dwave superconductor in presence of paramagnetic effects^{43}, SC topological insulators in presence of surface Dirac fermions^{44}, and disordered monolayer transition metal dichalcogenides^{45}. 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 spintriplet and mixedparity 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 swave superconductors, \(d_{x^2y^2}\)wave superconductors exhibit gapless quasiparticle spectrum, which allows even a weak amount of Zeeman field to spinpolarize the quasiparticle excitations, resulting in a stronger response to an external magnetic field^{46,47,48}. For a two dimensional system it is quite interesting to study the effect of an inplane (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 inplane magnetic field has been studied in various contexts, including quasi2D 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 dichalcogenides^{49,50}.
In this article we provide a detailed study of the competition between different superconducting pairings and their response to nonzero temperature and Zeeman field. We consider an extended attractive Hubbard model, with onsite and nearest neighbour attractive interaction on a square lattice, and treat the manybody interaction terms using meanfield 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 dvector 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 selfconsistently 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 nonzero magnetic field to lift the spin degeneracy. In “Ground State phase diagrams” section, we present the ground state phase diagram of nearestneighbour 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 nonzero 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 dvector 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 dvector formalism and its usefulness in characterizing triplet order parameters.
Results and discussion
For superconductors like hightemperature Cuprates, the electronic band is quite narrow i.e., the orbitals have a small overlap between adjacent atoms. We therefore make use of the tightbinding 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 twodimensional square lattice can be written as,
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 onsite attractive interaction between two different spin projections, \(\mathcal {H}_{\text {int}}^{\text {nn}}\) is the nearestneighbor (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,
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 nearestneighbor sites of the square lattice. The operator \(n_{i\sigma }\) is the electron occupation number operator at \(i{\text {th}}\) site for \(\sigma\) spinprojection. 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 onsite attractive interaction strength, whereas V is the intersite attractive interaction strength. The zcomponent of the magnetic field is given by B, which lies in the lattice plane. t is the usual nearestneighbor 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:
where we rearranged the Fermionic operators in particle–particle channel or so called pairing channel, using anticommutation rules. Similarly for the intersite term:
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 manyparticle interaction which in itself is very hard to analyse. We therefore make the meanfield approximation and via the Bogoliubovde Gennes framework, we solve the Hamiltonian computationally, and present the results in the following sections. The details to the BdG framework are given in “Bogoliubovde 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.
\(p_{x}\) state (and correspondingly \(p_y\) state): \(\Delta _{x}^{+}=\Delta _{x}^{} \in {\mathbb {C}};~~~~~~\Delta _{y}^{+}= \Delta _{y}^{}=0\).

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 selfconsistency 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.
\(\Delta ^{\uparrow \uparrow }(k)=2iV\left[ \Delta _{x}^{\uparrow }\sin (kx)+ \Delta _{y}^{\uparrow }\sin (ky)\right]\).

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 selfconsistent 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\).
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 inplane 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 nearestneighbor 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 puretriplet (\(p_{x}\pm ip_{y}\) type) state being stabilized in our previous study^{51}. As the interaction strength V is kept fixed while increasing nearestneighbor hopping amplitude t, effectively we are moving from a strongcoupling limit to a weakcoupling limit.
In the strongcoupling limit, where the electronpairs are tightly bound, we can appropriately apply the wellknown physics of spinsinglet and spintriplet pairing to the spin degrees of freedom associated with the electronpairs. 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 weakcoupling limit, the system gains energy via delocalization of electrons, and characterization of different triplet states in terms of local spin operators is not valid. Note that we use the phaselock 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 nearestneighbour hopping amplitude t we note that the OSP state (\(\uparrow \downarrow\) in figure) becomes energetically favourable. Furthermore, at any finite value of Zeemancoupling, \(\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 Zeemancoupling, 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 Zeemancoupling, OSP state is completely disfavoured in the weakcoupling limit (Fig.1d). One should take a note of the fact that this particular observation is validated by selfconsistent 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 Zeemancoupling in the two limiting cases (\(t = 0\) and \(t=(1/1.8) V\)). At strongcoupling (\(t=0\)) limit, OSP state is of lower energy compared to \(\downarrow \downarrow\)type ESP state below a critical Zeemancoupling (\(B\approx (2/1.8) V\)). Beyond such critical value of Zeemancoupling OSP state becomes completely disfavoured energetically with respect to the ESP states. On the other hand, at weakcoupling limit (\(t=(1/1.8) V\)), energy of OSP state remains comparable with \(\downarrow \downarrow\)type ESP state at moderate values of Zeemancoupling, and eventually becomes energetically disfavoured at higher values of Zeemancoupling. 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 Zeemancoupling.
Ground state phase diagrams
In the previous work^{51} 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^2y^2}} \Delta _{p_{x}},\Delta _{p_{y}}\) }, which were obtained using the standard definition of superconducting OPs, defined in literature^{51}, from the original set of pairing correlations { \(\Delta , \Delta _{\delta }^{+}, \Delta _{\delta }^{}\) }. Explicitly writing {\(\Delta\)} we have,
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 ESPtype superconducting solutions along with OSPtype 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
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 OSPtype triplet states to exist, while the ESP group contains pairing correlations (\(\Delta _{\delta }^{\uparrow }\),\(\Delta _{\delta }^{\downarrow }\)) that give rise to ESPtype 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 OSPtype 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 ESPtype triplet states to be stabilized because of their symmetrizedspin wavefunction (and therefore the orbital part of the wavefunction which is the superconducting gap function, should be antisymmetric) 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 selfconsistent 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 ESPtype 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 dvector formalism, defined in the “dvector 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:
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 intersite attraction V. We keep onsite attraction at a fixed value (\(U=t\)), as it is the intersite attraction V that plays a role in stabilizing different unconventional superconducting orders as deduced from our previous works^{51}. When intersite attraction V is small (\(V=0.4t\)) compared to onsite 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 halffilling (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 halffilling is \(d_{x^2y^2}\) type, while the singlet phase near lowdensity region is of \(s+s^{*}\) type based off our previous work^{51}. Pure triplet phase occurs near quarterfilling, while singlet–triplet mixed parity phase occurs between pure triplet phase and singlet phase of \(d_{x^2y^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^2y^2}\) type singlet phase occurs at lowV region, when V is still larger than onsite interaction U. At higher value of intersite attraction (\(V=3t\)), pure\(d_{x^2y^2}\) type singlet phase vanishes (Fig. 2c), while singlet–triplet mixed phase occupies the corresponding density region (near halffilling). 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 intersite interaction V is less than the onsite 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 onsite attraction is known to stabilize singlet swave SC order. While U is stabilizing singlet swave SC order, V is also playing an important role in stabilizing the singlet \(s^{*}\) and the unconventional \(d_{x^2y^2}\) SC orders in this region along with swave. When the lower density region prefers singlet s and \(s^{*}\) SC orders to other phases, near halffilling it is the singlet \(d_{x^2y^2}\)wave that prevails. Notice that these inferences come from a combination of the present work and the previous work^{51}. As intersite interaction V becomes larger compared to onsite interaction U, it not only stabilizes singlet \(s^{*}\) and \(d_{x^2y^2}\) SC orders, it also helps pure triplet and mixedparity phases to get stabilized at different density regimes. The pure triplet phase (green in color) occurs in a deltashaped region near quarter filling (\(0.4 \le \langle n\rangle \le 0.6\)), when the intersite attraction V is higher than the onsite attraction U but not as large as \(V=4t\). With V getting larger, the mixedparity (blue in color) phase becomes more stable, specially in the higher density region. It’s interesting to note that \(d_{x^2y^2}\)type singlet phase is prone to get stabilized near halffilling region and pure triplet phase has a tendency to get stabilized near quarterfilling. But depending on the strength of intersite interaction the system gains energy by stabilizing a new mixedparity (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 study^{51}. 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 intersite attractive interactions occurs via \(tJ\) 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 meanfield approach.
Finite temperature study
After exploring the stability of singlet, triplet, and mixedparity 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 wellknown 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 meanfield 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 mixedparity state exactly at halffilling \(\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 nonSC state. At \(T\approx 0.38\) the system makes a transition from a mixedparity SC state to a pure \(d_{x^2y^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 mixedparity superconducting state to a pure singlet superconducting state before making a transition to a nonsuperconducting 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 nonsuperconducting 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 nonsuperconducting phase when temperature is increased. While the transition temperature for singlet phase takes a domeshaped feature as a function of density, for triplet and mixedparity 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 nonsuperconducting state retaining the specific form of pairing symmetry, while the mixedparity 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 nonsuperconducting 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 Zeemancoupling 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 inplane 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 nonsuperconducting region (gray in color), which is quite expected as the magnetic field destroys superconductivity. It is interesting to note that the nonsuperconducting region appears at lower values of V suggesting that for the superconductivity to survive it really is a competition between the intersite interaction strength and the magnetic field strength.
From the ground state V\(\langle n\rangle\) phase diagram in the absence of a magnetic field Fig. 2e, we notice that when intersite interaction V is less than the onsite 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^2y^2}\) type singlet phase, which is stable near halffilling. It appears that presence of an external magnetic field of strength \(B = 0.1t\), destroys such pure \(d_{x^2y^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 Zeemancoupling 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 Zeemancoupling, (a) \(B=0\), and (b) \(B=0.1t\) presented in Fig. 4. For the purpose of the abovementioned illustration we keep the onsite interaction at \(U=t\), while the intersite interaction is kept fixed at \(V=4t\). While in the absence of a Zeemancoupling we see a continuous variation of \(\langle n\rangle\) with \(\mu\), in the presence of a Zeemancoupling, 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 Zeemancoupling 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 mixedparity 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 Zeemancoupling 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 mixedparity 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 nonsuperconducting 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 Zeemancoupling to an external magnetic field at low strength favours pure triplet superconductivity. Secondly, the disappearance of pure singlet and mixedparity 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^2y^2}\) etc., and ESP type states give rise to pure triplet superconducting phases, the abovementioned behavior is wellexplained.
We conclude this section by looking at the effect of Zeemancoupling at finite temperatures. Figure 3d shows the BT 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 distinctsymmetry 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 mixedparity 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 swave, \(s^{*}\)wave, \(p_{x}+ip_{y}\)wave and \(d_{x^2y^2}\)wave to define the nature of the orbital part of the gap function^{51,52}. In the case of pure singlet, or pure triplet phase the nature of the orbital part can be easily understood, while for mixedparity 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 onsite attraction potential U and the intersite attraction potential V are kept fixed at \(U = t\) and \(V = 2t\), for all the plots in Fig. 5.
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 mixedparity solution around half filling. Further characterization of the singlet phases, based on the definitions of swave, \(s*\)wave, and dwave are shown in the inset of Fig. 5b. Note that near half filling, where the singlet phase shows dtype pairing (inset of Fig. 5b), it is a mixed parity state with \(d_{x^2y^2}\) as a singlet component. The presence of swave and \(s^{*}\)wave character in pure singlet phase and \(d_{x^{2}y^{2}}\)wave character in the mixedparity phase is consistent with our previous findings when the Zeeman field is absent^{51}. 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 puresinglet, puretriplet, and mixedparity 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 (dvector formalism also allows us to probe the time reversal symmetry in a superconducting phase, refer “dvector formalism” section). Notice that even though Fig. 5a is generated for nonzero magnetic field \(B=0.12t\), we see the existence of time reversal invariant (TRPSC) superconducting phase appearing for lower densities. This is the enduring swave 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 nonsuperconducting (NS) region. Fig. 5a,b shows that, in the presence of a magnetic field, it is the ESPtype 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 ESPtype 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^2y^2}\) singlet component becomes dominant near half filling. The character of \(d_{x^2y^2}\) breaks the symmetry between the x and y components of the correlation functions, and as \(d_{x^2y^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 particlehole symmetric, at halffilling 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 spinresolved 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})\)vector^{13} 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 nonunitarity of a SC order is usually associated with the breaking of time reversal symmetry (TRS)^{53}. Nonunitarity 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 spinresolved 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 nonunitary (unitary when \(\Delta _\mathbf{q}=0\)) in nature when the triplet order parameter is finite. Note that nonunitarity of triplet gap function is defined in “dvector formalism” section. From Fig. 5e,f, we see that nonunitarity 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 nonunitary phase to exist and as we go close to halffilling, the nonunitarity is destroyed, consequently, the triplet superconducting phase becomes unitary exactly at halffilling 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 pseudospin 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 noninteracting 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^2y^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 Binduced difference in the shapes of the up and the down spin Fermi surfaces.
Conclusions
Using the unrestricted BdG meanfield 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, nonunitary, 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 antialigned, 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 Zeemancoupling, OSP state is completely disfavoured in the weakcoupling limit. The \(\uparrow \uparrow\) ESP state is favoured in the presence of Zeeman field. We use the dvector 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 halffilling. 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^2y^2}\)type, and finally to a nonsuperconducting 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 twodimensional systems, such as atomically thin conducting interfaces between insulating oxides.
Methods
Meanfield decoupling in pairing channel
Both the interaction terms, the onsite interaction and the intersite interaction term, are represented by twobody operator terms in the second quantization formalism. To reduce the complexity of these terms, we consider the meanfield treatment, where the manybody interaction effects are mimicked via the interaction of a single electron system with a meanfield created by the aggregated effect of rest of the electrons. Mathematically speaking, such a twobody operator term \({\hat{X}}{\hat{Y}}\) can be approximated as
where we have neglected the second order correction contribution. In practice, we can apply this idea of meanfield 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 meanfields of the form \(\langle c^{\dagger }c\rangle\), and pairing channel decoupling leads to meanfields 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 nonzero 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:
where, \(\Delta _{i} = \langle c_{i\downarrow }c_{i\uparrow }\rangle\) is defined as the onsite 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:
where, we have defined the intersite pairing correlations as:
After combining all the terms together, our effective meanfield Hamiltonian in real space looks like:
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 meanfield Hamiltonian in real space Eq. (16) it is apparent that the effective Hamiltonian is a function of

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

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:
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 meanfield Hamiltonian looks like:
where,
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 spindegeneracy in the system via Zeeman coupling. It’s worth mentioning that we work in the regime where magnetic field couples only with spindegree 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 onsite attraction strength U controls OSP states only, which is evident from the definitions of these pairing correlations Eq. (19), intersite attraction V controls both OSP and ESP states. This very fact suggests that V is meant to play an important role in stabilizing nontrivial superconducting orders with variety of pairing symmetries.
Bogoliubovde Gennes method
The mean field Hamiltonian \(\mathcal {H}_{\text {MF}}\) Eq. (18) can be written in terms of Nambu spinors \({\Psi }_{\mathbf {k}}\) as:
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:
\(\mathcal {H}_{\text {BdG}}(\mathbf {k})\) is known as the mean field Bogoliubovde Gennes Hamiltonian. We can diagonalize the BdG Hamiltonian by defining new Fermionic quasiparticle operators that mix the electronic operators as,
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,
here \(\{E_{\mathbf {k}\alpha }\}\) is the excitation spectrum for the Bogoliubov quasiparticles. The physical constraint of nonnegativity 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 pseudospin \(\alpha\). The BdG quasiparticles describe elementary excitations of the condensate, \(E_{g}\) being the ground state energy of the condensate.
dvector formalism
In the absence of magnetic field, all the triplet superconducting solutions should be degenerate. However, while performing selfconsistency, 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 selfconsistency 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 dvector formalism.
The most general superconducting gap function can be represented in a matrix form as,
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,
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 threecomponents 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 spinquantization axis, we may use this quantity to track down the overall triplet component of superconductivity. We therefore use,
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 nonunitary. With this definition, it is clear that only triplet pairing matrices can be nonunitary since,
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,
Therefore, using the dvector approach we can infer many properties of the superconducting wavefunction.
Selfconsistency and minimization
In this section we lay out how the method of selfconsistency 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 selfconsistent 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 rediagonalized. 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 selfconsistently is actually equivalent to the energy minimization. We use \(64 \times 64\) kgrid for \(V \ge 2t\), and a \(512 \times 512\) kgrid for the cases \(V < 2t\).
Code availability
All the codes are written from scratch, in Python or Fortran. Plots are generated using Python scripts.
References
 1.
Cooper, L. N. Bound electron pairs in a degenerate fermi gas. Phys. Rev. 104, 1189 (1956).
 2.
Bardeen, J., Cooper, L. N. & Schrieffer, J. R. Theory of superconductivity. Phys. Rev. 108, 1175 (1957).
 3.
Bardeen, J., Cooper, L. N. & Schrieffer, J. R. Microscopic theory of superconductivity. Phys. Rev. 106, 162 (1957).
 4.
Migdal, A. Interaction between electrons and lattice vibrations in a normal metal. Sov. Phys. JETP 7, 996–1001 (1958).
 5.
Eliashberg, G. Interactions between electrons and lattice vibrations in a superconductor. Sov. Phys. JETP 11, 696–702 (1960).
 6.
Eliashberg, G. Temperature green function for electrons in a superconductor. Sov. Phys. JETP 12, 1000–1002 (1961).
 7.
Tinkham, M. Introduction to Superconductivity (Courier Corporation, 2004).
 8.
Bennemann, K. & Ketterson, J. B. History of superconductivity: Conventional, hightransition temperature and novel superconductors. In Superconductivity (eds Bennemann, K. & Ketterson, J. B.) 3–26 (Springer, 2008).
 9.
Coffey, D. & Coffey, L. Quasiparticle decay effects in the superconducting density of states: Evidence for dwave pairing in the cuprates. Phys. Rev. Lett. 70, 1529 (1993).
 10.
Tsuei, C. & Kirtley, J. Phasesensitive evidence for dwave pairing symmetry in electrondoped cuprate superconductors. Phys. Rev. Lett. 85, 182 (2000).
 11.
Oates, D., Park, S.H. & Koren, G. Observation of the nonlinear Meissner effect in ybco thin films: Evidence for a dwave order parameter in the bulk of the cuprate superconductors. Phys. Rev. Lett. 93, 197001 (2004).
 12.
Maki, K., Kee, H. Y. & Morita, Y. Triplet superconductivity in a nutshell. J. Supercond. Novel Magn. 22, 71–74 (2009).
 13.
Mackenzie, A. P. & Maeno, Y. The superconductivity of sr 2 ruo 4 and the physics of spintriplet pairing. Rev. Mod. Phys. 75, 657 (2003).
 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).
 15.
Volovik, G. Fermion zero modes on vortices in chiral superconductors. J. Exp. Theor. Phys. Lett. 70, 609–614 (1999).
 16.
Sarma, S. D., Nayak, C. & Tewari, S. Proposal to stabilize and detect halfquantum vortices in strontium ruthenate thin films: Nonabelian braiding statistics of vortices in a p x+ i p y superconductor. Phys. Rev. B 73, 220502 (2006).
 17.
Beenakker, C. Search for majorana fermions in superconductors. Annu. Rev. Condens. Matter Phys. 4, 113–136 (2013).
 18.
Thouless, D. J. Perturbation theory in statistical mechanics and the theory of superconductivity. Ann. Phys. 10, 553–588 (1960).
 19.
Brueckner, K., Soda, T., Anderson, P. W. & Morel, P. Level structure of nuclear matter and liquid he 3. Phys. Rev. 118, 1442 (1960).
 20.
Anderson, P. W. & Morel, P. Generalized Bardeen–Cooper–Schrieffer states and the proposed lowtemperature phase of liquid he3. Phys. Rev. 123, 1911 (1961).
 21.
Leggett, A. J. A theoretical description of the new phases of liquid he3. Rev. Mod. Phys. 47, 331 (1975).
 22.
Tou, H. et al. Oddparity superconductivity with parallel spin pairing in upt 3: Evidence from 195 pt knight shift study. Phys. Rev. Lett. 77, 1374 (1996).
 23.
Tou, H. et al. Nonunitary spintriplet superconductivity in upt 3: Evidence from p 195 t knight shift study. Phys. Rev. Lett. 80, 3129 (1998).
 24.
Aoki, D., Ishida, K. & Flouquet, J. Review of ubased ferromagnetic superconductors: Comparison between uge2, urhge, and ucoge. J. Phys. Soc. Jpn. 88, 022001 (2019).
 25.
Lee, I. et al. Triplet superconductivity in an organic superconductor probed by nmr knight shift. Phys. Rev. Lett. 88, 017004 (2001).
 26.
Ran, S. et al. Nearly ferromagnetic spintriplet superconductivity. Science 365, 684–687 (2019).
 27.
Metz, T. et al. Pointnode gap structure of the spintriplet superconductor ute 2. Phys. Rev. B 100, 220504 (2019).
 28.
Sundar, S. et al. Coexistence of ferromagnetic fluctuations and superconductivity in the actinide superconductor ute 2. Phys. Rev. B 100, 140502 (2019).
 29.
Chronister, A. et al. Evidence for even parity unconventional superconductivity in sr2ruo4. Proc. Natl. Acad. Sci. 118, e2025313118 (2021).
 30.
Grinenko, V. et al. Split superconducting and timereversal symmetrybreaking transitions in sr2ruo4 under stress. Nat. Phys. 17, 748–754 (2021).
 31.
Benhabib, S. et al. Ultrasound evidence for a twocomponent superconducting order parameter in sr 2 ruo 4. Nat. Phys. 17, 194–198 (2021).
 32.
Sigrist, M. & Rice, T. Symmetry classification of states in high temperature superconductors. Z. Phys. B Condens. Matter 68, 9–14 (1987).
 33.
Tsuei, C. & Kirtley, J. Pairing symmetry in cuprate superconductors. Rev. Mod. Phys. 72, 969 (2000).
 34.
Annett, J. F. Unconventional superconductivity. Contemp. Phys. 36, 423–437 (1995).
 35.
Annett, J. F. Symmetry of the order parameter for hightemperature superconductivity. Adv. Phys. 39, 83–126 (1990).
 36.
Gorkov, L. P. & Rashba, E. I. Superconducting 2d system with lifted spin degeneracy: Mixed singlettriplet state. Phys. Rev. Lett. 87, 037004 (2001).
 37.
Lepori, L., Giuliano, D., Nava, A. & Perroni, C. Interplay between singlet and triplet pairings in multiband twodimensional oxide superconductors. Preprint at http://arxiv.org/abs/2107.01100 (2021).
 38.
Salomaa, M. & Volovik, G. Quantized vortices in superfluid he3. Rev. Mod. Phys. 59, 533 (1987).
 39.
Mineev, K. & Samokhin, K. Helical phases in superconductors. Zh. Eksp. Teor. Fiz 105, 747–763 (1994).
 40.
Romano, A., Cuoco, M., Noce, C., Gentile, P. & Annunziata, G. Fieldinduced transition from chiral spintriplet to mixedparity Fulde–Ferrell–Larkin–Ovchinnikov superconductivity. Phys. Rev. B 81, 064513 (2010).
 41.
Matsuo, S., Shimahara, H. & Nagai, K. Order parameter mixing effecting the Fulde–Ferrell state. J. Phys. Soc. Jpn. 63, 2499–2502 (1994).
 42.
Fujimoto, S. Unambiguous probe of parity mixing of cooper pairs in noncentrosymmetric superconductors. Phys. Rev. B 79, 220506 (2009).
 43.
Lebed, A. G. Cooper pairs with broken parity and spinrotational symmetries in dwave superconductors. Phys. Rev. Lett. 96, 037002 (2006).
 44.
Mizushima, T., Yamakage, A., Sato, M. & Tanaka, Y. Diracfermioninduced parity mixing in superconducting topological insulators. Phys. Rev. B 90, 184516 (2014).
 45.
Möckli, D. & Khodas, M. Robust paritymixed superconductivity in disordered monolayer transition metal dichalcogenides. Phys. Rev. B 98, 144518 (2018).
 46.
Ekin, J. et al. Correlation between dwave pairing behavior and magneticfielddependent zerobias conductance peak. Phys. Rev. B 56, 13746 (1997).
 47.
Kashiwaya, H. et al. Zeeman magnetic field responses of dwave superconductors observed by tunneling spectroscopy. J. Phys. Chem. Solids 67, 350–352 (2006).
 48.
Prohammer, M. & Carbotte, J. Upper critical field of sand dwave superconductors with anisotropic effective mass. Phys. Rev. B 42, 2032 (1990).
 49.
Wang, Y.L. et al. Parallel magnetic field suppresses dissipation in superconducting nanostrips. Proc. Natl. Acad. Sci. 114, E10274–E10280 (2017).
 50.
He, W.Y. et al. Magnetic field driven nodal topological superconductivity in monolayer transition metal dichalcogenides. Commun. Phys. 1, 1–7 (2018).
 51.
Nayak, S. & Kumar, S. Exotic superconducting states in the extended attractive Hubbard model. J. Phys. Condens. Matter 30, 135601 (2018).
 52.
Batra, N., Nayak, S. & Kumar, S. Topological transitions in a model for proximityinduced superconductivity. Phys. Rev. B 100, 214517 (2019).
 53.
Sigrist, M. & Ueda, K. Phenomenological theory of unconventional superconductivity. Rev. Mod. Phys. 63, 239 (1991).
Acknowledgements
We acknowledge the use of Computing Facility at IISER Mohali.
Author information
Affiliations
Contributions
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
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 http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Nayak, S., Batra, N. & Kumar, S. Pairing symmetries in the Zeemancoupled extended attractive Hubbard model. Sci Rep 11, 22724 (2021). https://doi.org/10.1038/s41598021021755
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41598021021755
Comments
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.