Giant magnetoresistance of Dirac plasma in high-mobility graphene

The most recognizable feature of graphene’s electronic spectrum is its Dirac point, around which interesting phenomena tend to cluster. At low temperatures, the intrinsic behaviour in this regime is often obscured by charge inhomogeneity1,2 but thermal excitations can overcome the disorder at elevated temperatures and create an electron–hole plasma of Dirac fermions. The Dirac plasma has been found to exhibit unusual properties, including quantum-critical scattering3–5 and hydrodynamic flow6–8. However, little is known about the plasma’s behaviour in magnetic fields. Here we report magnetotransport in this quantum-critical regime. In low fields, the plasma exhibits giant parabolic magnetoresistivity reaching more than 100 per cent in a magnetic field of 0.1 tesla at room temperature. This is orders-of-magnitude higher than magnetoresistivity found in any other system at such temperatures. We show that this behaviour is unique to monolayer graphene, being underpinned by its massless spectrum and ultrahigh mobility, despite frequent (Planckian limit) scattering3–5,9–14. With the onset of Landau quantization in a magnetic field of a few tesla, where the electron–hole plasma resides entirely on the zeroth Landau level, giant linear magnetoresistivity emerges. It is nearly independent of temperature and can be suppressed by proximity screening15, indicating a many-body origin. Clear parallels with magnetotransport in strange metals12–14 and so-called quantum linear magnetoresistance predicted for Weyl metals16 offer an interesting opportunity to further explore relevant physics using this well defined quantum-critical two-dimensional system.

A variety of mechanisms-both intrinsic and extrinsic-can lead to large magnetoresistance (MR) in metallic systems. The quest to understand these mechanisms has continued for longer than a century but many gaps still remain, which is especially obvious for the MR reported in newcomer materials such as various Dirac and Weyl systems 17-25 , strange metals [12][13][14] and so on. The history and current status of the research field are briefly reviewed in Methods. Whichever mechanism is behind a particular MR behaviour, it always relies on bending of electron trajectories by a magnetic field B and, accordingly, high carrier mobility μ is an essential prerequisite for the observation of large MR. Colossal MR (reaching about 10 6 % in a magnetic field of 10 T) was observed in a number of high-μ systems at liquid-helium temperatures 17-25 . However, because mobility decreases with increasing temperature T, this usually results in only a tiny MR above liquid-nitrogen temperatures. Those few materials in which carriers remain highly mobile at room temperature (such as doped graphene and indium antimonide) [26][27][28] are all non-compensated systems and, in agreement with the classical theory of normal metals 29 , their longitudinal resistivity ρ xx saturates in high B, leading again to little MR. Only the presence of extended defects 30-32 or a special design of four-probe devices 26,33 that creates strongly non-uniform current flows can lead to large-but extrinsic-MR (Methods).
As shown below, thermally excited charge carriers in monolayer graphene (MLG) at the neutrality point (NP) exhibit an anomalously high μ exceeding 100,000 cm 2 V −1 s −1 at room temperature, despite the fact that the system is strongly interacting [3][4][5][6][7][8] and the electronhole (e-h) scattering time τ P is ultimately short, being limited by the uncertainty principle τ P −1 ≈ Ck B T/h where k B and h are the Boltzmann and Planck constants, respectively, and C ≈ 1 is the interaction constant [3][4][5][9][10][11][12] . Importantly, unlike any known system with high μ at room temperature, the Dirac plasma is compensated (charge neutral) so that its zero Hall response allows non-saturating MR 29 whereas the high μ makes it colossal. To emphasize how unique magnetoresistivity ρ xx (B) of the Dirac plasma is, we provide its comparison with graphite (multilayer graphene 34 ) and charge-neutral bilayer graphene (BLG), another quantum-critical system exhibiting Planckian scattering but having massive charge carriers with modest mobilities 9,10 .

Giant MR in non-quantizing fields
Our primary devices were multiterminal Hall bars made from MLG encapsulated in hexagonal boron nitride (hBN; Fig. 1a). We have studied several such devices and focus here on two of them (devices D1 and D2) showing representative behaviour. At low T, their mobilities exceed 10 6 cm 2 V −1 s −1 at characteristic carrier densities of about 10 11 cm −2 , being limited by edge scattering despite the devices' size being more than 10 μm. The typical behaviour of ρ xx as a function of the gate-induced density n is shown in Fig. 1b. If the same curves are replotted on a log scale (Fig. 1c), it becomes clear that ρ xx responds to gate voltage only above a certain n dependent on T. This behaviour is commonly quantified as shown in Fig. 1c,d where δn marks the gate-induced density that leads to notable changes in ρ xx . At high T, the peak in ρ xx (n) broadens because of thermally excited electrons and holes in concentrations n th = (2π 3 /3)(k B T/hv F ) 2 , where v F is the Fermi velocity (Methods). The extracted δn evolves proportionally to T 2 as expected (Fig. 1d) and its absolute value is about 0.5n th , which means that to make changes in ρ xx visible on such log plots, gate-induced carriers are required in concentrations of about 50% of the thermal concentration. At low T, δn saturates typically at about 10 10 cm −2 because of residual charge inhomogeneity (e-h puddles of submicrometre scale) 1,2 . Below we focus on T > 100 K where thermal excitations totally overwhelm the residual δn.
The Dirac plasma's response to small fields is shown in Fig. 1e. One can see that the longitudinal resistivity at the NP ρ NP ≡ ρ xx (n = 0) increases proportionally to B 2 , as expected from the classical Drude model 29 . However, the changes in ρ NP are unexpectedly large for this T range. Indeed, if we consider 0.1 T as a characteristic field relevant for magnetic-sensor applications, then the relative magnetoresistivity Δ = [ρ xx (B) - ρ xx (0)]/ρ xx (0) reaches about 110% at 300 K near the NP (Fig. 1e) and increases by a factor of 3-4 at 200 K. For comparison, Δ in normal metals rarely exceeds a small fraction of 1% above liquid-nitrogen temperatures. Even high-quality encapsulated bilayer, few-layer and multilayer graphene exhibit Δ(0.1 T) reaching only about 1% at room temperature (Methods). Also, the renowned giant MR based on spin flipping in ferromagnetic multilayers yields changes in resistance that are one to two orders of magnitude smaller 35,36 than those observed here for the Dirac plasma.
Further characterization of the e-h plasma is provided in Fig. 2. It shows that Δ rapidly diminishes away from the NP at characteristic densities n ≈ n th (Fig. 2a). This is expected 29 because, for non-compensated systems, changes in ρ xx (B) should be small and saturate, if Hall resistivity ρ xy > ρ xx (Methods). In contrast, for charge-neutral systems (zero ρ xy ), the Drude model predicts non-saturating magnetoresistivity such that Δ = μ B 2 B 2 where μ B is the mobility in non-quantizing magnetic fields (Methods). The latter expression describes well the behaviour observed in small B (Fig. 2b). Figure 2c shows the extracted μ B as a function of T. The mobility exceeds 100,000 cm 2 V −1 s −1 at room temperature and grows above 300,000 cm 2 V −1 s −1 below 150 K. Although high μ values are well known for the Fermi-liquid regime in doped graphene, it is unexpected that the mobility remains high in the presence of Planckian scattering, characteristic of the quantum-critical regime in neutral graphene 5,6 . For comparison, bilayer and multilayer graphene also exhibit very high mobilities at liquid-helium temperature, but their ρ NP (B) are practically flat at elevated T (Fig. 1e), yielding μ B of only about 10,000 cm 2 V −1 s −1 at 300 K (Extended Data Figs. 2 and 3). The marked difference in electronic quality between the e-h plasmas of relativistic and non-relativistic fermions (in MLG and BLG, respectively) stems from the small effective mass m characteristic of the Dirac spectrum (μ ∝ m −1 ) and its low density of states, which reduces the efficiency of electron scattering (Methods). It is noted, however, that the Dirac spectrum on its own is insufficient for achieving giant values of Δ, and the high quality of MLG devices is paramount. This is emphasized by Extended Data Fig. 8, which shows magnetotransport for non-encapsulated graphene on a silicon oxide substrate. Such low-quality MLG exhibits three orders of magnitude smaller MR.
It is instructive to compare the found μ B with the zero-field mobility μ 0 . The latter can be evaluated using the standard Drude formula ρ NP −1 = 2n th eμ 0 , where e is the electron charge and the factor 2 accounts for equal concentrations of electrons and holes at the NP. Figure 2d

Article
shows that ρ NP quickly decreases with increasing T from liquid-helium temperatures to about 100 K but, as the Dirac plasma gets established (n th >> residual δn), ρ NP becomes almost T independent with a constant value of about 1 kΩ above 150 K (also, inset of Fig. 3b and Extended Data Fig. 1). The saturating behaviour of ρ NP is attributed to the onset of the quantum-critical regime in which the scattering is dominated by the Planckian frequency, τ P −1 . Indeed, ρ NP ≈ 1 kΩ yields C ≈ 0.7 close to unity, as expected [3][4][5][9][10][11][12] . This analysis also agrees with that of the quantum-critical behaviour reported for BLG 9,10 (Methods) and conclusions about MLG from other measurements 5 . Figure 2c shows that μ 0 evolves proportionally to 1/T 2 , as expected for Planckian systems with a Dirac spectrum (Methods). Surprisingly, μ 0 is two to three times smaller than μ B . As shown in Methods, this happens because μ B is less sensitive than μ 0 to the dominating e-h scattering. Qualitatively, the difference can be understood as arising from different relative motions of electrons and holes in zero and finite B. In zero B, the electric field forces electrons and holes to move in opposite directions so that e-h collisions are efficient in impeding a current flow. In contrast, cyclotron motion causes a drift of both electrons and holes in the same direction. Therefore, e-h collisions do not affect the Hall currents responsible for magnetoresistivity. This explanation is further substantiated by our measurements using screened graphene devices (encapsulated MLG with metallic gates placed at a distance of about 1-3 nm) 15 . The screening is found to suppress Coulomb scattering, which results in smaller C and, therefore, higher μ 0 (Extended Data Fig. 4a). However, the same screening has little effect on ρ NP (B) and hence μ B (Extended Data Fig. 4b), in agreement with theory. This consideration is equally applicable for e-h plasma of massive fermions and, indeed, a similar difference between μ 0 and μ B is observed for neutral BLG (Extended Data Fig. 2). The above analysis allows us to conclude that the anomalously large MR in low B arises owing to ultrahigh mobility of Dirac fermions, combined with ineffectiveness of e-h scattering in suppressing Hall currents.

Strange linear MR in the extreme quantum limit
In high B, magnetotransport in the Dirac plasma exhibits profound changes such that, above a few tesla, ρ NP (B) evolves from being parabolic to linear ( Fig. 3 and Extended Data Fig. 5). Slopes of this linear MR are found to be similar for all the studied devices (Fig. 3b, inset) and almost independent of T. The crossover between parabolic and linear dependences is marked by a flattened section on the curves which appears at T below 200 K. We attribute the flattening to the onset of Landau quantization (Extended Data Fig. 6). This attribution also agrees with the fact that at B ≈ 3-5 T, the main cyclotron gap between the zeroth and first Landau levels (LLs) reaches about 800 K, notably exceeding the thermal smearing k B T. As for the MR magnitude, Δ reaches about 10 4 % at 10 T and, despite the linear (slower than parabolic) dependence in quantizing fields, this is again among the highest for room-temperature experiments 32 . Comparison with multilayer and low-quality graphene   (Extended Data Figs. 3 and 8) shows the importance of both the Dirac spectrum and the electronic quality for such a giant MR response. Another notable feature of magnetotransport in the plasma of the zeroth LL is that ρ NP at a given B increases with increasing T (Fig. 3a and Extended Data Fig. 5). This contradicts the orthodox MR behaviour observed in all other systems, which results in lower Δ at higher T because of increased scattering 29 (Methods). To shed light on strange magnetotransport, we have also tested how ρ NP (B) is affected by proximity screening. Although the parabolic dependence of ρ NP in low B was practically unaffected (as discussed above), the screening greatly suppressed MR in quantizing B (Fig. 3b). The linear slopes of ρ NP (B) decrease from 5-8 kΩ T −1 in our primary devices to 1-3 kΩ T −1 in those with screening (inset of Fig. 3b), implying that magnetotransport in the zeroth LL depends on Coulomb interactions.
In discussing the high-B behaviour, we first note that the previously reported linear MR can in most cases be attributed to complex current flows that become increasingly non-uniform as ρ xy ∝ B increases (Methods). The involved mechanisms are based on either the spatial inhomogeneity or the presence of edges. To check for possible edge effects in our case, we have studied Corbino disks fabricated from encapsulated MLG and found very similar ρ NP (B) (Extended Data Fig. 7). Thus, for our zeroth-LL plasma with zero ρ xy , those extrinsic mechanisms can be ruled out (Methods). It may also be tempting to evoke Abrikosov's linear MR 16 predicted to occur in three-dimensional (3D) semimetals with Dirac-like spectra in the extreme quantum limit. However, the Born approximation used in the 3D model cannot be justified for two-dimensional (2D) transport in a smooth background potential 37 because charge carriers remain localized within electron and hole puddles.
For the lack of a theory suitable to describe the observed linear MR, we employ the simple Drude model by considering cyclotron-orbit centres as quasiparticles that circle along equipotential contours and, also, diffuse between them owing to electron scattering. The density of such quasiparticles is determined by the capacity of the LL, n LL = 2B/ϕ 0 >> n th , where ϕ 0 = h/e. For charge neutrality, the Drude model yields ρ NP (B) ≈ ρμ 2 B 2 (Methods), where n and μ in the standard expression ρ = 1/neμ should be substituted with n LL and μ Q , respectively, to reflect the density and mobility for the plasma of the zeroth LL. This leads to The linearity in B arises from the fact that the B 2 dependence inherent for compensated semimetals is moderated by the linear increase in the carrier density in the zeroth LL. Next, to estimate μ Q , we assume that Planckian scattering moves quasiparticles by a typical distance ℓ between equipotentials, resulting in the diffusion coefficient

Diffusion within individual puddles leaves carriers localized inside.
Only if a quasiparticle covers a distance of approximately ξ between neighbouring puddles, those processes contribute to macroscopic currents along the electric field and, hence, global conductivity. Accordingly, the timescale relevant for electron transport in the zeroth LL is given by τ tr ≈ ξ 2 /D >> τ p and the corresponding diffusion coefficient can be written as D tr = v T 2 τ tr ≈ v T 2 ξ 2 /D = ξ 2 /τ p . Then, using the Einstein-Smoluchowski equation, we find the transport mobility μ Q = eD tr /k B T ≈ eξ 2 /k B Tτ p = Ceξ 2 /h, where both v T and ℓ fell out from the final expression. This result suggests that the MR of the zero LL is linear in B and independent of T, as observed experimentally, and may also explain the suppression by proximity screening as smaller C result in lower μ Q . Furthermore, equation (1) can be rewritten as where ℓ B is the magnetic length. Equation (2) closely resembles the result of a formal extension of Abrikosov's model into the 2D case 37 . Although the above consideration catches the main physics and qualitatively agrees with our observations, further work is required to develop a microscopic theory of magnetotransport in the 2D Boltzmann plasma at the zeroth LL.

Outlook
The Dirac plasma in graphene exhibits the one of the highest MRs observed above liquid-nitrogen temperatures in both low and high fields. In low B, only ferromagnetic devices employing spin tunnelling 38 or the use of four-probe geometry 26,33 allow stronger electronic response to magnetic fields. In contrast to the latter phenomena, the giant MR of graphene stems from its magnetoresistivity ρ xx (B).
In quantizing fields, graphene experiences a system transformation, becoming an e-h plasma residing in the zeroth LL. Our observations are also relevant to the physics of strange metals that exhibit Planckian scattering. Strange metals display the renowned linear T dependence of their resistivity, in obvious contrast to our case. However, this difference arises only because strange metals have a fixed carrier density whereas the carrier density and effective mass in the Dirac plasma increase with T, leading to the constant ρ NP . Moreover, strange metals also exhibit linear MR that is weakly T dependent. This MR remains unexplained, although a recent ansatz 13,14 suggests that in Planckian systems, τ −1 should be defined by the largest relevant energy scale, either k B T or some magnetic-field-induced energy proportional to B.
The ansatz does not seem work for the Dirac plasma because the only relevant and sufficiently large magnetic energy is the cyclotron gap. It evolves as B 1/2 rather than linearly in B. Notwithstanding any differences, Planckian systems in high fields remain poorly studied, and graphene's Dirac plasma offers a model system to understand the relevant physics. The possibility to modify magnetotransport by tuning electron-electron interactions using proximity screening is especially appealing in this context.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-023-05807-0. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Brief history of linear MR
Studies of the electrical response of metals to magnetic fields go back to experiments by Lord Kelvin and Edwin Hall over one-and-a-half centuries ago 39,40 . Although the subject continued to attract sporadic attention during the following decades (see, for example, ref. 41 ), the first systematic study of MR phenomena is usually credited to Pyotr Kapitsa. In 1928-1929, he reported high-field studies of MR in 37 different materials 42,43 . This research brought up two major findings. First, some materials (for example, bismuth, arsenic, antimony and graphite) were found to exhibit MR exceeding 100% in a magnetic field of 30 T at room temperature, much higher than the others in that study. So large MR could not be explained by contemporary theories. Second, despite different absolute values of MR, all the studied materials followed a universal B dependence. In small fields, it was always parabolic, in agreement with the already accepted understanding that cyclotron motion of current-carrying electrons should bend their trajectories and, hence, increase resistivity. However, in fields above several tesla, MR was found to increase linearly, which was unexpected.
The first puzzle of large MR values was solved relatively quickly, owing to the development of the band theory. Most of the materials exhibiting large room-temperature MR in Kapitsa's experiments appeared to be semimetals so that the electric current was carried by both electrons and holes. It is now well understood that the reduced Hall effect in this case leads to ρ xx evolving in high fields approximately as 1/σ xx , where σ xx is the longitudinal conductivity. This is in contrast to the case of one type of charge carriers where ρ xx ≈ σ xx ρ xy 2 ('Drude model for charge-neutral graphene' below). The second puzzle of linear MR has attracted numerous theories and explanations. In general, there are several mechanisms that can cause linear MR and, even today, its observation often leads to controversies because it is difficult to pinpoint the exact origin.
One of the first mechanisms causing linear MR was proposed by Lifshitz and Peschanskii 44 . In 1959, they considered magnetotransport in polycrystalline metals with open Fermi surfaces. For certain orientations of the magnetic field with respect to crystallographic axes, such metals flaunt open cyclotron orbits that result in non-saturating MR proportional to B 2 (refs. 45,46 ). However, this quadratic behaviour occurs within only a narrow interval of angles, which decreases proportionally to B −1 . For the other angles, cyclotron orbits remain closed, and MR attributable to them saturates in high B. Averaging over all angles for polycrystalline samples resulted in linear MR, and this result helped to explain many-but not all-observations in the literature. Those ideas were further developed by Dreizin and Dykhne 47 who obtained MR proportional to B 4/3 and MR proportional to B 2/3 , depending on whether a metal with an open Fermi surface was compensated or not, respectively. Moreover, the authors presented a magnetotransport theory for not only polycrystalline but also inhomogeneous conducting media. Depending on the Fermi surface and compensation between charge carriers, various powers of B could be obtained including, for example, linear MR in compensated semimetals with 2D disorder 47 .
The MR theory relying on materials' inhomogeneity was expanded both theoretically and experimentally in the 1970s and 1980s. It was shown that macroscopic strain 48 , voids [49][50][51][52] and thickness variations 53,54 could lead to linear MR in high B (μB >> 1). The next step was taken in 2003 by Parish and Littlewood who considered the case of very strong inhomogeneity that could not be described by the earlier theories 55 . Using a random 2D resistance network, they obtained linear MR that starts from small magnetic fields (μB < 1) and could explain the behaviour observed in some disordered semiconductors 55 . The fundamental reason for MR in all the cases involving inhomogeneous media is the following. In the presence of regions with different magnetotransport coefficients, the arising Hall voltages (large for μB >> 1) necessitate substantial changes in the electric current distribution to satisfy boundary conditions at interfaces between different regions. As a result, the electric current becomes increasingly inhomogeneous, being squeezed into narrow streams near the interfaces. This current inhomogeneity increases the effective resistance of the medium 53,54 .
A different mechanism was suggested by Abrikosov 16, 56,57 . He pointed out that some materials exhibiting linear MR were neither polycrystalline nor inhomogeneous but single crystals with closed Fermi surfaces including graphite, bismuth and other materials 58,59 . To explain these observations, Abrikosov considered a Weyl (3D Dirac-like) spectrum so that, in quantizing B, all charge carriers collapsed onto the lowest (zeroth) LL. Assuming a scattering potential caused by screened charged impurities, linear MR was predicted in this case. Because of the essential role played by Landau quantization, the effect was called quantum linear MR 16,56,57 . The Abrikosov mechanism attracted considerable interest and was invoked as an explanation for many experiments 60,61 , even though the concerned materials often poorly matched the assumptions required by the theory (including being 2D rather than 3D systems). Unfortunately, Abrikosov provided no explanation for the physics behind his theory and only recently 37 it has been shown that his analysis is equivalent to calculations of diffusion of cyclotron-orbit centres in an electrostatic potential. This conceptual overlap requires mentioning of the earlier theories by Kubo and Ando for diffusion of cyclotron-orbit centres 62,63 . Furthermore, within the self-consistent Bohr approximation, linear MR was shown to appear in the 2D case for strongly screened charged impurities whereas, for short-range scattering, MR becomes sublinear 64 . The formal extension of Abrikosov's theory into two dimensions also leads to linear MR 37 .
Finally, two other mechanisms that result in linear MR have to be mentioned. First, Alekseev with colleagues showed that e-h annihilation at the edges of 2D semimetals could lead to linear MR 65,66 . This mechanism can be ruled in or out by comparing magnetotransport in Hall-bar and Corbino-disk devices, as done in our work. Second, so-called strange metals often exhibit resistivity that increases linearly not only with temperature but also with magnetic field 13,14 . Although such linear MR does not follow from the so-called holographic approach 11 , it was suggested 13,14,67 that the quantum-critical scattering rate τ −1 ≈ E/h could be controlled by the maximum relevant energy E in the uncertainty equation, that is, by either k B T or μ Bohr B where μ Bohr is the Bohr magneton. This could then explain both (T and B) linear dependences in strange metals. It is also worth mentioning that linear MR was recently reported in two other 2D strongly interacting systems, namely, twisted tungsten diselenide 68 and magic-angle graphene 69 . It was suggested that the MR had the same origin as in strange metals.

Earlier studies of MR in graphene and Dirac-type materials
Over the past decade, there have been numerous studies of magnetotransport using newly available materials such as graphene (see, for example, refs. 30-34,61,70-74 ), topological insulators (see, for example, refs. [75][76][77] ) and high-mobility Dirac and Weyl semimetals (see, for example, refs. 17-25 ). Graphene has attracted particular attention as a promising material for magnetic-field sensors owing to its high μ at room temperature. The first generation of graphene devices (graphene placed on oxidized silicon and so-called epitaxial graphene) exhibited relatively low μ, and their MR was also relatively low, reaching only about 100% in fields above 10 T (refs. 30-32,61,70-72 ; 'MR of low-mobility graphene' below). The MR typically originated from charge inhomogeneity and other disorder, although some reports suggested 61 the observation of Abrikosov's linear MR in doped multilayer graphene. Later research ruled out this explanation, arguing that the observed linear MR originated from a polycrystalline disorder 31, 55 .
The next generation of graphene devices using encapsulation with hBN exhibited exceptional electronic quality 2,78 . So far, magnetotransport in graphene-on-hBN devices has been studied at elevated T only for few-layer graphene 34 and MLG away from the NP 33 . Fewand multilayer graphene (graphite) exhibit a relatively low μ at elevated temperatures. This results in small quadratic MR in low B and also limits MR in high fields 32,34 , in agreement with our results in 'MR in multilayer graphene' below. Doped high-μ graphene exhibits saturating magnetoresistivity and its magnitude is small, as expected. It is noted, however, that if one uses a geometry that instigates a non-uniform current flow, it is possible to enhance the apparent MR in four-probe measurements, for example, using the so-called extraordinary MR configuration 26,33 . In such a geometry, the central part of a MLG device is replaced with a highly conducting metal (for example, gold). In zero B, the current mainly flows through the metal, despite being injected into graphene. The magnetic field curves the current trajectories and forces charge carriers to move through graphene, which is much more resistive than gold films. Accordingly, the apparent four-probe MR could reach extraordinary values of about 10 7 % at 9 T and room temperature. This is comparable to typical changes in Hall voltage that also require a four-probe geometry. It is noted that this extraordinary MR is not an intrinsic property of a material and, accordingly, translates into only modest changes for any two-probe measurements. Until now, no studies of magnetotransport at elevated T have been reported for charge-neutral MLG with high μ.
For completeness, let us mention extensive magnetotransport studies of 3D counterparts of graphene, which are different topological insulators, Dirac and Weyl semimetals, and other clean semimetals such as tungsten telluride (also suggested to be a Weyl semimetal 79 ). Many of them showed huge MR, which in some cases exceeded 10 6 % at liquid-helium temperatures [17][18][19][20][21][22]25 . Such colossal values were attributed to the high mobility of charge carriers in these materials (μ reaching above 10 6 cm 2 V −1 s −1 at 4 K, similar to encapsulated graphene). However, the mobility rapidly decayed with increasing T, which resulted only in a tiny low-B MR at elevated T. This is not the case for MLG that exhibits high μ at room temperature even at the NP, which results in the colossal quadratic MR in low B, as reported in this work.

Drude model for charge-neutral graphene
To evaluate the magnetotransport properties of our devices, we have used the standard two-carrier model for electrons and holes, which allows the longitudinal and Hall conductivities to be written as 80 where n e(h) is the carrier density of electrons (holes) and μ e(h) is the corresponding mobility. The relative magnetoresistivity is defined as . For the case of a compensated semimetal with n e = n h and equal mobilities for electrons and holes (μ e = μ h = μ), the above equations yield 2 2 This expression was used in this work to extract the magnetotransport mobility μ B from parabolic dependences of ρ NP (B) in small B.
Our analysis of the ρ xx (n)-peak broadening and the zero-field mobility μ 0 (see main text) have relied on theoretical expressions for the density of thermally excited electrons at the NP, n th . For MLG, this electron density is given by This expression can be obtained from the Boltzmann equations calculating the response of charge-neutral graphene to electric field and enforcing the resulting conductivity into a Drude-like form. It is noted that the above mass is proportional to the typical energy (thermal energy k B T of electrons and holes in the Dirac plasma) divided by their velocity squared, as expected for ultrarelativistic particles.
Using the same approach for BLG, we obtain its density of thermally exited electrons The above expressions for n th and m* have been used to evaluate conductivities of both charge-neutral MLG and BLG based on the standard Drude-like expression Additional examples of magnetotransport measurements for MLG Several (more than ten) monolayer devices (Hall bars and Corbino disks) were studied during the course of this work. To indicate variations in their magnetotransport behaviour, below we present measurements for another Hall bar (device D2) exhibiting notably higher remnant δn at low T. Its resistivity ρ NP (T ) at the NP is plotted in Extended Data Fig. 1a. Similar to device D1 (Fig. 2d), ρ NP of device D2 decreases with T and saturates above 200 K. In this device, the saturation occurs at higher T than in device D1 because of stronger inhomogeneity (compare Fig. 1d and Extended Data Fig. 1b). Despite an order-of-magnitude different inhomogeneities, both devices exhibit practically the same saturation value, ρ NP ≈ 1 kΩ. The same was valid for the other MLG devices. As discussed in the main text, we attribute the T-independent ρ NP in MLG to the entry of the Dirac plasma into the quantum-critical regime [3][4][5][9][10][11][12] . In this regime, the electron scattering time is determined by Heisenberg's uncertainty principle, where C is the interaction constant of about unity and depends on screening 3,4,11,12 . By plugging this scattering rate into equation (10) and using the effective mass from equation (8) and the carrier density given by equation (7), we obtain the quantum-critical resistivity NP 2 which is independent of T. The observed ρ NP ≈ 1 kΩ yields the interaction constant C ≈ 0.7, close to unity, as expected for Planckian-limit scattering 3-5,9-12 .
As for the MR behaviour of device D2, Extended Data Fig. 1c shows that ρ NP is parabolic in low B, similar to the case of device D1 in Fig. 1. The absolute value of Δ for device D2 is also similar, albeit slightly smaller, reaching 90% at 0.1 T and room temperature. The 20% reduction can be attributed to the lower electronic quality and homogeneity of device D2. Extended Data Fig. 1d shows zero-field and magnetotransport mobilities for device D2, which were extracted using the same approach as described in the main text. Both mobilities are slightly lower than those in Fig. 2. Nonetheless, at room temperature, μ B in device D2 still exceeds 100,000 cm 2 V −1 s −1 . Overall, the results presented in Extended Data Fig. 1 corroborate our conclusions that the Dirac plasma flaunts exceptionally high carrier mobility at elevated T, with no analogues among compensated metallic systems. The figure also reiterates the considerable differences between μ B and μ 0 , which were discussed in the main text and explained in 'Difference between zero-field and magnetotransport mobilities' below.

Electron-hole plasma in BLG
To emphasize how unique the Dirac plasma in MLG is, let us compare its magnetotransport properties with those of the closest electronic analogue, an e-h plasma at the NP in BLG. To this end, we fabricated and studied BLG devices that were also encapsulated in hBN to achieve high μ. They were double-gated and shaped into the standard Hall bars. At liquid-helium temperatures and away from the NP, the devices exhibited ballistic transport across their entire widths of about 10 μm. This was observed directly using bend resistance measurements 81 . The double-gating was required to tune the carrier density to the NP while maintaining zero bias between the two graphene layers. The latter ensured that no gap opened at the NP 82 , which otherwise would complicate the comparison 10 .
The typical behaviour of BLG's resistivity in zero B is shown in Extended Data Fig. 2a. Similar to MLG (Fig. 2d and Extended Data Fig. 1a), ρ NP (B = 0) of charge-neutral BLG reaches a few kiloohms at liquid-helium temperatures, but rapidly decreases to about 1 kΩ at higher T and becomes T independent above 50 K (Extended Data Fig. 2b). Such behaviour of high-quality BLG has already been reported recently, and constant ρ NP was attributed to the e-h plasma entering the quantum-critical regime 9,10 . Indeed, plugging the quantum-critical scattering rate τ C = k T h p −1 B into equation (10) and using the thermally excited density from equation (9), we obtain the resistivity for the e-h plasma in BLG The T independent value of ρ NP stems from the fact that both n th and scattering frequency τ p −1 evolve linearly with T. Equation (12) differs from equation (11) for MLG by only a factor of 2. From the data in Extended Data Fig. 2, we obtain C ≈ 1.4, close to unity as expected and in agreement with the previous reports 9, 10 . This value is two times larger than C for the Dirac plasma in MLG. We are unaware of any theory that would allow quantitative comparison between C in the two graphene systems. Nonetheless, the smaller value of the interaction constant in MLG compared with BLG could probably be understood as owing to the lower density of states in the Dirac spectrum.
In addition, we analysed δn(T ) for our BLG devices using the same approach as described for MLG in the main text. Above 50 K, δn in Extended Data Fig. 2c exceeds the remnant charge inhomogeneity (in the limit of low T ) by a few times, which ensures that the smearing of the peak in ρ xx at T > 100 K was dominated by e-h excitations. Extended Data Fig. 2c also shows that δn in BLG increased linearly with T, in agreement with equation (9) and qualitatively different from the quadratic behaviour of δn(T ) in MLG (equation (7) and Fig. 1d). Using the usually assumed value m* ≈ 0.03 m e for BLG (where m e is the free electron mass), we find δn ≈ 0.5 n th , similar to the case of MLG as discussed in the main text.
The response of BLG's e-h plasma to small B is shown in Extended Data Fig. 2d. Similar to the case of MLG, Δ evolves proportionally to B 2 but its absolute value is two orders of magnitude smaller than that in MLG, reaching only 1.5% at 0.1 T at room temperature. For completeness, we have evaluated the mobilities for the compensated e-h plasma in BLG, using the same approach as in the main text. Both magnetotransport and zero-field mobilities (μ B and μ 0 , respectively) are plotted in Extended Data Fig. 2e. They are found to be an order of magnitude lower than those for the Dirac plasma, which is the underlying reason behind the two-orders-of-magnitude smaller low-B MR in BLG compared with MLG (Δ ∝ μ 2 ). It is noted that μ 0 for BLG is approximately two times lower than μ B (Extended Data Fig. 2e), similar to the case of MLG in Fig. 2c. The difference between μ 0 and μ B is again attributed to electrons and holes moving against and along each other for longitudinal and Hall flows, respectively, as discussed in the main text and detailed in 'Difference between zero-field and magnetotransport mobilities' below.
Our experiments show that charge carriers in the Dirac plasma are several times more mobile than electrons and holes at the NP in BLG. The reason for the exceptionally high mobility in the Dirac plasma is twofold. First, the scattering rate τ C ∝ p −1 is approximately two times lower in MLG compared with BLG, as discussed above. Second, the effective mass for Dirac fermions at room temperature can be estimated from equation (8) as m* ≈ 0.01 m e , which is three times lower than the effective mass of charge carriers in BLG. Taken together, this suggests that the zero-field mobility µ eτ m = / * 0 for the Dirac plasma should be a factor of 6 higher than that for BLG's e-h plasma, in qualitative agreement with the experiment (compare Extended Data Figs. 1d and 2e).

Magnetotransport in multilayer graphene
Another electronic system with high-mobility charge carriers at room temperature is multilayer graphene (thin films of graphite). The material is an intrinsic semimetal with electrons and holes being present in approximately the same concentrations 83 . It is instructive to compare the magnetotransport properties of this nearly compensated semimetal with those of the Dirac plasma.
Our graphite devices were several nanometres thick (10-20 graphene layers) and shaped into Hall bars. To preserve the high electronic quality, the multilayer films were again encapsulated with hBN. Measurements for one of the devices are shown in Extended Data Fig. 3. Graphite's magnetoresistivity was found to increase quadratically in fields below 1 T. At room temperature, Δ was about 1.4% at 0.1 T, similar to the case of BLG and two orders of magnitude smaller than the MR of the Dirac plasma. Above 1 T, graphite exhibited notable deviations from the parabolic dependence bending towards a lower power and becoming practically linear in B at low T and above a few tesla. Room-temperature Δ reaches 80% and 3,500% at 1 T and 9 T, respectively, in agreement with a previous report for few-layer graphene 34 . Although MLG exhibits a few times larger Δ in high B, it is possible that the linear MR in graphite (first reported a century ago 42,43 and still not fully understood; see 'Brief history of linear MR') has the same origin as in MLG. This possibility requires further investigation because graphite's electronic spectrum is complicated and, also, strongly evolves with magnetic field 83 .
To evaluate the magnetotransport mobility μ B in graphite, we used the same approach as for MLG and BLG. The results are plotted in Extended Data Fig. 3b. At room temperature, μ B for the e-h system in graphite was found to be about 10,000 cm 2 V −1 s −1 , that is, a factor of more than 10 lower than that for the Dirac plasma in MLG (Fig. 2с) but close to μ B found for the e-h plasma in BLG (Extended Data Fig. 2e). This is perhaps not surprising as electronically, graphite is often considered as a stack of graphene bilayers. The provided comparison of graphene with its bilayers and multilayers highlights the unique nature of the Dirac plasma and its anomalously high mobility that results in the giant MR response, especially in low B. It is noted that μ B for multilayer graphene can be extracted more accurately, using both Hall and longitudinal measurements, which does not require the used assumption of e-h symmetry at the NP. The latter analysis 83 yields practically the same μ B as our intentionally simplified approach.

Difference between zero-field and magnetotransport mobilities
Magnetotransport in graphene's Dirac plasma was first analysed by Müller and Sachdev 84 and later by Narozhny with colleagues 85,86 . Below we provide analogous calculations, for completeness and to simplify our evaluation of the magnetoresistivity observed experimentally.
In the presence of electric E and magnetic B fields, the Boltzmann equations for electrons and holes at the NP can be written as where u e and u h are the drift velocities of electrons and holes, respectively, τ eh is the e-h scattering time and τ is the electron-impurity and/ or electron-phonon scattering times. The effective mass m* for the Dirac plasma is given by equation (8).
Taking the sum and difference between the top and bottom expressions in equation (13), we obtain As shown below, these coefficients determine the magnetotransport and zero-field mobilities. If equation (15) is placed into the left-hand side of the first line of equation (14) The first term defines the zero-B resistivity of the Dirac plasma and, as expected, depends on the total scattering rate 1/τ 0 . However, the second term is proportional to µ µ eτ m / = / * , that is, the absolute value of MR ρ xx (B) − ρ xx (0) is independent of e-h collisions and depends on only impurity and/or phonon scattering.
As for relative MR, we obtain The above analysis suggests different zero-field and magnetotransport mobilities, and their ratio is given by Our experiments found typical μ B /μ 0 of about 3, in agreement with the expectation that e-h scattering in the Dirac plasma should be the dominant scattering mechanism at room temperature.

Effect of proximity screening on mobility and MR
The observed difference between mobilities extracted from zero-field and magnetotransport measurements implies that μ 0 and μ B should be affected differently by screening. The latter mobility should be less sensitive to screening because e-h scattering does not contribute to Hall currents, as discussed above.
We have verified these expectations using MLG devices with proximity screening 15 . Such devices have previously been studied in the doped regime where electron scattering was found to be notably reduced by the screening 15 . Electron-hole interactions in charge-neutral graphene can also be expected to be modified by such proximity screening. We studied three MLG devices in which the graphite gate served as a metallic screening plate and was separated from graphene by a thin hBN layer (thicknesses of about 0.9 nm, 1.2 nm and 2.4 nm; inset of Extended Data Fig. 4a). In the particular case of the 2.4-nm device shown in Extended Data Fig. 4, we have found the screening to reduce ρ NP by a factor of about 2 below 250 K compared with similar-quality MLG devices without screening. The reduction in ρ NP yields a smaller interaction constant (about 0.4) and translates into higher μ 0 . It is noted that the difference between ρ NP observed for screened and unscreened devices reduces at higher T (Extended Data Fig. 4a). This can be attributed to the fact that the screening is sensitive to the average separation between charge carriers, which is proportional to n −1/2 . As the density of thermally excited carriers increases with T, the screening efficiency is reduced 15 .
The influence of proximity screening on magnetotransport in the Dirac plasma is found to be notably different from the case of zero B. Extended Data Fig. 4b shows that changes in ρ NP as a function of B remained practically the same for devices with and without screening. This agrees with the results in 'Difference between zero-field and magnetotransport mobilities', which predict that changes in ρ NP (B) should be insensitive to e-h scattering and, therefore, unaffected by proximity screening, in contrast to ρ NP (B = 0) that is dominated by this scattering mechanism.
For quantitative analysis of the observed screening effects, we have extracted e-h and electron-impurity (inelastic) scattering times (τ eh and τ, respectively) for the devices with and without proximity screening. To this end, we used the fact that the first (zero B) term in equation (18) depends on both τ eh and τ whereas the second term is determined only by τ. The results are plotted in Extended Data Fig. 4c. Both screened and unscreened devices exhibit similar τ that is, several times longer than τ eh . As expected, the proximity screening significantly suppresses electron interactions so that at about 150 K, τ eh is twice longer in the devices with proximity screening than for the standard encapsulated graphene. The difference is reduced at higher T, with possible reasons for this being mentioned earlier in this section.

Linear magnetoresistivity in high fields
As discussed in the main text, the parabolic MR is observed only in small magnetic fields up to about 0.1 T. In higher B, a linear MR behaviour emerges. We observed the linear dependence over a wide range of magnetic fields up to 18 T, the highest B available in our experiments. This is shown in Extended Data Fig. 5 for device D2. Again, the slope of ρ NP (B) depends weakly on T, and its absolute value is close to that exhibited by device D1 (within 20%), as shown in Fig. 3a. Overall, the described high-B behaviour was very similar for all five such MLG devices that we studied (inset of Fig. 3b). It is noted that the absence of T dependence for high-B MR indicates that many-body gaps caused by lifting of spin and valley degeneracies play little role within the discussed range of T and B. Otherwise, the gaps' smearing should have led to a strong T dependence.

Landau quantization at room temperature
We have attributed the observed linear MR in high B to the transition of the Dirac plasma into the quantized regime where the linear spectrum of MLG splits into dispersionless LLs. This condition is an essential prerequisite for discussing magnetotransport for the compensated Boltzmann gas in the zeroth LL.
In MLG, the main cyclotron gap at the filling factor ν = 2 in units of the kelvin (K) is given by 87 ℏ . The gap's size notably exceeds the thermal energy at room temperature already in fields of a few tesla. Previously, the Landau quantization has been reported for ultrahigh magnetic fields of 30-40 T where even the quantum Hall effect was observed at room temperature 87 . To demonstrate that Landau quantization in our devices becomes important at room temperature already in moderate B, Extended Data Fig. 6a shows the fan diagram measured for one of our Corbino devices at room temperature. The found peaks in inverse conductivity follow the main gaps at ν = ±2, as expected, and become clearly visible at B above 6 T. The Landau quantization is also visible in ρ xx measured in the standard Hall-bar geometry (Extended Data Fig. 6b). These observations support the description of high-B transport in neutral MLG in terms of the zeroth LL for the discussed temperature range up to 300 K.

Linear MR in Corbino devices
We have also used our Corbino devices to rule out edge effects in the appearance of strange linear MR. Extended Data Fig. 7 shows that the linear dependence ρ NP (B) was also observed in this geometry, exhibiting little difference with respect to the behaviour reported for the four-probe Hall-bar devices. Indeed, MR of Corbino disks is found to be weakly dependent on T and exhibit slopes with values close to those observed in the Hall-bar geometry (compare Fig. 3a and Extended Data Fig. 5). This proves that the linear MR is an intrinsic (bulk) effect and, for example, it is not related to e-h annihilation at graphene edges 66 or to spin and valley Hall currents reported for neutral graphene 88 .

MR of low-mobility graphene
To illustrate the importance of high quality for the reported MR behaviour of MLG in both low and high B, we have measured low-mobility devices obtained by exfoliation of graphene onto an oxidized silicon wafer (inset of Extended Data Fig. 8a). At liquid-helium temperatures, such devices exhibited strong charge inhomogeneity with δn ≈ 10 11 cm −2 (Extended Data Fig. 8a), which was nearly two orders of magnitude higher than that for hBN-encapsulated graphene (Fig. 1c). Even at 300 K, thermally excited density n th remained smaller than the residual δn, which means that electron transport near the NP in such devices was dominated by charge inhomogeneity (e-h puddles) at all T in the experiment. Accordingly, although ρ NP decreased with increasing T (Extended Data Fig. 8), similar to the case of our high-mobility devices, it only reached about 4 kΩ at room temperature, significantly away from the intrinsic value of about 1 kΩ for the Dirac plasma in the quantum-critical regime.
In small magnetic fields, ρ NP for MLG on silicon dioxide evolved quadratically with B (Extended Data Fig. 8b). The measured Δ was found to be more than two orders of magnitude smaller than in high-quality MLG (<1% at 0.1 T), which corresponds to about 8,500 cm 2 V −1 s −1 at the NP. With increasing B above 1 T, the MR of graphene on silicon dioxide deviated from the parabolic dependence and became sublinear at high T (Extended Data Fig. 8c), in agreement with the previous reports 30, 70 . Such sublinear behaviour may be attributed to short-range scattering 73 , which is present in graphene on silicon dioxide 89 , but further research is required to unambiguously identify the origins of high-B MR in low-mobility MLG. Nonetheless, our observations clearly show the importance of electronic quality for the observation of the linear magnetoresistivity.

Data availability
All relevant data are available from the corresponding authors. Source data are provided with this paper.