Ground states of atomic Fermi gases in a two-dimensional optical lattice with and without population imbalance

Here we consider a two-component ultracold Fermi gas with a short-range pairing interaction, $V_{\mathbf{k,k^{\prime}}}=U<0$, in 2DOL. The dispersion of noninteracting atoms without population imbalance is given by $\xi_{\mathbf{k}}=\epsilon_{\mathbf{k}}-\mu\equiv{k_{z}}^{2}/2m+2t[2-\cos(k_{x}d)-\cos(k_{y}d)]-\mu$, where $k_{z}$ is the momentum in the $z$ direction in the continuum dimension, $k_{x}$ and $k_{y}$ are the momenta in the lattice plane, $t$ and $d$ are the hopping integral and lattice constant in the $xy$ plane, respectively, and $\mu$ is the chemical potential. Following our recent works [39, 35, 15, 40], we take $t$ to be physically accessible, under the constraint $2mtd^{2}<1$ in our calculation. And the critical coupling for forming a two-body bound state of zero binding energy is given by $U_{\text{c}}=-1/\sum_{\textbf{k}}1/2\epsilon_{\textbf{k}}=-0.16072\sqrt{2m}/\sqrt{t}d^{2}$. Here and throughout we take the natural units and set $\hbar=k_{\text{B}}=1$.

At zero temperature, the mean-field BCS-Leggett ground state follows the gap and number equations [37]

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\frac{1}{U}+\sum_{\textbf{k}}\frac{1}{2E_{\textbf{k}}}\,,$		(1)
	$\displaystyle n$	$\displaystyle=$	$\displaystyle\sum_{\textbf{k}}\Big{(}1-\frac{\xi_{\textbf{k}}}{E_{\textbf{k}}}\Big{)}\,,$		(2)

where $E_{\textbf{k}}=\sqrt{\xi_{\textbf{k}}^{2}+\Delta^{2}}$ is the Bogoliubov quasiparticle dispersion, with an energy gap $\Delta$.

To make sure the mean-field solution is stable, we impose the requirement that the dispersion of the Cooper pairs be non-negative, both in the lattice plane and along the $z$ direction. To this end, we extract the inverse pair mass (tensor) using the fluctuating pair propagator, as given in the pairing fluctuation theory which was previously developed for the pseudogap physics in the cuprates [41] and extended to address the BCS-BEC crossover in ultracold atomic Fermi gases [1]. In particular, we mention that, compared to rival $T$-matrix approximations for the pairing physics, the pair dispersion as extracted from this theory is gapless below $T_{\text{c}}$, fully compatible with the mean-field gap equation. Here the pairing $T$ matrix is given by $t_{\text{pg}}(Q)=U/[1+U\chi(Q)]$, with the pair susceptibility $\chi(Q)=\sum_{K}G_{0}(Q-K)G(K)$, the bare Green’s function $G_{0}(K)=(\omega-\xi_{\mathbf{k}})^{-1}$, and the full Green’s function $G(K)=\frac{u_{\textbf{k}}^{2}}{\omega-E_{\textbf{k}}}+\frac{v_{\textbf{k}}^{2}}{\omega+E_{\textbf{k}}}$, where $u_{\textbf{k}}^{2}=(1+\xi_{\textbf{k}}/E_{\textbf{k}})/2$ and $v_{\textbf{k}}^{2}=(1-\xi_{\textbf{k}}/E_{\textbf{k}})/2$ are the BCS coherence factors, and $K\equiv(\omega,\mathbf{k})$, $Q\equiv(\Omega,\mathbf{q})$ are four momenta.

The inverse $T$-matrix $t_{\text{pg}}^{-1}(Q)$ can be expanded for small $Q$, given by $t_{\text{pg}}^{-1}(\Omega,\textbf{q})\approx a_{1}\Omega^{2}+a_{0}(\Omega-\Omega_{\textbf{q}}+\mu_{p})$, with $\Omega_{\textbf{q}}=B{q_{z}}^{2}+2t_{B}[2-\cos(q_{x}d)-\cos(q_{y}d)]$, and $\mu_{p}=0$ in the superfluid phase. Then we extract $B=1/2M$, with $M$ being the effective pair mass in the $z$ direction, and $t_{B}$ the effective pair hopping integral in the $xy$ plane. The sign of $a_{0}$ determines whether the fermion pairs are particle-like or hole-like, with positive $a_{0}$ for particle-like pairing and negative $a_{0}$ for hole-like pairing. For example, in a 3D lattice, in general one finds $a_{0}>0$ for fermion density below half filling, $a_{0}=0$ at half filling due to particle-hole symmetry, and $a_{0}<0$ above half filling. The sign of $a_{0}$ is controlled by the average of the inverse band mass ¹¹1For the latter case, one can perform a particle-hole transfermation so that it becomes $a_{0}>0$ and below half filling for holes.. While one could perform a particle-hole transformation for a pure lattice case, it does not seem to be feasible in our case since both lattice and continuum dimensions are present. The expressions for the coefficients $a_{1}$, $a_{0}$, $B$ and $t_{B}$ can be readily derived during the Taylor expansion. In this way, using the solution for ($\mu$, $\Delta$) from Eqs. (1)-(2), we can extract the pair dispersion $\tilde{\Omega}_{\textbf{q}}=(\sqrt{a_{0}^{2}+4a_{1}a_{0}\Omega_{\textbf{q}}}-a_{0})/2a_{1}$. The non-negativeness of the pair dispersion implies that the pairing correlation length (squared) $\xi^{2}=a_{0}B$ and $\xi_{xy}^{2}=a_{0}t_{\text{B}}d^{2}$ must be positive.

For the population imbalanced case, the spin polarization is defined via $p=(n_{\uparrow}-n_{\downarrow})/(n_{\uparrow}+n_{\downarrow})$, where spin index $\sigma=\uparrow,\downarrow$ refers to the majority and minority components, respectively. Then the dispersion of noninteracting atoms is modified as $\xi_{\mathbf{k}\sigma}=\epsilon_{\mathbf{k}}-\mu_{\sigma}$ , with $\mu_{\sigma}$ the chemical potential for spin $\sigma$.

Now the bare and full Green’s functions are given by

	$\displaystyle G_{0\sigma}(K)$	$\displaystyle=$	$\displaystyle\frac{1}{\omega-\xi_{\mathbf{k}\sigma}}\,,\quad\text{and}$
	$\displaystyle G_{\sigma}(K)$	$\displaystyle=$	$\displaystyle\frac{u_{\textbf{k}}^{2}}{\omega-E_{\textbf{k}\sigma}}+\frac{v_{\textbf{k}}^{2}}{\omega+E_{\textbf{k}\bar{\sigma}}},$

respectively, where $\bar{\sigma}$ is the opposite spin of $\sigma$, $E_{\textbf{k}\uparrow}=E_{\textbf{k}}-h$, and $E_{\textbf{k}\downarrow}=E_{\textbf{k}}+h$, with $\mu=(\mu_{\uparrow}+\mu_{\downarrow})/2$, and $h=(\mu_{\uparrow}-\mu_{\downarrow})/2$. Thus $E_{\textbf{k}\uparrow}$ becomes gapless, as it should, in order to accommodate the excessive majority fermions [See Eq. (5) below]. These gapless fermions will contribute in both the gap and number equations.

Following the BCS self-consistency condition and the number constraint, we arrive at the gap and number equations at zero $T$ in the presence of population imbalance:

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\frac{1}{U}+\sum_{\textbf{k}}\frac{\Theta(E_{\textbf{k}\uparrow})}{2E_{\textbf{k}}}\,,$		(3)
	$\displaystyle n$	$\displaystyle=$	$\displaystyle\sum_{\mathbf{k}}\Big{[}\Big{(}1-\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\Big{)}+\Theta(-E_{\textbf{k}\uparrow})\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\Big{]}\,,$		(4)
	$\displaystyle pn$	$\displaystyle=$	$\displaystyle\sum_{\mathbf{k}}\Theta(-E_{\textbf{k}\uparrow})\,,$		(5)

where $\Theta(x)$ is the Heaviside step function, and $n=n_{\uparrow}+n_{\downarrow}$ and $\delta n=n_{\uparrow}-n_{\downarrow}=pn$ are the total and the difference of fermion densities, respectively.

In the imbalanced case, the pair susceptibility is modified as $\chi(Q)=\sum_{K,\sigma}G_{0\sigma}(Q-K)G_{\bar{\sigma}}(K)/2$, which is consistent with the BCS self-consistency condition so that the pair dispersion remains gapless at $q=0$. Then we follow the same procedure as in the balanced case, and extract the inverse pair mass tensor along with coefficients $a_{0}$ and $a_{1}$ via the Taylor expansion of the inverse $T$ matrix, $t_{\text{pg}}^{-1}(Q)$.

Equations (3)-(5) form a closed set of self-consistent equations, and can be used to solve for ($\mu$, $h$, $\Delta$) as a function of ($U$, $t$, $d$, $p$), which is then further constrained by various stability conditions.

As shown in the 3D continuum and 1DOL cases, in the presence of population imbalance, not all solutions of Eqs. (3)-(5) are stable [19, 9, 43].

Following the stability analysis of Refs. [19, 9], the stability condition for the superfluid phase requires that for fixed $\mu$ and $h$, the solution for the excitation gap $\Delta$ is a minimum of the thermodynamic potential $\Omega_{\text{S}}$, which is demonstrated to be equivalent to the positive definiteness of the generalized compressibility matrix [9, 44]. Thus we have

$$\frac{\partial^{2}\Omega_{\text{S}}}{\partial\Delta^{2}}=\sum_{\mathbf{k}}\frac{\Delta^{2}}{E_{\mathbf{k}}^{2}}\Big{[}\frac{\Theta(E_{\textbf{k}\uparrow})}{E_{\mathbf{k}}}-\delta(E_{\textbf{k}\uparrow})\Big{]}>0\,,$$

(6)

where $\delta(x)$ is the delta function.

In addition, the positivity of the pair dispersion in the entire momentum space imposes another strong stability condition. Illustrated in Fig. 1 are the qualitative behaviors of the pair dispersion, for different signs of $a_{0}$ and $\xi^{2}$. For illustration purpose, a simple isotropic quadratic dispersion is assumed. In general, there are two branches of the dispersion, from the inverse $T$-matrix expansion up to the $\Omega^{2}$ order. The positive branch represents a propagating mode, while the negative branch represents a hole-like mode which contributes to quantum fluctuations. The case of $a_{0}>0$ and $\xi^{2}>0$ (Fig. 1(a)) corresponds to particle-like pairing, with a monotonically increasing energy and a positive effective pair mass, $B>0$ and $t_{\text{B}}>0$, so that $q=0$ is the bottom of the pair energy. For the $a_{0}<0$ case (Fig. 1(c)), this dispersion flips upside down into the blue-dashed hole mode. This corresponds to hole-like pairing, for which $q=0$ becomes a local maximum, with $B<0$ and $t_{\text{B}}<0$, similar to the hole band in a semiconductor. In case of a pure lattice, one could flip the sign of $a_{0}$ via a particle-hole transformation so that this blue-dashed line is flipped back to become positive as the dispersion for hole pairs. However, for our present case, due to the presence of the continuum dimension, there is no easy way to do a particle-hole transformation so that we have to stay with the (black solid) gapped positive branch, which is a flip of the hole branch in Fig. 1(a), as the dispersion of particle-like Cooper pairs. When $a_{0}=0$, the two branches become symmetric, without a gap. For all three cases, the coefficients of the $q^{2}$ terms in the inverse $T$ matrix expansion, $\xi^{2}$ and $\xi_{xy}^{2}$, must be positive. (Note that $a_{1}$ is always positive.) Indeed, as shown in Figs. 1(d-f), for a negative $\xi^{2}$, the dispersion $\tilde{\Omega}_{\mathbf{q}}$ of both particle-like (Figs. 1(d)) and hole-like (Figs. 1(f)) pairs quickly become diffusive and thus cease to exist, unless higher order terms, e.g., the $q^{4}$ terms, are included. In that case, the pair dispersion will reach a minimum at a non-zero $q$. Our numerics shows that in 2DOL, $\xi^{2}$ in the continuum dimension remains positive in general but $\xi_{xy}^{2}\propto a_{0}t_{\text{B}}$ in the lattice plane may indeed change sign so that $\xi_{xy}^{2}>0$ will constitute another stability requirement for the superfluid phase.

Finally, the superfluid density must also be positive definite in a stable superfluid. This, however, has been found to be a relatively weaker constraint in the cases of 3D continuum [19, 9].

As a representative transport property, superfluid density is an important quantity in the superfluid phase. While it is always given by $n/m$ at zero $T$ for the balanced case in 3D continuum, it will take the average of the inverse band mass in the presence of a lattice. Furthermore, in the presence of population imbalance, it may become negative [19, 9, 45], signaling an instability of the superfluid state. Here we shall also investigate the behavior of the anisotropic superfluid density $(n_{\text{s}}/m)$, and pay close attention to the population imbalanced case and the situations where it becomes negative.

The expression for superfluid density can be derived using the linear response theory. Following Refs. [19, 9, 41, 45, 46], we obtain for zero $T$

$$\left(\frac{n_{\text{s}}}{m}\right)_{i}=\sum_{\mathbf{k}}\frac{\Delta^{2}_{\text{sc}}}{E_{\mathbf{k}}^{2}}\left[\frac{\Theta(E_{\textbf{k}\uparrow})}{E_{\mathbf{k}}}-\delta(E_{\textbf{k}\uparrow})\right]\left(\frac{\partial\xi_{\mathbf{k}}}{\partial{k}_{i}}\right)^{2}\,,$$

(7)

where $i=x,y$ and $z$ for the lattice and the continuum directions, respectively.

Fermi surface plays an important role in the superfluid and pairing behavior of atomic Fermi gases. For 2DOL, it is very different from the 3D continuum or 3D lattice case, so is it from 1DOL [15, 39, 35]. This will lead to different physics. Here we first present the shape and topology of the Fermi surface for a series of representative sets of lattice parameters ($t,d$). Shown in Fig. 2 is the typical evolution behavior of the Fermi surface, calculated self-consistently in the noninteracting limit at zero temperature. The top row shows the evolution with the lattice constant, for $k_{\text{F}}d=1,2,3$ and $4$ at fixed hopping integral $t/E_{\text{F}}=0.05$. Then the bottom row shows the effect of hopping integral, with $t/E_{\text{F}}=0.01$, 0.04, 0.07, and $0.1$ and fixed $k_{\text{F}}d=3$.

The lattice constant $d$ provides a confinement in the momentum space; the larger $d$ the stronger confinement. The top row in Fig. 2 suggests that the Fermi surface becomes thicker along the $z$ direction as $d$ increases for fixed $t$. Indeed, fermions feel a stronger confinement in the lattice dimensions with a shrinking first Brillouin zone (BZ), as $k_{\text{F}}d$ increases from 1 to 4, and thus need to occupy higher $k_{z}$ states to keep the Fermi volume unchanged, so that the noninteracting fermionic chemical potential is pushed up. As a rough estimate, the maximum occupied $k_{z}$ increases by a factor of 16 from left to right. For relatively small $t/E_{\text{F}}=0.05$, the shape and topology of the Fermi surface evolve from a closed plate for $k_{\text{F}}d=1$ into one with only the top and bottom faces while completely open on the four sides at the BZ boundary of the lattice dimensions for $k_{\text{F}}d=3$ and 4. For the intermediate $k_{\text{F}}d=2$, the Fermi surface is open only at the center of the four side faces at the BZ boundary. At the same time, the effective filling factor in the lattice dimensions increases to nearly unity as $k_{\text{F}}d$ increases from 1 to 4. In this way, for large $d$, fermion dispersion on the Fermi surface on average becomes hole-like in the lattice plane, while it always remains particle-like in the continuum dimension.

On the other hand, a smaller $t$ makes the fermion energy less dispersive in the lattice dimensions, and thus the lattice band becomes narrower and more fully filled. In other words, fermions will tend not to go to higher $k_{z}$ states until the BZ at lower $k_{z}$ is fully occupied, leading to a flatter top and bottom of the Fermi surface. This will also pull down the noninteracting fermionic chemical potential. As shown in the bottom row in Fig. 2, the Fermi surface becomes thinner and flatter in the $z$ direction as $t/E_{\text{F}}$ decreases from $0.1$ to $0.01$ for fixed $k_{\text{F}}d=3$. In contrast, the $t/E_{\text{F}}=0.07$ and $0.1$ cases have a much more dispersive Fermi surface as a function of the in-plane momentum $(k_{x},k_{y})$. Fermions at high $(k_{x},k_{y})$ states are removed for relatively large hopping integral $t/E_{\text{F}}=0.07$ and $0.1$.

The evolution of the Fermi surface reveals that the in-plane fermion motion on the Fermi surface becomes hole-like for relatively small $t$ and large $d$. As a result, the nature of the in-plane and overall pairing in this case will also change from particle-like to hole-like when the contributions from lattice dimensions are dominant in the BCS and unitary regimes [36].

It should be mentioned that in the strong pairing regime, the detailed shape of the Fermi surface is no longer relevant, as pairing extends essentially to the entire momentum space. However, the confinement in the momentum space imposed by the lattice periodicity is always present and will govern the physical behavior in the BEC regime.

It is known from the 3D continuum case that the balanced case and the imbalanced case with $p\rightarrow 0^{+}$ are not continuously connected in the BCS and unitary regimes at $T=0$ [19, 47]. Population imbalance leads to very distinct behaviors. Therefore, we present in this section the balanced results only.

In Fig. 3, we present a typical phase diagram (a) in the $d$ – $U$ plane, for fixed relatively small $t/E_{\text{F}}=0.05$, and (b) in the $t$ – $U$ plane, for relatively large $k_{\text{F}}d=3$, corresponding to the cases of the top and bottom rows in Fig. 2, respectively. The lattice constant in panel (a) ranges from relatively small $k_{\text{F}}d=1$ with $2mtd^{2}=0.05$ to the upper limit $k_{\text{F}}d=2\sqrt{5}$ with $2mtd^{2}=1$ denoted by the horizontal (cyan) dotted line, and the hopping integral in panel (b) ranges from relatively small $t/E_{\text{F}}=0.01$ with $2mtd^{2}=0.09$ to the upper limit $t/E_{\text{F}}=1/9$ with $2mtd^{2}=1$ denoted by the horizontal (cyan) dotted line. In either panel, the (black dot-dashed) $\mu=0$ curve defines the boundary between the fermionic and the bosonic regimes. The (yellow) shaded region on the left of the (orange) dashed $a_{0}=0$ curve is a hole-like pairing regime with $a_{0}<0$, whereas the overall pairing evolves from hole-like into particle-like with $a_{0}>0$ across the $a_{0}=0$ curve. A PDW ground state with $t_{\text{B}}<0$ emerges within the grey shaded region, enclosed within the (green) $t_{\text{B}}=0$ curve. The entire phase space is a superfluid except for the PDW phase. Note that the PDW phase usually starts immediately before $\mu$ decreases down to zero, as the pairing strength increases. The fact that there are two branches of the $t_{\text{B}}=0$ curve indicates that there is an reentrant behavior of $T_{\text{c}}$ as a function of pairing strength. In the absence of population imbalance, similar reentrant behavior of superfluidity and associated PDW ground state have not been found in any other balanced systems with a short range pairing interaction, except in a very narrow range of density slightly above 0.53 in the attractive Hubbard model [46, 48, 49]. With a long-range anisotropic dipole-dipole interaction, however, such a reentrant behavior and PDW state have been predicted in the $p$-wave superfluid in dipolar Fermi gases [38].

As shown in Fig. 3, the interaction range for hole-like pairing extends toward stronger pairing regime with (a) increasing $d$ or (b) decreasing $t$. This can be explained by the evolution of the shape and topology of the Fermi surface, as shown in Fig. 2. As $d$ increases or $t$ decreases, the Fermi surface gradually opens up at the four $X$ or $Y$ points located at $(k_{x},k_{y})=(\pm\pi/d,0)$ and $(0,\pm\pi/d)$, and becomes fully open at the first BZ boundary for large $d$ small $t$, leading to an effective filling factor above $1/2$ in the lattice dimensions. In contrast to the 1DOL case, the existence of two lattice dimensions is enough to dominate the contributions of the remaining one continuum dimension (which is always particle-like due to its parabolic fermion dispersion), so that both the in-plane and the overall pairing becomes hole-like when $d$ is large or $t$ is small, with $a_{0}<0$ in the linear frequency term of the inverse $T$ matrix expansion. This is especially true in the weak coupling regime, where the superfluidity is more sensitive to the underlying Fermi surface. As the interaction goes stronger toward the BEC regime, the gap becomes large and the Fermi level (i.e., chemical potential $\mu$) decreases and then becomes negative, hence the shape of the non-interacting Fermi surface is no longer important. In this case, the contributions from the lattice dimensions will spread evenly across the entire BZ, so that the continuum dimension will become dominant, and the overall pairing eventually changes from hole-like to particle-like (with $a_{0}>0$). As shown in Fig. 2, within the occupied range of $k_{z}$, the average (or effective) filling factor within the first BZ in the $xy$ plane increases with increasing $d$ and/or decreasing $t$. Therefore, as $d$ increases, or $t$ decreases, the effect of the above-half-filling status persists into stronger pairing regime, and thus the hole-like pairing region in Fig. 3 extends toward right.

Shown in Fig. 4 is the behavior of (a) $\mu$ as a function of $U$, along with (b) $2n_{\text{p}}/n$, where $n_{\text{p}}\equiv a_{0}\Delta^{2}$, for $t/E_{\text{F}}=0.05$ and $k_{\text{F}}d=3$. Also plotted are $a_{0}$ and $\Delta$. This corresponds to a horizontal cut at $k_{\text{F}}d=3$ in Fig. 3(a) or at $t/E_{\text{F}}=0.05$ in Fig. 3(b). Inside the hole-like pairing regime, $a_{0}<0$ and thus the chemical potential $\mu$ goes above its noninteracting value. This can be seen from the expression [46, 36]

$$\sum_{\mathbf{k}}\Theta(-\xi_{\mathbf{k}})=n/2-a_{0}\Delta^{2}\,.$$

(8)

The chemical potential $\mu$ increases with the pairing strength, until it reaches a maximum where $n_{\text{p}}$ reaches a minimum. Here $2n_{\text{p}}/n$ is roughly the pair fraction, which reaches unity in the BEC regime. This plot is very close to its counterpart at $T_{\text{c}}$, which can be found in Ref. [36], since the temperature dependencies of both $\mu$ and $a_{0}$ are weak, except that here $a_{0}$ changes sign at a slightly larger $U/U_{\text{c}}$. As usual, the excitation gap $\Delta$ increases with $U/U_{\text{c}}$.

The PDW ground state in Fig. 3 with $t_{\text{B}}<0$ at a intermediate coupling strength for (a) relatively large $k_{\text{F}}d$ with fixed $t/E_{\text{F}}=0.05$ or (b) small $t$ with fixed $k_{\text{F}}d=3$ is associated with the strong inter-pair repulsive interaction, relatively large pair size and high pair density. Close to $\mu=0$, nearly all fermions have paired up with a relatively large pair size and a heavy effective pair mass, and the inter-pair repulsive interaction becomes strong. A large $d$ or small $t$ strongly suppresses the pairing hopping kinetic energy, and the large pair size and high pair density strongly reduce the pair mobility. All these factors lead to Wigner crystallization and hence PDW in the $xy$ plane, which can also be called a Cooper pair insulator. The negative sign of $t_{\text{B}}$ within the grey shaded region indicates that the minimum of the pair dispersion $\tilde{\Omega}_{\textbf{q}}$ has shifted from $\mathbf{q}=0$ to $\mathbf{q}=(\pi/d,\pi/d,0)$, with crystallization wave vector $(q_{x},q_{y})$ in the $xy$ plane. As the pairing interaction increases in the BEC regimes, the pair size shrinks and inter-pair repulsive interaction becomes weak; hence $t_{\text{B}}$ changes from negative back to positive, corresponding to a quantum phase transition from a PDW insulator to a superfluid.

Combining Figs. 2 and 3, we find that the emergence of hole-like pairing and the PDW phase is associated with the open Fermi surface topology. Once the Fermi surface is closed, both hole-like pairing and the PDW phase disappear.

In case of a closed Fermi surface, typical behaviors of the chemical potential $\mu$ and the excitation gap $\Delta$ for the balanced case can be seen from the $p=0$ lines in Fig. 5, calculated for $t/E_{\text{F}}=0.05$ and $k_{\text{F}}d=1$. Here $\mu$ decreases monotonically with $U/U_{\text{c}}$. Without a hole-like pairing regime, these solutions look qualitatively similar to other cases, e.g., in 3D continuum or 3D lattice, except that they follow a different asymptotic behavior in the BEC limit [36].

We now proceed and present our results for the population imbalanced case. With the added parameter $p$, the phase diagram becomes much more complicated. It renders the otherwise superfluid state unstable in the vast areas in the phase space.

To make the comparison easier, we begin by presenting phase diagrams in Fig. 6 in the same (a) $d$ – $U$ and $t$ – $U$ planes, as in Fig. 3, but with a tiny nonzero $p=0.001$. Here a normal gas phase (grey shaded) emerges in the weak coupling regime, delineated by the (black solid) $T_{\text{c}}^{\text{MF}}=0$ line, which is given by Eqs. (3)-(5) with $\Delta=0$. Indeed, in the presence of an imbalance, pairing cannot take place for an arbitrarily weak interaction. There exists a stable pSF phase (yellow shaded), defined by the (green solid) $t_{\text{B}}=0$ line and further confined by the stability condition (red solid line). The pSF phase resides in the low $d$ and large $t$ regime. A PDW ground state emerges in the dot shaded region, enclosed by the $t_{\text{B}}=0$ line and the dashed part of the (red) stability line. Then the rest unshaded space allows for an unstable mean-field superfluid solution, which may yield to phase separation. Now that the underlying lattice in the $xy$ plane breaks the continuous translational symmetry, the exotic Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states may possibly exist in part of the unstable region [50, *LO, 52, 53].

One can immediately tell that the vertical axes in Fig. 6 take different parameter ranges from those in Fig. 3, even though the imbalance $p=0.001$ is very small. While the $d$ – $U$ phase diagram in Fig. 6(a) is still calculated with $t/E_{\text{F}}=0.05$, the stable pSF phase is now restricted to relatively small $d$ (yellow shaded area). However, the $t$ – $U$ phase diagram has to be calculated at a much smaller $d$, with $k_{\text{F}}d=1.5$, as there is no stable pSF phase for $k_{\text{F}}d=3$ within the constraint $2mtd^{2}\leq 1$ (i.e., $t/E_{\text{F}}\leq 1/9$). In both cases in Fig. 6, the Fermi surface is closed. Unlike the balanced cases, one cannot find a stable superfluid solution with an open Fermi surface. For this reason, one does not find a hole-like pairing region in the weak coupling regime, but rather one in the strong coupling regime, on the right of the (blue dashed) $a_{0}=0$ line. Note that in the superfluid phase of hole-like pairing (on the right of the blue dashed line), both $a_{0}$ and $t_{\text{B}}$ are negative but the product $\xi_{xy}^{2}$ is positive. Outside the $t_{\text{B}}=0$ curve, we have $\xi_{xy}^{2}<0$, so that the mean-field superfluid solution becomes unstable, yielding to the PDW phase. The smallness of $p$ suggests that the ground state of $p\rightarrow 0^{+}$ is not continuously connected to the $p=0$ case, consistent with that in 3D continuum [19]. In comparison with Fig. 3, the current large PDW phase in the bosonic regime is totally a consequence of population imbalance.

Now we take $p$ as a varying parameter and explore phase diagrams in the $p$ – $U$ plane. Shown in Fig. 7 are the phase diagrams for (a) $(t/E_{\text{F}},k_{\text{F}}d)=(0.15,1)$, (b) (0.05, 1), and (c) (0.15, 1.5). Panels (b) and (c) show the effect of changing $t$ and $d$, respectively. In all three cases, there are three different phases, delineated by solid lines, as well as a PDW phase. A normal gas phase (grey shaded) takes the weaker coupling and larger $p$ area, on the left of the $T_{\text{c}}^{\text{MF}}=0$ curve. The vast majority is an unstable mean-field superfluid (unshaded), which should yield to phase separation or FFLO solutions. The stable pSF phase (yellow shaded) occupies only a small area. Finally, the PDW phase (dot shaded) takes the small region next to the pSF phase, bounded by the (red dashed) stability $\partial^{2}\Omega_{\text{S}}/\partial\Delta^{2}=0$ line and (green sold) $t_{\text{B}}=0$ line. When compared with panel (a), one readily sees that the pSF phase shrinks as $t$ decreases (panel (b)) and/or as $d$ increases (panel (c)). This is because both increasing $d$ and reducing $t$ lead to stronger momentum confinement in the lattice dimensions. In agreement with Fig. 6, the Fermi surface for all these three cases are closed. Note that the (red) stability line and the (green) $t_{\text{B}}=0$ line cross into each other, and the pSF phase is bounded by the stronger of these two conditions. Here also plotted are the lines along which the superfluid density vanishes. As found in 3D continuum, the positivity of superfluid density constitutes a much weaker stability constraint, as both lines of $(n_{\text{s}}/m)_{x}=0$ in the lattice dimension and of $(n_{\text{s}}/m)_{z}=0$ in the continuum dimension lie completely within the unstable area. Note that while the $(n_{\text{s}}/m)_{z}=0$ line looks very similar to its 3D continuum counterpart, the $(n_{\text{s}}/m)_{x}=0$ line exhibits an unusual nonmonotonic behavior, caused by the lattice effect. From the (violet dotted) $\mu=0$ curve, one readily sees that, as in Fig. 6, the pSF phase resides completely within the bosonic regime.

The fact that the pSF phase exists only in a small bosonic region (in both Fig. 6 and Fig. 7) is in stark contrast with the 3D continuum case, for which the stability line $\partial^{2}\Omega_{\text{S}}/\partial\Delta^{2}=0$ extends monotonically up to $p=1$, and a polarized superfluid exists for arbitrary imbalance $p$ in the BEC regime [19]. Apparently, this difference can be attributed to the presence of two lattice dimensions. Indeed, for 1DOL, with only one lattice dimension, the stability line already cannot extend to $p=1$. However, the pSF phase in 1DOL can extend all the way to the deep BEC limit [43]. This is also supported by the fact that with three lattice dimensions in a 3D attractive Hubbard model, one can barely find a pSF state except at very low density and extremely low $p$ [54]. Therefore, one can conclude that more lattice dimensions make it more difficult to have a stable pSF ground state.

This phenomena can be easily understood from the momentum distribution of paired fermions, which would be given by $v_{\mathbf{k}}^{2}$ had there been no imbalance. In 3D continuum, $v_{\mathbf{k}}^{2}$ in the deep BEC regime extends to the entire infinitely large momentum space in all directions, leading to a vanishingly small occupation for paired fermions. Therefore, the excessive majority fermions can readily occupy the low momentum states, with essentially no Pauli blocking from paired fermions. However, when one or more lattice dimensions are present, the momentum in these dimensions is restricted to the first BZ, so that $v_{\mathbf{k}}^{2}$ in these dimensions cannot be infinitesimally small even in the extreme BEC limit, which will cause a repulsion to excessive majority fermions. This repulsion increases with $p$, and may become costly enough so as to render the mean-field superfluid solution unstable. As a result, the distribution of paired fermions is now roughly given by that of the minority fermions, $n_{\mathbf{k}\downarrow}=\Theta(E_{\mathbf{k}\uparrow})v_{\mathbf{k}}^{2}$, which reduces to $v_{\mathbf{k}}^{2}$ for $p=0$.

Unlike the $p=0$ case, for which hole-like pairing takes place in the weaker coupling regime when $t$ is small and/or $d$ is large, here hole-like pairing occurs in the BEC regime via a completely different mechanism. As mentioned above, all the three cases shown in Fig. 7 have a closed noninteracting Fermi surface. As the pairing becomes stronger, the momentum distribution of $v_{\mathbf{k}}^{2}$ in the $xy$ plane extends to the entire first BZ, and becomes roughly a constant at strong coupling; in the absence of population imbalance, this would lead to a rough cancellation (via averaging over the inverse fermion band mass) due to the particle-hole symmetry of the lattice band. However, for any finite $p$, the excessive majority fermions will tend to occupy the low $(k_{x},k_{y})$ states, and thus expel paired fermions toward higher $(k_{x},k_{y})$ states, which have a negative (i.e., hole-like) band mass, leading to a net hole-like contribution to $a_{0}$ in the pair propagator, when integrated over the entire BZ. This also explains why the $a_{0}=0$ line leans toward weaker coupling with increasing $p$.

Shown in Fig. 8 is an example of the momentum distributions of $v_{\mathbf{k}}^{2}$ (left), $n_{\mathbf{k}\downarrow}$ (middle) and $\delta n_{\mathbf{k}}$ (right column) in the $(k_{x},k_{y})$ plane at different $k_{z}/k_{\text{F}}=0$ (top), 0.2 (middle) and 0.4 (bottom row), with $U/U_{\text{c}}=4$ and $p=0.05$, for $t/E_{\text{F}}=0.15$ and $k_{\text{F}}d=1.5$. This corresponds to a PDW state in Fig. 7(c). Indeed, the excessive fermion distribution, $\delta n_{\mathbf{k}}=\Theta(-E_{\textbf{k}\uparrow})$, occupies the low in-plane momentum part and below $k_{z}/k_{\text{F}}=0.4$ (right column). In addition, $v_{\mathbf{k}}^{2}$ (left column) remains roughly constant in the entire BZ and for $|k_{z}/k_{\text{F}}|\leq 0.4$. Most interestingly, the minority fermion distribution $n_{\mathbf{k}\downarrow}$ (middle column) is given by $v_{\mathbf{k}}^{2}$ but with a hole dug out at the center, due to the Pauli repulsion with the excessive fermions.

As a representative example, we show in Fig. 5 the behavior of (a) $\mu_{\sigma}$ and (b) the gap $\Delta$ for $p=0.05$ (red) and 0.1 (blue) with fixed $t/E_{\text{F}}=0.05$ and $k_{\text{F}}d=1$, as a function of $U$. They correspond to horizontal cuts at $p=0.05$ and 0.1 in Fig. 7(b), and should be compared with the $p=0$ case (black solid curves). The solid part of these lines are stable pSF solutions, while the dashed lines are unstable mean-field solutions. There are a few remarkable features. Firstly, the excitation gap changes only slowly with imbalance $p$, except that it does not have a solution below certain threshold of interaction strength. Secondly, at given pairing strength, $\mu_{\sigma}$ for $p=0.05$ and $p=0.1$ are very close to each other, but both far separated from the $\mu$ curve for $p=0$. This again indicates that the $p\rightarrow 0^{+}$ case is not continuously connected to the $p=0$ case; with a tiny bit of imbalance, $\mu_{\uparrow}$ and $\mu_{\downarrow}$ immediately split up. Lastly, $\mu_{\uparrow}$ increases slowly with pairing strength in the BEC regime. This is different from its counterpart in 3D continuum and 1DOL; for the former, $\mu_{\uparrow}$ decreases while for the latter $\mu_{\uparrow}$ approaches a $p$-dependent constant asymptote, as the pairing strength increases toward the BEC limit. This can be attributed to the emergence of hole-like pairing (with $a_{0}<0$) in the strong pairing regime as the number of lattice dimensions increases. To verify this idea, we have also checked the mean-field solution for imbalanced 3DOL, and found that, indeed, $\mu_{\uparrow}$ also increases with the pairing strength in the BEC regime at $T=0$, along with a negative $a_{0}$.

Finally, we present the typical behavior of the superfluid density in the imbalanced case. Shown in Fig. 9 are (a) $(n_{\text{s}}/m)_{z}$ and (b) $(n_{\text{s}}/m)_{x}$ in the continuum and lattice dimensions, respectively, as a function of $U/U_{\text{c}}$ for $p=0$, 0.05 and 0.1 at fixed $t/E_{\text{F}}=0.05$ and $k_{\text{F}}d=1$. Here solid and dashed lines are stable and unstable solutions, respectively. As expected, both are always positive for the balanced case. In addition, $(n_{\text{s}}/m)_{x}$ is much smaller than $(n_{\text{s}}/m)_{z}$, because it involves the average of the inverse band mass. For the imbalanced case, the superfluid density deviates continuously from its positive $p=0$ value as $p$ increases from 0. However, in the unitary and weak coupling regimes, both continuum and lattice components will become negative for $p\neq 0$. Furthermore, the superfluid density is more negative for smaller (but finite) $p$. This implies an immediate discontinuous jump from the $p=0$ value to a large negative value for $p=0^{+}$ in this regime. Note that for strong enough interaction, $(n_{\text{s}}/m)_{x}$ will again change sign to negative, but gradually rather than abruptly, as can already be seen from the $p=0.1$ curve. This has to do with the lattice induced confinement in the momentum space and the Pauli exclusion between paired and excessive fermions.

So far, it is not yet clear whether the PDW state can sustain a superfluid order, with and without an imbalance. If the answer is yes, then it will become a supersolid state rather than a Cooper pair insulator. We leave this to a future study.

It should be noted that we have worked with a system with homogeneous fixed densities. For this reason, we have not chosen to use $\mu$ and $h$ as control variables, which are more appropriate for systems connected with a large reservoir so that the chemical potentials are fixed or can be tuned separately. In such a case, all $h<\sqrt{\min(0,\mu)^{2}+\Delta^{2}}$ corresponds to the population balanced state. One can, however, convert between these two approaches, by calculating corresponding densities (and Fermi energy) for given $\mu$ and $h$, and performing a rescaling.