FFLO order in ultra-cold atoms in three-dimensional optical lattices

In the past two decades remarkable progress in cold atom physics has opened a new frontier in the construction and precise control of quantum systems. Following the development of a number of important experimental techniques, including Feshbach resonances and optical lattices, it was quickly suggested that ultra-cold atomic gases provide an ideal setting for the realization and investigation of a variety of exotic physical phenomena [1]. These systems provide experimental analogues to many condensed matter systems, but are also highly tunable and free of disorder. These experiments represent an exciting opportunity to simulate the fundamental mechanisms and models of condensed matter physics, for instance Cooper pairing of fermions and the Hubbard model, without the additional complexities presented by real materials. A number of experiments have already demonstrated the possibilities for ultra-cold atomic gases, including inducing superfluidity in fermionic systems and probing the BEC-BCS crossover [2–5].

In light of these advances, one system that has attracted considerable interest is an ultra-cold atomic gas in an optical lattice with unequal populations of two hyperfine states. The hyperfine states can be seen as two distinct spin species, and an attractive interaction can be induced between them, with its strength tunable, using a Feshbach resonance. This system represents an experimental simulation of the attractive fermionic Hubbard model. It was first suggested by Fulde and Ferrell (FF) [6], and separately by Larkin and Ovchinnikov (LO) [7], that the mismatched Fermi surfaces in a polarized system of this type could result in an instability to the formation of a condensate of finite-momentum pairs. However, the FFLO phase has eluded conclusive detection for nearly fifty years. Considering how challenging the observation of this phase has proven to be, reliable determination of the parameter domain in which this phase might exist, and its properties, remains an important goal.

Many efforts have been made, using a variety of theoretical and numerical techniques, to achieve this goal and to characterize the properties of the FFLO phase. However, in most cases these studies were limited to targeted states, fixed size simulation cells or to one- and two-dimensional lattices [8–14]. Three-dimensional lattices are in many ways the most direct and natural for optical lattice experiments with ultra-cold atomic gases, so these systems offer the most realistic possibility of observing FFLO states. With this in mind, we map the density-polarization phase diagram for spin-imbalanced fermions with attractive interactions in a 3D optical lattice in the present study.

While 3D systems may present great opportunities to observe the FFLO state experimentally, they present a considerable computational challenge. We carry out detailed calculations using the Hartree–Fock–Bogoliubov theory, which is the simplest quantitative approach. At the minimum, results from these mean-field calculations provide a qualitative understanding of the nature of the phases in a large region of the parameter space, and propose candidate phases for more elaborate (and computationally intensive) many-body approaches. In fact, comparisons [15] with quantum Monte Carlo (QMC) results [16] indicate that mean-field results provide not only qualitative but quantitatively useful information in the doped repulsive Hubbard model, which is closely related to the present model. Thus, a second goal of our study is to obtain numerically robust mean-field results to provide a foundation for QMC calculations on the same system.

Despite the simple nature of the mean-field approach, the determination of the correct ground state in the 3D lattice is far from straightforward [17]. To determine the stability of states that have 3D spatial dependence of the order parameter requires the use of cubic simulation cells, which quickly become computationally expensive as the system size increases. Additionally, 3D systems permit a wider range of potential ground-states, meaning the energy landscape will have more local minima and ground-state searches need to be increasingly thorough. We focus on moderate interaction strengths (U/t ⩽ 5), where this approach is most reliable. Several strategies are employed, using large scale computations, to validate the solutions and the extrapolation to the thermodynamic limit.

We find that, at high to intermediate densities, the ground state phase is of the canonical LO form independently of interaction strength, with counter-propagating pairs and order parameter going to zero on a regularly spaced array of parallel planes. This is the domain in which the effect of the optical lattice is most apparent on the shape of the Fermi surfaces, and consequently on the ground state phases. At low density, the Fermi surfaces become more spherical, as they would be in the continuum, and we find that the ground state is characterized by a superposition of pairs with non parallel momenta. In this region, where the impact of the optical lattice is less significant and these higher-dimensional states emerge, a larger interaction is required to induce pair ordering. Systematic information is obtained on the ground-state properties, especially in the first parameter regime. The physical origin of the phases and their connection to the Fermi surface topology and pairing are discussed.

Below we first describe our computational approach in section 2. In section 3 the results for the first parameter regime, namely at high to intermediate densities, are presented, with discussions of the effects of density and polarization, and of the interaction strength. Results more relevant to the continuum limit, i.e. at low densities, are then discussed in section 4. We conclude with a summary in section 5.

The starting Hamiltonian we study is,

$$\begin{eqnarray} &&\fl H= - \sum_{(ij)\sigma} t_{ij} c_{i\sigma}^\dagger c_{j\sigma} - \sum_i \left( Un_{i\uparrow} n_{i\downarrow} + \mu n_i + \frac{h}{2} m_i\right), \end{eqnarray} \tag{ 1 }$$

where c_iσ is a fermionic annihilation operator of spin σ on site i, $n_{i\sigma}=c^{\dagger}_{i\sigma} c_{i\sigma}$, n_i = n_i↑+n_i↓ and m_i = n_i↑ − n_i↓. In this paper we will only consider the Hubbard dispersion, i.e. t_ij = t if (ij) = 〈ij〉(i and j are near-neighbors) and t_ij = 0 otherwise. The interaction will be attractive, so U > 0. Further, we will be in the regime of negative scattering length, since we will be concerned with U/t ⩽ 5, as mentioned earlier. (A two-body bound state first appears at U/t = 7.913 55 for the Hubbard dispersion. See, e.g. [18]) The chemical potential μ and the 'magnetic field' h in the Hamiltonian control the density, n, and the polarization, p. Given a supercell of N lattice sites, these are defined by n_σ ≡ ∑_i〈n_iσ〉/N: n = n_↑ + n_↓, m = n_↑ − n_↓, and p ≡ m/n. The system is completely specified by the three parameters U/t, n, and p.

Our analysis of this Hamiltonian was performed using Hartree–Fock–Bogoliubov theory. We transform the Hamiltonian into a diagonalizable form by employing a standard mean-field approximation,

$$\begin{eqnarray} &&\fl \sum_i U n_{i\uparrow}n_{i\downarrow} = \sum_i Uc^{\dagger}_{i\uparrow} c_{i\uparrow} c^{\dagger}_{i\downarrow} c_{i\downarrow}\nonumber\\ &&\fl \rightarrow \sum_i U\left\{\langle c^{\dagger}_{i\uparrow} c^{\dagger}_{i\downarrow}\rangle c_{i\downarrow} c_{i\uparrow} + \langle c_{i\downarrow} c_{i\uparrow}\rangle c^{\dagger}_{i\uparrow} c^{\dagger}_{i\downarrow}\right.\nonumber\\ &&\qquad \left. + \langle c^{\dagger}_{i\uparrow} c_{i\uparrow}\rangle c^{\dagger}_{i\downarrow} c_{i\downarrow} + \langle c^{\dagger}_{i\downarrow} c_{i\downarrow}\rangle c^{\dagger}_{i\uparrow} c_{i\uparrow} \right\}, \end{eqnarray} \tag{ 2 }$$

with constant terms omitted from the expression above for brevity.

The FFLO phase is most distinctly characterized by a spatially modulated pairing order parameter. In order to accurately determine the relative stability of FFLO states with different real-space structures, we perform our calculations on simulation cells whose shapes accommodate those structures. The simulation cells are characterized by three basis vectors, L₁, L₂ and L₃, whose components are integers. Once the cell shape is chosen, we introduce Bloch states, defined as c_j(k) ∝ ∑_Lc_j+Lexp[ik · L] where L is a vector on the Bravais lattice having L₁, L₂ and L₃ as basis vectors, i.e. L = l₁L₁ + l₂L₂ + l₃L₃, and k is a vector that varies freely within the first Brillouin zone of the simulation cell reciprocal lattice. We use twist-averaged boundary conditions in all of our calculations.

Having applied the mean-field approximation, we can use the Bloch states described above to write the Hamiltonian as a sum of decoupled k-dependent pieces, H = ∑_kH(k), of the form

$$\begin{eqnarray} \fl H({\bf k})= [ {\bf c}^{\dagger}_\uparrow {\bf c}_\downarrow] \left[ \begin{array}{cc} {\bf H}_\uparrow({\bf k}) & \boldsymbol{\Delta} \\ \boldsymbol{\Delta}^\dagger & - {\bf H}_\downarrow^T({\bf G}-{\bf k}) \end{array} \right] [ {\bf c}_\uparrow {\bf c}^{\dagger}_\downarrow]^T \end{eqnarray} \tag{ 3 }$$

where c_↑ and c_↓ represent an array (row) of operators, {c_i↑(k)} and {c_i↓(G − k)} with the index i running over the N sites of the simulation cell. The vector G is defined so that θ = G · L is the twist angle of the pairing order parameter after a translation by L. H and Δ are N × N matrices with elements

$$\begin{eqnarray} [ {\bf H}_\sigma({\bf k})]_{ij} & = -t_{ij}({\bf k}) + \delta_{ij} ( D_{i\sigma} - \mu -s_\sigma h/2) \\ {[\boldsymbol{\Delta}]}_{ij} & = \delta_{ij} \Delta_i .\nonumber \end{eqnarray} \tag{ 4 }$$

In the above equation t_ij(k) = ∑_L exp(ik · L) t_i,j+L, s_↑/↓ = ±1 and D_iσ, Δ_i, μ and h are determined by the requirement that the Free energy F = 〈H〉 − TS is a minimum for the target average densities n_σ. (All of our calculations are performed at T/t = 0.01.) This amounts to the following self-consistency equations

$$\begin{eqnarray} D_{i\sigma} &= -U \int {\rm d}{\bf k} \langle c^{\dagger}_{i\sigma'}({\bf k})c_{i\sigma'}({\bf k}) \rangle \nonumber\\ \Delta_i &= -U \int {\rm d}{\bf k} \langle c_{i\,\downarrow}({\bf k})c_{i\,\uparrow}({\bf k}) \rangle\\ n_\sigma &= N^{-1} \sum_i \int {\rm d}{\bf k} \langle c^{\dagger}_{i\sigma}({\bf k})c_{i\sigma}({\bf k})\rangle.\nonumber \end{eqnarray} \tag{ 5 }$$

where in the first equation σ' is the opposite of σ.

We make the following initial ansatz for the spatial form of the order parameter,

$$\begin{equation} \Delta_i^{(0)} = \sum_ {\bf q} \Delta^{(0)}_ {\bf q} {\rm e}^{{\rm i} {\bf q}\cdot {\bf r}_i}. \end{equation} \tag{ 6 }$$

This represents a summation of plane wave modes characterized by a set of symmetry-related pairing vectors q. The spiral (FF) phase corresponds to a single $\Delta^{(0)}_ {\bf q}\ne0$ or, in real space, to $\Delta^{(0)}_i\propto {\rm e}^{{\rm i} {\bf q}\cdot {\bf r}_i}$. The linear (LO) phase has $\Delta^{(0)}_{\pm {\bf q}}\ne0$ with q ∝ (0, 0, 1), (0, 1, 1), or (1, 1, 1) and $\Delta^{(0)}_i\propto \cos( {\bf q}\cdot {\bf r}_i)$. In addition, we consider 2D structures of the form Δ_i ∝ cos (q_y · r_i) + cos (q_z · r_i), with q_y = |q|(0, 1, 0) and q_z = |q|(0, 0, 1), and 3D structures of the form Δ_i ∝ cos (q_x · r_i) + cos (q_y · r_i) + cos (q_z · r_i), with q_x = |q|(1, 0, 0) and q_y, q_z as before.

Our procedure allows us an unbiased search of the ground state within the general form of the candidate orders which are tested. Different choices of $\Delta_i^{(0)}$ determine different shapes of the simulation cell which, in turn, constrain the form of the self-consistent Δ_i. A typical example, for a linear phase, might have L₁ = (1, 0, 0), L₂ = (0, 1, 0), and L₃ = (0, 0, 50). After the shape of the simulation cell has been selected, we perform a scan over |q| to determine the optimal |q| corresponding to the minimum energy ground state for each q-direction (or for the higher-dimensional structures, the minimum energy for each set of q's). For each calculation in the scan, we sum over a sufficiently dense k-grid to remove all finite-size effect except for the constraint on the form of the order from the shape of the simulation cell. In the case above, for example, our calculation would use a k-point grid which has dimensions of a few in the L₃ direction and a few hundred in the L₁ and L₂ directions. This technique allows the calculation to accommodate the spatial modulation of the phase and approach the thermodynamic limit.

This procedure is sketched schematically for linear phases in figure 1. The calculations are to determine the true ground state among pair-ordered states with pairing vector q directed along either (1, 0, 0), (1, 1, 0), or (1, 1, 1). For each q-direction we perform a scan to determine the optimal |q|, varying the simulation cell size to ensure that it is commensurate with the targeted value of |q|. To rule out orders other than linear, we carry out searches for the 2D and 3D structures described above. Further, we increase the simulation cell size in directions other than q to verify the stability of the solution.

Zoom In Zoom Out Reset image size

We first consider the region of high to intermediate densities and low polarizations, where the characteristics of the ground-state phases of the system are significantly impacted by the presence of the optical lattice. This effect is most clearly reflected in the shape of the Fermi surface. At high density the Fermi surfaces of both spin species are very distinct from their spherical counterparts in the continuum. The nature of the pairing mechanism and its connection to the shape of the Fermi surfaces is further discussed below.

As described in section 2, the set of pairing wave vectors that leads to the minimum energy state determines the spatial structure of the pairing order parameter of that state. We find that in the optical lattice regime the system favors states with two q vectors, which results in an order parameter that is a linear pairing wave. The spiral state is energetically less favorable and never found to be the ground state in the regime we have investigated. This is similar to the situation in 2D [10] and is consistent with the results from the 3D repulsive Hubbard model [17] after particle-hole mapping. The properties of the linear phases, including the direction of the q vectors, exhibit dependence on density and polarization, and will be discussed in detail in section 3.1.

In figure 2, we present a characteristic example of the linear LO phase, in order to illustrate some of its real space properties. The ground state at these parameters is found to have q ∈ {|q|(0, 0, 1), |q|(0, 0, −1)}. At small polarizations and high densities such as this particular case, the domain wall nature of the pairing wave is evident. The densities of both spin species exhibit spatial modulation, with the density of the majority equal to the density of the minority at the peak of the order parameter. The greatest difference between the minority and majority density occurs at the nodes of the order parameter. This results in a peak of the spin density, which can be understood as the localization of excess spin at the nodes of the order parameter. The quantity α ≡ m π/|q| characterizes the total density of the excess spin within each nodal region (a stack of planes perpendicular to q. The overall charge density of the system shows no significant modulation in this case.

Zoom In Zoom Out Reset image size

The momentum-space properties of the same state are plotted in figure 3 using the gradient of the momentum distribution. This quantity gives the position of the underlying Fermi surface which, as shown later (figure 5), need not coincide with the non-interacting one. Illustrated on the plot is the pairing construction, k → −k + q, by which electrons (particles) near the Fermi surfaces of the two different spin species form pairs with finite momentum q. In this case, a slight modification of the shape of the interacting Fermi surfaces from the non-interacting ones allows ↑- and ↓- electrons along large sections of both Fermi surfaces to form pairs with a single pair of q's with common magnitude |q|. The resulting order parameter is a sum of plane waves, whose collective interference serves to lower the energy of the state and produce the standing wave structure visualized in figure 2. For the set of parameters corresponding to the state in the figure, and the slice of momentum-space plotted, a large fraction of the Fermi surface is smeared as a consequence of pair formation. The sharp features at the bottom of the minority Fermi surface identify a region where the Fermi surface is still intact and remains un-gapped. This is consistent with α ≠ 1 and a metallic nodal region [10]. In this case, the intact portion of the minority Fermi surface is small, indicating that most of the electrons near the Fermi surface have paired.

Zoom In Zoom Out Reset image size

Having highlighted the important features of the FFLO phase in the optical lattice regime, in both real and momentum space, we will now discuss in more detail the effect of density, polarization, and interaction strength on these features. A final phase diagram summarizing all our calculations is then presented in section 3.2.

3.1. Density and polarization

In this section we examine in further detail the characteristics of the ground-state phases as they depend on density and polarization. At each selected interaction strength U/t, we map out the complete n–p phase diagram. The behavior of the linear phase as a function of polarization, for n = 0.96, n = 0.60 and n = 0.24 at U/t = 4 is illustrated in figure 4. At large polarizations, near the onset of pairing order, the order parameter is small, large portions of the Fermi surfaces of the two spin species remain ungapped, and those that are gapped remain sufficiently sharp to be precisely located. As the polarization decreases, pairing is enhanced and the pairing order parameter increases as expected. Lower polarization is also where it is more likely to have (1, 1, 1) order, and a transition to it from (0, 0, 1) can be seen in the figure where the value of α increases significantly (going from high to low polarization), for n = 0.96 and n = 0.60. The appearance of (1, 1, 1) order involves larger Fermi surface reconstructions, in a way similar to the nesting mechanism for the formation of spin-density waves in the 3D repulsive case [17].

Zoom In Zoom Out Reset image size

Figures 5 and 6 visualize and compare the momentum- and real-space properties, respectively, for different values of the polarization. As already discussed, the underlying Fermi surface of the LO phase can deviate from the non-interacting one. The numerical solution can be understood by the momentum space nesting caused by the surface reconstruction and the pairing mechanism that ensues. At large polarizations, a larger q is required, and smaller portions of the Fermi surface can support pairing, hence weaker order parameter. Eventually, as one moves farther from the transition and deeper into the LO phase, the Fermi surface is heavily smeared, the order parameter comprises many (collinear) momenta. Correspondingly, in real space the order parameter remains purely sinusoidal, the density modulation is weak, and the excess spin is not localized at large polarization. As the polarization decreases, the physics is better understood in the language of weakly interacting domain walls, with the excess spin more localized at the nodes of the order parameter, and strong density modulation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 also captures the behavior of the ground state properties as a function of density. At high densities the presence of the underlying lattice has a significant effect on the shape of the Fermi surface. For states at high density |q| is large compared to states at the same polarization but lower density. Additionally, the effect of polarization on |q| is more prominent at higher density, where a larger spin imbalance is required to achieve the same polarization than is required at a lower density. This effect can be seen, for example, by comparing the slopes of |q|/π versus p for n = 0.96 and n = 0.24. The slope of the n = 0.96 curve is significantly steeper than the slope of the n = 0.24 curve. Also, all values of |q| are smaller for n = 0.24 compared to n = 0.60 and n = 0.96, which reflects the smaller mismatch between Fermi surfaces at lower density.

3.2. Interaction strength

In figure 7 we summarize the phase diagrams for three values of interaction strength. The interaction strength plays a significant role in determining the stability of pair ordered ground-state phases. The LO ground state becomes more stable as the polarization decreases and the interaction strength increases. This behavior is evident in the phase diagrams for U/t = 3, U/t = 4 and U/t = 5, which show that the area of phase space occupied by an ordered state grows larger with increasing interaction strength. The trend suggests that as the Fermi surfaces of the two spin species become closer and more similar in shape, pairing order becomes increasingly energetically favorable. This is especially true at higher interaction strengths, which allow a more significant reshaping of the Fermi surface to improve nesting and permit a larger number of electrons to participate in pairing.

Zoom In Zoom Out Reset image size

In each of the phase diagrams in figure 7, we show estimates of the parameter regions in which a local minimum exists for a pair-ordered state with pairing vector q directed along either (1, 0, 0), (1, 1, 0), or (1, 1, 1). These estimates are obtained by calculating U/t from the gap equation for Δ = 0 at fixed n and p with q ∝ (1, 0, 0), (1, 1, 0), or (1, 1, 1). For each q-direction we perform a scan in |q| to determine the minimum U/t required to induce pairing at the chosen n and p. We repeat this procedure for several hundred sets of n and p, which provides a map of the critical U across the phase-space. For a given U and q-direction this defines a curve in n and p outside of which the system will not have a pair-ordered solution with a pairing vector in the given q-direction. These curves are indicated for the different q-directions on the phase diagrams. They help guide our survey of the density-polarization phase space by indicating which states (defined by the direction of q) to consider in the fully self-consistent calculations. We then perform the numerical procedure outlined in section 2, and sketched in figure 1, which determines the true ground-state from the stable pair-ordered states. It is the full numerical search, the results of which are represented by the symbols in the phase diagrams, that provides the actual form of the order at each point.

In addition to affecting the overall stability of pair-ordered states relative to uniform states, the interaction strength also affects the density and polarization dependence of the transitions between the ordered phases, which are characterized by different sets of q vectors. At U/t = 3, we find that linear pairing order with the pairing-wave vector q directed along the (0, 0, 1)-direction is the ground-state for all values of density and polarization. We found no region of the U/t = 3 phase diagram in which the commensurate phase, defined by a density of one excess particle per node of the order parameter, is stable. This is seen in the U/t = 3 phase diagram, where no symbol reaches the color for α = 1. Instead, at low polarization and near half-filling α approaches 2/3. This behavior is caused by the nature of the LO ground-state at U/t = 3, which has q directed along the (0, 0, 1)-direction with |q| ≠ mπ. We observe that the commensurate phase has q ∝ (1, 1, 1), which does occur for U/t = 4 and U/t = 5.

At U/t = 4, a transition occurs between the linear phases with q ∝ (0, 0, 1) and q ∝ (1, 1, 1). The diagonal phase (q ∝ (1, 1, 1)) occupies the high to intermediate density and low polarization region of the phase space. In a portion of this region the commensurate phase is stable. At intermediate to high polarization, or for sufficiently low density, the pairing wave is directed along q ∝ (0, 0, 1) and the state is no longer commensurate.

The behavior at U/t = 5 is similar to that at U/t = 4, but with a larger region of stability for the diagonal phase. Again, in a portion of this region the commensurate phase is stable. As with U/t = 4, at large polarizations or low densities, the pairing wave is directed along (0, 0, 1), occupying a large portion of the phase space. The (0, 1, 1)-order is predicted by the gap equation to be stable in a large region but is never the true ground state.

The effect of increasing interaction strength is also apparent in the real-space character of the phases. This effect is visualized in figure 8. As interaction strength increases the pairing wave develops domain walls and the amplitude of the pairing wave and the density modulations grow. The larger density modulations cause the peaks of the spin density to become sharper, making the excess spin more localized.

Zoom In Zoom Out Reset image size

At low density the effect of the lattice on the shape of the Fermi surface is less significant and the properties of the system begin to resemble those of fermions in the continuum. In order to describe the experimental situation of trapped atomic gases, the Hamiltonian we have been using can be thought of as a discretized representation of the continuum [18]. The calculations must then be at the extremely dilute limit, with large supercells, to obtain realistic results in the thermodynamic limit in this situation. This is not the focus of the present study. However, we do extend our optical lattice studies above to selected lower densities. The results shed light on the approach to the continuum limit, which we discuss briefly here. A variety of work exists on mean-field and Hartree–Fock–Bogoliubov calculations in the continuum, see for example [19, 20] and references therein.

In this region, at large polarizations, we find that phases with a larger set of q's become energetically favorable relative to linear phases, which have just a single pair of q's. An example is illustrated in figure 9, which plots slices of the Fermi surfaces and spectral functions, and sketches the pairing construction, for a state with 2D modulations of the order parameter. The system forms pairs with q = ±|q|(0, 0, 1) and q = ±|q|(0, 1, 0), as compared to the case of linear order where pairs can form only with q = ±|q|(0, 0, 1). The additional q's allow for more pairing, again at little cost in kinetic energy, which lowers the total energy of the state.

Zoom In Zoom Out Reset image size

As depicted in the right panel, favorable nesting is achieved with four pairing wave vectors, which allows nearly every section of the majority spin surface to be covered by the minority surface. The sections that are not covered remain as bright spots, because the electrons in those regions have not paired and the Fermi surface remains intact.

In the dilute Fermi gas limit, the Fermi surfaces will be spherical and will not retain the features in the example above which made a 3D structure more favorable. However, more wave vectors can be involved which can create a more complicated structure of modulation to lower the interaction energy. This situation is seen in the electron gas, in which complex structures of spin-density waves are the true ground state in Hartree–Fock theory [21, 22]. Here we show one example, in figure 10, of the emergence of phases with higher-dimensional spatial variation of the order parameter as the lowest energy ground states of the system. At n = 0.24, p = 0.5833 the linear solution has the lowest energy. However, moving to lower density and polarization, but still near the onset of pairing order, the states with 2D and 3D spatial variation of the order parameter begin to have lower energy than the linear state. Finally, at n = 0.18, p = 0.444 the state with 3D structure emerges as the lowest energy ground state.

Zoom In Zoom Out Reset image size

We were able to perform calculations on moderately sized 3D simulation cells, up to 15³ sites, using GPUs to dramatically speed up the diagonalization. Even with these speed-ups, our search was somewhat limited by the rapidly increasing computational cost. To identify a genuine ground state with 3D structure, care was taken to ensure that the energy difference between the states with 2D and 3D structures was larger than any potential finite-size effect. For the case depicted in the rightmost panel of figure 10, this energy difference was ${\mathcal O}(10^{-6})$, whereas the convergence of the energy of both the 2D and 3D states to the thermodynamic limit was ${\mathcal O}(10^{-7})$, and the energy tolerance on the self-consistency loop was also ${\mathcal O}(10^{-7})$. The convergence to the thermodynamic limit was determined by comparing the energies from calculations using 100 k-points in each direction to calculations using 200 k-points in each direction for both states with 2D structures and states with 3D structures.

This result demonstrates the existence of a ground state with LO order of a 3D structure. The overall trend suggested by our results is that higher-dimensional ground states become increasingly stable, at relatively large polarization, with decreasing density near the onset of pairing order, and for n ≲ 0.18 the lowest energy ground state is likely to have a pairing order parameter with 3D spatial variation.

The 3D structure we observe corresponds to an order parameter that is the sum of six plane waves, as described in section 2. The set of q vectors is |q|{(±1, 0, 0), (0, ±1, 0), (0, 0, ±1)}. It has been suggested [23] that in the Fermi gas regime the most energetically favorable structure is a sum of eight plane waves of the (1, 1, 1) variety. As discussed in section 2 and indicated by figure 7, we expect solutions with q ∝ (1, 1, 1) to be stable only at small polarizations. However smaller polarizations result in a smaller |q|, as pictured in figure 4, and a smaller |q| corresponds to a longer wavelength pairing wave. This would require even larger 3D simulation cells, and thus lies outside the parameter region in which we have explored possible 3D structures.

We have carried out a systematic study of the phase diagram of spin-imbalanced fermions with attractive interactions in a 3D lattice. The phase space can be divided into two qualitatively distinct regimes, the optical lattice regime at high density and the Fermi gas regime at low density. In the optical lattice regime our survey involves detailed, fully self-consistent HFB calculations in which great care is taken to reach the true ground state at thermodynamic limit. The phase diagram in this regime was determined for up to intermediate interaction strengths. We find that the system favors linear pairing order of the LO type. At U/t = 3 the pairing vector q is directed along (0, 0, 1), and at U/t = 4 and U/t = 5 there is a transition from states with q along (0, 0, 1) at low polarizations to q along (1, 1, 1) at intermediate to high polarizations. The real and momentum space properties of these phases are determined. At low polarizations and high to intermediate densities the pairing wave is characterized by the presence of domain walls that become sharper with increasing interaction strength, and the localization of excess spin at the nodes. With increasing polarization and decreasing density the pairing wave becomes more sinusoidal and the excess spin is less strongly localized. Additionally, pairing becomes more stable with increasing interaction strength, as evidenced by the growing region of phase space occupied by pair-ordered phases. As we mentioned, the results are in the framework of mean-field theory. In the parameter regime studied here, they are likely to provide useful guides, as experience in the related repulsive doped Hubbard model has indicated. Future many-body calculations, for example QMC calculations, will build on these and be able to examine the effect of the mean-field approximation.

In the Fermi gas regime we searched for evidence of states with two and three-dimensional spatial modulation of the order parameter. With the use of GPUs to speed up the computation, we performed calculations on simulation cells large enough to accommodate both 2D and 3D structures. Our results provide evidence of the emergence of higher-dimensional states, which are most stable for low densities and high polarizations, near the onset of pairing order. These states occur as it becomes energetically favorable for the system to form pairs with a larger set of pairing vectors. Though our search was limited by the computational costs of large cubic simulation cells, our results suggest that for densities below n ≲ 0.18 the system supports FFLO states with 2D and 3D spatial variation of the order parameter, which makes this an interesting region for future theoretical and experimental exploration.

We acknowledge support from DOE (Grant no. DE-SC0008627) and NSF (Grant no. DMR-1409510). Computational support was provided by DOE leadership computing through an INCITE grant, and by the William and Mary SciClone cluster. We thank E Walter for help with computing.