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.


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, where, H kin represents the kinetic energy of the electrons, H µ is the chemical potential term, H onsite int is the on-site attractive interaction between two different spin projections, H nn int is the nearest-neighbor (nn) attractive interaction which would be responsible for all the unconventional superconducting phases, and finally H 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σ (c iσ ) creates (annihilates) an electron with spin σ at ith site of the lattice. The sum over ij represents the sum over the nearest-neighbor sites of the square lattice. The operator n iσ is the electron occupation number operator at ith site for σ spin-projection. The total electron occupation number at ith site is, therefore, defined as n i = n i↑ + n i↓ . The parameter U is the on-site attractive interaction strength, whereas V is the intersite 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, µ 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 anti-commutation rules. Similarly for the inter-site term: where δ is an index which denotes unit vectors along the direction of two nearest neighbors, namely δ = +x and +ŷ . 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 = ±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 � ↑↓ (k). 1. p x state (and correspondingly p y state): In general, p x and p y can be of the form: p x + e iθ p y . With the help of energetics and self-consistency we have checked that θ = π 2 − π 2 corresponds to the most stable solution, i.e. it takes the form p x ± 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 π/2 . This allows us to reduce our exploration space { } . Similarly, we fix the phase relation among, + x , − x , + y and − y .
Triplet combinations from � ↑↑ (k) and � ↓↓ (k).   (n i↑ n i+δ↑ + n i↓ n i+δ↓ + n i↑ n i+δ↓ + n i↓ n i+δ↑ ) . We find that � 1 = π 2 minimizes the average energy E (Fig. 1e) independent of the external parameters. Note that we also allowed ↑ and ↓ parring correlation functions to have a relative phase of 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 { ↑ x , ↑ y } (and { ↓ x , ↓ y } ) and we cannot, at this stage, comment on the phase difference between � ↑↑ (k) and � ↓↓ (k) because these results are independent of 2 .
Competition among the different triplet states. In this section we compare the energies of the three triplet states: We expect all three states |S z = 0� , |S z = 1� , and |S z = −1� 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� or |S z = −1� 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� triplet component enter into the Hamiltonian through the superconducting pairing correlation � ↑↓ (k) ; and |S z = 1� , |S z = −1� triplet contributions enter through the triplet superconducting correlations � ↑↑ (k) and � ↓↓ (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 , µ = −(1.5/1.8)V ) where we found pure-triplet ( p x ± ip y type) state being stabilized in our previous study 51 . 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.
. The variation shown here is independent of the choice of external parameters like the chemical potential µ , hopping parameter t, magnetic field B , and the relative phase between � ↑ δ and � ↓ δ i.e., 2 . This suggests that, independent of the external parameters, the most stable relative phase between � σ x and � σ y is � 1 = π/2 . (f,g) The plots shows the variation of average energy E for the three triplet states with the magnetic field B, for (f) t = 0 and (g) t = (1/1.8)V. www.nature.com/scientificreports/ 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 ( ↑↓ in figure) becomes energetically favourable. Furthermore, at any finite value of Zeeman-coupling, ↑↑-type ESP state remain energetically favourable as we expected, while the energy of ↓↓-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 ↓↓-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 ↓↓-type ESP state below a critical Zeeman-coupling ( B ≈ (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 ↓↓-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 work 51  Note that δ is an index for the unit vectors δ =x,ŷ along two independent directions on a square lattice. So { �, � + δ , � − δ } is actually a shorthand notation for { , + x , + y , − x , − 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 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 ( ,� + δ ,� − δ ) that allow the OSP-type triplet states to exist, while the ESP group contains pairing correlations ( � ↑ δ ,� ↓ δ ) 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 µ . Average electron density per site n is calculated for each value of µ . This allows one to plot magnitudes of the order parameters, defined in Eq. (10), with varying n . In the self-consistent approach, if we start with a random initial configuration of { } 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 n . 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 µ . In this case the d-vector formalism, defined in the "d-vector formalism" section, comes to our rescue. We find magnitudes of d 0 (k) and (9) www.nature.com/scientificreports/ d(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 n . In Fig. 2a-d, we plot the variation of singlet order parameter SP and triplet order parameter TP with n 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 works 51 . 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 work 51 . 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-n 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 www.nature.com/scientificreports/ 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 work 51 . 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 ≤ �n� ≤ 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 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 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.
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 SP and triplet TP order parameter with temperature for two different points in parameter space taken from the ground state V-n phase diagram. In Fig. 2f, we start with a mixed-parity state exactly at half-filling �n� = 1 . While the red dashed line represents the singlet order parameter SP , the green dotted line represents triplet order parameter TP . As temperature is increased, the triplet order parameter starts to fall off and around T ≈ 0.38 it vanishes completely. Interestingly, the singlet order parameter SP seems to remain constant, and begins to reduce only when 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 ≈ 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 SP starts to fall off and it completely vanishes at T ≈ 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-n phase diagram. In this case our zero temperature initial choice is at density �n� = 0.44 , where we find pure triplet state being stable. Notably with the increase of temperature the triplet order parameter TP falls off exactly the way it did for the previous case and the temperature where TP completely vanishes is exactly same as the previous case, i.e., T ≈ 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-n phase diagram to complete the analysis of the finite temperature effects. In the T-n 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 ≤ �n� ≤ 1.0 makes a transition to singlet superconducting state before making a transition to non-superconducting state. www.nature.com/scientificreports/ 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 inplane magnetic field serves our purpose. Motivated by the zero temperature V-n 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-n 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.
From the ground state V-n 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 ≈ 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-n 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 n with the chemical potential µ 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 n with µ , in the presence of a Zeeman-coupling, we clearly notice a discontinuity in the electron density at µ ≈ −2.5t . The discontinuity in the density is of the order δ �n� ≈ 0.1 , at V = 4t (Fig. 3b). www.nature.com/scientificreports/ To completely understand the effect of Zeeman-coupling on the superconducting phases, we present, in Fig. 3c, a zero temperature B-n 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-n 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 ∼ 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 µ = 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 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 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 n . 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.
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 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 SP and triplet order parameter 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  www.nature.com/scientificreports/ 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, Fig. 5a,b respectively. we observe that in the pure triplet state, near quarter filling,  www.nature.com/scientificreports/ 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 ↑ and ↓ 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 q vector averaged over the Brillouin zone, q ≡ k |q| . Here q-vector is defined in terms of the d-vector as q = i d × d * . The q(k)-vector 13 physically represents the net spin average present in the pairing state with momentum 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 ↑-spin electrons is different than that of the ↓-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 q vector, q tells us that the gap function is non-unitary (unitary when 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 �n� ≈ 0.6 . In the presence of the Zeeman field, the non-interacting Fermi surface splits into two, resulting in separate Fermi surfaces for ↑-spin (see Fig. 6c) and ↓-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 × 4 matrix at a given k F -point leads to two pairs of ±E(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 ↑-spin Fermi surface. Similarly, in Fig. 6d,e we plot the magnitude of gaps at ↓-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 ± ip 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 ± ip 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.

Conclusions
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, ↓↓ 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 ↑↑ 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.

Methods
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 XŶ can be approximated as www.nature.com/scientificreports/ 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 mean-fields of the form c † c , and pairing channel decoupling leads to mean-fields of the form c † c † . 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 X = c † i↑ c † i↓ and Ŷ = c i↓ c i↑ , we get: After combining all the terms together, our effective mean-field 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. Similarly ǫ ↓ (k) describes the kinetic energy of ↓-spin electrons with respect to an effective chemical potential µ ↓ eff = µ − 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) � ↑↑ (k) and � ↓↓ (k) represent superconducting correlations corresponding to equal spin pairing (ESP) states, while � ↑↓ (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 H MF Eq. (18) can be written in terms of
Nambu spinors k as: where, Nambu spinor k is a column vector of the form, � k = (c k↑ c † −k↓ c k↓ c † −k↑ ) and H BdG (k) is a 4 × 4 Hamiltonian matrix of the form: H BdG (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, 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 kα } 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 kα } to be a set of positive eigenvalues, and γ † kα (γ kα ) creates (annihilates) a Bogoliubov quasiparticle with momentum k and pseudo-spin α . 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.  where, ψ σ σ ′ (k) is the Cooper pair wavefunction with σ σ ′ pairing. The first term, � 0 (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 d(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 d(k) vector. A rotation of d(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 d(k) vector (and averaged over 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, 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, �(k) is called unitary if the product �(k)� † (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, where q = i(d × d * ) . The measure of q ≡ k |q| tells us if the superconducting phase is unitary. Also, the time reversal symmetry can easily be checked by, 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., µ, U, V , B . Corresponding to these, we will obtain a set of self-consistent SC pairing correlations, { } that defines the SC OP. To solve the BdG equations numerically, an initial guess of { } 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 { } 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 × 64 k-grid for V ≥ 2t , and a 512 × 512 k-grid for the cases V < 2t. (25) ψ σ σ ′ (k) = [� 0 (k) + d(k) · σ ]i(σ y ) σ σ ′ ,