When does a Fermi puddle become a Fermi sea? Emergence of Pairing in Two-Dimensional Trapped Mesoscopic Fermi Gases

We thank the referee for their useful and insightful feedback, which we believe has helped us to improve our manuscript. Below, please find our responses to the specific points raised. We number the referee’s comments from top to bottom as 1, 2, 3. In our resubmitted manuscript, we have written the associated changes in blue-coloured text to make them immediately obvious.

Point 1 (top):

We agree with the referee that it would be helpful to show the ground- and excited-state energies used to calculate the excitation energy separately. In Fig. 2 on page 8, we have created a new middle panel (b) where the ground- and excited-state energies from panel (a) [i.e., for 1+1, 2+2, and 3+3 fermions at zero effective range] are plotted separately as a function of the two-body binding energy. This makes it clear how the different dependencies of the excitation energies on the interaction strength arise in panel (a). To accommodate this change, we have reworked the caption of Fig. 2 on page 8 and edited/added two sentences in the first paragraph of Sec. 3.1 "Excitation Spectrum" on page 7 [lines 218–221 and 224–226].

Point 2 (middle):

We thank the referee for pointing out our lack of clarity regarding the "number of pairs" which saturates at 3 (for 3+3 fermions) in the strong-interaction limit, versus the "paired fraction" which saturates at 1; we always mean the number of pairs. To correct this, we have renamed Sec. 3.3 from "Paired Fraction" to "Number of Pairs" and replaced the phrase "paired fraction" in the Abstract with "number of pairs". We have also added a new sentence to the caption of Fig. 6 to remove any remaining ambiguity: "The maximum possible number of pairs is N_pair = 3.". (In Sec. 3.3 the text itself is clear on the fact that N_pair is the number of pairs with a maximal value of 3.)

The referee is correct to expect that for much stronger interactions than those shown, epsilon_b > 2 h-bar omega_r, all the fermions would form tightly bound bosonic dimers, reminiscent of the deep BEC regime in macroscopic systems. We have added a new sentence to this effect in Sec. 3.3 [lines 413–416]. (We had also mentioned this earlier on page 12 [footnote 8] where we first discuss our considered range of binding energies.) The reason why Fig. 6 does not show any red data points for epsilon_b > 2 h-bar omega_r is that for such strong binding energies, it is challenging to properly model the tight composite bosonic wave functions, and thus, to obtain fully numerically converged energies and structural properties within a reasonable time frame. We now explain this point in a new footnote [11] in Sec. 3.3.

Point 3 (bottom):

Indeed, our mention of the many-body limit was vague and hence confusing. We are very grateful to the referee for bringing this important point to our attention. In free space where BCS theory applies, the many-body limit is typically approached by increasing both the number of particles, N -> infinity, and the system volume, V -> infinity, in such a way that the density n = N/V remains constant. For our scenario where omega_r is fixed, we instead have n -> infinity when N -> infinity. In this case, the many-body two-component Fermi gas remains in the normal state for epsilon_b << h-bar omega_r and undergoes a quantum phase transition to a superfluid state at epsilon_b^c. If we wish to make our system amenable to BCS theory, then we could keep n constant by reducing the trap strength omega_r while increasing N. In that case, the trapping frequency would vanish (omega_r -> 0) for N -> infinity, which means the condition epsilon_b << h-bar omega_r would never be satisfied in the many-body limit. Consequently, a superfluid state with a finite gap would always exist at zero temperature for any non-vanishing interaction strength. We have clarified this point about the many-body limit in both sections of the text highlighted by the referee: in Sec. 3.1 "Excitation Spectrum" on page 9 [lines 263–269 and footnote 7] and in Sec. 3.3 "Number of Pairs" [lines 423–428].

We are very grateful for the referee’s careful reading of the manuscript and their insightful feedback. Below, please find our responses to the specific points raised. The numbers of our responses coincide with the numbers in the referee’s report. In our resubmitted manuscript, we have written the associated changes in orange-coloured text to make them immediately obvious.

Point 1:

The referee is correct. According to conventional standards, a "large" effective range (r_2D) would be one which satisfies: h-bar^2/(2m|r_2D|) > h-bar omega_r (noting that our definition of r_2D has units of squared length and is negative in value). This inequality is not satisfied if we substitute in the experimental values of {m, omega_r, r_2D} provided in Ref. [24], and thus, while the effective range in our study is finite and non-negligible, it is not conventionally "large". In consequence, we have altered the concerned sentence to read [page 6, lines 189−190]: "In practical computations, we tune the effective range to non-negligible negative values through a shape resonance...".

Point 2:

We agree with the referee that introducing the value of the critical binding energy earlier (than in Sec. 3.3) would enhance the coherency and comprehension of the text. In order to introduce the concept and value of the critical binding energy in Sec. 3.1, we have altered 5 sentences on pages 9−10 which now read as follows:

Lines 252−253: "At a critical binding energy (denoted by epsilon_b^c) the excitation energy Delta-E reaches a minimum."

Lines 264−267: "In this case, if the ground state has a closed-shell configuration, then the minimum value of Delta-E at the critical binding energy epsilon_b^c will decrease as N increases, eventually reducing to zero in the many-body limit so that pairs are coherently excited without any energy cost [25, 26]."

Lines 267−269: "In this limit, if epsilon_b is increased from zero to epsilon_b^c, then the many-body two-component Fermi gas will become unstable and undergo a second-order phase transition into a superfluid state."

Lines 277−279: "Notably, we find that as the magnitude of the negative effective range increases, the minimum value of Delta-E decreases and shifts to smaller binding energies, i.e., epsilon_b^c decreases."

Lines 284−285: "The value of the critical binding energy for this line is epsilon_b^c = 0.953 h-bar omega_r."

Point 3:

We thank the referee very much for their suggestion. We have conducted an initial analysis in this direction and obtained some preliminary results. In two dimensions, the two-body density matrix in Eq. (6), rho(r_1, r’_1; r_2, r’_2), is an eight-dimensional array. Therefore, we need to reduce the number of degrees of freedom before diagonalising it to obtain its eigenvalues (occupation numbers). In the first instance, we have done this by following Ref. [39] which treated a very similar problem in three dimensions. The idea is to first transform from the co-ordinates of the individual particles to the centre-of-mass (R, R’) and relative (r, r’) co-ordinates of the two spin-up-spin-down pairs, and to then set the relative distance vectors equal to each other, r = r’. Setting r = r’ is somewhat artificial; however, one then obtains a reduced two-body density matrix, rho_red(R, R’), which has the same number of degrees of freedom as the one-body density matrix, and thus, is diagonalisable once it has been decomposed into partial-wave contributions. Due to setting r = r’, the largest eigenvalue of rho_red(R, R’) does not directly indicate the onset of superfluidity. Nevertheless, the "condensate fraction" may be estimated by comparing the magnitudes of the eigenvalues of rho_red(R, R’) — e.g., see Eq. (16) in Ref. [39]. By analysing the reduced two-body density matrix in two dimensions, we have obtained plots of the occupation numbers and "condensate fraction" as a function of the two-body binding energy, akin to Fig. 11 in Ref. [39]. However, since rho_red(R, R’) only measures non-local correlations between spin-up-spin-down pairs that are characterised by the same relative distance vector, we believe this approach should be more carefully examined and compared with alternative methods. Hence, while the referee has given us a great idea, we have reached the conclusion that it is beyond the scope of the current work. We believe it deserves a more thorough and detailed investigation which should form the focus of a future paper. In the current paper, on the other hand, we are instead more focused on observables that are measurable in experiments.