Mixed temperature-dependent order parameters in the extended Hubbard model

TL;DR

This paper investigates how the symmetry of superconducting order parameters in the extended Hubbard model varies with temperature, revealing symmetry transitions below the critical temperature that impact the understanding of unconventional superconductors.

Contribution

It provides a detailed analysis of the temperature dependence of order parameter symmetries and introduces conditions on the temperature profile of the gap function in the extended Hubbard model.

Findings

Symmetry of the order parameter depends on temperature.

Existence of symmetry transitions below T_c.

Formulation of conditions on the temperature profile of the gap.

Abstract

The extended Hubbard model can host s-wave, d-wave and p-wave superconducting phases depending on the values of the on-site and nearest-neighbour interactions. Upon detailed examination of the free energy functional of the gap in this model, we show that these symmetries are often dependent on temperature. The critical points of this functional are constrained by symmetry and allow us to formulate stringent conditions on the temperature profile of the gap function, applicable to other models as well. We discuss the finite temperature phase diagram of the extended Hubbard model, and point out the existence of symmetry transitions below $T_c$. Understanding the nature of these transitions is crucial to assessing the symmetry of unconventional superconductors.

Equations

$$H$$

$${\Delta }_{\mathbf{k}}^{\mathrm{sing}}=\sum _{i\in \{0,s,d\}}{\Delta }_i{x}_{\mathbf{k}}^i;{\Delta }_{\mathbf{k}}^{\mathrm{trip}}=\sum _{i\in \{x,y\}}{\Delta }_i^{\mathrm{trip}}{x}_{\mathbf{k}}^i.$$

$$x_{\mathbf{k}}^0$$

$$x_{\mathbf{k}}^s$$

$$x_{\mathbf{k}}^d$$

$$x_{\mathbf{k}}^x$$

$$x_{\mathbf{k}}^y$$

$$H$$

$$V_{\uparrow \downarrow \downarrow \uparrow }=V_{\downarrow \uparrow \uparrow \downarrow }$$

$$V_{\uparrow \downarrow \uparrow \downarrow }=V_{\downarrow \uparrow \downarrow \uparrow }$$

$$V_{\uparrow \uparrow \uparrow \uparrow }=V_{\downarrow \downarrow \downarrow \downarrow }$$

$$c_{\nu }^{†}{c}_{\mu }^{†}{c}_{{\mu }^{\prime }}{c}_{{\nu }^{\prime }}\approx ⟨c_{\nu }^{†}{c}_{\mu }^{†}⟩c_{{\mu }^{\prime }}{c}_{{\nu }^{\prime }}+c_{\nu }^{†}{c}_{\mu }^{†}⟨c_{{\mu }^{\prime }}{c}_{{\nu }^{\prime }}⟩-⟨c_{\nu }^{†}{c}_{\mu }^{†}⟩⟨c_{{\mu }^{\prime }}{c}_{{\nu }^{\prime }}⟩,$$

$$H_{\mathrm{MF}}$$

$${\Delta }_{\mathbf{k}\sigma {\sigma }^{\prime }}\equiv -\frac{1}{N}\sum _{{\mathbf{k}}^{\prime }\tau {\tau }^{\prime }}{V}_{\sigma {\sigma }^{\prime }\tau {\tau }^{\prime }}(\mathbf{k},{\mathbf{k}}^{\prime })⟨c_{-{\mathbf{k}}^{\prime }\tau }{c}_{{\mathbf{k}}^{\prime }{\tau }^{\prime }}⟩,$$

$$E_{\mathrm{MF}}\equiv \frac{-1}{2N}\sum _{\mathbf{k}{\mathbf{k}}^{\prime }}\sum _{\sigma {\sigma }^{\prime }\tau {\tau }^{\prime }}{V}_{\sigma {\sigma }^{\prime }\tau {\tau }^{\prime }}(\mathbf{k},{\mathbf{k}}^{\prime })⟨c_{\mathbf{k}\sigma }^{†}{c}_{-\mathbf{k}{\sigma }^{\prime }}^{†}⟩⟨c_{-{\mathbf{k}}^{\prime }\tau }{c}_{{\mathbf{k}}^{\prime }{\tau }^{\prime }}⟩.$$

$$H_{\mathrm{MF}}=\sum _{\mathbf{k}}({\psi }_{\mathbf{k}}^{†}{h}_{\mathbf{k}}{\psi }_{\mathbf{k}}+{\xi }_{\mathbf{k}})+E_{\mathrm{MF}},$$

$${\psi }_{\mathbf{k}}\equiv \left(\begin{array}{c}{c}_{\mathbf{k}\uparrow }\\ c_{\mathbf{k}\downarrow }\\ c_{-\mathbf{k}\uparrow }^{†}\\ c_{-\mathbf{k}\downarrow }^{†}\end{array}\right);h_{\mathbf{k}}=\frac{1}{2}\left(\begin{array}{cc}{\xi }_{\mathbf{k}}\mathbb{I}& {\Delta }_{\mathbf{k}}\\ {\Delta }_{\mathbf{k}}^{†}& -{\xi }_{\mathbf{k}}\mathbb{I}\end{array}\right).$$

$$u_{\mathbf{k}}$$

$$v_{\mathbf{k}}$$

$$E_{\mathbf{k}}$$

$${\Gamma }_{\mathbf{k}}=\left(\begin{array}{c}{\gamma }_{\mathbf{k}\uparrow }\\ {\gamma }_{\mathbf{k}\downarrow }\\ {\gamma }_{-\mathbf{k}\uparrow }^{†}\\ {\gamma }_{-\mathbf{k}\downarrow }^{†}\end{array}\right)$$

$$E_g=E_{\mathrm{MF}}+\sum _{\mathbf{k}}({\xi }_{\mathbf{k}}-E_{\mathbf{k}}).$$

$$c_{\mathbf{k}\sigma }$$

$$c_{-\mathbf{k}\sigma }$$

$${\Delta }_{\mathbf{k}\sigma {\sigma }^{\prime }}=-\frac{1}{N}\sum _{\mathbf{k}\tau {\tau }^{\prime }}{V}_{\sigma {\sigma }^{\prime }\tau {\tau }^{\prime }}(\mathbf{k},{\mathbf{k}}^{\prime }){\Delta }_{{\mathbf{k}}^{\prime }{\tau }^{\prime }\tau }{g}_{{\mathbf{k}}^{\prime }},$$

$${\Delta }_{\mathbf{k}}^{\mathrm{sing}}$$

$${\Delta }_{\mathbf{k}}^x$$

$${\Delta }_{\mathbf{k}}^y$$

$${\Delta }_{\mathbf{k}}^z$$

$${\Delta }_{\mathbf{k}}^{\mathrm{sing}}$$

$${\Delta }_{\mathbf{k}}^{\lambda }=-\frac{1}{N}\sum _{{\mathbf{k}}^{\prime }}2V(\sin k_x\sin k_x^{\prime }+\sin k_y\sin k_y^{\prime }){\Delta }_{{\mathbf{k}}^{\prime }}^{\lambda }{g}_{{\mathbf{k}}^{\prime }},$$

$${\Delta }_{\mathbf{k}}^{\mathrm{sing}}$$

$${\Delta }_{\mathbf{k}}^{\lambda }$$

$${\Delta }_0$$

$${\Delta }_s$$

$${\Delta }_d$$

$${\Delta }_x^{\lambda }$$

$${\Delta }_y^{\lambda }$$

$$u$$

Full text

I Introduction

As more and more unconventional superconductors are discovered, the question of the symmetry of the superconducting order parameter has moved to the forefront as one of the most immediate and important questions to answer about any new material, especially given its close connection to the (often unknown) pairing mechanism. Experimentally, this question is difficult to answer, as the arduous history of the cuprates provides testament for Annett et al. (1996). The resolution of this question has been aided, in part, by phase-sensitive tunnelling measurements, which have been particularly useful in uncovering the gap symmetry of the heavy-fermion compound UPt3 Strand et al. (2009), shown to have different symmetries of the order parameter over different temperature ranges Norman (2011). Similarly, recent observations in LaAlO3/SrTiO3 Stornaiuolo et al. (2017) and underdoped La2-xSrxCuO4 (LSCO) Razzoli et al. (2013) suggest that a second component of the gap function develops below $T_{c}$. The purpose of this paper is to point out that one should generically expect to find a rich phase diagram below $T_{c}$, particularly at strong coupling. In fact we illustrate that symmetry transitions occur as a function of temperature in one of the simplest and most studied models for superconductivity, the extended Hubbard model on a square lattice. More particularly, this occurs in the regime of attractive nearest-neighbour interactions ($V<0$).

In one dimension (1D), the phase diagram of the extended Hubbard model with $V<0$ has been well established. In the presence of on-site attraction ($U<0$), the system hosts singlet superconductivity, while above a critical $U$, this transitions to triplet superconductivity; spin-density-wave order (SDW) and phase separation are observed as well Kuroki et al. (1994); Zhi-Guo et al. (2002); Sengupta et al. (2002). A definitive phase diagram in two dimensions (2D) is more difficult to produce owing in part to the presence of a sign problem Santos (2003). On a bipartite lattice such as the square lattice, sign-free quantum Monte Carlo away from half-filling requires attractive on-site, and repulsive nearest-neighbour interactions Wei et al. (2016). As a result, the $V<0$ regime is largely unexplored. This is unfortunate because an attractive inter-site interaction can arise, for example, due to polarons, as in the iron-pnictides Sawatzky et al. (2009). Thus, there are only a few benchmarks we use to frame the phase diagram of this region within. The ($V=0,\;U<0$) region is dominated by $s$-wave superconductivity Singer et al. (1998), and at half filling, the $U>0$ line contains an SDW transition Wu and Tremblay (2014); Buividovich et al. (2018). Doping slightly away from half filling yields $d$-wave superconductivity according to dynamical mean field theory Jarrell et al. (2001); Civelli et al. (2008), which makes this a good candidate model for the high-$T_{c}$ cuprates.

If we are to take any lessons from the exact results in 1D, it should be that both singlet and triplet superconductivity are present in the extended Hubbard model. However, in 2D we have the additional angular freedom of the gap profile in momentum-space. This makes the list of possible angular profiles much larger. In fact, as we show in this paper, one should not generically expect any rotation symmetries of the lattice to be retained by the gap function well below $T_{c}$.

Even when numerically exact methods like quantum Monte Carlo can be used, the question of the symmetry of the gap is far from answered. Such methods look for instabilities in specific angular channels and can easily miss mixed phases. Indeed as we will see, the energetic difference between mixed and non-mixed phases can be very small. In this paper we provide a systematic approach to find the exact profile of the gap function in momentum-space within the context of mean-field theory.

Constraints on the symmetry of the gap function are partially established within Landau-Ginzburg theory Annett (1990). In this paradigm, the gap function is segmented into pieces that transform under irreducible representations of the normal state symmetry group. Any of these individual pieces could form the superconducting state at $T_{c}$, but cannot be mixed at this temperature due to invariance of the free energy. However, at lower temperatures, when the magnitude of the order parameter is large, higher order terms in the Landau free energy become important. At these lower temperatures, the free energy becomes a complicated function of the gap, so invariance of the free energy under lattice symmetries does not impose such stringent conditions on the gap. As a result, mixing can occur. As we will see, the phases below $T_{c}$ are simply described by pitchfork bifurcations of critical points of the free energy. Such mixing and bifurcations have been predicted before in the context of anisotropic tight binding models Sörensen et al. (1991); Angilella et al. (1999). In this paper, we illustrate these ideas within a case study of the 2D extended Hubbard model. The zero-temperature superconducting phase diagram for this model is quite rich Nayak and Kumar (2018). We provide an efficient numerical method for determining the mean-field phase diagram, that is able to capture a phase that was missed with previous techniques. Our results indicate that many superconductors described by the extended Hubbard model are likely not strictly “s-wave” or “d-wave” but instead consist of gap functions that acquire multiple components below $T_{c}$.

II Model and Mean Field Solution

We consider the extended Hubbard model on a square lattice with unit lattice spacing, nearest neighbour hopping $t$, on-site interaction $U$, nearest-neighbour interaction $V$, and chemical potential $\mu$

[TABLE]

where $\langle i,j\rangle$ denotes nearest neighbour sites and $\sigma$ is a spin index. The standard pairing interaction, treated in all spin channels, results in a $2\times 2$ matrix $\Delta_{\alpha\beta}$, whose components satisfy a set of self-consistent gap equations (see appendix A for details). Since the purpose of this paper is to differentiate the competing superconducting symmetries, we focus on the attractive nearest-neighbour regime ($V<0$). In the half-filled positive $U$ regime, there is competing spin-density wave order with higher $T_{c}$ Micnas et al. (1990); Robaszkiewicz et al. (1981), but renormalization group results at zero temperature indicate that this phase does not survive much below $n=1$ Huang et al. (2013). Based on this evidence it seems unlikely that charge or spin-density waves play a significant role for $V<0,n<1$. We can therefore take our mean fields to include only superconducting order parameters $\Delta_{\alpha\beta}$.

First we consider the situation near $T_{c}$ where all gap components are small. The phase diagram at $T_{c}$ has been studied before Micnas et al. (1990). The gap equations decouple into three different sets of equations that correspond to the $s$-wave, $d$-wave and triplet $p$-wave phases. In general, each irreducible representation of the lattice point group will have a corresponding $T_{c}$ equation provided there is a part of the interaction that transforms under this representation. Unlike the $d$-wave case, the linearized $s$-wave gap equations can have two eigenvalues in some parts of the parameter space and therefore two critical temperatures. In Fig. 1, we show the phase diagram in the $U$-$V$ plane at quarter filling based solely on the critical temperatures determined by solving the linearized gap equations. Often, this diagram is taken to give the symmetry of the order parameter below $T_{c}$. However, we will see that the true phase diagram is much richer when we take into account the temperature dependence of the gap.

For the sake of generality, we may rewrite the gap function as

[TABLE]

Here the basis functions $x^{i}_{\bm{k}}$ are

[TABLE]

where $s_{\bm{k}}\equiv\frac{1}{2}(\cos k_{x}+\cos k_{y})$ and $d_{\bm{k}}\equiv\frac{1}{2}(\cos k_{x}-\cos k_{y})$. Such a decomposition can be made for any separable interaction and makes the expressions we derive valid for many models beyond the extended Hubbard model. Note that $x^{0}_{\bm{k}}$ and $x^{s}_{\bm{k}}$ belong to one irreducible representation of the square lattice point group ($A_{1}$), while $x^{d}_{\bm{k}}$ belongs to $B_{1}$ and ($x^{x}_{\bm{k}}$, $x^{y}_{\bm{k}}$) belong to another ($E$). For the singlet part, this means that the gap equations decouple whenever $\Delta_{d}=0$ or $\Delta_{0}=\Delta_{s}=0$. We refer to such solutions as pure solutions, and reserve “mixed solutions” for any case where components from multiple different irreducible representations are non-zero.

The natural question to ask is which symmetry of the gap would be observed for a given value of the interactions $U$, $V$, temperature $T$ and electron density $n$? The answer is whatever minimizes the free energy density. Minimization of the free energy is an important problem that requires strict numerical control, and we have developed a new method to deal it.

In appendix B and C, we show that the global minimum of the free energy proceeds through a series of bifurcations in parameter space as the temperature is lowered. Knowledge of these bifurcations provides an important numerical advantage. Identification of the symmetry of the order parameter at a given temperature requires the global minimization of a multidimensional function that is quite nonlinear. Without making use of these bifurcations, one would have to proceed by brute force search of a discretized parameter space in order to find all local minima Nayak and Kumar (2018). Instead, we leverage the bifurcations by doing a minimization within the lower-dimensional pure subspaces of the parameter space. The set of local minima found in these subspaces might not contain the global minimum, but one of them will lie within the basin of attraction of the global minimum. We can therefore initiate a local minimization scheme at each of these points and one of them will converge to the global minimum.

This procedure provides a significant reduction in computational time, scaling as the dimension of the largest pure subspace (2 in this model) rather than the dimension of the full paramter space (5 in this model). It may also be applied to any free energy minimization problem in which the critical points are continuous functions of some parameter (in this case temperature). For a single order parameter, this is equivalent to saying that the system only has a second-order phase transition. However, this does not preclude the existence of first-order transitions in the case of multiple competing components of the order parameter as we will see in Sec.III.

III Phase diagram

An example of a point in the phase diagram where a mixed solution takes over below $T_{c}$ is shown in Fig. 2. For these parameters, the system starts as a $s$-wave superconductor at $T_{c}^{s}$, but attains a lower free energy upon the emergence of $d$-wave components at a lower temperature forming an $s+id$ phase.

Fig. 3(a) shows a more complicated situation that is prevalent near the $T_{c}^{d}=T_{c}^{p}$ line. Here we see the development of a free energy minimum that mixes all symmetry sectors $s$-wave, $d$-wave and $p$-wave at a second-order transition below $T_{c}$. At an even lower temperature, this minimum gives way to the pure $p$-wave solution. Remarkably, both the ${\bm{k}}$-space symmetry and spin pairing change twice within the superconducting phase. Note that even the $p$-wave component is not continuous across this latter transition because the solutions on either side of the transition descend from different bifurcations of the free energy critical points. The $p$-wave component above the transition emerges from the $d$-wave phase, while the one below the transition emerges from the normal state in accordance with the bifurcation rules discussed in appendix C. This is more easily seen by comparing the free energies as in Fig. 3(b). Temperature profiles such as this would be accompanied by a jump in the specific heat at the onset of the mixed phase, as well as a second anomaly at lower temperature due to the reappearance of a pure phase. However, we note that the derivative of the free energy across the latter transition is remarkably smooth, as seen in Fig. 3(b), so the corresponding anomaly may be difficult to observe.

The zero-temperature phase diagrams are shown in Fig. 4 for different densities. Near half-filling, the phase diagram is dominated by a large mixed phase containing $s$, $p$ and $d$-wave components. Within this phase, the $s$ and $p$-wave components decrease as $V$ increases towards zero, first approaching a $d+ip_{x}$ phase, followed by a pure $d$-wave phase at weak coupling (not shown here). At lower fillings, a $p+ip$ phase is stabilized in agreement with the $T_{c}$ results, and an $s+id$ phase exists only near the border between $s$ and $s+d+ip_{x}$ phases, for all fillings. In Fig. 4 (c), we overlay the $T_{c}$ phase boundaries at quarter filling. This illustrates the dramatic distinction between the high and low temperature superconducting phases. The large region of the $d$-wave phase near $T_{c}$ is completely replaced at low temperatures by a phase with more broken symmetries. Even near half-filling, where the extended Hubbard model serves as a paradigmatic model for high-$T_{c}$ $d$-wave superconductors, the $d$-wave phase is mostly replaced by a mixed phase at low temperatures. Moreover, this rules out the possibility that such high-$T_{c}$ superconductors have a low temperature $s+id$ phase, since this $s+id$ phase occurs entirely to the left of the $T_{c}^{s}=T_{c}^{d}$ line.

Fig. 5 shows the angular dependence of the magnitude of the gap in the mixed phase. In this example, the $s$-wave component is very small so that the gap is effectively $d+ip_{y}$ (which is degenerate with $d+ip_{x}$). We see that the $p$-wave component opens a gap at the $d$-wave nodes and introduces a slight asymmetry between $\phi=0$ and $\phi=\pi$. In general, mixed gap functions can have zero, two or four nodes depending on the fermi level and the relative magnitude of the $\Delta_{i}$’s. This figure provides a prediction for photoemission experiments in some cuprate superconductors. While many cuprates display a node for all temperatures (e.g. Bi2Sr2CaCu2O8+δ Hashimoto et al. (2014)), this is not true in the deeply underdoped regime, particularly for underdoped LSCO Razzoli et al. (2013), which shows similar temperature dependence to Fig. 5. In this case, the authors of Razzoli et al. Razzoli et al. (2013) go beyond $T_{c}$ into the pseudogap phase. However, because mean field theory does not account for phase fluctuations, our results really probe the spectral gap, as discussed in Emery et al. Emery and Kivelson (1995). The nature of the pseudogap is of course widely debated, but we see that the changing symmetry of the energy gap in this region should not be construed as evidence for some competing order. On the other hand, our results show that one should generically expect changes to the gap symmetry at intermediate to strong coupling, arising solely from the pairing order parameter. This is consistent with the idea of preformed pairs, which is supported by recent experiments Zhou et al. (2019); Yuli et al. (2009). Fig. 6 shows that a transition from a $d$-wave to a $d+ip$ gap at low temperatures indeed fits well with the experimental data for LSCO.

With these examples, we see that the full temperature-dependent phase diagram is quite complicated. In general, there is a significant portion of the phase diagram where $d$-wave gives way to a mixed state below $T_{c}$ which may be followed by a $p+ip$ phase at even lower temperatures.

IV Conclusion

We have studied the temperature dependence of the symmetry of the superconducting gap in the 2D extended Hubbard model, focusing on the regime of attractive nearest neighbour interaction ($V<0$). We used a mean field approach to determine the gap equations that set the allowed symmetries of the gap.

The precise symmetry of the gap for given parameters ($U,V,n,T$) is determined by the minimum of the free energy, which can be considered as a critical point of a Morse function in a multi-dimensional parameter space. This perspective, outlined in the appendix B and C, has conceptual and technical advantages that may translate to other models as well.

We discovered that even within this simple model, there is a rich variety of symmetry phase transitions that occur as a function of temperature, observable through the specific heat or photoemission. In particular, we note that there are portions of the phase diagram where there is a different pure symmetry at $T=T_{c}$ and $T=0$ (e.g $d$-wave and $p+ip$), with multiple symmetry transitions at intermediate temperatures. It is interesting that even deep within the parameter range where $d$-wave has the highest $T_{c}$ (i.e. far from the $T_{c}$ phase boundaries), where we expect that different symmetries are least competitive, the mixed solution is dominant at low temperatures. At zero temperature where the magnitude of the gap is large, it appears that $d$-wave symmetry is never preferred at strong coupling, and the energy can always be lowered by the addition of $s$ or $p$-wave components. Provided these new phases are stable to quantum fluctuations, there are two possibilities; either the extended Hubbard model is not sufficient to capture the symmetry of $d$-wave superconductors with large $V$, or experiments to determine the gap symmetry should include fits to more complicated gap functions containing small $s$ or $p$-wave components. Indeed, we have compared the temperature-dependence found in the mixed phase with experiments in underdoped LSCO and found a striking similarity. It should be noted that the presence of p-wave in our result produces an asymmetry between the antinodes of the gap that was not tested for experimentally. Of course, it is still possible that the symmetry is given by $d_{x^{2}-y^{2}}+d_{xy}$, as suggested in Razzoli et al. Razzoli et al. (2013). More low-temperature data at both antinodes is needed to differentiate these two possibilities. $d_{xy}$ components do not appear in the extended Hubbard model, but may be included with e.g. a next-nearest neighbour hopping term. In that case one might find transitions from pure $d_{x^{2}-y^{2}}$ to mixed $d_{x^{2}-y^{2}}+d_{xy}$.

A word about the interaction parameter ranges in this paper is warranted. We focused mainly on strong nearest-neighbour coupling, and found that in this regime, the range of symmetry transitions is significant. Such transitions may also occur at weak coupling, but in this regime the free energies of the different phases are so similar as to make it difficult to discern the transition. Therefore it would be easiest to observe this phenomenon for large $|V|/t$. One promising avenue to explore this is in the context of ultra-cold atoms, where $U/t$ is controlled via the optical lattice potential, $n$ and $T$ are controlled via evaporative cooling, and $V$ can be included through dipolar interactions Tarruell and Sanchez-Palencia (2018); Menotti et al. (2008).

Lastly, we point out that the phase diagram is highly dependent on the symmetries of the normal state and therefore may be further enriched by the addition of Rashba or Zeeman terms.

Acknowledgements.

We wish to thank Elia Razzoli for sharing photoemission data with us. This work was supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC) as well as Alberta Innovates - Technology Futures (AITF). J.H. was supported by fellowships from NSERC, AITF, and the French National Research Agency (ANR).

Appendix A Mean Field Theory

After Fourier transforming the Hamiltonian we restrict our consideration to interactions between pairs of zero total momentum; the absence of finite-momentum pairing is expected in an unpolarized system. The resulting Hamiltonian reads

[TABLE]

Here $\xi_{\mathbf{k}}=-2t(\cos k_{x}+\cos k_{y})-\mu$, $N$ is the total number of sites in the lattice and the non-zero components of the interaction are

[TABLE]

where we have used the $s$-wave and $d$-wave basis functions: $s_{\bm{k}}\equiv\frac{1}{2}(\cos k_{x}+\cos k_{y})$, and $d_{\bm{k}}\equiv\frac{1}{2}(\cos k_{x}-\cos k_{y})$.

We apply mean-field theory to Eq. (LABEL:eq:Hk) following Ref. Sigrist (2005). Including all anomalous pairing mean fields 111We ignore Hartree-Fock terms here. These terms will renormalize the bandwidth and the chemical potential, but will not qualitatively change our results Micnas et al. (1990)., the generic four-fermion interaction (suppressing momentum labels for now) becomes

[TABLE]

where we have removed terms of order $(c_{\sigma}c_{\sigma^{\prime}}-\langle c_{\sigma}c_{\sigma^{\prime}}\rangle)^{2}$ and their hermitian conjugates.

The result is a Hamiltonian with only fermion bilinears:

[TABLE]

where we have defined the gap function

[TABLE]

and the energy due to the product of mean fields

[TABLE]

This may be compactly written as

[TABLE]

where

[TABLE]

Here $\mathbb{I}$ is the $2\times 2$ identity matrix, and $\Delta_{\bm{k}}$ is the $2\times 2$ matrix with components $\Delta_{{\bm{k}}^{\prime}\sigma\sigma^{\prime}}$. We will exclusively consider gap functions that are unitary: $\Delta^{\dagger}_{\bm{k}}\Delta_{\bm{k}}=|\Delta_{\bm{k}}|^{2}\mathbb{I}$, where $|\Delta_{\bm{k}}|^{2}\equiv\frac{1}{2}\mathop{\mathrm{Tr}}\Delta^{\dagger}_{\bm{k}}\Delta_{\bm{k}}$. A non-unitary gap allows for spin-polarized mean fields ($\Delta_{\uparrow\uparrow}$ and $\Delta_{\downarrow\downarrow}$) to be arbitrary, and are only warranted in proximity to a ferromagnetic instability Sigrist (2005). In particular, all superconductors that preserve time-reversal symmetry are unitary, as seen by the operation of time-reversal on the gap function $\Delta_{\bm{k}}\rightarrow\sigma_{y}\Delta^{*}_{\bm{k}}\sigma_{y}$. Moreover, some time-reversal symmetry-broken phases are unitary, including the $s+id$ phase.

$h_{{\bm{k}}}$ is diagonalized by the Bogoliubov transformation where

[TABLE]

Here we have used the fact that antisymmetry of $\Delta_{\bm{k}}$ requires that $v^{T}_{-{\bm{k}}}=-v_{\bm{k}}$ and $v^{*}_{-{\bm{k}}}=-v^{\dagger}_{\bm{k}}$ and we have defined the quasi-particle spinors

[TABLE]

via $\psi_{\bm{k}}=U_{\bm{k}}\Gamma_{\bm{k}}$. We have also identified the ground state energy

[TABLE]

In terms of these quasi-particle operators, we have

[TABLE]

This allows us to evaluate the mean-field expectation value $\langle c_{-{\bm{k}}\tau}c_{{\bm{k}}\tau}\rangle$ since the $\gamma_{\bm{k}}$ operators satisfy the fermion anticommutation relations and therefore have the occupation number distribution $\langle\gamma^{\dagger}_{{\bm{k}}\sigma}\gamma_{{\bm{k}}\sigma}\rangle=f(E_{\bm{k}})$, where $f(E)$ is the Fermi function. The defining gap equation (14) then becomes the self-consistency condition

[TABLE]

where we have defined $g_{{\bm{k}}}\equiv\frac{1}{2E_{\bm{k}}}(1-2f(E_{{\bm{k}}}))$.

It is conventional to separate the gap matrix into its singlet and triplet contributions:

[TABLE]

Since $\Delta_{-{\bm{k}}}=-\Delta_{\bm{k}}^{T}$, and $g_{{\bm{k}}}=g_{-{\bm{k}}}$, the singlet gap equation only contains contributions from the part of the interaction that is parity-even.

[TABLE]

Likewise, the triplet parts only contain contributions from the part of the interaction that is parity-odd. In fact, in this model all triplet components satisfy the same gap equation

[TABLE]

for $\lambda=x,y,z$.

The ${\bm{k}}$-dependence of Eqs. (A) and (32) determines the appropriate ansatz for the gap functions:

[TABLE]

and $|\Delta_{\bm{k}}|^{2}=|\Delta_{\bm{k}}^{\rm sing}|^{2}+\sum_{\lambda}|\Delta_{\bm{k}}^{\lambda}|^{2}$.

These components satisfy the following gap equations:

[TABLE]

Note that the triplet channels are degenerate since we may write Eqs. (38) and (39) as an eigenvector equation ${\bm{\Delta}}^{\lambda}=A{\bm{\Delta}}^{\lambda}$, where $A$ is a $2\times 2$ matrix independent of $\lambda$ and ${\bm{\Delta}}^{\lambda}=(\Delta^{\lambda}_{x},\Delta^{\lambda}_{x})^{T}$. In other words, for all three values of $\lambda$, ${\bm{\Delta}}^{\lambda}$ is an eigenvector of $A$ with eigenvalue $1$, but this can only be the case if all the ${\bm{\Delta}}^{\lambda}$ are equal to one another up to multiplicative constants which can be absorbed into the definition (38). Specifically, if ${\bm{\Delta}}^{x}/c_{1}={\bm{\Delta}}^{y}/c_{2}={\bm{\Delta}}^{z}/c_{3}$, then the rescaling ${\bm{\Delta}}^{\lambda}\rightarrow{\bm{\Delta}}^{\lambda}/\sqrt{|c_{1}|^{2}+|c_{2}|^{2}+|c_{3}|^{2}}$ allows us to write the spectrum (and the free energy) in terms of a single triplet component $\Delta^{\rm trip}_{\bm{k}}\equiv\Delta_{x}\sin k_{x}+\Delta_{y}\sin k_{y}$ so that $E_{{\bm{k}}}=\sqrt{\xi_{\bm{k}}^{2}+|\Delta^{\rm sing}_{\bm{k}}|^{2}+|\Delta^{\rm trip}_{\bm{k}}|^{2}}$.

Appendix B Free Energy

Combining the quasi-particle energy density with the ground state energy density yields the internal energy density:

[TABLE]

The entropy density is that of a free Fermi gas:

[TABLE]

Combining these produces the mean-field free energy density $f_{\rm MF}=u-Ts$:

[TABLE]

$f_{\rm MF}$ gives the correct value of the free energy at the solutions to the gap equation, i.e. its critical points, but it is not the correct functional to minimize in order to obtain these critical points. That functional comes from the finite temperature variational theorem Feynman (1955)

[TABLE]

where the mean-field Hamiltonian $H_{\rm MF}$ is given in Eq. (13). Of course, if we know the solutions to the gap equation, we can simply plug them into $f_{\rm MF}$ and compare the resulting free energies, but it will prove fruitful to work with the variational free energy $f[\Delta_{\bm{k}}]$. Here, the mean-field expectation of an operator $X$ is given by

[TABLE]

We must compute

[TABLE]

The trace of the quartic term is readily evaluated in the eigenbasis of $H_{\rm MF}$. Using Wick’s theorem and the fact that the occupation number distribution of the quasi-particles $\langle\gamma^{\dagger}_{{\bm{k}}\alpha}\gamma_{{\bm{k}}\alpha}\rangle$ is the Fermi function $f(E_{\bm{k}})$, we obtain the identity

[TABLE]

Using this identity and applying the Bogoliubov transformation (23), (24), we can evaluate the quartic term

[TABLE]

The free energy functional follows from this,

[TABLE]

where we have separated the interaction into its parity-even and odd parts:

[TABLE]

Minimization of the free energy is an important problem that requires strict numerical control. Ensuring that one has obtained the global minimum of this function is numerically difficult in a brute-force approach, which amounts to extensively searching for all local minima in the five-dimensional parameter space spanned by the $\Delta_{i}$ components. In this regard, we have made use of the theory of Morse functions Palais (1963); Morse (1934, 1934) to aid our numerical search for the global minimum, which allows us to effectively reduce the dimensionality of the parameter space one needs to explore from five to two in this model.

Appendix C Critical Points

The gap equations for the extended Hubbard model (35)-(39) can be reproduced from the solutions to the functional minimization problem:

[TABLE]

while the chemical potential is fixed by the number equation

[TABLE]

which upon using (51) becomes,

[TABLE]

The functional minimization of $f$ can be thought of as a minimization with respect to the $N$ terms of $\Delta_{\bm{k}}$ (one for each value of ${\bm{k}}$). However, we can alternatively think of this as a minimization over a lower-dimensional parameter space spanned by ($\Delta_{0}$, $\Delta_{s}$, $\Delta_{d}$, $\Delta_{x}^{\lambda}$, $\Delta_{y}^{\lambda}$), so that we need only solve

[TABLE]

for the all the components $\Delta_{i}\in\{\Delta_{0},\Delta_{s},\Delta_{d}\}$ and $\Delta^{\lambda}_{i}\in\{\Delta_{x}^{\lambda},\Delta_{y}^{\lambda}\}$. The first derivative is given by

[TABLE]

where we have defined $G_{\bm{k}}\equiv g_{\bm{k}}+\frac{|\Delta_{\bm{k}}|^{2}}{E_{\bm{k}}}\frac{\partial g_{\bm{k}}}{\partial E_{\bm{k}}}$ and used the basis functions (3)-(7). Since $G_{\bm{k}}$ is positive definite, these equations simply restate what we already know; the critical points of $f$ occur at the values of $\Delta_{i}$ that satisfy the gap equation. One conceptual advantage of this formulation is that we can now compute a simple Hessian in a low-dimensional parameter space. For now, let us focus on real singlet order parameters to better understand these critical points. In that case, the Hessian elements evaluated at the critical points are

[TABLE]

where $G_{\bm{k}}$ is evaluated at the solutions to the gap equation.

In this case the critical points have two symmetries apparent from the gap equation. First, there is inversion symmetry; if $(\Delta_{0},\Delta_{s},\Delta_{d})$ is a critical point, then so is $(-\Delta_{0},-\Delta_{s},-\Delta_{d})$. This is required by antisymmetry of the pair wavefunction. Second, there is mirror symmetry under reflection in the $\Delta_{0}$-$\Delta_{s}$ plane, since any solution $(\Delta_{0},\Delta_{s},\Delta_{d})$ has a corresponding solution $(\Delta_{0},\Delta_{s},-\Delta_{d})$ upon replacing $k_{x}\leftrightarrow k_{y}$, $k_{x}^{\prime}\leftrightarrow k_{y}^{\prime}$. Moreover, the curvature (57) is a rank-2 tensor under these transformations. That is, under inversion,

[TABLE]

and under reflection in the $\Delta_{0}$-$\Delta_{s}$ plane it follows the Householder transformation Householder (1958)

[TABLE]

This means that the gradient flow of $f$ near the critical points respects these symmetries, and the index of any critical points related by these symmetries is the same. The index $\gamma_{i}$ of a critical point $i$ is defined as the number of orthogonal directions along which $i$ is a maxima of the corresponding function (in this case the free energy). The index plays a large role in the theory of Morse functions Palais (1963). A Morse function is a smooth real function that has no degenerate critical points. For our purposes, the free energy is a Morse function, since the curvature can only vanish on a set of measure zero in the $\{T,U,V,n\}$ parameter space. As a consequence, the free energy satisfies the following Morse condition (valid for functions $f(\{\Delta_{i}\})$ that increase without bound as each $|\Delta_{i}|$ goes to infinity):

[TABLE]

where $\chi(M)$ is the Euler characteristic of the domain manifold of $f$ Morse (1934). We may choose $M$ to be the three-dimensional space $\{\Delta_{0},\Delta_{s},\Delta_{d}\}$, or any of the pure subsets $\{\Delta_{0},\Delta_{s}\}$ or $\{\Delta_{d}\}$, because any extrema of $f$ on a pure subset is guaranteed to be a solution to the full gap equation. All of these cases correspond to flat manifolds with $\chi(M)=1$. The combination of symmetries of $f$ and the Morse condition restrict the possible minima of the free energy in this parameter space.

C.1 Corollaries of the Morse condition

Continuing with the case of a real singlet gap function, the simplest question to ask is what can happen at $T_{c}$ according to these restrictions. Above $T_{c}$, the system is in the normal state, there is one critical point, and on any of the above manifolds, Eq. (60) reads

[TABLE]

We immediately see that no mixed solution can emerge from this state. This is because mixed solutions are four-fold degenerate by inversion and mirror symmetry, so that the Morse condition reads

[TABLE]

denoting the normal index by $\gamma_{N}$ and the mixed index by $\gamma_{m}$. This equation has no integer solutions. Thus the condition expressed in Ref. Annett (1990) that mixed solutions cannot emerge at $T_{c}$ is an immediate consequence of the symmetry of the gap equation.

We may also prove the uniqueness (modulo sign) of the pure $d$-wave solution. In fact, the following will hold for any single-parameter gap function, including the Hubbard model with only on-site attraction. Existing proofs of the uniqueness of the BCS gap solution are quite nontrivial Vansevenant (1985). Here we will see that it is a simple consequence of symmetry and the Morse condition.

Proof.

For a single component gap $\Delta_{\bm{k}}=\Delta x_{{\bm{k}}}$, with a separable interaction $V_{{\bm{kk}}^{\prime}}=Vx_{\bm{k}}x_{{\bm{k}}^{\prime}}$, the Hessian reads

[TABLE]

The single-component gap equation reads

[TABLE]

so that at a critical point

[TABLE]

$\frac{\partial g_{{\bm{k}}^{\prime}}}{\partial E_{{\bm{k}}^{\prime}}}$ is negative definite, while $G_{{\bm{k}}}$ is positive definite, so for any attractive interaction $V<0$, the curvature is positive and all solutions must be minima of the free energy. Since the solutions are symmetric under inversion $\Delta\rightarrow-\Delta$, they must come in pairs. For $n$ such pairs, the Morse condition reads

[TABLE]

whose only solutions are $\{\gamma_{N}=0,n=0\}$ and $\{\gamma_{N}=1,n=1\}$ 222Since this is a one-dimensional problem, this last statement is just Rolle’s theorem.. Thus, the single-component BCS gap equation admits one solution modulo sign. ∎

Returning to the extended Hubbard model, this result applies to solutions that are pure $d$-wave (in this case, $x_{\bm{k}}\rightarrow d_{\bm{k}}$ and $V\rightarrow 4V$ in the arguments above), guaranteeing their uniqueness as well as the instability of the normal state. At a critical $d$-wave temperature $T_{c}^{d}$, two $d$-wave critical points emerge from the normal state critical point as the temperature is lowered. In fact, this kind of bifurcation of critical points is generic.

The solutions to the gap equation are continuous functions of $T$, so critical points cannot appear in pairs at arbitrary points in the parameter space but must grow from the normal state or an existing superconducting state as the temperature is lowered. The Morse condition then provides a conservation of indices. For example, a pure state with index $\gamma_{p}$ that grows from the normal state at a temperature $T_{c}^{p}$ is two-fold degenerate (barring pathological accidental degeneracies). This degeneracy comes from the inversion symmetry of the zeros of Eqs. (C), (56) in the $\Delta_{i}$ parameter space. The corresponding critical points therefore must satisfy

[TABLE]

where $\gamma_{N_{1}}$ is the index of the normal state at $T>T_{c}^{p}$, and $\gamma_{N_{2}}$ is the index of the normal state at $T<T_{c}^{p}$. The solution is

[TABLE]

So the normal state bifurcates into two pure states passing its index to the new solutions. Mixed states can grow out of pure states via the same bifurcation mechanism. An example of the evolution of critical points as the temperature is lowered is shown in Fig. 7.