Pairing transition in a double layer with interlayer Coulomb repulsion

We study the effect of interlayer Coulomb interaction in an electronic double layer. Assuming that each of the layers consists of a bipartite lattice, a sufﬁciently strong interlayer interaction leads to an interlayer pairing of electrons with a staggered order parameter. We show that the correlated pairing state is dual to the excitonic pairing state with a uniform order parameter in an electron-hole double layer. The interlayer pairing of electrons leads to strong current-current correlations between the layers. We also analyze the interlayer conductivity and the ﬂuctuations of the order parameter, which consists of a gapped and a gapless mode.

Layered electronic systems have attracted substantial attention over several decades because new physical effects can be observed, which neither exist in a single layer nor in an isotropic three-dimensional systems.The role of the layered structures might be important for physical systems ranging from high-T c superconductors [1] to new quantum phases in twisted bilayers with a "magic angle" [2,3].Another direction of recent research is associated with multilayer graphene [4] and transition metal dichalcogenide multilayers [5], where an anomalous giant magnetoresistance [6], superconductivity [7], and the formation of exciton condensates [8][9][10][11][12][13] have been discussed and observed.Other interesting effects in layered systems are based on the application of an external magnetic field.In recent experiments the emergence of novel correlated many-body states between quasiparticles from different layers was observed [14,15].In this article we propose electron interlayer pairing caused by a repulsive interlayer interaction.As we show in a mean-field calculation, there is a second-order phase transition from two uncorrelated Fermi liquids in the two layers to a correlated pairing state if a critical interaction strength is exceeded.Moreover, we discuss a duality transformation between the pairing states of an electronic double layer and the excitonic pairing state in an electron-hole double layer.

I. MODEL
In an isolated (e.g., graphenelike) two-dimensional layer the electrons are subject to hopping when we assume that the Coulomb interaction within the layer is screened and renormalizes only the hopping parameters [16][17][18].In the case of two parallel layers there is also a Coulomb interaction which acts between the electrons of the two layers.This interlayer Coulomb interaction can be adjusted by inserting a dielectric material between the layers.Interlayer tunneling is ignored here to demonstrate only the effect of the interlayer Coulomb interaction.It may have some effect on the form of the order parameter though, as it was observed in the case of an attractive interlayer Coulomb interaction [19].This model has some similarity with exciton models in double layers [20], where we have electrons in one and holes in the other layer.While the Coulomb interaction between electrons and holes is attractive, the Coulomb interaction is repulsive in the electronic case.This implies that we do not expect the formation of excitonic Cooper pairs but a collective state built by electron pairs in the two layers, provided the Coulomb interaction is strong enough and the thermal fluctuations are weak at sufficiently low temperatures.Then the electronic Hamiltonian consists of the two independent hopping terms H ↑ (upper layer) and H ↓ (lower layer) and the repulsive interaction of the electrons between the two layers H I with where r and r are lattice sites of each layer and n sr is the local electron density operator in the layer s = ↑, ↓.The interlayer Coulomb interaction U r,r is short ranged in the lateral direction due to screening [21,22] and can be approximated by the diagonal matrix U r,r ≈ gδ r,r /2, where g is an effective parameter that decreases with the distance between the layers and depends on the interlayer dielectric material.This approximation is valid for interlayer distance that is large in comparison with the lattice constant of the layers, as described in Appendix A, where also the dependence of g on the various physical parameters is given.The resulting total Hamiltonian H = H ↑ + H ↓ + H I has the structure of the Hubbard Hamiltonian, provided that we interpret the layer index as the z component of the electronic spin.For a bipartite lattice (e.g., a honeycomb lattice or simple square lattice) at half filling we introduce a sublattice representation of the coordinates r = (R, j), with j = 1 and 2 for sublattices A and B, respectively.The coordinates R refer to sublattice A only.Then the hopping Hamiltonians H s (s = ↑, ↓), in terms of fermionic creation (annihilation) operators ĉ † s;R (ĉ s;R ), read as where σ j=1,2 are Pauli matrices.The fermionic annihilation operators are written as column vectors ĉs;R = (c s;R , d s;R ) T , whose upper (lower) component refers to sublattice A (B).The interaction term (1), together with the diagonal approximation, then becomes

A. Mean-field analysis
As a possible ansatz for the local-order parameter we consider a staggered order parameter with opposite signs on the two sublattices: where a global phase of reflects the global U(1) invariance of the model.A special mean-field solution fixes this phase.
For simplicity, we can choose a real positive for the subsequent calculation.The staggered order parameter is similar to the order parameter of an antiferromagnetic phase in the Hubbard model.The staggering sign translates into a σ 3 Pauli matrix in the corresponding quadratic mean-field Hamiltonian (cf.Appendix B): In the staggered order parameter we can replace (−1) j by a general phase factor exp(iϕ j ) with some phases ϕ j .It turns out, though, that only ϕ j = π j gives a stable mean-field solution.The Fourier-transformed kernel matrix in Eq. ( 5) is a 4 × 4 matrix with the twofold degenerate eigenvalues ±E q , with E q = √ 2 + h 2 1;q + h 2 2;q (cf.Appendix C).With this result we can treat the grand-canonical free-energy in the mean-field approximation.It is plotted for the subcritical interaction strength g < g c in Fig. 1(a) and for g > g c in Fig. 1(b).Thus, depending on the interaction strength g, the free energy has either a single minimum at = 0 for interaction g g c or a continuously degenerate minimum for g > g c .The continuous degeneracy implies a gapless Goldstone mode.The corresponding phase transition and its consequences for transport properties and the gapless fluctuations of the order parameter are discussed subsequently.
To obtain the value of we solve the mean-field equation δ F = 0.This reads for solutions = 0 as where ρ(E ) is the density of states (DOS).The solution of this equation provides us the value of as a function of the temperature and the coupling strength g.This is visualized in Fig. 2(a).Moreover, the critical temperature T c for the pairing transition is plotted as a function of g in Fig. 2(b).
Here we compare the electron dispersion of a honeycomb lattice (with a linearly vanishing DOS at E = 0) and of the square lattice (with a logarithmically divergent DOS at E = 0).The honeycomb lattice has a nonzero threshold of g for a pairing transition, whereas a square lattice has always a nonzero transition temperature for any nonzero interaction.
In contrast to the constant DOS ρ, which gives the critical temperature T c ∝ t exp(−1/gρ ), the logarithmic DOS causes a renormalization of the effective coupling gρ to larger values.It is also possible to use other dispersions to enhance the critical temperature, e.g., with a van Hove singularity at the Fermi level [20].

B. Duality transformation
It should be noticed that we can replace the electrons in the lower layer by holes and the repulsive interaction by an attractive interaction g → −g.The resulting hopping Hamiltonian H ↑ − H ↓ − H I describes electrons in the upper layer and holes in the lower layer which attract each other through an attractive interlayer Coulomb interaction.This duality transformation can also be applied to the mean-field Hamiltonian H MF .It amounts to the replacement of the staggered order parameter of the electronic system by a uniform order parameter for the electron-hole system.The latter means formally replacing the Pauli matrix σ 3 by the 2 × 2 unit matrix σ 0 in Eq. ( 5).This leads to the following mean-field Hamiltonian of the exciton gas [22]: where ĉ † ↓;R is the creation operator for a hole (i.e., an annihilation operator of an electron) in the lower layer.This Hamiltonian has the same dispersion E q as the Hamiltonian H MF of the electronic double layer, implying that the two problems can be mapped formally onto each other.In other words, both systems describe pairing of particles, where the main difference is that the order parameter is uniform in the electron-hole double layer and staggered with respect to the two sublattices in the electron-electron double layer.Here it should be pointed out that a modulated order parameter is a common phenomena in many condensed-matter systems, which can be of different origins.Typical cases are the Fulde-Ferrell-Larkin-Ovchinnikov phase in superconductors [23,24] or the formation and melting of the Wigner crystal phase in bilayer systems, which has been intensively studied in Refs.[25][26][27][28].

C. Current-current correlations
To analyze the effect of the interlayer Coulomb interaction we create a local current density at time t = 0 at site R on the lower layer which oscillates with frequency ω.Then we measure the local current at time t > 0 at site R on the upper layer, in analogy to the drag effect [21,22].This measurement probes how much a local current in the lower layer will generate currents in the upper layer.It can be described in terms of linear response theory, which is associated with the current-current correlator between the two layers and reads with a positive parameter α.For low temperatures the Fourier components of the correlator become the following (cf.Appendix D): with Its real part is plotted in Fig. 3 and exhibits a singularity at q sing (ω) with The singular wave vector is visualized in Fig. 3(b).It corresponds with an inverse length scale (wavelength) for the spatial distribution of the current in the upper layer, caused by the local current in the lower layer with frequency ω.As long as the frequency is less than the gap 2 , there is a mode with a finite wavelength.When ω is equal to the gap, this wavelength diverges and the strongest correlation appears with a zero mode q sing (2 ) = 0.And finally, when ω exceeds the gap, q sing (ω) becomes imaginary, which indicates an exponential spatial decay of the strongest correlation.The ω dependence of the q = 0 component of the currentcurrent correlator also has a characteristic behavior at the gap:

D. Interlayer conductivity
From the current-current correlator we obtain the interlayer conductivity tensor as [29] This gives the following for the diagonal conductivity: with the Fermi-Dirac function f β (E k ) at inverse temperature β.The real and imaginary parts read at low temperatures as Im The real part of the conductivity is zero when the frequency is less than the gap.But when ω reaches the gap, it jumps to a finite value, indicating that there is a current between the two layers when the photon energy hω is equal or larger than the gap.This jump is similar to the jump observed for pairing in electron-hole layers at ω = 0 and T > 0 [30].The imaginary part, on the other hand, is always nonzero in the pairing phase and has a logarithmic divergence when the real part jumps (cf.Fig. 4).

E. Quantum fluctuations of the order parameter
Quantum fluctuations around the mean-field solution consist of two spectral branches in the low-energy sector: a gapless branch due the degenerate ring of Fig. 1 and a gapped branch.The dispersion of both branches can be evaluated within the low-energy approximation and a gradient expansion p ∼ 0 (cf.Appendix E): indicating a stable mean-field solution.The result for a honeycomb lattice is visualized in Fig. 4(b).

F. Discussion and summary
In an electronic double layer there is an interlayer pairing transition due to a strong Coulomb interaction.Its critical temperature T c increases almost linearly with large coupling strength (cf.Fig. 2).For weak coupling the behavior depends on the details of the DOS: On the honeycomb lattice T c = 0 if g < g c and for the square lattice T c vanishes exponentially with 1/g.The pairing phase is accompanied by strong interlayer current-current correlations.These correlations diverge if the relation q 2 + ω 2 = 2 is satisfied (cf.Fig. 3).This reflects a long-range current-current correlation if ω 2 = 2 holds.The real part of the interlayer conductivity is nonzero only for frequencies ω larger than the gap of the pairing phase.Although there is a divergent current-current correlation for ω less than the gap, the real part of the interlayer conductivity σ μμ is zero.Due to this characteristic behavior the drag effect would be a good candidate to identify experimentally the pairing phase in electronic double layers.The effective coupling strength can be tuned by the interlayer distance and by the coefficient of the interlayer dielectric material, which enables the experiment to drive through the pairing transition.Moreover, the lattice structure of the layers determines the density of states.This could be used, for instance, by applying a strain field, to change the pairing transition temperature.And finally, the duality relation provides a basis to compare the large number of experiments on excitonic double layers with  17) and (18).(b) Gapless and gapped spectral branches of the order parameter fluctuations, calculated for = 0.125t (blue lines), = 0.25t (black lines), = 0.375t (red lines), and = 0.5t (green lines).Only positive energies at positive momenta are shown.The momentum K BZ is the square root of Brillouin zone area.(For paper printed black-white version: smallest has been used to calculate the lowest gapped and upper higher gapless curve.)more recent experiments on electronic double layers in order to clarify some of the recent experimental results.
where s and s refer to the layer index and ĉ1,s = c s , ĉ2,s = d s , and correspondingly for ĉ † i,s as defined in the main text, it is easy to show the following identity which is valid at every lattice site R (from right to left): where σ μ=1,2,3 are Pauli matrices and σ 0 is a two-dimensional unity matrix.In the mean-field approximation we put densities in both layers at the same value (zero at half filling) and introduce the following order parameters: The only order parameter which has a stable mean-field solution is 3 , which is exclusively considered in this paper.Because of the U(1) symmetry of the problem, there are infinitely many solutions which describe thermodynamically one and the same system.In the special solution which is discussed in the main text, we consider a purely real order parameter and call it .Equation (B2) is easily generalized to generic nonlocal interactions on bipartite lattices: Here we can prepare the product of two densities as which is valid for all combinations of indices, i.e., for s = s and R = R too.

APPENDIX C: EIGENBASIS OF THE MEAN-FIELD HAMILTONIAN
In general, the Fourier-transformed kernel matrix of the mean-field Hamiltonian Eq. ( 5) reads as K q = j h j;q σ j e −iφ σ 3 e iφ σ 3 j h j;q σ j .

APPENDIX D: CURRENT-CURRENT CORRELATIONS
We expand the correlator of Eq. ( 9) in terms of the eigenstates {|n } with the corresponding eigenvalues {E n } of the mean-field Hamiltonian Eq. ( 5) (cf.Appendix C).This gives which can be expanded for small q to give Cμμ (q; ω) ≈ This correlation is nonzero only in the gapped phase, while it vanishes in the gapless phase.The angular integration is performed using the relation to yield at low temperatures Eq. (10).|BZ| n,m [ f β (E n,q ) − f β (E m,q )] n, q + p|K (1)  p |m, q m, q|K (1)  −p |n, q + p E n,q+p − E m,q − iα , (E2) Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s)and the published article's title, journal citation, and DOI.Open access publication funded by the Max Planck Society.

FIG. 1 .
FIG. 1. Mean-field potential.(a) Below critical interaction strength g c .(b) Above critical interaction strength g c .

FIG. 2 .
FIG. 2. (a) Solution of the gap equation for different temperature values from left to right: T /t = 0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, and 0.5.(b) Comparison of the critical temperatures as functions of the interaction strength for the square lattice [blue (gray in printed version) curve] and the honeycomb lattice (black curve) in units of the hopping amplitude t.

FIG. 3 .
FIG. 3. (a)The real part of the current-current correlator of Eq. (11) as a function of the momentum for different frequencies.(b) Position of the singularity in the correlator at the edge of the gap for different frequencies and momenta according to Eq. (12).

FIG. 4 .
FIG. 4. (a) Real and imaginary part of the conductivity according to Eqs. (17) and(18).(b) Gapless and gapped spectral branches of the order parameter fluctuations, calculated for = 0.125t (blue lines), = 0.25t (black lines), = 0.375t (red lines), and = 0.5t (green lines).Only positive energies at positive momenta are shown.The momentum K BZ is the square root of Brillouin zone area.(For paper printed black-white version: smallest has been used to calculate the lowest gapped and upper higher gapless curve.) )where f β (E n ) is the Fermi-Dirac distribution at inverse temperature β.Then the Fourier-transformed correlator becomesCμν (q; ω) = d 2 k (2π ) 2 n,m [ f β (E m,k ) − f β (E n,k )] n, k + q| ĵμ↑ |m, k m, k| ĵν↓ |n, k + q ω − iα − E m,k + E n,k+q , (D2) To obtain the spectra of order parameter fluctuations we assume small quantum deviations from the mean-field values * → * δ RR + B † RR and → δ RR + B RR .The effective Hamiltonian for the fluctuating order parameter B can be found from the second order of the Schrödinger perturbation theory if we treat the termK (1) RR = 0 B † RR σ 3 B RR σ 3 0 , i.e., K (1) p = 0 B † p σ 3 B −p σ 3 0 , (E1) as a small perturbation.(The first-order term is zero because of the mean-field condition.)In Fourier representation we have δH = 1 2 d 2 p |BZ| d 2 q