Multi-band pairing of ultrarelativistic electrons and holes in graphene bilayer

Physics Letters A 374 (2009) 326–330

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Multi-band pairing of ultrarelativistic electrons and holes in graphene bilayer Yu.E. Lozovik ∗ , A.A. Sokolik Institute of Spectroscopy, Russian Academy of Sciences, 142190 Troitsk, Moscow Region, Russian Federation

a r t i c l e

i n f o

Article history: Received 3 October 2009 Accepted 9 October 2009 Available online 22 October 2009 Communicated by V.M. Agranovich PACS: 74.78.Na 74.20.-z 81.05.Uw Keywords: Graphene Excitons BCS Massless fermions

a b s t r a c t Multi-band pairing of effectively ultrarelativistic electrons and holes in asymmetrically biased graphene bilayer in strong coupling regime is considered. In this regime, the pairing affects both conduction and valence bands of the both graphene layers, and the order parameter is a matrix, which indices correspond to the bands. For band-diagonal s-wave pairing, we derive the system of multi-band gap equations for the gaps in the valence and conduction bands and solve it in the approximation of constant gaps and in the approximation of separable pairing potential. For a characteristic width of the pairing region of order of magnitude of the chemical potential, the gap values are not much different from single-band BCS estimations. However, if the pairing region is wider, then the gaps can be much larger and depend exponentially on its energy width. We also predict gapped and soliton-like oscillations of a relative phase of the gaps and unpairing of quarter-vortices at Kosterlitz–Thouless transition. © 2009 Published by Elsevier B.V.

1. Introduction Among various generalizations of the original Bardeen–Cooper– Schrieffer (BCS) theory of superconductivity [1] is its extension on multi-band case, when particles in several different bands are simultaneously involved into the pairing. The problems of multiband pairing can be divided into two physically distinct groups. In the first one, there are more than one Fermi surfaces in different bands, and the pairing involves particles from vicinities of these surfaces. Such a model of pairing was first considered in the seminal paper [2] and later was applied to describe superconductivity in a number of materials including MgB2 [3] and novel iron pnictide high-temperature superconductors [4,5]. In the second group, the Fermi surface occupies only one band, but the pairing affects also remote bands, not containing the Fermi surfaces. This is mainly the case of superconducting pairing of relativistic particles [6,7] (see also [8] for review), where antiparticle bands provide additional phase space for Cooper pair constituents and thus lead to increase of the gap with respect to a single-band case. Fabrication of graphene, carbon monolayer with two-dimensional honeycomb lattice [9], allows to study effectively ultrarelativistic fermions in condensed-matter systems [10]. This possibility is provided by its unique band structure: in the vicinity of the Fermi level of undoped graphene, electrons obey two-dimensional

*

Corresponding author. Tel.: +7 496 7510881; fax: +7 496 7510886. E-mail address: [email protected] (Yu.E. Lozovik).

0375-9601/$ – see front matter © 2009 Published by Elsevier B.V. doi:10.1016/j.physleta.2009.10.045

Dirac equation for massless particles and move with the Fermi velocity v F ≈ 106 m/s (see, e.g., the reviews [11]). Ultrarelativistic electron dynamics in graphene lead to multi-band picture of superconducting pairing [12,13]. In the present Letter, we consider multi-band pairing of spatially separated electrons and holes in graphene bilayer, a system consisting of two parallel asymmetrically gated graphene layers [14]. Nanometer-scale interlayer separation D can prevent recombination of electrons and holes, while their Coulomb attraction leads to BCS-like condensation and superfluidity of electron–hole pairs (similarly to electron–hole condensation and superfluidity in coupled semiconductor quantum wells, predicted in [15]). In [14], electron–hole pairing in graphene bilayer was considered in the weak coupling regime when the pairing is of single-band BCS type, whereas a possibility to achieve also strong-coupling regime experimentally was argued. However, on increase of the coupling strength, there is no crossover from BCS state to Bose–Einstein condensation (BEC) in a gas of local pairs (BCS-BEC crossover, typical to systems of attracting fermions [16]), since ultrarelativistic particles cannot form true bound states due to absence of gap in the energy spectrum [17] (although in magnetic field the crossover is possible even in graphene [18]). At strong coupling, the pairing involves both the bands of paired particles, containing their Fermi surfaces (conduction band of electron-doped layer and valence band of hole-doped one), and the remote bands (valence band of electron-doped layer and conduction band of hole-doped one), which is a manifestation of ultrarelativistic nature of the particles (Fig. 1). At the same time, the paired particles are condensed

Yu.E. Lozovik, A.A. Sokolik / Physics Letters A 374 (2009) 326–330

327

are non-zero, similarly to [6] (Fig. 1). A spin-valley structure of (1) is assumed to be factorized in a form of (4 × 4) unitary matrix P . One-particle excitations in the system under band-diagonal pairing have Bogolyubov-like dispersion laws and are represented by two branches, corresponding to particle and “antiparticle” excitations (see [6–8,12])

E + ( p) = Fig. 1. Band picture of multi-band pairing of ultrarelativistic electrons and holes. Thick horizontal lines show chemical potential positions μ and −μ in electron(left) and hole-doped (right) layers, and the region of pairing is marked by hatching. Under band-diagonal pairing, the electron–hole condensate is described by two anomalous Green functions F ++ and F −− .

in strongly overlapping Cooper pairs, as in usual BCS state, so we call the state formed in graphene bilayer at strong coupling as “ultrarelativistic BCS” state. In contrast to BEC of electron–hole pairs, ultrarelativistic BCS state can demonstrate rather high critical temperature and anomalous superfluid properties. Subsequent study of electron–hole pairing in graphene bilayer, presented in [19–23], has given various estimations of a superfluid transition temperature, however a possibility of multi-band character of the pairing was not considered in these works. Simple multiband pairing models in graphene with a contact interaction were considered in [24] for a problem of charge- and spin-density wave formation in single-layer graphene in parallel magnetic field, in the works [12,13] for a superconductivity in graphene and in [25] for electron–hole graphene bilayer. Here we perform multi-band pairing analysis taking into account a finite range of screened pairing interaction and using physically realistic estimations for an ultraviolet momentum cutoff in gap equations. Recent success in a fabrication of graphene bilayer [26] makes the problem under consideration to be of high importance. The Letter is organized as follows. We formulate the multi-band pairing model in Section 2 and then solve gap equations at zero temperature in two approximations: in the constant gap approximation (Section 3) and in the approximation of separable potential (Section 4). In Section 5 we discuss a Kosterlitz–Thouless transition temperature in the system and Section 6 is devoted to conclusions.

E − ( p) =

   

2

v F | p| − μ

2

v F | p| + μ

+ Δ++ (p ), + Δ−− (p ),

where Δ++ and Δ−− are the band-diagonal anomalous selfenergies (gap functions). Here we consider, for simplicity, equal concentrations of electrons and holes. Note that large gap values allows the pairing picture under consideration to hold even at rather high imbalances of electrons and holes chemical potentials of order of the gap in the spectrum Δ++ ( p F ). Here we do not take into account frequency dependence of selfenergies, assuming that the role of dynamical effects is reduced to some effective restriction of the pairing region in the momentum space, analogously to restriction of superconducting pairing in metals to the Debye frequency-sized region near the Fermi surface [1]. The anomalous Green functions (1) are expressed through the gaps Δγ γ as

F γ γ ( p , i εn ) = −

Δγ γ (p )

εn2 + E γ2 (p )

(2)

,

while the gaps are expressed through F γ γ by the self-consistent gap equations:



d p 

 p γ p γ  2 2 ( 2 π ) εn γ      × V p − p F γ  γ  p , i εn ,

Δγ γ (p ) = − T

(3)

where V ( q) is Fourier-transform of a potential of interlayer electron–electron interaction. The overlap factors  p γ | p γ   of inip γ  states in a scattering of ultrarelativistic tial | p γ   and final | particle result from folding over components of its spinor wave function, and their squares has the form     p γ p γ  2 = 1 + γ γ cos(ϕ − ϕ ) ,

2. Multi-band gap equations

(4)

2

The multi-band pairing can be described by the anomalous Matsubara Green functions



(1 )

(2)+



F γ1 γ2 ( p , τ ) = − T c p γ (τ )c p γ (0) , 1 2

(1)

where T is an ordering operator in the imaginary time τ and . . . means averaging over thermodynamic ensemble. The oper(1) (1) (2) (2) ators c p γ = a p γ , c p γ = a p −γ are connected with the destruction (i )

 located in the layer i operators a p γ of electrons with momentum p (i = 1 for electron-doped layer, i = 2 for hole-doped layer) and populating the band γ (γ = +1 for conduction band, γ = −1 for valence band) via electron–hole transformation in the hole-doped layer. Such a transformation allows the bands of both kinds of pairing particles, containing their Fermi surfaces, to be denoted by γ = +1 and named as “conduction bands”, and the bands remote from the Fermi surfaces can be denoted by γ = −1 and named as “valence bands”. The anomalous Green functions (1) are matrices with the indices γ1 and γ2 corresponding to conduction and valence bands, so various structures of these matrices are possible. Leaving this issue to further study (see also [7,27]), here we assume the simplest case of band-diagonal pairing, when only F γ1 γ2 and F γ+1 γ2 with γ1 = γ2

 and p rewhere ϕ and ϕ  are azimuthal angles of the vectors p spectively (see also more detailed consideration of the pairing in the multi-band regime, which will be published separately [27]). At s-wave pairing, when F γ γ ( p , i εn ) and Δγ γ ( p ) are indepen , the self-consistent gap equations (3), with dent of direction of p taking into account (2), (4), after integrating over ϕ  and summing over εn = π T (2n + 1), take the form:

  p dp V 0 ( p , p  ) + γ γ  V 1 ( p , p  )

Δγ γ ( p ) =

γ  =±

×

2π

2

E γ  ( p ) Δγ  γ  ( p  ) tanh , 2E γ  ( p  ) 2T

(5)

where the following angular harmonics of the pairing interaction are introduced:

V l ( p, p ) =

2π

dϕ 2π

V





p 2 + p  2 − 2pp  cos ϕ cos(lϕ ).

(6)

0

As long as V 0 = V 1 , Eqs. (5) for the two gap functions Δ++ ( p ) and Δ−− ( p ) are mixed and their relative phase is therefore fixed.

328

Yu.E. Lozovik, A.A. Sokolik / Physics Letters A 374 (2009) 326–330

However, when this fixation is weak, gapped and soliton-like excitations, corresponding to oscillations of the relative phase and being solutions of the sine-Gordon equation, can arise (similarly to [28]). The pairing interaction can be approximated by the Coulomb potential, statically screened in the random-phase approximation, which is well justified in graphene bilayer due to large number of fermionic flavors [21,23,29]. However, to study dynamical effects, we consider more general case of dynamically screened interlayer interaction [15,30]:

V (q, ω) =

v q e −qD 1 − 2v q Π(q, ω) + v q2 Π 2 (q, ω)[1 − e −2qD ]

,

(7)

where v q = 2π e 2 /εq is a bare Coulomb interaction, ε is a dielectric permittivity of a surrounding medium and Π(q, ω) is a polarization operator of graphene (it is the same for both layers due to the particle-hole symmetry). Hereafter we consider the most favorable case for pairing, when p F D 1 (p F = μ/ v F is the Fermi momentum) and the interlayer distance D is the smallest characteristic distance in the system. In this case, a behavior of the system is governed by only one dimensionless parameter, the ratio of characteristic Coulomb to quantum kinetic energy:

rs =

e2

εh¯ v F

2.19

≈

ε

(the maximal value rs = 2.19 is achieved when graphene layers are suspended in vacuum). The cutoff energy w, specifying half-width of a pairing region around the Fermi surface, can be estimated as a characteristic frequency of the lower branch of plasma oscillations in the bilayer and found as the first zero of denominator of (7) [14]:



w=

if rs 1, if rs ∼ 1.

8μrs , 2μ,

(8)

The cases rs 1 and rs ∼ 1 correspond to the weak and strong coupling regimes respectively. In the static limit, the long-wavelength asymptotic of the polarization operator is [31] Π(q, 0) ≈ −4N , where N = μ/2π h¯ v 2F is the density of states at the Fermi level, and (7) takes a simple form:

N V (q, 0) =

rs q/ p F + 8rs

(9)

.

The integral equations (5) along with (6) and (9) are the starting point for our further consideration. We use two approximations to solve these equations at T = 0: a constant gap approximation and a separable potential approximation.

Fig. 2. The dimensionless coupling constants at p F D 1 as functions of rs : λ0 (solid line), λ1 (dotted line), λa (dashed line), λb (dash-dotted line).

λ0 = N V 0 ( p F , p F ) and λ1 = N V 1 ( p F , p F ) are dimensionless s- and p-wave harmonics (6) of the pairing potential at the Fermi level. Using (6) and (9), the coupling constants can be calculated explicitly:

λa = λb =

rs

π rs

π



2π r s − 1 −



16rs2 − 1 arccos

−2π rs + 1 + 

16rs2 16rs2 − 1



1

arccos

(12)

,

4rs 1

4rs

.

(13)

The coupling constants λa,b as well as λ0,1 are presented in Fig. 2. At rs ∼ 1, when the pairing interaction is effectively shortrange, its p-wave harmonic tends to zero, λ0 ≈ 1/8, λ1 ≈ 1/48π rs , and thus λa and λb approach each other. In the opposite case rs 1, the pairing interaction becomes long-lange, s- and p-wave harmonics are close by the order of magnitude, and one cannot neglect difference between λa ≈ (rs /π )(−1 − ln 2rs ) and λb ≈ rs /π . In the asymptotical case of small gaps Δ± μ, (10) and (11) give the following solution for Δ± :



Δ+ = 2μ exp − Δ− =

1

Λ

+

w

μ

λb λa − (λa2

− λb2 ) A˜ −



−1 ,

(14)

Δ+ ,

(15)

where we have introduced the effective coupling constant 3. Approximation of constant gap Suppose that the gap functions Δγ γ ( p ) are constant in the pairing region with the energy half-width w, partly occupying the remote bands: Δ++ ( p ) = Δ+ Θ( w + μ − v F p ), Δ−− ( p ) = Δ− Θ( w − μ − v F p ) (Δ± are real and positive, Θ(x) is a unit step function). Replacing, as usual, the potentials V l by their values at the Fermi level, we reduce the system (5) to



Δ+ = λa A + Δ+ + λb A − Δ− , Δ− = λb A + Δ+ + λa A − Δ− , w /μ±1



A± = 0

x dx

2 (x ∓ 1)2 + (Δ± /μ)2

(10)

,

(11)

where λa = (λ0 + λ1 )/2 and λb = (λ0 − λ1 )/2 are the dimensionless intra- and interband coupling constants respectively, while

Λ=

λa − (λa2 − λb2 ) A˜ − ˜ 2− 1 − (λa2 − λ2 ) A

(16)

b

and the positive increasing function of w /μ, being an asymptotic of A − (11) at Δ± → 0 and characterizing a contribution of the ˜ − = [ w /μ − 1 − ln( w /μ)]/2. The forvalence band into pairing: A mula (14) differs from the similar formula (6) in [13] by taking into account finite range of the pairing interaction expressed in replacement of λ0 /2 by (16). Fig. 3 shows the gap Δ+ as a function of rs , obtained by numerical solving of (10), for w = 8μrs and w = 2μ, in accordance to (8). At large w the gap Δ+ depend on w exponentially, which is in agreement to (14) and in contrast to usual BCS-like one-band result [14]

1 , Δ+ = 2w exp − λa

(17)

Yu.E. Lozovik, A.A. Sokolik / Physics Letters A 374 (2009) 326–330

Fig. 3. The conduction-band gap Δ+ in the model Δ± = const, normalized on the chemical potential μ, as a function of rs at T = 0, p F D 1: at w = 8μrs (solid line) and at w = 2μ (dashed line). Dotted line: the gap parameter Δ˜ + , normalized on μ and calculated in the approximation of separable potential.

where the pairing half-width w appears only as a pre-exponential factor. The ratio Δ− /Δ+ is always smaller than unity and is as smaller, as considerably λa and λb differ, according to (10); on increase of rs , the ratio Δ− /Δ+ increases and tends to unity. Note that the approximate expressions (14)–(15) provide rather accurate estimations of Δ± , deviating from the numerical results only at Δ± > μ. According to (8), the estimations w = 8μrs and w = 2μ are valid only at rs 1 and rs ∼ 1 respectively, so the appropriate values of Δ+ are given in Fig. 3 by the solid line at small rs and by the dashed line at large enough rs . Thus, in the model used the maximum value of Δ+ is (10−7 –10−6 ) × μ, which is unessentially different from the weak-coupling one-band BCS-like result (17) [14,21,23]. However, as will be shown in the next section, going beyond the approximation Δ± = const can rise the gaps considerably owing to integration over wide region of momenta in the gap equations. 4. Approximation of separable potential In separable potential approximation we take into account a dependence of the gap functions Δ± ( p ) and potentials V 0,1 ( p , p  ) on momenta in (5). Similarly to [32], we single out a separable part of the harmonics (6) of the statically screened poten(sep) tial (9) as follows: N V l ( p , p  ) = λl v l ( p ) v l ( p  ), where v l ( p ) = V l ( p , p F )/ V l ( p F , p F ). The replacement of the original potentials (sep) V l ( p , p  ) by their separable parts V l ( p , p  ) in the gap equations (5) lead to slight underestimation of the gap Δ+ ( p F ) with (sep)

respect to exact solution of (5), since V l

( p , p  ) ≈ V l ( p , p  ) when

(sep) p or p  is close to p F and V l ( p , p  ) < V l ( p , p  ) when the both momenta p and p  are far from the Fermi surface.

With the separable potential in (5), the gap functions Δγ γ ( p )

take the form Δγ γ ( p ) = Δ˜ γ v a ( p ) + Δ˜ −γ v b ( p ), where v a,b ( p ) = [ v 0 ( p ) ± v 1 ( p )]/2, and the gap equations (5) become algebraic equations with respect to two gap parameters Δ˜ ± :



Δ˜ + = λa A Δ˜ + + λb C Δ˜ + + λa C Δ˜ − + λb B Δ˜ − , Δ˜ − = λa C Δ˜ + + λb A Δ˜ + + λa B Δ˜ − + λb C Δ˜ − ,

(18)

+ + U− , B = U+ + while the coefficients of these equations A = U aa bb bb + − − U aa , C = U ab + U ab ,

329

Fig. 4. The exponent ln(Δ˜ + /2μ) as a function of rs . Solid line: the result in the separable potential approximation, according to (22). Dotted line: the same, but with neglecting a contribution of the valence band. Dashed line: the curve ln(Δ˜ + /2μ) = −2/λ0 , corresponding to the BCS limit.

γ

Uij =

vF

∞

pF

0

v i ( p ) v j ( p ) p dp  2 ( v F p − γ μ)2 + Δ2γ γ ( p )

(19)

are functions of Δ˜ ± themselves. Since the integrals (19) converge, we do not need any momentum cutoffs. + diverges logarithIn the asymptotical case Δ˜ ± → 0 only U aa mically, whereas all the other integrals (19) converge to finite + = ln(2μ/Δ +, A = ˜ + ) + U˜ aa limits. Extracting this singularity as U aa ˜ and denoting the limits of A, ˜ B and C at Δ˜ ± → 0 ln(2μ/Δ˜ + ) + A ˜ 0 , B 0 and C 0 respectively, we obtain the asymptotical solution as A of (18):

 Δ˜ + = 2μ exp − 1 − λa ( A˜ 0 + B 0 ) − 2λb C 0        + λa2 − λb2 A˜ 0 B 0 − C 02 / λa − λa2 − λb2 B 0 , Δ˜ − = Δ˜ +

λb + (λa2 − λb2 )C 0 λa − (λa2 − λb2 ) B 0

.

(20) (21)

The results of numerical solving of (18) deviate from (20)–(21) only at rs > 2. Assuming, for simplicity, that λa ≈ λb ≈ λ0 /2 (which is well justified at strong coupling, according to Fig. 2), we get

2 Δ˜ + = 2μ exp − + A˜ 0 + B 0 + 2C 0 . λ0

(22)

This result is similar to (14) and differs from the BCS result by a presence of additional terms in the exponent; these terms are positive anywhere except a small region at rs 1. Fig. 4 shows the exponents a in the approximate formula Δ˜ + = 2μ exp(a) for the gap parameter Δ˜ + in three cases: according to the formula (22), taking into account only conduction band contributions (i.e. neglecting the integrals U i−j ), and in the BCS limit when the exponent is a = −2/λ0 . We see that, at large enough rs , the additional (with respect to the BCS result) contributions to the exponent, resulting from the valence band and from the integration on momentum in wide region in the conduction band, are numerically large, approximately equal to each other and grow linearly with rs . The curve for Δ˜ + , obtained by numerical solving of (18), is drawn as a dotted line in Fig. 3. In the separable potential approximation, Δ˜ + can reach large values comparable to μ at large rs .

330

Yu.E. Lozovik, A.A. Sokolik / Physics Letters A 374 (2009) 326–330

This result demonstrates a sharp contrast to the BCS estimation (17), although is close to the curve of Δ+ at w = 8μrs in Fig. 3, since the characteristic size of a momentum region, providing a major contribution to the gap equations and being 8p F rs , acts similarly to restricting the pairing region to an energy half-width of w ∼ 8μrs , according to (9). 5. Kosterlitz–Thouless transition temperature We have performed analysis of an instability of the unpaired state of the system in the approximation of Bethe–Salpeter equation, which will be published in detail elsewhere [27]. This analysis allows to determine the mean-field transition temperature T c , playing the role of an upper limit for a Kosterlitz–Thouless transition temperature T KT in a two-dimensional system. Our estimations show that the critical temperature to zero-temperature gap ratio T c /Δ+ ( T = 0) in the multi-band regime is the same as in the BCS theory, i.e. exp(γ )/π (γ ≈ 0.577 is the Euler constant). However, the consideration of the Kosterlitz–Thouless transition temperature deserve several remarks. This temperature is determined by a condition of vortex unpairing

πρs ( T KT ) 2

= N KT T KT ,

⎜ 0 0 0

1 0 0 1 0 0

Acknowledgements The work is supported by the Russian Foundation for Basic Research and by the Program of the Russian Academy of Sciences. A.A.S. acknowledges support given by the Dynasty Foundation and by the Russian Science Support Foundation. References

(23)

where ρs ( T ) is superfluid density (phase stiffness) at finite temperature, equal to μ/4π at T = 0 and vanishing at T c [23]. The factor N KT refers to vortex fractionalization; it is 1 for usual vortices and 4 for half-vortices [24] having twice smaller energy than the former ones. In graphene bilayer, quarter-vortices are possible as topological excitations of the form ⎛ iϕ ⎞ e 0 0 0

P = S⎝

thus gapped and soliton-like excitations, corresponding to oscillations of this phase, are possible. The self-consistent gap equations are solved in two approximations: in the approximation of constant gaps and in the approximation of separable potential. Owing to integration over wide region of momenta in the both valence and conduction bands in the gap equations, taken into account in the latter approximation, the gap values can be significantly larger than in the conventional constant gap approximation. Determination of Kosterlitz–Thouless transition temperature in the strong coupling regime should take into account a possibility of quarter-vortex unpairing and corrections for a superfluid density due to participation of the remote bands in the pairing.

0⎟ + ⎠S 0 1

in U (4) order parameter, where S is some unitary matrix, smoothly varying in space, and ϕ is an azimuthal angle measured around the quarter-vortex core. These quarter-vortices have four times smaller energy than usual vortices, as can be shown by analogy with [24], and so are the most favorable to be created; this leads to N KT = 42 = 16 for our system (which is different from N KT = 4, supposed in [23]). At weak coupling, when one-band regime effectively takes place, T KT differs insignificantly from the mean-field temperature T c . However, at strong coupling, a participation of the remote bands in the pairing should lead to corrections of an expression for ρs ( T ) as compared to one-band model; detailed analysis of superfluid properties of graphene bilayer in the multi-band regime at strong coupling and determination of corresponding T KT will be addressed to future work. 6. Conclusions We have considered the pairing in electron–hole graphene bilayer in the multi-band regime at strong coupling, when both valence and conduction bands of both graphene layers participate coherently in the pairing. The ground state of the system is described by a system of multi-band Gor’kov equations, which involve two order parameters, corresponding to valence and conduction bands, and are governed by intra- and interband coupling constants, determined by finite-range properties of the pairing interaction. A relative phase of two gaps is fixed in the ground state,

[1] J. Bardeen, L.N. Cooper, J.R. Schrieffer, Phys. Rev. 108 (1957) 1175. [2] H. Suhl, B.T. Matthias, L.R. Walker, Phys. Rev. Lett. 3 (1959) 552. [3] H. Schmidt, J.F. Zasadzinski, K.E. Gray, D.G. Hinks, Phys. Rev. Lett. 88 (2000) 127002. [4] V. Barzykin, L.P. Gorkov, JETP Lett. 88 (2008) 131. [5] M.V. Sadovskii, Phys.-Usp. 51 (2008) 1201. [6] R.D. Pisarski, D.H. Rischke, Phys. Rev. D 60 (1999) 094013. [7] T. Ohsaku, Phys. Rev. B 65 (2001) 024512; T. Ohsaku, Phys. Rev. B 66 (2002) 054518. [8] M.G. Alford, A. Schmitt, K. Rajagopal, T. Schäfer, Rev. Mod. Phys. 80 (2008) 1455. [9] K.S. Novoselov, et al., Science 306 (2004) 666; K.S. Novoselov, et al., Nature 438 (2005) 197. [10] M.I. Katsnelson, K.S. Novoselov, Solid State Commun. 143 (2007) 3; M.I. Katsnelson, K.S. Novoselov, A.K. Geim, Nature Phys. 2 (2006) 620. [11] A.K. Geim, K.S. Novoselov, Nature Mater. 6 (2007) 183; M.I. Katsnelson, Mater. Today 10 (1–2) (2007) 20; A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, A.K. Geim, Rev. Mod. Phys. 81 (2009) 109. [12] T. Ohsaku, Int. J. Mod. Phys. B 18 (2004) 1771. [13] N.B. Kopnin, E.B. Sonin, Phys. Rev. Lett. 100 (2008) 246808. [14] Yu.E. Lozovik, A.A. Sokolik, JETP Lett. 87 (2008) 55; Yu.E. Lozovik, S.P. Merkulova, A.A. Sokolik, Phys.-Usp. 51 (2008) 727. [15] Yu.E. Lozovik, V.I. Yudson, JETP Lett. 22 (1975) 274; Yu.E. Lozovik, V.I. Yudson, JETP 44 (1976) 389; Yu.E. Lozovik, V.I. Yudson, Solid State Commun. 19 (1976) 391. [16] P. Nozières, S. Schmitt-Rink, J. Low Temp. Phys. 59 (1985) 195. [17] Yu.E. Lozovik, A.A. Sokolik, J. Phys.: Conf. Ser. 129 (2008) 012003. [18] O.L. Berman, Yu.E. Lozovik, G. Gumbs, Phys. Rev. B 77 (2008) 155433. [19] H. Min, R. Bistritzer, J.-J. Su, A.H. MacDonald, Phys. Rev. B 78 (2008) 121401(R); R. Bistritzer, A.H. MacDonald, Phys. Rev. Lett. 101 (2008) 256406. [20] C.-H. Zhang, Y.N. Joglekar, Phys. Rev. B 77 (2008) 233405. [21] M.Yu. Kharitonov, K.B. Efetov, Phys. Rev. B 78 (2008) 241401(R). [22] R. Bistritzer, H. Min, J.-J. Su, A.H. MacDonald, arXiv:0810.0331v1 [cond-mat]. [23] M.Yu. Kharitonov, K.B. Efetov, arXiv:0903.4445v1 [cond-mat]. [24] I.L. Aleiner, D.E. Kharzeev, A.M. Tsvelik, Phys. Rev. B 76 (2007) 195415. [25] B. Seradjeh, H. Weber, M. Franz, Phys. Rev. Lett. 101 (2008) 264404. [26] H. Schmidt, T. Lüdtke, P. Barthold, et al., Appl. Phys. Lett. 93 (2008) 172108. [27] Yu.E. Lozovik, A.A. Sokolik, arXiv:0910.4366v1 [cond-mat]. [28] Yu.E. Lozovik, A.V. Poushnov, Phys. Lett. A 228 (1997) 399. [29] S.M. Apenko, D.A. Kirzhnits, Yu.E. Lozovik, Phys. Lett. A 92 (1982) 107. [30] S. Das Sarma, A. Madhukar, Phys. Rev. B 23 (1981) 805. [31] B. Wunsch, T. Stauber, F. Sols, F. Guinea, New J. Phys. 8 (2006) 318; E.H. Hwang, S. Das Sarma, Phys. Rev. B 75 (2007) 205418. [32] V.A. Khodel, V.V. Khodel, J.W. Clark, Nucl. Phys. A 598 (1996) 390.