Strong-coupling theory of condensate-mediated superconductivity in 2D materials

We develop a strong-coupling theory of Bose-Einstein condensate-mediated superconductivity in a hybrid system, which consists of a two-dimensional electron gas with either (i) parabolic spectrum or (ii) relativistic Dirac spectrum in the vicinity of a two-dimensional solid-state condensate of indirect excitons. The Eliashberg equations are derived and the expressions for the electron pairing self-energy due to the exchange interaction between electrons mediated by a single Bogoliubov excitation (a bogolon) and the bogolon pairs are found. Furthermore, we find the superconducting order parameter and estimate the critical temperature of the superconducting transition. The critical temperature reveals its linear dependence on the dimensionless coupling constant. It is shown, that the bogolon-pair-mediated interaction is the dominant mechanism of electron pairing in hybrid systems in both the weak and strong coupling regimes. We calculate the effective bogolon-electron interaction constant for both parabolic and linear electron dispersions and examine the dependence of the critical temperature of electron gas superconducting transition on exciton condensate density.


I. INTRODUCTION
The first microscopic description of the superconducting (SC) state belongs to the celebrated Bardeen-Copper-Schrieffer (BCS) theory [1][2][3][4]. Using this theory, one can explain the emergence of the SC gap and estimate the critical temperature of SC transition T c . Later, Green's functions formalism-based Migdal-Eliashberg (later Eliashberg) theory was developed [5][6][7]. It provides a more rigorous and accurate basis for the estimation of T c and allows to account for the Coulomb repulsion between electrons. Generally, by employing the McMillan's approach [8], the formation of Cooper pairs by electron-phonon interaction can be solely determined by the effective electron-phonon coupling strength λ. As compared with the BCS theory, which only works in the weak coupling regime (of λ 1), the Eliashberg theory is more general. In particular, it avoids the Debye frequency cutoff, thus allowing to extend the estimations to the strong coupling regime (of λ close to or larger than unity).
The price to pay is that the analytical and numerical calculations by the Eliashberg theory are usually more complicated than the ones by the BCS. It requires the solution of multiple coupled equations, which is often tricky, and some assumptions are in order. Since 1960s, various numerical approaches implying different special assumptions and simplifications have been developed in order to find the approximate solutions of the Eliashberg equations. Very recently, there have been suggested powerful density functional-based techniques (called the density functional theory for superconductors) [9][10][11]. In general, the phonon-mediated superconductivity is well studied, though there still exist many open questions, in particular, regarding the SC transi-tion in novel two-dimensional (2D) materials and their stacks due to phonon-mediated pairing.
However, electron-phonon interaction is not the only possible route to the formation of Cooper pairs. Among others, there have been reported several proposals for exciton [12][13][14], exciton-polariton [15][16][17], and cavityphoton-mediated [18] coupling mechanisms between electrons. In our recent works [19,20], we have developed a BCS-like theory for the indirect exciton condensatemediated superconductivity in hybrid systems consisting of layers of a (Bose-condensed) exciton gas and a two-dimensional electron gas or graphene. Employing the Bogoliubov theory, we considered quasiparticle excitations above the exciton condensate called bogolons. They interact with the electron gas via Coulomb forces, which distinguishes this problem from the phononassisted electron-electron interaction in conventional superconductors.
However, all the proposals listed above are based on the BCS-like approach. Thus, the weak coupling between electrons and phonons (or other particles like bogolons) is assumed, which means that λ should, strictly speaking, be small. In this paper, we build a strong-coupling theory for bogolon-mediated superconducting in hybrid Bose-Fermi systems. By employing the Green's functions technique, we calculate the electron-bogolon coupling strength λ and give an estimation for the critical temperature of SC transition. The bosonic subsystem which we consider here is indirect exciton gas. It can be substituted by other quasi-condensates like direct exciton gas or exciton-polariton condensates, which have been predicted [21][22][23][24][25][26] and studied experimentally [27][28][29][30]. The condensation in these systems has been reported at relatively high temperatures, even approaching the room temperature. The paper is organized as follows. In Sec. II, we introduce the Hamiltonian of bogolon-electron interaction. In Sec. III, we build the bogolon-mediated Eliashberg theory and derive the Eliashberg spectral function. Furthermore, in Sec. IV, we calculate the electron-bogolon coupling constant λ and estimate the critical temperature by considering electrons with parabolic and linear dispersions. Finally, in Sec. IV we summarize the results.

II. SINGLE-BOGOLON AND BOGOLON-PAIR-MEDIATED INTERACTION
Let us consider a hybrid system, in which a layer of electron gas (2DEG) and a layer of a Bose-Einstein condensate are in the vicinity of each other in z-direction ( Fig. 1). For the BEC part, we consider indirect excitons, where the formation of a BEC has been recently reported [30,31]. Indirect excitons consist of electrons and holes residing in n-and p-doped layers which are separated in z direction. These layers can be made of, e.g., MoS 2 or WSe 2 materials separated by several layers of hexagonal boron nitride (hBN) [30,32]. The 2D electron gas (described by the field operator Ψ σ (r) with the position vector r and the spin σ = {↑, ↓}) and the exciton gas (described by the field operator Φ(r)) layers are also spatially separated by a hBN layer. The two types of particles, electrons and indirect excitons, are coupled through the Coulomb interaction (in standard form) [33,34], (1) Here, R is the in-plane position vector of the exciton center-of-mass position. We consider the case when the most of excitons are being in the ground state (in BEC). Then, in the weakly-interacting regime, we can write the exciton field operator as Φ(R) = √ n c + ϕ(R), where n c is the condensate density and ϕ(R) is the field operator for non-condensed excitons. Applying the Fourier and the Bogoliubov transformations, from (1), we can find the Hamiltonian for the one-bogolon and bogolon-pairmediated interactions [19,20,35,36] (putting = k B = 1 below for simplicity and restoring these constants later) where g p is the Fourier image of electron-exciton interaction; c p,σ and b p are the annihilation operators for the electrons and bogolons, respectively. The Bogoliubov coefficients are defined as [37] where M is the exciton effective mass, s = κn c /M is the sound velocity, κ = e 2 0 d/ 0 is the exciton-exciton interaction strength in the reciprocal space, e 0 is electron charge, is the dielectric constant, 0 is the dielectric permittivity, and ω p = sp(1 + p 2 ξ 2 h ) 1/2 is the spectrum of bogolons with the definition of healing length ξ h = 1/2M s.
To get an analytic form of the electron-exciton interaction, we disregard the peculiarities of the exciton internal motion (relative motion of the electron and hole in the exciton). In monolayers of transition metal dichalcogenides, the exciton binding energy is usually very large: it might even exceed the room temperature. Thus, the assumption that an individual exciton is in its ground state with respect to its relative electron-hole motion is reasonable, and only the exciton center-of-mass motion plays an important role. Then, the electron-exciton interaction in direct space reads with r ee = l 2 + (r − R) 2 and r eh = (l + d) 2 + (r − R) 2 ; here, d is an effective size of the boson, which is equal to the distance between the n-and p-doped layers in the case of indirect exciton condensate, and l is the separation between the 2DEG and the BEC as shown in Fig. 1. Performing the Fourier transformation, we find Using the conventional perturbation theory, one may argue that H 1 H 2 since the density of non-condensed particles is a small quantity as compared with the condensate density, i.e.
However, with the Bogoliubov coefficients, the summation such as v p + u −p in (2) can drastically change the situation. In the long wavelength limit, ξ h p 1, where we are interested in this work, we have ω p ≈ sp and u p ≈ −v p ≈ M s/p. Then, in this limit, one has v p + u p → 0 whereas u p v p → ∞. This suggest the contribution of H 2 may surpass the one in H 1 as will be shown in Sec. III.

III. THE ELIASHBERG EQUATIONS
We will work in the framework of the Nambu-Gor'kov formalism [38,39]. First, let us introduce the twocomponent field operator, where c kσ is the annihilation operator of an electron with the momentum k and spin σ. Then, the generalized Green's function in the (2 × 2) matrix form readŝ or using (8), we can write it aŝ (10) The diagonal terms in Eq. (10) represent the standard Green's functions of electron quasiparticles, whereas the off-diagonal terms are the Gor'kov's anomalous Green's functions. Performing the Fourier transform from the imaginary time domain to the Matsubara frequency representation, we havê To define the bogolons' Green's function, we first introduce the notation Then, the normal and anomalous free Green's functions of the bo- respectively. Switching to the Matsubara frequency domain, we find In the long-wavelength limit pξ h 1, for the Bogoliubov coefficients (4) and (5) we have u p ≈ −v p , and consequently, D (p, iω n ) ≈ −A (p, iω n ). Using the operators A p and A † p , the bogolon-mediated interaction with the electronic subsystem given in Eqs. (2) and (3) reads Henceforth, we perform a perturbation theory expansion over the weak interacting terms Eqs. (14) and (15) in order to obtain the Dyson equation for the electron Green's function. The general solution of the Dyson equation readsĜ is the unperturbed Green's function, σ ν=0,1,2,3 are the Pauli matrices, and ξ k is the electron dispersion measured with respect to the chemical potential, ξ k = ε k − µ.
Following the steps of the standard Eliashberg theory, we separate the electron self-energy into two terms: the usual Coulomb contributionΣ c and the electron-bogolon contributionΣ eb (in full analogy with the electronphonon contribution).
The Coulomb contribution is given bŷ where V (k−p) are the matrix elements of static screened Coulomb interaction between the electronic states k and p, andĜ od is the off-diagonal component of the Green's function [4]. The self-energies of electron-bogolon (Σ 1b ) and electron-bogolon-pair (Σ 2b ) interaction can be calculated by the Dyson equation up to first order, in accordance with the diagrams presented in Fig. 2, Electron self-energy diagrams. The double solid lines stand for the Green's functions of an electronĜ in (10). The zigzag lines denote the condensate particles. (Thus, each zigzag line gives the factor √ nc). The dotted lines represent bogolons (the particles excited from the condensate). The wiggly lines stand for the electron-exciton interaction gp. Panels (a-d) correspond to 1b process in (19), and panels (e-h) correspond to 2b process in (21). Physically, the diagrams in (a-d) describe the excitation of condesate particle to the non-condesed state by a moving electron, whereas (e-h) describe the condensate polarization due to the moving electron.
where p k and ω n are the Matsubara frequencies of the fermions and bosons, respectively. Then, in the longwavelength limit, we see thatΣ 1b → 0, because the nor-mal and anomalous Green's functions of the bogolon cancel each other out. For the 2b processes, we introduce the polarization operator P (p, iω n ), which reads In this work, we will restrict ourselves to the case when the contribution of N q terms is negligible. As we have shown in earlier works [19,20], the N q -containing correction only results in quantitative difference (more precisely, it results in an increase of the critical temperature of SC transition), and it is negligible when the size of the condensate is small. Then, the polarization operator simplifies to and the self-energy can be rewritten in the form, We have already demonstrated, that the self-energy contribution due to 2b process is dominant,Σ 2b Σ 1b . We want to note, that all other terms like three-and four-(and more) bogolon-mediated processes give smaller contribution since all these terms emerge in the perturbative expansion, where the small parameter is electron-exciton interaction strength g p . The termsΣ 1b(2b) are of the same order g 2 p , which is the leading order of the expansion. To proceed further, it is a common practice to decompose the matrix self-energy and rewrite it as a linear combination of Pauli matrices with scalar functions as coefficients [40], where Z(k, ip k ) is the mass renormalization function, χ(k, ip k ) is the energy shift, φ(k, ip k ) andφ(k, ip k ) are the order parameters. By the gauge transformation [4], we can set the order parameterφ to zero. Then, us-ing (16) and (17), the Green's function readŝ Carefully replacing (24) and (25) by the exact forms of self-energy in (18) and (23), we come up with the Eliashberg equations, where N F is the density of states per spin at the Fermi level and is the Eliashberg electron-bogolon spectral function [4] with The superconducting gap can now be found as For the 2b process, we have where H(ω) is the Heaviside step function (see the details of the calculation in Appendix A).
The equation ∆(k, ip k ) = 0 gives the normal-state solutions. The critical temperature T c can be defined as the highest temperature, for which ∆(k, ip k ) = 0. We will use it in what follows.

IV. RESULTS
In this section, we apply the Eliashberg theory to different hybrid systems. In particular, we consider 2DEG with parabolic dispersion case and the electron gas with the linear dispersion case.
First of all, let us simplify the Eliashberg equations (27)-(29) using the following approximations: (i) since the superconducting pairing mainly occurs within a narrow energy window around the Fermi surface, we restrict the consideration to the electrons with k f [4,40,41]. Then, we can put χ(k f , ip k ) = 0 and only solve Eqs. (27) and (29); (ii) we assume the anisotropy of the Fermi-surface is weak, and thus use the isotropic formulation of the Eliashberg equations. Then, we can relabel the scalar functions (for brevety): In the parabolic dispersion case (and under the approximations discussed above), Eqs. (27) and (29) read [42] (see the details of derivations in the Appendix B) Here, we introduce the dimensionless Coulomb interaction µ * c . By definition, it represents a double-average over the Fermi-surface (FS) of the term V (k − p) in (29), For large class of superconductors [4,40,43], µ * c is in the range 0.1 ∼ 0.2. We will use this value instead of calculating the double-average term.
The λ-function in Eq. (34) reads where m is a an integer number, which indicates the difference between two Matsubara frequencies; F(φ, m) is the Elliptic integral of the first kind; the density of state per spin is N F = m e /(2π) with the effective mass of electron m e ; and L is the effective size of the condensate introduced as the necessary cut-off. In order to calculate T c , we will employ the technique discussed in [4,42]. Since the transition temperature is defined as the point, when the energy gap is infinitesimally small, the value of T c can be found by setting ∆ = 0 in all of the denominators of Eqs. (34) and (35). This gives where n and v are the Matsubara frequencies' indices, ip n = π(2n + 1)T , and sgn(x) is the signum function. Now, the equations for Z 2n+1 is independent of the gap function ∆ 2n+1 . The critical temperature can be found by putting the determinant of the matrix to zero, det{M nv } = 0, where Furthermore, we can calculate the T c numerically.
Let us now consider a 2DEG with a linear dispersion in, e.g., a graphene layer [44][45][46][47]. For simplicity, we will disregard the spinor structure of the wave function assuming that we deal with the doped graphene with the Fermi energy sufficiently away from the Dirac point. Furthermore, we will disregard the contribution from different valleys [48,49]. This approximation is valid since the large non-zero wavevector (∼ k f ) strongly depresses the interaction g p in Eq. (7).
Given the assumptions discussed above, we come up with a similar system of equations as the one for the parabolic case: Eqs. (34) and (35). However, we find different formula for the coupling constant (38)   Let us, first, understand the principle dependence of the critical temperature on the coupling constant. Comparing Eqs. (37), (38), (39) and (43), we see that it is convenient to denote the following dimensionless parameters, Then, we can investigate the critical temperature in terms of b E = πT /sk f and λ p(l) 0 (Fig. 3). Let us compare the critical temperatures calculated using the Eliashberg and BCS theories. In BCS, it reads where ω D = M s 2 /2 is the frequency cut-off, and γ = exp C 0 with C 0 ≈ 0.577 the Euler's constant. The parameter χ in parabolic [19] and linear [20] cases reads Figure 3 shows that the Eliashberg and BCS theory curves converge well when λ p(l) 0 is small, as expected; with the increase of λ p(l) 0 , then discrepancies between the two theories become more and more pronounced. As λ p(l) 0 → ∞, the asymptotic limit gives a nearly linear dependence of the Eliashberg critical temperature on the coupling coefficient λ p(l) 0 . Such dependence is not typical for the conventional superconductors discussed in the early works [4,42], in which the dependence T c ∼ √ λ in the framework of the Einstein model was found.
The nearly-linear dependence which we find can be also checked analytically if we only consider the leading order of (40) and (41) (see the details in the Appendix D). The asymptotic behaviour reads Since the Elliptic integral converges faster as one increases b E , we can treat F(arccos φ 0 , b E ) as a constant. The inset plot in Fig. 3 shows an estimation of the asymptotic behavior of the dimensionless critical temperature with fixed F(arccos φ 0 , b E = 0.5). Figure 4 and Fig. 5 show the critical temperature of SC transition as a function of condensation density. For the parabolic dispersion, we use the parameters typical for MoS 2 , and for the linear dispersion case we use the parameters typical for graphene. Another interesting question is the dependence of electron-bogolon interaction constant on n c and n e . In the zeroth order of the electron-bogolon interaction constant, from Eqs. (44) and (45) we find This discrepancy in the parametric dependencies can be understood by recalling the difference between the densities of states in linear and parabolic cases. Let us address the remaining assumptions, which we used in the calculations. First, to derive Eqs. (38) and (43), in addition to the standard approximations discussed in Sec. IV, we also approximated the excitonelectron interaction by g p ≈ e 2 0 d 2 0 . Such a simplification is valid if k f d and k f l 1. This assumption imposes a restriction on the maximal allowed value of n e for considered distances d and l. If the separations d and l are in a nanometre scale, the density of electrons should be n e ≤ 10 13 cm −2 .
Second, the polarization operator which we consider in Eq. (22) does not include the N q terms. These terms are relatively small when the size of the condensate L BEC is small, as it has been shown in the analysis in the framework of the BCS theory [19,20], since these contributions give an integral truncation from L −1 BEC to k f /2. In other words, we find the lower bound for the SC gap and the T c . For electrons, we considered the Green's function with account of the finite temperature since we are interested in the critical temperature of the SC transition. As the earlier works on bogolon [19,20] and phonon [51]mediated SC point out, in the BCS theory the non-zero terms in (21) contribute to a non-monotonous temperature dependence, and they only increase the critical temperature. Thus, the simplifications we use are of the technical nature, and they do not qualitatively change our results or conclusions.

CONCLUSIONS
We developed the Eliashberg theory of the Bose condensate-mediated superconductivity in hybrid twodimensional Bose-Fermi systems, and showed that bogolon-pair-mediated pairing of electrons represents the dominant mechanism not only in weak but also in the strong-coupling regime, while single-bogolon pairing is suppressed and thus, does not play any significant role. We started with an analytic expression for the self-energy of electrons in a two-dimensional material due to their interaction with bogolons; then, we calculated the electronbogolon coupling constant in the cases of parabolic and linear electron dispersions, and presented the corresponding estimations of the critical temperature of superconducting transition, which turns out relatively high. It was demonstrated, that the critical temperature of the superconducting transition depends on the dimensionless coupling constant linearly, which is not common for the conventional superconductors. We expect our theory and the estimations to impact such research areas as low-dimensional superconductors, novel two-dimensional Dirac materials, and the density functional theory for mesoscopic superconductivity. Again, using the notation q = k − p and noticing that the integral over q only depends on the absolute value of q, we find λ (k, p, n, v) = N F κ 2 n 2 c |g q | This integral is analytical (for ab > 0) [52], Then, denoting α nv ≡ p n − p v we find λ (k, p, n, v) = N F κ 2 n 2 c g 2 k−p 4s 2 L 2 α 2 nν + (s|k − p|) where the approximation where F is the elliptic function of the first kind, and Finally, Performing a similar procedure, we obtain Eq. (35). Let us separately address the asymptotic limit of large λ p(l) 0 . Following the arguments discussed in [42], we can assume b E to be large in this case. Then, the leading order contribution of (37) or (43) is with the zero Matsubara frequency index, and we have neglected the terms λ(n) with |n| > 1. Then, Eqs. (40) and (41) are truncated. For Eq. (40), we find Z ±1 = 1 + λ(0), and for Eq. (41) we have Z ±1 ∆ ±1 = ∆ ±1 λ(0) + ∆ ∓1 λ(±1). This results in a self-consistent equation, Furthermore, we can assume F(φ 0 , 1 √ 1+b 2 E ) to be constant for simplicity since it converges to a constant when the temperature increases, as it is shown in Fig. 6. Then, we find the following estimation for T c , where the symbol F c means we have fixed b E as a constant inside F φ 0 , 1 √