Anisotropic collisions of dipolar Bose–Einstein condensates in the universal regime

In general, the description of quantum mechanical scattering becomes simpler when the relative kinetic energy between the collision partners is negligibly small. Specifically, in this limit the de Broglie wavelength of relative motion greatly exceeds the length scale over which the two particles exert forces on one another, rendering details of their interaction, if not insignificant, at least fairly simple to account for [1]. This concept was first articulated by Enrico Fermi in the context of collisional broadening of spectral lines in a gas of Rydberg atoms [2]; was instrumental in the understanding of low-energy scattering of neutrons from nuclei [3]; and has found expression in modern times as the bedrock upon which our understanding of ultracold gases is based [4]. Ordinarily, the key aspect of this simplicity is that the scattered state contains only a single (or very few) eigenstate(s) of angular momentum, i.e., one or only a few partial waves. Ultracold identical bosons have a differential cross section which is isotropic and independent of the collision energy such that a single length scale, referred to as the scattering length, suffices for its description. The situation is only marginally more complicated for ultracold identical fermions in which the differential cross section is $\propto {k}^{4}| \hat{k}\cdot \hat{k}^{\prime} {| }^{2}$, where k is the magnitude of the relative wave vector and $\hat{k}$ ($\hat{k}^{\prime}$) is a unit vector along the direction of relative incoming (outgoing) momentum [5]. Once again, a single scalar quantity, the scattering volume, sufficiently characterises the underlying potential. These simple results dictate many important dynamical properties of quantum gases such as the efficiency of evaporative cooling [6–9] and the speed of sound [10].

However, this situation is radically altered when the potential energy between the collision partners does not decay suitably fast. If the interaction potential decays as $1/{r}^{n}$, where r is the distance separating the two particles, then for $n\leqslant 3$ a different scenario emerges [11, 12]. Ultracold gases containing atoms or molecules which possess magnetic [13–18] or electric [19–22] dipole moments therefore provide a magnificent opportunity to observe a novel regime of low-energy quantum scattering. When the dipoles are aligned along a chosen axis by an externally applied field, the overall magnitude of the cross section is still determined by a characteristic length scale, the dipole length, which depends on the dipole moment [23]. However, the dipole alignment direction and the $1/{r}^{3}$ asymptotic decay of the interaction potential conspire to create a differential cross section which now involves a very large set of partial waves in an essential way due to the coupling between different angular momentum states. Moreover, it presents a novelty in that such a differential cross section depends explicitly on the relative momenta before the collision (as well as after) with respect to the polarization axis of the dipoles [24]. The consequence of this dependence has been measured indirectly through its effect on the equilibration rate of a dipolar gas, which is taken out of equilibrium by diabatically changing the trap along a certain direction [25, 26]. It was seen that the equilibration rate can vary by as much as a factor of four, depending on the angle between the dipole alignment direction and the dynamic axis of the trap.

We report a direct measurement of the differential scattering cross section of ¹⁶²Dy, which has an exceptionally large magnetic moment of $9.93\,{\mu }_{{\rm{B}}}$, where ${\mu }_{{\rm{B}}}$ is the Bohr magneton. We obtain this measurement by colliding two Bose–Einstein condensates (BECs) at a relative velocity significantly greater than the width of their respective internal velocity distributions. Furthermore, most of these measurements are made in the dilute limit where the majority of particles pass through the opposing condensate without experiencing a collision. In this situation, scattered particles emerge on a spherical halo defined by the constraints of conserved energy and momentum. The angular distribution of the particles on the sphere is almost entirely determined by the differential cross section. We vary the dipole alignment direction relative to the collision axis and take absorption images of the post-collision number density after a long time-of-flight (TOF). To obtain a detailed theoretical understanding of these images, we perform Monte-Carlo simulations, which capture the essential features of the experiment. Such simulations yield an approximate understanding of how multiple-collisions (defined when an individual atom collides again after it has been ejected from its original condensate) can corrupt the direct correspondence between the observed halo and the differential cross section.

Related experiments on bosonic alkali gases which also observed halo-like structures due to atomic collisions were reported in [27–30] and recently in fermionic ⁴⁰K [31, 32]. Remarkably beautiful experiments were reported in [33] in which metastable helium was employed. In addition to the halo structure, the internal energy of metastable helium allowed for pair correlations to be measured via spatially and temporally resolved single-atom counting. This stimulated theoretical interest in the strength of such correlations [34–36]. All effects from dipolar interactions are negligible in these experiments involving alkali and He gases, and the reported halo is well described by a very small set of partial waves. A unique experiment reported in [37] studied the creation of artificial partial waves in ⁸⁷Rb BEC collisions mediated by an optical potential which modifies the constraints due to energy and momentum conservation. Using this approach they were able to engineer differential scattering with several partial waves, in spite of the fact that a bare collision would have been entirely s-wave. Our work presents a natural and intriguing extension of these experiments, demonstrating the nature of differential scattering in the presence of dipolar interactions which are anisotropic and unavoidably couple to a large number of angular momentum states during the collision.

We produce a single BEC of ¹⁶²Dy spin polarized in the J = 8, ${m}_{J}=-8$ absolute ground state with $6.3(4)\times {10}^{4}$ atoms as described in [38]. An optical lattice is used to diffract the condensate into the $2n{\hslash }{{\bf{k}}}_{L}$ momentum states, where the lattice wave vector ${{\bf{k}}}_{L}=2\pi \sin (\theta /2)/\lambda$ depends on the lattice lasers wavelength λ and alignment angle θ, and n is an integer which labels the different diffraction orders. The optical lattice is formed by two fiber-coupled beams derived from the same Ti:Sapphire laser with $\lambda =741\,\mathrm{nm}$ that propagate along $\hat{x}+\hat{y}$ and $\hat{y}-\hat{x}$ such that $\theta =\pi /2$, as shown in figure 1(a). Both beams are linearly polarized along $\hat{z}$.

Zoom In Zoom Out Reset image size

The condensate is split into the $\pm 2{\hslash }{{\bf{k}}}_{L}$ diffraction orders with high efficiency by using a precisely timed two-pulse sequence [39]. Thus, two spatially overlapping condensates are produced and collide with relative momentum ${{\bf{p}}}_{\mathrm{rel}}=4{\hslash }{k}_{L}\hat{x}$ as they separate. Immediately after applying the lattice grating, all trapping potentials are removed, and the two condensates collide and expand. After a 22 ms TOF, a column-integrated density profile of the atoms in the $\hat{x}$–$\hat{z}$ plane is measured via absorption imaging along $\hat{y}$. We manipulate the dipole alignment direction by applying a bias magnetic field relative to both the collision axis and the imaging axis. In all configurations, the magnetic field magnitude is held at 1.58 G and away from any Feshbach resonance [40]. Three instructive cases are shown in figure 1, which provides a schematic of our experiment. The differential cross section is isotropic when the dipole alignment is parallel to the relative momentum, as in figure 1(a). However, with the dipole alignment perpendicular to the relative momentum, as in figures 1(b) and (c), the differential cross section is clearly anisotropic. While the cross section is the same in panels (b) and (c), the rotation with respect to the imaging direction allows different projections of the 3D halo to be imaged.

3.1. Two-body scattering theory

The central theoretical element is the differential scattering cross section for identical bosons interacting via a combination of dipolar and short-range interactions. This is found by solving the two-body Schrödinger equation in the center-of-mass rest frame,

$$\begin{eqnarray}&&\left[-\displaystyle \frac{{{\hslash }}^{2}}{m}{{\rm{\nabla }}}_{{\bf{r}}}^{2}+{V}_{{\rm{d}}}(\hat{\epsilon },{\bf{r}})+{V}_{\mathrm{sr}}({\bf{r}})\right]\psi ({\bf{r}})=\displaystyle \frac{{{\hslash }}^{2}{k}^{2}}{m}\psi ({\bf{r}}),\end{eqnarray} \tag{ 1 }$$

where m is the mass of a single particle and ℏ is the reduced Planck's constant. The dipole–dipole potential energy for aligned dipoles is ${V}_{{\rm{d}}}(\hat{\epsilon },{\bf{r}})=(2{{\hslash }}^{2}{a}_{d}/m)[1-3{(\hat{\epsilon }\cdot \hat{r})}^{2}]/{r}^{3}$, where $\hat{\epsilon }$ is the direction of alignment. The potential energy from short-range van der Waals interactions can be approximated by ${V}_{\mathrm{sr}}({\bf{r}})=(4\pi {{\hslash }}^{2}a/m){\delta }^{(3)}({\bf{r}})$, provided one stays within the first-order Born approximation. The strength of the dipole interaction is determined by the length scale ${a}_{d}={\mu }_{0}{\mu }^{2}\,m/8\pi {{\hslash }}^{2}$, where ${\mu }_{0}$ is the vacuum permeability and μ is the magnetic moment of a single particle [41]. Similarly, the strength of the short-range interaction is determined by the s-wave scattering length a. This problem is solved in detail in [24], see also [12]. Briefly, in the limit where $\max (| a| ,{a}_{d})\times k\ll 1$, a solution can be found which takes the asymptotic form $\psi ({\bf{r}})={{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}+f(\hat{\epsilon },\hat{k},\hat{k}^{\prime} ){{\rm{e}}}^{{\rm{i}}{k}{r}}/r$, where ${{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}$ is an incoming wave and f is the scattering amplitude that depends on the dipole alignment direction, the incoming wave vector ${\bf{k}}=k\hat{k}$, and the outgoing wave vector ${\bf{k}}^{\prime} =k\hat{k}^{\prime}$. The differential scattering cross section for identical bosons is given by $\tfrac{{\rm{d}}\sigma }{{\rm{d}}{\rm{\Omega }}}=\tfrac{1}{2}| f(\hat{\epsilon },\hat{k},\hat{k}^{\prime} )+f(\hat{\epsilon },\hat{k},-\hat{k}^{\prime} ){| }^{2}$, which in the low energy limit and within the first-order Born approximation, is found to be

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\sigma }{{\rm{d}}{\rm{\Omega }}}=2{a}_{d}^{2}{\left[\displaystyle \frac{{(\hat{k}\cdot \hat{\epsilon })}^{2}+{(\hat{k}^{\prime} \cdot \hat{\epsilon })}^{2}-2(\hat{k}\cdot \hat{\epsilon })(\hat{k}^{\prime} \cdot \hat{\epsilon })(\hat{k}\cdot \hat{k}^{\prime} )}{1-{(\hat{k}\cdot \hat{k}^{\prime} )}^{2}}-\displaystyle \frac{2}{3}+\displaystyle \frac{a}{{a}_{d}}\right]}^{2},\end{eqnarray} \tag{ 2 }$$

where $\hat{k}$ ($\hat{k}^{\prime}$) is a unit vector along the direction of incoming (outgoing) relative momentum. Note that our formula for $\tfrac{{\rm{d}}\sigma }{{\rm{d}}{\rm{\Omega }}}$ accounts for bosonic indistinguishability, relevant for our experiment. Alternatively, the differential cross section for fermions would be found from $\tfrac{1}{2}| f(\hat{\epsilon },\hat{k},\hat{k}^{\prime} )-f(\hat{\epsilon },\hat{k},-\hat{k}^{\prime} ){| }^{2}$. Plots of the differential cross section for three different dipole alignments are shown in figures 1(d)–(f). For our experimental parameters the dimensionless parameter ${{ka}}_{d}/(2\pi )\approx 0.04$ (corresponding to a de Broglie wavelength much larger than the dipole length). In [66] we calculated the correction due to the second-order term in the Born approximation, from which we can estimate the effects to lie within the $0.6 \%$ level at our current collision energy. The total cross section, found by integrating equation (2) over all possible outgoing directions, is given by [24]

$$\begin{eqnarray}&&\sigma (\eta )={a}_{d}^{2}\displaystyle \frac{\pi }{9}\left\{72\displaystyle \frac{{a}^{2}}{{a}_{d}^{2}}-24\displaystyle \frac{a}{{a}_{d}}[1-3{\cos }^{2}(\eta )]+11-30{\cos }^{2}(\eta )+27{\cos }^{4}(\eta )\right\},\end{eqnarray} \tag{ 3 }$$

where η is the angle between $\hat{k}$ and $\hat{\epsilon }$. We note that the spatial anisotropy of this cross section is caused by $\hat{\epsilon }$. Equations (2) and (3) are universal in the sense that they are insensitive to any details of the short-range potential.

3.2. Many-body considerations: Monte-Carlo simulation

When two atoms from different condensates collide, their relative velocity is much larger than the width of each condensate's internal velocity distribution. Therefore, these primary collisions effectively occur at a fixed angle between $\hat{k}$ and $\hat{\epsilon }$. However, each atom then has a finite probability of a secondary scattering event. These secondary scattering events occur with an essentially random angle between $\hat{k}$ and $\hat{\epsilon }$ and corrupt the direct correspondence between the observed halo and the differential scattering cross section. For this reason we use a direct-simulation Monte Carlo (DSMC) algorithm, which is able to keep track of such multiple collision events, in order to quantitatively understand the experimental data as accurately as possible.

The system dynamics can be separated into a low energy part (relevant to the two condensates) and a high energy part (relevant to the halo). These two parts have high and low phase space densities, respectively. Because we are primarily interested in the halo, we focus on the classical kinetic equation (the Boltzmann equation) for the phase space distribution function, which is given by

$$\begin{eqnarray}&&{\partial }_{t}f({\bf{r}},{\bf{p}},t)+\displaystyle \frac{{\bf{p}}}{m}\cdot {{\rm{\nabla }}}_{{\bf{r}}}f({\bf{r}},{\bf{p}},t)=C[f],\end{eqnarray} \tag{ 4 }$$

where the left-hand side contains free-streaming terms and the right-hand side includes effects from two-body collisions. This is given by

$$\begin{eqnarray}&&C[f]=\int \displaystyle \frac{{{\rm{d}}}^{3}{{\bf{p}}}_{1}}{{h}^{3}}\int {{\rm{d}}}^{2}\hat{{\rm{\Omega }}}\displaystyle \frac{{\rm{d}}\sigma }{{\rm{d}}{\rm{\Omega }}}{v}_{r}[f^{\prime} {f}_{1}^{\prime }-{{ff}}_{1}],\end{eqnarray} \tag{ 5 }$$

where $f=f({\bf{r}},{\bf{p}})$, ${f}_{1}=f({\bf{r}},{{\bf{p}}}_{1})$, $f^{\prime} =f({\bf{r}},{\bf{p}}^{\prime} )$, and ${f}_{1}^{\prime }=f({\bf{r}},{{\bf{p}}}_{1}^{\prime })$ account for the four momenta (two incoming, ${\bf{p}}$, ${{\bf{p}}}_{1}$, and two outgoing, ${\bf{p}}^{\prime}$, ${{\bf{p}}}_{1}^{\prime }$) associated with a two-body collision, and ${v}_{r}=| {\bf{p}}-{{\bf{p}}}_{1}| /m$ is the relative velocity. Note that ${\bf{p}}={\hslash }{\bf{k}}$ connects the momentum to the wavevector, which was used in section 3.1. The momenta in equation (5) are related by energy and momentum conservation such that ${\bf{p}}+{{\bf{p}}}_{1}={\bf{p}}^{\prime} +{{\bf{p}}}_{1}^{\prime }$ and $| {\bf{p}}-{{\bf{p}}}_{1}| =| {\bf{p}}^{\prime} -{{\bf{p}}}_{1}^{\prime }|$. The integration variable is $\hat{{\rm{\Omega }}}=({\bf{p}}^{\prime} -{{\bf{p}}}_{1}^{\prime })/| {\bf{p}}^{\prime} -{{\bf{p}}}_{1}^{\prime }|$. Equation (4) does not provide an accurate description of the atoms within the high-phase-space-density regions, i.e., the condensates. These would presumably be described within the context of a Gross-Pitaevskii (GP)-type theory of BEC evolution [42–44]. However, the goal of simulating both the high and low-phase-space-density components of the gas is beyond the scope of our current work. We simply wish to approximately capture the density evolution of the condensates. This in turn yields a reasonable prediction for the scattering rates and therefore the formation and deformation of the scattering halo.

We solve equation (4) numerically using a DSMC algorithm, details of which have been published in [45]. We use an initial condition corresponding to both spatial and momentum densities proportional to the probability distribution function,

$$\begin{eqnarray}&&P({\bf{R}};{\bf{x}})=\max \left\{\displaystyle \frac{15}{8\pi \prod _{\alpha }{R}_{\alpha }}\left[1-\displaystyle \sum _{\alpha }\displaystyle \frac{{x}_{\alpha }^{2}}{{R}_{\alpha }^{2}}\right],0\right\},\qquad \alpha =\{1,2,3\},\end{eqnarray} \tag{ 6 }$$

where ${\bf{x}}$ (R_α) can denote either spatial or momentum coordinates (widths). The spatial part corresponds to a Thomas−Fermi condensate density, with the widths calculated in the manner prescribed by [46, 47]. The momentum widths are found by fitting to the experimental image of the expanded condensate after a 22 ms TOF. The system is then divided into two halves, which propagate along the positive and negative x-axis, respectively.

For a detailed account of the simulation method, we refer the reader to [45], which adjusts the original DSMC approach [48] into a version appropriate to dipolar gases. Similar methods have been employed in the study of ultracold gases, see for instance [49–51]. In particular, [49] also studied halo formation in the case of s- and d-wave collisional cross sections. Briefly, our computational algorithm uses test-particles with phase-space coordinates which are sampled from the initial distribution in equation (6). The test particles move classically from one time-step to the next. At each time-step particles are binned in position space. The bin-size represents the finite resolution of the delta-function within the numerics. Within each bin, the collision probability for each pair of particles is evaluated according to equation (3) using the correct value of $\hat{k}$ corresponding to each pair. Collisions are then chosen to occur stochastically in accordance with these probabilities. The post-collision velocities are also chosen stochastically, in accordance with the differential scattering cross section in equation (2).

3.3. Validity and limitations of the theoretical model

The experimental absorption images are presented in the first column of figure 2, i.e., in panels (a1), (b1), ... (e1). Each image corresponds to the average of forty experimental runs at a fixed dipole alignment direction, as stated in the figure caption. Imperfections in the diffraction leaves residual atoms with a momentum near $0{\hslash }{{\bf{k}}}_{L}$. This was also observed in [62]. These atoms contribute noise to the measurement both near $0{\hslash }{{\bf{k}}}_{L}$ and in the $| 2{\hslash }{{\bf{k}}}_{L}|$ halos from collisions with the $\pm 2{\hslash }{{\bf{k}}}_{L}$ BECs. However, the effect of the dipole alignment direction is clear from the remarkable variation between images (a1), (b1), ... (e1). We run a family of simulations for each value of the dipole alignment direction, but with varying s-wave scattering length⁵ . These families of simulation results are shown as the smaller images in figure 2 panels (a2)–(a9), (b2)–(b9), ... (e2)–(e9). In addition to the simulations which we present here, we have also looked at initializing the test particles using Gaussian distribution functions in place of equation (6), and also small variations in momentum and spatial widths R_α. However, we found that such variations have very little effect on our conclusions.

Zoom In Zoom Out Reset image size

We perform a pixel-by-pixel comparison between the experimental and simulated scattering images. Each pixel corresponds to an area of 2.6 μm × 2.6 μm. For each dipole alignment direction, the average experimental image ${ \mathcal E }$ is obtained as well as the standard error for each pixel ${{ \mathcal E }}^{\sigma }$. The experimental image and each simulation image ${ \mathcal S }$ (including unscattered atoms) is normalized to suppress the effects of atom number variations, and a mask is then applied to the images to exclude the unscattered BECs and atoms not diffracted by the lattice. An example of a masked image is shown in figure 3(f). The weighted sum of the squared residuals is calculated for each simulated image $S={\sum }_{{ij}}{({{ \mathcal E }}_{{ij}}-{{ \mathcal S }}_{{ij}})}^{2}/{({{ \mathcal E }}_{{ij}}^{\sigma })}^{2}$, where the subscript ij denotes the pixel in the ith row and jth column of an image. This corresponds to the cost function for a weighted least squares regression. The results of this analysis are shown in figures 3(a)–(e). Though a particular choice of scattering length, a, minimizes S for each dipole alignment direction, the analysis is not sensitive enough to define a 1σ confidence interval. However, as seen in figure 3, all the minima fall within the range of 120a₀–160a₀, which is consistent with measurements reported in our previous work [26, 66]. We believe that the fluctuations in the location of the minima in figure 3 are primarily due to an incomplete knowledge of the initial density of the BEC. Recent results suggest that dipolar interactions in strongly dipolar BECs can result in the formation of droplets which require beyond-mean-field effects to be quantitatively understood, see [67–69]. While we expect an initial peak density of $2.1(2)\times {10}^{14}$ cm⁻³ within the mean-field theory, we have run additional simulations with slight variations in the spatial widths of equation (6) which determines initial density. These variations lead to the density increasing (decreasing) by as much a factor of 1.37 (0.67). The results (for dipoles aligned along the $\hat{x}$-axis) are shown in figure 3(a). As can be observed by the green line in figure 3(a), running a simulation at a higher density generates a small but noticeable shift in the minimum towards a lower scattering length. This occurs because the increased density generates a larger number of scattered atoms in the halo, and therefore a shift toward a lower scattering length compensates against this change. Similar statements can be made about the simulations at a reduced density.

Zoom In Zoom Out Reset image size

We also analyze the collision frequency and, through a detailed comparison between simulation and experiment, establish the relevance of multiple collisions. These are collisions involving atoms that have already been scattered out of the original condensates. One might expect that a simple prediction for the mean number of collisions might be ${{\rm{\Gamma }}}_{c}=\langle n\rangle \sigma {p}_{\mathrm{rel}}/m$, where $\langle n\rangle$ is the mean density. However, because the expansion occurs over a time scale comparable with the time scale over which the clouds overlap, this formula significantly over-predicts the number of collisions. Unfortunately, a simple analytic theory of expansion in dipolar condensates is not yet available, so we must turn to the simulations which account empirically for condensate expansion. The nature of the DSMC method allows one to label particles and keep track of collisions as they occur. The results are shown in figure 4 for dipoles aligned along $\hat{x}$, a case in which the presence of multiple collisions is strong. Figure 5 shows results where the dipoles are aligned along both $(\hat{x}+\hat{z})/\sqrt{2}$ and $\hat{z}$, and the presence of multiple collisions is less relevant. The data can be reasonably well understood starting from equation (3) and noting that, when $a=140{a}_{0}$, the total cross section has a maximum at $\eta =0^\circ$ (corresponding to dipoles aligned along $\hat{x}$) and a minimum at $\eta =90^\circ$ (corresponding to dipoles aligned anywhere in the $\hat{y}$–$\hat{z}$ plane). Specifically, we have $\sigma (0^\circ )/\sigma (90^\circ )\approx 2.6$ and⁶ $\sigma (0^\circ )/\sigma (45^\circ )\approx 1.6$. For this reason, we observe the unscattered fraction of the gas increases from ∼40% in figure 4(b), to ∼60% in figure 5(b), and to ∼80% in figure 5(f). From this, we conclude that the experimental absorption image with dipoles aligned along the $\hat{x}$ axis is strongly affected by atoms which have undergone multiple collisions. Indeed, these atoms appear to constitute the majority of the absorption image in this case. Reducing the initial density to suppress multiple collisions is not practicable due to the sharp loss in signal-to-noise. However, for all other alignment directions the experimental absorption image (away from the condensate regions) is dominated by a single scattering event, and is therefore in close correspondence with the differential scattering cross section of equation (2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In conclusion, we have measured the differential scattering between identical dipolar bosons in the low-energy (universal) regime by imaging the halos of atoms scattering from two colliding Dy BECs. The results depend strongly on the angle between the dipole alignment direction and the collision axis and are well described by the analytic formula in equation (2), derived under the first-order Born approximation. Although it is not intended to be entirely quantitative, a classical Monte-Carlo simulation of the many-body dynamics provides a reasonable higher-order approximation to account for finite-density effects and finite-momentum distribution widths. We allow the s-wave scattering length to vary within the simulations, and note that discrepancies between simulation and experiment are minimized at a scattering length which remains consistent with previous measurements [26, 66].

The measurements provide a beautiful demonstration of the theoretical prediction for a low energy scattering amplitude, which dates back to O'Malley in 1964 [12]. The envisaged physical system under consideration at that time was not dipolar collisions, but rather the scattering of an electron by a non-spherical atom, such as atomic oxygen. However, we see that the interaction potential is asymptotically equivalent in the two cases. Future work may use such experiments to probe the complex collisional physics of dipolar condensates and degenerate Fermi gases near Feshbach resonances in and among the dense and ultradense spectra observed in the dysprosium and erbium systems [17, 40, 65, 70–72].

We gratefully acknowledge John Bohn for early contributions to this work. We thank Jack DiSciacca and Wil Kao for experimental assistance and Ian Spielman for a stimulating discussion. We are grateful to the NSF and AFOSR for funding support. YT acknowledges partial support from a Stanford Graduate Fellowship. The research leading to these results received funding from the European Unions Horizon 2020 research and innovation programme under grant agreement No 658311.