# Ballistic transport, chiral anomaly and emergence of the neutral electron - hole plasma in graphene

###### Abstract

The process of coherent creation of particle - hole excitations by an electric field in graphene is quantitatively described using a dynamic ”first quantized” approach. We calculate the evolution of current density, number of pairs and energy in ballistic regime using the tight binding model. The series in electric field strength up to third order in both DC and AC are calculated. We show how the physics far from the two Dirac points enters various physical quantities in linear response and how it is related to the chiral anomaly. The third harmonic generation and the imaginary part of conductivity are obtained. It is shown that at certain time scale the physical behaviour dramatically changes and the perturbation theory breaks down. Beyond the linear response physics is explored using an exact solution of the first quantized equations. While for small electric fields the I-V curve is linear characterized by the universal minimal resistivity , at the conductivity grows fast. The copious pair creation (with rate ), analogous to Schwinger’s electron - positron pair creation from vacuum in QED, leads to creation of the electron - hole plasma at ballistic times of order . This process is terminated by a relaxational recombination.

###### pacs:

81.05.Uw 73.20.Mf 73.23.Ad22.3.10

## I I. Introduction

It has been demonstrated recently that a graphene sheet, especially one suspended on leads, is one of the purest electronic systems. It became increasingly evident that electronic mobility in graphene is extremely large exceeding that in best semiconductor 2D systems GeimPRL08 . The scattering of charge carriers in suspended graphene samples of submicron length is so negligible that the transport is ballisticAndreiNN08 ; BolotinPRL08 . The ballistic flight time can be estimated as where is the graphene velocity characterizing the ”relativistic” spectrum of graphene near Dirac points and is the length of the sample. During the ballistic flight conductivity calculated neglecting interactions with phonons, ripplons, disorder and between electrons, etc., are consistent with experiments at least near the Dirac point at which no carriers are present Lewkowicz .

In a simplified model of a single graphene sheet (neglecting scattering processes and electron interactions) the chemical potential is located right between the valence and conductance bands and the Fermi ”surface” consists of two Dirac points of the Brillouin zone Castro . A lot of effort has been devoted theoretically to the question of dc and ac transport in pure graphene due to the surprising fact that the conductivity is finite without any dissipation process present. A widely accepted ”standard” value of the ”minimal dc conductivity” at zero temperature,

(1) |

was calculated using a simplified Dirac model Castro ; GusyninPRB06 ; Ando02 ; Ziegler06 ; Katsnelson06 ; Beenaker by a variety of methods. Direct application of the the Kubo formula at zero frequency, disorder strength, temperature and chemical potential GusyninPRB06 ; Ziegler06 utilizes certain ”regularizations” like the regulator which is interpreted as an ”infinitesimal disorder”. The regulator is removed at the end of the calculation. A more customary application of the Kubo formula starts with finite frequency. As noted by Ziegler and others Ziegler06 the order of limits makes a difference and several other values different from were provided for the same system. The standard value is obtained using a rather unorthodox procedure when the dc limit is made before the zero disorder strength limit is taken. If the order of limits is reversed, one obtains Ziegler06 :

(2) |

When the limit is taken holding one can even obtain a value of Ziegler06 , thus solving the ”missing ” problem (the same value was obtained very recently using yet another regularization in Beneventano09 ). Indeed, at least early experiments on graphene sheets on substrates provided values roughly 3 times larger than Novoselov05 . Recent experiments on suspended graphene AndreiNN08 demonstrated that the dc conductivity is lower, , as temperature is reduced to . Hence, while seems to be inconsistent with experiment, one still faces the question of what is the proper theoretical value. Since the conductivity of clean graphene in the infinite sample is a well defined physical quantity, there cannot be any ambiguity as to its theoretical value. Another theoretical approach to the problem was the use of Landauer formalism to the graphene sheet conductance. The conductivity appears as a limiting value of infinite width in this approach Katsnelson06 ; Beenaker and one obtains . Also resumming the series in disorder by a self consistent Born approximation gives Ando02 , while other resummations can lead to a different result FradkinPRB86 .

In contrast to the controversy with respect to the dc conductivity, both the experimental and the theoretical situation for the ac conductivity in the high frequency limit is quite different. The theoretically predicted value in the Dirac model is independent of frequency under condition Ando02 ; Varlamov07 . The Dirac model becomes inapplicable when is of order of or larger, where is the hopping energy of graphene. It was shown theoretically using the tight binding model and experimentally in GeimScience08 that the optical conductivity at frequencies higher or of order becomes slightly larger than . Moreover, in light transmittance measurements at frequencies down to it was found equal to within 4%. The tight binding model does not contain any other time scales capable of changing the limiting value of the ac conductivity all the way to . Therefore one would expect that the dc conductivity is rather than .

As is shown in the present paper, using the dynamical approach (a brief account of some results was published in Lewkowicz ; Rosenstein10 ), this is indeed the case. The basic physical effect of the electric field is a coherent creation of electron - hole pairs, mostly, but not exclusively, near the two Dirac points. To effectively describe this process we develop a dynamical approach to charge transport in clean graphene using the ”first quantized” approach to pair creation physics similar to that developed in relativistic physics Gitman to describe electron - positron pairs creation. To better visualize the phenomenon of resistivity without dissipation directly at zero temperature, doping, frequency and disorder we describe an experimental situation as closely as possible by calculating directly the time evolution of the electric current after switching on an electric field. In this way the use of a rather formal Kubo or Landauer formalism is avoided and as a result no regularizations are needed. Although we consider an infinite sample, the dynamical approach allows us to obtain qualitative results for finite samples by introducing time cutoffs like ballistic flight time. Various other factors determining transport can be conveniently characterized by time scales like the relaxation time for scattering of phonons or impurities.

We show in detail that some aspects the linear response physics are not dominated by the two Dirac points of the Brillouin zone at which the spectrum is gapless. For example, large contributions (infinite, when the size of the Brillouin zone ( is the lattice spacing), is being considered infinite) to the conductivity from the vicinity of the Dirac points are canceled by contributions from the region between them. This phenomenon is related to the chiral anomaly in field theory Smit . We analyse the use of massless Dirac (Weyl) model approximation to the tight binding model and propose an effective chirally invariant regularization for it.

In addition to the analysis of the linear response, we determin the range of applicability of the linear response approximation by both calculation of higher orders in (dc and ac) electric field and solving the model exactly in the nonlinear electromagnetic response regime. Only the zero chemical potential case (no net charge) is considered. In this respect the work is complimentary to an earlier study by Mikhailov Mikhailov in which ballistic nonlinear electromagnetic response to an ac field was calculated at finite chemical potential using the Boltzmann equation within a semi-classical approach and major effects we study were omitted. We first obtain perturbative corrections up to the third order in the electric field . At this order qualitatively new phenomena like the third harmonic generation and the imaginary part of conductivity appear. The calculation of the corrections allows us to estimate the time scale at which perturbation theory breaks down. At this scale, , the physical behaviour is expected to change qualitatively. Therefore one has to resort to nonperturbative methods. Certain aspects of nonlinear ac response at zero chemical potential were also studied recently by Mishchenko MishchenkoPRL09 . In his work disorder (taken into account phenomenologically) dominates the purely ballistic effects by cutting off the ballistic times at relaxation time scale before is reached.

Physics of the simplest tight binding model beyond the linear response is explored via an exact solution of the first quantized equations. It is a peculiar property of the tight binding model in the dynamical approach, that the exact solution of the equations for arbitrary momentum can be reduced to that for ; the constant electric field being in the direction. Moreover, the remaining equations, using Floquet theory, can then be reduced to a recursion relation. We use this property to effectively calculate the long ballistic flight evolution of various physical quantities like the current density and the number of created pairs. While for small ballistic times, , the conductivity settles at , at the current grows linearly. This increase can be explained using Schwinger’s electron - positron pair creation mechanism Schwinger . The pair creation is a two dimensional version of that in QED with the creation rate proportional to Cohen . This, in absence of relaxation channels (for times below ), leads to the creation of a neutral electron - hole plasma at ballistic times of order . This process cannot continue forever and is terminated by a relaxational recombination. The applicability of Schwinger’s formula for the electron - hole pairs creation rate was debated over a long time and we set the limitation on the applicability of this exact formula to graphene.

The rest of the paper is organized as follows. The tight binding model and the dynamical approach are described in section II. The perturbation theory in electric field is described in section III. The minimal conductivity is obtained and the role of states beyond the neighbourhoods of the Dirac points (and their relation to ”axial anomaly” and the Nielsen-Ninomya theorem) is clarified. The dynamical linear response approach to the ac field is considered. In Section IV perturbative corrections beyond linear response (up to third order in ) are calculated . The third harmonic and inductive contributions to conductivity are discussed. The exact solution for arbitrary field is constructed using the Floquet theory in Section V. It is used in Section VI to discuss the pair creation rate, conductivity and speculate about the electron - hole plasma. Finally Section VII contains discussion and conclusions.

## Ii II. The model and the dynamic approach to its physical properties

### ii.1 The tight binding model in an electric field

Electrons in graphene are described sufficiently accurately for our purposes by the 2D tight binding model of nearest neighbour interactions Castro . The Hamiltonian in space is

(3) |

where

(4) |

with being the hopping energy; and

(5) |

with vector potential for (). Later it is generalized to more general fields including the ac field.

### ii.2 The dynamical ”first quantized” approach.

Since the crucial physical effect of the field is mainly a coherent creation of electron - hole pairs (though not always near the Dirac points, see below), a convenient formalism to describe the pair creation is the ”first quantized” formulation described in detail in Gitman . The spectrum before the electric field is switched on is divided into positive and negative energy parts describing the valence and conduction band:

(6) | |||||

(7) |

where . A second quantized state is uniquely characterized by the first quantized amplitude,

(8) |

which is a ”spinor” in the sublattice space. It obeys the matrix Schroedinger equation in sublattice space

(9) |

where is defined in Eq.(5). The initial condition corresponding to a second quantized state at in which all the negative energy states are occupied and all the positive energy states are empty is

(10) |

A physical quantity is usually conveniently written in terms of . We will be interested mostly in the current density (multiplied by factor due to spin)

(11) |

as well as the energy and the density of the electron - hole pairs. The summation is over the honeycomb-lattice Brillouin zone, see Fig. 1, in which two Brillouin zones are outlined. The energy of the electrons is changing due to the applied electric field (with no dissipation in the ballistic regime, see however Discussion) and is given by

(12) |

The amplitude of lifting an electron into the conduction band (defined for the Hamiltonian without the electric field) is , where the scalar product is defined by , and consequently the density of pairs reads

(13) |

Since the Hamiltonian depends on time via the vector potential that shifts the momentum , a more useful definition of the density of pairs taking into account the shift of the momentum will be given in Section III.

### ii.3 Units and conventions

The units of energy, time and distance are defined by the microscopic values and correspondingly. The unit of the electric field will be , so that the dimensionless electric field is and the unit of conductivity will be . Effectively one can set . In these units the first quantized equation reads

(14) | |||||

where the tight binding Hamiltonian matrix element takes the form

(15) |

where and for the dc field. We will use extensively its expansion in powers of :

(16) |

For example, the off - diagonal matrix element for the component of the current is , so that the current density, using the first quantized approach, is

(17) |

The next two sections deal with the perturbative treatment of the electric field and later (after having shown that it fails at certain time scale) we will switch to a nonperturbative method.

## Iii III. Linear response, the pseudo - diffusive behaviour and the parity anomaly

### iii.1 Expansion in electric field

Expanding in the dimensionless electric field strength (assumed for simplicity homogeneous and constant after switching on the field at ), , one obtains for the first correction the following equation:

(18) |

Writing the correction as , it takes the form

(19) | |||||

with initial conditions , . We use the abbreviations

(20) |

and

(21) |

The coefficient denotes the amplitude of accumulation of electrons in the conduction band, whereas denotes the amplitude of remaining in the valence band. Solving the equations one obtains

(22) |

The expansion for the current, Eq.(11), in the electric field up to the first order contains the following momentum contributions, :

(23) |

The zero order contribution is

(24) |

where

(25) |

The correction gives the ’paramagnetic’ and ’diamagnetic’ contributions to the current densities:

(27) |

Both corrections for a specific momentum diverge at large ballistic times as , however one still has to integrate over the valence band momenta.

### iii.2 Integration over momenta and physical interpretation of the ”quasi -Ohmic” behaviour.

To first order in electric field the current density is

(28) |

where the zero order contribution can be written as a derivative with respect to of a periodic function

(29) |

It vanishes upon integration, in accordance with the Bloch theorem, since the Brillouin zone can be chosen in such a way that it exhibits periodicity in the field () direction. For example, we can integrate over the rectangular area shown in Fig.1 containing two Brillouin zones:

(30) |

It should be noted that at every and due to the continuity of even at the Dirac points. Therefore one is left with the linear response.

The conductivity can be divided into an apparently linearly divergent part and a finite part, . As will be explained shortly, the contributions to the first term at large stem mostly from the immediate neighbourhoods of the two Dirac points, while contributions to the second come from whole Brillouin zone. The part diverging linearly with time can be written (after some algebra) as a full derivative with respect to :

(31) | |||||

This too vanishes upon integration over the Brillouin zone, since it is again an integral of a derivative of a periodic function, albeit this cancellation is nontrivial. In Fig. 2a the integrand of Eq.(31) is plotted within the integration area specified above, Eq.(30).The integrand is negative along the line connecting two Dirac points along the field direction, for example

(32) |

(recall that in our units), and positive elsewhere. In Fig. 2b several cross sections for various values are shown. One observes, that as approaches the integrand goes to at corresponding to Dirac point . Consequently for all the integrand is a regular function, and thus the integral over vanishes. At the integral over is finite, yet does not influence the two-dimensional integral.

The remaining part

(33) |

gives a finite result. Unlike the divergent part of the conductivity, the integral for (in units of ) is dominated by contributions from the vicinity of the two Dirac points (circles in Fig.2a of radius (in units of )). Indeed the prefactor is bounded from above by and integral over the area outside the circles vanishes at . The contribution of a single Dirac point is obtained by integrating to infinity (in polar coordinates centered at the Dirac point),

At long time the upper limit can be replaced by infinity which results

(35) |

Here one can safely multiply the result by the valley degeneracy . Returning to physical units, one obtains , Eq.(2).

The dependence on time was calculated numerically. After an initial increase on the natural time scale the conductivity approaches via oscillations, see Fig. 1 in ref. Lewkowicz . The amplitude of oscillations decays as a power

A physical picture of this resistivity without dissipation is as follows. The electric field creates electron - hole excitations in the vicinity of the Dirac points in which electrons behave as massless relativistic fermions with the graphene velocity playing the role of the light velocity. For such particles the absolute value of the velocity is and cannot be altered by the electric field and is not related to the wave vector . On the other hand, the orientation of the velocity is influenced by the applied field. The electric current is , thus depending on orientation, so that its projection on the field direction is increased by the field. An important observation, made for example in ref.Sachdev , is that the electric field, while creating charge transport, does not change the overall momentum. As a consequence the effects of impurities do not affect the pair creation in the same way they affect ”free” carriers. These two sources of current, namely creation and reorientation are roughly of the same size in the linear response. As we will see below, the situation changes at ballistic times when the linear response breaks down.

### iii.3 Pair creation rate

In analogy to the linear response in the current, the pair creation rate per unit area, as defined in Eq.(13), can be calculated and to leading order in is:

(36) |

The result of the numerical integration over the Brillouin zone was given in Fig. 2 of ref. Rosenstein10 . It is dominated by the two Dirac points and at large times (compared to ) behaves as

(37) |

It is however well known that a dc electric field in QED (even with massive relativistic electrons) renders the vacuum unstable Schwinger ; Gitman , with the ”renormalized” number of pairs depending on . Therefore nonperturbative effects are expected and will be studied below .

The notion of ”renormalized” number of pairs is a consequence of the fact that for unstable systems the Fermi level should be continuously ”renormalized” as function of time. In the first quantized formalism this corresponds to a continuous modification of the ”initial” state defining holes. This leads to the following definition of the pair density Nussinov ,

(38) |

It is used in the framework of the relativistic Dirac model and was recently extensively discussed in ref.Cohen in connection with graphene (using the instanton approximation).

### iii.4 ac field

A similar calculation within the dynamic first quantization formalism for the evolution of the current density can be performed for an arbitrary time dependence of the homogeneous electric field. Let us consider an arbitrary time dependence of the component of the vector potential , subject to a ”switching on” condition, , for . Fourier transforms are defined by

(39) |

and similarly for the current density and other quantities.The next to leading order in first quantized tight binding equations are:

(40) |

The switching on condition, , can be taken into account by the substitution regularizing the way the Fourier transform is defined for zero frequency. From Eq.(40) one obtains

(41) |

In the same order the ac conductivity is an integral over the sum of the ”paramagnetic” and the ”diamagnetic” contributions given by

(42) |

Subtracting a full derivative (independent of frequency), one can rewrite this as

(43) |

Real and imaginary parts, taking the definition into account, are

(44) | |||||

Therefore ac and dc conductivity are equal, as was also shown in recent calculationsMacDonald . This is consistent with both the Kubo formula derivations Varlamov07 and optical experiments GeimScience08 . Similar results were obtained in multiple layered graphene MacDonald .

A similar calculation can be performed directly in the space (similar to the dc case) and shows that after a short transient one obtains the value for the ac conductivity independent of frequency.

## Iv IV. Beyond linear response. Weyl fermions, parity anomaly and the Nielsen - Ninomiya theorem

As will be discussed below, an electric field should not necessarily be very small in graphene experiments, so that corrections to linear response are of interest. In addition it is important to determine at what ballistic time scale perturbation theory breaks down. In this section we present a perturbative calculation of the leading nonlinear effect in both dc and ac, and discuss the ways to regularize the Weyl model ”correctly” in order to calculate these corrections.

### iv.1 Ultrarelativistic fermions near Dirac points

The tight binding model employed here has two Dirac points around which the spectrum becomes ultra - relativistic , , see Fig.1. In our units the graphene velocity is . The matrix element of the tight binding Hamiltonian can be expanded around as

(45) |

The Weyl field describing the ”left handed”Weinberg ; Smit fermions , namely a field satisfying the relativistic Dirac equation with zero fermion mass

(46) | |||||

can be constructed by the unitary transformation

(47) |

The matrix element around the second Dirac point is different,

(48) |

and the Weyl fermions are ”right handed” particles that obey

(49) | |||||

without any unitary transformation required. They belong to a different representation of the 2+1 dimensional Lorentz group Luscher .

This is not just a peculiarity of the model, but a rather general feature GusyninPRB06 of massless fermions on a lattice Smit ; Nielsen . It is well known that in order to get a massless spectrum of fermionic excitations with any ultraviolet cutoff (hexagonal lattice is an example), they come in multiple locations on the Brillouin zone (species ”doubling”). In the Hamiltonian formulation the multiplicity is , where is the dimensionality of space, in graphene. In addition, the graphene fermions are ”staggered”, meaning the spinor components ”live” on different sublattices of the honeycomb lattice. This reduces the doubling to . The doubling is intimately linked to the parity (discrete chiral symmetry) Luscher ; Smit . The two Dirac points have opposite chirality so that there is no ”parity anomaly”.

It is sometimes claimed in condensed matter literature that, at least while doing linear response, one can concentrate on the two neighbourhoods of the Dirac points and neglect the rest of the Brillouin zone. Even more, often the calculation’s result is just multiplied by the factor (the valley degeneracy). Below we show that this is generally not an accurate description of what happens. This is not an academic question, since the low energy Weyl theory is simpler than the tight binding model and will be used in the next section to extend the calculations into higher orders in the electric field. One therefore needs a proper ultraviolet regularization of the Weyl model. Similar regularization issues are well known in field theory whenever chiral anomaly is encountered. Roughly, a ”correct” regularization should respect the chiral symmetry leading to important constraints on the number and charges of massless fermions. Otherwise the model is called ”anomalous” and results of perturbative calculations become arbitrary. The most famous example is the requirement that in each generation of elementary particles (leptons and quarks) the sum of charges should be zero Peskin .

We therefore discuss in some detail the cancellation of infrared divergencies and the correct application of the ”ultra - relativistic” approximation to the tight binding model.

### iv.2 Cancellation of ultraviolet divergencies and the approximate chiral symmetry.

The ultra - relativistic approximation, Eq.(45), fails when applied to the first term in the expression for conductivity given in Eq.(31). At first glance the integral in Eq.(31) is dominated by the two Dirac points since the integrand diverges there, see Fig. 2. For the both (widely separated) regions of the Brillouin zone, and , it has the same asymptotic form

(50) |

The integral over the neighbourhood of each Dirac point converges in the ”infrared” (here meaning ) due to the Jacobian , but is linearly divergent in the ”ultraviolet” and both have the same sign, see Fig.2. Hence the integration cannot be extended to infinity and the size of the Brillouin zone serves as a natural ultraviolet cutoff.

It is important to note that there is no cancellation of the divergence between the Dirac points since both have the same sign! The divergencies are however canceled by the contributions from regions of the Brillouin zone between the Dirac points in which the ultra - relativistic approximation is not valid. Therefore in the ”ohmic” regime during ballistic times one is not allowed to neglect states far from the Dirac points and replacement by the Weyl equation is incorrect. Due to cancellation of the whole ”divergent” term, Eq.(31), one can devise an appropriate regularization in the UV in which these contributions are cancelled and even generalize the procedure to higher orders in the electric field. We propose and use such a scheme below in section IVC. Simple recipes like the momentum cutoff regularization (circles in Fig.1) or giving the fermions an infinitesimal mass GusyninPRB06 , commonly used, may lead to unphysical results. As a consequence more sophisticated regularizations like the -function regularization Beneventano09 , ”dimensional regularization” Peskin and ”Pauli - Villars” were developed for continuous field theories like the Weyl model. The tight binding model is very similar in this respect to the ”lattice Hamiltonian” regularizations of the field theory (Hamiltonian meant here with time kept continuous) and also satisfies the chiral invariance criterion Smit .

### iv.3 Nonlinear response in dc. Where does the linear response break down?

Since the current density is an odd function of the electric field, the leading nonlinear correction to it appears in the third order in the field. The calculation along the lines of subsection IIIA is quite straightforward albeit tedious. One can use the Weyl model instead of the tight binding model due to the following reason. It is well known in field theory that generally chiral anomaly effects, including the ultraviolet divergencies discussed in section IIID, appear only in one loop calculations Peskin . We checked and found that indeed up to the third order the expression is finite in the ultraviolet. Within the dynamical approach it is natural to perform the calculations without Fourier transforming into the space.

Up to the third order the current density is

(51) |

The correction therefore is growing as a large power of the ballistic time. It becomes as large as the leading order for . Hence the perturbation theory breaks down on the time scale of

(52) |

As will be seen in Section V, this agrees well with the crossover time obtained from a nonperturbative calculation. Of course this time should be larger than any other possible relaxation time (due to impurities, phonons etc.) present in the system. A similar expansion was obtained for the number of electron-hole pairs Rosenstein10 .

The result is the same for the tight binding and the corresponding Weyl model. The only issue would be the first order calculation within the Weyl scheme since it is divergent in the ultraviolet. In Weyl model we used the momentum cutoff that will be described in section IVB.

### iv.4 ac response and the third harmonic generation

An analogous calculation for the ac field switched on at results in

(53) | |||||

The first term is just the reflection of the initial conditions and is non - universal. In the limit one recovers the dc result. The corrected value of the real (dissipative) part of the ac conductivity is

(54) |

One observes that the perturbation theory is inapplicable for

(55) |

This is less restrictive compared to the dc condition, Eq.(52) for .

The present formulas can be compared with the results of the dynamical calculation beyond linear response by Mishchenko MishchenkoPRL09 in which a phenomenological model of relaxation was employed. For the inverse relaxation time , his nonperturbative result is

This is consistent with the second correction terms in Eq.(54).

The inductive part was absent in linear response and now appears for the first time:

(56) |

Due to the two conditions, Eq.(55), it is much smaller than . The third harmonic is generated with the real part

while the inductive part is absent at the present order.

## V V. The exact solution of the first quantized tight binding equations

### v.1 Application of the Floquet theory and reduction to 1D

It is a peculiar property of the tight binding matrix Eq.(4) that the solution for arbitrary can be reduced to that for . Shifting the time variable to , one can define an ’universal wave function’

(57) |

The Schroedinger equation (14) is now void of any dependence,

(58) | |||||

where the dimensionless frequency is in units of . The component of the momentum enters the solution via the initial conditions Eq.(10) only. In terms of the universal function it takes the form

(59) |

Taking the Fourier transform of Eqs. (58), the equations in frequency space turn into recursion relations:

(60) | |||||

(61) |

Floquet theory (and the recursion relations) assure that the frequencies are discrete and form two series both indexed by an integer

(62) |

The ”central” frequency will take generally two incommensurate values. Writing the Fourier amplitude as , one obtains, combining the two equations, Eq.(60) and Eq.(61), the forward and the backward recursions,

(63) | |||||

(64) | |||||

with normalization . The and have to be determined by the consistency between the forward and the backward equations.

Generally there are two independent Floquet frequencies, however our system is special with the exact relation

(65) |

which can be proven rigorously. This greatly simplifies the solution. To obtain a simple approximation, we expand around the easily solvable case of corresponding to . In this case the solution consists of two frequencies

(66) |

where the superscript denotes zeroth order in . For experimentally accessible cases , the Floquet frequencies are close to . We expand the Fourier coefficients in powers of :

(67) |

where we do not distinguish for the time being between the two Floquet branches. Let us look for a solution with the initial choice

(68) |

which in particular means that at this order. It can be shown that the ”central” frequency does not have a first order correction in . Using this fact the forward recursion can be rewritten as:

Inserting the expansions for and into this relation yields the expression for ,

(70) | |||||

so that the positive Floquet frequency (up to second order) is . Moreover, it turns out that the expansion in converges in the whole relevant range, , and greatly assists the numerical solution of the recursion relation.

After the coefficients and Floquet frequencies are found, the solution of Eq. (58) is written in terms of the Fourier series:

(71) | |||||

The second line follows from Eq.(61), are defined in Eqs.(62) and (LABEL:pert) and we utilized Eq.(65) for simplification. The two complex coefficients are fixed by the initial conditions Eq.(59). This solution is used to calculate the evolution of the current density, the energy and the number of electron - hole pairs.

### v.2 Results. Nonlinear regime and Bloch oscillations.

The current density divided by the electric field (in physical units), , is shown in Fig. 3 in ref. Rosenstein10 for various values of the dimensionless electric field . The (microscopic) unit of the electric field is very large , so that at realistic fields . After an initial fast increase on the microscopic time scale , approaches the universal value the conductivity rises linearly above the constant ”universal” value : and settles there, consistent with the linear response, calculated in Section IIIB. Beyond the crossover time

(72) |

where . This is consistent with the Landau - Zener calculation (instantons) within the Weyl model GavrilovPRD96 ; Dora10 , see below. The crossover time, defined by , is therefore

(73) |

and is consistent with the perturbation theory breakdown ballistic time, Eq.(52). This linear increase regime can be considered as a precursor of the Bloch oscillations, but is still universal in the sense that it depends solely on neighborhood of the Dirac points.

The current on the time scale of order

(74) |

exhibits Bloch oscillations, as is seen in Fig.4 of ref. Rosenstein10 , where larger fields in the range are shown. It turns out that the current vanishes at points given exactly at multiples of with being the period of the Bloch oscillations. One observes a peculiar feature that (apart from the ”relativistic” initial constant segment) the time dependence of is similar for different electric fields. Indeed, if one plots versus , all the curves nearly coincide. Moreover, one obtains the following excellent fit for the current

(75) |

The Bloch time is approximately the time required for the electric field to shift the momentum across the Brillouin zone . This time scale is very long for experimentally achieved fields, much longer than the ballistic flight time. For example in a sample of submicron length, , . If one assumes that the electric current can reach , so for a typical width one has an electric field corresponding to (the voltage in such case would be quite large ). The first maximum of the Bloch oscillation will be seen at flight time of , which is of the same order as . If one uses a value of the current typical to transport measurements , the electric field is just corresponding (voltage ), and is therefore out of reach. See, however a recent proposal Dragoman08 .

### v.3 Crossover at in the Weyl model

It is expected that the transition to the nonlinear regime is dominated by the neighborhoods of the Dirac points. Therefore it can be obtained also within the Weyl model provided it is properly regularized in the ultraviolet region consistent with chiral symmetry, as was discussed in section IV. For the calculation of the current density it suffices to renormalize the electric current by subtracting the UV divergent term Eq.(50):

(76) |

The results are practically indistinguishable compared to the tight binding model for times larger than a microscopic time scale and much smaller than the Bloch time. Some examples of the tight binding model were presented in Fig. 3 of ref. Rosenstein10 . In Fig. 3 the conductivity in units of