Local excitation and valley polarization in graphene with multi-harmonic pulses

We elucidate the mechanism of strong laser pulse excitation in pristine graphene with multi-harmonic pulses, linearly polarized parallel to the line connecting the two di ﬀ erent Dirac points in the Brillouin zone and with a maximal vector potential given by the distance of those points. The latter two conditions have emerged from our previous work [Kelardeh et al. , Phys. Rev. Res. , 2022, 4 , L022014] as favorable for large valley polarization. We introduce a novel compacted representation for excitation, locally resolved in the initial conditions for the crystal momenta. These maps are our main tool to gain insight into the excitation dynamics. They also help with understanding the e ﬀ ect of dephasing. We work out why a long wavelength and a moderate number of overtones in the harmonic pulse generate the largest valley polarizations.


Introduction
Illuminating condensed-matter targets with intense short pulses affects nonlinear photo-absorption which depends crucially on the band structure of the system and is commonly probed by high-harmonic generation. 1,2 Pioneering work has shown that the basics of the ensuing electron dynamics follows a similar twostep process in reciprocal space, 3,4 as is familiar from atoms and molecules in real space. Firstly, a valence electron tunnels at the gap to the energetically closest conduction band appearing there with wave vector k. Subsequently, the wave vectork moves on the conduction band driven by the lasers's vector potential A according tokðtÞ ¼ k þ AðtÞ. Owing to the gap, electronic transitions back to the valence band are exponentially suppressed. This is not the case for gapless 2D systems such as graphene, rendering the ensuing (coherent) excited electron dynamics potentially more complicated. Moreover, dephasing may alter the nal excitation for laser pulses longer than the coherence time. Physical properties of interest are the excitation probability and valley polarization. [5][6][7][8] The latter quanties the difference of excitation originating from initial wave-vectors in the rst Brillouin zone (BZ) belonging to the domains of the two Dirac points K and K 0 of the vanishing gap. These two domains of Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany. E-mail: us@pks.mpg.de; rost@pks.mpg.de triangular shape form a unit cell of the BZ, well adapted to describe valley polarization, see Fig. 1. Valley polarization (VP), detectable via even harmonics, 6 has generated considerable interest, since it promises a means of storing and manipulating quantum information by phase changes only, although this is still quite a distant goal. 8 In recent work, 9 we have shown that laser elds of sufficient strength to affect excitation, can be chosen such that transitions between bands occur in graphene essentially only close to the Dirac points. Hence, band transitions happen at instants of time when the intense pulse driven electron wave vectorq(t) passes by K or K 0 leading to a sequence of (coherent) electron amplitudes in the conduction band.
As an immediate consequence and in contrast to common belief, linearly polarized pulses can generate substantial valley polarization for a gapless system 9 under two conditions: (i) the linear polarization is parallel to the difference vector D ¼ K À K 0 of the two Dirac points and the vector potential A 0 is of the order of the distance D ¼ 4p=ð ffiffiffiffiffi 27 p aÞz0:90 a.u. between the two Dirac points, where a ¼ 1.43 A is the nearest-neighbor C-C spacing in graphene. (ii) The pulse couples differently to electrons at K and K 0 . This requires an asymmetry † in the driving pulse which is, from a theoretical point of view, most simply achieved with a half cycle pulse. 9 A sufficiently short pulse also ensures that dephasing does not play a role. Experimentally, such pulses can be generated thanks to recent experimental advances, 10,11 requiring, however, special equipment. With one notable exception, 12 to the best of our knowledge pristine (gapless) graphene exposed to light polarized in the graphene plane has only been investigated with linear polarization perpendicular 13,14 to D or with ultrashort circular pulses, where the vector potential's component parallel to D was symmetric, 14,15 with jmax{A D (t)}j ¼ jmin{A D (t)}j. These kind of pulses cannot generate VP.
In the following, we investigate the interplay of coherent dynamics and dephasing in graphene as a gapless system with a multi-harmonic pulse whose characteristics in terms of pulse shape, pulse length and wavelength can be Fig. 1 Sketch for the compact representation of the triangular q-domains around K and K 0 , which we will use below. We fold the Voronoi triangles into a rectangle (a / b). Then the two outer rectangles, confined by the dotted lines, are overlaid with the inner ones (b / c) by a p-rotation around the center between K and K 0 , marked with crosses. Since the dynamics are symmetric with respect to the dashed line, we use the left part for presenting the local excitation c + (q) and the right part for the local polarization c À (q), respectively. The probabilities c AE are defined below in eqn (5), their corresponding color scales are given in the right panel of d. The momenta are given in terms of D ¼ jK À K 0 j, the distance between the K-points. systematically changed, while maintaining a basic asymmetry of the pulses to enable valley polarization. The main quantities which will provide detailed insight are local probabilities, resolved w. r. t. the initial conditions in the Brillouin zone. We will analyze excitation and valley polarization. To this end we calculate the time-dependent density matrix in a two-band approximation with a phenomenological dephasing rate g and we introduce a compact graphical representation for the local probabilities. If the laser pulses are intense as in our case, the two-band approximation is conrmed to provide an excellent description. 14 The model is briey described in Section 2, where we also explain the multi-harmonic pulse and illustrate with an example, the basic light-driven electronic excitation dynamics. Section 3 presents the results for different numbers of overtones in the pulses and uncovers how excitation and valley polarization depend on the parameters of the pulse, elucidating how intense pulses steer electrons in graphene. Section 4 uses excitation distributions in connection with the compacted local probabilities to work out the differences between coherent and dephased dynamics and to identify those initial conditions in the unit cell which mainly contribute to excitation and/or VP. The paper ends with a short summary in Section 5. We use atomic units, unless stated otherwise.

Graphene dynamics including dephasing
The time-dependent Hamilton matrix for the well known two-band approximation for graphene reads in reciprocal space with the Pauli matrices s x and s z . The Hamilton matrix eqn (1), follows from an instantaneous diagonalization of the standard hopping matrix, 16 which denes also the so-called Houston basis. 17 In this basis we obtain the (instantaneous) band gap E q+A(t) and the (instantaneous) eigenvectors V q+A(t) and W q+A(t) for the valence and conduction band, respectively. Find a derivation for H in appendix A.
Note that both, the vector potential A(t) and the electric eld affect the time-dependent Hamilton matrix H(q,t) in eqn (1) and inuence the time evolution of the density matrix through v vt r q ðtÞ ¼ Ài Â ℍðq; tÞr q ðtÞ À r q ðtÞℍðq; tÞ where g is the dephasing rate. Products with r, H, and s x,z in eqn (2) are matrixmatrix multiplications.

Compact visualization of differential excitation and valley polarization
Excitation in the triangular q-domains VK and VK 0 around K and K 0 of the initial conditions, representing together a unit cell of the Brillouin zone (see Fig. 1a), are dened as averaged local excitations c(q) in the respective domains of area V in qspace, with conduction-band occupation c as a function of initial conditions q given by the density matrix element c(q) ¼ r q,cc (t / N). In terms of P K and P K 0 , the total excitation p averaged over the entire unit cell 2V, and the total valley polarization h are dened by We obtain c(q) by numerically solving eqn (2) for a dense set of initial conditions q, sampling each of the two triangular domains with 512 Â 222 ¼ 113 664 initial conditions. Therefore we can dene excitation p and valley polarization h locally in the initial conditions q with a resolution inversely proportional to the sampling. A natural graphical representation consists of color maps (scales are given in Fig. 1d) lling the area of the two triangles in Fig. 1a. Due to the underlying hexagonal symmetry, we can fold the triangles into a rectangle, cf. the transition from Fig. 1a and b. 9 Since in the end, p and h emerge from sums and differences over local conduction-band occupations c(q), one can perform part of this operation graphically by overlaying equivalent areas in the K and K 0 triangles to arrive at Fig. 1c, which produces an even more compact representation. ‡ This second step consists essentially of a p-rotation around the center of the line segment K-K 0 , which superimposes each point in VK with one point in VK 0 : q / q 0 ¼ f(q). By means of this mapping f we can dene the local conduction-band occupations for p and h, respectively, They enable a detailed analysis, either in contour plots (Fig. 5,7,9) or distribution functions (Fig. 8). Of course, p and h, dened traditionally in eqn (4), follow equivalently from their local contributions eqn (5) For convenience, below we will use the terms "local excitation" and "local polarization" for c + and c À , respectively.

Asymmetric pulse shapes from multi-harmonic pulses
Large ionization probability all the way to damaging the sample is easy to achieve by increasing the eld strength, in close analogy to atoms or molecules, where damaging corresponds to depletion of the ionized electron orbital. This, however, implies that excitation also occurs further away from the Dirac points at nite band gaps, and eventually, for large enough eld strength throughout the unit cell with the consequence of a vanishing VP.
The generation of signicant VP hinges upon preserving the vicinity of the Dirac points as they are dominantly responsible for the excitation. This provides the means for selective excitation of the two valley domains upon breaking the equivalence in the excitation via the two Dirac points through suitable pulse forms. These requirements select combinations of appropriate eld strength, vector potential and frequencies of the light pulse. An obvious choice fullling these requirements is a pulse close to a half-cycle, of suitable strength and linearly polarized along D, as demonstrated before. 9 To make the interplay of laser parameters for excitation and VP more apparent and also, to assess the effect of dephasing, we will use multi-harmonic pulses. Asymmetric by construction, they look like pulse trains and provide a means to study the interplay of coherent multi-excitation electron dynamics and decoherence. We dene these pulse trains as superpositions of harmonic overtones with an envelope determining the overall pulse length T, The relevant parameters in the pulse eqn (6) are the fundamental frequency u (dened by the wavelength l), the amplitude A 0 , the (full-width-at-halfmaximum) pulse duration T, the number of harmonics n and the harmonic phase f. For n > 1 and f ¼ 0 the periodic part in eqn (6a), i.e. the term given by the sum cosðkwtÞ=k, is highly asymmetric. The contrast between maximum and minimum within one period T u ¼ 2p/u is given by C n ¼ ÀS n (t max )/S n (t min ) with t max ¼ T u . For odd n, t min ¼ T u /2, for even n there is no explicit analytical expression of t min . Hence, for odd n whereas we have to calculate C n for even n numerically. Typical values are C 1 ¼ 1, C 2 ¼ 2, C 8 z 3.69, and C 32 z 5.74. Apparently, the contrast increases with n, see Fig. 2.

A paradigmatic example of excitation by a multi-harmonic pulse
Before we discuss excitation and VP resulting from different multi-harmonic pulses, we briey illustrate with Fig. 3 how the strong eld excitation mechanism proceeds through the Dirac points with the specic harmonic pulse shown in Fig. 3c. The time-dependent crystal momentumq(t) ¼ q + A(t), with the ycomponent e y $q(t) shown in Fig. 3b, is mapped into the unit cell in Fig. 3a with the green line and the initial condition q as a white-lled circle. Maximal and minimal excursion occur close to the maximum of the envelope as marked with arrows. During its evolution,q(t) oscillates back and forth on the green line with the amplitude shown in Fig. 3b. Once this amplitude is large enough,q(t) passes repeatedly by the Dirac point K 0 (cross in the lower triangle and dashed line in Fig. 3b) as indicated by thin black lines in Fig. 3b-d. At these times, a sequence of transitions occurs between the valence and the conduction band as stated before and is clearly visible in Fig. 3d. Whether they lead to net excitation or de-excitation at a specic encounter of K 0 depends on the relative phases of the electron excitation amplitudes.
The nal excitation c(q) for each initial condition is quantitatively indicated by the heat map in Fig. 3a revealing a quite complicated pattern across the unit cell but clearly conned laterally (in the q x -direction) around the straight line passing through the Dirac points. Overall, the asymmetric pulse affects theq(t) with initial condition in the K 0 triangle not equivalently to those starting in the K triangle regarding excitation. This holds true despite the fact that the Dirac points, responsible for the transitions, are visited by equal areas of initial conditions from both triangles of the unit cell, irrespective of the pulse shape. Differences in the excitation-the prerequisite for VP-emerge from the intricate coherent dynamics ofq(t), sensitively depending on the respective initial conditions (see  Fig. 3d). How these VP generate differences and excitation in the rst place, depends on global features of the pulse and will be elucidated below.

Graphene in two-color pulses
For a rst impression, we investigate how the averaged excitation p and VP h depend on wavelength and dephasing, we present in Fig. 4 results for harmonic pulses eqn (6) with n ¼ 2, containing the fundamental frequency corresponding to the wavelength, and the rst overtone. One can see that p decreases while h increases almost linearly and changes sign. Moreover, dephasing reduces p by roughly 1/3, while it hardly affects h. Interestingly, the pulse with the highest h at l ¼ 8 mm is quite long with T ¼ 160 fs. One can understand these trends better from the local excitation and polarization c AE (q), shown in Fig. 5. Not surprisingly,  Shown is c(q), for l ¼ 2 mm only, and in compact form, cf. sketch in Fig. 1, the excitation c + (q) and the polarization c À (q). We show both results with (lower row) and without (upper row) damping. Color scales are given in Fig. 1d.
dephasing blurs the interference pattern leaving less unexcited area (white) in the unit cell, cf. upper and lower row in Fig. 5. This is particularly obvious for the green color maps showing c + in compact representation, as all the upper panels have more white space than the lower ones.
One notes also, that c + is more spread over the unit cell and is overall larger for smaller wavelengths (more green areas). For increasing wavelength, the excitation shrinks laterally towards the K 0 point (black circle) and offers the selectivity which is necessary for large c À , and indeed, valley polarization increases towards larger wavelengths with the warm colors dominating the cold colors. Since we keep the vector potential xed at the distance between the two Dirac points, A 0 ¼ D, decreasing the wavelength implies lowering the eld strength. Therefore, transitions between the bands happen closer to the Dirac points. Next, we investigate what happens if n is increased systematically beyond 2 in multi-harmonic pulses.

Multi-harmonic pulses: dependence on number of harmonics
The arguments about the wavelength dependence of p and h given above for n ¼ 2 do not explicitly depend on the number of harmonics in the pulse. Therefore, we expect the same qualitative behavior when increasing n. Fig. 6 reveals that this is only the case for the effect of dephasing, which remains generally larger for p than for h and hardly changes for p when n is increased. Surprisingly at rst glance, neither p nor h depend monotonically on n. Rather p and h exhibit in most cases a maximum for a specic n max (l)˛ [2,32], which depends on the wavelength.
Indeed, the VP approaches zero for small, and towards large, n. The limit n ¼ 1 corresponds to an ordinary pulse which is symmetric in the vector potential and therefore cannot generate VP. The reason for vanishing VP towards large n is more subtle: large n implies that with large time derivatives of the vector potential's harmonic components cos(kut) there are large electric elds, cf. upper row in Fig. 2. These pulse parts impair the selectivity of the excitation and therefore the ability of the pulse to generate large h in the rst place, despite the fact that increasing n produces increasingly asymmetric pulses in the vector potential, cf. lower row in Fig. 2. The even more complicated pattern, with two extrema of smaller absolute values of h for l ¼ 2 mm in Fig. 6, can be also attributed to lower selectivity, since for smaller wavelength, excitation is in general less conned in the unit cell of initial conditions.
The spread of the excitation over the unit cell with increasing n and for small l can be directly read off the compact color maps for c + in Fig. 7, where one sees that the panels for l ¼ 2 mm contain more green than those for l ¼ 8 mm and that n ¼ 32 has in general the most dark green areas. Of course, this dependence on n is only pronounced if the wavelength is long enough to not cover the unit cell already for n ¼ 2 in the excitation, which is the case for l ¼ 2 mm: the corresponding p is only weakly dependent on n (see Fig. 6).

Multi-harmonic pulses: dependence on the harmonic phase
For completeness we also show in Fig. 6 different harmonic phases f whose effect is obvious: they hardly inuence p, but with increasing f the vector potential within a single period becomes more symmetric (see lower row in Fig. 2) decreasing VP and for f ¼ p/2 with fully symmetric A(t), it is impossible to generate VP. Interestingly, dephasing has less inuence for more symmetric vector potentials (the difference between p with and without dephasing is much less for f ¼ p/2 than for f ¼ 0). Regarding VP, there is in general very little inuence of dephasing regardless of the harmonic phase f. Note that the pulse length T varies with the wavelengths between 40 and 160 fs, i.e., ultrashort pulses are not required for VP.   Most striking is the qualitative difference of coherent dynamics with smooth distributions (grey areas) and dephased dynamics exhibiting pronounced peaks (colored lines); note, however, the logarithmic scale. Not surprisingly, c À ¼ 0% dominates as VP is difficult to achieve. How is this central peak and the (smaller) peaks at c À ¼ AE50% related to the dominant peaks in the distributions of c + ? This can be answered by looking at the differential probabilities in Fig. 9. Trivially, all initial conditions which do not lead to excitation (not Fig. 8 Distribution of local excitation c + and polarization c À for three wavelengths (l ¼ 2, 4, 8 mm) and a two-color pulse (n ¼ 2) with maximal contrast (f ¼ 0). Results for a finite damping time g ¼ 1/10 fs (colored lines) are compared to those without damping (gray shaded). shown) contribute to the peak at c À ¼ 0% in Fig. 8, as can be seen from Fig. 9g and also from the fact that all ED have their maximum at c + ¼ 0. Note, however, that also the highest abundance of excitation (peak at c + ¼ 50% in Fig. 8d-f) does not contribute to c À , which becomes clear when looking at Fig. 9a, b and g.

The effects of dephasing and excitation distributions
The peak at c + ¼ 25%, on the other hand, is the reason for the peak at c À ¼ 50%. This becomes apparent, when comparing Fig. 9c with Fig. 9h. Although this peak constitutes the main contribution to VP, it does not come from the initial conditions with the highest excitation, which correspond to the green area in Fig. 9a. Rather, the strongest VP comes from excitation at the tip of the triangle in Fig. 9h and the red area in Fig. 9b. The principles that the systematics for c À > 0 follows are conrmed for c À < 0 in Fig. 9, where the corresponding differential probabilities are shown for consistency. Since negative c À is only dominant for l ¼ 2 mm, we do not discuss it here further.

Conclusions
We have given a detailed picture of how intense laser pulses affect excitation and valley polarization in graphene, as a representative gapless 2D material. Guided by insight from our previous work, 9 we have focused on pulse trains with asymmetric shapes in each optical period generated by multi-harmonic components, linearly polarized along D ¼ K À K 0 , the line connecting the Dirac points. We have constructed a compact representation of conduction-band occupations in terms of initial conditions in the unit cell, which has given us means to understand the electron dynamics from the local excitation and valley polarization.
The guiding principles to achieve large valley polarization are two-fold: rstly, a pulse, linearly polarized along D, and with a maximal vector potential of A 0 ¼ D. Secondly, (for xed A 0 ) a moderate eld strength to avoid uncontrolled excitation throughout the unit cell. The latter can be achieved by long wavelengths, and a nite number of overtones n in the harmonic pulse, where n is limited by the requirement of moderate eld strength, which-the derivative of the vector potential-becomes large for large n. Although interferences from different excitation amplitudes are omnipresent, they do not inuence excitation and valley polarization strongly for the scenarios investigated. This has become clear from the qualitative agreement of results with and without dephasing. It also means, that here dephasing does not pose an obstacle to achieve large valley polarization in gapless 2D materials.

Appendix: Hamilton matrix and density matrix equation
For completeness we derive here briey the Hamilton matrix eqn (1) in the density matrix eqn (2). The starting point is the usual tight-binding two-band Hamilton matrix 16 where g 0 ¼ 3.03 eV and a ¼ 1.43Å. Quantities in the tight-binding basis get an overbar. In a laser pulse the crystal momentum becomes q /q(t) ¼ q + A(t). We can write both density matrix and Hamilton matrix by means of the transformation TðqÞ, that makes eqn (8) The unitary matrix TðqÞ, with , is oen referred to as the Houston basis. 17 The equation of motion for the Houston basis density matrix r q follows from the tight-binding version by means of eqn (10) v vt r ¼ Ài with the Hamilton matrix TðqÞ T † ðqÞ: The time the derivative of eqn (11) can be rewritten . In order to arrive at eqn (1) we note that Eq ¼ E c (q) À E v(q ) and that the Houston basis TðqÞ consists of two eigenvectors Vq and Wq, for valence and conduction band, respectively. Whereas these vectors are complex, all matrix elements of the 2nd term in eqn (11) are real. The diagonal elements are identical and can thus be omitted. There we have the expression given in eqn (1). We introduce a dephasing term with a dephasing rate g in order to allow for long driving pulses. It will dampen exclusively coherences in the Houston basis density matrix À g 2 Â r q ðtÞ À s z r q ðtÞs z Ã ¼ Àg 0 r vc r cv 0 ! :

Conflicts of interest
There are no conicts to declare.