Wave packet dynamics in a monolayer graphene

The dynamics of charge particles described by Gaussian wave packet in monolayer graphene is studied analytically and numerically. We demonstrate that the shape of wave packet at arbitrary time depends on correlation between the initial electron amplitudes $\psi_1(\vec r,0)$ and $\psi_2(\vec r,0)$ on the sublattices $A$ and $B$ correspondingly (i.e. pseudospin polarization). For the transverse pseudospin polarization the motion of the center of wave packet occurs in the direction perpendicular to the average momentum $ {\vec p_0}=\hbar \vec{k_0}$. Moreover, in this case the initial wave packet splits into two parts moving with opposite velocities along $ {\vec p_0}$. If the initial direction of pseudospin coincides with average momentum the splitting is absent and the center of wave packet is displaced at $t>0$ along the same direction. The results of our calculations show that all types of motion experience {\it zitterbewegung}. Besides, depending on initial polarization the velocity of the packet center may have the constant component $v_c=uf(a)$, where $u\approx 10^8 cm/s$ is the Fermi velocity and $f(a)$ is a function of the parameter $a=k_0d$ ($d$ is the initial width of wave packet). As a result, the direction of the packet motion is determined not only by the orientation of the average momentum, but mainly by the phase difference between the up- and low- components of the wave functions. Similar peculiarities of the dynamics of 2D electron wave packet connected with initial spin polarization should take place in the semiconductor quantum well under the influence of the Rashba spin-orbit coupling.


I. INTRODUCTION
At last years the dynamics of wave packets in 2D electron gas and other systems in solids including the phenomenon of zitterbewegung (ZB) or trembling motion has been the subject of numerous studies. 1,2,3,4,5,6,7,8,9,10,11,12,13 Firstly the oscillatory motion analogous to the relativistic zitterbewegung in twodimensional systems with the structural and bulk inversion asymmetry was investigated by Schliemann et al. 3 In the recent work by authors the detailed studying of the electron wave packet dynamics in the semiconductor quantum well under the influence of the Rashba spinorbit coupling was performed. 14 It was shown (analytically and numerically) that the initial wave packet splits into two parts with different spin polarizations propagating with unequal group velocities. It was demonstrated also that the splitting and overlapping of wave packets leads to the damping of zitterbewegung.
As well known, the electron zitterbewegung in relativistic physics at first time was predicted by Schrödinger 15 (see also 16 ). This phenomenon is caused by the interference between positive and negative energy states in the wave packet. The frequency of ZB motion is determined by the gap between these two states and when the amplitude of oscillations in a particle position is of the order of the Compton wave length. This phenomenon was discussed also in Refs. [17][18][19].
In the papers by Rusin and Zawadzki 12,13 the evolution of the wave packet in a monolayer and bilayer graphene as well as in carbon nanotubes was analyzed. The exact analytical expressions for two components of wave function and average value of position operator were found for bilayer graphene, which allowed to obtain analytical results for the ZB of Gaussian wave packet. It was shown that the transient character ZB in bilayer graphene is due to the fact that wave subpackets related to positive and negative electron energies move in opposite directions, so their overlap diminishes with time. At the same time the dynamics of the wave packets in a monolayer graphene in Ref. [12] was not investigated fully.
In this work we present the detailed description of wave packet evolution including the phenomenon of ZB of the packet center in a monolayer graphene. We investigate the influence of the initial pseudospin polarization on the dynamical characteristics of the moving electron wave packet. Our analytically results are illustrated by a graphic presentation.
The outline of this paper is as follows. The expression for the wave function at t > 0 is found in Sec.II for the arbitrary initial pseudospin polarization. The time dynamics of Gaussian wave packet and, in particular, the ZB phenomenon are described in Sec.III for different correlations between the initial electron amplitudes on the sublattices A and B. We conclude with remarks in Sec.IV.

II. BASIC EQUATIONS
Graphene is a single layer of carbon atom densely packed in a honeycomb lattice. The two-dimensional Hamiltonian describing its band structure has the form 20,21,22,23Ĥ where u ≈ 10 8 cm/s,ˆ p = (p x ,p y ) is the momentum operator defined with respect to the centre of the valley centered at the corner of the Brillouin zone with wave vector K. Pauli matrices σ i operate in the space of the electron amplitude on two sites (A and B) in the unit cell of a hexagonal crystal. This internal degree of freedom plays a role of a pseudospin. The Dirac-like HamiltonĤ determines the linear dispersion relation Here p = p 2 x + p 2 y , s = 1 for the electron in the conduction band and s = −1 for the valence band ("hole" branch of quasiparticles). The corresponding eigenfunctions are given by with e iϕ = px+ipy p . The time-evolution of an arbitrary initial state ψ( r, 0) in Shrödinger representation can be found with the help of Green's function G µν ( r, r ′ ) where µ, ν = 1, 2 are matrix indices, corresponding to the upper and lower components of ψ( r, t). These components are related to the probability of finding electron at the sites of the sublattices A and B correspondingly. The standard expression for Green's function is G µν ( r, r ′ , t) = s=±1 d p ϕ p,s;µ ( r, t)ϕ * p,s;ν ( r ′ , 0). (5) Using Eq.(3) for ϕ p,s;µ ( r, t) we find Let us represent the initial wave function by Gaussian wave packet having the width d and nonvanishing average momentum p 0y =hk 0 where coefficients c 1 and c 2 determine the initial pseudospin polarization. We suppose that the packet width d is much greater than the lattice period and consequently ψ( r, 0) is smooth enveloping function. We suppose also that the most of the states in valence band are unfilled, that corresponds to negative Fermi level located far from Dirac point (see also 24 ). Substituting Eqs.(8a, 8b) in Eq.(4) and using the expressions (6) and (7) we obtain where, for notational convenience, φ 1,2 ( r, t) denote the functions Using the cylindrical coordinates in Eqs. (11), (12) and integrating over the angular variable, we have x + a + iy where J 0 (z), J 1 (z) are Bessel functions. For the sake of convenience we introduce in Eqs. (13), (14) and everywhere below the dimensionless variables, measuring the distance in the units of initial width of wave packet d and time in d/u units. Besides, instead of the wave vector k 0 we consider the parameter a = k 0 d.

III. ZITTERBEWEGUNG OF GAUSSIAN WAVE PACKET WITH DIFFERENT PSEUDOSPIN POLARIZATION
Now we describe the time dynamics of Gaussian wave packets, in particular, the ZB phenomenon and the influence of the initial pseudospin polarization on the characteristics of trembling motion. i). Following Ref. [12] let us firstly consider the model problem when the lower component of initial wave function is equal to zero, i.e. the parameters c 1 = 1, c 2 = 0 in Eq.(8a). That means that at the initial moment of time the electron probability is located at the sites of the sublattice A. It is not difficult to show that this packet is formed by the states with positive and negative energies. The relative weight of these states is equal to one. The wave function for t > 0 can be found using Eqs.(9), (10): where the functions φ 1 ( r, t), φ 2 ( r, t) are defined by Eqs. (13), (14). In Fig.1 we represent the full electron density at the moment t = 7 for initial wave packet, Eq.(8b) with width d = 2 nm and k 0 = 0.6 nm −1 . As one can see, at t > 0 this packet splits in two parts moving along y axis with opposite velocities so that the electron probability density is symmetrical with respect to y: ρ(x, y, t) = ρ(x, −y, t) (note that at the case k 0 = 0 the electron probability density has a cylindrical symmetry at all time). For enough large time the width of both parts of the packet increases with time due to effect of dispersion. One can check that in this situation the contributions of two components of wave functions ψ 1 ( r, t) and ψ 2 ( r, t) in full electron density are equal. In other words the electron probability distributes with the time on the sides of sublattice A and B. Note at the same time ρ(x, y, t) = ρ(−x, y, t) and the packet center oscillates along x direction (zitterbewegung).
To analyze this motion we find the average value of position operator. To do it, we use the momentum representation. The upper (C 1 ( p, t)) and lower (C 2 ( p, t)) components of wave function (15) in this representation can be easily obtained from Eqs. (11), (12). After that the usual definition readily yieldsȳ where I 1 (z) is a modified Bessel function of the first order. In Eq.(17b) the integral term represents the zitterbewegung. Note that average valuex(t) depends on only one parameter a (in the dimensionless variables). The obtained functionsx(t) which describes the typical transient zitterbewegung are plotted in Fig.2. After the oscillation disappears the center of the packet is displaced by amount which equals to the first term Eq. (15). In the case when the wave packet width is large enough and the inequality a = dk 0 ≫ 1 takes place, Eq.(17b) reduces to 25x As it follows from Eqs. (17), (18) for given initial polarization of wave packet the ZB occurs in the direction perpendicular to the initial momentum p 0y =hk 0 , just as for bilayer graphene 12 and for the semiconductor quantum well in the presence of the Rashba spin-orbit coupling 3,14 . One can see from Eq.(18) that the trembling motion has a transient character as it was described in Refs. [12,14] and at t ≫ 1 x(t) → 1/2a. We should notice that Eqs.(17b), (18) coincide with corresponding formulas of Ref. [14]. This is because the Hamiltonian for the system with Rashba-coupling where α is a Rashba coupling constant, transforms into Hamiltonian for monolayer graphene, Eq.(1), if we make the replacement in Eq. (19) x ii). Let us consider now the case when c 1 = c 2 = 1, that is pseudospin is directed along x axis at t = 0. Then from Eqs.(9), (10) Fig.3 illustrates the corespondent electron probability density at the time moment t = 7 for initial wave packet, Eq.(8b), for the same parameters as in Fig.1. One can see that the initial wave packet at t > 0, as in previous case, splits into two parts propagating along y in opposite directions so that the symmetry concerning this axis, i.e. ρ(x, y, t) = ρ(x, −y, t), retain during the time (as the case i)). The distribution of the probability density along x axis clearly demonstrates that its maximum is displaced in the positive direction that corresponds to the motion of the packet centre along x axis. The velocity of such motionv x = dx dt consists of both constant as well as oscillatory parts. Really, a straightforward calculation yields the average value of position operator x andȳ(t) = 0 like for the case i). In Fig.4 we demonstrate the dependencex(t) for various values of parameter a. When the parameter a increases, the amplitude of ZB and the period of oscillations decrease. At a ≫ 1 we have from Eq.(22) We see that the character of motion of wave packet is changed. Now the center of wave packet moves along x direction with constant velocity, which is determined by the first term in Eqs. (22), (23) and performs the damping oscillations. This constant velocity has a maximum value which is equal to 1/2 (in the units of u) for the narrow -width wave packet, or for the case k 0 = 0. When the width of packet is increased the velocity of motion of its centre is decreased. The frequency and amplitude of the zitterbewegung for a ≫ 1 are the same as in the case i). However, the first term in Eq. (22) corresponding to the motion of wave packet with constant velocity reduces the effect of ZB at least for a < ∼ 1 (Fig.4). It is not difficult to show that as in the other two-band systems the phenomenon of ZB in graphene is a result of an interference of states corresponding to positive and negative eigenenergies of Hamiltonian, Eq.(1). For wide enough packet a = k 0 d ≫ 1 and at time t > 1 when the ZB disappears two split parts of initial wave packet (see Fig. 3) move along y axis with opposite velocities u and −u. In this situation the subpackets moving in the positive and negative directions consist of the states with positive and negative energies correspondingly.
iii). When the initial pseudospin is along y axis the wave function at t > 0 has the form In Fig.5 the full electron density for the same moment of time and for the same parameters as in previous cases is shown. As one can see, the initial wave packet does not split into two parts at t > 0 unlike in the cases i) and ii). This result is confirmed by the straightforward calculations. Indeed, one can show that the eigenenergy states corresponding to propagation in the positive direction along y axis give the dominant contribution in total wave function, Eq. (24). For wide packets (a ≫ 1) almost all of these states belong to the positive branch of energy.
The results of calculations of average values of x and y for this polarization lead tō Thus in the considering case the wave packet propagates along y axis and the zitterbewegung has the "longitudinal" character 26 . Note that in a numerical work 18 by Thaller the authors have observed similar oscillatory behavior of the expectation value of the position operator for one -dimensional relativistic electron in vacuum. In Ref. [18] it was also shown that apart from the rapid oscillation, the wave packet drifts slowly even when its average momentum is zero.
In Fig.6 we represent the dependenceȳ(t) for different values of parameter a. As one can see, even at zero value of a the oscillations are absent. In this case, as it follows from Eq.(26), the drift velocity is equal to 1/2 (in the units of u). As above, Eq.(26), takes more simple form at a ≫ 1ȳ Comparing Eqs. (18), (23), (27), we see that the amplitude for the "longitudinal" zitterbewegung is much smaller than the amplitude of "transverse" zitterbewegung. This fact can bee seen as a consequence of special form of the initial wave function, which in the given case consists of (at a ≫ 1) the states with positive energy mostly. That makes the interference between the positive and negative components difficult, i.e. decreases the ZB. Moreover, at any values of the parameter a the integral term in Eq.(26), corresponding to the oscillating motion, is negligible in comparison with the first term, and one may neglect the effect of the "longitudinal" ZB.
As was demonstrated above, the direction of the average velocity depends not only on module of the components ψ 1 ( r, 0) and ψ 2 ( r, 0), but also on their phases. In general case for the initial Gaussian packet the probability density becomes asymmetric and the average position operator has two components where ϕ is an arbitrary phase difference between the up and low components of wave function andx(t),ȳ(t) are determined by Eqs. (22), (26). For illustration we show in Fig.7 the electron probability density obtained for the initial packet, Eq.(28), with ϕ = π/4. It is clear that the phase ϕ determines the direction of the average velocity of the packet center. Using the expression for velocity operatorˆ v = u σ and Eq.(28) we obtain (in the dimensionless variables) at t = 0: v x (0) = cos ϕ, v y (0) = sin ϕ. (30) The components of the velocity for a large enough time, when the trembling motion stops, can be found from Eqs. (22),(26) and (29) for arbitrary parameter a (31) In particular, as it follows from Eq.(31) for a ≪ 1, the direction of the motion of wave packet center at large time coincides with the initial one, Eq(30). In other limiting case a ≫ 1 (and for not too small ϕ) asymptotic direction of the average velocity is along 0Y axis, i.e. along the average momentum of wave packet p y =hk 0 . Thus, by changing the initial phase ϕ, one can govern the packet motion and consequently the direction of dc current. To illustrate this, we plot in Fig.8 the packet center trajectories for two initial phases: ϕ = π/4 and 3π/4. Note that the packet motion with the constant velocity predicted above (see Eqs. (22), (26)) should lead to the existence of the direct current in the corresponding direction.
To check our formalism let us consider the plane wave as the starting point. In this case it is easy to obtain the expression for the average value of electron velocity¯ v(t). Really, in the Heisenberg picture the kinetic velocity is (in the dimensional variables) where In these equations p(t) = p(0). Let the initial momentum p oy =hk 0 . Then, using the solutions of Eqs.(32), where ω = 2uk 0 and σ i (0) = σ i -Pauli matrixes (i = 1, 2, 3). So, if in the initial state pseudospin is along z direction, i.e.σ z (0) = 1 (case i)) we obtain from Eq.(34a) thatv x (t) = u sin ωt which leads tō Returning to the original variables in Eq. (18) and setting d = ∞ we see that this expression coincides with Eq.(35). We get similar results also for other initial polarizations.

IV. CONCLUDING REMARKS
We have studied the quantum dynamics of charge particles represented by Gaussian wave packets in twodimensional single layer of carbon atoms (graphene). We investigated numerically also the spatial evolution of the initial wave packet and demonstrated the effect of the packet splitting for the pseudospin polarization perpendicular to the average momentum. The analytical expressions for the average values of position operators were obtained. These expressions clearly demonstrate that the evolution of wave function is accompanied by the zitterbewegung and strongly depends on the initial pseudospin polarization. In particular, if the initial pseudospin polarization coincides with initial average momentum, the packet center moves and exhibits the ZB along the same direction. In this case the second term in Eq.(26) describing the longitudinal oscillations (the "longitudinal" ZB) is essentially smaller than the first one connected with the motion with constant velocity. As for other systems with two-band structure 12,14,18 , the ZB in monolayer graphene has a transient character.
It was also predicted that apart from the packet center exhibits the trembling motion it can move with constant velocity (for the pseudospin polarization along x and y axis). The direction of this velocity depends on not only the orientation of average momentum p 0 , but also on the module of the components ψ 1 ( r, 0), ψ 2 ( r, 0) and the differences of their phases (see Eqs. (28),(31)).
All above calculations have been done for the K point of the Brillouin zone in graphene. Similar results can be found for initial wave packet with wave vector k in the valley centered in inequivalent point K ′ . The Dirac Hamiltonian around K ′ point can be written as 27 This expression can be obtained from Hamiltronian around K point given by Eq.(1) by replacementp x → −p x . Thus valuesx(t) for the wave packet of different polarizations (and corresponding components of velocity) change sign whileȳ(t) remain unchanged (see also 24 ).
In conclusion we would like to stress that the packet motion with the constant velocity (see Eqs. (22), (26)) leads to the appearance of the dc current. For the experimental detection of this current one needs sensitive current meters. Experimental observation of trembling motion is currently a more difficult task since it is necessary to use femtosecond techniques. 6,12 If new methods of formation of wave packets with different pseudospin polarizations will be elaborated then their trajectories and spatial separations can be observed probably with the help of devices with quantum point contacts and gates (see for example 28 ). i , where F (y) is arbitrary function, propagates along y direction without changing of its form ψ(y, t) = F (y − ut) 1 i .