{'fraction_non_alphanumeric': 0.038840819442677185, 'fraction_numerical': 0.018036645883700218, 'mean_word_length': 4.693045997826874, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thermopower of multilayer graphene

We systematically calculate thermopower of biased and unbiased multilayer grphene systems. The effect of screening to a bias field perpendicular to the graphene planes is taken into account self-consistently under the Hartree approximation. The model including nearest neighbor hopping and the more complete Slonczewski-Weiss-McClure (SWMcC) model are both considered for a comparison. The effect of impurity scattering is studied for monolayer and unbiased bilayer graphene and is treated in terms of the self-consistent Born approximation. For a monolayer graphene, only when the effect of impurity scattering is taken into account, could all the qualitative aspects of the experimental results be correctly reproduced. A small energy gap opens for the biased trilayer graphene. When the carrier density is high enough, a gap is also found for unbiased trilayer and quad-layer graphene. The biased bilayer graphene shows the largest thermopower among all the systems studied.


I. INTRODUCTION
Recently, thermopower of graphene systems have attracted much attention. [1][2][3][4][5][6][7][8][9][10][11][12][13] In experiments on monolayer graphene (MLG), while the high carrier density region is well accounted for in terms of the Mott's relationship, the low carrier density region shows deviation away from this behavior. [2][3][4] This is explained in terms of either an electron-hole puddle model [1] or the coherence effect between the conduction band and the valence band mediated by impurity scattering [7]. Besides the MLG, a bilayer graphene (BLG) system is also very interesting. The opening of a gap by an external electric field applied perpendicular to the layer planes [14][15][16][17][18][19] introduces a new degree of tunability to the system and is recently shown by the authors to enhance greatly the thermopower. [11] Despite the above progresses, several questions remain to be answered. First, the dependence of thermopower on carrier density as a function of temperature for a MLG has two features which are not fully explained in previous works. The first one is the shift of thermopower peak with respect to the charge neutral point as temperature decreases. In addition, thermopower shows a monotonic dependence on temperature for carrier densities smaller than that of the thermopower peak. [1,7,9] The second question is the systematic dependence of thermopower on the number of layers. Further more we are interested in finding out whether or not an applied biased potential will enhance thermopower for the multilayer graphene as the same as the bilayer system. Below we shall address these questions.
Full microscopic calculations of thermopower of MLG in the presence of charged impurity scattering are performed in terms of the self-consistent Born approximation. Better agreement with experiments [2][3][4] as compared with previously reported results are achieved. [1,7,9] We also study thermopower of impure unbiased bilayer graphene. Results which are qualitatively similar to those in monolayer graphene are obtained, but for much higher impurity concentrations.
For a multilayer graphene with layer number from two to six, we consider both the nearest neighbor tight binding model and the more general Slonczewski-Weiss-McClure (SWMcC) model. [20,21] For a multilayer system the effect of screening is important in determining the band structure, as it produces non-uniform charge densities across the layers. [14,15,[22][23][24][25][26][27][28] In this work, by using the self-consistent Hartree approximation we obtain potential energies of different layers for unbiased systems with three or more layers and biased systems with a perpendicular electric field for two or more layers. For unbiased graphene multilayers, peak values of the thermopower are close to each other irrespective of the model used. For a biased multilayer graphene, only thermopower of the bilayer shows significant enhancement, thermopower of the trilayers also increases slightly, but thermopower of systems with more layers decreases under a bias. The huge enhancement of thermopower in a bilayer graphene is a direct result of gap opening under a perpendicular electric bias. [11] In a trilayer graphene, a much smaller gap opens and hence the thermopower enhancement is very small. While for layer number larger than three, no gap opens under a bias. Thus no thermopower enhancement is expected for a multilayer graphene with layer number larger than three.

II. MODEL AND METHOD
A. The SWMcC model We consider multilayer graphene systems stacked in the standard AB (Bernal) stacking between consecutive layers. [29][30][31][32] In order to correctly reproduce the band structure, a tight binding type of model up to the fifth nearest neighbor hopping is employed. [20-22, 29, 31, 32] For an N -layer graphene, the Hamiltonian is transformed in the xy-plane to wave vector space, whereas the direction perpendicular to the plane is kept in the real space.

B. Nearest neighbor tight binding model
Apart from the full SWMcC model, a simplified model up to nearest neighbor intralayer and interlayer hopping is widely used to study many problems. In this case, the Hamiltonian matrix is written as [29,30,[33][34][35] in which In the absence of any external field or gate voltage, it is known that the above Hamiltonian could be decomposed into subsystems of bilayer and monolayer graphene. [30,33,36,37] For an odd-layered multilayer graphene, there are one MLG subsystem and (N − 1)/2 different BLG subsystems. While for an even-layered multilayer, there are N/2 different bilayers. Label the different BLG subsystems with an index m, and take the basis vector as , the Hamiltonian matrix of the m-th BLG subsystem is written as where γ 1m =2γ 1 sin mπ 2(N +1) . m takes the value of N − 1, N − 3, · · · , 2 (for odd N ) or 1 (for even N ).

C. Screening effects
In our former work, a bias voltage is shown to greatly enhance the thermopower of BLG. [11] It is unclear whether this is true also for other multilayer graphene.
In the presence of a bias voltage or a charged gate, the potential energy difference would be induced between different layers. To determine the distribution of potential energies, self-consistent calculations are required to take into account of the Coulomb screening arising from charge redistribution among the different layers. [14,22,24,38] For unbiased multilayer graphene with layer number larger than two, the charge and potential energy distribution should also be obtained selfconsistently.
There are two schemes to account for the charge redistribution. One is to consider the screening of a charged gate with fixed carrier density by the multilayer graphene system [22]. The other approach is to consider the charge redistribution problem of the multilayer graphene system in the presence of a perpendicular electric field, while the carrier density is controlled by another gate voltage applied equally to all layers [38]. Here we follow the second approach.
Suppose an external electric field E 0 is applied to the free-standing multilayer graphene system along the direction perpendicular to the graphene planes from layer 1 to layer N (E 0 =0 for unbiased systems). After establishing equilibrium, the i-th layer has a total excess electron density of n i . According to Gauss's law, the screening electric field pointing from the j-th to the (j+1)- means screening, e is the absolute value of the electron charge, and ǫ r =3 is the static dielectric constant of the multilayer graphene [7]. The total electric field is thus 4 a 2 is the effective area per carbon site), which means the number of excess electrons per site residing on the i-th layer. Define two parameters V 0 =eE 0 d 0 and γ= e 2 d0 2ǫrΩ0 . d 0 ≃3.5Å is the vertical distance between adjacent layers. The potential energy difference between the j-th and the (j+1)-th layer is The total difference of on-site potential energy is thus Since only differences in the potential energies are relevant, we would take the potential energy of the first and N -th layer symmetric with respect to zero by shifting the energy reference. Thus we takeṼ N =(V N − V 1 )/2 and V 1 =-Ṽ N . Potential energy differencesṼ j+1 −Ṽ j between consecutive layers are not influenced by this shift of energy reference, and are still determined by Eq. (6). For a certain V 0 and a fixed average doping, we calculate the set of doping {x i } self-consistently. Label the 2N ×2N matrix which diagonalizes the Hamiltonian matrix as U (k), that is Here, we have added the on-site potential energies ψ k to the model and have kept the notation of the Hamiltonian matrix unchanged as H(k). The diagonal matrix H V (k) is k independent and is defined as diag[Ṽ 1 ,Ṽ 1 , ... ,Ṽ N ,Ṽ N ]. The α-th column of U (k) stores the α-th eigenvector of H(k) corresponding to an eigenenergy of ǫ α (k)=[H d (k)] αα . Suppose the number of wave vectors considered in the Brillouin zone is N k , the self-consistency condition for determining x i is thus (8) µ is the chemical potential. f (x)=1/(e βx + 1) is the Fermi distribution function, in which β represents the inverse temperature 1/k B T with k B denoting the Boltzmann constant. After a set of convergent results for {x i ,Ṽ i ; i = 1, ..., N } is obtained, the corresponding on site potential energiesṼ i (i = 1, ..., N ) are substituted into the model to calculate the thermopower of interest in this work.

D. Scattering by charged impurities
Transport properties of graphene systems sensitively depend on the nature and concentration of impurities. In particular, charged impurities are shown to dominate the various transport properties of MLG [39][40][41][42] and also play a very important role in the transport of BLG [43]. In the following, we would consider the effect of charged impurity scattering in the MLG and unbiased BLG systems. The treatment would be at the level of self-consistent Born approximation (SCBA). [7,33,[44][45][46][47][48] Due to the limitation in computation time and the fact that calculations performed for clean systems are enough to answer if a bias could enhance thermopower of a multilayer graphene, in this work we would concentrate on clean systems for layer number larger than two.
The formalism is presented in detail in our former work on thermopower of gapped BLG [11]. The thermopower is represented in terms of the linear response coefficients as [1,7,11,49] where T is the absolute temperature and e is the magnitude of electron charge. The linear response coefficients are obtained as In the Matsubara notation, the correlation function reads [49] The two subindices 'i' and 'j' both run over 1 and 2. Ω is the total effective area of the system defined as the area per layer multiplied by N , the number of layers. d=2 is the dimensionality. j 1 and j 2 denote particle current and heat current operators, respectively. [11] The two linear response coefficients in terms of which the thermopower is obtained are written as where the subindex j is either 1 or 2. The kernels are defined as [11] (13) with Γ 1 (k, z, z ′ ) as the vertex function corresponding to the wave vector k. Here, z=ǫ ± i0 + and z ′ =ǫ + i0 + . j k 1 is the matrix for the charge current at wave vector k and is obtained from the Hamiltonian matrix H(k) of the multilayer graphene as [11,50,51] In H(k) the potential energy distribution H V (k) with the on site energies determined self-consistently are incorporated, which affects the energy spectrum but does not affect the current operator.
For an N -layer multilayer graphene, the Green's function is represented in terms of a 2N ×2N matrix as The free Green's function is obtained as G 0 where I 2N is a 2N -dimensional unit matrix and µ is the chemical potential determined from H(k) for a certain carrier density and temperature. The self energies and the vertex functions are determined self-consistently in terms of the following iteration functions [7,11] n i is the impurity concentration averaged to per layer. v i (q) is the electron-impurity scattering potential. For charged impurity scattering, v i (q) is taken as of the Thomas-Fermi type [1,7,46] v ǫ r is the effective dielectric constant from lattice and substrate, ǫ r =3 is adopted in this work [7,41,46]. d i is the distance between the impurities and the graphene plane and would be set as zero in the present work [1,7]. q T F is the Thomas-Fermi wave number and is obtained from the long-wavelength-limit static polarizability of the corresponding noninteracting electron system [1,46] as with the static polarizability The particle number operator is defined as A factor of '2' comes from the two fold degeneracy in spin. The subindex 'c' means retaining only connected Feynman diagrams in evaluating the expectation value. [11] The self-consistent Born approximation (SCBA) enters when making averages over impurity configurations as [52] where N i = n i Ω is number of impurities in the system under consideration. For clean systems, G k (z) reduces to the free Green's function and the vertex function Γ 1 (k, z, z ′ ) reduces to the bare charge current matrix j k 1 . [11] III. RESULT AND DISCUSSION

A. Monolayer graphene
The experimental results [2][3][4] for thermopower of a MLG show two interesting features. First, for a certain temperature, thermopower follows 1/ |x| (x is the average number of excess electrons per site) for high carrier densities but then the magnitude decreases and changes sign close to the charge neutrality point. Theoretical calculations in terms of both the Boltzmann transport equation [1] and the microscopic Kubo's formula [7] have successfully reproduced the above features. The deviation of thermopower from the 1/ |x| behavior close to the charge neutrality point is ascribed to electron-puddle formation [1] or coherence between the conduction band and valence band [7], which both imply the coexistence of electron-like and hole-like characters of the carriers. The second feature is that, as temperature decreases the thermopower decreases monotonously for all carrier densities and the peak position also shifts toward lower carrier densities. [2,3] The above behaviors are very similar to the results presented in Figs. 1(c) and 1(d).
Results in Fig. 1 are obtained in terms of the nearest neighbor hopping model with hopping integral as 2.7 eV. As shown in Fig. 1(a) and 1(b), calculations on clean systems are unable to reproduce the second feature mentioned above. [9] Previous works on transport properties of a MLG had confirmed the importance of charged impurities. [39][40][41][42] Here, we treat the impurity scattering in terms of SCBA. [7,44] We would neglect the shift of chemical potential by the impurity potential. [7,11] The results for a series of temperatures with impurity concentration n i = 10 12 cm −2 , which corresponds to a moderately disordered sample [40], are shown in Fig. 1(c). For the electron-hole symmetric band as considered here, the relationship S(−x) = −S(x) generally holds. The full curve for 300 K in Fig. 1(c) is explicitly calculated, the above relationship is seen to hold very well. For all other curves, only the hole doping part with x ≤ 0 is calculated explicitly and the part with x > 0 is obtained in terms of S(x) = −S(−x). It is clear that, Fig. 1(c) agrees qualitatively very well with experiment by Zuev et al. [2] As shown in Fig. 1(d), increasing n i to 5 × 10 12 cm −2 , which corresponds to dirty graphene samples, the peak positions are shifted to higher carrier densities and the maximum values of thermopower are further suppressed for all temperatures, which are qualitatively similar to the experimental results by Wei et al. [3] In the degenerate region for relatively high doping concentrations, the semiclassical Mott's relationship generally holds. [1,53] In these cases, for a certain carrier density (or, gate voltage), the low temperature thermopower scales linearly with temperature. Experimentally, the Mott's relationship holds very well for graphene at high carrier densities but breaks down close to the charge neutrality point. Fig. 2 shows our results for three impurity concentrations. The linear temperature dependence holds even for doping close to (but still larger than) the peak position. For carrier densities smaller than at the peak position (e.g., x≃±1.4 × 10 −4 for T=300 K and n i =10 12 cm −2 ) where electron-like and hole-like carriers coexist, the linear Mott's relationship is not valid. As impurity concentration increases, the peak position shifts to higher carrier densities, region of carrier density for the violation of linear temperature dependence is correspondingly enlarged.
We would like to compare our results with previous calculations taking into account of scattering by charged impurities. Hwang et al. [1] study the problem in terms of semiclassical Boltzmann equation approach. There, the high and low carrier density regions are treated separately. Hence it is difficult to provide the dependence of peak position with temperature. Another work by Yan et al. [7] starts from the microscopic Kubo's formula but has only focused on the low temperature limit and is also unable to reproduce the temperature dependence of the peak position observed experimentally [2][3][4]. Compared to the above two works, we have shown that in terms of fully microscopic calculations incorporating the effect of scattering by charged impurities, both the two features mentioned in the beginning of this section could be reproduced successfully. Our results also confirm the picture that the peak position shifts to higher carrier density with the increase of both temperature and impurity concentration. [7,11]. Having achieved a reasonable success by using the SCBA to treat scatterings by charged impurities in a MLG, we shall proceed to study the bilayer system. We have previously studied the effect of charged impurity scattering on the thermopower of gapped (biased) BLG [11]. In a gapped BLG, localization effect becomes important in the low carrier density region, which is out of the reach of SCBA. [54][55][56] There we restricted our attention to the systems with dilute impurity concentrations. [11] However, in an unbiased BLG with a finite carrier density at the Fermi surface we expect to have a much stronger screening which allows us to consider much larger impurity concentrations.
As a comparison to the results for MLG, we present in Fig. 3 thermopower of clean BLG and impure BLG for two typical impurity concentrations: 10 12 cm −2 and 5×10 12 cm −2 , obtained in terms of the nearest neighbor hopping model with γ 0 =3 eV and γ 1 =0.3 eV. [11] For each impurity concentration, the carrier density dependence of thermopower for two typical temperatures, 300 K, 200 K are considered. Lower temperature requires much larger number of wave vectors to converge, which are difficult to approach. In contrast to the MLG, impurity concentration up to 10 12 cm −2 does not change the qualitative aspects of the results. Only when the impurity concentration is increased to be as large as 5×10 12 cm −2 , the result becomes similar to that observed in experiments. [12,13] In both experiments, the samples show semiconducting like transport behaviors with a low mobility [12,13], indicating appreciable impurity concentrations and are thus in agreement with our results.
The result of Fig. 3(c) is very similar to that of the Fig. 1(c) for a MLG except with an impurity concentration five times larger. This is due to a much better screening in a bilayer system. We present in Fig. 4 the Thomas-Fermi screening wave vectors of both systems at two different temperatures: 300 K and 50 K. At 300 K, q T F of BLG is more than five times that of MLG around the charge neutral point. This difference becomes even larger at 50 K. For the carrier density where the thermopower peaks (x ∼±2×10 −4 ) for BLG, q T F of BLG is still more than two times that of MLG. According to Eq. (16), Eq. (17), and Eq. (18), the effective strength of impurity scattering is approximately proportional to n i q −2 T F . Hence a BLG with a much larger concentration of charged impurity shows the same qualitative behavior as observed in MLG with a smaller impurity concentration.
The stronger screening in BLG as compared to that in MLG has been used to make the conjecture that short-range scatterers may also play an important role in the transport of gapless BLG in addition to the charged impurities. [57] It would be interesting to see how short-range scatterers would change the results presented above.

C. Unbiased multilayer graphene
Now, we study the evolution of thermopower as a function of layer number for an unbiased multilayer graphene. From now on we would only consider clean graphene multilayers. First, we consider the simple nearest-neighborhopping model, ignoring temporarily the screening effect. As shown in previous works, this simplified multilayer graphene model could be decomposed into independent subsystems of MLG and BLG. [30,33,36,37] The two linear response coefficients L 11 and L 12 are then obtained as summations of L (m) 11 and L (m) 12 , which are the corresponding values for the various subsystems labeled by m. [30] Thermopower of multilayers is then obtained from Eq. (9). The room temperature thermopowers obtained in this way are shown in Fig. 5(a) with layer number up to six. Peak value of the thermopower does not show very large variation with layer number. Different from MLG, thermopower for larger layer number samples usually show secondary peaks associated with the onset of higher conduction or valence bands contributing to transport. [11] Results presented in Fig. 5(a) does not take the effect of screening into account, so electrons on different layers feel the same electrostatic potential. However, as discussed in Sec. II C, multilayer graphene with layer number larger than two should have different carrier densities for different layers. When this effect is treated self-consistently, results for the multilayer graphene described in terms of the simplified model are presented in Fig. 5(b). Except for some quantitative difference, the qualitative behavior is the same as in Fig. 5(a). The positions and magnitudes of various peaks are separated into two groups depending on whether N is even or odd. For both cases, the peak value of thermopower increases and the peak position continuously shifts to higher carrier density with increasing layer number. In the presence of screening effect, site energies of the various layers become unequal, hence the subsystem decomposition is not valid and Fig. 5(b) is obtained starting from the full model.  [20-22, 29, 31, 32], with the screening effect incorporated. Since the particlehole symmetry is now explicitly broken, the curve becomes asymmetric with respect to the point of zero carrier density. Specifically, the positive and negative peak value becomes unequal and the position of zero thermopower is shifted away from zero doping. Another qualitative difference as compared to the results obtained from the simplified model is that, the peak value of the even-layer (N =2, 4, 6) and odd-layer (N =3, 5) systems show opposite behaviors compared to Fig. 5(b). This nontrivial result arising from further neighbor hoppings contained in the SWMcC model needs experimental confirmation.
The reason for the above grouping between even and odd layer systems could be understood qualitatively from the subsystem decomposition in Sec. IIB. Thermopower of various bilayer like graphene subsystems form a series which show monotonous dependency on the effective interlayer hopping γ 1m in Eq. (5). Thermopower of monolayer like subsystem does not fall into this series. Since monolayer like subsystems are contained only in odd-layer graphene, it makes sense that the even layer and odd layer multilayer graphene should group separately. Though self-consistent calculation and inclusion of further neighbor hopping both make the subsystem decomposition invalid, numerical results contained in Figs. 5(b) and 5(c) are still consistent with the above picture.

D. Biased multilayer graphene: The electronic structure
Transport properties, like conductivity and thermopower, sensitively depend on the underlying electronic structure. In order to understand the thermopower of biased multilayer graphene to be presented in the next section, we first discuss in this section the electronic structure of a biased multilayer graphene, taking trilayer and quad-layer graphene systems as two examples. In the presence of an electric field perpendicular to the layer plane, inversion symmetry of a multilayer graphene is broken explicitly. In this case, even for the bilayer system, charge redistributes to screen the electric field. As in the above section, we take this screening into account in terms of the Hartree approximation by treating separate layers as parallel plates with certain density of net charge which are to be determined selfconsistently. [14, 22-24, 38, 58] We focus on two related questions. One is the robustness of the Dirac-like linear dispersive band in the presence of an electric field for a multilayer graphene with an odd number of layers. The other is the possibility of opening a full energy gap by an external electric field, similar to the BLG. [11,14,15] We find that the effect of temperature on the electronic structure is very tiny. Thus, only zero temperature results are presented here. For the biased systems, without loss of generality, the bare interlayer potential energy dif- ference induced by the electric field is taken as V 0 =1 eV. For an unbiased charge neutral trilayer graphene, as shown in Figs. 6(a), the system is gapless. But when the carrier density is increased to be as high as 2×10 13 cm −2 , a small gap about 5 meV is induced below the chemical potential. Two Dirac cone like structures still exist in the system, but they do not touch at the cone vertex and are usually immersed in bilayer like bands. Applying a perpendicular electric field as high as V 0 =1 eV, a small gap about 17 meV is induced even for the neutral system x=0, as shown in Figs. 6(c). The magnitude of the gap enhances slightly with the increase of x as shown in Fig.  6(d).
For the biased quad-layer graphene, as shown in Figs. 7(c) and 7(d), a full gap is not observed. For the unbiased quad-layer graphene, though the neutral system is still gapless (Fig. 7(a)), a full gap below chemical potential is observed when the carrier density is high enough, as shown in Fig. 7(b). The latter feature, together with that in Fig. 6(b), arises from the screening effect. [31,32,36,58] Screening effect for the unbiased multilayer graphene systems for N >2 matters because, for layer number larger than two, N/2 (for even N ) or (N + 1)/2 (for odd N ) non-equivalent layer sets appear which usually have different on site energies and electron densities.
Finally, we present in Fig. 8 the density of states (DOS) of the trilayer and quad-layer graphene for several typical parameter sets. The corresponding room temperature chemical potentials are marked by the dotted vertical lines. As compared to the unbiased systems, external bias introduces some Van Hove singularities in the density of states.

E. Biased multilayer graphene: The thermopower
Previously, we have shown that the opening of a gap in biased BLG greatly enhances thermopower of the system [11,59]. Calculations in the previous section show that a small gap opens also in other multilayer graphene systems, such as the trilayer graphene. It is thus interesting to see whether or not the thermopower is likewise enhanced in these systems.
Results for multilayer graphene systems at room temperature with layer number up to 5 are presented in Fig.  9. Strength of the external electric field is taken as V 0 =1 eV. Fig. 9(a) contains results obtained from the nearest neighbor hopping model. Fig. 9(b) is for the SWMcC model.
Compared with zero bias results, the large enhancement of thermopower under a bias is specific to the BLG system. Though peak value of thermopower is slightly enhanced in trilayer graphene (Fig. 9(c)), it decreases under a bias for both the quad-layer and quintuple-layer systems.
In a biased BLG, peak value of room temperature thermopower is roughly proportional to the gap size. An external field as large as V 0 =1 eV opens a gap of about 288 meV, which yields a peak room temperature thermopower of about 412 µV/K. [11] From the previous subsection, V 0 =1 eV gives a gap of about 20 meV for trilayer graphene. Following this thread of argument, room temperature thermopower of the biased (V 0 =1 eV) trilayer graphene would be around 50 µV/K, which is in qualitative agreement with Fig. 9(b).
For layer number larger than three (e.g., Fig. 7(c) and 7(d) for the quad-layer system), usually no full gap opens in the biased multilayer graphene system since the bottom and top layers almost decouple for so thick samples. In addition, as could be seen from Figs. (6), (7) and (8), the chemical potential is usually not situated inside the gap, even though a gap opens for certain doping and external electric field. Both of these effects suppress the thermopower for biased multilayer graphene systems with N ≥3.
To confirm the self-consistency of the results, we analyze the applicability of the Mott's formula, which relates the thermopower to the longitudinal conductivity as [1,5,7,53] S = − π 2 k 2 B T 3e ∂lnσ(µ) ∂µ . The conductivity is obtained from the linear response coefficient as σ=e 2 L 11 /T . We present in Fig. 10 results obtained by a full microscopic calculation and by fitting the Mott's formula, taking the room temperature thermopower of biased (V 0 =1 eV) trilayer and quad-layer systems as two examples. The Mott's formula holds well at temperatures smaller compared to the Fermi energy. As shown in Fig. 10, the qualitative behavior of thermopower is reproduced very well by the Mott's formula even at the room temperature. Only close to the peak position, is the difference appreciable. From experiences gained on the monolayer and bilayer graphene, it is safe to say that impurity scattering would not enhance thermopower much. Thus, we conclude that biased BLG systems has the largest room temperature thermopower among all the multilayer graphene systems.

IV. SUMMARY
In this work, we systematically calculate thermopower of biased and unbiased multilayer grphene systems. The effect of screening of an external electric field is taken into account self-consistently under the Hartree approximation so that charge densities are different between inner and outer layers. Both the model with only nearest neighbor hopping and the more general SWMcC model with further neighbor hoppings are considered. The effect of impurity scattering is considered for monolayer and unbiased bilayer graphene in terms of the selfconsistent Born approximation. For monolayer graphene, only when the effect of impurity scattering is taken into account, we could obtain results consistent with experiments. The electronic structure and thermopower of biased multilayer graphene systems are calculated in which the screening effect is self-consistently incorporated. The biased bilayer graphene shows the largest room temperature thermopower.