A Tagliacozzo

G Campagnano

D Giuliano

P Lucignano

Benoit Jouault

Thermal transport driven by charge imbalance in graphene in magnetic field, close to the charge neutrality point at low temperature: Non local resistance

Thermal transport driven by charge imbalance in graphene in a magnetic field close to the charge neutrality point at low temperature: Nonlocal resistance

I. INTRODUCTION

Taming quantum transport at the edges of high mobility graphene Hall bars provides control of the quantization of the fractional Hall effect and is the prerequisite for the implementation of non-Abelian braiding statistics of excitations, which has been proposed as a tool for alternative quantum information processing 1 . In a two-dimensional electron gas, Quantum Hall Effect is a nowadays paradigmatic example of a "bulk" incompressible insulating phase, with compressible chiral edge states, carrying a charge either integer (Integer Quantum Hall Effect) 2 , or fractional (Fractional Quantum Hall Effect) 3 . In addition to chiral current-carrying states, at so-called "non-Laughlin fillings", such as ν = 2/3, additional counterpropagating neutral boundary states have been predicted [4][5][6] , which can in principle be detected by e.g looking at the extra charge noise they generate, once excited 7,8 , by means of momentum-resolved tunneling 9 , or measuring current-current correlations in a pertinent generalization of the quantum point contact scatterer between edge states proposed in 10 . Evidence for the neutral edge modes has recently been provided in a simultaneous measurement of both the chemical potential and temperature at a quantum Hall edge heated by means of a quantum point contact. As a result, it has been found that, while the charge is exclusively transported downstream, when the edge is expected to have additional structures such as the neutral branch of counterpropagating modes at ν = 2/3, heat can be transported upstream. In addition, an unexpected bulk contribution to heat transport was also found at Integer Quantum Hall fillings, in particular at ν = 1 11 . In graphene Hall bars in an external magnetic field, the observation of large nonlocal resistances near the Dirac point adds fuel to the mystery about the role of the external confining potential that defines the edge.

Clean graphene at the Charge Neutrality Point ( CNP) attracts a lot of attention because it can provide an example of a strongly-interacting quasi-relativistic electronhole plasma, known as a Dirac fluid 12 . Coulomb interaction is controlled by the dimensionless parameter α f and the only relevant energy scale is temperature T which determines the inelastic scattering rate τ c ,

τ -1 c ∼ α 2 f k B T with α f = e 2 r v F ∼ 2 r ∼ 1. (I.1)
Eq.(I.1) is a hallmark of many quantum-critical systems [12][13][14] . While in a Fermi liquid the heat current is strictly related to the mass current, close to the CNP the thermal conductivity is expected to be enhanced, because particles and holes move both in the direction of the thermal gradient with a strongly reduced relaxation rate, due to the linearity of the energy dispersion relation 15 .

Backscattering of Dirac fermions is also suppressed if intervalley scattering, which comes into play only in presence of strong disorder, is neglected 16 . Violations of the Wiedemann Franz law [START_REF] Ashcroft | Solid State Physics[END_REF] have been reported, as well, as a consequence of the strong Coulomb interactions between thermally excited charge carriers at the CNP [START_REF] Crossno | [END_REF] . Noticeably, a strong enhancement of the Lorentz number by a factor of about 2.5 its Wiedemann-Franz value has also been derived in a ballistic bilayer graphene in the presence of trigonal warping term and with the electrochemical potential close to the Lifshitz energy, as a consequence of the van Hove singularities in the single particle density of states 19 . In a Dirac fluid, the violations of the Wiedemann Franz law have been reported in a temperature window k B τimp < T < k B τ ph , where 1/τ ph refers to the relaxation mechanism introduced by inelastic phonon scattering at higher temperatures. The reason for this is that charged impurities in the substrate generate a local finite carrier concentration (particle or hole puddles 20 ), which can give rise to a position dependent chemical potential µ(r), so that both electric and thermal transport may be dominated by an elastic scattering rate,

τ -1 imp ∼ Ze 2 r 2 ρ imp max[µ, k B T ] .
τ -1 imp is naturally proportional to the impurity density and is responsible for restoring a Fermi liquid-like behavior. At larger doping, when the chemical potential exceeds k B T , the inelastic-scattering rate tends to the familiar Fermi-liquid form ∼ T 2 /µ if the interactions are screened. Favourable experimental bath temperatures to monitor the Dirac fluid properties in boron nitride (hBN) encapsulated graphene has been found to be [START_REF] Crossno | [END_REF] T > 40 o K.

Regarding charge transport, the longitudinal conductance is known to be finite at the CNP, in the absence of any disorder: σ Q = 4e 2 /πh 21,22 . Indeed, the Zitterbewegung, with the creation of virtual zero total momentum electron hole pairs, appears to be responsible for the evanescent zero energy Dirac modes, which provide the finite conductivity. 23,[START_REF] Katsnelson | Graphene: Carbon in Two Dimensions[END_REF] In this paper, we are mainly interested in giant nonlocal voltages appearing in graphene in the QHE regime close to the CNP at low temperature and in the presence of a magnetic field. These nonlocal voltages have been observed by various authors [START_REF] Abanin | [END_REF][26][27][28] , including some of us. 29 , though under operating conditions rather far from the two limits (Dirac and Fermi liquid behavior) mentioned above. However, here we claim that one cannot explain the unusual electric properties found close to the CNP without accounting for thermal diffusion, even at such low temperatures.

For the sake of clarity, we now present some typical experimental results, complementary of those shown in Ref. [29]. Figure 1a shows the longitudinal and transverse magnetoresistances observed in a graphene Hall bar epitaxially grown onto SiC. The magnetic field is orthogonal to the sample plane, the temperature is T = 1.7 K and the Fermi energy is tuned close the CNP (Hall concentration p 10 10 cm -2 , mobility 10, 000 cm 2 /Vs). The graphene Hall bar is macroscopic and has a length 400 µm and a width w 100 µm. The peak in the longitudinal resistance R xx at B = 0 is due to weak or strong localization and its discussion is out of the scope of this paper. When the magnetic field is turned on, this peak disappears, while another peak can be resolved around B = ±2.5 T. At approximately the same value of B, the transverse resistance has an ambipolar behavior and goes to zero. We conclude that for some reasons (intrinsic to graphene on SiC), the Fermi energy smoothly increases with B and crosses the CNP around B = 2.5T. More re- markably, Figure 1b shows the nonlocal resistance measured across the bar at two different places, when the Hall bar is current biased at one of its extremities, along its width (transverse direction), while it is kept as an open circuit along its length (longitudinal direction). Here we define the resistance R ij,kl = V kl /I ij , where V kl is the voltage drop between contacts k and l, and I ij is the current biased between contacts i and j (see inset in Fig. 1b). Remarkably, the Ohmic resistance R Ohmic decreases with the distance from the applied current much more slowly than one would as one would naively expect from the classical spreading of the charge flow inside the Hall bar, which would yield

R Ohmic = 4 π ρ xx e -πD/w (I.2)
where ρ xx is the resistivity and D is the distance separating the current injection point from the voltage detection (see Fig. 1b). The measured nonlocal resistances observed in Fig. 1b are much larger than the ones predicted by Eq. I.2: at B = 3 T, Eq. I.2 predicts R 28,37 300 Ω and R 28, 46 10 Ω while experimentally one finds R 28,37 10 kΩ and R 28,46 2 kΩ. Therefore, some other explanation must be found.

If the channel length matches the spin diffusion length in graphene, the nonlocal configuration makes it possible to detect spin-related signals in the spin Hall regime. However, the distance between probes 2 and 4 for the Hall bar presented in Fig. 1 is 200 µm, which is much larger than the usual spin relaxation lengths reported in graphene. Nevertheless, the nonlocality was first attributed to the existence of a Zeeman spin Hall effect (ZSHE), where the combination of both Lorentz force and Zeeman spin splitting gives rise to a spin imbalance which propagates along the Hall bar. 30 Later on, additional experimental works demonstrated that a large part of the nonlocal voltage is insensitive to the magnetic field parallel to the graphene plane -suggesting the predominance of orbital effects over ZSHE.

Thermal effects can also give rise to nonlocal voltages. At the excitation point, heat flow is induced by Ettingshausen and Joule effects, perpendicular to the charge flow. This heat flow diffuses and produces a nonlocal voltage far from the charge current, via Nernst effect. 27,31 However, the observed nonlocal resistances yield Nernst coefficients which seem to be unrealistically large (700 µV/K in Ref. [27], more than 10 mV/K in Ref. [29]).

In two dimensional electron gases, nonlocal resistances in the QHE regime were first observed by McEuen et al. 32 and unambiguously attributed to a coupling of the edge and bulk conducting pathways. Later, the same model was used to explain the appearance of nonlocal voltages in other 2D topological insulators 33 . This model was also extended to graphene near the CNP, to explain dissipative QHE [START_REF] Abanin | [END_REF] and nonlocal voltages 28,29 . In Ref. [25], a plateau for the Hall conduction at filling factor ν = 0 has been found, accompanied by a peak ρ xx > 0 centered at ν = 0. It was argued that this peak is originated by Zeeman spin splitting of the former ν = 0 Landau level centered at fillings ν = ±1/2. The resistance ρ xx > 0 accompanying the unconventional ν = 0 plateaus could be the result of interference of edge states located at opposite boundaries via the bulk sandwiched in between, having diffusive conductance. A simple semiclassical description of the interference provides a convincing model (henceforth named "Abanin et al. model") 34 .

The same approach can be used to explain nonlocality. Indeed, the peak in R xx around B = ±2.5 T in Fig. 1b corresponds to the diffusive longitudinal peak observed at the CNP in Ref. [25]. In Ref. [29], a numerical model, again based on the coupling of both edge and bulk pathways, was used. Quantitatively, an appropriate choice of the conductivities of the carriers (particles and holes close to the CNP) could explain most of the non-local resistances.

However, remarkably, the nonlocal resistances are also asymmetric under flipping B 28 , and the closer the measuring contacts are to the applied bias, the larger is the asymmetry. This can be readily observed in Fig. 1b. The model for the NLR in terms of dissipative scattering of the edge states 29 accounts for the magnitude of the effect but does not catch the origin of the asymmetry when the B field is reversed. In Ref. [28], an explanation for the asymmetry has been proposed, assuming that grain boundaries in the graphene sheet have different transmission either for spin up or down. This explanation is valid for the experiments of Ref. [28], which were performed on polycristalline graphene. However, the +B/ -B asymmetry is also clearly visible in Fig. 1b and in Ref. [29], where graphene has been obtained by epitaxy on SiC and does not contain grain boundary. Also, edge dislocations could play a role, as well 35 .

The purpose of our work is to propose another and more generic approach to justify further contributions to the nonlocal voltages and their asymmetry. In line with the approach proposed by Ref. [34], we consider transport in the interior (henceforth named "bulk") of the graphene Hall bar away from the edges by reconsidering thermal effects. Close to particle-hole symmetry point, both particles (e -) and holes (h + ) contribute to transport. Here we explore the possibility that the non-equilibrium conditions induced by the electric field E applied at the one end of the Hall bar and Lorenz force could give rise to local charge imbalance between particles and holes which could have long relaxation times 15,36 at low temperature and could diffuse under the action of Joule heating and the thermal gradient generated by the injected current (Ettingshausen effect), far away from the source. The contribution to the non-local voltage produced in this way depends on the orientation of the magnetic field because particles and holes exchange their role when the field is flipped. This could be the origin of the anisotropy when the B field is reversed. If this interpretation holds, the magnitude of the anisotropy is a measure of the relative weight between the contribution to dissipation in the transport at the edges and the thermal diffusion of the charge imbalance in the bulk.

The bulk carriers may be created by leakage from the edge states into the bulk as assumed in the Abanin et al. model. Charged impurities in the substrate which may have a preferred charge and weak disorder could also locally enhance the charge imbalance 37 . Or intrinsically, by adding next-nearest-neighbor hopping terms in the tightbinding calculation for graphene bands, the Fermi velocity at the Dirac cone can be different between particles and holes. Elastic scattering is ineffective in reducing the charge imbalance. On the other hand, thermal inelastic scattering may produce particle hole pairs with different relaxation times, τ c → τ ee , τ hh , τ eh , so that the charge imbalance does not relax. e-ph scattering, as a source of imbalance relaxation is ruled out at low temperatures.

Close to the quantum critical point ( µ, T ∼ 0), there is an emergent relativistic invariance of the interacting electronic Hamiltonian for the clean sample. In this regime, a relativistic hydrodynamical approach is expected to apply as long as ωτ c << 1 (where τ c is the inverse cyclotron frequency).

The paper is organised as follows:

In Section II we first review standard results about the linear response to external perturbations including electric and thermal gradient. Therefore, we extend those results to systems with two types of charge carriers, e -and h + . The Boltzmann transport theory, which applies well to the Fermi Liquid regime µ >> T , cannot be safely applied here, because µ ∼ T . Nervertheless, assuming the perturbation to be small, the relativistic Dirac electron picture and the Lorenz invariant relativistic equations of motion are still valid in the linear approximation, because the constraints posed by the conservation laws on the energy-momentum tensor including dissipative processes ( viscosity and thermal conduction) [START_REF] Landau | Fluid Mechanics[END_REF] . In the limit in which the disordered impurity distribution is smooth on the scale of the lattice constant, one can assume that the charge imbalance can be described by two different chemical potentials µ e , µ h , respectively for e -and h + . At equilibrium, one has µ e = -µ h = µ and, accordingly, in this case one can define the electro-chemical potential µ = (µ e -µ h )/2 and the imbalance chemical potential

µ I = (µ e + µ h )/2.
Section III reviews the model of Ref. [34] for the leakage of edge current into the bulk of the Hall bar. In particular, while the original version of the model focuses onto the linear response to an applied longitudinal field E, parallel to the edges (x direction), here we add a leakage of the edge states towards a resistive propagation away from the edges, to justify a nonzero ρ xx = 0 coexisting with the plateaux of σ xy at ν = 0. An expression for the Hall voltage V H , neglecting thermal effects is derived. Eventually, the model is indeed shown to reproduce some of the conductance features of both HgTe QW and graphene Hall bar. However, NLR across the Hall bar is probed in a different setup, by injecting current across the Hall bar (ŷ direction) at the origin of our x-axis and measuring the voltage at a point x = D > 0.

In Section IV the model of the previous Section is extended, to include the thermal current and a transverse applied electric field, together with the Peltier thermal gradient which induces diffusion of the imbalance carriers in the bulk along the x direction, with particles spatially separated from holes, away from the applied bias. The diffusive longitudinal carrier transport with velocity v F is connected to the charge imbalance via an additional constitutive equation, Eq.(II.12). Phenomenologically, it is possible to include the charge imbalance relaxation rate in this equation. The diffusion equations involving the temperature and the charge imbalance chemical potential on one side and the imbalance charge density and the non local voltage (NLV) on the other, are derived and approximately solved for the 'near region', close to the origin. Nernst effect fixes the boundary condition at the origin. Details of the derivation are given in Appendix A,B,D,E and F. The consistency of the equations is discussed, as well as their physical content. The Fermi liquid transport parameters are discussed in Appendix C.

In Section V an estimate of the magnitude of the Ettingshausen parameter, P E , is derived for a magnetic field ∼ 2 T esla. The thermal conduction and the Peltier and Nernst coefficients are chosen in the range of experimental values quoted in the graphene literature. Thermal transport necessarily involves also some thickness of the graphene Hall bar which cannot be a priori determined. The delicate point here is the estimate of the three dimensional (3d-) thermal energy per carrier eQ and the related 3d-effective carrier density n 3d . The bulk conductivities for particles and holes are extracted from the model of Ref. [29]. The parameter choice is discussed in Appendix C. The result for P E is about one or even two order of magnitude larger than the one for bismuth.

In Section VI the non local voltage is presented, which gives rise to the NLR. The NLV is plotted in Fig. (5), as a function of the magnetic field at various distances D from the origin and the physical picture that emerges from the plots is presented and discussed.

Section VII includes the Summary of the crucial points in the derivation of the results and the Conclusions.

II. LINEAR RESPONSE IN HYDRODYNAMICS

In this Section the hydrodynamical picture for matter and energy transport is reviewed, by also taking into account the additional features of the Dirac fluid. In particular, the total average charge density, expressed in terms of the average number densities of the particle and hole fluids is given by ρ = -e(n e -n h ) ( -e < 0 is the electron charge in the following). The charge current density is J = ρ u -e (ν e -ν h ), where u is the center of mass velocity of the total fluid and ν e,h are the fluctuations produced by the applied electric field E and thermal gradient (-∇T ). The particle current density is J n = n u + ν e + ν h , while the energy current density is given by

J Q = T s u Q + µ e ν e + µ h ν h = ( + P) u Q + µ I J n - µ e J.
(II.1)

In Eq.(II.1) u Q is the center of mass velocity of the J Q current, s is the entropy per unit volume and + P = T s + µ e n e + µ h n h is the enthalpy per unit volume. P = P -M B includes the thermodynamic pressure of the carrier fluid and the work done by B onto the currents induced by the magnetization M( B) at the boundary [START_REF] Cooper | [END_REF] . Within linear response theory, one obtains the relations:

J J Q = σ α T α κ E -∇T , (II.2)
with σ,α, α and κ being 2 × 2 matrices which depend on the coordinate, as well as on the flavour (e -or h + ) label.

To set up the notation, we start by considering a onecomponent fluid. By definition, we get for the con-ductivity matrix σ = σ xx 1 + σ xy ˆ , where ˆ is the twodimensional antisymmetric tensor: xx = yy = 0, xy = yx = 1. α, α and κ are the thermoelectric conductivities which determine the Peltier, Seebeck, and Nernst effects. We will assume that, on the scale of the measured samples, the response kernels can be taken in the uniform and static limit (by also assuming that, in performing the calculations, the q → 0 limit is taken first, followed by the ω → 0 limit). Due to Onsager reciprocity, there is no difference between α and α as long as currents are uniform ( q = 0 limit ).

The thermal conductivity, κ, is defined as the heat current response to -∇T , in the absence of electric current (J = 0, i.e. electrically isolated boundaries) and is given by

κ = κ -T ασ -1 α (II.3)
Instead, κ applies to samples connected to conducting leads, allowing for a stationary current flow.

In our case, the relativistic stress energy tensor provides the energy flux T 0 α from the four-momentum conservation. We write the energy current for a Fermi liquid as follows, to match with the expected form from Eq.(II.2) 40 :

J Q = + P n J - σ e 2 -T ∇ µ T + e c v × B → + P n J - + P ne 2 σ 1 T (-∇T ) + + P n e σ 0 σ B σ E - ( + P) n 2 σ e 2 ∇P, (II.4) σ = σ 0 1 + (ω c τ ) 2 1 -ω c τ ω c τ 1 , σ-1 = σ -1 0 1 ω c τ -ω c τ 1 , σ 0 ∼ σ 3d = ρ v F τ k F . (II.5)
σ 3d is an appropriate scalar reference conductivity of the 3d system assumed to be uniform and ω c = eB/( ck F /v F ) is the cyclotron frequency for the linear energy dispersion. (We have interpreted ρv y as corresponding to σ yx E x for the x component of the energy current e v × B, where σ yx = σ B ≡ ρ c/B, with reference to a Drude metal in the Hall configuration.)

The first term describes a convective flux which corresponds to the last term in Eq.(II.1). In the case of a two-component fluid the first term at the right hand side of Eq.(II.1), referring to the center of mass convective motion of the two components, should also appear. The energy flux also includes a term depending on the derivatives of P that is absent in the non relativistic result [START_REF] Landau | Fluid Mechanics[END_REF] . However, we recognize here that the x component ∂ x P ∝ ∇ × M × B is a term arising from the energy contribution of the magnetization currents ∇ × M flowing at the boundaries of the sample [START_REF] Cooper | [END_REF] . In the following, we will take care of the edge-bulk interaction in the Hall bar in a different way, which eventually enables us to get rid of this term. = κ xx and -+P n e 3d σ B ˆ = T α. In the case of a one-component plasma for a Drude metal, n ≡ n 3d and ρ are the carrier and charge density, respectively. The transverse dc response yields the standard Hall effect, with σ xy = ρc/B, and the transverse Peltier effect, α xy = ( + P) 3d c/(T B) = s c/B, which can be interpreted as charge ρ and entropy s density [START_REF] Cooper | [END_REF] , drifting with the velocity v D = c E × B/B 2 . Hence, the Nernst coefficient is: which diverges for n 3d → 0, but goes as 1/B at large B (ω c τ >> 1) 41 . When T >> µ, by posing (( + P)/n) 3d = k B T in the definition of κ xx , we get a rewriting of the Wiedemann-Franz law, to be compared with the Fermi liquid result:

Consistency

e N = E y -∂ x T = -σ -1 α xy = σ -1 0 s c/B -s c/B 0 xy = 1 σ 0 ( + P) 3d c T B = + P n 3d 1 T 1 ω c τ , (II.6)
κ xx σT = π 2 3 k B e 2 .
(II.7)

In a clean relativistic system at B = 0, Lorenz invariance implies σ xx (ω → 0) = α xx (ω → 0) = 0 for a one-component plasma 14 . In a reference frame moving at the constant velocity v D with respect to the laboratory frame, the observed electric field vanishes and, hence, in that frame the charge currents vanish. As v D ⊥ E, transforming back to the laboratory frame, the longitudinal field is still vanishing. These results hold beyond the hydrodynamic description even when ω c τ ee >> 1, as long as Lorentz invariance holds.

In graphene, electric conductivity reaches the minimum value 21,42 σ Q = 4e 2 /πh at the CNP and κ xx is finite. Using Eq.(II.3), we get 14 :

κ F L xx (ω = 0, µ) = + P n e 2 3d σ 3d T 1 - σ B σ 3d 2 .
(II.8)

Similarly,

α F L xy = ( + P) 3d T B c.
(II.9)

Moving now to the two component fluid for a uniform isotropic system, with average 3d carrier number density n 3d = n e + n h , a derivation similar to the one given in Eq.(II.4) provides contributions to the fluctuations ν e,h appearing in the charge and particle current densities J and J n , in terms of the gradient of the temperature, of µ, and of µ I , by means of the conductivities σ ee , σ hh and σ eh . At zero applied magnetic field, the conductivities are constructed from the Drude relaxation times τ ee , τ hh and τ eh as 15 :

e (J -ρu) = -e 2 (ν e -ν h ) = [σ ee + σ hh -2σ eh ] eE -T ∇ µ T + [σ ee -σ hh ] -T ∇ µ I T -e 2 (J n -nu) = -e 2 (ν e + ν h ) = (σ ee -σ hh ) eE -T ∇ µ T + [σ ee + σ hh + 2σ eh ] -T ∇ µ I T . (II.10)
While σee and σhh in Eqs.(II.10) are easily interpreted as particle and hole conductivities, respectively (in the next Sections they will be generically referred to as σ ab,i , with a, b = x, y and i = 1, 2, respectively for e and h), σeh = σhe refers to a drag conductivity between particles and holes. σ ee + σ hh -2σ eh plays the role of σ 3d 15 . We consider the limit in which both x-component terms of J are negligibly small (that is E x ≈ ∂ x µ ≈ 0 and to linear order in σee -σhh , as µ I ∝ σee -σhh , as well) and the only contribution relevant to our derivation is the x-component of the second term at the r.h.s. of the particle current J n , of dimension [energy/t 2 ]. As this is the most important contribution appearing in our model, in the following we extensively discuss about it.

While Hall transport is essentially a 2d phenomenon, particle current and thermal transport are essentially 3d. Accordingly in the following we and we will have to carefully account for that difference. Also, defining the conductivities requires here some care. Here, σ i ≡ σ xx,i and σ xy,i (i = 1, 2) are 3d conductivities and have dimension [1/t]. Also, an aspect ratio must be be introduced to take into account the effective thickness d ≈ 0.3 Å of the graphene sheet grown on top of the SiC surface in the Hall bar, so that the σ ab,i conductivities are related to the longitudinal resistance R according to:

1 σ 3d = R d w . (II.11)
This implies that a volume unity V ∝ [ wd] can be introduced, which will appear in the rest of the paper.

As mentioned in the Introduction, the particle and energy currents along the x direction are related to the imbalance between electrons and holes induced by the applied electric field E ŷ and to the corresponding electric current J y . We are mostly interested in the diffusion of the x-component of J n along the x direction.

Here we introduce δJ n ≡ J n,A -J n,B , where A,B denote areas close to each of the two edges, which is not a locally conserved current. In a steady state, the diffusion ∂ x δJ n of the particle current density is induced by the particle/hole imbalance charge ρ I , which we will appropriately define in Section IV, Eq.(IV.6), in terms of the chemical potential imbalance. Accordingly, we write:

∂ ∂x (J n • x) = -v F 1 n V ρ I e .
(II.12)

Here n V plays the role of a phenomenological effective relaxation length for the thermal diffusion λ Q , which is ∼ µm. (Phonon scattering is expected to play no role as it is virtually frozen out at low temperature.) On the basis of the previous remarks, we are eventually able to build up a possible interpretations of the non local resistance of Fig. 1, by assume an open circuit in the x direction parallel to the Hall bar edges. On injecting a current across the Hall bar at x = 0, a nonequilibrium charge imbalance, as well as a temperature gradient -∂ x T (x), are generated in the x direction, due to the Ettingshausen effect. Enhanced thermal conduction κ xx (x) along the Hall bar drives the particle current, with space separated particle and holes close to edges A and B so that the imbalance reaches the point x = x b where a transverse voltage v y is measured. Charge conduction in the edge channels is assumed not to be influenced by thermal processes, but they provide some dissipation and a leakage into "bulk" states of carriers of opposite charge at opposite edge states.

III. THE ABANIN ET AL. MODEL

In this Section we recall the main features of the model of Ref. [34]. While the model includes the longitudinal dissipation coexisting with the leakage of carriers from the edges into the bulk of the Hall bar, it does not include thermal effects. However, its careful discussion is an important preliminary step for us, in order to extend it by including thermal effects, as well, which will be the subject of Section IV.

Under the action of an orthogonal magnetic field B, close to the CNP (µ ≈ 0), particle and holes states counterpropagate at the edges of the graphene Hall bar. Let x be the direction along the edges and ŷ the direction orthogonal to the edges and w be the width of the Hall bar. We denote by A the upper edge at y = 0 and by B the lower edge at y = w. In the following, labels A,B refer to the region close to the edge A and B, respectively, while the interior of the Hall bar will be denoted as the "bulk" (see Fig. (3)). Charge transport along the edges is described by the electrochemical potentials c 1,2 for particles and holes, respectively. (Here, we do not account for the spin, which merely provides a factor of 2 in the final results). As the Hall bar is not fully insulating in the bulk, we assume nonzero isotropic bulk conductivities σ xxi = σ yyi = σ i (the label i = 1, 2 is for e and h respectively) so that, if ψ 1,2 (y) are the electrochemical potentials for the bulk, σ yyi ŷ • ∇ψ i (y) is the bulk current orthogonal to the edges. The terms describing leakage of current from the edges into the bulk are respectively given by g (ψ i -c i ) A and g (c i -ψ i ) B . At edge A the linear response to an electric field E in the x direction is

E = γ (c 1 -c 2 ) A + (-) i g (ψ i (y = 0) -c i ) A , while, at edge B, E = -γ (c 1 -c 2 ) B -(-) i g (ψ i (y = w) -c i ) B .
The electrochemical potentials c 1,2 and ψ 1,2 (y) are complementary to the density fluctuations ν e,h of the previous Section. In the absence of thermal effects, we assume mirror symmetry w.r.to the longitudinal axis lying halfway between edges A and B, that is, we set ψ i (0) = -ψ i (w) = s i . The derivative of the bulk electrochemical potential in the ŷ-direction is linearized as a finite difference ŷ

• ∇ψ i ∼ (ψ i (0) -ψ i (w))/w = 2 s i /w.
The electric field E is in the x-direction and appears in the four equations for the edges (1e, 2e, 3e, 4e):

1e : -E = -γ (c 1 -c 2 ) A + g (s 1 -c 1 ) A 2e : E = γ (c 1 -c 2 ) A + g (s 2 -c 2 ) A 3e : E = -γ (c 1 -c 2 ) B -g (s 1 + c 1 ) B 4e : -E = γ (c 1 -c 2 ) B -g (s 2 + c 2 ) B (III.1)
The bulk currents explicitly appear in the four equations for the bulk (1b, 2b, 3b, 4b):

1b : 0 = -2 σ 1 w s 1 + σ xy,1 E + g(c 1 -s 1 ) A 2b : 0 = 2 σ 2 w s 2 -σ xy,2 E -g(c 2 -s 2 ) A 3b : 0 = 2 σ 1 w s 1 -σ xy,1 E + g(s 1 + c 1 ) B 4b : 0 = -2 σ 2 w s 2 + σ xy,2 E -g(s 2 + c 2 ) B (III.2)
The conductivities γ , γ and g , g differ just by unities, as the former (primed) ones have dimension [1/L] while the latter ones have dimension [(L t) -1 ]. Accordingly, the chemical potentials have dimension [e/L] and the the electric field E has dimension [e/L 2 ]. Note that, the sum of Eq.s(III.1:1e,2e) minus the sum of Eq.s(III.1:3e,4e) yields:

1 2 [(c 1 + c 2 ) A -(c 1 + c 2 ) B ] = (s 1 + s 2 ), (III.3)
where, at the right hand side, we write [(s 1 + s 2 ) A + (s 1 + s 2 ) B ] /2 = (s 1 + s 2 ) ≡ s + . This is consistent with isotropy in the bulk. This approximation holds approximately also when thermal effects are taken into account, as long as the relaxation length for the carrier imbalance λ Q is large enough with respect to the width w and the length of the Hall bar. The left hand side is the definition of the Hall voltage V H . In the absence of thermal effects the Hall voltage can also be defined from the bulk potentials as V H = s + . The corresponding electrochemical potentials µ, being charge dependent, involve the differences instead of the sums:

c -,A = (c 1 -c 2 ) A , c -,B = (c 1 -c 2 ) B , s -= (s 1 -s 2 ).
The equality in Eq.(III.3) states the absence of the carrier imbalance between edges and bulk, which is a typical situation in the absence of thermal effects. Indeed, in the presence of thermal effects (which we take into account in the next Section by assuming that they play a role in the 3d-like bulk but not at the edges), it is indeed violated and, as it will clearly appear in the following, the non-equilibrium difference

(c 1 + c 2 ) A -(c 1 + c 2 ) B -2(s 1 + s 2 ) ≡ c +A -c +B -2 s +
provides the imbalance chemical potential µ I close to the edges A and B.

We now choose as independent variables: c -,A , c -,B , s 1 , s 2 and drop the apex in the definitions of the conductivities, with an appropriate choice of the units for E. From Eq.s (III.2) we get ( σ xy,i ≡ η i ):

   -(2γ + g) 0 g -g 0 -(2γ + g) -g g γ 0 0 -2 σ2 w 0 -γ 2 σ1 w 0       c -,A c -,B s 1 s 2    =    -2 E 2 E -(η 2 -1) E (η 1 + 1) E    (III.4)
The solution is:

(c 1 -c 2 ) A = -(c 1 -c 2 ) B = 4 + wg 1 σ 1 + 1 σ 2 -wg η 2 σ 2 - η 1 σ 1 • E/D D = 2(2γ + g) + gwγ 1 σ 1 + 1 σ 2 s 1 = [γg(η 1 + η 2 ) + 2(2γ + g) η 1 σ 2 + 2gσ 2 ] • E/(2σ 1 σ 2 D) s 2 = [γg(η 1 + η 2 ) + 2(2γ + g) η 2 σ 1 -2gσ 1 ] • E/(2σ 1 σ 2 D) (III.5)
The Hall Voltage of Eq.(III.3) is readily obtained as:

V H = -2w(2γ + g) [η 1 (σ 2 + λ) + η2 (σ 1 + λ)] • E/(2σ 1 σ 2 D), (III.6)
where λ = gγw/(2γ + g), η1 = η 1 + g/(2γ + g) and η2 = η 2 -g/(2γ + g). The ratio g/(2γ + g), plays the role of a contribution to the Hall conductivity determined by current leaking from the edges into the bulk. Eq.s(III.6) were originally derived in Ref. [34]. For completeness, here we report the plots of the longitudinal and transverse resistivity ρ xx , ρ xy , together with the corresponding conductances G xx and G xy as functions of the filling ν (Fig.( 4)) 34 . A kind of plateau at ν ≈ 0 is recognizable in G xy accompanied by a large peak in ρ xx .

IV. THERMAL RELAXATION OF CHARGE IMBALANCE

Here we extend the model of Section III to include energy transport and thermal effects on the conductance.

To keep in touch with the experimental setup we discuss, we consider the Hall bar of length as an open circuit in the longitudinal (x) direction. A current is injected in the ŷ direction, orthogonal to the edges A and B, by applying an electric field E at the contacts at x = 0. We assume that no thermal effects involve the edges. Therefore Eq.s(III.1,1e-4e) for the edge propagation in the x-direction do not change, except for the fact that there is no driving electric field at open circuit. However, the electrochemical potentials are expected to depend on the x coordinate along the edges and their derivative replaces the external electric field that was applied in the model of Section III, when the circuit was closed. Therefore, the equations for the edges are now given by : 1e :

∂ x c 1A = -γ (c 1 -c 2 ) A + g (s 1 -c 1 ) A 2e : -∂ x c 2A = γ (c 1 -c 2 ) A + g (s 2 -c 2 ) A 3e : -∂ x c 1B = -γ (c 1 -c 2 ) B -g (s 1 + c 1 ) B 4e : ∂ x c 2B = γ (c 1 -c 2 ) B -g (s 2 + c 2 ) B (IV.1)
We now look for the components of J and of J Q in Eq.s (II.1) respectively along the ŷ direction and along the x direction. Eventually, on inserting them in Eq.(II.2), we trade it for a system of differential equations for the temperatures and the chemical potentials only.

Close to the edge A, or B, the components of J in the ŷ direction, coming from opposite edges, can be recovered by respectively summing Eq.s(III.2: 1b,2b ) and Eq.s(III.2: 3b,4b ). The terms σ yx,i E x appearing in Eq.s(III.2: 1b-4b) vanish. The portion of charge current lost by the edge potential (c 1 -c 2 ) A is, to first order in g, the fraction of current emerging from the imbalance close to the edge, g 2γ+g t -1 0 ∂ x (c 1 +c 2 ) A , as derived from the difference between Eq.s(IV.1, 1e,2e). t 0 is an appropriate time scale which we will not have to specify in the following. The terms ∝ g should not be counted twice. The corresponding manipulations apply in the region close to edge B. We therefore set:

1b : J A • ŷ -ρ u = -(ν e -ν h ) A = -2 σ1 w s 1 -σ2 w s 2 A + g(c 1 -c 2 ) A -g (s 1 -s 2 ) A + g 2γ+g t -1 0 ∂ x (c 1 + c 2 ) A 3b : J B • ŷ + ρ u = (ν e -ν h ) B = 2 σ1 w s 1 -σ2 w s 2 B -g(c 1 -c 2 ) B -g (s 1 -s 2 ) B + g 2γ+g t -1 0 ∂ x (c 1 + c 2 ) B . (IV.2)
Due to the open circuit condition along x, u is in the ŷ direction, as a consequence of the applied electric field E ŷ.

On the contrary, the components of the current J n , which give the fluctuations in the carrier transport along the x direction, can be induced from the difference between Eq.s(III.2: 1b,2b) for area A, or Eq.s(III.2, 3b,4b)for area B. The extra term to be considered coming from Eq.s(IV.1:1e,2e) for edge A is given by t -1 0 ∂ x (c 1c 2 ) A . The same happens for region B. Therefore, we get:

2b : (J n -n u) A • x = -(ν e + ν h ) A = -2 σ1 w s 1 + σ2 w s 2 A + g(c 1 + c 2 ) A -g(s 1 + s 2 ) A + t -1 0 ∂ x (c 1 -c 2 ) A 4b : (J n -n u) B • x = (ν e + ν h ) B = 2 σ1 w s 1 + σ2 w s 2 B -g(c 1 + c 2 ) B -g(s 1 + s 2 ) B + t -1 0 ∂ x (c 1 -c 2 ) B (IV.3)
(we have reported the labels 1b -4b on the very left, in correspondence with the ones of Eq.s(III.2)).

The gradients of the electrochemical and imbalance chemical potentials, ∇µ, ∇µ I appear here explicitly, according to the identification

∂ y µ = ∂ y (µ e -µ h )/2 ∼ [(c 1 -s 1 ) A -(c 2 -s 2 ) A ]/w and ∂ y µ I = ∂ y (µ e + µ h )/2 ∼ [(c 1 + c 2 ) A -(s 1 + s 2 ) A ]
/w, so that Eq.(IV.2,IV.3) are a special case of Eq.s(II.10), adapted to our case. Although Eq.s(IV.3) refers to the fluctuating part of the particle current J n only, they are sufficient to complete the characterization of the thermal effects from Eq.(II.2). This is because, according to Eq.(II.4), just linear terms in the gradients should be added to ( + P) u Q + ( + P) J to obtain J Q (see also Eq.s(II.10)). According to Eq.(II.1), J Q -( +P)u Q -( +P)J ≈ µ I J n , if we approximate µ/e with ( + P)/n. This implies that it is enough to know the contributions -e [ν e -ν h ] and -[ν e + ν h ] to the charge current density and to the number current density, respectively, to recover informations about the fluctuating components, which are proportional to the gradients of the temperature and of the chemical potential.

A. Response to thermal and potential gradients Following our derivation above, we now derive the equations describing the bulk charge and particle cur-rents arising in response to applied electric field and thermal gradients, E and ∇T , via the electrical and thermal conductivities corresponding to Eq.s(II.2), by means of Eq.s(IV.2, IV.3). In our geometric, the electric field is orthogonal to the edges, while the thermal gradient is parallel. In the reference frame in which u Q = 0, we get the following Eq.s(IV.4: 1b',3b') for the charge current oriented toward the bulk and orthogonal to the edges, and the following Eq.s(IV.4: 2b',4b') for the energy current in the x direction, produced by particle-hole imbalance: 1b :

-σE

+ α yx (-∂ x T ) A = -ρ u -2 σ 1 w s 1 - σ 2 w s 2 A + g(c 1 -c 2 ) A -g (s 1 -s 2 ) A + g 2γ + g t -1 0 ∂ x (c 1 + c 2 ) A 2b : T A α xy E + κ xx (-∂ x T ) A = + P n e 3d -2 σ 1 w s 1 + σ 2 w s 2 A + g(c 1 + c 2 ) A -g(s 1 + s 2 ) + t -1 0 ∂ x (c 1 -c 2 ) B 3b : σE -α yx (-∂ x T ) B = ρ u + 2 σ 1 w s 1 - σ 2 w s 2 B -g(c 1 -c 2 ) B -g (s 1 -s 2 ) B + g 2γ + g t -1 0 ∂ x (c 1 + c 2 ) B 4b : -T B α xy E -κ xx (-∂ x T ) B = + P n e 3d 2 σ 1 w s 1 + σ 2 w s 2 B -g(c 1 + c 2 ) B -g(s 1 + s 2 ) + t -1 0 ∂ x (c 1 -c 2 ) B .
(IV.4) (Note that, from now on, we refer all the physical quantities to a 3-d system so that, e.g., the conductivities have dimension sec -1 , according to Eq.(II.11)). When comparing Eq.s(IV.4) with Eq.s(III.2), we see that now the terms at the left hand side of Eq.s(IV.4) are no longer equal to 0, due to the electric field E = 0 (in the ŷ direction) and to the coupling to the thermal gradient, according to Eq.(II.2). Also, note that, in Eq.s(IV.4: 2b',4b'), we have introduced the prefactor +P n e 3d , to account for the proportionality in Eq.(II.4) and in Eq.(IV.3).

Within linear approximation, we assume that Eq.s(IV.1:1e-4e) still hold to lowest order and we use them to substitute the derivative term ∂ x (c 1 + c 2 ) A and ∂ x (c 1 -c 2 ) A from Eq.s(IV.1:1e, 2e) into Eq.s(IV.4: 1b') and Eq.s(IV.4: 2b') respectively, and similar derivative terms from Eq.s(IV.1:3e, 4e) into Eq.s(IV.4: 3b') and Eq.s(IV.4: 4b') respectively, to get:

1b : -σ E + α yx (-∂ x T ) A = -ρ u -2 σ1 w s 1 -σ2 w s 2 A -g 2γ 2γ+g s -,A 2b : T A α xy E + κ xx (-∂ x T ) A = -2 +P ne 3d σ1 w s 1 + σ2 w s 2 A 3b : σ E -α yx (-∂ x T ) B = ρ u + 2 σ1 w s 1 -σ2 w s 2 B -g 2γ 2γ+g s -,B 4b : -T B α xy E -κ xx (-∂ x T ) B = 2 +P ne 3d σ1 w s 1 + σ2 w s 2 B . (IV.5)
Eq.s(IV.5) are pretty remarkable, as they express the stationary linear response to electrical and thermal per-turbations just in terms of the bull electrochemical potentials, with the coupling of the bulk to the edges encoded in Eq.s(IV.1).

We now use ∆T = T A -T B to denote the difference in the fluctuations of the temperature in the ŷ direction. The average temperature T (x)/2 = (T A + T B )/2, fluctuates around the temperature of the bath T = 1 o K. The fluctuation within T (x) is a small fraction of the bath temerature. Eq.s(IV.1,IV.5) determine set of 8 equations in the 9 unknowns T (x), ∆T, s -,A ± s -,B , c -,A ± c -,B , c +,A , ±c +,B , s + , which eventually yields the space dependent relaxation along the x direction, when a steady state perturbation acts at x = 0.

To close the corresponding set of linear differential equations, we add to it the constitutive equation in Eq.(II.12), which quantifies the energy flux in the x-direction when the charge and the thermal imbalance diffuse along the Hall bar. In the presence of charge imbalance, Eq.(III.3) is violated, so that the charge imbalance in the bulk has now be defined as

ρ I = n (c +,A -c +,B -2 s + ) , (IV.6)
with reference to the remarks after Eq.(II.12) [START_REF]referring to the Poisson electrostatic problem N ∂ 2 δc/∂y 2 = 4πρI ( with clear meaning of the symbols)[END_REF] . It follows that constitutive equation, Eq.(II.12), can be rephrased in the present scheme. Extracting from the difference of Eq.s(IV.3) an expression for δJ n • x, which appears on the left hand side of Eq.(II.12), we set

-2∂ x 1 2w (σ 1 -σ 2 ) [(s 1 -s 2 ) A + (s 1 -s 2 ) B ] + 1 2w (σ 1 + σ 2 ) 2 s + = -2 v F V (c +,A -c +,B -2 s + ). (IV.7)
We remind that we assume for simplicity that our bulk potentials satisfy the mirror symmetry between opposite sides of the edges, s +A ≈ -s +B = s + . This is a weak restriction that could be lifted at the cost of clarity, but does not invalidate the core of our arguments and our results. The space derivative on the left hand side of Eq.(IV.7) defines the length scale for diffusion (in units of L), when the space dependence of the unknown potentials c's and s's has been determined.

The set of nine differential equations can be simplified by neglecting the dependence on ∆T , which corresponds to perform the derivation at ∆T = 0. A straightforward but boring derivation ( see Appendices) leads, to first order in the dimensionless model parameter g = g (2γ + g)L 2 t 2 0 , which quantifies the coupling between bulk and edge, to a two-equation set in the unknowns T and ρ I . Assuming that -σ E = ρ u, we get:

κ xx (-∇ 2 xT ) + α xy ∇ x T EL = -4 + P n e 3d v F n V σ 1 σ 2 (σ 2 1 + σ 2 2 ) ρ I L 2 .
(IV.8)

∂ 2 xρ I + 2 w L v F (σ 1 + σ 2 )V ∂ xρ I -g ρ I = - e n L n + P 3d σ 1 -σ 2 σ 1 + σ 2 w 4 ∂ 2 x 1 σ 1 - 1 σ 2 κ xx + 1 σ 1 + 1 σ 2 α yx + P n e 3d (-∂ xT ) + 1 σ 1 - 1 σ 2 κ xx α xy κ xx E L T . (IV.9)
The system is written in terms of the dimensionless space coordinate x = x/L. L is a length unit at B = 2 T esla that will be found in Section V. It is assumed to depend on B -1 , as the classical cyclotron radius r c = p F c/eB. The length scale of variation of ρ I in Eq.(IV.9) is a -1 L, where

a ≡ 2 wL (σ 1 + σ 2 )V v F , (IV.10)
fixes the length scale for the particle/hole imbalance relaxation. Note that Eq.(II.11) makes the parameter a independent of d. ( + P) 3d is the enthalpy per unit volume. An important parameter is the effective carrier density n 3d ∼ N 3d /V, whose determination we discuss at length in Appendix C. Here, we choose to define the thermal energy per particle (and unit charge) as:

Q ≡ ( + P) 3d V eN 3d ≡ + P n e 3d .
(IV.11)

Choosing Q = 0.12 × 10 -3 V for the enthalpy per unit particle, we are able to match reasonably well the val-ues expected from a Fermi Liquid approach and those experimentally found for α xy , κ xx and e N (see Appendix C). (n ≡ n 3d ∼ N 3d /V will have to be self-consistently determined (see Eq.(V.10))). E appearing in Eq.(IV.8) is the transverse applied electric field close to x = 0: it depends on ρ I away from x = 0, and should be self-consistency recovered from the transverse potential generated by ρ I at any distance. For the sake of simplicity, here we linearize the corresponding equations, by considering two main x-regions: one close to the contacts where the external current J y = σE y is applied ("near zone") and one "far" from x = 0 ("far zone"), where the influence of the externally applied field vanishes. E ≡ E y is the applied field at x ≈ 0. We keep E = 0 constant within each zone. In the far zone E = E nl is self-consistently determined, giving rise to the non local potential v nl .

At x = 0, the boundary condition for ∂ x T (x = 0) is provided by the Nernst coefficient:

∂ x T x=0 = sign(B) (σ -1 α) xy E y (x = 0). (IV.
12)

It will be shown in Section that this boundary condition together with Eq.(V.3) implies δµ I ∼ k B T (see Eq.(V.7)).

Our derivation has ignored the relaxation of ∂T due to inelastic scattering processes including acoustic phonon excitation Away from the origin, the length scale for these relaxation processes is ∼ 10 µm 44 .

B. Solution of Eq.s(IV.8,IV.9) for x ∼ 0 To solve Eq.s(IV.8,IV.9) for x ∼ 0, we set x = x/L and introduce a dimensionless temperature τ = (T / o K), as well as a dimensionless imbalance energy δµ I = eρ I L 2 /(k B o K). As a result, we can write the system of differential equations for τ and δµ I as

∂ 2 x τ -∂ x (θ τ ) = β δµ I , (a) ∂ 2 x δµ I + a ∂ x δµ I -g δµ I = ∂ 2 x {k ∂ x τ -θ u τ } . (b) (IV.13)
Eqs.(IV.13) depend on the dimensionless parameters

θ ≡ sign(B) α xy E L κ xx , β ≡ + P n e 2 3d k B κ xx v F n V 4 σ 1 σ 2 (σ 2 1 + σ 2 2 ) , (IV.14)
as well as on the parameters k and u, whose full expression is given in Eq.s(D.1,D.2), which depend on the imbalance ratio q,

q = σ 1 -σ 2 [σ 2 1 + σ 2 2 ] 1/2 , (IV.15)
that is linear in the difference between the particle and the hole conductivities. The longitudinal conductivities σ i 's are fairly isotropic and, therefore, they are expected to be pretty insensitive to the orientation of the magnetic field B. However, we expect that the difference in chemical potential between region A and region B changes sign when B is flipped, as flipping B implies exchanging particles and holes with each other. This effect, which we discuss in the following, corresponds to breaking the symmetry under B → -B. For x ∼ 0, it is E ∝ θ that fully determines the current in the ŷ direction. Consistently with the experimental data, we assume E w = 10 -4 V takes the value given by the experiment. Moreover, we also assume B > 0 (the case B < 0 will be analyzed in Section VI, when discussing the non local voltage). It is important to remark that changing the sign of B corresponds to having q → -q, that is, to exchange with each other the Hall dynamics of particles and holes.

A more quantitative analysis will be given in the Appendices, by establishing the numerical estimate of the parameters involved. Here we discuss the general features and our specific approach to the solution of the problem.

The bath temperature of 1 K corresponds to the value in the experiment. The longitudinal conductivities σ 1,2 are fixed by assuming the reference resistance R ∼ 10 kΩ and the inverse conductivities require an aspect ratio 1/σ 3d = R d/w, according to Eq.(II.11). Inserting ∼ 4 w ∼ 400 µm, d ∼ 1 Å, L = 0.1 µm and the Fermi velocity v F = 10 8 cm/sec in the parameter a defined in Eq.(IV.10), gives a ∼ 10 -3 , so that the scale of variation for ρ I is a -1 L ∼ 100 µm. The appropriate order of magnitude for the length scale L is self-consistently determined in the next Section.

If the edge/bulk leakage parameter g ∝ 10 -6 , it can be shown (see Eq.s(D.6) in the Appendix D) that, employing the parameters defined in Eq.(D.8), the system Eq.(IV.13) can be cast in a form in which all known quantities are O(1), which is particularly amenable for drawing plots from the numerical data. Nevertheless, for the general discussion here we will keep using Eq.s(IV.13) as our reference, as they appear to be more appealing for the sake of the physical interpretation of the results.

We now discuss the general features of the system by using an approximate analytical solution in the region x ≈ 0 and accounting for the boundary conditions of Eq.(IV.12).

We show in Appendix E that, inserting Eq.s(IV.13,a) for T into Eq.s(IV.13,b), one recovers a higher order equation for δµ I . In dimensionless units, this can be cast in the form of set of first order differential equations as: 16) with the additional equation for τ = T / o K:

   ∂ x δµ I = ξ ∂ x ν = δµ I ∂ x ξ = -(a -θ -k β) ξ + [(a -u β) θ + g] δµ I -g θ ν, , (IV.
∂ 2 x τ -∂ x (θ τ ) = β δµ I . (IV.17)
It is useful to shift a -u β → a and to define F = (aθ -κ β) and G = a θ + g, where κ = k -u. The boundary conditions are satisfied by F > 0 and G, θ, g < 0. The system of Eq.(IV.16) provides three possible solutions. The corresponding eigenvalues characterize their decay rate in real space when moving away from the applied perturbation. One eigenvalue is real and the other two may be still real or complex conjugate. It is remarkable that in both cases the eigenvalues λ 1,2 are independent of θ, if and only if κ ∝ q = 0. Indeed, it is shown in Appendix E that, for κ = 0, the three eigenvalues are given by:

λ 1,2 = - a 2 ± a 2 4 -|g|, λ 3 = -θ sign[g]. (IV.18)
If a 2 4 -|g| > 0, the three eigenvalues are real, otherwise λ 1,2 are complex conjugate. Let us consider the case κ = q = 0 for a while. It follows that, for κ = 0, only the solution corresponding to λ 3 is relevant. It takes the form

δµ I = ηθ 2 e θx (IV .19)
(with η to be determined in the following). Here, we assume θ < 0 and x > 0. In any case, we will always choose the decaying solution for x > 0, which fixes the sign of E for a given orientation of B orthogonal to the strip. From Eq.(IV.17), we obtain:

τ = 1 + η β θ x e θ x.
(IV.20)

Imposing the boundary condition at x = 0 as per Eq.(IV.12), to O(θ 2 ) we obtain:

∂ xτ | x=0 = η β θ = κ xx o Kα xy (σ -1 α) xy θ (IV.21)
which gives βη ∼ 40, once the parameters are chosen as described in Appendix C. From the definition of β, the boundary condition Eq.(IV.21) can be written as:

η Q e o Kα xy (σ -1 α) yx κ xx k B κ xx 2 v F n V 2 σ 1 σ 2 (σ 2 1 + σ 2 2 ) ∼ 1. (IV.22)
This provides an estimate for η ≈ 2.35 × 10 11 . Our definitions of the physical quantities are fully consistent, as we show in the following.

Let us now add some charge imbalance at time t = 0. For E = 0, a thermal flux -∂ xδJ Q = -e Q ∂ xδJ n moves from x = 0 into the bulk, parallel to the x axis. From Eq.(II.5), by differentiating Eq.(IV.20) and using the definition of θ in Eq.(IV.14), one obtains

-∂ xJ Q | x=0 = κ xx ∂ 2 xT x=0 -α xy EL ∂ xT /2L = κ xx η βθ 2 × o K/2L = (∂ xT )α xy E/2.
(IV.23) (We have divided by 2 because we only consider the flux in the direction x > 0.) ∂ xJ Q is no longer vanishing as from Eq.(IV.13 (a)) for µ I = 0. Instead, it now generates a thermal gradient.

On the other hand, by substituting v F /(n V) from Eq.(II.12) in Eq.(IV.22), we get:

-e Q ∂ x δJ n L 2 α xy (σ -1 α) yx κ xx 2 κ xx 2 σ 1 σ 2 (σ 2 1 + σ 2 2 ) = θ 2 , (IV.24)
where we have identified the fluctuation of the chemical potential, ρ I L 2 /(k B o K) ≡ δµ I = ηθ 2 , on the r.h.s. from Eq.(IV.20). Using Eq.(IV.23), Eq.(IV.24) reads:

L ∂ x T α xy (σ -1 α) yx κ xx 2 σ 1 σ 2 (σ 2 1 + σ 2 2 ) ≈ θ (IV.25)
which gives ∂ xT ∼ 40 θ o K once more, as in Eq.(IV.21) because the ratio 2 σ 1 σ 2 /(σ 2 1 + σ 2 2 ) ≈ 1.

V. ETTINGSHAUSEN PARAMETER

The first step is to estimate the leakage factor g and, consequently, the actual number of free carriers, N 3d , that have been redistributed between the areas A and B of the Hall bar. In dimensionless units, the energy associated to the transverse voltage across the Hall bar (in the ŷ-direction) is given by:

v y = e (c +,A -c +,B )/(k B o K). (V.1)
Following the derivation of Appendix F, we get, to O[q 2 ] and O[g] :

∂ 2 x v y -g v y = g e L 1 σ 1 + 1 σ 2 w 4 k B Q α xy EL T o K -κ xx - σ 1 -σ 2 σ 1 + σ 2 α xy Q ∂ xτ . (V.2)
In order to establish the consistency, in this Section we concentrate on the area around x ≈ 0 where the external bias is applied. We will fix N 3d , together with the scale L, by assuming an applied electric field at the origin E ∼ V /m. The r.h.s. includes a term ∝ [α xy EL T / o K -κ xx ∂ xτ ] which describes the particle flux J n flowing away from x ≈ 0. As metallic contacts are applied at x ≈ 0, we expect a flow in the contacts of the charge carriers lost because the first two terms. The last term ∝ q ∂ xτ relates the transverse voltage to the thermal gradient in the presence of the orthogonal magnetic field. A thermal gradient in presence of a current with a magnetic field is named Ettingshausen effect. We now focus on it as, in our case, it is proportional to the charge imbalance. If we assume that, close to the origin, the external source sustains the particle flux ∂ xJ n given by Eq.(IV.23) into the contacts, so that ∂ x [∂ xτ -θ τ ] ≈ 0, according to the definition of θ given by Eq.(IV.14), an estimate of the Ettingshausen parameter can be obtained. Self-consistency also allows to determine the effective 3d carrier density N 3d , as we show in the following.

For x ∼ 0, we can neglect the term ∂ 2

x v y at the l.h.s of Eq.(V.2). Within this approximation, g drops out for x ∼ 0, where the charge distribution is fixed by the applied electric field. Of course, g will play an important role in the next Section, when estimating the non local transverse voltage. By keeping just the last term in Eq.(V.2), which we denote by v ⊥ , we get

v ⊥ ≈ e L q 1 σ 1 + 1 σ 2 w 4 k B α xy ∂ xτ. (V.3)
It is important to realize that, for x ∼ 0, where the current is fed in, v ⊥ does not flip its sign when B changes sign. This is because α xy is an odd function of B and, when B changes sign, particles and holes exchange their position between area A and B, so that q changes sign, as well. Eq.(V.3) can be usefully rewritten, according to the definition of β given by Eq.(IV.14), as

v ⊥ ≈ β-1 a q n L 2 Q (σ 1 + σ 2 ) 2 2(σ 2 1 + σ 2 2 ) α xy κ xx ∂ xτ. (V.4)
where a is defined in Eq.(IV.10). To set up the selfconsistency, we insert v ⊥ = e qEw/(k

B o K), EL = (κ xx /α xy )θ
in Eq.(V.4). From Eq.(IV.25), observing that

σ 2 1 +σ 2 2
2σ1σ2 ∼ 1, we get:

κ xx α xy θ 1 k B w L ≈ β-1 a n L 2 Q e θ (σ -1 α) xy , (V.
1 = n L 2 = 10 10 cm 2 (0.1 µm) 2 . (V.9)
Here further consistency requires that

ρ I e = 14 × 10 12 cm -3 = q n 3d . (V.10)
With a density n 3d = 1.3 × 10 15 cm -3 , at B = 2 T esla, the bulk Hall conductance is

σ B = n 3d ec B = 10 12 sec -1 , (V.11)
a value roughly consistent with the ratio κ F L xx /α F L xy given in Eq.(C.6), which provides σ B ≈ 5 × 10 12 sec -1 . According to Eq.(V.10), this density requires q ∼ 10 -2 . The inequality σ B /σ 3d << 1 is therefore confirmed.

The electrostatic energy per unit volume due to the charge fluctuation is:

e 2 n 3d ρ I L 2 = δµ I n 3d = 1 2 ηθ 2 n 3d ≈ 46.7 m eV cm 2 Å . (V.12)
It is remarkable that Eq.(V.3) entails the definition of the Ettingshausen parameter, as we show in the following, by resorting to the Fermi Liquid forms of the transport parameters α F L xy and κ F L xx . The Ettingshausen ratio is:

P E = ∂ x T |B|J y (V.13)
Now, substituting the parameters given in Eq.(II.9,II.8) into Eq.(V.4), with

α F L xy κ F L xx = σ B σ 3d 1 Q 1 1 -σ B σ 3d 2 ,
we get:

v ⊥ = β-1 q a n L 2 (σ 1 + σ 2 ) 2 2(σ 2 1 + σ 2 2 ) σ B σ 3d 1 1 -σ B σ 3d 2 ∂ xτ Posing e qEw/k B o K = v ⊥ again, and dropping 2 (σ 2 1 + σ 2 2 ) (σ 1 + σ 2 ) 2 1 - σ B σ 3d 2 ≈ 1, we are left with e k B n L 2 βa -1 w L = σ B ∂ x τ o K σ 3d E , (V.14)
or

P F L E = 1 B ∂ x τ o K σ 3d E = 1 k B c βa -1 n L 2 1 n 3d w L . (V.15)
With B = 1 T esla and βa -1 ∼ 10 -6 ,

P F L E = 1.1 × 10 -14 cm 3 eV c o K = 1.7 × 10 -2 meter o K Amp T esla , (V.16)
as Amp/meter = 8.97 × 10 9 V /sec. This is our first result, to be compared with the case of Bismuth: P E Bismuth = 7.5 × 10 -4 meter o K Amp T esla . At the bottom of Section VI we argue that consistency with the measured NLR magnitude ( R N L ∼ kΩ) suggests that the actual value of P E far from the source should be lowered of about two orders of magnitude with respect to the one reported in Eq.(V.16). Given the uncertainities in the effective bulk density, our result cannot be sharper.

VI. THE TRANSVERSE NON LOCAL VOLTAGE

We now turn back to Eq.(V.2) and examine it for x ∼ D, far away from the origin. Defining

r = sign(B) e L 1 σ 1 + 1 σ 2 w 4 k B Q κ xx , (VI.1)
Eq.(V.2) reads:

∂ 2 x v y -g v y = g r α xy κ xx E nl L T o K -1 - σ 1 -σ 2 σ 1 + σ 2 α xy κ xx Q ∂ xτ , (VI.2)
where T = 1 o K. For x ∼ D, there is no applied electric field. However, as we explain in the Introduction, the thermal gradient at the origin determines a charge imbalance that propagates in the x direction with a very low relaxation rate. Once transported at x ∼ D, the charge imbalance generates the electric field E nl , which is the source of the B-dependent contribution to the non local resistance. Note that the product α xy E nl does not depend on sign(B), although both factors do. A consistency condition can be set in Eq.(VI.2) by posing v y = eE nl w/(2k o B K) (the factor 1/2 arises from the observation that the electric field is not given by an external source but created when the charge difference accumulates close to the edges). Switching to the variables x = x/a -1 L and τ = aτ , we get:

∂ 2 x v y + g a -2 2 r α xy κ xx L w k B T e -1 - 1 τ 2 Q v y = -g a -2 r 1 -q α xy κ xx Q ∂ xτ . (VI.3)
In Eq.(VI.3), we have added an extra term -v y /τ 2 Q , to account for the intrinsic relaxation. In any case, even in the absence of this term, v y is a decaying function at large values of x, where the thermal gradient goes to 0, as a consequence of the fact that g a -2 < 0 and that the quantity within the square brackets is positive. The l.h.s. of Eq.(VI.3) implies a decay length λ of the non local voltage. For τ Q → ∞, one obtains

λ ∼ σ 1 + σ 2 2 σ 1 σ 2 T α xy Q -1 -1/2 L/ g. (VI.4)
We argue that this length is well defined, because the difference in the square bracket is, on very general grounds, greater than zero. Recalling that the Nernst coefficient e N = σ -1 3d α xy ∼ E y /∂ x T and that the orbital motion of carriers driven by the magnetic field is roughly circular, we recognize a torque acting on one carrier which generates a work L = δx eE y per unit angle. On the other hand, in the same wedge the energy per particle provides a thermal contribution to the change in the free energy ∼ s δT , where s = e Q/T is the entropy per unit particle. Ultimately, the decay rate is related to a change of the free energy per particle δf where f = u -T s is the drift due to the electric field. In fact, δf = -L -s δT . As δf has to be negative in order for the system to evolve towards equilibrium, one obtains L + sδT > 0, with L < 0 (as it is performed by the external source) and sδT > 0 as δT 0 is induced by the source.

Indeed, one obtains

2r α xy κ xx L w k B T sign(B) e = σ 1 + σ 2 2 σ 1 σ 2 T |α xy | Q > 1 (VI.5)
for T = 1 o K and σ 1 = 7.19×10 14 , σ 2 = 4.79×10 14 , while, on the r.h.s. of Eq.(VI.3), Q α xy /κ xx ≈ 0.12 × 10 -3 .

It is important to note that the value chosen for the energy per particle, e Q, is ∼ k B T , but, in presence of magnetic field, the magnetization energy -M B has to be added to it, which makes the argument of the square root in Eq.(VI.4) even more positive, when B increases.

It is also remarkable that the closer σ 1 is to σ 2 , the longer the decay length λ is.

The space scale is a -1 L ∼ 100 µm and, in this scale, λ ∼ 1, or even larger. Indeed in our case the enthalpy per particle, including the magnetization work, as well, is given by eQ = 0.12 10 -3 -4 × 10 -4 |b| eV , where b is the magnetic field in units of 2 T esla.

The inhomogenous term in Eq.(VI.3) dominates at distances D ∼ 3 a -1 L from the origin. In the plots, we have chosen a phenomenological B dependence τ Q ∝ B -1/2 and we have considered a Hall bar which is infinitely long in the x > 0 direction.

In deriving Eq.(VI.3), B was assumed positive. For negative B values, α xy ∝ 1/B changes sign, while κ xx is, in general, an even function of B and, in particular, here it is taken independent of B. No other functional dependence on B is introduced, except for τ Q ∝ 1/ |B|. However convergency of the solution of Eq.(IV.19) requires that θ of Eq.(IV.14) is an even function of B. The prefactor r as well as the imbalance q, and σ B of Eq.(V.11) are odd in B, because particle and hole exchange their role by flipping the magnetic field. In solving the differential system, the initial condition for the integral of the imbalance chemical potential, ν, is also odd: ν(0) = sign(B).

Fig. (5) displays the main result of this work. We describe here qualitatively the picture that emerges from the model with the help of the plots of the relevant quantities.

According to the sketch in Fig. (2), the applied electric field E y and the corresponding current J y are oriented from edge B (the lower edge of the picture) to edge A (the top edge of the picture), independently of the orientation of the magnetic field B. When B > 0, the carriers leaking in the bulk from the edges moving toward x > 0 are particles close to edge A and holes close to edge B. The carriers moving in the other (opposite) direction, impinging on the Hall bar boundary, have opposite sign, but are assumed to be absorbed by the boundary and, so, they do not enter our discussion. For x ∼ 0, the system is assumed to be thermalized by the boundary, but an increase of temperature with respect to the thermal bath is expected, due to the Joule heat accompanying the applied current. Let us consider the case B > 0 first. Nernst effect provides a thermal gradient which moves the carriers away from the applied field region. The increase of temperature drops relatively fast away from the origin. This is reported in Fig. (6) showing the thermal gradient. At large distances the temperature decreases at a rate decreasing with the distance, till it becomes constant. Under the effect of the thermal gradient, the carriers diffuse in the Hall bar as proved by the space dependence of the difference in chemical potential between edge A and edge B, δµ I , which is plotted in Fig. (7) at increasing magnetic fields. As the relaxation time across the bar is rather long, the carriers can reach regions of the Hall bar where the effect of the applied field has vanished, so that δµ I keeps finite also at distances D ≈ 6 a -1 L. This implies that a voltage difference develops between edges which is at even with the applied voltage. For B > 0, this gives a positive non local voltage, which is reported in Fig. (8).

Let us now assume B < 0. In this case the positive and negative carriers leaked in the bulk of the Hall bar and drifting at x > 0 have opposite sign with respect to sketch i), so that a positive δµ I develops between the two boundaries (see Fig. (9)) and the corresponding voltage difference is at odd with respect to the applied one. This induces heat diffusion away from the injected current and the thermal gradient is opposite to the one of Fig. (6), as plotted in Fig. (10). This fact increases the distances at which the perturbation diffuses. The result is that not only the sign of the voltage correction is opposite (see Fig. (5) for negative magnetic fields), but there is a marked difference in amplitude between the contributions coming from opposite orientations of the magnetic field .

While the non local voltage has an exponential decrease not far from the origin (see Fig.( 8)), it is power law far from the origin as shown in the inset of Fig. (8). The decay toward zero of the imbalance chemical potential for larger magnetic fields is slower the larger the magnetic field is, as shown in Fig. (11) and in Fig. (12). They are opposite in sign, but the first one adds up to the applied voltage, the second one is subtracted. It is remarkable that the relative damping ratio with distance is clearly weaker for B < 0.

The contribution to the non local voltage of Fig. (5), appropriately scaled, according to the experiment, by fixing the constant energy scale k B T 0 , adds up to the non local voltage difference derived in Ref. [29]. The corresponding transverse resistance arises from edge states only, which, along their path, suffer some dissipation at the various contacts of the Hall bar. While the contribution coming from the edges is fully symmetric with reversing the orientation of the magnetic field, the contribution derived here introduces an asymmetry which is found in the experiment. The asymmetry tends to reduce with increasing D, what is fully consistent with our plots.

The differential system is fully homogeneous. Unknown functions depend on an overall scale normalization. Close to the origin, the magnitude of ∂ x τ is determined by E, field according to the Nernst boundary condition of Eq.(IV.12). This is the reason why we could extract the Ettingshausen parameter in the near zone. To extract the magnitude of the NLV away from the origin, the differential equation system should have been solved allowing for a selfconsistent E field at any distance (θ → θ(x)) . This has not been done. Therefore, in the far zone the magnitude of the NLV is strictly speaking undetermined. We could try to use our estimate of the Ettingshausen parameter of Eq.(V.16) to determine the normalization, by using nAmp as a probe current and 1 T esla as the reference magnetic field, and extracting a temperature derivative ∂T /∂(a -1 L) ≈ 10 -2 o K/(a -1 L). However, when we plug this in Eq.(VI.3) and we evaluate the NLR by means of the same current, we find that our estimate is two-three orders of magnitude larger than the measured one, amounting to kΩs.

VII. SUMMARY AND FINAL REMARKS

Near the Dirac point, graphene Hall bars in an applied orthogonal magnetic field of the order of a few Tesla, exhibit large nonlocal resistances which cannot be fully explained within the framework of the "ideal" Quantum Hall Effect (QHE) in graphene. Here, we discuss nonlocal resistances in a Hall bar configuration corresponding to an open circuit in the longitudinal direction, parallel to the length of the Hall bar. In these conditions, one would expect zero voltage between contacts in the transverse direction, and that a transverse current injected at one boundary of the Hall bar, driven by an external electric field E gives rise to a local resistance that exponentially vanishes far from that boundary. Instead, the decay is power-law and there is a marked asymmetry between positive and negative values of the B field, which is reduced at increasing distance from the injection point. Weak spin orbit coupling in graphene excludes that the NLR is due to Spin Hall Effect 30,45,46 . The model of McEuen et al. 32 has been used by one of us and collaborators to describe qualitatively and semi-quantitatively the experimental results 29 . The model rests on the assumption that a non insulating bulk plays a role in transport close to the CNP. This is actually consistent with the model by Abanin et al. 30 , which we review in Section III and accounts for some dissipative effects measured in graphene at the CNP by assuming leakage of edge carriers into the bulk. The ν = 0 Landau level would be splitted by Zeeman interaction and a sort of plateaux in the xyconductance would survive at the CNP accompanied by a sizeable longitudinal bulk resistance. Indeed, graphene on SiC, in which QHE was measured with a metrological precision, is relatively disordered, with strong intervalley scattering. However, these models are unable to catch the asymmetry under a change in sign of the applied the magnetic field.

In this work we have reconsidered the possibility that dissipation in the bulk can be accompanied by local charge imbalance and thermal gradients diffusing across the Hall bar, thus being responsible for the NLR. We can expect that even in the presence of relatively strong magnetic field, close to the CNP, in the high-mobility bulk graphene, the Fermi surface vanishes, leading to ineffective screening and linear energy dispersion, which implies the formation of a strongly-interacting quasi-relativistic electron-hole Dirac fluid 12 . In the relativistic approach, which is reviewed in Section II, particle current with energy transfer can be a neutral excitation mode independent of the charge transfer. The linear energy dispersion is also responsible for long recombination relaxation times of the quasiparticles, so that local charge imbalance can persist during the diffusion. We propose that the McEuen model accounts for the large part of the NLR, but long lived thermal effects occurring at low temperature in the bulk of the Hall bar may add a relevant contribution which, in particular, could be the source of the asymmetry when flipping B.

In fact the McEuen model requires bulk conductivities σ 1 , σ 2 for particles and holes respectively, which are not the same.

On these side effects we focus on. At low temperature (T ∼ 1 o K), with local non equilibrium induced by a locally applied transverse electric field, a few-Tesla B induces space separation of oppositely charged carriers leaking from opposite edges, with little recombination relaxation rate. The charge imbalance generated in the bulk of the Hall bar by Joule heating close to the source is drifted far from the injection point by a Nernst thermal gradient at open circuit. In our framework, this is monitored by the space dependence of the imbalance chemical potential which does not decay exponentially along the Hall bar over distances comparable with the sample length. The bulk charge imbalance accumulated at opposite edges turns into a non local voltage that can be probed by an external transverse current which produces the NLR. The contribution discussed here appears in Fig. (5) and is strongly asymmetric under changing the sign of B.

Various remarks are in order. As mentioned at the end of Section VI, we assume that the applied perturbation is localized close to the origin and only the drift of charges induced by the thermal gradient is accounted for in solving the diffusion equations derived in the hydrodynamical limit, far away from the injection point. Self-consistency is invoked to fix the non local voltage in the absence of sources. This implies that, far from the source, the set of differential equations is homogeneous and the solution is up to an unknown constant energy scale k B T 0 . We have tried to estimate this constant in conjunction with the Ettingshausen parameter, but it can be of course fitted from experiment. Once the polarization of the electric field E is fixed, the response of the system to opposite directions of the applied B field is certainly asymmetric, as particles are exchanged with holes when flipping the orientation of the magnetic field. The asymmetry is stronger close to the injection point. The asymmetry with B found experimentally can help in fixing the ratio between the contribution coming from Eq.(VI.3) and the one coming from the McEuen model.

The charge imbalance arises from the spatial separation of charges, as it can be seen from the sketch of Fig. (2). In fact, the magnetic field tends to keep carriers of opposite sign far away from one another. This reduces the chances of recombination. However some charge imbalance also originates from the difference between the bulk conductivities σ 1 , σ 2 of the opposite charge carriers, which could be the source of a reduction of the Coulomb interaction between charges. This difference parametrized by the variable q of Eq.(IV.15) increases the asymmetry with the magnetic field orientation as q is odd with B and appears explicitely on the r.h.s. of Eq.(VI.3). It is remarkable that, according to our model, the closer the bulk conductivities are, the longer is the decay length λ of the NLV, as can be seen from Eq.(VI.4).

We assume that the Dirac model, in presence of two carriers of opposite charge (at the CNP), is the appropriate low energy framework to discuss the thermal, as well as the charge transport. Hence we only consider the lowest bath temperatures attained in the experiment which are T ∼ 1 ÷ 1.7 o K. In Eq.(V.7) an estimate is given of the charge imbalance chemical potential δµ I in units k B o K and is comparable with the bath temperature. Being estimated close to the origin, where the perturbation acts, δµ I ∼ 2 appears as the maximum attainable value. Hence we are at the limit of validity for our approach, because, in this range of values, the Dirac model gives up by turning into the more appropriate Fermi Liquid model. In passing, we note that, at these temperatures, similar, but p doped, graphene/SiC samples have been used as barrier for Al/graphene/Al superconducting Josephson junctions, showing a remarkable hysteretic loop of the critical current when the magnetic field is cycled and an Aslamazov-Larkin (AL) paraconductivity temperature dependence consistent with a Fermi Liquid picture 47,48 .

In the Dirac model the relativistic stress energy tensor provides the energy flux T 0 α from the four-momentum conservation. The full approach does not contain more general dissipative terms except in Eq.(VI.3) for v y , in which a phenomenological dissipative term has been added. The recombination time τ Q ∼ 1/ √ B has been measured and found rather long 36,49 . This extra term has been added to account for the relatively fast decay of v y with increasing B. If we neglect this term, the decay is determined by λ given by Eq.(VI.4) which gives a weaker dependence on B.

The numerical values for the electrical bulk conductivities are derived from Ref. [29], while the numerical values chosen for the thermal transport 3-d parameters (Nernst and Peltier parameters, thermal conductivity) have been found in the literature. In various occasions we have compared those values with the ones derived from 3-d hydrodynamics of a Fermi Liquid 14,41 . This has given us the chance of fixing the order of magnitude of two quantities that are relevant in thermal transport, but are difficult to extract from experiment: the enthalpy per particle eQ and the particle density in the effective 3d system n 3d . Thermal transport is an essentially 3-d phenomenon and the knowledge of the sheet carrier density of graphene n , as determined by the QHE, is not enough. The value we have chosen for eQ is not far from k B T plus the magnetization work -M B which freezes the kinetic energy term. This reduces λ with increasing B.

Our results cannot be valid in the limit of small magnetic field (B 1 T esla). Outside of the QHE regime the model, in the absence of edge modes, is invalid from the outset. However the trend when the magnetic field is reduced is qualitatively well represented by our results of Fig. (5). Of course the transport parameters should be recalculated in a Dirac frame and in presence of QHE quantization of the conductivity. Nevertheless, hydrodynamics is still valid at very low energy (very low temperatures), lower that the Zeeman splitting at ν = 0. The starting guess is that the bulk of the Hall bar is not totally insulating even at such a small working temperature.

We conclude by noting that further studies of the thermal transport in chemical vapour deposited graphene can help by making the fractional quantum Hall effect an efficient and controlled playground for applications in alternative quantum information processors.

1 σ 1 - 1 σ 2 α xy ∆T E Q - w 4 1 σ 1 + 1 σ 2 g 2γ 2γ + g (s -,A + s -,B ) (A.1) s -,A + s -,B = -K -(-∇ x T ) - w 4 1 σ 1 - 1 σ 2 α xy T E Q -g 2γ 2γ + g w 4 1 σ 1 + 1 σ 2 (s -,A -s -,B ) (A.2) 2s + = -K + (-∇ x T ) - w 4 1 σ 1 + 1 σ 2 α xy T E Q - w 4 1 σ 1 - 1 σ 2 g 2γ 2γ + g (s -,A -s -,B ) , (A.3)
where

K ± = w 4 1 σ 1 ± 1 σ 2 κ xx Q + 1 σ 1 ∓ 1 σ 2 α yx . (A.4)
Singling out the -2∂ x s + term from Eq.(IV.7) and substituting it into the derivative of Eq.(A.3), one gets:

K + (-∇ 2 x T ) + w 4 1 σ 1 + 1 σ 2 α xy ∇ x T E Q + w 4 1 σ 1 - 1 σ 2 g 2γ 2γ + g ∂ x (s -,A -s -,B ) = σ 1 -σ 2 σ 1 + σ 2 ∂ x (s -,A + s -,B ) - 2 v F w (σ 1 + σ 2 ) 1 n V ρ I .
In the comparison with the derivative of Eq.(A.2):

∂ x (s -,A + s -,B ) = -K -(-∇ 2 x T ) - w 4 1 σ 1 - 1 σ 2 α xy ∇ x T E Q -g 2γ 2γ + g w 4 1 σ 1 + 1 σ 2 ∂ x (s -,A -s -,B ) ,(A.5)
the contribution from ∂ x (s -,A -s -,B ) cancels out, giving eventually the final form of the equation we are looking after:

κ xx (-∇ 2 x T ) + α xy ∇ x T E = -4 + P n e 2 3d 1 n V v F σ 1 σ 2 (σ 2 1 + σ 2 2 ) ρ I , (A.6)
which can be cast in the form (θ ≡ α xy E L/κ xx ):

∂ 2 xτ -θ ∂ xτ = βδµ I , (A.7) β = + P n e 2 3d v F n V k B κ xx 4 σ 1 σ 2 (σ 2 1 + σ 2 2 ) × 1 - σ B σ 3d 2 -1 , (A.8)
where

δµ I = eρ I L 2 k B o K and τ = T / o K. ∂ x (c +,A -c +,B ) = -(2 γ + g )(c -,A + c -,B ) + g (s -,A -s -,B ) (B.1) ∂ x (c -,A -c -,B ) = -g (c +,A + c +,B ) (B.2) ∂ x (c +,A + c +,B ) = -(2 γ + g )(c -,A -c -,B ) + g (s -,A + s -,B ) (B.3) ∂ x (c -,A + c -,B ) = -g [c +,A -c +,B -2 s + ] (B.4) ∂ 2 x (c -,A -c -,B ) = g (2γ + g )(c -,A -c -,B ) -g 2 (s -,A + s -,B ) ∂ 2 x (c +,A -c +,B ) = g (2γ + g ) [c +,A -c +,B -2s + ] + g ∂ x (s -,A -s -,B ) (B.5)
Using Eq.(B.5) and Eq.(IV.7), we also obtain:

∂ x (c +,A -c +,B -2 s + ) = -(2γ + g) t 0 (c -,A + c -,B ) + g t 0 (s -,A -s -,B ) - 2 w v F (σ 1 + σ 2 )V (c +,A -c +,B -2 s + ) + σ 1 -σ 2 σ 1 + σ 2 ∂ x (s -,A + s -,B ) (B.6)
The parameters g and γ have dimensions [1/(Lt 0 )] where t 0 is the relaxation time scale that we have introduced in the definition of the dimensionless leakage parameter from the edge states to the bulk, g = g(2 γ + g)L 2 t 2 0 to properly match the dimensions ( see below for details). By differentiating Eq.(B.6) and by plugging into it Eq.s(A.1), keeping terms up to O(g 2 ), one recovers the desired equation, connecting the imbalance potential with the thermal excess temperatures, given by:

∂ 2 x (c +,A -c +,B -2 s + ) + 2 w v F (σ 1 + σ 2 )V ∂ x (c +,A -c +,B -2 s + ) -g (2γ + g) τ 2 Q (c +,A -c +,B -2 s + ) + σ 1 -σ 2 σ 1 + σ 2 ∂ 2 x K -(-∂ x T ) + w 4 1 σ 1 - 1 σ 2 α xy E Q T = -g τ Q 1 - 2γ τ -1 Q 2γ + g w 4 1 σ 1 + 1 σ 2 σ 1 -σ 2 σ 1 + σ 2 ∂ x K -(-∂ x ∆T ) + w 4 1 σ 1 - 1 σ 2 α xy E Q ∆T , (B.7)
Note that the expression in curly brackets is the same, both for T and for ∆T . One can assume that one of the two is constant. We neglect any ŷ-dependence of ρ I except for the mirror symmetry, so to make it independent of ∆T = T A -T B . Hence the right hand side is put equal to zero:

K -(-∂ x ∆T ) + w 4 1 σ 1 - 1 σ 2 α xy E Q ∆T = cnst, (B.8)
where the constant can be put to zero. Our choice is consistent with the initial choice s +,A = -s +,B = s + . Finally, multiplying both sides by e n L 2 /k B o K and resorting to the dimensionless space coordinate x = x/L, we get:

∂ 2 xδµ I + 2 w L v F (σ 1 + σ 2 )V ∂ xδµ I -g δµ I = - e k B n L 2 Q σ 1 -σ 2 σ 1 + σ 2 w 4 L ∂ 2 x 1 σ 1 - 1 σ 2 κ xx + 1 σ 1 + 1 σ 2 α yx Q (-∂ xτ ) + 1 σ 1 - 1 σ 2 κ xx α xy κ xx E L τ . (B.9)
is lower than the chemical potential, the expressions for the transport coefficients α F L xy and κ F L xx are:

α F L xy = + P n e 3d σ B T . κ F L xx ≈ + P n e 2 3d σ 3d T 1 - σ B σ 3d 2 . (C.1)
These are quantities defined in 3-d, where their 3 dimensions are respectively given by [α] ∼ e/(Lt T ) and by [κ] ∼ e 2 /(L 2 t T ). On the other hand, in the range T = 1 ÷ 10 o K the value reported for graphene is 50 α xy ≈ 75 nAmp/ o K, which is a 2-d quantity. As we consider thermal effects in the graphene grown onto SiC as a "bulk " phenomenon, we extend the definition of α xy to 3d by dividing it by ∼ 0.3 nm, so that we eventually obtain:

α xy ≈ 2.24 × 10 11 V /( o K sec). (C.2)
Ref. [41] provides a numerical estimate ( + P) ∼ 0.12 × 10 3 (T / o K) 3 eV /cm 2 , also given in Ref. [18], when the chemical potential exceeds thermal energy, even at temperature 40 o K. The question is how to relate this quantity to the corresponding 3-d quantity( + P) 3d , that is required for the enthalpy per unit (3-d) volume.

We make the following choice. We define a total energy per particle by multiplying ( + P) 3d by the volume V = w × 1 Å ≈ 10 -11 cm 3 and dividing it by the effective number of 3-d carriers N 3d ∼ 10 3 and we pose:

Q = ( + P) 3d V N 3d e ≈ 0.12 × 10 -3 V (C .3)
(here and in the following V ≡ V olt). We will show that this value for N 3d is also suggested by our self consistent procedure. N 3d can be interpreted as the effective number for carriers contributing to the thermal processes in the As σ 3d L d-2 = R -1 , an aspect ratio is required, so that σ -1 3d = R d/w ≈ 0.5 × 10 -15 sec, with d = 1 Å and R = 10 kΩ. Here R is the longitudinal resistance and w the area.

Here σ B /σ 3d = e 2 R/h • [hc/(eB w)] • /w(note that, as it happens in the hydrodynamic approach, the limits B → 0 and ω → 0 do not commute). We find Let us rearrange the system of differential equations involving T and δµ I , Eq.(A.7) and Eq.(B.9), in a more practical form so to be able to discuss the relaxation length scale of the two unknown functions in the x-direction.

σ B = N 3d V e c B ≈ 2
We write the prefactor of (-∂ xτ ) as k ≡ u k u where k u is given by:

k u = 1 + σ 1 + σ 2 σ 2 -σ 1 α yx κ xx Q, (D.1)
and u is given by:

u = e k B n L 2 Q σ 1 -σ 2 σ 1 + σ 2 w 4 L 1 σ 1 - 1 σ 2 κ xx (D.2)
Eq.s(IV.13) follow: The Fermi velocity v F = c/300 = 10 8 cm/sec. We take the parameter that describes the leakage from edge to bulk g ∼ 10 -6 , so that ga -2 ∼ 10 0 . For β ∼ 10 -9 and a ∼ 10 -3 Eq.(D.1) and Eq.(D.2) give ( see in the next, Subsection D) k ≡ βk a -1 = q 2 + 3 × 10 -3 q u ≡ βu a -1 = q 2 , q = σ 1 -σ 2 [σ 2 1 + σ 2 2 ] 1/2 . (D.8)

The order of magnitude of θ is fixed by the applied voltage (expressed in V olts ≡ V ), as derived from the experiment:

Ew ∼ 10 -4 V → R ∼ 10 -4 V 10 n Amp ∼ 10 kΩ

→ θ ≡ θ a -1 = α xy κ xx Ew a -1 L w ≈ 3 × 10 -2 . (D.9)
The final equation is:

∂ 2 xµ + [1 -(k -u ) -θ ] ∂ xµ + (g + θ ) µ -g θ µ d x = 0 (D.10)
Appendix E: Properties of the imbalance charge energy solution

Having discussed the order of magnitudes and the general features of the solutions of the system Eq.(D.3), we now turn to a numerical solution to highlight the dependence of δµ I (x) and τ (x) on the magnetic field B, both far from x = 0 and close to the origin.

We have shown in Appendix D that the system of Eq.(D.3) allows for a rewriting with parameters and unknown functions of O(1). The system of Eq.(D.6) has been used in the numerical results which are discussed in Section V of the text. with g = ga -2 and primed parameters given in Eq.(D.8, D.9). As discussed in the text, the solution of the homogeneous equation Eq.(E.1, (b)) has to be added to the solution of the system of coupled equations.

Here we keep the form of the system given in Eq.(E.1), but we drop all the primes for simplicity and we use ∂ instead of ∂ x.

For g, θ = 0, a numerical solution is obtained by deriving an equation for µ only and, by inserting the solution into Eq(D.3,(a)). Multiplying Eq.(D.3,(a)) by k and substituting it in Eq.(D.3,(b)) , we get:

(k -u) ∂ 2 (θ τ ) = ∂ 2 µ + [1 -(k -u)] ∂ µ -g µ (E.2)
Let us consider ∂ (θ τ ) = θ ∂ τ + τ ∂ θ. The second term on the r.h.s. is reminiscent of ∂ x E y -∂ ỹ E x where E x = 0 because, in the bulk, the circuit is open ( no current along x ) and assumed to be equipotential. This z-component of ∇ × E vanishes, because there is no current in the z direction. We conclude that only the first term survives with a θ which is practically constant in x, in an interval of size θ -1 L close to the origin, or zero for any x far from the origin. Following Eq.(E.1, E.2), we define

ν(x) = ν 0 + x 0 µ(x ) dx (E.3)
and we get:

   ∂ 2 τ -θ ∂ τ = µ, (k -u) θ ∂ τ = ∂ µ + [1 -(k -u)]µ -g ν ∂ ν = µ
or, substituting the second into the first:

   (k -u) ∂ 2 τ = ∂ µ + µ -g ν, (k -u) θ ∂ τ = ∂ µ + [1 -(k -u)] µ -g ν ∂ ν = µ
We now multiply the first equation by θ and differentiate the second one. Equating the two of them to each other, we obtain:

∂ 2
x µ + (a -k β) ∂ xµ -g µ = θ ∂ x µ + θ a -u β µ -g θ ν.

When µ has been found, one can plug it into the first one of Eq.s(E.4) to find ∂ x τ . The equation can be written as a first order differential system:

   ∂ x µ = ξ ∂ x ν = µ ∂ x ξ = -[1 -(k -u) -θ] ξ + [g + θ] µ -g θ ν,
or, in the matrix form, For θ = 0, solutions converging at infinity require g, θ < 0.

The vanishing of the imbalance parameter q implies the vanishing of (k -u). Here we show that, for θ = 0, the vanishing of q, provides two real or complex conjugate solutions for λ which are independent of θ and a real solution z 3 = -θ sign[g].

If z 1 , z 2 , z 3 are the roots of the cubic equation the following relations hold:

z 1 z 2 + z 1 z 3 + z 2 z 3 = -G, z 1 + z 2 + z 3 = -F, z 1 z 2 z 3 = -gθ. (E.6)
Given z 1,2 = -a/2 ± (1 -4|g|/2, it is easily found that z 3 = -θ sign[g] satisfy Eq.s(E.6). When q = 0, one solution will be real and the other two will be complex or all of them will be real.

Appendix F: The non local voltage

The equation for the non local voltage, in dimensionless units v y = e (c +,A -c +,B )/(k B o K) is derived as follows.

We go back to Eq.(B.5): where K ± is defined in Eq.(A.4). With Eq.s(A.1, A.2), by keeping only first order in g, we get:

∂ 2 x (c +,A -c +,B ) -g (2γ + g )(c +,A -c +,B ) = g (2γ + g ) K + -∂ x T + w 4 α xy E Q 1 σ 1 + 1 σ 2 T -g ∂ x K -(-∂ x ∆T ) + w 4 α xy E Q 1 σ 1 - 1 σ 2 ∆T
According to Eq.(B.8), we have chosen to impose the vanishing of the second square bracket on the right hand side, so that the final equation for the Hall voltage of Eq.(V.1), reads

FIG. 1 :

 1 FIG. 1: (a) Longitudinal (black line) and transverse (red line) magnetoresitances Rxx and Rxy as measured on a graphene Hall bar close to the CNP at T = 1.7 K and I = 10 nA. (b) Non local magnetoresistances R28,37 (close to the current source, blue line) and R28,46 (far from the current source, green line). The inset is a sketch of the Hall bar which labels the contacts. (c) Temperature dependence of the nonlocal resistances R28,37 and R28,46 at a fixed magnetic field B = 3 T.

FIG. 2 :

 2 FIG. 2: Sketch of the nonlocal experiment. A charge flow is applied across the Hall bar on its left part (yellow colour). The magnetic field induces space separation of the charges of opposite sign, preferentially leaking from the edges (blue and read circles) in the bulk of the Hall bar. Thermal diffusion along the Hall bar drives the charge imbalance of the chemical potential between electrons (blue line) and holes (red line) away from the source, which relaxes poorly over microscopic distances. The thermal and Lorentz force acting on the imbalance flow are compensated by a transverse electric field which can be detected by nonlocal probes. Flipping the direction of the magnetic field requires exchanging the read and blue circles and lines with each other.

FIG. 3 :

 3 FIG. 3: (a) Transport in the Hall bar geometry. The edge states are denoted by red and blue lines and carry opposite charge carriers. The various couplings between the edge states and the bulk are also indicated. The notation used in the text for the edge chemical potentials ci (i = 1, 2) and bulk chemical potentials si is reported.

FIG. 4 :

 4 FIG. 4: Results of the Abanin model, reproduced from Ref.[34]. Density dependence of the transport coefficient ρxx, ρxy, Gxx= ρxx/(ρ 2 xx + ρ 2 xy ), Gyx = ρxy/(ρ 2 xx + ρ 2 xy ). The conductivities are modeled by Gaussians centered at ν ± 1, σ (1,2) xx = exp(-A(ν ± 1) 2 ) and σ (1,2) xy are given by the semicircle relation. A = 5, gw = 1 and γw = 6 are used.

FIG. 5 :

 5 FIG. 5: Normalized [e× non local voltage] vs b = B/(2T esla), at various distances from the origin: D/a -1 L = 2., 4., 6., 8.

FIG. 7 :

 7 FIG. 7: Normalized imbalance chemical potential δµI vs distance D/a -1 L from the origin, at different magnetic fields ( b = B/(2 T esla) = 2., 4., 6., 8.)

FIG. 9 :

 9 FIG. 9: Normalized imbalance chemical potential δµI vs distance D/a -1 L from the origin, at different magnetic fields (b = B/(2 T esla) = -2., -4., -6., -8.)

FIG. 11 :

 11 FIG. 11: Normalized imbalance chemical potential δµI vs b = B/(2 T esla) at various distances D/a -1 L = 2., 4., 6., 8. from the origin

≈ 4 .

 4 84 × 10 12 sec -1 (C.5) With these conductances, the ratio of Eq.(C.4) is estimated as:× 10 -2 V (C.6)The Nernst coefficient in the limit of large fields and relatively clean samples 50 is e 3d N = (σ -1 α) xy ≈ 1 mV / o K.Appendix D: Coupled differential equation and relaxational scales

∂ 2 x 2 x

 22 τ -∂ x (θ τ ) = β δµ I , (a) ∂ 2 x δµ I + a ∂ x δµ I -g δµ I = ∂ 2 x {k ∂ x τ -θ u τ } . (b) (D.3) A numerical solution is obtained by deriving an equation for δµ I only and by inserting the solution into Eq(D.3,(a)). Multiplying Eq.(D.3,(a)) by k and substituting it in Eq.(D.3,(b)) , we get:(k -u) ∂ 2 x (θ τ ) = ∂ 2 x δµ I + (a -k β) ∂ x δµ I -g δµ I (D.4)From now on we shift a → a -u β by summing and subtracting u β to a, and ridefine it as a thus obtaining:(k -u) ∂ 2 x (θ τ ) = ∂ 2 x δµ I + [a -(k -u) β] ∂ x δµ I -g δµ I . (D.5)We cast the system in length units a -1 L: ( x = x/(a -1 L) and introduce the new variables µ ≡ βa -1 δµ I , τ ≡ a τ to get:∂ 2 x τ -∂ x θa -1 τ = µ , (a) ∂ 2 x µ + ∂ x µ -ga -2 µ (D.6) = ∂ ( βka -1 ) ∂ x τ -θa -1 ( βua -1 ) τ . (b)where a ≡ 2 w (σ 1 + σ 2 )V v F , g ≡ g(2γ + g) L 2 t 2 0 . (D.7)

∂ 2 x

 2 τ -∂ x (θ τ ) = µ , (a) ∂ 2 x µ + ∂ x µ -g µ (E.1) = ∂ 2 x {k ∂ x τ -θ u τ } . (b)

F

  = 1-(k-u)-θ and G = g+θ. As discussed above, a → aa -1 ∼ 1, k β → k βa -2 , g → ga -2 , u β → u βa -2and θ are all O(1). The characteristic equation giving the space dependence of the functions here e λx is λ3 + F λ2 -G λ + g θ = 0 (E.5)

∂ 2 x

 2 (c +,A -c +,B ) = g (2γ + g ) [c +,A -c +,B -2s + ] + g ∂ x (s -,A -s -,B ).Inserting Eq.(A.1) and Eq.(A.3) in it, we get ( T ≈ 1 o K):∂ 2 x (c +,A -c +,B ) -g (2γ + g )(c +,A -c +,B ) = g -K --∂

-κ xx - σ 1 -σ 2 σ 1

 11 + σ 2 α xy Q ∂ xτ .(F.1)with x = x/L, τ = T / o K, to O[(σ 1 -σ 2 ) 2 ].

  This is fully consistent with the initial definition δµ I = e ρ I L 2 /k B

				Comparison with Eq.(IV.6) requires that
				ρ I e	= n	δµ I e	k B e	o K = 2.3 ×	10 10 cm 2	86.17 × 10 -6 27.16 × 0.529 Å
								= 14 × 10 4 cm 2 Å -1 .	(V.8)
				5)
	consistent with the definition of β.
	With E ∼ V /m we get
	θ =	α xy κ xx	EL ∼ 10 -5	(V.6)
	and
	δµ I ∼ η θ 2 ∼ 2.3 .	(V.7)

o K, provided we set L ∼ 0.1 µm:

  3-d strip. This value for N 3d allows us to match reasonably well the values from a Fermi Liquid approach of Eq.(C.1) and those experimentally found for α xy , κ xx and e N . The value reported by experimentalists for κ xx ,adopted by us, is κ xx ≈ 1÷10W att/( o K ×meter) .

	The ratio κ F L xx /α F L xy is given by:
	κ F L xx α F L xy	=	+ P n e 3d	σ 3d σ B	1 -	σ B σ 3d	2	,	(C.4)

Acknowledgments

We thank P. Brouwer, A. De Martino and R. Egger for fruitful discussions. This work was founded by project PLASMJAC Univ. Napoli "Federico II".

Appendix A: Diffusion equation for the average temperature T Here we derive the diffusion equation for the temperature T (x), induced by the imbalance charge ρ I (x) defined in Eq.(IV.6), ρ I = n (c +,A -c +,B -2 s + ).

Using the notation introduced in Section III, We assume the mirror symmetry between the side A and B of the bar, so that s + = s +,A = -s +,B in the following. Solving Eq.s(IV.5) with respect to the potentials "s", we get (Q

Appendix B: Equation for the imbalance energy δµI

We now derive the equation that relates the imbalance energy δµ I to the average temperature T . From Eq.s(IV.1), we get:

The scattering time is the e -e inelastic one, because this is the time scale which qualifies the imbalance relaxation rate: τ -1 ee = α 2 f k B T / . The elastic scattering is unable to relax the charge imbalance. At B = 2 T esla, one gets:

with k F = (π n ) 1/2 and n = 10 10 /cm 2 , v F = c/300, B = 2 T esla.

In terms of the enthalpy, the thermal response functions in the limit ω c τ >> 1 are given by Eq.s(II.6,II.8). In the Fermi liquid regime, in which the thermal energy