High thermoelectric performance in excitonic bilayer graphene

We consider the excitonic effects on the thermal properties in the AB-stacked bilayer graphene. The calculations are based on the bilayer generalization of the usual Hubbard model at the half-filling. The full interaction bandwidth is used without any low-energy assumption. We obtain the unusually high values for the electronic figure of merit even at the room-temperatures which is very promising for the thermoelectric applications of the AB-bilayer structure. We discuss the effects of the interlayer Coulomb interaction and temperature on different thermal parameters in the bilayer graphene and we emphasize the role of the charge neutrality point in the thermal properties and within the excitonic insulator transition scenario. The calculated values of the rate of thermoelectric conversion efficiency suggest the possibility of high-performance device applications of AB-bilayer graphene.


INTRODUCTION
The discovery of the very high heat conductivity in graphene has led to a great scientific interest and surges many experimental and theoretical studies [1][2][3][4][5][6][7][8][9] on the thermal and thermoelectric properties of this unusual material. The anomalous high values of thermopower have been observed also in graphene monolayer in contact with a semiconductor substrate [10,11]. Recently, the thermal properties in bilayer graphene (BLG) systems have been also under intense studies [12][13][14][15][16][17][18][19][20] owing to the plausible device applications of the bilayer graphene as well as fundamental interest in this material due to the large band gap opening comparable to those in semiconducting systems [21]. The fundamental effect which observed experimentally was the electric field tunability of the thermoelectric power (TEP) [14] under the applied perpendicular electric field which opened new possibilities for the thermoelectric applications using bilayer graphene-based devices. Recently, it has been predicted, theoretically, [19] that the electric field effect could enhance the room-temperature thermopower in bilayer graphene more than four times that of the monolayer graphene and unbiased bilayer graphene. This effect has been confirmed by the experimental results at the low-temperature and with the dual gated bilayer graphene [22]. It has been shown that one can realize even higher thermoelectric power in a clean dualgated bilayer graphene device at the large electrical displacement field. The frequency dependence of the thermal conductivity in the AA and AB bilayer graphene and also the bias voltage dependence have been intensively studied in a series of theoretical works [23][24][25][26]. The thermoelectric efficiency of a high-performance thermoelectric material, for both power generation and cooling, is determined by its own enhancing (dimensionless) figure of merit (FOM) ZT , discussed well in Ref. [27], where a * e-mail:v.apinyan@intibs.pl δ-shaped transport distribution is found to maximize the thermoelectric properties. It has been shown in Ref. [28] that the thermoelectric power in bilayer graphene can be drastically increased by introducing the nanopores and by shifting the Fermi level. In this way, the authors obtained for the electronic part of the thermoelectric figure of merit the room temperature value ZT = 2.45 which renders the bilayer graphene as to be competitive with other power generation systems. It has been shown recently [29] that the thermoelectric figure of merit is considerably enhanced in bilayer graphene at the lowtemperature limit (ZT > 3) when a rigid substrate suppresses out-of-plane phonons by reducing their contribution to the thermal conductivity. By using an atomic tight-binding Hamiltonian, the authors in Refs. [30][31][32] have shown that one can achieve very high values of ZT in graphene devices composed of two partially overlapped graphene sheets coupled to each other in the cross-plane direction. The authors have shown that the figure of merit can reach the values higher than one in the junctions consisting especially of gapped graphene materials. Particularly, they have revealed that the phonon conductance in these layered structures gets strongly reduced compared to that in single layer graphene. The excitons can also be involved in thermal conductivity. The excitonic thermal conductivity has been reported in a few works [33][34][35]. For the first time, the excitonic thermal conductivity was calculated in Ref. [33]. The thermal conductivity of the excitonic insulator has been calculated in Ref. [36], in the semimetallic limit. Experimentally, the excitonic thermal conductivity was measured only in 1988 and reported in Ref. [37], where the provided data evidently shows the possibility of excitonic thermal transport in solids. Moreover, it appears that the bilayer exciton system displays a remarkable similarity to a thermoelectric couple as it was suggested in Ref. [38], where the authors suggested the experimental realization of a thermocouple device based on the bilayerexcitons, which has the potential to increase considerably the figure of merit ZT. Recently, the thermal transport in the exciton-condensate Josephson junction, composed of two graphene bilayers, has been studied [39] and the quasiparticle and interference heat currents through an insulating barrier have been calculated. It was shown in Ref. [39] that the phase difference of the condensates has an important impact on the thermal rectification effect in the junction.
Notwithstanding the numerous studies on the thermal properties in bilayer graphene heterostructures, there is no information about the influence and the role of the electron-hole bonding and electron-electron interaction on the thermal properties in the single bilayer graphene. In the present paper, we decided to accomplish this lack in the literature and we considered the excitonic effects in the thermal properties of the bilayer graphene structure. Basing on our previous results, concerning the excitonic gap parameter and chemical potential [38] in the suspended AB stacked bilayer graphene with excitons, we considered here the interacting bilayer Hubbard model with the full bandwidth of the interaction strength and without any low-energy assumptions. By using the equations of motion for the electric and heat current densities in the system and by employing the Kubo linear response theory, we calculated the in-plane components of the heat transport coefficients and we found their temperature dependence. The calculations show that the excitonic effects strongly enhance the thermoelectric figure of merit (2.5 < ZT e < 3 at the room temperatures) for a given range of frequencies of the external photonic field and for the selected values of the local interlayer interaction parameter. Moreover, the high values of the thermoelectric conversion efficiency parameter indicate about the important role that play the excitons in AB-BLG which opens the possibility for the direct device application of the bilayer graphene as a promising thermoelectric material. Particularly, the thermoelectric energy conversion efficiency is strongly enhanced (> 20%) at the high photon's frequencies and at the large Carnot efficiency. A detailed discussion about all obtained thermoelectric effects is given and the conditions for the high thermoelectric performance in AB-BLG system are highlighted. The presented work could have the important impact in the technological applications of the bilayer graphene as an efficient thermoelectric generator, which operates at the room temperatures.
The paper is organized as follows: in the Section II, we introduce the considered model and we discuss the singleparticle and excitonic Green's functions. In the Section III, we calculate the electronic and thermal current densities within the Mahan's formalism and in Section IV we present the thermal kinetic coefficients within the Kubo theory. Furthermore, in the Section V, we give the numerical calculation results and in Section VI we give a conclusion to our paper. In the A, the electronic and heat current density operators are calculated in details and B is devoted to the calculation of heat-heat response function.

II. THE BILAYER HUBBARD MODEL
Here, we consider the AB Bernal stacked bilayer graphene (BLG) and we describe it with the help of the interacting bilayer Hubbard Hamiltonian. The AB-BLG system is composed of two coupled honeycomb layers with sublattices A 1 , B 1 and A 2 , B 2 , in the bottom and top layers respectively. In the z-direction the layers are arranged according to Bernal Stacking (BS) order [39], i.e. the atoms on sites A 2 in the top layer lie just above the atoms on sites B 1 in the bottom layer graphene, and each layer is composed of two inter-penetrating triangular lattices. The schematic representation of the AB-stacked bilayer graphene is shown in Fig. 1. The Hamiltonian of the interacting BLG system at the half-filling is where the parameters γ 0 and γ 1 are the intralayer and interlayer hopping amplitudes. The summation mm ′ , in the first term, in Eq.(1), denotes the sum over the nearest neighbors lattice sites m and m ′ in different honeycomb layers in the bilayer graphene. We keep the small letters a 1 , b 1 and a 2 , b 2 for the electron operators on the lattice sites A 1 , B 1 and A 2 , B 2 respectively. The summation index ℓ = 1, 2 denotes the layers in the BLG. Particularly, we use ℓ = 1 for the bottom layer, and ℓ = 2 for the top layer. Furthermore, we have η = a 1 , b 1 for ℓ = 1, and η = a 2 , b 2 , for ℓ = 2. The spin variable σ takes two possible directions (σ =↑, ↓). Next, n ℓ,ηm in Eq.(1) is the total electron density operator for the layer ℓ and n ℓ,mησ = η ℓ,mσ η ℓ,mσ is the electron density operator for the η-type fermion in the layer ℓ and with the attached spin σ. We consider the suspended BLG structure with pure electronic layers and without initial doping in the system. The condition of half-filling in each layer reads as n ℓ = 1, for ℓ = 1, 2, where n ℓ is the total electron density operator for the layer ℓ. Furthermore, U , in the Hubbard term in Eq.(1), parametrizes the intralayer Coulomb interaction, while the parameter W denotes the local interlayer Coulomb repulsion between the electrons located on sites B 1 and A 2 , in different layers. In Fig. 1, here, we have shown different physical parameters entering in the Hamiltonian, in Eq. (1). The excitonic formation between the layers is also shown in Fig. 1. We will put γ 0 = 1, as the unit of energy in the system, and we set k B = 1, = 1 and e = 1 throughout the paper. For a simple treatment at equilibrium, we initially suppose the balanced BLG structure, i.e., with the equal chemical potentials in the both layers µ 1 = µ 2 . The Hamiltonian, given in Eq.(1) could be reduced to an effective Hamiltonian which is linear in the electron density terms. This procedure is discussed in details in Ref. [38] and we give in A the linearized form of it after the Hubbard-Stratanovich decouplings of the interaction terms.
A. The partition function and the Green's functions in the BLG We write the total action of the BLG system in the following form Here, the first two terms are the Berry terms for the layers with the indices ℓ = 1, 2 where we have introduced the imaginary time τ at each lattice site position m. The upper limit of integration over τ , in all terms in Eq. (3), is given as β = 1/T where T is the temperature. The Hamiltonian H (τ ) of the BLG system, in the last term in Eq.(3), is described in Eq.(1). Furthermore, after some linearization procedure [38,40] of the nonlinear interaction terms in Eq.(1) and transforming the fermionic operators into the Fourier space representation we can rewrite the fermionic action in the following form where we have introduced the four component fermionic Nambu-spinors Ψ k,σ (ν n ) = [a 1k,σ , b 1k,σ , a 2k,σ , b 2k,σ ] T at each discrete state k in the reciprocal space and for each spin direction σ =↑, ↓. The fermionic Matsubara frequencies ν n are ν n = π(2n + 1)/β with n = 0, ±1, ±2.... Next, G −1 k,σ (ν n ) stands for a 4 × 4 inverse Green's function matrix. It is defined as follows Here, the diagonal elements in the matrix G −1 k,σ are defined as E 1 (ν n ) = −iν n − µ eff 1 and E 2 (ν n ) = −iν n − µ eff 2 , where the effective chemical potentials µ eff 1 and µ eff 2 , have been introduced in the form µ eff The parameter ∆ in Eq.(6) represents the excitonic order parameter which is defined as ∆ = W b † 1,σ (r, τ )a 2,σ (r, τ ) . We considered here the case of the homogeneous BLG structure where the pairing is between the particles with the same spin orientations, i.e. ∆ σσ ′ = ∆ σ δ σσ ′ ≡ ∆. The parameter γ ℓk , in Eq. (6), is the renormalized energy dispersion in the layer ℓ andγ ℓk = γ ℓk t, where γ ℓk = δ e −ikδ ℓ . The vectors δ ℓ are the nearest neighbor vectors in different layers ℓ = 1, 2. The components of δ ℓ , for the layer with ℓ = 1, are δ = −a 0 / √ 3, 0 , and a 0 = √ 3a is the sublattice constant, while a is the carbon-carbon length in the graphene sheets). In the layer with ℓ = 2 (the top layer in the stacking), we have δ By the convention, we put a ≡ 1, for both layers. It is not difficult to realize that γ 2k = γ * 1k ≡ γ * k and it follows that γ 2k =γ * 1k ≡γ * k , where we have omitted the layer index ℓ. The partition function of the system takes the following form where the action in the exponent is given in Eq. (5). Furthermore, by assuming the half-filling condition in each layer of BLG, we obtain a set of self-consistent equations for the excitonic order parameter ∆ and the chemical potential µ in the BLG where the dimensionless coefficients α ik , in Eq.(8) with i = 1, ..4, are given as Here, The coefficients β ik in Eq.(9), with i = 1, ..4 are given by the relations The function n F (x), in Eqs. (8) and (9), is the Fermi-Dirac distribution function n F (x) = 1/ e β(x−µ) + 1 and the energy parameters ε ik define the excitonic band structure in the interacting BLG. They are given as We have introduced the effective Fermi levelμ in Eq. (13) in terms of the effective chemical potentials µ 1eff and µ 2eff , i.e.,μ = (µ 1eff + µ 1eff ) /2. The solution of the system of equations given in Eqs. (8) and (9) and the following discussion about the behavior of the excitonic order parameter and chemical potential are well discussed in Refs. [38,41]. In the present work, we will use the data for ∆ and µ, obtained in Ref. [38], in order to calculate the thermal properties in the suspended bilayer graphene. Just for the completeness, and the logical continuation, we will follow the results given in the work in Ref. [41] and we rewrite here the explicit expressions of the fermionic Green's functions in the interacting bilayer graphene system, under consideration. For the sublattices A and B in the bottom layer in BLG we have for the normal single-particle Green functions [42] where the k-dependent coefficients γ ik are defined as where We see that a ′ 1k , a ′ 2k and a ′ 3k could be obtained from the coefficients a 1k , a 2k and a 3k , given in Eq. (11), by interchanging the effective chemical potentials µ and by setting simultaneously ∆ + γ 1 = 0. The anomalous or the excitonic Green function in bilayer graphene is defined as F a2b1k (ν n ) = − a 2k (ν n )b † 1k (ν n ) and we obtain where the coefficients β ik are given in Eq. (12) above. For the layer ℓ = 2, it is not difficult to show that G a2k (ν n ) = G b1k (ν n ) and G b2k (ν n ) = G a1k (ν n ). Furthermore, the Green functions in Eqs. (14) and (17) will be used for the calculation of the thermal kinetic coefficients in bilayer graphene.

III. THE ELECTRONIC AND HEAT CURRENT DENSITY OPERATORS
The thermal conductivity is principal for the study of thermal properties in solid state materials. In semiconducting materials with the band gap in the excitation spectrum, the electrons can have a predominant contribution to the heat conduction. It has been shown recently [43] that the thermal transport properties in bilayer graphene could be manipulated by the sp 3 interlayer bonding and the role of the interlayer phonon coupling has a small impact on the in-plane thermal conductivity in the graphene structures. Here, we consider the role of the excitons on the thermal transport properties in the suspended AB-BLG. We suppose that the temperature gradient along x-axis causes an electric field along that axis. In Fig. 2, we have presented the electron and hole currents through the stable excitonic condensate state in the AB-BLG structure. The directions of the temperature gradient ∇T , thermoelectric field E, total heat current j H and electric (j e ) and hole (j h ) currents are shown in the picture. The AB-BLG structure is maintained at the constant temperature gradient, between the hot (in red) and cold (in blue) sides of the structure. The charge and heat currents are related within the linear response equations since the temperature gradient ∇T leads to the electric potential gradient ∇V in the system. We write here the linear response equations for the electric current density (j e ) and total heat current density (j H ). Then the total heat current density is related to the energy current density j E via the relation where µ is the chemical potential in the system. With the definition of thermal current density in Eq.(18), we have a set of coupled linear equations for the electric and heat current density. After defining the appropriate forces for the generation of currents we have in general Here, ∇µ is the gradient of the chemical potential in the system, which is the consequence of the electron concentration gradient. Next, F e (ω) is the force acting under effective electric field E(ω) produced by the temperature gradient ∇T . Thus, F e (ω) = −eE(ω). The frequency dependent transport coefficients L ij (ω) with i, j = 1, 2 are the imaginary parts of retarded response functions [44] which will be calculated later on. With these definitions, the Onsager relation between non-diagonal transport coefficients is L 12 (ω) = L 21 (ω). Furthermore, we assume that there are no concentration gradients in the considered problem and we put ∇µ = 0. We can obtain the analytical expressions of important thermal parameters in the system just by reverting the linear response equations in Eqs. (19) and (20) by writing them in the following equivalent form where σ e (ω) is the electrical conductivity, caused by the temperature gradient, S(T ) is the Seebeck coefficient, Π is the Peltier coefficient and κ is the thermal conductivity which is defined as the coefficient of proportionality (when the electric current density is zero) between the total heat current density and the temperature gradient, i.e., j H (ω) = −κ∇T . We have In the definitions, above, we have putted e ≡ 1. The thermoelectric conversion efficiency is evaluated with the help of the thermoelectric figure of merit ZT which is the regrouped combination of the principal thermoelectric parameters S(T ), σ e and κ ZT = S 2 σT /κ.
The thermal conductivity κ is given in general as κ = κ el + κ ph with the electronic (κ el ) and phononic (κ ph ) contributions. The thermoelectric efficiency represents the ratio between the thermal power input and the electrical power output [45], for a given material, and the maximum value of it is defined as [46] where T h and T c are the temperatures in the hot and cold sides of the material and the temperature difference is ∆T = T h − T c . The average T av is the mean operating temperature in the sample and we have T av = (T h +T c )/2. Thus the thermoelectric conversion efficiency is the product of the Carnot efficiency η c = (∆T /T h ) and a reduction factor which is the function of the material's figure of merit parameter, defined as Z = S 2 σ/κ. It is important to notice here that 50% of the waste heat and the solar thermal energy lie in the range 300-500 K, while no efficient thermoelectric materials have been found for low-temperature applications (< 500 K). For example, a thermoelectric generator, composed of a material with the FOM ZT = 1, is expected to reach only 5% energy conversion efficiency when ∆T = 100 K. Most research efforts are focused on materials development with the increasing FOM ZT and the maximum ∆T which is highlighted by the Carnot efficiency η c . Recently, it has been shown that the thermoelectric power in monolayer and bilayer graphene follows the semiclassical Mott formula [15,47] at the low-temperature limit, while at the high carrier density the temperature dependence was shown to be linear, which demonstrates the absence of phonon drag modes in the bilayer graphene structure. Regarding this, we will neglect the phononic part of thermal conductivity κ ph and further derivations concern only to the electrons. As it has been mentioned earlier, we consider the longitudinal thermal transport in bilayer graphene by supposing the temperature gradient along x-direction in the planes of the BLG. In order to evaluate the thermal parameters given in Eq. (23) we should calculate the linear response coefficients L ij (ω) i, j = 1, 2 which are diagonal in i, j for the isotropic systems. They are given with the help of the Kubo formula, and, in the Matsubara notations, we have for L ij (ω m ), (ω m here are the bosonic Matsubara frequencies ω m = 2πm/β) along x-direction The dynamical transport coefficients in Eqs. (23), (24) and Eq.(26) are obtained after the analytical continuation in Eq.(26) ıω m → ω + iδ, where ω are the real frequencies and δ is an infinitesimal constant [44]. We have Here, in this paper, we will omit the notation of the components along x for the principal thermoelectric parameters in Eq. (29) and the response coefficients L ij (ω). After Eqs. (19) and (20), the current density operators in the time-ordered correlation function in Eq. (26) are defined by the convention that j 1x = j ex and j 2x = j Hx , where j ex and j Hx are the electronic and heat current density operators. To calculate the current-current correlation functions in Eq. (26), for a given response function, we should know the explicit expressions of the electronic and heat current density operators in our bilayer graphene system. The calculation of those operators could be done within the Mahan's formalism evaluated in Ref. [48], for the electric and heat current density operators. Thus, for the electric current density operator, we introduce the electron polarization operator P e , which is the sum over the position R m (of a lattice site m) and the electron number operator n m at that site, i.e., P e = m R m n m , while for the energy current density operator j E we construct an operator R E which is the sum over position R m and effective Hamiltonian density h effm (see in A), i.e., R E = m R m h effm . Furthermore, by using the continuity equations and the Heisenberg equation of motion for the current density operators it is easy to show that the electron and heat current density operators satisfy the following equations Here, R mx means the x-component of the position vector R m at the lattice site m. Furthermore, after some calculations (see in A) we can write the electronic and heat current density operators in the following forms where the fermionic spinors Ψ † kσ (ν n ) and Ψ kσ (ν n ) have been introduced in Eq.(5), in the Section II A, while the electron velocity matrix operator υ kx is defined as (see in A, for details) Here v kx and v * kx are the x-components of the kdependent electron velocity operator v k and its complex conjugate. We have defined here The total heat current density operator along the x direction could be calculated as j Hx = j Ex − µj ex . After calculating the commutator [H, R E ] (see in A 2, for details), we get where the total velocity matrix operator υ H kx is defined as Each element of the matrix υ H kx in Eq.(34) represents a 2 × 2 matrix, and we have The matrix υ † tk in Eq. (34) is the Hermitian conjugate of the matrix υ tk . The elements of the matrices υ gk and υ tk are introduced in the following form The function Λ k in the expression of v 1k is where the summation vector u takes the values 0, ±a 1 , ±a 2 with the vectors a 1 and a 2 being the unit cell vectors in graphene lattice structure IV.

THERMAL LINEAR RESPONSE FUNCTIONS
Here, we give the results of calculations of the response functions using the definitions in Eqs. (26) and (27). For this, we put the expressions for the electronic and heat current densities in Eqs. (30) and (33) into the expression of the correlation function in Eq.(26) and we apply the Wick's theorem [42] for the statistical average of the product of four fermionic operators. Then, after using the definitions of the normal and excitonic Green's functions in Eqs. (14) and (17), we obtain Furthermore, we perform the summation over fermionic Matsubara frequencies ν n in Eq.(41) and the analytical continuation, given in Eq. (27), in order to obtain the expression of the retarded real-frequency response function. After taking the imaginary part of the retarded function, we get For the non-diagonal element L 12 (ω m ) we have Again, after performing the Matsubara summation in Eq.(43) we get for the real-frequency dependent coefficient L 12 (ω) The coefficient L 21 (ω) is the same as L 12 (ω). The second diagonal coefficient L 22 (ω) could be also calculated using Eq.(33) for the total heat current density operator. The analytical expression of the coefficient L 22 (ω m ), in terms of the normal and excitonic Green functions is given in B. Here, we present only the real frequency value of it, after analytical continuation into upper half complex plane iω m → ω + iδ The Thus, for the k-dependent coefficient |v 1kx | 2 , in the first term in Eq. (45), we obtain Next, the x component of the k-dependent velocity op- and we have Furthermore, when calculating numerically the response coefficients L ij (ω), we will replace δ-Dirac functions, figuring in Eqs. (42), (44) and (45), by a δ-Lorentzian with the appropriate broadening parameter.

V. THE RESULTS AND DISCUSSION
A. The thermal conductivity In this section, we discuss the numerical results concerning the thermal conductivity in the bilayer graphene. Particularly, we will show the interaction induced excitonic effects in the behavior of this function. In Fig. 3, we have shown the low-frequency behavior of the thermal conductivity function κ(ω). Different values of the interaction parameter W have been considered and the zero temperature limit is considered. We see in Fig. 3 that there are double peaks, appearing at the low-frequency region. When increasing the interlayer Coulomb interaction parameter W , the positions of those peaks are displacing to the right (the "blue"-shift effect in the spectrum) on the ω-axis and also the thermal conductivity is decreasing when augmenting W . The low frequency peaks disappear when W is crossing the critical value at the charge neutrality point (CNP) W c = 1.4899γ 0 (for details see in Refs. [38] and [41]). The similar single peaks, in the low-frequency region, have been observed in Ref. [23] concerning the biased AA-stacked bilayer graphene, studied within the tight-binding model Hamiltonian. Particularly, it has been shown that the intensities of peaks decrease when the temperature increased. The doubling of peaks in the case of our ABstacked bilayer graphene could be related to the excitonic effects considered here. The opposite behavior was shown for κ(ω) when varying the applied bias voltage (see in Ref. [23]). It is important to note that at the value W = 1.48γ 0 , close to the CNP value W CNP = 1.4899γ 0 , the thermal conductivity is very small (see the dark-red curve in Fig. 3). At the CNP point, the exact solution for the chemical potential in the AB-BLG, discussed in Ref. [38], jumps into its upper bound solution branch µ LB → µ UB (such a behavior of the chemical potential has been observed also experimentally in Ref. [49]), and there are two solutions for µ and the excitonic order parameter ∆ at W = W CNP (this becomes clear when solving the system of coupled equations in Eqs. (8) and (9)). One of those solutions is the lower bound solution (LB): µ LB = −1.871γ 0 and ∆ LB = 0.0364γ 0 , and the another one is the upper bound solution with µ UB = −0.499γ 0 and ∆ UB = 0.0249γ 0 . Furthermore, at the CNP, κ(ω) becomes very small, especially at the upper bound solution for the chemical potential µ UB , shown in the picture, in Fig. 4. We see in Fig. 4 that the double-peak structure in the thermal conductivity function practically disappears at the value W = W CNP = 1.4899γ 0 . For the values of W higher than W CNP the double-peak structure is completely absent. The spectrum of the electronic thermal conductivity function at the high photon's frequencies is shown in Fig. 5, where we considered the full bandwidth of the interlayer Coulomb interaction parameter, passing through the low and strong interaction limits. We observe in Fig. 5 the additional and very broad peaks at the high frequency region. In this case, there is a large"blue"-shift broadening of thermal conductivity when increasing W at W > W CNP while at the values W < W CNP a small "red"-shift effect is also visible (see the curves corresponding to W = 0, W = 0.5γ 0 , W = γ 0 and W = 1.48γ 0 ). At the CNP point, we have shown only the curve corresponding to µ UB , because the other plot for µ LB is very similar. The existence of thermal conductivity peaks in the high frequency domain could be explained as the result of the high-frequency excitonic effects recently observed in bilayer graphene when consid- ering the optical conductivity in this system [40,50,51]. The small values of κ(ω) imply the large values for the thermoelectric figure of merit, thus, the relatively small values of the thermal conductivity in the high frequency region, as compared with those in the low frequency regime, could be interesting for the thermoelectrical device applications of the BLG.

B. The figure of merit and electrical conductivity
The thermoelectric FOM, given in Eq.(24) can be rewritten as a function of the response functions as The ω dependence of FOM is important to reveal out the optimal frequencies of the incident photons at which the FOM attains its maximal values. Theoretically, the thermoelectric figure of merit in bilayer graphene system has been discussed recently in Refs. [26,29]. Particularly, the ω-dependence of FOM has been considered in Ref. [26] concerning the AA-stacked bilayer graphene and its temperature dependence has been found for the biased structure. Here, we consider the effects of interactions and the excitonic pairing on the behavior of FOM. Our treatment goes far beyond the standard tightbinding approaches discussed in Refs. [26,29] and permits to consider properly the realistic thermoelectric application mechanisms for the bilayer graphene. In Fig. 6, we have presented the results for FOM for the same values of the interaction parameter W as given in Fig. 3. It is clear that at the certain values of photon's frequencies the values of FOM are very high. The highest FOM value ZT e = 2.592 is attained for W = γ 0 and at ω = 2.05γ 0 .
The ω-dependence of FOM, for W ≥ W CNP , is given in Fig. 7. It is clear from the results in Figs. 6 and 7 that the highest FOM values are attained in the region W < W CNP . The positions of maxima of curves in Figs. 6 and 7 don't change when increasing the temperature up to room temperature region 0 ≤ T ≤ 300 K (not presented here). The frequency dependence of the Seebeck coefficient is presented in Fig. 8. We have shown the curves which correspond to (from bottom to top) W = 0, 0.5γ 0 , γ 0 , 1.48γ 0 , 1.4899γ 0 , 2γ 0 , 3γ 0 , 4γ 0 . There are two curves in Fig. 8 corresponding to the critical value W CNP = 1.4899γ 0 . One of them (with the lower bound solution for µ: µ = µ LB ) is in the negative region S(ω) < 0, while the other one (with the upper bound solution for µ: µ = µ UB ) is in the positive region S(ω) > 0. They define, indeed, the outermost boarders of the coefficient S(ω) at which the it changes its sign. We see also in Fig. 8 that the Seebeck coefficient does not vary considerably with the photon's frequencies ω. We see clearly in that the Seebeck coefficient changes its sign when crossing the CNP point (see the curves in dark yellow and orange, in Fig. 8, corresponding the lower and upper bounds of the chemical potential µ). It has been shown experimentally by Yuri M. Zuev, et al., [52], that the sign of the TEP changes across the charge neutrality point as the majority carrier density switches from electron to hole. This result is consistent with our behavior of the Seebeck coefficient, given in Fig. 8. The similar effects in thermopower have been shown also in Refs. [17,18]. The temperature dependence of the FOM and the electrical conductivity function σ e (T ) are shown in Figs. 9, 10 and 11. Different interaction limits are considered in the pictures. In Fig. 10, the electrical conductivity function σ e is normalized to the minimum conductivity quanta σ 0 = e 2 / . We see in Fig. 9 that at W < W CNP the FOM values are very close to the value ZT e = 1, for a very large interval of temperature, while for W > W CNP FOM is decreasing rapidly when increasing the temperature (see the small-dashed curve in green and the solid curve in orange). The highest values of the electrical conductivity function are obtained at W = γ 0 , which corresponds to the behavior of FOM presented in Figs. 6 and 9. It is worth to mention here that the external field frequency ω, for each curve in Fig. 9, is fixed to the value at which the FOM attains its maximal values in Figs. 6 and 7 (for the case T = 0). The behavior of the thermoelectric figure of merit as a function of temperature, starting from zero up to room temperature value T = 300 eV (corresponding to the normalized temperature T /γ 0 = 0.0086), is given in Fig. 11 below. We see that the FOM is very stable to the variation of temperature and the considerable deviation from this behavior takes place at the value W = 2γ 0 = 6 eV (see the small-dashed curve in green). The perfect thermoelectric FOM is achieved at W = γ 0 = 3 eV, at which very high values of FOM (given in the interval 2.4 ≤ ZT e ≤ 2.592) are obtained in the interacting BLG with excitons in the temperature interval T ∈ [0, 300]K. This result suggests that the interacting bilayer graphene could be a high-FOM performance thermoelectric device, operating at the room temperatures. The electrical conductivity function σ e is continuously increasing with T in this range of temperature (see in Fig. 10) and contributes in the high values of ZT e (see in Eq. (24)) and we will see that the principal shape of thermoelectric conductivity κ is due to the electrical conductivity function σ e . The difficult experimental task, related to the observation of our results is related to the finding of appropriate values of external field frequencies at which the BLG will be calibrated to the suitable interaction value W . This could be achieved with the help of the adiabatic variation of the external field values (see in Ref. [53], where such a situation is described in the context of the excitonic Josephson tunnel junction based on AB-BLG) until the suitable effect takes place in the BLG. It is also interesting to discuss the explicit W -dependence of thermal parameters in the interacting AB-BLG in order to observe the changes of those parameters at the CNP point and finding the similarities in their behavior. In Figs. 12, 13, 14, and 15, we shown the behavior of principal thermoelectric quantities defined in Eq. (29). In Fig. 12, the FOM is given for W in the interval W ∈ [0, 4γ 0 ]. This is a large enough interval that covers the small and strong interaction limits in the AB-BLG. Then we have considered several values of frequencies ω, from low (ω = 0.1γ 0 = 0.3 eV) up to very high frequencies (ω ∼ 3γ 0 = 9 eV). We see in Fig. 12 that for all fixed values of the photon's frequencies ω the thermoelectric FOM parameter has a drastic decrease at the value W = W CNP , corresponding to the charge neutrality point. It is clear in Fig. 12 that the largest FOM values are given in the interval of Coulomb interaction W ∈ (0, W CNP ), i.e., for W below the CNP value and for which the AB-BLG system is in the mixed state composed of the excitonic insulator and condensate states (see in Ref. [38,41] Fig. 12. The behavior of thermal conductivity κ as a function of interaction parameter W is given in Fig. 13. Contrary, the parameter κ takes the small values in the interval where ZT e is high (this is  conductivity vanishes at the CNP point for the photon's frequencies ω ∈ [0.2γ 0 , γ 0 ], while it is finite at the CNP for the frequencies in the interval ω ∈ [1.8γ 0 , 3γ 0 ]. In Fig. 14, we have shown the critical behavior of the Seebeck coefficient S as a function of W . It is clear in Fig. 14 that the coefficient S reflects well the particle-hole symmetry in our BLG system and changes the sign exactly at the CNP point given by W = W CNP . Moreover, the Seebeck coefficient does not vary when changing the frequencies and only small changes happen in the electronregion S < 0, near the CNP point. This is due to the fluctuations of the chemical potential in the deep vicinity of the CNP point. The behavior of the electrical conductivity σ e as a function of W is given in Fig. 15 and it is very similar to that of the electronic thermal conductivity function κ, given in Fig. 13. As it was mentioned earlier, it repeats the shapes of the curves of κ in Fig. 13 but with different scales. This suggests that in the limit T → 0 the Wiedemann-Franz law should be satisfied. It is worth to mention that the electrical conductivity function presented here (as a consequence of temperature gradient in the system) differs much from the optical conductivity function discussed in Ref. [40], where a threshold behavior is observed in the conductivity spectrum. We see in Fig. 15 that the electrical conductivity in the AB-BLG system is not a threshold process and the function σ e is non zero for all values of W , except the charge neutrality point, i.e., W = W CNP . Another interesting behavior of the electrical conductivity function σ e (ω) is related to the case of the non-interacting AB-BLG system. This case is presented here, in Fig. 16. The conductivity limit corresponding to the non-interacting monolayer graphene σ e = σ MG = e 2 /4 could be achieved for three different values of the external photon's frequencies: low-frequency limit ω 01 = 0.617γ 0 , medium frequency limit ω 02 = 2.618γ 0 and high-frequency limit with ω 03 = 4.254γ 0 . The electrical conductivity corresponding to the non-interacting bilayer graphene with the minimal electrical conductivity σ e = σ Bi = e 2 /2 FIG. 16. (Color online) The electrical conductivity, normalized to the elementary conductivity quanta σ0 = e 2 / , as a function of the normalized photon's frequencies ω/γ0. The zero interaction limit is considered in the picture. Indicated red points in the picture correspond to different photon's energies at which the conductivity function gets the values equal to the minimal conductivities in the non-interacting monolayer (σe = e 2 /4 ) and bilayer (σe = e 2 /2 ) graphene systems. The temperature is set at zero. [54][55][56] could be achieved also at three different values of the photon's energies. Namely, at ω 01 = 0.452γ 0 (low-frequency), ω 02 = 2.742γ 0 (medium-frequency) and ω 03 = 3.404γ 0 (high-frequency). Those points are indicated on the curve of the normalized electrical conductivity function in Fig. 16. The inset, in Fig. 16 shows the interaction dependence of the electrical conductivity function. The frequency ω/γ 0 = 0.3266 (or ω = 0.97 eV) is chosen such that σ e = σ 0 in the limit of the noninteracting bilayer W/γ 0 = 0.

C. Thermoelectric efficiency
In general, it is known that higher average FOM ZT av and larger temperature difference ∆T will produce the higher energy conversion efficiency η eff . For the traditional heat engines which rival the Carnot efficiency η c = ∆T /T h with ∆T = 400 K and ZT av ≈ 3 we have the thermoelectric conversion efficiency η eff ∼ 25% [57,58]. Here, we calculate the thermoelectric conversion efficiency coefficient η eff for two different regimes of the Carnot efficiency. Namely, η c = 0.4 which corresponds to temperature difference ∆T = 200 K with T c = 300 K and T h = 500 K and η c = 0.6 which corresponds to temperature difference ∆T = 300 K with T c = 200 K and T h = 500 K. In Fig. 17, we have shown the ω-dependence of the coefficient η eff (in %) for the first case. Different values of W are considered in the picture. We see that for W ≤ W CNP the coefficient η eff is very high in a given interval of the photon's frequencies, situated in the deep ultraviolet and vacuum ultraviolet range of the spectrum. Particularly, when ω ∈ (1.84γ 0 , 2.6γ 0 ) the values of the coefficient η eff are in the range η eff ∈ (7, 14.14)%. Contrary, at W = 2γ 0 (see the large dot-dashed curve blue, in Fig. 17) the coefficient η eff is reduced considerably. Thus, at W = γ 0 and for ω ∈ (5.52, 7.8) eV the thermal performance of the interacting AB-BLG is very high (the solid black curve in Fig. 17). In Fig. 18), we have presented the W -dependence of the coefficient η eff with η c = 0.4 and for different values of the photon's frequencies ω. The general behavior of η eff in Fig. 18 is governed by the thermoelectric FOM, presented in Fig. 12, in the Section V B. Namely, there is a drastic decrease of η eff at the charge neutrality point W = W CNP and the high values of η eff correspond to the high energy photons (see the solid green line, in Fig. 18) with ω = 2.4γ 0 = 7.2 eV). The thermoelectric coefficient enhances considerably when increasing the Carnot efficiency to the value η c = 0.6 (see in Figs. 19 and 20 here). We observe in Fig. 19 that for W < W CNP (see the curves at W = 0 and W = γ 0 ) the coefficient η eff is higher than 5% practically for all values of the external photon's energies. In the narrow regions on ω-axis (ω ∈ (0, 0.495γ 0 ) = (0, 1.485) eV) and ω ∈ (2.3γ 0 , 2.83γ 0 ) = (6.9, 8.49) eV, we have η eff 10% for the case W = 0 (see the dashed red curve in Fig. 19). For W = γ 0 (see the solid black curve in Fig. 19) we have η eff 10% for the photon's energies in the intervals ω ∈ (0.161γ 0 , 0.545γ 0 ) = (0.483, 1.635)eV (low-energy photons) and ω ∈ (1.787γ 0 , 2.87γ 0 ) = (5.361, 8.61) eV (high-energy photons). There is also a very narrow and high energy interval ω ∈ (1.95γ 0 , 2.23γ 0 ) = (5.85, 6.69) eV, where we have η eff 20%. For the W -dependence of η eff , in this case, we have the same behavior as for the previous case with η c = 0.4, but the values of η eff are much higher. The obtained here values for the energy conversion efficiency coefficient permit to qualify definitively the bilayer graphene among the high performance thermoelectric materials and open a new prospect towards the usage of bilayer graphene as a material with outstanding thermoelectric properties.

VI. CONCLUSION
We considered here the thermoelectric properties in the excitonic BLG system with the Coulomb interaction effects, included in the band structure dispersion relation. By supposing the temperature gradient in the individual graphene-layers in the AB-stacked BLG, we considered the electrical and heat current components along the x-axis in the two-dimensional layers. Then after using the Mahan's method [44] for the electronic and heat currents density operators, we evaluated the explicit analytical forms of the linear response coefficients within the Kubo-Green formalism. The principal thermoelectric parameters have been calculated in the excitonic BLG and the role of the charge neutrality point has been discussed. We obtained the additional single peaks (due to the interlayer excitons) in the spectrum of the thermal conductivity function κ, for the values of the interaction parameter W in the region W < W CNP and in the low- frequency regime ω ∈ [0, γ 0 ]. While for W ≥ W CNP we get the usual behavior for the conductivity function in the low-frequency regime (see for example in Ref. [23]). The observed double-peak structure in the low-frequency κ spectrum and the additional peaks in the high photon's frequency region are attributed to the existence of the excitonic insulator state in the BLG [38] and the existence of the excitonic condensate state in the AB-BLG, even for the case W = 0 [41]. We have calculated the Seebeck coefficient and electrical conductivity function σ normalized to the minimal conductivity quanta e 2 / . We have shown that the thermoelectric figure of merit (FOM) could hit very high values (ZT 2.5) in the excitonic BLG, for a given interaction parameter regime and a given interval of the photon's energy. Additionally, we have calculated the temperature dependence of the FOM for different Coulomb interaction regimes and we found nearly a constant picture up to room temperatures suggesting the high stability of FOM with respect to temperature changes and we have shown the possibility of room temperature applications of our BLG as an efficient thermoelectric engine. Another remarkable effect is related to the electrical conductivity mediated by the temperature gradient in the BLG. Indeed, we have shown that the limits of the minimal electrical conductivities e 2 /4 and e 2 /2 typical for the non-interacting monolayer and bilayer graphene systems could be reached at three different operating frequencies for our BLG system (at W = 0 or ∆ = 0). Moreover, we found that thermopower changes the sign at the charge neutrality point W = W CNP by defining the branches for electrons (S < 0, at W < W CNP ) and holes (S > 0, at W > W CNP ), thus, behaving in accordance with the definition of the excitonic insulator state, given in Ref. [38]. The very high values of the thermoelectric conversion efficiency coefficient obtained here (η eff 25%, for not very large temperature differences ∆T ), qualifies the excitonic BLG as a high-performance thermoelectric energy conversion material and open a new prospect for the practical technological applications of the excitonic AB-BLG such as in the modern solar cell energy converters, thermoelectric nano-converters for the effective heating or cooling the microprocessors, the heat converter engines for the ecological locomotions, and etc..
At the end, we have to distinguish, among many cited extraordinary thermal properties in the AB-BLG, two principal advantages for which the thermoelectric applications of such a system are strongly recommended. First one, is related to the high values of ZT (ZT ∼ 2.5) at the room temperatures which suggest the possibility of direct device applications of the excitonic BLG system as a high-performance thermoelectric material operating at the room temperatures. The high values of the thermal conversion efficiency (η eff 25%), from another hand, imply that the AB bilayer graphene system is a new efficient thermoelectric energy conversion device with the unique excitonic band structure. In our suggestion, the excitonic bilayer graphene could be reliable also as a model system for the upcoming thermometry devices at the nanoscale and the new generation laser sensors.
The authors have no competing interests to declare Appendix A: Mahan's formalism

Electric current density operator je
The equation of continuity for the electric current density operator in the Matsubara time formalism reads as where n(τ ) is the total electron density operator in the bilayer graphene n(τ ) = n 1 (τ ) + n 2 (τ ). The electric polarization operator, introduced in the Section III is related to the electron density operator n by P Thus, the calculation of the electric current density operator j e reduces to the calculation of the standard commutator in the right-hand side in Eq.(A2). The effective Hamiltonian density operator h effn after the linearization procedure described in Ref. [38] is Here, the parameter ∆ appears after the linearization of the interlayer interaction term in Eq.(1), while the effective chemical potentials µ 1eff and µ 2eff , defined in the Section II A, result from the linearization of the intralayer Hubbard interaction term U and the interlayer interaction term W . The total electron density operator n m , in the commutator in Eq.(A2), is n m = ℓ,σ n ℓ,mσ , where the summation index ℓ = 1, 2 denotes the layers in the bilayer graphene. Then, for the layer ℓ = 1 we get from Eq.(A2) and for the layer ℓ = 2 we have Combining the results in Eqs.(A4) and (A5) we obtain for the total electron current density operator j e = j ℓ=1 e +j ℓ=2 e the following result It is not difficult to show that the commutator of the total electron density operator with the interlayer part of the Hamiltonian density, given in Eq.(A3) is zero. Indeed, we have The vanishing of this expression becomes apparent when passing into the Fourier space representation for the fermionic operators b † 1mσ and a 2mσ . We have and we obtain where ∆k is the wave vector difference ∆k = k ′ − k.
The sum over the lattice vectors m in the expression in Here, a 1 and a 2 are the unit cell vectors defined in the Section III and m 1 and m 2 are integer numbers. We have used, in Eq.(A10), the fact that m e i∆kRm = N δ (∆k) and mi m i = 0, because we consider the infinite lattice structure. Furthermore, we will use the expression of the electric current density operator, given in Eq.(A6) when evaluating the linear response coefficients, presented in the Section IV. The Fourier transformed form of the expression in Eq.(A6) could be done easier by using the rules in Eq. (8). After a simple calculation we write where the k-dependent velocity operator v k is defined in Eq. (32). Furthermore, we can use the fermionic Nambu spinors notations, introduced in Section II A, and we write the electric current density operator in the more symmetric form. Namely, we have where the electron velocity matrix υ k is At the end, we remark that the commutation of the electronic density operator n m with the other terms containing the single electron density operators is obvious and we don't give here those derivations.

Total heat current density operator jH
We define an operator which is the sum over the position R m and the effective Hamiltonian density operator h effm given in Eq.(A3), i.e., R E = m R m h effm . Next we use the equation of continuity ∂H/∂τ − i∇j E = 0 which leads to the following equation for the energy current density operator j E = i∂R E /∂τ , i.e., j E = i [R E , H], where we have used the Heisenberg equation of motion for the operator R E : ∂R E /∂τ = − [R E , H]. Then we have Therefore, the calculation of the energy current density operator is reduced to the calculation of the commutator between the effective Hamiltonian densities at different site positions h effn and h effm . When calculating this commutator with the effective Hamiltonian density given in Eq.(A3), we should take into account the commutator between (a) intralayer-intralayer hopping terms, (b) intralayer and interlayer hopping terms, (c) intralayer and effective chemical potential terms and finally (d) interlayer and chemical potential terms. The last commutator is zero as it was discussed above. We present here the calculations of the other commutators, cited here.

a. Contribution of the intralayer hopping terms
The commutator between the intralayer terms with the Hamiltonian densities h effn (γ 0 ) and h effm (γ 0 ) has the Furthermore, we can use the transformations in Eq.(A8) and we have It is easy to realize that = iγ 0 (γ 1 + ∆) mδ1 σ R m+δ1 a † 1,mσ a 2,m+δ1σ − h.c.
Here, we have used the fact that δ 2 = −δ 1 . Let's remark also that δ1 δ 1 e ikδ1 = a ′ =0,a1,a2 Thus the summation over δ 1 is equivalent to the summation over the lattice vectors a ′ = 0, a 1 , a 2 . Therefore, we Here, we present the calculation of the commutator between the intralayer hopping term and the last two terms in the effective Hamiltonian density given in Eq.(A3). The commutator which contains the effective chemical potential µ 1eff is ue −iku a † ℓ,kσ a ℓ,kσ +i γ 0 (γ 1 + ∆) N kσ a ′ a ′ a † 1,kσ a 2,kσ e ika ′ + b † 2,kσ b 1,kσ e ika ′ −h.c.) Next, the total heat current density operator j H is defined in Eq. (22) in the Section III. Taking into account the expression of the electric current operator, given in Eq.(A11) we can write for j H j H = j E − µj e = 1 N kσ ℓ a † ℓ,kσ a ℓ,kσ v 1k + 1 N kσ a † 1,kσ a 2,kσ + b † 2,kσ b 1,kσ v 2k +b † 2,kσ a 2,kσ v 3k + a † 1,kσ b 1,kσ v 4k + h.c. , (A26) where the velocity operators v ik with i = 1, ...4 are defined in Eq. (37), in the Section III. After using the Nambu spinor representation for the fermionic operators, the result in Eq.(A26) can be rewritten in the form given in Eq. (33), in the Section III. For calculating the correlator under the integral in the right hand side in Eq.(B1) we use the form of the total heat current density operator, given in Eq.(A26), in the Section A 2 c. When expanding the correlator we get 100 terms from which survive only 18, because many of them vanish when expanding the four fermion correlators with the help of the Wick theorem. As an example, we consider a term which gives the product of the normal Green's function with the single excitonic correlator The first two terms in Eq.(B2) vanish because of the symmetry of the action of the system. Then we have F a2b1 = F b1a2 , and we get Here, we give the final result for the function in Eq.(B1) Furthermore, we use the explicit forms of the normal and excitonic Green's functions given in Eqs. (14) and (17) and we perform the summation over the fermionic Matsubara frequencies ν n . Then we use the Eq.(27) and we get the expression of the response coefficient, given in Eq. (45), in the Section IV.