How do pairing order and condensate fraction behave across the BCS-BEC crossover?

A benchmark map of the attractive Hubbard ground state

The paper asks how the superconducting ground state of the two-dimensional attractive Hubbard model evolves as pairing moves from weakly bound Cooper pairs toward tightly bound bosonic pairs. Its answer is a systematic, numerically exact set of ground-state data rather than a single number: the authors characterize the interaction-driven and density-induced BCS-BEC crossovers using pair wave functions and the fermion momentum distribution, and report high-precision results for the on-site spin-singlet pairing order, the condensate fraction, the double occupancy, and the total energy density.

The parameter window is broad and explicitly stated: interaction strength U/t = 2 to 12, where t is the nearest-neighbor hopping amplitude, and fermion density n = 0.25 to 1.0. System sizes run from L = 4 to 24, and the authors perform finite-size extrapolations to obtain thermodynamic-limit values. The stated purpose is to establish quantitative benchmarks for understanding superconductivity across the crossover and for validating future theoretical and numerical approaches.

Why auxiliary-field quantum Monte Carlo anchors the numbers

The results come from the numerically exact zero-temperature auxiliary-field quantum Monte Carlo method, also called projector quantum Monte Carlo. The ground state is reached by imaginary-time projection from a trial wave function, with the projection parameter Θ controlling convergence; observables are evaluated as expectation values in the projected state. The authors note that their simulations are free of the minus sign problem for arbitrary fermion density, which is what allows large-scale calculations and reliable finite-size extrapolation to the thermodynamic limit.

Several implementation choices matter for interpreting the numbers. The trial wave function is the ground-state Slater determinant of the kinetic part, with a tiny twisted boundary phase added only when constructing it to lift degeneracies at the Fermi level. A third-order Trotter-Suzuki decomposition is used, and residual Trotter error is removed by convergence tests in the imaginary-time step Δτt until results are statistically indistinguishable within uncertainties; similar tests are done for the projection parameter Θt. Measurements are taken over an imaginary-time interval 0.45Θ ≤ τ ≤ 0.55Θ rather than at a single midpoint, which the authors justify by the remaining projection time on both sides being long enough to ensure ground-state convergence.

Pairing order and condensate fraction are not the same quantity

The two headline quantities answer different questions and should not be substituted for one another. The pairing order is characterized through the mean-squared pairing order ⟨Δ²⟩ and the pair-pair correlation function, including the value at the separation r = (L/2, L/2); the condensate fraction n_c is obtained from the momentum-space pairing matrix, which the fast Fourier transform implementation gives direct access to. One describes the strength and spatial structure of pairing correlations, the other the occupation of the zero-momentum pairing channel.

This distinction shapes how the paper's benchmark data should be read. The authors report both quantities across the same interaction and density ranges, and both are extrapolated to the thermodynamic limit, so a reader comparing theories or numerics against this dataset needs to match the observable, not just the parameter point. The paper also provides high-precision double occupancy and total energy density as complementary quantities, which probe local interaction energy and overall energetics rather than pairing order.

Finite-size corrections decide the thermodynamic-limit values

Extrapolation is not a formality here. The authors observe oscillations in ⟨Δ²⟩, P[r = (L/2, L/2)], and n_c for systems with L < 12, and a simple linear fit to periodic-boundary-condition data with L ≥ 12 yields evidently inconsistent thermodynamic-limit results for ⟨Δ²⟩ and P[r = (L/2, L/2)]. To resolve this, they run simulations under twist-averaged boundary conditions and apply an additional finite-size correction, using the noninteracting system as the exactly solvable reference and averaging over twist angles.

With that combined procedure, the corrected ⟨Δ²⟩ and P[r = (L/2, L/2)] extrapolate to consistent thermodynamic-limit values, and the same calculation for n_c also reaches reliable thermodynamic-limit data. The authors apply this combination mainly at U/t = 2 with n = 0.25 to 0.875, primarily for ⟨Δ²⟩ and n_c. For double occupancy and total energy density, twist-averaged boundary conditions typically give well-converged results up to L = 24, so finite-size extrapolation is unnecessary there. This means the correction burden is not uniform across observables or interaction strengths.

Next-nearest-neighbor hopping adds a second crossover

Beyond the standard model, the authors investigate the effects of next-nearest-neighbor hopping t₂ on the superconducting properties. In its presence they identify a distinct t₂-induced BCS-BEC-like crossover and explore how the energetics, pairing order, and condensate fraction evolve with t₂ at U/t = 6 and n = 0.875.

This is a separate scenario from the interaction-driven and density-induced crossovers mapped in the standard model, and it is reported at a single interaction strength and density rather than across the full U/t = 2 to 12 and n = 0.25 to 1.0 grid. Readers should therefore treat the t₂ results as a targeted probe of band-structure effects, not as a second complete phase map.

What these benchmarks can and cannot settle

The paper's own framing is that the qualitative structure of the ground state was already relatively well established, while a comprehensive quantitative characterization remained limited. The contribution here is therefore the quantitative dataset: unbiased, high-precision values covering the pair-pair correlation function, mean-squared pairing order, condensate fraction, double occupancy, and energy densities, which the authors describe as an extensive numerical dictionary for the two-dimensional attractive Hubbard model. They state that these data can serve as benchmarks for various theoretical and numerical approaches and can offer useful guidance for interpreting computed superfluid transition temperatures.

The scope conditions are set by the model and method. This is the two-dimensional attractive Hubbard model on the stated lattice with on-site attraction, solved at zero temperature for the ground state; the reported quantities are ground-state pairing order, condensate fraction, double occupancy, and energy density, not the superconducting transition temperature itself. The t₂ investigation is confined to U/t = 6 and n = 0.875. The authors also position the work as a foundation for further investigations of more complex systems, such as the attractive Hubbard model with spin-orbit coupling and on composite lattices, which are outside the present evidence.

Start with the observable, then the parameter point

For anyone using this dataset, the decisive step is to fix which quantity is being compared before fixing the parameter point. Pairing order and condensate fraction are measured differently, extrapolated under different correction needs, and answer different physical questions, so a mismatch on the observable invalidates a comparison even when U/t and n agree. The second decisive step is the interaction strength: the twist-averaged boundary conditions plus finite-size correction are applied mainly at U/t = 2 for ⟨Δ²⟩ and n_c, while double occupancy and total energy density converge without extrapolation up to L = 24.

The most useful priority this paper establishes is therefore a comparison discipline rather than a single headline number: match the observable, match the boundary-condition and finite-size treatment, and only then compare against the thermodynamic-limit values. The t₂-induced crossover at U/t = 6 and n = 0.875 is the clearest case where a reader should not assume the standard-model behavior carries over, because it is reported as a distinct crossover driven by band structure rather than by interaction strength or density alone.

Sources used in this answer