Fate of a strongly correlated $d$-wave superconductor in a Zeeman field: The Fulde-Ferrel-Larkin-Ovchinnikov perspective

The Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase is an unconventional superconducting state found under the influence of a strong Zeeman field. This phase is identified by finite center-of-mass momenta in the Cooper pairs, causing the pairing amplitude to oscillate in real space. Repulsive correlations, on the other hand, smear out spatial inhomogeneities in $d$-wave superconductors. We investigate the FFLO state in a strongly correlated $d$-wave superconductor within a consolidated framework of Hartree-Fock-Bogoliubov theory and the Gutzwiller approximation. We find that the crucial effects of strong correlations lie in shifting the BCS-FFLO phase boundary towards a lower Zeeman field and thereby enlarging the window of the FFLO phase. In the FFLO state, our calculation features a sharp midgap peak in the density of states, indicating the formation of strongly localized Andreev bound states. We also find that the signatures of the FFLO phase survive even in the presence of an additional translational-symmetry-breaking competing order in the ground state. This is demonstrated by considering a broken-symmetry ground state with a simultaneous presence of the $d$-wave superconducting order and a spin-density wave order, often found in unconventional superconductors.

Figure 1

Energetics of ground states with $h$ for different broken symmetries, calculated within RIMT (a) and IMT (b) methods (see Sec. 2 for model parameters). A $d$-wave BCS state (blue curve), an FFLO state (red curve) and a normal state (green curve) are considered for comparison. Though the FFLO ground state is found sandwiched between the $d$-wave BCS state (low $h$) and the normal state (large $h$) both in RIMT and IMT, it is realized in a wider range of $0.18\le h\le 0.38$ in RIMT findings (a) than that from IMT (b) window of $0.22\le h\le 0.38$. We mark the FFLO and BCS regions by pink and blue shades, respectively. The inset in each panel shows that the nature of modulations of the superconducting pairing amplitude within the FFLO phase changes from a stripe to a checkerboard form (see Sec. 3a for details). Such a changeover occurs at $h\approx 0.25$ in the RIMT calculation (a) and at $h\approx 0.33$ in the IMT scheme (b).

Figure 2

Determination of $q^{*}$, the optimal modulation wave vector of the pairing amplitude, versus $h$, from RIMT (a) and IMT (c) calculations. The $h$ dependence of $q^{*}$ from RIMT and IMT findings is shown in panels (b) and (d), respectively. For small $h, E_t\left(q\right)$ achieves the minimum for $q=0$, indicating a uniform BCS $d$-wave ground state. On the other hand, $E_t\left(q\right)$ features a minimum at finite $q^{*}$ for $0.18\le h\le 0.38$ in RIMT calculations in panel (a), and for $0.22\le h\le 0.38$ in plain IMT calculations in panel (c), signaling an FFLO ground state. For large $h>h_2$, the energy minimum is lost and the $d$-wave pairing amplitude becomes feeble, indicating the collapse of FFLO phase. The evolution of $q^{*}$ with $h$ is strongly sublinear, with apparent saturation at large $h$ within the RIMT scheme (b), whereas its dependence is approximately linear in IMT calculations (d). The insets in the panels (b) and (d) focus on the energy landscapes with respect to $q$ for specific $h=0.25$ in panel (b) and $h=0.26$ in panel (d) within the FFLO region in RIMT and IMT, respectively.

Figure 3

Evolution of different order parameters as a function of $h$, obtained from (a) RIMT and (b) IMT calculations. The primary orders shown here are root-mean-square pairing amplitude ${\mathrm{\Delta }}_{\mathrm{rms}}$ (blue curves), average magnetization $m_0$ (red curves), and the intertwined SDW $m_{2q}$ (magenta curves) and CDW $n_{2q}$ (green curves) orders (defined in Sec. 3b), all presented at $q=q^{*}\left(h\right)$. In $d$-wave BCS regime ${\mathrm{\Delta }}_{\mathrm{rms}}$ remains nearly constant and $m_0$ is small. The onset of FFLO regime at $h_1$ [$h_1=0.18$ and 0.25 in panels (a) and (b)] is signaled by a sharp fall of ${\mathrm{\Delta }}_{\mathrm{rms}}$ and a steep rise in $m$. These two orders keep decreasing and increasing in the FFLO regime, with a much slower rate in RIMT results than in IMT. The exit from FFLO regime to NS at $h_2$ ($h_2=0.38$ in both panels) is signaled by a near vanishing of ${\mathrm{\Delta }}_{\mathrm{rms}}$ while $m_0$ reaches its normal-state value. The self-generated orders $m_{2q}$ and $n_{2q}$ survive only in the FFLO regime and are stronger in RIMT results than in IMT findings.

Figure 4

Evolution of the spatial profiles of order parameters [see Eqs. (4) and (6) and the text following Eq. (7) for definitions] in the FFLO phase from IMT calculations. The spatial profiles are featured in a $3\times 3$ panel on the left side using color-density plots. The evolution of $\mathrm{\Delta }$ is shown in the left column with $h$ increasing from top to bottom. Similarly, magnetization $m$ and accompanying CDW order are shown in the middle and right columns, respectively. The wavelength of the stripe modulation decreases with increasing $h$ from 0.26 to 0.32, and finally leads to a checkerboard pattern at a higher field ($h=0.34$). The line plots on the right side provide a cross-sectional view of these profiles along the $\widehat{x}$ direction for stripe ($h=0.26,0.32$) and along the $\widehat{x}+\widehat{y}$ direction for checkerboard ($h=0.34$) modulation.

Figure 5

Evolution of spatial profiles of order parameters [see Eqs. (4) and (6) and the text following Eq. (7) for definitions] in the FFLO phase, similar to Fig. 4 but from RIMT calculations, is shown to emphasize the role of strong repulsive correlations here. The line plots on the right provide cross-sectional view of the order parameter profiles along the $\widehat{x}$ direction for stripe ($h=0.18,h=0.26$) and along the $\widehat{x}+\widehat{y}$ direction for checkerboard ($h=0.35$) modulation. The spatial profiles here are obtained from our real-space calculations in a $40\times 40$ lattice system. Note the sharp rise and fall of orders across the line of zero of $\mathrm{\Delta }$ (discussed in Sec. 3d) and the robustness of the wavelength of modulation of the order parameters with increasing $h$, consistent with our momentum-space findings in Fig. 2.

Figure 6

Behavior of the effective hopping parameter $t_{\mathrm{eff}}$ [subtracting $t_{\mathrm{eff}}(h=0)$] versus $h$ from RIMT (a) and IMT (b) findings. The $h$ dependence of the effective magnetic field $h_{\mathrm{eff}}$ from RIMT (c) and IMT (d) outcomes. The strong correlations in panel (a) cause up- and down-spin hopping (denoted by red and blue curves, respectively) to branch out weakly in the BCS regime, while their difference rises sharply in the FFLO phase. The difference increases strongly again upon exiting FFLO phase into NS. In contrast, the up- and down-spin hoppings maintain nearly identical value across the entire range of $h$ within the IMT method (the weak difference is only due to Fock shifts). The behavior of $h_{\mathrm{eff}}$ follows a linear trend with bare $h$ in both RIMT (c) and IMT (d) findings for all $h$. The RIMT calculation results in a somewhat lower $h_{\mathrm{eff}}$ in the BCS region and significantly lower $h_{\mathrm{eff}}$ in the FFLO and NS region compared to IMT results. $h_{\mathrm{eff}}$ faces jumps at $h_1$ in both RIMT and IMT calculations—the jump being more significant in RIMT. $h_{\mathrm{eff}}$ also undergoes a weak jump near $h_2$ in RIMT, unlike in IMT, where only a continuous change in $h_{\mathrm{eff}}$ across $h_2$ is found.

Figure 7

Spin-resolved average DOS from RIMT (left column) and IMT methods (right column). Panels (a) and (b) show DOS for $d$-wave BCS state, featuring standard profile where up- and down-spin DOSs are oppositely shifted in energy from Fermi level by $h_{\mathrm{eff}}$. Panels (c) and (d) present DOS in the FFLO phase close to corresponding $h_1$ highlighting the signature of a bound state through the midgap peak. Notice that the presence of strong correlations makes the midgap peak much sharper in panel (c) compared to IMT outcome in panel (d). Further increase of $h$ deep inside FFLO regime in panels (e) and (f) begins to broaden and subsequently split the midgap peak. Finally, for $h>h_2$ pairing amplitude collapses altogether as seen in panels (g) and (h), and the resulting DOS features standard profile of tight-binding electrons in the presence of a magnetic field in the normal state.

Figure 8

Spatial distribution of the low-lying wave function $|{\varphi }_n\left({\mathbf{r}}_i\right){|}^2\equiv |u_{i,n\uparrow }{|}^2+{\left|v_{i,n\uparrow }\right|}^2$ for $E_n\approx -h_{\mathrm{eff}}$ to illustrate its boundedness. Panel (a) demonstrates that $|{\varphi }_n\left({\mathbf{r}}_i\right){|}^2$ at $h=0.18$ is a more tightly bound state in RIMT calculations than in plain IMT result for $|{\varphi }_n\left({\mathbf{r}}_i\right){|}^2$ in panel (b) for $h=0.26$. Note that the lack of tight boundedness in the IMT wave function leads to a relatively broader midgap peak in DOS, e.g., in Fig. 7.

Figure 9

Energetics of the proposed competing phases with different broken symmetries, calculated within RIMT (a) and IMT (b) methods. We consider the $d$-wave BCS state (blue dashed curves), the FFLO phase (red dashed curves), the normal state (green dashed curves), a state with both the $d$-wave BCS order and a competing SDW order with $(\pi ,\pi )$ modulation (dark blue curves), a state with FFLO modulation $+$ competing SDW order with $(\pi ,\pi )$ modulation (pink curves), and a normal state with $(\pi ,\pi )$ SDW order (dark green curves) as the possible candidates for comparison, from which the true GS emerges as the one with minimum energy at the given $h$. The state with FFLO modulation with SDW order (pink shade) is realized as the GS between $d$-wave $\mathrm{BCS}+\mathrm{SDW}$ state for low $h$ and the normal state at large $h$ with both RIMT and IMT schemes. This phase is realized for a range of $h_1\approx 0.23$ to $h_2\approx 0.45$ in RIMT findings (a), and from IMT (b) window of $h_1\approx 0.2$ to $h_2\approx 0.43$.

Figure 10

Evolution of the spatial profiles of different order parameters [see Eqs. (4) and (6) and the text following Eq. (7) for definitions] with $h$, with a ground state that features a competing SDW order along with the $d$-wave superconducting order and at $h=0$, for the chosen set of model parameters, as mentioned in Sec. 3e. The $2\times 3$ panels with color-density plots show the spatial profiles of superconducting pairing amplitude (left panels), magnetization (middle panels), and charge density (right panels) on the left side for $h=0.2$ (top panels) and $h=0.36$ (bottom panels) obtained from plain IMT calculations. The $2\times 3$ right-side panels are similar depictions from RIMT calculations for $h=0.25$ and $h=0.35$, respectively.

Figure 11

Spatial profiles of the magnitude of the pairing amplitude with $d$-wave [panels (a) and (b)] and extended $s$-wave [panels (c) and (d)] symmetries, shown in the FFLO phase for $h=0.18$ (left panels) and $h=0.35$ (right panels) within RIMT calculations. In panels (e) and (f) corresponding average DOSs for up-spin electrons are shown, which are obtained by fully suppressing the $d$-wave pairing and considering only the self-consistent extended $s$-wave pairing. The extended $s$-wave pairing has very small amplitudes in regions where the $d$-wave pairing in suppressed, and elsewhere it is negligible. In fact, the resulting DOS for both values of $h$ resembles that of the underlying NS of a nearest-neighbor tight-binding model, confirming the irrelevance of the extended $s$-wave pairing in the FFLO phase.

Figure 12

Differences of the components of $E_t$ near the phase boundaries ($h_1$ and $h_2$) as a function of $h$ in RIMT (a) and IMT (b) methods. The shaded regions in pink represent the FFLO state. Here, $E_t$ gain of the FFLO equals $E_t$ of the competing phase, such as BCS or NS, minus $E_t$ of the FFLO phase, and $E_t$ loss of the FFLO equals $E_t$ of the FFLO phase minus $E_t$ of the competing phase, such as BCS or NS. $E_t$ gain of the FFLO (red curve near $h_1$) phase through $E_K$ and $E_m$ crosses the $E_t$ loss of the FFLO (blue curve near $h_1$) due to $E_p$ with respect to the BCS state at $h_1$, which occurs early in RIMT. At $h_2$ the difference of $E_t$ components between FFLO and NS approaches zero. Consequently, $E_t$ gain and $E_t$ loss of FFLO with respect to NS merge with each other at the FFLO-NS boundary.

Figure 13

Differences of the components of $E_t$ near the phase boundaries as a function of $h$ in RIMT (a) and IMT (b) methods. The shaded regions in pink represent $\mathrm{FFLO}+\mathrm{SDW}$ state. Here, $E_t$ gain of the $\mathrm{FFLO}+\mathrm{SDW}$ equals $E_t$ of the competing phase, such as $d$-wave $\mathrm{BCS}+\mathrm{SDW}$ or NS-SDW, minus $E_t$ of the $\mathrm{FFLO}+\mathrm{SDW}$ phase, and $E_t$ loss of the $\mathrm{FFLO}+\mathrm{SDW}$ equals $E_t$ of the $\mathrm{FFLO}+\mathrm{SDW}$ phase minus $E_t$ of the competing phase, such as $d$-wave $\mathrm{BCS}+\mathrm{SDW}$ or NS-SDW. $E_t$ gain of the $\mathrm{FFLO}+\mathrm{SDW}$ (black curve near $h_1$) through $E_m$ crosses the $E_t$ loss of the $\mathrm{FFLO}+\mathrm{SDW}$ (magenta curve) due to $E_p$ and $E_K$ at $h_1$, which occurs at 0.23 and 0.2 in RIMT and IMT, respectively. At $h_2$ the difference of $E_t$ components between $\mathrm{FFLO}+\mathrm{SDW}$ and NS-SDW approaches zero. Consequently, $E_t$ gain and $E_t$ loss of $\mathrm{FFLO}+\mathrm{SDW}$ with respect to NS merge with each other at the $h_2$ boundary, which is $\approx 0.43$ and $\approx 0.45$ in RIMT and IMT.