Visualizing the BEC-BCS crossover in the two-dimensional Fermi gas: pairing gaps and dynamical response functions from ab initio computations

TL;DR

This paper presents ab initio Monte Carlo calculations of the pairing gap, spectral function, and structure factors in a 2D Fermi gas, offering detailed insights into the BEC-BCS crossover relevant for experiments.

Contribution

The study introduces first-principles numerical methods to compute dynamical properties of the 2D Fermi gas, enabling direct experimental comparison.

Findings

Computed pairing gap across the BEC-BCS crossover

Obtained spectral functions using analytic continuation

Visualized density and spin structure factors

Abstract

Experiments with ultracold atoms provide a highly controllable laboratory setting with many unique opportunities for precision exploration of quantum many-body phenomena. The nature of such systems, with strong interaction and quantum entanglement, makes reliable theoretical calculations challenging. Especially difficult are excitation and dynamical properties, which are often the most directly relevant to experiment. We carry out exact numerical calculations, by Monte Carlo sampling of imaginary-time propagation of Slater determinants, to compute the pairing gap in the two-dimensional Fermi gas from first principles. Applying state-of-art analytic continuation techniques, we obtain the spectral function, and the density and spin structure factors providing unique tools to visualize the BEC-BCS crossover. These quantities will allow for a direct comparison with experiments.

Figures

Equations

$$\hat{H}=t\sum _{\vec{k},\sigma }\epsilon (\vec{k}){\hat{c}}_{\vec{k},\sigma }^{†}{\hat{c}}_{\vec{k},\sigma }+U\sum _i{\hat{n}}_{i,\uparrow }{\hat{n}}_{i,\downarrow },$$

$$G^p(\vec{k},\tau )$$

$$G^h(\vec{k},\tau )$$

$$A(\vec{k},\omega )=⟨{\hat{c}}_{\vec{k}}\delta (\omega -\hat{H}){\hat{c}}_{\vec{k}}^{†}⟩+⟨{\hat{c}}_{\vec{k}}^{†}\delta (\omega -\hat{H}){\hat{c}}_{\vec{k}}⟩$$

$$S^{\hat{O}}(\vec{k},\omega )=⟨{\hat{O}}_{\vec{k}}\delta (\omega -\hat{H}){\hat{O}}_{-\vec{k}}⟩,$$

Full text

It is truly unusual when, starting from a microscopic Hamiltonian, theory can achieve an exact description of a strongly correlated fermionic system which, at the same time, can be realized in a laboratory with great precision and control. Experiments with ultracold atoms Bloch et al. (2008); Giorgini et al. (2008) have provided a possibility to realize such a scenario. The accuracy that can be reached in experiments with Fermi atomic gases and optical lattices is exceptional, thus offering a unique setting to explore highly correlated quantum fermion systems. In this paper, we demonstrate that, from the theoretical side, advances in computational methods now make it feasible to obtain numerically exact results for not only equilibrium properties, but also excited states. We compute the pairing gap, spectral functions and dynamical response functions in the two-dimensional Fermi gas across the range of interactions, which will allow direct comparisons with spectroscopy or scattering experiments. The dynamical properties provide a powerful tool to probe the behavior of the system and to visualize the crossover from a gas of molecules to a BCS superfluid.

We study the Fermi gas with a zero-range attractive interaction, which has generated a great deal of research activity Bloch et al. (2008); Giorgini et al. (2008); Carlson et al. (2003); Martiyanov et al. (2010); Fröhlich et al. (2011); Magierski et al. (2011); Wlazłowski et al. (2012); Makhalov et al. (2014); Ong et al. (2015); Murthy et al. (2015); Ries et al. (2015); Fenech et al. (2016); Boettcher et al. (2016); Carlson and Reddy (2005); Luo et al. (2016). The interest of the system is very wide, ranging from condensed matter physics Bauer et al. (2014); Magierski et al. (2011) to nuclear physics, with possible important applications also in the study of neutron stars Gezerlis and Carlson (2008); Lacroix (2016). This system describes experiments with a collection of atoms, for example 6Li, which are cooled to degeneracy in an equal mixture of two hyperfine ground states, labeled $|\uparrow\rangle$ and $|\downarrow\rangle$. Feshbach resonances allow the tuning of the interactions by varying an external magnetic field, making the system a unique laboratory to explore many-body physics Chin et al. (2010); Bloch et al. (2012). Starting from a weakly interacting BCS regime, where the attraction between particles induces a pairing similar to the one observed in ordinary superconductors, a crossover is observed as the interaction strength is increased, leading to a BEC regime where the Cooper pairs are tightly bound such that the system behaves as a gas of bosonic molecules. While both the BCS and the BEC regimes are well understood, the crossover regime provides an excellent example of a strongly interacting quantum many-body system Bloch et al. (2008); Giorgini et al. (2008).

We focus in particular on the two-dimensional (2D) Fermi gas, which has recently been realized experimentally using an highly anisotropic trapping potential Martiyanov et al. (2010); Fröhlich et al. (2011). The 2D system is important, since some of the most interesting physical phenomena, such as high temperature superconductivity Lee et al. (2006), Dirac fermions in graphene Castro Neto et al. (2009) and topological superconductors Qi and Zhang (2011), nuclear “pasta” phases Yu and Bulgac (2003) in neutron stars are two-dimensional in nature. Quantum fluctuations are known to be enhanced in 2D, making it even more important to have quantitative results beyond mean-field approaches.

Experiments are just beginning to measure properties in the 2D gas Feld et al. (2011); Makhalov et al. (2014); Ries et al. (2015); Boettcher et al. (2016); Cheng et al. (2016). An array of calculations have been performed Bauer et al. (2014); He et al. (2015); Klawunn (2016); Shi et al. (2015); Galea et al. (2016); Bertaina and Giorgini (2011), although much less is available in 2D compared to three-dimensional systems. The most direct connection with experiments is through response functions and accurate many-body data on dynamical response at low temperatures would provide crucial and fundamental missing link. However, these are much more challenging theoretically and computationally He (2016).

In this paper, we develop the capabilities to obtain unbiased results for imaginary-time correlation functions in spin-balanced Fermi gas systems, using first principles auxiliary-field quantum Monte Carlo (AFQMC) Zhang ; LeBlanc et al. (2015); Qin et al. (2016); Shi and Zhang (2013) methods. This provides a unique approach to excitations and dynamical response functions. Focusing on the BEC-BCS crossover regime, we compute the pairing gap as a function of the interaction strength, the spectral function, which can be measured experimentally in photoemission spectroscopy Feld et al. (2011), and the density and spin structure factors, which can be measured in two-photon scattering experiments Hoinka et al. (2012).

As the range of the interaction in the Fermi gas system of cold atoms is much smaller than the average inter-particle distance, the system can be modeled using a lattice Hamiltonian Werner and Castin (2012):

[TABLE]

where the label $i$ runs over a square lattice with $\mathcal{N}_{s}=L\times L$ sites hosting a total of $\mathcal{N}_{p}$ fermions, half with each spin $\sigma$ ($=\uparrow$ or $\downarrow$). The momentum ${\vec{k}}=(k_{x},k_{y})$ is defined on the reciprocal lattice with units $2\pi/L$ and $k_{x}$,$k_{y}\in[-\pi,\pi)$. The dispersion is $\varepsilon({\bf{k}})=k_{x}^{2}+k_{y}^{2}$ and $t=\hbar^{2}/(2mb^{2})$, with $b$ the lattice parameter. The attractive on-site interaction $U/t$ is tuned Werner and Castin (2012); Shi et al. (2015) for each lattice density $n=\mathcal{N}_{p}/\mathcal{N}_{s}$ and Fermi momentum $k_{F}=\sqrt{2\pi n}/b$ to produce the desired scattering length $a$, defined as the position of the node of the zero-energy $s$-wave solution of the two-body problem.

The ground state wave function $|\Psi_{0}\,\rangle$ of $\hat{H}$ is sampled using the AFQMC method Shi et al. (2015). For hamiltonian (1) with $U/t<0$, the sampling is not affected by the sign problem, so that numerical results can be obtained free of any bias for each set of parameters $\{\mathcal{N}_{s},\mathcal{N}_{p},U/t\}$. Accelerated sampling techniques with force bias are used, together with other technical improvements Shi and Zhang (2016), which greatly improves the efficiency of our calculations. This allows us to reach large system sizes in order to reliably extrapolate to the continuum and then to the thermodynamic limit Shi et al. (2015).

Our computation of the dynamical correlation functions here relies on a new algorithm which improved the computational scaling in the calculation of imaginary-time correlation functions Vitali et al. (2016) from $\mathcal{O}(\mathcal{N}_{s}^{3})$ in standard algorithms Hirsch (1985); Feldbacher and Assaad (2001); Assaad and Imada (1996); Motta et al. (2014, 2015); Alhassid, Y. et al. (2014) to $\mathcal{O}(\mathcal{N}_{s}\,\mathcal{N}_{p}^{2})$. The algorithm lets fluctuations related to creation/destruction operators or density/spin operators propagate in imaginary time, coupled to the stochastic evolution of the underlying AFQMC random walk or path-integral Vitali et al. (2016). The dynamical correlation functions are obtained as suitable combinations of matrix elements involving the Slater determinants Vitali et al. (2016). In the Fermi gas systems, the calculation is at the dilute limit, with $\mathcal{N}_{s}\gg\mathcal{N}_{p}$, so that a drastic speedup is achieved. This allows us to study lattices of $\mathcal{N}_{s}\sim 2000$ sites in order to, as illustrated below, reach proper convergence of the results to the realistic limit.

The exact imaginary-time correlation functions allow one to access a number of important physical quantities. We compute the pairing gap $\Delta$ from the large imaginary-time behavior of the dynamical Green’s functions:

[TABLE]

where the superscripts $p$ and $h$ indicate particle and hole, and $E_{0}$ is the ground-state energy. Moreover, we estimate the spectral function

[TABLE]

and the density and spin dynamical structure factors:

[TABLE]

where the operator $\hat{O}$ is $\rho_{\vec{k}}=\hat{n}_{\vec{k},\uparrow}+\hat{n}_{\vec{k},\downarrow}$ for density and $S_{\vec{k}}=(\hat{n}_{\vec{k},\uparrow}-\hat{n}_{\vec{k},\downarrow})/2$ for spin, and the brackets indicate ground-state expectations. These functions are obtained from analytic continuation of the imaginary-time Green’s functions and density-density or spin-spin correlation functions, using the Genetic Inversion via Falsification of Theories (GIFT) method Vitali et al. (2010); Bertaina et al. (2017, 2016); Motta et al. (2016); Arrigoni et al. (2013); Nava et al. (2013); Saccani et al. (2012).

Figure 1 shows the computed pairing gap across different interaction strengths. The most standard approach to determine the pairing gap requires separate calculations corresponding to different particle numbers. While the spin-balanced calculation is free of the sign problem, the $({\mathcal{N}_{p}\pm 1})$ calculations are not. This makes our approach through imaginary-time Green’s functions in Eq. (2) especially advantageous, since the Monte Carlo sampling remains at ${\mathcal{N}_{p}}$ and thus sign-problem-free Vitali et al. (2016). We are able to determine the Green’s function with high statistical accuracy in the asymptotic regime in imaginary time $\tau$, so that the behavior is dominated by an exponential whose exponent gives the quasi-particle peak at the targeted momentum ${\vec{k}}$ (see Supplementary Materials). A double exponential function is used in the fit to account for any residual effects. The imaginary-time interval on which the fit is performed is stochastically varied and sampled, and multiple data sets are generated to remove statistical correlations in imaginary time within each run. The final statistical uncertainty reflects the combined effects from the AFQMC error bars and the fitting procedure. We then scan ${\vec{k}}$ to locate the minimum/maximum (for the particle/hole Green’s functions) for the pairing gap, as further illustrated in Fig. 2.

In Fig. 1 we also show the BCS mean-field prediction, as well as the current best many-body results, from recent diffusion Monte Carlo (DMC) calculations Galea et al. (2016); Bertaina and Giorgini (2011). It is seen that our pairing gap is compatible with the DMC results on the BEC side of the crossover, but is consistently smaller for larger values of $\ln(k_{F}a)$. The smaller gap value is not surprising, since the DMC contains a fixed-node (FN) approximation which gives an upper bound on the computed energy. It is reasonable to expect that the trial wave function used for FN is of higher quality for the spin-balanced system compared than for the $(\mathcal{N}_{p}\pm 1)$ systems, which would lead to an overestimation of the pairing gap. Our results on the BCS side are consistent with the rescaled BCS results $\Delta_{BCS}/e$ from the theory by Gorkov and Melik-Barkhudarov, which is expected to be exact in the BCS limit ($\log(k_{F}a)>>1$) Petrov et al. (2003); Gorkov and K. (1961).

Figure 2 plots the computed quasi-particle peaks as a function of $k\equiv|{\vec{k}}|$, together with the spectral function, for four values of the interaction parameter. The zero of the energy is set equal to the chemical potential, which we can compute exactly Shi et al. (2015). We will refer to the function $A(\vec{k},\omega)$ as particle and hole spectral function respectively for $\omega>\mu$ and $\omega<\mu$. The particle spectral function originates from the first term on the right in Eq. (3), physically representing states available for additional particles injected into the system, while the hole spectral function, originating from the second term, contains information about states occupied by the particles in the system, which are thus accessible by the creation of holes. In each panel, we show also mean-field prediction for the quasi-particle energies Randeria et al. (1989): $E_{\pm}(\vec{k})=\pm\sqrt{\left(\hbar^{2}k^{2}/2m-\mu_{\rm BCS}\right)^{2}+\Delta_{\rm BCS}^{2}}$, where $\Delta_{\rm BCS}$ is the gap and $\mu_{\rm BCS}$ the chemical potential in BCS theory. The non-interacting spectral function, $A^{0}(\vec{k},\omega)=\delta\left(\omega-\left(\hbar^{2}k^{2}/2m-\varepsilon_{F}\right)\right)$, is also shown for reference. In the AFQMC spectral functions obtained from the GIFT analysis, shown in the color plot, quasi-particles peaks are still visible, which are broadened from many-body correlations, resulting in a non-zero imaginary part of the self-energy, and are renormalized with respect to the BCS dispersion relations. The quasi-particle peaks computed directly from AFQMC are shown by symbols. These were obtained following the procedure described in Fig. 1. Results from different system sizes are shown, which indicate convergence to the bulk limit within numerical resolution. (Separate calculations were also carried out to verify that these densities are indistinguishable from the dilute limit Shi et al. (2015).)

The behavior of the spectral function provides a clear visualization of the BEC-BCS crossover. In the BEC regime at $\ln(k_{F}a)=0$, a large gap, of the order of the energy needed to break a molecule, separates the two branches, which are roughly momentum-independent for $k\leq k_{F}$. A smooth evolution of the spectral function is observed. In the BCS regime at $\ln(k_{F}a)=1.5$, it starts to resemble the non-interacting behavior, where a gap is still present at the Fermi momentum, as in conventional superconductors. The intermediate values of the interaction show a smooth crossover between the two regimes. Viewed in the reverse direction, gradual and significant departures from the BCS results are seen as the interaction strength is increased.

We also compute two-body dynamical correlations in imaginary time, which can again be obtained using our method with computational cost linear in ${\mathcal{N}}_{s}$ Vitali et al. (2016). From these, we apply analytic continuation to obtain the density and spin dynamical structure factors, $S^{\rho}(\vec{k},\omega)$ and $S^{S}(\vec{k},\omega)$, which can be measured experimentally using two-photons Bragg spectroscopy Hoinka et al. (2012). In particular, the high momentum behavior is very interesting as it provides a highly sensitive probe of the BEC-BCS crossover. We focus our attention on $k=4k_{F}$, close to the value recently investigated experimentally in three-dimensions Hoinka et al. (2012).

The results are plotted in Fig. 3 as functions of the frequency $\omega$ for the four values of the interaction parameter. In addition to AFQMC, we have also performed self-consistent dynamical BCS theory calculations for the same system, following the approach in Ref. Combescot et al. (2006) which studied the three-dimensional Fermi gas. The results are shown in the figure for comparison. Because the theory yields the response functions directly, it helps to provide an additional gauge on the reliability of analytic continuation analysis. We observe that the dynamical BCS theory gives results on the response functions that are qualitatively reasonable. Significant differences arise from the AFQMC results, however, for example in the peak position in the spin structure factor for strong interactions. Direct comparisons of the imaginary-time correlation functions show significant differences between AFQMC and dynamical BCS theory as well, manifesting particle correlation effects absent in the latter.

In the density response, a large peak is seen at $\omega\simeq\omega_{R}/2$ in the deep BEC regime. Since the particles are tightly paired to form molecules in this regime, the response of the system at high momentum is dominated by the recoil of the molecules themselves, whose mass is twice the atomic mass. In contrast, the response on the BCS side is simply a free particle recoil with the bare mass of the atoms. The behavior of the density response in the crossover regime interpolates between the two physical pictures, as is evident from Fig. 3. Starting from $\ln(k_{F}a)=0$, we observe a gradual shift of the spectral weight from the dominant molecular contribution towards the second peak at $\omega\simeq\omega_{R}$. At $\ln(k_{F}a)=1$, the second peak dominates, and the molecule peak almost disappears. By $\ln(k_{F}a)=1.5$, the response becomes qualitatively similar to the non-interacting one.

The spin response, on the other hand, is not sensitive to the molecular mode at $\omega\simeq\omega_{R}/2$, since the positive and negative fluctuations on the spin-$\uparrow$ and spin-$\downarrow$ particles cancel each other. However, we observe that, as it happens in three-dimensions Hoinka et al. (2012), the intensity of the peak is smaller on the BEC side of the crossover, and the position of the peak is shifted towards higher energies. This corresponds to a suppression of the spin susceptibility, related to the increased energy required to remove atoms from the molecules.

In summary, we have performed ab initio calculations of the pairing gaps and dynamical correlation functions for the two-dimensional interacting Fermi atomic gas. Numerically exact AFQMC predictions are provided for the pairing gap. From unbiased imaginary-time correlation functions computed by AFQMC for the many-body ground state, the spectral function and the density and spin dynamical structure factors are obtained, via analytic continuation, across the BEC-BCS crossover. Much larger system sizes are reached in our simulations by the development and implementation of several technical advances. Many internal validations and self-consistency checks are performed and careful error quantifications are carried out to maximize the robustness and reliability of the results. The results will allow benchmarks of further theoretical and computational developments, and direct comparisons with experiments. The exact pairing gaps will also be crucial as an input for formulating density-functional theory in 2D Bulgac (2013); Zou et al. (2016). Excitations and dynamical correlation functions provide excellent tools for visualizing the BEC-BCS crossover. In interacting many-fermion systems in general, they connect directly with experimentally accessible measurements. Our approach opens up many new possibilities for the computational studies of strongly interacting fermionic cold atomic systems. It is hoped that the results presented here will also serve as an illustration of state-of-the-art computational capabilities, and will stimulate additional theoretical and experimental activities. The feedback from such activities will in turn spur further computations and additional developments.

We thank J. Carlson, A. Gezerlis, and Lianyi He for helpful discussions. This work was supported by NSF (Grant No. DMR-1409510). E.V. and S. Z. were also supported by the Simons Foundation. Computing was carried out at the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575, and the computational facilities at William and Mary.