Superfluid-insulator transition and the BEC-BCS crossover in the Rashba moat band

TL;DR

This paper investigates the superconducting transition in a 2D electron gas with Rashba spin-orbit coupling, revealing a tunable BEC-BCS crossover driven by chemical potential or gate distance.

Contribution

It demonstrates how the BEC-BCS crossover can be controlled in a Rashba system through chemical potential and Coulomb interactions.

Findings

Superconducting transition follows BEC or BCS depending on chemical potential.

Crossover can be tuned by gate distance or chemical potential.

Coulomb repulsion influences the transition regime.

Abstract

We study the superconducting transition in a two-dimensional electron gas with strong Rashba spin-orbit coupling. We assume low electron density, such that only the majority spin band participates in the transition. We show that the superconducting transition follows either the Bose-Einstein condensation (BEC), or the Bardeen-Cooper-Schrieffer (BCS) scenarios, depending on the position of the chemical potential with respect to the bottom of the majority band, and the strength of the Coulomb repulsion between electrons. Hence, the BEC-BCS crossover in this system can be driven either by the change in the chemical potential, or the distance to a gate.

Equations

$$H_0=\kappa p^2+\alpha \hat{\mathbf{z}}\cdot (\mathbf{p}\times \mathbf{\sigma }),$$

$${\epsilon }_{\pm }=\kappa p^2\pm \alpha p=\kappa (p\pm p_0{)}^2-E_R,$$

$${\xi }_{\mathbf{p}}=\kappa (p-p_0{)}^2+\Lambda ,$$

$${\chi }_{+}=\frac{1}{\sqrt{2}}[1,i{e}^{i{\theta }_{\mathbf{p}}}]$$

$$H_{int}^{\prime }=\frac{1}{2V}\sum _{\sigma {\sigma }^{\prime }}\sum _{\mathbf{p}{\mathbf{p}}^{\prime }\mathbf{q}}{U}_{\mathbf{q}}{c}_{\mathbf{p}+\mathbf{q},\sigma }^{†}{c}_{{\mathbf{p}}^{\prime }-\mathbf{q},{\sigma }^{\prime }}^{†}{c}_{{\mathbf{p}}^{\prime },{\sigma }^{\prime }}{c}_{\mathbf{p},\sigma },$$

$$H_{int}=\frac{1}{8V}\sum _{\mathbf{p}{\mathbf{p}}^{\prime }\mathbf{q}}{U}_{\mathrm{eff}}(\mathbf{p},{\mathbf{p}}^{\prime },\mathbf{q})c_{\mathbf{p}+\mathbf{q}}^{†}{c}_{{\mathbf{p}}^{\prime }-\mathbf{q}}^{†}{c}_{{\mathbf{p}}^{\prime }}{c}_{\mathbf{p}},$$

$$U_{\mathrm{eff}}(\mathbf{p},{\mathbf{p}}^{\prime },\mathbf{q})=U_{\mathbf{q}}(1+e^{-i{\theta }_{{\mathbf{p}}^{\prime }}+i{\theta }_{{\mathbf{p}}^{\prime }-\mathbf{q}}})(1+e^{-i{\theta }_{\mathbf{p}}+i{\theta }_{\mathbf{p}+\mathbf{q}}}).$$

$$H_p=-\frac{g}{2V}\sum _{\mathbf{p}{p}^{\prime }}{e}^{-i{\theta }_{\mathbf{p}}+i{\theta }_{{\mathbf{p}}^{\prime }}}{c}_{{\mathbf{p}}^{\prime }}^{†}{c}_{-{\mathbf{p}}^{\prime }}^{†}{c}_{-\mathbf{p}}{c}_{\mathbf{p}}.$$

$$H_c=\frac{1}{8V}\sum _{\mathbf{p}{\mathbf{p}}^{\prime }\mathbf{q}}{U}_{\mathrm{eff}}(\mathbf{p},{\mathbf{p}}^{\prime },\mathbf{q})c_{\mathbf{p}+\mathbf{q}}^{†}{c}_{{\mathbf{p}}^{\prime }-\mathbf{q}}^{†}{c}_{{\mathbf{p}}^{\prime }}{c}_{\mathbf{p}}.$$

$$U_{\mathbf{q}}=\frac{e^2}{2ϵq}(1-e^{-qd}),$$

$$|\Psi ⟩=\sum _{\mathbf{p}}{a}_{\mathbf{p}}{c}_{\mathbf{p}}^{†}{c}_{-\mathbf{p}}^{†}|0⟩.$$

$${\alpha }_{\mathbf{p}}=\frac{1}{2{\xi }_p-E}\frac{g}{V}\sum _{\mathbf{k}}{\alpha }_{\mathbf{k}}\equiv \frac{c}{2{\xi }_p-E},$$

$$1=\frac{g}{V}\sum _{\mathbf{p}}\frac{1}{2{\xi }_p-E},$$

$$E=2\Lambda -\frac{g^2{p}_0^2}{8\kappa }.$$

$$g_{c0}=\frac{4\sqrt{\kappa \Lambda }}{p_0}.$$

$$\zeta =\sqrt{\frac{\kappa }{\Lambda }}.$$

$$(2{\xi }_q-E){\alpha }_{\mathbf{q}}+\frac{1}{V}\sum _{\mathbf{k}}{\tilde{U}}_{\mathbf{q},\mathbf{k}}{\alpha }_{\mathbf{k}}=\frac{g}{V}\sum _{\mathbf{p}}{\alpha }_{\mathbf{p}},$$

$$g_{\mathrm{eff}}=g-\frac{e^{2}d}{4ϵ}.$$

$$g_{cb}=\frac{e^{2}d}{4ϵ}.$$

$$g_c=g_{c0}+\frac{E_{c}d}{2p_0},$$

$$2{\xi }_q{\alpha }_q-g{p}_0{\int }_p{\alpha }_p+{\int }_k\mathcal{K}(q,k){\alpha }_k=E{\alpha }_q,$$

$$\mathcal{K}(q,k)=\frac{E_c}{2\pi }{\int }_0^{2\pi }\frac{{\cos }^2({\theta }_{\mathbf{qk}}/2)(1-e^{-d|\mathbf{q}-\mathbf{k}|})}{|\mathbf{q}-\mathbf{k}|}d\theta ,$$

$${\int }_p\equiv {\int }_{-\infty }^{\infty }\frac{dp}{2\pi }.$$

$$\mathcal{K}(q,k)\to \mathcal{K}\left(p_0+\frac{q-k}{2},p_0-\frac{q-k}{2}\right)\equiv \mathcal{K}(q-k).$$

$$\mathcal{K}(q)\approx \frac{E_c}{2\pi p_0}\ln \frac{4{\beta }^2{d}^2{p}_0^2}{1+\frac{{\beta }^2}{4}{d}^2{q}^2},$$

$$U(x)=\frac{E_c}{2\pi p_0|x|}\exp \left(-\frac{2|x|}{\beta d}\right).$$

$$-2{\alpha }^{\prime \prime }(x)-g\delta (x)\alpha (x)+\frac{E_c}{2\pi |x|}{e}^{-\frac{2|x|}{\beta d}}\alpha (x)=(E-2\Lambda )\alpha (x).$$

$${\alpha }_0(x)=A\sqrt{|x|}{K}_1\left(\sqrt{\frac{E_c|x|}{\pi }}\right),$$

$$K_1(z)\sim \sqrt{\frac{\pi }{2z}}{e}^{-z}.$$

$${\zeta }_0=\frac{\pi }{E_c}\gg 1.$$

Full text

I Introduction

BEC-BCS crossover between the limits of Bose-Einstein condensation of tightly bound two-electron molecules and superconducting transition of widely overlapping pairs lies at the heart of the current understanding of fermionic superfluidity. The subject has received an extensive theoretical and experimental attention in variety of systems including 2D systems with spin-orbit-coupling Leggett (1980); Chen et al. (2012); Ries et al. (2015); Valtolina et al. (2015); He and Huang (2012), culminating in the observation of the crossover phenomenon in dilute alkali gasesAltmeyer et al. (2007); Papp et al. (2008). In this work, we study the BEC-BCS in a two-dimensional (2D) electron gas with low carrier density, realized in a 2D heterostructure with structural inversion asymmetry, and the ensuing Rashba spin-orbit coupling. Systems with similar conditions have been previously reported Gui et al. (2004); Koroteev et al. (2004); Singh and Romero (2017); Ast et al. (2008, 2007); Jo et al. (2017).

There have been a number of theoretical studies of fermionic superfluidity in systems with a single particle spectrum having a minimum on a line in 2D, or a surface in 3D – the “moat” Gor’kov and Rashba (2001); Cappelluti et al. (2007); Goldstein et al. (2015); Vyasanakere and Shenoy (2011); Chaplik and Magarill (2006); Takei et al. (2012). The former option is realized in the lower spin subband of a Rashba 2DEG, and is considered in this work. The main attractive feature of such systems is their effectively one-dimensional (1D) thermodynamics, brought about by the diverging 1D density of states (DoS) near the bottom of the moat band. The 1D DoS of moat bands makes it easier to form two-particle bound states, and to form Bose-Einstein condensates. For instance, in 2D systems, which are the primary interest in the present work, the critical temperature shows a power-law dependence on the BCS coupling constant, as opposed to the exponentially small result for the usual parabolic dispersion in 2DEGs with constant DoS Cappelluti et al. (2007). This gives the proposal a higher chance of experimental observation.

In this paper, we study the competition between BCS attraction and Coulomb repulsion in the phase diagram of a 2DEG with a moat band. The repulsion is necessarily present in an electronic system, and the lack of screening in a 2D system with low carrier density makes Coulomb effects strong. We show that electrostatic repulsion between electrons, even when partially screened by an external gate, is of qualitative importance at the superconducting transition. Strong electrostatic repulsion induces an interesting BCS-BEC crossover in the presence of a shallow Fermi sea in the moat band.

The rest of the paper is organized as follows. In section II we construct the model of a 2DEG with Rashba SOC. Assuming a low density limit next we introduce an attractive and a repulsive interaction term projected on the lower moat-like band. In section III we treat the problem by considering a single pair wavefunction and showing that it can describe many features of superfluid-insulator transition including the effects of Coulomb repulsion. Section IV extends the problem to the case of a shallow ring-like Fermi sea. There we study the BEC-BCS crossover and construct the phase diagram of the system in the absence and presence of Coulomb repulsion. Finally we provide a brief summary of the results in Section V. Our main findings are summarized in phase diagrams in Figures 5 and 7.

II Model of a 2DEG with a moat band

For a system with the Rashba spin-orbit coupling that stems from structural inversion asymmetry, the single particle Hamiltonian in the momentum space has the following form:

[TABLE]

where $p$ is the particle’s momentum, and $\alpha$ is the characteristic velocity that determines the strength of the spin-orbit coupling. The two eigen-bands of $H_{0}$ are given by

[TABLE]

where we defined the Rashba energy, $E_{R}=\alpha^{2}/(4\kappa)$. The corresponding dispersions are plotted in Fig 1.

The lower (“moat”) band, $\varepsilon_{-}$, has a minimum at a circle in momentum space, which is centered at the origin, and has the radius of $p_{0}=\alpha/2\kappa$. The energy minimum is given by $\varepsilon_{-}(p_{0})=-E_{R}$. Therefore, the single-particle energy of the moat-shaped band, measured from the chemical potential, can be written as

[TABLE]

where $\Lambda$ is the minimum energy of the moat band measured from the chemical potential (see Fig. 1).

The single-particle wave functions that correspond to the bands of Hamiltonian (1) are represented by spinors, in which the spin orientation is determined by a particle’s momentum due to the presence of the spin-orbit coupling. The spinors for the $``\pm"$ bands have the following form in the basis of $\sigma_{z}$ eigenstates:

[TABLE]

where $\theta_{\bm{p}}$ is the angle momentum ${\bm{p}}$ makes with the $\hat{x}$-axis, hence $e^{i\theta_{\bm{p}}}=\hat{p}_{x}+i\hat{p}_{y}$.

Now we consider the general form of density-density interaction

[TABLE]

where $U_{\bm{q}}$ is the Fourier transform of the interaction potential, and $c_{{\bm{p}},\sigma}$ is the annihilation operator of a particle with momentum ${\bm{p}}$ and spin $\sigma$. We are interested in the low-density case, in which the chemical potential lies far below $\varepsilon_{+}$ band bottom, such that this band’s influence on the particle dynamics is insignificant. This allows to project the interaction Hamiltonian on the lower ($``-"$) band with the help of Eq. (4):

[TABLE]

with

[TABLE]

In what follows we are interested in the BEC-BCS crossover in the system of Fermions described by Hamiltonians (1) and (6). Specifically, we restrict ourselves to the case where the interaction Hamiltonian, that is the amplitudes $U_{\bm{q}}$, contain attraction in the Cooper channel, and repulsion in the Coulomb channel.

The Cooper channel of the problem is described by a contact attractive interaction of $U_{\bm{q}}\to-2g$ only between electrons with opposite momenta $\bm{p^{\prime}}=-\bm{p}$. We assume that this attractive interaction $g$ stems from the usual electron-phonon interaction, its specific value being determined by the details of the material hosting two-dimensional electron gas. Note that spins of these electrons are also opposite to each other. In the Cooper channel, the phase factor in (7) reduces to $(1+e^{-i\theta_{{\bm{p}}}+i\theta_{\bm{p+q}}})^{2}$; however since $c_{\bm{p}}c_{-\bm{p}}=-c_{-\bm{p}}c_{\bm{p}}$ even in momentum terms vanish and the pairing interaction can be written as

[TABLE]

For the Coulomb channel however, we work with full interaction form presented in (7). The significance of the phase factor, as we show in the rest of the paper, shall depend on the external screening mechanism that we introduce below. Therefore, the repulsive part of the Hamiltonian maintains the general form in (6) and (7).

[TABLE]

The long-range nature of the Coulomb interaction brings up the issue of screening in a system of charged Fermions. In the two-dimensional (2D) system with low-electron density, screening of Coulomb interaction via polarization is not effective. Instead we assume that our 2D layer of Rashba material is screened by a proximate conducting gate, see Fig. 2, which also provides a simple way to control the (electro-)chemical potential of the system by varying the gate voltage. The Coulomb potential, describing instantaneous screening by the gate’s image charges, takes the following form:

[TABLE]

where $d/2$ is separation between the gate and the layer and $\epsilon$ is dielectric constant of the material in which two-dimensional conducting layer is formed (semiconducting heterostructure).

III BEC regime: pair wave function approximation

In this Section, we analyze the ground states of the low-density Fermionic gas described by Hamiltonians (1), (8) and (9). Specifically, we study the evolution of the ground state of the system as the gate voltage is adjusted to control the position of the electrochemical potential with respect to the bottom of the moat band, assuming it remains below the moat band bottom, $\Lambda\geq 0$ in Fig. 1: we refer to this situation as the BEC regime, see below.

For non-interacting particles, $\Lambda\geq 0$ would mean that particles do not enter the two-dimensional layer, which remains insulating. However, in the presence of particle attraction (the Cooper channel), there can be bound states of pairs of particles. For strong enough attraction, their energies can lie below the chemical potential for pairs (twice the chemical potential for single particles). In that case, pre-formed pairs would actually enter the system, and form a Bose-Einstein condensate (BEC) of pairs. We refer to this kind of situation as the BEC limit.

It is clear that since there is a finite binding energy required for a pair to enter the system, the spatial extent of such a pair is also finite, being determined by roughly the inverse binding energy. On the other hand, for the range of parameters where the pairs are just (energetically) able to enter the system, their density is vanishingly small, and one indeed deals with the case of a dilute Bose gas of pairs.

The diluteness of the pair gas implies that in order to find the conditions at which the BEC occurrs, one only needs to determine the conditions required for the first pair to enter the system. Therefore, in this Section we solve the single-pair Schrödinger equation to determine when the pair bound state goes under the chemical potential for pairs. Since the chemical potential plays only an auxiliary role in such considerations, setting the value below which the bound state energy needs to go, one can formally solve the single-pair problem even for $\Lambda<0$. This is especially relevant in the case with Coulomb repulsion, which effectively increases the pair fugacity above $2\Lambda$, hence $\Lambda=0$ is not special within the single-pair problem in that case. However, the physical starting point - a single pair in the system - is not correct in the $\Lambda<0$ case, since one must start with a Fermi sea, hence we do not consider it.

III.1 Neutral Fermions with short-range attraction

We start with the case without the Coulomb interaction, pertaining to neutral Fermions with short-range attraction. In order to model and study this limit we introduce a simple ansatz wave-function made of a superposition of Cooper pairs at different possible momenta, which is an exact eigenstate of the pairing Hamiltonian:

[TABLE]

The form of this wave function is dictated by the form of the pairing Hamiltonian, Eq. (8), and in the absence of the Coulomb repulsion (when $H=H_{0}+H_{p}$) the wave function (11) is an exact eigenstate of the total Hamiltonian. The Schrödinger equation $H|\Psi\rangle=E|\Psi\rangle$ gives us the following form for the bound state wave-function in the momentum-space

[TABLE]

where $\alpha_{\bm{p}}=a_{\bm{p}}e^{-i\theta_{\bm{p}}}$, $\xi_{p}$ is given by Eq. (3), and $c$ is a constant given by the integral over ${\bm{k}}$. Solving Eq. (12) self-consistently leads to a simple integral equation for $E$ as follows:

[TABLE]

Approximating the integral is straightforward and gives us

[TABLE]

A negative value of $E$ implies that it is energetically favorable for a pair to enter the system. Therefore, the condition $E=0$ determines the critical value of the coupling constant, $g$, at which a condensate forms in the system (at a given value of $\Lambda$). It yields the critical interaction strength for the superconductor-insulator transition in the system of neutral fermions, $g_{c0}$

[TABLE]

This results coincides with that of Ref. Goldstein et al., 2015, which was obtained using field-theoretical methods. Here, we just emphasize the very simple nature of this problem, which reduces to calculating a pair bound state.

Using the treatment presented above, one can also easily confirm that at the critical value of $g$, when the boundstate goes under the chemical potential for pairs, the corresponding wave-function has finite spatial extension – the pair size. It is easy to check that this characteristic length scale is given by

[TABLE]

III.2 Effects of Coulomb Repulsion

Now if we include the Coulomb repulsion, $H_{c}$, in the Hamiltonian, the Schrödinger equation becomes

[TABLE]

where $\tilde{U}_{{\bm{q},\bm{k}}}=U_{{\bm{q}-\bm{k}}}(1+\cos(\theta_{\bm{q}}-\theta_{\bm{k}}))/2$.

The key difference brought about by the Coulomb repulsion is the fact that now the bound state formation requires a finite strength of the attraction. In the case of pure attraction, the bound state exists at any positive value of $g$, see (13), while a finite critical value $g_{c0}$ (15) follows from the condition that the bound state must be $2\Lambda$ deep in order for pairs to enter the system. Below, we will refer to the critical value of $g$ needed for the bound state formation as $g_{cb}$, keeping the notation $g_{c}$ for the value of $g$ that marks the onset of superconductivity.

The form of solutions to this equation depends on the strength of screening. One should distinguish the limits of strong, $dp_{0}\ll 1$, and weak, $dp_{0}\gg 1$, screening, which we discuss in detail below.

III.2.1 Strong screening

When the screening is strong, i.e. when $dp_{0}\ll 1$, $U_{q}$ (defined in (10)) can be approximated by the constant $e^{2}d/2\epsilon$, reflecting the fact that it becomes a short-range contact-like interaction. Therefore, in that case the Coulomb term in (17) can be absorbed in the constant attraction term on the right hand side. Absorbing the Coulomb term, one should note that due to the presence of the angular factor, $(1+\cos(\theta_{\bm{q}}-\theta_{\bm{k}}))/2$, it would diminish by a factor of $1/2$ over the angular integral. With this consideration we obtain a renormalized effective attraction

[TABLE]

It is obvious that the value of $g$ required for the formation of a bound state is

[TABLE]

In turn, just like in the case of pure attraction, the onset of superconductivity is marked by the value of $g$ required for the bound state to cross the value of the chemical potential for pairs, $2\Lambda$. As a result, strong screening of the Coulomb interaction by the conducting gate results in a weak enhancement of the critical strength of attraction $g_{c}$,

[TABLE]

where $g_{c0}$ is given by Eq. (15), and we defined $E_{c}=e^{2}p_{0}/2\epsilon$ as the characteristic Coulomb energy scale in the moat band.

In Fig. 3 the analytic approximation (20) is compared with the results found by solving the Schrödinger equation (17) numerically. One can see that the strong-screening results - the linear in $d$ enhancement of the critical attraction strength - persists to $dp_{0}\sim 0.1$.

III.2.2 Weak screening

In the remainder of this Section, we focus on determining the condition for the appearance of the bound state in the two-particle spectrum in the limit of weak screening, $dp_{0}\gg 1$, for which the Coulomb repulsion of Eq. (10) is strong and $q$-dependent. We will also focus on describing qualitative features of this bound state. Full numerical solution of Eq. (17) is deferred till Section IV, where it will be used to determine the condensation line on the BEC side of the phase diagram in the $(g,\Lambda)$ plane.

While the Schrödinger equation (17) can indeed be solved numerically, considerable insight into the the problem can be gained by transforming it into an effective Schrödinger equation in auxiliary one-dimensional “real” space, which is reciprocal to the radial direction in the momentum space.

The transformation proceeds by reducing Eq. (17) to a radial equation (for a rotationally-symmetric $\alpha_{\bm{q}}=\alpha_{q}$):

[TABLE]

where the Coulomb kernel ${\cal K}(q,k)$ is defined as

[TABLE]

$\theta_{\bm{qk}}$ being the angle between ${\bm{q}}$ and ${\bm{k}}$, and we introduced a shorthand notation

[TABLE]

Furthermore, as explained above, we are interested in finding two-electron bound states with vanishing energy, $E\to 0$, for the purpose of finding the boundary of the normal-superconducting transition. For $\Lambda\ll E_{R}$, only momenta in the close vicinity of $p_{0}$ participate in the formation of such shallow bound states, which is enforced by the term with $\xi_{q}$ in Eq. (21). Combined with the fact that the Coulomb kernel $\cal{K}$ is sharply peaked around $q=k$, this allows us to count all momenta from $p_{0}$ and extend all integrations over the magnitude of such reduced momenta to infinite limits, as well as substitute

[TABLE]

This implies that Eq. (21) maps onto a 1D Schrödinger equation, in which $2\xi_{q}=2\kappa(q-p_{0})^{2}+2\Lambda\to-2\kappa\partial_{x}^{2}+2\Lambda$ plays the role of an effective 1D kinetic energy (shifted by $2\Lambda$, as we count the energeis from the chemical potential); $-gp_{0}\delta(x)$ represents the short-range attractive potential; ${\cal K}(x)$ - the Fourier transform of ${\cal K}(q)$, is a long-range repulsive potential.

As demonstrated in Appendix A, for relatively small $q\ll p_{0}$, which participate in the pairing problem at realistic values of parameters, one can approximate ${\cal K}(q)$ with

[TABLE]

where $\beta=e^{\gamma}\approx 1.78$, with $\gamma\approx 0.58$ being the Euler-Mascheroni constant. Expression (25) is of interpolation nature, and connects the $q\to 0$ limit of ${\cal K}(q)$ with the $q\gg 1/d$ limit.

The “real”-space potential that corresponds to the kernel (24) is derived in Appendix A, and is given by Eq. (68), which we repeat here for convenience:

[TABLE]

This potential defines the behavior of the effective 1D quantum-mechanical problem down to distances $|x|\gtrsim 1/p_{0}$. At smaller scales, the large-momentum behavior of ${\cal{K}}(q,p)$ becomes important, and Eq. (24) ceases to hold.

From now and until the end of this Section, we will measure all energies in units of $E_{R}$, all momenta in units of $p_{0}$, all coordinates in units of $1/p_{0}$, and $g$ in units of $\kappa$. Then the effective Schrödinger equation in real space is written as

[TABLE]

To qualitatively understand the solutions of Eq. (27), we start with the $d\to\infty$ limit when $e^{-\frac{2|x|}{\beta d}}\to 1$. First, we focus on finding the critical value of $g$, $g_{cb}$, at which the first bound state emerges below the continuum – this signals the appearance of a two-particle bound state in the original problem. To find $g_{cb}$, we set $E=2\Lambda$, and retain only the solution that does not grow indefinitely for $x\to\infty$. At distances $|x|\gtrsim 1$, such solution is

[TABLE]

where $A$ is an arbitrary constant, and $K$ is the modified Bessel function of the second kind. Note that the asymptotic behavior of this function at large values of the argument is given by

[TABLE]

This implies that even though the state $\alpha_{0}$ has energy infinitesimally below the continuum, it has a finite spatial extent $\alpha(x)\sim\exp(-\sqrt{|x|/\zeta_{0}})$, determined by

[TABLE]

(Note that in dimensionful units $\zeta_{0}=\pi E_{R}/(E_{c}p_{0})$.) The critical value of $g$ is found from the condition that $\alpha_{0}(x)$ has a jump in derivative at the origin determined by $g$ – the strength of the attractive $\delta$-function – in the standard way. We note that for $x\to 0$, the solution presented in Eq. (28) has a logarithmically divergent derivative. This divergence is unphysical, since the form of the Coulomb potential is not valid near the origin. Therefore, we cut-off the divergence at $x\sim 1$, to obtain

[TABLE]

It is clear that the critical $g_{cb}$ given by Eq. (31) also persists for finite $d\gtrsim\zeta_{0}$.

In turn, for $1\ll d\lesssim\zeta_{0}$, the repulsive potential behaves qualitatively like a $\delta$-functional one, simply renormalizing the strength of the attractive potential. That is,

[TABLE]

Performing the integral, and taking the $d\gg 1$ asymptotics of the result, we observe that in this limit the repulsive interaction renormalizes $g$ by

[TABLE]

The bound state appears when the renormalized strength of the $\delta$-function potential is positive, hence in this regime

[TABLE]

Note that for $d\sim\zeta_{0}$, Eqs. (31) and (34) are parametrically the same.

Taken together, Eqs. (19), (31), and  (34) determine the dependence of the critical attraction strength required to bind a two-particle bound state for all values of $d$. This dependence is schematically illustrated in Fig. 4.

The effective Schrödinger equation (27) can also help us determine the dependence of the bound state energy, $E_{b}$, on $g>g_{cb}$. To this end, we write this equation for $\alpha_{0}$, as well as for $\alpha_{g}$, the latter being assumed to solve the equation for $g$ above the critical value, and take into account that $2\Lambda-E_{b}>0$ :

[TABLE]

We now multiply the first equation with $\alpha_{g}$, and the second one with $\alpha_{0}$, subtract one from the other, and integrate over the entire axis. We immediately get

[TABLE]

The value of the prefactor in front of $g-g_{cb}$ depends on $d$ and the value of $g$ as compared to $g_{cb}$. In two cases, either for $d\lesssim\zeta_{0}$ and any $g$, or for $g\gg g_{cb}$, the size of the bound state in the effective $\delta$-potential is determined by its energy, $\zeta_{g}\propto 1/\sqrt{\varepsilon_{cb}(g)}$. Thus we have

[TABLE]

and hence

[TABLE]

.

In the limit of large $d\gtrsim\zeta_{0}$, and for $g-g_{cb}\ll g_{cb}$, we can estimate

[TABLE]

and conclude that

[TABLE]

In principle, the boundary of the insulator-superconductor transition can be obtained by solving Eq. (36), and setting $E_{b}(g)=0$. We will not pursue the task of its determination, deferring it until Section IV.2, where many-body aspects of the condensation problem are discussed. The intuition used in that Section, however, heavily rests on the results of this one.

IV Variational BCS wave-function

We now improve our treatment by using a BCS-type variational wave function. We restrict ourselves to the study of the many-body ground state of the system at $T=0$, and aim to investigate the situation with small but finite electron density in the moat band. Later on, we use the zero-temperature results to estimate the critical temperature of the system.

It should be clarified here that by considering a BCS variational wave function, we are assuming that the normal state of the model at $\Lambda<0$, i.e. when the chemical potential is inside the moat band, is a Fermi liquid one. This assumption is based on previous studies Berg et al. (2012); Ruhman and Berg (2014) of the phase diagram of the repulsively-interacting electrons which found that other, more exotic, electron phases in this interesting regime - such as Wigner crystal, nematic and ferromagnetic phases - require significant Coulomb repulsion (see, for example, Fig.1 of Ref.Ruhman and Berg, 2014). Such a strong repulsion is not consistent with the superconducting state we are considering here.

As shown below, the BCS wave function allows us to study the behavior of the superconducting order parameter $\Delta$ across the crossover region between the BEC and BCS limits, in the absence and presence of Coulomb interaction.

The variational wave function that we consider is

[TABLE]

where $u_{\bm{p}}$ and $v_{\bm{p}}$ are the variational functions; $v_{\bm{p}}$ is multiplied by the phase factor which is chosen to compensate for the phase factors in the attractive channel, Eq. (8). Similar to Section III we first study the model without repulsion and then will investigate the effects of Coulomb repulsion.

IV.1 Phase diagram for neutral fermions

We start with summarizing the phase diagram in $(g,\Lambda)$ coordinates for the case of neutral fermions with attraction. It is depicted in Fig. 5.

Qualitatively, for $g>0$, the phase diagram contains two phases: an insulator, and a superconductor. The superconducting phase is stable for all values of $\Lambda\leq 0$ (chemical potential is inside, or touches the majority band), and for $g>g_{c}(\Lambda)$ for positive $\Lambda$ (when the chemical potential is below the majority band), see Fig. 5. For $\Lambda>0$, the transition from the insulator into a superconductor upon increasing $g$ occurs into a BEC-type superconducting phase, which is represented by a dilute gas of bound pairs of electrons. For $\Lambda\leq 0$, the superconducting phase has BCS character, in which it is no longer permissible to disregard the fermionic nature of the electrons comprising the pairs in the condensate, and think of it as a dilute Bose gas.

There are two notable lines in the phase diagram of Fig. 5: the line of superconductor-insulator transition, and the line of BEC-BCS crossover. The superconductor-insulator transition line is determined from the condition that the bound state energy for electronic pairs lies at the chemical potential for pairs (chosen to be zero in our analysis). In other words, counted from the bottom of the two-electron continuum, the binding energy of the pair must be equal to $2\Lambda$. The corresponding critical value of $g$ was calculated in Section III.1, and is given by Eq. (15).

The BEC-BCS crossover line has the inherent vagueness in its definition, which depends on the criterion for the crossover. Here we adopt the following logic for drawing such lines: the BEC-BCS crossover is a crossover between regimes in which the Fermi statistics of pair constituents either does not, or does matter. The importance of carrier statistics is dictated by occupations of states in momentum space. So we define the crossover as the situation in which the maximum occupation of states in the momentum space is 1/2. It will become apparent below that for the case of neutral fermions, and within the model adopted in the present paper, the BEC-BCS line is then the vertical $\Lambda=0$ one - the $g$-axis.

To get further quantitative insight into the behavior of the superconducting order parameter, we study its variation along a horizontal cut of the phase diagram in Fig. 5, that is, we vary $\Lambda$ at a given value of $g$. We will show below that the superconducting gap, $\Delta(\Lambda)$, is represented by a universal curve, which we determine below.

Minimization of (1) and (8) over the BCS wave function (41) leads to the standard self-consistent equation Jishi (2013) for superconducting order parameter $\Delta$

[TABLE]

where $\xi_{\bm{p}}$ is the single electron dispersion (3). Integrating over the angle and noting that the most important part of the radial integral comes from the vicinity of $p_{0}$, we approximate $pdp\sim p_{0}dp$ and extend the $p$-integral to $\pm\infty$. This gives

[TABLE]

where $\tilde{\Lambda}=\Lambda/\Lambda_{c0}$ and $\tilde{\Delta}=\Delta/\Lambda_{c0}$ are made dimensionless with the help of $\Lambda_{c0}=g^{2}p_{0}^{4}/16E_{R}$. This integral equation holds for both positive (chemical potential below the band) and negative (chemical potential inside the moat band) $\tilde{\Lambda}$. Note that the characteristic energy $\Lambda_{c0}$ scales up like $g^{2}$. It is clear that for reasonably small $g$ $\Lambda_{c0}$ stays much smaller than the Rashba energy $E_{R}$. Therefore $\tilde{\Lambda}$ can be large even in the ‘deep moat’ limit $|\Lambda|\ll E_{R}$.

Fig. 6 shows numerical solution $\tilde{\Delta}(\tilde{\Lambda})$ of Eq. (43). Finite $\tilde{\Delta}$ for positive $\tilde{\Lambda}$ describes BEC regime - condensed state of two-electron pairs. The order parameter vanishes at $\Lambda=\Lambda_{c0}$ with a discontinuity in its derivative. It can be shown that asymptotic behavior of $\tilde{\Delta}$ near the transition point is $\tilde{\Delta}\sim\sqrt{\frac{8}{3}(1-\tilde{\Lambda})}$.

The other feature of function $\tilde{\Delta}(\tilde{\Lambda})$ (Fig. 6) is that it has a maximum on the negative side, when the chemical potential is inside the moat band. That implies that the highest critical temperature happens when the chemical potential is slightly above the band edge.

Moreover, on the BCS side, when electrons fill up a shallow ring-shaped Fermi sea, the order parameter, $\tilde{\Delta}$, has two distinct behaviors: when $|\Lambda|$ is of order of $\Lambda_{c0}$, the order parameter is of the same order and is proportional to $g^{2}$. For larger $|\Lambda|$, when $|\tilde{\Lambda}|\gg 1$, the order parameter decays exponentially as $\Delta\propto|\Lambda|\exp\left(-\frac{\pi}{2}\sqrt{|\Lambda|/\Lambda_{c0}}\right)$. To see this, expand $(x^{2}-|\Lambda|)^{2}\approx 4|\Lambda|(x-|\Lambda|^{1/2})^{2}$ and approximate the integral in (43) by the region around $x\approx|\Lambda|^{1/2}$. The large-$|\Lambda|$ has a familiar form $\exp(-{\text{const}}/g)$ for the order parameter in a standard BCS theory suggesting a superconducting state with strongly overlapping electron pairs at large enough $|\Lambda|$.

Finally, the BEC-BCS crossover line can be deduced from the usual BCS coherence factor, $v_{p}$, whose square defines the occupation of states in the momentum space:

[TABLE]

It is clear that for $\Lambda>0$, and hence $\xi_{p}>0$, all occupation numbers are smaller than 1/2. In turn, for $\Lambda<0$, when $\xi_{p}$ changes sign at the Fermi momenta, there is a region in momentum space – where there is “water” (i.e. electron liquid) in the moat – where the occupation numbers exceed 1/2. For $\Lambda=0$, regardless of $g$, the maximum occupation is exactly 1/2, and is reached at $p=p_{0}$, at the bottom of the moat. Hence the $\Lambda=0$ is the BEC-BCS crossover line according to the criterion we chose, as shown in Figure 5.

IV.2 Phase diagram with Coulomb repulsion

Inclusion of Coulomb repulsion between electrons changes the phase diagram, which is now shown in Fig. 7. It is convenient to describe phases appearing in this diagram starting with the behavior of the system along $\Lambda=0$ line. As we know from Section III.2, in the presence of the Coulomb repulsion, there is a finite value of $g$, which we denoted with $g_{cb}$, needed to achieve binding of pairs even at $\Lambda=0$. Therefore, for $0\leq g<g_{cb}$, the $\Lambda=0$ line is the boundary between the insulating phase on the right, and the normal metal phase on the left. In other words, the system behaves as a semi-metal with diverging one-dimensional density of states for $g<g_{cb}$, and $\Lambda=0$. To the left of this semimetallic phase, there lies a normal metal phase. The doping level of this metal is lower than the nominal one prescribed by $|\Lambda|$, due to electrostatic effects.

Superconducting phases can be reached from both metallic and insulating phases. As in the case of neutral fermions, one reaches BEC- or BCS-like states increasing $g$ on the $\Lambda\gtrless 0$, respectively. Exactly at $\Lambda=0$, one goes into a BEC superconductor, since the BEC-BCS line lies on the left of the $g$-axis in this case, see below.

Having described the qualitative features of the phase diagram for charged fermions, we elaborate now on its quantitative aspects. The Coulomb repulsion brings two new terms to the total energy, one through the Cooper channel (third term in the equation below) and the Hartree-Fock term (forth term below), to which the attractive interaction gives no contribution,

[TABLE]

where $\tilde{U}_{{\bm{k,k^{\prime}}}}=U_{{\bm{k-k^{\prime}}}}(1+\cos(\theta_{\bm{k}}-\theta_{\bm{k^{\prime}}}))/2$, as before. Minimizing the energy over $v_{\bm{k}}=\sin\alpha_{\bm{k}}$ (with $u_{\bm{k}}=\cos\alpha_{\bm{k}}$) and introducing fields

[TABLE]

we obtain desired coupled self-consistent equations for the superconducting and Hartree-Fock (HF) fields

[TABLE]

It is clear from Eqs. (47) that the Coulomb repulsion has two-fold effect on the pairing problem. First, it changes the coupling constant in the equation for the order parameter $\Delta_{\bm{p}}$, $g\to g-\tilde{U}_{\bm{p,p^{\prime}}}$, making it a coupling matrix. Physically, this describes interaction of electrons within a pair, and the effects of this interaction have been described in Section III.2.

Second, the Coulomb repulsion leads to the appearance of the Hartree-Fock renormalization of the single-particle spectrum, described by $\Phi_{\bm{p}}$. The first bracket in the equation for $\Phi_{\bm{p}}$ is just $2v_{\bm{p^{\prime}}}^{2}$, twice the particle number at the momentum $\bm{p}^{\prime}$. Moreover, the second bracket in the same equation is always positive because Coulomb interaction (10) is maximum at $\bm{p}=0$, $\tilde{U}_{0}=e^{2}d/2\epsilon$. Therefore the new collective field $\Phi_{\bm{k}}$ is a kind of ${\bm{k}}$-dependent chemical potential which provides a positive shift of the single particle spectrum $\xi_{\bm{k}}$. This is the same as downward renormalization of the chemical potential and has the effect of ‘pushing out’ electrons from the system. The renormalization is proportional to the Coulomb repulsion experienced by a particle at momentum ${\bm{k}}$ due to all other particles (at momenta ${\bm{k}^{\prime}}$).

In a usual good conductor such a renormalization is neglected due to the strong screening of Coulomb interaction. This point is made explicit in classic papers on strongly coupled superconductivity, see for example Scalapino et al. (1966); Schrieffer (1999) (note that $\Phi$ is denoted as $\chi$ there). The reason is that renormalization of the spectrum, which represents a “simple scale change” (Ref.Scalapino et al., 1966) of electron dispersion, is a small effect on the scale of Fermi energy $E_{F}$ which is present in both normal and superconducting states. In our case, the normal state is either a band insulator, or a metal with a very small carrier density, hence the scale change due to $\Phi$ is crucial in determining whether or not electrons can populate the bottom of the moat band, either in the form of interacting electron gas or bound two-electron molecules.

Therefore the combination of low-density and weak screening (which is determined by the distance due to the external gate) results in the HF field $\Phi_{\bm{p}}$ playing an unusually important role. The basic physics is that of an electrostatic repulsion between two electrons attempting to form a bound state. In the absence of screening and retardation (our attractive potential $g$ is frequency independent constant) the electrostatics of the pair becomes crucial. As discussion in Sec.IV.2.2 below shows, weakly screened Coulomb interaction strongly renormalizes bare chemical potential (see Eq.(61) below) and squeezes electrons out, providing a kind of ‘Coulomb blockade’ phenomenon in a macroscopic setting. This crucial physics is missing in Ref.Goldstein et al., 2015.

It is clear from Eq. (47) that the presence of the Coulomb interaction leads to $\bm{p}$-dependence of both $\Delta_{\bm{p}}$ and $\Phi_{\bm{p}}$. However, it is not so clear if the possible solutions are radially symmetric or not. We have carefully investigated this question numerically and concluded that the lowest energy solution of the coupled equations (47) is s-wave symmetric. We would like to add here that frequently evoked Kohn-Luttinger physics Vafek and Wang (2011), which generally promotes pairing in the finite angular momentum channel and is mediated by repulsive inter-electron interactions via Friedel-like polarization of the fermionic background Kagan et al. (2015), is not expected to contribute in the present situation of a very dilute two-dimensional electron gas. In this case, one needs to include diagrams up to the third-order in repulsive interaction Chubukov (1993). This leads to an exponentially small in the density $n$ dependence of the resulting critical temperature, $T_{c}\sim\exp[-1/n^{3}]$, which is not competitive with the estimate of the critical temperature based on our consideration, see discussion at the end of this Section.

Therefore below we investigate solution with s-wave symmetry and perform angular integration to obtain one dimensional integral equations as follows:

[TABLE]

where the Coulomb kernel ${\cal K}(p,p^{\prime})$ is given by Eq. (22). The upper momentum cut-off, $p_{\rm max}$, is taken to be $2p_{0}$ in numerical solutions of the above equations.

Armed with Eqs. (LABEL:eq:self-consis-colmb-1d), we are going to consider the limits of strong ($dp_{0}\ll 1$) and weak-screening ($dp_{0}\gg 1$), similar to Section III.

IV.2.1 strong screening: $dp_{0}\ll 1$

We start with the limit of strongly screened Coulomb interaction. In this limit

[TABLE]

With this approximation, the $p$-dependence of $\Delta$ and $\Phi$ disappears, and the effect of Coulomb interaction is to reduce $g$ by $E_{c}d/2p_{0}$ (compare with (20)) resulting in an effective attraction

[TABLE]

Following the notation used in Section  IV.1, we define $\Lambda_{c,\text{eff}}=g_{\text{eff}}^{2}p_{0}^{4}/16E_{R}$, introduce dimensionless combinations $\tilde{\Lambda}=\Lambda/\Lambda_{c,{\rm eff}}$ and $\tilde{\Phi}=\Phi/\Lambda_{c,{\rm eff}}$ and extend integration limits due to the fast convergence of the integral. The first equation in (LABEL:eq:self-consis-colmb-1d) becomes

[TABLE]

which is just the same as Eq. (43) upon switching $\tilde{\Lambda}+\tilde{\Phi}$ to $\tilde{\Lambda}$. The second equation in Eq. (LABEL:eq:self-consis-colmb-1d) can be written as

[TABLE]

Eq.(51) shows that $\tilde{\Delta}$ in terms of $\tilde{\Lambda}+\tilde{\Phi}$ is given by the same universal function that is plotted in Fig. 6. Therefore, with appropriate rescaling $\tilde{\Lambda}\to\tilde{\Lambda}+\tilde{\Phi}$ one can find dependence of $\Delta$ on the bare $\Lambda$. For instance, from Fig. 6 we know that the order parameter $\tilde{\Delta}$ turns to zero at the critical point $\tilde{\Lambda}+\tilde{\Phi}=1$. Eq. (52) shows that at this point $\tilde{\Phi}$ is also zero, and therefore we conclude that $\tilde{\Lambda}=1$, that is $\Lambda=\Lambda_{c,\text{eff}}$, is the bare critical $\Lambda$ in this case too.

The same logic shows that the maximum value of $\Delta$ does not depend on the renormalization of $\Lambda$ and is given by $\Delta_{\rm m}\approx 1.65\Lambda_{c,\text{eff}}$.

Finally, similarly to Sec.IV.1 there too is a region where $|\tilde{\Lambda}+\tilde{\Phi}|\gg\Lambda_{c,\text{eff}}$ and the order parameter is exponentially small. If $g/g_{\text{eff}}-1$ is very small the asymptote in terms of bare $\Lambda$ has the same form as in the case with no Coulomb interaction, because then $\tilde{\Lambda}+\tilde{\Phi}\sim\tilde{\Lambda}$. In the other limit, when $g_{\text{eff}}\ll g$, although we get a different asymptote as far as $\Lambda$-dependence is concerned, we still get a similar $g$-dependence in $\Delta$ that goes as $\Delta\sim\exp(-\text{const}/g_{\text{eff}})$.

To summarize, the case of strongly screened Coulomb interaction is very similar to that with no Coulomb interaction. $\Lambda+\Phi$ plays the role of effective $\Lambda$ where, thanks to the strong screening, $\Phi$ slightly modifies the bare $\Lambda$.

IV.2.2 Weak screening: $dp_{0}\gg 1$

Now we consider the weak screening limit where distance to the screening gate $d$ is large in comparison with $1/p_{0}$, Coulomb interaction is only weakly screened and therefore repulsion between electrons is strong.

From the results of Section (III.2), it is apparent that even in the $dp_{0}\gg 1$ case one should distinguish two regimes: $d\lesssim\zeta_{0}\equiv\pi E_{R}/E_{c}p_{0}$, and $d\gtrsim\zeta_{0}$. As we saw in Section (III.2), $\zeta_{0}$ plays the role of the size of the two-particle bound state on its appearance for $d\to\infty$. In the $d\lesssim\zeta_{0}$ regime, which can be called the regime of intermediate screening, one can still introduce an effective local coupling constant, see Eq. (33). While the electrostatic effects, described by $\Phi_{\bm{p}}$, are strong in this case, the basic physics is very similar to the strong screening case. Therefore, in what follows we concentrate on the $d\gtrsim\zeta_{0}$ limit of really weak screening.

We start the analysis of Eqs (LABEL:eq:self-consis-colmb-1d) with noting that the $p$-dependence of the $\Phi$ field can be neglected. To see that, note that for large $dp_{0}$ the maximum of function ${\cal K}(p,p^{\prime})$ goes like $\ln(d)$, see Eq. (64) of Appendix A. Therefore the $\cal K$-term in the second bracket in the equation for $\Phi$ (LABEL:eq:self-consis-colmb-1d) can be dropped in comparison with the first, $d$ term, and the $\Phi$-equation simplifies to

[TABLE]

where the particle density $n$ is given by

[TABLE]

Next, let us assume that $\Delta_{p}\ll\xi_{p}+\Phi$. We are going to show that this assumption holds as long as the pair condensate is dilute. Under this assumption, the first of Eqs. (47) can be linearized as

[TABLE]

We then notice that the above equation is equivalent to the single pair Schrödinger equation (17) upon substituting the following ansatz for $\Delta_{p}$:

[TABLE]

where $\delta$ is a normalization constant to be found from Eq. (53). The proposed above ansatz works if one identifies

[TABLE]

The condition (57) requires $E_{b}<0$, since $\Phi$ always has a positive value proportional to the density (see Eq. (53)). In other words, we can say that if the single pair Scrödinger equation (17) has a solution with $E<0$, then through (56) and (57) one can construct a solution for the linearized BCS equation (55) out of it. Note that for the case of $\Lambda>0$, such a solution for (17) is automatically a bound state. That is because for $\Lambda>0$, the bottom of the continuum is always above the chemical potential. We therefore conclude that $\Phi=-E_{b}/2$, where, as before, $E_{b}$ is the two-particle bound state energy counted from the chemical potential for pairs. Since on the $\Lambda>0$ side $\Phi\neq 0$ corresponds to the presence of a superconducting condensate, we conclude that

[TABLE]

is the equation for the superconductor-insulator transition line, $g_{c}(\Lambda)$. Therefore, it can be obtained from the solution of the two-particle problem, described in Section III.2. The result of numerical solution of this problem is shown in Fig. 7, on the $\Lambda_{\rm eff}>0$.

Now we show that our approximation, $\Delta_{p}\ll\xi_{p}+\Phi$, equivalent to $\delta\alpha_{p}\ll 1$, is consistent as long as the condensate is dilute. To see that, consider the ansatz (56) and use the approximation, $\Delta_{p}\ll\xi_{p}+\Phi$, to write the density (54) as

[TABLE]

where $\delta\ll 1$, as appropriate for the onset of superconductivity, and $\alpha_{p}$ is a normalizable bound state wavefunction of the single pair Schrödinger equation (17). This implies that $\alpha_{p}$ is peaked around $p_{0}$, and has a finite width proportional to the inverse of the pair-size, $\zeta$. Hence, one can estimate $\int\alpha_{p}^{2}dp\sim\alpha_{p_{0}}^{2}/\zeta$. Applying this to Eq. (59), we conclude that in order for the ansatz Eq. (56) to be valid, we must have

[TABLE]

At the point of transition, $g=g_{c}$, the bound state energy of a single pair, as well as density-dependent quantities $\Phi$ and $n$, are zero, while the pair has a finite size for $\Lambda>0$, which means (60) is valid at the transition point. As $g$ increases above $g_{c}$ the HF field $\Phi$ and the density $n$ both increase with $g-g_{c}$, while the pair-size slowly shrinks. Therefore, we conclude that in the case of $\Lambda>0$, one can think of the superconducting state as a dilute condensate of well-defined pairs as long as the system is in the vicinity of the transition point. Fig. 8 shows that the solution constructed by the ansatz (56) matches the numerical solution of Eq. LABEL:eq:self-consis-colmb-1d very well.

Some insight into the shape of the $g_{c}(\Lambda)$ line for the superconductor-insulator transition on the $\Lambda>0$ side of the phase diagram can be deduced from Eqs. (38) and (40). Setting $E_{b}=0$ in those equations, we see that in the vicinity of $\Lambda=0$ the transition line has a finite slope, $g_{c}-g_{cb}\propto\Lambda$, while for $g_{c}-g_{cb}\gtrsim g_{cb}$ we have $g_{c}-g_{cb}\propto\sqrt{\Lambda}$. Both of these observations are confirmed numerically in Fig. 7.

Next we consider the case $\Lambda<0$, when the bare chemical potential is inside the band. This situation is illustrated in the insert of Figure 1. For $g=0$, one has a gas of interacting electrons, with density strongly reduced from the nominal value for the non-interacting case due to the strong electrostatic effects. We can determine the normal-state density by noting that since $\xi_{p}+\Phi=\kappa(p-p_{0})^{2}-|\Lambda+\Phi|$ in Eq. (54) changes sign as a function of $p$, there is a finite electron density even in the absence of superconductivity, $\Delta_{p}=0$. For $d\gg\zeta_{0}$, solving Eqs. (53) and (54) together leads to

[TABLE]

This shows that negative bare $\Lambda$ causes positive $\Phi$, the magnitude of which is comparable to $|\Lambda|$, so that the renormalized Fermi energy of the electrons counted from the bottom of the band, $|\Lambda+\Phi|$, is much smaller than the bare one, $|\Lambda|$. The Coulomb repulsion pulls the effective chemical potential down (decreases $|\Lambda+\Phi|$) so as to reduce the density of electrons. We can define a normal-state Fermi momentum, counted from $p_{0}$, associated with the normal state Fermi energy, $|\Lambda+\Phi|$:

[TABLE]

Qualitative features of the superconducting transition on the $\Lambda<0$ side depend on the relation between $p_{F}$ and $1/\zeta_{0}$, since $1/\zeta_{0}$ determines the spread of momenta around $p_{0}$ that participate in the formation of the two-particle bound state. For $p_{F}\lesssim 1/\zeta_{0}$, the transition is strongly affected by the formation of the bound state on the BEC side. In that case, one expects $\Delta_{p}$ with a strong momentum dependence, whose form is largely similar to the one obtained for $\Lambda=0$ at large momenta, $|p-p_{0}|\gtrsim p_{F}$ (see Fig. 9), with some modifications around the normal-state Fermi surface. In this regime, we expect a value slightly less than $g_{cb}$ for $g_{c}$, the reduction being due to the increased right hand side of the self-consistency equation – the first of Eqs. (LABEL:eq:self-consis-colmb-1d) – in the presence of zeros of $\xi_{p}$. For $p_{0}\gg p_{F}\gtrsim 1/\zeta_{0}$, the states that participated in the formation of the bound state are completely Pauli-blocked, and BEC physics is irrelevant. In this limit, one expects the pairing physics to be dominated by the Fermi surface, with only the value of $\Delta_{p}$ at the Fermi momentum being important. Since we are interested in the BEC-BCS crossover physics, we do not pursue the problem of determining the normal-superconducting phase boundary in full rigor, and present only numerical results here, see Fig. 7. Of importance to us is the undoubtful fact that the BCS phase does exist on the $\Lambda<0$ side, and the critical $g_{c}(\Lambda)$ separating the normal and superconducting phases lies below the BEC-BCS crossover line. This latter line will be described below, and is precisely determined numerically.

To see how the $\Delta_{p}$ solutions of Eq. LABEL:eq:self-consis-colmb-1d behave on the $\Lambda<0$ side, we present a typical numerical solution for $\Lambda=-E_{R}$ in Fig. 9. The lowest panel corresponds to the smallest $g$ where $\Lambda+\Phi<0$ and we are in the BCS region; the middle panel belongs to the boundary case where $\Lambda+\Phi=0$, which corresponds to the BEC-BCS crossover, see below; and finally for an even larger $g$ in the top panel $\Lambda+\Phi$ becomes positive and we are in the BEC region. Notice how in the BCS region $\Delta_{p}$ has two very sharp dips on the Fermi momenta and as we go into the BEC region it becomes more like the solutions on the positive side of $\Lambda$ presented in Fig. 8.

It is also noteworthy that we are still within the validity of the moat-band model even for $\Lambda=-E_{R}$, since the normal-state doping level is much smaller than $|\Lambda|$, see Eq. (61). This tendency is seen in Fig. 9: while the bare value of $\Lambda=-E_{R}$ in all three cases shown there, the renormalized values of $|\Lambda+\Phi|$ are indeed much smaller than $|\Lambda|$, and are in agreement with Eq.(61).

The BEC-BCS crossover line can be deduced from the same considerations as in the case of neutral fermions, except the condition $\Lambda=0$ should be replaced with $\Lambda+\Phi(g,\Lambda)=0$, which implicitly defines a line in the $(g,\Lambda)$ plane. It is clear that this line exists only on the $\Lambda\leq 0$ side of the phase diagram, and starts from point $(g_{cb},0)$. The rest of the line can be obtained numerically in the following way: first, one sets $\Lambda+\Phi=0$ in the equation for $\Delta_{p}$, and thus obtains a solution for every $g>g_{cb}$. Having obtained a solution for $\Delta_{p}$ for a particular $g$, one then substitutes this solution into the equation for $\Phi$, Eq. (53), in which $\Phi$ itself in the left hand side is replaced with $|\Lambda|$. As a result, one obtains $|\Lambda(g)|$ that defines the BEC-BCS crossover line. The result of implementing this program is shown in Fig, 7 as a dashed line on the $\Lambda_{\rm eff}<0$ side. This concludes the construction of the phase diagram for the case of charged fermions.

We leave the comprehensive study of the model at finite temperature for future work. Here we would like to provide a simple estimate of the magnitude of the temperature at which superconductivity develops. The critical temperature of a two-dimensional superconductor is controlled by the Berezinskii-Kosterlitz-Thouless (BKT) physics which describes binding of vortex-antivortex pairs of the superconducting phase. Typically, the BKT transition temperature $T_{\rm BKT}$, which describes the onset of the two-dimensional phase coherence, is much lower than the mean-field theory result $T_{\rm mf}$, roughly set by the magnitude of the superconducting order parameter. In case of neutral fermions with moat dispersion and no Coulomb repulsion Ref. He and Huang, 2012 finds that $k_{B}T_{\rm BKT}=\hbar^{2}\pi n/(4m_{\rm B})$. Here $n$ is the superconducting density, which controls superfluid stiffness, and $m_{\rm B}$ is the effective mass of the pairs, which approaches to $4m$ in the limit of large SOC, where $m=1/2\kappa$ in our notation. For an estimate, we assume a favorable set of parameters: material with large SOC, $E_{R}=200{\rm meV}$ Ast et al. (2007), small effective mass $m=0.05m_{e}$, and moderate attractive interaction $g/\kappa=0.1$. For $\Lambda=0$ numerical solution of (LABEL:eq:self-consis-colmb-1d) and (54) gives $n\approx 0.004k_{0}^{2}$, where the moat wave vector $k_{0}=\sqrt{2mE_{R}}/\hbar\approx 5.1\times 10^{8}{\rm m}^{-1}$. This leads to a large critical temperature $T_{\rm BKT}=\pi nE_{R}/(8k_{B}k_{0}^{2})\approx 3.5$ Kelvin. As discussed in the bulk of the paper, Coulomb interaction between electrons strongly suppresses superfluid density $n$. Indeed, turning on finite $E_{c}=e^{2}k_{0}/(2\epsilon)$, even with a large $\epsilon=100\epsilon_{0}$ Goldman (2014) and intermediate distance to the gate $dk_{0}=1$, strongly reduces the density. We find $n\approx 9\times 10^{-5}k_{0}^{2}$ and $T_{\rm BKT}\approx 0.08$ Kelvin. Increasing distance to the gate to $dk_{0}=10$ suppresses superfluid density completely, $n=0$, and removes the superconducting state altogether. This simple estimate shows once again the importance of the proper treatment of the weakly screened electrostatic repulsion in the problem of dilute two-dimensional electron gas with moat dispersion.

V Conclusion

Two-dimensional electron gas with a moat dispersion represents convenient model system for studying a superconductor - band insulator transition in the limit of low density of charged carriers. Moat dispersion leads to the enhanced density of states which strongly amplifies effects of pairing and repulsive interactions.

BEC-BCS crossover as a function of the strength of attractive interaction $g$ on the negative-$\Lambda$ axis represents one of the key results of our study. Another qualitatively novel feature is the appearance of the bound electron pairs of finite size, Eq.(30), driven by the competition between local attractive $g$ and long-ranged repulsive Coulomb interaction.

While we model the attraction phenomenologically as a contact-like interaction of BCS type, several recent studies have pointed out that in the metals with strong spin-orbit coupling the attraction may be induced by the spin-orbit interaction itself Mahmoodian and Chaplik (2018); Gindikin and Sablikov (2018). We leave investigation of this interesting possibility for the future.

Our study highlights importance of electrostatic contributions to the electron dispersion in the low density limit. Our findings are of relevance to two-dimensional superconducting heterostructures and surface states where spin-orbit interaction is particularly prominent.

Acknowledgements.

We thank Eugene Mishchenko for numerous insightful discussions at the beginning of this investigation. This work is supported by the NSF grants DMR-1507054 (H.A. and O.A.S.) and DMR-1853048 (DAP).

Appendix A Coulomb kernel, and the effective potential in “real” space

In this Appendix, we derive the approximate expression for the Coulomb kernel, Eq. (24), and the effective 1D potential energy for the effective 1D motion that describes the radial problem in the SPWF problem, Eq. (26).

The full form of the Coulomb kernel is given in Eq. (22). For convenience, we repeat it here:

[TABLE]

Since we are interested in ${\cal K}(p)\equiv{\cal K}(p_{0}+p/2,p_{0}-p/2)$ with $p\ll p_{0}$, we can neglect the $\cos^{2}(\theta_{\bm{q}k}/2)$ term in the definition of ${\cal K}(q,k)$.

We now calculate ${\cal K}(p)$ for several important values ot its argument, starting with the value of ${\cal K}(0)$. To this end, we set $q,k\to p_{0}$ in Eq. (63), introduce a new variable $t=\sin(\theta/2)$, such that $d\theta\approx 2dt$, to obtain

[TABLE]

where $\beta=e^{\gamma}$ with $\gamma\approx 0.58$ being the Euler-Mascheroni constant, and to write the second approximate equality, we neglected $O(e^{-2dp_{0}})$ terms, which are negligible in the weak screening limit, $dp_{0}\gg 1$.

We now turn to the behavior of ${\cal K}(p)$ for $1/d\ll p\ll p_{0}$. In this limit, the exponential term in the definition of ${\cal K}(q,k)$, Eq. (63), can be neglected, and we arrive at

[TABLE]

The accuracy of the final expression here is $O(p^{2}/p_{0}^{2})$.

Given the two limits, Eqs. (64) and (65), we arrive at Eq. (25) of the main text as an interpolating function that connects them, and serves as a good approximation to the full Coulomb kernel.

We now discuss the effective 1D potential that corresponds to as explained in Section III.2.2, the effective potential energy corresponds to the Fourier transform of the Coulomb kernel, ${\cal K}(q)$, see Eq. (24). It is thus given by

[TABLE]

By integrating by parts, this integral is reduced to

[TABLE]

At this point one can substitute ${\cal K}(p_{0}+q/2,p_{0}-q/2)\to{\cal K}(q)$, which is now appropriate, since it yields a convergent result. This way we obtain

[TABLE]