Theory of pump-probe photoemission in graphene: Ultrafast tuning of Floquet bands and local pseudospin textures

The control of physical properties of solids with short laser pulses is an intriguing prospect of ultrafast materials science. Continuous-wave high-frequency laser driving with circular polarization was predicted to induce a light-matter coupled new state possessing a quasi-static band structure with an energy gap and a quantum Hall effect, coined"Floquet topological insulator". Whereas the envisioned Floquet topological insulator requires well separated Floquet bands and therefore high-frequency pumping, a natural follow-up question regards the creation of Floquet-like states in graphene with realistic pump laser pulses. Here we predict that with short low-frequency laser pulses attainable in pump-probe experiments, states with local spectral gaps at the Dirac points and novel pseudospin textures can be achieved in graphene using circular light polarization. We demonstrate that time- and angle-resolved photoemission spectroscopy can track these states by measuring sizeable energy gaps and quasi-Floquet energy bands that form on femtosecond time scales. By analyzing Floquet energy level crossings and snapshots of pseudospin textures near the Dirac points, we identify transitions to new states with optically induced nontrivial changes of sublattice mixing that can lead to Berry curvature corrections of electrical transport and magnetization.

The ultrafast optical manipulation of materials by femtosecond laser pulses is rapidly becoming a major guiding theme in condensed matter physics [9][10][11] .At the same time, the quest for novel topological states of matter has triggered enormous research activity since the discovery of topological insulators 12 .Merging both of these vibrant fields, a recent work reports the coupling of short laser pulses to surface Dirac fermions in the topological insulator Bi 2 Se 3 13 .On the theory side, the investigation of quasistatic Floquet eigenvalues suggests that topologically trivial bands can be tuned by nonlinearly polarized, continuous driving fields to become topologically nontrivial [5][6][7][8][14][15][16][17][18][19][20] , coined "Floquet topological insulator". The iplied topologically protected edge states were indeed observed in photonic crystals with helical waveguides 21 .These developments pose the natural questions how transient nontrivial states of matter form in real time from equilibrium trivial bands, whether this transition can be tracked on femtosecond time scales using pump-probe spectroscopy, and which nonequilibrium states of matter can be created which have no equilibrium counterparts.In this work, we answer these questions by simulating the real-time development of single-particle energy gaps in graphene coupled to short laser pulses, using realistic parameters for timeresolved, angle-resolved photoemission spectroscopy (tr-ARPES).The subsequent analysis of snapshots of pseudospin textures around the Dirac points allows us to identify characteristics of dynamically generated local Berry curvatures and nontrivial quantum geometries, both in accord with Haldane's original model and beyond.
Guided by Haldane's model of a time-reversal-symmetry-breaking induced Chern insulator on the honeycomb lattice 4 , we start from the minimal honeycomb-lattice tight-binding model of graphene.We drive this system by coupling to a time-dependent, spatially homogeneous electric field with tunable temporal width σ, photon frequency Ω, and linear or circular light polarization, called the "pump pulse".This setup is sketched in Fig. 1a.
We track the time-and momentum-resolved single-particle spectrum of the pump-driven electrons using a short "probe pulse" that emits photoelectrons and thereby generates a photocurrent 22,23 , as measured experimentally with tr-ARPES (see Methods).
To set the stage for our results, we briefly outline the basic ingredients for the lowenergy physics of Haldane's equilibrium model.We start from two Dirac cones with effective Hamiltonian v D (q x σ x ⊗ τ z + q y σ y ⊗ I).Here v D is the Dirac point velocity, the Pauli matrices σ label pseudospin arising from the graphene sublattices A and B, and τ labels the valley degree of freedom corresponding to the Dirac cones around K and K .Momentum q is measured from the respective Dirac points.In Haldane's model, an effective mass term mσ z ⊗ τ κ leads to an energy gap ∆ = 2m at the Dirac points (Fig. 1b).Its relative sign between K and K depends on its origin: If the gap is induced by introducing a staggered sublattice potential breaking inversion symmetry, τ κ = τ 0 , implying that the effective mass term has the same sign at K and K , and the out-of-plane pseudospin component P z is the same at both Dirac points (Fig. 1c).By contrast, if the gap originates from breaking timereversal symmetry, τ κ = τ z , hence P z points in opposite directions (Fig. 1d).Computation of the corresponding Chern number shows that the former state is a trivial band insulator with canceling contributions from K and K , while the latter is a nontrivial Chern insulator with additive Berry phases leading to a nonzero Chern number.
We now come to the discussion of our nonequilibrium results.We first characterize the nonequilibrium band structures using tr-ARPES spectra.Our method allows us to use arbitrary photon frequencies and driving fields A(t), with typical electric field profiles as shown in Fig. 3f and tunable pump-probe delay time t.We choose realistic optical pump pulses with photons at Ω = 1.0 and 1.5 eV, well below the electronic bandwidth of 17 eV.
Fig. 1 shows the tr-ARPES spectra on a momentum cut along K − K near K for before (t = −92 fs) and at peak field (t = 0 fs).We first perform a calculation using a pump pulse with linear polarization along the k y direction (Fig. 2a).This pump preserves time-reversal symmetry, and the spectrum remains gapless at the Dirac point energy.The only effect of the linearly polarized pump is a slight shift of the Dirac point location to the left, with a maximal shift occurring at peak field.
Next, we turn to circular light polarization, thereby breaking time-reversal symmetry.In Floquet theory, the quasistatic eigenvalue spectrum at finite driving field A shows copies of the original bands shifted by integer multiples of Ω, the so-called Floquet sidebands.Energy gaps of n-th order in the field open at avoided level crossings of sidebands which differ by n photon energies.For circular light, an energy gap of second order in the field opens at the Dirac point (see Methods).In our tr-ARPES simulation, for a moderate pump-field strength and 1.5 eV photons, an energy gap at K is induced, accompanied by avoided level crossing gaps nearby (Fig. 2b).Doubling the pump-field strength, the energy gaps become more pronounced, and clear Floquet-like sidebands appear (Fig. 2c).Lowering the photon frequency to 1.0 eV, the energy gap is smaller at peak field (Fig. 2d2) than for smaller field (Fig. 2d1), indicating a nonmonotonic field dependence of the gap.
In order to analyze the gap behavior in more detail, we now fix the momentum at K and analyze the tr-ARPES energy distribution curves (EDCs) near the Dirac point energy as a function of delay time (Fig. 3).For moderate field strength A max = 0.125, the gap continuously increases towards peak field for Ω = 1.5 eV (Fig. 3a), with some saturation already visible for lower photon frequency Ω = 1.0 eV (Fig. 3d).Increasing the maximum field strength to A max = 0.250, a slight decrease in gap size before peak field is visible for Ω = 1.5 eV (Fig. 3b).For Ω = 1.0 eV, the gap strongly decreases and even closes around peak field (Fig. 3e).Finally, gap closing and reopening is found for A max = 0.500 and Ω = 1.5 eV (Fig. 3c).
We observe that the maximal gap size that can be reached by starting from zero field and continuously turning on the pump is limited by the photon frequency Ω.On closer inspection, we find that the maximal gap scales linearly with the driving frequency.The The fact that a gap opens under circularly polarized light, which breaks time-reversal symmetry, and its absence under linearly polarized light, which preserves time-reversal symmetry, implies a close analogy with Haldane's model.This analogy suggests a nontrivial structure of the pseudospin texture for circular light polarization, as sketched in Fig. 1d.In order to demonstrate this, we now analyze the time-and momentum-resolved, band-specific pseudospin content for a particular gauge choice (see Methods).In Fig. 4b, we present pseudospin textures below the Dirac point energy near both K and K for the driven system at photon frequency Ω = 1.5 eV and at delay times corresponding to states p and d, with effective field strengths x eff (t) ≈ 0.1 and 1.4, respectively.For state p, the pseudospin com-ponents P x and P y are polarized along the k y and k x axes with a p-wave structure.P z peaks at K and K and, importantly, changes its sign between K and K .This behavior suggests a nontrivial quantum geometry with a mass term changing sign between valleys according to Fig. 1d.
After tuning the pump-field strength through the critical value x eff ≈ 1, the pseudospin textures shown on the right in Fig. 4b emerge.To start with, the P z components both at K and K change sign compared to state p, indicating that state d is also nontrivial.Intriguingly, the P x and P y textures develop additional nodes forming a d-wave pattern around the Dirac points, implying that the pseudospin vector has doubled its winding number when going around a Dirac point.This d-wave pseudospin texture matches the one expected for effective low-energy Hamiltonians which are quadratic in momentum measured from the respective Dirac point, (q x + iq y ) 2 , as obtained for bilayer graphene 24 in a perpendicular electric field 25 , and reminiscent of spin-orbital textures in topological insulators 26 .In this analogy, we stress that the effective "perpendicular electric field" generated by the circularly polarized light pulse points in opposite directions at K and K , in contrast to the static electric field applied in Ref. 25.This suggests the possibility to dynamically engineer effective Haldane multilayer models with complex pseudospin textures, as sketched in Fig. 4c (also see Methods).
Our combined results show that band gaps induced by breaking time-reversal symmetry in graphene are within reach under realistic experimental conditions.In particular, the achievable energy resolution for probe-photon energies which are sufficiently high to reach the Dirac points 27,28 should allow for the detection of photoemission gap sizes exceeding 100 meV.The change in pseudospin texture near critical driving offers the exciting opportunity of optical manipulation of quantum geometries and local Berry curvatures near Dirac points on ultrafast time scales.In a ribbon geometry, edge states might be induced 14 in addition to edge states already present for zigzag edges in equilibrium 29 , such that optical switching of nontrivial edge states appears as a possibility.Finally, the spectroscopic detection of pseudospin textures requires access to orbital band content.To this end, hexagonal structures with inequivalent orbitals on the A and B sublattices having different photoemission probeenergy cross sections could be examined.A candidate material for this purpose is MoS 2 .
Given its direct band gap of 1.8 eV at the Dirac point 30 , it seems unlikely that laser-induced band inversion can be achieved in MoS 2 , but the demonstration of pseudospin imbalance at the two Dirac points by circularly polarized light would already by intriguing.Alternatively, artificial hexagonal lattices with sublattice potentials have already been demonstrated with cold atoms. 31Thus the proposed pseudospin textures could in principle also be realized in driven ultracold quantum gases. 32

METHODS SUMMARY
The simulations presented here start from a minimal tight-binding model of spinless electrons with nearest-neighbor hopping on the honeycomb lattice 1,33 .The pump pulse drives the electrons via minimal coupling to a gauge field A(t) = A max p σ (t) (sin(Ωt)e x + cos(Ωt)e y ), with a Gaussian shape function p σ (t) = exp(−(t−t p ) 2 /(2σ 2 )) for a pulse of width σ centered around time t p , and photon frequency Ω.The phase shift of π/2 between the x and y components, represented by the unit vectors e x and e y , describes circular light polarization.
For comparison we also study linearly polarized light with vanishing x component.The The time-and angle-resolved photoemission spectroscopy (tr-ARPES) is computed from the trace of the nonequilibrium lesser Green function by a postprocessing step involving the probe-laser-pulse shape s W (t) with time resolution W , which leads to an effective tr-ARPES energy resolution ∝ 1/W 22,23 (see Methods for details).The delay time of a pump peaked at t p and probe peaked at t pr is given by t = t pr − t p .
The simulation parameters for the pump-probe setup are as follows: The pump laser field has frequencies ω = 1.5 (1.0) eV, implying oscillation periods of 2.58 (3.86) fs.Its temporal width is σ = 66 fs, unless indicated otherwise.We vary the peak vector potential A max = 0.125 . . .0.500, corresponding to peak electric field strengths 130 . . .530 mV/ Å for ω = 1.5 eV, and to 90 . . .350 mV/ Å for ω = 1.0 eV.The photoemission probe pulse has a width W = 20 fs.This choice of parameters is motivated by the hierarchy of time scales in the system: The oscillation period for the pump laser light, the temporal width of the probe pulse which controls the time and energy resolution for the tr-ARPES signal, and the temporal width of the pump pulse which controls the nonequilibrium state and ensures a well-defined center frequency for the pump pulse (see Fig. 3f).
Simulations are performed for an initial equilibrium sample temperature T = 232 K.We typically use 150,000−500,000 time steps for the computation of the time evolution operator, depending on the maximal simulation time and field parameters.This typically corresponds to a time step size of 0.0017 fs.Green functions are evaluated on a grid with 1,500−5,000 real-time steps and a step size of 0.17 fs.The Dirac point velocity is given by v D = 4.2 eV a C−C , where a C−C = 1.42 Å is the carbon-carbon distance 1 .We choose a chemical potential µ = 0.5 eV, which sets the Dirac point energy E D = −0.5 eV relative to µ.This choice is motivated by the fact that typical graphene samples on substrates are doped, and that states in both the lower and upper Dirac points are occupied in the initial equilibrium states and therefore nicely visible in the tr-ARPES spectra.

Keldysh Green functions. We define the Keldysh-contour time-ordered (T ) Keldysh
Green functions where the α and β operators are fermionic annihilation and creation operators in orbital or band basis, see below.Including the field via Peierls substitution, the time-dependent Hamiltonian for A and B sublattice orbitals reads with the Hamiltonian matrix elements g(k We define the Hamiltonian matrix for momentum k in orbital basis, In the absence of a driving field, this Hamiltonian is diagonalized by a rotation R(k) at t = 0, which is a time before the pump pulse is turned on.For later times, the given rotation does not diagonalize the Hamiltonian except for accidental cases where the gauge field is an integer multiple of a reciprocal lattice vector.The time-dependent Hamiltonian matrix in band representation is then given by such that h band (k, 0) is diagonal by construction.Note that t = 0 is used here as a notation for the earliest real time we consider on the Keldysh contour, not to be confused with delay time 0, which refers to the time where the pump-pulse envelope is maximal.
The computation of double-time propagators requires the evaluation of the time evolution operators Since H(k, t) at different times do not commute with each other, the time ordering T in U(k, t, t ) is taken into account by discretization of the real time axis and multiplication of the resulting time-step evolution operators.We then obtain the time evolution operator as 2 × 2 matrices in band basis, where N t,t is the number of fine time steps of size ∆t between t and t , and the product is understood as time-ordered with later times t − j∆t to the left.Computation of the lesser Green function matrix elements in band basis results from 34 where time "0" refers to an initial time where the system is in equilibrium before the pump pulse is turned on, and the corresponding time-independent matrix of band occupation represented by a † k a k is diagonal with Fermi functions for the corresponding equilibrium eigenvalues, f ( (k)).These Green functions are then rotated back to orbital space by the inverse of the transformation R(k) in order to compute the photocurrent and pseudospin content: tr-ARPES formalism.The computation of the time-resolved photocurrent involves normalized Gaussian probe-pulse shape functions s W (t) of width W centered around time t.The photocurrent is then obtained from 22 Here is the gauge shift for a double-time single-particle propagator 35 .In the main text, we also use the effective probe-averaged field strength given by the convolution of probe (width W ) and pump (width σ).For Gaussian shapes, this expression simplifies to where t is measured relative to the center of the pump, and σ2 = σ 2 +W 2 .Thus the effective field "seen" by the probe at zero delay time is reduced by the ratio σ/σ.
Pseudospin content.Pseudospin content is extracted by expanding the Green function matrices in orbital representation in Pauli matrices, The respective pseudospin content P x,y,z (k, t) specifically for states below the Dirac point energy E D is obtained by computing the analogue of the tr-ARPES response (9) for G < x,y,z , with j ∈ {x, y, z}, integrating (13) in a frequency window between E D −Ω/2 and E D , chosen to represent states in the n = 0 manifold near the Dirac point, and normalizing the resulting vector for each momentum point.
We note that the continuity of pseudospin in momentum space requires some care with the choice of eigenvectors for the orbital-to-band-basis transformation R(k) around Dirac points, since the phase of each eigenvector is arbitrary momentum by momentum.In practice, we fix this phase by making the first element of each numerically obtained eigenvector real and positive.Moreover, in contrast to the Berry curvature obtained from pseudospin representations, the orientation of pseudospin vectors is basis-dependent.For the lattice problem, in-plane pseudospins can be rotated by moving the origin in the real-space unit cell, whereby the Hamiltonian matrix elements g(k) acquire a momentum-dependent phase factor.However, our conclusions about the relative orientations of pseudospins between Dirac points, or the flipping of pseudospins and lowered point-group symmetries upon increasing the driving-field strength, are independent of this basis choice.
Floquet eigenvalues.The observation of gap closing triggers the question whether this effect is ubiquitous to the full nonequilibrium simulation on a lattice, or whether its precursors are already present in Floquet theory for linear bands.Floquet theory assumes a continuous ac driving field with constant envelope and frequency Ω.If Ω is the largest energy scale in the problem, Floquet theory provides a compact description in terms of a small number of time-independent basis states.The matrix elements of the static Floquet Hamiltonian for a driven system with period τ = 2π/Ω are given by Here integers m and n label the number of photons absorbed in the m-th and n-th Floquet manifold.Specifically for Dirac fermions, we represent momentum q as a complex number q = q x + iq y and use minimal coupling to a field A(t) with circular polarization represented as A exp(iΩt).The effective Floquet Hamiltonian is a Floquet matrix of 2 × 2 matrices in orbital space.For a Dirac cone with ideal linear bands and velocity v D one obtains This structure leads to avoided level crossings of order m − n in the field at the intersection of sidebands derived from the m-th and n-th manifold.Specifically at the Dirac point, where q = 0, the dimensionless Floquet eigenvalues as a function of dimensionless field x = v D A/Ω are analytically obtained as The initial quadratic band splitting is followed by a series of level crossings, which eventually lead to gap closings accompanied by changes in the pseudospin textures as shown below.
Pseudospin textures with lower point-group symmetry.The pseudospin textures that emerge after gap closings can be qualitatively understood when applying a perturbative "k • p" expansion for the dressed states in a Floquet picture for linear Dirac bands.hybridize only via processes involving two or more photons.To second order in q, we obtain an effective low-energy Hamiltonian The associated pseudospin texture displays a d-wave pattern, indicative of the analogy of H eff to bilayer graphene with an applied perpendicular electric field. 25Further gap closings occurring at higher fields -though experimentally challenging -paint an even more intriguing picture: As the relevant low-energy dressed states require increasingly higher-order photon processes to hybridize, the low-energy k • p dynamics progressively maps to trilayer graphene (with f-wave pseudospin texture), quadruple-layer graphene (g-wave pseudospin texture), and further multilayer realizations of the Haldane model, as sketched in Fig. 4c.
field strength A max required to close the gap is smaller for lower photon frequency, also with a linear relation.Converting the scale of A into an energy scale by multiplying with the Dirac point velocity v D , we investigate the size of the dimensionless gap ∆/Ω at peak field (t = 0 fs) as a function of the dimensionless effective field x eff = v D A eff /Ω, where A eff is the probe-averaged pump pulse envelope (see Methods) at zero delay time.The resulting curve is shown in Fig. 4a.The gap closing occurs for a dimensionless critical field x eff ≈ 1. Remarkably, the data for Ω = 1.0 eV and 1.5 eV, and for two different pump pulses with widths σ = 66 fs and 99 fs, collapse onto a single curve in the frequency regime under investigation, showing the same scaling relation as the Floquet eigenvalues (see Methods).The critical field distinguishes between states p and d, before gap closing and after gap reopening.
))), and V = 2.8 eV is the nearest-neighbor hopping matrix element matching the graphene bandwidth and Dirac point velocity.In equilibrium, the Hamiltonian has two Dirac points at momenta K and K given by 2π/3 1, ±1/ √ 3 ≈ (2.0944, ±1.2092),where momenta are measured in multiples of the inverse of the carbon-carbon distance a C−C . 1 At the K point, the Floquet Hamiltonian H 0 = m Ωm(|A, m A, m| + |B, m B, m|) + v D A(|A, m + 1 B, m| + h.c.) is readily diagonalized by a set of dressed states |±, n = cos θ |A, n + 1 + sin θ |B, n that mix n-and (n + 1)-photon Floquet states.Here, θ = arccos Ω+Ω 2 is the mixing angle, with Ω = √ Ω 2 + 4A 2 the dressed-state splitting.In this basis the momentum-dependent perturbation V q = v D (q x + iq y ) m |A, m B, m| + h.c.couples dressed states |±, n and |±, n + 1 .To discuss the pseudospin texture, we derive an effective low-energy Hamiltonian in dressed-state basis.Around the trivial gap closing A → 0, the k • p term couples the conduction and valence states to linear order and gives rise to a conventional massive Dirac fermion with p-type pseudospin texture.Approaching the first nontrivial gap closing at finite A however, the relevant low-energy dressed states

FIG. 1 .FIG. 4 .
FIG. 1. Graphene with band gap.(a) Honeycomb lattice irradiated with circularly polarized light.(b) Band structure with a band gap at the Dirac points.(c) Trivial gap structure in Haldane's model: The pseudospin P z (k) points in the same direction at K and K .The winding number contributions cancel.(d) Nontrivial gap structure: P z (k) points in opposite directions at K and K , leading to a nonzero Chern number +1 or −1.