Electronic response of graphene to an ultrashort intense terahertz radiation pulse

We have recently reported a study (Ishikawa 2010 Phys. Rev. B 82 201402) on a nonlinear optical response of graphene to a normally incident terahertz radiation pulse within the massless Dirac fermion (MDF) picture, where we have derived physically transparent graphene Bloch equations (GBE). Here we extend it to the tight-binding (TB) model and oblique incidence. The derived equations indicate that interband transitions are governed by the temporal variation of the spinor phase along the electron path in the momentum space and predominantly take place when the electron passes near the Dirac point. At normal incidence, the equations for electron dynamics within the TB model can be cast into the same form of GBE as for the MDF model. At oblique incidence, the equations automatically incorporate photon drag and satisfy the continuity equation for electron density. Single-electron dynamics strongly depend on the model and pulse parameters, but the rapid variations are averaged out after momentum-space integration. Direct current remaining after the pulse is generated in graphene irradiated by an intense monocycle terahertz pulse, even if it is linearly polarized and normally incident. The generated current depends on the carrier-envelope phase, pulse intensity and Fermi energy in a complex manner.


Introduction
There has been rising interest in graphene [1,2] over a broad spectrum of fields because of its potential application in carbon-based electronics as well as the possibility of mimicking and testing quantum relativistic phenomena [3][4][5][6]. While unique properties including finite conductivity at zero carrier concentration [7] and ac and dc universal conductance [8][9][10] are predicted and observed, interest in the optical response of graphene is even further boosted by the advent of terahertz (THz) radiation technology [11,12]. Especially, remarkable progress in the development of intense THz sources has enabled the generation of phasestable, monocycle to a few cycles, transients with field amplitudes exceeding 100 kV cm −1 and even 1 MV cm −1 [12][13][14][15][16]. This has opened up a new field of high-field physics in condensed matter, especially graphene [17,18]. Along these lines, a theoretical description of nonlinear optical responses of graphene has become one of the key issues [19][20][21][22][23][24][25][26][27][28][29][30]. Floquet analysis has provided valuable information, e.g. on dressed band structure, gap opening [25,27,28] and frequency-dependence of the nonlinear response [20,22,23]. However, this approach, assuming a continuous wave, cannot be directly applied to the case of ultrashort pulse irradiation. As for time-domain microscopic theories, kinetic approaches [18,19] based on the Boltzmann equation and second-quantization approaches [29,30] have been reported.
One distinct aspect of THz-pulse interaction with graphene, compared with the case of near-infrared and visible pulses, is that the photon energy is much smaller than that required for vertical transition in most of the momentum space. On the other hand, a THz pulse can induce a large displacement of electrons in the momentum space: transition from one valley to another is even possible. We have recently studied [26] nonlinear optical responses of graphene, starting from the time-dependent Dirac equation (TDDE) and focusing on the electronic motion in the momentum space (intraband dynamics) and the interband transition along its trajectory. We have shown that the TDDE can be cast into a form of generalized optical Bloch equations, referred to as graphene Bloch equations (GBE) hereafter, which describe the interplay between intraband and interband dynamics in a physically transparent fashion. The previous study [26] was, however, limited to the massless Dirac fermion (MDF) picture and normal incidence.
In this study, we extend it to the tight-binding (TB) model with nearest-neighbor interactions and oblique incidence. Starting from the time-dependent Schrödinger equation (TDSE), we formulate an exact, fully nonperturbative, nonlinear theory of single-electron dynamics and derive formulas for the induced current, applicable to an arbitrary waveform and polarization of THz pulse, both for normal (section 2) and oblique incidence (section 3). In the case of normal incidence, the derived equations can be further cast into the same form of GBE as within the MDF approximation (section 2). Then, applying them to the case of linear polarization, we show that interband transitions take place predominantly when the electron passes near the Dirac point along its trajectory in the momentum space. We examine how the results depend on the model, pulse parameters and initial electron momentum (section 4). Also, we briefly analyze the single electron dynamics when it takes a circular path around the Dirac point from the adiabatic limit involving Berry's phase [31,32] to the nonadiabatic limit leading to full population oscillation (section 5). Finally, we study the macroscopic electric current induced by an intense monocycle THz pulse and show that the carrier generation through the enhanced interband transition near the Dirac point leads to the generation of a direct current (dc) that remains after the end of the pulse (section 6).

Normal incidence
In this section, we treat normal incidence of an optical pulse whose electric field E(t) and vector potential A(t) = − E(t) dt is in the graphene plane (x y plane). The pulse may be of arbitrary time-dependent polarization, while [26] only considered linear polarization.

Temporal evolution of the electronic wave function
The TDSE for the two-component wave function ψ(t) of the electron with an initial wave vector of k (p =hk is the canonical momentum) is given by where the time-dependent Hamiltonian H (t) is of a form where (t) = |h(t)| and θ (t) = − arg h(t).
where π(t) = p + eA(t) is the kinetic momentum with e (> 0) being the elementary charge. is the magnitude of the energy eigenvalue for the field-free case; its value k is given by which we depict in figure 1(a). Sufficiently near the Dirac point, the MDF picture is valid [5]; equation (3) can be approximated as where v F = 3γ a 2h ≈ c/300 and the Dirac point is taken as the origin of π. Then, we obtain and θ corresponds to the directional angle of κ satisfying κ x = |κ|cos θ and κ y = |κ|sin θ (around K), −|κ|sin θ (around K ).
If h(t) is constant independent of time, for both the TB and MDF models, the TDSE equation (1) has the following two solutions: whose energy eigenvalues are ± , where the upper and lower signs correspond to the upper band (electron band) and the lower band (hole band), respectively.
For the general case of time-dependent h(t), let us make the following ansatz: i.e. a linear combination of the instantaneous upper and lower band states (Volkov states) where the instantaneous temporal phase or dynamical phase (t) is defined as Substituting equations (8) and (9) into equation (1), we obtain, as a temporal variation of the expansion coefficients c ± (t), Introducing the interband coherence ρ = c + c * − and population difference n = |c + | 2 − |c − | 2 , one can rewrite equation (11) into the form of GBE [26]: We previously derived equations (11)-(13) for the MDF model [26], but they are equally valid for the TB model. While equations (11)-(13) basically describe interband transitions, they incorporate field-induced intraband dynamics through (t) and θ(t) and are physically more transparent than the TDSE. At the same time, they are totally equivalent with the TDSE involving no approximation. The electron dynamics described by the GBEs may be experimentally probed, e.g. by time-resolved pump-probe photoemission spectroscopy [34].
It should be noted that if θ (t) were defined by the principal value θ(t) = −Arg h(t), with Arg z defined in the interval (−π, π], as plotted in figure 1(b), θ(t) would undergo 2π jumps on the white lines linking the Dirac points in figure 1(b). Here, instead, we are to define θ(t) = − arg h(t), with arg z = Arg z + 2πn (n is an integer), in such a way that it varies continuously along the path of κ(t). Then, if κ(t) eventually returns to its original position in the k-space along a trajectory surrounding a Dirac point, ψ ± (t) acquires a geometrical phase of π, in addition to the dynamical phase (t). Berry's phase [31,32] is, hence, incorporated in the GBEs (see section 5).
It is noteworthy that GBEs do not explicitly contain optical matrix elements but that the interband transitions are governed byθ , which is, especially within the MDF picture, the temporal variation of the polar angle along the path that κ(t) takes on the k plane. This indicates that an optical-pulse irradiation can be viewed as a means to move an electron along a path in the k-space. The coupling elementθ (t) is, in the MDF picture, explicitly given bẏ which is nonlinear in E(t).

Electric current
The electric current by a single electron is given by −ej e , where j e is defined as With the solutions ρ(t) and n(t) of the GBEs (12) and (13), equation (15) can be rewritten as, for the case of the TB model, where the first and second terms correspond to the contribution from the intraband current and interband polarization, respectively. To take into account the Fermi distribution, we solve equations (12) and (13) with initial Fermi-Dirac function with µ, k B and T being the chemical potential, Boltzmann constant and temperature, respectively. Then, we can obtain the macroscopic electric current J(t) by integration over the honey-comb lattice Brillouin zone depicted in figure 1(a) as where g s = 2 denotes the spin-degeneracy factor. The carrier current j c is calculated by replacing n with the carrier occupation n + 1 in equation (16).
In the MDF picture, each component of equation (15) can be explicitly written as The macroscopic electric current J(t) is given by where g v = 2 denotes the valley-degeneracy factor. The carrier current j c is calculated by replacing n with the carrier occupation n + 1 in equations (18) and (19), i.e.
If we use the MDF model, the whole dynamics is invariant under the multiplication of quantities of energy dimension,hω, v F p, v F e A, µ, k B T andh/t, by a common factor. In the weakfield limit, one can derive the universal conductivity e 2 /4h from equations (20)-(22) (see the appendix).

Oblique incidence
In this section, we extend the formulation developed in the previous section to the case of oblique incidence. Let us consider an oblique incidence of a pulse whose phase velocity vector in the x y plane is u = (u x , u y ), related to the carrier frequency ω and the wave vector q in the x y plane as q = (ω/u x , ω/u y ). Note that |u| = c/sinχ, with c and χ being the vacuum velocity of light and the angle of incidence, respectively. The limit u x , u y → ∞, i.e. q → 0 corresponds to normal incidence. Figure 2. Schematic representation of the relation between p and ± for the case of p and q parallel to the y-axis.

Massless Dirac fermion picture
Let us begin with the MDF picture. The TDSE reads Through variable transformation (t, x, y) → (τ, x, y) with τ ≡ t − x/u x − y/u y , i.e. in a frame of reference comoving with the pulse, equation (23) can be rewritten as where σ and I denote the Pauli and identity matrix, respectively. This equation can also be obtained by replacing π(t) in equation (5) with π(τ ) + ih q ω ∂ ∂τ . To obtain insight into equation (24), let us consider the field-free case, i.e. A(τ ) = 0 and π(τ ) = p. The eigenvalue equation for ψ of temporal and spatial dependence ∝ exp which can be solved analytically. It should be noticed that replacingp = p + ω q with p would lead to the regular eigenvalue equation (6) and that the momentum p in the τ -frame corresponds to p + ω q in the t-frame. These circumstances are schematically depicted in figure 2. Let us denote the two solutions of equation (25) as +,p (> 0) and −,p (< 0). These two states would be coupled by a photon that has the same in-plane phase velocity u as the incident pulse. This observation indicates that equation (24) automatically incorporates photon drag [35].
We now consider the solution of equation (24) in the presence of an external field A(τ ). Again, let us write ψ(τ ) as a linear combination of the instantaneous upper and lower band states: where and θ ± (τ ) denote the directional angle of vector π + ±,π ω q. Unlike for normal incidence, the dynamical phase ± and directional angle θ ± are distinct for each band. Then, we obtain, as the temporal evolution of the expansion coefficients c ± , where = + − − and θ = θ + − θ − . In contrast to equation (11), these coupled equations contain diagonal terms (the second terms) and cannot be cast into a simple form of GBE. One recovers equation (11) for the limit of normal incidence, i.e. q → 0, for which θ → 0. It should be noted that the dynamics described by these equations are not unitary, i.e. a(τ ) = |c + (τ )| 2 + |c − (τ )| 2 is not constant and a(τ ) = 1 even after the pulse in general, whose physical implication will be discussed below in section 4.3.
The single electron current can be evaluated by of which each component can be explicitly written as whereθ = (θ + + θ − )/2 and ρ = c + c * − . Again, the first and second terms correspond to intraband currents and the third terms to the interband polarization oscillating in time. Let us denote the electric current for a single electron initially in the upper (lower) band as j + e (j − e ) and also denote the current by ψ − (τ ) as j 0 e (τ ). Then, the macroscopic current density can be calculated as

Tight-binding model
Here we follow an approach similar to the one in the previous subsection. The TDSE in the τ -frame is where π = p + eA(t − q · r/ω) = p + eA(τ ). In the τ -frame, h π in H π is given by For the field-free case, this can be rewritten as then the energy eigenvalue satisfies Hence, as in the case of Dirac fermion, equation (34) incorporates photon drag. Let us rewrite equation (34) as denote the two solutions of equation (37) as + (> 0) and − (< 0) and make ansatz equations (26) and (27) with θ ± being the phase angle of h * π . Here, h π is given by equation (36) with p replaced by π. Then, equation (29) holds as coupled differential equations describing the temporal evolution of the expansion coefficients c ± (τ ).
The single electron current equation (15) can be evaluated as which tends to equation (16) if q → 0; the macroscopic current density can be calculated by

Comparison with other approaches
It may be worth, at this stage, briefly comparing our approach with some others. Common approaches consider a distribution function, say f k , in the k-space and look at its temporal evolution. This view may be said to be analogous in spirit to the Eulerian specification in hydrodynamics, i.e. a way of looking at fluid motion that focuses on specific locations in the space. References [18,19,29,30] belong to this class. In contrast, our approach, which may be called a dynamic first-quantized approach [22], follows the motion of each electron in the k-space starting from the TDSE, thus, analogous to the Lagrangian specification to follow an individual fluid parcel as it moves. While the two views can, in principle, be transformed to each other, our approach helps us to intuitively visualize the large displacement of each electron (in the k-space) induced by an intense THz radiation field and allows us, for instance, to obtain simple solutions for a circular path around the Dirac point (see section 5). Moreover, whereas [18,19,29,30] apparently assume normal incidence, here we have also presented a treatment of oblique incidence and can account for the effect of photon momentum (see figure 11 below), which have been neglected in those studies. Kao et al [22] and López-Rodríguez and Naumis [27,28] took a Lagrangian approach and analyzed the cases of constant electric fields [22] and continuous waves [22,27,28]. The latter authors also considered oblique incidence. Our approach basically includes these as special cases. Moreover, we have derived equations (11) and (29) of motion for the expansion coefficients c ± , which provide us with a clear physical view.
The microscopic approaches in [18,29,30] include carrier-carrier and/or carrier-phonon scattering processes, which have been neglected in this study. Although further theoretical developments are required, it is expected that our approach may, in principle, treat these effects microscopically as jumps in the k-space. Alternatively, the mean-field effect may be taken into account through a self-consistent treatment as in [19].

Linear polarization
In this section, we examine the electron dynamics and generated current predicted by our theory developed in the previous sections, for the case of linear polarization.

Initial momentum parallel to the field
Let us begin with the simplest case where the canonical momentum is parallel to the field, i.e. p//A(t), for the case of normal incidence.
In the MDF picture,θ (t) is nonzero only at the Dirac point where θ(t) changes stepwise. As can be easily verified, equation (1) has the following two analytical solutions: where θ k denotes the directional angle of k, i.e. the initial value of θ andn ≡ (cos θ k , sin θ k ).
One can also verify that these solutions satisfy equations (8)- (11). It can be easily shown that, for each branch, the current j e = v F ψ † σ ψ is given by respectively. j e remains constant, i.e. shows no response for the case of p//A(t). This holds true even when the instantaneous kinetic momentum π = p + eA(t) passes by the  For example, for k = (k, 0) and an incident light linearly polarized along the x axis, (t) is given by

Electron dynamics and current at normal incidence
In order to gain insight into how the electron behaves when π (t) passes near the Dirac point, let us consider a flat-top pulse with a half-cycle ramp-on, normally incident and linearly polarized along the x-axis: where A 0 > | p x |. We numerically solve the GBEs for an electron initially in the lower band using the Bulirsch-Stoer method with adaptive stepsize control [36] for parameters p x /e A 0 = −3/4,hω/v F e A 0 = 9.46 × 10 −3 and different values of p y /e A 0 , where p x and p y are measured from the Dirac point K = ( 2π 3 √ 3a , 2π 3a ) (figure 1). We set γ = 2.8 eV in this section. These parameters correspond, e.g. to radiation of 0.25 (1,4) THz frequency with a field strength of 2.0 (31.5, 504) kV cm −1 . Figure 3 plots the calculated temporal evolution of the population difference n(t) between the upper (conduction) and lower (valence) states. As we have seen in equation (14), the coupling elementθ (t) that governs interband transition is nonlinear in E(t) and, when p y is much smaller than e A 0 , strongly peaks within a narrow time window around π x (t) = p x + A(t) = 0. This leads to a quasi-step-like population transfer (enhanced interband transitions) when π x (t) changes its sign, e.g. when the electron momentum passes near a Dirac point, while We also compare the results of the MDF and TB models in figure 3. For fixed values of p x /e A 0 , p y /e A 0 andhω/v F e A 0 , the MDF model gives the same result, as we have already mentioned. The results of the TB model, on the other hand, depend on the pulse frequency (or field strength). We can see that all the four solid lines overlap with each other till ≈ 1.3 cycles (before the second population jump), since the first jump takes place near the Dirac point, where the MDF model is a good approximation. After the first jump, the electron is in a superposition of the upper and lower states. Then, the phase difference accumulated in ρ (equation (12)) depends on whether one treats the dynamics within the MDF or TB model and further on field strength in the latter, throughθ and . Therefore, the population dynamics sensitively depends on the simulation conditions. Figure 4 illustrates the dependence on pulse polarization up to four optical cycles. We denote the polarization direction by the angle ξ measured counterclockwise from the x direction and use the initial momentum of the electron rotated by the same angle ξ around the Dirac point K , to make a consistent comparison. Again, although the result is virtually independent of polarization before the second interband jump, the dynamics is very sensitive to ξ thereafter. One also notices that, at the lower intensity (panel (a)), the MDF is still a relatively good approximation before the third jump ( 2.1 cycles), compared with the case of the higher intensity (panel (b)). In figure 5, we show the results for longer times up to sixty cycles. Remarkably, the population dynamics exhibits quasi (but not completely)-periodic behavior. Its period varies with the polarization angle. For example, it is about 22 cycles at ξ = 0 • , while about six cycles at 45 • .  (20)) and TB (equation (17)) models at different polarization angles ξ in the latter. The y component vanishes for the MDF model and the TB model at ξ = 0 • due to symmetry and virtually vanishes also for the other values of ξ . We see that the current does not linearly follow the ramp-up of the pulse shown in figure 3, indicating a nonlinear nature of the electronic response of graphene, which will be discussed further in section 6. The MDF model overestimates the current only slightly; the variation with the polarization direction of the pulse within the TB model is negligibly small. This observation may come as a surprise, especially after having seen the sensitive dependence in figures 4 and 5, but can be understood as follows. The sensitiveness of the single-electron response to the model and pulse conditions is due to that of the phase of ρ (equation (12)) accumulated throughθ and along the electron path in the momentum space. The latter is expectedly affected also by the initial momentum p. Indeed, in figure 7, which illustrates the dependence of the x component of single-electron fluence | j c,x (t)| 2 dt integrated up to six cycles on initial momentum, we see that the contribution to the current rapidly changes with the initial momentum and forms a complicated interference pattern, which also depends on whether we calculate using the MDF or TB model. Integrated over the momentum space, such fine interference as well as the sensitive dependence in figures 4 and 5 are mostly averaged out, resulting in the weak dependence seen in figure 6.

Single-electron dynamics at oblique incidence
We have examined normal incidence so far. Let us now turn to oblique incidence. Here, we consider a pulse propagating in the y direction (u x = 0 and u y = c/sinχ) and polarized in the x direction whose vector potential is given by where A 0 (τ ) denotes the envelope function. Specifically, we assume a central frequency of 10 THz, a peak intensity of 10 7 W cm −2 corresponding to a 87 kV cm −1 field strength and a Gaussian temporal profile with a full-width at half-maximum (FWHM) of 300 fs. Figure 8 plots the population dynamics, calculated within the MDF model, of an electron initially in the valence band with an initial t-frame momentump = p + ω q ofp x = 0 andp y = A max /5, where A max denotes the peak vector potential. We compare the results for χ = π 4 and 0. Both, composed of enhanced interband transitions near the Dirac point, are virtually indistinguishable from each other. In figure 9 we show the single-electron intraband current j intra e , i.e. the first and second terms of equations (31) and (32). Again, the result for oblique incidence (panel (a)) is indistinguishable from that for normal incidence (panel (b)) on the scale of the figures. These observations are due to the fact that the carrier velocity v F is so much smaller than the propagation velocity c/sinχ of the pulse that the field seen by the electron is practically unchanged by the spatial displacement.
A closer look into the total population a(τ ), however, reveals a fine but important difference, as can be seen in figure 10. It remains at unity for the case of normal incidence, since equation (11) describes a unitary temporal evolution. As we have mentioned in section 3, on the other hand, a(τ ) is not constant at oblique incidence and does not return to unity even   where ∇ is defined on the x y plane. This equation can be transformed into in the τ -frame, from which we obtain When the pulse propagates along the y direction, this equation can be rewritten as Another effect can be found if one looks at the sum of j intra e,y forp = (0, A max /5) and (0, −A max /5), which is plotted in figure 11. While the sum vanishes due to symmetry at normal incidence, it is finite at oblique incidence due to the net momentum exchange with the pulse upon interband transition (figure 2). We also notice that both the deviation of the total population from unity in figure 10 and the sum of j intra e,y in figure 11 are nearly exactly proportional to sinχ.

Circular path surrounding the Dirac point
In this short section, let us examine an important, analytically solvable problem, namely, how the electron behaves when its kinetic momentum vector π (t) varies along a path circulating around the Dirac point in such a way that Within the MDF picture, the electron takes a circular path around the Dirac point and such a situation may typically be realized under a normal incidence of a circularly polarized pulse. In this case, equation (11) can be transformed into a form of the Rabi oscillatioṅ where η ≡ 2 /hω, which can be analytically solved. For an initial condition c + (0) = 1 and c − (0) = 0, we obtain and then In terms of an analogy to Rabi oscillations, η plays the role of detuning. Figure 12 shows the population |c + (θ)| 2 of the conduction band.
In the opposite limit η → 0, corresponding to a quick rotation near the Dirac point, the solutions read as The electron population is completely transferred to the lower band at θ(t) = π . Moreover, whether the electron arrives at the other side of the Dirac point by a counter-clockwise [θ(t) = π] or clockwise [θ (t) = −π ] half rotation, it ends up with the same wave function, in contrast to the adiabatic limit. These observations are consistent with our discussion in subsection 4.1.
After a complete rotation [θ (t) = 2π ], the electron population is transferred back to the upper band, but, since c + = −1 instead of unity, Berry's phase is cancelled by a phase acquired through the interband dynamics, again in contrast to the adiabatic limit. It should be noticed that equation (56) also holds for any path of π(t) sufficiently close to the Dirac point. Indeed, if we approximate (t) as vanishing, it is easy to see that the following solutions: satisfy equation (11). The results in figure 12 can be used to roughly estimate the vertical width of the distribution presented in figure 7. The population transfer is significant, say, for η 1. Using the definition of η ≡ 2 /hω, the first part of equations (51), (6) and (14) with E y = 0, A y = 0 and π x = 0, one can transform this into where p y is here measured from the Dirac point and E 0 denotes the field amplitude. This criterion corresponds to a full-width of ∼0.06 for the case of figure 7, which is indeed found to be a reasonable estimate.

Direct current generation by a monocycle terahertz radiation pulse
In this section we investigate the response of graphene to an ultrashort (typically monocycle) intense THz pulse and predict that a direct current is generated that remains after the pulse, even when the pulse is linearly polarized and normally incident. Let us consider a pulse whose vector potential is given by where A 0 (t) and φ CE denote the envelope function and the carrier-envelope phase (CEP), respectively. Figure 13(a) shows the temporal evolution of the current J x (t) for the case of a monocycle Gaussian 10 THz, 100 fs pulse with a field amplitude E 0 of 316 kV cm −1 and φ CE = π 2 whose vector potential is shown in figure 13(b). We set T = 0, µ = 0 and v F = c/300 (γ = 3.09 eV) in this section, unless otherwise stated. Since the MDF and TB models give virtually the same results (see figure 16(a) below), the former model has been used for figures 13-15. In the beginning of the pulse, the generated current basically follows the change of the vector potential with some nonlinearity. On the other hand, in the latter part of the pulse (t 0), the curve for the current is largely distorted, compared with the vector potential and the current remains finite and virtually constant after the pulse. This indicates that direct current can be generated in the bulk part of graphene by monocycle terahertz pulses. The contributions from the intraband current and interband polarization are also plotted in figure 13(a). One can see that the dc component is almost entirely due to the former, which implies that carriers remain after the pulse, asymmetrically distributed in momentum space. Since the interband polarization gives only a very small fluctuation, we plot only the intraband current in figures 14 and 15 below.
In order to clarify the mechanism of the carrier generation, let us examine how the variation of Fermi energy µ and pulse intensity affect dc generation. Figure 14(a) presents the CEP dependence of the generated current that remains after the pulse with a field amplitude of 316 kV cm −1 and for two different values of µ = 0 and 0.3 eV. Since the photon energy (0.041 eV) is much smaller than 0.3 eV, the interband transition through a single-or few-photons absorption is suppressed for the case of µ = 0.3 eV. Nevertheless, basically a similar amount of current is generated for the two cases, though the generated current oscillates with CEP more largely for µ = 0 eV than for 0.3 eV. Therefore, the carrier generation is not due to the singleor few-photons jumping between the bands. Another possible mechanism of carrier generation is the enhanced interband transition near the Dirac point. For this to take place, the displacement |eA(t)| in the momentum space must be greater than the momentum p µ ≡ |µ/v F | corresponding to the Fermi energy. In figure 14(b) we show a comparison for the case of a field amplitude of 100 kV cm −1 . One finds that dc generation is almost completely suppressed for µ = 0.3 eV. It should be noticed that p µ = 2.41 × 10 −2 a.u. for µ = 0.3 eV is smaller than e A 0 = 4.04 × 10 −2 a.u. for E 0 = 316 kV cm −1 but larger than 1.28 × 10 −2 a.u. for 100 kV cm −1 . Thus, we can identify the enhanced interband transition as the origin of the carrier generation leading to the dc current. Indeed, the direction of the remaining current can be understood as follows. Carrier electrons in the conduction band are most efficiently generated when the displacement in the momentum space reaches its maximum, i.e. around t = 0. Some of the valence electrons originally with negative p x are lifted to the conduction band in the positive p x region at this time and then brought back to the negative p x region still in the conduction band, which leads to a positive remaining current, as can be confirmed in figure 13(a).
In figure 14(a) we also plot the result for a three-cycle pulse (blue line). One can see that the remaining current is substantially reduced compared with the case of a monocycle pulse. Hence, the dc generation is a feature unique to ultrashort pulses in the monocycle regime.
Let us now look closely at CEP dependence. In figure 14 we have already seen that the dependence of the dc current J dc remaining after the end of the pulse on φ CE is mainly composed of a large, slow oscillation with a period of 2π , superposed by a faster oscillation. Figure 15(a) shows how this oscillation varies with pulse intensity. Interestingly, the more intense the pulse, the higher the modulation frequency. As a consequence, remarkably, the value of J dc at a fixed carrier envelope phase is a complicated function of intensity (or equivalently, field amplitude). The dependence is clearly nonlinear, but not of simple integer orders. At φ CE = 45 • , for example, J dc jumps from 1.3 × 10 −6 to 1.2 × 10 −5 when the peak field amplitude increases from 100 to 200 kV cm −1 , while it is nearly unchanged even when the field amplitude further increases to 316 kV cm −1 . Moreover, at φ CE = 10 • , for example, J dc does not vary monotonically with the intensity, even changing its sign.
It follows from the space inversion symmetry and the symmetry between c + and c − in equation (11) that the generated direct current J dc 1 and J dc 2 for the vector potential A 1 (t) and A 2 (t), respectively, satisfies For the vector potential given by equation (59) with a Gaussian envelope A 0 (t), these properties translate into which can be confirmed in figures 14 and 15(a). As a consequence, the CEP dependence J dc (φ CE ) can be expanded in a Fourier sine series as where the spectral amplitudes f 2n+1 are real. We show the spectrum ( f 2n+1 as a function of spectral order 2n + 1) in figure 15(b). The spectrum is mainly composed of the fundamental component (first order) and the modulation component (approximately 3rd, 7th and 11th/13th orders, at 100, 200 and 316 kV cm −1 , respectively) that shifts to higher orders with an increase in intensity. The behavior in figure 15(b) emphasizes the complex nature of the intensity and CEP dependence of dc generation. J dc (φ CE ) vanishes at φ CE = 0 and π. In these cases, as plotted in figure 13(b), the temporal profile of the vector potential A(t) is symmetric to the positive and negative directions. This indicates that one of the origins of the dc generation is a positive-negative asymmetry of the vector potential, i.e. that of the displacement of electrons in the momentum space. Indeed, the asymmetry is largest at φ CE = π 2 and 3 2 π , for which |J dc (φ CE )| reaches maximum. It should be, however, emphasized that the interband transition enhanced in the vicinity of the Dirac points is another essential factor, since, if this is suppressed by a large Fermi energy, there is no appreciable remaining current, as we have seen in figure 14(b).
All the figures presented in this section so far have been obtained within the MDF picture. The results are virtually unchanged even if we use the TB model. Let us now examine how MDF results deviate from TB results if we further increase the pulse intensity. Figure 16 shows the comparison for a field amplitude E 0 from 316 kV cm −1 to 10 MV cm −1 . It is to reemphasize that the pulse acts as an agent to displace electrons in the momentum space and that the interband transition, leading to carrier generation, predominantly takes place when electrons pass sufficiently close to Dirac points. Hence, we should analyze this figure in terms of peak vector potential amplitude A max rather than intensity and field strength. The maximum momentum displacement multiplied with a is e A max a = 0.108, 0.343, 1.08 and 3.43 for E 0 = 0.316, 1, 3.16 and 10 MeV, respectively, which should be compared with the value 4π 3 √ 3 = 2.42 corresponding to the distance between the Dirac points K and K . Figure 17 illustrates the dependence of ∞ −∞ | j c,x (t)| 2 dt on the initial momentum. In the case of E 0 = 316 kV cm −1 and 1 MV cm −1 (figures 17(a) and (b), respectively), carrier and current generation is well localized near the Dirac points and the dynamics in each valley are separated from each other. Thus, the MDF picture is a good approximation to the TB model, as can be seen in figures 16(a) and (b). For E 0 = 3.16 MV cm −1 ( figure 17(c)), the electron dynamics begin to be delocalized, leading to a quantitative difference between the TB and MDF results ( figure 16(c)). At E 0 = 3.16 MV cm −1 , the vector potential has so large an amplitude that electrons can move between the two Dirac cones ( figure 17(d)). The MDF model fails to describe such a situation, leading to even a qualitative discrepancy between the two models, as is clearly seen in figure 16(d). We can also see that dc generation is reduced. Surprisingly, however, even in this case, the current depends on pulse polarization ξ only moderately ( figure 16(d)). The vertical width of each distribution in figure 17 can reasonably be predicted by the criterion equation (58) as ∼0.044, 0.078, 0.14 and 0.25, respectively.

Conclusions
Based on both the MDF and the TB models, we have derived equations that describe the coherent population dynamics of electrons in graphene irradiated by a terahertz radiation pulse of arbitrary wave form, angle of incidence and polarization. These equations are equivalent to the time-dependent Schrödinger equation and, at the same time, provide a physically transparent description of intraband dynamics and interband transitions and polarization. For the case of normal incidence, these can be further cast into the form of GBE (12) and (13). We have also derived the formula for a single-electron current, the integration of which over the momentum space leads to a macroscopic current.
The calculations using the derived equations have revealed a sensitive dependence of the single-electron dynamics on the model, pulse polarization with respect to the lattice and the electron's initial momentum. This dependence is, however, mostly averaged out after the momentum-space integration and thus the macroscopic current exhibits only moderate dependence.
With the help of our formulation, we have obtained exact analytical solutions for the dynamics of a single electron whose kinetic momentum takes a circular path around the Dirac point. The solutions account for Berry's phase in the adiabatic limit and predict full population oscillation in the limit of rapid circulation.
Furthermore, we have shown that a direct current is generated when graphene is subject to a monocycle intense THz pulse, due to enhanced carrier generation taking place when electrons pass near the Dirac points. from equation (14), where φ denotes the initial value of θ(t) satisfying p x = p cos φ and p y = p sin φ and also (t) ≈ ω 0 t, (A.2) with ω 0 = v F p/h. By substituting equations (A.1) and (A.2) to equations (12) and (13) we obtain, using equation (21), (A.7) Now that we have a linear relation between j x and A x , it is convenient to redefine E x (t) and A x (t) as Then, equation (A.7) is rewritten as (A.9) Let us now integrate this with respect to p, assuming intrinsic graphene (µ = 0). The total integrated current equation (20) is determined by the pole at ω = 2ω 0 , i.e. resonant interband transition, leading to J x (t) = e 2 4h E x (t), (A.10) which indicates the universal conductivity.