Signatures of the BCS-BEC crossover in the yrast spectra of Fermi quantum rings

Ulrich Ebling, 2, ∗ Ali Alavi, 4 and Joachim Brand 2, 3 Centre for Theoretical Chemistry and Physics, New Zealand Institute for Advanced Study, Massey University, Private Bag 902104, North Shore, Auckland 0745, New Zealand Dodd-Walls Centre for Photonic and Quantum Technologies, PO Box 56, Dunedin 9056, New Zealand Max Planck Institute for Solid State Research, Heisenbergstraße 1, 70569 Stuttgart, Germany Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW, United Kingdom (Dated: December 1, 2020)

Of particular interest are ring configurations, where translational invariance along one spatial dimension makes (angular) momentum a good quantum number. This allows for the study of yrast states, which are defined as the lowest energy state at given value of the total momentum. Yrast states in a bosonic superfluid are intimately connected [24][25][26][27][28][29][30][31] to localized nonlinear waves known as dark solitons [32]. It was shown that measuring the position of all bosons in an yrast state reveals a darksoliton-like particle depletion [30], and that wave-packetlike superpositions of yrast states emulate the behavior of classical dark solitons [31]. While the yrast states are fragmented quantum condensates, breaking the translational symmetry restores single condensation and classical soliton features [28]. Dark solitons in Fermi superfluids have been identified in experiments [33][34][35], but many predictions from mean-field theory have not yet been tested [36][37][38][39][40][41][42]. Moreover, there is an intriguing connection [43,44] between dark solitons and the predicted Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase * uebling97@gmail.com of imbalanced superfluids [45,46]. Dispersion relations of yrast states were analyzed in the context of dark soliton physics in the Yang-Gaudin model, a one-dimensional (1D) Bethe-ansatz solvable model of a fermionic superfluid in Ref. [47]. An overview of computational studies of quantum rings can be found in the recent review literature [48][49][50].
Beyond the purely one-dimensional models of quantum rings, a second spatial dimension can be added by considering stripe-or ladder-type lattice configurations as it is done in this work. In Refs. [48,49] these are referred to as quasi-1D geometries. While the Bethe ansatz is unavailable for such models and mean-field theory is not valid in the strongly-correlated regime, direct numerical simulation is very challenging due to the fact that Hilbert space size increases exponentially both with particle number and the number of lattice sites. Quantum Monte Carlo (QMC) simulations are still possible, although the fermion sign problem [51,52] provides a challenge for the simulation of fermionic many-body problems.
QMC methods relevant in the field of ultracold quantum gases each have their own strengths and weaknesses. Diffusion Monte Carlo has no basis set dependence but either converges to a bosonic ground state or requires node-fixing, which introduces an approximation [53]. Auxiliary-field QMC is sign-problem free for the attractive balanced Hubbard model [54], but the Hubbard-Stratonovich transformation involved in this method breaks symmetries of the Hamiltonian and thus does not allow for the study of yrast states. Recent work has suggested a solution [55] but it has yet to be seen whether the method can be implemented efficiently. Determinant Monte Carlo and related methods can handle finite temperature and extrapolate to zero temperature, but they fix the chemical potential instead of particle number [56][57][58][59]. All these existing methods have in common that they can study overall ground-state properties while it is not possible to study yrast states, because the total momentum cannot be easily constrained.
For this work we use Full-Configuration-Interaction Quantum Monte Carlo (FCIQMC), a method originally developed for strongly correlated electrons in the context of quantum chemistry [60,61]. FCIQMC has been applied with great success to a large number of problems in this field [62,63], and recently to ultracold atoms [64,65]. This method can find the ground state energy and many-body wave function in a Fermi system by expanding the wave function into a set of Slater determinants. A stochastic version of exact diagonalization of the Hamiltonian in this basis is achieved by simulating the dynamics of a walker population in Slater determinant space. FCIQMC mitigates the sign problem by walker anihilation to a certain degree but does not eliminate it [66]. By performing a stochastic projection to the ground state of a Hamiltonian in a given Fock basis directly, i.e. without resorting to a Hubbard-Stratonovich transformation, it is easy to respect symmetries of the Hamiltonian. In particular, it is possible to obtain energies and observables from yrast states by ground state projection in a plane-wave basis because the FCIQMC algorithm conserves total momentum if the Hamiltonian does. With the FCI method taking into account all correlations in the system, we probe the BCS-BEC crossover from the non-interacting to the strongly attractive regime in the Hubbard model.
In this paper, we present a QMC study of yrast states in the Hubbard model for spin-1 2 fermions. Using a filling factor much smaller than one, this system resembles a continuum superfluid with the difference that momentum is replaced by lattice momentum. We study the crossover from 1D to 2D geometry in the case of attractive on-site interactions. For 1D Hubbard chains, we obtain exact results using the Bethe ansatz and compare them to QMC results. Then we use FCIQMC to investigate the crossover into the 2D geometry by increasing the number of sites in the transverse direction. We find that in 2D, the shape of the yrast dispersion changes considerably, but that for increasing interaction strength, the more typical shape of the 1D spectrum is restored, which indicates soliton-like physics and is a signature of the transition into the superfluid regime. We investigate in more detail the behavior of the local minima of the yrast spectrum at the so-called "umklapp" points where sufficient quasi-momentum is added to boost either all or half of the constituent fermions by one unit in order to form a ring current. We find signatures of the transition from non-interacting Fermi gas to a paired superfluid and of the BCS to BEC crossover in the excitation energy and in the pair correlation functions for the first half and full umklapp points. Last, we focus on the yrast states around the maxima of the dispersion, which are related to dark solitons. We calculate the inertial mass of possible solitons for different system sizes and interaction strengths. We find an increase of the inertial mass by a factor of 2 when the transverse dimension is large enough for the system to be considered truly 2D, which is indica-tive of a transition from dark soliton to a solitonic vortex [67,68]. From mean-field and basic hydrodynamic theory, in the 1D to 2D crossover dark solitons are replaced as stable yrast excitations by solitonic vortices [68,69], or vortex pairs [41,70], which have larger inertial mass [34,71,72]. By looking at the pair densities of these yrast states, we find that fragmented condensation into more than one momentum state takes place, as expected for superfluid yrast states [28].
This paper is organized as follows: After introducing the model in Sec. II and the FCIQMC approach in Sec. III, we discuss yrast spectra of a 1D Hubbard chain obtained by the Bethe ansatz in Sec. IV. The energies and spin correlation functions for the umklapp points in the 1D to 2D crossover are discussed in Sec. V before analysing the physics of the maxima of the yrast dispersion by computing their effective mass and momentumspace pair densities in Sec. VI, and drawing conclusions in Sec. VII.

II. SYSTEM
To study yrast states in an ultracold fermionic superfluid, we use the Hubbard model in the regime of low densities. In this regime, the Hubbard model approximates a discretized free space. The correspondence to free space becomes exact in the low-density limit. We focus on the crossover between the 1D and 2D geometry, therefore our Hubbard model corresponds to a rectangular lattice with L × W sites, where 1 ≤ W < L, but with the same lattice spacing α in both dimensions. We are using periodic boundary conditions in both directions, and thus our systems has the topology of a torus. The Hubbard Hamiltonian in momentum representation is where the operators c ( †) k,s create (annihilate) a fermion with lattice momentum k and spin s. Wave vectors can take on values k = (k x , k y ) = 2π(n x /αL, n y /αW ), where n x , n y are integers. The single-particle dispersion for the Hubbard model is given by where t is the hopping amplitude and U < 0 is the interaction parameter. In the low density regime, mostly the low-lying momentum states are occupied where the dispersion relation (2) is approximately parabolic. Thus the Hubbard model approximates a continuum Fermi gas. Throughout this paper we present results obtained for a particle number of N = 10, with 5 fermions in each spin state, and a lattice length of L = 21, while the width W varies from 1 to 11. Energies will be given in units of hopping amplitude t and momenta in longitudinal lattice units P 0 = 2π αL .

III. FCIQMC
FCIQMC is a numerical method originally created in the context of strongly correlated electron systems and quantum chemistry. It can find ground states of fermionic many-body systems by expanding the many-body wave function in terms of Slater determinants which in our case are of the form FCIQMC then obtains the expansion coefficients C i by using a stochastic population dynamics approach to solve the imaginary-time Schrödinger equation. The automatic antisymmetrization of the wave function by expanding it in Slater determinants ensures that unlike Diffusion Monte-Carlo, FCIQMC always finds a fermionic wave function. FCIQMC's capacity for overcoming the so-called "fermion sign problem", which here manifests in fluctuations of the sign of each of the coefficients C i and which cannot be pre-determined, depends on the importance of annihilation events among the walkers in establishing the sign structure of the sampled wavefunction. When this effect is strong, the full FCIQMC method requires a number of walkers which scales with the size of the Hilbert space, making it impractical for large spaces.
To counter this, we use a range of modifications to FCIQMC, which facilitates calculations when the sign problem is strong. We make use of a similaritytransformed Hamiltonian which makes the many-body wave function more compact in Hilbert space [73]. We also use the initiator approximation [61], which can introduce a bias into the energy that disappears in the limit of large walker number. In order to control this undesirable bias, we first compare QMC results with exact results in the 1D case and then adjust the walker number until the initiator bias is eliminated. For the 2D systems, we successively increase the walker number with increasing W . Walker numbers used in this paper range from N W = 1 × 10 6 for weakly-interacting 1D chains to N W = 2×10 8 for 2D systems at U/t = −5. The most demanding computations were run on up to 400 processor cores using up to 2 Gigabyte memory per core.

IV. BETHE ANSATZ RESULTS
For a one-dimensional Hubbard chain (W = 1), the system is integrable and the Hamiltonian (1) can be diagonalized using the Bethe ansatz [74]. For a balanced Fermi system with N fermions, energy and total momen-  tum are given by with N dimensionless quasi-momenta κ j which must be obtained, alongside N/2 rapidities Λ α , by solving the Lieb-Wu equations Solving the Lieb-Wu equations via root finding can be done with great accuracy and polynomial effort with particle number. The Bethe ansatz thus provides us with an exact reference for the one-dimensional Hubbard chain. This allows us to compare QMC results with exact results to determine the parameter range where we can consider FCIQMC to be reliable in the sense that a possible initiator bias is smaller than the statistical uncertainty inherent in the Monte Carlo approach. In general, larger values of |U/t| lead to stronger correlations in the many-body wave function, which then requires a larger number of Slater determinants to be accurately represented. Also, the fermionic sign problem becomes more severe, which tends to increase the initiator bias in the calculated energy. In FIG. 1, we compare results obtained using the Bethe ansatz and FCIQMC results for a 1D chain (L = 21, W = 1). We see that for |U/t| ≤ 5 the agreement is very good. We therefore mainly use interaction strengths of |U/t| ≤ 5 in this paper, which covers the entire BEC-BCS crossover and typical values achievable in experiments [59].
V. UMKLAPP POINTS Figure 1 illustrates the characteristic shape of the yrast dispersion, which is concave downward resembling inverted parabolas in the intervals 0 < P < P 0 N/2 and P 0 N/2 < P < P 0 N , and with local minima at integer multiples of P 0 N/2. It can be understood by looking at the non-interacting case: To increase the total momentum of the system by a single unit of P 0 , first a fermion at the Fermi surface is excited. The resulting hole can then be filled by another fermion to increase the momentum again, and so forth, with the energy tracing the inverted parabolic part of the Hubbard lattice dispersion relation of Eq. (2). The first local minimum of the yrast spectrum, called the half-umklapp point, in our case (N = 10) at P/P 0 = 5 is reached when all particles of one spin component have each been boosted by P 0 . At that point, the Fermi surfaces of both components are shifted with respect to each other, but there are no holes in the Fermi seas of either spin species.
The full umklapp point at P/P 0 = 10 is reached when both spin components, or all particles are boosted. In the continuum limit, where full Galilean invariance is restored, the excitation energy of the umklapp point is determined by the boost only and is independent of interactions, since the state is strictly a boosted ground state. In the lattice system, where Galilean invariance is broken, a weak interaction dependence at the umklapp points is nevertheless observed as seen in Fig. 1.
For a non-interacting 2D system, constructing the yrast dispersion from hole excitations leads to a different shape, which in the thermodynamic limit in an isotropic 2D system is linear. In our case, as is shown in FIG. 2 for W = 11, due to finite-size effects in our mesoscopic system, the spectrum for U = 0 has linear segments but is not perfectly linear. It is remarkable that for increasing interaction strength, the parabolic shape of the 1D spectrum with the umklapp points at P/P 0 = 5, 10 is restored.
We can take a look at how the pairing in the system changes the characteristics of the umklapp points by calculating two-body correlation functions. It is possible to obtain the reduced two-body density matrix Γ s1,s2,s3,s4 ( k 1 , k 2 , k 3 , k 4 ) = Ψ 1 |c † k1,s1 c † k2,s2 c k3,s3 c k4,s4 |Ψ 2 (9) from FCIQMC by simultaneously running two statistically independent QMC simulations with solutions Opposite-spin pair correlation function g ↑↓ (k1, k2) = c † k 1 ,↑ c † k 2 ,↓ c k 2 ,↓ c k 1 ,↑ for the ground state (left) and the halfumklapp point (right) at U/t = −1 (top) and U/t = −5 (bottom). The crossover from a weakly interacting Fermi gas to a BEC of bound pairs is clearly visible, as the half-umklapp point changes its characteristics from one spin component shifted in momentum space with respect to the other to a translation of the whole system. The half-umklapp point of the Fermi system becomes the first full umklapp point of a fully paired superfluid with 5 bosonic pairs. Ψ 1 , Ψ 2 , to avoid biases [75]. This is valuable even in the 1D case as obtaining the same quantity from the Bethe-ansatz solution is not feasible. To illustrate the BEC-BCS crossover, we show the opposite-spin pair correlation function g ↑↓ ( k 1 , k 2 ) = Γ ↑↓↓↑ ( k 1 , k 2 , k 2 , k 1 ) for the ground state and half-umklapp point in the 1D case in FIG. 3. Strong pair correlations with k 1 +k 2 = 0 clearly emerge as interaction strength is increased from U/t = −1 to U/t = −5. The P/P 0 = 5 half-umklapp point at weak interactions exhibits mostly the physics of a noninteracting system, with two Fermi seas displaced with respect to each other. This is in stark contrast to the situation at U/t = −5, where the correlation function is the same as for P = 0 but shifted in momentum space. The mere translation of the pair correlation function is consistent with interpreting the system as superfluid of 5 bosonic pairs, where total momentum P/P 0 = 5 correponds to a full umklapp (i.e. Galilean boost of the ground state) in contrast to the weakly-interacting Fermi gas, which only reaches a half umklapp point at this momentum.
The yrast excitation energy at the half umklapp point P h = P 0 N/2 and at the full umklapp point P u = P 0 N are shown in Fig. 4. The lattice results can be compared to the expected excitation energies in an equivalent freespace system. The behavior of the umklapp points of an attractive 1D Fermi gas in free space has been studied using the Yang-Gaudin (YG) model [47]. There, the energy of the full umklapp point is independent of the interaction strength, as it represents simply a translation of the entire Fermi sea in momentum space. In the YG model this energy is given by E − E 0 = P 2 u /2N m, where N is the total particle number and P u = P 0 N the full umklapp momentum. In units of the Hubbard model parameters, this excitation energy is E − E 0 = 40π 2 t/L 2 and is depicted as the upper dashed horizontal line in Fig. 4 (a). The half-umklapp point however drops by factor of 2 from E − E 0 = P 2 h /2mN ↑ to E − E 0 = P 2 h /2mN , as it changes from being the half-umklapp point of a system of N fermions to being the full umklapp point of a gas of N/2 bosons. This is shown as the two lower dashed horizontal lines in Fig. 4 (a), with energies E−E 0 = 20π 2 t/L 2 and E − E 0 = 10π 2 t/L 2 .
For the non-interacting case we observe expected behavior with energy values slightly lower than the YG model. This is because the YG model uses a quadratic dispersion written in parameters of the Hubbard model as and k ≤ YG k . However, for finite interactions the energy values we calculate for the Hubbard model do not behave like for the YG model. In our system, lattice effects dominate once the interaction is strong enough. It is known that for U/t → −∞, the asymptotic effective Hamiltonian of the Hubbard model corresponds to a bosonic system with a one-boson-per-site hard-core condition and repulsive next-nearest neighbor interactions [76]. The effective Hamiltonian also has a global pre-factor of U 2 /t, meaning that the entire spectrum has the same scaling in the asymptotic regime. For the (half-)umklapp points, the asymptotes 80P 2 0 /m × t 2 /U and 20P 2 0 /m × t 2 /U are presented in Fig. 4 (a) as the green and purple dashed lines, respectively. We see that finite-size effects are reduced as the umklapp energies approach these asymptotic lines.

VI. MAXIMA OF THE YRAST SPECTRUM
The point P = P h /2 near the first local maximum of the yrast spectrum (such as at P/P 0 = 2, 3 in Fig. 2) is of particular interest as a point where in the 1D homogenous case dark solitons appear [27,47,77] that are stationary with respect to background and with phase step π across the soliton. Dark solitons in a Fermi superfluid are characterized by a localized density depression and a phase jump in the superfluid order parameter around this depression. In a system with periodic boundary conditions, this phase jump must be compensated by a phase gradient along the system, which corresponds to a constant counterflow velocity v cf . In addition, the soliton can be associated with an inertial mass m I , related to the curvature of the yrast dispersion. We extract these parameters from our calculated dispersions by fitting the quadratic function around the first local maximum at momenta P/P 0 = 1, 2, 3, 4, where P h = N P 0 /2 = 5P 0 .
FIG. 5. Inertial mass of the first yrast maximum vs. interaction strength. Unlike in free space, in the lattice, the inertial mass asymptotically approaches −∞ linearly with the interaction strength U/t. In the 2D regime, where a parabolic yrast dispersion reappears for strong interactions, the inertial mass is larger by a factor of 2.
The results for the inertial mass are shown in Fig. 5, where we only show data points for the parameters where the yrast dispersion closely resembles a parabolic shape. This would correspond to a regime where a superfluid is present and the particles are strongly paired. This is mostly the case for W ≤ 7, where the system is still effectively one-dimensional and transverse momentum states are sparsely populated. For the larger and more 2D systems with W = 7, 9, we find that for |U/t| = 4, 5, a parabolic yrast spectrum reappears (see also Fig. 2). For the effective mass, there is an increase in magnitude by a factor of 2 as we increase the system size to W = 11. This indicates a further change to the properties of the system, likely a transition from a soliton state to a vortex pair. This scenario is closely related to the snaking instability of a planar dark soliton in a two-dimensional superfluid, where the soliton decays into pairs of oppositely charged vortices as the system becomes wide enough [41,68,70]. However, we cannot directly show a potential density depression caused by the soliton. The reason is that, unlike in mean-field theory, our technique provides the many-body wave function of a (translationally invariant) eigenstate of total momentum and thus a superposition of solitons at all possible positions. The real-space singleparticle density we can calculate is flat. To map out the dark soliton as described in [77], we would need access to higher-order density matrices beyond the two-body density matrix. Therefore, it remains to be seen if the increase in inertial mass really corresponds to a solitonvortex pair transition. Instead, we investigate more closely the pair condensation, for which a relevant quantity is the pair Green's function where ψ and ψ † denote creation and annihilation operators, respectively, in position space. The Fourier transform of this Green's function is the momentum-space pair density. For a system with a homogeneous density, it can be directly obtained from the momentum representation of the two-body density matrix This quantity indicates whether Bose-Einstein condensation of pairs occurs. While we find evidence of Bose-Einstein condensation of pairs by a peak in the pair density that grows with increasing interaction strength for the ground state and umklapp point, the situation is more complex for general yrast states. In Fig. 6 we show the pair density for a 2D system with 21 × 11 sites for interaction parameters U/t = −1 and U/t = −5, for the ground state and the P/P 0 = 2 yrast state. For weak interactions, the structure of the pair density is determined mostly by the structure of the Fermi sea, or the noninteracting yrast state. For stronger interactions, the ground state exhibits one sharp peak at zero momentum indicating strong pairing correlations, as expected for crossover to a BEC of pairs. However, for P/P 0 = 2, there are actually two peaks, for longitudinal momenta 0 and 2P 0 . Similar features appear in the smaller 2D and 1D systems. We now show that this double-peak feature signifies the presence of fragmented condensation. Fragmentation occurs when during the transition to a Bose-Einstein condensate, more than one state becomes macroscopically occupied [78]. In Fig. 6, the momentum states (0, 0) and (2P 0 , 0) dominate the pair density. For the 1D system, where obtaining the reduced density matrices is easier, we plot the pair densities of several momenta for the yrast state with P = 3P 0 in Fig. 7 as a function of interaction strength. We see that in addition to P = 0, the pair density at P = P 0 strongly increases as well. Similarly, for other yrast states we see the same phenomenon, a strong signature of fragmented condensation. It is worth noting that the pair density with two peaks obtained here is similar to the case of an imbalanced Fermi gas, where Fermi surfaces of different size lead to FFLO pairing with non-zero total momentum and the signature is a twopeaked pair-density. This has been studied for 1D and 2D Hubbard models [44,57,58,79,80]. Yrast states in our balanced system start with holes in one of the Fermi seas for weak interaction, which also leads to pairing with non-zero total momentum.

VII. CONCLUSIONS
In this paper, we have used the FCIQMC method to study the crossover from weakly interacting fermions to a condensate of bosonic pairs for yrast states in mesoscopic Fermi systems. With this method we can treat larger systems than are accessible to the previously used deterministic CI or exact diagonalization methods [48,49] and can probe the full transition from a 1D chain not just to quasi-1D rings, but also to full 2D systems. We obtain energy spectra and reduced two-body density matrices The ground state shows a rapid growth of the pair density at zero momentum, as expected for Bose condensation of pairs, while the half-umklapp point is identical to the ground state but shifted by one lattice momentum unit. For the yrast states however, we observe near equal growth of both pair momenta P = 0 and P = P0, meaning that fragmented condensation into zero and non-zero momentum states is taking place.
for yrast states in the attractive Hubbard model in these geometries.
Comparisons with exact Bethe ansatz results in 1D show very good agreement and demonstrate that FCIQMC can accurately provide yrast states for mesoscopic systems with 10 fermions and 21 × W sites where W ranges from 1 to 11. We find that the shape of the yrast spectrum changes from the typical inverted parabolas in 1D to a quasi-linear spectrum in 2D. However, as interaction strength is increased, the parabolic dispersion including the half umklapp point is restored. While the quasi-linear spectrum is a consequence of the 2D geometry and the Pauli exclusion principle, with stronger interactions the exact shape of the non-interacting Fermi sea plays less of a role until we see the expected universal concave downward dispersion of a spinless superfluid. This indicates a transition to a fully paired Fermi superfluid.
We further find that mesoscopic effects can cause significant deviations for certain geometries from the general behavior of the umklapp energies, which otherwise does not differ much between 1D and 2D systems. Specifically we find that the half-umklapp excitation energy for a lattice of 21 × 11 points increases with interaction strength at intermediate values of U/t contrary to the general trend displayed by all other systems under study. This originates in a rearrangement of the fermions in momentum space, where for this particular geometry, the half-umklapp point is a different configuration than the ground state with one spin component shifted in momen-tum space. This can be of importance for experiments on mesoscopic Fermi systems.
By calculating the pair correlation function, we can follow the pairing process by means of which the fermionic half-umklapp point becomes the first full umklapp point of the Bose condensate of pairs.
In the fully paired regime, we calculated the inertial mass of the dark-soliton-like local maximum of the yrast dispersion. We found a sudden increase by a factor of 2 between the narrow stripe geometry and our largest system in a lattice of 21 × 11 sites. We interpret this as a possible change in the system geometry where a soliton is no longer stable and the yrast state is instead provided by pair of oppositely charged vortices.
The most striking feature of yrast states in the Hubbard model is revealed to be fragmented condensation, which occurs around the maxima of the yrast dispersion, away from the umklapp points. We find that here, multipeaked pair densities appear, where in addition to pairing with zero total momentum, the amplitude of other total momentum pairs becomes large for strong interactions. These pair densities of yrast states share some similarities with FFLO states, which are characterized by a doublepeaked momentum pair density [44]. In the FFLO case the origin of this is a mismatch of Fermi surfaces, which have different size due to the spin imbalance. In our case of the yrast states, the Fermi surfaces are shifted in momentum space. In addition, both yrast and FFLO states are related to solitons in the real space density [44].