Temperature-driven BCS-BEC crossover in a coupled boson-fermion system

TL;DR

This paper introduces a 2D bose-fermi model demonstrating a temperature-driven BCS-BEC crossover, with implications for high-temperature superconductors, graphene, and cold atomic systems, using Monte Carlo simulations.

Contribution

It presents a simple coupled boson-fermion model showing temperature-dependent interactions and a novel BCS-BEC crossover observed through spectral changes.

Findings

Temperature-dependent long-range interactions between bosons.

Observation of a temperature-induced BCS-BEC crossover.

Spectral shift from Fermi wave vector to Gamma point.

Abstract

We propose a simple bose-fermi model in two dimensions, with a coupling that converts pairs of opposite spin fermions into localized bosons and vice versa. We show that tracing out one of the degrees, either the bosons or fermions, generates temperature-dependent long range effective interactions between bosons as well as effective attractive interactions between fermions. Using Monte Carlo techniques we obtain the thermodynamic properties and phase stiffness as a function of temperature, dominated by vortex-antivortex unbinding of the bosons. Remarkably in the fermion sector we observe a temperature-induced BCS-BEC crossover signaled by a distinct change of their spectral properties: the minimum gap locus moves from the Fermi wave vector to the $\Gamma$ point. Such a model is relevant for describing aspects of high $T_c$ superconductivity in cuprates and pnictides, superconducting islands on graphene, and bose-fermi mixtures in cold atomic systems.

Equations

$$H_{\mathrm{BF}}$$

$$H_{\mathrm{PF}}=-t\sum _{⟨ij⟩,\sigma }{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+\tilde{g}\sum _i\left(e^{i{\theta }_i}{\hat{c}}_{i\uparrow }{\hat{c}}_{i\downarrow }+H.c.\right),$$

$$Z=\sum _{\{{\theta }_i\}}\mathrm{Tr}{e}^{-\beta H_{\mathrm{PF}}(\vec{\theta })}=\sum _{\{{\theta }_i\}}{e}^{-\beta \mathcal{H}(\vec{\theta })},$$

$$\mathcal{H}(\vec{\theta })=-\frac{1}{\beta }\sum _n\ln \left[1+e^{-\beta E_n(\vec{\theta })}\right],$$

$$V(\Delta\theta,T)\to\left\{\begin{array}[]{ll}A(T)(\Delta\theta)^{2}&\mbox{ for }0\leq\Delta\theta<\pi,\\[8.61108pt] A(T)(\Delta\theta-2\pi)^{2}&\mbox{ for }\pi\leq\Delta\theta<2\pi,\end{array}\right.$$

$$V(\Delta \theta ,T)\to A(T)\cos (\Delta \theta )$$

$${\rm Y}(\beta )=\frac{2}{\beta {\pi }^2}{\int }_0^{\beta }\left[⟨\bar{\mathcal{E}}(\pi ){⟩}_{{\beta }^{\prime }}-⟨\bar{\mathcal{E}}(0){⟩}_{{\beta }^{\prime }}\right]d{\beta }^{\prime },$$

$$⟨\bar{\mathcal{E}}(\varphi ){⟩}_{\beta }\equiv \frac{1}{Z}\sum _{\{{\theta }_i\}}\bar{\mathcal{E}}\left(\vec{\theta },\varphi \right)e^{-\beta \bar{\mathcal{E}}\left(\vec{\theta },\varphi \right)}$$

$$\bar{\mathcal{E}}\left(\vec{\theta },\varphi \right)\equiv \sum _n\left[E_n(\vec{\theta },\varphi )-\mu \right]f\left[E_n(\vec{\theta },\varphi )-\mu \right].$$

$${{\rm Y}}_L={{\rm Y}}_{L\to \infty }\left[1+\frac{1}{2}\frac{1}{\ln (L)+C}\right],$$

$${\left\{\frac{2}{\pi }{\int }_0^{\beta }\left[⟨\bar{\mathcal{H}}(\pi ){⟩}_{{\beta }^{\prime }}-⟨\bar{\mathcal{H}}(0){⟩}_{{\beta }^{\prime }}\right]d{\beta }^{\prime }-2\right\}}^{-1}-\ln (L)=C.$$

$$\mathcal{G}({\mathbf{R}}_i,{\mathbf{R}}_j,z)={\left\{z-\mathcal{H}(\vec{\theta })\right\}}_{ij}^{-1},$$

$$A(\mathbf{k},\omega )=-\frac{1}{\pi }\mathrm{Im}G\left(\mathbf{k},\omega +i{0}^{+}\right),$$

$$\sum _{{\mathbf{R}}_i}\sum _{{\mathbf{R}}_j}{e}^{i\left(\mathbf{k}\cdot {\mathbf{R}}_i-{\mathbf{k}}^{\prime }\cdot {\mathbf{R}}_j\right)}⟨\mathcal{G}({\mathbf{R}}_i,{\mathbf{R}}_j,z)⟩=G(\mathbf{k},z)\delta (\mathbf{k}-{\mathbf{k}}^{\prime }).$$

$$\Delta \equiv \frac{1}{N}\sum _i⟨c_{i\uparrow }{c}_{i\downarrow }⟩$$

$$\mathrm{IPR}=\frac{{\sum }_{ni}|{\Psi }_{ni}{|}^4}{{\sum }_{ni}|{\Psi }_{ni}{|}^2},$$

$$H({\theta }_1,{\theta }_2)=$$

$$E_1$$

$$E_2$$

$$E_3$$

$$E_4$$

$$Z=\int d{\theta }_2\int d{\theta }_2\mathrm{Tr}{e}^{-\beta H({\theta }_1,{\theta }_2)}=2\pi \int d\Delta \theta e^{\beta V(\Delta \theta ,\beta )},$$

$$V(\Delta \theta ,\beta )=-\frac{1}{\beta }\sum _{i=1}^4\ln \left[1+e^{-\beta E_i(\Delta \theta )}\right]$$

Full text

I Introduction

Strongly interacting systems lead to emergent new phases with spontaneously broken symmetries, such as magnets with broken time reversal symmetry and superconductors with broken gauge symmetry. Even within the superconducting phase the system can show effects of interactions evolving from a BCS regime with large Cooper pairs compared to inter-particle spacing to a strongly coupled regime where the Cooper pairs are tightly bound Randeria and Taylor (2014). Such an evolution of a system from a BCS to BEC regime has been predicted theoretically as a function of increasing pairing interaction and observed in experiments Ketterle and Zwierlein (2008). In this paper we propose such a BCS to BEC crossover at a fixed coupling strength but with increasing temperature. We believe this is the first such prediction originating from an explicit calculation of a model.

The original motivation for our model comes from recognizing that strongly interacting systems, such as fractional quantum Hall effect, frustrated magnets and high–$T_{c}$ superconductors, can be described in terms of emergent degrees of freedom that interact via fluctuating gauge fields Fradkin (2013). In high temperature superconductors, for example, strong on–site Mott interactions generate effectively a two-component system in which the fermionic holes become superconducting in a matrix of bosonic fluctuations of the spin singlets. While many issues are still hotly debated, such as the role of intertwined order, it is nevertheless remarkable that the broad phenomenology can be understood in terms of a two-component response: (a) the hole density determining superconducting transition temperature $T_{c}$ and (b) the spin singlets generating the pseudogap scale $T^{\ast}$ below which a soft gap opens up in the density of states Anderson et al. (2004).

Below, we investigate the spectral properties of a simple boson-fermion (BF) model as a function of coupling and temperature. The BF model, proposed in the 1950s in the context of superconductivity Schafroth (1954), describes itinerant fermions hybridizing with bosons composed of pairs of tightly bound opposite–spin fermions. The idea reoccurred in 1980s in the context of electrons interacting with local lattice deformations Ranninger and Robaszkiewicz (1985) and high temperature superconductivity Ranninger et al. (1995); Friedberg and Lee (1989); Geshkenbein et al. (1997); Domański et al. (2003); Maśka and Mierzejewski (2001). Recently, the boson–fermion model has been adopted to describe resonance superfluids in the BCS–BEC crossover regime Shin et al. (2008). The model is rather general and can be applicable to a host of other realizations such as graphene with superconducting islands, and bose-fermi mixtures in cold atoms.

Model: The Hamiltonian of the BF model is given by

[TABLE]

where $\hat{c}_{i\sigma}$ ($\hat{b}_{i}$) are fermionic (bosonic) annihilation operators, $n^{f}$ ($n^{b}$) is the concentration of fermions (bosons), $\mu$ is the chemical potential, and $E_{B}$ is the bosonic level. In the following we tune the model parameters $\mu$ and $E_{B}$ to have the average value $n^{f}=1$. We set the hopping integral $t$ as the energy unit ($t=1$). The bosons are assumed to be localized, so their kinetic energy is neglected. We assume that the number of bosons per lattice site is large which allows us to approximate $\hat{b}_{i}\>\rightarrow\>\sqrt{n_{i}^{\rm b}}e^{i\theta_{i}}$. We then fix the number of bosons, but fermions are in the grand canonical ensemble with the chemical potential used to control their concentration.

The BF model does not include a direct interaction between bosons. However, through the interaction with mobile fermions, effective boson–boson interaction gets mediated. Correspondingly, the properties of the fermions are also affected by their coupling to the bosons. Since the state of the bosons is strongly temperature dependent, the fermion–boson interaction leads to nontrivial behavior of the fermions, most strikingly, a temperature-induced BCS–BEC crossover. Fig. 1 shows the main idea of the proposed approach.

To perform a numerical study of the model in Eq. (1) we trace out the fermionic degrees of freedom, and obtain an effective boson–boson interacting Hamiltonian. We next determine the thermodynamics of this model using Monte Carlo (MC) simulations. What distinguishes this problem from standard MC applications is that each MC step for the bosons requires diagonalizing the fermionic Hamiltonian. Such a fully quantum problem is restricted to small systems and does not allow us to access the phase transitions in the bosonic sector. This explains why we approximate the bosonic operators by $c$-numbers $\hat{b}_{i}\>\rightarrow\>\sqrt{n_{i}^{\rm b}}e^{i\theta_{i}}$ so that quantum fluctuations of the bosonic fields are neglected while still retaining the classical phase fluctuations.

For the sake of simplicity, we assume that the boson density is uniform $n_{i}^{\rm b}=n^{\rm b}$, however, this assumption is not crucial and density fluctuations can be taken into account within the proposed approach (using, e.g., grand–canonical MC methods Adams (1974, 1975)). The resulting Hamiltonian describes itinerant fermions interacting with classical phases $\theta_{i}$ [phase–fermion (PF) model]:

[TABLE]

where the effective coupling includes the density of bosons $\tilde{g}\equiv g\sqrt{n^{\rm b}}$.

Experimental Systems: We propose systems that can be directly described by the effective Hamiltonian (2).

(a) The first one is graphene decorated with an array of superconducting islands Kessler et al. (2010); Han et al. (2014); Eley et al. (2012). For islands sufficiently close to each other, Cooper pairs can tunnel directly between them and the system can be described as an array of Josephson junctions. In such a system as the temperature drops below a Berezinsky-Kosterlitz-Thouless (BKT) transition temperature Berezinskii (1971); Kosterlitz and Thouless (1973), the phases on the islands get locked and become phase coherent. The BF model suggests that even if the distance between the islands is larger than the range of the superconducting proximity effect, phase ordering in the entire array could occur mediated by the (normal) carriers in graphene that are sensitive to the phases of the superconducting islands. Because of the relatively large size of the superconducting islands, the fluctuations of the number of Cooper pairs can be neglected and the PF model (2) is clearly applicable.

(b) Another example is to take a 3D Bose-Einstein condensate of tightly bound Cooper pairs (molecules) and break it up into a 2D array of 1D tubes, by using an optical lattice applied along two dimensions. For a sufficiently deep lattice potential, the phases between the tube-BECs are uncorrelated, hence the system as a whole is insulating. As the lattice depth is decreased, tunneling between the tubes can drive long range phase coherence across the entire system. In addition to the bosons, if unpaired fermions are present, they will experience the phase fluctuations arising from the tube-BECs. If the lattice potential can localize the molecules but allow the unpaired fermions to tunnel between the BECs, the latter can mediate the interaction between the bosons leading to long range coherence of the isolated bosons.

II Bosonic Sector:

Effective classical Hamiltonian

The Hamiltonian (2) represents a system described by both quantum (itinerant fermions) and classical (phases) degrees of freedom. We generalize the approach in Ref. Maśka and Czajka (2006) proposed for the Falicov-Kimball model to study the thermodynamics of the system here. The partition function of the PF model is a sum over all possible configurations of phases and a trace over the fermionic degrees of freedom

[TABLE]

where the second term is obtained after numerically tracing out the fermions exactly for a given configuration of phases $\vec{\theta}\equiv(\theta_{1},\ldots,\theta_{N})$. Here ${\cal H}$ is a classical temperature–dependent “Hamiltonian”

[TABLE]

and $E_{n}(\vec{\theta})$ are its single-particle eigenvalues for a given set of $\theta$’s. Note that ${\cal H}(\vec{\theta})$ is the free energy of the fermionic subsystem. The model explicitly takes into account the time scales separating fast fermion dynamics, treated in a fully quantum way, from the slow phase dynamics described by the classical degrees of freedom ${\theta_{i}}$.

Some of the questions we address about the PF model in the bosonic sector are: (1) What are the effective interactions between bosons generated by tracing out the fermions? (2) Is there a phase transition from a low temperature ordered state of boson phases to a high temperature disordered state? (3) Is the transition described by vortex-antivortex unbinding and is it in the BKT universality class? (4) Are there any differences in the phases or the transitions in the PF model and the standard BKT transition in the XY model with nearest neighbour interactions? In the discussion below, we answer these questions.

A simple analysis of two sites (see Appendix A) shows that in the weak coupling regime the effective interaction potential between bosons or phases deviates from the expected $\cos(\Delta\theta)$ of the usual XY model not only in the functional form but also in that the interaction $V$ is strongly temperature dependent. We present the limiting behavior for $V$ here (the full interaction potential is given in the Appendix A), given by

[TABLE]

at low temperature and

[TABLE]

at high temperature. The parameter $A(T)\to 0$ for $T/t\to\infty$ indicating that the interaction strength decreases with temperature (see Fig. 8). For stronger coupling $\tilde{g}$ the angle dependence is very well described by $\cos(\Delta\theta)$ though the coupling strength is strongly temperature dependent.

To identify the phase transition we follow the original approach by Fisher Fisher et al. (1973) and calculate the response of the system to a phase twist. Specifically, we fix $\theta_{i}=0$ at the left edge of the system and $\theta_{i}=\pi$ at the right one and obtain the superfluid stiffness $\Upsilon$ from

[TABLE]

where $\langle\ldots\rangle$ denotes an average over the phase configurations $\vec{\theta}$ generated in MC sampling

[TABLE]

and $\bar{\mathscr{E}}(\vec{\theta},\phi)$ is the usual quantum mean value of the Hamiltonian for a given configuration $\vec{\theta}$

[TABLE]

$E_{n}(\vec{\theta},\phi)$ are eigenvalues for a difference of phases on opposite edges equal to $\phi$ and $f(\ldots)$ is the Fermi–Dirac distribution function.

Fig. 2a shows the increase of phase fluctuations with increasing temperature, the formation of vortex-antivortex bound pairs (Fig. 2b), and their unbinding at a higher temperature (Fig. 2c).

Figure 3a shows the behavior of the superfluid stiffness $\Upsilon$. With increasing system size $\Upsilon$ gets steeper, but the systems are still too small to display a signature of the universal jump in the stiffness. On the other hand, in the XY model $\Upsilon$ converges rapidly at low temperature even for very small systems Hsieh et al. (2013). The strong size dependence of $\Upsilon$ in the PF model at low temperatures can be attributed to the fact that the inter–boson coupling is long–ranged as it is generated by the itinerant fermions that travel through the entire system.

Given the stronger size dependence and the inability to observe a jump of the stiffness, we explore a different method to locate $T_{c}$. The stiffness at the BKT transition $T_{c}$ in the thermodynamic limit is given by Kosterlitz (1974) $\Upsilon_{L\to\infty}=\frac{2}{\pi}T_{c}$. From the knowledge of the Kosterlitz renormalization-group scaling at $T_{c}$ Weber and Minnhagen (1988),

[TABLE]

where $C$ is a constant given by:

[TABLE]

Since $C$ does not depend on the system size, if one plots the l.h.s. of Eq. (11) [denoted as ${\cal A}_{L}(T)$] as a function of temperature, different sizes cross at the same point for $T=T_{c}$, as shown in Fig. 3b. The vanishing of the fitting errors and crossings ${\cal A}_{L}(T)$ for all $L$ at $T=T_{c}$ indicates that the itinerant–fermions-mediated interaction between the classical phases indeed drives a phase transition in the bosonic sector that is in the BKT universality class.

The inset in Fig. 3a demonstrates that the stiffness reduction close to the BKT temperature is due to phase fluctuations and is not because of the fermions. It shows $\Upsilon$ calculated from Eq. (7), but under the assumption that the phases are fully ordered and frozen and temperature affects the system only through the Fermi-Dirac distribution function [see Eq. (9)]. One can see that in this case the stiffness remains finite to temperatures more than an order of magnitude higher than when phase fluctuations are allowed.

By finding the critical temperature for different values of $\tilde{g}$ we construct the phase diagram shown in Fig. 4a, where the solid red and dashed blue lines represent $T_{c}$ and the position of the specific heat maximum, respectively. For weak coupling $T_{c}$ increases with $\tilde{g}$, reaches a maximum at $\tilde{g}\approx 2$ and decreases in the strong-coupling regime.

With increasing temperature, the rapid increase of the number of vortices is accompanied by a maximum of the specific heat that becomes more and more pronounced as the system size increases. The maximum, however, does not necessarily indicate a phase transition. It is known that in the XY model the position of the specific heat maximum is about 10% above the actual critical temperature. Fig. 4a shows that for the PF model the separation in temperature between the specific heat maximum and the BKT transition is much larger.

III fermion sector

It has been shown in the previous section that upon tracing out the fermionic degrees of freedom, the PF Hamiltonian leads to a model with effective interactions between phases on different sites. At sufficiently low temperatures the system develops long range phase coherence in the bosonic subsystem. As the temperature is increased the phases fluctuate and undergo a BKT transition at $T_{c}$. In this section, we analyze how the development of phase coherence in the bosonic sector affects the spectral properties of fermions.

Spectral function: For a given configuration of the bosonic phases $\vec{\theta}$ the real-space fermionic Green function is given by

[TABLE]

where ${\cal H}(\vec{\theta})$ is given by Eq. (4). The spectral function is defined as

[TABLE]

where

[TABLE]

$\langle\ldots\rangle$ denotes averaging over configurations $\vec{\theta}$ generated in the MC sampling. By this procedure the classical bosonic degrees of freedom are integrated out and we can analyze the fermion spectral functions and density of states (Figs. 5) renormalized by their coupling to the phases.

The case of $\tilde{g}=1$ ($\tilde{g}=4$) represents a situation where the interaction is weaker (stronger) than the interaction leading to maximum $T_{c}$ in Fig. 4a at intermediate coupling.

We find that for weak coupling, $\tilde{g}=1$, the system behaves like a BCS superconductor as expected with a gap that closes at $T_{c}$. However, for strong coupling, $\tilde{g}=4$, the spectral function at low temperatures is BCS–like with a minimum gap at the Fermi energy. However upon increasing the temperature there are two interesting observations: (a) the gap does not close even for arbitrarily high temperatures; and (b) the spectral function looks BEC–like with a minimum gap moving to the $\Gamma$ point. A similar change in the spectral function has been reported in the attractive Hubbard model on two coupled triangular lattices Loh et al. (2016) at $T=0$. We report here the first evidence for a temperature–driven phase transition from a BCS superconductor to a disordered bosonic insulator. The mechanism responsible for this transition can be illustrated with the PF model on only two sites, and provides useful insight; see Appendix A.

Essentially, temperature affects the fermionic spectrum primarily via the bosonic phases. At low temperatures, fermions interact with ordered phases and the spectral lines are relatively narrow. At high temperatures, the fermions are scattered by the disordered phases that changes their momenta and energies and leads to broadening of the spectral lines. Note here that the broadening due to the Fermi–Dirac distribution function is negligible and does not affect their shape.

Pairing amplitude: The difference between the weak and strong coupling regimes is seen also in the magnitude of induced pairing between the fermions. Fig. 6a shows the dependence of the pairing amplitude

[TABLE]

on temperature.

In weak coupling $\Delta$ decreases with increasing temperature, whereas in strong coupling regime it remains almost temperature independent, slightly increasing for $\tilde{g}$ around $4$. In both cases at some characteristic temperature $T_{\Delta}$ (indicated in Fig. 6a), the pairing amplitude $\Delta$ reaches a constant value and does not change with further increase of temperature. For strong coupling $T_{\Delta}$ is difficult to determine because of its weak temperature dependence. In weak coupling $T_{\Delta}$ tracks the specific heat maximum, as shown in Fig. 4a. Additionally, in this limit fluctuations of $\Delta$ are strongly pronounced at $T_{\Delta}$, as shown in Fig. 4b.

We attribute the difference in the temperature dependence of $\Delta$ to the balance of different contributions to the total energy: in weak coupling the fermionic pairs are weakly bound and with the increase of temperature, it is energetically favorable to break up the pairs and gain the kinetic energy for the unpaired fermions. As a result, $\Delta$ decreases with increasing temperature.

On the other hand, for strongly bound pairs it is preferable to stay in the paired state even at the cost of localization of the pair wave functions. This makes the pairing amplitude almost temperature independent. The system nevertheless undergoes a temperature–driven superconductor–insulator transition due to the vanishing of the superfluid stiffness.

Analogous to the spectral functions, the temperature dependence of $\Delta$ is a result of the phase disorder and not merely the consequences of broadening of the Fermi-Dirac distribution function.

Mobility of fermions: The properties of fermions are connected to their mobility. In order to find the regime where the single fermion wave functions are localized we calculate the inverse participation ratio (IPR), a quantity that in the thermodynamic limit tends to zero for delocalized states and remains finite for localized ones, defined by:

[TABLE]

where $\Psi_{ni}$ is amplitude of the $n$-th state at site $i$. In order to perform finite size scaling we simulate large systems up to $48\times 48$ using the Travelling Cluster Approximation Kumar and Majumdar (2006); Mukherjee et al. (2015).

Our results are presented in Fig. 6c. Above the BKT transition, fermions are delocalized only for very weak coupling. The transition between localized and delocalized states is almost independent of temperature. Our calculations of the density of states at the Fermi level shows that the fermions get localized while the energy gap is still closed. We attribute the localization in this regime to interactions with disordered phases. There exists a unitary transformation of the fermionic operators that transfers the bosonic phases to the hopping integral $t\to t\exp\left[i(\theta_{i}-\theta_{j})/2\right]$ which shows that the localization of fermions can be understood also as the localization due to a random magnetic field Furusaki (1999).

With increasing coupling a gap opens in the spectrum. For $\tilde{g}\leq 1$ the gap opens at the BKT temperature; for larger $\tilde{g}$, however, the temperature at which the gap opens increases, whereas the $T_{\rm BKT}$ remains constant up to $\tilde{g}=2$ and then decreases. The overall phase diagram is presented in Fig. 6d.

IV Summary

We have proposed an effective boson–fermion model where the interaction between the classical and quantum degrees of freedom leads to highly nontrivial properties. In particular, we have shown that fermions mediate an effective interaction between bosonic phases that is long–ranged and temperature–dependent. This interaction is significantly different from that in the 2D XY model, but nevertheless also leads to a phase transition in the BKT universality class. The bosonic phase configurations fluctuate with temperature and affect the properties of fermions. At low temperatures the phases are almost aligned uniformly and their interaction with fermions drives a BCS–type pairing in the fermionic subsystem. However, when temperature increases, the bosnic phases become disordered and leads to localization of the fermions. When the coupling of the fermions with the bosonic phases is strong, the fermions are still paired, but in this regime they form a BEC state with different spectral properties: the minimum gap locus shifts from the Fermi surface to the $\Gamma$ point. we thus find a temperature–driven BCS–BEC transition.

We have so far assumed in the PF Hamiltonian that the bosons are localized. By allowing them to hop between lattice sites one can introduce a direct interaction that competes or cooperates with the interaction mediated by the fermions.

Acknowledgements.

M.M.M. acknowledges support by National Science Centre (Poland) under Grant No. DEC-2013/11/B/ST3/0082 and N.T. acknowledges DOE-Basic Energy Sciences grant DE-FG02-07ER46423.

Appendix A Two–site model

The PF model can be analytically solved for a system composed of two sites.

The Hamiltonian is given by

[TABLE]

with eigenenergies

[TABLE]

where $\Delta\theta\equiv\theta_{1}-\theta_{2}$. The partition function

[TABLE]

where the effective temperature–dependent interaction between phases $\theta_{1}$ and $\theta_{2}$ is given by

[TABLE]

is equal to the free energy of the fermionic subsystem. In the weak coupling regime the potential has very different dependence on $\Delta\theta$ at low and high temperature (Fig. 8a). The limits of the angle-dependent part of $V$ are given by Eqs. (5) and (6). For stronger coupling the difference is much smaller and for $\tilde{g}=4$ (Fig. 8b) the angle dependence is very well described by $\cos(\Delta\theta)$ at arbitrary temperatures. Note, that in many cases $V(\Delta\theta)$ has the same form as in the XY model.

It can be seen in Fig. 8 that the magnitude of the interaction strongly decreases when temperature increases.

Since in the 2D PF model fermions travel through the entire system, the phase–phase interaction that is mediated by them has long–range character. The knowledge of the interaction in a two-site model allows one to determine the importance of the long–range part of this effective interaction. To this end we have performed MC simulations for a classical model in which phases on a 2D square lattice were interacting only with their nearest neighbors with the interaction potential given by Eq. (23). In this way the microscopic mechanism that leads to the effective phase–phase interaction is the same as in the original PF model, but the range of the interaction is limited to nearest neighbors only. Since for a wide range of the model parameters the interaction in the two–site model is well described by $\cos(\Delta\theta)$, one can expect that the BKT transition should be observed also in the classical model with only nearest–neighbor interaction. MC simulations confirms this assumptions, but it turns out that the critical temperature is higher than in the original PF model with long–range interaction. The convergence with increasing system size is similar to that presented in Fig. 3 and the difference in the critical temperature is, for example, around 20% for $\tilde{g}=4$.

A.1 BCS–BEC Crossover

While there is no collective behavior in the two–site model, it still can be useful to get an insight into the mechanism of the temperature–driven BCS–BEC crossover.

At low temperature $\Delta\theta$ is small and there are two doubly-degenerate BCS–like eigenvalues $\pm\sqrt{\varepsilon_{k}^{2}+\tilde{g}^{2}}$ (the Brillouin zone includes only points 0 and $\pi$ and $\varepsilon_{0}=t$ and $\varepsilon_{\pi}=-t$). However, at high temperature $\Delta\theta$ has a uniform distribution between 0 and $2\pi$. It can be easily shown that then the distribution of the eigenvalues is strongly peaked at $\pm|\varepsilon_{k}\pm\tilde{g}|$. The transition is illustrated in Fig. 9.

Appendix B Vortex-antivortex interaction

The BKT transition results from a logarithmic dependence of the vortex-antivortex interaction on the distance between them. The same behavior can be demonstrated for the PH model. Namely, we fix phases at corners of an elementary plaquette in such way that they form a vortex. Another fixed phases form an antivortex at some distance $d$. Then we start MC simulations for all other phases and after the system gets thermalized we measure its energy $E$ with respect to the energy of a similar system without any fixed phases $E_{0}$. The simulations were performed at low temperature applying the simulated annealing method, so the final configuration is close to the ground state. Fig. 10 shows $E-E_{0}$ as a function of $d$.

For small $d$ the dependence is perfectly described by the logarithmic function. For larger $d$, when is becomes comparable to the linear size of the system, the phases are affected by the presence of the system boundaries (we used periodic boundary conditions) and the energy deviates.

Appendix C Methods

Eqs. (3) and (4) allow us to perform the classical MC simulation for the PF model: in each MC step a phase in a randomly chosen lattice site is changed, the Hamiltonian matrix is diagonalized and the new configuration is accepted or rejected according to the Metropolis criterion. The difference between MC simulations for, e.g., the Ising model and the present one is that here we use the free energy of the fermionic subsystem (4) instead of the internal energy in the acceptance criterion. It is also much more time-consuming, because $\theta$’s are continuous variables, not just $\pm 1$ like in the Ising model, and for a 2D $L\times L$ system in each MC step $2L\times 2L$ matrix has to be diagonalized. On the other hand, the size of the Hilbert space is still much smaller than in the case of systems with full quantum correlations. Additionally, it is possible to significantly speed up the calculations using, e.g., Chebyshev expansion of the fermion density of states Motome and Furukawa (1999), an algorithm that directly updates the spectrum of a successive Hamiltonian matrix based on the spectrum of previous Hamiltonian matrix Alvarez et al. (2007), or Restricted Boltzmann Machines Huang and Wang (2017). Some of the present calculations were performed with the help of the Traveling Cluster Approximation (TCA) Kumar and Majumdar (2006); Mukherjee et al. (2015).