Revealing Ultrafast Phonon Mediated Inter-Valley Scattering through Transient Absorption and High Harmonic Spectroscopies

Processes involving ultrafast laser driven electron-phonon dynamics play a fundamental role in the response of quantum systems in a growing number of situations of interest, as evidenced by phenomena such as strongly driven phase transitions and light driven engineering of material properties. To show how these processes can be captured from a computational perspective, we simulate the transient absorption spectra and high harmonic generation signals associated with valley selective excitation and intra-band charge carrier relaxation in monolayer hexagonal boron nitride. We show that the multi-trajectory Ehrenfest dynamics approach, implemented in combination with real-time time-dependent density functional theory and tight-binding models, offers a simple, accurate and efficient method to study ultrafast electron-phonon coupled phenomena in solids under diverse pump-probe regimes which can be easily incorporated into the majority of real-time ab-initio software packages.


I. INTRODUCTION
Time resolved spectroscopies, such as time and angle resolved photoemission, time resolved photoluminesence and transient absorption spectroscopy (TAS) constitute fundamental tools to study the flow of energy in materials following excitation by light.Understanding the microscopic details of the excitation and relaxation pathways can serve as the basis for deterministic manipulation of material properties for technological applications such as enhanced photodetectors [1,2], long lived optically controlled qubit registers [3,4] and attosecond control of magnetic ordering for ultrafast spintronics [5][6][7].In parallel to developments pushing the time, energy, and momentum resolution of these spectroscopic techniques, there has been a plethora of phenomena studied under novel conditions such as transient phases and Floquet renormalization under strong parametric driving [8][9][10][11][12], chemical reaction rate modification under exposure to cavity confined fields [13,14], and exotic quantum phases when layering 2D materials [15,16].The study of the microscopic origins of these phenomena pushes the boundaries of theoretical tools which are useful near equilibrium conditions, in particular for one of the most fundamental processes for understanding the behavior of materials: the electron-phonon interaction.
Treating coupled electron-phonon dynamics in simula-arXiv:2306.16010v1[cond-mat.mtrl-sci]28 Jun 2023 tions of periodic systems is typically limited to either phenomenological coupling and decay terms, or coarse approximations such as the two-temperature model [17][18][19].Going beyond these options leads one to consider an explicit treatment of the phonon degrees of freedom, which can be achieved through the time-dependent Boltzmann equation (TDBE) [19][20][21][22][23][24], which is derived in the perturbative limit of the electron-phonon interaction.Attempts to move beyond some of the constraints of the TDBE include approaches based on density matrix formalism, such as the Bloch equation [25][26][27], the Hierarchical Equations of Motion [28] or the time-dependent Density Matrix Renormalization Group [29,30] for 1D systems, and non-equilibrium Green's function approaches based on solving the Kadanoff-Baym equations [31][32][33].
Conversely, the electronic system can be treated in a fully ab initio manner and one can still capture the effect of phonon fluctuations, even in cases with very strong coupling, via static displacement approaches which sample phonon distortions in large supercells [34,35].In the adiabatic limit, this approach has recently been shown to be equivalent to the Feynman expansion to all orders of the electron-phonon interaction perturbation [36], and is analogous to the nuclear ensemble average technique from molecular physics, where electronic properties are obtained by averaging over the nuclear coordinate distribution on the ground Born-Oppenheimer state [37].Breaking the high symmetry equilibrium lattice structure of periodic systems has been found to significantly alter the results of ab-initio calculations; disorder has been argued as being responsible for capturing a significant portion of what is usually attributed to electron correlation when using ostensibly uncorrelated DFT methods [38], as well as being the dominant factor in phase transitions typically argued to be rooted in electronic structure [39].
A natural step beyond including the static disorder of the nuclei is to also consider their dynamics, which opens the possibility to treat time-dependent response properties where the nuclear forces are dependent on nonequilibrium electronic configurations.A mean field treatment that can be applied in this context is the multi-trajectory Ehrenfest approach (MTEF), which allows one to recover the quantum statistics of the equilibrium nuclear subsystem [40] whilst approximating the time-evolution using Ehrenfest trajectories.This approach has been shown to capture Franck-Condon physics [41] as well as time-resolved out of equilibrium system dynamics [42][43][44].Propagating the system in real time allows for coherent electronic evolution at short time scales, while accounting for all orders of interaction with both external fields and the phonon system (at the mean field level).
Although the basic ingredients required for MTEF are already available in most real-time ab initio codes, and despite some existing formulations of reciprocal space semi-classical dynamics in the literature [45][46][47], to the authors' best knowledge, application of this method has not been widely explored in periodic systems, nor has the static displacement approach been widely applied to time-resolved phenomena.Rather, most uses of ab initio semi-classical dynamics in periodic systems involve only a single trajectory [48][49][50], typically using a primitive unit cell or initializing the phonon distortion with classical molecular dynamics simulations which often fail to capture the exact quantum statistics of the initial state [51].Nonetheless, inclusion of more phonon modes through supercell dynamics allows for detailed study of fundamental processes such as the relaxation of excited electrons through phonon emission, anharmonic phononphonon scattering and phase transitions, even at the single trajectory level [52,53].
Therefore, in this work we formulate the Wigner representation for the phonon subsystem in a generic manner using ab initio dispersion relations of real materials, yielding a framework that systematically captures the equilibrium properties of the phonon system.While MTEF is well known to suffer from zero point energy (ZPE) leakage and incorrect thermalization between the quantum and classical systems at long time scales, there are several schemes that can systematically correct these failures with added computational cost [54,55].However, as our focus here is limited to studies of the short time-scale dynamics of strongly driven systems, we will not pursue such corrections and simply aim to discover what can be achieved at the mean field level.
As an illustrative example we study ultrafast phonon mediated electronic reorganization following valley selective laser excitation in monolayer hexagonal boron nitride (hBN), demonstrating the ability of MTEF to capture many of the essential features of the process which has been argued to be a significant driver of valley selective relaxation in Transition Metal Dichalcogenides (TMDs) [31,[56][57][58].Due to negligible spin-orbit coupling in hBN, the two K valleys in the Brillouin Zone (BZ) are degenerate, while valley specific selection rules are preserved upon interaction with circularly polarized light [59][60][61][62], making this a prototypical material to study the relaxation of excited charge carriers theoretically.Using a real space ab initio supercell approach with time-dependent density functional theory (TDDFT), and a reciprocal space tight binding (TB) model, we find that including phonon fluctuations leads to a rapid redistribution of excited charge carriers within a characteristic time scale of less than 30 fs, and that these results converge with as little as two trajectories.
We show that in the static limit, the harmonic Wigner distribution of phonon momenta and coordinates reduces to the phonon distribution obtained in the reciprocal space picture of Williams-Lax theory [35].We compare the results of propagating the electronic system with MTEF to the limit of frozen phonons, as well as the TDBE approach, finding broad qualitative agreement between these approaches across temperatures from 0−2000 K while reproducing the experimentally observed lowtemperature behaviour of analogous TMD systems.Finally, we demonstrate the flexibility of our method to predict and recreate spectroscopic experiments by simu-lating two different transient absorption measurements.First we propose a circularly polarized TAS measurement and demonstrate the signal corresponds directly to the valley asymmetry decay.Second, we replicate a recently performed study [63] using the ellipticity of harmonics generated under extreme laser driving to track tunable valley selective excitation in hBN using bichromatic counter-rotating 'trefoil' pumps, also showing a rapid decay of the observable signal.Both cases demonstrate the ease with which our approach can be applied within tight binding and fully ab initio real time electronic dynamics software packages under arbitrary pump-probe setups.

II. THEORY
Here we briefly summarize how the description of the phonon subsystem can be formulated based on the phonon dispersion of real materials within MTEF by appealing to the Wigner distribution of the phonon system.We then elaborate on the connections between this approach and the static displacement formalism, and we briefly outline the TDBE that will be used for later comparison.We generally use atomic units throughout, though sometimes for clarity is written explictly.

A. Multitrajectory Ehrenfest dynamics
While a variety of routes to derive the Ehrenfest equations of motion are available, we focus here on how the MTEF approach can be derived as the uncorrelated solution to the quantum-classical Liouville equation (QCLE) [64,65], which describes the approximate time-evolution of the total density matrix.The partial Wigner transform is employed, which for an arbitrary operator is defined as: Here the list of variables (R, P) represents the full set of nuclear position R = (R 1 , . . ., R N ) and momentum P = (P 1 , . . ., P N ) variables which are vectors in N × d Cartesian dimensions, and we note that the operator character of objects in the electronic Hilbert space is unchanged by the partial Wigner transform.
In the mean field limit the density operator can be factorized ρW = ρ n,W (R, P, t)ρ e (t), and the Wigner function of the nuclear degrees of freedom can be represented by an ensemble of N t independent trajectories, ρ n,W (R, P, t) = Nt i=1 w i δ(R − R i (t))δ(P − P i (t)), with weights w i .The time evolution of the electronic density and the phase space coordinates is given by the Ehrenfest equations of motion [40]: where M are the nuclear masses.These equations are solved for the independent trajectories with initial conditions (R i (0), P i (0)) sampled from ρ n,W .Observables are constructed by averaging over the ensemble of trajectories; in the case of equal trajectory weights,

B. The Phonon Subsystem
For completeness we restate some textbook definitions of the phonon coordinates, in particular drawing from the work of Brüesch [66] and Giustino [67], with special emphasis on the often less well described conjugate momenta, which play an important role in the MTEF method (however, for a notable exception see [68]).
For a supercell composed of N p = N 1 × . . .× N l primitive cells in l periodic dimensions, we utilize the Born von-Karman (BvK) boundary conditions (see SI.1 A).Each primitive cell contains N c unique atoms which for a given lattice configuration have equilibrium positions in the primitive cell of R 0 α , for α = 1, . . ., N c .The equilibrium position of a given atom α within the supercell is specified by the primitive cell position R p ; for primitve cell index p = 1, . . ., N p , and the primitive cell equilibrium position, R 0 αp = R p + R 0 α .We further denote small displacements from these positions via δR αp = R αp − R 0 αp .Canonically conjugate momenta, P αp , can then be defined by the commutation relation Using the textbook definition of the interatomic force constant matrix and it's Fourier transform, the dynamical matrix (see SI.1 B), we define the complex normal coordinates and momentum for a given phonon quasimomentum q and branch ν as the following linear transformation : Here, e αν (q) ∈ C d is the normal mode of vibration describing the displacement of atom α in the primitive cell with mass M α .M 0 is a reference mass, taken to be the mass of the proton.The inverse of this transformation reads as p qν e iq•Rp (M 0 /M α ) 1/2 e αν (q)z qν It is easily shown that the complex normal position and momenta obey the normal commutation relations for phonons [66][67][68][69]: [z qν , P −q ν ] = i δ qq δ νν . (5) The redundancy in the complex normal coordinates, seen by z −qν = z * qν , can be removed by introducing the so called real normal coordinates [66,67].To show how this is done, we start by partitioning the q grid in the first BZ into three sets, in the manner of Giustino and Brüesch.Call set A the set of vectors invariant under inversion modulo addition with a reciprocal lattice vector G, i.e. q = −q + G .Set B and C are partitioned in such a way that all vectors q ∈ C are obtained from q ∈ B via inversion, i.e. q = −q + G.This leads to the following definitions whereby the following properties hold true Now we can rewrite (4) as where the sum over C has been included by taking twice the real part of the sum over B. With this removal of the redundancy, one can see that this linear transformation contains dN c N p independent variables corresponding to x qν for q ∈ A and x qν , y qν ∈ B, with canonical conjugates r qν for q ∈ A and r qν , s qν for q ∈ B, defining canonical commutation relations Following Giustino we define the characteristic length scale of the phonon frequencies as for q ∈ B, C and rescale the coordinates as zqν = z qν /l qν and Pqν = P qν l qν / .With this linear transformation we use the properties of the normal modes subject to the BvK boundary conditions to rewrite the real space nuclear Hamiltonian in the harmonic limit as: 1. Phonon Wigner function The Wigner transform of the density matrix of a set of uncoupled phonons in the canonical ensemble can be written in terms of the reduced coordinates as [70]: for β = 1/k b T , where the reduced coordinates are now treated as continuous degrees of freedom.We use this distribution to sample the phonon modes when performing real space supercell calculations; the nuclear configuration associated with a particular phonon coordinate configuration is obtained by simply using Eq. ( 9) after sampling the reduced coordinates from Eq. ( 13).
The required inputs are e αν (q) and ω qν , which are easily obtained via Density Functional Perturbation Theory (DFPT) calculations such as the implementation in Quantum Espresso [71,72].
We note here that while sampling the equilibrium phonon distribution in Eq. ( 13) exactly captures the quantum statistics of the phonon system only in the harmonic limit, the time evolution is not constrained to this limit as the nuclei are subject to the (Ehrenfest) forces coming from the driven electronic system.In this sense, the phonon coordinate picture can be viewed as a convenient basis in which the nuclear system can be initialized, rather than a limitation of MTEF to the harmonic approximation.

C. Connection to Static Displacement Methods
In the special displacement method (SDM) introduced by Zacharias and Giustino [34,35], static thermodynamic properties of periodic systems (within the Born-Oppenheimer approximation) are calculated by taking specific phonon coordinate configurations from the Williams-Lax nuclear coordinate distribution, written in reciprocal space [35].For an arbitrary operator Ô the expectation value at a given temperature is given as with σ 2 qν = l 2 qν (2n qν,T + 1), for the Bose-Einstein occupation of the mode, n qν,T = [exp (β ω qν ) − 1] −1 , and O {xqν ,yqν } (T ) referring to the expectation value of Ô at temperature T with respect to the electronic system, evaluated at phonon configuration {x qν , y qν }.
The connection to the Wigner transformation of the phonon coordinates can been seen by inserting σ qν and n qν,T for the distribution functions of Eq. (14).Keeping in mind our definition of l qν in Eq. (11), one immediately arrives at Eq. ( 13) for the position coordinates.Therefore, including the phonon momenta is in some sense an extension of this method to dynamic ions.Thus, given the success of the SDM in capturing phonon renormalization using small numbers of displacement samples, one can consider sampling approaches inspired by this method which incorporate momenta.We explore this possibility in SI.6, but find no significant advantage for the observables studied here.
In what follows, we refer to a simulation protocol as 'static' when we sample exclusively positions from Eq. (13) (i.e. from the distribution function in Eq. ( 14)), constrain the phonon coordinates at this initial configuration and evolve only the electronic system in time.This approximation is sometimes referred to as the 'clamped ion', 'frozen phonon', or 'static disorder' approximation.In contrast we refer to a protocol as 'dynamic' when both phonon coordinate and momenta initial conditions are sampled from Eq. ( 13) and the phonon coordinates evolve along with the electronic system according to Eq. (2).

D. Time-dependent Boltzmann equation
In order to make connections with statistical mechanics approaches based on perturbation theory, we also make comparisons with a time-dependent Boltzmann equation treatment of the problem.Starting from Fermi's golden rule for the first-order rate equation for electron occupation in band n at crystal momentum k, f nk and phonon occupation n qν due to electron-phonon scattering processes: [19,21] and Here g ν mn (k, q) are the electron-phonon matrix elements in the band basis.Since we focus mainly on the short time scale dynamics, we do not include phonon-phonon scattering in our TDBE analysis.

III. TIGHT BINDING MODEL
We treat the electronic structure of hBN in a widely used DFT-based tight binding model with nearest neighbor hopping between the two inequivalent sublattice sites.While, in principle, the electron-phonon coupling matrix elements can be extracted from existing DFPT packages, these quantities are almost always given as absolute values which can be used in the TDBE.As the time evolution of the electronic dynamics in MTEF also requires the complex phase of the coupling matrix elements, we derive a BZ extensive expression for the coupling.We treat the nuclear dependence of the electronic Hamiltonian by modeling the hopping term to be exponentially dependent on interatomic distance: where t 0 is the equilibrium hopping term, b is the electron-phonon coupling factor, and d 0 is the equilibrium distance between sublattice sites i and j.This fitting is common in literature for graphene tight binding models, and b can be related to experimental observables [73][74][75].
For simplicity we restrict our study to two dimensions by only considering phonon branches with no out of plane component.We perform a second order expansion of ionic displacement from equilibrium while restricting the hopping term expansion to first order, ignoring the second order electron-phonon coupling coefficients, or Debye-Waller terms, which is quite often done in literature [67,75].By rewriting the electronic operators in terms of plane waves and the nuclear displacements in terms of phonon complex normal coordinates, we obtain the following reciprocal space Hamiltonian (see SI.2 for details): Here X = (z, P) is the collection of phonon coordinates, αk are the electronic site operators for sites α = {a, b} with onsite energies ∆ α = ±∆ responsible for opening the gap E g = 2|∆|.Nearest neighbor sites are connected by vectors δ.The matrices M (q, ν) describe the coupling of the phonon normal coordinates to the electronic system, resulting in the scattering of electronic planewaves from k → k ± q, and depend on coupling terms g δ ν (q), which themselves depend on e αν (q).The required inputs for this model are the lattice constant a 0 , onsite energies ∆ α , hopping term t 0 , phonon dispersion ω qν , polarization e ν (q) and electron-phonon coupling factor b. In this work, we obtained the lattice constant and phonon information from Quantum Espresso cell relaxation and phonon dispersion calculations, and fit the hopping term, onsite energies and coupling factor from a symmetric two band approximation to the conduction band calculated using an uncorrected LDA xc functional.See the Computational Methods section for further details.
The present tight-binding model clearly involves a major simplification, as it treats the electronic system at the independent particle level.However, this simplified treatment is not entirely unreasonable as one can eliminate excitonic effects in experimental setups by placing the monolayer onto a conductive substrate, which screens the electric field and allows the study of free carriers [76,77].Furthermore, the purpose of the TB model is mainly in developing comparisons with ab initio simulations, which go beyond this limitation.Nevertheless, it would be interesting to go beyond this limit in future applications.

Implementing the MTEF Method with the TB Model
For each initial condition X 0 that is sampled from ρ ph,W (X), one can find the single particle orbitals for ĤW (X 0 ), and construct the density operator, where f ( 0 l , T ) is the Fermi occupation at temperature T evaluated for the orbital energy l (X 0 ) = 0 l at this initial configuration.In the above expression for the electronic density we have suppressed the dependence of the orbital on the phonon coordinates.One may then propagate a set of trajectories, associated with the set of initial conditions, according to the MTEF equations of motion (Eq.( 2)), which in the case of the present TB model are more conveniently expressed as follows, We have validated the accuracy of the MTEF treatment of this form of tight binding model by comparing with the nonequilibrium Green's function (NEGF) results of Nery and Mauri [36] for the phonon renormalization of the electronic spectral function in graphene, and we find that the spectral function generated from MTEF simulations is in good agreement with the NEGF results.See the SI.4 for details.

A. Excited Carrier Relaxation in hBN
Due to the lack of inversion symmetry, hBN is a strong insulator with a gap around the K/K points that is depicted in panel (b) of Fig. 1.
This also leads to the potential for selective excitation in the BZ around the K/K points upon irradiation with circularly polarized light.
For simplicity, we did not include direct coupling between the laser field and the phonon modes, and as such the laser field is only indirectly coupled to the phonons via the modification of the electronic occupations.To couple the electronic system to the laser field in the longwave approximation we use a Peierls substitution to modify the electronic wavevector k: We use a circularly polarized laser pulse A(t) using a cos 2 envelope with polarization defined by: We utilize a pump pulse duration of ten optical cycles at carrier frequency ω = 5 eV, such that T pump ≈ 8.3 fs (FWHM ≈ 4.1 fs), with amplitude A 0 = 5 a.u.corresponding to a peak intensity of 7.88 × 10 11 W/cm 2 .The handedness of the laser is controlled by s = {1, −1}, creating left and right circularly polarized light, and we distinguish the K points by which sign of s excites them as K + /K − .
After exposure to the pump pulse we track the occupation of the conduction band states, f ck = |c k | 2 , by taking the expectation value of the projector: Strictly speaking this projector is only valid when the laser field is turned off, otherwise one must project onto the Houston states |ck(t) [78,79].For simplicity we project onto the equilibrium band states, and indicate when the laser field is on (and therefore where this measure is only approximate) using grey shading in the background.In panels (c-f) of Fig. 1 we show a snapshot of the conduction band occupations taken 200 fs after the circularly polarized pump pulse, using the different approaches mentioned above.In all cases the system is initialized at 300 K.In Fig. 1 (c) we initialize the ionic system at the equilibrium geometry with zero velocity.The excited population established after the pulse has a large imbalance between the K − and K + valley carrier distributions and remains completely static throughout the propagation in this case due to a lack of decay channels.
The marginal population around K + , mostly around the K ± M lines, and the dip at K − is due to pumping slightly above the gap, as seen in Fig. 1 (b), though there is no population at the symmetry forbidden K + point itself.
Often in literature, a dynamics simulation will be referred to as Ehrenfest if the ions are allowed to move according to mean field forces, without a unique specification of their initial conditions.We show that performing the simulation with fixed ions or with dynamic ions starting with zero initial velocity, but in either case starting from the equilibrium geometry, results in no qualitative change in the excited electron occupation.By breaking the symmetries of the equilibrium geometry and including static disorder in the phonon system, there is a fundamental difference in the evolution of the excited electronic system.
The results of TDBE dynamics are shown in panel (d) of Fig. 1.The initial carrier distribution rapidly equalizes between the K − /K + valleys, while simultaneously contracting towards the lowest energy states available in the valley bottoms.These dynamics are restricted exclusively to relaxation of the electronic occupations because, as seen in panel (b), the optical phonons initially have no energy available to contribute to scattering the electronic system uphill in energy.
For the case of static disorder, shown in Fig. 1(e), the valley occupations are effectively equalized by 200 fs, however there is no change in occupied energy levels, as seen by the ring around K ± which is present in the occupation immediately after excitation as in panel (c).This can be explained in the framework of band folding: displacing the ionic positions in a supercell is equivalent to folding over the primitive cell BZ bands.This allows the excited charge carriers at K − to have a decay channel into the energetically degenerate folded k points around K + .However, this is a strictly elastic scattering process and therefore is effectively restricted to states inside the initially excited energy window.
The results of allowing for dynamic ion motion at the MTEF level are shown in Fig. 1 (f).Here we see that in addition to homogenization of valley population, there is also some scattering of carriers out of the energetic window in which they were initially excited.One can see that in addition to scattering downwards in energy towards the K ± points, there is also excitation of the electrons upward in the valleys along the K ± M lines.Without exact numerical results it is difficult to address the accuracy of this phenomena, however it could be related to the ZPE leakage of the mean-field dynamics: although the optical phonon occupation is approximately zero initially, over time the ZPE of phonons near the gamma point is drained into the electronic system, incrementally exciting the charge carriers beyond what one sees in the TDBE results at this temperature.However, note that when propagating without pumping the electronic system, there is no loss of phonon ZPE and no heating of the electronic system as the band gap forbids the promotion of electrons due to the order of magnitude smaller phonon energies.See SI.7 for further details.
To assess the relaxation timescales in more detail, we can track the flow of excited charge carriers throughout the BZ by integration of the conduction band occupation within the populated region A K ± around the K ± points: In the valleytronics literature this measure has been used to define the 'valley polarization' or 'valley asymmetry' [4,58,62] and (for TMDs) has been shown by perturbation theory to be strongly influenced by electron-phonon coupling [80].For our purposes we use this term to refer strictly to the population imbalance between the K ± valleys, without reference to the spin resolved bands typically associated.We choose the integration regions to be |k − K ± | < 0.36 Å−1 , roughly corresponding to the populated areas Fig. 1(c).In Fig. 2(a), we see that static displacement, MTEF dynamics and TDBE all capture an extremely rapid population rearrangement within the first 50 fs.In the static displacement method there is a slight inversion of polarization which slowly decays over a 50 fs-2 ps timescale, and the total population outside the valley regions in Fig. 2(b) remains fixed due to carrier energies being confined to their initial excitation energy window.Small changes to the radius of integration do not change the characteristics of these plots due to the normalization with respect to region population in Eq. ( 25), and simply changes the quantitative values in Fig.  25) over the first 50 fs, with an inverse scattering rate γ fit to these data shown as the dotted line.

2(b).
The TDBE and MTEF dynamics clearly display a second time scale, starting from about 50 fs in the case of MTEF and about 100 fs for TDBE, where the valley equalization slows into an asymptotic-like behavior when the population imbalance dips below 0.15.Simultaneously, in the TDBE results there is a down-scattering of carriers which were initially outside of the valley regions as the excited electronic population emits energy into the phonon bath and assumes a more thermal distribution.In contrast, MTEF shows the presence of up-scattering processes as the ZPE from the phonon system is absorbed by the electronic system due to the classical nature of the approximation.See SI.7 for further details.
We gain further insight to the effect of the phonon bath on the ultra-fast excited carrier relaxation or 'valley depolarization' by varying the initial temperature of the system.In Fig. 3 we show the characteristic time scales of depolarization obtained from fitting the valley asymmetry in the first 50 fs to an exponential decay function, f (t) = a 0 +a 1 exp (−t/τ dp ), across multiple temperatures.The most apparent feature is an effective independence of the depolarization rate on temperature until about 300K which can be easily attributed to the high phonon energies in hBN, meaning that the optical phonon branches only begin to have significant occupation at higher temperatures.This behavior has been experimentally observed in the steady-state photoluminescence polarization of MoS 2 , which is directly proportional to the valley lifetime [56], as well as in valley asymmetry decay times in WSe 2 measured by time-dependent Kerr rotation measurements [31].
The depolarization timescales τ dp can be related to the average phonon occupation via a linear function of the scattering rate: γ = 1/τ dp where γ = γ 0 + α n ph .Using the fit τ dp values, we further fit the scattering rates by taking at each temperature the expected oc-cupation value of the Bose-Einstein distribution at the average optical phonon energy, 0.16 eV, showing very good agreement with the data.The fit scattering parameters are α Dynamic = 0.072 fs −1 , α Static = 0.081 fs −1 , and α TDBE = 0.065 fs −1 .Although the ZPE leakage present in the MTEF dynamics can affect electronic populations at long time scales, these results show that for short time scales the MTEF decay rates broadly agree with the TDBE and static displacement approaches; all methods predict this timescale to depend on thermal activation of the phonon optical modes in a consistent manner.
In the context of the MTEF and static disorder simulations, it is worth pointing out that these results converge with an extremely small number of trajectories.For example, while the data shown in Fig. 2 was obtained with N t = 380 trajectories, we can reproduce a graphically converged result with high probability through as little as two random samples.See SI.5 for further analysis.

B. Tracking Valley Asymmetry with Transient Circular Dichroism
Direct experimental observations of time dependent valley populations has been done in TMDs by performing time resolved measurements, such as time and angle resolved photoemission [76,77], time dependent Kerr rotation spectroscopy [31], helicity resolved two-dimensional electronic spectroscopy [81], and TAS [27,82].In this section we focus on TAS using circularly polarized light, Transient Circular Dichroism (TCD).
Pumping the system with circularly polarized light produces an electronic current, j pump (t), that depends on the polarization of the pump.Next, one probes the system with a much weaker circularly polarized probe pulse with a time delay between it's envelope center and the pump center of τ and a duration of T pr to generate the pump-probe current j pump-probe (t, τ ), which encodes the excitation of the system induced by the pump at this delay, and depends on the polarization of the probe.Taking the difference between these currents allows one to calculate the intermediary transient optical conductivity (TOC) [83]: Here i, j are the cartesian directions of the pump and probe, and we have inserted the mask function W but the same probe, σ no-pump (ω) we obtain the true TOC: As with the pump pulse, we also use a cos 2 envelope for the probe, with a strength of A 0 = 0.01 a.u. at the same frequency as the pump for a single optical cycle, T pr ≈ 0.8 fs.The carrier envelope phase (CEP) between the pump and probe for each delay is fixed to zero.See SI.3 for the explicit formulas for the electronic current operator in the TB model and section VI E for the MTEF simulation protocol.
The results for the real component of the TOC calculated in the TB model are shown in Fig. 4, with a delay spacing of ∆τ = 0.5fs and a propagation time of T f = 30fs.The pump chirality in all cases is s pump = −1, the same used in the results of Fig. 1, and the chirality of the probe is s probe = −1 in the left column and s probe = +1 in the right column.The periodic fluctuation of the signal is due to the CEP locking, and can in principle be removed if desired by averaging the results over several CEPs [84,85].Starting with the static equilibrium geometry results on the top row, in panel (a), we see that after the pump, at around τ = 4fs, the signal is saturated around the pumping frequency, whereas in panel (b) following the same pump, but with a probe pulse of the opposite chirality, there is a minimal response indicating a very small population slightly above 5eV.By directly comparing to Fig. 1(c) we see that the signal in Fig. 4(b) corresponds to the small population around (and energetically above) K + induced by the pump, while the saturated signal in Fig. 4(a) corresponds clearly to the large population in the K − valley.
The TOC signal is of course unchanging after the pulse due to the lack of decay channels, again corresponding to the equilibrium results shown in Fig. 2. Turning to the MTEF results, we see a sharp attenuation of the signal in panel (c) following the pump, as well as a broadening of the range of the signal with respect to panel (a) due to scattering within and outside the K − valley.Concurrent with the attenuation of the K − valley signal around 5 eV there is a buildup of signal in panel (d) indicating a buildup of carrier density in the K + region.
We can generalize these predictions of the TOC to the calculated conduction band populations without reliance on defining valley integration areas by taking the difference of the columns shown in in Fig. 4(d) and (c), and integrating over the energy axis: where σ same indicates the TOC obtained by pumping and probing with the same chirality of light as in the left columns of Fig. 4, and σ cross refers to pumping and probing with opposite circular polarizations, as in the right columns of Fig. 4. In Fig. 5 we perform a TCD calculation using frozen phonons and compare the measure in Eq. ( 29) coming from the TCD directly to the valley polarization calculated via Eq.( 25) using the occupation data from the pump only part of the TAS calculation.It is clear that there is good agreement between these two measures.
Finally we perform a TCD calculation using the TDDFT program Octopus [86] 29).The results of the same calculation done through TDDFT are also included.The grey shaded region is the time period in which the laser field is active.The characteristic decay times fit to the TAS data via an exponential starting from time 0 are τ dp = 32 ± 2fs for the tight binding calculation and τ dp = 13.5±0.8fsfor the TDDFT calculation.
The exponential fit to the TDDFT data is also plotted, with the grey line being a guide to the eye.
a delay spacing ∆τ of 5 fs.Due to the size of this calculation, 1800 ions in periodic boundary conditions, it was computationally infeasible to do the full calculation with dynamical ions.However we can compare to a full MTEF calculation sampling the Γ point optical phonons with ten dynamical calculations in a primitive cell in gray.The fully ab initio calculations show a smaller degree of valley polarization at the peak of the pulse, and the supercell calculations show an extremely rapid depolarization, on a timescale roughly two to three times faster than the tight binding model.As expected, the displacements in the primitive cell calculation, while demonstrating some fluctuation of the signal, do not depolarize due to being restricted to Γ point phonons, incapable of scattering electronic states from K − to K + .

C. Tracking Valley Asymmetry with Optical Harmonic Polarimetry
While testing the predictions of the above TCD calculation experimentally is in principle possible -as the generation of circularly polarized light in the range of the experimentally measured hBN gap of 6eV can be achieved via high harmonic generation (HHG) [87] -a much easier measure of the valley asymmetry has been proposed by Jiménez-Galán and colleagues [62,88] and recently experimentally tested by Mitra et.al [63] which is based off the ellipticity of harmonics generated by a linearly polarized, high-intensity, off-resonant probe.In this section we recreate this experiment in silico and test the robustness of this measure when including phonon induced scattering in the simulation.
The basic idea was explained succinctly in the SI of [62]: when driving the system with an off-resonant probe aligned in the Γ − M direction of the BZ, the current response induced parallel to the driving field (j ) will not be affected by the population of the K + /K − valleys, however the current response in the perpendicular direction (j ⊥ ) must be.The contribution to j ⊥ is proportional to the anomalous Hall conductivity arising from conduction band population in regions of non-zero Berry curvature.For equal K + and K − valley populations the anomalous Hall current will have exactly counteracting contributions from these two regions.For unequal population, a non-zero j ⊥ will emerge, which has been used as an experimental measure for valley polarization in TMDs for nearly a decade [89].
Furthermore, since the Berry curvature around the two valleys have opposite sign, the anomalous current arising from a population imbalance around either valley will always be completely out of phase (π) with respect to one another, and both π/2 out of phase with respect to j .This means that the radiation emitted as a result of this current will have an elliptical polarization, with ellipticity ∈ [−1, 1] corresponding to the K − /K + population imbalance, and gives an all-optical measure of the valley asymmetry.Mitra et.al [63] use this phenomenon as a measure of the degree of valley asymmetry in hBN following a bichromatic 'trefoil' light pulse -which has been proposed as a tunable driver of valley selective excitation in monolayer hBN and graphene, as well as bulk TMDs and twisted bilayer graphene [90,91] -finding a small but apparently valley selective signal about 100 fs after excitation.
While the ellipticity of emitted harmonics has been proposed to detect other system properties such as topologically insulating phases, the utility of this measure has been called into question by ab initio simulations due in part to the many technical difficulties related to the generation and interpretation of HHG signals, even in theoretically ideal conditions [51,[92][93][94][95]. Furthermore the HHG spectrum has been found to be very sensitive to Γ point phonon distortions [50].Given these open questions in the literature, along with the results in Sections IV A and IV B indicating that phonons should induce ultrafast sub-30 fs homogenization of any valley asymmetry in monolayer hBN, in this section we test the robustness of this measure of valley asymmetry upon inclusion of phonon degrees of freedom and the time dependence of the asymmetric valley population signal following a trefoil pump in our TB and TDDFT approaches.

Robustness of Harmonic Ellipticity under Phonon-Induced Valley Equilibration
We define the intensity of light emitted due to the nonequilibrium current in a given Cartesian direction x, y, as: where m(t, x) is a mask function.For clean high harmonic generation, we utilize a probe pulse which is far from resonance with the energy gap in hBN, ω probe << E g , allowing us to drive the system at very high intensities, and utilize a 'super-sine' envelope for both the probe pulse and the Fourier transform [96]: where w = 0.75, τ is the center of the mask, and m(t, τ ) is defined to be zero when |t − τ | ≥ T probe /2.This mask is useful as it begins and ends exactly at 0 while having a short and smooth ramp time, maximizing the number of optical cycles present at full intensity.We drive the system along the mirror axis, parallel to the B-N bond, and define this to be the y direction corresponding to pumping along the Γ − M line in the BZ.The probe is depolyed at multiple delays τ relative to the center of a pump envelope: Where I probe is the probe intensity, c is the speed of light, and φ is the CEP, which is always fixed from the start time, i.e., φ = ω probe (τ − T probe /2).The ellipticity of emitted light calculated by Eq. ( 30) at a given energy ω can be determined via the Stokes parameters ([97] Eq.SI.2): where h is the helicity of the signal corresponding to the handedness: -1 for right, 1 for left, 0 for linear, and φ x,y is the phase of the signal.The ellipticity ranges continuously from [−1, 1] defining fully right circularly polarized light to fully left circularly polarized light.One can define the ellipticity of a given harmonic by taking the normalized harmonic yield for each harmonic n, weighted by the ellipticity: where the integral goes from the energy at (n−1/2)ω probe to (n+1/2)ω probe .With these definitions in hand we can test the ellipticity of the harmonics generated in hBN after irradiation by the on-resonant pump from sections IV A and IV B which, as already demonstrated, produces very strong valley asymmetry.
In Fig. 6, we see the HHG spectrum for a fixedequilibrium geometry calculation, 25 fs after being pumped with left circulalry, right circularly or linearly polarized light at ω pump = 5 eV and T pump ≈ 8.3 fs, generated by a probe with ω probe = 1 eV, and T probe = 30 fs.As one would expect, the fundamental matches the polarization of the driving probe field, while there is flipping of the ellipticity at each subsequent harmonic, until reaching the band edge at 4.43 eV, whereupon the conduction band electron response dominates the signal.These results hold for both the TB and the TDDFT spectra, though there are naturally some quantitative differences, in particular concerning the intensity of the emitted light.Under linearly polarized pumping, the ellipticity of the spectrum is effectively flat for the TB model, with a small signal at the 3 rd harmonic, however the TDDFT results show an elliptical response at both the fundamental and 3 rd harmonic, which nonetheless is washed out in the yield calculation.Although the intensity weighted ellip- ticity can be quite large, due to the normalization in Eq. ( 34) the calculated yield does not necessarily reflect this.
Although the 3 rd harmonic signal is small, we find that it provides the cleanest time resolved information when scanning over probe delays, shown in Fig. 7.Here we see precisely the same behavior as in the previous results.When the ions are fixed at the equilibrium geometry, there is a fluctuation in the signal for both the TB and TDDFT results, but it remains non-zero due to a lack of valley asymmetry decay channels.In contrast when including phonon dynamics either via MTEF or static disorder there is a suppression of the initial signal followed by a rapid decay corresponding to the equilibration of valley population seen in Fig. 2.

Robustness of Trefoil Valley Selectivity under Phonon-Induced Valley Homogenization
By combining two counter-rotating circularly polarized fields with a fundamental frequency ω tr and it's second harmonic 2ω tr with a relative amplitude of 2:1, and a delay t d of the 2ω tr field relative to the fundamental, one obtains a laser with a triangular 'trefoil' shape in the plane which can be arbitrarily rotated by tuning t d .The orientation of the trefoil was shown numerically via the semi-conductor Bloch equation to preferentially excite electrons into the K + or K − valleys in monolayer hBN [62,63].The explicit form of the trefoil gauge field we use is: where again m(t) is an envelope function.We chose ω tr = 0.6 eV for a duration T pump = 30 fs as in the experimental paper, with a super-sine envelope and an intensity of 1.67 × 10 12 W/cm 2 .The effect of driving the TB system with this laser is shown via the conduction band occupations plotted in Fig 8.
In Fig. 8, panels (a-c) show the results for the static equilibrium geometry with the trefoil pulse oriented at −30 • , 0 • , and 30 • rotations relative to the x axis in the BZ.The excitation at K + for the 30 • rotation in Fig. 8(c) is quite strong and agrees qualitatively with Fig. 2e of Ref. [62], which shows a similar pumping frequency.At the opposite tuning in panel (a) when K − should be more populated, there is indeed some excitation, and looking closely one can see that it's texture also resembles the 'grape-cluster' structure of the K + valley in panel (c), although the magnitude of excitation is lower.Qualitatively this reproduces the relative difference in excitation density reported in Fig. 3d and 3e in [62], confirming that our model also captures this phenomenon.Halfway between these two tunings, in panel (b), one still sees a strong preferential excitement in the K + valley over the K − valley.
Incorporating phonon dynamics via MTEF in panels (d), (e), and (f), these trends broadly remain 50 fs after pumping, though the fine structure of the excitation seen away from the K valleys is washed out by electronphonon scattering, causing a radially symmetric distribution away from the BZ boundaries that decreases towards Γ.Inclusion of the phonon system via static disorder produces BZ occupations virtually identical to the MTEF results at this timescale.
In all cases, looking within a small region around K + /K − when pumping at −30 • versus 30 • one sees there is indeed an asymmetry in the excitation at the extrema of trefoil orientation.To see how this manifests in the HHG signal, we recreate the experimental setup by utilizing a linearly polarized probe of the same intensity and fundamental frequency as the trefoil pump (ω probe = ω tr ) to create an HHG spectrum.As the HHG from the pump alone has no contribution to the 3ω tr harmonic channel it serves as a natural signal region to study with the probe.In Fig. 9 we compare the ellipticity of the 3 rd harmonic to the valley asymmetry calculated directly from the conduction band occupation induced exclusively by the trefoil pump for fixed equilibrium geometry in panels (a) and (b), and MTEF/static displacement in panels (c) and (d).Because of the diffuse excitation, we choose integration regions for Eq. ( 25) in the BZ formed from vertices at the high symmetry points bisecting the lines The triangular regions inscribed by the dotted lines are formed by vertices at the high symmetry points between the K + /K − points, and are used to define the integration regions for the valley asymmetry calculation in Fig. 9.
between adjacent K + /K − valleys, which neatly encapsulates the excitation around K + for a 30 • trefoil pump.The equilibrium geometry results in panel (b) show that while the valley asymmetry integrated in the conduction band generally shows a tuning of the valley occupations, beginning preferentially in the K − valley at −30 • and in the K + valley at 30 • , it is not exact, as at 0 • there should in principle be more equal population, yet there is clearly a preferential excitation.This is partially reflected in the ellipticity in panel (a), where there is a strong polarization signal at 30 • and 0 • , however, the small degree of K − population under −30 • pumping fails to contribute to a significant negative elliptical emission.
These signals become cleaner when looking at the MTEF/static results in panel (c).There is a much smaller degree of polarization at 0 • , while −30 • is negative throughout.In all pump cases, there is a decay in the signal over the 70 fs depicted which is commensurate with the decay of valley asymmetry seen in panel (d).
Comparing the strength of the signal when using oppositely tuned pumps, these simulations indicate that the relative signal difference is very weak after 100 fs.

V. CONCLUSION
We have derived a general expression for the Wigner transformation of the phonon thermal density matrix applicable to the phonon dispersion for real materials taken from ab initio calculations, and utilized it to perform MTEF calculations in both a reciprocal space TB model and a real space TDDFT supercell approach.This methodology allowed us to simulate the relaxation of asymmetrically excited charge carriers in hBN, to track this through the transient absorption spectra using circularly polarized probes, and to recreate an experimental observation of harmonic ellipticity as a measure of valley population imbalance.The inclusion of phonon degrees of freedom, either dynamically (including anharmonic effects) or statically, produced fundamentally different spectroscopic simulation results due to the complete exclusion of a critical decay channel in the common equilibrium geometry framework.
We discussed connections between the MTEF method and static displacement approaches, finding that MTEF is analogous to an extension of the reciprocal space William-Lax / Zacharias-Giustino coordinate distribution to include nuclear velocities/momenta.We have also demonstrated that MTEF and static displacement methods can capture the phonon driven sub-30 fs valley depolarization on time scales commensurate with the TDBE results.Further, we showed that while static displacement methods are restricted to elastic scattering, MTEF suffers from unrealistic electron heating at long time scales due to ZPE leakage.Despite this, we showed that the temperature dependence of the short time scale relaxation of excited charge carriers agrees well across these simulation methods, exhibiting a plateau in intervalley relaxation lifetimes at low temperature that agrees relatively well with experimental measurements in similar systems [31,56].Furthermore we found that these results converge with a very small number of samples.This work extends MTEF into the domain of the ab initio treatment of periodic systems, providing a starting point for the hierarchy of semiclassical dynamics approaches which build on from the mean field limit, correcting for some of it's most serious shortcomings.However, even with issues in the long time-scale behavior of MTEF, this simulation method provides a unique tool to study the ultrafast response behavior of systems with strong electron-phonon coupling in laser driven regimes far from thermal equilibrium.Inclusion of phonon dynamics has the potential to be utilized to incorporate phonon dynamics into fully ab initio simulations of parametric driving, Floquet engineering and driven phase transitions, while already inclusion of static disorder can substantially change predictions of spectroscopic measurements via inclusion of elastic electron-phonon scattering under arbitrarily complex pump-probe configura-tions.
The rapidity with which our results converge for the results presented here is highly encouraging for other far from equilibrium observables, and given the simplicity of incorporating our method into existing real time simulation protocols of non-equilibrium driven phenomena in quantum materials, as well as the fundamental qualitative differences in the resulting temperature dependent dynamics resolved in a systematic framework, we think that this will be a valuable tool going forward.

VI. COMPUTATIONAL DETAILS A. Tight Binding Model
The lattice parameter of a 0 = 4.734 Bohr was obtained from cell relaxation in Quantum Espresso.The phonon frequencies and displacements were calculated on a Monckhorst-Pack (MP) grid of 30 × 30 using Quantum Espresso, with a 64 × 64 MP electronic k point grid.We utilize a 30 2 k and q grid in the BZ for our Equilibrium, MTEF and static TB simulations in sections IV A and IV B which is sufficiently dense to converge our results compared to a 36 2 grid.In section IV C we used a 60 × 60 k-grid.
The tight binding parameters t 0 and ∆ α = ±|∆| were fit from the conduction band calculated over a 30 × 30 MP grid in the TDDFT code Octopus according to the analytical dispersion relation of the equilibrium geometry tight binding model: An LDA functional was used with Hartwigsen-Goedeker-Hutter (HGH) norm-conserving pseudopotentials with the simulation box having a real space grid spacing of 0.35 Bohr, and being periodic only in two dimensions with a vacuum of 1a 0 on either side of the monolayer.The band gap of E g = 2|∆| corresponds to the uncorrected LDA direct gap of 4.43 eV, and t 0 = 2.68 eV.Although hBN has rather high energy phonons compared to many other materials, the electronic gap remains an order of magnitude larger, meaning that none of the physics of the electron dynamics within the upper band presented here is affected by not correcting this gap.The hopping parameter b was calculated by fitting the bands to Eq. ( 36) for a range of ionic configurations with maximum displacement of 3% of the lattice constant, about 0.14 Bohr from equilibrium in steps of 0.005 Bohr, using a first order expansion of Eq.( 17).The resulting value of b = 2.87 is quite close to the value suggested in [74] (b = 3.3), where it was approximated from Slater-Koster parameters between nearest and next nearest neighbors, and is also in a similar range to the value of 3.37 which has been estimated for graphene [75].
The tight binding model is written in Python and C++/CUDA, and is available at gitlab.comunder kevin.lively_mpsd/graphene-tight-binding

B. MTEF Simulations
We integrated the MTEF equations of motion with an fourth order Runge-Kutta algorithm for the electronic portion simultaneously with a velocity-Verlet type scheme for the phonon configuration, using a time step of 0.1 a.u..

C. TDDFT Simulations
The TDDFT simulation was done with a time step of 5 × 10 −3 fs using a 16 th order Lanczos expansion of the exponential propagator in an approximated enforced time reversal symmetry framework.Only the Γ point was included in the 30 × 30 supercell BZ, i.e. equivalent to a 30 × 30 k grid in the primitive cell BZ under the BvK boundary conditions.We used the same grid spacing, pseudopotentials, and simulation box geometry that were used to fit the tight binding model.

D. TDBE Simulations
All of our TDBE results are generated using input from the tight-binding model Eq. ( 18) in the main text, with band energies nk = ±| l (X = 0)|, and electron-phonon matrix elements taken via projection of the coupling terms onto the band states g ν mn (k, q) = ψ nk |l qν M (qν)|ψ mk .
The delta functions used in the TDBE scattering rate equations are approximated by Gaussian distributions δ(x) ≈ 1 η √ π e −( x η ) 2 with η = 10 meV, and the resultant scattering rates converged with respect to number of k/q points and η [24,98].The equations of motion themselves are integrated using a fourth order Runge-Kutta algorithm with f nk (t), n qν (t) inserted into each derivative and a timestep of 1fs.The initial excited electronic carrier population is set via interpolation of the carrier occupations from the equilibrium geometry using a 72 2 k MP grid just after the cessation of the pulse onto a denser 90 2 k grid.Because we are working with highly non-thermal electronic distributions in a system with extremely strong electron-phonon coupling on top of this gaussian approximation, we find there is a slight drift in the total excited population, and that the ratio of population drift to total population decreases as the initial total excited population decreases.Therefore with the exception of Fig. 1(d), which is included on the same scale for illustrative purposes, the initial TDBE electronic population distribution is set as 1/100 th of the MTEF occupation at the end of the pump, i.e. f TDBE ck (T pump ) = 0.01f MTEF ck (T pump ).

E. Transient Optical Conductivity
The simulation protocol for calculating the transient optical conductivity is as follows: 1. Sample a phonon system initial condition and initialize the electronic system according to Eq. ( 20), giving a total system initial condition IC 0 .
2. Expose the system to a probe of a given chirality s probe from IC 0 , and propagate for a duration of T f .
3. Reset to IC 0 and propagate the system under the influence of the pump to the maximum delay time desired τ max +T f .At each time τ that one wants to have data for, save the instantaneous state of the system IC τ .
4. Load each IC τ and propagate under the same pump, but with an added probe of chirality s probe centered at τ + T probe /2.
7. Repeat steps (1-6) for N t samples and average the results

F. Graphical Representation of Data
The data in figures 1, 4 and 8 were interpolated with the bicubic interpolation function of matplotlib.
Supplemental Information: Revealing Ultrafast Phonon Mediated Inter-Valley Scattering through Transient Absorption and High Harmonic Generation Spectroscopies

SI.1. DETAILS OF THE PHONON NORMAL MODE DESCRIPTION
We follow primarily the convention from Feliciano Giustino's Rev. Mod.Phys.[67], restating many definitions here to make the text self-contained.These definitions are used in the derivation of the phonon wigner distribution and the electron-phonon coupling matrix elements.

A. The Born-Von Kármán Boundary Conditions
The crystalline primitive unit cell is defined by the primitive lattice vectors a i , for i = {1, . . ., l} for the number of periodic dimensions l and the p th unit cell is specified by R p = i n i a i with integers n i ∈ [0, N i − 1].The BvK supercell contains N p = i N i primitive unit cells.The primitive vectors of the reciprocal lattice are denoted by b j , fulfilling the duality condition a i • b j = 2πδ ij .Consider Bloch wave vectors q i defined on a uniform grid in the first Brillouin zone, q = i (m i /N i )b i with integers m i ∈ [0, N i − 1].From these definitions we have the following sum rules: (SI.1)

B. Normal Mode Coordinates
Within the BvK boundary conditions atoms are identified by their positions w/rt the primitive cell, R 0 αp = R p +R 0 α , where p = 1, . . ., N c identifies the primitive cell, α indicates the specific atom within the primitive cell, N c is the number of atoms within each primitive unit cell, and R 0 α denotes the equilibrium position within the primitive unit cell defined by being a minimum energy configuration for a given lattice configuration.We can further identify small displacements from these positions via δR αp = R αp − R 0 αp .For small displacements from the minimum energy configuration, we can write the potential as Since the IFC must be invariant under any operation which maps between supercells, it obeys certain symmetry operations.We can encode this via the Fourier transform of the IFC, defined as the dynamical matrix [99] where M α is the mass of the α th ion.The dynamical matrix is hermitian and positive definite allowing real eigenvalues, denoted as ω 2 qν α D α,α (q)e α ν (q) = ω 2 qν e αν (q).(SI.5) with corresponding unit vectors δ0 p .We take the Fourier transform by replacing the real space lattice site operators with their planewave counterparts: Note that we make the distinction between the primitive cell lattice sites R p and sublattice sites R 0 αp belonging to the α ∈ {A, B} sublattices: R 0 αp .Although we have ionic displacements δR αp , these are fundamentally defined w/rt to the sublattice sites R 0 αp which define the periodicity of the crystal, and therefore how the Fourier transform is defined.
Inserting equations (SI.13) and (4) into equation (SI.11), we have the following: Where µ δ = {0, −a 2 , −a 2 − a 1 } are the primitive lattice vectors connecting primitive cell R p to p = {p − 1, p, p + 1}, and in the second line we have summed through p, applying the BvK boundary conditions Eq. (SI.2), and defined the electron-phonon coupling term as e aν (q) .(SI.15) We can group g δ ν (q) and everything in parentheses together and call it the electron-phonon coupling operator M (q, ν).By summing q through C, we obtain the actually implemented form of the tight binding Hamiltonian in Eq. ( 18): ĤW (X) = H ph (X) + Ĥe + ν,q∈A x qν M (q, ν) + ν,q∈B x qν M (q, ν) + M † (q, ν) + iy qν M (q, ν) − M † (q, ν) , (SI.16) where H e gathers the bare hopping and onsite energy terms.This means that the equations of motion for each x i qν , y i qν trajectory are explicitly (SI.17) Looking at the structure of M (q, ν) one can see that M (q, ν) + M † (q, ν) will be hermitian while M (q, ν) − M † (q, ν) will be anti-hermitian.Therefore expectation values of the former will always be real and imaginary for the latter.

SI.3. ELECTRONIC CURRENT OPERATOR
The current density operator can be defined as for the super cell surface area Ω and the electron position operator Where i is summed over primitive cells and the number operator ni+δ is understood to be a b sublattice number operator to reduce notational clutter.We can utilize the definition of the sublattice anticommutation, where α, β can be either â or b site operators to derive the identity: Taking the commutator and utilizing the identity (SI.21) we obtain Expanding in plane waves and phonons as before we have: We take Ω to be N p times primitive cell area Ω 0 = 19.5bohr 2

SI.4. ELECTRONIC STATE RENORMALIZATION IN GRAPHENE
To provide a sanity check of our derivation of g δ ν (q) in the tight binding model, we first apply our tight binding Hamiltonian in Eq. ( 18) to graphene and compare with recent results from Nery and Mauri [36] using static displacement averaging to calculate the electron-phonon interaction renormalized spectral function A k (ω).
Reference [36] utilizes an electron-phonon coupling parameter η, which in our model corresponds to η = bt 0 /d 0 .We reproduced their Figure 8 results using the same parameter set of ∆ α = 0, t 0 = 2.58eV and d 0 = 1.413Å for two different coupling strengths, η expt = 4.42 eV/ Å(a value extracted from experiment), and η inter = 2.5η expt by running MTEF dynamics on a 36 × 36 k and q grid for 12fs.We calculate the retarded Green's function for band n via: By taking the Fourier transform of the mean field propagated Green's function, using the mask function W (x) = 1 − 3x 2 + 2x 3 , we obtain G W,nk (ω).The spectral function is of course defined to be A The results shown in Fig. SI.1, have good agreement with the Nery and Mauri results, namely a significant broadening of the signal with increasing temperature and coupling strength, as well as a decrease in signal intensity away from K and Γ for both cases.The small negativities in the spectral function for high energy k points near Γ in the lower panel are due to the extremely anharmonic dynamics arising from the artificially strong electron-phonon coupling constant.Such spectral negativity is a known feature of MTEF calculations with anharmonic forces [41].
We further compare our model to a specific peak at k = 0.75K in Fig. SI.2.The various sampling approaches that are compared are explained in section SI.6.In panels (a) and (b) we see that the sampling approaches used in the text, dynamic and static, agree very well with data extracted from [36], with the primary difference being slightly broader tails.With this check we are confident that simply replacing the parameters of our tight binding model with those fit to hBN is a reasonable approach to the electron-phonon coupling in this system.

SI.5. CONVERGENCE OF VALLEY HOMOGENIZATION
For this analysis we use the Normalized Root Mean Squared Displacement (NRMSD) to quantify how much a given time dependent signal varies from another: We take f (t) = VA Nt (t) calculated with the largest number of trajectories available, in the case of the data in Fig. 2, N t = 380.We let g(t) = VA N t (t) be the same signal calculated with a random selection of N t << N t trajectories.The NRMSD value for this collection of N t trajectories tells us how much this signal deviates from the more converged value calculated with N t trajectories.By pulling different random collections of size N t from the data we can histogram the resulting NRMSD values, giving us a probability distribution of what errors we can expect with this number of trajectories.The results of this analysis are plotted in Fig. SI.3.Plotted in panel (a) in black is the same MTEF dynamic phonon data plotted in Fig. 2(b), alongside an example of a signal with an NRMSD of 0.02 in red.It's clear that this signal already appears qualitatively converged, yet we see from the error probability distributions in panel (b) that this a high error outlier of a result, even when using only two trajectories!Instead, with quite reasonable numbers of samples, the results rapidly converge towards signals which are graphically indistinguishable from results using over an order of magnitude more trajectories.

SI.6. ALTERNATIVE SAMPLING APPROACHES
In the special displacement method (SDM) of Zacharias and Giustino [34,35], the convergence of sampling the configuration distribution in Eq. ( 14) can be accelerated for small supercell sizes (i.e.sparse q grids in the BZ) by taking a handful of specific configurations.These involve taking positions of magnitude |x qν |, |y qν | = σ qν and carefully A summary of the various choices of initial conditions for the phonon reduced coordinates and momenta.Although the electronic system properties can be obtained via time evolution for all these methods, the phonon system is dynamic only when rqν , sqν = 0.The symbol "→ nqν,T " indicates convergence with increasing Nt.conditions and dynamics choices is shown in Table SI.6.The principle question is whether these different sampling techniques lead to altered observables, or faster convergence.We start by looking at a phonon renormalized property of the electronic system in a thermal equilibrium state.Taking the graphene spectral function in Fig. SI.2(a), we see that using 'dynamic' sampling, that is drawing phonon momenta and coordinates in a straight-forward manner from Eq. 13 and n T -dynamic already lead to slightly different spectral function intensities.In Fig. SI.2(c), we track the convergence of this signal by taking the average error as a function of time.Given that we are Monte Carlo sampling an initial distribution this error is defined as the standard deviation of the observable calculated with N t trajectories, over √ N t .Panel (c) shows that the signal produced by n T sampling converges marginally faster than that produced by dynamic sampling, but that at large configurations, there is virtually no difference in the two approaches.
We see something similar when comparing the 'static' and 'ZG-static' sampling approaches in , with the primary difference of this comparison being that the spectral function magnitudes for these two methods agree nearly exactly.The outlier however appears to be n T -static, which has a significantly different spectral intensity than any static or dynamic approach.This indicates that the majority of the effect responsible for renormalization of the electronic spectral function at this k point can be accounted for via displacement within the first standard deviation of the phonon distribution.Of course in the harmonic limit, setting phonon momentum and displacement absolute values to ±|σ qν | will result in motion in phase space entirely on the circumference of the phase space circle of radius √ 2|σ qν |. obtained by full Monte Carlo sampling, or even Monte Carlo sampling on the ± √ 2|σ qν | radial points.
However, this is a value calculated at thermodynamic equilibrium.To see the effect for far from equilibrium properties, we turn to the temperature dependence of the valley homogenization timescale in the TB hBN model, seen in Fig. SI.4.The n T -dynamic and ZG-static sampling approaches sit almost precisely on top of their coresponding dynamic and static results.Quite notably the n T -static results again constitute an outlier to the other methods, with a scattering fit parameter of α n T −static = 0.129fs −1 , compared to α Static = 0.081fs −1 , α ZG−static = 0.082fs −1 , α Dynamic = 0.072fs −1 , α n T −Dynamic = 0.068fs −1 and α TDBE = 0.065fs −1 .Clearly this sampling method, while producing the correct thermal occupation value through strong static displacement, appears to overestimate the effect of phonon renormalization, when compared to Monte Carlo sampling of the phonon coordinate / momenta distribution.
Finally we can investigate whether there is any significant advantage in the dynamical case to sampling with the n T -dynamic approach.We repeat the NRMSD probability analysis of section SI.5 for the time dependent valley asymmetry in an MTEF calculation and report the results in Fig. SI.5.We find that the probability of obtaining a signal using a small number of trajectories, which has a small deviation from the signal one would obtain using a much larger number, looks effectively identical between the two sampling approaches, thus in this case this sampling method confers no appreciable advantage.

SI.7. ZPE LEAKAGE
We track the effects of ZPE leakage when propagating with dynamic MTEF by looking at the phonon occupation number n qν calculated via Eq.(SI.26).Given sufficient samples, n qν will converge towards the expected thermal occupation n qν,T , which for the optical phonons at 300K means approximately 0. Therefore as the ZPE drains out of a given phonon mode, the occupation will go towards −1/2.
In Fig. SI.6 we show a selection of the occupation numbers of the highest energy optical branch phonon modes over time.In panel (a) when the system is exposed to a laser pulse, we see that on timescales commensurate with the long time scale excitation outside the valleys in Fig. 2(b), there is a loss of energy in the phonon modes.The initial occupation value of 0 requires a large number of trajectories to converge to exactly, but as discussed throughout the text, the dynamical electronic observables of interest converge rapidly.Panel (b) shows the phonon occupation when the electronic system is not pumped.In this case, due to the large electronic gap, there is nowhere for the phonon energy to go, and instead one just sees oscillation of the occupation numbers over time.This oscillation is due to anharmonic forces arising from exposure to the electronic system, and can in principle be analyzed to capture the renormalized phonon frequencies.
In Fig. SI.6(c) we take the phonon occupations from the last time step of panel (a) at 2ps and plot them against the distance of their q vector from Γ.There is very clearly a direct correlation between the how close a phonon mode is to Γ and the amount of energy it loses to the electronic system.This may be related to the fact that the phonon modes are coupled directly to the electronic system via a linear dependence through the nearest neighbor hopping term.Therefore the modes closer to Γ which correspond to a coherent reduction in the nearest neighbor distance throughout the system, most strongly excite the electrons and allow a conduit for vibrational energy to go into the electronic system.

FIG. 1 .
FIG. 1. Panel (a) shows a schematic of circularly polarized light exciting a sheet of hBN.The upper half of panel (b) shows the electron bands calculated using DFT with an LDA xc functional, alongside the tight binding (TB) bands fit to match the band onset.The Fermi level has been centered between the conduction and valence bands.The arrow indicates the approximate energy of excitation around a specific K valley, and the shading indicates the occupation of states in the TB model just after the laser pulse, interpolated from the MTEF results.The lower half of panel (b) shows the phonon dispersion calculated using DFPT, with shading indicating the initial phonon occupation number at 300 K, with the translational modes at Γ set to zero.The bottom four panels show the occupation of the conduction band 200 fs after the laser pulse under different approximations to the electron-phonon dynamics: (c) initially equilibrium geometry with zero velocity, (d) TDBE, (e) static displacement, and (f) MTEF dynamics.The phonon system in panels (d-f) is sampled at/set to 300 K.The labeling of the otherwise degenerate K + /K − points is based off the polarization of the pump and the scale is capped at 0.8, to emphasize small differences.
FIG. 2. (a):Valley population asymmetry at T = 300 K defined by Equation(25).Panel (b) shows the fraction of conduction band population outside the K ± regions.The grey region corresponds to when the pump is turned on and the dotted line at zero is a guide to the eye.

FIG. 3 .
FIG.3.The timescale of an exponential fit of the valley depolarization rates calculated with Eq. (25) over the first 50 fs, with an inverse scattering rate γ fit to these data shown as the dotted line.
FIG. 4. The real part of the TOC: Re [σ(ω, τ )] calculated in the TB model.The arrows indicate the chirality of the pump and probe, either pumping and probing with the the same handedness or opposite handedness.Panels (a) and (b) are the equilibrium geometry results and the MTEF results at 300 K are shown in panels (c) and (d).The left most panels show the spectral weight of the pump, |E(ω)| which is active during the grey highlighted time span.The logarithmic scale has been set to emphasize the signal in the cross valley data, while minimizing the fluctuations arising from CEP locking.

50 FIG. 6 .
FIG. 6.The intensity weighted ellipticity and corresponding elliptical yield of the harmonics generated in a fixedequilibrium geometry calculation 25 fs after excitation with left circularly, right circularly or linearly polarized light using the same pump parameters as in sections IV A and IV B. The vertical dashed line indicates the conduction band edge at 4.43 eV.(a): the TB model results, (b): TDDFT results, (c): the elliptical yield of panel (a) calculated with Eq. (34), (d): the elliptical yield of panel (b).

FIG. 8 .
FIG. 8.The occupation of the conduction band 50 fs after irradiation with the trefoil pump in the TB model.The path of the trefoil pump in the plane is arbitrarily rescaled and drawn within the BZ as the blue line.Panels (a), (b) and (c) show the results for the static equilibrium geometry when the degree of rotation of the trefoil with respect to the BZ x axis is at −30 • , 0 • and 30 • respectively.Panels (d), (e) and (f) show the MTEF results at the same rotations -Note the smaller colormap scale.The triangular regions inscribed by the dotted lines are formed by vertices at the high symmetry points between the K + /K − points, and are used to define the integration regions for the valley asymmetry calculation in Fig.9.

FIG. 9 .
FIG.9.The left column shows the ellipticity yield of the 3 rd harmonic generated after irradiation with the three trefoil pulse rotations seen in Fig.8, compared to the right column showing valley asymmetry calculated by applying Eq. (25) to the triangular dotted regions in Fig.8.Panels (a) and (b) show the static equilibrium geometry results while panels (c) and (d) show the MTEF results in bold lines and static results in dashed lines.Note the order of magnitude difference between the y-axes scales for the two approaches.
FIG. SI.1.The spectral functionA k (ω) of graphene calculated in the tight binding model through MTEF, with electron-phonon coupling constants η and temperatures matching the simulations from Ref.[36], Figure8.See text for details.

FIG. SI. 4 .
FIG. SI.4.The temperature dependence of the short time scale valley population homogenization, as seen in Fig.2, calculated with all of the sampling techniques in Table SI.6.

2 FIG. SI. 5 .
FIG. SI.5.The probability distribution of getting a particular NRMSD of the valley polarization seen in Fig. 2 when selecting an arbitrary Nt number samples taken via n T -dynamic sampling versus straight forward Monte Carlo sampling.

FIG. SI. 6 .
FIG. SI.6.The phonon occupation number in a selection of the highest energy optical modes over time, when the electronic system is either (a) exposed to a pump or (b) allowed to propagate without a pump.Panel (c) shows the phonon occupation numbers from panel (a) at the final time, plotted against the distance of their q vector from Γ.
. 7. The elliptical yield of the 3 rd harmonic calculated with Eq. (34), following excitation with the same circularly polarized resonant pump from sections IV A and IV B. The grey region corresponds to the duration of the pump, while the vertical dotted line indicates the final time of pump-probe overlap. FIG 2)p δR αp δR α p , (SI.2)where we have defined the so called Interatomic Force Constant matrix (IFC) as C, with C αp,α p ∈ R d×d .Treating the nuclear positions as operators, we define canonical momenta P αp by the canonical commutation relation [R αp , P αp ] = i δ α,α δ p,p and write the Hamiltonian operator for the nuclei in real space as