# Nature of the mixed-parity pairing of attractive fermions with spin-orbit
coupling in optical lattice

Yu-Feng Song  Hefei National Laboratory for Physical Sciences at Microscale
and Department of Modern Physics, University of Science and Technology of
China, Hefei, Anhui 230026, China  Institute of Modern Physics, Northwest
University, Xi’an 710127, China  Youjin Deng<EMAIL_ADDRESS>Hefei
National Laboratory for Physical Sciences at Microscale and Department of
Modern Physics, University of Science and Technology of China, Hefei, Anhui
230026, China  Hefei National Laboratory, University of Science and Technology
of China, Hefei 230088, China  Yuan-Yao He<EMAIL_ADDRESS>Institute of
Modern Physics, Northwest University, Xi’an 710127, China  Shaanxi Key
Laboratory for Theoretical Physics Frontiers, Xi’an 710127, China  Hefei
National Laboratory, University of Science and Technology of China, Hefei
230088, China

###### Abstract

The admixture of spin-singlet and spin-triplet pairing states in
superconductors can be typically induced by breaking spatial inversion
symmetry. Employing the numerically exact auxiliary-field Quantum Monte Carlo
method, we study such mixed-parity pairing phenomena of attractive fermions
with Rashba spin-orbit coupling (SOC) in two-dimensional optical lattice at
finite temperature. We systematically explore the evolution of the essential
pairing structure in both singlet and triplet channels versus the temperature,
fermion filling, SOC and interaction strengths, via computing the finite-size
results of condensate fraction and pair wave function. Our numerical results
reveal that the singlet channel dominates in the fermion pairing and the
triplet pairing has relatively small contribution to the superfluidity for
physically relevant parameters. In contrast to the singlet channel mainly
consisted of the on-site Cooper pairs, the triplet pairing has plentiful
patterns in real space with the largest contributions from several nearest
neighbors. As the SOC strengh increases, the pairing correlation is firstly
enhanced and then suppressed for triplet pairing while it’s simply weakened in
singlet channel. We have also obtained the Berezinskii-Kosterlitz-Thouless
transition temperatures through the finite-size analysis of condensate
fraction. Our results can serve as quantitative guide for future optical
lattice experiments as well as accurate benchmarks for theories and other
numerical methods.

##  I Introduction

The fermion paring and corresponding superconductivity and superfluidity [1]
are of great interest in condensed matter physics. The fundamental ingredient
is the Cooper pair consisting of two spin-$1/2$ electrons [2]. Given the
spatial inversion symmetry, the pair wave function can be decoupled into
orbital and spin channels resulting in two states of Cooper pairs, even parity
with spin-singlet and odd parity with spin-triplet [3]. Majority of known
superconductors (SCs) fall into the spin-singlet case, such as simple metals
[4] and high-$T_{c}$ cuprates [5, *Lee2006, *Armitage2010]. Nevertheless, the
triplet paring has been observed or suggested to exist in far fewer realistic
systems, e.g., superfluid 3He [8], $\text{UPt}_{3}$ and
$\text{Sr}_{2}\text{RuO}_{4}$ [9]. Without inversion symmetry, the parity
conservation is broken and thus the mixing of singlet and triplet paring
states can emerge [10, 11, 12, 13]. Such mixed-parity pairing state has been
experimentally verified in various three-dimensional (3D) noncentrosymmetric
SCs [14, 15, 16, 17, 18, 19], which induces intensive interests due to many
exotic properties [12] including fertile superconducting gap structures [15,
20], anisotropic magnetic response [21, 22] and topological superconductivity
[23].

The appearance of the mixed-parity pairing in noncentrosymmetric systems can
be attributed to the arise of the antisymmetric spin-orbit coupling (SOC)
[24], which has become one of the key elements for condensed matter physics
[25]. For correlated fermion systems, SOC acts as another dimension and
induces many exotic states of matter, including spintronics [26], topological
phases [27, *Qi2011] and unusual superconductivity [13]. Specifically, it
typically breaks the spatial inversion symmetry and mixes the spin species,
and thus renders the coexistence of spin-singlet even-parity and spin-triplet
odd-parity pairing. Moreover, it was shown [16] that tunning the SOC strength
can even change the dominant component of the mixed-pairity pairing from the
singlet in $\text{Li}_{2}\text{Pd}_{3}\text{B}$ to the triplet in
$\text{Li}_{2}\text{Pt}_{3}\text{B}$, as replacing the Pd atom by Pt atom. To
date, most of the study for the SOC induced singlet-triplet mixed pairing
phenomena concentrates on the 3D systems including the noncentrosymmetric SCs
[13] and interacting Fermi gas [29, 30, 31].

In physically more relevant two-dimensional (2D) systems, the interplay
between the reduced dimensionality and enhanced quantum fluctuations can
induce fascinating and unique quantum phenomena [32, 33, *Giorgini2008, 35,
*Kotov2012]. A typical representative is the Berezinskii-Kosterlitz-Thouless
(BKT) transition [37, 38, 39, 40] of superconductivity and superfluidity.
Similar to the 3D analog, inclusion of SOC to 2D attractive fermions can also
induce mixed-pairity pairing [10], which has been relatively much less
studied. Experimentally, the recently elegant realization of synthetic SOC for
fermions [41, 42] with ultracold atoms, especially in 2D optical lattice [43,
44, 45], substantively paves the way for exploring novel quantum phenomena
related with SOC. Thus, a systematically theoretical study with high precision
on the mixed-pairity pairing in 2D is highly demanded to shed light on
problems closely related to ultracold atom experiments. For example, finding
the best condition to observe the spin-triplet pairing in optical lattice, in
comparison to the achieved singlet pairing [46], should be a useful guide for
experiments.

To date, most theoretical work on 2D attractive fermions with SOC falls into
the Fermi gas and approximate theories [47, *Anna2012, 49, 50, 51, 52].
Interesting results such as singlet and triplet contributions to the
condensation [47, *Anna2012] are presented in these studies, but still need
careful verifications from unbiased approaches. Nevertheless, numerically
exact calculations for such systems are rare. Auxiliary-field Quantum Monte
Carlo (AFQMC) simulations have been performed for the ground state of 2D Fermi
gas [53], and for the lattice system at finite temperatures [54] as well as
its ground state [55]. The authors in Ref. [54] focused on the properties of
BKT transition temperatures and anisotropic spin susceptibility without
touching the pairing structure, which were limited to 12$\times$12 finite
lattices. The pairing structure were discussed in Ref. [55] only for the half-
filling case, for which the BKT transition disappears and thus it was of less
interest to experiments.

In this paper, we study the mixed-parity pairing of attractive fermions with
Rashba SOC in 2D optical lattice, applying finite-temperature AFQMC algorithm
[56, 57, 58, 59]. We mainly concentrate on the condensate fraction and pair
wave functions to reveal the pairing structures of both singlet and triplet
channels for physically relevant regimes of the temperature, fermion filling,
SOC and interaction strengths. We also present the determination of BKT
transition temperature from the finite-size analysis of condensate fraction
results. The rest of the paper is organized as follows. In Sec. II, we
introduce the lattice model that we use to describe the 2D attractive fermions
with Rashba SOC in optical lattice, and the AFQMC method. In Sec. III, we
present our numerical results, including the pairing structures, the pairing
correlations and calculations of the BKT transition temperature. Finally, Sec.
IV summarizes this work, and discusses its connections with the optical
lattice experiments.

##  II Model and method

We describe the 2D attractive fermions with Rashba SOC in optical lattice
using the following square lattice model Hamiltonian [54, 55] as

$\displaystyle\hat{H}=$
$\displaystyle\sum_{\mathbf{k}\sigma}\varepsilon_{\mathbf{k}}c_{\mathbf{k}\sigma}^{{\dagger}}c_{\mathbf{k}\sigma}+\sum_{\mathbf{k}}(\mathcal{L}_{\mathbf{k}}c_{\mathbf{k}\downarrow}^{{\dagger}}c_{\mathbf{k}\uparrow}+\mathrm{H.c.})$
(1) $\displaystyle\hskip
8.5359pt+U\sum_{\mathbf{i}}\Big{(}\hat{n}_{\mathbf{i}\uparrow}\hat{n}_{\mathbf{i}\downarrow}-\frac{\hat{n}_{\mathbf{i}\uparrow}+\hat{n}_{\mathbf{i}\downarrow}}{2}\Big{)}+\mu\sum_{\mathbf{i}\sigma}\hat{n}_{\mathbf{i}\sigma},$

with $\varepsilon_{\mathbf{k}}=-2t(\cos k_{x}+\cos k_{y})$,
$\mathcal{L}_{\mathbf{k}}=2\lambda(\sin k_{y}-i\sin k_{x})$, and
$\hat{n}_{\mathbf{i}\sigma}=c_{\mathbf{i}\sigma}^{{\dagger}}c_{\mathbf{i}\sigma}$
representing the density operator with spin $\sigma=\uparrow,\downarrow$ on
the lattice site $\mathbf{i}=(i_{x},i_{y})$. The momentum $k_{x}$ and $k_{y}$
are defined in units of $2\pi/L$ with the system size $N_{s}=L^{2}$. We denote
the fermion filling as $n=N/N_{s}$ with $N$ as the total number of fermions in
the system. The nearest-neighbor hopping $t$, on-site Coulomb interaction $U$
($<0$), the SOC strength $\lambda$, and chemical potential $\mu$ are model
parameters. Within the above formulation, the system is at half filling with
$n=1$ for $\mu=0$ due to the particle-hole symmetry [54], and it is hole doped
for $\mu>0$. Throughout this work, we set $t$ as the energy scale, and we
focus mostly on the doped systems with the fermion filling $n<1$.

The previous study [55] showed that the model in Eq. (1) has a supersolid
ground state with coexisting charge and superfluid long-range orders at half
filling. Away from this special point, the superfluidity survives for
arbitrary filling with arbitrary interaction strength [53]. Since the SOC term
breaks the spin SU(2) symmetry and results in two helical bands for
noninteracting case, the corresponding superfluid state with interaction is
composed of both spin-singlet and triplet Cooper pairs, whose pairing
properties are the main content of this work.

We then apply the finite-temperature AFQMC algorithm [56, 57, 58, 59] to
numerically solve the lattice model in Eq. (1). It is free of fermion sign
problem at arbitrary filling due to the time-reversal symmetry [60]. The
scheme of the AFQMC method is first to decouple the two-body interactions into
free fermions coupled with auxiliary fields and then to calculate the
fermionic observables through importance sampling of the field configurations.
Practically, the imaginary-time discretization of the inverse temperature as
$\beta=M\Delta\tau$, the symmetric Trotter-Suzuki decomposition
$e^{-\Delta\tau\hat{H}}=e^{-\Delta\tau\hat{H}_{0}/2}e^{-\Delta\tau\hat{H}_{I}}e^{-\Delta\tau\hat{H}_{0}/2}+\mathcal{O}[(\Delta\tau)^{3}]$
(with $\hat{H}_{0}$ and $\hat{H}_{I}$ as the free and interaction parts of the
Hamiltonian), and the Hubbard-Stratonovich (HS) transformation are
successively implemented. The discrete HS transformation with the
spin-$\hat{s}_{z}$ decomposition rather than the usual charge channel [57] is
adopted for the attractive $U$ interaction to suppress the fluctuations of
pairing related observables. Other algorithmic advances and techniques applied
here include the Fast Fourier Transform (FFT) between the real and momentum
space [59], the delayed version of local update [61], and the $\tau$-line type
of global update [62], which together improve the efficiency of the numerical
simulations. For further details of the AFQMC algorithm, we refer to the
reviews in Ref. [63, 64].

##  III Numerical results

In this section, we present the AFQMC simulation results of the lattice model
in Eq. (1), including the pairing structure, the pairing correlation functions
and BKT transition. Our AFQMC calculations reach the linear system size $L=20$
with the temperature as low as $T/t=0.025$ to sufficiently access the
superfluidity (quasi-long-range ordered or quasi-condensate) regime. We mainly
concentrate on the pairing properties away from half filling, for which the
charge density wave does not have long-range order (see details in Appendix
A). The parameter $\Delta\tau t=0.05$ is chosen mostly in this work, which has
been tested to safely eliminate the Trotter error, except for the strong
interactions where smaller $\Delta\tau$ is applied. Periodic boundary
conditions in both directions are applied for all the calculations.

###  III.1 Condensate fractions and pair wave functions

The contributions of spin-singlet and triplet channels to the fermion pairing
can be quantified by the corresponding condensate fractions [47, *Anna2012].
On the other hand, properties of the Cooper pairs, including their sizes and
the fermion momentum, can be obtained from the pair wave functions [53, 55].
The computation of these quantities involves the following pairing matrix in
momentum space [53, 55] as

$\displaystyle
M(\mathbf{k},\ell;\mathbf{k}^{\prime},\ell^{\prime})=\langle\Delta_{\ell}^{{\dagger}}(\mathbf{k})\Delta_{\ell^{\prime}}(\mathbf{k}^{\prime})\rangle,$
(2)

with $\ell=s$ or $t_{\uparrow}$ or $t_{\downarrow}$, and
$\Delta_{\ell}^{{\dagger}}(\mathbf{k})$ as spin-singlet and triplet pairing
operators with zero center-of-mass momentum as

$\displaystyle\hskip
17.07182pt\Delta_{s}^{{\dagger}}(\mathbf{k})=\frac{1}{\sqrt{2}}(c_{\mathbf{k}\uparrow}^{{\dagger}}c_{-\mathbf{k}\downarrow}^{{\dagger}}-c_{\mathbf{k}\downarrow}^{{\dagger}}c_{-\mathbf{k}\uparrow}^{{\dagger}})$
(3)
$\displaystyle\Delta_{t_{\uparrow}}^{{\dagger}}(\mathbf{k})=c_{\mathbf{k}\uparrow}^{{\dagger}}c_{-\mathbf{k}\uparrow}^{{\dagger}}\hskip
22.76228pt\Delta_{t_{\downarrow}}^{{\dagger}}(\mathbf{k})=c_{\mathbf{k}\downarrow}^{{\dagger}}c_{-\mathbf{k}\downarrow}^{{\dagger}}.$

Note that the third component of the triplet pairing eliminates by symmetry.
In $N_{s}=L^{2}$ finite system, the pairing matrix $\bf{M}$ is a $3N_{s}\times
3N_{s}$ matrix. Attributed to the FFT algorithm applied in our numerical
calculations, we can compute the equal-time, momentum-space single-particle
Green’s function matrix
$\mathbf{G}=\\{G_{\mathbf{k}\sigma,\mathbf{k}^{\prime}\sigma^{\prime}}=\langle
c_{\mathbf{k}\sigma}c_{\mathbf{k}^{\prime}\sigma^{\prime}}^{\dagger}\rangle_{\tau}\\}$,
and thus directly measure the pairing matrix defined in Eq. (2) for a single
auxiliary-field configuration via Wick decomposition.

Figure 1:  Illustration of the condensate fractions and pair wave functions in
momentum space. Plotted in panel (a) are the total, singlet and triplet
condensate fractions versus the temperature. The inset is the corresponding
fermion filling. Panel (b) presents the magnitudes of momentum-space pair wave
functions $|\Psi_{c}(\mathbf{k})|$ in both singlet and triplet channels for
two temperatures $T/t=0.025$ and $0.30$. The inset plots the Fermi surfaces
(red and blue dashed lines) for the fermion filling of $T/t=0.025$ case and
the high-symmetry path (black solid lines with arrows) with
$\boldsymbol{\Gamma}$, $\mathbf{X}$ and $\mathbf{M}$ points in Brillouin zone.
These calculations are performed for $L=20$ system with $U/t=-4$,
$\lambda/t=0.5$ and $\mu/t=0.5$.

Then from the leading eigenvalue $N_{c}$ of the pairing matrix, we can obtain
the total condensate fraction as $n_{c}=N_{c}/(N/2)$ [65, 66]. The
corresponding eigenstate of $N_{c}$ is the momentum-space pair wave function
$\boldsymbol{\Psi}_{c}$, which consists of the singlet and triplet components
as
$\boldsymbol{\Psi}_{c}=(\boldsymbol{\Psi}_{c,s},\boldsymbol{\Psi}_{c,t_{\uparrow}},\boldsymbol{\Psi}_{c,t_{\downarrow}})^{\rm
T}$ with every component as a $N_{s}$-dimensional vector. For the lattice
model in Eq. (1), the two triplet channels are degenerate as
$\boldsymbol{\Psi}_{c,t_{\uparrow}}=\boldsymbol{\Psi}_{c,t_{\downarrow}}$ due
to the spin-inversion symmetry. Thus, we define the overall triplet pair wave
function as
$\boldsymbol{\Psi}_{c,t}=\sqrt{2}\boldsymbol{\Psi}_{c,t_{\uparrow}}$. Then
within normalized $\boldsymbol{\Psi}_{c}$, we assign the condensate fractions
of spin-singlet and triplet pairing as
$n_{c,s}=n_{c}\times(\boldsymbol{\Psi}_{c,s}\boldsymbol{\Psi}_{c,s}^{\rm T})$
and $n_{c,s}=n_{c}\times(\boldsymbol{\Psi}_{c,t}\boldsymbol{\Psi}_{c,t}^{\rm
T})$, respectively. Thus, the relation $n_{c}=n_{c,s}+n_{c,t}$ obviously
holds, with $n_{c,s}/n_{c}$ and $n_{c,t}/n_{c}$ as the contributions of
singlet and triplet channels to the pairing. The square
$|\boldsymbol{\Psi}_{c,\ell}({\bf{k}})|^{2}$ ($\ell=s,t$) stands for the
probability of fermions with momentum $\bf{k}$ participating the pairing. We
can further obtain the corresponding real-space pair wave functions
$\psi_{c,s}(\bf{r})$ and $\psi_{c,t}(\bf{r})$ by Fourier transform of
$\boldsymbol{\Psi}_{c,s}$ and $\boldsymbol{\Psi}_{c,t}$. Similarly,
$|\psi_{s}({\bf{r}})|^{2}$ and $|\psi_{t}({\bf{r}})|^{2}$ represent
probabilities of spin-singlet and triplet Cooper pairs with distance $\bf{r}$
of the two fermions, and they actually reflect the size of the pairs.

We note that the finite-temperature condensate fraction should vanish in the
thermodynamic limit for both the quasi-condensate and disordered phases.
Alternatively, its finite-size results can present plentiful information about
the pairing properties [53, 55, 67]. Attributed to the distinct behaviors of
pairing correlation accross the BKT transition, the finite-size condensate
fraction should also undergo a qualitative change around the transition, which
serves as an effective tool to locate the BKT transition temperature [67].
Thus, throughout this work, we always deal with finite-size AFQMC results of
condensate fraction at finite temperatures.

As a demonstration, we show the typical results of condensate fractions and
pair wave functions for a specific group of parameters in Fig. 1. With
lowering temperature, both the spin-singlet and triplet condensate fractions
monotonically increases from the high-temperature normal state to the low-
temperature superfluid phase, and then saturates to the ground-state values as
indicated by the plateau achieved with $T/t\leq 0.06$ results. They also
exhibit a rapid increase at a specific temperature, for which the BKT
transition should be responsible (see Sec. III.4). Remarkably, the triplet
channel has a rather small contribution to the total condensate fraction (less
than $15\%$ approaching $T=0$). As shown in Fig. 1(b), the momentum-space pair
wave functions of both singlet and triplet channels have peaks around the
Fermi surfaces of the two helical bands for the intermediate interaction
$U/t=-4$ at low temperature, which resembles the results without SOC [59, 67].
This is also consistent with the fundamental picture of BCS theory that the
fermions around the Fermi surfaces dominates the pairing in the weakly
interacting regime [2]. In contrast, the high-temperature results of the pair
wave functions seem featureless as the system is in the normal state. Note
that the node at $\boldsymbol{\Gamma}$ point in triplet pair wave function
indicates its antisymmetry, while the singlet component is symmetric without
node.

###  III.2 The mixed-parity pairing structure

Based on the results and discussions in Sec. III.1, we then concentrate on the
mixed-parity pairing structure for physically relevant parameter regimes,
revealed by the numerical results of condensate fractions and pair wave
functions. tunning parameters including the temperature, interaction strength,
SOC and the chemical potential are accounted for in our AFQMC simulations.

The condensate fractions versus the temperature typically shares similar
behavior as the results shown in Fig. 1(a), with differences lying in the
specific numbers and BKT transition temperatures. With varying interaction
strength, SOC and chemical potential, we perform the AFQMC simulations with
$T/t=0.10$, at which most systems studied as follows fall into the superfluid
phase and the condensate fractions are close to the corresponding $T=0$
results [similar to Fig. 1(a)].

Figure 2:  The condensate fractions $n_{c}$, $n_{c,s}$ and $n_{c,t}$ versus
various tunning parameters. (a) Tune interaction strength $U/t$ with
$\lambda/t=0.5$ and $\mu/t=0.5$; (b) Tune SOC $\lambda/t$ with $U/t=-4$ and
$\mu/t=0.5$; (c) Tune the fermion filling $n$ by chemical potential $\mu/t$
with $U/t=-4$ and $\lambda/t=0.5$. The insets in Panel (a) and (b) plot the
results of corresponding filling. These calculations are performed for $L=20$
system with temperature $T/t=0.1$. Figure 3:  The singlet (top) and triplet
(bottom) pair wave functions in momentum space,
$|\boldsymbol{\Psi}_{c,s}({\bf{k}})|$ and
$|\boldsymbol{\Psi}_{c,t}({\bf{k}})|$, versus chemical potential $\mu/t$ with
corresponding fermion filling $n$ (error bars are in the fourth/fifth digits
and are thus neglected) shown on top of the plots. The error bars of $n$ are
in the fourth or fifth digits and thus are neglected. The noninteracting Fermi
surfaces at $T=0$ of the two helical bands are also plotted with the green and
blue solid lines. These calculations are performed for $L=20$ system with
$T/t=0.10$ and $U/t=-4$, $\lambda/t=0.5$. Figure 4:  The singlet (top) and
triplet (bottom) pair wave functions in real space, $|\psi_{s}({\bf{r}})|$ and
$|\psi_{t}({\bf{r}})|$, versus chemical potential $\mu/t$. Simulation
parameters are the same as Fig. 3.

In Fig. 2, we present one of the key results in this work, as condensate
fractions with tunning parameters other than the temperature. First, with
increasing on-site interaction, the system hosts the BCS-BEC crossover from
extended Cooper pairs to tightly bounded molecules [68]. Fig. 2(a) shows that
the singlet condensate fraction simply increases during the crossover, while
the triplet contribution has a peak around $U/t=-4$ (for $\lambda/t=0.5$).
This difference can be understood as follows. Turning on the interaction can
first enhance the pairing in both channels as well as the condensate
fractions. Cross some intermediate $U$, the interaction begins to frustrate
the triplet pair formation and continues to increase the singlet pairs, due to
the nature of the attraction between fermions with unlike spins. The results
suggest that the triplet contribution to the pairing is most significant
[$\sim 21\%$ as in Fig. 2(a)] in the intermediate interaction regime, whose
specific value of $U/t$ should depend on SOC strength. These results are
qualitatively consistent with those from the ground-state calculations of 2D
spin-orbit-coupled Fermi gas [53]. Then, with tunning SOC strength, the
condensate fractions for $U/t=-4$ are plotted in Fig. 2(b). The decreasing of
singlet condensate fraction with $\lambda/t$ can be explained by the enlarged
bandwidth $W(t,\lambda)$ [54] and the reduced effective interaction $U/W$.
However, the triplet condensate fraction is first enhanced by SOC, due to the
fact that SOC is the essential source of triplet pairing in presence of
Hubbard interaction. Then the effect from reduced $U/W$ sets in, and the
competition results in a broaden peak around $\lambda/t=0.5\sim 1.0$. The
biggest contribution from the triplet channel to the pairing is $\sim 30\%$
around $\lambda/t=1.3$, where nevertheless the total condensate fraction is
only $0.043$. Finally, in Fig. 2(c), we show the numerical results versus the
fermion filling (by tunning the chemical potential). Both the singlet and
total condensate fractions reach the maximum around $n=0.80$, while the
triplet one possesses a wide plateau regarding the filling. The triplet
contribution saturates to largest value $\sim 20\%$ towards the low filling
regime. As discussed above, for the simulation temperature $T/t=0.10$, the
system evolves from the normal state at half filling (with $\mu=0$) to the
superfluid phase with increasing doping. Thus, the results in Fig. 2(c) might
indicate that the maximal BKT transition temperature is achieved around the
filling $n=0.80$ [54]. Combining all the results in Fig. 2, we can conclude
that the spin-singlet pairing always has the predominant contribution than the
triplet channel to the mixed-parity pairing in the system.

Figure 5:  The triplet pair wave function in real space $|\psi_{t}({\bf{r}})|$
versus (top row) interaction strength $-U/t=2,4,6,8$ with
$T/t=0.10,\lambda/t=0.5,\mu/t=0.5$; (middle row) the temperature
$T/t=0.025,0.05,0.14,0.20$ with $U/t=-4,\lambda/t=0.5,\mu/t=0.5$; (bottom row)
SOC strength $\lambda/t=0.3,0.6,1.2,1.5$ with $T/t=0.10,U/t=-4,\mu/t=0.5$.
These calculations are performed for $L=20$ system.

We then turn to the results of the pair wave functions. First, their
evolutions versus the chemical potential in momentum space and real space are
illustrated in Fig. 3 and Fig. 4, respectively. For half filling, our results
are quantitatively consistent with the $T=0$ results in Ref. [55]. Increasing
the chemical potential results in smaller fermion filling, and the
corresponding noninteracting Fermi surfaces at $T=0$ of the two helical bands
(dashed lines in Fig. 3), which are determined from the corresponding fermion
filling from finite-T AFQMC calculation, shrinks towards circles. It’s clear
that for the intermediate interaction $U/t=-4$ the pair wave functions in both
channels show sharp peaks in the vicinity of the Fermi surfaces, regardless of
the filling. With the increasing interaction, the results should gradually
become smooth in the whole Brillouin zone (not shown) without apparent peaks
[53] indicating the deviation from BCS theory. In contrast, the singlet and
triplet pair wave functions in real space show significant difference, as
shown in Fig. 4. The localized peaks in singlet pair wave function
$|\psi_{s}({\bf{r}})|$ clearly shows that the singlet pairing mainly has a
local origin with on-site pairs. However, the Pauli principle prohibits such
on-site triplet pair formation, resulting in zero value at $\bf{r}=0$.
Instead, the triplet pair wave function $|\psi_{t}({\bf{r}})|$ is more
extended and has very rich patterns and evolutions along with decreasing
fermion filling. Multi-peak structures appear in $|\psi_{t}({\bf{r}})|$, with
the largest amplitude locations changing from the next-nearest-neighbor (NNN)
sites at half filling, to intermediate fourth-nearest-neighbor ($4^{\rm
th}$-NN), and finally to the nearest-neighbor (NN) sites at low filling, as
shown in Fig. 6(a). These finite-size triplet Cooper pairs within several NN
sites can be explained by the real-space nature of Rashba SOC, which is
actually NN spin-flip hopping. As a result, successive SOC hops can enhance
the possibility to find another fermion at neighboring sites with the same
spin as the one located at origin. Moreover, towards smaller filling, both of
$|\psi_{s}({\bf{r}})|$ and $|\psi_{t}({\bf{r}})|$ show more extended behaviors
due to the enlarged wavelength $\sim 2\pi/k_{F}$.

Figure 6:  The amplitudes of real-space triplet pair wave function
$|\psi_{c,t}(\mathbf{r})|$ with the distance $r$ equal to NN, NNN, $3^{\rm
rd}$-NN and $4^{\rm th}$-NN sites, with (a) tunning the fermion filling $n$
and (b) varying the SOC strength. The other simulation parameters for Panel
(a) and (b) are the same as Fig. 4 and the bottom row of Fig. 5, respectively.
These calculations are performed for $L=20$ system.

As for the other tunning parameters, the pair wave functions in momentum space
show similar behaviors as illustrated in Fig. 3, and in real space for singlet
channel as $|\psi_{s}({\bf{r}})|$ are also dominated by the center peak as
Fig. 4. Thus, we now concentrate on triplet pair wave function in real space
$|\psi_{t}({\bf{r}})|$ shown in Fig. 5. With increasing interaction strength
[top row of Fig. 5], $|\psi_{t}({\bf{r}})|$ gradually evolves from rather
extended pattern with multipeaks along diagonals, to local peaks at NN lattice
sites, which illustrates the BCS-BEC crossover behavior in triplet pairing
channel. With decreasing temperature [middle row of Fig. 5], the peaks in
$|\psi_{t}({\bf{r}})|$ (located at NNN and $4^{\rm th}$-NN sites) simply
become more significant and eventually stabilize, indicating entering the
superfluid phase from the normal state. Increasing SOC strength, the $T=0$
AFQMC simulations at half filling [55] showed a diamond pattern of
$|\psi_{t}({\bf{r}})|$ with enhanced peak values at both NNN and $3^{\rm
rd}$-NN sites. It behaves differently away from half filling [bottom row of
Fig. 5]. As shown in Fig. 6(b), SOC first enhances all the finite-range
triplet pairing for $\lambda/t<0.75$, where NNN and $4^{\rm th}$-NN pairing
play the leading role. The NNN component is then further promoted by SOC, and
NN and NNN pairing gradually becomes comparable towards large SOC, resulting
in instead a square pattern as illustrated in Fig. 5. All the qualitative
behaviors of pair wave functions results in Fig. 3, Fig. 4 and Fig. 5 do not
change with the system size.

###  III.3 The pairing correlation functions

In Sec. III.2, the results in Fig. 2 clearly present an optimal SOC strength
and fermion filling regime where the triplet condensate fraction reaches the
maximum. In this section, we pursue to understand this point from the aspect
of the pairing correlation functions.

Figure 7:  Vertex contribution of on-site singlet (top), NN triplet (middle)
and NNN triplet (bottom) pairing correlation functions
$\bar{P}_{\ell}(\mathbf{r})$ (with $\ell=s,t$) versus SOC strength. The
correlations with distance $r=5\sim 10$ (along the $x$ axis) are plotted.
These calculations are performed for $L=20$ system with $T/t=0.10$ and
$U/t=-4$, $\mu/t=0.5$.

We define the real-space singlet and triplet pairing operators as

$\displaystyle\hat{\Delta}_{s,\bf{i}}=(c_{\bf{i}\uparrow}^{\dagger}c_{\bf{i}\downarrow}^{\dagger}+c_{\bf{i}\downarrow}c_{\bf{i}\uparrow})/2,$
(4)
$\displaystyle\hat{\Delta}_{t,\bf{i}}=(c_{\bf{i}\uparrow}^{\dagger}c_{\bf{i}+\boldsymbol{\delta}\uparrow}^{\dagger}+c_{\bf{i}+\boldsymbol{\delta}\uparrow}c_{\bf{i}\uparrow})/2.$

with $s$ and $t$ as singlet and triplet. For the triplet, we concentrate on NN
and NNN pairing with $\boldsymbol{\delta}=(1,0)$ and
$\boldsymbol{\delta}=(1,1)$ denoting the corresponding lattice vectors. We
then measure the real space correlation functions
$P_{s}(\mathbf{r})=\langle\hat{\Delta}_{s,\mathbf{i}}\hat{\Delta}_{s,\mathbf{i}+\mathbf{r}}\rangle$
and
$P_{t}(\mathbf{r})=\langle\hat{\Delta}_{t,\mathbf{i}}\hat{\Delta}_{t,\mathbf{i}+\mathbf{r}}\rangle$,
and the structure factors as their Fourier transformation
$S_{\ell}(\mathbf{q})=\sum_{\mathbf{r}}P_{\ell}(\mathbf{r})e^{i\mathbf{q}\cdot\mathbf{r}}$
with $\ell=s$ or $t$. To directly evaluate the pure interaction contribution,
we have also obtained the vertex contribution for the pairing correlations and
structure factors, as $\bar{P}_{s}(\mathbf{r})$ and
$\bar{S}_{\ell}(\mathbf{q})$, by subtracting the uncorrelated part [69].

Figure 8:  Vertex contribution of on-site singlet (top), NN triplet (middle)
and NNN triplet (bottom) pairing structure factors
$\bar{S}_{\ell}(\mathbf{q}=\boldsymbol{\Gamma})$ (with $\ell=s,t$) versus SOC
strength. These calculations are performed for $L=16,18,20$ systems with
$T/t=0.10$ and $U/t=-4$, $\mu/t=0.5$.

In Fig. 7, we present the vertex contributions to the pairing correlation
functions of on-site singlet, NN and NNN triplet channels, with tunning SOC
strength. All the positive vertex contributions to the correlations in Fig. 7
reveal that the on-site attractive interaction enhances the singlet and
triplet pairing correlations with specific distances as $L/4\leq r\leq L/2$
(as $L=20$). These results contribute more than $90\%$ of the corresponding
bare correlation functions (not shown). It’s clear that SOC simply suppresses
the singlet pairing correlation, while the NN and NNN triplet correlations
show broaden peaks around $\lambda/t=0.8$ and $\lambda/t=1.0$. Moreover, the
singlet correlation is stronger than the triplet ones around two orders of
magnitude, indicating the dominant role of singlet channel. These results are
in accordance with the behaviors of the corresponding condensate fractions
shown in Fig. 2(b). The almost collapsed numerical data for different
distances in Fig. 7 also highlight the superfluid phase of the system for the
chosen parameters. Then, the vertex contributions of the pairing structure
factors $\bar{S}_{\ell}(\mathbf{q}=\boldsymbol{\Gamma})$ (with $\ell=s,t$)
with increasing SOC strength are illustrated in Fig. 8. They show similar
behaviors as the real-space correlation functions. The negative vertex of
$\bar{S}_{t}(\boldsymbol{\Gamma})$ with NN triplet for $\lambda/t<0.3$ and
$\lambda/t>1.2$ means that the NN triplet pairing is not favored in these
regimes. The growing numbers of $\bar{S}_{\ell}(\boldsymbol{\Gamma})$ in Fig.
8 for all three quantities (especially in the intermediate SOC regime) with
increasing system size also suggest quasi-long-range pairing orders.

Figure 9:  Vertex contribution of (a) on-site singlet, (b) NN triplet and NNN
triplet pairing structure factors $\bar{S}_{\ell}(\boldsymbol{\Gamma})$ versus
fermion filling. These calculations are performed for $L=20$ systems with
$T/t=0.10$ and $U/t=-4$, $\lambda/t=0.5$.

Similarly, the results of condensate fractions versus the fermion filling in
Fig. 2(c) can also be alternatively understood from the pairing correlations.
Fig. 9 plots the vertex contributions of the pairing structure factors
$\bar{S}_{\ell}(\boldsymbol{\Gamma})$ versus fermion filling. The results of
on-site singlet structure factor has the same nonmonotonic behavior as the
singlet condensate fraction in Fig. 2(c). Instead, the NN and NNN triplet
structure factors show more interesting signatures with different peak
locations, revealing that the triplet channel is first governed by the NNN and
then by the NN pairing from half filling to low filling regime ($n<0.5$).
These validate the results of condensate fractions in Fig. 2(c) and pair wave
functions in Fig. 3. Moreover, the wide plateau of the triplet condensate
fraction in Fig. 2(c) can be explained by the accumulated results of the NN
and NNN triplet pairing correlations in Fig. 9(b).

As for the temperature and interaction strength, we have also obtained the
vertex contributions of both singlet and triplet pairing correlations. In
Appendix B, we have presented the results of vertex contributions
$\bar{P}_{\ell}(\mathbf{r})$ versus temperature.

###  III.4 BKT transition temperature from condensate fractions

In previous studies of 2D attractive Hubbard model, the BKT transition
temperature was usually determined by the finite-size scaling of the pairing
structure factor or from the universal jump property of the superfluid density
[70, 71, 54]. However, these quantities can become significantly small for low
filling system, which makes the finite-size scaling even harder. In Ref. [54],
the BKT transition temperatures for the same system as we study were
calculated from superfluid density with systems up to $L=12$. Such AFQMC
simulations need to compute the superfluid density from dynamical current-
current correlation functions, which definitely cost much more computational
effort to reach high-precision results.

Alternatively, numerical studies in 2D XY models solidly confirm that the
finite-size BKT transition temperature $T_{\rm BKT}(L)$ has a following form
[72, 73] as

$\displaystyle T_{\rm BKT}(L)=T_{\rm BKT}(L=\infty)+\frac{a}{(\ln L+b)^{2}},$
(5)

with $a,b$ as coefficients related to the specific problem, and $T_{\rm
BKT}(L=\infty)$ as the final answer under thermodynamic limit. The second term
in Eq. (5) containing the logarithm of linear system size $L$ already
indicates the strong finite-size effect. As a result, biased result of $T_{\rm
BKT}(L=\infty)$ might be obtained if only a small group of systems with not
large enough sizes are accessed in the calculations. Based on Eq. (5), we
could extrapolate the precise $T_{\rm BKT}(L=\infty)$ from the finite-size
$T_{\rm BKT}(L)$ results. In the previous study of the 2D interacting Fermi
gas without SOC [67], it was found that the first-order derivative of
condensate fraction over the temperature shows a peak and its location can be
identified as $T_{\rm BKT}(L)$. Such calculations do not involve dynamical
measurements and are thus computationally much cheaper and high-precision
results of $T_{\rm BKT}(L)$ can be yielded. Similar formula as Eq. (5) was
also confirmed dealing with the convergence of $T_{\rm BKT}$ with number of
fermions for 2D Fermi gas in Ref. [67]. Thus, in the following, we also
concentrate on calculating $T_{\rm BKT}(L)$ from the condensate fractions and
reaching $T_{\rm BKT}(L=\infty)$ using Eq. (5), for the system with SOC.

Figure 10:  Determination of the BKT transition temperatures from correlation
ratio and total condensate fraction. (a)(b) are the singlet pairing structure
factor $S_{s}(\boldsymbol{\Gamma})$ and the corresponding correlation ratio.
The inset in panel (b) plots the cross points of finite-size correlation
ratios. (c)(d) are the total condensate fraction and its first order
derivative (after polynomial fitting). The inset in panel (d) plots $T_{\rm
BKT}(L)$ after the best fitting using Eq. (5), reaching the final result as
$T_{\rm BKT}(L=\infty)/t=0.135(4)$. These calculations are performed for
$L=8\sim 20$ systems with $U/t=-4,\lambda/t=0.5,\mu/t=0.5$.

In Fig. 10, we first demonstrate the determination of BKT transition
temperature from the spin-singlet pairing structure factor
$S_{s}(\boldsymbol{\Gamma})$ (defined in Sec. III.3) and the total condensate
fraction as a comparison. The correlation ratio for
$S_{s}(\boldsymbol{\Gamma})$ is defined as $R_{\rm
corr}=1-S_{s}(\boldsymbol{\Gamma}+\mathbf{q})/S_{s}(\boldsymbol{\Gamma})$ with
$\mathbf{q}$ as the smallest momentum on the lattice, i.e., $(2\pi/L,0)$ or
$(0,2\pi/L)$. It resembles the Binder cumulant which converges to unity in
ordered phase and vanishes in the disordered phase in thermodynamic limit.
Then the cross points of the finite-size $R_{\rm corr}$ results can be
approximately viewed as the transition temperature. As shown in Fig. 10(b),
the cross points of $R_{\rm corr}$ indeed move to the lower temperature with
system size but not with a well defined behavior. Instead, for the total
condensate fraction in Fig. 10(c), we first perform a polynomial fitting to
the numerical data and then compute its first-order derivative and get the
location of the peak as $T_{\rm BKT}(L)$ [shown in Fig. 10(d)], which avoids
the step error involved in the numerical derivative. We have further
calculated the error bars of $T_{\rm BKT}(L)$ applying the standard
bootstrapping technique. Finally, we use Eq. (5) to extrapolate the final
result of BKT transition temperature $T_{\rm BKT}(L=\infty)=0.135(4)$, as
plotted in inset of Fig. 10(d). The details of the bootstrapping calculations
of $T_{\rm BKT}(L)$ are presented in Appendix C. These results also indicates
large finite-size effect in $R_{\rm corr}$ as the cross point of $L=18$ and
$L=20$ is $T/t=0.158$, which strongly deviates from $T_{\rm BKT}(L=\infty)$.

Figure 11:  Determination of the BKT transition temperatures from singlet and
triplet condensate fractions. (a)(b) are the singlet condensate fraction and
its first-order derivative (after polynomial fitting). (c)(d) are the triplet
condensate fraction and its first-order derivative. The inset in panel (b)(d)
plots $T_{\rm BKT}(L)$ after the best fitting using Eq. (5), reaching the
final result as $T^{s}_{\rm BKT}(L=\infty)/t=0.135(4)$ and $T^{t}_{\rm
BKT}(L=\infty)/t=0.135(4)$, respectively. Simulation parameters are the same
as Fig. 10.

For the mixed-parity pairing we stuty, we also have the numerical data of
singlet and triplet condensate fractions. From them, we can separately
extrapolate the BKT transition temperatures for the spin-singlet and triplet
superfluidity as $T_{\rm BKT}^{s}(L=\infty)$ and $T_{\rm BKT}^{t}(L=\infty)$,
which are expected be the same. In Fig. 11, we illustrate the determination of
$T_{\rm BKT}$ from both singlet and triplet condensate fractions. The
procedure is exactly the same as that in Fig. 10(c) and (d), and the details
can also be refered in Appendix C. These calculations produce the final
results of $T_{\rm BKT}^{s}(L=\infty)=0.135(4)$ and $T_{\rm
BKT}^{t}(L=\infty)=0.135(4)$. These results are indeed consistent as
expectation, meaning the BKT transition for the quasi-long-range mixed-parity
pairing order happens simultaneously for singlet and triplet channels.

Based on the results in Fig. 10 and Fig. 11, we have obtained the BKT
transition temperature $T_{\rm BKT}(L=\infty)=0.135(4)$ for the parameter
$U/t=-4,\lambda/t=0.5,\mu/t=0.5$ (with fermion filling $n=0.6795$ at the
transition point). This result is consistent with the $T_{\rm BKT}$ computed
for filling $n=0.7$ in Ref. [54]. We can then conclude that, similar to
previous studies [67], it’s also an efficient way to determine BKT transition
temperature from (total, singlet and triplet) condensate fractions for
attractive fermion systems with SOC.

##  IV Summary and Discussion

The mixed-parity pairing phenomena is theoretically a natural result for
fermionic systems with broken inversion symmetry [24, 10], and it has been
experimentally observed in various three-dimensional superconductors with SOC
[13]. In addition, the experimental realization of SOC with an artificial
gauge field in optical lattice by ultracold atoms [43, 44, 45] provides the
opportunity to perform more systematic and deeper study of the mixed-parity
pairing in a more controlled manner. Our AFQMC numerical results in this work
can not only serve as quantitative guide for such 2D optical lattice
experiments, but also present some new physical results on the essential
pairing structure of the corresponding mixed-parity pairing.

In summary, we have applied the numerically exact finite-temperature AFQMC
method to study the pairing properties of attractive fermions with Rashba SOC
in 2D optical lattice. We evaluate the contributions of the spin-singlet and
triplet channels to the mixed-parity pairing. With the scanning of
temperature, fermion filling, SOC and interaction strengths, we find that the
singlet pairing plays a dominant role with relatively small triplet
contribution in most relevant parameter regimes. From the pair wave functions,
we find that, for intermediate interaction ($U/t=-4$), the singlet pairing
mainly consists of local Cooper pairs while the triplet channel is rather
extended with major contributions from several nearest neighbors. Especially,
in low filling regime ($n<0.5$), the triplet pairing is dominated by NN
fermion pairs, in contrast with the NNN ones around half filling. Via the
vertex contribution of pairing correlations, we have shown that the triplet
pairing is first enhanced and then suppressed with increasing SOC, and there
exists an optimal SOC strength for observing the triplet pairing. Finally, we
have demonstrated the computation of the BKT transition temperature from the
finite-size results of total, singlet and triplet condensate fractions,
suggesting it also as an efficient method for systems with SOC. Our numerical
results will surely offer useful benchmarks for future optical lattice
experiments as well as theories and other numerical methods.

Our work also has implications for achieving the spin-triplet
superconductivity and superfluidity. Considering the fact that the triplet
pairing is only confirmed to exist in very rare systems, it might be a way out
to pay more attention to the systems with mixed-parity pairing. Specifically,
if one can control the triplet contribution to the pairing by tunning physical
parameters (for example, the SOC strength) in such systems, we might access
the special case in which the triplet channel dominates, similar to
$\text{Li}_{2}\text{Pt}_{3}\text{B}$ [16]. Unfortunately, our work shows it’s
very unlikely to realize such special case for the system described by the
lattice model in Eq. (1). Instead, there are actually other possibilities,
such as further including the Dresselhaus SOC, and NN or NNN attractive
interactions. The former was found to be useful in promotion of the triplet
contribution in interacting Fermi gas within the mean-field theory [47,
*Anna2012]. The latter is apparently supported by our numerical results as the
triplet pairing is mainly contributed by NN and NNN Cooper pairs. We leave
these open possibilities to future work.

###### Acknowledgements.

Y.Y.H. acknowledges Peter Rosenberg and Shiwei Zhang for valuable discussions.
This work was supported by the National Natural Science Foundation of China
(under Grant No. 12047502, 12204377 and 12275263), the Innovation Program for
Quantum Science and Technology (under Grant No. 2021ZD0301900), and the Youth
Innovation Team of Shaanxi Universities.

##  Appendix A Structure factor of the density-density correlation function

In Ref. [55], it was found that, at half filling, the long-range charge
density wave (CDW) order with checkerboard pattern coexist with the pairing
order in ground state for the lattice model in Eq. (1). We have also checked
this and our numerical results suggest that the long-range CDW order should
not exist for the case away from half filling.

We compute the density-density correlation function defined as
$D(\mathbf{r})=\frac{1}{4}(\langle\hat{n}_{\mathbf{i}}\hat{n}_{\mathbf{i}+\mathbf{r}}\rangle-\langle\hat{n}_{\mathbf{i}}\rangle\langle\hat{n}_{\mathbf{i}+\mathbf{r}}\rangle)$
(with
$\hat{n}_{\mathbf{i}}=\hat{n}_{\mathbf{i}\uparrow}+\hat{n}_{\mathbf{i}\downarrow}$),
and the corresponding momentum-space structure factor as $S_{\rm
CDW}(\mathbf{q})=\sum_{\mathbf{r}}D(\mathbf{r})e^{i\mathbf{q}\cdot\mathbf{r}}$.
The leading component of $S_{\rm CDW}(\mathbf{q})$ appears at
$\mathbf{q}=\mathbf{M}=(\pi,\pi)$ point, consistent with the CDW order with
the checkerboard pattern.

Figure 12:  The density-density structure factor at $\Gamma$ point for
different quantities. The detailed parameters are shown in the figure.
Compared with half-filling case ($\mu=0$), the structure factor has an obvious
decrease, which indicates that the superconducting order and CDW do not
coexist in the doped system.

In Fig. 12, we illustrate the results of the CDW structure factor $S_{\rm
CDW}(\mathbf{M})$ with various tunning parameters. First, with doping as
increasing the chemical potential, $S_{\rm CDW}(\mathbf{M})$ immediately
decrease from the half-filling result by approximately an order of magnitude
for $n=0.94$, which suggests the significant suppression of CDW order away
from half filling. Second, the results with lowering temperature with
$\mu/t=0.5$ (around $n=0.68$) explicitly shows that $S_{\rm CDW}(\mathbf{M})$
first decreases, then reaches a minimum and gradually saturates towards $T=0$,
indicating the absence of long-range order. The results with varying SOC and
interaction strengths show some enhancements of $S_{\rm CDW}(\mathbf{M})$ for
specific regimes, but its largest values are still much smaller than the half-
filling results, which also suggests the short-range correlations.

##  Appendix B Vertex contributions of the pairing correlation functions
versus temperature

In Sec. III.3, we have shown the numerical results of vertex contributions for
the pairing correlation functions with tunning SOC strength and chemical
potential. Here, we present more results with varying temperature.

Figure 13:  Vertex contribution of on-site singlet (top), NN triplet (middle)
and NNN triplet (bottom) pairing correlation functions $P_{\ell}(\mathbf{r})$
(with $\ell=s,t$) versus temperature, with $r=\sqrt{2}L/2=10\sqrt{2}$ [as
$\mathbf{r}=(L/2,L/2)$] as the largest distance. Results with several SOC
strengths $\lambda/t=0\sim 1.5$ are plotted. These calculations are performed
for $L=20$ system with $T/t=0.10$ and $U/t=-4$, $\mu/t=0.5$.

In Fig. 13, we present the vertex of reals-space pairing correlations
$P_{\ell}(\mathbf{r})$ (with $\ell=s,t$) with the largest distance as
$r=\sqrt{2}L/2$ [as $\mathbf{r}=(L/2,L/2)$] on the lattice with on-site
singlet, NN and NNN triplet channels versus the temperature, for several SOC
strengths as $\lambda/t=0\sim 1.5$. All the results show enhancements with
decreasing temperature and plateaus appear saturating to the $T=0$ results,
indicating the quasi-long-range pairing order at low temperature regime. At
low temperature regime, the triplet correlations reach the maximum at
$\lambda/t=0.9$ for both NN and NNN pairing, consistent with the results shown
in Fig. 7. Besides, these results also illustrate that the results at
$T/t=0.10$ is very close to the $T=0$ correspondences.

##  Appendix C The determination of the BKT transition temperature $T_{\rm
BKT}(L)$ and $T_{\rm BKT}(L=\infty)$

In this section, we present the details for the determination of the BKT
transition temperature $T_{\rm BKT}(L)$ and $T_{\rm BKT}(L=\infty)$.

Figure 14: Condensate fraction versus temperature. The points are QMC data for
system size $L=18$ and $20$, where the error bars denote the standard error.
The curves are the fitting results of different random
$n_{c,i}(L,T,q)=\bar{n}_{c}(L,T)+N(0,q\sigma(L,T))$. For simplicity, we plot
five random curves for each system. The parameters are $U/t=-4$,
$\lambda/t=0.5$ and $\mu/t=0.5$.

Based on the numerical data of condensate fraction $\bar{n}_{c}(L,T)$
(including the total, singlet and triplet) and the corresponding standard
error $\sigma(L,T)$, we apply the bootstrapping technique by first generating
a set of random data by

$\displaystyle n_{c,i}(L,T,q)=\bar{n}_{c}(L,T)+N(0,q\sigma(L,T)),$ (6)

where $i$ denotes the $i$-th random data with $q=1,2,3$ for different range of
deviation, and $N(0,q\sigma(L,T))$ stands for the normalized Gaussian
distribution with expectation and standard deviation as $0$ and
$q\sigma(L,T)$. The whole process follows the Gauss analysis and can quickly
construct a large number of $n_{c,i}(L,T,q)$. Then we fit $n_{c,i}(L,T,q)$ for
every set of random data with a fourth-order polynomial of temperature around
the transition point, and compute the peak location of its first-order
derivative, and then we take it as $T_{\rm BKT}(i,L,q)$. With the full set of
$T_{\rm BKT}(i,L,q)$, one can perform data analysis and obtain
$\overline{T}_{\rm BKT}(L,q)$ with the stand deviation as its error bar.
Compared with the method to only fit $n_{c}(L,T)$ with the original data, this
bootstrapping method can additional present a reasonable error bar for
$\overline{T}_{\rm BKT}(L,q)$.

Figure 15: The distribution of $T_{\rm BKT}(L,q)$ based on the bootstrapping
calculations. We have generated 150000 random data for each system size and
$q$. It is well illustrated that the average values $\overline{T}_{\rm
BKT}(L)$ for different $q$ are identical. Simulation parameters are the same
as Fig. 14. Figure 16: BKT transition temperature $T_{\rm BKT}(L,q)$ versus
$(\ln L-0.7)^{-2}$ for the total, singlet and triplet condensate fractions.
The fittings are based on the correction formula Eq. (5), where $-0.7$ is
determined from the fitting. The singlet and triplet channels give very
similar results. Simulation parameters are the same as Fig. 14.

In order to show the process, we take $L=18$ and $20$ for example. Fig. 14
shows the original data and fitting polynomials of five random sets of data.
It is shown that the fourth-order polynomials can capture the essential
behavior of $n_{c}$ around the transition point. By generating 150000 samples,
Fig. 15 shows the histogram of results for $\overline{T}_{\rm BKT}(L,q)$,
which are fairly consistent with Gaussian distributions. It is also well
illustrated that, for different $q$, the average $\overline{T}_{\rm BKT}(L,q)$
are obviously identical for both system sizes. The only difference of these
results the data set generated by different Gaussian noise (different $q$)
shows in the standard deviations of these distributions. As espectation, the
distribution is wider (with larger standard deviation) for larger $q$.

Finally, to obtain $T_{\rm BKT}$ in the thermodynamic limit, we perform
fittings of $T_{\rm BKT}(L,q)$ with the formula in Eq. (5). Fig. 16 shows the
fitting results with different $q$. The total condensate fraction as well as
its two channels give similar $T_{\rm BKT}$, as shown in Fig. 16. To achieve a
confident estimate of $T_{\rm BKT}(L=\infty)$, we adopt the results of $q=3$
as the final result as presented in Sec. III.4 of the main text.

## References

  * Ginzburg [2004] V. L. Ginzburg, Nobel lecture: On superconductivity and superfluidity (what i have and have not managed to do) as well as on the “physical minimum” at the beginning of the xxi century, Rev. Mod. Phys. 76, 981 (2004).
  * Cooper [1956] L. N. Cooper, Bound electron pairs in a degenerate fermi gas, Phys. Rev. 104, 1189 (1956).
  * Anderson [1984] P. W. Anderson, Structure of ”triplet” superconducting energy gaps, Phys. Rev. B 30, 4000 (1984).
  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theory of superconductivity, Phys. Rev. 108, 1175 (1957).
  * Tsuei and Kirtley [2000] C. C. Tsuei and J. R. Kirtley, Pairing symmetry in cuprate superconductors, Rev. Mod. Phys. 72, 969 (2000).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a mott insulator: Physics of high-temperature superconductivity, Rev. Mod. Phys. 78, 17 (2006).
  * Armitage  _et al._ [2010] N. P. Armitage, P. Fournier, and R. L. Greene, Progress and perspectives on electron-doped cuprates, Rev. Mod. Phys. 82, 2421 (2010).
  * Lee [1997] D. M. Lee, The extraordinary phases of liquid 3he, Rev. Mod. Phys. 69, 645 (1997).
  * Mackenzie and Maeno [2003] A. P. Mackenzie and Y. Maeno, The superconductivity of ${\mathrm{sr}}_{2}{\mathrm{ruo}}_{4}$ and the physics of spin-triplet pairing, Rev. Mod. Phys. 75, 657 (2003).
  * Gor’kov and Rashba [2001] L. P. Gor’kov and E. I. Rashba, Superconducting 2d system with lifted spin degeneracy: Mixed singlet-triplet state, Phys. Rev. Lett. 87, 037004 (2001).
  * Frigeri  _et al._ [2004a] P. A. Frigeri, D. F. Agterberg, A. Koga, and M. Sigrist, Superconductivity without inversion symmetry: Mnsi versus ${\mathrm{c}\mathrm{e}\mathrm{p}\mathrm{t}}_{3}\mathrm{S}\mathrm{i}$, Phys. Rev. Lett. 92, 097001 (2004a).
  * Yip [2014] S. Yip, Noncentrosymmetric superconductors, Annual Review of Condensed Matter Physics 5, 15 (2014).
  * Smidman  _et al._ [2017] M. Smidman, M. B. Salamon, H. Q. Yuan, and D. F. Agterberg, Superconductivity and spin–orbit coupling in non-centrosymmetric materials: a review, Reports on Progress in Physics 80, 036501 (2017).
  * Bauer  _et al._ [2004] E. Bauer, G. Hilscher, H. Michor, C. Paul, E. W. Scheidt, A. Gribanov, Y. Seropegin, H. Noël, M. Sigrist, and P. Rogl, Heavy fermion superconductivity and magnetic order in noncentrosymmetric ${\mathrm{c}\mathrm{e}\mathrm{p}\mathrm{t}}_{3}\mathrm{S}\mathrm{i}$, Phys. Rev. Lett. 92, 027003 (2004).
  * Yuan  _et al._ [2006] H. Q. Yuan, D. F. Agterberg, N. Hayashi, P. Badica, D. Vandervelde, K. Togano, M. Sigrist, and M. B. Salamon, $s$-wave spin-triplet order in superconductors without inversion symmetry: ${\mathrm{li}}_{2}{\mathrm{pd}}_{3}\mathrm{B}$ and ${\mathrm{li}}_{2}{\mathrm{pt}}_{3}\mathrm{B}$, Phys. Rev. Lett. 97, 017006 (2006).
  * Nishiyama  _et al._ [2007] M. Nishiyama, Y. Inada, and G.-q. Zheng, Spin triplet superconducting state due to broken inversion symmetry in ${\mathrm{li}}_{2}{\mathrm{pt}}_{3}\mathrm{B}$, Phys. Rev. Lett. 98, 047002 (2007).
  * Xu  _et al._ [2020] X. Xu, Y. Li, and C. L. Chien, Spin-triplet pairing state evidenced by half-quantum flux in a noncentrosymmetric superconductor, Phys. Rev. Lett. 124, 167001 (2020).
  * Ishihara  _et al._ [2021] K. Ishihara, T. Takenaka, Y. Miao, Y. Mizukami, K. Hashimoto, M. Yamashita, M. Konczykowski, R. Masuki, M. Hirayama, T. Nomoto, R. Arita, O. Pavlosiuk, P. Wiśniewski, D. Kaczorowski, and T. Shibauchi, Tuning the parity mixing of singlet-septet pairing in a half-heusler superconductor, Phys. Rev. X 11, 041048 (2021).
  * Cameron  _et al._ [2022] A. S. Cameron, Y. S. Yerin, Y. V. Tymoshenko, P. Y. Portnichenko, A. S. Sukhanov, M. C. Hatnean, D. M. Paul, G. Balakrishnan, R. Cubitt, A. Heinemann, and D. S. Inosov, Singlet-triplet mixing in the order parameter of the noncentrosymmetric superconductor ${\mathrm{ru}}_{7}{\mathrm{b}}_{3}$, Phys. Rev. B 105, 094519 (2022).
  * Bonalde  _et al._ [2005] I. Bonalde, W. Brämer-Escamilla, and E. Bauer, Evidence for line nodes in the superconducting energy gap of noncentrosymmetric ${\mathrm{cept}}_{3}\mathrm{Si}$ from magnetic penetration depth measurements, Phys. Rev. Lett. 94, 207002 (2005).
  * Frigeri  _et al._ [2004b] P. A. Frigeri, D. F. Agterberg, and M. Sigrist, Spin susceptibility in superconductors without inversion symmetry, New Journal of Physics 6, 115 (2004b).
  * Samokhin [2005] K. V. Samokhin, Paramagnetic properties of noncentrosymmetric superconductors: Application to ${\mathrm{c}\mathrm{e}\mathrm{p}\mathrm{t}}_{3}\mathrm{S}\mathrm{i}$, Phys. Rev. Lett. 94, 027004 (2005).
  * Gao  _et al._ [2022] Z. S. Gao, X.-J. Gao, W.-Y. He, X. Y. Xu, T. K. Ng, and K. T. Law, Topological superconductivity in multifold fermion metals, Quantum Frontiers 1, 3 (2022).
  * Bychkov and Rashba [1984] Y. A. Bychkov and É. I. Rashba, Properties of a 2D electron gas with lifted spectral degeneracy, Soviet Journal of Experimental and Theoretical Physics Letters 39, 78 (1984).
  * Manchon  _et al._ [2015] A. Manchon, H. C. Koo, J. Nitta, S. M. Frolov, and R. A. Duine, New perspectives for rashba spin–orbit coupling, Nature Materials 14, 871 (2015).
  * Hirohata  _et al._ [2020] A. Hirohata, K. Yamada, Y. Nakatani, I.-L. Prejbeanu, B. Diény, P. Pirro, and B. Hillebrands, Review on spintronics: Principles and device applications, Journal of Magnetism and Magnetic Materials 509, 166711 (2020).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * Hu  _et al._ [2011] H. Hu, L. Jiang, X.-J. Liu, and H. Pu, Probing anisotropic superfluidity in atomic fermi gases with rashba spin-orbit coupling, Phys. Rev. Lett. 107, 195304 (2011).
  * Zhou and Zhang [2012] K. Zhou and Z. Zhang, Opposite effect of spin-orbit coupling on condensation and superfluidity, Phys. Rev. Lett. 108, 025301 (2012).
  * Powell  _et al._ [2022] P. D. Powell, G. Baym, and C. A. R. Sá de Melo, Superfluid transition temperature and fluctuation theory of spin-orbit- and rabi-coupled fermions with tunable interactions, Phys. Rev. A 105, 063304 (2022).
  * Arovas  _et al._ [2021] D. P. Arovas, E. Berg, S. Kivelson, and S. Raghu, Many-body physics with ultracold gases, arXiv: 2103.12097 (2021).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Giorgini  _et al._ [2008] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of ultracold atomic fermi gases, Rev. Mod. Phys. 80, 1215 (2008).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * Kotov  _et al._ [2012] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Castro Neto, Electron-electron interactions in graphene: Current status and perspectives, Rev. Mod. Phys. 84, 1067 (2012).
  * Berezinsky [1972] V. L. Berezinsky, Destruction of Long-range Order in One-dimensional and Two-dimensional Systems Possessing a Continuous Symmetry Group. II. Quantum Systems., Sov. Phys. JETP 34, 610 (1972).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * José [2013] J. V. José,  _40 Years of Berezinskii–Kosterlitz–Thouless Theory_ (WORLD SCIENTIFIC, 2013).
  * Ryzhov  _et al._ [2017] V. N. Ryzhov, E. E. Tareyeva, Y. D. Fomin, and E. N. Tsiok, Berezinskii-Kosterlitz-Thouless transition and two-dimensional melting, Physics Uspekhi 60, 857 (2017).
  * Cheuk  _et al._ [2012] L. W. Cheuk, A. T. Sommer, Z. Hadzibabic, T. Yefsah, W. S. Bakr, and M. W. Zwierlein, Spin-injection spectroscopy of a spin-orbit coupled fermi gas, Phys. Rev. Lett. 109, 095302 (2012).
  * Wang  _et al._ [2012] P. Wang, Z.-Q. Yu, Z. Fu, J. Miao, L. Huang, S. Chai, H. Zhai, and J. Zhang, Spin-orbit coupled degenerate fermi gases, Phys. Rev. Lett. 109, 095301 (2012).
  * Huang  _et al._ [2016] L. Huang, Z. Meng, P. Wang, P. Peng, S.-L. Zhang, L. Chen, D. Li, Q. Zhou, and J. Zhang, Experimental realization of two-dimensional synthetic spin-orbit coupling in ultracold fermi gases, Nature Physics 12, 540 (2016).
  * Meng  _et al._ [2016] Z. Meng, L. Huang, P. Peng, D. Li, L. Chen, Y. Xu, C. Zhang, P. Wang, and J. Zhang, Experimental observation of a topological band gap opening in ultracold fermi gases with two-dimensional spin-orbit coupling, Phys. Rev. Lett. 117, 235304 (2016).
  * Sun  _et al._ [2018] W. Sun, B.-Z. Wang, X.-T. Xu, C.-R. Yi, L. Zhang, Z. Wu, Y. Deng, X.-J. Liu, S. Chen, and J.-W. Pan, Highly controllable and robust 2d spin-orbit coupling for quantum gases, Phys. Rev. Lett. 121, 150401 (2018).
  * Hartke  _et al._ [2023] T. Hartke, B. Oreg, C. Turnbaugh, N. Jia, and M. Zwierlein, Direct observation of nonlocal fermion pairing in an attractive fermi-hubbard gas, Science 381, 82 (2023).
  * Dell’Anna  _et al._ [2011] L. Dell’Anna, G. Mazzarella, and L. Salasnich, Condensate fraction of a resonant fermi gas with spin-orbit coupling in three and two dimensions, Phys. Rev. A 84, 033633 (2011).
  * Dell’Anna  _et al._ [2012] L. Dell’Anna, G. Mazzarella, and L. Salasnich, Tuning rashba and dresselhaus spin-orbit couplings: Effects on singlet and triplet condensation with fermi atoms, Phys. Rev. A 86, 053632 (2012).
  * Sun  _et al._ [2013] Q. Sun, G.-B. Zhu, W.-M. Liu, and A.-C. Ji, Spin-orbit coupling effects on the superfluidity of a fermi gas in an optical lattice, Phys. Rev. A 88, 063637 (2013).
  * Devreese  _et al._ [2014] J. P. A. Devreese, J. Tempere, and C. A. R. Sá de Melo, Effects of spin-orbit coupling on the berezinskii-kosterlitz-thouless transition and the vortex-antivortex structure in two-dimensional fermi gases, Phys. Rev. Lett. 113, 165304 (2014).
  * Devreese  _et al._ [2015] J. P. A. Devreese, J. Tempere, and C. A. R. Sá de Melo, Quantum phase transitions and berezinskii-kosterlitz-thouless temperature in a two-dimensional spin-orbit-coupled fermi gas, Phys. Rev. A 92, 043618 (2015).
  * Lee and Kim [2017] J. Lee and D.-H. Kim, Induced interactions in the bcs-bec crossover of two-dimensional fermi gases with rashba spin-orbit coupling, Phys. Rev. A 95, 033609 (2017).
  * Shi  _et al._ [2016] H. Shi, P. Rosenberg, S. Chiesa, and S. Zhang, Rashba spin-orbit coupling, strong interactions, and the bcs-bec crossover in the ground state of the two-dimensional fermi gas, Phys. Rev. Lett. 117, 040401 (2016).
  * Tang, Ho-Kin  _et al._ [2014] Tang, Ho-Kin, Yang, Xiaosen, Sun, Jinhua, and Lin, Hai-Qing, Berezinskii-kosterlitz-thoules phase transition of spin-orbit coupled fermi gas in optical lattice, EPL 107, 40003 (2014).
  * Rosenberg  _et al._ [2017] P. Rosenberg, H. Shi, and S. Zhang, Ultracold atoms in a square lattice with spin-orbit coupling: Charge order, superfluidity, and topological signatures, Phys. Rev. Lett. 119, 265301 (2017).
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Monte carlo calculations of coupled boson-fermion systems. i, Phys. Rev. D 24, 2278 (1981).
  * Hirsch [1983] J. E. Hirsch, Discrete hubbard-stratonovich transformation for fermion lattice models, Phys. Rev. B 28, 4059 (1983).
  * White  _et al._ [1989a] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis, and R. T. Scalettar, Numerical study of the two-dimensional hubbard model, Phys. Rev. B 40, 506 (1989a).
  * He  _et al._ [2019] Y.-Y. He, M. Qin, H. Shi, Z.-Y. Lu, and S. Zhang, Finite-temperature auxiliary-field quantum monte carlo: Self-consistent constraint and systematic approach to low temperatures, Phys. Rev. B 99, 045108 (2019).
  * Wu and Zhang [2005] C. Wu and S.-C. Zhang, Sufficient condition for absence of the sign problem in the fermionic quantum monte carlo algorithm, Phys. Rev. B 71, 155115 (2005).
  * McDaniel  _et al._ [2017] T. McDaniel, E. F. D’Azevedo, Y. W. Li, K. Wong, and P. R. C. Kent, Delayed Slater determinant update algorithms for high efficiency quantum Monte Carlo, The Journal of Chemical Physics 147, 174107 (2017).
  * Scalettar  _et al._ [1991] R. T. Scalettar, R. M. Noack, and R. R. P. Singh, Ergodicity at large couplings with the determinant monte carlo algorithm, Phys. Rev. B 44, 10502 (1991).
  * Assaad and Evertz [2008] F. Assaad and H. Evertz, World-line and determinantal quantum monte carlo methods for spins, phonons and electrons, in  _Computational Many-Particle Physics_, Lecture Notes in Physics, Vol. 739, edited by H. Fehske, R. Schneider, and A. Weiße (Springer Berlin Heidelberg, 2008) pp. 277–356.
  * Chang  _et al._ [2015] C.-C. Chang, S. Gogolenko, J. Perez, Z. Bai, and R. T. Scalettar, Recent advances in determinant quantum monte carlo, Philosophical Magazine 95, 1260 (2015).
  * Yang [1962] C. N. Yang, Concept of off-diagonal long-range order and the quantum phases of liquid he and of superconductors, Rev. Mod. Phys. 34, 694 (1962).
  * Astrakharchik  _et al._ [2005] G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Momentum distribution and condensate fraction of a fermion gas in the bcs-bec crossover, Phys. Rev. Lett. 95, 230405 (2005).
  * He  _et al._ [2022] Y.-Y. He, H. Shi, and S. Zhang, Precision many-body study of the berezinskii-kosterlitz-thouless transition and temperature-dependent properties in the two-dimensional fermi gas, Phys. Rev. Lett. 129, 076403 (2022).
  * Strinati  _et al._ [2018] G. C. Strinati, P. Pieri, G. Röpke, P. Schuck, and M. Urban, The bcs–bec crossover: From ultra-cold fermi gases to nuclear systems, Physics Reports 738, 1 (2018), the BCS–BEC crossover: From ultra-cold Fermi gases to nuclear systems.
  * White  _et al._ [1989b] S. R. White, D. J. Scalapino, R. L. Sugar, N. E. Bickers, and R. T. Scalettar, Attractive and repulsive pairing interaction vertices for the two-dimensional hubbard model, Phys. Rev. B 39, 839 (1989b).
  * Paiva  _et al._ [2004] T. Paiva, R. R. dos Santos, R. T. Scalettar, and P. J. H. Denteneer, Critical temperature for the two-dimensional attractive hubbard model, Phys. Rev. B 69, 184501 (2004).
  * Fontenele  _et al._ [2022] R. A. Fontenele, N. C. Costa, R. R. dos Santos, and T. Paiva, Two-dimensional attractive hubbard model and the bcs-bec crossover, Phys. Rev. B 105, 184502 (2022).
  * Tomita and Okabe [2002] Y. Tomita and Y. Okabe, Probability-changing cluster algorithm for two-dimensional $\mathrm{XY}$ and clock models, Phys. Rev. B 65, 184405 (2002).
  * Nguyen and Boninsegni [2021] P. H. Nguyen and M. Boninsegni, Superfluid transition and specific heat of the 2d x-y model: Monte carlo simulation, Applied Sciences 11, 10.3390/app11114931 (2021).