Relativistic charge and energy
transport phenomena in graphene
nanostructures 

Dominik Suszalski 

Rozprawa doktorska 

Promotor: prof. dr hab. Adam Rycerz
Promotor pomocniczy: dr Grzegorz Rut 



Uniwersytet Jagielloński
Instytut Fizyki im. Mariana Smoluchowskiego
Zakład Teorii Materii Skondensowanej i Nanofizyki 


Kraków, lipiec 2020 

Wydział Fizyki, Astronomii i Informatyki Stosowanej Uniwersytet Jagielloński 
Oświadczenie 

Ja niżej podpisany Dominik Suszalski (nr indeksu: 1079343) doktorant Wydziału Fizyki, Astronomii i Informatyki Stosowanej Uniwersytetu Jagiellońskiego oświadczam, że przedłożona przeze mnie rozprawa doktorska pt. „Relativistic quasiparticles in charge and energy transport through graphene nanostructures” jest oryginalna i przedstawia wyniki badań wykonanych przeze mnie osobiście, pod kierunkiem prof. dr hab. Adama Rycerza. Pracę napisałem samodzielnie. 
Oświadczam, że moja rozprawa doktorska została opracowana zgodnie z Ustawą o prawie autorskim i prawach pokrewnych z dnia 4 lutego 1994 r. (Dziennik Ustaw 1994 nr 24 poz. 83 wraz z późniejszymi zmianami). 
Jestem świadom, że niezgodność niniejszego oświadczenia z prawdą ujawniona w dowolnym czasie, niezależnie od skutków prawnych wynikających z ww. ustawy, może spowodować unieważnienie stopnia nabytego na podstawie tej rozprawy. 
Kraków,dnia29.07.2020 ................................. 

podpis doktoranta 
Contents 

1 Introduction 11 
1.1 Aimandscopeofthisthesis......................... 11 

2 Graphene structure 13 
2.1 Crystalstructureofgraphene ........................ 13 

2.2 Electronicstructureofgraphene....................... 13 

2.3 Crystalstructureofbilayergraphene . . . . . . . . . . . . . . . . . . . . 14 
2.4 Electronicstructureofbilayergraphene................... 15 


2.5 The 2-band approximation for bilayer graphene . . . . . . . . . . . . . . 15 
2.6 ThePeierlssubstitution ........................... 16 

2.7 Localdensityofstates ............................ 16 

3 Scattering matrix formalism 18 
3.1 Theconceptofleads ............................. 18 

3.2 Scatteringmatrix............................... 18 

3.3 Transport characteristics at constant temperature . . . . . . . . . . . . . 19 
3.4 Thermoelectricproperties .......................... 20 

3.5 Adiabaticpumping .............................. 20 

4 Investigating graphene 21 
4.1 Pursuingthetrigonalwarpingstrength . . . . . . . . . . . . . . . . . . . 21 
4.1.1 Motivation ............................... 21 

4.1.2 Seebeck coefficient: the simple approach . . . . . . . . . . . . . . 21 
4.1.3 Temperaturedependentanalysis. . . . . . . . . . . . . . . . . . . 24 

4.2 Thermoelectric properties in the presence of the bandgap . . . . . . . . . 27 
4.2.1 Introduction.............................. 27 

4.2.2 Seebeck coefficient and T <BR>(α) <BR>models................. 27 

4.2.3 Thermoelectricfigureofmerit(ZT) . . . . . . . . . . . . . . . . . 30 


4.3 Conductivity of the large sample of bilayer graphene . . . . . . . . . . . . 31 
4.3.1 Introduction.............................. 31 

4.3.2 Zero-temperatureresults....................... 31 

4.3.3 Finitetemperatureeffects ...................... 34 

4.3.4 Spontaneous symmetry breaking in bilayer graphene . . . . . . . . 36 
5 Selected devices made out of graphene 37 
5.1 Aharonov-Bohmeffectwithouttwoslits . . . . . . . . . . . . . . . . . . 37 
5.1.1 Introduction.............................. 37 

5.1.2 Earlierwork.............................. 37 

5.1.3 Magnetoconductance at non-zero chemical potential . . . . . . . . 38 
5.1.4 Magnetoconductance in the absence of cylindrical symmetry . . . 38 
5.2 Mesoscopicvalleyfilter............................ 41 

5.2.1 Conceptofvalleytronics ....................... 41 

5.2.2 Thesetup ............................... 41 

5.2.3 Roleofelectrostaticpotential .................... 42 

5.2.4 Roleofthemassterm ........................ 42 

5.2.5 Operatingthevalleyfilter ...................... 45 



5.3 Adiabatic quantum pumping with mechanical kinks in graphene . . . . . 49 
5.3.1 Aimofthesection .......................... 49 

5.3.2 Simplekinkmodel .......................... 50 

5.4 Dynamicsinthesimplekinkmodel ..................... 50 

5.4.1 TheSSHkinkmodel ......................... 54 

5.4.2 ResultsfortheSSHkinkmodel ................... 55 

6 Conclusions 60 








Publications constituting the thesis 
1. 
D. Suszalski, G. Rut and A. Rycerz, Lifshitz transition and thermoelectric properties of bilayer graphene, Physical Review B 97, 125403 (2018). 

2. 
D. Suszalski, G. Rut and A. Rycerz, Thermoelectric properties of gapped bilayer graphene, Journal of Physics: Condensed Matter 31, 415501 (2019). 

3. 
D. Suszalski, G. Rut and A. Rycerz, Mesoscopic valley filter in graphene Corbino disk containing a p–n junction, Journal of Physics: Materials 3, 015006 (2019). 

4. 
D. Suszalski, G. Rut and A. Rycerz, 	Conductivity scaling and the effects of symmetry-breaking terms in bilayer graphene Hamiltonian, Physical Review B 101, 125425 (2020). 

5. 
R. Adam and D. Suszalski, Graphene disk in a solenoid magnetic potential: Aharonov-Bohm effect without a two-slit-like setup, Physical Review B 101, 245429 (2020). 

6. 
D. Suszalski and A. Rycerz, Adiabatic quantum pumping in buckled graphene nanoribbon driven by a kink, http://arxiv.org/abs/2002.08507, accepted to Acta Physica Polonica B. 


Supplement: 
7. D. Suszalski and A. Rycerz, Adiabatic pumping driven by moving kink and quantum standard ampere in buckled graphene nanoribbon, http://arxiv.org/abs/2007.11145, (unpublished). 
Acknowledgments 
I would like to express my deep gratitude to my supervisor, Prof. Adam Rycerz for his time and guidance. I am also very grateful for his critical reading of this Thesis. 
My deep thanks are extended to my supporting supervisor, Dr. Grzegorz Rut, for his presence and all the friendly advice given. 
I thank Prof. Anton Akhmerov for hospitality in his research group and introduction to Kwant package. 
Many thanks go to Dr. Maciej Fidrysiak, Dr. Danuta Goc-Jagło, Dr. Andrzej Kapanowski,Dr. AndrzejKądzielawa,Prof. JózefSpałek,KrzysztofNowakowski,Wojciech Tarnowski, Franciszek Sobczuk, Dr Bas Nijholt, Dr Kim Pöyhönen, Dr Piotr Rożek, Dr Daniel Varjas, Dr Adriaan Vuik for fruitful conversations and all help received. 
Special thanks go to my wife for her love, support and ongoing presence, and my parents for all the effort they put in raising me. 
Last but not least, I would like to thank the Lord for life, love and all the beautiful world to investigate. 
The work was supported by the National Science Centre of Poland (NCN) via Grant No. 2014/14/E/ST3/00256. The financial support from dotation KNOW from Krakowskie Konsorcjum “Materia-EnergiaPrzyszłość” im. Mariana Smoluchowskiego is acknowledged. 
Abstract 

This thesis is concentrated on quantum transport properties of mono-and bilayer graphene systems in ballistic regime. The Landauer–Büttiker formalism is utilised to analyse transport properties of ultraclean samples. The aims of this work are twofold. The first one is concerned in analysing the properties of graphene itself, both mono-and bilayer. The second one is concentrated on novel electronic devices that might be constructed thanks to extraordinary properties of graphene. The particular systems are chosen weighting the theoretical interest, experimental feasibility and possible future applications. 
The first part addresses three main questions. The first one How one may reliably measure the trigonal warping strength in bilayer graphene? is replied with vast analysis of thermoelectric properties of bilayer graphene. Two experimentally accessible ways has been proposed. The first one requires measurement of carrier concentration corresponding to the secondary maximum of the Seebeck coefficient, while the other requires measurement of temperature at which the Seebeck coefficient reaches its global maximum on the doping – temperature plane. 
The second question: How appearance of a bandgap will modify the thermoelectric properties? is followed by analysis of two system geometries. The traditional lead-sample-lead geometry and the sharp potential step geometry allowing to pinpoint the crucial properties of the system in idealised case. The thermoelectric properties occur to be defined dominantly by the presence of the bandgap. The special attention is put to unusual behaviour of figure of merit (ZT). 
The third question: What is a conductivity of large piece of bilayer graphene? considers the impact of non-leading terms in effective bilayer graphene Hamiltonian on the conductanceofthelargesystematzerodoping. Thethreedifferentbehavioursofconductivity (diverging, levelling off or decaying) with growing system length are recognized and interpreted. The analysis includes the finite temperature effects and comparison with available experimental results. 
The second part is also naturally divided into three parts. The first one concentrates on possibility of pinpointing the Aharonov-Bohm effect without two slit setup in the Corbino disk made out of graphene. It occurred that pseudodiffusive transport in bilayer grapheneisaffectedbymagneticfluxpiercingthecentralelectrodeandsufficientlystrong to allow its experimental measurement. The robustness of the effect is analysed in terms of non-zero doping and breaking of rotational symmetry of the system. 
The second system is a proposition of a mesoscopic valley filter – the first device in valleytronicsbeingthemodificationofelectronicsthattakesadvantageofvalleypolarised currents in graphene. As proposition does not require use of experimentally problematic single atom precise sample forming or application of strain-induced pseudomagnetic fields, we believe it allows experimentalconstructionof the first valley filter – a milestone in development of valleytronics. 
The last system consists of buckled nanoribbon with mechanical kink generator. Such a system has been analysed before, by other authors, and proved to be controllable via gentle shaking of one of its ends, taking advantage of unusual negative radiation pressure effect. Our central question: Does buckled graphene nanoribbon might be used as an electric pump? is addressed and analysed in terms of SSH model. The surprising result show it should be possible not only to construct such a pump, but also to use it as a standard current device (ampere prototype). 
Streszczenie 

Dyssertacja jest skoncentrowana na kwantowych cechach przewodnictwa w ballistycznych układach grafenu i dwuwarstwy grafenowej. Do analizy przewodnictwa wykorzystano formalizm Landauera–Büttikera. Cele niniejszej pracy dzielą się na dwie ogólne grupy. Pierwsza grupa koncentruje się na analizie własności monowarstwy i dwuwarstwy grafenowej. Druga natomiast rozważa własności konkretnych układów wykonanych z grafenu. Wybór układów jest podyktowany interesującą fizyką, realizowalnością oraz możliwymi, przyszłymi zastosowaniami. 
Pierwsza część pracy zadaje trzy główne pytania. Po pierwsze, Jak wiarygodnie zmierzyć wartość tunelowania skośnego w dwuwarstwie? pytanie to jest analizowane w świetle własności termoelektrycznych dwuwarstwy grafenowej. Dwie dostępne eksperymentalnieprocedurysązaproponowane. Pierwszawymagapomiarukoncentracjinośników ładunku odpowiadającej pobocznemu maksimum współczynnika Seebecka, podczas gdy drugawymagapomiarutemperaturyodpowiadającejnajwyższejwartościwspółczynnika Seebecka w przestrzeni parametrów domieszkowanie – temperatura. 
Drugie pytanie Jak obecność przerwy energetycznej wpływa na własności termoelektryczne? jest analizowane dla dwóch różnych geometrii układu. Tradycyjnej geometrii złącze-próbka-złącze oraz geometrii z gwałtownym skokiem potencjału, która umożliwia analizę kluczowych własności dwuwarstwy w sytuacji uproszczonej. Analizowane własności termoelektryczne okazują się być w dominującym stopniu zdeterminowane przez obecność przerwy energetycznej. Szczególny nacisk jest położony na analizę współczynnika wydajności (ZT). 
Trzeciepytanie Jaka jest przewodność dużego fragmentu dwuwarstwy grafenowej? dotyczy kwestii wpływu niewiodących wyrazów Hamiltonianu na przewodnictwo dużych układów niedomieszkowanej dwuwarstwy grafenowej. Trzy istotnie różne zachowania przewodnictwa (rozbieżne, stabilizujące się oraz zanikające) zostały zaobserwowane przy wzroście długości układu w zależności od obecności/braku poszczególnych wyrazów. Analiza uwzględnia również efekty skończonej temperatury oraz porównanie z wynikami eksperymentalnymi. 
Drugaczęśćpracyrównieżnaturalniedzielisięnatrzyczęści. Pierwszakoncentrujesie na możliwości obserwacji efektu Aharonova-Bohma bez tradycyjnej geometrii dwuszczelinowej w dysku Corbino wykonanym z grafenu. Okazuje się, że strumień pola magnetycznego przechodzący przez centralną elektrodę wpływa na pseudodyfuzyjny transfer ładunku w stropniu umożliwiającym eksperymentalną obserwację. Trwałość efektu jest analizowana ze względu zmienne domieszkowanie i złamanie symetrii obrotowej układu. 
Drugi z rozważanych systemów jest propozycją konstrukcji mezoskopowego filtra dolinowego – pierwszego z urządzeń postulowanych w ramach dolinotroniki, będącej wersją elektroniki wykorzystującą możliwość dolinowego polaryzowania prądu w grafenie. Jako że proponowany układ nie wymaga przygotowania z precyzją sięgającą pojedynczych atomów, ani wykorzystania pochodzących z naprężeń pól pseudomagnetycznych, mamy nadzieję, że pozwoli on na realistyczną eksperymentalnie konstrukcję pierwszego filtru dolinowego – kamienia milowego rozwoju dolinotroniki. 
Ostatni z rozważanych systemów składa się ze ściśniętego paska grafenu połączonego z generatorem mechanicznych kinków. Taki układ był niedawno analizowany przez innych autorów ze względu na możliwość kontroli ruchu solitonów (kinków) poprzez delikatne potrząsanie jednym z końców paska. Mechanizm kontroli wykorzystuje niecodzienne zjawisko ujemnego ciśnienia promieniowania. Nasze pytanie Czy ściśnięty pasek grafenowy może posłużyć jako pompa elektronowa? jest analizowane w ramach modelu SSH (Su-Schrieffer-Heeger). Dość nieoczekiwanie wyniki wskazują na możliwość wykorzystania takiego układu jako wzorca ampera. 
1 Introduction 

Mastering new materials has been one of key ingredients driving human development. Starting from ancient times and mastering the stone and brass, through mastering the ironandsteeluptonowadays. Understandingtheelectronbehaviourgunfiredelectricity revolution. Last century resulted in mastering semiconductors, liquid crystals, optical fibres and many other materials. The possibilities opened by those breakthroughs do not need to be enlisted for a modern man. To say that we achieved a lot would be an understatement, yet it seems we are still at the beginning of the journey. 
Many researchers struggle in order to understand high temperature superconductors, that are already used in medical tomography. The fascinating world of topological insulators and Majorana edge modes is now opening, holding the hope for realisation of scalable quantum computers. Single atom manipulation via Scanning Tunnelling Microscope is already achievable in almost every nanoscale laboratory. There’s Plenty of Room at the Bottom[1] and there are probably wonders that not a man has ever dreamed about. Down there, among many others, there is graphene – a unique material holding very special place. 
Graphene is a stable, self-standing, two-dimensional material of thickness of a single carbon atom. Such materials were long thought to be non-existing and even impossible to exist due to reasons connected with mechanical stability of in-plane vibartions [2]. In finite temperatures its mechanical properties are exceptional making it one of the strongest material ever tested [3]. Graphene is a great and promising thermal conduc<BR>tor [4]. It has unusually high carrier mobility making it a great current conductor [5]. Its phenomenal electronic feasibility has drawn interest of many researchers investigating graphene devices in areas like unconventional superconductivity [6], spintronics [7], valleytronics or physics of Quantum Hall Effect [8]. Despite purely laboratory charac<BR>ter of first graphene production method (2004) it has already found first applications and is available in commercial products including cars [9], earphones [10], helmets [11], tires [12], lubricants [13], conducting ink [14], conducting packing foil for electronics [15] and photodetectors [16]. Graphene allowed also to relax operating conditions of Quan<BR>tum Hall resistance standards from Tmax ≈ 1.3 K, Bmin ≈ 10 T, and Imax ≈ 30 µA to Tmax ≈ 10 K, Bmin ≈ 3.5 T, and Imax ≈ 0.5 mA significantly reducing their size, cost and complexity [17]. 
1.1 Aim and scope of this thesis 
Among different interests and hopes concerning graphene this thesis is devoted to its electronic properties in ballistic regime. The Landauer–Büttiker formalism is used throughout the thesis. Despite all the research done, there are still some questions concerning properties of graphene itself. This thesis addresses some of them 
• 
What is a value of trigonal warping strength in bilayer graphene? 

• 
How does it affect thermoelectric properties? 

• 
How thermoelectric properties are going to change if the bandgap appear? 

• 
What is the conductivity of infinite sample of ultraclean bilayer graphene? 


The first question is one of few elementary questions about electronic structure of bilayer graphene that is still open. The second one is a consequence of approach taken to the first one. The third one is a natural extension and a complementary research. The fourth question, probably the most complex one, is devoted to the analysis of effective low energy Hamiltonian and significance of terms appearing in it for conductivity of large samples. 
To understand graphene well may yet be only a temporary goal. The real value of graphene lies in possibilities opened by its understanding and mastering of its fabrication. Among enormous number of graphene applications considered in contemporary discussion this thesis concentrates on three special cases 
• 
realisation of Aharonov-Bohm effect in Corbino geometry (without two slit setup), 

• 
construction of a mesoscopic valley filter, 

• 
construction of an electric pump in graphene nanoribbon. 


The first task is an interesting possibility to realise Aharonov-Bohm effect without standard two slit setup. Such a realisation would also allow to show the magnetoconductance oscillations in graphene Corbino disk that are not experimentally available in constant magnetic field situation [18]. The second task is a new proposition of construc<BR>tion of mesoscopic valley filter -the device that might be use to valley polarise current flowing through graphene. The basic idea of valley polarised current has given rise to the valleytronics -the modern electronics, yet to be implemented due to experimental difficulties. The mesoscopic nature of proposed filter should help to overcome some of them. The third task has brought interesting result, showing possibility of constructing mechanically driven electric pump that transports a quantum value of charge per cycle. Such a pump might be then used to construct a device serving as a current standard. The author believes that two last system might lead to construction of devices that push further limits of the human development. 
2 Graphene structure 
2.1 Crystal structure of graphene 
Grapheneisoneofcrystallinecarbonallotropes. Itis2-dimensionalhoneycombcrystal of carbon atoms. The primitive lattice vectors are 
√ 
a1 = a, 3a, a2 =(a, 0), (1)
22 
with lattice constant of graphene a =0.246 nm [19]. Within the unit cell there are two atoms located at positions 
a
dA = (0, 0), dB = √ , 0 , (2)
3 
leading to presence of two triangular sublattices, often referred to as A and B sublattices (or atoms). 
The reciprocal lattice of a triangular lattice is also a triangular lattice leading to hexagon shape of Wigner-Seitz cell being the first Brillouin zone. Its primitive vectors 
are defined as  
�  �  �  �  
b1 =  2π a , 2π √ 3a  ,  b2 =  2π a , − 2π √ 3a  .  (3)  

It is important to note that out of six corners of the Brillouin three may be connected by inverse lattice vectors. Thus there exists two inequivalent corners called K and K’ points. Hereinafter, the following convention is chosen 
4π 4π
K = , 0 and K = − , 0 . (4)
3a 3a 
2.2 Electronic structure of graphene 
The band structure of graphene crosses the zero-energy level at K and K’ points [20]. In the vicinity of each of them (called also valleys) there is a conical structure that resembles dispersion relation of relativistic, massless particles. Around the K and K’ points the effective Hamiltonian may be written as: 
Hξ = v0p 
· (ξσx ,σy)+ V (x
)+ M(x
)σz , (5) 
where ξ = 1(−1) corresponds to K(K’) valley, v0 ≈ 106 m/s is the Fermi velocity, p 
= −in' is a 2-dimensional momentum operator, V (x
) is a potential that has the same value on both sublattices and M(x
) is a staggered potential. The above Hamiltonian is equivalent to the 2-dimensional Dirac Hamiltonian, where V (x
) plays the role of external potential and M(x
) plays a role of a mass term. It is also worth to note that second nearest neighbour tunnelling within a graphene lattice does not introduce terms linear 

Figure 1: Structure of tight-binding parameters in bilayer graphene. The matrix elements γ3 and γ4 are displayed only once for clarity of the figure. The B(A) sublattice in upper(lower) lattice are marked with black dots. They are vertically aligned and connected with the biggest interlayer matrix element γ1. The other sublattices are marked by coloured dots. 
in momentum to the Hamiltonian. The peculiar Hamiltonian structure is the source of bizarre behaviour of electrons in graphene, such as spinor transformation rules for the wavefunction, Klein tunnelling or lack of backscattering for normal incidence. 
2.3 Crystal structure of bilayer graphene 
Bilayer graphene consists of two graphene layers based on top of one another. In order to minimize energy of such system the second lattices needs to be shifted with respect to the first one forming so-called Bernard stacking [21]. In this stacking atoms from one sublattice in the lower sheet are positioned directly under the atoms from the second sublattice in the sheet above. Such pairs of vertically aligned atoms are called dimers, while the other atoms will be referred to as non-dimer ones. Both real and reciprocal primitive cells are the same as in the case of graphene. The K(K’) points are inherited from their definition in monolayer graphene. Infinite self-standing bilayer graphene is an energetically stable system. 
2.4 Electronic structure of bilayer graphene 
Theelectronicpropertiesofbilayergraphenearemorecomplicatedthanthoseofsingle layer graphene. This is not surprising as in addition to intra-layer tunnellings one need to include interlayer ones. The tight binding Hamiltonian including three dominant interlayer tunnellings takes a form: 
⎞
⎛ 
Hb = 
⎜⎜⎝ 

δAB/2 − U/2 v0πγ1 −v4π† v0π† −δAB/2 − U/2 −v4π† v3π 
γ1 −v4π −δAB/2+ U/2 v0π† −v4πv3π† v0πδAB/2+ U/2 
⎟⎟⎠

, 

(6) 

where π = e−iθ(ξpx +ipy), π† = e−iθ(ξpx −ipy), θ is the angle between armchair direction and the x-axis, ξ equals +1(−1) for K(K’) valley, v3 = v0γ3/γ0, v4 = v0γ4/γ0, U is the electrostatic bias between layers and δAB is the staggered potential describing irreducible part of spontaneous bandgap in bilayer graphene. One might also add additional potential (Δ ) introducing energy difference between dimer and non-dimer sites. Such a term introduces electron-hole asymmetry to the system, similarly to the γ4 tunnelling. As γ4 is present in model considered Δ may be skipped as redundant clarifying interpretation of γ4. 
Thevaluesof γ0 =3.16 eV and γ1 =0.38 eVarewidelyaccepted[22]. Thevalueof γ4 ≈ 
0.15 eV is known with reasonable precision, while γ3 remains somehow a mystery with values reported from 0.1 eV to 0.4 eV. Experiments report no presence of spontaneous bandgap in bilayer graphene based on semiconductor. The value of interlayer energy offset U is defined by external out of plane electric field. Provided values are used throughout the thesis unless stated otherwise. 
The difference in band structure of bilayer and monolayer graphene comes predominantly from γ1 term. Its presence changes the graphene band structure around K(K’) points significantly. It is no longer linear and become quadratic in momentum. Moreover, bands split leading to appearance of low energy conductance band (that touches highest valence band) and high energy conductance band (and corresponding lower valence band). It may be also interpreted as giving mass to electron pseudoparticles. 
2.5 The 2-band approximation for bilayer graphene 
In the low energy limit the electronic properties of the bilayer graphene origins almost entirely from the lowest conductance band and the highest valence band. It is thus possible to catch them in the 2-band approximation [21] 
Hˆ2 = hˆ0 + hˆ3 + hˆ4 + hˆU + hˆAB, (7) 
1 0(π†)2 
hˆ0 = − , (8)
2mπ2 0 
ˆ0 πv3a 0(π†)2 
h3 = v3 −√ , (9)
π† 0 π2 0
43n 2vv4 π†π 0 
hˆ4 = , (10)
γ1 0 ππ† 
U 10 2v2 π†π 0 
hˆU = −− , (11)
20 −1 γ12 0 −ππ† 
δAB 10 
hˆAB = , (12)
20 −1 
with Hˆ2 being the 2-band approximate Hamiltonian. The hˆ0 term dominates in low energy limit, thus other terms might be considered a correction to it. It is worth noting that second part of hˆ3 term is only a renormalization of effective electron mass and the unique role of γ3 comes from the first part. 
2.6 The Peierls substitution 
The inclusion of magnetic field is somehow unique among many different physical effects. Usually inclusion of additional effects correspond to addition of a new term to effective Hamiltonian. In the case of magnetic field it is rather modification of existing terms. The exact formula is given by the Peierls substitution 
p 
→ p 
+ eA, (13) 
where e is electron charge and A is the vector potential corresponding to magnetic field considered B = '× A. The vector potential is a gauge field giving one some freedom in parametrization of the problem. Throughout this thesis, the symmetric gauges will be chosen. 
2.7 Local density of states 
Density of states is a useful quantity giving insight into the quantum system considered. Itdescribesthenumberofquantumstatespresentingivenenergyrange. Formally, it is defined as 
 
ρ(E)=δ(E − En), (14) 
n 
where δ(x) is Dirac delta, En is energy of state n and sum is performed over all eigenstates of Hamiltonian. Analysis of density of states reveals presence or lack of valence/conduction band and the bandgap, that is crucial for electronic properties of the system. In some cases additional, useful information may be obtained from its local version. It is unsurprisingly called the local density of states and in the lattice models is defined as 
ρloc(Rj,E)= |ψ(j)|2δ(E − En), (15)
n n 
where Rj is a position of site j and ψn (j) is the value of n-th Hamiltonian eigenstate wavefunction at site j. Local density of states gives information where are localized states corresponding to given energy. For practical reasons, usually in order to obtain better looking plots, both density of states and local density of states is smeared by substituting 
δ(x) → 1 �, (16)
πx2 + �2 where � defines (small) energy scale over which the densities are smeared. 
3 Scattering matrix formalism 
3.1 The concept of leads 
One of the main difficulties in the domain of quantum transport is connecting the macroscopic circuits (batteries, amperometers, voltmeters or probes in general) with the investigated nanostructure. Due to the quantum nature of such systems the measurement rules known from macroscopic world do not longer apply. The theoretical object allowing modelling of macroscopic probes connected to the true quantum system is a lead. The lead is a thermodynamic reservoir of electrons occupying delocalized states. Connecting two or more of such leads to the quantum nanostructure allows one to characterise transport properties of the structure. The electrons coming from an incoming lead may be reflected back from a nanostructure or transmitted to a second lead giving rise to the flow of the electric current through the nanostructure. If a number of open channels -quantum delocalized states allowing to deliver or pick up electrons -in the leads is much bigger than the number of open channels in the nanostructure the current flowing through the system does not usually depend on the microscopic details of the leads and may be solely predicted by the details of nanostructure and thermodynamic properties of the leads (the temperature and chemical potential for electrons). 
It is striking that transport properties of the nanostructure are described most conveniently solely in terms of states in the leads being at once independent on the specific choice of the leads1 . 
3.2 Scattering matrix 
The quantum system may be then described by so called scattering matrix (Sˆ(E)) connecting the incoming electron states in the leads (Ψin(E)) with outgoing electron states in the leads (Ψout(E)): 
Ψout(E)= Sˆ(E)Ψin(E), (17) 
or explicitly specifying states from left (L) and right (R) leads (the dependence on chemical potential is omitted for clarity): 
ΨL,out sˆLL sˆLR ΨL,in
= . (18)
ΨR,out sˆRL sˆRR ΨR,in 
Traditionally, the blocks building the scattering matrix are renamed into transmission (t and t ) and reflection (r and r ) parts of scattering matrix: 
ˆ
ˆˆrˆt
sLL sLR
ˆ
S = ≡ ˆ. (19)
sˆRL sˆRR trˆ
1The transmission amplitudes and form of scattering matrix blocks may (and will) of course depend on the exact choice of the leads, yet the set of transmission eigenvalues will approximately not. 
Of great theoretical interest are probabilities of electron incoming in channel n being transmitted through (or reflected by) the nanostructure. These probabilities are given by diagonal values formed by transmission and reflection matrices: 
 n 
Tn =tˆ†tˆnn , (20)
 n 
Rn =rˆ† rˆnn . (21) 
3.3 Transport characteristics at constant temperature 
The most straightforward transport characteristic of a quantum nanosystem is conductance (G) being the inverse of the electric resistance. To obtain a formula for the conductance of a system it is convenient to start with formula for the current flowing through the system [23]: 
⎧  ∞⎨ 
dkx
I =2sevx(kx)fL(E) (22)
⎩2π 
n 0 
⎫ 
 0 ⎬
dkx 
+vx(kx)[Rn(E)fL(E) + (1 − Rn(E))fR(E)] (23)
2π ⎭ 
−∞ 
 ∞

2se 
= (1 − Rn(E)) [fL(E) − fR(E)] dE (24)
h 
n 0 
 ∞ 
2se 
= Tn(E)[fL(E) − fR(E)] dE (25)
h
0 n 
2se2 ≈ UTn(E) , (26)
h 
n 
obtaining desire formula for the system conductance: 
G(E)=2sG0 Tn(E), (27) 
n 
where G0 = e2/h is conductance quantum, fL(R)(E) is the Fermi-Dirac distribution of electrons in the left (right) lead, U is the (small) voltage and 2s stands for spin 
t tt 
degeneracy. In the derivation, identities 1 − Rn =1 − Rn = Tn following 
nnn from the unitarity of the scattering matrix have been used. In the final approximation the assumption that transmission coefficients does not significantly depend on chemical potential in the scales defined by voltage U and temperature kbT has been made. 
3.4 Thermoelectric properties 
Despite a fact that electric current is the current that drives electronic and informatic revolutions the thermoelectric properties of materials are also of significant importance. The electronic part of thermal current (and thermal properties of a material) may be easily represented in Landauer–Büttiker formalism. Starting from formulas for electrical and thermal currents [24]: 
ge 
I = −T (E)[fL(E) − fR(E)]dE, 
h 
g

IQ = T (E)[fL(E) − fR(E)](E − µ)dE, 
h 
where g is the degeneracy factor, µ is the (average) chemical potential in the leads and T (E) is shorthand notation for Tn(E). In the linear response regime (infinitesi
n 
mal potential and temperature differences) the Seebeck (S) coefficient and the thermal conductance Kel are given by: 
t 
V 

ΔT 

L1 
, (28)
S = − 

= 

I=0 eT L0 
= 

L0L2 − L21 
I=0 TL0 
IQ 
, 

(29)

Kel = 
ΔT 

where  
Ln  =  g h  T (E)  − ∂f ∂E  (E − µ)n ,  (30)  

with f being the Fermi Dirac distribution with thermodynamic parameters (chemical potential and temperature) taken as mean values between the leads. 
3.5 Adiabatic pumping 
The analysis of transport properties of a nanostructure up to now assumed that the nanostructure itself does not change during the time of experiment. When it does, this process may generate the current between the leads even in the absence of voltage and temperature difference between the leads. Assuming that typical time corresponding to change of the scatterer is much larger than time of flight of electron through the scatterer one may write adiabatic formula for a charge transferred [23] 
ie dSˆSˆ†
ΔQα = − Trα dt (31)
2π dt 
where ΔQα is charge transferred to a lead α and trace is taken over elements corresponding to states originating from the lead α. 
4 Investigating graphene 
4.1 Pursuing the trigonal warping strength 
4.1.1 Motivation 
Trigonal warping is the common name for γ3 tunnelling -arguably the second most important interlayer tunnelling in bilayer graphene. This term breaks the rotational symmetry of the approximate band structure in the vicinity of the K and K points. Moreover it leads to the splitting of each of two Dirac cones into 4 new cones (see Fig. 2 (a)). A single cone out of 4 new cones is centred exactly at the K (K ) point, while 3 remaining cones are centred at momentum p = γ1v3/v2 at three different directions. The central cone is approximately circular, while the satellite cones are approximately elliptical with ratio of major-an minor-semi axes equal 3. 
The splitting of Dirac cones leads to enhancement of minimal conductivity of bilayer graphene by (approximately) factor of 3 [25] being the possible explanation for exper<BR>imentally measured value [26]. The splitting of the Fermi surface does not survive to higher chemical potentials. This has obvious consequence of appearance of the Lifshitz phase transition reconnecting the splitted Fermi surface. The Lifshitz transition appears 
at Lifshitz energy  
EL =  γ1 4  γ3 γ0  2  (32)  

and is accompanied by the presence of van Hove singularity in bulk density of states (see Fig. 2 (b)) [27]. 
Despite being profound in its consequences the Lifshitz transition has till now been avoiding precise experimental determination. The available experimental results cover the range of EL ∼ 0.1 − 1 meV [29, 30, 22]. The motivation of this work was to propose an experimentally convenient way to measure trigonal warping strength with satisfying precision. 
4.1.2 Seebeck coefficient: the simple approach 
For the clarity of reasoning and in order to concentrate on the influence of trigonal warping strength the simplest possible model has been chosen. Explicitly, our model is equivalent to setting v4 = U = δAB =0 in Hamiltonian (6). 
The Lifshitz transition is accompanied by van Hove singularity in density of states. One may thus expect abrupt behaviour of conductance of the system in the vicinity of the Lifshitz energy. Unfortunately there are too small to smoke gun the transition in experiment. The well-known Mott formula for metals2 [31] states that Seebeck coeffi<BR>cient is proportional to logarithmic derivative of transmission with respect to chemical potential. One may thus expect that due to low conductance at the vicinity of Lifshitz 
2S =(π2/3)e−1k2 T [∂lnT (E)/∂E]E=µ
B

Figure 2: (a) The Fermi surface for E =0.5EL, E = EL, E =1.5EL depicted by grey solid lines, red solid lines and blue dashed lines respectively. (b) Density of states (purple solid lines) with the case γ3 =0 (grey solid line) and linear approximation for γ3 
=0 (black dash dotted line, see eq. (33)). Adapted from [28]. 
transition and the presence of derivative in the Mott formula the abrupt changes originating from van Hove singularity will be magnified in dependence of Seebeck coefficient on chemical potential. 
Theanalysisstartswithsimplified2-bandHamiltonianforbilayergraphene(7). Treat<BR>ing each Dirac cone separately and expanding the Hamiltonian to linear terms in momentum leads to simple formula for bulk density of states: 
|E|
ρ(E) ≈ ρ0 , (33)
EL where ρ0 is the bulk density of states of bilayer graphene in the absence of trigonal warping. One may thus expect linear dependence of conductance of the system as the function of chemical potential. Mode matching numerical analysis confirms these suspicion for chemical potential E � EL/2 (see Fig. 3). Unfortunately for a system with linear dependence of conductance on chemical potential the Seebeck coefficient is uniquely defined by temperature and holds no information about the slope of conductance vs chemical potential dependence. Moreover one may not simplyapply Mott formula to numerical analysis due to significant contribution from Fabry-Pérot resonances. One is forced to use non-simplified equation (28). Fortunatelyenough,thefullnumericalcalculationofSeebeckcoefficientshowsanomaly -the additional maximum at the vicinity of Lifshitz energy (EL) for appropriate temperature range (see Fig. 4). When the temperature is too large, it average out all details connected to Lifshitz transition. When it is to small, the Fabry-Pérot oscillations give 

Figure 3: The zero temperature conductivity of bilayer graphene for 3 distinct values of trigonal warping (marked on picture). The vertical dashed lines corresponds to Lifshitz energies for considered values of trigonal warping (compare eq. 32). The inset is a zoom-in close to charge neutrality point. The horizontal lines corresponds to conductivities σ =2 σMLG = (8/π)e2/h and σ =6 σMLG. Reprinted from [28]. 

Figure 4: Left: The Seebeck coefficient for γ3 =0.3 eV for 4 different temperatures (see the picture). The dashed vertical line marks the Lifshitz energy. Right: The Lorentz number with respect to the value provided by Wiedemann-Franz law. The temperatures are the same as in the left panel. Inset shows the crossover temperature corresponding to disappearance of second maximum of Seebeck coefficient (triangles) and second minimum of Lorentz number (triangles) with respect to Lifshitz energy. Reprintedfrom [28]. 
higher contribution than the anomaly. 
Similar anomaly is present in the Lorentz number being the ratio of electronic part of thermal conductance to the conductance L = Kel/T G though this time it represents itself as an additional minimum (see Fig. 4) below the value given by Wiedemann-Franz law (LWF =(π2/3)(kB/e)2), obeyed by the standard Fermi gas (with constant DOS). 
MeasuringthustheSeebeckcoefficientortheLorentznumberasafunctionofchemical potential allows to find a value of Lifshitz energy. One may be worried about possibility of measuring the electronic part of thermal conductance, yet such worries would be misplaced. The electronic part of thermal conductance may be measured as a difference betweenthetotalthermalconductanceandthephononthermalconductance. Measuring the phonon thermal conductance is possible by electrostatically opening the energy gap, reducing the electronic contribution of thermal conductance to negligible value. 
The real inconvenience of such measurement comes from the need to measure values considered as a function of chemical potential, while most of experiments provides them rather as the function of carrier concentration. This difficulty might however be overcomed. 
4.1.3 Temperature dependent analysis 
In order to avoid measurement of chemical potential one may measure global properties of the dependence of Seebeck coefficient on the chemical potential. The procedure proposed is to find (for a given temperature) the global maximum of Seebeck coefficient Figure 5: The carrier concentration corresponding to the maximal Seebeck coefficient 

as a function of temperature for three different values of trigonal warping. 
Reprinted from [28]. 
as a function of chemical potential (controlled by the bias gate voltage). Than for the bias gate voltage corresponding to the maximal value of Seebeck coefficient one should measure the carrier concentration. Such a measurement should be repeated for different values of temperature. From measured temperature dependence of carrier concentration corresponding to the maximal Seebeck coefficient one may obtain the value of trigonal warping strength (compare Fig. 5). It should be stressed, that despite of finding the maximum of Seebeck coefficient as a function of chemical potential, the chemical potential need not to be measured. The corresponding bias gate potential is fully capable of controlling the measurement process. 
During the investigation of system properties another method has been found. It is based on the observation that the maximal value of Seebeck coefficient (as the function of chemical potential) grows up to some temperature and than starts to fall down. It occurred that the temperature corresponding to the maximal value of Seebeck coefficient dependslinearlyontheLifschitzenergy. Itisthusenoughtomeasuretemperaturecorresponding to global maximum of the Seebeck coefficient as the function of both chemical potential and temperature (compare Fig. 6). Once again, the chemical potential does not need to be measure, even once. 

Figure 6: Temperature corresponding to global maximum of the Seebeck coefficient (on the chemical potential vs temperature plane) as a function of Lifshitz energy. Adapted from [28]. 
4.2 Thermoelectric properties in the presence of the bandgap 
4.2.1 Introduction 
Previous analysis has been restricted to situation, where there is no bandgap in the system. It is valuable to include in the analysis also the case with the bandgap present. Such gap may be easily realised experimentally by applying out of plane electric field.3 In order to do so, one needs to include U term from Hamiltonian (6) on top of all parameters considered till now. In the absence of trigonal warping, the bandgap is equal 
s 
to Δ= |U|t⊥/U2 + t2 . The presence of trigonal warping leads to minor correction, 
⊥
that usually might be ignored. 
As the presence of the bandgap suppress the minimal conductance via evanescent modes and the bandgap influence is of primary interest one may get rid of the Fabry-Pérot oscillations by considering the sharp potential step geometry instead of usual lead-sample-lead geometry. Such sharp potential step geometry might be understood as a geometry joining two regions (weakly-and heavily-doped) at an abrupt line. Unfortunately despite producing nice looking, qualitatively correct predictions it also overestimates the conductance by the approximate factor of 2 (due to lack of backscattering at the second interface). All the analysis is than performed for two system geometries the abrupt potential step one and the traditional lead-sample-lead one. The first gives clear view, while the second gives results easier to realize in typical experiments. 
4.2.2 Seebeck coefficient and T (α) models 
As one might expect opening the bandgap leads to strong suppression of the conductance in the bandgap regime followed by rapid increase just after the bottom of the conduction band (see Fig. 7). Following the Mott formula one would expect the appearance of the maximum of Seebeck coefficient at the bottom of conduction band. That is not the case for finite temperatures. Instead it appears at the vicinity of the Goldsmid-Sharp value |S|GS =Δ/(2eT ) derived for wide semiconductors. 
max 
In order to understand that behaviour one might apply the family of simple models quantifying the transmission of the system 
�n n 
E − 11
δ 2 Δ+ δE + 2 Δ if α =0 
T (α) 
= C(Δ) × nnα−1 (34)
Θ |E|− 1 Δ |E|− 1 Δ if α> 0
22 
where C(Δ) is parameter controlling the value of transmission, δ(x) is the Dirac delta function and Θ(x) is Heaviside step function. 
Applying the above family of models to the calculation of Seebeck coefficient leads to theasymptotic(u =Δ/(2kBT ) » 1)formulaforthemaximumof the Seebeckcoefficient 
kB 11 
|S|(α) ≈ u − ln(u)+ α − + O(u −1) , (35)
max e 22 
3Some experiments report appearance of spontaneous bandgap in the suspended graphene of the order 
of 2.5 meV <BR>[32, <BR>33]. <BR>Such <BR>a <BR>bandgap <BR>is <BR>much <BR>smaller <BR>than <BR>the <BR>gaps <BR>analysed <BR>in <BR>this <BR>work <BR>and <BR>is <BR>not <BR>
found in the graphene placed on semiconductor, thus it is of little importance for this analysis. 

Figure 7: Conductance and Seebeck coefficient as a function of chemical potential for different values of the bandgap (from 0 meV to 20 meV every 5 meV). Solid lines correspond to potential step geometry, while dashed lines correspond to lead-sample-lead geometry. The inset shows zoom-in of conductance around the bottom of the conductance band (the curves where appropriately shifted). Reprinted from [34]. 

Figure 8: (a) Maximal value of the Seebeck coefficient as a function of Δ/(2kBT ) for sharp potential step geometry. (b,c) Deviation from Goldsmid-Sharp value in the presence (absence) of trigonal warping. Points show numerical calculations for T =1K (open circles), T =5 K (full circles) and T = 10 K (triangles). Lines shows predictions of T (α) models for α =0 (blue dashed), α =1 (red solid) and α =2 (magenta dash dotted). Reprinted from [34]. 
whichclearlyreproducesGoldsmid-Sharpvalueintheleadingterm. Moreoveritpredicts this maximum to appear at the chemical potential 
kBT 2Δ
|S|
µ≈ ln (36)
max 2 kBT 
that is much lower value than the one expected by naive implementation of Mott formula (µexpected =Δ/2). This lead to conclusion that the properties of Seebeck coefficient (at sufficiently low temperatures) are predominantly determined by the presence of bandgap and the details of band structure or transmission give just minor corrections. Indeed, the comparison of numerical calculations with predictions of T (α) models gives great agreement (see Fig. 8). One should though stress that influence of trigonal warping strength is still visible in low temperature data (see Fig. 8 (b,c)). 

Figure 9: The figure of merit (ZT) as a function of chemical potential for different values 
ofthegap. Thetemperatureissetto 5 K.The colorand solid/dashed encoding 
of lines is the same as in Fig. 7. Reprinted from [34]. 
4.2.3 Thermoelectric figure of merit (ZT) 
Having analysed the Seebeck coefficient one may now wonder how the thermoelectric figure of merit (ZT) will behave in such a system. Usually, the maximum of ZT corresponds to the maximum of the Seebeck coefficient due to its square dependence on it. However, in the system considered that is no longer the truth. The reason for that lies in a strong suppression of the figure of merit in the gap regime due to almost negligible conductance. One might have expected that the ratio of electronic to thermal electron conductance would still remain almost constant and thus the suppression of conductance should not diminish the figure of merit. Such reasoning would yet be misplaced, as the figure of merit does depend on total thermal conductivity and the phonon thermal conductance is almost unaffected by electric field opening the gap. 
In gapped bilayer graphene the maximum of figure of merit corresponds to the bottom of the conduction band (and top of the valence band), see Fig. 9. This is the result of interplay betweenenergies withdisappearing conductance and regions with disappearing Seebeck coefficient. Again the results for finite sample system are qualitatively the same as the results for potential step geometry, being just diminished by some numerical factor. 
Once again, the T α models (see Eq.(34)) are capable of both describing the figure of merit and predicting the maximum of it with satisfying precision (especially for α =1, see Fig. 10). Just like in the case of the Seebeck coefficient one may conclude that the properties of the figure of merit are determined by the presence of the bandgap, while the other properties of band structure play at most secondary role. 
Phononsplaydominantroleinthermalconductanceduringthedeterminationoffigure of merit for bilayer graphene. Among different degrees of freedom, the dominant contribution to phonon thermal conductance comes from out of plane modes (ZA phonons). Thus suppressing them by placing bilayer graphene on rigid substrate may open practical possibility to enhance the maximal value of the figure of merit, approximated by a factor of 3. 
4.3 Conductivity of the large sample of bilayer graphene 
4.3.1 Introduction 
Bilayer graphene is a crystal characterized by several non-equivalent parameters. Researchers included different sets of them, appropriately to the goal chosen. Here the following question is addressed 
What is the D.C. conductance of large piece of ultra-clean bilayer graphene and which microscopic parameters are crucial for appropriate description of it? 
In order to pursue crucial model parameters the symmetry driven approach, ennobled in long history of physics, has been taken. Among different parameters in Eq.(6) only U is omitted, as the interest lies in intrinsic properties of bilayer graphene, not the ones originating from external fields. Parameters γ0 and γ1 describes the nearest neighbour hoppings within the layers and between them respectively. There are thus unavoidable and form a minimal model of bilayer graphene. The trigonal warping strength γ3 introduces the breaking of rotational symmetry, while γ4 breaks the electron-hole symmetry. The δAB term represents the irreducible part of bandgap that may not be closed by external electric field. Such a gap may spontaneously appear due to electron-electron interaction [35], that are otherwise ignored in this work. 
For the simplicity and clarity of discussion unique, non-zero values has been chosen for each parameter according to reported experimental results [21]. Namely γ3 =0.3 eV, γ4 =0.15 eV and δAB =1.5 meV has been chosen. 
4.3.2 Zero-temperature results 
The zero temperature results in the presence of a gap are not surprising. The bilayer graphene becomes an insulator, while for short samples (L< 100 l⊥ with l⊥ = 
√ 
3aγ0/(2γ1) ≈ 1.77 nm being the effective length appearing in bilayer graphene transport analysis that originates from coupling between the layers) it preserves some transmission via evanescent modes (see Fig. 11 (c,d)). 
The physics becomes more rich in gapless scenarios. In absence of both skew interlayer tunnellings the conductance of the system is just twice the conductance of single 

|S|
(µmax) and figure of merit (µZT ) as a function of Goldsmid-Sharp value. 
max
Point symbol encoding is the same as in Fig. 8. Lines correspond to T (α) model predictions for α =1 (solid) and α =2 (dash dotted). Reprinted from [34]. 

layer graphene, despite presence of the coupling (it is a non-trivial result, see Ref.[37]). Introduction of trigonal warping (and thus breaking of the rotational symmetry) slowly increases the conductance up to the value of 3 σ0, with σ0 being the conductance of two decoupled graphene layers. Introduction of γ4 term in place of trigonal warping leads in opposite to slow, power law decrease of conductance with the system length (σ ∝ L−2.0). One would expect the situation with presence of both trigonal warping and γ4 term to be a competition between two opposite forces. Instead, presence of both skew tunnelling terms leads to unbounded increase of conductance with the system length (see Fig. 11 (a)). The reason of this approximately linear increase is the presence of propagating modes at zero energy (or any other energy). The propagating modes originate from elevation of three of secondary Dirac cones. The secondary cones are introduced by trigonal warping, yet elevated by γ4 term. Thus only presence of both terms may introduce conducting states at zero energy to the bilayer graphene. This situation where interplay of two factors leads to appearance of qualitatively new phenomenon is a beautiful example of the rich structure of quantum mechanics and physics in general. 
It is worth to note, that Fano factor does not significantly differ from the pseudodiffusive value (F =1/3), except for the regimes of vanishing conductance, when unsurprisingly it approaches value of 1. 
4.3.3 Finite temperature effects 
At the beginning of discussion of influence of temperature on conductance properties of the system one should note that interaction induced gap depends on temperature and disappears for a critical value of temperature. In order to reproduce reported temperature dependence [38] the gap will be than modelled as 
⎧  �  
δAB(T ) = δAB(0) ⎨ ⎩ tanh 0  1.74  TC T  − 1  if T TC if T > TC ,  (37)  
with TC = 12 K and δAB(0) = 1.5 meV.  

One of the most striking features of finite temperature conductance is that it almost does not depend on γ4 term, that becomes irrelevant. One may thus conclude that electron-hole symmetry does not play important role in finite temperature conductance of bilayer graphene (see Fig. 12). 
Unsurprisingly, T> 0 leads to increase of conductance with higher temperature leading to higher increase of conductance. The increase seams unlimited with system length for any finite temperature considered. The presence of thermically reduced gap leads to diminished increase of conductance, with effect disappearing altogether with gap closing. 
The most interesting part of conductance is probably the regime of sufficiently short, gapless systems where the effects of finite temperature are almost negligible. In this regime the quantum properties of conductance via evanescent states are the most important and dominate the finite temperature effects. 
The electronic noise for finite temperature is dominated by Nyquist-Johnson noise, making the Fano factor irrelevant. 

4.3.4 Spontaneous symmetry breaking in bilayer graphene 
Another interesting property of conductance of bilayer graphene is the fact that it may be used as an example of spontaneous symmetry breaking system with three possible scenarios. Considering σ(L) being the order parameter the spontaneous symmetry breaking is gun-smoked by non-commuting order of limits 
lim L→∞  [...] σ = σ0,  (38)  
lim d→∞ lim L→∞ lim δAB→0 σ = ∞,  (39)  
lim d→∞ lim δAB→0 lim L→∞ σ = lim δAB→0 σ = 0,  (40)  

where limit of infinite interlayer distance d (d →∞) correspond to simultaneous limits γ3 → 0, γ4 → 0 and γ1 → 0 (the last limit being irrelevant for the final result). One may than conclude that in bilayer graphene breaking of the rotational symmetry, breaking of electron-hole symmetry and breaking of sublattice equivalence (due to γ1 =0) may appear spontaneously. 
5 Selected devices made out of graphene 
5.1 Aharonov-Bohm effect without two slits 
5.1.1 Introduction 
Aharonov-Bohm effect is a bizarre quantum effect that has no analogue in classical mechanics. It is unique in the sense that it allows to influence time evolution of charged particles by the magnetic field that is never met by particles wavefunctions. Despite many experimental realizations [39, 40] its existence in orthodox version presented here has been long questioned and the effect was attributed to residual Lorentz force (originating from e.g. leakage fields) acting on the electron wavefunction [41, 42]. Nowadays the discussion is settled, by beautiful experimental realization [43] exploiting shielding of the magnetic field with superconductor. 
In vast majority of experimental realizations of Aharonov-Bohm effect the trajectories encircling the magnetic field flux are considered as they provide most straightforward interpretations. Inthisworkanalternativeapproachistaken. Insteadofelectroncurrent encircling the magnetic field, the electron wavefunctions themselves encircle magnetic field flux at each moment of experiment and electron current flows radially according to rotational symmetry of the system. 
5.1.2 Earlier work 
The functional form of transmission coefficient for Corbino disk in uniform magnetic field at charge neutrality point is long since known [44] 
Tj =1 , (41)
cosh2 [ln(Ro/Ri)(j + φd/φ0)] 
where j corresponds to total angular momentum (Jz = nj, j = ±1/2, ±3/2, ...), φd = π (Ro 2 − Ri 2) B is the magnetic flux piercing the disk, φ0 = 2(h/e) ln(Ro/Ri) being period of conductance oscillations, B standing for magnetic field, Ro(Ri) for outer (inner) disk radius. Unfortunately, the uniform magnetic field setup occurred to be unsuitable to measure oscillations of magnetoconductance, leaving space for more subtle approach involving Aharonov-Bohm effect. 
The picture is even more simple in the case of Aharonov-Bohm effect, when all the magnetic field flows through the hole in the inner lead and does not flow through the disk itself. In such a situation above formula holds, but with simplified substitution [45] 
φd ≡ φi, (42) φ0 ≡ φAB = h/e, (43) 
where φi stands for total magnetic flux piercing the inner lead. Nonetheless, the analysis for non-zero chemical potentials were missing. 
5.1.3 Magnetoconductance at non-zero chemical potential 
The magnetoconductance of the system has been derived using mode matching analysis. The closed form solution for the transmission coefficients has been derived [18] and used to analyse the magnetoconductance. The total conductance of the system grows almost linearly starting from charge neutrality conductance via evanescent modes (see Fig. 13 (a)). The magnetoconductance is perfectly periodic with a period given by Aharonov-Bohm flux (φAB) (see Fig. 13 (b)). Notably enough, the magnetoconductance oscillations are quite large (of the order of 0.15 φAB for chosen Ro/Ri =5). Moreover the oscillations appear at almost all energies (see Fig. 13 (d)) with contraposition to the uniform magnetic field case, where they appear only in the vicinity of Landau levels [44]. The further analysis shows that similarly high magnetoconductance oscillations appear also for different ratios of disk radii (see Fig. 13 (c)), yet they might be shifted toward higher chemical potential. 
5.1.4 Magnetoconductance in the absence of cylindrical symmetry 
In experimental reality the perfect cylindrical symmetry might be hard to achieve. Almost all imperfections leads to mode mixing and many to local charge density fluctuations (i.e. p-n puddles). The natural question arises, whether results of the previous section survive in lower symmetry reality. In order to investigate this issue lets consider application of in plane electric field leading to addition of linear potential 
U0r 
U(r, ϕ)= − sin(ϕ), (44)
Ro where U0 stands for maximal value of potential at the edge of the system. Such a term introduces mode mixing to the system removing possibility of finding compact closed formsolutionfortransmissioncoefficients. Thenumericalanalysishasbeenthenapplied, following the standard procedure. The results are compactly summarised in Fig. 14. As expected, presence of in-plane electric field enhanced the zero-energy transmission (via conducting states formed in heavily doped regions). The change in total conductance is almost negligible for the higher chemical potentials, as one may consider the uniformly doped sample as being ’averaged’ version of a sample with linear potential. The interest lies though in the magnetoconductance. At the charge neutrality point the amplitude of magnetoconductance oscillations is strongly diminished. This may be understood, provided the additional conducting states are roughly localized either in upper or lower half circle. Such states does not encircle the magnetic field flux and there is no reason to believe they should be significantly influenced by its change. Fortunately, the amplitude of magnetoconductance oscillations is preserved even for U0 values being higher than p-n puddle potential variations (being of the order of 10 meV [47]) provided than one not restrict the analysis to charge neutrality point (see Fig. 14). This proves both the robustness of results obtained and importance of non-zero chemical potential analysis. 

Figure 13: (a) The conductance of the system as a function of doping for the ratio Ro/Ri =5. Lines depict φi =0 (blue solid) and φi = φAB/2 (red dashed). Blackdashdottedlineintheinsetistheaverageover φi. (b)Themagnetoconductancefordifferentdoping. Dashedlinemarkspseudodiffusiveconductance Gdiff =2g0/ ln(Ro/Ri) with g0 =4e2/h. (c) Magnitude of magnetoconductance oscillations as a function of doping for different radii ratio. The lines for Ro/Ri =5 (Ro/Ri = 10) are shifted by 0.25 g0 (0.5 g0). (d) 2-dimensional map of nodal lines marking conductance unaffected by magnetic field on the radii ratio – doping plane. ΔG = G(φAB/2) − G(0). Reprinted from [46]. 

5.2 Mesoscopic valley filter 
5.2.1 Concept of valleytronics 
The presence of two Dirac cones (valleys) in a dispersion relation of graphene does usually manifest itself as an additional (to spin) two-fold degeneracy of electron states. This follows from equivalence of valleys -there are exchanged under combined action of time reversal and charge conjugation operators. As long as this symmetry holds and there is no coupling between layers the degeneracy of states is be preserved. Breaking those symmetries may lead however to more complex behaviour. 
Valleytronic is probably most far going idea of utilizing the valley inequivalence [48]. It postulates to physically separate the currents originating from both graphene valleys in analogy to spintronics postulating separation of currents with respect to spin polarization. Such valley polarized current might be than used to store information and give rise to new type of electronics — the valleytronics itself. 
The valley polarization is defined as 
TrTξ=1 − TrTξ=−1
P = , (45)
TrTξ=1 + TrTξ=−1 
with Tξ being the transmission matrix for valley ξ (ξ =1, −1 corresponds to K, K valley). 
The first valleytronic device, crucial for the whole idea, would be then the valley filter (or polariser) of the current. There are many propositions of constructing such devices, all employing either system preparation with single-atom precision, or presence of strained induced pseudomagnetic fields, both being problematic for experimental realisation. In opposition to them this work tries to propose construction of mesoscopic valley filter -the valley filter that operation will not depend on single-atom details, but rather on mesoscopic properties without experimental difficulties of strain induced pseudomagnetic fields. 
5.2.2 The setup 
In order to achieve the goal lets concentrate on effective graphene Hamiltonian with two additional potentials 
(ξπxσx + πyσy) ψ(r, φ)=[E −V(r, φ) −M(r, φ)σz] ψ(r, φ) (46) 
with V(r, φ) being the electrostatic potential and M(r, φ) being the position dependent mass term. The system in Corbino geometry has been chosen due to clarity of picture following from the lack of edge states. The electrostatic potential is assumed to origin from in-plane electric field leading to form of 
V(r, φ)= −eEr sin(φ − φV ), (47) 
where E stands for electric field and φV defines orientation of direction of electric field. For further use define V = eERo with Ri, Ro being respectively the internal and external Corbino disk radii. The ratio Ro/Ri =4 has been chosen throughout all section. For the sake of simplicity, the step profile of mass term has been incorporated 
M(r, φ)= MΘ(φ − φM )Θ(π + φM − φ), (48) 
with M quantifying the mass potential strength and φM defining the mass orientation. Such a mass term may originate from chemical functionalization [49] or the adsorption of hexagonal boron nitride [50]. Throughout the whole paper φV = π/2 and φM =0 has been preserved. The filter is designed in such a way, that after system fabrication it might be fully controlled by gate voltages adjusting the in-plane electric field and chemical potential. It is worth noting that fully functioning filters has been modelled also for continuous change in spatial mass distribution following very similar procedure to the one described here. 
5.2.3 Role of electrostatic potential 
The analysis of system behaviour starts with the case of zero mass (M =0). In the case of lack of external electric field and zero chemical potential the system is in pseudodiffusive limit and the conductance of the system slightly oscillates around the approximate value of 6g0, with g0 = e2/h (see Fig. 16). For small electric field the behaviour of the system does not differ significantly up to some threshold value of magnetic field, where the conductance drops to the value of 4g0. The stronger the electric field, the lower the threshold magnetic field. Above that magnetic field the current flows along the p-n interface through so-called snake states, with single snake state per valley explaining the additional factor of 2 (on top on spin degeneracy). The snake states are most intuitively described via the semi-classical picture involving Klein tunnelling. In that picture travelling electron is subjected to classical cyclotron trajectories and consecutive changes into and from a hole at the moment of traversing the p-n interface (see Fig. 15). Despite easily understandable picture the quantum mechanical states them<BR>selves does not exhibit ’snake shape’ structure and usually are evenly distributed around the interface preserving the translational symmetry along the (straight) p-n interface. Nonetheless, the formation of snake states, being just special kind of Quantum Hall states, is well understood and will be of primary interest of this work. 
The crucial property of this snake states used in construction of valley filter is their appearance on only single side of the system for a given voltage polarization between the leads (see Fig. 17). It is enforced by their Quantum Hall nature, that allows them to traverse the interface in exactly one direction. At the end of this paragraph I would like to note that up to know the currents originating from both valleys behave exactly in the same way. In order to construct a working valley filter one needs additional tool. 
5.2.4 Role of the mass term 
The primary role of the mass term is to break the valley symmetry. It is possible, because mass term breaks the symplectic symmetry that gives rise to effective timereversal symmetry in each of the valleys. Combined with magnetic field that breaks 




the true time-reversal symmetry (involved in symmetry between valleys) it should be capable of allowing the existence of valley filter. 
Lets start with the case of zero magnetic and electric fields. Oddly enough, in configuration considered the valley currents are already well separated (see Fig. 18). One may relate it to forming of the zero doping edge states formation, with mass step playing the role of effective edge of the system. Unlike in true edge case, in mass step setup there is no reason to enforce sharpness of the step. In addition to the cases considered here, the systems with regularised step profile has been considered, all of them leading to same qualitatively results. 
5.2.5 Operating the valley filter 
Having described influence of both electric field and mass term it is time to combine theminavalleyfilter. Thespatialorientationofpotentialsisthesameasintheexamples above. The electric potential V is set to 1 meV and magnetic field is set to large value of 
s 
B =1 T. Large value means here that the magnetic length lB = n/eB and cyclotron radius are significantly smaller than distance between leads. Note however that such value describes moderate magnetic fields and is easily achievable in laboratory with use of permanent magnets. The main result of this part of work – the working diagram of the filter – is presented at Fig. 19. 

The diagram may be easily understood in idealized case (B →∞, ERi « E, M) depicted on Fig. 20. The analysis is based on assumption that effective potentials experiencedbycurrentsoriginatingfrombothvalleysmightbeapproximatedbyeffective potentials for lower and upper component of pseudospin wavefunction (see Fig. 21). Such an approximation is valid for the case of strong electric potential [52] and seems to survive to the range of parameters considered here. 
Having made the assumption about effective potentials one should only add that the eacheffectivep-njunctionwithinsinglevalleywillprovideoneconductingstate. Therole of the chemical potential is restricted to shifting the p-n junctions. As long as effective p-n junction crosses the internal lead the conducting state provided by it increases the conductance by g0. The situation is especially clear for the M =0 case, where the conductance is equal to 2g0 as long as p-n junction crosses the internal lead and 0 when it does not. The appearance of small mass leads in the first place to the splitting of the p-n junction as long as its original position stay within the upper halfcircle (being 


the region of M =0). As long as exactly one effective p-n junction crosses the internal lead the conduction is equal g0 and all the current is fully valley polarized. In the high mass limit the upper half circle effective potential is dominated by the mass effective potential (but of opposite signs for different valleys). The appearance of the effective p-n junction will thus depend on the polarization of lower half circle (ruled by the chemical potential). If the polarizations are the same there is no p-n junction and no current, if there are opposite there is the single conducting state. Following the fact of dominant contribution of mass term being opposite for different valleys, there is always just single conducting state and the current is always fully polarized. 
In the case of finite fields there are two main corrections. The first one is coming from finite width of the conducting states, leading to increase of effective radius of the inner lead by approximately lB. The second one is connected to finite value of V , which leads to formation of additional, not fully developed conducting states at the edges of P = −1 domain. Those additional states may be easily pinpointed by more detailed analysis of exact form of effective potentials (see Fig. 21). As they brought no additional insight into the physics of the system nor influence the interior of well developed domains of fully polarized current they detailed analysis will be omitted here. 
Having understood the physics of working valley filter two possible ways of switching the polarization of the current may be proposed. The first one assumes changing the chemical potential (via the bias gate potential) at given value of in-plane electric field, while the second assumes changing the in-plane electric field at constant and close to zero chemical potential. Both propositions require only changing the gate potentials, while keeping the mass term and magnetic field constant. 
The author believe that proposed here scheme of valley filter is experimentally realizable and may lead to breakthrough in realization of valleytronic devices. 
5.3 Adiabatic quantum pumping with mechanical kinks in graphene 
5.3.1 Aim of the section 
Recently, very interesting study of buckled graphene nanoribbon as a mechanical realisation of φ4 model has appeared [53]. The formation, reflection and annihilation of solitons (kinks and antikinks) has been predicted. Moreover, it has been shown that one may control the kink movement by periodically shaking one of the edges of the nanoribbon. Such shaking produces phonons that may either push the kink away from this end either pull it toward the source of phonons (being the realization of radiation within the model) resulting in straightforward realization of striking negative radiation pressure effect. 
Having such mechanical system described one may ask about characterization of its electronic properties, crowned with question 
May such a system be used as an electric pump? 
What really is a pump? Pump is a device working in a periodic cycle that produces the net flow. Usually, we encounter pumps that enforce flow of liquid or gas using electrical or combustion engines (with human heart being noticeable exception). Here, electric pump means the pump that produces flow of electrons (the electric current) out of mechanical motion of the nanoribbon (ignoring the source of mechanical force). It is worth to note, that usually electricity is produced using Faraday’s law or by chemical reaction. Here the interest lies in none of them. Electric pump understood in this way has been already proposed taking advantage of laser irradiation [54], strain fields [55], tunable magnetic barriers [56], Landau quantization [57], or sliding the Moiré pattern in twisted bilayer graphene [58]. 
Despitedifferencesindetails, thegeneralsystemgeometrywillbethesamethroughout the whole section. The nanoribbon considered is long (along y-dimension) and narrow (along x-dimension). Near both ends of the nanoribbon there are 2 electrodes connected to the opposite edges of the nanoribbon (see Fig. 22). The leads are made out of flat graphene. The energy offset in the leads is assumed to be U∞ = −0.5 t0 throughout all section. 
5.3.2 Simple kink model 
Lets consider simple model of kink, by defining elevation of atoms as 
  
y − y0 πx 
z = H tanh sin2, (49)
λW
here W isthe nanoribbonwidth, y0 isa position ofa center ofthe kink, H is the maximal elevation of atoms in graphene buckled nanoribbon in the absence of kink and λ is the parameter describing kink profile. The λ =3a has been set throughout the work in order to resemble the kink profile obtained from molecular dynamics simulations [53]. 
Having defined atom positions the standard modification of tunnellings is applied to tight binding model of electrons in graphene [60] 
δdij
tij = −t0 exp −β, (50)
d0 
where tij is the tunnelling between the atoms i and j, t0 is the nearest neighbour tunnelling value for the relaxed graphene, β =3 is the dimensionless electron-phonon cou
√ 
pling, δdij is the elongation of the bond length and d0 = a/ 3 is the equilibrium bond length. 
Suchasystem obviously lacks relaxation, as itmayonlyextend interlayerbondlength. It may be then at most the first step in this journey, yet it is capable of catching few key concepts. 
5.4 Dynamics in the simple kink model 
We start the analysis of the system dynamics by calculating the conductance of the nanoribbon in the absence of the kink, in the buckled nanoribbon (being the case of 

√ 
y0 →∞) and in the presence of the kink at the center of the system. The results are depicted in the Fig. 23. The conducting nanoribbon (W = 10a) has been chosen. As it may be easily seen the conductance of buckled nanoribbon does not significantly differ from the case of flat nanoribbon with exception of the charge neutrality point. There is a wideregion(µ0 ∈ [−0.2 t0, 0.2 t0])where conductanceisclose to 2e2/h. The 2 degeneracy factor origins from the spin degeneracy. Outside of that region the conductance becomes heavily noisy, thus the analysis will be restricted to the region itself. The behaviour of the system becomes very different in the presence of the kink. The nanoribbon stops to conduct, with the exception of two very narrow resonances, that disappear, when the kink is shifted from the center of the nanoribbon (y0 = L/2 with L being the nanoribbon length). 
Described behaviour of the system is very promising. The kink seems to play the role of impenetrable barrier for electrons. One may thus hope, that moving it through the system will sweep the electrons from the ribbon to the (second) lead. The working mechanism would be then very simple and similar to the working mechanism of the snow plow. 
The calculations seem to confirm the thesis above (see Fig. 24). For purpose of calculating the charge transferred the kink has been moved from effective −∞ position tothe effective +∞. Thechargehasbeencalculatedwithuseoftheformula(31). Asone may see the pumping efficiency is approximately linear with the chemical potential and grow with kink height. The results may be understood in terms of mechanism proposed. In the regime considered the band structure of unbuckled nanoribbon consists of two states with linear dispersion relation E = ±nvF ky. Thus, the total pumpable charge present in the nanoribbon at given chemical potential is 
Qtot |µ0| Leff |µ0|
≈ = 35.3(51)
e nvF πt0 where Leff = L − Wlead. As one may see, the moderate kink height H =2a already corresponds to almost maximal theoretical pumping efficiency. It suggests that buckling of the nanoribbon has not significantly changed corresponding density of electrons. It is worth to note that system has Z2 symmetry corresponding to the reflection of z-axis. Due to that the scattering matrix does not change under z →−z transformation. Thus full pumping cycle, corresponding to transport of two kinks through the system, that restores the original mechanical state of the system corresponds to two identical cycles of its electronic properties. Having that in mind, the term cycle will correspond to the full mechanical cycle of the pump, while all the charge will be referenced for half-cycle or in other words per kink. Summing up the results, buckling nanoribbon has little influence on its electronic properties, the kink may play a role of impenetrable barrier for electrons and it may sweep them through the system. The problem posses the Z2 symmetry corresponding to halving the electronic cycle length with respect to mechanical cycle length. Now it is time to verify, if those findings appear in more realistic model of graphene nanoribbon. 


5.4.1 The SSH kink model 
In order to validate results obtained with simple kink model the Su-Shrieffer-Heeger (SSH) Hamiltonian has been adopted [60]. It takes a form 
δdij ††
HSSH = − t0 exp −βci,scj,s + cj,sci,sd0
<ij>,s 
1 
+ Kd (dij − d0)2 
2 
<ij> 
1 n2 (52) 
+ Kθ θ�(j) − θ0
2 
j {�(j)}
⎛⎞ 
+ Vδ ⎝2π − θ�(j)⎠, 
j {�(j)} 
with a constrain 
(dij − d0)=0. (53) 
<ij> 
The modified tunnellings are the same as in the equation 50, Kd quantizes the energy corresponding to change of the bond length, Kθ parametrizes the resistance to angle bending, Vδ describes the influence of out of plane deformations and θ0 = 120◦ . The sum over set of angles {�(j)} is sum of 3 angles with the common vertex (j). Its value is equal to full angle (2π) in flat graphene. 
The minimization procedure requires simultaneous minimization of both electron and phonon parts of the Hamiltonian. The minimum is obtained by standard iterative procedure. Within the model considered the configuration of atoms and bond lengths depends on the chemical potential. However for the range of chemical potentials analysed here the difference is negligible and will be omitted. The results are provided for µ =0 being the charge neutrality point with number of electrons equal to number of atoms. Two 
√ 
systems are considered, the metallic armchair nanoribbon (W = 10 a, L =90 3 a) and 
√ 
for insulating armchair nanoribbon (W = 11 a, L =90 3 a). 
5.4.2 Results for the SSH kink model 
The case of metallic armchair is quite similar to the one described previously. The charge pumped per kink grows approximately linearly with the chemical potential and for W /W =0.90 is close to the total charge available in the conduction band (see Fig. 25). The hypothesis about sweeping the electrons by kink is supported by analysis of local density of states. For illustrative purposes, we take µ =0.04 t0, the local density of states is calculated according to formula (15) with smearing (16), � =5 · 10−3 t0. The results are presented on the figure 26. As one may see, there is no electron density in the vicinity of the kink. The kink forms an impenetrable barrier for electrons, what explains the pumping mechanism. 
The case of insulating armchair is more interesting. For sufficiently strong buckling, already of the order of W =0.96 W , the charge pumped per kink levels at the value of 2e for significant range of chemical potential (see Fig. 25). This is the pumped charge quantization. In order to understand its origin once again the local density of states will be proved useful. As one may see on Figure 26 the local density of states is concentrated in the vicinity of the kink itself. It represents electron bound state that accompanies the mechanical kink. The single bound state (though with a twofold spin degeneracy) located in the bangdap is easily visible in the density of states, see figure 27 d). Such bound state is absent in the metallic nanoribbon, as there is no bandgap and all states are delocalized. 
In the case of insulating armchair the pumping mechanism is thus very different. The kink and accompanying electron bound state are formed at the beginning of the nanoribbon. As the nanoribbon is insulating the bound state is initially unoccupied. As kink adiabatically traverses the system, the bound state is following it. Passing the first electrode the bound state is filled by electron (2 electrons counting the spin degeneracy). Later the connection to the first electrode is broken by insulating nanoribbon. At the end of the system the kink is anihilated and there is no place for the electron occupying till now the bound state. The electron(s) is thus forced to leave the system by second electrode located near the system end. In opposition to the case of metallic armchair the kink does not sweep electrons in front of itself, it rather loads them on itself and carry them between the electrodes. 


� =5 · 10−4 t0. The electron phonon coupling β =3 was present during the bond length optimization for buckling W /W =0.90. Reprinted from [61]. 

The important note is needed to be made. The density of states corresponding to systems optimized with and without electron phonon coupling (corresponding to β =3 or β =0) are quite similar, yet without the electron phonon coupling the bound state localized crucial for quantization of charge pumped does not appear (see Fig. 27). Thus in order to properly describe pumping in the insulating nanoribbon one has to include electron density impact on the bond lengths. 
Concluding results of this section, the metallic armchair nanoribbon behaves similarly to the system described by simple kink model, efficiently sweeping electron density on front of itself. On the other hand the description of insulating armchair nanoribbon is moreinvolvedandrequiresfullSSHmodel. Thechargepumpedperkinkintheinsulating armchair is quantized, opening the possibility of using it as a quantum standard of ampere. 
6 Conclusions 
The transport properties of both mono-and bilayer graphene has been studied. The special attention has been put on the analysis of thermoelectric properties and conductivity scaling for bilayer graphene systems. Moreover, three graphene devices has been proposed and analysed. The purposes of those devices are: visualisation of Aharonov-Bohm effect, valley polarization of electric current, and serving as a current standard. 
Our initial motivation of analysing the thermoelectric properties of bilayer graphene was to find a suitable way to determine the value of trigonal warping strength. This task results with three propositions: (i) to measure the anomaly of thermoelectric properties in sub-Kelvin temperatures, (ii) to measure the carrier concentration corresponding to the maximal Seebeck coefficient, and (iii) to measure temperature corresponding to global maximum of the Seebeck coefficient on the doping-temperature plane. The last two propositions do not need to be performed at sub-Kelwin temperatures, the approximate temperature of 1 K should be sufficient. It is worth to stress that proposed measurements does not require neither strong magnetic fields nor measurement of conductivity scaling with the systems length. All of them might be performed at a single sample. 
The thermoelectric properties of the bilayer graphene has been analysed with respect to varying doping, temperature and external electric field opening the bandgap. The Seebeck coefficient, Lorentz number and figure of merit has been studied numerically. Observed behaviour of discussed quantities has been rationalised in term of simple effective models. The violation of Wiedemann-Franz law in the pseudodiffusive regime has been found and reported. Among different results, the dominant role of bandgap in forming the thermoelectric characteristic should be marked. 
The influence of trigonal warping strength γ3, γ4 term and the bandgap on conductivity of long sample at zero doping point has been studied. The four different behaviour has been identified. The semiconductor behaviour characterises all systems with the bandgap. For the systems without a gap three regimes appear, the well known pseudodiffusive regime for γ3 = γ4 =0, the asymptotic pseudodiffusive regime γ3 =0,γ4 =0 and divergent pseudodiffusive regime for γ3 =0 and γ4 =0. The available experimental data allow to conclude that the conductivity of bilayer graphene may not be properly described without γ3. The role of γ4 is revealed only in samples much longer than the ones experimentally achievable so far. 
Use of graphene Corbino disk to measure the Aharonov-Bohm effect has been proposed. ThispropositionavoidsuseoftwoslitsetuptypicalformagneticAharonov-Bohm effect. Instead it presents that properties of the wavefunction propagating away from central solenoid in any direction on xy plane might be significantly affected by magnetic flux flowing through the center of the system. As in standard Aharonov-Bohm setup the wavefunction is affected, despite never seeing the magnetic field itself. The amplitude of the effect has been investigated with respect to system geometry, doping, magnetic field and in-plane electric field being a measure of system imperfections. The effect has been found to be robust and possibly measurable over wide range of parameters considered. 
Taking advantage of the presence of two Dirac cones in graphene dispersion relation the construction of mesoscopic valley filter has been proposed. The construction has been based on functionalization of graphene Corbino disk in a preparation phase, constantmagneticfield, andgatesprovidingcontroloverthedopingandthein-planeelectric field. The system has been designed to be controllable and switchable between current polarisation solely by means of gate potentials. Almost perfect filtering efficiency exceeding 99% has been observed even in the scenarios with lowered assumptions about system manufacturing precision. The device is believed to be scalable and robust against long range manufacturing imperfections. 
Turning recently considered mechanical kink representation of 1 − D field theory in grapheneintoelectricpumpservingasacurrentstandarddevicehasbeenproposed. The construction is based on clamped and buckled graphene nanoribbon and the control is performed by mechanical shaking of one of the system edges. The mechanical solitons are found to be efficient pumps pushing out almost all available electrons in their way. Thespecialcaseofsolitonisaccompaniedbytheelectronboundstatethatmightbeused to deliver quantum charge per pump cycle, potentially serving as the current standard device. 
References 
[1] Miniaturization. Reinhold Publishing Corporation, 1961. 
[2] N. D. Mermin. Crystalline order in two dimensions. Phys. Rev., 176:250–254, Dec 1968. 
[3] Ke Cao, Shizhe Feng, Ying Han, Libo Gao, Thuc Hue Ly, Zhiping Xu, and Yang Lu. Elasticstrainingoffree-standingmonolayergraphene. Nature Communications, 11(1):284, Jan 2020. 
[4] Jae-Ung Lee, Duhee Yoon, Hakseong Kim, Sang Wook Lee, and Hyeonsik Cheong. Thermal conductivity of suspended pristine graphene measured by raman spectroscopy. Phys. Rev. B, 83:081419, Feb 2011. 
[5] Luca Banszerus, Michael Schmitz, Stephan Engels, Jan Dauber, Martin Oellers, Federica Haupt, Kenji Watanabe, Takashi Taniguchi, Bernd Beschoten, and Christoph Stampfer. Ultrahigh-mobility graphene devices from chemical vapor de-position on reusable copper. Science Advances, 1(6), 2015. 
[6] Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero. Unconventional superconductivity in magic-angle graphene superlattices. Nature, 556(7699):43–50, Apr 2018. 
[7] Wei Han, Roland K. Kawakami, Martin Gmitra, and Jaroslav Fabian. Graphene spintronics. Nature Nanotechnology, 9(10):794–807, Oct 2014. 
[8] Yuanbo Zhang, Yan-Wen Tan, Horst L. Stormer, and Philip Kim. Experimental observation of the quantum hall effect and berry’s phase in graphene. Nature, 438(7065):201–204, Nov 2005. 
[9] Graphene Flagship. 	https://graphene-flagship.eu/news/Pages/Graphene-for<BR>Composites-Applications2.aspx. Accessed: 18.06.2020. 
[10] MediaDevil. 	https://mediadevil.com/products/artisanphonics-cb-01-nanene<BR>graphene-wood-earphones. Accessed: 18.06.2020. 
[11] MomoDesign. https://en.momodesign.com/products/fgtr-graphene-1-0. Accessed: 18.06.2020. 
[12] Vittoria. 	https://www.vittoria.com/eu/corsa-competition-race.html. Accessed: 18.06.2020. 
[13] Graphene-XT. 	https://www.graphene-xt.com/en/gxt-lubricant-2/. Accessed: 18.06.2020. 
[14] S. GrapheneTech. 	http://www.graphene-tech.net/product/conductive-ink/. Accessed: 18.06.2020. 
[15] Graphene-XT. 	https://www.graphene-xt.com/en/gxt-esd-2/. Accessed: 18.06.2020. 
[16] Emberion. https://www.emberion.com/products/vis-swir/. Accessed: 18.06.2020. 
[17] 
R. Ribeiro-Palau, F. Lafont, J. Brun-Picard, D. Kazazis, A. Michon, F. Cheynis, 

O. 
Couturaud, C. Consejo, B. Jouault, W. Poirier, and F. Schopfer. Quantum hall resistance standard in graphene devices under relaxed experimental conditions. Nature Nanotechnology, 10(11):965–971, Nov 2015. 


[18] Dominik Suszalski and Adam Rycerz. 	Adiabatic quantum pumping in buckled graphene nanoribbon driven by a kink, 2020. 
[19] R Saito, G Dresselhaus, and M S Dresselhaus. Physical Properties of Carbon Nanotubes. published by Imperial College Press and Distributed by World Scientific Publishing Co., 1998. 
[20] Michael Wimmer. Quantum transport in nanostructures: From computational concepts to spintronics in graphene and magnetic tunnel junctions, December 2009. 
[21] EdwardMcCannandMikitoKoshino.Theelectronicpropertiesofbilayergraphene. Reports on Progress in Physics, 76(5):056503, Apr 2013. 
[22] A. B. Kuzmenko, I. Crassee, D. van der Marel, P. Blake, and K. S. Novoselov. Determinationof the gate-tunableband gapand tight-binding parametersinbilayer graphene using infrared spectroscopy. Phys. Rev. B, 80:165406, Oct 2009. 
[23] Yuli V. Nazarov and Yaroslav M. Blanter. 	Quantum Transport: Introduction to Nanoscience. Cambridge University Press, 2009. 
[24] Keivan Esfarjani, Mona Zebarjadi, and Yoshiyuki Kawazoe. Thermoelectric properties of a nanocontact made of two-capped single-wall carbon nanotubes calculated within the tight-binding approximation. Phys. Rev. B, 73:085406, Feb 2006. 
[25] József Cserti, András Csordás, and Gyula Dávid. Role of the trigonal warping on the minimal conductivity of bilayer graphene. Phys. Rev. Lett., 99:066802, Aug 2007. 
[26] A. S. Mayorov, D. C. Elias, M. Mucha-Kruczynski, R. V. Gorbachev, T. Tudorovskiy, A. Zhukov, S. V. Morozov, M. I. Katsnelson, V.I. Fal ´ 
ko, A. K. Geim, and 
K. S. Novoselov. Interaction-driven spectrum reconstruction in bilayer graphene. Science, 333(6044):860–863, 2011. 
[27] Francisco Guinea, Antonio Castro Neto, and Nuno Peres. 	Electronic states and landau levels in graphene stacks. Physical Review B, 73, 04 2006. 
[28] Dominik Suszalski, Grzegorz Rut, and Adam Rycerz. Lifshitz transition and thermoelectric properties of bilayer graphene. Phys. Rev. B, 97:125403, Mar 2018. 
[29] L. M. Malard, J. Nilsson, D. C. Elias, J. C. Brant, F. Plentz, E. S. Alves, A. H. Castro Neto, and M. A. Pimenta. Probing the electronic structure of bilayer graphene by raman scattering. Phys. Rev. B, 76:201401, Nov 2007. 
[30] L. M. Zhang, Z. Q. Li, D. N. Basov, M. M. Fogler, Z. Hao, and M. C. Martin. Determination of the electronic structure of bilayer graphene from infrared spectroscopy. Phys. Rev. B, 78:235408, Dec 2008. 
[31] Melvin Cutler and N. F. Mott. Observation of anderson localization in an electron gas. Phys. Rev., 181:1336–1340, May 1969. 
[32] F. Freitag, J. Trbovic, M. Weiss, and C. Schönenberger. Spontaneously gapped ground state in suspended bilayer graphene. Phys. Rev. Lett., 108:076602, Feb 2012. 
[33] Wenzhong Bao, Jairo Velasco, Fan Zhang, Lei Jing, Brian Standley, Dmitry Smirnov, Marc Bockrath, Allan H. MacDonald, and Chun Ning Lau. Evidence for a spontaneous gapped state in ultraclean bilayer graphene. Proceedings of the National Academy of Sciences, 109(27):10802–10805, 2012. 
[34] Dominik Suszalski, Grzegorz Rut, and Adam Rycerz. Thermoelectric properties of gapped bilayer graphene. Journal of Physics: Condensed Matter, 31(41):415501, jul 2019. 
[35] Valeri N. Kotov, Bruno Uchoa, Vitor M. Pereira, F. Guinea, and A. H. Castro Neto. Electron-electron interactions in graphene: Current status and perspectives. Rev. Mod. Phys., 84:1067–1125, Jul 2012. 
[36] Dominik Suszalski, Grzegorz Rut, and Adam Rycerz. Conductivity scaling and the effects of symmetry-breaking terms in bilayer graphene hamiltonian. Phys. Rev. B, 101:125425, Mar 2020. 
[37] M. I. Katsnelson. Minimal conductivity in bilayer graphene. The European Physical Journal B -Condensed Matter and Complex Systems, 52(2):151–153, Jul 2006. 
[38] Youngwoo Nam, Dong-Keun Ki, Mikito Koshino, Edward McCann, and Alberto F Morpurgo. Interaction-induced insulating state in thick multilayer graphene. 2D Materials, 3(4):045014, oct 2016. 
[39] R. G. Chambers. Shift of an electron interference pattern by enclosed magnetic flux. Phys. Rev. Lett., 5:3–5, Jul 1960. 
[40] Akira Tonomura, Tsuyoshi Matsuda, Ryo Suzuki, Akira Fukuhara, Nobuyuki Osakabe, Hiroshi Umezaki, Junji Endo, Kohsei Shinagawa, Yutaka Sugita, and Hideo Fujiwara. Observation of aharonov-bohm effect by electron holography. Phys. Rev. Lett., 48:1443–1446, May 1982. 
[41] S. M. Roy. Condition for nonexistence of aharonov-bohm effect. Phys. Rev. Lett., 44:111–114, Jan 1980. 
[42] P. Bocchieri and A. Loinger. Charges in multiply connected spaces. 	Il Nuovo Cimento A (1965-1970), 66(2):164–172, Nov 1981. 
[43] Nobuyuki Osakabe, Tsuyoshi Matsuda, Takeshi Kawasaki, Junji Endo, Akira Tonomura, ShinichiroYano, andHirojiYamada. Experimentalconfirmationofaharonovbohm effect using a toroidal magnetic field confined by a superconductor. Phys. Rev. A, 34:815–822, Aug 1986. 
[44] Adam Rycerz. Magnetoconductance of the corbino disk in graphene. Phys. Rev. B, 81:121404, Mar 2010. 
[45] M. I. Katsnelson. Quantum transport via evanescent waves in undoped graphene. Journal of Computational and Theoretical Nanoscience, 8(6):912–918, Jun 2011. 
[46] Adam Rycerz and Dominik Suszalski. Graphene disk in a solenoid magnetic potential: Aharonov-bohm effect without a two-slit-like setup. Phys. Rev. B, 101:245429, Jun 2020. 
[47] 
C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, 

T. 
Taniguchi, P. Kim, K. L. Shepard, and et al. Boron nitride substrates for high-quality graphene electronics. Nature Nanotechnology, 5(10):722–726, Aug 2010. 


[48] A. Rycerz, J. Tworzydło, and C. W. J. Beenakker. Valley filter and valley valve in graphene. Nature Physics, 3(3):172–175, Mar 2007. 
[49] D W Boukhvalov and M I Katsnelson. Chemical functionalization of graphene. Journal of Physics: Condensed Matter, 21(34):344205, jul 2009. 
[50] B. Sachs, T. O. Wehling, M. I. Katsnelson, and A. I. Lichtenstein. Adhesion and electronic structure of graphene on hexagonal boron nitride substrates. Phys. Rev. B, 84:195414, Nov 2011. 
[51] Dominik Suszalski, Grzegorz Rut, and Adam Rycerz. Mesoscopic valley filter in graphene corbino disk containing a p–n junction. Journal of Physics: Materials, 3(1):015006, nov 2019. 
[52] Dominik Suszalski, Grzegorz Rut, and Adam Rycerz. Mesoscopic valley filter in graphene corbino disk containing a p–n junction. Journal of Physics: Materials, 3(1):015006, Nov 2019. 
[53] R.D.Yamaletdinov,V.A.Slipko,andY.V.Pershin.Kinksandantikinksofbuckled graphene: A testing ground for the ϕ4 field model. Phys. Rev. B, 96:094306, Sep 2017. 
[54] Pablo San-Jose, Elsa Prada, Henning Schomerus, and Sigmund Kohler. 	Laserinduced quantum pumping in graphene. Applied Physics Letters, 101(15):153506, 2012. 
[55] Yongjin Jiang, Tony Low, Kai Chang, Mikhail I. Katsnelson, and Francisco Guinea. Generation of pure bulk valley current in graphene. Phys. Rev. Lett., 110:046601, Jan 2013. 
[56] Evgeny Grichuk and E. Manykin. Adiabatic quantum pumping in graphene with magnetic barriers. The European Physical Journal B, 86(5):210, May 2013. 
[57] Babak Abdollahipour and Elham Moomivand. 	Magnetopumping current in graphene corbino pump. Physica E: Low-dimensional Systems and Nanostructures, 
86:204 – 209, 2017. 
[58] ManatoFujimoto,HenriKoschke,andMikitoKoshino.Topologicalchargepumping by a sliding moiré pattern. Phys. Rev. B, 101:041112, Jan 2020. 
[59] Dominik Suszalski and Adam Rycerz. 	Adiabatic quantum pumping in buckled graphene nanoribbon driven by a kink, 2020. 
[60] 	Properties of Carbon Nanotubes. Imperial College Press, 1998. 
[61] Dominik Suszalski and Adam Rycerz. Adiabatic pumping driven by moving kink and quantum standard ampere in buckled graphene nanoribbon, 2020. 
PHYSICAL REVIEW B 97, 125403 (2018) 

Lifshitz transition and thermoelectric properties of bilayer graphene 
Dominik Suszalski, Grzegorz Rut, and Adam Rycerz 
Marian Smoluchowski Institute of Physics, Jagiellonian University, Łojasiewicza 11, PL–30348 Kraków, Poland 
(Received 19 December 2017; published 5 March 2018) 
This is a numerical study of thermoelectric properties of ballistic bilayer graphene in the presence of a trigonal warping term in the effective Hamiltonian. We find, in the mesoscopic samples of the length L>10 μmat sub-Kelvin temperatures, that both the Seebeck coefficient and the Lorentz number show anomalies (the additional maximum and minimum, respectively) when the electrochemical potential is close to the Lifshitz energy, which can be attributed to the presence of the van Hove singularity in a bulk density of states. At higher temperatures the anomalies vanish, but measurable quantities characterizing the remaining maximum of the Seebeck coefficient still unveil the presence of massless Dirac fermions and make it possible to determine the trigonal warping strength. Behavior of the thermoelectric figure of merit (ZT) is also discussed. 
DOI: 10.1103/PhysRevB.97.125403 <BR>
I. INTRODUCTION 
It is known that thermoelectric phenomena provide valuable insight into the details of the electronic structure of graphene and other relativistic condensed-matter systems that cannot be <BR>solely <BR>determined <BR>by <BR>conductance <BR>measurements <BR>[1]. <BR>Such <BR>a fundamental perspective has inspired numerous studies on Seebeck and Nernst effects in mono-(MLG) and bilayer (BLG) graphenes <BR>[2–9] as well as in other two-dimensional systems <BR>[10–13]. The exceptionally high thermal conductivity of graphenes <BR>has <BR>also <BR>drawn <BR>significant <BR>attention <BR>[14–19]after a seminal work by Balandin et al. [20]. <BR>A <BR>separate <BR>issue <BR>concerns <BR>thermal and thermoelectric properties of tailor-made graphene systems <BR>[1,21–26], <BR>including <BR>superlattices <BR>[21], nanoribbons [22–25], <BR>or <BR>defected <BR>graphenes <BR>[25,26], for which peculiar electronic structures may result in high thermoelectric figures of merit ZT >2 <BR>at <BR>room <BR>temperature <BR>[24,25]. 
Unlike in conventional metals or semiconductors, thermoelectric power in graphenes can change a sign upon varying the gate <BR>bias <BR>[2–4], making it possible to design thermoelectronic devices <BR>that <BR>have <BR>no <BR>analogues <BR>in <BR>other <BR>materials <BR>[27]. <BR>In <BR>BLG <BR>the <BR>additional <BR>band-gap <BR>tunability <BR>[28–30] was utilized to noticeably enhance the thermoelectric power in a dual-gated setup <BR>[8]. <BR>
At sufficiently low temperatures, one can expect thermoelectric properties of BLG to reflect most peculiar features of its electronic structure. These features include the presence (in the gapless case) of three additional Dirac points in the vicinity of each of the primary Dirac points K and K1 [31–34]. In turn, when varying chemical potential the system is expected to undergo the Lifshitz transition at μ=±EL (the Lifshitz energy) <BR>[34]. <BR>What <BR>is <BR>more, <BR>electronic <BR>density <BR>of <BR>states <BR>(DOS) <BR>shows van Hove singularities at μ=±EL. Unlike in systems with Mexican-hat band dispersion, for which diverging DOS appears at the bottom of the conduction band and at the top of the <BR>valence <BR>band <BR>[12,13], in BLG each van Hove singularity separates populations of massless Dirac-Weyl quasiparticles (|μ|<EL) with approximately conical dispersion relation, and massive chiral quasiparticles (|μ|>EL) characterized by the effective mass meff ≈0.033 me, with me being the free-electron mass. Although the value of EL is related to several directly-measurable quantities, such as the minimal conductivity <BR>[35–37], available experimental results cover the full range of EL ∼0.1–1 <BR>meV <BR>[34]. <BR>
The purpose of this paper is to show that thermoelectric measurements in ballistic BLG (see Fig. 1) <BR>can <BR>provide <BR>insights <BR>into the nature of quasiparticles near the charge-neutrality point and allow one to estimate the Lifshitz energy. We consider a relatively large, rectangular sample of ballistic BLG (with the length L=17.7 μm and the width W=20 L) and calculate its basic thermoelectric properties (including the Seebeck coefficient Sand the Lorentz number L) within the Landauer-Büttiker <BR>formalism <BR>[38,39]. Our main findings are outlined in Fig. 1, <BR>where <BR>NS min)—the number of maxima 
(NL 
max 
(minima) of S(L) appearing for μ>0 is indicated in the EL -T parameter plane. For instance, a handbook value of EL/kB ≈ 10 <BR>K <BR>[40] <BR>leads <BR>to <BR>the <BR>anomalies, <BR>including <BR>additional <BR>ex-trema at μ≈EL, at sub-Kelvin temperatures. We further show that even for T 2 1 K (at which NS =NL =1)
max min 
the value of EL determines the carrier concentration corresponding to the remaining maximum of S(or the minimum of L). 
The paper is organized as follows: The model and theory are described in Sec. II, <BR>followed <BR>by <BR>the <BR>numerical <BR>results <BR>and discussions on the conductance, thermopower, validity of the Wiedemann-Franz law, the role of phononic thermal conductivity, and the figure of merit (Sec. III). <BR>A <BR>comparison <BR>with the linear model for transmission-energy dependence (see Appendix) is also included. The conclusions are given in Sec. IV. <BR>
II. MODEL AND THEORY 
A. The Hamiltonian 
We start our analysis from the four-band effective Hamiltonian <BR>for <BR>low-energy <BR>excitations <BR>[34], <BR>which <BR>can <BR>be <BR>
2469-9950/2018/97(12)/125403(10) 125403-1 ©2018 American Physical Society DOMINIK SUSZALSKI, GRZEGORZ RUT, AND ADAM RYCERZ PHYSICAL REVIEW B 97, 125403 (2018) 
300 
200 


t=0 0.2eV 0.3eV  (d) 0.2eV t=0 0.3eV (c)  

(2π/W)×Nopen [1/nm] 
CARRIER DENSITY [ μm −2 ] 
100 
0 
0.12 
0.08 
0.04 
-1.0 -0.5 0 0.5 1.0 1.5 
kx/kL 
1.0 
(b)  

0.5 
0 
0 
0 0.5 1.0 1.5 2.0 0 0.4 0.8 1.2 

E/EL E [meV]
FIG. 1. The system studied numerically in the paper (inset) and the thermoelectric phase diagram (main) for Bernal-stacked bilayer graphene. Thermal and electric currents (IQ,I) flow between the leads, modeled as infinitely-doped graphene regions with electrochemical potentials μL(R) →∞(at a fixed μR −μL ≡eV, with the electron charge −e and the voltage V), and temperatures TL(R), attached to the rectangular sample. Additional gate electrodes (not shown) are used to tune the chemical potential in the sample area μ at zero bias between the layers. The number of distinct maxima of the Seebeck coefficient (NS ) and minima of the Lorentz factor 
max
(NL 
min), occurring for 0 <μ<∞, are indicated in the Lifshitz energy– temperature parameter plane. The border of the Fabry-Pérot transport regime T<TF−P (in which NS min »1), corresponding to the 
max,NL system length L =104 l⊥=17.7 μm, is also depicted. 
written as 
⎛⎞ 
0 vFπ ξt⊥ 0 
⎜⎟
vFπ† 00 v3π 
⎜⎟
H =ξ, (1) 
⎝ ξt⊥ 00 vFπ† ⎠ 0 v3π† vFπ 0 
where the valley index ξ =1(−1) for K (K1) valley, vF =
√ 3 t0a/(2¯h) '106 m/s is the asymptotic Fermi velocity de
fined via the intralayer hopping t0 =3.16 eV and the lattice paiθ(−
rameter a =0.246 nm, π =¯i∂x ), θ denotes the 
he−+∂yangle between the main system axis and the armchair direction. (For the numerical calculations, we set ¯hvF =0.673 eV nm.) The nearest-neighbor interlayer hopping is t⊥=0.381 <BR>eV <BR>[41] <BR>
1
defining l⊥=hvF/t⊥ v3 1/t0
¯=1.77 nm, and =vFtwith tbeing the next-nearest neighbor (or skew) interlayer hopping. The <BR>Hamiltonian <BR>(1) <BR>leads <BR>to <BR>the <BR>bulk <BR>dispersion <BR>relation <BR>for <BR>electrons <BR>[31,34] 
[fr
()1/2
(e) 12 212
E± (k) =2 t⊥+vF +2 v3k2 ±r(k), 
FIG. 2. Physical consequences of the dispersion relation given by <BR>Eq. <BR>(2). <BR>(a) <BR>Equienergy <BR>surfaces <BR>for <BR>E =0.5 EL (gray solid lines), E =EL (red solid line), and E =1.5 EL (blue dashed line) for the crystallographic orientation θ =0. (b) Density of states (purple solid line) <BR>and <BR>the <BR>approximating <BR>expressions <BR>given <BR>by <BR>Eqs. <BR>(7), <BR>(8), <BR>and <BR>(9) <BR>(gray solid, black dashed, and black dotted-dashed line, respectively). 
(c) 
Carrier density and (d) the number of open channels for different values of t1 (specified for each line). Solid and dashed lines in panel 

(d) 
corresponds to θ =0and θ =π/6, dotted-dashed lines represent the <BR>approximating <BR>Eq. <BR>(12). <BR>


where k ≡(kx,ky) is the in-plane wave vector (with k =0 referring to the K or K1 point), k =|k|, and the angle 0 ≤ ϕ< 2π can be defined as the argument arg z of a complex number 
−iθ(kx
z =e +iky). (3) 
(h)(e)
For holes, we have E± (k) =−E± (k)[42]. <BR>
B. Low-energy electronic structure 
Basic <BR>consequences <BR>of <BR>Eq. <BR>(2) <BR>are <BR>illustrated <BR>in <BR>Fig. <BR>2.In <BR>the energy range |E|<EL, with the Lifshitz energy 
1  v3  2 
EL = t⊥, (4)
4 vF
there are four distinct parts of the Fermi surface [see Fig. 2(a)], <BR>centered <BR>at <BR>z =z0,...,z3, where z0 =0, zj = kL exp(2πij/3),j =1,2,3, and 
kL = t⊥v3 . (5)
¯
hv2 
F 
For |E|≥EL the Fermi surface becomes connected, and the transition at E =±EL is accompanied by the van Hove singularity in the density of states ρ(E)[seeFig. 2(b)], <BR>which <BR>can be defined (for electrons) via
(2)22(2)
1222 22 222
r(k) =t⊥−h¯v3 k+h¯vkt⊥+h¯v3 k 
4FE A(E) 
32 dE1ρ(E1) ≡n(E) = , (6)
+2ξt⊥h¯v3vk3 cos 3ϕ, (2) π2
F 
0 
125403-2 LIFSHITZ TRANSITION AND THERMOELECTRIC … PHYSICAL REVIEW B 97, 125403 (2018) 
where n(E) is the physical carrier density (taking into ac-count spin and valley degeneracies gs =gv =2) depicted in Fig. 2(c), <BR>and <BR>A(E) denotes the area bounded by the Fermi surface in the (kx,ky) <BR>plane <BR>[43]. <BR>In <BR>particular, <BR>taking <BR>the <BR>limit <BR>of t1→0wehave 
2meff t⊥ 
ρt1→0(|E|«t⊥) ≈= ≡ρ0, (7)
2
πh¯π(¯hvF)2 
where we have introduced the effective mass relevant in the absence of trigonal warping (EL =0).Atfinite t1 (EL > 0) the value of ρ0 defined <BR>in <BR>Eq. <BR>(7) <BR>is <BR>approached <BR>by <BR>the <BR>actual <BR>
1
ρ(E)for |E|2 EL [see Fig. 2(b)]. <BR>Also, <BR>in <BR>the <BR>t= 0 case, we find that the approximating formula 
ρ0|E|
ρ(E) ≈ [1 +0.33(E/EL)2], (8)
EL reproduces the actual ρ(E) with 1% accuracy for |E|≤EL/2 (being the energy interval most relevant for discussion presented in the remaining parts of this paper). For |E|«EL, leaving <BR>only <BR>the <BR>leading <BR>term <BR>on <BR>the <BR>right-hand <BR>side <BR>of <BR>Eq. <BR>(8) <BR>brought us to 
4|E|
ρ(|E|«EL) ≈ , (9)
π(¯hv3)2 
which can be interpreted as a double-monolayer DOS with the Fermi velocity replaced by v3. 
Although the trigonal-warping effects become hardly visible in ρ(E)for E 2 EL, characteristic deformations of the Fermi surface can be noticed also for E »EL. We point out that a compact quantity taking this fact into account, which can <BR>be <BR>determined <BR>directly <BR>from <BR>Eq. <BR>(2) <BR>without <BR>resorting <BR>to <BR>quantum transport simulations, is the number of propagating modes (open channels) Nopen(θ,E) presented in Fig. 2(d).It <BR>can be defined as a total number of solutions, with real kx,of equations 
(p)(p)
E(kx,q1ky) =E, E− (kx,q1ky) =E, (10)
+ 
where p =e for electrons (E> 0) or p =h for holes (E< 0), q =0, ±1,±2,..., and 1ky =2π/W (we suppose the periodic boundary conditions along the y axis) that correspond to a chosen sign of the group velocity, e.g., 
)(p)(p)
(vg±,q =∂E± (kx,q1ky)/∂ky > 0. Apart from the t1→0 limit, for which 
 
|E|t⊥ 1 
Nopen(t1→0,|E|«t⊥) ≈2, (11)
(¯hvF)2 1ky 
the number of open channels is anisotropic and shows the periodicity with a period π/3. In the low-energy limit 
|E| 1 
Nopen(θ,|E|«EL) ≈2F(θ) , (12)
¯1ky
hv3 
where 
   1/2 
L 
82π 
F(θ) =1 + 1 − cos2 θ + j,
93 
j=1,2,3
≈3.126 +0.029 cos 6θ. (13) 
The anisotropy is even more apparent for |E|2 EL. In particular, Nopen(θ =0,E) grows monotonically with increasing E, whereas Nopen(θ =π/6,E) has a shallow minimum at E ≈1.11EL. 
It is also visible in Fig. 2(d) <BR>that the effects of increasing tare essentially opposite at different energy ranges: For |E|2 EL, Nopen grows systematically with t1;for |E|«EL we have Nopen ∝1/t1 following <BR>from <BR>Eq. <BR>(12). <BR>Such <BR>a <BR>feature <BR>has <BR>no <BR>analogues in behaviors of other characteristics presented in Fig. 2. <BR>
C. Thermoelectric properties 
In the linear-response regime, thermoelectric properties of a generic nanosystem in graphene are determined via Landauer-Büttiker expressions for the electrical and thermal currents [44,45] 
gsgve 
I =−dET(E)[fL(E)−fR(E)], (14)
h 
IQ = gsgv dET(E)[fL(E)−fR(E)](E−μ), (15)
h 
where gs =gv =2 are spin and valley degeneracies, T(E) ≡ Tr(tt†) with t being <BR>the <BR>transmission <BR>matrix <BR>[36], <BR>fL(R) is the distribution functions for the left (right) lead with electrochemical potential μL(R) and temperature TL(R). Assuming that μL −μR ≡−eV and TL −TR ≡1T are infinitesimally small [hereinafter, we refer to the averages μ =(μL +μR)/2 and T =(TL +TR)/2], we obtain the conductance G,the Seebeck coefficient S, and the electronic part of the thermal conductance Kel, <BR>as <BR>follows <BR>[39] <BR>
 
I  
 
G ==e 2L0, (16) 
 
V
1T=0 
 
V  L1 
 
S =− = , (17) 
 
1TI=0 eTL0 
 
IQ  L0L2 −L2 
 1
Kel == , (18) 
 
1TI=0 TL0 
where Ln (with n =0,1,2) is given by 
gsgv ∂fFD 
Ln = dET(E) − (E −μ)n, (19)
h ∂E 
with fFD(μ,T,E) =1/[exp ((E−μ)/kBT) +1 ] the Fermi-Dirac distribution function. 
By definition, the Lorentz number accounts only the electronic part of the thermal conductance, 
Kel L0L2 −L2 
L == 1 . (20)
TG e2T2L2 
0 
The thermoelectric figure of merit accounts the total thermal conductance (Ktot =Kel +Kph) 
TGS2 Kel L2 
ZT == 1 , (21)Ktot Kel +Kph L0L2 −L2 
1 
where the phononic part can be calculated using 
1 ∂fBE 
Kph = hω (22)
dω ¯Tph(ω),
2π ∂T 
with fBE(T,ω) =1/[exp (¯hω/kBT) −1 ] the Bose-Einstein distribution function and Tph(ω) the phononic transmission 
125403-3 DOMINIK SUSZALSKI, GRZEGORZ RUT, AND ADAM RYCERZ PHYSICAL REVIEW B 97, 125403 (2018) 
spectrum. For BLG in a gapless case considered in this work, we typically have Kph ∼Kel (see Sec. III <BR>D)[46]. As Tph(ω)inEq. (22) <BR>is <BR>generally <BR>much <BR>less <BR>sensitive <BR>to <BR>external electrostatic fields than T (E)inEq. <BR>(18) <BR>it <BR>should <BR>be possible—at least in principle—to independently determine Kph and Kel in the experiment. 
It can be noticed that ultraclean ballistic graphene shows approximately linear transmission to Fermi-energy dependence T (E) ∝|E|(where E =0 corresponds to the chargeneutrality <BR>point) <BR>[47–50]. Straightforward analysis (see Appendix A) <BR>leads <BR>to <BR>extremal <BR>values <BR>of <BR>the <BR>Seebeck <BR>coefficient <BR>as a function of the chemical potential 
Smax =−Smin ≈kB/e =86 μV/K, (23) 
providing yet another example of a material characteristic given <BR>solely <BR>by <BR>fundamental <BR>constants <BR>[47,51]. Similarly, the Lorentz number reaches, at μ =0, the maximal value given by 
9 ζ(3) kB 2 
Lmax ==2.37 LWF, (24)
2ln 2 e 
1
with LWF =π2(kB/e)2 being the familiar Wiedemann-Franz 
3 
constant. Although the disorder and electron-phonon coupling may affect the above-mentioned values, existing experimental works report Smax and Lmax close <BR>to <BR>those <BR>given <BR>by <BR>Eqs. <BR>(23) <BR>and <BR>(24) <BR>for <BR>both <BR>MLG <BR>and <BR>BLG, <BR>provided <BR>the <BR>temperature <BR>is <BR>not <BR>too <BR>low <BR>[2–7,18]. <BR>
At low temperatures, the linear model no longer applies, partly <BR>due <BR>to <BR>the <BR>contribution <BR>from <BR>evanescent <BR>modes <BR>[47,52] and partly due to direct trigonal-warping effects on the electronic structure (see Sec. II <BR>B). <BR>For <BR>this <BR>reason, <BR>thermoelectric <BR>properties <BR>calculated <BR>numerically <BR>from <BR>Eqs. <BR>(16)–(21)are discussed next. 
III. 
RESULTS AND DISCUSSION 

A. 
Zero-temperature conductivity 


For T →0Eq. <BR>(16) <BR>leads <BR>to <BR>the <BR>conductivity <BR>
G(T →0)Lg0L 
σ(T →0) == Tr(tt†), (25)
WW 
with the conductance quantum g0 =4e 2/h. As the right-hand side <BR>of <BR>Eq. <BR>(25) <BR>is <BR>equal <BR>to <BR>Tr(tt†) with a constant prefactor, σ(T →0) gives a direct insight into the transmission-energy dependence that defines all the thermoelectric properties [see Eqs. <BR>(16)–(21)]. 
In order to determine the transmission matrix t for a given electrochemical potential μ we employ the computational scheme <BR>similar <BR>to <BR>that <BR>presented <BR>in <BR>Ref. <BR>[36]. <BR>However, <BR>at <BR>finite<BR>precision arithmetics, the mode-matching equations become ill defined for sufficiently large L and μ, as they contain both exponentially growing and exponentially decaying coefficients. This difficulty can be overcome by dividing the sample area into Ndiv consecutive, equally-long parts, and matching wave functions for all Ndiv +1 <BR>interfaces <BR>[53]. <BR>
Numerical results are presented in Fig. 3. <BR>A <BR>striking <BR>feature <BR>of all datasets is the presence of quasiperiodic oscillations of the Fabry-Pérot type. Although such oscillations can be regarded as artifacts originating from a perfect, rectangular 

FIG. <BR>3. <BR>Zero-temperature <BR>conductivity <BR>[see <BR>Eq. <BR>(25) <BR>in <BR>the <BR>main <BR>text] for L =W/20 =104 l⊥=17.7 μmand θ =0[54] <BR>plotted <BR>as <BR>a <BR>function of the chemical potential. The value of skew-interlayer hop-ping t1is specified for each line. Remaining tight-binding parameters are <BR>given <BR>below <BR>Eq. <BR>(1) <BR>in <BR>the <BR>main <BR>text. <BR>Vertical <BR>lines <BR>mark <BR>values <BR>of <BR>the <BR>Lifshitz <BR>energy, <BR>given <BR>by <BR>Eq. <BR>(4), <BR>for <BR>t1=0.2eV and t1=0.3eV. Inset is a zoom in, for low chemical potentials, with horizontal lines depicting σ =2 σMLG =(8/π) e 2/h and σ =6 σMLG. 
shape of the sample area (vanishing immediately when, e.g., samples <BR>with <BR>nonparallel <BR>edges <BR>are <BR>considered, <BR>see <BR>Ref. <BR>[55]) <BR>their periodic features are useful to benchmark the numerical procedure applied. 
1
In particular, for t=0, the conductivity shows abrupt features at energies associated with resonances at normal incidence (ky =0) <BR>[52], <BR>namely <BR>
��2
πn
1
En(thvFl⊥ (26)
=0) ≈±¯,n =1,2,3,..., 
L 
where the approximation refers to the parabolic dispersion relation applying for |En|«t⊥, or equivalently for n « L/(πl⊥) ≈3180 in our numerical example. In turn, the separation between consecutive resonances is 
1 2n+1 π ¯
hvF 
1En(t=0) =|En+1 −En|≈ 
t⊥ L 
π ¯||

hvF En
≈2 , (27)
Lt⊥ 
with the last approximation corresponding to n »1. 
1
For t=0 the analysis is much more cumbersome even at low energies, as we have resonances associated with four distinct Dirac cones. However, resonances at normal incidence associated with the central cone, occurring at En ≈π ¯
hv3n/L (n =±1, ±2,... ), allow us to estimate the order of magnitude of the relevant separation as 
1 πhv¯3 π ¯EL
hvF 
1En(t=0) ∼=2 ≡kBTF−P, (28)
L Lt⊥ 
finding that the period of Fabry-Pérot oscillations is now energy independent and should be comparable with 1En(0) 
125403-4 LIFSHITZ TRANSITION AND THERMOELECTRIC … PHYSICAL REVIEW B 97, 125403 (2018) 
given <BR>by <BR>Eq. <BR>(27)for <BR>μ =En ≈EL(t1). The data displayed in Fig. 3 <BR>show that the oscillation period is actually energy independent in the surprisingly-wide interval of μ � 1.5EL, with the multiplicative factor 1En(t1)/1En(0)|EL(t1) ≈3.
μ=
The oscillation amplitude is also enhanced, in comparison to the t=0 case, for μ � 1.5EL.For μ 2 1.5EL, both the oscillation period and amplitude are noticeably reduced, resembling 
1
the oscillation pattern observed for the t=0 case. It is also visible in Fig. 3, <BR>that <BR>the <BR>mean <BR>conductivity <BR>(averaged <BR>over <BR>the oscillation period) linearly increase with μ for μ � EL, with a slope weakly dependent on t1. Such a behavior indicates that Tr(tt†) <Nopen [see Fig. 2(d)] <BR>which <BR>can <BR>be <BR>interpreted <BR>as <BR>a backscattering (or transmission reduction) appearing when different classes of quasiparticles are present in the leads and in the sample area. For larger μ, the transmission reduction is still significant, but its dependence on t1 is weakened, and the sequence of lines from Fig. 2(d) <BR>is reproduced. 
A detailed explanation of the above-reported observations, in terms of simplified models relevant for |μ|«EL and for |μ|2 EL, will be presented elsewhere. Here we only notice that the linear model for transmission-energy dependence is justified, for |μ|� EL, with the numerical results presented in Fig. 3. <BR>
The <BR>rightmost <BR>equality <BR>in <BR>Eq. <BR>(28) <BR>defines <BR>the <BR>Fabry-Pérot <BR>temperature, which can be written as 
πt⊥t1l⊥ t1l⊥ 
TF−P ==13890 K × . (29)
kBt0Lt0L 
1
For L =104 l⊥, we obtain TF−P =88mK if t=0.2eV, 1
or TF−P =132 mK if t=0.3 eV. For higher temperatures, Fabry-Pérot oscillations are smeared out due to thermal excitations involving transmission processes from a wider energy window <BR>[see <BR>Eqs. <BR>(16) <BR>and <BR>(19)]. 
B. Thermopower and Wiedemann-Franz law 
As the finite-T conductivity is simply given by a convolution of T(E) =Tr(tt†) with the derivative of the Fermi-Dirac function, we proceed directly to the numerical analysis of the Seebeck coefficient and the Lorentz number given by Eqs. <BR>(17)–(20)[56]. <BR>In <BR>Fig. <BR>4, <BR>these <BR>thermoelectric <BR>properties <BR>
1
are displayed as functions of μ,for afixed t=0.3eV (corresponding to EL/kB ≈10 K) and varying temperature. Quasiperiodic oscillations are still prominent in datasets for the lowest presented temperature, T =80 mK ≈0.6 TF−P, although it is rather close to TF−P. This is because all the abrupt features of T(E) are magnified when calculating S,or L, since they affect the nominator and the denominator in the corresponding <BR>Eq. <BR>(17), <BR>or <BR>Eq. <BR>(20), in a different manner. For T =0.2K ≈1.5 TF−P the oscillations vanish for S and are strongly suppressed for L; instead, we observe the anomalies: the secondary maximum of S and minimum of L, located near μ =EL. The secondary maximum of S vanishes for T =T/S=0.515 K, but L still shows the two shallow minima at this temperature. (We find that the minima of L merge at T/L=1.20 K =2.33 TS, the corresponding dataset is omitted 
/ 
for clarity.) For T =2 K, each of S and L shows a single extremum for μ> 0. 
The crossover temperatures TS and T/Las functions of EL,
/ 
varied in the range corresponding to 0.1eV ≤t1≤0.35 eV, 

FIG. 4. Seebeck coefficients S and Lorentz number L for t1= 
0.3 eV, as a function of the chemical potential. The temperature is specified for each line in the top panel and is the same in both panels. Vertical lines mark the Lifshitz energy; horizontal line in bottom panel corresponds to the Wiedemann-Franz value L =LWF =31 π2(kB/e)2. Inset shows crossover temperatures, corresponding to vanishing of the secondary maximum of S (triangles) and minimum of L (circles), plotted as functions of the Lifshitz energy, together with the best-fitted linear <BR>functions <BR>[see <BR>Eqs. <BR>(30)and <BR>(31)]. 
are also plotted in Fig. 4 <BR>(see the inset). The least-squares fitted lines are given by 
T/,Sfit =0.0504(5) ×EL/kB, (30) 
L
T/,fit =0.1176(3) ×EL/kB, (31) 
with standard deviations of the last digit specified by numbers in parentheses. 
These findings can be rationalized by referring to the onset on low-energy characteristics given in Sec. II <BR>B <BR>(see Fig. 2). <BR>In <BR>particular, the abrupt features of T(E) near E =EL, attributed to the van Hove singularity of ρ(E) shown in Fig. 2(b),ortothe <BR>anisotropy of Nopen(θ,E)inFig. 2(d), <BR>are <BR>smeared <BR>out <BR>when <BR>calculating thermoelectric properties for energies of thermal excitations 
kBT 2 0.1 EL. (32) 
However, some other features, related to trigonal-warping effects on Nopen(θ,E)or n(E)[seeFig. 2(c)]awayfrom <BR>
125403-5 
DOMINIK SUSZALSKI, GRZEGORZ RUT, AND ADAM RYCERZ 	PHYSICAL REVIEW B 97, 125403 (2018)


0 2 4 6 8 10 12 

/kBT 
FIG. 5. Same as Fig. 4 <BR>but plotted versus the dimensionless variable μ/(kBT ). Solid lines in both panels represent the datasets for T =0.2K and T =2 K. Dashed-dotted lines correspond to Eqs. <BR>(A4)and <BR>(A5) in Appendix, following from the linear model for transmission-energy dependence. Dashed line in bottom panel marks L =LWF. Inset shows the finite-T conductivity (solid lines) σ =GL/W [see <BR>Eqs. <BR>(16)and <BR>(19) in the main text] and the linear fit (dash-dot line) to the corresponding T =0 dataset in Fig. 3. <BR>
E =EL, visible in thermoelectric properties, may even be observable at higher temperatures. 
C. Comparison with the linear model for transmission-energy dependence 
In Fig. 5 <BR>we display the selected numerical data from Fig. 4,for <BR>T =0.2 K and T =2 K, as functions of μ/(kBT ) [solid lines] in order to compare them with predictions of the linear model for transmission-energy dependence T (E) ∝|E|[dashed-dotted lines] elaborated in Appendix. For T =2K, both S and L show an agreement better than 10% with the linear model for μ EL ≈5 kBT .For T =0.2 K, larger deviations appear for low chemical potentials due to the influence of transport via evanescent waves, which are significant for μ< ¯For larger μ, a few-percent agreement 
hvF/L ≈2–3 kBT . 
with the linear model is restored and sustained as long as 

FIG. 6. Relative electronic contribution to the thermal conductance for t1=0.3 eV as a function of the chemical potential. The temperatures are (top to bottom) T =80mK, 0.2K, 0.515 K, and 2 K. Inset shows phononic [blue line] and electronic [red line] thermal conductivities (κ =KL/(2dW), with d =0.335 mn the separation between layers) as functions of temperature, with the chemical potential fixed at μ =μmax corresponding to the maximal Seebeck coefficient. 
Another remarkable feature of the results presented in Fig. 5 <BR>becomes apparent when determining the extrema: The maximal thermopower corresponds to μ(S) /(kBT ) =2.0at 
max
T =2K, or to μ(S) /(kBT ) =2.2at T =0.2 K; the minimal 
max(L)
Lorentz number corresponds to μmin/(kBT ) =5.4at T =2K, 
(L)
or to μmin/(kBT ) =4.6at T =0.2 K. In other words, an almost perfect agreement with the linear model [see, respectively, <BR>the <BR>second <BR>equality <BR>in <BR>Eq. <BR>(A6), <BR>or <BR>the <BR>second <BR>equality <BR>in <BR>Eq. <BR>(A7) <BR>in <BR>Appendix] <BR>is <BR>observed <BR>provided <BR>that <BR>
¯
hvF 
(S)(L)
« kBT ∼μmax ∼μmin « EL. (33)L 
In consequence, the effects that we describe may be observable for the sample length L> 10 μm. 
D. Electronic and phononic parts of the thermal conductance 
Before discussing the thermoelectric figure of merit ZT we first display, in Fig. 6, <BR>values <BR>of <BR>the <BR>dimensionless <BR>prefactor <BR>in <BR>the <BR>last <BR>expression <BR>of <BR>Eq. <BR>(21), <BR>quantifying <BR>relative <BR>electronic <BR>contribution to the thermal conductance. The phononic transmission <BR>spectrum <BR>[see <BR>Eq. <BR>(22)] <BR>was <BR>calculated <BR>numerically <BR>by <BR>employing, for the sample length L =17.7 μm, the procedure presented <BR>by <BR>Alofi <BR>and <BR>Srivastava <BR>[57] <BR>adapting <BR>the <BR>Callaway <BR>theory <BR>[58] <BR>for <BR>mono-and <BR>few-layer <BR>graphenes <BR>[59]. <BR>The <BR>results show that in sub-Kelvin temperatures the electronic contribution usually prevails, even if the system is quite close to the charge-neutrality point, as one can expect for a gapless conductor. For T> 1 K, however, the phononic contribution overrules the electronic one in the full range of chemical 
μ EL ≈40 kBT . 	potential considered. 
125403-6 LIFSHITZ TRANSITION AND THERMOELECTRIC … PHYSICAL REVIEW B 97, 125403 (2018) 

FIG. 7. Effective carrier concentration (n−p)max (a) and chemical potential μmax (b) corresponding to the maximal Seebeck coefficient Smax (c) as functions of temperature. The figure of merit ZT(μmax) (d) is also displayed. Solid lines represent the numerical results for different values of t1 [specified in panel (a)]. Dashed lines mark predictions of the linear model for transmission-energy dependence T(E) ∝|E|. 
A direct comparison of the phononic and the electronic and thermal conductivities calculated in the physical units (see inset in Fig. 6) <BR>further <BR>shows <BR>that, <BR>if <BR>the <BR>chemical <BR>potential <BR>is adjusted to μmax ≡ μ(S) for a given temperature, both 
max 
properties are of the same order of magnitude up to T = 2K. Also, for μ = μmax, we find that Kel = Kph at the temperature 
1
Tel−ph ≈ 0.3 K, which is almost insensitive to the value of t. 
E. Maximal performance versus temperature 
In Fig. 7 <BR>we present parameters characterizing the maximal thermoelectric performance for a given temperature (0 < T ≤ 2 K). As the existing experimental works refer to the carrier concentration rather than to the corresponding chemical potential, we focus now on the functional dependence of the former on T (and t1). 
Taking into account that the maximal performance is expected for μ ∼ kBT (see previous subsection), and that a gapless system is under consideration, one cannot simply neglect the influence of minority carriers. For the conduction band (μ> 0), the effective carrier concentration can be written as 
∞ 
n−p =  dE ρ(E)f(μ,E)  
0  
0  
−  dE ρ(E)[1 − f(μ,E)],  (34)  
−∞  

where we have supposed the particle-hole symmetry ρ(E) = ρ(−E). [For the valence band (μ< 0), the effective concentration p-n is simply given by the formula on the right-hand side of <BR>Eq. <BR>(34) <BR>with <BR>an <BR>opposite <BR>sign.] <BR>Next, <BR>the <BR>approximating <BR>Eqs. <BR>(7) <BR>and <BR>(8) for the density of states lead to 
n−p ≈ ρ0kBT 
 
1 
y, if t= 0, 
× 2  (35)
1
τLI1(y) + 0.33 τI3(y), if t= 0,
L 
where y = μ/kBT, τL = kBT/EL, and we have defined 
∞∞
(x + y)n (x − y)n 
In(y) = dx − dx. (36) 
ex + 1 ex + 1
−yy 
(In particular, I0(y) = y.) Numerical evaluation of the integrals <BR>in <BR>Eq. <BR>(35)for <BR>y = ymax given <BR>by <BR>Eq. <BR>(A6) <BR>in <BR>Appendix <BR>brought us to 
(n−p)max ≈ ρ0kBT 
 
1.949 if t1= 0, 
×() (37)3.269 τL1+3.23 τ2 if t1= 0.
L
In turn, the carrier concentration corresponding to the maximum of S for a given T is determined by the value of EL. [A similar expression for the minimum of L,see <BR>Eq.(A7)in <BR>Appendix, is omitted here.] 
Solid lines in Fig. 7(a) <BR>show the values of (n−p)max calculated <BR>from <BR>Eq. <BR>(34) <BR>for <BR>the <BR>actual <BR>density <BR>of <BR>states <BR>and <BR>the chemical potential μ = μmax [displayed with solid lines in Fig. 7(b)] <BR>adjusted <BR>such <BR>that <BR>the <BR>Seebeck <BR>coefficient, <BR>obtained <BR>numerically <BR>from <BR>Eq. <BR>(17), <BR>reaches <BR>the <BR>conditional maximum (Smax)[seeFig. 7(c)] <BR>at <BR>a <BR>given <BR>temperature <BR>T (and one of the 
1
selected values of t= 0, 0.2eV, or 0.3 eV). The numerical results are compared with the linear-model predictions (dashed lines <BR>in <BR>all <BR>panels), <BR>given <BR>explicitly <BR>by <BR>Eq. <BR>(37) <BR>[Fig. <BR>7(a)]or <BR>Eq. <BR>(A6) <BR>in <BR>Appendix <BR>[Figs. <BR>7(b) <BR>and 7(c)]. <BR>Again, <BR>the <BR>linear <BR>model shows a relatively good agreement with corresponding data obtained via the mode-matching method; moderate devia
1
tions are visible for t= 0 when μmax 2 EL/2. In such a range, both ρ(E) <BR>no <BR>longer <BR>follows <BR>the <BR>approximating <BR>Eq. <BR>(9), <BR>and <BR>the <BR>sudden rise of T(E) near E ≈ EL starts to affect thermoelectric properties. 
Figures 7(c) <BR>and Fig. 7(d) <BR>display, respectively, the maximal Seebeck coefficient (Smax) and figure of merit [ZT(μmax)] 
1
as functions of temperature. For t= 0, the former shows broad peaks, centered near temperatures corresponding to μmax ≈ 0.4 EL, for which the prediction of the linear model [see <BR>Eq. <BR>(A6) <BR>in <BR>Appendix] <BR>is <BR>slightly <BR>exceeded <BR>(by <BR>less <BR>then <BR>
1
10%), whereas for t= 0 a monotonic temperature dependence, approaching the linear-model value, is observed. The figure of merit (calculated for μ = μmax) shows relatively fast temperature decay due to the role of phononic thermal conductivity (see Sec. III <BR>D). <BR>We <BR>find <BR>that <BR>ZT(μmax), although being relatively small, is noticeably elevated in the presence of 
1
trigonal warping in comparison to the t= 0 case. The behavior of Smax presented in Fig. 7(c) <BR>suggests a procedure, allowing one to determine the trigonal-warping strength via directly measurable quantities. 
1
For any t= 0, one can determine a unique global maximum of S = S(μ,T), which is reached at μ = μS and T = TS .
max max
Our numerical findings for 0.1eV ≤ t1≤ 0.35 eV are presented in Fig. 8, <BR>where <BR>we <BR>have <BR>plotted <BR>(instead <BR>of <BR>μS ), the 
max
optimal effective concentration (n−p)S [see the inset]. The 
max 
125403-7 DOMINIK SUSZALSKI, GRZEGORZ RUT, AND ADAM RYCERZ PHYSICAL REVIEW B 97, 125403 (2018) 


FIG. 8. Temperature corresponding to the maximal thermopower as a function of the Lifshitz energy (data points). Least-squares fitted linear <BR>dependence <BR>[see <BR>Eq. <BR>(38)] <BR>is <BR>also <BR>displayed <BR>(line). <BR>Inset <BR>shows <BR>the carrier concentration versus ρ0kBT S , together with the model 
max
prediction <BR>(dashed <BR>line) <BR>and <BR>the <BR>linear <BR>fit <BR>(solid <BR>line) <BR>[see <BR>Eqs. <BR>(37) <BR>and <BR>(39)]. <BR>
best-fitted lines displayed in Fig. 8 <BR>are given by 
TS = 0.203(2) × EL/kB, (38)
max, fit 
(n − p)S = 0.816(5) × ρ0kBTS , (39)
max, fit max
where numbers in parentheses are standard deviations for the last digit. A few-percent deviation of the actual (n−p)S 
max 
from predictions of the linear model [see dashed line in the inset, <BR>obtained <BR>from <BR>Eq. <BR>(37) <BR>by <BR>setting <BR>τL = kBT S 
max, fit/EL ≈ 0.20] is relatively small taking into account that the existence of a global maximum of S(μ,T ) is directly linked to the breakdown of the linear model occurring for μ ∼ kBT 2 EL 
1
(and therefore is not observed in the t= 0 case). We further notice <BR>that <BR>Eqs. <BR>(38) <BR>and <BR>(39) provide direct relations between the two independent driving parameters corresponding to the optimal thermopower, TS and (n−p)S , and the trigonal 
max max
warping strength quantified by EL. 
IV. CONCLUSIONS 
We have investigated the thermopower, violation of the Wiedemann-Franz law, and the thermoelectric figure of merit, for large ballistic samples of bilayer graphene in the absence of electrostatic bias between the layers (a gapless case) and close to the charge-neutrality point. Although the thermoelectric performance is not high in such a parameter range, we find that low-temperature behavior of thermoelectric properties is determined by microscopic parameters of the tight-binding Hamiltonian, including the skew-interlayer hopping integral responsible for the trigonal warping, and by the relativistic nature of effective quasiparticles (manifesting itself in linear energy dependence of both the density of states and the electrical conductivity). 
In particular, at sub-Kelvin temperatures, clear signatures of the Lifshitz transition, having forms of anomalies in chemicalpotential dependences of the Seebeck coefficient and the Lorentz number, occurs in a vicinity of the Lifshitz energy (defined by the microscopic parameters and quantifying the trigonal-warping strength). The anomalies are blurred out by thermal excitations above the crossover temperatures (different for the two thermoelectric properties) that are directly proportional to the Lifshitz energy. 
At higher temperatures (of the order of 1 K) the trigonalwarping strength can be determined from thermoelectric measurements following one of the two different approaches: (i) finding the carrier concentration corresponding to the maximal thermopower as a function of temperature, or (ii) finding the optimal temperature, i.e., such that the thermopower reaches its global maximum. The first possibility is linked to the properties of massless quasiparticles, due to which the carrier concentration corresponding to the maximal thermopower depends approximately quadratically on temperature and reciprocally on the Lifshitz energy. On the other hand, existence of unique optimal temperature (equal to 2 K if the handbook value of the Lifshitz energy EL/kB ≈ 10 K is supposed) is related to the gradual conductivity enhancement, and subsequent suppression of the thermopower, with increasing population of thermallyexcited massive quasiparticles above the Lifshitz energy. 
To conclude, we have shown that thermoelectric measurements may complement the list of techniques allowing one to determine tight-binding parameters of bilayer-graphene Hamiltonian. <BR>Unlike <BR>the <BR>well-established <BR>techniques <BR>[34](or <BR>the other <BR>recently <BR>proposed <BR>[36,37]), they neither require high-magnetic-field measurements nor refer to conductivity scaling with the system size. Instead, the proposed singledevice thermoelectric measurements must be performed on large ballistic samples (with the length exceeding 10 μm), such that quantum-size effects define the energy scale much smaller then the Lifshitz energy. 
As we have focused on clean ballistic systems, several fac-tors which may modify thermoelectric properties of graphenebased <BR>devices, <BR>including <BR>the <BR>disorder <BR>[34], <BR>lattice <BR>defects <BR>[60], or <BR>magnetic <BR>impurities <BR>[61], <BR>are <BR>beyond <BR>the <BR>scope <BR>of <BR>this <BR>study. <BR>However, recent progress in quantum-transport experiments on ultraclean freestanding monolayer samples exceeding 1 μm size <BR>[49,62] allows us to expect that similar measurements would become possible in bilayer graphene soon. Also, as the effects we describe are predicted to appear away from the charge-neutrality point, the role of above-mentioned factors should be less significant than for phenomena appearing precisely at the charge-neutrality point, such as the minimal conductivity <BR>[63,64]. Similar reasoning may apply to the role of <BR>interaction-induced <BR>spontaneous <BR>energy <BR>gap <BR>[65–67](we notice that experimental values coincide with energy scales defined by quantum-size effects, e.g., ¯
hvF/L ≈ 3meV for L = 250 <BR>nm <BR>in <BR>Ref. <BR>[67]). <BR>
Note added. Recently, we become aware of theoretical works on strained monolayer graphene reporting quite similar, double-peak spectra of the Seebeck coefficient for sufficiently high <BR>uniaxial <BR>strains <BR>[68]. <BR>
ACKNOWLEDGMENTS 
We thank Colin Benjamin and Francesco Pellegrino for the correspondence and one of the referees for pointing out the role of the phononic part of the thermal conductance. The work was supported by the National Science Centre of Poland 
125403-8 
LIFSHITZ TRANSITION AND THERMOELECTRIC … 	PHYSICAL REVIEW B 97, 125403 (2018) 
(NCN) via Grant No. 2014/14/E/ST3/00256. Computations were partly performed using the PL-Grid infrastructure. D.S. acknowledges the financial support from dotation KNOW from Krakowskie Konsorcjum “Materia-Energia-Przyszłość” im. Mariana Smoluchowskiego. 
APPENDIX: LINEAR MODEL FOR 
TRANSMISSION-ENERGY DEPENDENCE 

At sufficiently high temperatures, thermoelectric properties given <BR>by <BR>Eqs. <BR>(16)–(21) become insensitive to the detailed functional form of T(E), and simplified models can be considered. Here we assume T(E) =C|E|, with C being a dimensionless <BR>parameter. <BR>In <BR>turn, <BR>Eq. <BR>(19) <BR>leads <BR>to <BR>
y ∞
D dx xdx 
L0 = y + 
β cosh x +1 cosh x +1
0 y 
D 
= ln (2 cosh y +2), 	(A1)
β 
y 	∞
D x 2dx xdx 
L1 =	+y
β2 0 cosh x +1 y cosh x +1 
D π2 
−y)
=+y 2 −y ln (2 cosh y +2) +4Li2(−e, 
β23 (A2) 
y ∞
D x 2dx x 3dx 
L2 = y + 
β3 0 cosh x +1 y cosh x +1 
D π23 
= y −y +y 2 ln(2 cosh y +2)
β33 
−y)
−8yLi2(−e −y) −12Li3(−e, (A3) where D =(gsgv/h) C, β =1/kBT, y =βμ, and Lis(z)is the <BR>polylogarithm <BR>function <BR>[69]. <BR>Subsequently, <BR>the <BR>Seebeck <BR>coefficient <BR>and <BR>the <BR>Lorentz <BR>number <BR>[see <BR>Eqs. <BR>(17) <BR>and <BR>(20)in the main text] are given by 
 	 
π2 −y)
kB +y 2 +4Li2(−e 
S =−y + 3 	, (A4) 
eln (2 cosh y +2)
2 	π2 −y)
Kel kB y +y 3 −12Li3(−e 
= 
TG e ln (2 cosh y +2) 
 
π2 +y 2 +4Li2(−e −y)2
− 3 	. (A5)
ln (2 cosh y +2) 
As <BR>the <BR>right-hand <BR>sides <BR>in <BR>Eqs. <BR>(A4) <BR>and <BR>(A5) depend only on a single dimensionless variable (y) they are convenient to be compared with thermoelectric properties obtained numerically via the mode-matching method (see Sec. III <BR>for details). In particular, <BR>the <BR>function <BR>of <BR>Eq. <BR>(A4) <BR>is <BR>odd <BR>and <BR>has <BR>a <BR>single <BR>maximum for y> 0, i.e., 
(S)
Smax =1.0023 kB/e for y=1.9488, (A6)
max 
which <BR>is <BR>approximated <BR>by <BR>Eq. <BR>(23) <BR>in <BR>the <BR>main <BR>text. <BR>Analo<BR>gously, <BR>the <BR>function <BR>of <BR>Eq. <BR>(A5) <BR>is <BR>even, <BR>and <BR>has <BR>a <BR>maximum <BR>at y =0, <BR>that <BR>brought <BR>us <BR>to <BR>Eq. <BR>(24) <BR>in <BR>the <BR>main <BR>text. <BR>It <BR>also <BR>reaches a minimum 
Lmin =3.0060 (kB/e)2 ≈0.91 LWF for y(L) =4.5895,
min 
(A7) 
1
with the Wiedemann-Franz constant the LWF =π2(kB/e)2.
3 
For y →∞we have L →LWF. 
[1] P. Dollfus, V. H. Nguyen, and J. Saint-Martin, J. Phys.: Condens. Matter <BR>27, <BR>133204 <BR>(2015). <BR>
[2] Y. M. Zuev, W. Chang, and P. Kim, Phys. <BR>Rev. <BR>Lett. <BR>102, <BR>096807 <BR>(2009). <BR>
[3] P. Wei, W. Bao, Y. Pu, C. N. Lau, and J. Shi, Phys. <BR>Rev. <BR>Lett. <BR>102, <BR>166808 <BR>(2009). <BR>
[4] J. G. Checkelsky and N. P. Ong, Phys. <BR>Rev. <BR>B <BR>80, <BR>081413(R) <BR>(2009). <BR>
[5] E. H. Hwang, E. Rossi, and S. Das Sarma, Phys. <BR>Rev. <BR>B <BR>80, <BR>235415 <BR>(2009). <BR>
[6] S. G. Nam, D. K. Ki, and H. J. Lee, Phys.Rev.B <BR>82, <BR>245416 <BR>(2010). <BR>
[7] C.-R. Wang, W.-S. Lu, and W.-L. Lee, 	Phys. <BR>Rev. <BR>B <BR>82, <BR>121406(R) <BR>(2010). <BR>
[8] C.-R. Wang, W.-S. Lu, L.-Hao, W.-L. Lee, T.-K. Lee, F. Lin, I.-C. Cheng, and J. Z. Chen, Phys. <BR>Rev. <BR>Lett. <BR>107, <BR>186602 <BR>(2011). <BR>
[9] M. M. Wysokiński and J. Spałek, J. <BR>Appl. <BR>Phys. <BR>113, <BR>163905 <BR>(2013). <BR>
[10] D. M. Newns, C. C. Tsuei, R. P. Huebener, P. J. M. van Bentum, 
P. C. Pattnaik, and C. C. Chi, Phys. <BR>Rev. <BR>Lett. <BR>73, <BR>1695 <BR>(1994). <BR>
[11] T. Cao, Z. Li, and S. G. Louie, Phys. <BR>Rev. <BR>Lett. <BR>114, <BR>236602 <BR>(2015). <BR>
[12] L. Seixas, A. S. Rodin, A. Carvalho, and A. H. Castro Neto, Phys. <BR>Rev. <BR>Lett. <BR>116, <BR>206803 <BR>(2016). <BR>
[13] H. Sevinçli, Nano <BR>Lett. <BR>17, <BR>2589 <BR>(2017). <BR>
[14] M. T. Pettes, I. Jo, Z. Yao, and L. Shi, Nano <BR>Lett. <BR>11, <BR>1195 <BR>(2011). <BR>
[15] A. A. Balandin, Nat. <BR>Mater. <BR>10, <BR>569 <BR>(2011). <BR>
[16] Y. Xu, Z. Li, and W. Duan, Small <BR>10, <BR>2182 <BR>(2014). <BR>
[17] W. Zhao, Y. Wang, Z. Wu, W. Wang, K. Bi, Z. Liang, J. Yang, 
Y. Chen,Z.Xu, andZ.Ni, Sci. <BR>Rep. <BR>5, <BR>11962 <BR>(2015). <BR>
[18] J. Crossno, J. K. Shi, K. Wang, X. Liu, A. Harzheim, A. Lucas, 
S. Sachdev, P. Kim, T. Taniguchi, K. Watanabe, T. A. Ohki, and 
K. C. Fong, Science <BR>351, <BR>1058 <BR>(2016). <BR>
[19] X. Zhang, Y. Gao, Y. Chen, and M. Hu, Sci. <BR>Rep. <BR>6, <BR>22011 <BR>(2016). <BR>
[20] A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, 
F. Miao, and C. Lau, Nano <BR>Lett. <BR>8, <BR>902 <BR>(2008). <BR>
[21] Y. Chen, T. Jayasekera, A. Calzolari, K. W. Kim, and M. B. Nardelli, J. <BR>Phys.: <BR>Condens. <BR>Matter <BR>22, <BR>372202 <BR>(2010). <BR>
[22] W. Huang, J.-S. Wang, and G. Liang, Phys. <BR>Rev. <BR>B <BR>84, <BR>045410 <BR>(2011). <BR>
[23] L. Liang, E. Cruz-Silva, E. C. Girão, and V. Meunier, Phys. Rev. B <BR>86, <BR>115438 <BR>(2012). <BR>
[24] H. Sevinçli, C. Sevik, T. Çağin, and G. Cuniberti, Sci. <BR>Rep. <BR>3, <BR>1228 <BR>(2013). <BR>
125403-9 
DOMINIK SUSZALSKI, GRZEGORZ RUT, AND ADAM RYCERZ 	PHYSICAL REVIEW B 97, 125403 (2018) 
[25] M. S. Hossain, F. Al-Dirini, F. M. Hossain, and E. Skafidas, Sci. Rep. <BR>5, <BR>11297 <BR>(2015). <BR>
[26] Y. Anno, Y. Imakita, K. Takei, S. Akita, and T. Arie, 2D <BR>Mater. <BR>4, <BR>025019 <BR>(2017). <BR>
[27] X. Chen, L. Zhang, and H. Guo, Phys.Rev.B <BR>92, <BR>155427 <BR>(2015). <BR>
[28] E. McCann, Phys.Rev.B <BR>74, <BR>161403(R) <BR>(2006). <BR>
[29] H. K. Min, B. Sahu, S. K. Banerjee, and A. H. MacDonald, Phys. Rev. <BR>B <BR>75, <BR>155115 <BR>(2007). <BR>
[30] Y. Zhang, T.-T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, 
M. F. Crommie, Y. R. Shen, and F. Wang, Nature <BR>(London) <BR>459, <BR>820 <BR>(2009). <BR>
[31] E. McCann and V. 	I. Fal’ko, Phys. <BR>Rev. <BR>Lett. <BR>96, <BR>086805 <BR>(2006). <BR>
[32] M. Orlita, P. Neugebauer, C. Faugeras, A.-L. Barra, M. Potemski, 
F. M. D. Pellegrino, and D. M. Basko, Phys.Rev.Lett. <BR>108, <BR>017602 <BR>(2012). <BR>
[33] M. I. Katsnelson, Graphene: Carbon in Two Dimensions (Cambridge University Press, Cambridge, 2012). 
[34] E. McCann and M. Koshino, Rep. <BR>Prog. <BR>Phys. <BR>76, <BR>056503 <BR>(2013). <BR>
[35] A. G. Moghaddam and M. Zareyan, Phys. <BR>Rev. <BR>B <BR>79, <BR>073401 <BR>(2009). <BR>
[36] G. Rut and A. Rycerz, Europhys. <BR>Lett. <BR>107, <BR>47005 <BR>(2014). <BR>
[37] G. Rut and A. Rycerz, Phys. <BR>Rev. <BR>B <BR>93, <BR>075419 <BR>(2016). <BR>
[38] M. 	Paulsson and S. Datta, Phys.Rev.B <BR>67, <BR>241403(R) <BR>(2003). <BR>
[39] K. Esfarjani, M. Zebarjadi, and Y. Kawazoe, Phys. <BR>Rev. <BR>B <BR>73, <BR>085406 <BR>(2006). <BR>
[40] See e.g. M. S. Dresselhaus and G. Dresselhaus, Adv. <BR>Phys. <BR>30, <BR>139 <BR>(1981); <BR>reprinted <BR>in <BR>51, <BR>1 <BR>(2002). <BR>
[41] The values of t0 and t⊥ are taken from: A. B. Kuzmenko, I. Crassee, D. van der Marel, P. Blake, and K. S. Novoselov, Phys. Rev. <BR>B <BR>80, <BR>165406 <BR>(2009). <BR>
[42] The combined particle-hole-reflection symmetry applies as we have neglected the next-nearest-neighbor intralayer hopping t2, which <BR>is <BR>not <BR>determined <BR>for <BR>BLG <BR>[34]. <BR>Effects <BR>of t2 ∼ t⊥ are insignificant when discussing the band structure for |E|∼EL «t⊥. 
[43] We notice a generic relation between ρ(E) and the cyclotronic mass mC(E) =(πh¯2/2)ρ(E). 
[44] R. Landauer, IBM <BR>J. <BR>Res. <BR>Dev. <BR>1, <BR>223 <BR>(1957). <BR>
[45] M. Buttiker, Y. Imry, R. Landauer, and S. Pinhas, Phys.Rev.B <BR>31, <BR>6207 <BR>(1985); <BR>M. <BR>Buttiker, <BR>Phys. <BR>Rev. <BR>Lett. <BR>57, <BR>1761 <BR>(1986); <BR>IBM <BR>J. <BR>Res. <BR>Dev. <BR>32, <BR>317 <BR>(1988). <BR>
[46] For BLG, contributions of out-of-plane phonons, governing the thermal properties at low-temperatures, are additionally reduced in comparison to MLG, see: A. I. Cocemasov, D. L. Nika, and 
A. A. Balandin, Nanoscale <BR>7, <BR>12851 <BR>(2015); <BR>D. <BR>L. <BR>Nika <BR>and <BR>A. <BR>
A. Balandin, J. <BR>Phys.: <BR>Condens. <BR>Matter <BR>24, <BR>233203 <BR>(2012); <BR>N. <BR>Mingo and D. A. Broido, Phys.Rev.Lett. <BR>95, <BR>096105 <BR>(2005). <BR>
[47] J. Tworzydło, B. Trauzettel, M. Titov, A. Rycerz, and C. W. J. Beenakker, Phys.Rev.Lett. <BR>96, <BR>246802 <BR>(2006). <BR>
[48] R. Danneau, F. Wu, M. F. Craciun, S. Russo, M. Y. Tomi, J. Salmilehto, A. F. Morpurgo, and P. J. Hakonen, Phys. <BR>Rev. <BR>Lett. <BR>100, <BR>196802 <BR>(2008). <BR>
[49] M. Kumar, A. Laitinen, and P. J. Hakonen, arXiv:1611.02742. <BR>
[50] In some geometries, the transmission is suppressed by approximately constant factor for E< 0 due to contact effects; see supplemental <BR>information <BR>in <BR>Ref. <BR>[49]. <BR>
[51] R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, Science <BR>320, <BR>1308 <BR>(2008). <BR>
[52] I. Snyman and C. W. J. Beenakker, Phys.Rev.B <BR>75, <BR>045322 <BR>(2007). <BR>
[53] Typically, using double-precision arithmetic, we put Ndiv =20 for L =W/20 =104l⊥ and |μ|≤2 meV. Numbers of different transverse momenta ky =2πq/W,where q =0,±1,± 2,...,±qmax, necessary to determine Tr(tt†)inEq. <BR>(25) <BR>with <BR>a 10-digit accuracy were found to be 1955 ≤2qmax +1 ≤8149, with the lower (upper) value corresponding to μ =0and t1= 0.1eV (μ =2 meV and t1=0.35 eV). 
[54] For L =W/20 =104l⊥, effects of the crystallographic orientation are negligible starting from |μ|2 ¯
hvF/L ≈38μeV, as the summation <BR>over <BR>normal <BR>modes <BR>in <BR>Eq. <BR>(25) <BR>involves <BR>numerous <BR>contributions, corresponding to different angles of incident, that are comparable. 
[55] G. Rut and A. Rycerz, Acta <BR>Phys. <BR>Polon. <BR>A <BR>126, <BR>A-114 <BR>(2014). <BR>
[56] Due to the presence of Fabry-Pérot oscillations, the well-known 
−1k2
Mott formula S =(π2/3)e BT ×[∂ln T(E)/∂E]E=μ cannot be directly applied. 
[57] A. Alofi and G. 	P. Srivastava, Phys.Rev.B <BR>87, <BR>115421 <BR>(2013); <BR>A. <BR>Alofi, <BR>Theory <BR>of Phonon Thermal Transport in Graphene and Graphite, Ph.D. Thesis, University of Exeter, 2014; http://hdl.handle.net/10871/15687. <BR>
[58] J. Callaway, Phys. <BR>Rev. <BR>113, <BR>1046 <BR>(1959). <BR>
[59] The model parameters are choosen such that the experimental results <BR>of <BR>Ref. <BR>[14] <BR>are <BR>reproduced <BR>if <BR>L =0.25 μm. Namely, the Callaway model parameters are BU =4.77 ×10−25 sK−3 and BN =3.18 ×10−25 sK−3, the defect concentration is cd = 
1.6 ×10−4, and the concentration of C13 isotope is 10−2.Remaining parameters of the phonon dispersion relations are same as <BR>in <BR>the <BR>first <BR>paper <BR>of <BR>Ref. <BR>[57]. <BR>
[60] M. S. Dresselhaus, A. Jorio, A. G. Souza Filho, and R. Saito, Phil. <BR>Trans. <BR>R. <BR>Soc. <BR>A <BR>368, <BR>5355 <BR>(2010). <BR>
[61] B. Uchoa, T. G. Rappoport, and A. H. Castro Neto, Phys. Rev. Lett. <BR>106, <BR>016801 <BR>(2011); <BR>J. <BR>Hong, <BR>E. <BR>Bekyarova, <BR>W. <BR>A. <BR>de <BR>Heer, <BR>
R. C. Haddon, and S. Khizroev, ACS <BR>Nano <BR>7, <BR>10011 <BR>(2013). <BR>
[62] P. 	Rickhaus, P. Makk, M.-H. Liu, E. Tóvári, M. Weiss, R. Maurand, K. Richter, and C. Schönenberger, Nat. <BR>Commun. <BR>6, <BR>6470 <BR>(2015). <BR>
[63] A. S. Mayorov, 	D. C. Elias, M. Mucha-Kruczynski, R. V. Gorbachev, T. Tudorovskiy, A. Zhukov, S. V. Morozov, M. I. Katsnelson, V. I. Fal’ko, A. K. Geim, and K. S. Novoselov, Science <BR>333, <BR>860 <BR>(2011). <BR>
[64] S. Samaddar, I. Yudhistira, S. Adam, H. Courtois, and C. B. Winkelmann, Phys. <BR>Rev. <BR>Lett. <BR>116, <BR>126804 <BR>(2016). <BR>
[65] G. M. Rutter, S. 	Jung, N. N. Klimov, D. B. Newell, N. B. Zhitenev, and J. A. Stroscio, Nat. <BR>Phys. <BR>7, <BR>649 <BR>(2011). <BR>
[66] F. Freitag, J. Trbovic, M. Weiss, and C. Schönenberger, Phys. Rev. <BR>Lett. <BR>108, <BR>076602 <BR>(2012). <BR>
[67] W. 	Bao, J. Velasco, Jr., F. Zhang, L. Jing, B. Standley, D. Smirnov, M. Bockrath, A. H. MacDonald, and C. N. Lau, Proc. <BR>Natl. <BR>Acad. <BR>Sci. <BR>USA <BR>109, <BR>10802 <BR>(2012). <BR>
[68] A. Mani and C. Benjamin, Phys. <BR>Rev. <BR>E <BR>96, <BR>032118 <BR>(2017); <BR>arXiv:1707.07159. <BR>
[69] K. Oldham, J. Myland, and J. Spanier, An Atlas of Functions (Springer-Verlag, New York, 2009), Chap. 25. 
125403-10 


Journal of Physics: Condensed Matter 
J. Phys.: Condens. Matter 31 (2019) 415501 (8pp) https://doi.org/10.1088/1361-648X/ab2d0c 
Thermoelectric properties of gapped bilayer graphene 
Dominik Suszalski, Grzegorz Rut and Adam Rycerz 
Marian Smoluchowski Institute of Physics, Jagiellonian University, Łojasiewicza 11, PL–30348 Kraków, Poland 
E-mail: rycerz@th.if.uj.edu.pl 
Received 15 April 2019, revised 14 June 2019 Accepted for publication 26 June 2019 Published 15 July 2019 
Abstract 
Unlike in conventional semiconductors, both the chemical potential and the band gap in bilayer graphene (BLG) can be tuned via application of external electric field. Among numerous device implications, this property also designates BLG as a candidate for high-performance thermoelectric material. In this theoretical study we have calculated the Seebeck coefficients for abrupt interface separating weakly- and heavily-doped areas in BLG, and for a more realistic rectangular sample of mesoscopic size, contacted by two electrodes. For a given band gap (Δ) and temperature (T) the maximal Seebeck coefficient is close to the Goldsmid– Sharp value |S|GS =Δ/(2eT), the deviations can be approximated by the asymptotic 
max 
expression |S|maxGS −|S|max =(kB/e) ×[12 ln u + ln 2 −12 + O(u−1)], with the electron charge −e, the Boltzmann constant kB, and u =Δ/(2kBT)» 1. Surprisingly, the effects of trigonal warping term in the BLG low-energy Hamiltonian are clearly visible at few-Kelvin temperatures, for all accessible values of Δ : 300 meV. We also show that thermoelectric figure of merit is noticeably enhanced (ZT > 3) when a rigid substrate suppresses out-of-plane vibrations, reducing the contribution from ZA phonons to the thermal conductivity. 
Keywords: graphene, bilayer graphene, thermoelectrics, Seebeck effect, ballistic transport 

S Supplementary material for this article is available online 
(Some figures may appear in colour only in the online journal) 
1. Introduction 
In recent years, bilayer graphene (BLG) devices made it possible to demonstrate several intriguing physical phenomena, including the emergence of quantum spin Hall phase [1, 2], the fractal energy spectrum known as Hofstadter’s butterfly [3, 4], or unconventional superconductivity [5, 6], just to mention a few. From a bit more practical perspective, a number of plasmonic and photonic instruments were designed and build [7–9] constituting platforms for application considerations. BLG-based thermoelectric devices have also attracted a significant attention [10–13], next to the devices based on other two-dimensional (2D) materials [14–16]. 
In a search for high-performance thermoelectric material, one’s attention usually focusses on enhancing the dimensionless figure of merit [17, 18] 
GS2T
ZT = , (1)
K 
where G, S and K are (respectively): the electrical conductance, the Seebeck coefficient quantifying the thermopower, and the thermal conductance; the last characteristic can be represented as K = Kel + Kph, with Kel (Kph) being the electronic (phononic) part. This is because the maximal energy conversion efficiency is related to ZT via [19] 
√ 
ΔT 1 + ZTav − 1 
ηmax = √ ,
Tc (2)
Th1 + ZTav + 
Th 
where Tc (Th) is the hot- (or cool) side temperature, ΔT = Th − Tc , and Tav =(Tc +Th)/2. In particular, for 
1
ZTav = 3 we have ηmax > 3 ΔT/Th (with ΔT/Th the Carnot efficiency), and therefore ZT > 3 is usually regarded as a 
1361-648X/19/415501+8$33.00 1 © 2019 IOP Publishing Ltd Printed in the UK J. Phys.: Condens. Matter 31 (2019) 415501 
condition for thermoelectric device to be competitive with other power generation systems. 
As the Seebeck coefficient is squared in equation (1), a maximum of ZT—considered as a function of the driving parameters specified below—is commonly expected to appear close to a maximum of |S|. Here we show this is not always the case in gapped BLG: when discussing in-plane thermoelectric transport in a presence of perpendicular electric field such that the band gap is much greater than the energy of thermal excitations (Δ » kBT), the maximal absolute thermopower |S|max corresponds to the electrochemical potential relatively close 
|S|
to the center of a gap, namely µmax ≈±1 kBT ln(2Δ/kBT),
2 
whereas the maximal figure of merit ZTmax appears near the maximum of the valence band (or the minimum of the con-
ZT
duction band), i.e. µ≈±Δ/2. In contrast to ZTmax, |S|max
max 
is not directly related to the value of the transmission probability near the band boundary, and these two quantities show strikingly different behaviors with increasing Δ for a given T. 
Qualitatively, one can expect that thermoelectric performance of BLG is enhanced with increasing Δ, since abrupt switching behavior is predicted for the conductance G when passing µ = ±Δ/2 for Δ » kBT [20, 21]. (Moreover, the gap opening suppresses κel, reducing the denominator in equation (1).) Such a common intuition is build on the Mott formula, according to which S is proportional to the logarithmic derivative of G as a function of the Fermi energy EF. One cannot, however, directly apply the Mott formula for gapped systems at nonzero temperatures, and the link between a rapid increase of G(EF) for EF ≈ Δ/2 and high |S| is thus not 
|S|
direct in the case of gapped BLG, resulting in µmax « Δ/2. 
The results of earlier numerical work [22] suggest that |S|max, obtained by adjusting µ for a given Δ and T, is close to 
Δ 
|S|GS = , (3)
max 
2eT 
being the Goldsmid–Sharp value for wide-gap semiconductors [23]. In this paper, we employ the Landauer–Büttiker approach for relatively large ballistic BLG samples, finding that equation (3) provides a reasonable approximation of the actual |S|max for Δ ∼ kBT only. For larger Δ, a logarithmic correction becomes significant, and the deviation exceeds kB/e (≈86 µV K−1) for Δ 2 10 kBT. We further find that— although |S|max grows monotonically when increasing Δ at fixed T and may reach (in principle) arbitrarily large value— ZTmax shows a conditional maximum at Δ*(T) ∼ 102 kBT (for T : 10 K). An explanation of these findings in terms of a simplified model for transmission-energy dependence is provided. 
2. Model and methods 
The two systems considered are shown schematically in figures 1(a) and (b). The first system (hereinafter called an abrupt interface) clearly represents an idealized case, as we have supposed that both the doping and temperature change rapidly on the length-scale much smaller than the Fermi wavelength for an electron. Therefore, a comparison with the second system, in which chemical potentials and temperatures are attributed to macroscopic reservoirs (the two leads), separated by a sample area of a finite length L, is essential to validate the applicability of our precictions for real experiments. 
We take the four-band Hamiltonian for BLG [24] 
−U/2 vFπξt⊥ 0
⎛⎞ 
vFπ† −U/20 v3π
H = ξ ,
⎜
ξt⊥ 0 U/2 vFπ†⎟
(4)
⎝⎠ 
0 v3π† vFπ U/2 
where the valley index ξ = 1 for K valley or ξ = −1 for K' valley, π = px + ipy, π† = px − ipy, with p =(px, py) the 
√ 
carrier momentum, vF = 3at0/(2!) is the Fermi velocity, v3 =(t'/t0)vF , U is the electrostatic bias between the layers, and a = 0.246 nm is the lattice parameter. Following [25], we set t0 = 3.16 eV—the nearest-neighbor in-plane hopping energy, t⊥ = 0.381 eV—the direct interlayer hopping energy; the skew interlayer hopping energy is set as t' = 0 or t' = 0.3 eV in order to discuss the role of trigonal warping1. The band gap Δ ≈|U|for |U|«t⊥ and t' = 0 (remaining details are given in supplementary information, section I (stacks.iop.org/ JPhysCM/31/415501/mmedia)). Solutions of the subsequent Dirac equation, HΨ=EΨ, with Ψ=(ΨA1, ΨB1, ΨB2, ΨA2)T the probability amplitudes, are matched for the interfaces separating weakly- and heavily-doped regions allowing us to determine the energy-dependent transmission probability T(E) (see supplementary information, section II). 
At zero temperature, the relation between physical carrier concentration (the doping) and the Fermi energy EF can be approximated by a close-form expression for t' = 0 and |EF|−Δ/2 «t⊥, namely 
t⊥ 
n(EF) ≈ max (0, |EF|−Δ/2) , (5)
π(lvF)2 
following from a piecewise-constant density of states in such a parameter range [24]. The prefactor in equation (5) 
−2
(1/π) t⊥/(!vF)2 ≈ 0.268 nm eV−1. In a general situation, it is necessary to perform the numerical integration of the density of states following from the Hamiltonian H (4) (see supplementary information, section I); however, equation (5) gives a correct order of magnitude for Δ » kBT. In the remaining part of the paper, we discuss thermoelectric characteristics as functions of the electrochemical potential µ, keeping in mind that the doping n = n(|µ|) is a monotonically increasing function of |µ|. 
Next, we employ the Landauer–Büttiker expressions for the electrical and thermal currents [26, 27] 
gsgve
I = −J dET(E)[ fL(E) − fR(E)] , (6)
h 
gsgv
IQ = J dET(E)[ fL(E) − fR(E)] (E− µ), (7)
h 
1 Remaining parameters of the systems studied numerically are, for an abrupt interface of figure 1(a): W = 103 l⊥ = 1.77 µm (with l⊥ = nvF/t⊥); and for a rectangular setup of figure 1(b): L = W/20 = 104 l⊥. 
2 J. Phys.: Condens. Matter 31 (2019) 415501 

Figure 1. Systems studied numerically in the paper (left) and their basic thermoelectric characteristics (right). (a) BLG strip of width W with a voltage source driving electric current through the strip. Two gate electrodes, with different temperatures T1 and T2, induce an abrupt interface (dash-dot line) between the weakly-doped (light area) and the heavily-doped (dark area) regions. (b) A finite, weaklydoped section of the strip (of length L) with leads driving both thermal and electric currents. In both cases, additional gate electrodes (not shown) allow to tune the electrostatic bias between the layers (See footnote 1.). The coordinate system is also shown. (c) Zero-temperature 
conductance, specified in the units of g0 = 4e2/h, and (d) the Seebeck coefficient at T = 5 K as functions of the chemical potential. Solid lines correspond to the system of panel (a), dashed lines correspond to the system of panel (b). (e) Zoom-in of the conductance as a function of the chemical potential measures from the conduction band minimum (µ =Δ/2) for the system of panel (a). Skew-interlayer hopping is fixed at t1 = 0.3 eV. The electrostatic bias between the layers is varied from U = 0 to U = 20 meV in steps of 5 meV (c) and (d) or from U = 10 meV to U = 150 meV in steps of 20 meV (e). (The extreme values of U are specified for corresponding lines.) 
where gs = gv = 2 are spin and valley degeneracies, fL and fR are the distribution functions for the left and right reservoirs, with their electrochemical potentials µL and µR, and temperatures TL and TR. We further suppose that µL −µR ≡−eV and TL − TR ≡ ΔT are infinitesimally small (the linear-response regime) and define µ =(µL +µR)/2 and T =(TL +TR)/2. The conductance G, the Seebeck coefficient S, and the electronic part of the thermal conductance Kel, are given by [28] 
G = I = e2L0, (8)
V    ΔT=0 
VL1
S = − = , (9)
ΔT    I=0 eTL0 
IQ L0L2 − L2 
Kel == 1, (10)
ΔT    I=0 TL0 
where 
Ln = gsgv J dET(E) (− ∂fFD ) (E− µ)n 
h ∂E (11) (n= 0, 1, 2) 
with fFD(µ, T, E)= 1/ [ exp ((E − µ)/kBT)+ 1 ] being the Fermi–Dirac distribution function. In particular, for T → 0, equation (8) reduces to G =(gsgve2/h)T(µ), the well-known zero-temperature Landauer conductance. 
The phononic part of the thermal conductance, occuring in equation (1), can be calculated using 
1 ∂fBE
Kph = J dω lω Tph(ω), (12)
2π∂T 
with fBE(T, ω)= 1/ [ exp (!ω/kBT) − 1 ] the Bose–Einstein distribution function and Tph(ω) the phononic transmission spectrum. We calculate Tph(ω) by adopting the procedure developed by Alofi and Srivastava [29] to the two systems considered in this work (see supplementary information, section IV). 
3. Numerical results 
Before discussing the thermoelectric properties in details, we present zero-temperature conductance spectra, which represent the input data to calculate thermoelectric properties (see section 2). Since the Hamiltonian given by equation (4) is 
3 J. Phys.: Condens. Matter 31 (2019) 415501 

Figure 2. (a) Maximal absolute value of the Seebeck coefficient and (b), (c) its deviation from the Goldsmid–Sharp value 
|S|GS 
=Δ/(2eT)calculated numerically for the system of
max 
figure 1(a) (datapoints) as functions of the energy gap for different temperatures T = 1 K (open circles), T = 5 K (full circles), and T = 10 K (triangles). Skew-interlayer hopping integral is t1 = 0.3 
= SGS
eV (a) and (b) or t1 = 0 (c). Lines depict Smax (black
max
dotted line), and predictions of the model for transmission-energy dependence given by equation (14) with α = 0 (blue dashed line), α = 1 (red solid line), and α = 2 (magenta dash-dot line—omitted in panel (a) for clarity). 
particle-hole symmetric, it is sufficient to limit the discussion to µ � 0. 
Typically, the conductance G(µ) of the finite-strip section, compared with the case of an abrupt interface, is reduced by approximately 50% near the conduction band minimum (µ =Δ/2) due to backscattering on the second interface, and slowly approaches the abrupt-interface limit for µ » Δ/2; see figure 1(c). Additionally, for a rectangular setup of figure 1(b) oscillations of the Fabry–Pérot type are well-pronounces, see [13]. In contrast, the Seebeck coefficient is almost identical for both systems, see figure 1(d). We further notice that the conductance near µ =Δ/2, displayed in figure 1(e), is gradually suppressed with increasing U. (Hereinafter, the bandgap Δ is determined numerically for the dispersion relation following from equation (4), see supplementary information, section I. In general, Δ < |U| [24]). 

1 
T(α)  δ(E −12 Δ) + δ(E + 2 Δ) if α= 0 
= C(Δ) ×Θ(|E|−12 Δ)(|E|−12 Δ)α−1 if α>0 
(13) with the prefactor C(Δ) quantifying the transmission propability near the band boundary |E|≈12 Δ, δ(x) being the Dirac delta function, and Θ(x) being the Heaviside step function. The analytic expressions presented here, and later in section 4, are (unless otherwise specified) valid for any α � 0, althought when comparing model predictions with the numerical data we limit our considerations to integer α. 
In particular, equation (13) leads to the maximal absolute value of the Seebeck coefficient 
−1 √ 
|S|(α) × (kB/e)≈ J(u + α)(u + α − 1) − ln (√ u + α + u + α − 1)
max 
1 13 − 7α + α2 
= u − ln (4u)+ α − ++ O(u−2) for α � 0, (14)
22 
The close overlap of the thermopower spectra presented in figure 1(d) allows us to limit the discussion of |S|max to the case of an abrupt interface, see figures 2(a)–(c). In order to rationalize the deviations of the numerical data from |S|GS 
max 
given by equation (3), we propose a family of models for the transmission-energy dependence, namely 
8u 
where the first asymptotic equality corresponds to u =Δ/(2kBT)» 1 (see supplementary information). It is 
(α)
clear from figure 2(a) that |S|max with α = 1 (red solid line) reproduces the actual numerical results (datapoints) notice
(α)
ably better than |S|max with α = 0 (blue dashed line) or |S|GS 
max 
(black dotted line). What is more, the deviations |S|GS −|S|
max max, 
4 J. Phys.: Condens. Matter 31 (2019) 415501 

Figure 4. Left: thermoelectric figure of merit displayed as a function of the chemical potential for temperatures (a) T = 1 K and (b) T = 5 
K. Solid lines correspond for the system of figure 1(a), dashed lines correspond for the system of figure 1(b). (Remaining parameters are same as in figures 1(c) and (d).) Right: maximal value of the figure of merit for the system of figure 1(a), versus the energy gap for temperatures (c) T = 1 K and (d) T = 5 K. Inset in panel (c) shows the optimal gap Δ* as a function of temperature. Different lines in panels (c) and (d) correspond to the limit of rigid substrate eliminating ZA phonons (solid lines) or the free-standing sample (dashed-dotted 
lines). 
displayed as functions of Δ, allows one to easily identify the effects of trigonal warping, see figures 2(b) for t' = 0.3 eV and 2(c) for t' = 0. 
Next, we investigate the figure of merit (ZT) given by equation (1). For this purpose, it is necessary to calculate both the electronic part of thermal conductance (Kel), that is determined by the energy-dependent transmission T(E) (see section 2), as well as the phononic part (Kph), presented in figure 3. We find that for T < 10 K the two systems considered show almost equal Kph (∝ T3), and different values of ZT (see figures 4(a) and (b)) follow predominantly from the G reduction discussed above. Unlike |S|max, a value of which (for a given T) is limited only by the largest experimentally-accessible Δ ≈ 300 meV [24, 30], ZTmax shows well-defined conditional maximum for the optimal bandgap Δ*/kBT ≈ 100–150 (at the temperature range 10 K� T � 1 K) and decreases for Δ > Δ*; see figures 4(c) and (d). 
Also in figures 4(c) and (d), we compare ZTmax for free-standing BLG, in which all polarizations of phonons (LA, TA, and ZA) contribute to the thermal conductance (see supplementary information, section IV) with an idealized case of BLG on a rigid substrate, eliminating out-of-plane (ZA) phonons. In the latter case, ZTmax is amplified, approximately by a factor of 3 (for any Δ), exceeding ZT = 3 for T = 1 K and Δ ≈ Δ* = 10 meV. 
4. Discussion 
Let us now discuss here why we have identified apparently different behaviors of |S|max and ZTmax with increasing Δ. To understand these observations, we refer to the model T(α)(E) 
(α)
given by equation (13) with α � 0, for which |S|max, approximated by equation (14), corresponds to 
max √
µ|S|≈ ln (√ u + α + u + α − 1)
kBT 
(15)
1 ( 2Δ )
≈ ln for Δ » kBT.
2 kBT 
In contrast, the chemical potential corresponding the the maximal ZT is much higher and can be approximated (in the Δ » kBT limit) by 
or 
Δ 
µ ZT ≈− 1.145 kBT for α = 1, (16)
max 
2 
Δ 
µ ZT ≈ + 0.668 kBT for α = 2, (17)
max 
2 
where we have further supposed that Kph » Kel, being equivalent to 
5 J. Phys.: Condens. Matter 31 (2019) 415501 

Figure 5. Chemical potential corresponding to the maximal 
|S|
absolute thermopower presented in figure 2 (µmax) and the maximal figure of merit presented in figure 4 (µZT ) as function of the energy 
max
gap. Points correspond to same datasets as in figures 2(a) and (b) (skew-interlayer hopping integral is t1 = 0.3 eV). Lines depict approximating equation (15) with α = 1 (red solid line) and α = 2 (red dash-dot line), equations (16) (blue solid line), and (17) (blue dash-dot line). 
ZT ≈ GS2 T . (18)Kph(T) 
As the last term in equation (18) depends only on T we can focus now on the power factor (GS2), a maximal value of which can be approximated by 
gsgvkB2 α−1
(GS2)≈M(α) (kBT)C(Δ), (19)
maxmax 
h 
(α)
where the prefactor Mmax depends only on α and is equal to 
1.27 for α = 1 or to 4.06 for α = 2 (see also supplementary information, section III). 
The positions of maxima visible in figures 4(a) and (b) are numerically close to the approximation given by equation (16). Also, the data visualized in figure 5, together with the lines corresponding to equations (16) and (17), further 
support our conjucture that the model T(α)(E), given by equation (13) with α = 1 (the step-function model), is capable of reproducing basic themoelectric characteristics of gapped BLG with a reasonable accuracy. 
Although it might seem surprising at first glance that the 
step-function model (α = 1) reproduces the actual numerical results much better than T(α)(E) with α = 2, as the conductance spectra in figure 1(e) exhibit linear, rather than step-line, energy dependence for µ � Δ/2. However, for the temperature range of 1 K : T : 10 K the T(E) behavior for E − Δ/2 » 0.1 eV, visualized with the data in figure 1(c), preveils (notice that the full width at half maximum for −∂fFD/∂E in equation (11) is ≈3.53 kBT). For this reason, a simple step-function model grasps the essential features of the actual T(E). 
|S|
ZT
Apart from pointing out that µmax « µmax ≈ Δ/2 for Δ » kBT (some further implications of this fact are discussed below), the analysis starting from T(α)(E) models also leads to the conclusion that—unlike |S|max that is not directly related to C(Δ)—for the figure of merit we have: ZTmax ∝(GS2)max ∝C(Δ) (see equations (18) and (19)). It becomes clear now that a striking ZTmax suppression for large Δ is directly link to the local G suppression for large U, illustrated in figure 1(e). Power-law fits to the datasets presented in figures 4(c) and (d), of the form ZTmax ∝ Δ−γ for Δ > Δ*, lead to γ ≈ 0.5. It is worth to stress here that the dispersion relation, and also the number of open channels as a function of energy above the band boundary, Nopen(|E|−Δ/2) (see supplementary information, section II), is virtually unaffected by the increasing Δ. Therefore, the average transmission for an open channel near the band boundary decreases with Δ. This observation can be qualitatively understood by pointing out a peculiar (Mexican hat-like) shape of the dispersion relation for Δ > 0 [31]. In the energy range Δ/2 < |E| < |U|/2 there is a continuous crossover from zero transmission (occuring for |E| < Δ/2) to a high-transmission range (|E| > |U|/2). As the width of such a crossover energy range, (|U|−Δ)/2, increases monotonically with Δ, the continuity of T(E) implies that the average transmission near |E|≈Δ/2 decreases with Δ. 
For a bit more formal explanation, we need to refer the total transmission probability through an abrupt interface (see supplementary information, section II). For the incident wavefunction with the momentum parallel to the barier (conserved during the scattering) Mky (where ky = 2πq/W and q = 0, ±1, ±2, ... assuming the periodic boundary conditions) we have 
Tky = |tnm|2 jx (ψIIn ) /jx (ψIm) , 
(20) 
m,n 
where {tm} is the 2 × 2 transmission matrix (m, n = 1, 2 are
n 
the subband indices) to be determined via the mode-matching, and jx(ψXn ) is the x-component of electric current for the wave function propagating in the direction of incidence, with X = I, II indicating the side of a barrier. For E =E +Δ/2 (with 0 <E« Δ/2) there are propagating modes in a weaklydoped area (X = II) with n = 1 (the lower subband) only; in the simplest case without trigonal warping (t1 = 0) we find that the relevant current occuring in equation (20) scales as 
jx (ψII 1 ) ∝ v;/Δ. (21) 
The above allows us to expect that the full transmission scales as T(E)=Lky Tky ∝ Vr/Δ also for t1
= 0, provided that the band gap is sufficiently large (Δ » EL , with 
1
EL = 4 t⊥ (t1/t0) 2 the Lifshitz energy). As the number of propagating modes is approximately Δ-independent, scaling 
√
roughly as Nopen ∝ E, we can further predict that zero- (or low-) temperature conductance should follow the approximate scaling law 
6 J. Phys.: Condens. Matter 31 (2019) 415501 
√ 
G(µ)∝ (µ − Δ/2)/ Δ(for µ � Δ/2). (22) 
This expectation is further supported with the numerical data presented in figure 1(e). 
An initial increase of ZTmax with Δ for 0 < Δ � Δ�, also apparent in figures 4(c) and (d), can be understood by pointing out that the electronic and phononic parts of the thermal conductance are of the same order of magnitude (Kel ∼ Kph ) in such a range. In consequence, an upper bound to ZT can be written as (up to the order of magnitude) ZT � TGS2/Kel, allowing a rapid increase of ZT with Δ (see supplementary information, section III), until Kel (decreasing with Δ) is overruled by Kph (Δ-independent). 
In the remaining part of this section, we briefly discuss the possible influence of electron–electron interactions, neglected in our numerical analysis. 
Several experimental works on free-standing BLG report an intrinsic (or spontaneous) band gap of Δint(T = 0) ≈ 1.5 meV vanishing above the critical temperature Tcrit ≈ 12 K [32–34]. To the contrary, no signatures of an intrinsic band gap are reported for BLG in van der Waals heterostructures (VDWHs) [35], in which thermoelectric properties may be significantly enhanced due to the suppression of out-of-plane (ZA) vibrations. Possibly, the above-mentioned difference could be attributed to a modification of the effective dielectric constant due to the materials surrounding a BLG sample in VDWHs. In fact, a basic mean-field description, relating Δint > 0 to the alternating spin order, allows one to expect that Δint ∼ t0 exp(−const. × t0/Ueff) (where const. ∼1 is determined by the bandwith), and thus a moderate decrease of the effective Hubbard repulsion (Ueff ) strongly suppresses Δint; see [36]. 
Although temperatures considered in this paper (0 < T � 10 K) are essentially lower then Tcrit, we focus on the case with a bias between the layers |U|≈Δ ∼10–100 meV, and much smaller Δint should not affect the physical properties under consideration. 
Additionally, the maximal ZT appears near the bottom of the conduction band or the top of the valence band 
ZT
(|µ|≈Δ/2), where one of the layers is close to the charge
max
neutrality, and thus one can expect the Coulomb-drag effects to be insignificant [37]. 
5. Concluding remarks 
We have numerically investigated thermoelectric propeties of large ballistic samples of electrostatically-gapped BLG. A logarithmic deviation of the maximal absolute thermopower from the Goldsmid–Sharp relation is identified and rationalized with the help of the step-function model for the transmission-energy dependence. In addition to the earlier findings that the trigonal warping term modifies the density of states 
[24] and transport properties [13] also for Fermi energies � 1 meV, we show here that signatures of trigonal warping may still be visible in thermoelectric characteristics for the band gaps as large as Δ ∼ 100 meV. 
Next, the analysis is supplemented by determining the total 
(i.e. electronic and phononic) thermal conductance, making it possible to calculate the dimensionless figure of merit (ZT). The behavior of maximal ZT with the increasing gap can also be interpreted in terms of the step-function model, provided that we supplement the model with the scaling rule for the typical transmision probability for an open channel near the minimum of the conduction band (or the maximum of the valence band), which is ∝ Δ−0.5 (for large Δ). This can be attributed to the Mexican-hat like shape of the dispersion relation. 
Although some other two-dimensional systems with the Mexican-hat like (or quartic) [15] dispersion also show enhanced thermoelectric properties, two unique features of BLG are worth to stress: (i) the possibility of tuning both the chemical potential and the band gap in a wide range, and 
(ii) the ballistic scaling behavior of transport characteristics with a barrier width. It is also worth to stress that the value of ZT > 1 (or ZT > 3 in the absence of out-of-plane vibrations) is reached at T = 1 K for a moderate bandgap Δ ≈ 10 meV, 
−1
Acorreponding to the electric field of ≈3 mV ˚, making it 
possible to consider an experimental setup in which the spontaneous electric field from a specially-choosen substrate (e.g. ferroelectric) is employed to reduce the energy-consumption in comparison to a standard dual-gated setup [30]. 
High values of ZT at a 1 K temperature may not have practical device implications per se, however, we hope that scaling mechanisms identified in our work will help to find the best thermoelectric among graphene-based (and related) systems. The necessity to reduce the phononic part of the thermal conductance with a simultaneous increase of the maximal accesible band gap (possibly in a setup not involving an extra power supply to sustain the perpendicular electric field), further accompanied by some magnification of the electric-conductance step on the band boundary, strongly suggests to focus future studies on graphene-based van der Waals heterostructures. 
Acknowledgments 
We thank Romain Danneau and Bartłomiej Wiendlocha for discussions. The work was supported by the National Science Centre of Poland (NCN) via Grant No. 2014/14/E/ST3/00256. Computations were partly performed using the PL-Grid infrastructure. 
Author contributions 
DS and AR developed the code for a rectangular BLG sample, DS performed the computations for a rectangular BLG sample; 
G.R. developed the code for the system with abrupt interface and performed the computations; all authors were involved in analyzing the numerical data and writing the manuscript. 
Additional information 
Supplementary information accompanies this paper. 
7 J. Phys.: Condens. Matter 31 (2019) 415501 
ORCID iDs 
Adam Rycerz 
https://orcid.org/0000-0001-7873-0183 
References 
[1] Qiao Z, Tse W K, Jiang H, Yao Y and Niu Q 2011 Twodimensional topological insulator state and topological 
phase transition in bilayer graphene Phys. Rev. Lett. 

107 <BR>256801 
[2] Maher P et al 2013 Evidence for a spin phase transition at 
charge neutrality in bilayer graphene Nat. Phys. 9 <BR>154 <BR><BR>
[3] Dean C T et al 2013 Hofstadter’s butterfly and the fractal quantum Hall effect in Moiré superlattices Nature 497 <BR>598 
[4] Ponomarenko L A et al 2013 Cloning of Dirac fermions in 
graphene superlattices Nature 497 <BR>594 <BR><BR>
[5] Cao Y et al 2018 Unconventional superconductivity in magicangle graphene superlattices Nature 556 <BR>43 
[6] Cao Y et al 2018 Correlated insulator behaviour at half-filling in magic-angle graphene superlattices Nature 556 <BR>80 
[7] Yang L, Deslippe J, Park C H, Cohen M L and Louie S G 2009 Excitonic effects on the optical response of graphene and bilayer graphene Phys. Rev. Lett 103 <BR>186802 
[8] Bonaccorso F, Sun Z, Hasan T and Ferrari A C 2010 Graphene photonics and optoelectronics Nat. Photon. 4 <BR>611 
[9] Grigorenko A N, Polini M and Novoselov K S 2012 Graphene plasmonics Nat. Photon. 6 <BR>749 
[10] Wang C R et al 2011 Enhanced thermoelectric power in dualgated bilayer graphene Phys. Rev. Lett 107 <BR>186602 
[11] Chien Y Y, Yuan H, Wang C R and Lee W L 2016 Thermoelectric power in bilayer graphene device with ionic liquid gating Sci. Rep. 6 <BR>20402 
[12] Mahapatra P S, Sarkar K, Krishnamurthy H R, Mukerjee S and Ghosh A 2017 Seebeck coefficient of a single van der Waals junction in twisted bilayer graphene Nano Lett. 17 <BR>6822 
[13] Suszalski D, Rut G and Rycerz A 2018 Lifshitz transition and thermoelectric properties of bilayer graphene Phys. Rev. B 
97 <BR>125403 
[14] Lee M J et al 2016 Thermoelectric materials by using twodimensional materials with negative correlation between electrical and thermal conductivity Nat. Commun. 7 <BR>12011 
[15] Sevinçli H 2017 Quartic dispersion, strong singularity, magnetic instability, and unique thermoelectric properties in two-dimensional hexagonal lattices of group-VA elements Nano Lett. 17 <BR>2589 
[16] Qin D, Yan P, Ding G, Ge X, Song H and Gao G 2018 Monolayer PdSe2: a promising two-dimensional thermoelectric material Sci. Rep. 8 <BR>2764 
[17] Mahan G D and Sofo J O 1995 The best thermoelectric Proc. Natl Acad. Sci. USA 93 <BR>7436 
[18] Kim H S, Liu W, Chen G, Chu C W and Ren Z 2015 Relationship between thermoelectric figure of merit and energy conversion efficiency Proc. Natl Acad. Sci. USA 
112 <BR>8205 
[19] Ioffe A F 1957 Semiconductor Thermoelements and Thermoelectric Cooling (London: Infosearch) 
[20] Cutler M and Mott N F 1969 Observation of Anderson Localization in an electron gas Phys. Rev. 181 <BR>1336 
[21] Rut G and Rycerz A 2014 Minimal conductivity and signatures of quantum criticality in ballistic graphene bilayer Europhys. Lett. 107 <BR>47005 
[22] Hao L and Lee T K 2010 Thermopower of gapped bilayer graphene Phys. Rev. B 81 <BR>165445 
[23] Goldsmid H J and Sharp J W 1999 Estimation of the thermal band gap of a semiconductor from Seebeck measurements 
J. Electron. Mater. 28 <BR>869 
[24] McCann E and Koshino M 2013 The electronic properties of bilayer graphene Rep. Prog. Phys. 76 <BR>056503 
[25] Kuzmenko A B, Crassee I, van der Marel D, Blake P and Novoselov K S 2009 Determination of the gate-tunable band gap and tight-binding parameters in bilayer graphene using infrared spectroscopy Phys. Rev. B 80 <BR>165406 
[26] Landauer R 1957 Spatial variation of currents and fields due to localized scatterers in metallic conduction IBM J. Res. Dev. 
1 <BR>233 
[27] Büttiker M, Imry Y, Landauer R and Pinhas S 1985 Generalized many-channel conductance formula with application to small rings Phys. Rev. B 31 <BR>6207 
[28] Esfarjani K and Zebarjadi M 2006 Thermoelectric properties of a nanocontact made of two-capped single-wall carbon nanotubes calculated within the tight-binding approximation Phys. Rev. B 73 <BR>085406 
[29] Alofi A and Srivastava G P 2013 Thermal conductivity of graphene and graphite Phys. Rev. B 87 <BR>115421 
[30] Zhang Y, Tang T-T, Girit C, Hao Z, Martin M C, Zettl A, Crommie M F, Shen Y R and Wang F 2009 Direct observation of a widely tunable bandgap in bilayer graphene Nature 429 <BR>820 
[31] McCann E 2006 Asymmetry gap in the electronic band structure of bilayer graphene Phys. Rev. B 74 <BR>161403 
[32] Yankowitz M et al 2014 Band structure mapping of bilayer graphene via quasiparticle scattering APL Mater. 2 <BR>092503 
[33] Grushina A L, Ki D-K, Koshino M, Nicolet A L A, Faugeras C, McCann E, Potemski M and Morpurgo A F 2015 Insulating state in tetralayers reveals an even–odd interaction effect in multilayer graphene Nat. Commun. 
6 <BR>6419 
[34] Nam Y, Ki D-K, Koshino M, McCann E and Morpurgo A F 2016 Interaction-induced insulating state in thick multilayer graphene 2D Mater. 3 <BR>045014 
[35] Kraft R et al 2018 Tailoring supercurrent confinement in graphene bilayer weak links Nat. Commun. 9 <BR>1722 
[36] Hirsch J E 1985 Two-dimensional Hubbard model: numerical simulation study Phys. Rev. B 31 <BR>4403 
[37] Gorbachev R V et al 2012 Strong Coulomb drag and broken symmetry in double-layer graphene Nat. Phys. 8 <BR>896 
8 
Supplementary Information 

Thermoelectric properties of gapped bilayer graphene 
Dominik Suszalski, Grzegorz Rut, and Adam Rycerz 
Marian Smoluchowski Institute of Physics, Jagiellonian University, Łojasiewicza 11, PL–30348 Kraków, Poland 
(Dated: June 15, 2019) 
I. DISPERSION RELATION, THE ACTUAL BANDGAP, AND DOPING 
The Hamiltonian H given by Eq. (4) in the main text leads to the dispersion relation for electrons [1] 
 
 
  1/2 
11 1
(e)2 22
E(k)=tU2 + !2v vk2 ±
±⊥ + F + 3
Γ(k)
,
2 4 2 
 
 2
  
1 
22 
⊥ −!2 v3k2
2
+ !2 vF k2
t
22 
⊥ +U2 +!2 v3k2
2
+2ξt⊥!3 v3vF k3 cos 3ϕ,
with Γ(k)= (S1)
t
4
where ± refers to the upper/lower band, k ≡ (kx,ky) is the wavevector, with k =0 referring to the K or K’ point (marked by ξ =1 or ξ = −1, respectively), k ≡|k|, and the angle 0 : ϕ< 2π is the argument argz of a complex 
(h)(e)
number z = kx + iky. For holes, we simply have E(k)= −E(k). This is a consequence of the combined particle
±± 
hole-reflection symmetry, which is preserved as we have neglected the next-nearest neighbor intralayer hopping (t2). The precise value of t2 (∼t⊥) is difficult to determine for BLG [2]; however, its effects are insignificant when discussing the band structure for |E|« t⊥. 
The value of a band gap Δ can be determined numerically via 
Δ 
(e)(h)
= minE(k)= −maxE(k). (S2)
−−
2 
For v3 =0 the band gap can be approximated as follows 
 
Δ ≈ Δ01 − Δ0 
 
√ 
2v32 
|U|t⊥
5v
3
− Δ2 
0 
with Δ0 (S3)
=
 
,4U2v2 FU2 + t2 
⊥ 
t⊥vF 
(Δ0 is the exact value of a gap in the absence of trigonal warping). Although its accuracy is better then 0.5% for t' =0.3 eV and |U| : 300 meV, the approximating Eq. (S3) is insufficient to present the numerical data in the energy scale used in Fig. 1(e) in the main text; instead, we have determined Δ for a given U directly by performing the minimization in Eq. (S2) numerically. 
In general, Δ < |U| for |U| > 0, leading to a peculiar (the Mexican hat-like) profile of the dispersion relation for 
(e)(h)
E(k) and E(k) bands in the presence of a gap (see the second paper in Ref. [1]). 
−− 
The dispersion relation given by Eq. (S1) also allows us to define the density of states at the Fermi energy (EF ) in a compact form 
ρ(EF ) ≡ ρ(e)(EF )+ ρ(h)(EF )= 1 
π2
JJ
JJ 
 
c=e,h 
∂ 
∂E 
A(c) 
m 
(E) 
JJJJ 
E=EF 
, 
(S4) 
m=± (e,h)
where we took into account spin and valley degeneracies (gs = gv =2) and A(EF ) denotes the area bounded 
± 
by the Fermi surface in the (kx,ky) plane. (In case E = EF is beyond the energy range of a given band, we put 
(c)
Am (E)=0.) In particular, for v3 =0 and max(0, |EF |− Δ/2) « t⊥, we have 
t⊥
ρ(EF ) ≈ Θ(|EF |− Δ/2) , (S5)
π(!vF )2 
with Θ(E) being the Heaviside step function. The prefactor in Eq. (S5) is often written in a form 2meff /(π2!), with the effective mass meff ≈ 0.033 me (where me is the free-electron mass) coinciding with the standard cyclotronic mass mC (EF )=(π!2/2)ρ(EF ). Integrating Eq. (S5) over the energy, we obtain the physical carrier concentration (or doping) given by Eq. (5) in the main text. 
At finite temperatures (T> 0), the effective carrier concentration can be defined via its relation to the Hall resistance, namely RH = −1/(eneff ) (with the electron charge −e), leading to 
  ∞   0 
neff = ne − nh = dEρ(E)fFD(µ, T, E) − dEρ(E) [1 − fFD(µ, T, E)] , (S6) 0 −∞ 
where ne (nh) is the concentration of electrons (holes) and fFD(µ, T, E)=1/ [ exp ((E−µ)/kBT )+1] is the Fermi-Dirac distribution function for a given electrochemical potential µ. 
II. TRANSMISSION PROBABILITY 
A. A semi-inifinite system in BLG 
In the case of an abrupt potential barrier separating weakly-and heavily doped regions (see Fig. 1(a) in the main text) we start from the substitution in the Hamiltonian (see Eq. (4) in the main text) p → p + eA, with a vector potential A =(Ax,Ay) (not further specified, since it is disregarded in the forthcoming calculations). The i-th component of the current (with i = x or y) corresponding to the Hamiltonian H reads, for the K valey, 
⎛⎞
0100 
ji (ψ)= ψ† 
· 
∂H 
∂Ai 
JJJJ JJJJ 
∂H 
⎜⎝ 
100 ν 
0001 
⎟⎠
· ψ 
with 
, 
(S7)
= evF
∂Ai
Ai→0 Ai→0 
0 ν 10 
where we have defined a dimensionless parameter ν = v3/vF . 
Next, the transmission probability can be found by simple mode-matching of wavefunctions in the heavily-and weakly-doped regions. Although deriving the corresponding functions is not a challenging task, they cannot be presented in a closed form. Below, we briefly present a procedure for finding the solutions corresponding to waves moving in desired directions (i.e., to the left or to the right): 
(i) The wavefunctions satisfying the Dirac equation Hψ (x, y)= Eψ (x, y) can be found algebraically in 
to the momentum representation, πi = !ki. The general solution has a form of a spinor ψ (x, y)= exp [i (kxx + kyy)] (ψ1,ψ2,ψ3,ψ4)T . There are four independent solutions, each corresponding to a different value of kx (for a system with periodic boundary conditions the transverse wave number ky =2πn/W , with n =0, ±1, ±2, ..., and W being the width of the system). 
(ii) In order to determine the corresponding direction of propagation for a given solution of the Dirac equation, one needs to follow the subsequent two-step procedure. First, one has to check a sign of the current employing Eq. (S7). The positive (or negative) sign indicates the transmitted/incoming (or reflected) wave. Second, when dealing with the wavefunction carrying zero current, it is important to check the sign of the imaginary part of kx. Positive (or negative) Im(kx) corresponds to a solution decaying exponentially to the right (or to the left). When discussing the transmission through a single barrier separating the heavily-and weakly doped regions we can, however, limit ourselves to the solutions with real kx =0 (which are normalizable with respect to the carried current). Out of four solutions of the Dirac equation, there are at most two carrying a positive or negative current. 
(iii) Subsequent calculation of the transmission probability becomes similar to the one usually performed for systems containing a sample area between two highly-conducting contacts. In the contact region, here modelled as heavilydoped bilayer graphene (for x< 0) the wavefunction takes a form 
1(2) 1(2) 1(2)
ψ1 ψ2
ψI (x)= ψ(x)+ rL,I (x)+ rL,I (x) , (S8)
R,I 12 
while the wavefunction in the weakly-doped sample area (x> 0) reads 
1(2) 1(2)
ψ1 ψ2
ψII (x)= t1 R,II (x)+ t2 R,II (x) . (S9) 
The lower indexes, R and L, refer to the waves moving to the right (jx > 0) and to the left (jx < 0), respectively. The 
1(2) 1(2) 1(2) 1(2)
upper indexes, 1(2), correspond to the two subbands. The parameters r, r, t, tare closely related to the 
1 212 
reflection and transmission probabilities and can be calculated via mode-matching at x =0, namely: ψI (0) = ψII (0). The final transmission probability is given by a sum over the all possible current ratios 
m 2 
ψn 
= |t|/jx ψm . (S10)
Tky n jx R,II R,I m,n 
The comparison of the Landauer-Büttiker conductance, following from Eq. (S10), with the number of open channels, which is determined solely by the dispersion relation given by Eq. (S1), is presented in Fig. 1. 

FIG.1: Theconductance(top)andthenumber ofopen channels(bottom)foran abruptpotentialbarrierinBLG,corresponding to the width W = 103 l⊥ =1.77 µm and the trigonal-warping strength t� =0.3 eV, displayed as functions of the chemical potential. The electrostatic bias between the layers (U ) is specified for each line. Right panels are zoom-ins, for low values of µ−Δ/2, for the data shown in left panels. 
B. Mode-matching for a rectangular sample 
In order to determine the transmission probability at a given energy: T (E), for a rectangular sample attached to the two heavily-doped regions (the leads), we employ the computational scheme similar to the presented in Ref. [3]. However, at finite-precision arithmetics, the mode-matching equations become ill-defined for sufficiently large L and µ, since they contain both exponentially-growing and exponentially-decaying coefficients. This difficulty can be overcome by dividing the sample area (0 <x<L) into Ndiv consecutive, equally-long parts, and then matching the wave functions for all (Ndiv+1) interfaces. Typically, using the double-precision arithmetic, we put Ndiv = 20 for L = W/20 = 104 l⊥. The necessary number of different transverse momenta ky =2πq/W (where q =0, ±1, ±2,..., ±qmax) 
varies with the energy E, scaling roughly as 2qmax +1 ∝|E|− Δ/2 (for |E| > Δ/2). For instance, to determine T (E) with a 10-digit accuracy for Δ= U =0 we took: 1955 2qmax +1 8149, with the lower (upper) value corresponding to E =0 and t ' =0.1 eV (E =2 meV and t ' =0.35 eV). 
III. SIMPLIFIED MODELS FOR TRANSMISSION-ENERGY DEPENDENCE 
A. Basic definitions 
For sufficiently large Δ, transmission spectra T (E) for either the single-barrier or the rectangular-sample case show essentially abrupt switching near E ≈ Δ/2, with some secondary details becoming irrelevant when calculating thermoelectric properties at nonzero temperature. Therefore, one can consider a family of simplified models for transmission-energy dependence, as given by Eq. (13) in the main text 
1
δ(E − 1 Δ) + δ(E + Δ) for α =0 
22
T (α)(E)= C(α)(Δ) × α−1 , (S11)
Θ(|E|− 1 Δ) |E|− 1 Δ for α> 0
22 
with δ(x) being the Dirac delta function, and Θ(x) being the Heaviside step function. A compact form of Eq. (S11) 
implies that the prefactor C(α)(Δ) is dimensionless for α =1 only; in general, we have C(α)(Δ) = eV−α+1 (for α � 0). 
The cummulants Ln, determining the thermoelectric properties via Eqs. (8)–(10) in the main text, are well-defined for the transmision-energy dependence of the form T (α)(E) (S11) with arbitrary α � 0. Substituting the derivative of the Fermi-Dirac distribution function 
∂fFD 11 
− = · , (S12)
∂E 4 cosh2 [(E − µ)/2kBT ] kBT 
and introducing the dimensionless variables 
Eµ Δ 
t = ,v = ,u = , (S13)
kBTkBT 2kBT 
we can write 
gsgv (kBT )n   t − vgsgvn+α−1
L(α) 
= dt (t−v)n T (α)(t · kBT ) cosh−2 ≡ (kBT )C(α)(Δ) L(α), (S14)
nn
h 42 h 
where the last factor (Ln (α)) is dimensionless and C–independent. Explicitely, we obtain 1 u − v 1 −u − v 
L(0) 
=(u − v)n cosh−2 +(−u − v)n cosh−2 , (S15)
n 
424 2 
  ∞   −y
1 t − v 1 t − v 
L(α) 
= dt (t − v)n(t − u)α−1 cosh−2 + dt (t − v)n(−t − u)α−1 cosh−2 (α> 0).
n 
424 2
y −∞ 
(S16) 
B. Maximal absolute thermopower 
Subsequent approximations, performed when calculating Ln (α) for y � 1, depend on mutual relation between x and 
y. In particular, for the Seebeck coefficient, we first rewrite Eq. (9) from the main text as follows 
kB (α) � (α)
S = LL. (S17)
10 
e Since the maximal |S| is expected for |v|« u, we can employ the approximation 
−t+v
1 t − vefor u<t< ∞,
cosh−2 ≈ (S18)
t−v
42 efor −∞ <t< −u, 
(valid for u ± v � 1), in order to get closed-form approximating expressions for the first three cumulants 
(α) v−u −v−u
L≈ e + e, (S19)
0 
(α) v−u(−v + u + α) − e −v−u(v + u + α),
L≈ e (S20)
1 
L(α) ≈ e v−u [(v − u)(v − u − 2α)+ α(α + 1)] + e −v−u [(v + u)(v + u +2α)+ α(α + 1)] . (S21)
2 
The right-hand side of Eq. (S17) can now be approximated by 
kB
S ≈ [(u + α) tanh v − v ] , (S22)
e 
with the maximal absolute value 
√√
−1
|S|max (kB/e)≈ (u + α)(u + α − 1) − ln u + α + u + α − 1 , (S23) 
appearing for 
√√
|S|
v≈± ln( u + α − 1+ u + α). (S24)
max 
These are Eqs. (14) and (15) in the main text.