Gate-tunable direct and inverse electrocaloric effect in trilayer graphene
Natalia Cortés, Oscar Negrete, Francisco J. Peña, Patricio Vargas
GGate-tunable direct and inverse electrocaloric effect in trilayer graphene
Natalia Cort´es, ∗ Oscar Negrete,
1, 2
Francisco J. Pe˜na, and Patricio Vargas
1, 2 Departamento de F´ısica, Universidad T´ecnica Federico Santa Mar´ıa, Casilla 110V, Valpara´ıso, Chile Centro para el Desarrollo de la Nanociencia y la Nanotecnolog´ıa, 8320000 Santiago, Chile (Dated: January 11, 2021)The electrocaloric (EC) effect is the reversible change in temperature and/or entropy of a materialwhen it is subjected to an adiabatic electric field change. Our tight-binding calculations linked toFermi statistics, show that the EC effect is sensitive to the stacking arrangement in trilayer graphene(TLG) structures connected to a heat source, and is produced by changes of the electronic density ofstates (DOS) near the Fermi level when external gate fields are applied on the outer graphene layers.We demonstrate the AAA-stacked TLG presents an inverse EC response (cooling), whereas the ECeffect in ABC-stacked TLG remains direct (heating) regardless of the applied gate field potentialstrength. We reveal otherwise the TLG with Bernal-ABA stacking geometry generates both theinverse and direct EC response in the same sample, associated with a gate-dependent electronicentropy transition at finite temperature. By varying the chemical potential to different Fermi levels,we find maxima and minima of the DOS are located near the extremes of the electronic entropy,which are correlated with sign changes in the differential entropy per particle, giving a particularexperimentally measurable electronic entropy spectrum for each TLG geometry. The EC effect inquantum two-dimensional layered systems may bring a wide variety of prototype van der Waalsmaterials that could be used as versatile platforms to controlling the temperature in nanoscaleelectronic devices required in modern portable on-chip technologies.
I. INTRODUCTION
The rapid increase in the need for efficient, environ-mentally friendly and low-cost materials with the capa-bility of cooling, has opened interesting possibilities toexplore material properties beyond the traditional vapor-compression method applied in household and industrialrefrigeration [1, 2]. As a clean and efficient alternative totuning the systems temperature, can be the implementa-tion of caloric effects on solid-state materials, driven byapplied pressure, electrical, magnetic and/or mechanicalfields, conducting to the barocaloric, electrocaloric (EC),magnetocaloric and elastocaloric effect, respectively [2–4]. Thermal response on solid systems is produced byeither of those external field acting on and/or removedfrom the sample under adiabatic conditions [5], leadingto large temperature changes ∆ T associated with isother-mal entropy changes ∆ S T , generally near phase transi-tions [3, 6].In the last decades, the EC effect—the change in tem-perature and entropy as an electric field is applied—has been widely studied in three-dimensional dielectricmultilayer capacitors (MLC), by considering a differentnumber of layers with diverse thicknesses in the samples,typically of the order of micrometers [7, 8], and find-ing thickness-dependent thermal responses [8]. In single-crystal MLC, the temperature change reaches ∼ . ∗ [email protected] nanometers, a giant EC response of ∼
12 K is manifestedclose phase transitions [11]. Theory predicts larger cool-ing power in thin-film MLC of ceramics and polymers byvarying the number of layers of the systems [12]. The fea-sible scalability of MLC systems to reduced dimensional-ity may allow expanding the EC effect to novel materialsfor electrical refrigeration at the atomistic level [13, 14].Nanoscale systems including van der Waals multilayersand heterostructures, comprise a wide variety of quantumsolids [15] that could be used as prototype miniaturizedmaterials to get caloric effects. A giant EC response of ∼
23 K has experimentally been achieved in a ferroelec-tric heterobilayer because of interface-induced interac-tions [16]. Theoretical studies on graphene nanoribbonsunder an electric and magnetic field show controllable en-tropy changes magnitude and entropy sign changes dueto the applied magnetic field [17]. First-principle andthermodynamic calculations in a two-dimensional (2D)monolayer of MoTe , have reported a temperature changeof 10 −
15 K near the structural phase transition that oc-curs when the monolayer is subject to electrostatic gating[18].We explore here the thermal response of atomicallythin layered nanostructures, focusing on gated trilayergraphene (TLG). Both experimental and theoreticalstudies have revealed the electronic band structure anddensity of states (DOS) of TLG samples strongly dependon the stacking pattern they possess, showing significantchanges under applied electrical potentials [19–27]. Thismakes TLG materials suitable to inspect the EC effectand novel thermodynamic quantities such as the differen-tial entropy per particle (DEP) [28–33]. Through Blochelectrons in parametrized tight-binding models linked toFermi statistics, we show how the thermal response varieswithin each TLG stacking as a result of the quantum- a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n thermodynamic processes involved at low-energies. Aseach system is gated, we find entropy changes for AAA-stacked TLG linearly increase with temperature, show-ing an inverse EC effect (cooling) up to room temper-ature. In contrast, the entropy changes for ABC TLGnon-linearly increase as temperature increases, giving adirect EC effect (heating). Remarkably, the ABA-TLGarrangement presents a dual EC response as cools andheats within the same sample, displaying gate-dependentmaxima for entropy changes. We also analyze the num-ber of electrons thermally excited of ABA-TLG at lowand high temperatures, unveiling its direct connection tothe gate-tunable electronic entropy and dual EC effect. II. QUANTUM-THERMODYNAMIC MODEL
We simulate TLG structures with simple hexagonalAAA, Bernal ABA, and rhombohedral ABC stackingorder by using a π -orbital tight-binding (TB) model,that captures the main low-energy symmetry propertiesaround the Fermi level for each system [20–22, 34–38].To construct the TB Hamiltonians, we consider the unitcells of the TLG with three pairs of atoms { A B } , { A B } and { A B } , which are located respectively in thebottom, central and top graphene layers, as shown in Fig.1. We use different parametrization (hoppings) depend-ing on the stacking type of the TLG. For each systemwe have first taken a minimal set of parameters fromfitting photoemission spectroscopy spectra to TB mod-els [39]. These parameters include a strong intralayerhopping γ between Ai and Bi nearest neighbors withineach graphene layer ( i = 1 , , γ , as used forAAA-stacked TLG in Fig. 1(a). For ABC and ABA-stacked TLG, this parametrization also incorporates aweaker nearest-layer coupling γ , Fig. 1(b) and 1(c) re-spectively. In particular, for ABA-TLG, we also havetaken a full set of parameters adopting the Slonczewski-Weiss-McClure parametrization [40], which contains oneextra nearest-interlayer hopping γ as well as two next-nearest interlayer hoppings γ and γ . An asymmetriconsite-potential energy difference with magnitude 2 V g ,induced by a bias voltage between the external graphenelayers, can act as an external gate potential on each car-bon atom of top + V g and bottom − V g graphene layersof the TLG structures, similar to bilayer graphene un-der an external voltage [41]. This asymmetric potentialis a crucial quantity to tune the electronic DOS and toproduce a caloric response near the Fermi level of TLGsystems.The Hamiltonians for the TLG samples ˆ H α ( kkk ) =ˆ H α ( kkk ) + ˆ V ( α = AAA , ABA , ABC) include a diago-nal matrix ˆ V of the gate potential proportional to V g .The matrices of hopping ˆ H α ( kkk ) can be obtained withthe matrix elements H αjj (cid:48) ( kkk ) = (cid:80) RRR e ikkk · RRR E αjj (cid:48) ( RRR ), where E αjj (cid:48) ( RRR ) = (cid:104) φ j ( rrr ) | ˆ H α | φ j (cid:48) ( r − Rr − Rr − R ) (cid:105) is the hopping integralbetween the atomic orbitals | φ j (cid:105) at and | φ j (cid:48) (cid:105) at latticevector RRR with j ( j (cid:48) ) = Ai, Bi . RRR = a C-C (1 , RRR = a C-C (cid:0) − / , +( − ) √ / (cid:1) are in-plane nearest-neighborvectors with a C-C = 1 .
42 ˚A the carbon-carbon distancewithin a graphene layer, and kkk = ( k x , k y , k z ) is the mo-mentum. The TLG Hamiltonians are represented in thebasis with components { ψ A , ψ B , ψ A , ψ B , ψ A , ψ B } .The density of states D ( E, V g ) is a fundamental gate-tunable electronic quantity, which allows us to calculatethe electronic thermodynamics properties of the TLGsystems. The DOS depends of both, the electronic levelwith energy E and the gate potential V g applied on theouter carbon atoms of each TLG. As the TLG structurescan be considered as quasi 2D systems, we numericallycalculate the DOS using a 2D Brillouin zone (BZ) in re-ciprocal kkk -space (i. e., k z = 0). We use a fine mesh ofabout ten million of kkk points in the area enclosed by thegreen triangle of Fig. 1(d), and for every kkk -state we eval-uate the energy levels coming from each band of ˆ H α ( kkk ).By calculating the DOS, we can obtain the number ofelectrons N in each unit cell of the trilayers N ( V g , T, µ ) = (cid:90) E h E l D ( E, V g ) n F ( E, T, µ ) dE, (1)where E l ( h ) is the lowest (highest) electronic energyeigenvalue of the considered Hamiltonian, n F ( E, T, µ ) =1 / [ e β ( E − µ ) + 1] is the Fermi-Dirac function distributionwith β = 1 /k B T , k B is the Boltzmann constant, µ is thechemical potential and T temperature. At T = 0 K, theelectronic levels are occupied up to the Fermi energy E F ,and n F it converts in the Heaviside function, so that weobtain N = (cid:82) E F E l D ( E, V g ) dE = 6 for the unit cell of eachTLG. For finite temperatures T >
0, we can calculatethe chemical potential µ ( V g , T ) by inversion of Eq. 1.The total entropy of the system S = S latt ( T ) + S e ( V g , T ) includes two terms, the entropy of the lattice S latt ( T ), giving account of the phonon contribution, inwhich we assume it is only dependent on temperatureand not on V g [42]. In our approximation we neglect S latt ( T ) and use S ∼ S e ( V g , T ), where the electronic en-tropy S e ( V g , T ) = − k B (cid:90) E h E l D ( E, V g ) F ( n F ) dE, (2)depends on the gate potential and temperature, and iscalculated in the triangular area of the BZ of Fig. 1(d),where the number of electrons is fixed to N = 6. In Eq.2, F ( n F ) = n F ln n F + (1 − n F ) ln(1 − n F ) , (3)its approximated by a Lorentzian-like function L ( E, T, µ ) = C/ [ e ( | E − µ | / k B T ) / + 1]. By consid-ering C = 1 .
4, low and high T values, we obtain A1 B1A2 B2A3 B3 A1 B1A2 B2A3 B3 γ γ γ γ γ γ γ γ γ k x k y k z K’K +V gate A1 B1 A3 B3 A2 B2 γ γ γ γ γ γ γ γ γ γ - V gate - V gate T +V gate +V gate - V gate γ TT (a) AAA (b) ABC (c) ABA (d) M 𝚪 γ FIG. 1. Schematic representation of trilayers graphene with diverse stacking patterns, hexagonal AAA (a), rhombohedral ABC(b) and Bernal ABA (c). Blue (red) sphere show A ( B ) carbon sublattice connected through intralayer γ and interlayerhoppings γ , γ , γ , γ and γ as indicated; ± V gate is the gate potential applied to the outer graphene layers on each system,the orange square enclosing each structure represent the thermal source at temperature T . (d) Shows the Brillouin zone inmomentum space for all trilayers graphene, where green lines highlight the triangular k -path with corners located at highsymmetry points Γ = (0 , K = πa (cid:0) , (cid:1) and K (cid:48) = πa (cid:0) √ , (cid:1) , whit a = √ a C-C = 2 .
46 ˚A the 2D graphene lattice constant. excellent agreement between Eq. 3 and L ( E, T, µ ) with − F ( n F ) ≈ L ( E, T, µ ), so that Eq. 2 transforms as S e ( V g , T ) (cid:39) k B (cid:90) E h E l D ( E, V g ) L ( E, T, µ ) dE. (4)The main contribution of L ( E, T, µ ) to the electronicentropy is given by their temperature-dependent width,which increases as temperature increases. Within therange of temperature we work, k B T (cid:28) γ , the chemicalpotential remains constant at µ = 0 eV in order to fulfill N = 6 for the unit cell of each TLG.The EC effect is in general quantified by either thetemperature change that occurs when the electric fieldchanges adiabatically, or the entropy change induced byisothermal application or removal of the electric field.Typically there are two ways to measure the caloric ef-fects in solid materials: ( i ) through the variation of tem-perature in an adiabatic thermodynamic path ∆ T (di-rect measurement), and ( ii ) by means of the entropychange in an isothermal path ∆ S T (indirect measure-ment), where the subindex T indicates constant tem-perature. It is important to note that is experimentallychallenging to obtain ∆ T as compared to ∆ S T , as directmeasurements generally require precision calorimetry [1].The EC effect in our TLG structures is then obtainedthrough the indirect way, by means of electronic entropychanges calculations as a function of T through Eq. 4 at µ = E F = 0, − ∆ S e,T ( T ) = S e ( V g = 0 , T ) − S e ( V g , T ) . (5)In case to obtain − ∆ S e,T >
0, we are in presence ofthe direct EC effect, that is the TLG system is capa-ble to heat as V g increases. In the opposite case when − ∆ S e,T <
0, the system present an inverse EC effectand so the sample cools down.
A. Differential electronic entropy per particle
To get essential insights into the microscopic originof the EC effect in TLG structures, we have calculateda fundamental thermodynamic quantity, the differentialelectronic entropy per particle s (DEP). The DEP is anexperimentally measurable quantity in gated 2D electronsystems, allowing superior sensitivity in comparison toa.c. calorimetry [28], and can be directly linked to thethermoelectric power or alternatively to the Seebeck co-efficient [43]. The DEP, s = ∂S e /∂N is found from theMaxwell relation as [28–33] s ( µ, V g , T ) = (cid:18) ∂S e ∂N (cid:19) T = − (cid:18) ∂µ∂T (cid:19) N , (6)and its related to the DOS through N in Eq. 1 by means s ( µ, V g , T ) = ( ∂N/∂T ) µ ( ∂N/∂µ ) − T .To interpret our results, we make a useful trans-formation of the form s ( µ, V g , T ) = ( ∂S e /∂N ) T =( ∂S e /∂µ ) T ( ∂µ/∂N ) T . The function that predominantlygoverns the DEP corresponds to the entropy’s derivativewith respect to µ . Consequently, we expect that at max-ima and minima of the electronic entropy as a functionof µ , the DEP vanishes, giving account of the thermody-namic fingerprints for each TLG arrangement.As the number of particles N can vary in adoped/gated system, we obtain a caloric effect directlyconnected to the DEP, which could be used as a feasi-ble procedure to determine the entropy changes given byfluctuating charge carriers, see Appendix for more de-tails. III. THE ELECTROCALORIC RESPONSE INTRILAYERS GRAPHENEA. AAA-stacked TLG
AAA-stacked TLG has been recently exfoliated [39],preserves metallic character in the presence of externalelectric fields [44] as in the monolayer graphene case, andpossesses mirror reflection symmetry with respect to thecentral layer. In this structure the { Ai , Bi } sublatticesmatch to the nearest-neighbor layer carbon sites { Ai +1, Bi + 1 } with vertical hopping γ . The whole systemeffectively can be seen as a direct superposition of threegraphene monolayers as shown in Fig. 1(a). Because ofthe crystal symmetry, the AAA-TLG Hamiltonianˆ H AAA ( kkk ) = − V g γ f γ f γ f ∗ − V g γ f γ f ∗ γ f γ f γ f ∗ γ f ∗ γ f γ f ∗ V g γ f γ f ∗ γ f ∗ V g , (7)allows analytical eigenvalues connecting monolayer-graphene modes ε , and AA-bilayer graphenelike modes ε , , , ε , = ± γ | f | , (8a) ε , = ± γ | f | − (cid:113) V g + 2 γ | f | , (8b) ε , = ± γ | f | + (cid:113) V g + 2 γ | f | . (8c) f ( k x , k y ) = e − ik x a C-C + 2 cos( √ k y a C-C ) e ik x a C-C and f ( k z ) = e ik z c are respectively an in-plane and out-of-plane momentum-dependent functions with c = 3 . K point for different values of the gate potential0 ≤ V g ∼ γ . Due to the structural symmetry, theDirac points K , K (cid:48) are degenerate, where around eachone of them there are three pairs of linear branches cre-ating a diamond structure with Dirac cones throughoutenergy-momentum space [44]. At the Fermi level E F = 0and momentum K , or charge neutrality point (CNP),one of the Dirac cones remains constant regardless of the V g strength, resembling the monolayer graphene modes[Eq. 8(a)]. When V g = 0 (thin dark blue solid lines)and because of interlayer interaction γ , a pair of Diraccones shift away from the CNP in a symmetric way by E = ±√ γ (cid:39) .
25 eV. As V g is turned on, a pair of conesdisperse to E = ± (cid:112) V g + 2 γ , Eq. 8(b) and 8(c). Theeffect of the gate potential therefore can be represented asa renormalization of the interlayer hopping energy [45],preserving metallic character for the bands as in the un-perturbed case ( V g = 0). At low-energies, the DOS pre-serves electron-hole symmetry [ D ( E, V g ) = D ( − E, V g )]in Fig. 2(b), and shows a global minimum at the CNP, G M ( a ) K Energy [eV] V g = 0 e V V g = 0 . 0 5 e V V g = 0 . 2 0 e V V g = 0 . 4 0 e V k [ 1 / (cid:1) ] ( b ) D O S [ 1 / e V ]
FIG. 2. Low-energy electronic spectra for AAA-stacked tri-layer graphene. (a) Band structure in the vicinity of the K point. (b) Density of states obtained for the triangular areaenclosed by green lines in Fig. 1(d). Vertical dashed line in(a) highlight the K point, horizontal dashed lines in (a)-(b)are fixed at E F = 0 eV. The hopping parameters are γ = 3 . γ = 0 .
18 eV [39]. which causes the major contribution to the EC responseat low and high temperatures. At the energies given bythe shifted Dirac cones from the Fermi level, the DOSshows two symmetric discontinuities, not providing statesto the EC effect as the L function in Eq. 4 does not cap-ture those states up to room temperature.The EC effect (or entropy changes) for AAA-stackedTLG is presented in Fig. 3(a). The − ∆ S e,T response islinear as a function of T , following the linear behavior ofthe electronic spectra near E F of Fig. 2. The EC effect isinverse for all V g and the entire temperature range from0 to 300 K, and reaches a high value (cid:39) − . µ eV/Kat room temperature for large gate voltage ( V g = 0 . T = 30 K and Fig.4(b) T = 300 K. The electronic entropy at T = 30 K al-most preserves the DOS shape of Fig. 2(b) as the L func-tion slowly smears the DOS at low temperature, while at T = 300 K the entropies are largely smoothed with onemagnitude order higher than T = 30 K because of higher T . The electronic entropy at µ = 0 and V g = 0 for bothtemperatures T = 30 K and T = 300 K is the smallestentropy with respect to the other entropies with larger V g . This means S e ( V g = 0 , T ) < S e ( V g (cid:54) = 0 , T ), hence wealways obtain − ∆ S e,T < V g we have consideredhere.The DEP as a function of µ is shown in Fig. 4(c) and4(d) for the same temperatures as for the electronic en-tropy calculations. The DEP for AAA TLG its an oddfunction of µ , that is s ( − µ ) = − s ( µ ), and at µ = 0 van-ishes for all gate fields. At low temperature in Fig. 4(c),the DEP shows a stairslike shape, mainly determined bythe slope ( ∂S e /∂µ ) T , with maxima (minima) localizednear the discontinuities of the electron (hole) regime ofthe electronic entropy. The states of the maxima andminima of the DEP are around the shifted Dirac conesof the electronic dispersion of Fig. 2(a), which do not con-tribute to the EC effect as we mention above. At roomtemperature in Fig. 4(d), the DEP is smoothed and onemagnitude order higher than T = 30 K, as also seen forthe electronic entropy, demonstrating a general linear be-havior for the electronic and thermodynamics quantitiesin AAA-stacked TLG. V g = 0 . 0 5 e V V g = 0 . 1 0 e V V g = 0 . 2 0 e V V g = 0 . 4 0 e V -D S e ,T [ (cid:2) (cid:1) eV/ K] A B C s t a c k i n g( b ) -D S e ,T [ (cid:2) (cid:1) eV/ K] T [ K ] ( a ) A A A s t a c k i n g V g = 0 . 0 5 e V V g = 0 . 1 0 e V V g = 0 . 2 0 e V V g = 0 . 4 0 e V T [ K ]
FIG. 3. The isothermal entropy change − ∆ S e,T as a func-tion of T for different gate potentials as indicated. (a) AAA-stacked and (b) ABC-stacked trilayer graphene. The chemicalpotential is set to µ = 0. Dashed vertical lines indicate T = 30K and T = 300 K, where we calculate the electronic entropyand DEP in Fig. 4 and Fig. 6. - 0 . 4 - 0 . 2 0 . 0 0 . 2 0 . 40 . 0 20 . 0 40 . 0 60 . 0 8 T = 3 0 0 KT = 3 0 0 K ( a ) T = 3 0 K V g = 0 e V V g = 0 . 0 5 e V V g = 0 . 2 0 e V V g = 0 . 4 0 e V S e ,T [ (cid:2) (cid:1) eV/K] - 0 . 4 - 0 . 2 0 . 0 0 . 2 0 . 40 . 20 . 40 . 60 . 8 - 0 . 4 - 0 . 2 0 . 0 0 . 2 0 . 4- 3- 2- 10123 T = 3 0 K DEP [ (cid:1) eV/K] m [ e V ] - 0 . 4 - 0 . 2 0 . 0 0 . 2 0 . 4- 3 0- 2 0- 1 001 02 03 0 ( d )( c ) ( b ) m [ e V ] FIG. 4. Electronic entropies (a)-(b) and differential entropyper particle (c)-(d) as a function of the chemical potential forAAA-stacked trilayer graphene. Left panels stand for T = 30K, right panels T = 300 K. Vertical dashed lines indicate µ = 0. B. ABC-stacked TLG
Rhombohedral ABC-stacked trilayer graphene is oneof the commonly stable crystal obtained in experimen-tal procedures [24, 39, 46], which can be considered asa zero-gap semiconductor material in the unperturbedcase, and as a semiconductor when external electric-fieldpotentials are applied [19]. This stacking geometry pos-sesses inversion symmetry, but lacks mirror symmetry[24] as is schematized in Fig. 1(b). The TB Hamiltonianconsidering the hopping interactions in Fig. 1(b) is givenbyˆ H ABC ( kkk ) = − V g γ f γ f γ f ∗ − V g γ f γ f ∗ γ f γ f γ f ∗ γ f ∗ γ f
00 0 0 γ f ∗ V g γ f γ f ∗ γ f ∗ V g , (9)where f ( kkk ) = f ( k x , k y ) e ick z . The band structure andDOS from Eq. 9 considering 0 ≤ V g ≤ γ are plottedin Fig. 5(a) and Fig. 5(b) respectively. The electronicspectra preserve electron-hole symmetry around E F = 0for all V g . In the electronic structure there are valenceand conduction bands dispersing according to the hop-pings between the different atoms. For V g = 0 (thinnestgray lines), the lower conduction band and higher valenceband touch at E F near the K point, where the DOS showsa local maximum. Furthermore when V g = 0, two other ( b )( a ) Energy [eV] V g = 0 e V V g = 0 . 0 5 e V V g = 0 . 2 0 e V V g = 0 . 4 0 e V k [ 1 / (cid:1) ] D O S [ 1 / e V ] G M K FIG. 5. (a) Band structure near the K point, (b) DOS fromthe green triangle area in Fig. 1(d) for ABC-stacked trilayergraphene with diverse gate potential V g as indicated. Hori-zontal (vertical) dashed lines indicate E F = 0 eV ( K point).The hopping parameters are γ = 3 .
10 eV, γ = 0 . γ = 0 . pairs of bands shift away from E F and are degenerate,crossing near E (cid:39) ± . γ between the direct bonded B A B A T >
900 K. When V g is turned on, an energy gap between the lower elec-tron (conduction band) and higher hole branch (valenceband) opens at the K point as well as in their vicinitybecause of the potential energy difference between thebottom and top graphene layers[19, 22, 36, 44, 47, 48].The gaps non-linearly increase as V g increases [27, 44],the DOS vanishes for the gap energies as expected, andshows two high symmetric peaks (van Hove singularities)near the valence and conduction band edges, whose statesare mainly contributed by the carbon atoms of the topand bottom graphene layers [27]. The bandgap edges forlow-gate potential V g = 0 .
05 eV ( V g = 0 . ∼
80 K (140 K), while the states contribute for
T > V g = 0 . , . S e ( V g = 0)] is higher than the other gate-dependententropies, because the states of the DOS evaluated at D ( E F = 0 , V g = 0), are fully capturated by the L function at low and high temperatures. We also havechecked that there is not a direct-inverse transition[ − ∆ S e ( V g , T ) = 0] up to T = 500 K. Moreover, an in-teresting effect occurs for temperatures below 80 K forall V g , where the entropy changes are constant and donot present differences in the caloric response because ofthe presence of the band gap. As the temperature in-creases in Fig. 3(b), the maxima of the entropy changescorrespond to − ∆ S e ( V g = 0 .
05 eV , T = 137 K) (cid:39) . µ eV/K and − ∆ S e ( V g = 0 . , T = 236 K) (cid:39) . µ eV/K. For high V g and T > − ∆ S e equally increasesfor V g = γ / . V g = γ = 0 . (cid:39) . µ eV/K at room temperature.These behaviors are consistent with the electronic en-tropy as a function of µ in Fig. 6(a) at T = 30 K and Fig.6(b) T = 300 K. For both temperatures at µ = 0, theelectronic entropy fulfills the condition S e ( V g = 0 , T ) >S e ( V g (cid:54) = 0 , T ), giving the direct EC response, then theABC-TLG sample heats, opposite to the AAA-TLG case.At low temperature in the vecinity of µ = 0 in Fig. 6(a),the only entropy contributing to the EC effect is thatfor V g = 0, while the entropies when V g (cid:54) = 0 vanishand do not provide states to the EC effect as expected.For temperatures lower than 80 K, we have verified that − ∆ S e ( V g (cid:54) = 0) (cid:39) S e ( V g = 0) because of the band gap. Atroom temperature in Fig. 6(b), the electronic entropiesare nearly five times larger than entropies at T = 30 K,showing non null values for all V g in the gap region. Theentropy is almost flat for V g = 0, while the other entropieshave a parabolic shape as V g increases, with almost thesame values for high gate-potential fields V g = 0 . , . µ = 0, as expected from the EC response in Fig.3(b).The DEP of ABC-stacked TLG also is an odd functionof µ at low and high temperatures, as shown in Fig. 6(c)and Fig. 6(d) respectively. At T = 30 K the peak value - 0 . 2 - 0 . 1 0 . 0 0 . 1 0 . 20 . 0 00 . 0 80 . 1 6 ( d )( c ) ( b )( a ) T = 3 0 0 KT = 3 0 0 KT = 3 0 K V g = 0 e V V g = 0 . 0 5 e V V g = 0 . 2 0 e V V g = 0 . 4 0 e V S e ,T [ (cid:1) eV/K] - 0 . 2 - 0 . 1 0 . 0 0 . 1 0 . 20 . 00 . 40 . 8 - 0 . 2 - 0 . 1 0 . 0 0 . 1 0 . 2- 4- 2024 T = 3 0 K
DEP [meV/K] m [ e V ] - 0 . 2 - 0 . 1 0 . 0 0 . 1 0 . 2- 0 . 4- 0 . 20 . 00 . 20 . 4 m [ e V ] FIG. 6. Electronic entropies (a)-(b) and differential entropyper particle (c)-(d) as a function of µ for ABC-stacked trilayergraphene. Left panels stand for T = 30 K, right panels T =300 K. Vertical dashed lines indicate µ = 0, and vertical dash-dotted lines in (a)-(b) at µ = ± .
127 ( µ = ± .
18) eV highlightthe maxima at T = 30 ( T = 300) K for S e,T ( V g = 0 . of | s ( µ (cid:39) | is three magnitude orders lower than that ofAAA TLG, and increases as the gate field increases. Asthe EC effect in Fig. 3(b) at T = 30 K has not dependenceon V g , the DEP presents sizable dip-peak structures inthe band gap zone, with large slopes in the vicinity of µ = 0 for V g (cid:54) = 0, while s ( µ ) = 0 for V g = 0 except near E F . The chemical potential dependence for s ( µ ) giveszeros when S e,T ( µ ) reach maxima, as indicated with ver-tical dash-dotted lines particularly for V g = 0 . s ( µ ) and maxima of S e,T ( µ )are located near the band edge Van Hove singularities ofthe electronic spectra in Fig. 5. This behavior is similarto the thermopower of biased bilayer graphene [41] andother gapped 2D materials such as germanene under anexternal electric field [30] and semiconducting dichalco-genides [31]. At room temperature, the DEP in Fig. 6(d)decreases one magnitude order along with smoother reso-nances, showing a sinelike shape about E F for all V g (cid:54) = 0.Near µ = 0 the slopes decrease as the gate potential de-creases. The slope reverses for V g = 0 because of themaximum of S e,T ( µ = 0) instead a minimum. These re-sults are in agreement with the EC response in Fig. 3(b). C. ABA-stacked TLG
Trilayer graphene with Bernal-ABA stacking is themost stable geometry as resembles a graphitelike atomicstructure [19]. The ABA TLG possesses mirror symme-try with respect to the middle graphene layer. This struc-tural symmetry is broken in the presence of electrostaticpotentials [21, 47], driving to semimetallic character withtunable overlap of the electronic bands [19, 23, 27], as wellas controllable gaps at K point and near it [20, 27, 44].The TB Hamiltonian considering all hoppings in Fig. 1(c)is given byˆ H ABA ( kkk ) = − V g γ f γ f γ f γ f γ f ∗ − V g γ f γ f γ f γ f ∗ γ f ∗ γ f γ f ∗ γ f γ f ∗ γ f ∗ γ f ∗ γ f ∗ γ f ∗ γ f ∗ γ f γ f V g γ f γ f ∗ γ f ∗ γ f γ f ∗ V g . (10)We analyze first a simple case for ABA TLG, labelednearest-layer (NL) TB model, including hopping param-eters γ , γ and γ , with γ = γ = γ = 0 from spec-troscopy experiments [39]. We take different values ofthe external gate potential in the limit 0 ≤ V g (cid:38) γ forthe electronic spectra and EC calculations.Figure 7(a) show the electronic band structure around K point, Fig.7(b) correspond to the DOS. Both quan-tities are shown in the vicinity of E F = 0. The energyspectra are symmetric about zero energy with and with-out the presence of the gate field, preserving electron-hole symmetry [21]. For V g = 0 (solid black lines), twobands mimic the k -linear electronic dispersion of mono-layer graphene at the K point, and one pair of parabolicbands overlap along KM k -axis near E F . The whole un-perturbed bands ( V g = 0) give a superposition of statesof monolayerlike and Bernal-bilayerlike graphene states[34, 35, 47]. The DOS for V g = 0 shows a sharp mini-mum at E F with two symmetric peaks surrounded it [seeinset of Fig.7(b)]. The major contribution of the edgepeaks its from the AB-bilayerlike electron and hole bandedges, whose states are fully captured by the L functionfor T >
10 K.The gate field breaks mirror reflection symmetry withrespect to the central layer and mixes the linear andparabolic bands in the region of E F [21, 47], so that when V g = 0 .
05 eV (red dotted lines) and V g = 0 . E F , and the bilayerlike bands remain almost unper-turbed. The DOS for these cases presents a minimumat E F , and shows large electron and hole edge peaks asmore states are overlapped near E F , where the electronicstates contribute to the EC effect for T >
30 K. Whenthe gate potential increases comparable to γ ∼ V g = 0 . γ (cid:39) V g = 0 . K point as wellas along the KM direction near E F , while a gap opensalong K Γ as γ is non zero [20, 39]. For these large-gate field potentials, the DOS almost vanishes at E F andshows one pair of symmetric peaks centered near the gapedges, contributing to the EC response when T ≥
70 K.Other pair of peaks are centered in the vicinity of theminima of the highest valence band and maxima of thelowest conduction band, not providing states to the ECeffect for
T <
170 K.When all the interactions between the atoms are takeninto account in the ABA TLG system, labeled as next-
FIG. 7. Electronic spectra for ABA-stacked trilayer graphenenear the Fermi level using nearest layer NL (a)-(b) and next-nearest layer NNL (c)-(d) TB models. Left panels show theelectronic band structure about K point in the Brilloin zone.Right panels ilustrate the density of states for the triangulararea enclosed by green lines in Fig. 1(d). (a)-(b) γ = 3 . γ = 0 .
39 eV and γ = 0 .
25 eV [39]. (c)-(d) γ = 3 . γ = 0 .
39 eV, γ = − .
020 eV, γ = 0 .
315 eV, γ = 0 . γ = − .
04 eV [36]. Horizontal (vertical) dashed lineindicates E F ( K point). Inset in (b) shows a zoomed area forthe spectra of the NL TB model around E F = 0. nearest layer (NNL) TB model, the energy spectra areno longer electron-hole symmetric about E F as in the NLcase because of the additional electron hoppings. How-ever, the full band structure for V g = 0 can still be con-sidered as a combination of the band diagram of a sin-gle graphene layer and that of AB-bilayer graphene asshown in Fig. 7(c). At K point, the valence and con-duction parabolic bands open a gap near E F for V g = 0,where the DOS shows a flat minimum in the hole zone( E F (cid:46) T >
60 K. When V g = 0 .
05 eV (dotted red lines),the gate potential causes anticrossing of the parabolicbands, similar as seen when negative and positive low-charged gates act on a suspended ABA TLG [20]. At E F = 0 the DOS increases, showing nearly a maximuminstead of a minimum as in the NL case. The electronicstates belonging to this maximum are fully captured bythe L function for T >
30 K. For V g = 0 . K point, and showing a highly antisymmetric DOS, whosestates near E F contribute to the EC effect for T > V g = 0 . V g = 0 . K point as well as along both K Γ and
KM k -axis, similar to the effect of gate-inducedhigh charge density on ABA TLG [20]. The DOS forthese high V g values show minima at zero energy, as inthe band structure there are allowed states only near the K point, contributing to the EC response for T ≥
30 K.In Fig. 8(a) and 8(b), we present respectively theEC response for both the NL-ABA and NNL-ABA TBmodel. Both cases are qualitatively similar as they showa dual behavior seen as a combination of the directand inverse EC effect, whose gate-dependent inversion − ∆ S e,T ( V g ) = 0 occurs at relatively low temperatures.Moreover, noticeable differences are observed for the en-tropy changes depending on the magnitude of the appliedgate field in each ABA-TB case. We observe a small di-rect response zone (heating) for low gate-fields V g = 0 . V g = 0 . − ∆ S e ≤ . µ eV/K at T ≤
12 –see insets of Fig. 8. As tempera-ture increases for these low-gate field values, the entropychanges become inverse for T ≤
21 and remain constantas temperature increases, cooling down up to T = 300 K.The EC effect for high-gate fields V g = 0 . , . , . − ∆ S e (cid:39) . µ eV/K at T (cid:39)
23 K. The entropy inversion tempera-ture also is nearly the same for all high V g , occurringat T ∼
45. The entropy changes maxima for high-gatefields increase almost twice in the ABA-NNL model, withhigher temperatures in the range 26 K ≤ T ≤
30 K. Thisis because the L function captures a maximum insteadof a minimum of the reference DOS [ D ( V g = 0)] in theNNL-ABA model. Remarkably, the inversion tempera-ture in the ABA-NNL model varies with high-gate fields,56 K ≤ T ≤
69 K. As the temperature increases, the en-tropy changes increase up to room temperature, coolingmore for high gate field values. The inversion tempera-tures for the ABA-NNL system also can be obtained inzone II of Fig. 10(a), where the reference entropy (blackline) crosses with all other entropies.The direct-inverse response in ABA-TLG is capturedby the electronic entropy as a function of µ in Fig. 9(a)and 9(b). We present here only the ABA-NNL plots asboth ABA system EC results are equivalent. Figure 9(a)shows S e,T at T = 20 K resembles the DOS shape of Fig.7(d) with broken electron-hole symmetry. At µ = 0, thereference entropy for V g = 0 is lower than the entropy for V g = 0 .
05 eV, and larger than the entropies for all othergate field values, see Fig. 10(a). Consequently, the samplecools for V g = 0 .
05 eV, and heats for higher gate fields at T = 20 K, as expected from the NNL-EC response in Fig.8(b). When the temperature increases up to T = 150 Kin Fig. 9(b), the electronic entropies are nearly four timesthe entropies at T = 20 K, and are smoothed becauseof higher temperature. The reference entropy shows aflat line shape near µ = 0 and is smaller than all otherentropies, indicating the system cools for all gate fieldsat T = 150, as also seen from the EC response.Figure 9(c) and 9(d) show the DEP for the ABA-NNL V g = 0 . 0 5 e V V g = 0 . 1 0 e V V g = 0 . 2 5 e V V g = 0 . 3 0 e V V g = 0 . 4 0 e V -D S e ,T [ (cid:2) (cid:1) eV/ K] T [ K ] ( b )( a )
A B A - s t a c k i n g N N LA B A - s t a c k i n g N L
T [ K ] -D S e ,T [ (cid:2) (cid:1) eV/ K] T [K ] -D S e ,T [ (cid:2) (cid:1) eV/ K] T [K ]
FIG. 8. The isothermal entropy changes − ∆ S e,T as a functionof temperature T for different gate potentials as indicated.(a) ABA-stacked trilayer graphene considering NL hoppings,(b) ABA-stacked trilayer graphene with NNL hoppings. Thechemical potential is set to µ = 0 eV. Horizontal dashed lineindicates − ∆ S e,T = 0, while vertical dashed lines stand for T = 20 K and T = 150 K, where the entropies are calculatedas a function of µ in Fig. 9. Insets show a zoom for eachABA-EC response at relatively low temperatures. - 0 . 2 - 0 . 1 0 . 0 0 . 1 0 . 20 . 0 00 . 0 40 . 0 8 V g = 0 e V V g = 0 . 0 5 e V V g = 0 . 1 0 e V V g = 0 . 3 0 e V V g = 0 . 4 0 e V T = 1 5 0 KT = 1 5 0 K ( b )( a )
T = 2 0 KT = 2 0 K S e ,T [ (cid:2) eV/K] - 0 . 2 - 0 . 1 0 . 0 0 . 1 0 . 20 . 00 . 20 . 4 - 0 . 2 - 0 . 1 0 . 0 0 . 1 0 . 2- 1 6 0- 8 008 01 6 0 ( d ) m [ e V ] DEP [ (cid:2) (cid:1) eV/K] m [ e V ] - 0 . 2 - 0 . 1 0 . 0 0 . 1 0 . 2- 8 0- 4 004 08 0 ( c ) FIG. 9. Electronic entropies (a)-(b) and differential entropyper particle DEP (c)-(d) as a function of µ for ABA-stackedtrilayer graphene considering the NNL-TB model. Left panelsstand for T = 20 K, right panels T = 150 K. Vertical dashedlines indicate µ = 0 for all plots. Vertical dash-dotted lines at µ = − .
083 eV and µ = 0 . µ = − .
065 eV, µ = 0 .
083 eV in (b) highlight the maxima for S e,T ( V g = 0 . model at T = 20 K and T = 150 K respectively. s ( µ ) isnot a symmetric function of µ because of broken electron-hole symmetry for all gate-potential fields. Near zerochemical potential and V g = 0 , . , . ( b )( a ) V g = 0 e V V g = 0 . 0 5 e V V g = 0 . 1 0 e V V g = 0 . 2 5 e V V g = 0 . 3 0 e V V g = 0 . 4 0 e V S e ,T [ (cid:1) eV/ K] T [ K ]
I I I I I I n th[1012cm-2] T [ K ]
FIG. 10. (a) Electronic entropy S e ( T ), (b) density of ther-mally excited electrons n th ( T ) versus temperature calculatedat E F = 0 in ABA-stacked TLG within the NNL-TB modelfor different gate potentials V g . In (a), I: S e ( V g = 0) ≥ S e ( V g (cid:54) = 0) , II: S e ( V g = 0 . , . ≥ S e ( V g = 0) ≥ S e ( V g = 0 . , . , . S e ( V g = 0) ≤ S e ( V g (cid:54) = 0).Zone I, II and III in (b) corresponds to the same temperaturerange as in (a). There are two other smaller resonances when µ lies inthe vicinity of the discontinuities of the electronic en-tropy, shifted to larger µ as V g increases. These entropydiscontinuities at 20 K are located near the minima ofthe valence band and maxima of the conduction band ofthe electron dispersion. s ( µ ) decreases for low gate fields V g = 0 . , . E F , near where theEC response displays a sign change at ∼
21 K [see insetof Fig. 8(b)]. As the ABA TLG gets cool for all V g at150 K [Fig. 8(b)], the DEP in Fig. 9(d) reduces by almosthalf near µ = 0, and the resonances are now spread out,mixed with the shifted ones close to S ( µ ) discontinuities.Because the ABA-stacked TLG is capable to heats andcools within the same sample as a direct and inverse ECresponse are obtained in Fig. 8, it is interesting to knowhow the electronic entropy S e is related to the chargedensity at the microscopic level. We calculate the num-ber of carriers thermally excited of ABA-stacked TLGthrough the NNL-TB model. We determine the densityof excited electrons n th from the valence band to the con-duction band of Fig. 7(c) as a function of temperature forselected values of the gate potential. In a similar way toEq. 1, n th above the Fermi level at a given temperature isgiven by n th ( V g , T ) = A − (cid:82) E h µ D ( E, V g ) n F ( E, T, µ ) dE ,where A is the unit cell area of the TLG, and we take µ = 0 as the Fermi energy at T = 0.The temperature dependence of S e ( T ) and n th ( T ) isshown in Fig. 10(a) and Fig. 10(b) respectively. The frac-tion of excited electrons follows a very similar pattern asthe electronic entropy because charge carriers contributeto the entropy by filling electronic states as temperatureincreases. In Fig. 10(b), electrons are excited at lowtemperature for low gate potentials V g = 0 . , . E F = 0, see Fig.7(c) and 7(d). As the gate potential increases, the en-ergy gaps of ABA-TLG increase, and a higher tempera- ture is needed to excite electrons from the valence bandto the conduction band. This produces a suppression(non-linear behavior) of n th at low temperatures for highelectric potentials. From ∼
70 K, n th linearly increaseswith temperature at high temperatures, which may corre-spond to a monotonic DOS as a function of energy, thatis gapless parabolic electronic bands near E F . Similarbehavior for the density of excited electrons has been ex-perimentally observed in graphene multilayers includingAB-bilayer and ABA-TLG, correlated with a theoreticaltemperature-dependent gap approximation, where bothprocedures give account of a finite-temperature electronicphase transition [49]. In our case, large gate-field values V g = 0 . , . , . E F . This causes thedirect EC effect − ∆ S e > − ∆ S e < S e ( T ), whereelectrons are thermally excited in a linear way as a func-tion of temperature, indicating gapless bands near E F . IV. CONCLUSIONS
Gated trilayer graphene structures show an elec-trocaloric response that depends on the stacking patternand the electronic character they possess near the Fermilevel. Trilayer graphene with AAA stacking remainsmetallic under the applied gate field, presenting largeentropy changes as the gate potential increases, linearlycooling up to room temperature. In striking contrast,ABC-stacked trilayer graphene converts into a semicon-ductor when the gate field is turned on, non-linearlyheating as temperature increases. The energy spectraof ABA-stacked TLG geometry present semimetallic aswell as gapped character with gate voltage, cooling andheating in the caloric response in a linear and non-linearmixed way, giving a combination of AAA and ABC elec-trocaloric effect. We verified that gate-dependent ther-mally excited charge at E F = 0 for ABA-TLG can be di-rectly linked to the electronic entropy as well as to the ECresponse, where charge suppression is related to the di-rect EC effect, and linearly excited electrons correspondto the inverse EC effect.The differential entropy per particle shows particu-lar signatures for each TLG structure at low and hightemperatures. Measuring this thermodynamic quan-tity could be used as a cheaper and novel experimen-tal technique to identify the stacking pattern in trilayersgraphene and-or similar layered 2D materials. It is worthmentioning that a new caloric effect based on the entropyper particle and EC effect connection can be obtained,see Appendix. We refer to this mechanism as charge-carrier caloric effect , as the variable quantity are chargecarriers instead of an external field. We propose this ef-0fect may be measured by varying the charge density andkeeping constant the electric field in the sample, similarto a 2D electron gas experimental setup [28]. V. ACKNOWLEDGMENTS
N.C. acknowledges support from ANID Fondecyt Post-doctoral Grant No. 3200658, P.V. and O.N. acknowledgesupport from ANID PIA/Basal AFB18000. F.J.P. ac-knowledges support from ANID Fondecyt, Iniciaci´on enInvestigaci´on 2020 grant No. 11200032, and the financialsupport of USM-DGIIE.
APPENDIX: CHARGE-CARRIER CALORICEFFECT
The connection between the electrocaloric effect andthe differential entropy per particle s defined in Sec. II A,can be obtained by analyzing the Helmholtz free energy F and internal energy U of a general system that includesa variable number of particles N in its formulation.First, we define the differential internal energy of thesystem as dU = T dS + E · d P + µdN, (11)where the first term corresponds to heat δQ = T dS , with S the total entropy; the second term is the definitionof electric work W = (cid:82) E · d P , where P is the electricpolarization and E is the electric field. The last term µdN , where µ is the chemical potential, is the energyneeded to add or remove one particle of the system. Onthe other hand, F is defined as F = U − T S − E · P , (12)and its differential is dF = − SdT − P · d E + µdN, (13)where we have used Eq. 11. From Eq. 13, we can obtainthe following Maxwell relations: (cid:18) ∂S∂ E (cid:19) T,N = (cid:18) ∂ P ∂T (cid:19) E ,N , (14)and (cid:18) ∂S∂N (cid:19) T, E = − (cid:18) ∂µ∂T (cid:19) E ,N , (15) where the last equation is the same as presented in Eq.6, but including the thermodynamic variable associatedwith the electric field, and the subscripts refer to vari-ables held constant during partial differentiation.To take into account the electronic entropy variationsin our system, it is convenient to specify the total differ-ential entropy, whose expression is given by dS = (cid:18) ∂S∂T (cid:19) E ,N dT + (cid:18) ∂S∂ E (cid:19) T,N d E + (cid:18) ∂S∂N (cid:19) T, E dN. (16)If we calculate the entropy change in an isothermal pro-cess ( dT = 0) along a path of the form S ( E , N ) → S ( E , N ) → S ( E , N ), we obtain∆ S T = (cid:90) E E (cid:18) ∂S∂ E (cid:19) T,N = N d E + (cid:90) N N (cid:18) ∂S∂N (cid:19) T, E = E dN. (17)Employing the Maxwell relations of Eq. 14 and Eq. 15,we finally obtain∆ S T = (cid:90) E E (cid:18) ∂ P ∂T (cid:19) E ,N d E − (cid:90) N N (cid:18) ∂µ∂T (cid:19) N, E = E dN. (18)This last equation occurs for an isothermal process inwhich the electric field and the number of particles varyin the process. The terms associated with the DEP corre-sponds to the integrand in the second term of Eq. 17 andEq. 18. This means that the DEP is directly connected tothe electrocaloric phenomena. The charge-carrier caloriceffect will correspond to the measurement of a process inwhich the electric field remains fixed, but a variation inthe number of particles is present in the initial and fi-nal stages of the process. Hence, the pure caloric effectobtained from the DEP is given by∆ S DEP T = − (cid:90) N N (cid:18) ∂µ∂T (cid:19) N, E dN. (19)It is essential to highlight that when the number of par-ticles remains fixed, the usual expression of the elec-trocaloric effect is recovered from Eq. 18∆ S T = (cid:90) E E (cid:18) ∂ P ∂T (cid:19) E ,N d E . (20)As the differential entropy per particle is an experimen-tally measurable quantity, we believe that the associated charge-carrier caloric effect can be performed through in-direct measurements, not requiring precision calorimetry.A similar procedure has been recently reported by a di-rect measurement of ∂µ/∂T , obtaining very small electriccurrents, caused by a periodic temperature modulationin a 2D diluted electronic system [28]. [1] J. Scott, Annu. Rev. Mater. Res. , 229 (2011). [2] J. Shi, D. Han, Z. Li, L. Yang, S.-G. Lu, Z. Zhong, J. Chen, Q. Zhang, and X. Qian, Joule , 1200 (2019).[3] X. Moya, S. Kar-Narayan, and N. D. Mathur, Nat.Mater. , 439 (2014).[4] B. Li, Y. Kawakita, S. Ohira-Kawamura, T. Suga-hara, H. Wang, J. Wang, Y. Chen, S. I. Kawaguchi,S. Kawaguchi, K. Ohara, et al. , Nature , 506 (2019).[5] M. Valant, Prog. Mater. Sci. , 980 (2012).[6] M. Otoniˇcar and B. Dkhil, Nat. Mater. , 9 (2020).[7] S. Kar-Narayan and N. Mathur, J. Phys. D: Appl. Phys. , 032002 (2010).[8] X. Moya, E. Defay, N. D. Mathur, and S. Hirose, MRSBull. , 291 (2018).[9] X. Moya, E. Stern-Taulats, S. Crossley, D. Gonz´alez-Alonso, S. Kar-Narayan, A. Planes, L. Ma˜nosa, andN. D. Mathur, Adv. Mater. , 1360 (2013).[10] B. Nair, T. Usui, S. Crossley, S. Kurdi, G. Guzm´an-Verri,X. Moya, S. Hirose, and N. Mathur, Nature , 468(2019).[11] A. Mischenko, Q. Zhang, J. Scott, R. Whatmore, andN. Mathur, Science , 1270 (2006).[12] S. Kar-Narayan and N. Mathur, Appl. Phys. Lett. ,242903 (2009).[13] S. Lisenkov and I. Ponomareva, Phys. Rev. B , 140102(2009).[14] I. Ponomareva and S. Lisenkov, Phys. Rev. Lett. ,167604 (2012).[15] K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H.Castro Neto, Science , aac9439 (2016).[16] S. E. Shirsath, C. Cazorla, T. Lu, L. Zhang, Y. Y. Tay,X. Lou, Y. Liu, S. Li, and D. Wang, Nano Lett. , 1262(2019).[17] M. Reis and S. Soriano, Appl. Phys. Lett. , 112903(2013).[18] D. A. Rehn, Y. Li, and E. J. Reed, Phys. Rev. Mater. , 114004 (2018).[19] M. Aoki and H. Amawashi, Solid State Commun. ,123 (2007).[20] A. Avetisyan, B. Partoens, and F. Peeters, Phys. Rev.B , 195401 (2009).[21] M. Koshino and E. McCann, Phys. Rev. B , 125443(2009).[22] A. Avetisyan, B. Partoens, and F. Peeters, Phys. Rev.B , 115432 (2010).[23] M. Craciun, S. Russo, M. Yamamoto, J. B. Oostinga,A. Morpurgo, and S. Tarucha, Nat. Nanotechnol. , 383(2009).[24] C. H. Lui, Z. Li, K. F. Mak, E. Cappelluti, and T. F.Heinz, Nat. Phys. , 944 (2011).[25] M. Yankowitz, F. Wang, C. N. Lau, and B. J. LeRoy, Phys. Rev. B , 165102 (2013).[26] K. Zou, F. Zhang, C. Clapp, A. MacDonald, and J. Zhu,Nano Lett. , 369 (2013).[27] Y.-P. Wang, X.-G. Li, J. N. Fry, and H.-P. Cheng, Phys.Rev. B , 165428 (2016).[28] A. Y. Kuntsevich, Y. Tupikov, V. Pudalov, and I. Bur-mistrov, Nat. Commun. , 1 (2015).[29] V. Y. Tsaran, A. Kavokin, S. Sharapov, A. Varlamov,and V. Gusynin, Sci. Rep. , 1 (2017).[30] D. Grassano, O. Pulci, V. Shubnyi, S. Sharapov,V. Gusynin, A. Kavokin, and A. Varlamov, Phys. Rev.B , 205442 (2018).[31] V. Shubnyi, V. Gusynin, S. Sharapov, and A. Varlamov,Low Temp. Phys. , 561 (2018).[32] Y. Galperin, D. Grassano, V. Gusynin, A. Kavokin,O. Pulci, S. Sharapov, V. Shubnyi, and A. Varlamov,J. Exp. Theor. Phys. , 958 (2018).[33] I. Sukhenko, S. Sharapov, and V. Gusynin, Low Temp.Phys. , 264 (2020).[34] B. Partoens and F. Peeters, Phys. Rev. B , 075404(2006).[35] B. Partoens and F. Peeters, Phys. Rev. B , 193402(2007).[36] M. Koshino and E. McCann, Phys. Rev. B , 165409(2009).[37] S. Yuan, H. De Raedt, and M. I. Katsnelson, Phys. Rev.B , 235409 (2010).[38] Y. Mohammadi, R. Moradian, and F. S. Tabar, SolidState Commun. , 1 (2014).[39] C. Bao, W. Yao, E. Wang, C. Chen, J. Avila, M. C.Asensio, and S. Zhou, Nano Lett. , 1564 (2017).[40] M. S. Dresselhaus and G. Dresselhaus, Adv. Phys. , 1(2002).[41] L. Hao and T. Lee, Phys. Rev. B , 165445 (2010).[42] Z. Kutnjak, B. Roˇziˇc, and R. Pirc, Wiley Encyclopediaof Electrical and Electronics Engineering , 1 (1999).[43] C. Goupil, W. Seifert, K. Zabrocki, E. M¨uller, and G. J.Snyder, Entropy , 1481 (2011).[44] B.-R. Wu, Appl. Phys. Lett. , 263107 (2011).[45] I. Redouani, A. Jellal, A. Bahaoui, and H. Bahlouli,Superlattice Microstruct. , 44 (2018).[46] W. Bao, L. Jing, J. Velasco, Y. Lee, G. Liu, D. Tran,B. Standley, M. Aykol, S. Cronin, D. Smirnov, et al. ,Nat. Phys. , 948 (2011).[47] F. Guinea, A. C. Neto, and N. Peres, Phys. Rev. B ,245426 (2006).[48] S. Latil and L. Henrard, Phys. Rev. Lett. , 036803(2006).[49] Y. Nam, D.-K. Ki, D. Soler-Delgado, and A. F. Mor-purgo, Science362