Light-Induced Valleytronics in Pristine Graphene

It is often assumed that achieving valley selective excitation in pristine graphene with all-optical means is not possible due to the inversion symmetry of the system. Here we demonstrate that both valley-selective excitation and valley-selective high-harmonic generation can be achieved in pristine graphene by using the combination of two counter-rotating circularly polarized fields, the fundamental and its second harmonic. Controlling the relative phase between the two colours allows us to select the valleys where the electron-hole pairs and higher-order harmonics are generated. Present work offers a robust recipe to realise valley-selective electron excitations in materials with zero bandgaps and zero Berry curvature.

The realisation of atomically-thin monolayer graphene has led to breakthroughs in fundamental and applied sciences [1,2]. Charge carriers in graphene are described by the massless Dirac equation and exhibit exceptional transport properties [3], making graphene very attractive for novel electronics applications. One of the most interesting features of graphene and gapped graphene materials is the electron's extra degree of freedom, the valley pseudospin, associated with populating the local minima K and K in the lowest conduction band of the Brillouin zone. This extra degree of freedom has the potential to encode, process and store quantum information, opening the field of valleytronics [4].
The monolayer graphene, as opposed to gapped graphene materials, presents a fundamental challenge for valleytronics: it has zero bandgap and zero Berry curvature. These aspects are generally considered to be a major impediment for valleytronics. In gapped graphene materials, valley selectivity is achieved by matching the helicity of a circularly polarized pump pulse, resonant to the bandgap, to the sign of the Berry curvature [5][6][7][8][9].
Recently demonstrated [10] sub-cycle manipulation of electron population in K and K valleys of tungsten diselenide, achieved with the combination of a resonant pump pulse locked to the oscillations of the THz control pulse, represents a major milestone. Precise sub-cycle control over the driving light fields opens new opportunities for valleytronics, such as those offered by the new concept of a topological resonance, discovered and analysed in Refs. [11][12][13]. Single-cycle pulses with the controlled phase of carrier oscillations under the envelope offer a route to valleytronics in gapped graphene-type materials even when such pulses are linearly, not circularly, polarized [14]. It is also possible to avoid the reliance on resonances in gapped graphene-type materials, breaking the symmetry between the K and K valleys via a light-induced topological phase transition, closing the gap in the desired valley [15].
Thus, with its zero bandgap, zero Berry curvature, and identical dispersion near the bottom of the valleys, pristine graphene appears unsuitable for valleytronics -a disappointing conclusion in view of its exceptional transport properties. We show that this generally accepted conclusion is not correct, and that the preferential population of a desired valley can be achieved by tailoring the polarization state of the driving light pulse to the symmetry of the lattice. Our proposal offers an all-optical route to valleytronics in pristine graphene, complementing approaches based on creating a gap by using a heterostructure of graphene with hexagonal boron nitride [16][17][18][19], or by adding strain and/or defect engineering [20][21][22][23].
The key idea of our approach is illustrated in Fig. 1, which shows graphene in real (a) and reciprocal (b) space, together with the structure of the incident electric field (a) and the corresponding vector potential (b). The field is made by superimposing two counterrotating circularly polarized colours at the frequencies ω and 2ω. The Lissajous figure for the total electric field is a trefoil, and its orientation is controlled by the relative two-colour phase φ, i.e., the two-colour delay. In the absence of the field, the two carbon atoms A and B, in real space, are related by the inversion symmetry. When the field is turned on, this inversion symmetry is broken: the electric field in panel (a) always points from the atom A to the one of the three neighbouring atoms B during the full laser cycle, but not the other way around. If the centre of the Lissajous figure is placed on the atom B, the field points in the middle between its neighbours. One can control this symmetry breaking by rotating the trefoil, interchanging the roles of the atoms A and B. Thus, the bi-circular field offers a simple, all-optical, ultrafast tool to break the inversion symmetry of the graphene lattice in a controlled way. Such controlled symmetry breaking naturally controls the relative excitation probabilities induced by the same laser field in the adjacent K and K valleys of the Brillouin zone.
Panel (b) provides the complementary reciprocal-space perspective. In the laser field, the electron crystal momentum follows k(t) = k i + A(t), where k i is the initial crystal momentum and A(t) is the laser vector potential, shown in (b) for the electric field shown in (a). The asymmetry between the two valleys with respect to the vector potential is immediately visible. Panels (c) and (d) provide additional support to this qualitative picture. One can observe that the two field-free valleys, K and K , are only identical near their very bottoms. As soon as one moves away from the bottom, the valleys start to develop the trefoil structure, with K and K being the mirror images of each other. How the symmetry of the vector potential fits into the symmetry of the valley away from their bottoms will control the dynamics and the excitation probability.
In sufficiently strong fields, excitation happens not only at the Dirac point where the gap is equal to zero, but also in its vicinity. For the vector potential in panel (c), the average gap seen by the electron when following the vector potential from the Dirac point in the K valley, i.e., moving along the trajectory K + A(t), is less than in the K valley, i.e., when following the trajectory K + A(t). In sufficiently strong and low-frequency fields, such that the bandgap along the trajectories K + A(t) and K + A(t) quickly exceeds the photon energy, the excitation probability should therefore be higher in panel (c). For the same reason, rotating the vector potential as shown in panel (d) should favour population of the K valley.
These expectations are confirmed by our numerical simulations. In the simulations, graphene is exposed to the bicircular field with the vector potential Here, A 0 = F ω /ω is the amplitude of the vector potential for the fundamental field, F ω is its strength, f (t) is the temporal envelope of the driving field, φ is the sub-cycle phase difference between the two fields, and R is the ratio of the electric field strengths for the two fields, leading to R/2 ratio for the amplitudes of the vector potentials. The amplitude of the fundamental field was varied up to F ω = 15 MV/cm, leading to the maximum fundamental intensity 3×10 11 W/cm 2 , with the fundamental wavelength λ = 6µm. This laser intensity is lower than the one used in HHG in graphene experimentally [47] and is below the damage threshold. We have also varied R, using R = 1 and R = 2. A pulse duration of 145 femtosecond with sin-square envelope is employed in this work. Our findings are valid for a broad range of wavelengths and field intensities. To obtain the total population of the different valleys, we have integrated the momentum-resolved population over the sections shown in Fig. 2(a).
To quantify the amount of the valley polarisation, we used the valley asymmetry parameter defined as where n K c and n K c are electron populations at the end of the laser pulse in the conduction band around K and K valleys, respectively.  Fig. 2(c, d). The asymmetry in the valley population is negligible up to an intensity of 2 × 10 11 W/cm 2 and gradually increases with intensity [ Fig. 2 (b)]. This shows that the valley asymmetry is observed only when the laser pulse is able to drive the electrons to the anisotropic part of the conduction band.
The current generated by the electron injection into the conduction band valleys is accompanied by harmonic radiation and makes substantial contribution to the lower order harmonics, such as H4 and H5 for our two-colour driver. These harmonics are stronger whenever the dispersion (k) is more nonlinear. In this respect, for the electrons following the trajectories K + A(t) and K + A(t), the current-driven high harmonic generation from 6 the K valley is preferred for the vector potential in Fig. 1(c). Conversely, for the vector potential in Fig. 1(d) low-order, current-driven harmonic generation should be preferred from the K valley. Indeed, following the vector potential, the electron is driven against the steeper walls in the K valley in panel (c) and against the steeper walls in the K valley in panel (d). These qualitative expectations are also confirmed by our numerical simulations, shown in Fig. 3. The same laser parameters used in Fig. 2 are used in this simulation.
In general, the interband and intraband harmonic emission mechanisms in graphene are coupled [48][49][50], except at low electric fields [50], leading to a complex interplay between the interband and intraband emission mechanisms. The former should be more sensitive to the dephasing than the latter. To this end, we have calculated the the harmonic spectrum for 7 different dephasing times T 2 , see Fig. 3(a). We find that the harmonic emission is essentially T 2 -independent, until at least H13, suggesting that the intraband emission mechanism is generally dominant. Fig. 3(b) shows polarisation-resolved high-harmonic spectrum. The (3n + 1) harmonics follow the polarisation of the ω pulse (left-handed circular polarisation), whereas (3n + 2) harmonics follow the polarisation of the 2ω pulse (right-handed circular polarisation), while 3n harmonics are missing, just like in atomic media [51]. In this context, we note that while harmonic generation with single-colour circularly polarized drivers is forbidden in atoms, such selection rules do not generally arise in solids [31]. However, the ellipticity-dependence studies on graphene show very weak harmonic yield for drivers with higher ellipticity [47][48][49], as we have also observed. As in atoms, application of bicircular fields allows for efficient generation of circularly polarized harmonics in graphene.  Fig. 1(c), while the maximum contribution of the K valley corresponds to the vector potential orientations such as shown in Fig. 1(d). Therefore, we are able to control the valley-polarisation of the harmonics by controlling the two-colour phase φ. We have also checked that these results do not depend on the specific directions of rotation of the two driving fields. That is, we find the same results when simultaneously changing the helicities of both driving fields.

METHODS
In the simulations, we used the nearest-neighbour tight-binding approximation to obtain the ground state of graphene with a hopping-energy of 2.7 eV [52]. The lattice parameter of graphene is chosen to be 2.46Å. The resultant band-structure has zero band-gap with linear dispersion near the two points in the Brillouin zone known as K and K points.
The density matrix approach was used to follow the electron dynamics in graphene.
Time-evolution of density matrix element, ρ k mn , was performed using semiconductor Bloch equations within the Houston basis |n, k + A(t) as [53,54] Here, k mn and d k mn are, respectively, the band-gap energy and the dipole-matrix elements between m and n energy-bands at k. Dipole matrix elements were calculated as d k mn = -i u k m |∇ k |u k n , where u k n is the Bloch function corresponding to the total wavefunction. F(t) and A(t) are, respectively, the electric field and vector potential of the laser field and are related as F(t) = -∂A(t)/∂t. A phenomenological term accounting for the decoherence is added with a constant dephasing time T 2 . Conduction band population relaxation was neglected [55].
The total current was calculated as Here, p k nm are the momentum matrix elements, obtained as p k nm = n, k ∇ kĤk m, k . The off-diagonal elements of momentum and dipole-matrix elements are related as d k mn = ip k mn / k mn . Finally, the harmonic spectrum was determined from the Fourier-transform of the timederivative of the total current as Here, integration is performed over the entire Brillouin zone.

DATA AVAILABILITY
The datasets that support the findings of this study are available from the corresponding authors on reasonable request.