Revivals of quantum wave-packets in graphene

We investigate the propagation of wave-packets on graphene in a perpendicular magnetic field and the appearance of collapses and revivals in the time-evolution of an initially localised wave-packet. The wave-packet evolution in graphene differs drastically from the one in an electron gas and shows a rich revival structure similar to the dynamics of highly excited Rydberg states. We present a novel numerical wave-packet propagation scheme in order to solve the effective single-particle Dirac-Hamiltonian of graphene and show how the collapse and revival dynamics is affected by the presence of disorder. Our effective numerical method is of general interest for the solution of the Dirac equation in the presence of potentials and magnetic fields.


Introduction
Since its experimental discovery, graphene has become a hot topic in solid state physics [1,2]. On the theoretical side, the similarities between the effective singleparticle Hamiltonian for graphene in the low-energy regime and the Dirac equation for massless particles have been noted [3,4,5] and lead to the prediction of "relativistic effects" in graphene. One example is the non-transient zitterbewegung in the presence of an external magnetic field [6,7]. In order to study the zitterbewegung and other transport effects in graphene one needs to solve the time-dependent Dirac equation with high accuracy and efficiency. Previous time-dependent propagation methods have often been restricted to eigenstate decompositions. However, the evaluation of integrals involving a finite set of eigenfunctions is cumbersome and also limited to energetically low lying eigenstates. Additionally the eigenstate approach can only be used if the analytic solutions are known, since the numerical determination of a full set of energy eigenstates is in general not feasible.
In section 2 we show that the prior knowledge of the eigenstates is not required in order to obtain a reliable solution of the time-dependent Dirac-Hamiltonian of graphene. Our algorithm is accurate up to the machine error and stable for any propagation time. In addition it is very flexible and thus allows us to calculate the time-evolution in inhomogeneous magnetic fields and arbitrary potential landscapes. In section 3 we study quantitatively the scattering in realistic potentials and obtain the fidelity (which is related to the mobility) of graphene in a magnetic field. For wave-packets with large initial momenta we reveal a rich collapse and revival structure, which has a strong similarity to the dynamics of highly excited states of Rydberg atoms [8,9]. We derive analytic results for the semiclassical expectation value of the position operator and the autocorrelation function which are in very good agreement with the numerical results. Subsequently, we study the effect of an impurity potential and analyse the suppression of the autocorrelation function corresponding to a rapid fidelity decay. Interestingly, revivals in the centre-of-mass coordinate survive the disorder-perturbations if the correlation length of the potential is on the order of the cyclotron radii of the occupied states. In general, our algorithm provides a basic method to model time-dependent transport through mesoscopic systems [10,11,12].

Wave-packet propagation for the graphene-system
The electronic structure of graphene is obtained from the tight binding model of the two-dimensional honeycomb lattice. In the resulting band structure two inequivalent vanishing bandgaps (often referred to as different valleys or Dirac points) appear, which are located at the corners of the Brillouin zone [13,14]. The band structure around each Dirac point is conveniently described in the framework of an effective k · p perturbation theory and the resulting Hamiltonian [15] has the same structure as the 2 + 1 dimensional Dirac Hamiltonian, where the velocity is given by v ≈ 10 6 m/s, and σ x , σ y denote the corresponding Pauli matrices [16]. The linear approximation is valid for energies close to the Dirac points. For higher energies the trigonal warping gets important [17]. We neglect this effect here, since we only consider wave packets which can be decomposed into eigenstates from −0.7 eV to 0.7 eV. The Hamiltonians (1) at the two Dirac points are obtained by switching the sign of τ = ±1. Here, we study the bulk region of a sheet of graphene and therefore we can neglect effects of the boundaries which may couple the valleys and induce symmetry breaking [18,19]. In the following we present a new algorithm for the propagation of a wave-packet on a sheet of graphene in an anisotropic magnetic field B z (r) oriented perpendicularly to the sheet. The Hamiltonian can also contain a position-dependent potential landscape V (r). The resulting system is described by the following Hamiltonian For simplicity we only use the Hamiltonian with τ = 1, but it is straightforward to apply the same algorithm to the second Dirac point. The vector potential of the magnetic field A(r) = (A x (r), A y (r), 0) is connected to the magnetic field B(r) = ∇ × A(r) and enters the Hamiltonian asÃ(r) = A x (r) − iA y (r). As initial state of the system we choose an arbitrary two spinor wavefunction ψ(t = 0) = {ψ A (t = 0), ψ B (t = 0)}. The action of the time-evolution operator on the wavefunction is given by In principle, there exist several possibilities to calculate the action of the time-evolution operator. For a system with a discrete set of eigenstates defined byĤ|φ n = E n |φ n a favourable way is given by the expansion in eigenstates. In this case, the time-evolution becomes Such a propagation algorithm was already used in references [6,7] to evaluate the zitterbewegung in graphene and also for other semiconductor materials [20]. Due to the lack of analytical known solutions of the Dirac equation, this scheme can only be applied to a limited set of problems. Thus it is favourable to develop an algorithm which does not require any knowledge of the eigenstates and can be applied to arbitrary mesoscopic systems. A similar starting point exists in quantum chemistry, where no analytic solutions are known for the molecular potential energy surfaces. In all these cases it is favourable to use an algorithm which directly gives the time-evolution of the system without evaluating all eigenstates first [21]. All relevant information about stationary properties, i.e. the energy spectrum or the Green's function and the local density of states, is encoded in time-dependent observables, which are extracted afterwards from the time-evolved wave-packet [22,10].
For the time-dependent solution of the Dirac Hamiltonian, we use a reliable and efficient algorithm based upon the polynomial expansion of the propagator in Chebyshev polynomials [23]. The main idea of the algorithm is to expand the time-evolution operator on the interval [−1, 1] using a set of Chebyshev polynomials, which are particularly suitable since they distribute the numerical error of the expansion equally on the whole interval. As a first step the eigenvalues of the HamiltonianĤ have to be projected into the interval of the expansion by introducing a normalised Hamiltonian The scaling-energy ∆E has to be big enough to cover all eigenenergies contributing to the initial wave-packet. With the aid of (5) the time-evolution is given by where T n stands for the series of Chebyshev polynomials of the first kind. The expansion coefficients are given by a n ∆Et 2h with the Bessel function J n . Only the action of the polynomial ofĤ norm on the initial state is needed to calculate the time-evolution. Thus it is sufficient to recursively calculate the vectors P n by the following modified recursion relations for Chebyshev polynomials: P n = − 2iĤ norm P n−1 + P n−2 .
In order to avoid the fermion doubling problem of a simple first-order nearest-neighbour derivative [24] we use a Fourier representation of the spinors and apply the kinetic energy operator in momentum space [21]. Another possibility is to introduce certain non-local first-order derivatives [25]. In our case, the action of the normalised Hamiltonian acting on a two spinor wavefunction (ψ A , ψ B ) is given bŷ where F −1 and F stand for the two-dimensional Fast Fourier Transform. This type of propagation algorithm allows us to evaluate the time-evolution of the wave-packet for arbitrary long propagation times, as long as the wavefunction does not leave the region of interest. A good measure of the propagation quality is given by the comparison of the analytic autocorrelation function for the homogeneous magnetic field with the numerical result. Our results show that the maximum error is indeed on the order of the unavoidable finite machine precision. In contrast to all other propagation algorithms for graphene we are able to incorporate arbitrarily shaped potentials and also inhomogeneous magnetic fields in our configurations.

Model and observables
In the following we apply the propagation algorithm to an infinite sheet of graphene put in a perpendicular magnetic field B(r) = (0, 0, B) and set the potential V (r) to zero. This system was proposed by Rusin et. al. [6] and Schliemann [7] to study the non-transient zitterbewegung of a wave-packet on graphene. The eigenstates of the Dirac Hamiltonian in a perpendicular magnetic field are known and given by the harmonic oscillator eigenfunctions [26]. In the lowest Landau level, the ground state wavefunction has Gaussian shape and is located on one sublattice. As initial condition we consider a kicked Gaussian of the form The width is taken as a 0 =h/(eB) and the initial momentum k can point along any direction in the x − y plane, but is set to k = k 0 e y in the following. The pseudospin σk of the wave-packet is determined by the phase-relation between the two spinor entries. For ϕ = π/2, the wave-packet is mostly electron-like since the momentum is pointing along the y-direction. Consequently the wave-packet is mostly hole-like if ϕ = 3π/2 and a mixture of both states for intermediate angles.
During the propagation we track several observables, including the expectation value of the position operator We refer to r(t) as the "centre-of-mass" although the calculated motion describes massless particles. The trembling motion of the centre-of-mass is visible in the videos linked from figure 1(a) for different initial setups. Another important property, the local pseudospin, is given by ψ(r, t)| σk|ψ(r, t) at each gridpoint r and colour coded in the video. The local pseudospin reveals the local electron or hole character of the wave-packet and shows an anti-symmetry with respect to the x-axis if the pseudospin of the initial state is oriented orthogonal to the initial momentum. An advantage of our method is the extraction of stationary properties, like the local density of states (LDOS) and the Green's function, out of the time-evolution [10]. For example the relative strength of each energy eigenstate with respect to the initial state is given by the Fourier transform of the autocorrelation function The same concept allows us to calculate the E × B drift motion and eigenstates in the quantum Hall regime of graphene for high bias currents [26]. The LDOS is given by the Laplace transformation of the autocorrelation function As shown in figure 1(b) the LDOS of the trembling wave-packet reveals pronounced peaks at the energies of the Landau levels E n = v sign(n) 2eBh|n|. The peak height (which is proportional to the enclosed area in the LDOS) indicates the overlap of the initial Gaussian wave-packet with the different Hermite polynomials of the Landau levels. For ϕ = π/2, mostly electrons are occupied whereas for ϕ = π holes and electrons are equally occupied. Our computational scheme derives the overlap information without calculating the numerically unstable integrals over Hermite polynomials used in previous methods. Therefore we can increase the initial momentum k 0 of the wavepacket and study the quantum mechanical propagation of high Landau levels.
The trembling motion on graphene is present since both, electron and hole-like states, contribute to the initial wave-packet. For high initial momenta the index of the average Landau level is given bȳ which follows from the semiclassical quantisation condition (A.16). Our quantum mechanical calculations support this assertion as shown in the LDOS of figure 2. A strongly kicked Gaussian wave-packet mostly occupies one energy branch if the initial pseudospin is parallel to the initial momentum. For an initial momentum which puts the Gaussian wave-packet into the 24.5 th Landau level, the hole contributions are 400 times smaller compared to the electron parts. In such a scenario the trembling motion may be neglected and other interesting revival phenomena appear. In atomic physics revival effects are commonly observed for valence electrons in the Coulomb potential of Rydberg atoms and have been studied in theory [8] as well as in experiments [9,27]. Reference [28] contains a recent review about revivals of quantum wave-packets. In order to get an idea about the rich revival structure of a quantum wave-packet on graphene, it is instructing to look at the numerical results in figure 3 and the corresponding video. Wave-packet revivals are typically not present in semiconducting material in a perpendicular magnetic field since there the dynamics is governed by a purely quadratic Hamiltonian. Graphene and toplogical insulators form their own class due to their non-quadratic Hamilton operators.

Revivals in the centre-of-mass position
As first phenomena we study the collapse and the revival of the centre-of-mass of a graphene wave-packet. The quantum mechanical position of the centre-of-mass cannot be found analytically. Thus we introduce a classical picture where the centre-of-mass is given by a weighted sum over the centre-of-masses moving along quantised cyclotron orbits. The centre-of-mass of the classical cyclotron motion in the n-th Landau level is defined by as derived in Appendix A. The occupation P n depends on the initial momentum k 0 , the initial pseudospin and the width a 0 of the wave-packet. As an approximation we find . Revivals and fractional revivals of a cyclotron wave-packet in a magnetic field of B = 1 T. The red pin stands for the quantum mechanical centre-of-mass, whereas the green pin is calculated with the semiclassical expression (18). The initial wave-packet in (a) was chosen such that the contributing eigenstates follow the Gaussian distribution (17) withn = 60, σ = 1, leading to a Poincaré cyclotron time of t0 = T cl . During the time-evolution the centre-of-mass collapses to the centre of a circular wavefunction and later shows pronounced revivals. For example (b) displays a quarter revival at t = 1/2n T cl , (c) a third revival at t = 2/3n T cl and (d) a half revival at t =n T cl . All revivals have specific symmetries which also occur in the subsequent revivals. After a certain time a localised wave-packet emerges in (e) as a mirror revival at t = 2n T cl . In this case, the quantum mechanical centre-of-mass is at the opposite side compared to the semiclassical centre-of-mass. This phenomenon has already been seen in Rydberg revivals [29]. (f) The full revival where the semiclassical and the quantum mechanical centre-of-mass coincide occurs for t = 4n T cl . The related movie (movie2.mp4, 3.0 MB) which illustrates the dynamics is available from http://www.quantumdynamics.de/graphene. from numerical calculations for high average Landau levels a Gaussian distribution of the form shown in figure 2. Heren stands for the average occupied Landau level and σ denotes the Landau level spread. The sum over the Landau levels is normalised ifn ≫ σ is fulfilled.   (20). A small deviation is visible which stems from the fact that the centre-of-mass of the initial wave-packet has already a strong angular dispersion (see figure 3) and thus the average cyclotron radius l is overestimated.
In this limit the cyclotron radii l n do not differ substantially from the cyclotron radius of the average occupied Landau level l := ln. Thus, the time-evolution of the average position is approximated by The last sum contains both, the centre-of-mass collapse as well as the revivals. Numerical results of this classical centre-of-mass motion are shown in figure 3. In the corresponding video the classical approximation is marked by the green pin, which almost perfectly coincides with the exact quantum mechanical result indicated by the red pin.

Collapse:
Next, we extract the times of the initial collapse of the centre-of-mass and the revivals out of equation (18). We expand the cyclotron frequency ω n around the average Landau leveln: For high Landau levels the energy spectrum becomes denser and the sum can be approximated by an integral which we evaluate analytically after the series is reduced to the first two terms: Equation (20) shows that the centre-of-mass motion completes one cyclotron orbit during the cyclotron time of the average Landau leveln and also a Gaussian decay on the timescale Our numerical calculations in figure 4 show the same behaviour. Very small deviations from the exact result stem from the initial radial dispersion of the Gaussian wave-packet, which is not included in the classical theory.

Revivals:
The revival of the classical centre-of-mass takes place if all particles orbiting along the contributing Landau orbits reunite at the initial position. We calculate the revival time from the Taylor expansion of the classical centre-of-mass motion inserted in equation (18). In order to meet the conditions for a revival |z(T rev(CM) )| = l, all elements of the sum must have the same phase. Thus, the n-dependent addend of the exponent inω ′n t must be a multiple of 2πi. This leads to a revival time of the classical centre-of-mass given by In this case the classical centre-of-mass of the wave-packet is on the opposite side compared to the centre-of-mass of the quantum mechanical calculation. This phenomenon is visible in the comparison between the classically and the quantum mechanically calculated expectation value of the x coordinate for the first revival time shown in figure 5(a). Exactly the same effect occurs in the mirror revivals of Rydberg atoms [29,30]. For the second revival the classical and the quantum mechanical centreof-mass coincide again as is evident from figure 5(b).

Revivals in the autocorrelation function
In contrast to the centre-of-mass revivals, the revivals in the autocorrelation function can be understood without a semiclassical approximation. By definition the time-evolution of the autocorrelation function is given by  Again, we expand the square root dependence in a Taylor series leading to an approximate autocorrelation function where E ′ n stands for the first derivative of the dispersion and E ′′ n , E ′′′ n denote the higher order derivatives. It is astonishing that all orders are prominently encoded in the quantum mechanical motion and influence the picture at the relevant revival times. From the term iEnt/h in the exponent of equation (27) we obtain the phase oscillation of the autocorrelation function on the timescale This oscillation takes place when the wave-packet moves over the initial position. At the classical cyclotron time the wave-packet comes back to the initial position and thus the autocorrelation function approximately reaches its initial value, as shown in figure 6(a). For a revival the terms  proportional to the second derivative have to be a multiple of 2π. This leads to a revival for which is best seen in the time-evolution of the absolute value of the autocorrelation function ( figure 6(b)). Of course, the revival hierarchy continues to higher orders. In atomic physics the next time is the so called "super revival time" The huge range of timescales requires a highly accurate wave-packet propagation, since the phase of the propagated wave-packet needs to be accurate for at least 4|n| 2 phase oscillations in the autocorrelation function. The used polynomial propagation is well suited for this task since the accumulated error is reduced by using only very few but long timesteps.

Fractional revivals
Interestingly, we can also distinguish fractional revivals occurring between the full revivals [31] which are again encoded in the autocorrelation function [32]. Fractional revivals have been experimentally measured in Rydberg atoms for the first time [33]. The underlying phase effect can be also employed for prime number factorisation [34,35]. Recently several experiments have implemented this scheme into a NMR setup [36], cold atoms [37] or optically [38] and successfully factorised numbers. Note that the phase properties for wave-packets in graphene differ from those of the aforementioned experiments.
In order to extract the fractional revivals one has to calculate the autocorrelation function in a Poincaré map for classical cyclotron time differences. Figure 7 displays the fractional revivals for several t = (n + p/q) T cl where n is an integer and p/q represents an irreducible fraction. Suitable fractions for different fractional revivals are found by slicing the circle depicted in figure 7(a). For a full revival the autocorrelation function for p/q = 0 needs to be studied, whereas from p/q = 1/2 mirror revivals are identified. The next higher fractional revivals with a triangular shape have the time shifts p/q = 0, p/q = 1/3, p/q = 2/3 and p/q = 1/2, p/q = 1/6, p/q = 5/6. The absolute value of the autocorrelation for these times is shown in figure 7(b) and indicates the times when revivals occur. A full revival takes place if the absolute value for p/q = 0 becomes unity again, which is the case for t = T rev . The revival for t = 1/2 T rev has a minimum for p/q = 0 and a maximum for p/q = 1/2, explaining the discrepancy between the classical centre-of-mass and the quantum mechanical result. Fractional revivals with two peaks are visible at times t = 1/4 T rev and t = 3/4 T rev . In this case the absolute value of the autocorrelation function has to reach almost 1/2. Fractional revivals with a triangular shape can occur in two different ways. For the first one a maximum is located at the opposite side of the initial starting point. In this case the autocorrelation function for p/q = 1/2, p/q = 1/6, p/q = 5/6 is 1/3 as shown in the numerical data. As second possibility for a triangular revival the same result occurs for p/q = 0, p/q = 1/3, p/q = 2/3. Also the next higher revivals can be seen in the numerical data but are not presented here to keep clarity.

Effect of impurities
For graphene based devices, a good understanding of the effect of impurities on the electronic motion is required. Suspended graphene has a high mobility leading to mean free path lengths of about 2 µm [39]. However, the mobility is reduced by the formation of intrinsic ripples which stabilise the two-dimensional system [40]. The lattice deformations involved can be added to the Dirac Hamiltonian by a position dependent gauge field [41,42], which can be handled by the propagation algorithm described in section 2. Additionally the curvature leads to a spin orbit interaction [43]. Here, we study the effect of a scalar potential which is commonly used to model impurities [44].
We use a simple model with a Gaussian distributed random potential and analyse the propagation of a cyclotron wave-packet through this perturbed environment. The potential was generated by the convolution of a random potential U n at the grid points r n of a square lattice: The potential has a vanishing average expectation value U imp (r) = 0, a variance of U imp (r) 2 = U 2 0 , and is correlated by U imp (r)U imp (r ′ ) = U 2 0 exp (−|r − r ′ | 2 /(2ξ 2 )) with the correlation length ξ. A similar potential has already been used in conductance calculations for graphene on a honeycomb lattice [44]. We focus on cases where the correlation length ξ is longer than the grid spacing leading to a potential which is a smooth function with respect to the graphene unit cell. In the following, we study several impurity configurations and calculate the propagation of a wave-packet which has the same initial shape as a Landau level wave-packet of the clean system. As first observable we investigate the autocorrelation function. All the calculated results have in common that the autocorrelation function vanishes much faster compared to the clean system. In figure 8(a), the grey line shows the autocorrelation function of the clean system with prominent revivals whereas the autocorrelation function in the presence of an impurity potential decays after a few picoseconds. This phenomenon happens for many different setups since the wave-packet simply drifts away. Some impurity configurations lead to a wave-packet motion, where the wave-packet returns to its initial position after some time. The return results in an increased autocorrelation function but differs from the revivals occurring in the clean systems, which reflect the initial occupation of the Landau levels. As a measure for the deviation of a disordered system from the clean system we calculate the fidelity We track the overlap of the wave-packet propagated by the unperturbed HamiltonianĤ and the perturbed HamiltonianĤ ′ =Ĥ+U imp . For the free electron gas, expressions are known which describe the fidelity decay in disordered systems for differently correlated impurity potentials [45]. We obtain similar numerical results, shown in figure 8(b). For short times the fidelity is approximately unity and then decays exponentially. After longer times, the fidelity does not approach zero since there is a non-vanishing probability to return back to the initial position. However we cannot infer from the autocorrelation function whether wave-packet revivals occur or not, since the exponential fidelity decay removes all relevant information from the calculated data.
Consequently we have to find other observables which signify revivals even if an impurity potential is present. One particularly suitable observable is the centre-of-mass motion of the wave-packet. In figure 9(a-c) the centre-of-mass motion of a Dirac wavepacket is shown for different impurity potentials. In figure 9(a) the wave-packet does not drift away and thus the revivals lie on top of each other. This leads to a revival in the autocorrelation function which has nothing in common with the wave-packet revivals of the clean system since the revival time depends crucially on the exact choice of the impurity potential. However the majority of the disorder configurations lead to a drift of the wave-packet as shown in figure 9(b,c). In all cases the wave-packet moves at the beginning on a shrinking spiral since the centre-of-mass is located at the position with the highest probability density (figure 10(a)). In this state the centre of the spiral already moves. After a few picoseconds the centre-of-mass collapses to the centre of a rotating circular wave-packet as shown in figure 10(b,c). In the following the wavepacket is moving without oscillations on a smooth line because the density is equally distributed along a circular orbit. As in the clean system, the first revival occurs in the form of a mirror revival at time 1/2 T rev . Then the centre-of-mass collapses again and produces a full revival at T rev . A way to remove the structure of the impurity potential is extracting the corresponding velocity of the wave-packet out of the centre-of-mass motion which is shown in figure 9(d). The initial wave-packet moves on the circular orbit with a velocity slightly smaller than the "speed of light" of graphene. When the centre-ofmass collapses this velocity drops down to the drift velocity of the wave-packet which is defined by the potential gradient and at least two orders of magnitude smaller. The revivals occur at the same time for all the impurity configurations. The substructure of the revivals cannot be seen in the centre-of-mass motion of the wave-packet. However the full wavefunction shows fractional revivals ( figure 10).
In reference [46] the centre-of-mass propagation was used to establish a dipole moment. In our case this dipole moment is maximal if the wave-packet revival takes place and almost zero if the wave-packet is collapsed similar to the velocity of the wave-packet in figure 9(d). If such a dipole moment is experimentally measurable, the wave-packet revivals can be observed even if disorder displaces the wave-packet far away from its initial position.

Conclusion
In this work we presented a new algorithm for the propagation of wave-packets on a sheet of graphene which allows us to calculate the time-evolution in arbitrary shaped potentials and inhomogeneous magnetic fields. We studied time-dependent effects like the zitterbewegung, which is very pronounced in graphene. We focussed on wave-packet revivals, which occur in graphene due to the special structure of the Landau levels. We gave a detailed explanation of the centre-of-mass motion and the autocorrelation function. The applied semiclassical description quantitatively revealed the collapse of the centre of mass as well as the revivals. From the autocorrelation function we extracted several different timescales which all signify various revivals. A detailed analysis of the autocorrelation function for different Poincaré times confirmed fractional revivals which have also been seen in our numerical calculations which supported the results throughout the whole work. Finally we analysed the effect of impurities on the wave-packet revivals by a Gaussian correlated model potential. Our calculations showed that revivals in the autocorrelation are strongly suppressed due to the fidelity decay. In contrast to that the centre-of-mass revivals are still present if the correlation length of the impurity potential is long enough with respect to the average cyclotron radius. The here demonstrated longtime accuracy of the propagation algorithm opens the window towards realistic device simulations on the micrometre range, also including time-dependent external fields. and hole-like states (λ = −1) with the corresponding energy dispersion E λ = λ vh|k|. Then it is possible to construct a wave-packet which is restricted to one band. The time-evolution of the momentum and the centre of this wave-packet are defined by two coupled equations of motioṅ where B = Be z stands for an isotropic magnetic field perpendicular to the plane and E(r) is a position dependent electric potential. The Berry curvature Ω λ (k) enters the equations of motion in momentum space but is zero in the case of a gap-less perfect sheet of graphene. As a result of the linear dispersion, relation (A.2) reduces toṙ = λvk/|k| and therefore describes a centre-of-mass which always travels with the speed v parallel or antiparallel to the wave vector k. For a vanishing electric field the two coupled equations of motion can be solved. The momentum of the wave-packet is controlled by the differential equation where the initial momentum of the wave-packet is given by k 0 in direction ϕ. The cyclotron frequency depends on the momentum and is given by The resulting real space propagation (A.2) is given by r = r 0 + λ l c sin(ω c t + ϕ) − cos(ω c t + ϕ) . (A.7) Thus electron-like and hole-like states propagate in opposite directions on orbits with the radius In order to obtain the quantisation of the Landau levels from the closed cyclotron orbits we have to study the phase changes along the trajectories. The main contribution is given by the classical action known from the electron case, while an additional contribution comes from Berry's phase. Although the Berry curvature is vanishing, Berry's phase [48] is still present and defined by the vector potential where ϕ is the direction of the vector k. Hence we get a vector potential of which is used to calculate the additional phase-change Γ λ (ω) = ω A λ (k)dk (A. 12) for an arbitrary adiabatic transition on the path ω(t) in momentum space. The result is Berry's phase Γ ± (ω) = ±π (A.13) for a single closed orbit on graphene [2]. The adiabatic transition in momentum space is automatically fulfilled in a magnetic field, since all trajectories are cyclotron orbits and thus sufficiently smooth. The classical phase change and Berry's phase have to be added up to build a new EBK quantisation condition which reads [49] where ν is the Maslov index depending on the caustics traversed on one cyclotron orbit (here ν = 2). Plugging the formula for the cyclotron motion in the quantisation condition the resulting Landau levels are given by