{'fraction_non_alphanumeric': 0.033161752496619776, 'fraction_numerical': 0.007400906136584671, 'mean_word_length': 4.231819310002482, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dirac cones, Floquet side bands and theory of time resolved ARPES

Pump-probe techniques with high temporal resolution allow one to drive a system of interest out of equilibrium and at the same time, probe its properties. Recent advances in these techniques open the door to studying new, non-equilibrium phenomena such as Floquet topological insulators and superconductors. These advances also necessitate the development of theoretical tools for understanding the experimental findings and predicting new ones. In the present work, we provide a theoretical foundation to understand the non-equilibrium behavior of a Dirac system. We present detailed numerical calculations and simple analytic results for the evolution of a Dirac system irradiated by light. These results are framed intuitively by appealing to the recently revitalized notion of side-bands. We find that, under the application of circularly polarized light, a Dirac point only ever splits into two copies of sidebands. Meanwhile, the application of linearly polarized light leaves the Dirac point intact while producing side bands. Our immediate interest in this work is in connection to time and angle resolved photoemission experiments, where we find excellent qualitative agreement between our results and those in the literature. However, our results are general and may prove useful beyond this particular application and should be relevant to other pump-probe experiments.


I. INTRODUCTION
One of the greatest triumphs in the last decade of condensed matter research has been the theoretical prediction 4-6 and subsequent experimental realization [7][8][9] of the topological insulator (TI). These materials are insulating in the bulk, while their edge plays host to topologically protected metallic modes with energies lying in the band gap of bulk states. The existence of these edgestates makes TIs of great fundamental and practical interest with applications ranging from quantum computation to spintronics. Moreover, the discovery of topological systems leads to a new classification of possible states of matter.
While many of the topological systems can be understood by non-interacting, clean systems at equilibrium the study of topological states is not limited to those. The effects of disorder, for example may drive a system in and out of a topological state. 10,11 It is therefore interesting to ask whether there is a knob that can be tuned to alter the topological properties of a system. One auspicious route towards the generation of a TI comes from considering time-periodic perturbations [1][2][3] . In these systems, a time-periodic perturbation, is applied to a topologically trivial system and drives into a non-equilibrium topological state. As continuous time translational invariance is broken, it is no longer appropriate to discuss energy eigenstates. One must instead talk about their quasi-energy spectrum, which is the closest analogue to an energy spectrum for a system with discrete timetranslational invariance 33 . The topological state created with an external, time-periodic perturbation is called a Floquet topological insulator (FTI) and it exhibits edge-states in the gap of its quasi-energy spectrum 12 .
The notion of a FTI has garnered much attention lately, and has enjoyed experimental validation in the field of photonic crystals, where Floquet states can be simulated in the laboratory 34 . However, a solid state verification of a Floquet topological state and several issues regarding feasibility have been raised 35 . The first of these is that most available periodic perturbations are not perfectly periodic, but have an envelope function in addition to the periodic signal. The second, perhaps more pressing, issue has to do with the experimentally available frequencies. Present discussions in the literature are valid in the large frequency limit Ω Ω BW , Ω being the applied frequency and Ω BW being the frequency of the band-width of the system. However, available technology in terahertz is sub-bandwidth. This small frequency is believed to be problematic as it will lead to a complicated quasi-energy structure which may obscure any potential topological effects.
Given the above complications, our goal in the present Paper is to understand the behaviour of a topological system in the presence of a non-periodic and sub-bandwidth external perturbation, while probing the system continuously over time. We will work with Dirac cone dispersion, typical for a three dimensional topological insulator surface, in order to develop a fundamental understanding and will not discuss a Floquet topological insulator at this point. Our work is inspired by measurements of Wang et al in Ref. [3]. This group used time resolvedangle resolve photoemission spectroscopy (TR-ARPES) to view the evolution of surface states of Bi 2 Se 3 , a three dimensional topological insulator. We find excellent qualitative agreement with these results . We employ the language of sidebands recently used arXiv:1603.09718v2 [cond-mat.str-el] 8 Nov 2016 in Refs. [1,2]. This language allows us to develop the following physical picture of the time-dependent system. The pump pulse excites the system out of equilibrium. Its time dependence is generally composed of fast oscillations modified by a slow envelope function. The fast oscillations normalize the band dispersion and produces side band copies. The side bands are populated statistically with weights which depend approximately as Bessel functions on the ratio pump pulse amplitude to its frequency. Since the pump amplitude is time dependent through its envelope the side band weights are also time dependent. This time dependence allows the system to interpolate between its equilibrium state at the distant past to the Floquet/side band picture when the pump is applied.
The above picture leads to the three main results of this work. The first of these is that even though applying a sub-bandwidth perturbation to a system may "fold" many states into the Floquet zone, only a few of these states have any spectral weight and contribute to physical processes. In the present example, we consider a Dirac cone, which has effectively an infinite band-width. We find that only states within a couple Ω from the Fermi surface have any statistical weight in our side-band picture. Second, we work in a regime where the time scale over which the pump pulse envelope is changing is much longer than the period time of the drive. In this regime we develop simple, analytic expressions. In other cases (such as a quench) the same formalism can be applied and solved numerically. It should be noted that when the pump electric field is turned off, the system does not necessarily relax immediately to its equilibrium state. This is particularly clear in the case of the physical gauge we adopt, as explained in Appendix C and may lead to interesting effects like persistent Hall response 36 . Finally, we show that the structure of a Dirac cone colludes with circularly polarized light to produce only two sidebands for momenta near the Dirac point. This is quite remarkable; the spectral weight of the equilibrium Dirac point states is entirely shared between two sidebands. It therefore behaves as two massive Dirac points, with different masses. These two cones share the spectral weight of the original Dirac cone, and the weights are found analytically as a function of time.
The intuition developed here, as well as the satisfactory results in view of recent experiments 3 , will add to a current ongoing discussion in the literature regarding the stability of Floquet-states [37][38][39][40][41] . Our side-band interpretation in concert with an understanding of Floquet states and the results of Ref. [3] provide an intuitive physical picture of the side bands and their probabilities.
The rest of this paper is organized as follows. In the following section we discuss some fundamentals of the Floquet formalism in order to introduce the side-band intuition of Refs. [1,2]. We move on to present our model and methods. In Section III we present our results and discussion for two polarizations of light. The appendices detail various technical aspects of the work.

A. Sidebands
We begin with a brief discussion of Floquet theory as it pertains to the language of sidebands. Consider the time-dependent Schrödinger equation where H(t + T ) = H(t) is a Hamiltonian with period T . Defining Ω = 2π/T the principle result of Floquet theory is that the steady states of the above system can be written as 33 where (H(t) − i ∂ t )|φ(t) = η|φ(t) and |φ(t + T ) = |φ(t) . The eigenvalues η are typically called the quasienergies. The quasi-energies are only unique up to integer multiples of Ω, as can be seen by noting that e inΩt |φ(t) is an eigenvalue of (H(t) − i ∂ t ) with quasi-energy η + n Ω and also meets the boundary condition |φ(t + T ) = |φ(t) . Thus all quasi-energies are defined within a first "Floquet zone", an interval of energies of width Ω. The center of this zone is, of course, arbitrary. The quasienergy spectrum in the first Floquet zone can be copied at integervals of Ω above and below to generate the full quasi-energy spectrum. We now introduce the side-bands. Since |φ(t) is periodic, we are free to express it as a discrete Fourier series |φ(t) = n e −inΩt |n . The full wave function reads The states |n are determined by solving the eigenvalue equation m (H n−m − n Ωδ n,m ) |m = η|n where H n = T 0 dt T e inΩt H(t). The intuitive picture we wish to take away from Eq. (3) is the following. In a time periodic system the steady states are a linear combination of definite energy states with energies η + n Ω and probability n|n . This follows from either inspecting Eq. (3) or by noting that the average energy over one cycle of the period reads For pedagogical reasons we no consider the application of the above theory to a time-independent system. To be more concrete, let's say we have an applied perturbation with frequency Ω but a vanishingly small amplitude. In this limit H n = δ n,0 H and the eigenvalue equation becomes (H − n Ω)|n = η|n . The solution to this system is |n = δ N,n |ζ with η = E − N Ω, where H|ζ = E|ζ and N is an integer that takes E and moves it into the first Floquet zone we have chosen for our problem. Thus when the time periodic fields are turned off the system is, of course, found in eigenstates of the static Hamiltonian. These eigenstates can, of course, be defined in a first Floquet zone, but if this first Floquet zone does not contain E (i.e. if N = 0 in the language above) then there exists a quasi-energy E − N Ω, but this state has zero probability of being occupied in the first Floquet zone because n|n = δ n,N . One must move to the N th Floquet zone where this state is occupied with unit probability. The morale of this exercise is the following. When working in the Floquet zone the quasienergy spectrum might be dense with folded bands. However, the 'occupation' of a given state (its weight in the time dependent wavefunction) maybe zero in the first Floquet zone, leaving only a few relevant states.
Starting from the above limit, as we turn on the time dependence there are two effects that take place. First, electrons beginning in the original eigenstates develop some probability to absorb or emit photons and thus their unit probability of being found in one Floquet zone gets smeared into other, adjacent Floquet zones. This creates "copies" of the original band structure analogous to those proposed by Tien and Gordon several decades ago 42 . Second, unlike the physics of Ref. [42], in our present system these side-bands can also be modified in a non-trivial way. This occurs when states corresponding to absorbing/emitting different number of photons hybridize. This can lead to important effects such as gaps in these side-bands opening.
The above interpretation is important when applying a probe of well defines energy to a time periodic system. An example for this situation is studies in Refs. [1,2] where we calculate the transport properties of periodically driven quantum well heterostructures. Namely, the edge-states in these systems, whether naturally occurring or driven, are split into side-bands. As a result, certain transport signatures of these edge states, for example 2e 2 /h conductance, are fragmented. In Refs. [1,2] we have discussed how a sum rule 24 can be used to salvage these transport signatures. This sum rule is rooted in the understanding that systems in a time-periodic field have their energy bands modified by the time-periodic perturbation and also that these bands are split into side-bands. Crucially, these side-bands are only occupied with a certain probability, and, for reasonable field strengths, this probability decreases with the separation in energy between the original energy eigenvalue and the side-band eigenvalue that we're interested in. Thus it is usually appropriate to treat only eigenstates within several multiples of Ω from the Fermi level.
The above observations are important to keep in mind when applying Floquet theory to look at the quasienergies by themselves. When the energy scale Ω is small compared to the band-width of the equilibrium model, the quasi-energy spectrum becomes very convoluted as many eigenstates are "folded" back into the Floquet zone. Making predictions based on this spectrum alone then becomes an arduous task. The discussion above, and the results to follow, illustrate that one must keep in mind that even though the quasi-energy spectrum may become complicated in this limit, only quasi-energies resulting from folding of energies within a few Ω of the Fermi energy contribute significantly to observables. The information about these probabilities is contained in the often ignored side-band states |n and their statistical weight.
Our model is a generic Dirac cone and no cut-off is considered, thus our effective band-width is infinite. We subject this system to terahertz frequency light Ω ∼ 30meV. Looking only at the quasi-energy spectrum of this system the Dirac cone will be folded back into the Floquet zone infinitely many times and would thus be meaningless. We therefore approach the system in a slightly different manner, while keeping in mind the sideband language discussed above. Provided that the field is turned on slowly compared to the frequency of the light, the system evolves into a state described by a splitting of its original bands into side-bands. In cases where the operator describing the external field commutes with the static Hamiltonian at all times, this side-band splitting is the only effect of the light, i.e. we see no hybridization and no gap opening. In all other cases there are additional modifications of the side bands. In either case, we see that for physical field strengths only the first couple of side-bands carry any spectral weight in these simulations, in spite of the fact that the system is subjected to low-frequency light. These central results of our work are summarized in the schematic in Fig. 1. This intuition should be relevant to related experiments on timedependent systems and will be crucial in driving a topological state with externally applied light.

B. Model Hamiltonian
We begin with the following Dirac Hamiltonian where v F is the Fermi velocity, σ i = σ i is a vector of Pauli matrices and µ the Fermi energy. The above Hamiltonian is immediately applicable to the surface of a three dimensional topological insulator (TI) and should also be relevant to graphene in the limit where any applied field doesn't induce intervalley scattering. We now envisage the above system irradiated by an electromagnetic field. To keep our theoretical model simple we assume this field is spatially constant over the sample size. This should be approximately true for the terahertz type radiation considered here where the wavelength of the light should be tens of microns 3 . We model this electromagnetic field as follows where E 0 is the amplitude of this pump pulse, T pump is the width of the pulse and E Ω (t) is the monochromatic Schematic picture of the main results of this paper. The original Dirac cone is split into side bands, with side bands further away from the original cone receiving less "weight". In the figure this is signified using lighter colours for less probable sidebands. Now, if the operator describing the time periodic field, V (t), commutes with the original Hamiltonian, H k (t), then this splitting is all that happens. If these two operators do not commute sidebands hybridize and the band structure becomes modified by, e.g., having gaps opened. component of the field. In this paper we consider two scenario's: (1) Linearly polarized light, with E Ω (t) = sin Ωtx and (2) circularly polarized light in which case E Ω (t) = sin Ωtx − cos Ωtŷ. We introduce the above field via minimal coupling, ignoring the Zeeman effect, as we expect the dominant contribution to come from the electron's orbital motion. We choose a Gauge such that the electric scalar po- where we have chosen in initial condition such that A pump (t) → 0 for t → −∞. This choice of initial condition is, of course, immaterial and represents the gauge freedom of the problem. We show in the appendix that within the formalism we use in this paper 35,43,44 this choice of initial condition does not change any of our observations.
Let us define the frequency scale associated with the pump pulse envelope Ω pump = 2π/T pump . We work in the limit Ω pump Ω in which case it is appropriate to write (see Appendix A for more details) whereẼ . The evolution of our time-dependent system is now described through a minimal coupling of the above pump field to our Dirac Hamiltonian via k → k − eA pump (t). Thus sidebandsthe (time-dependent) Hamiltonian we work with is as follows To complete our discussion of the models we must de-fine the probe pulse profile. For this we take the envelop probe where T probe is the width of the probe, assumed to be much shorter than the width of the pump, T probe T pump , and t O is the delay time between the pump and probe peaks. t O is effectively the time at which we are "viewing" the system. In the above model we have (arbitrarily) assigned t = 0 to be the time at which the pump pulse is maximal.
In our simulation we take experimentally relevant values for the parameters from Ref. [3]. Namely, we estimate v F 3.6 eVÅ, µ 300meV, Ω 120meV. For convenience we define ω F = µ/ . We take a pump-pulse with a full width half-max (FWHM) of 250fs (T pump 106.16fs) and a probe-pulse with 35 T probe = 26fs. Finally, to fully illustrate the conceptual power of our findings we take E 0 7.5 × 10 −3 V/Å, slightly exaggerated from the estimates of Ref. [3].

C. Photocurrent
A simplified picture of the technology involved in ARPES is to think of the experimental set-up as measuring the particle current of electrons ejected from the sample at a wave vector k, energy ω and time t O (relative to the pump maximum time). This measurement is called the photocurrent, I(k x , k y , ω, t O ). Typically this quantity involves complicated momentum, orbital, and time dependent matrix elements. To develop a solid understanding for this problem we will work under the assumption that these matrix elements are the same for all orbitals, momenta and times. Under this approximation the relevant quantity to calculate is 35,43,44 In the above G < k (t 1 , t 2 ) is the 2×2 lesser Green's function matrix of the system in spin space. It is obtained by evolving the equilibrium states of the original Dirac cone from the distant past to the present. It is defined as follows where c † kα creates and electron with momenta k and spin α.
Our theory relies on knowing the solutions to the Dirac equation at all times, as these states can be used to construct the above Green's function. We therefore define the states subject to the initial condition |ψ k,α (t → −∞) = |φ k,α where |φ k,α are the eigenstates of the equilibrium system labeling the chirality of the state. Once these wave functions are known the lesser Green's function of the system can be constructed (see Appendix B) where f (E kα ) is a Fermi function and |ψ σ k,α (t) is the spin σ component of the state |ψ k,α (t) .
The theory described in the rest of this paper involves determining G < k (t, t ) either analytically or numerically and then making use of Eq. (9) to estimate the results of a TR-ARPES experiment.

III. RESULTS AND DISCUSSION
Using the methodology outlined above, we present our results and interpretation of calculations relevant to TR-ARPES measurements. For the sake of clarity, we divide our discussion into two categories. First, we consider light polarized along the x direction of the sample. Next, we allow for circularly polarized light. Certain limits of these two set-ups can be solved analytically and crucial insight can be gained into the distribution of states in a non-equilibrium system. We begin with linearly polarized light.

A. Linearly Polarized Light
We consider an electric field along the x direction only. In this case a closed form solution to the Dirac equation can be found along the k y = 0 cut of the Brillouin zone: In this case the time dependent field commutes with the Hamiltonian for the chosen momenta and the wave functions can be written as × e ieαv F t tr dt Ax(t )/ |φ kx,0,α where t r → −∞ is the "turn-on" time for the field. This gives the Green's function × |φ kx,0,α φ kx,0,α | note that the above is independent of t r . We eventually need Tr G < kx,0 (t, t ) , where the trace is over spin degrees of freedom. This trace can be performed in any complete basis and becomes particularly simple when we choose the states |φ kx,0,α , which leaves Our discussion has been exact until this point. We would now like to make an approximation to simplify the above trace. We recall that A x (t) = E0 Ω e −t 2 /2T 2 pump cos Ωt and expand it in the limit T pump 2π/Ω. Integration by parts may be used to show that to leading order in E 0 Ω 2 e −t 2 /2T 2 pump sin Ωt − e −t 2 /2T 2 pump sin Ωt Using the above, the identity e ix sin Ωt = m J m (x)e imΩt and assuming the probe pulse is much shorter than the pump pulse gives the following result for the photocurrent (for technical details see Appendix D) where with A(t) = eE0v F Ω 2 e −t 2 /2T 2 pump . The above formula is our main analytic result for this part of the paper. It provides a nice picture of the side-band splitting that occurs in the presence of a periodic field. Owing to the nature of the applied field, which commutes with the Hamiltonian, none of the original bands are dressed. The exponent describes peaks not just at energy eigenvalues v F k x − µ, but also at integer values of Ω above and below this value. This indicates that there are copies of the original  In the distant past we again only see the Dirac cone, as the pump field starts to turn on we see copies of this cone (sidebands) begin to develop and these copies develop avoided crossings. As the field becomes full turned on the weight of the original Dirac cone is shifted into other sidebands and these gaps become more evident.
band structure at multiples of Ω above and below the original pattern.
The Bessel function pre-factor gives the weights of these side-band peaks. These weights depend on the probe time due to the time dependence of the driven system. The physical picture which emerges here is as follows. While the periodic part of the pump pulse is responsible for the existence of side bands and their dispersion, their relative contribution to the photocurrent is given by the above Bessel functions. Therefore, the ration of the pump envelope to its frequency determine the side band weights at any given time. In addition, it is evident from Eq. 18 that the probe pulse determines the time resolution, as the effective gauge field is a weighted average of the pump over the probe duration.
Owing to the simple structure at k y = 0 there is no interference/avoided crossing of sidebands. Thus in the limit of a wide pump pulse the system is split into sidebands and the population of these side bands is given by the instantaneous weighted average of the pump envelope function.
We now turn to numerics in order to test the validity of our analytic results and to extend our analysis to finite k y . For this we integrate the Dirac equation numerically. We begin by fixing k y = 0 and comparing our analytic treatment to exact numerics. Fig. 2 shows I(k x , 0, ω, t O ) for k x = 0 and k x = 0.05Å −1 for several values of t O . As can be seen in the figure, there is excellent agreement between our approximate formula above and the numerics. Fig. 2 also nicely illustrates the side-band interpretation discussed above. We see that all of the spectral weight associated with the original peaks in the distant past (before the pump pulse hits the system) gets redistributed into sidebands separated by Ω.
Next we move on to present results going beyond the scope of the analytic results. We plot I(k x , 0, ω, t O ) and I(0, k y , ω, t O ) in Figs. 3 and 4. First, the results for I(k x , 0, ω, t O ) (within the purview of the analytic approach above) nicely confirm the intuition developed above; we see no renormalization of the energy bands and a simple development of of sidebands. These sidebands are evident by the copies of the Dirac cone seen in the above plots. Second, I(0, k y , ω, t O ) goes beyond our analytic approach above. We see a twofold effect as the pump-pulse hits the system. The primary effect is a splitting of the system into side bands. Secondary, we see that the light renormalizes the side-band structure opening gaps in energies where level crossing occurs at equilibrium.

B. Circularly Polarized Light
beyond We now shift our focus to the more involved problem of circularly polarized light. Circularly polarized light makes even the k y = 0 cut along momentum space intractable analytically. We can, however, make progress right at the equilibrium Dirac point k x = k y = 0, the Γ-point. Here we have The above can be written as To solve for the evolution under this Hamiltonian we transform to a rotating frame by letting |ψ Γα (t) = e iµ(t−tr)/ e iΩtσz/2 |ψ Γα (t) . Our equation of motion for the wave function then reads Our purpose in finding |ψ Γα (t) is to build the Green's function G < Γ (t, t ) and ultimately convolve this Green's function with the probe pulse envelope. Therefore, a good first approximation would be to find the wave function in the vicinity of t O , the peak time of the probepulse. We therefore make the somewhat crude approximation A(t) → A ef f (t O ) in the above equation of motion, Eq. 22. This yields an effective (rotating frame) Hamiltonian which is time independent. The above equation of motion can be therefore solved to give: is an effective, time-independent, Hamiltonian in the vicinity of t O . |ψ Γα (t r ) is the rotating frame wave function in the distant past. By requiring that the wavefunction |ψ Γα (t) at t → t r be a simple spinor (whose direction can be arbitrarily chosen due to the degeneracy at the Dirac point) we find |ψ Γα (t r ) = e −iΩtrσz/2 |φ Γα (25) where {|φ Γα } are the eigenstates at the Dirac point in the distant past. Using the above, noting that the eigenvalues of , and performing some additional manipulations which are left for Appendix E we arrive at the following approximation for the photocurrent where α, β, s, s all run over ±1 and A s The gap betweem states is 2E ef f (t O )/ , which is consistent with the gap found in [45] in a purely Floquet system (T probe → ∞). Importantly, the distance separating the peaks in the Gassians above is much larger than the width of the peaks. We therefore discard terms where s = s . Further, we note that the eigenvalues in the distant past Γα = −µ are independent of α (as we're at the Dirac point). These two observations along with  27) for I(0, 0, ω, tO). The left shows I(0, 0, ω, tO) for different delay times tO while the right shows I(0, 0, ω, tO) also for different delay times but this time with a pulse FWHM of 500fs instead of 250fs. In these plots the bottom plot is for tO = −500fs, the middle for tO = −100fs and the top is for tO = 0fs. In all plots the solid line is the approximation in Eq. (27) while the circles are numerical results. On the top row we see reasonable agreement between numerics and our approximation for all three delay times. When we turn up the pulse width, which effectively makes the "turn-on" time slower, we see that the agreement becomes excellent. some additional straightforward, but tedious, algebra lead to the simplified result Examining the above shows that the ARPES spectrum from the Γ point shows the following features at energies E (measured from µ) with weights P . It is obvious from the above the there's no additional spectral weight in any other energy. When the amplitude A(t) is shut offâ z → 1 and we see E 1 , E 2 → 0 with weights going to unity. At the same time E 3 , E 4 → ± Ω, albeit with zero weight.
Our interpretation of the above is as follows. As the pump probe is turned on, the original two-fold degeneracy at the Dirac point is lifted and a gap is opened up with width the weight of these states is 1+âz 2 which decreases with field strength. The peaks at the other two energies correspond to single sidebands of the states E 1 and E 2 . Put another way, E 3 = E 1 − Ω while E 4 = E 2 + Ω. The weights of these side bands increase with field strength. Interestingly, unlike our treatment of the linearly polarized light, there is no statistical weight given to any other side bands, all of the spectral weight is found within two sidebands. Note that the same approximations were made in both cases. We can trace this phenomenon back to the transformation we made to the rotating frame. While in general, this transformation leaves the Hamiltonian time dependent, here it does not since we have used an effective field strength. In the rotating frame we find two solutions to our Hamiltonian and transforming back to the original frame can split each one of these into two side band due to the dimensionality of the transformation operator. It is interesting to note that the same behavior was found by Dehghani et al. in Ref. [30].
With the above analytic analysis let us move on to numerical methods in an effort to validate the above description and further explore momenta where an approximate solution is not tractable. We do this with the sideband language discussed above in mind.
We begin with a simulation at the Gamma point. Fig. 5 shows I(0, 0, ω, t O ) as a function of ω for various different values of t O . Both our approximate analytic expression as well as our numerics are displayed in this plot. We see that the approximation provided above is in good agreement with the numerics with respect to both the size of the gap and the position of the sidebands, it also shows that this approximation becomes better when the width of the pump-pulse gets larger.
Before continuing we would like to highlight the fact that the discussion here appears to be more general than the Dirac model we have used. To check that our conclusions are not simply a coincidence of this model we have gone beyond the Dirac cone model by including higher order corrections 46 and also by studying a lattice model for TIs. Our conclusions of only two side-bands do not change. In order to keep the discussion of the main text simple we have included the details of calculations on these models in the appendix.
Let us now move on to explore a wider range of momentum using our numerical protocol. Fig. 6 shows the time-evolution of the ARPES spectrum for a cut such that k y = 0. A cut along k x = 0 looks very similar and such plots would not add to the present discussion. In the figure we see effects common to all results in this work. As the field strength is turned on the original Dirac cone is copied into sidebands, each of which is populated only with a certain weight. States in these sidebands then hybridize with each other leading to gaps. Most notable, our analytic result for the Γ point is verified at the center of the momentum cut.

C. Timescales
We close this work with a short discussion on the timescales required to see the side-band physics that we have discusses here. In the majority of this paper we have focused on a hirarchy T pump T probe 2π/Ω. Focusing on these timescales was not only relevant from an experimental viewpoint 3 but also aiding in our derivation of approximate analytic results. Here we will briefly explore what happens when these conditions are relaxed. First, we have explored the effect of varying T pump on the development of side-bands. We have found that varying T pump down to even half of T Ω = 2π/Ω one can still see the development of side-bands. Results of this can be seen in the top plot of Fig. 7. We see that the major effect of varying T pump is that the effective field strength that the system sees is decreased (visible by noting the smaller gap and side-band amplitudes). We can understand this heuristically as the system having less time to see the field while the envelope is at large values. Second, we have studied how changing T probe can change our observations. Not surprisingly, 1/T probe sets our energy resolution in the photocurrent. A very narrow probe width leads to very broadened side-band peaks that can overlap with each other. For sharp side-band peaks the width of the probe field should be made as large as possible. The intuition behind this appears to be that the probe field needs to observe the system for at least a few periods, T Ω , to properly observe the sidebands. Our results are summarized in the bottom panel of Fig. 7. In the bottom plot we see that the probe pulsewidth sets our ability to resolve the side-band peaks.

IV. CONCLUSIONS
We have provided both simple analytic results and rigorous numerical simulations of TR-ARPES in a Dirac system. Our results show that the time-evolution of an ARPES spectrum can be understood using the language of probabilistic occupation of sidebands in a time periodic system 1,2 . The population of these side bands depends non-linearly on the envelope of the pump field. As expected, when the pump is turned off gradually the side bands are not populated while the original band is populated with probability 1. The shape of the probe pulse determines the time resolution and therefore, when the probe is sharp the TR-ARPES photocurrent follows the time evolution of the side band population. Our results are in qualitative agreement with those of the experimental work in Ref. [3].
Our work also highlights the fact that not all side bands are equally important. We showed that despite the repeated folding of the Dirac cone into the Floquet zone only a few side bands, which are displaced by a few Ωs from the equilibrium energy, contribute to the time resolved ARPES signal. This point pertains not only to the results of Ref. [3] but also to other measurements on Floquet topological systems such as transport 1,2 .
Finally, we have explored an interesting interplay between a Dirac point and circularly polarized light. Our results suggest that under the application of circularly polarized light the Dirac point is gapped and copied into two sidebands only. This is in contrast to other systems, and other points in the Brillouin zone, where many sidebands can be seen.

V. ACKNOWLEDGMENTS
The authors are thankful for useful discussions with Nuh Gedik. Financial support for this work was provided by the NSERC and FQRNT (TPB) and the Vanier Canada Graduate Scholarship (AF). Numerical calculations for this work were performed using McGill HPC supercomputing resources.
As discussed above, we choose to describe this electric field in a gauge where the scalar potential is zero. Thus we have where we have chosen in initial condition such that Let us define the frequency scale associated with the pump pulse Ω pump = 2π/T pump . We work in the limit such that there are many oscillations within the pump field envelope. We now define E Ω (t) = −Ẽ Ω (t) Ω and integrate by parts to obtain pumpẼ The above process could in principle be iterated to produce a perturbative expansion in Ωpump Ω , although we stop here for practicality. We could alternatively write Ωpump Ω = T Tpump , which tells us this expression is valid in the limit T T probe ; i.e. the amplitude changes on a much longer time scale than the period of oscillation. We neglect all but the leading order terms. Continuing the procedure above shows that the next to leading order term is of order Ωpump Ω 2 . Thus we work in a regime where where h k is the Dirac Hamiltonian. Therefore H k+eA0 (t) is independent of our arbitrary choice of A 0 , meaning that U k+eA0 (t, t r ) and thus G k+eA0 (t, t ) are gauge invariant as well.
Of course the most natural choice is to set A 0 = 0 so that the (equilibrium) Hamiltonian before the switch on time t r is simply h k , as one would like it to be. For this choice of initial condition the Gauge invariant Greens function and the traditional Greens function are identical.

Appendix D: Linearly Polarized Light
Here we go over the details leading to Eq. (18) in the main text. Assuming T pump 2π/ω, Ωt − e −t 2 /2T 2 pump sin Ωt and using the identity e ix sin Ωt = m J m (x)e imΩt in the wave functions for the linarly polarized light the photocurrent with k y = 0 reads In the above, J m (αA(t 1 )) describes a splitting of the eigenstates into sidebands, labeled by m, where the amplitudes of these sidebands depend on time. Meanwhile is the profile of the probe pulse and e −i(ω+αv F kx+ µ −mΩ)t1 describes having energies not just at ±v F k x − µ, but also at values m Ω above and below these values. The integral above is of course intractable to perform exactly. We can make progress using a series of appropriate approximations. The first is that the probe pulse is much shorter than the pump pulse. Therefore A(t) changes very slowly over the dura- probe . We can therefore simply replace the J m (αA(t 1 )) term with its value at the peak t 1 = t O . However we observe that a better approach is to replace A(t) = eE0v F Ω 2 e −t 2 /2T 2 pump with a weighted average over the probe pulse. Thus we define Once we have made this replacement the remaining integral can be done analytically and gives The width of the peaks described by the Gaussians above are set by the frequency scale 1/T probe . The two Gaussians describe peaks centred at αv F k x − µ + mΩ and αv F k x − µ + m Ω. Thus the separation between the two peaks is (m − m )Ω. If the decay scale 1/T probe is much smaller than the smallest separation Ω, i.e. 1/T probe Ω then the peaks do not overlap at the same frequency and so the major contribution to the double sum comes from the m = m terms. Working in this approximation gives Appendix E: Circularly Polarized Light

Analytic Work
We begin with the approximate wave function found in the main text, reproduced here for convenience |ψ Γα (t) = e iµ(t−tr)/ e iΩtσz/2 (E1) × e −iH ef f (t−tr)/ e −iΩtrσz/2 |φ Γα . We note that in order to find the trace of the Green's function we will require φ Γ,β |ψ Γα (t) (and its complex conjugate). In the distant past the field is turned off and so the Hamiltonian is h Γ (t → −∞) = −µσ 0 . Thus we are free to choose any "initial" set of states, provided they are orthonormal. For convenience we choose {|φ Γα } to be (1, 0) T and (0, 1) T , which we label |φ α with α = ±1.
From this point forward we will drop the Γ subscript in the interest of brevity. This leads to the following φ Γ,β |ψ Γα (t) = e i(µ+β Ω/2)t/ e −i(µ+α Ω/2)tr/ In order to calculate the matrix elements of e −iH ef f (t−tr)/ we write the argument of the exponential as follows is a unit vector and we remind the reader that Making use of the identity e −ixâ· σ = cos (x) σ 0 − i sin (x)â · σ one can show that where A s α,β = (δ α,β − sâ · ψ β | σ| ψ α ) /2 which leads to φ Γ,β |ψ Γα (t) = Recalling that we are working under the assumption ≥ Ω/2 1/T probe , the distance separating the peaks in the Gaussians above (which is 2E ef f (t O )/ ) is much large than the width of the peaks. We therefore discard terms where s = s which gives us

Beyond the Dirac Cone Model
One of the main results of this work is that the Dirac cone, under the application of circularly polarized light, admits only two side-bands at the Dirac point. Here we go beyond the Dirac cone description in order to argue that this behavior is universal and not a peculiarity of the linearized Dirac cone model. Towards this end we have added so-called hexgonal warping and velocity renormalization terms to our model 46 , our equilibrium system then reads h k = ( v F + Λk 2 )(k × σ) + λ 2 (k 3 + + k 3 − )σ z − µσ 0 (E9)  where Λ and λ characterize the hexagonal warping and velocity renormalization terms respectively. In the above we have defined k ± = k x ± ik y . We use as estimates for these parameter values estimated from experiment, λ = 50eVÅ 3 and Λ = 100eVÅ 3 . Our numerical results are shown in Fig 8. We see that the inclusion of these terms does not change the conclusion that only two sidebands exist. We have run simulations at λ and Λ values two orders of magnitude larger than the physical values above and still find the same behaviour. In addition to the above continuum model we have also used the lattice model of a TI developed by Marchand and Franz 47 to look at the behaviour at the Dirac point. The equilibrium Hamiltonian is given by: ξ k = 2λ (sin(k x )σ y − sin(k y )σ x ) M k = −2t (cos(k x ) + cos(k y )) − µ (E10) We calculate the photocurrent again, in the framework of this model. While the band curvature away from the Dirac point it definitely apparent, the behaviour at the Dirac point is essentially the same as in the linearized model. Fig. 9 is a snapshot of the time result ARPES spectrum at the center of the pump pulse, for k y = 0 and both circularly and linearly polarized light.