{'fraction_non_alphanumeric': 0.03445692883895131, 'fraction_numerical': 0.00850187265917603, 'mean_word_length': 4.22933803368586, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non linear magnetotransport theory and Hall induced resistance oscillations in graphene

The quantum oscillations of nonlinear magnetoresistance in graphene that occurs in response to a dc current bias are investigated. We present a theoretical model for the nonlinear magnetotransport of graphene carriers. The model is based on the exact solution of the effective Dirac equation in crossed electric and magnetic fields, while the effects of randomly distributed impurities are perturbatively added. To compute the nonlinear current we develop a covariant formulation of the migration center theory. The analysis of the differential resistivity in the large magnetic field region, shows that the extrema of the Shubnikov de Hass oscillations invert when the dc currents exceeds a threshold value. This results are in good agreement with the experimental observations. At small magnetic field, the existence of Hall induced resistance oscillations are predicted for ultra clean graphene samples. These oscillations originate from Landau-Zener transitions between Landau levels, that are tilted by the strong electric Hall field.

in crossed electric and magnetic fields, while the effects of randomly distributed impurities are perturbatively added. To compute the nonlinear current we develop a covariant formulation of the migration center theory. The analysis of the differential resistivity in the large magnetic field region, shows that the extrema of the Shubnikov de Hass oscillations invert when the dc currents exceeds a threshold value. These results are in good agreement with the experimental observations. In the small magnetic field regime, corresponding to large filling factors, the existence of Hall induced resistance oscillations are predicted for ultra clean graphene samples. These oscillations originate from Landau-Zener tunneling between Landau levels, that are tilted by the strong electric Hall field.

INTRODUCTION
Since the discovery of graphene 1 , a lot of excitement has been generated, due to its unique electronic properties and its potential applications in electronic devices 2,3 . The dynamic of charge carriers in graphene is described by a Dirac-like Hamiltonian, in the vicinity of the so-called Dirac points (where the conduction and valence band merge), leading to an energy spectrum that has linear dependence on momentum 4,5 . The linear dispersion relation, simulates the physics of quantum relativistic effects for chiral massless fermions, except for the fact that in graphene the Dirac fermions move with a speed v F , which is 300 times smaller than the speed of light c. Hence, graphene exhibits a variety of pseudo relativistic phenomena, providing a connection between condensed matter physics and quantum-relativistic phenomena 6 . Among others we can cite: the Zitterbewegung and its relation with the minimal electrical conductivity at vanishing carrier concentration, the unconventional quantum Hall effect, and the Klein tunneling 6 .
The study of magnetic field effects on the electron gas dynamics has attracted the attention of physicists for long time. Magnetic oscillations, such as the Shubnikov-de Haas (SdH) oscillations of the electric conductivity, and the de Haas-van Alphen oscillations of magnetization, are a useful tool to study the shape of the Fermi surface 7 . Furthermore, it has been found that magnetotransport phenomena are enhanced when the electron gas is confined to two dimensions, leading to one of the most remarkable phenomena in condensed matter: the quantum Hall effect (QHE) 8 . More recently, non-equilibrium and non-linear magnetotransport phenomena in high mobility two-dimensional electron gases (2DEG) have acquired great experimental and theoretical interest. Microwave-induced resistance oscillations (MIRO) and zero resistance states were discovered 9,10 in 2DEG subjected to microwave irradiation and moderate magnetic fields. An analogous effect, Hall field-induced resistance oscillations (HIRO) and zero differential resistance states, have been observed in high mobility samples in response to a dc-current excitation 11,12 . Analogous phenomena to those described for conventional 2DEG in the previous paragraph, have been discovered in graphene, including the integer 13,14 and fractional quantum Hall effect 15 . However the non-linear magnetotransport effects in graphene are less studied.
A recent experiment 16 has analyzed the effect of a strong dc current bias on SdH oscillations in graphene. A phase inversion on the SdH oscillations of the differential resistance is observed with increasing bias, in which the oscillation maxima evolve into minima and vice versa. On the work of Wang and Lei 17 the nonlinear effects on graphene magnetotransport are studied within the balance-equation approach.
In this paper we analyze a theoretical model for the nonlinear magnetotransport of graphene carriers. The model is based on the exact solution of the effective Dirac equation in crossed electric and magnetic fields, while the effects of randomly distributed impurities are perturbatively added. We take advantage of the relativistic-like structure of the model to compute the nonlinear current. With this purpose we develop a covariant formulation of the migration center theory. Based on this model the quantum oscillations of nonlinear magnetoresistance in graphene that occur in response to a dc current bias are investigated.
The analysis of the differential resistivity in the large magnetic field region shows, in good agreement with the experimental observations 16 , that the extrema of the Shubnikov de Hass oscillations invert when the dc currents exceeds a threshold value. At small magnetic field, the existence of Hall induced resistance oscillations are predicted for ultra clean graphene samples. These oscillations originate from Landau-Zener transitions between Landau levels, that are tilted by the strong electric Hall field.
The outline of the paper is presented as follows. In Sec. II the effective Dirac equation for charge carriers in graphene subject to crossed electric and magnetic fields, as well as randomly distributed impurities is described. We present the method, based in Lorentz-like transformations, that allow us to solve the problem. Section III describes the migration center theory, and the computation of the impurity assisted current. The details of the current controlled scheme that lead to the expression for the longitudinal resistivity are discussed in Sec. IV. In Sec. V we present the main results and discussion of the nonlinear effects in SdH oscillations as well as the Hall induced resistance oscillations in graphene.
The formulation of the covariant migration center theory is described in the appendix.

II. MODEL
We consider the motion of charge carriers in an imperfect graphene sample, in the presence of an in-plane electric field and a magnetic field perpendicular to the graphene plane. The graphene defects are modeled by the impurity scattering potential V . The goal of this work is to study the nonlinear response to a DC electrical current density (J dc ) of graphene placed in a magnetic field. In a typical experimental configuration the electric field is not explicitly controlled, instead the longitudinal current is fixed, J x = J dc , while the transverse current cancels, J y = 0. Then, in general, the electric field has both longitudinal (E x ) and transverse (E y ) components, that have to be determined in terms of J dc . The calculations are simplified if we perform a rotation by an angle φ = Arctan(E y /E x ) to a S L frame with coordinates r = (x , y ):, hence, the electric field points along the y axis: E = (0, E). In section III the current density J is calculated from the wave function solution that is obtained in the present section, an inverse rotation to the original S L frame will lead to a set of two implicit equation for the unknown E x and E y , that can be determined by a self-consistent iteration in terms of J dc .
The dynamic of the problem in the S L frame is governed a Dirac equation, that for convenience is written as where v F is the Fermi velocity and σ ≡ (σ x , σ y ) are the Pauli matrices that act on the two component spinor Ψ, the pseudo-spin degrees of freedom are associated to the two graphene sub-lattices. In the previous equation Π 0 = p 0 − e v F φ and Π = p − e c A. When necessary, we use a covariant notation with space-time 2 + 1-vectors, e.g: x µ = (v F t, r ); µ = 0, 1, 2. The momentum operator is p µ = i ∂/∂x 0 , −i ∇ , and Π µ = p µ − e c A µ . Notice that the gauge potential for constant electric and magnetic fields is given as: A µ = − c v F Ey , −By , 0 . The Dirac matrices are selected as γ ν = (σ z , iσ y , −iσ x ) 18 .
The impurity potential is decomposed in terms of its Fourier components where r i represents the impurity positions, that are chosen at random; N i is the number of impurities; and the δ parameter determines the rate at which the perturbation is turned on and off. The Fourier transform of the disorder potential, for screened Coulomb interactions, can be approximated by 19 V (q) = 2πe 2 / (q + 4α g k F ) where k F is the Fermi momentum, α g = e 2 / v F is the effective fine structure constant and is the material permittivity.
However, as far as we do not consider regions very close to the Dirac point 19 , the short range scattering approximation can be used, leading to a constant value V (q) This is possible as long as the condition cE < v F B is fulfilled.
We follow the procedure presented in the work of Lukose et. al. 20 and extend the method to take into account the impurity scattering effects. In order to cancel the electric field effect, a boost along the x axis is performed, with a boost parameter θ given by Under the Lorentz boost, coordinates transform as covariant three-vectors:x µ = Λ ν µ x ν ; where the explicit components of Λ ν µ can be read from the following relations: Whereas, the wave function transforms as a Dirac spinor By re expressing the Dirac equation (1) in terms ofΨ(x) with the aid of Eqs. (5) and (4) one finds thatΨ(x) satisfies an equation with the same structure withΠ µ =p µ − e cÃ µ , but now the gauge potential is given asÃ µ = 0, −Bỹ, 0 , wherẽ B = B 1 − β 2 . Hence in the boosted frame the charge carrier dynamics corresponds to a problem with vanishing electric field and a reduced magnetic fieldB. Whereas the transformed impurity potential takes the form here γ ≡ cosh θ = 1/ 1 − β 2 , and the matrix M is given by Notice that in the S A frame the impurity potential acquires an additional time-dependent term; the boost induces a superposition of oscillatory terms, each with a momentum depen- When the impurity potential is not included, Eq. (6) is exactly solvable. The eigenstates and spectrum are given by the well-known relativistic Landau levels (LL): here the Landau number n is an integer and the quantum numberk x corresponds to the translation symmetry along the x direction. The function φ n are the displaced harmonic oscillator eigenfunctions, centered atỸ kx = l 2 βk x ; and the magnetic length in the S A system, l β , is related to the magnetic length l c = c /eB by the relation l β = l c γ 1/2 . Notice that the energy spectrum is degenerate with respect tok x . Applying the inverse boost to the solution in Eq. (8), yield the eigenfunctions and spectrum in the original S L frame: notice that, as expected, the electric field breaks the LL degeneracy, the Landau levels are tilted with respect to k x with a slope given by v F β. The displaced harmonic oscillator wave function is now centered at Y n,kx , given as The wave function center position depends both on k x and on the LL number n, as a result the LL are mixed, in particular the electric field mixes the particle and hole solutions, consequently the eigenfunctions are not orthogonal.
In order to obtain the current, within the linear response formalism, it is convenient to use an orthogonal base, otherwise the algebra turns rather cumbersome. To avoid this problem and since the eigenfunctions in S A are orthogonal, it is convenient to work out the current densityJ µ in the S A frame, and then apply a Lorentz transformation to obtain the current density J µ in the original S L frame.

III. CURRENT DENSITY IN THE NONLINEAR REGIME
We now turn to the calculation of the current density in the boosted S A system where the dynamic is ruled by Eq. (6). To consider the dynamic of an electron in a magnetic field it is customary to use the migration center theory [21][22][23] , where the electron coordinate r splits into the guiding center X = (X, Y ) and the cyclotron or relative coordinate η, i.e.
With this decomposition the η x and η y coordinates become non-commutative, similarly for X and Y : It is possible to develop a covariant extension of migration center theory using the Schwinger proper time formalism. In what follows we present the general aspects of the covariant migration center theory, details are worked out in the appendix. The covariant expression for the relative coordinate η µ is selected as where the electromagnetic tensor F µν is constructed from the constant E and B fields as given in Eq.(36), and the quantity the Lorentz transformation. In the S L frame the 2 + 1-vector η µ is given by where we recall that the electric field points along the y−axis and β = cE/v F B. Instead in the S A frame, the relative coordinate is given byη µ = (c/eB)(0, −Π y ,Π x ). It is easily demonstrated that the relative coordinate transform as a covariant 2 + 1-vector under the Lorentz transformations Eqs. (4):η µ = Λ ν µ η ν . The three-vector guiding center coordinate is then defined as X µ = x µ − η µ . The guiding center velocity is obtained from the Heisenberg equation of motion (see the appendix), in the boosted S A frame it reads: It is a straightforward exercise to prove (see the appendix) that the current density probability J µ defined in terms of the guiding center velocity as J µ = Ψ † γ 0 dX µ /dτ Ψ, is also a Lorentz 2 + 1-vector.
We now turn to the calculation of the current density. The current density J µ is computed from the impurity, thermal, and time average of the guiding center velocity where S is the sample area. In the S A frame the density matrix ρ satisfies the Von Neumann equation: The Hamiltonian in Eq. (6) splits into the first part that is exactly solved, Eq.(8), and the impurity potentialṼ , the latter is treated as a perturbation. The density matrix is also decom- where f n is the Fermi distribution function. Since the eigenfunctions |ν = |n,k x are orthogonal and complete, the matrix elements of ∆ρ are readily obtained as: E n is the energy of the state |n,k x given in Eq. (8). If we consider randomly distributed impurities, the terms that have linear dependence on the impurity potential vanish once we take the average over the impurities. In the auxiliary reference frame S A the current space components get contributions that arise solely from impurity scattering, it is worked out as: where n i is the impurity density, the LL degeneracy per unit area is given by g k = 1/(2πl 2 c ) and e 3 is an unitary vector in the direction of the magnetic field. The transition matrix element |V (q)| 2 |Θ mn (ζ)| 2 includes the function Θ mn = γ D |m|,|n| + sgn(nm)D |m|−1,|n|−1 − βγ sgn(m)D |m|−1,|n| + sgn(n)D |m|,|n|−1 , (17) here the matrix elements of the displacement operator D nm , are given by the following where ζ = l β √ 2 (q x + iq y ). Finally, the function G nm displays a Lorentzian distribution the effect of the quasiparticle scattering is introduced by means of δ that produces the broadening of the Landau levels. In the δ → 0 limit, G nm embodies the energy conservation condition as discussed below. Due to the elastic scattering, the impurities transfer a momentum q x =k x −k x to the electrons. Recalling that the center of the oscillatory wave function is given asỸ = l 2 βk x , it follows that the momentum transfer is equivalent to a hopping or shifting of the guiding center ∆Ỹ =Ỹ −Ỹ = l 2 β q x . In the δ → 0 limit Eq. (19) gives the resonance conditionẼ n =Ẽ m + v F βq x . If interpreted in the S A system, the inter-Landau transitions are induced by the time dependent impurity potential Eq. ( 7) that oscillates with the frequency ω q = v F βq x . Alternatively, the results can be interpreted according to the displacement mechanism: using the momentum transfer q x =k x −k x and noticing that E n,kx = γ Ẽ n − v F βk x , the resonance condition can be rewritten as an energy conserving condition in the S L frame E n,kx = E m,k x , and the transition can be interpreted as a Landau-Zener tunneling between LL, tilted by the electric field Eq. (9) 11 .
When the Lorentzian distribution is considered in Eq. (16) the summation over the n and m indices does not converge. While it is customary to set a cut-off energy to avoid the convergence issue 24 , we substitute the Lorentzian distribution with a Gaussian distribution, leading to a natural cut-off: Finally, the thermal and impurity average of the time component of the current density in Eq. (14) is easily computed utilizing Eqs. (13) and (8), yieldingJ 0 = eñ e v F , whereñ e is the electron density in the boosted frame. IV.

CURRENT CONTROLLED SCHEME
We are now ready to obtain the current density J in the laboratory frame, from the corresponding current densityJ µ in the auxiliary frame.J µ was calculated in the previous section: the time component is simply given asJ 0 = eñ e v F ; whereas the corresponding space components are obtained by computing the expressions contained in Eq. (16). It is easy to verify thatJ imp x cancels because of the angular integrations; hence, we only need to compute the componentJ imp y . Applying the inverse Lorentz transformation the current J µ in the S L frame is obtained as: where the electron density n e is related to the electron density in the boosted frame according to the relation n e = γñ e . Recalling that the imposed dc-current points along the x− axis, we need to perform an additional rotation by an angle φ = Arctan(E y /E x ) (see the beginning of Sec. II). Consequently these leads to the following conditions: The last equalities apply for the Hall configuration: the longitudinal current is fixed at the value J x = J dc , while the transverse current cancels. Eqs. (21) represent two implicit equations for the unknown E x and E y . The equations can be solved following a self-consistent iteration method. However, it is easily verified that for a typical Hall experiment a correct approximate solution to the previous equation is obtained assuming that the electric field is almost transverse with respect to the dc bias current: i.e. E x E y or φ ≈ π/2. Retaining terms to first order in = π/2 − φ, the explicit solution is worked out as The first equality is obtained from the leading contribution, it sheds the well known result for the Hall electric field E H induced by the dc bias in the magnetic field. The second equation follows from the next order contribution, it gives the longitudinal electric field in terms of the impurity scattering current obtained in Eq. (16); the result allow us to calculate the longitudinal resistivity, ρ xx = E x /J dc , as well as the longitudinal resistivity differential It is interesting to notice that the β parameter, Eq. (3), takes the form: β H = J dc /en e v F , so β H can be understood as the rate between the current density and charge carriers density moving at a top speed v F . Clearly, the formalism breaks down when J dc > en e v F .

V. RESULTS AND DISCUSSION
Results are presented for the differential resistivity r xx both in the strong and weak magnetic field regimes. As explained in section III, the displacement mechanism is the physical origin of the nonlinear effects: impurity scattering induces a spatial hopping of the charge carriers between tilted Landau levels. However, there is an additional effect that has to be consider. Recent experimental results 16 suggest that in graphene a current dc bias produces a strong electron heating effect. It was found that the charge carriers effective temperature depends linearly on the bias current as T ∼ T 0 + αJ dc . The parameter α was In a typical graphene sample with a mobility µ ≈ 10 4 cm 2 V −1 s−1, the dimensionless parameter is estimated asδ ≈ 0.14. In Fig. (1) we present the longitudinal resistivity ρ xx (B) for various values of the dc current J dc . As observed in the experiments 16 , the amplitude of the SdH oscillations decreases with increasing temperature. Fig. (2) displays the SdH oscillation for the differential resistivity as function of B and for the same values of J dc .
Besides the damping effect on the oscillations amplitude, it is observed that as J dc increases the maxima of the oscillation evolve into minima, and vice versa. These nonlinear effects on SdH oscillations originate mainly from the current heating effect. If we compare the results of r xx (J dc ) without the heating effect (α = 0), with the result of the linear response approximation r xx (0) we find essentially the same results, confirming the fact that the displacement mechanism has a negligible effect on the SdH oscillations.
We now explore the possibility of detecting nonlinear oscillations in the low magnetic region. The study of ultra clean semiconductor samples 2DEG has shown that dc current excitation leads to Hall field-induced resistance oscillations 11,12,28,29 . In order to simulate a high mobility graphene sample the value of the LL width in Eq. (24) (16). For this purpose we consider the limit δ → 0, in Eq. (19) or Eq. (20), that leads to the resonance condition here we are considering levels well above the Dirac point, so all the Landau numbers n and m are positive. The momentum transfer ∆k is equivalent to a hopping or shifting of the guiding center ∆Y = Y − Y = l 2 β q x . An estimation of the momentum transfer ∆k can be obtained from the overlap of the wave functions between the n and m LL. The overlap, obtained in section III is proportional to the function D nm (ζ)in Eq. (18), where ζ = l β √ 2 (q x + iq y ). While this functions is highly oscillatory, the envelope of the function can be approximated by 30,31 n! m! ζ m−n e −|ζ| 2 /2 L m−n n (|ζ| 2 ) → ζ m+n e −|ζ| 2 /2 √ n! + m! , with a maxima at ζ max = √ m + n. Hence, the maximum overlap corresponds to a momen-tum transfer of ∆k = 2(m + n)/l β . Inserting this result in Eq. (25) the following selection rule is obtained In HIRO we are interested in large LL; so we consider contributions from states that are close to the Fermi level, that is far away from the Dirac point, so n F ≈ πl 2 c n e /2. As we assume that n F 1 we write m = n F + m and n = n F + n , with (m , n ) n F . So the resonance condition (26) becomes with N = m − n . We recall that, even though the linear spectrum in graphene implies zero rest mass, their cyclotron mass is not zero, but it is given as 13 m * = √ πne v F . Hence, using the definition of the cyclotron frequency, ω c ≡ eB m * c , the resonance condition can be written as The HIRO oscillations are periodic in 1/B, according to the condition (28)  It was mentioned in section (III) that according to the interpretation in the S A system, the inter LL transitions are induced by the time dependent potential Eq. (7) that oscillates with the frequency ω q = v F βq x . The transition matrix element peaks at the momentum transfer q x = q * = 2 √ n F /l c , hence we find that the leading oscillation frequency of the periodic potential, is precisely the Hall frequency ω q * ≡ ω H , given in Eq. (28).
It is also interesting to observe that from Eq. (25) the hopping distance for a resonant transition between LL is expressed as that depends both on the magnetic field and J dc . This equation can be interpreted as a selection rule that determines ∆Y for arbitrary transition between LL. If we now consider the limit of large LL levels, the selection rule reduces to HIRO may be observable in ultra clean graphene with mobilities above µ ≈ 10 5 cm 2 V −1 s−1.
Finally we also mention that it will be interesting to look after the possibility of producing zero-differential resistance states, similar to those observed in 2DEG 11 .

Appendix: COVARIANT MIGRATION CENTER THEORY
In this appendix, we work out the details of the covariant migration center formalism.
First it is convenient to write the equation of motion in a manifestly covariant form where we recall that the Dirac matrices are select as γ ν = (σ z , iσ y , −iσ x ) 18 The Green function G(x, x ) is obtained in terms of the unitary evolution operator U (x, x , τ ) by means of the relation The evolution operator satisfies the equation with the following boundary conditions: For the purpose of this paper, we are mainly interested in deducing the covariant guiding center velocities. First we consider the Heisenberg equation of motion for the x µ coordinates, it gives The Dirac probability current, is then written as j µ =ψγ µ ψ = ψ † γ 0 γ µ ψ. It is easily verified that the probability current transform as a covariant 2 + 1-vector under the Lorentz-like transformations Eqs. (4,5):j µ = Λ µ ν j ν . As explained in section III the electron coordinate x µ splits into the guiding center X µ and the cyclotron or relative coordinate η µ , i.e. x µ = X µ + η µ , where a manifestly covariant expression for the relative coordinate η µ is given in terms of the velocity operator Π µ and the the electromagnetic tensor F µν in Eq. (11); where F µν is constructed from the constant fields: c v F E and B as The quantity With the aid or the Heisenberg equations of motion and using the definition forη µ , the relative proper velocity dη µ /dτ can be obtained both in the S L and in the S A frames. In the boosted S A frame we obtain: where e 3 is an unitary vector in the direction of the magnetic field. On the other hand, dη µ /dτ evaluated in the S L sheds the following result notice that the temporal component dη 0 /dτ is obtained in terms of spatial x-component.
It is a lengthly but straightforward exercise to demonstrate that the probability current defined in terms of the relative velocityψ (dη µ /dτ ) ψ transform as a covariant three-vector.
Summarizing all these results, and taking into account that the relative coordinate is obtained as X µ = x µ − η µ , it immediately follows that the current density probability J µ defined in terms of the guiding center velocity as J µ = Ψ † γ 0 (dX µ /dτ ) Ψ, is also a Lorentz three-vector. The guiding center velocity dX µ /dτ in the S A frame as given in Eq. (13)