Effective p-wave Fermi-Fermi Interaction Induced by Bosonic Superfluids

We study the two-dimensional Bose-Fermi mixture on square lattice at finite temperature by using the determinant quantum Monte Carlo method within the weakly interacting regime. Here we consider the attractive Bose-Hubbard model and free spinless fermions. In the absence of boson-fermion interactions, we obtain the boundary of the collapsed state of the attractive bosons. In the presence of boson-fermion interactions, an effective p-wave interaction between fermions will be induced as far as the bosons are in a superfluid state. Moreover, we find the emergence of the composite fermion pairs at low temperatures.

bosons and fermions. Throughout this paper, we take t b = t f = 1 as the energy unit. We employ the finite temperature determinant Monte Carlo (DQMC) method [20][21][22][23][24][25] for numerical simulations. It is a field-theoretic method where many-body propagators resulting from two-body interactions can be transformed into one-body propagators by using the Hubbard-Stratonovich (HS) transformation. The resulting integrals can then be computed by Monte Carlo sampling. In ref. 26 , the DQMC method was first used to study Bose-Fermi mixtures via a combination of bosonic and fermionic Monte Carlo techniques. Though some exact results can be obtained in the small size lattices, the sign problem is still quite severe in most parameter ranges (strong repulsive interaction in the large size lattices at low temperature). Interestingly, we find that there are two cases where the sign of the determinant is always equal to one. One example is the attractive Bose-Hubbard model at low particle density, whose Hamiltonian can be acquired by ignoring the whole fermionic part of the Eq. (1). Another case is in the Bose-Fermi mixture with the interaction between bosons and fermions being weaker than that of the bosons, namely ≤ U U bf b and U b ≤ 0. On the other hand, a repulsion between bosons and fermions prefers a demixing to minimize the overlapping region, whereas in the case of an attraction the mixture can collapse, as long as the particle numbers are sufficiently large or interspecies interaction is strong 27,28 . In the following, we mainly focus on dilute Bose-Fermi mixture with weak interspecies interaction where the sign of the corresponding determinants in these two cases which are used to express the trace over one-body propagators is always positive in our calculation (Monte Carlo samples up to 10 4 ) so that the many-dimensional integrals can be performed accurately by the Monte Carlo sampling.

Results
Bose-Hubbard model with attractive interaction. In this section, we investigate the Bose-Hubbard model by considering only the bosonic part of the Hamiltonian in Eq. (1). According to the previous studies on the atomic Bose gas with negative scattering lengths, the number of the bosons would be quite limited due to the collapse of the BEC [29][30][31] . More specifically, the maximum of the boson number N c of a system under experiment is inversely proportional to the negative scattering length a s , that is N c · a s = constant. Here we choose the average particle density n b of a square lattice system with size 6 × 6 to be 0.02 ± 0.0006 which could be realized within experimental reach 12 . The on-site attractive interactions are set to be U b = −0.1, −0.2, −0.3.
Off-diagonal long-range order of bosons. Here we study the superfluidity of bosons at low temperature. The off-diagonal correlation function † g r bb ( ) i j = 32,33 , where = − r i j is the distance between two lattice site i and j, has been investigated. We rewrite the correlation function as g r for a better comparison. The calculations are performed for several different inverse temperatures with lattice size 4 × 4, 6 × 6, 8 × 8, and the corresponding results are presented in Fig. 1(a-c) respectively. Note that the bosons in our system are quite dilute. As a result, the on-site interactions have little effects on the density distribution or the off-diagonal correlation, and we only present results with U b = −0.2. From the off-diagonal correlation function shown in Fig. 1(a-c), we can observe that at the high temperature such as β = 0.5, the off-diagonal correlation goes down to zero quickly. While as the temperature is progressively reduced, the decay of the correlation become more and more slowly with a finite value at the longest distance r = r max which indicates that the finite size system has entered the superfluid state.
To show finite size effect, finite size extrapolations of the off-diagonal long-range order g(r max ) has been done for lattice size 4 × 4, 6 × 6, 8 × 8 at various inverse temperature. The results are shown in Fig. 1(d). For high temperature β = 0.5, 1.0, 2.0, the off-diagonal long-range order decay quickly to zero. While a finite value g(rmax ) ≈ 0.014 can be obtained by doing linear fitting at β = 5.0. To verify the superfluidity of bosons in the thermodynamic limit, larger system sizes are need to be considered.
Critical value of the collapsed state at finite temperature. The influence of a larger number of bosons, which can be tuned by increasing the chemical potential ε b , can also be analyzed by the numerical method. Interestingly, the particle densities at certain lattice sites increase abruptly while densities at other sites are decreased. This phenomenon is consistent with the collapsed state observed in previous experimental studies 31 . To better characterize the collapsed state, we introduce the inverse participation ratio (IPR) α 34 , which is usually adopted as the criteria for localization in numerical studies. The IPR is defined as: www.nature.com/scientificreports www.nature.com/scientificreports/ where α i is the particle density at site i, and N is the number of lattice sites. α has a small value of O(N −1 ) when the system is in a uniform state. However, it will reach 1 when all the particles are localized at a certain lattice site. In the inset of Fig. 2, the variation of α as the chemical potential ε b changed is presented at the inverse temperature of  The critical value of density vs the absolute value of attractive interaction U b for various temperatures with lattice size 6 × 6. Obviously, the region of the collapse state enlarges as the temperature goes down. The inset displays the participation ratio of density versus the chemical potential ε b for several attractive interaction at inverse temperature β = 0.5. Note that normally the density of boson increases linearly as the the chemical potential decreases.
www.nature.com/scientificreports www.nature.com/scientificreports/ β = 0.5. It can be found that when the system is in the normal state, the IPR is a constant value α = 1/36. While if the system enters the collapsed state in which the system becomes unstable, the IPR will suddenly jump to a much larger value closed to 1. For example, when U b = −0.3, as indicated by the blue dotted line in the inset, α is quite close to 1. As we suppress the attractive interaction between the bosons, the abrupt change of IPR still shows up as the system goes into the collapsed state (see the dotted lines corresponding to U b = −0.2 and −0.1, respectively). We have also calculated the IPR of the system at several other temperatures (not shown here), and find that when the system undergoes a phase transition between a normal state and the collapsed state, the IPR will always show a sudden jump at the critical value.
In Fig. 2, we also show the critical value of the particle number density n c which corresponds to the onset of the sudden jump of IPR for different temperatures β = 0.25, 0.5 and 3.0. It can be observed that the critical value will be reduced as we enhance the attractive interaction between the bosons. Besides, as the temperature is lowered down (i.e., with β increasing from 0.25 to 3.0), the reduction of the critical value will become sharper. This is consistent with the result from the previous study in which the equation describes the connection between the critical value and the interaction strength , with k being a constant which decreases for lower temperatures. A lattice of size 4 × 4 has also been investigated, and similar results are obtained (results not shown here). We can conclude that the regime of the collapsed state will be enlarged as the temperature goes down.
Bose-Fermi mixture with spinless fermions. Setup. For the sake of simplicity, we consider the spinless fermion only in the Bose-Fermi mixture system. To avoid the sign problem, we set the interaction strength between the bosons and fermions U bf to be equal to the interaction between bosons U b unless otherwise stated. Here we mainly concentrate on the square lattice system with size N = 6 × 6 in which densities of bosons and fermions are tuned to be = . ± . , −0.1 and −0.2. It is noteworthy that the boson-fermion interaction will change from U bf to − U bf under particle-hole transformation, which can be shown by changing the operator 16 .
Charge fluctuation and p-wave pair correlation. Due to the coherence of the superfluid bosons and interaction between boson and fermion, two fermions can affect each other via the bosons in the mixture. To study the influences of the superfluidity on the fermions, we calculate the equal time charge fluctuation of the fermions and the density fluctuation between bosons and fermions. The numerical results for systems with various coupling strengths are presented in Fig. 3. Generally speaking, the charge fluctuation of fermions is given by Note that the angle bracket denotes both the quantum mechanical average over the whole space-time configuration and the statistical average of the samples (HS field configurations) 24 . To reduce the statistical errors, we redefine charge fluctuation as: f More interestingly, due to the interaction between the bosonic and fermionic components in the Bose-Fermi mixture, the effective p-wave pairing interactions between fermions could be induced. The p-wave pair correlation is defined as corresponds to the lattice site involving the p-wave pairing. Similar to the definition of the charge fluctuation, we can also redefine the p-wave pair correlation as where P(r) 0 is the uncorrelated pair correlation 36,37 . Namely, for each g g f f 〈 〉 in P(r), there will be a term of the form 〈 〉〈 〉 g g f f . We have calculated the pair correlation along both the x and y-direction in our numerical calculations and the results are shown in Fig. 4. P r ( ) becomes weaker as the distance increases in both directions. The sign of U bf has no significant effect on the results. From Fig. 4(a), we can also find that at the same distance, a stronger coupling between bosons and fermions will lead to a stronger p-wave pairing and the sign of P r ( ) will not change. However, in the x-direction (see Fig. 4(b)), the correlation will change sign at the same distance while the strength of U bf changes. The conclusion that stronger coupling between bosons and fermions will result in stronger p-wave pair correlation still holds along the x-direction. So the pair correlation will always be zero if there is no coupling between the bosons and fermions but will become finite if U bf is turned on. We can conclude that the superfluid bosons mediate the effective p-wave pairing between the fermions.
At this stage, we need to address whether the superfluidity of the boson is a necessary condition for the induction of p-wave correlation of fermions. According to the results shown in Fig. 1, the bosons enter the superfluid state, which is demonstrated to be the second order phase transitions, with the temperature goes down. This www.nature.com/scientificreports www.nature.com/scientificreports/ superfluid state often can be described by the off-diagonal long-range order. Here we consider the g(r max ) as the long-range order parameter of the bosons with = r 3 2 max for system sites 6 × 6. Meanwhile, we take the p-wave correlation P r ( ) max as the off-diagonal long-range order parameter for the p-wave of the fermion with r max = 3 in the y-direction. For clarification, P r ( ) max of the different interaction strengths U bf are rescaled to the U bf = −0.1 by dividing | | U bf . And g(r max ) is divided by 1200. These results are shown in Fig. 5. Note that the weak Bose-Fermi couplings have little effect on g(r max ), so we only show the g(r) of U bf = −0.2 in the figure. It is obvious to find that the boson goes into the superfluid state as the temperature goes down shown by the real olive line. Due to the finite size effect, g(r max )/1200 turns into zero smoothly. It is important that the P r ( ) max of the different interaction strengths U bf (dashed line) show the almost same trend with g(r max ), which indicate that the superfluid state of the bosons is a prerequisite for the p-wave correlation of the fermions.
Both the charge fluctuation and the p-wave pair correlation indicate that an effective interaction 38 can be induced by the superfluid bosons. To further understand this effective interaction, we use the DQMC method to calculate a simplified spinless fermion model whose Hamiltonian is given by: The V eff here represents the effective interaction between the nearest neighboring fermions. Since the boson-fermion coupling is weak, we can neglect the interactions between fermions separated by larger distances and only consider the nearest-neighboring situation. Generally, the spinless fermion model is sensitive to the sign problem at low temperature except the positive definite case 39 . Fortunately, we just care about the weak interaction region (V eff ∈ [−0.5, 0]) at the low density = . m ( 0 2 ) where the sign is always positive with samples up to 10 4 . For a better comparison, we introduce the ratio r P r P r ( ) ( ) / ( ) where P r ( ) bf and P r ( ) eff are the p-wave pair correlation at distance r in the mixture and spinless fermion system respectively. It is similar to the ratio of charge fluctuation r C r C r ( ) The numerical results about the C r ( ) and P r ( ) for lattice sizes 6 × 6 and 8 × 8 are shown in Fig. 6. Surprisingly, both the ratios of charge fluctuation and p-wave pair correlation γ(r) p/c at | | = .
. . What's more, the strength of effective interaction is a half of the Bose-Fermi interaction, which is much different form the induced effective interaction in 3D Bose-Fermi mixture 12 .
Finite size extrapolations for P-wave pair correlation function. As the temperature goes down, the p-wave pair correlation decays slowly shown in the inset of Fig. 7. It is obvious that the lower of the temperature, the slower decay of the p-wave pair correlation. And a finite value at the longest distance r = r max = 3, which is served as the off-diagonal long range order of p-wave pair correlation function, can be observed at low temperature β = 2.0, 5.0 for lattice size 6 × 6.
To investigate the finite size effect, the p-wave pair correlation functions for different lattice sizes 4 × 4, 6 × 6, 8 × 8 have been calculated at various inverse temperatures.
The results of p-wave pair correlation function P r ( ) max have been shown in Fig. 7. It is observed that the p-wave pair correlation function decays quickly to zero as a function of lattice size at the high temperature β = 2.0, 1.0, 0.5, while appears a finite value at low temperature β = 5.0, 3.0. To verify the superfluidity of the fermion unquestionably in the thermodynamic limit, larger size and low temperature are still needed.
Composite pair correlation function. In this subsection, we concentrate on the correlations between the bosons and fermions which can be characterized by the composite pair Greens function. A general definition of the www.nature.com/scientificreports www.nature.com/scientificreports/

Greens function is
with Δ(i) = b i c i . Due to the particle-hole transformation, the sign of G bf changes when the interaction U bf is changed to −U bf . For the sake of clarity, we redefine the composite pair correlation as: bf bf j i j i By using this definition, the composite pair correlation is equal to zero when there is no interaction between the bosons and fermions.  www.nature.com/scientificreports www.nature.com/scientificreports/ The numerical results of G r ( ) bf is presented in Fig. 8. It is obvious that the stronger interaction between bosons and fermions may result in stronger composite pair correlation. However, the evolution of G r ( ) bf as a function of r shows a non-monotonic behavior. The values of G r ( ) bf for different U bf intersects and vanishes around r = 2.5 and reaches the negative maximum around r = 3 and then become much suppressed as r further increases. As the inset displayed, the correlation becomes smeared out by increasing temperatures.

Discussion
In conclusion, we have studied the 2D Bose-fermi mixture on square lattice by using the DQMC method at a finite temperature in the ranges where the sign problem can be ignored. Our results show the boundary of the collapsed state at various temperatures in the attractive Bose-Hubbard model. More interestingly, when the interaction between bosons and fermions turns on, the off-diagonal long-range order of composite fermions pair shows a fairly flat finite value at a distance r = 2.5 at low temperature, which denotes that the composite fermi pairs exist in the mixture. Meanwhile, we also have observed a slowing decay p-wave pair function due to the effective attractive interaction induced by the superfluid boson, which denotes that the p-wave superconductivity may be detected in the Bose-Fermi mixture. This super-pairing mechanism is different from the conventional BCS pairing mechanism in which effective interactions are induced by exchanging massless phonons. Moreover, to overcome finite size effect, finite size extrapolations have been done for p-wave pair correlation function at various inverse temperatures. It is obvious that a larger finite value can be observed as the temperature goes down. To verify the superfluid of the fermion unquestionably in the thermodynamic limit, larger size system and lower temperature are still needed to do finite size scaling. None zero off-diagonal long range order may be observed in larger system near zero temperature. While there is a notorious sign problem in the DQMC method when the particle number is away from half-filling in the large lattice size system at low temperature. In this case, the results will drastically deviate form the exact solutions. More effort is needed to obtain the results in the thermodynamic limit.
By the way, spinless non-interaction fermions embedded into superfluid boson have been considered in this work. Here the superfluidity of bosons play a vital role to induce the effective interaction between fermions and non-zero p-wave pair correlation have been observed. It is easy to imagine that the s-wave pairing or d-wave pairing may be observed in the two components non-interaction Fermi system which could be realized in the experiment. On the other way, bosons with mediated interaction embedded into the superfluid fermions is also worth to consider in the experiment. The strong correlated fermions may induce the effective interaction between bosons. Then both of them may enter into an exotic superfluid state.