Sub-cycle multidimensional spectroscopy of strongly correlated materials

Strongly correlated solids are extremely complex and fascinating quantum systems, where new states continue to emerge, especially when interaction with light triggers interplay between them. In this interplay, sub-laser-cycle electron response is particularly attractive as a tool for ultrafast manipulation of matter at PHz scale. Here we introduce a new type of non-linear multidimensional spectroscopy, which allows us to unravel the sub-cycle dynamics of strongly correlated systems interacting with few-cycle infrared pulses and the complex interplay between different correlated states evolving on the sub-femtosecond time-scale. We demonstrate that single particle sub-cycle electronic response is extremely sensitive to correlated many-body dynamics and provides direct access to many body response functions. For the two-dimensional Hubbard model under the influence of ultra-short, intense electric field transients, we demonstrate that our approach can resolve pathways of charge and energy flow between localized and delocalized many-body states on the sub-cycle time scale and follow the creation of a highly correlated state surviving after the end of the laser pulse. Our findings open a way towards a regime of imaging and manipulating strongly correlated materials at optical rates, beyond the multi-cycle approach employed in Floquet engineering, with the sub-cycle response being a key tool for accessing many body phenomena.

laser-cycle-averaged modifications of material properties.Yet, in strongly correlated systems the sub-laser-cycle time scale is also highly relevant: a typical electron-electron interaction parameter U ∼ 1eV corresponds to the time-scale ∆t ∼ 1/U ∼ 1 fs, well below the cycle of a standard infrared driver, with the respective dynamics potentially leading to such remarkable features as e.g. a transition from Coulomb repulsion to effective electron-electron attraction induced by half-cycle pulses [35,36].
One way to probe and control the sub-laser-cycle electronic response is to use few-cycle pulses with controlled carrier-envelope phase (CEP) [17,19,20,37].In solids, these pulses have been used to detect photoemission delays [38] and quantify the time-scale of non-linear response to light [39], image surface states in topological insulators [40], resolve and control highly non-linear electronic response in bulk dielectrics, 2D materials, and nano-structures [8,[41][42][43][44].Yet, the physical picture of electron-electron correlations evolving on the sub-cycle scale in strongly correlated systems remains elusive.
Here we introduce a sub-cycle multidimensional spectroscopy of electron dynamics in solids and apply it to a strongly correlated system.Our approach uses the CEP-dependence of the correlated multi-electron response to decode the complex interplay between different many-body states, triggered by the interaction with a few-cycle mid-IR control pulse.Ultimately, the analysis of the multi-dimensional spectra allows us to uncover the physical picture of the underlying correlated dynamics in this system, both in space and time.Fundamentally, we show that the subcycle one-particle response is able to track many-body dynamics by providing direct access to many-body response functions.The possibility to spectroscopically analyze the underlying excitation pathways, as introduced below, is key to understanding non-thermal materials control.
We consider a half-filled Hubbard model on the two-dimensional square lattice for fermions, supporting a realistic two-dimensional band dispersion with the characteristic van Hove singularity and sharp band edges.The lattice is driven by a strong field linearly-polarized along the lattice diagonal, triggering a fully two-dimensional response (in contrast to previously employed Bethe-or hypercube-lattices [45] or one-dimensional chains [12,46].)To treat the non-perturbative time-dependent problem, we employ the non-equilibrium extension [47,48] of the dynamical meanfield theory (DMFT) [49].The algorithm and its realization are described in Ref. [50] (see Methods).The method was benchmarked against exact diagonalization of one-dimensional finite chain [12], see Supplementary Information (SI).The implementation is based on the NESSi simulation package for non-equilibrium Green's functions [51].We employ a single-band Hubbard model with parameters for undoped La 2 CuO 4 (LCO): the lattice constant a 0 = 3.78 Å, the nearest-neighbour hopping T 1 = 0.43 eV [52], the Hubbard U = 2.5 eV [53].We use few-cycle pulses with central wavelength of λ = 1500 nm (with ω = 0.827 eV) and duration of 7.7 fs (full width at half-maximum, FWHM), with a total simulation time 32.8 fs.To demonstrate that our results are typical for the low-frequency regime ω U , we also present simulations for λ = 3000 nm (ω = 0.413 eV).
Figures 1(a,c,e) show the temporal profile of the occupied density of states (see Methods) for the field strengths F 0 from 0.1 to 2.0 V/A.The voltage across a unit cell approaches the hopping rate, a 0 F 0 ∼ T 1 , at F 0 ∼ 0.1 V/A (I 0 ∼ 1.6 × 10 11 W/cm 2 ).Thus, F 0 ∼ 0.1 V/A could modify the effective hopping rate within the laser cycle and alter the structure of the correlated system.
Indeed, the transfer of spectral weight from the quasi-particle peak (QP), located near zero energy to the Hubbard bands becomes prominent as soon as F 0 approaches 0.1 V/A, see Fig. 1(a): after the transition at ∼ 17.5 fs, the spectral density remains predominantly in the lower Hubbard band and does not return back to the QP after the pulse.Already for this field, Fig. 1(a) hints at the importance of the sub-cycle response: the cycle-averaged renormalized hopping T 1 → T 1 J 0 (F 0 a 0 /ω) = 0.97T 1 does not lead to any substantial changes in the spectral density, let alone to the major restructuring observed in Fig. 1(a) (see SI).
At higher fields (Fig. 1(c,e)), we see substantial transfer of the spectral density to the upper Hubbard band (situated at E = 1.25 eV), with the electron density peaked at the energies corresponding to the upper and lower Hubbard bands.Crucially, this dichotomic structure survives well after the end of the pulse.Figures 1(a,c,e) thus signify transition from a metallic to a highly correlated state in which the light-driven current is fully quenched (see Fig. 4(d) below.) To understand these complex many-body dynamics, we first look at the cuts (Figure 1(b,d,f)) of the electron density for specific energies corresponding to the lower Hubbard band (LHB, E = −1.25 eV), upper Hubbard band (UHB, E = 1.25 eV) and QP (maximizing at −0.215 eV).The exchange of population in Fig. 1(d) has three distinct regimes, marked as three shaded areas: around 4-11 fsec (red), 11-16 fsec (green), and beyond 16 fsec (blue).The first regime (red shading) shows decreasing electron density at the energy corresponding to the QP peak and increasing density at the energies corresponding to LHB and U HB (see also Figure 1(b)), with the populations at LHB and U HB energies oscillating out of phase (see also Figure 1(b)).In the second regime (green shaded area in Figure 1(d)), the density at U HB energy raises, while the density at LHB energy decreases.The third regime is most surprising: we observe in-phase oscillations of the electron density at U HB and LHB energies (blue shaded area in Figure 1(d)).At higher intensities the electron densities at U HB and LHB energies become locked: both populations are equal and oscillate exactly in phase (blue shaded area at t > 14 fs in Figure 1(f)).The maxima of these locked populations are synchronized with the minima in the density located at the QP .
The locking of populations at three key energies of the system is correlated with the onset of metal to insulator transition observed in Figs.1(c,e).Cartoon in Fig. 2(a) illustrates the three key field-free many-body states of our system using the language of DMFT [54].Characteristic many-body states contributing to the signal at LHB mainly involve electrons localized on singly occupied lattice sites, QP represents superposition of delocalized and localized electrons, while doubly occupied and unoccupied sites are the characteristic features of many-body states contributing to U HB. Analysis (see SI) of Fig. 1 (f) suggests that the rate of flow of electron density from LHB and U HB bands maximizes near zeroes of field oscillation (F (t) 0), while the rate of flow of electron density to LHB and U HB bands maximizes near instantaneous maxima of the field (|F (t)| F 0 ).The oscillations of the density at QP are in phase with the laser field: the minima coincide with F (t) 0, the maxima coincide with |F (t)| F 0 .FIG. 2. Spectroscopic nature of the Keldysh Green function G < .(a) Cartoon view of the key many-body states corresponding to the spectroscopic signal at the energies of LHB, QP , U HB. Orange circles stand for lattice sites, green circle symbolises electron localized on a lattice site, white arrows show orientation of electron spin, gray circle and a grey cloud symbolize delocolized electron, two oppositely oriented white arrows symbolize a doubly occupied site.Double headed arrows indicate possible subcycle transitions driven by the field in the phase locking regime (see Fig. 1 (f)).(b) Scanning delay τ between the pump and probe pulses (violet), and delay t between the pump-probe pair and the control pulse (red), yields two dimensional data set for the Green function G < (t, t−τ ) emulating photionization signal.Scanning the carrier-envelope phase (φCEP ) of the control pulse yields the third dimension of the spectroscopic signal G < (t, t − τ, tCEP ).(c) Fourier transform of G < (t, t − τ, tCEP ) with respect to all arguments yields G < (Ωτ , Ωt, ΩCEP ).Cartoon view of |G < (Ωτ , Ωt, ΩCEP )| for Ωτ fixed at the energy corresponding to LHB.Red-brown peaks illustrate the Floquet ladder associated representing laser-dressed LHB state.The appearance of green peaks at QP ± nω and U HB is due to non-adiabatic transitions between the laser dressed states LHB ↔ QP and U HB ↔ QP .The extension of green peaks in the ΩCEP dimension quantifies the sub-cycle response time.
To decode the underlying physics and establish the role of sub-cycle electron dynamics we need a spectroscopy sensitive to such dynamics, which we introduce below.To this end, we exploit the full spectroscopic nature of the one particle Keldysh Green's function G < (t, t − τ ), which can be retrieved from multi-pulse time and angular-resolved photoemission (trARPES) experiments (see e.g [55] and SI.) Figures 1(b-d) show the Fourier image of such signal with respect to τ , yielding the population of the occupied states (with energy Ω τ ) prior to photo-ionization.The intense low-frequency field plays the role of a control pulse, which modifies our system between these two events.Thus, formally, we have a sequence of three pulses (Fig. 2(b)) reminiscent, but not identical, to the set-up of non-linear 2D spectroscopy [56].Fourier transforming G < (t, t − τ ) with respect to τ , selecting a value of Ω τ , and Fourier transforming G < (t, Ω τ ) with respect to the delay t between the control pulse and the pump-probe pair, we obtain the spectrum of the states (transition frequencies) that the Floquet state associated with Ω τ couples to (see Methods).Indeed, for the dynamics described by a single Floquet state, G < (t, Ω τ ) should behave periodically as a function of t, and its Fourier transform will only show the sidebands at ±nω (n is integer, see the ladder of red-brown peaks in Fig. 2(c)).In contrast, in the presence of non-adiabatic transitions between the Floquet states, G < (t, Ω τ ) becomes aperiodic and its Fourier transform will show the new frequencies appearing due to non-adiabatic excitations (see green peaks in Fig. 2(c).) Next, we can reveal the underlying sub-cycle dynamics by scanning the CEP φ CEP of the IR (control) pulse and recording the resulting response G < (t, t − τ, t CEP ) as a function of t CEP = φ CEP /ω.The Fourier transform of G < (t, t − τ, t CEP ) with respect to t CEP shows the speed of response [57]: from instantaneous to cycle averaged.The broader the resulting spectrum is with respect to the CEP (see schematic in Fig. 2(c)), the stronger is the non-trivial CEP dependence, the stronger, faster, and more sensitive to the instantaneous electric field are the nonadiabatic transitions between the Floquet states.We now focus on non-adiabatic dynamics to highlight the part of the interaction, which is fundamentally different from the cycle-averaged response.To this end, we consider the difference between the derivatives of the Green's function with respect to t and t CEP , ∆G ) contains only non-adiabatic transitions and provides direct access to non-equilibrium two-particle dynamics via the respective Green's functions K pm ijσ (t, t − τ ) and Kpm ijσ (t, t − τ ) (Methods): where f (t) is the envelope of the short pulse, is an amplitude of non-adiabatic transition between time-dependent states evolving from the field-free eigenstates |Ψ m (t 0 ) , |Ψ p (t 0 ) , ∆G (1)< ij are the non-adiabatic terms of one particle nature, and Here U is the on-site Coulomb interaction, c † iσ (c jσ ) are the fermionic creation (annihilation) operators for site i (j) and spin σ, n iσ = c † iσ c iσ is the particle number operator.To visualize the non-adiabatic transitions encoded in ∆G < (Ω τ , Ω t , Ω CEP ) we fix Ω τ and plot the resulting 2D spectrum |∆G < (Ω τ , Ω t , Ω CEP )| as a function of Ω t (vertical axis) and Ω CEP (horizontal axis), as schematically shown in Fig. 2(c).
At the intermediate field F 0 = 0.5 V/A, (Fig 3(c,d)) the peaks at Ω CEP = 0 in Fig. 3(c) (for Ω τ = −1.25 eV, LHB) correspond to the Floquet ladders associated with both QP and LHB, with broad overlapping steps.The CEP-sensitive dynamics becomes more significant in the LHB-QP transitions.The dominant part of the upper peak at Ω CEP = 2ω is due to the QP lifted by one photon or U HB.
The spectroscopic portrait of the system in the high-field regime (F 0 = 2 V/A), Fig 3(e,f) contains several new features.First, we see broad spectrum along the vertical axis (Ω t ) for a wide range of horizontal axis (Ω CEP ).Second, the sub-cycle dynamics is very important and the spectrum has individual sub-cycle cut-offs: the highest positive and highest negative Ω t depend on Ω CEP .These cut-offs appear to be proportional to the instantaneous values of the laser field and increase with increasing Ω CEP > 0. Third, the overall 2D spectra are shifted towards positive values of Ω t for LHB, corresponding to absorption, and negative values of Ω t for UHB, corresponding to emission.c,d) correspond to the lower laser frequency ω = 0.413 eV, demonstrating that our results are not sensitive to the frequency of the non-resonant field.
This spectroscopic information points to the following simplified picture of many-body dynamics in the locking regime.The synchronized cycling of charge flow along the two "circuits," LHB → QP → LHB and U HB → QP → U HB, appears to dominate the dynamical photon-dressed many-body state just before it freezes into the final state after the end of the laser pulse.At the peak of the field, a localized electron in LHB is promoted to the QP, leaving an empty site behind.At the same time, the strong field also destroys doublons: a doublon from UHB loses one of its electrons into the QP.The total energy does not change during these two synchronized processes.Near the instantaneous zero of the field, the QP electrons localize on the empty or half-filled sites of the lattice completing the oscillation cycle: localization at an empty site contributes to the population in LHB, while localization at a half-filled site contributes to the population in UHB.However, in contrast to light-driven dynamics in dielectrics [11], the flow of charge between LHB and QP or UHB and QP cannot be understood as single electron oscillation.Indeed, the time-evolution of double occupancy (Fig. 4(b)) shows that in the locking regime the doublon production saturates at 0.25, the maximal value for an uncorrelated system.Thus, every produced doublon decays and its fragments randomly occupy lattice sites.This "randomness" limits the total energy increase to zero energy and precludes population inversion, i.e. achieving higher electron density at U HB compared to LHB.This apparent randomness could be a sign of entanglement destroyed by our observation.For example, calculating double occupancy we trace out a multitude of different ways in which doublons are crated and destroyed each laser cycle, i.e. we do not follow their individual Feynman paths when calculating this observable and therefore we partially destroy the entanglment in our system.
In correlated systems, large increase in electron temperature can transform a metal into a bad metallic or insulatinglike state.However, in our case, the opening of the gap and the peculiar dynamics observed already within a small fraction of the laser cycle is clearly not thermal.In contrast to phonon-driven transitions [58], our mechanism is purely electronic.
Our results show the power of sub-cycle response to record and address non-equilibrium many-body dynamics.Here we considered the first derivative with respect to ∆ ≡ ∂ ∂t − ∂ ∂t CEP , which gave us access to non-equilibrium two-particle correlations.Likewise, the n-th derivative, ∆ n , contains the non-equilibrium n-particle Green's functions.If the correlations are strong, as is the case here, odd and even multiphoton pathways contributing to the CEP-dependent response can involve photon absorption by different electrons, with correlated interaction establishing coherence between these interfering events.Note that in contrast to the standard equilibrium (diagonal in m) expression for the two-particle Green's functions K mm ijσ (t, t − τ ) and K mm ijσ (t, t − τ ), ∆G < ij features non-equilibrium two-particle correlations, which would not be recorded in long pulses, as ∆G < ij = 0 in this case, corresponding to the standard Floquet regime.
Our findings demonstrate the possibility of manipulating phases of correlated systems with strong non-resonant fields on the sub-cycle time-scale, in a manner that is robust against the frequency of the driving field.While in atoms or molecules such non-resonant light-induced modifications of the electron density vanish once the light is turned off, in a strongly correlated system light-induced changes can lead to persistent modifications surviving after the end of the pulse, as we see here.The pulse can thus transfer the system into a correlated state inaccessible under equilibrium conditions, with non-adiabatic transitions triggering non-equilibrium many-body correlations, see Eq.( 2).The achieved final state is controlled by charge density and currents shaped on the sub-laser-cycle time-scale; the spectroscopy introduced here can provide key insights in analysing and designing such excitation pathways.

A. Simulations
The Hamiltonian is where i, j label the lattice sites, U is the on-site Coulomb interaction, c † iσ (c jσ ) are the fermionic creation (annihilation) operators for site i (j) and spin σ, n iσ = c † iσ c iσ is the particle number operator.The hopping amplitudes T ij (t) between the sites i and j include the nearest-neighbor (T 1 ) and the next-neighbor (T 2 ) terms.The external low-frequency laser field (frequency ω < U, W = 8T 1 ) is included via the Peierls substitution, where A(t) is the field vector-potential, F(t) = −∂A(t)/∂t.The one-particle dispersion is The total energy E tot (t) = E kin (t) + E pot (t) includes potential and kinetic terms, where G< k (t, t) is the gauge-invariant [59] lesser Green function.The momentum distribution function is The population density is calculated as Due to the limitation in time data, the selected Fourier transform produces some blur on the graph of occupied states for the first 5 fs.The time-resolved photoemission intensity is given by where S(t − t p ) is the envelope of the probe pulse centered at t p [60].

B. Direct access to non-adiabatic many-body dynamics
Here we develop the strategy for isolating the non-adiabatic response in the one particle Green's function.In a long laser pulse, the dependence of the instantaneous electric field on the CEP amounts to the overall time shift t = t + t CEP .This trivial dependence is of no interest and should be removed when analysing the Green's function G < ij (t, t − τ, t CEP ).To see how this should be done, let us assume for the moment that its CEP dependence amounts only to the overall time shift of the argument t = t + t CEP : If this is the case, the Green's function should obey the following equation: Therefore, differentiating the Green's function G < ij (t, t − τ, t CEP ) with respect to t and t CEP and subtracting the resulting terms we obtain the differential contribution ∆G < ij (t, t − τ, t CEP ), which no longer contains the trivial dependence: We use ∆G < ij (t, t − τ, t CEP ) for building the 2D spectroscopy maps as a function of Ω t and Ω CEP .Now we can explicitly evaluate ∆G < ij (t, t − τ, t CEP ) (Eq.14 ) for an arbitrary Hamiltonian H(t).
Since we consider coherent dynamics, we can rewrite Eq. ( 16) for the Green's function as follows (see e.g.[61], [62]): where |Ψ m (t 0 ) are field free eigenstates of the system, t 0 is the initial moment before the laser pulse.Inserting the resolution of identity I = n |Ψ n (t 0 ) Ψ n (t 0 )| on the field free eigenstates and switching from Heisenberg to Schrödinger picture by transforming the temporal dependence from the operators to the wave-functions we obtain: where |Ψ m (t) are time-dependent basis states evolving from the field free states under the influence of the full propagator: |Ψ m (t) = T e −i t 0 H(t )dt |Ψ m (t 0 ) , and T is time-ordering operator.Since the evolution is unitary the time-dependent basis states remain orthogonal to each other at any time t within the pulse.
Since our Hamiltonian explicitly depends on t and t CEP , H(f (t), t + t CEP ), where f (t) is a short pulse envelop, the derivatives of the Hamiltonian with respect to each of this times can be explicitly calculated: where we have introduced an operator ∆ ≡ ∂ ∂t − ∂ ∂t CEP .To evaluate ∆G < ij (t, t − τ ) we take into account that one can rewrite the Schrödinger equation for ∆Ψ m in the equivalent form connecting ∆H and ∆Ψ m explicitly (|Ψ m (t 0 ) =0): Writing ∆G < ij explicitly we obtain: Substituting Eq. 20 into above equations and limiting ourselves to the terms of the leading order with respect to ∂ ∂t and ∂ ∂f , we find that ∆G < ij is proportional to the amplitudes of non-adiabatic transitions a mp (t) = Ψ m (t)| ∂ ∂f |Ψ p (t) and between the quasienergy states: where H H (t − τ ) is the Hamilton operator in the Heisenberg picture.Substituting the explicit expressions for the the commutators c † i (t − τ ), H H (t − τ ) and c j (t), H H (t) for the Hubbard model and focussing on the correlated part of the Hamiltonian we obtain that ∆G < ij (∆G < ijσ ) plotted in Figure 3 directly reflects the non-equilibrium two-body Green's functions +∆G is an amplitude of non-adiabatic transitions between time-dependent states evolving from the field-free eigenstates |Ψ m (t 0 ) , |Ψ p (t 0 ) , ∆G (1)< ij are the non-adiabatic terms of one particle nature and Here U is the on-site Coulomb interaction, c † iσ (c jσ ) are the fermionic creation (annihilation) operators for site i (j) and spin σ. n iσ = c † iσ c iσ is the particle number operator.Note that in contrast with the standard equilibrium (diagonal in m) expressions: which would appear in ∂G ijσ /∂t, off-diagonal K pm ijσ (t, t − τ ) and Kpm ijσ (t, t − τ ) in Eqs.(29,30) feature non-equilibrium two-particle correlations, which would not be recorded in the long pulses (∆G < ij = 0), corresponding to the standard regime of cycle averaged field driven dynamics, which lays at the foundations of the Floquet engineering.
C. Recovering full Green's function from photo-electron measurements.
While time and angular-resolved photoemission (trARPES) experiments are directly related to the Green's function, going back from trARPES to the Green's function is nontrivial.In particular, analyzing the photoemission from a single pulse is restricted by energy-time uncertainty [60].The multi-pulse spectroscopy does not suffer from this limitation, and allows for full retrieval of the Green function, see e.g [55] and discussion below.
The full information in the Green's function can be retrieved by suitable measurements, e.g., exploiting the dependence of the photoemission signal on the phase delay between interfering parts of the photoemission pulse [55].With this in mind, we can say that G < (t, t − τ ) emulates the photoionization signal arising from the interference of two photo-ionization events at t and t − τ .
To demonstrate how a time-resolved photo-emission experiment may in principle resolve the full Green's function, we start from the general expression given in Ref. [60], where orbital and momentum indices are omitted form simplicity, ω is the frequency of a probe pulse, and s(t) its envelope.It is easy to see that a single Gaussian probe pulse of width δt implies a measurement of G < (ω, t) with an uncertainty-limited filter in time and frequency.However, with suitable pulses, Eq. (33) shows that in principle the full time dependence can be retrieved from experiment.For example, to measure G < (t, t ) in a given time window, we choose an orthonormal basis φ n (t) for time-dependent functions in that interval, and expand −iG < (t, t ) = n,n φ * n (t)g n,n φ n (t ).The matrix g n,n is hermitian and positive definite.A probe pulse S(t) = φ n (t) then measures the diagonal components, I < = g n,n .A probe pulse S(t) = φ n (t) + e iϕ φ m (t) gives I < = g n,n + g m,m + e −iϕ g n,m + e iϕ g n,m , so that off-diagonal components g n,m can be obtained by scanning the phase difference ϕ.
D. Spectrocopic nature of the double-time lesser Green's function The double-time lesser Green's function provides information about the spectrum of occupied states of the system and is indispensable for visualizing electronic structure and dynamics.We review the emergence of laser-dressed states in its structure using the approach presented in section "Direct access to non-adiabatic many body dynamics" starting from Eq. (17).Consider the typical Floquet regime corresponding to CW pulse.The quasienergy states can be written as: where E m (E n ) is the quasienergy of a Floquet state m (n) and f m (t) (f n (t)) is a periodic function of time.Introducing auxiliary functions Φ nm (t), Φ (+) we can rewrite the expression for Green's function Since the functions f mn (t), f nm (t − τ ) are periodic, we can expand them in Fourier series: Thus, The situation changes dramatically in the presence of non-adiabatic transitions between the Floquet states.Suppose such non-adiabatic transition couples the Floquet state m to another Floquet state m .Then Φ nm (t) = λ mm (t)e i(En−Em)t f nm (t) + λ m m (t)e i(En−E m )t f nm (t), (41) and where the coefficients λ mm (t) represent the amplitudes of non-adiabatic transitions between the quasienergies E m and E m .In this case, the product Φ mn (t − τ )Φ nm (t), contributing to G < ij (t, t − τ ) in Eq.( 37), acquires three additional terms (Eqs.(44)(45)(46)): Non-adiabatic transitions from state with quasienergy E m to state with quasienergy E m lead to new features in the spectrum both along Ω τ and Ω t dimensions.Indeed, while the term Eq.( 43) is similar to Eq. ( 37), the term Eq.( 46) adds new frequency to the spectrum along Ω τ dimension.In addition, terms represented by Eqs.(44,45) oscillate (in time t) at "new" frequencies ±(E m − E m ), which will appear along Ω t direction.Thus, we see that non-adiabatic transitions lead to significant restructuring of the spectrum encoded in G < ij (t, t − τ ) and revealed after Fourier transforms wrt to τ and t.
To reveal the time-scale of non-adiabatic transitions we need to employ the additional dimension, sensitive to sublaser cycle features of electron dynamics.The carrier-envelope phase (CEP) is a natural choice.Scanning CEP we obtain G < ij (t, t − τ, t CEP ).Fourier transform wrt t CEP [57] allows us to tag each non-adiabatic transition and quantify the role of sub-laser cycle dynamics in restructuring the spectrum of the system and in the formation of the final insulating state.

II. SUPPLEMENTARY INFORMATION A. Cycle-averaged Floquet picture
Here we show that the metal-insulator transition we find in our simulations can not be understood in terms of the simple cycle-averaged picture, which underlies the standard mechanism of Floquet engineering.It becomes clear after the analysis of results of laser-driven dynamics at the lowest field strength F 0 = 0.1 V/A, shown in Fig. 5. Comparing the prediction based on the cycle averaged picture show in in Fig. 5(a), which uses renormalized hoppings, with the full simulation shown in Fig. 5(b), we see that panel (a) does not describe our observations.Specifically, the red curve in Fig. 5(a) presents equilibrium density of occupied states calculated with renormalized hoppings T ef f ij = T ij J 0 (AR ij ) and the same U (i.e.U = 2.5 eV).It is essentially identical to the electron density prior to the pulse, because for F 0 = 0.1 V/A the Bessel function J 0 1, meaning that the hoppings are hardly modified.In Fig. 5(b) (identical to Fig. 1a) we see that at about 17.5 fs the peak of the density is shifted towards the energy around -1.25 eV.Remarkably it stays there after the field is off.The cycle averaged picture ( panel (a), red curve) suggests that the peak of the density should be around zero energy at all times.The stark contrast between the density in Fig. 5a and the density on in Fig. 5b after 17.5 fs means that the sub-cycle modification of the electron density is crucial for establishing the state of the system observed after 17.5 fs and after the pulse is off.

B. Temporal dynamics in the locking regime
We have established three different regimes of electron dynamics, described in Fig. 1 (b,d,f).Here we focus on the most interesting regime prominent at F 0 = 2 V/A, which corresponds to pronounced locked oscillations of electron densities at LHB (red) and UHB (blue) energies: after 13 fs the populations at these energies are nearly equal, locked in phase, and oscillate out of phase with the electron density at the QP energy (green).Here we provide additional information to show that these oscillations are well synchronized with the instantaneous electric field.
Fig .6 (a) specifies the time instants at which the instantaneous field is F (t) = 0.It shows that the populations at LHB (red) and UHB (blue) reach maxima just before the instantaneous zero of the field.The instantaneous rate of population decay (time derivative of red (LHB) and blue (UHB) curves) appears to be maximized at F (t) = 0. Thus, the rate of flow of electron density from LHB and U HB bands maximizes near zeroes of field oscillation (F (t) 0).The population at the QP (green) reaches minima at the zeroes of the field, meaning that it is in phase with the oscillations of the laser field.
Fig. 6 (b) specifies time instants at which the instantaneous field reaches its maximal value |F (t)| = F 0 .It shows that the populations in LHB (red) and UHB (blue) reach minima just before the maximum of the field, when the electron density at the energy of QP maximises.The instantaneous rate of the population increase (time derivative of red (LHB) and blue (UHB) curves), appears to be maximized near the maxima of the laser field oscillation (|F (t)| F 0 ).
Thus, the charge density oscillates between QP and LHB, U HB and these oscillations are synchronized with the instantaneous field.Figs. 5 and 6 emphasize the importance of the sub-laser cycle dynamics and motivate the development of the sub-cycle multidimensional spectroscopy described in the main text.This spectroscopy allows one to zoom into the sub-cycle electron dynamics and identify dominant pathways of charge flow responsible for the metal-insulator transition in our system.

C. Benchmark simulations
We have benchmarked our IPT-DMFT on square lattice against the code described in Ref. [12], performing exact diagonalization for the finite 12-site one-dimensional chain.We set the hoppings T = 1 eV, on-site Coulomb repulsion U = 6 eV, pulse vector potential amplitude A 0 = 5, pulse FWHM is 3fs, pulse central frequency ω = 10 eV.In order to compare our two-dimensional lattice model to one-dimensional chain, we choose linear pulse polarization along [10] direction and relatively large field amplitude.
Although the physics is different between 1D and 2D systems due to the existence of closed loops and additional scattering channels in two dimensions, the resulting HHG spectra look qualitatively similar (see Figs. 7 and 8).

FIG. 1 .
FIG. 1. Temporal evolution of density of states showing light-induced transition from metallic to Mott-insulating states.(a) F0=0.1 V/A, (c) 0.5 V/A,(e) 2.0 V/A.Artifact of the Fourier transform which appears at 0-2.5 fs is covered by shadow.(b,d,f) Oscillations of electron density at energies corresponding to LHB (E = −1.25 eV, red), U HB (E = 1.25 eV, blue), and QP (E = −0.215eV, green) for F0 = 0.1 V/A (b), F0 = 0.5 V/A (d) and F0 = 2 V/A (f).Red, green and blue shading marks three different regimes of field-driven dynamics.Red shading: density at LHB and U HB oscillates out of phase; blue shading: locking regime, density at LHB and U HB oscillates in phase; green shading: intermediate regime.Vertical lines indicate the maxima in LHB and U HB populations, correlated to minima in QP populations in locking regime.
The respective peaks at Ω CEP = 4ω (Fig 3(e)) and Ω CEP = −4ω (Fig 3(f)) indicate their strongly sub-cycle nature.Forth, while the direct non-adiabatic transitions between LHB and UHB (orange circles in Fig 3(e)) become visible, the non-adiabatic transitions are still dominated by the couplings LHB ↔ QP and U HB ↔ QP .Indeed, the onset of locking is synchronized with the saturation of energy transfer from the field to the system (see Fig 4(a)).For the highest field (F 0 = 2V /A) the maximum energy saturates at zero.The onset of locking (see Fig 4(c)) is accompanied by the suppression of the current (see Fig 4(d)), which is fully quenched at ∼ 18 fs (see Fig 4(c,e)) when the insulating state is established.Note that Figs 4(

a mn k a nm
k +l e −ilωt e −i(Em−En+k ω)τ(40) and the Fourier transform wrt τ yields Ω τ = E m − E n + k ω, where E m − E n represent spectral energies on the vertical axis of Fig1(a,c,e).If we fix Ω τ = −1.25 eV (LHB), then k = 0. Fourier transforming Eq. 40 wrt t we obtain Ω t = lω and a nm l are the amplitudes of the Floquet ladder starting from zero energy, i.e. the Floquet ladder corresponding to LHB for Ω τ = −1.25 eV.Thus, in the standard Floquet picture, fixing Ω τ = E (as it is done in Fig3) leads to observation of standard Floquet ladder from the state E.

FIG. 6 .
FIG.6.Temporal oscillations of electron density at the key energies of the system for F0 = 2 V/A, shown vs oscillations of the laser electric field.(a) The vertical lines mark time instants, at which the instantaneous field is equal to zero, F (t) = 0. (b) Vertical lines mark time instants, at which the instantaneous field is maximal, |F (t)| = F0.