Fulde-Ferrell-Larkin-Ovchinnikov state in strongly correlated d-wave superconductors

TL;DR

This paper explores the FFLO state in strongly correlated d-wave superconductors, revealing how strong correlations influence phase boundaries, induce localized states, and persist despite competing orders.

Contribution

It introduces a comprehensive framework combining Hartree-Fock-Bogoliubov theory and Gutzwiller approximation to study FFLO states in correlated d-wave superconductors, highlighting new effects of correlations.

Findings

Strong correlations shift the BCS-FFLO boundary to lower Zeeman fields.

The FFLO state exhibits a sharp mid-gap peak indicating localized Andreev bound states.

FFLO signatures persist even with competing spin-density wave order.

Abstract

The Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase is an unconventional superconducting state found under the influence of 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 Gutzwiller approximation. We find that the profound 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 mid-gap 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.

Equations

$${\mathcal{H}}_{\mathrm{Hub}}=-t\sum _{⟨ij⟩\sigma }\left({\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+h.c\right)+U\sum _i{\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow }.$$

$${\mathcal{H}}_{\mathrm{t}-\mathrm{J}}=-t\sum _{⟨ij⟩\sigma }\left({\tilde{c}}_{i\sigma }^{†}{\tilde{c}}_{j\sigma }+h.c\right)+J\sum _{⟨ij⟩}\left({\tilde{\mathbf{S}}}_{\mathbf{i}}.{\tilde{\mathbf{S}}}_{\mathbf{j}}-\frac{{\tilde{n}}_i{\tilde{n}}_j}{4}\right).$$

$${\mathcal{H}}_{\mathrm{GA}}$$

$${\Delta }_{ij\sigma }=\sum _{\sigma }⟨{\psi }_0|{\hat{c}}_{i\sigma }{\hat{c}}_{j\overline{\sigma }}|{\psi }_0⟩,$$

$${\tau }_{ij\sigma }=⟨{\psi }_0|{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }|{\psi }_0⟩,$$

$$n_{i\sigma }=⟨{\psi }_0|{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{i\sigma }|{\psi }_0⟩;m_i=\frac{1}{2}\sum _{\sigma }\sigma n_{i\sigma }.$$

$${\mathcal{H}}_{\mathrm{MF}}$$

$$n_{i\sigma }=n_0^{\sigma }+2n_{2q}^{\sigma }\cos (2\mathbf{q}.{\mathbf{r}}_i)+2n_{4q}^{\sigma }\cos (4\mathbf{q}.{\mathbf{r}}_i)+..$$

$$g_{i,i+\delta }^{t\sigma }=g_0^{t\sigma }+2g_{2\mathbf{q}\delta }^{t\sigma }\cos (2\mathbf{q}.r_i)+2g_{4\mathbf{q}\delta }^{t\sigma }\cos (4\mathbf{q}.r_i)+..$$

$$g_{i,i+\delta }^{J,\nu }=g_0^J+2g_{2q\delta }^J\cos (2\mathbf{q}.{\mathbf{r}}_i)+2g_{4q\delta }^J\cos (4\mathbf{q}.{\mathbf{r}}_i)+..$$

$${\mathcal{H}}_{\mathrm{MF}}$$

$${\mathcal{H}}_{\mathrm{MF}}={\Psi }^{†}{\hat{\mathrm{H}}}_{\mathrm{MF}}\Psi .$$

$${\hat{\mathrm{H}}}_{\mathrm{MF}}=\left[\begin{array}{cc}{\hat{\xi }}_{kp\sigma }& {\hat{\Delta }}_{kp}\\ {\hat{\Delta }}_{kp}^{\ast }& -{\hat{\xi }}_{kp\overline{\sigma }}\end{array}\right],$$

$${\hat{\xi }}_{kp\sigma }={\delta }_{kp}{\xi }_{p\sigma }^{(r)}+{\delta }_{k\pm 2q,p}{\xi }_{p\sigma }^{(r)},$$

$${\hat{\Delta }}_{kp}={\delta }_{-k\pm q,p}{\Delta }_{kp}^{(r)},$$

$${\Psi }^{†}=\left[\begin{array}{ccccccc}{\hat{c}}_{k_1\sigma }^{†},& .& ,& {\hat{c}}_{k_N\sigma }^{†},{\hat{c}}_{k_1\overline{\sigma }},& .& ,& {\hat{c}}_{k_N\overline{\sigma }}\end{array}\right].$$

$$\sum _p\left[\begin{array}{cc}{\hat{\xi }}_{kp\sigma }& {\hat{\Delta }}_{kp}\\ {\hat{\Delta }}_{kp}^{\ast }& -{\hat{\xi }}_{kp\overline{\sigma }}^{\ast }\end{array}\right]\left[\begin{array}{c}{u}_{p,n\sigma }\\ v_{p,n\overline{\sigma }}\end{array}\right]=E_{n\sigma }\left[\begin{array}{c}{u}_{k,n\sigma }\\ v_{k,n\overline{\sigma }}\end{array}\right],$$

$$\sum _j\left[\begin{array}{cc}{\hat{\xi }}_{ij\sigma }& {\hat{\Delta }}_{ij}\\ {\hat{\Delta }}_{ij}^{\ast }& -{\hat{\xi }}_{ij\overline{\sigma }}^{\ast }\end{array}\right]\left[\begin{array}{c}{u}_{j,n\sigma }\\ v_{j,n\overline{\sigma }}\end{array}\right]=E_{n\sigma }\left[\begin{array}{c}{u}_{i,n\sigma }\\ v_{i,n\overline{\sigma }}\end{array}\right],$$

$${\hat{\xi }}_{ij\sigma }$$

$${\hat{\Delta }}_{i,j}=\left\{J{g}_{i,j}^{J,xy}{\Delta }_{ij\downarrow }-\frac{J}{2}\left(g_{i,j}^{J,z}+1\right){\Delta }_{ij\uparrow }\right\}{\delta }_{i+\delta ,j},$$

$${\Delta }_{Q^{\prime }}=\frac{1}{N}\sum _k\left[⟨c_{-k+Q^{\prime }\downarrow }{c}_{k\uparrow }{⟩}_0{\eta }_k+⟨c_{-k+Q^{\prime }\downarrow }{c}_{k\uparrow }{⟩}_0{\eta }_{-k+Q^{\prime }}\right],$$

$$m_Q=\sum _{k,\sigma }\frac{\sigma n_Q^{\sigma }}{2};n_Q=\sum _{k,\sigma }{n}_Q^{\sigma },$$

$${\mathcal{H}}_{\mathrm{MF}}$$

$$N_{\sigma }(\omega )$$

$$g_{ij}^{t\sigma }=\sqrt{g_i^{t\sigma }{g}_j^{t\sigma }}$$

$$g_i^{t\sigma }=\sqrt{\frac{2{\delta }_i\left(1-{\delta }_i\right)}{\left(1-{\delta }_i^2+4m_i^2\right)}\frac{1+{\delta }_i+\sigma 2m_i}{1+{\delta }_i-\sigma 2m_i}}$$

$$g_{ij}^{J,xy}=g_i^{J,xy}{g}_j^{J,xy}$$

$$g_{ij}^{J,xy}=\frac{2{\delta }_i\left(1-{\delta }_i\right)}{\left(1-{\delta }_i^2+4m_i^2\right)}$$

$$g_{ij}^{J,z}=g_{ij}^{J,xy}\frac{2\left({\Delta }_{ij}^2+{\tau }_{ij}^2\right)-4m_i{m}_j{X}_{ij}^2}{2\left({\Delta }_{ij}^2+{\tau }_{ij}^2\right)-4m_i{m}_j}$$

$$X_{ij}^2=1+\frac{12\left(1-{\delta }_i\right)\left(1-{\delta }_j\right)\left({\Delta }_{ij}^2+{\tau }_{ij}^2\right)}{\sqrt{\left(1-{\delta }_i^2+4m_i^2\right)\left(1-{\delta }_j^2+4m_j^2\right)}}$$

Full text

I Introduction

A magnetic field destroys superconductivity in two ways. One is through orbital effect, which couples the magnetic field to the orbital motion of the electrons. This creates vortices in a superconductor by puncturing holes in the superconducting pairing amplitude, through which the magnetic flux lines penetrate. With increasing magnetic field, the density of the vortices increases and the pairing amplitude fails to recover between them, causing superconductivity to collapse progressivelyTinkham (2004). The second one is the Zeeman effect, where the magnetic field couples to the spin degrees of freedom of the electrons. This strains the spin-singlet configuration of the Cooper pairs owing to the split Fermi surfaces of spin-up and spin-down electrons. As the Zeeman field $h$, increases, a superconductor to normal state (NS) transition occurs at the Clogston-ChandrashekharClogston (1962); Chandrasekhar (1962) limit, where the magnetization energy due to the Fermi surface splitting overcomes the condensation energy of the Cooper pairs. However, it was later shown that the Cooper pairs can survive beyond the Clogston-Chandrashekhar limit by having a finite pairing momentum, which can make the superconducting pairing amplitude spatially modulating, as proposed by Fulde and FerrellFulde and Ferrell (1964) and by Larkin and Ovchinnikov independentlylarkin and Ovchinnikov (1964). The advantage of having spatially modulating pairing amplitude in this Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state is to accommodate the local magnetization in regions where the pairing amplitude vanishes. Hence, both the superconducting pairing amplitude and the magnetization survive in the same system by periodically avoiding each other and the superconductivity is stabilized even at sufficiently large $h$. As a result, the system undergoes a transition from a usual BCS superconducting state at small $h$ to a FFLO state at an intermediate range of $h$ and finally to a NS at very high $h$.

Even though the existence of the FFLO state had long been a theoretical truismZhou and Ting (2009); Yang and Sondhi (1998); Wang et al. (2006); Jin et al. (2006); Cui and Yang (2008); Yanase and Sigrist (2009a); Cui et al. (2012); Shimahara (1994); Samokhvalov et al. (2010); Maśka et al. (2010); Loh et al. (2011); Yanase (2009),it evaded experiments for many years. Two primary reasons behind this are the presence of disorder and the dominance of orbital effect of the applied magnetic field. A modulating FFLO phase is extremely sensitive to disorder Aslamazov (1969); Takada (1970)(see however, Ref. Loh and Trivedi, 2010) and orbital effects, if dominant, often destroys superconductivity even before the appearance of the FFLO phase. The relative importance of the Zeeman effect and the orbital effect is characterized by the Maki parameterSaint-James et al. (1969) $\rm{\alpha}=\sqrt{2}\rm{H}^{\rm{orb}}_{c_{2}}(T=0)\left[\rm{H}_{\rm{p}}(T=0)\right]^{-1}$, where $\rm{H}^{\rm{orb}}_{c_{2}}$ is the orbital critical field and $\rm{H}_{p}$ is the Pauli limiting field. The orbital-limited materials, characterized by $\alpha<1$, are thus unfavorable to conceive the FFLO phase. Despite of these facts, recently this exotic phase has gained renewed interests, due to the experimental indications of its presence in some unconventional superconductors. This includes heavy fermion superconductors like $\mathrm{CeCoIn_{5}}$, $\mathrm{CeCu_{2}Si_{2}}$, some organic superconductorsMoll et al. (2011); Wright et al. (2011); Mayaffre et al. (2014); Uji et al. (2012); Agosta et al. (2012); Cho et al. (2009); Lortz et al. (2007); Beyer et al. (2012) and even some of the Fe-based superconductorsCho et al. (2011); Zocco et al. (2013). Extensive experimental evidences have been obtained in $\mathrm{CeCoIn_{5}}$ in favor of the existence of this inhomogeneous superconducting state (commonly coined as the Q phase in these materials) in the presence of magnetic field Tayama et al. (2002); Radovan et al. (2003); Gratens et al. (2012); Martin et al. (2005); Bianchi et al. (2003a); Correa et al. (2007); Spehling et al. (2009); Kumagai et al. (2006); Koutroulakis et al. (2010); Capan et al. (2004); Kenzelmann et al. (2008, 2010). This material is a Pauli-limitedMatsuda and Shimahara (2007) ($\alpha\approx 5$), quasi-2D, d-wave superconductorVorontsov and Vekhter (2006); Zhou et al. (2013) and can be obtained in a reasonably clean form. More recently, nuclear magnetic resonance experiments on $\mathrm{CeCu_{2}Si_{2}}$ Kitagawa et al. (2018) and organic superconductor $\rm{\kappa-(BEDT-TTF)_{2}Cu(NCS)_{2}}$ Mayaffre et al. (2014) revealed the presence of Andreev bound states when magnetic fields were applied parallel to their conduction planes, one characteristic of the FFLO phase.

Strong electronic interactions play a crucial role in most of these unconventional superconductors. $\rm{CeCoIn_{5}}$ is believed to be strongly correlated for having electrons with high effective massesPetrovic et al. (2001), for its proximity to antiferromagnetic instabilityBianchi et al. (2003b); Hu et al. (2012) and for spin-dependent quasiparticle mass enhancement observed in De Haas-van Alphen oscillation measurementsMcCollam et al. (2005). Other candidates like $\mathrm{CeCu_{2}Si_{2}}$Takenaka et al. (2017); Pang et al. (2018) and $\rm{\kappa-(BEDT-TTF)_{2}Cu(NCS)_{2}}$H. McKenzie (1998) are also prototypes of strongly correlated superconductors. The dominance of strong correlations among these materials leads to the natural question: what is the role strong electronic correlations in the FFLO phase? Earlier theoretical studies on d-wave superconductors suggest that strong interactions change the nature and the degree of inhomogeneities and smear out the small scale spatial charge fluctuationsGarg et al. (2008); Chakraborty and Ghosal (2014). Strong correlations are thus expected to modify the existence or the nature of the FFLO phase which exhibit periodic inhomogeneities in charge and spin densities. Moreover, the interplay of strong correlations with the magnetic field is also expected to uncover interesting physics. The effect of spin-dependent mass enhancement on the FF phase, a homogeneous counterpart of the FFLO phase, has already been studied recentlySpałek (2011).

In this work, we investigate the role of strong correlation in the FFLO state in a d-wave superconductor within an integrated framework of Hartree-Fock-Bogoliubov theory and Gutzwiller approximation. The crux of our findings are as follows: (i) At $T=0$, near the optimal doping for superconductivity, strong correlations renormalize the different energy scales of the system. As a result, their subtle balance shifts the boundaries of the FFLO phase and consequently, it increases the FFLO window of the Zeeman field. (ii) The behaviors of the order parameters and the pairing momenta in the presence and absence of strong correlations are contrasting in nature owing to the renormalizations of different parameters in the Hamiltonian. (iii) Strong correlations cause a sharper and narrower mid-gap peak appearing at the density of states in the FFLO phase. This is due to a strong localization of the Andreev bound states at regions having zero superconducting pairing amplitude owing to its sharper fall near those regions. (iv) The signature of the FFLO phase survives even when the $h=0$ ground state (GS) has competing orders. We show this by considering a GS that has a competing spin-density wave (SDW) order in addition to the d-wave BCS order.

The rest of the paper is organized as follows: In Sec. II we give the details of the model used in our calculations emphasizing the way the effect strong correlation is introduced through Gutzwiller factors. We also discuss the computational method of our study at $T=0$. In Sec. III.1-III.4 we present our results where we compare the phase diagrams, the behavior of the parameters and the observables with respect to the applied magnetic field in the presence and absence of strong correlations. In Sec. III.5 we discuss the phase diagrams in the presence of an additional competing SDW order. Finally, we conclude in Sec. IV.

II Model and methods

We describe our system by the microscopic Hubbard Hamiltonian:

[TABLE]

Here, $t$ is the hopping energy of the electrons to its nearest neighbor, denoted as $\langle ij\rangle$, on a 2D square lattice, and $U$ is the on-site repulsion energy between the electrons. In the strongly correlated limit $U\gg t$, an effective low energy Hamiltonian can be obtained from ${\cal H}_{\rm{Hub}}$ which lives in a restricted Hilbert space that prohibits double occupancy of any site due to strong on-site repulsions. The resulting Hamiltonian is known as the $t-J$ model, which can be considered as the perturbative expansion of ${\cal H}_{\rm{Hub}}$ up to the quadratic order in ${t}/{U}$:

[TABLE]

Here, the exchange interaction $J=4t^{2}/U$ emerges via Schrieffer-Wolff transformation. The renormalized creation operator ${\tilde{c}}^{\dagger}_{i\sigma}=\left(1-\hat{n}_{i\overline{\sigma}}\right)\hat{c}^{\dagger}_{i\sigma}$ is defined to operate on the Hilbert space that excludes all double occupancies. In order to probe the FFLO state, we introduce the Zeeman field, $h$ and redefine ${\cal H_{\rm{t-J}}}$ to ${\cal H_{\rm{t-J}}}-\sum_{i\sigma}\sigma h\hat{n}_{i\sigma}$. Analyzing Hamiltonian in Eq. (2) is challenging due to the projection operators which are usually dealt using Variational Monte Carlo methods. A simpler implementation of the projections is achieved through a Gutzwiller approximation (GA) Zhang et al. (1988), an approximation method which incorporates the Hilbert space restrictions in a spirit similar to that of a mean field theory. Within the framework of the GA, we write the GS wave-function as, $|\psi\rangle=\Pi_{i}(1-\hat{n}_{i\uparrow}\hat{n}_{i\downarrow})|\psi_{0}\rangle$, where $|\psi_{0}\rangle$ is the GS wave-function in the unrestricted Hilbert space. In GA, the effects of projection are mimicked in Gutzwiller renormalization factors (GRFs)Ko et al. (2007) represented as $g^{t,\sigma}_{ij}$, $g^{J,z}_{ij}$, $g^{J,xy}_{ij}$, which depend on the local densities, magnetization, pairing amplitude and kinetic energies. The explicit expressions of the Gutzwiller factors are provided in App. A.1. Finally, using GA, the Hamiltonian for strongly correlated d-wave superconductors in a Zeeman field is shaped as,

[TABLE]

We define our local order parameters as follows:

[TABLE]

[TABLE]

[TABLE]

Using the above order parameters, the mean-field decomposition of ${\cal H}_{\rm{GA}}$ in Hartree, Fock and Bogoliubov channels, whose details are provided in App. A.2, leads to,

[TABLE]

where, ${\delta}$ represents nearest-neighbor spacing from the $i^{\rm{th}}$ site and $g^{t\sigma}$, $g^{xy}$, $g^{z}$, $\mu_{\sigma}$,$\tau_{\sigma}$ and $\tau_{\overline{\sigma}}$ are abbreviated notations corresponding to $g^{t\sigma}_{i,i+\delta}$, $g^{J,xy}_{i,i+\delta}$, $g^{J,z}_{i,i+\delta}$, $\mu+\sigma h$, $\tau_{i,i+\delta\sigma}$, $\tau_{i,i+\delta\overline{\sigma}}$ respectively. We have used the notation, $\phi_{i\sigma}=\partial W/\partial n_{i\sigma}$, with $W=\langle\psi_{0}|{\cal H}_{\rm{GA}}|\psi_{0}\rangle-\lambda\left(\langle\psi_{0}|\psi_{0}\rangle-1\right)-\mu\left(\sum_{i}n_{i}-\langle n\rangle\right)$. Here, $\lambda$ is a Lagrange multiplier fixing the wave-function renormalization $\langle\psi_{0}|\psi_{0}\rangle$=1 and $\mu$ is the chemical potential which takes care of the average density $\langle n\rangle=N^{-1}\sum_{i}n_{i}$ of the system. In this work, we focus on the d-wave symmetry of the superconducting pairing amplitude, defined as, $\Delta_{i}=\sum_{\sigma}(\Delta^{+\hat{x}}_{i,\sigma}+\Delta^{-\hat{x}}_{i,\sigma}-\Delta^{+\hat{y}}_{i,\sigma}-\Delta^{-\hat{y}}_{i,\sigma})/4$. To gauge the effects of strong correlations, we compare two sets of results obtained in the presence and absence of strong correlations. In the presence of strong correlations, we rely on the framework based on Gutzwiller approximation augmented with Hartree-Fock-Bogoliubov theory, as discussed already, whereas, in the absence of strong correlations, we employ Hartree-Fock-Bogoliubov theory on the unrestricted Hilbert space, i.e, without any double occupancy prohibition. Operationally, this is equivalent to setting the Gutzwiller factors to unity. We will use the notation IMT (Inhomogeneous Mean-field Theory) to refer to the calculation based on Hartree-Fock-Bogoliubov theory and RIMT (Renormalized Inhomogeneous Mean-field Theory) to refer to the scheme based on Hartree-Fock-Bogoliubov theory and Gutzwiller approximation in the discussions from now on. We set $U=12t$, for RIMTNorman et al. (1995) and choose $U=3.077t$ and set the Gutzwiller factors to unity for IMT. This yields the same d-wave superconducting gap at $h=0$ in these two schemes. We place $\langle n\rangle$ at $0.84$, a value which ensures a pristine superconducting state away from the dominance of the competing orders. For example, such $\langle n\rangle$ value is known as the optimal doping for cuprate superconductors. This is convenient for our case, since the primary motivation here is to study the effect of strong correlation in the FFLO phase i.e, a superconducting phase with finite momentum Cooper-pairs. However, it is also interesting to study the effect in the presence of competing orders which we will discuss later in Sec. III.5. In our investigation, we will focus on the LO state where the pairing profile in the lattice is $\Delta_{i}\approx 2\Delta_{{q}}\cos(\mathbf{{q}}.\mathbf{{r}}_{i})$. Here, $\mathbf{{q}}$ is the pairing momentum of the Cooper pairs and $\mathbf{{r}}_{i}$ denotes the position of the $i^{th}$ lattice site. Such a behavior of the pairing amplitude arises due to the coupling of the single particle states $|\mathbf{{k}},\sigma\rangle$ with both $|-\mathbf{{k}}+\mathbf{{q}},\overline{\sigma}\rangle$ and $|-\mathbf{{k}}-\mathbf{{q}},\overline{\sigma}\rangle$. The states $|-\mathbf{{k}}+\mathbf{{q}},\overline{\sigma}\rangle$ and $|-\mathbf{{k}}-\mathbf{{q}},\overline{\sigma}\rangle$ also connect to $|\mathbf{{k}}+2\mathbf{{q}},\sigma\rangle$ and $|\mathbf{{k}}-2\mathbf{{q}},\sigma\rangle$ respectively in the Cooper channel. This links the single particle state $|\mathbf{{k}},\sigma\rangle$ with the states $|\mathbf{{k}}\pm 2\mathbf{{q}},\sigma\rangle$, $|\mathbf{{k}}\pm 4\mathbf{{q}},\sigma\rangle$ and so on with progressive weaker coupling. Subsequently, an intertwined spin-density wave (SDW) and a charge-density wave (CDW) order is generated with modulating wave-vectors $2\mathbf{{q}}$, $4\mathbf{{q}}$, etc, which seeds in many intriguing consequences in the presence of strong electronic correlations. Note that, allowing only FFFulde and Ferrell (1964) pairing does not generate any such SDW or CDW order. We finally solve the mean-field Hamiltonian in the momentum space, owing to the periodic inhomogeneity in the FFLO phase. Drawing from the above discussions, the spatial profile of spin-densities look like:

[TABLE]

where $n_{\rm{Q}}^{\sigma}=N^{-1}\sum_{k}\langle\hat{c}^{\dagger}_{k+Q\sigma}\hat{c}_{k\sigma}\rangle_{0}$; where $Q=0$, $\pm 2q$, $\pm 4q$ and so on. Here, $\langle..\rangle_{0}$ signify expectation value with respect to the unprojected wave-function $|\psi_{0}\rangle$. As the GRFs are functions of local spin-densities, we fuse the following ansatz for them:

[TABLE]

[TABLE]

where $\nu=`xy\textrm{'}$ or $`z\textrm{'}$. Using the periodic translational symmetry of this phase, we solve Eq. (7) in the momentum space. We take the spin-densities up to $n^{4{q}}_{\sigma}$ mode, pairing amplitude up to $\Delta_{3q}$ (which along with the higher order modes also arise naturally with a weaker coupling as a result of the connected chains of single particle states) and neglect the higher order modes to simplify our calculations. We have checked our results by considering the higher order modes in the spin-densities and the pairing amplitude in our Gutzwiller mean-field theory calculations, and our results indicate that their effects are minimal in determining the phase boundaries and the important features of the physical observables we have studied here.

Building on these ideas, ${\cal H}_{\rm MF}$ in the momentum space becomes,

[TABLE]

The explicit expressions of $\xi^{(\rm r)}_{k\sigma}$, $\xi^{(\rm r)}_{k\pm 2q\sigma}$ and $\Delta^{(\rm r)}_{k,-k\pm q}$ are given in App. A.3. The Hamiltonian in Eq. (11) can be written in Nambu space, as

[TABLE]

Here, $\hat{\rm{H}}_{\rm{MF}}$ is a $2\rm{N}\times 2\rm{N}$ matrix, which has the following form

[TABLE]

where $\hat{\xi}_{kp\sigma}$ ,$\hat{\Delta}_{kp}$ are expressed as follows

[TABLE]

[TABLE]

[TABLE]

Here, $\rm{N}=\rm{L}\times\rm{L}$, indicating the total number of sites in the lattice with $\rm{L}$ being the length of the same. We exploit the translational symmetry of the system and block diagonalize $\hat{\rm{H}}_{\rm{MF}}$ in Eq. (13) into smaller matrices. A typical size of the block is $2\rm{L}\times 2\rm{L}$, which can be further reduced depending on the periodicity of the order parameters in the lattice. Most of our calculations are for $\rm{L}=200$, except the results shown in real space are for $L=40$. We diagonalize the resulting Hamiltonian using the transformations, $\hat{c}_{k\sigma}=\sum_{n}\left(u_{k,n\sigma}\gamma_{n\sigma}-\sigma v^{\ast}_{k,n\sigma}\gamma^{\dagger}_{n\overline{\sigma}}\right)$ , where $\gamma_{n\sigma}$ and $\gamma^{\dagger}_{n\overline{\sigma}}$ are Bogoliubov quasiparticle operators and $u_{kn\sigma}$ and $v_{kn\sigma}$ satisfy the equation:

[TABLE]

The Hamiltonian in Eq. (11) in the presence of Zeeman field is expected to describe phases like d-wave BCS, FFLO, polarized NS, as well as CDW or SDW order and interplay among these in the GS.

As mentioned earlier, we will also present some of the results in real space in Sec. III.1, III.2 and III.5, for which we use Bogoliubov de-Gennes (BdG) transformations G. de Gennes , $c_{i\sigma}=\sum_{n}(\gamma_{n,\sigma}u_{i,n}-\sigma\gamma^{\dagger}_{n\overline{\sigma}}v^{\ast}_{i,n})$, for diagonalizing Eq. (7). $\gamma^{\dagger}_{n\sigma}$ and $\gamma_{n,\sigma}$ are the creation and annihilation operators of the Bogoliubov quasiparticles. This results in the following eigen-equation,

[TABLE]

which is self-consistently solved for all the local order parameters defined in Eq. (4), (5) and 6). The matrix equation in Eq. (18) leads to the following equations,

[TABLE]

[TABLE]

and similarly for $v_{i,n\sigma}$. Diagonalizing the BdG matrix in Eq. (18) is numerically expensive. So, we solve the BdG matrix to obtain the real space pictures only for $L=40$. Considering smaller systems comes with an energy cost for the FFLO state, as it puts constraint on the possible $q$ values for the variational determination of the lowest energy state.

III Results

In the following we discuss the fate of FFLO phase due to strong correlations at $T=0$ by contrasting the behavior of different observables obtained within RIMT and IMT calculations.

III.1 Phase diagram

Earlier theoretical studies on superconductors subjected to a Zeeman field, $h$, suggest that a BCS superconductor with spatially uniform pairing amplitude undergoes a phase transition to an FFLO phase with modulated order parameter at $h=h_{1}$. Upon increasing $h$ further, superconductivity gets fully suppressed for $h\geq h_{2}$ (here, $h_{2}>h_{1}$) leading to a spin polarized normal state (NS) Yang and Sondhi (1998); Zhou and Ting (2009); Cui and Yang (2008); Loh and Trivedi (2010). The FFLO phase is thus sandwiched between an uniform BCS (for low $h\leq h_{1}$), and a spin-polarized NS (at large $h\geq h_{2}$). These three phases are identified by blue, pink and white shades in Fig. 1 (and also later in Fig. 2, 3, 6, 11, 12). These phases are realized in our calculations as the GS for a given $h$ within the framework of a Hartree-Fock-Bogoliubov description of ${\cal H}_{\rm t-J}$, in both RIMT and IMT calculations. The location of the phase boundaries, i.e. the values of $h_{1}$ and $h_{2}$ differ from the two methods of calculations, as shown in Fig. 1. The differences arise because of the physics of strong electronic correlations captured by the RIMT method, as we proceed to discuss below.

We identify the GS by considering the possible broken symmetry solutions, and then choosing the one with the lowest energy as shown in the main panel of Fig. 1. We find that the introduction of the Gutzwiller factors lowers $h_{1}$ (while keeping $h_{2}$ more or less unaltered) and thus enhances the window of stability of the FFLO GS. In particular, for our model parameters i.e., $J/t=0.33$ and average density $\langle n\rangle=0.84$, we obtain $h_{1}\approx 0.18$ and $h_{2}\approx 0.38$ (expressed in the units of $t$) by RIMT method in Fig. 1(a). In contrast, the IMT calculation, based on weak-coupling description, results into $h_{1}\approx 0.24$ and $h_{2}\approx 0.38$, as shown in Fig. 1(b) for comparable model parameters (as discussed in Sec. II). This is a narrower window of $h$ than what is found from RIMT. The possibility of correlation induced enhancement of FFLO phase in the parameter space is exciting in the context of strongly correlated systems.

Our investigation on the origin of this enhancement suggests that a subtle balance among relevant energy scales plays the most crucial role. In particular, the intricate interplay between components of the total mean field energy $E_{t}$ ($=\langle\psi_{0}|{\cal H_{\rm{GA}}}|\psi_{0}\rangle$) corresponding to the the pairing amplitude (the pairing energy $E_{p}$), the spin-imbalance (the magnetic energy $E_{m}$) and the hopping of electrons (the kinetic energy $E_{K}$) decide the window of the FFLO phase. All of these components are renormalized by the GRFs (also see the discussions in the next paragraph). We note that, within the FFLO regime the nature of spatial modulation of the superconducting pairing amplitude and magnetization changes from a stripe modulation at lower range of $h$, to a checkerboard pattern at higher $h$ side of FFLO phase, as highlighted in the insets of Fig. 1. Here, stripe implies a unidirectional modulation with $\mathbf{q}=\pm q(1,0)$ whereas, the checkerboard modulation with wave vector $\mathbf{q}$ is identified with an equal superposition of modulations with $\pm q(1,0)$ and $\pm q(0,1)$. This transformation of the nature of modulation is consistent with the earlier calculations Zhou and Ting (2009) on a square lattice. We find this transition from stripe to checkerboard modulation to occur at $h=0.25$ in RIMT findings (inset of Fig. 1), and at $h=0.33$ in IMT method (inset of Fig. 1(b)). We also mention here that, though the main panels of Fig. 1 present results obtained on a much larger system of size $200\times 200$, the results in the insets highlighting the transition in the modulation pattern are solved on a system size of $40\times 40$. This results in a weaker energy resolution compared to the traces presented in the main panel. A smaller system size for the checkerboard pattern of modulation is needed due to its reduced translational symmetry. In fact, we will only consider the stripe modulation of the FFLO phase to illustrate several of the subsequent analysis, specifically, in Sec. III.2-III.4, unless otherwise mentioned. This is because, we obtain an enhanced clarity of the results, characterizing the FFLO phase, in larger systems.

In order to develop a deeper understanding of the change of phase boundaries between RIMT and IMT results, we note that both $E_{m}$, and $E_{K}$ of the total energy $E_{t}$ lead to an energy gain in the FFLO phase as $h$ is increased. In contrast the spatially modulated pairing amplitude results in an energy loss, when compared to a homogeneous BCS state. It is this fine balance between these gain and loss of energy which dictates the boundaries between different phases as $h$ increases. On the other hand, the prohibition of double occupancy through Gutzwiller factors in RIMT renormalizes the separate components of energy differently, and hence it is natural that the aforementioned balance will occur for different $h$ values in RIMT and IMT calculations. We have elaborated this aspect quantitatively in App. B to provide a comprehensive picture.

As discussed already, the FFLO phase is characterized by a spatially modulating pairing amplitude and magnetization. It is crucial to decide the correct modulation wave vector. Other than possibilities for the different natures of modulations (consistent with the square symmetry of the underlying lattice), e.g. stripe and checkerboard patterns as discussed already, we also determine the optimal magnitude of $q^{\ast}$ by variationally minimizing the total energy $E_{t}$ of the FFLO state over the entire range of $q$. This is illustrated in Figs. 2(a) and (b) where $E_{t}(q)$ traces are presented (considering stripe pattern of modulation) using RIMT and IMT calculations respectively. Representative $h$ values in Fig. 2 are chosen from each of BCS, FFLO and spin-polarized NS. We find for very small $h\ll h_{1}$ that $E_{t}(q)$ has minimum value at $q=0$. For $h\lesssim h_{1}$, a second minimum in $E_{t}(q)$ emerges at a finite $q^{\ast}$, though $E_{t}(q^{\ast})\gtrsim E_{t}(q=0)$. In the FFLO phase for $h_{1}\leq h\leq h_{2}$, however, $E_{t}(q)$ develops a global minimum at $q^{\ast}$.

The value of $q^{\ast}$ is expected to increase Yang and Sondhi (1998) with $h$. A larger $q^{\ast}$ ensures a larger number of nodes in the spatial profile of the pairing amplitude causing more domain walls. These domain walls supports the magnetization arising from $h$. Hence, with increasing $h$, FFLO state gains more energy by increasing domain wall density. As a result, $q^{\ast}$ increases with $h$. Our results, however, establishes that the nature of this rise is different in RIMT and IMT calculations. While $q^{\ast}$ follows an apparent sub-linear increase in RIMT method as in Fig. 2(c), it rises approximately linearly in IMT results as seen in Fig. 2(d). In RIMT, $q^{\ast}$ increases rapidly for $h\gtrsim h_{1}$ and then gets saturated in a large part of the FFLO window. This is because of an increased role of the effective repulsion between domain walls at higher $h$. In RIMT, $\Delta$ undergoes a rapid change of magnitude and sign across narrow domain walls, particularly at small $h$ (for reasons discussed in Sec. III.2). Hence, the effective repulsion between the domain walls has little role when their density is small near $h_{1}$. However, at large $h$, the increased effective repulsion between these domain walls does not allow their density to rise as much for large $h$, leading to a near-saturation of $q^{\ast}$. In contrast, the qualitative behavior of $q^{\ast}$ in IMT agrees well with the previous studies Yang and Sondhi (1998). Here, the profile of the pairing amplitude is sinusoidal in the entire FFLO window, and the role of effective repulsion remains weak in the entire FFLO regime which results into an approximately linear $q^{\ast}(h)$.

In addition to the variational determination of $q^{\ast}$ outlined above, the energetics of the FFLO, the BCS (i.e. $q=0$) and the underlying spin-polarized NS are compared to identify the true GS, for each value of $h$. The $q$-resolution of our calculation is enhanced as we have exploited the translational symmetry of the FFLO phase across the lattice as mentioned earlier in Sec. II and solve the eigen-system of Eq. (17) for ${\cal H}_{\rm MF}$ in the momentum space for a large system. This yields good precisions for the individual components of energy.

III.2 Order parameters

The phase diagram in Fig. 1 identifies the boundaries between distinct phases, whereas, the energy minimum at $q^{\ast}$ decides the stability of the broken symmetry GS. The different energy scales in each of these states and their behavior with $h$, is characterized by the $h$-dependence of various order parameters characterizing our system. With this motivation, we study in Fig. 3, the behavior of the root mean square pairing amplitude ($\Delta_{\rm{rms}}$), average magnetization ($m_{0}$), and the self-generated intertwined SDW ($m_{2q}$) and CDW ($n_{2q}$) orders as a function of $h$. The $\Delta_{\rm{rms}}$ takes a value $\Delta_{0}$ in the homogeneous BCS state and $\sqrt{2}\Delta_{q}$ in the FFLO state. Here,

[TABLE]

[TABLE]

where $Q^{\prime}=0$ or $q$, $Q=0$ or $2q$ and $\eta_{k}=2({\rm cos}k_{x}-{\rm cos}k_{y})$ is the d-wave form factor. The behavior of different order parameters is contrasted from the two calculations: RIMT and IMT. For all $h$, the average magnetization $m_{0}$ attains a higher value in RIMT, as seen in Fig. 3(a). This is due to the reduced bandwidth upon prohibition of double occupancy. The superconducting order parameter is also found stronger in RIMT, see Fig. 3. This is because of our choice of exchange coupling $J$ in the Hamiltonian in Eq. (2) in the two methods of calculations to obtain the same value of superconducting energy gap at $h=0$ in RIMT and IMT methods. The energy gap and the pairing amplitude differ in RIMT calculations as they obtain different Gutzwiller renormalization. The $m_{0}$ and $\Delta_{\rm{rms}}$ calculated within RIMT scheme experience little change with $h$, inside the FFLO region (up to $h\approx 0.36$, a value close to $h_{2}$), as depicted in Fig. 3(a). These order parameters, reach their NS values with further increase in $h$. Such behavior of the order parameters is related to the saturation of $q^{\ast}$ for a wide window of $h$, followed by a quick change of $q^{\ast}$ near $h_{2}$ within the FFLO region (see Fig. 2). The order parameter values largely depend on $q^{\ast}$. For a given $q^{\ast}$, increasing $h$ causes only little changes in the order parameters. In IMT, these orders change continuously across the FFLO regime, finally attaining their NS values beyond $h_{2}$. The magnitude of the coexisting SDW ($m_{2q}$) and CDW ($n_{2q}$) order, which are self-generated due to the modulated pairing amplitude, though have small values, also increases with strong correlations.

Having understood the behavior of the global order parameters, we next focus on their spatial profiles for different values of $h$. The self-consistent spatial structure of modulating pairing amplitude, magnetization and charge density, are depicted in color-density plots on the left side in Fig. 4 and Fig. 5 from IMT and RIMT schemes respectively, whereas, the cuts on the right side emphasize the one-dimensional modulations. The panels from top to bottom present results for increasing $h$. These results are obtained using real-space BdG simulations carried out on a system of size $40\times 40$. We find from Fig. 4 that the wavelength of the stripe modulations of all three order parameters from IMT calculations decrease as $h$ is increased from $h=0.26$ to $h=0.32$. The corresponding wavelength is inversely proportional to the pairing momentum of the Cooper-pairs. In the RIMT scenario, however, the wavelength of stripe modulation changes marginally by going from $h=0.18$ to $h=0.26$, which is roughly consistent with the weak dependence of $q^{*}$ on $h$ within FFLO regime in the RIMT scheme, as shown in Fig. 2, also discussed earlier in Sec. III.1. The magnetization $m$ nucleates near the location of nodes of the superconducting pairing amplitude $\Delta$, thus the modulating wavelength of magnetization becomes half of the superconducting pairing amplitude. The modulation in local density $n$ also has this same wavelength.

In the RIMT scenario, $\Delta$ goes through a sharp fall where it changes sign. This is because strong electronic repulsions act to suppress the nanoscale density fluctuations locally causing a relatively smooth variations in the spatial density which in turn flattens the small-scale variations in the superconducting pairing amplitude self-consistently in the lattice Chakraborty and Ghosal (2014); Garg et al. (2008). As a result, the magnitudes of the higher order modes in the superconducting pairing amplitude and the density modulation increase in a self-consistent manner. These factors make a steeper spatial variation of the superconducting pairing amplitude $\Delta$ where it changes the sign in the presence of strong correlations as depicted in Fig. 5(a1,a4) compared to that from the IMT outcomes as shown in Fig. 4(a1,a4). This leaves a narrower space for local $m$ to nucleate compared to what is found from IMT method. The sharpness of the fall of $\Delta$ in RIMT, however, reduces with increasing $h$, reflecting a reduction in the dominance of strong correlations in higher fields.

III.3 Renormalized parameters

The distinctive features associated with the global and local properties of the order parameters at different values of $h$ are dictated by the renormalized parameters in the Hamiltonian. The mean field Hamiltonian in Eq. (7) can be re-cast in terms of the renormalized parameters using the following form,

[TABLE]

The explicit expressions of $t^{\sigma}_{\rm eff}(i)$ and $\mu^{\sigma}_{\rm eff}(i)$ can be obtained by comparing Eq. (7) with Eq. (23). The major contribution in $t^{\sigma}_{\rm eff}$ in RIMT, comes from the GRF for hopping, $g^{t\sigma}$ ($=0.275$ at $\langle n\rangle=0.84$, when $h=0$), which restricts the hopping solely to the unoccupied sites. Strong correlations also induce spin-dependence in the hopping parameters at finite $h$. At finite $h$, the number of sites occupied by the down-spin species decreases in the lattice. Thus the up-spin electrons find it easier to hop around and vice versa (assuming up-spin is favored by $h$), also reflected in the expressions of $g^{t\sigma}$ in Eq. (26). We show in Fig. 6 the evolution of renormalized hopping parameter $t^{\sigma}_{\rm eff}=\sum_{i,\delta}t^{\sigma}_{\rm eff}(i,\delta)$ as a function of $h$. The reduced $t_{\rm eff}$ ($\approx 0.4$ at $h=0$), thereby the reduced bandwidth and the spin-dependence of $t^{\sigma}_{\rm eff}$ enhance the average magnetization in RIMT calculations. Within the IMT framework, on the other hand, the renormalization of hopping parameter ($t_{\rm eff}\approx 1.13$ at $h=0$) and its spin dependence is negligible, as seen in Fig. 6(b), arising only from the Fock-shifts in ${\cal H}_{\rm MF}$. The externally applied magnetic field $h$ also gets renormalized to an effective magnetic field $h_{\rm eff}$ ($=\sum_{i\sigma}\sigma{\mu}^{\sigma}_{\rm eff}(i)/2$) in both RIMT and IMT by Hartree shifts of ${\cal H}_{\rm{MF}}$, defined in Eq. 7. The variation of $h_{\rm eff}$ with respect to the external field $h$ is depicted in Figs. 6(c) and (d) from RIMT and IMT methods respectively. The suppression is significant in RIMT due to the additional action of $\phi_{i\sigma}$ i.e, the derivatives of GRFs, e.g. $(dg^{J,z}/dn)$, consumed within ${\mu}^{\sigma}_{\rm eff}$ in RIMT calculations. $t_{\rm eff}(i)$ locally plays a key role in homogenizing small scale inhomogeneities in the spatial profiles of the densities of the systemChakraborty and Ghosal (2014).

III.4 Density of states

Another key feature of a superconducting state is its single particle density of states (DOS), which carry specific signatures when $h$ is turned on. In order to explore the effects of correlations on the DOS, we study it at different values of $h$ within the RIMT and IMT schemes. The spin-resolved DOS are evaluated using BdG eigenvalues {$E_{n\sigma}$} and eigenvectors {$u^{n}_{i,\sigma}$,$v^{n}_{i,\sigma}$} as,

[TABLE]

DOS characterizing different phases obtained by tuning $h$ from RIMT and IMT methods are shown in Fig. 7 on the left and right columns respectively. For the d-wave BCS state at a low $h(<h_{1})$, the DOS of the two spin flavors split in which the up-spin DOS gets shifted towards the left and the down-spin towards the right with respect to the Fermi level by an amount of $h_{\rm eff}$ as shown in Figs. 7(a) and (b).

For intermediate magnetic fields $h_{1}\leq h<h_{2}$, the FFLO state is identified by a mid-gap peakLoh and Trivedi (2010); Zhou and Ting (2009) appearing at $\omega=\mp h_{\rm eff}$ for up- and down-spin DOS respectively. This is because, the paired states in the FFLO phase reside near the Zeeman-split Fermi surfaces and the single particle states, which cause finite magnetization, occupy the energies in between them. In real space, the single particle states are piled up at the zeros of the superconducting pairing amplitude $\Delta$ where it changes sign and form domain walls. The near ”square-wave” nature of the $\Delta$ modulation within RIMT, as seen in Fig. 5, supports domain walls within narrow regions in the real space. The resulting spatial profile of $m$ thus features strong peaks at these narrow domain walls. Such strong localization of the single-particle states conduce a rather sharp mid-gap feature in the resulting density of states from RIMT calculations, as seen in Fig. 7(c). These mid-gap states and the corresponding mid-gap peak in DOS are reminiscent of the bound states formed due to the Andreev reflections along the nodal lines of a superconductor. Note that the superconducting order parameter changes its sign on the nodal line. The nature of modulation of pairing amplitude in IMT scheme, however, maintains near-sinusoidal form (higher harmonics less relevant) and as a result, the mid-gap feature in corresponding density of states is much less sharp, as can be seen from Fig. 7(b).

As $h$ increases, the mid-gap peak in the FFLO state broadens with the increased separation of the two Fermi surfaces. In real space, the domain walls get closer to each other with increasing $h$. This facilitates stronger hybridization of the Andreev bound states, resulting larger bandwidth of the mid-gap peak. The low energy states with respect to the up-spin and down-spin Fermi surfaces fill up with increasing field, until the FFLO to NS transition occurs, where the gap disappears completely.

The sharpness of the mid-gap feature of DOS in the RIMT results, particularly at lower magnetic fields at $h=0.18$ in Fig. 7(c), is also due to the reduction of bandwidth to $8t_{\rm eff}$ ($t_{\rm eff}\approx 0.4$ at $h=0$). This reduction further becomes spin-dependent in the presence of $h$, as shown in Fig. 6(a). The reduced bandwidth also make the DOS better resolved with closely spaced energy levels in RIMT.

The sharp change of the pairing amplitude and sign near the zeros in RIMT results, makes the wave-function corresponding to the mid-gap energy far more localized near the domain walls compared to the IMT findings. This is shown for the lowest lying wave-function $|\phi_{n}(\bf{r})|^{2}$ (here, $E_{n}\approx-h_{\rm eff}$) in Figs. 8(a) and (b) from RIMT and IMT calculations respectively to highlight their contrast. The boundedness of the low-lying wave-functions reduce as the steepness of the $\Delta$ modulation decreases with the increase in $h$.

III.5 Fate of FFLO phase in the presence of competing order

Our results in the previous subsections illustrate how an application of Zeeman field $h$ generates an FFLO state from a pristine d-wave superconducting GS (for $h=0$). One of the hallmarks of most strongly correlated superconductors is that they often carry translational symmetry broken orders in their GS. For example, charge ordersCampi et al. (2015); Peng et al. (2016) and antiferromagnetismSachdev (2012); Lu et al. (2014) in cuprate superconductors, spin-density wave (SDW) order in heavy-fermion superconductorsHu et al. (2012); Bianchi et al. (2003b) and Fe-based superconductorsSchmiedt et al. (2014); Si et al. (2016). Motivated by this, we consider a broken symmetry GS with d-wave superconducting order and a commensurate SDW orderYanase and Sigrist (2009b) and scan the phase space traced by $h$. The SDW order is rife with the d-wave superconducting order and often found in the strongly correlated superconductors at small doping values. We begin with the Hamiltonian in Eq. $(\ref{hamil})$, and its spin rotational symmetry gets broken by the SDW order even at $h=0$.

For this calculation, we fix $\langle n\rangle=0.9$. With this average density, we focus on addressing the following question: can a translational symmetry breaking SDW order in the GS at $h=0$ perturb the modulating pairing amplitude in the FFLO phase at finite $h$? The phase diagram including SDW order using RIMT and IMT schemes are shown in Figs. 9(a) and (b) respectively. Here, we allowed the following competing states to emerge at different $h$ values: d-wave BCS state, FFLO state, d-wave BCS state coexisting with SDW order (d-wave BCS+SDW), a state with FFLO modulation coexisting with this SDW order (we coin this as FFLO+SDW state), normal state (NS) with itinerant magnetism and a NS with SDW order (NS-SDW). The final energetics at BdG self-consistency determines the true GS.

We compare the phase diagrams including the competing SDW order for both IMT and RIMT calculations. For a legitimate comparison, we fix $J\approx 1.47$ in IMT calculations. This fetches the same d-wave pairing gap in the presence of this SDW order at $h=0$ in RIMT and IMT. The phases in the GS and the phase boundaries, as obtained from the energetics, are shown in Fig. 9. The calculations for the case with strong correlation is done in a 40$\times$40 lattice, yielding weaker $q$ resolution than our results obtained in momentum space on systems of size 200$\times$200. Our crucial finding from this study is that, the GS for all $h$, shown in Fig. 9, are essentially those obtained in Fig. 1, but in addition accommodate the competing SDW oder, as we discuss below. Note that, while the SDW order is put in by hand at $h=0$, it survives in the self-consistent GS for all $h$ in both RIMT and IMT, as shown in panel (a) and (b) of Fig. 9.

We find that the d-wave BCS+SDW state (dark blue curves) and the NS-SDW state (dark green curves) are always energetically favorable over the FFLO state (red curves) generated from the pristine dSC state in both RIMT and IMT, as shown in Fig. 9. However, we further find that, the FFLO+SDW state energetically survives for a window of $h$ in the phase diagram. This state carries a modulated superconducting order parameter like in the FFLO state, with an additional $(\pi,\pi)$ modulation in the magnetization. The spatial modulation of different order parameters are shown in Fig. 10 on the left side with IMT results, and on the right side with RIMT results. The upper and lower panel figures correspond to two strengths of $h$. The FFLO+SDW state turns out to be the lowest energy state sandwiched between the d-wave BCS+SDW (at low $h\leq h_{1}$) and NS-SDW (at high $h\geq h_{2}$) states. In RIMT scheme, FFLO+SDW phase ranges from $h_{1}\approx 0.23$ to $h_{2}\approx 0.45$ and in IMT this region ranges from $h_{1}\approx 0.2$ to $h_{2}\approx 0.43$. The balances of energy gain and loss from the individual components of the total energy deciding the boundaries of the FFLO+SDW phase for RIMT and IMT schemes are shown in Fig. 12 in App. B. In fact, this window in $h$ where FFLO+SDW is the ground state appears wider compared to the window where FFLO was energetically favorable if SDW order was ignored (as shown in Fig. 1) for both the RIMT and IMT calculations.

Thus, our toy calculation indicates that the signatures of modulating pairing amplitude, i.e. the impression of FFLO survives with competing orders in the corresponding GS at $T=0$.

IV Conclusion

In conclusion, we have studied the effects of strong electronic correlations in the FFLO state of a d-wave superconductor. Thus our results are of great relevance for the search of FFLO signatures in strongly correlated d-wave superconductors, such as, $\rm{CeCoIn_{5}}$ and the cuprates. Our findings indicate that the strong correlations renormalize all relevant energy scales, whose intricate balance decides the phase space for the FFLO state. Consequently, we found an increased window of the magnetic field for the FFLO phase. We make definitive predictions for the behaviors of the order parameters, pairing momenta and the DOS – all feature interesting distinctions between RIMT and IMT findings. In RIMT, pairing momentum rises sharply from zero near the lower critical field and remains nearly saturated over a large part of the FFLO phase. Such near-saturation of modulating wave-vector has also been reported in scattering experiments done in the Q phase of $\rm{CeCoIn_{5}}$ Kenzelmann et al. (2010). Though, the nature of the Q phase is a subject of current debateRaymond et al. (2014); Martiny et al. (2015), it shares broad similarity with a putative FFLO state. Strong interactions were found to homogenize small-scale variations in the $\Delta$ landscape, which cause it to change sign rather sharply near $\Delta=0$. This, in turn, localizes a high density of Andreev bound states on these domain regions, leading to a narrow and sharp mid-gap peak in the density of states of FFLO phase within RIMT. This is consistent with the recent NMR experiment on $\rm{CeCu_{2}Si_{2}}$, which is suggestive for the presence of high density of Andreev bound states in its inhomogeneous superconducting state in Zeeman field Kitagawa et al. (2018).

A natural question might arise: Why do we not encounter FFLO phase in strongly correlated superconductors, e.g. high $T_{c}$-cuprate superconductors? We already mentioned that the Q phase, a ‘cousin’ of the FFLO phase, has already been observed in $\rm{CeCoIn_{5}}$ Kenzelmann et al. (2010). In addition, we have not considered the orbital effects of the applied field in our analysis, assuming only the Zeeman effects. In reality, it is challenging to disentangle the orbital and Zeeman effects of an applied field. In particular, for cuprate superconductors, the orbital effects produces vortices at weaker field strengths ($H^{\rm{orb}}_{c_{2}}\sim 100\rm{T}$ for YBCO near the optimal dopingGrissonnanche et al. (2014)) compared to the Clogstron-Chandrasekhar limit ($H_{p}\sim\rm{T}$ for YBCO (within BCS theory) near the optimal dopingO’Brien et al. (2000)), where Zeeman effects become crucial! Thus, it is quite possible that homogeneous superconductivity might become completely disordered pre-empting FFLO modulations due to proliferation of vortices!

On the other hand, Pauli-limited superconductors possessing large Maki parametersSaint-James et al. (1969) ($>1.8$), which also are strongly correlated, such as heavy-fermion superconductors $\rm{CeCoIn_{5}}$, $\rm{CeCu_{2}Si_{2}}$ and some of the organic superconductors show signatures of FFLO phase, when exposed to magnetic field. We also have not included the effects of quantum phase fluctuations in either of our recipe: IMT and RIMT. It will be interesting to explore the effects of quantum phase fluctuation on the FFLO physics, particularly at lower doping where strong correlation effects are significant.

V Acknowledgements

The authors acknowledge computational facilities at IISER Kolkata. AD acknowledges her fellowship from UGC (India).

Appendix A Obtaining renormalized mean-field Hamiltonian ${\cal H}_{\rm{MF}}$ of Eq. $(\ref{MF})$

A.1 Detailed expressions for GRFs used in Eq. $(\ref{hamil})$Christensen et al. (2011); Yang et al. (2009) as functions of the order parameters defined in Eq. $(\ref{op1},\ref{op3},\ref{op6})$

In the process of handling the constraint of strong on-site repulsion ($U\gg t$) in the Hubbard model at our starting point in Sec. II, we employed Gutzwiller projection operator to focus on a restricted Hilbert space, which prohibits all the double occupancies from the system. While this reduces Hamiltonian in Eq. $(\ref{Hubb})$ to that in Eq. $(\ref{tJ})$, standard manipulation demands further simplifications for the implementation of the constraints. One intuitive and elegant (through approximation) way for implementing the constraint is called Gutzwiller approximation, which actually gets rid of the constraints in expense of renormalization of the Hamiltonian parameters locally. The resulting Hamiltonian is given in Eq. $(\ref{hamil})$. Here, the Gutzwiller renormalization factors (GRFs), i.e, all g’s are local variables, to be determined self-consistently. The form of such GRFs can be derived from a phase space argumentEdegger et al. (2007) or from infinite dimensional calculationsGebhard (1990) and they also depend on the broken symmetry ground state that we wish our ground state would describe. For our case, we wish to accommodate d-wave superconducting order as well as magnetization in our ground state and within such premises, the explicit form of GRFs have been worked out in Ref. Ogata and Himeda, 2003. For completeness we list below the expressions of these GRFs in terms of all the order parameters.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, $\delta{i}=1-n_{i}$, where $n_{i}=\sum_{\sigma}n_{i\sigma}$; $\Delta_{ij}=\sum_{\sigma}\Delta_{ij\sigma}/2$; $\tau_{ij}=\sum_{\sigma}\tau_{ij\sigma}/2$.

A.2 Details of mean-field decomposition of ${\cal H}_{\rm{t-J}}$ defined in Eq. $(\ref{tJ})$

The method of derivation of the mean-field Hamiltonian of nature as given in Eq. (7) starting from the bare Hamiltonian like in Eq. (3) is standard and can be found in the literature Chakraborty and Ghosal (2014). Here, we proceed to describe the procedure specifically for our investigation of FFLO state, for the sake of completeness.

We minimize $\langle\psi_{0}|{\cal H}_{\rm{GA}}|\psi_{0}\rangle$ with respect to $|\psi_{0}\rangle$ (Sec. II), under the constraints of fixed total electron density $N^{-1}\sum_{i}n_{i}=\langle n\rangle$ and normalization of the wavefuntion $\langle\psi_{0}|\psi_{0}\rangle=1$ or equivalently we minimize the functional $W=\langle\psi_{0}|{\cal H}_{\rm{GA}}|\psi_{0}\rangle-\lambda\left(\langle\psi_{0}|\psi_{0}\rangle-1\right)-\mu\left(\sum_{i}n_{i}-\langle n\rangle\right)$ as follows:

[TABLE]

which leads to the renormalized mean-field Hamiltonian ${\cal H}_{\rm{MF}}$ of Eq. $(\ref{MF})$.

A.3 Expressions of $\xi^{(\rm r)}_{k+Q,\sigma}$ and ${\Delta}^{(r)}_{k,-k+Q}$ of Eq. $(\ref{kmf})$

For the periodic and clean system we have considered, ${\cal H}_{\rm{MF}}$ in Eq. (7) can be best solved in momentum space after performing Fourier transformation, as obtained in Eq. $(\ref{kmf})$. Here, we present the detailed expressions $\xi^{(\rm r)}_{k+Q\sigma}$, and ${\Delta}^{(r)}_{k,-k+Q}$ appearing in Eq. $(\ref{kmf})$.

[TABLE]

[TABLE]

Here, $\gamma_{k}=2(\cos(k_{x})+\cos(k_{y}))$; $\eta_{k}=2(\cos(k_{x})-\cos(k_{y}))$; $g^{t\sigma}_{0,\delta}=(1/4N)\sum_{i,\delta}g^{t\sigma}_{i\delta}$; $g^{J,xy}_{0}=(1/4N)\sum_{i,\delta}g^{J,xy}_{i\delta}$ ; $g^{J,z}_{0}=(1/4N)\sum_{i,\delta}g^{J,z}_{i\delta}$; $g^{t\sigma}_{Q\delta}=N^{-1}\sum_{i}g^{t\sigma}_{i\delta}\cos(\bf{Q}.\bf{r}_{i})$; $g^{J,xy}_{Q\delta}=N^{-1}\sum_{i}g^{J,xy}_{i\delta}\cos(\bf{Q}.\bf{r}_{i})$; $g^{J,z}_{Q\delta}=N^{-1}\sum_{i}g^{J,z}_{i\delta}\cos(\bf{Q}.\bf{r}_{i})$; $n^{\sigma}_{Q}=N^{-1}\sum_{k}\langle\hat{c}^{\dagger}_{k\sigma}\hat{c}_{k-Q,\sigma}\rangle_{0}$; $\tau^{\sigma}_{Q}=(1/4N)\sum_{k}\langle\hat{c}^{\dagger}_{k\sigma}\hat{c}_{k-Q,\sigma}\rangle_{0}\gamma_{k}$; $\tau^{\sigma}_{Q,\delta}=N^{-1}\sum_{k}\langle\hat{c}^{\dagger}_{k\sigma}\hat{c}_{k-Q,\sigma}\rangle_{0}\cos(\mathrm{k}.\mathrm{\delta})$; $\Delta^{\prime}_{Q}=(1/4N)\sum_{k}\langle\hat{c}^{\dagger}_{-k+Q\uparrow}\hat{c}^{\dagger}_{k\downarrow}\rangle_{0}\eta_{k}$ ; $\Delta^{\prime}_{Q,\delta}=N^{-1}\sum_{k}\langle\hat{c}^{\dagger}_{-k+Q\uparrow}\hat{c}^{\dagger}_{k\downarrow}\rangle_{0}\cos(\mathrm{k}.\mathrm{\delta})$. $\delta$ denotes nearest-neighbor spacing.

Appendix B Evolutions of the energy gain and loss with respect to the magnetic field at the phase boundaries

Our main result in Sec. III showed that the phase boundaries between BCS, FFLO and NS move around depending on the inclusion of strong correlations in the fold of the calculation. It was also argued that the signature of strong correlations renormalize different components of energy in a different manner, such that the subtle balance between these components are achieved at different strengths of Zeeman field $h$, causing the phase boundaries to be different for RIMT and IMT. Here, we illustrate the above statement in the following manner in terms of our results.

The differences of the energy components of the competing phases without the competing SDW order near their phase boundaries with respect to $h$ obtained from the RIMT and IMT scheme are shown in Fig. 11. The pairing energy ($E_{p}$) favors a uniform BCS phase than the FFLO phase, because the spatial modulation of order parameter comes for an energy cost, most easily seen from a Ginzburg-Landau expansion of free energyTinkham (2004). In contrast, the magnetization energy ($E_{m}$) and effective kinetic energy cost ($E_{K}$) is better accommodated in the FFLO phase, because, of the nucleation of the $E_{m}$ and $E_{K}$ at the domain walls formed along the nodes of $\Delta$ in real space in the FFLO phase. Because of the GRFs appearing in the Hamiltonian ${\cal H}_{\rm MF}$ in Eq. (11), which are self-consistently determined for specific location in the parameter space, independent components of energies evolve differently with $h$. As a result, the changeover from BCS to FFLO and subsequently from FFLO to normal state can occur at different $h_{1}$ and $h_{2}$ in principle, from RIMT and IMT calculations. The crossing of $(E_{m}+E_{K})$ gain and $E_{p}$ loss in the FFLO phase as found at $h_{1}$ reduces in RIMT, because, the rate of change of energy gain in FFLO due to $(E_{m}+E_{K})$ increases compared to the energy loss from $E_{p}$ near $h_{1}$. This is expected because strong correlations homogenize the small scale variations in $\Delta$ in the FFLO phase and makes the zero region of $\Delta$ narrower by steepening the fall of $\Delta$ near the same as depicted in Fig. 5. Therefore, even near a reduced $h_{1}$, FFLO state becomes energetically favorable compared to the BCS state as the loss of $E_{p}$ in the FFLO state is lesser compared to the gain in $E_{m}$ which mainly stems from the increased magnetization due to lowering of bandwidth in RIMT. On the other hand, inside the FFLO region the gap filling or in other words, the decay of $E_{p}$, ($E_{p}\approx 0$ determines $h_{2}$) occurs at a relatively slower rate in RIMT due to renormalization of the effective magnetic field $h_{\rm eff}$ as shown in Fig. 6(c) and therefore $h_{2}$ does not shift in a similar fashion as $h_{1}$. Fig. 12 shows the balances in the energy components at the phase boundaries of the competing phases with the competing SDW order from the RIMT and IMT calculations. The competing phases here are dwave BCS+SDW phase, FFLO+SDW phase and the underlying normal state. Near $h_{1}$, the FFLO+SDW phase is favored by $E_{m}$ and dwave BCS+SDW phase is favored by $E_{p}$ and $E_{K}$. Renormalization of the parameters due to Gutzwiller factors and as a result of that the renormalized energy components causes different $h_{1}$ values in RIMT and IMT. $h_{2}$ in the two cases are also different due to different renormalization of the parameters in RIMT and IMT.