Probing ultrafast electron correlations in high harmonic generation

We show here that electron-electron ( e - e ) interactions that evolve dynamically on ultrafast timescales can be imprinted onto the high harmonic generation (HHG) spectrum driven by intense laser ﬁelds, even far from any resonant multielectronic behavior of the medium and over a wide energy range. Speciﬁcally, we ﬁnd that in bi-elliptical HHG, the ellipticities of the bichromatic pumps that maximize/minimize the HHG yield are sensitive to the level of theory describing the process, and considerably shift upon inclusion of e - e interactions. We explore this phenomenon by performing ab initio calculations with varying levels of theory on noble gas (Ar, Kr, and Xe) and molecular (carbon monoxide) systems. Interestingly, the sensitivity to e - e interactions in atomic systems arises due to mean-ﬁeld Coulomb repulsion, while in molecules the strong laser ﬁeld also excites ultrafast dynamical correlations. Our approach can be used to benchmark theories for multielectron systems interacting with strong ﬁelds. For instance, it can be used to determine the validity of various approximations, such as, the adiabatic exchange-correlation approximation in time-dependent density functional theory.

Here we show that HHG in a bichromatic, bielliptical, driving pulse geometry [19,54] can be effectively utilized to probe ultrafast (dynamical) many-body interactions.In this approach, the HHG spectrum is measured vs the driving field's ellipticities and/or relative phases, yielding spectrograms that carry unique spectral signatures of the target medium.We find that the structure of the spectrograms is sensitive to the level of theory used in calculations, and even though corrections due to e-e interactions are generally small, they are qualitative and systematic, i.e., causing generic effects such as shifting of spectrogram peaks and changing their relative intensities.Importantly, these spectral signatures arise over wide energy ranges and in systems far from exhibiting multielectronic resonances.We demonstrate this phenomenon by performing ab initio calculations on several atomic (Ar, Kr, Xe) and molecular (carbon monoxide, CO) systems with varying levels of theory, showing that the spectrogram structure changes with different choices of exchange-correlation (XC) functionals in time-dependent density functional theory (TDDFT), and with the inclusion/exclusion of dynamical e-e interactions.The source of this sensitivity is investigated and found to be electrostatic e-e repulsion in atoms, while in molecules it also arises from laser excitations of ultrafast correlations.pulses are transversely elliptically polarized with an equal ellipticity (ε), and opposite handedness, resulting in the following electric field: where E 0 is the amplitude corresponding to the laser power (I 0 ), is the beam's amplitude ratios, φ is a relative phase, ω is the fundamental frequency corresponding to a wavelength of 800 nm, and A(t ) is an envelope function.A bielliptical spectrogram is generated by measuring the harmonic spectrum emitted from the medium as it interacts with E(t ) in Eq. ( 1) while scanning ε in the driving field from −1 to 1 [see Fig. 1(a) for an illustration].This type of approach was shown to be selective towards chirality [23] and atomic orbital sizes [55].Here we show that it can also be used to probe dynamical e-e interactions that evolve upon interaction with ultrashort driving pulses.
The interaction of the medium with E(t ) is described within the fixed-nuclei and electric-dipole approximations by the following Hamiltonian (atomic units are used throughout): where r j is the coordinate of the jth electron, ∇ 2 j is the Laplacian operator with respect to r j , R I is the position of the Ith nuclei with charge Z I , and we neglect spin-orbit interactions for simplicity.The HHG response can be conceptually found by solving the many-body time-dependent Schrödinger equation defined by H (t ) and Fourier-transforming the induced nonlinear polarization in the medium: where ρ(r, t ) is the electron density.Unfortunately, an exact numerical solution for Eq. ( 2) is currently feasible for systems with only few electrons; hence further approximations must be made.In this article, we utilize Kohn-Sham (KS) TDDFT to solve for ρ(r, t ) directly using the real-space grid-based code OCTOPUS [56][57][58].Here the electron density itself is propagated (instead of the full wave function) under an equation of motion that is formally equivalent to Eq. ( 2), but where the e-e interaction term is split to two terms: (i) a Hartree Coulomb repulsion term and (ii) a quantum-mechanical interaction that is approximated by an XC functional of the electron density (leading to the potential V XC [ρ]).The equations are cast into an auxiliary KS system where they reduce to single-particle equations for the KS orbitals, ϕ KS j (r, t ), governed by the one-body Hamiltonians: where V KS (r, t ) is given by the sum of interactions of each electron with the nuclei, the Hartree term, and the XC potential, respectively: and we have used the adiabatic approximation for the XC functional.Note that since h KS j (t ) itself depends on ρ(r, t ), the KS orbitals are propagated in tandem to construct the timedependent density -ρ(r, t ) = j |ϕ KS j (r, t )| 2 , from which the HHG emission can be recovered.
On the other hand, it is common practice in HHG calculations to simply neglect e-e interactions altogether.In this approximation, the Hamiltonian in Eq. ( 2) reduces to onebody terms: where V eff (r) is an effective potential describing the interaction of electrons with the nuclei and may also contain some information about e-e interactions and screening in the field-free system (i.e., prior to the laser being turned on).V eff can be determined empirically in order to fit some physical property of the ground-state system (e.g., ionization potential), or it may be obtained by other physical considerations.
In this work, we set V eff for a noninteracting electrons system using the latter approach-we use ab initio density-functional (DFT) calculations to find the ground-state electron density of the system, ρ 0 (r) = ρ(r, t = 0), and set V eff as equal to the KS potential that is associated with ρ 0 , such that V eff (r) = V KS (r, t = 0) (see the Appendix for details).Following this process, the individual orbitals ϕ j (r, t ) are propagated separately to yield the HHG response in the noninteracting electrons system.The numerical approach described above allows us to consistently compare the HHG response from an interacting electron system (using TDDFT where interactions evolve dynamically) to a noninteracting one (where e-e interactions are frozen at the ground-state configuration).We emphasize that the initial orbitals and potentials are identical in both approaches and so are their energy eigenvalues corresponding to the ionization potentials; therefore, any deviation in the HHG response between calculations is directly associated with dynamically evolving many-body interactions (and not due to interference of emission from multiple orbitals, as was explored for instance in Ref. [59]).It is worthwhile to point out that one intuitively expects the differential response from both systems to be small, because even in strong fields the electron density is only weakly modulated (ρ(t )-ρ 0 ρ 0 ), meaning that the variation in the dynamically evolving KS potential should be accordingly small.Nonetheless, even small changes in the potential energy may turn out to largely impact the HHG response, as it is an extremely nonlinear process.

III. RESULTS IN ATOMIC MEDIA
Having introduced the theoretical model, we first explore bielliptical HHG from noble gas atoms, where Fig. 1(b) presents the calculated bielliptical HHG spectrogram from atomic Ar (total intensity).It is instructive to discuss the general structure of the spectrogram prior to addressing the role of e-e interactions.As seen in Fig. 1(a), the driving laser field's time-dependent polarization rapidly changes with ε, which is the physical cause for the modulation in the HHG yield in Fig. 1(b).At the edges of the spectrogram, both components of the field are circularly polarized such that E(t ) exhibits a threefold rotational dynamical symmetry [see Fig. 1(a)].In this geometry the emission is determined by the interference of three recollision events per optical cycle [60], and the emission of 3n harmonics (for integer n) is symmetry forbidden [see Fig. 1(b)] [19,61].On the other hand, at the center of the spectrogram E(t ) is cross-linearly polarized and exhibits a reflection dynamical symmetry [62].In this geometry the emission is dominated by two recollision events per optical cycle [63,64], and all harmonic orders are symmetry allowed.In between, the spectral yield oscillates as the dominating HHG trajectories are modified [54,55].
We now analyze the effects of e-e interactions in the HHG bielliptical spectrograms from Ar. Figure 1(b) presents calculations using TDDFT and a corresponding noninteracting electrons model [with Eqs. ( 1)-( 6), see Appendix 1-3 for numerical details].While both models qualitatively capture the same features, clear deviations between them arise over a wide energy range.This is despite the fact that the presented spectra is well below argon's Cooper minima at 50 eV [32,36].In particular, we note that the maximizing/minimizing ellipticities (the position of peaks/valleys in the spectrogram) largely vary with the inclusion of e-e interactions for some harmonics (up to ε ≈ 0.1).For example, the peak of the 13th harmonic at ε = 0.42 in the TDDFT calculation is shifted to ε = 0.34 for the noninteracting electrons calculation.Similarly, the intensity ratios between different peaks are also highly sensitive to the level of theory.We further note that at ε = 0, the harmonic yield must either be at a local maxima or a local minimum, since the spectrogram is symmetric about its center (because the medium is reflection invariant).This constraint can be efficient for separating different levels of theory, because in a few cases, the noninteracting electrons model predicts a minima at ε = 0 while the TDDFT calculations predict a maxima, or vice versa [e.g., see harmonic 25 in Fig. 1(b)].Notably, these discrepancies between the models tend to vanish with increasing laser power (see Appendix 5).This result can be intuitively interpreted as a generalized form of the strong-field approximation (SFA)-under intense enough illumination, both the atomic potential and the e-e interaction terms become negligible compared to the lightmatter interaction term.We further note that this enhanced sensitivity to e-e interactions results from the bielliptical HHG scheme that allows measuring robust experimental features (in the spectrograms) rather than measuring absolute HHG intensities.This is very advantageous, since it allows probing even small effects that arise from correlations.
In order to shed further light on the effects of e-e interactions on bielliptical HHG spectrograms, we utilize the factorization picture [14], where P(t ) decouples to where a ion , a prop , and a rec denote the ionization, propagation, and recombination contributions to the HHG emission from a specific channel, respectively, and the summation is performed over all contributing channels that are calculated separately (e.g., different trajectories, different electrons, etc.).Even though Eq. ( 7) is formally correct only when applying the saddle-point method to solve the response of the noninteracting electrons system, it can still provide an intuitive interpretation for the response that does include e-e interactions (because corrections are small) [53].Within this picture, e-e interactions can be thought of as inducing slight perturbations to the electron trajectories that contribute to the HHG emission (i.e., changing their duration, physical trajectory, etc.), which translates to variations in phases and probabilities of different channels that comprise P(t ).Since P(t ) results from the interference of all channels, phase shifts can lead to shifting of spectrogram peak/valley positions, and changes to channel probabilities lead to spectrogram intensity variations.In particular, we note one mechanism that affects the probability coefficients of different channels, which is changes induced to a ion over time-in the TDDFT calculations, the Ar atoms ionize during their interaction with the laser field such that the overall positive charge on the ion increases between optical cycles.This in turn reduces the ionization probability from cycle to cycle, attenuating the photoelectron yield (see Appendix 5 and Fig. 3) and reregulating the a ion term for subsequent channels.This effect vanishes when neglecting dynamical e-e interactions.The discrepancy in ionization yields between the models suggests a breakdown of the factorization picture, as the contribution from different channels can no longer be fully separated (the ionization of one channel influences the propagation of another).We note that this effect could depend on the level of theory being used (e.g., choice of XC functional in TDDFT).Most importantly, the bielliptical HHG spectrogram probes all three processes in Eq. (7).That is, not only are the ionization and recombination processes probed at a variety of angles and intensities, but so is the propagation of multielectron wave packets in the continuum, including correlations between the ionized fragments and the ionic core.Consequently, the spectrograms are highly suitable for benchmarking theoretical approaches.
We also perform similar calculations on Kr and Xe atoms [Figs.1(c) and 1(d)].We find that each spectrogram contains fingerprints unique to the target medium.That is, even though Ar, Kr, and Xe have relatively similar valence electronic structure, the maximizing/minimizing ellipticities for each atom are slightly different.Additionally, in all of these systems distinct discrepancies between the TDDFT and the noninteracting electron calculations arise, even though we are well below the energies of Cooper minima or the giant multielectron resonance in Xe [39][40][41], further indicating the generality of the technique.We also note that Ref. [55] already contains experimental data from atomic Ar that can be compared to in order to test our results; however, the experiment in this case was too noisy to benchmark the theory.
Next we wish to investigate which term in the KS potential is responsible for the deviation between the TDDFT and noninteracting calculations, i.e., is it the mean-field Coulomb e-e repulsion term, the XC term, or both?Thus we recalculated the bielliptical HHG spectrograms within the TDDFT approach [Eqs.(3)- (5)] but where the XC potential in Eq. ( 5) was kept frozen during propagation (while the Hartree term was evolved in time).The results of this calculation are very close to the full TDDFT calculation for all three tested atomic systems [Ar, Kr, and Xe in Figs.1(b)-1(d)].Hence we conclude that for atomic noble gases, the dynamical evolution of the XC potential is mostly negligible, whereas mean-field e-e repulsion is the dominant source of deviations between the levels of theory.This result indicates that it's possible to model e-e interactions in HHG as perturbations to noninteracting electron models of noble gases (since the deviation due to e-e interaction does not originate from a quantum entanglement).

IV. RESULTS IN CARBON MONOXIDE
We next move to explore molecular systems.Figure 2 presents bielliptical HHG spectra from x-aligned CO molecules showing similar peak/valley structures but somewhat more complex behavior.Specifically, the HHG intensity oscillates more rapidly and nonuniformly compared to the atomic systems (e.g., most harmonics in Fig. 2 show an abrupt cusp at ε = 0).We attribute this to the angularly dependent molecular ionization rates that are effectively scanned as ε is varied (the field's peak intensity changes its spatial orientation with ε with respect to the molecular axis).Remarkably, the inclusion of dynamically evolving e-e interactions leads to deviations in the HHG spectrograms that are even more pronounced than in the atomic case.We verify that the deviations are observed using different XC functionals [Figs.2(a) and 2(b)].Interestingly, Figs.2(a) and 2(b) show that in the case of CO, the dynamical evolution of both the Hartree term and the XC term is significant (unlike in atoms where the Hartree term played the dominant role).This means that in the molecular case, quantum XC evolves dynamically on ultrafast timescales with the laser excitation, reminiscent of recent results in strongly correlated solids [66].Figure 2(c) further presents the calculated HHG spectrograms from aligned CO using different approximations for the XC functional, showing a different structure-these present exemplary theoretical predictions for experiments that can determine which functional is more well behaved, i.e., benchmark the theory (see Appendix 5 for similar results from Ar). Notably, local density approximation (LDA) and Perdew-Burke-Ernzerhof (PBE) XC functionals lead to very different HHG spectra from aligned CO, suggesting that there might be a breakdown of the LDA in strongly driven CO.This should be examined in future work, along with a systematic study of different XC functionals and molecules, and comparisons to available experiments.

V. EXTENSIONS TO MULTIDIMENSIONAL SPECTROSCOPY
Lastly, we discuss the possibility to extend this technique in two main avenues: (i) rather than scanning only ε in E(t ), one could scan other degrees of freedom such as phases, intensities, etc. to generate multidimensional spectrograms; (ii) rather than observing the harmonic yield, one may focus on the HHG emission ellipticity.In the Appendix we present FIG. 2. Bielliptical HHG spectrograms from aligned CO calculated with and without dynamical e-e interactions.(a) Integrated intensity per harmonic order vs ε with a full TDDFT model (using the LDA and a SIC) [65], compared to the corresponding noninteracting electrons model, and the corresponding calculation where the XC potential is frozen.The intensity of each harmonic is normalized to its maximum along the ε axis for clarity.Selected harmonics are shown.(b) Same as (a), but where PBE [67] XC is used with a SIC [65].(c) Comparison between the full TDDFT calculated spectra in (a) and (b).Calculations correspond to a laser power of 10 14 W/cm 2 and using 2 = 1/3, φ = 0 (see Appendix 1-3 for details).
data that addresses both of these extensions: (i) in Appendix 6 we show that the maximizing phase in the cross-linear HHG geometry is indeed sensitive to the level of theory in an analogous fashion to the laser ellipticities (though showing a much reduced sensitivity).This result conceptually proves that extensions to multidimensional spectroscopy might enhance sensitivity to probing correlations.(ii) In Appendix 7 we find that for aligned CO, the ellipticity response of bielliptical HHG is as useful for probing correlations as the intensity-based response, which can be further utilized for ultrafast spectroscopy.For atomic media on the other hand, we observe a peculiar result-the HHG emission ellipticity is largely insensitive to the level of theory (see Fig. 6).In fact, we find that the emitted harmonics' ellipticities are also insensitive to the target atomic species and to whether or not the harmonics are generated perturbatively (below the ionization threshold) or nonperturbatively (above the ionization threshold).This suggests that the ellipticity response in bielliptical HHG from atoms depends mainly on the driving laser parameters and not on the medium properties.The physical origin of this result is still not known (see further discussion in Appendix 7).

VI. CONCLUSIONS
In summary, we explored bielliptical HHG from atomic (Ar, Kr, Xe) and molecular (CO) systems using ab initio TDDFT calculations.We found that in this geometry, the spectral response of the system is affected by dynamically evolving e-e interactions that are imprinted onto the HHG emission.This occurs over a wide energy range and far from any multielectronic resonant behavior.The sensitivity is expressed as qualitative and systematic (easy to measure) modifications to the harmonic yield vs the driving field ellipticity, such as shifting the maximizing/minimizing ellipticities for the yield of a given harmonic order.We analyzed the origin of these effects in terms of Coulombic and quantum contributions, and found that in noble gases the mean-field Coulomb repulsion term dominates the total role of e-e interactions, while in molecules, mean-field effects and ultrafast quantum correlations are both important.Lastly, we demonstrated that the sensitivity to the level of theory can be used to benchmark time-dependent multielectron theories robustly and effectively (by comparing experimental results to predictions), including TDDFT [3], density matrix approaches [7], Green's function approaches [68,69], and wave-function-based methods [43,45,53].Extensions of this approach to multidimensional spectroscopy of electron correlation are also possible.We note that we have neglected here phase-matching and focal averaging effects, corresponding to experimental conditions of thin nonlinear media (e.g., as in Refs.[19,20,23]).Looking forward, our work paves the way for robust methods of testing theoretical approaches in strong-field physics using HHG, e.g., perturbative corrections to standardly used models [14,49], the adiabatic XC functional approximation in TDDFT [3], testing new XC functionals [70,71], or new theories.We also believe that the presented results will be useful for the development of future ultrafast spectroscopy techniques.
gratefully acknowledges the support of the Adams Fellowship Program of the Israel Academy of Sciences and Humanities.The authors thank Eliyahu Bordo, Gil Ilan-Haham, Michael Krüger, and Ido Kaminer from the Technion, Israel, for fruitful discussions.

APPENDIX 1. Ground-State DFT Calculations
All DFT calculations were carried out using the real-space grid-based code OCTOPUS [56][57][58].The equations were discretized on a cylindrical grid of radius 25.2 bohr (xy plane) and length 32.8 bohr (z axis).Atomic species were centered to the grid origin, and the CO molecule was aligned along the x axis, with both atoms symmetrically positioned with respect to the origin.Calculations were performed using two XC functionals: (i) the LDA with an added self-interaction correction (SIC) [65], implemented in an optimized effective potential (OEP) method (within the Krieger-Li-Iafrate (KLI) approximation [72]); (ii) using a PBE [67] XC functional with the same SIC [65].The frozen core approximation was used for inner orbitals, which were treated with appropriate norm-conserving pseudopotentials [73].The [He] inner shells of Ar, C, and O, the [Ar]3d 10 inner shell of Kr, and the [Kr]4d 10 inner shells of Xe were replaced with pseudopotentials.The Kohn-Sham (KS) equations were solved to selfconsistency with a tolerance <5 × 10 −7 Hartree, and the grid spacing was converged to x = y = z = 0.4 bohr, such that the total energy per electron was converged <10 −3  Hartree.
The molecular structure of CO was relaxed to <10 −4 Hartree/bohr in forces within the LDA, leading to a bond length of 1.136 Å, compared to the experimental value of 1.128 Å (0.7% error).For atomic species we found that the LDA with a SIC correctly reproduced the ionization potential of the valence p states, yielding 0.569, 0.507, and 0.443 Hartree for Ar, Kr, and Xe, respectively, comparable with the experimental values of 0.579, 0.514, and 0.445 Hartree, respectively (errors for each species are 1.7%, 1.4%, and 0.4%, respectively).In Ar, using PBE XC with a SIC lead to relatively similar results of 0.552 Hartree ionization potential (4.7% error from experimental value).For CO we found that the LDA with a SIC led to an ionization potential of 0.585 Hartree, comparable to the experimental value of 0.514 Hartree (13.8% error), while PBE XC with a SIC leads to better results of 0.524 Hartree (1.9% error).
For the noninteracting electron models, the exact same KS orbitals as found in the ground-state calculations were used as the initial states, and the effective potential was taken to be the KS potential of the ground-state DFT calculation (i.e., the potential of which the KS orbitals are eigenstates of).

TDDFT HHG Calculations
The time-dependent KS equations were propagated with a time step t = 0.05 a.u. and by adding an imaginary absorbing boundary potential of width 7.2 bohr.The initial state was taken to be the ground state of the system.The propagator was represented by a fourth-order Taylor expansion and was explicitly time-reversal symmetric.The grid size, absorbing potential, and time step were tested for convergence.The envelope function A(t ) was taken as a flat-top envelope with two-optical-cycle-long rise and drop sections and a fouroptical-cycle-long flat top (where the optical cycle is T = 2π/ω).The dipole acceleration [calculated from Eq. (3) in the main text] was filtered with a super-Gaussian mask.The HHG spectra were calculated by Fourier transforming the second derivative of P(t ), calculated numerically with an eighth-order finite-difference method.For each calculated HHG spectra, the integrated harmonic yield was calculated by integrating over each harmonic order from ω ∈ [nω − 0.5ω, nω + 0.5ω], where n is the harmonic order.The calculations were only performed for positive values of ε, since due to the reflection symmetry of all species along the x axis, the HHG response is fully symmetric for ε → −ε.In the spectrograms, steps of ε = 0.02 were used, and plots are presented with parabolic interpolation.For phase scans presented in Appendix 6 of this Appendix, a step of φ = 0.025π was used.Appendix 7 shows the ellipticity data of the HHG emission, calculated directly for each harmonic order at its peak using Stokes parameters.
For the frozen XC potential calculations, exactly the same procedure as described above was utilized except that the XC part in the time-dependent KS equations was kept frozen to its ground-state value during propagation.

Noninteracting Electron HHG Calculations
For calculations in the noninteracting electron models, exactly the same procedure as described above for TDDFT was used.The only difference is that both the Hartree term and the XC term in the KS potential were kept frozen to their ground-state value during propagation.In this case, the KS equations are fully decoupled as described in the main text, leading to a noninteracting electron model with an effective potential that is equal to the ground-state KS potential.

SFA Calculations
In later sections of this Appendix, SFA calculations of the bielliptical HHG response from Ar are presented for a few cases in order to compare to other methods.These calculations were performed following the technical approach presented in Refs.[74,75] and using transmission-matrix dipole elements of hydrogenic 1s orbitals (even though Ar has a valence p shell).

Additional Results From Argon
We present here complementary results for the bielliptical HHG spectrograms from atomic Ar.First we present the total ionization yield vs time in the HHG response of atomic Ar for two different levels of theory-TDDFT using the LDA with a SIC, and a corresponding noninteracting electrons model.This is shown in Figs.3(a) and 3(b) for two values of pump ellipticities, ε = 0, 1.In the TDDFT calculations, the ionization rate is much lower than in the noninteracting electrons calculation.This is a result of the increasing positive charge on the ionized system, which effectively raises the ionization potential for the valence states as the calculation progresses (since both calculations use the same initial orbitals and initial ionization energies).Figures 3(a) and 3(b) show that early in the calculation the total ionized yields in both models match.At the point in time where roughly 0.01e charge has been ionized, the ionization rates depart.Even though the total ionization yields are small (less than 0.2e units of charge are ionized in all cases, which corresponds to less than 3.5% of each orbital), there are still discrepancies between the calculations since there is exponential sensitivity in the ionization step to the effective potential felt by the electrons.The magnitude of this effect likely depends on the level of theory used in calculations.
Next we consider the bielliptical HHG spectrograms calculated from atomic Ar at a higher intensity of 4 × 10 14 W/cm 2 compared to results presented in the main text [which are presented for a power of 10 14 W/cm 2 in Fig. 1(b)] using several levels of theory.Figure 4(a) shows results from TDDFT calculations (using the LDA with a SIC) compared to a corresponding SFA calculation.(The noninteracting electrons calculations are not presented, since at this laser power the degree of ionization is too high, causing numerical artifacts in our limited size grid.)Notably, at this intensity the SFA calculations match very well with the full TDDFT calculation.Very small deviations between the levels of theory are observed in terms of the maximizing/minimizing ellipticities (smaller than ε = ±0.04),while some larger differences are seen in the relative peak intensities (mostly for higher pump ellipticities).The extremely good correspondence between SFA and TDDFT calculations in this case can be attributed to a type of generalized SFA, where for high enough laser powers the e-e interaction term becomes negligible during propagation compared to the laser-matter term, and thus does not affect the dynamical response.This effect is even more pronounced for higher harmonics, as expected.We note that SFA calculations performed at the lower laser power, on the other hand, produce results that are indeed slightly offset from the full TDDFT calculations, with the same order of deviations as seen from the noninteracting electrons model [see Fig.We further present bielliptical HHG responses from atomic Ar under the same conditions as in the main text Fig. 1(a) but using two different XC functionals: LDA and PBE [see Fig. 4(c)].We note that the long-range asymptotic part in both cases is corrected by the same SIC; hence any difference in the response is likely a result of the XC descriptions of the bound part of the electronic density.In this case the response is almost identical using the two XC functionals, even though there is a slight difference in the ground-state ionization potential in each case (8.3% difference in ionization potential).In fact, the only noticeable differences are seen in harmonics 20 and 22 in Fig. 4(a), showing slight variations in relative spectrogram peak intensities.These in turn can be used in comparison with experiments to determine which functional is more appropriate.

Phase-Resolved HHG Spectrograms
In this section we present results in a similar nature to those shown in Fig. 1(b) in the main text but where rather than scanning the ellipticity in the pump fields, we scan the relative two-color phase (φ) in order to produce spectrograms (for the cross-linear geometry, where ε = 0).Results are shown in Fig. 5 for Ar, comparing the TDDFT calculation to the noninteracting electrons model.As seen, a similar type of sensitivity to e-e interactions is observed (i.e., the maximizing phase is sensitive to the inclusion of e-e interactions in the calculation).However, in this case the sensitivity is reduced compared to the bielliptical spectrogram.Practically, the reduced sensitivity can be attributed to the simpler structure of the phase-resolved spectrogram that has only one dominant peak per each harmonic order.Physically, it originates from the fact that the phase scan alone does not cover a wide range of HHG responses as in the bielliptical case (e.g., for any value of phase there are only two interfering electron trajectories per optical cycle).Nonetheless, the fact that there is a sensitivity to the level of theory in this spectrogram indicates that the approach in the main text can be extended to multidimensional spectroscopy, where other parameters of the laser can be scanned: total laser power, relative power between the beams, relative of the beams, elliptical major axes of the beams, and their relative phase.

Bielliptical HHG Ellipticity Spectrograms
We present here results for the bielliptical HHG spectrograms from both atomic and molecular media but where the emission ellipticity is analyzed rather than the intensity.Results are shown for several different species and using various level of theory.First, Fig. 6(a) presents calculated HHG ellipticity spectrograms from atomic Ar, Kr, and Xe at FIG. 6. Bielliptical HHG ellipticity spectrograms from noble gases.(a) Ellipticity-resolved HHG response from several atomic species using TDDFT calculations (within the LDA with a SIC).(b) Same as (a) but from atomic Ar and using various levels of theory.Calculations in (a), (b) are performed under similar conditions to Fig. 1 in the main text.(c) Same as (a) but from atomic Ar compared to He at a higher laser power of 4 × 10 14 W/cm 2 .The ellipticity is directly plotted on a scale of −1 to 1 for each harmonic order (i.e., including helicity).
the same level of theory (using full TDDFT within the LDA, with a SIC).As seen in Fig. 6(a), practically the same response is seen in all atomic species, except for a few small deviations.This is despite the fact that the atoms have slightly different ionization potentials.Similarly, Fig. 6(b) presents results for Ar using different levels of theory (TDDFT, noninteracting electrons, and SFA), showing that typically the same response is observed in all cases, though slight variations can be seen in higher harmonics.Results of the same nature are obtained from Kr and Xe (not presented).
We further investigate this behavior by performing similar calculations on atomic helium that has a much larger ionization potential (at 0.903 Hartree).In this case, again the same ellipticity response is observed for all harmonic orders of both He and Ar [see Fig. 6(c)].Remarkably, this result is independent of the fact that many of the harmonics which are above the ionization potential for Ar (nonperturbative HHG, above the 11th harmonic) are below the ionization potential of He (perturbative regime, below the 15th harmonic).
We conclude that in atomic media, the ellipticity response from bielliptical HHG is apparently a quite pure function of the laser parameters and is largely independent of the atomic species and the level of theory (for the explored laser parameters, species, and theory levels).We note on the other hand, that the intensity of the HHG emission does show a relatively large dependence on these degrees of freedom, as presented throughout the main text and the Appendix.Furthermore, the ellipticity response itself is much more sensitive to changes in the pump ellipticities than the intensity response-the ellipticity spectrograms in Fig. 6 oscillate very rapidly along the ε axis.Therefore, it is overall a very peculiar result that on one hand the intensity of the emission can be used to probe a great deal of properties of the interacting medium, while the ellipticity of the emission cannot.This result is even more confusing considering that both the intensity and the ellipticity of the HHG emission arise from the same fundamental quantities-the emitted HHG field components polarized along the x and y axes.We do not yet understand FIG. 7. Bielliptical HHG ellipticity spectrograms from aligned CO with varying levels of theory (full TDDFT using the LDA with an added SIC and a corresponding noninteracting electrons model).Calculation is performed under similar conditions to that in Fig. 2 in the main text.The ellipticity is directly plotted on a scale of −1 to 1 for each harmonic order (i.e., including helicity).
the physical source of this behavior, but hypothesize that it has to do with the phase picked up by the electrons as they are accelerated in the continuum.
Lastly, Fig. 7 presents bielliptical HHG ellipticity spectrograms from aligned CO, showing that in the molecular case the ellipticity response can indeed be utilized to probe e-e interactions.Here, the ellipticity shows sensitivity on the same scale as the total intensity spectrograms, where in the noninteracting electrons model the ellipticity oscillates much more rapidly vs ε compared to the full TDDFT calculation.Thus, the response is highly sensitive to the inclusion of e-e interactions, especially in the sign of the helicity as ε is tuned away from zero, which may be positive for the TDDFT calculation but negative in the noninteracting electrons model (see harmonics 11-21 in Fig. 7).

FIG. 1 .
FIG. 1. Bielliptical HHG spectrograms from atomic Ar, Kr, and Xe, calculated with and without dynamical e-e interactions.(a) Illustration of the driving electric field's time-dependent polarization: by tuning ε, the polarizations of both beams are controlled, resulting in the presented Lissajou curves, where arrows indicate the direction of time.(b) Integrated intensity per harmonic order vs ε from Ar, with a full TDDFT model [using the local density approximation (LDA) and a self-interaction correction (SIC)] [65], compared to the corresponding noninteracting electrons model and the corresponding calculation where the XC potential is frozen.The intensity of each harmonic is normalized to its maximum along the ε axis for clarity.Selected harmonics are shown.(c), (d) Same as (b), but for atomic Kr and Xe, respectively.Calculations correspond to a laser power of 10 14 W/cm 2 , and 2 = 1/3, φ = 0 (see Appendix 1-3 for details).

FIG. 3 .
FIG. 3. (a), (b) Total electron ionization yield vs time from Ar in similar settings to Fig. 1(b) in the main text for two values of pump ellipticities, respectively, and comparing TDDFT calculations (using the LDA with a SIC) to a noninteracting electrons model. 4(b)].

FIG. 4 .
FIG.4.(a) Same as Fig.1(b) in the main text but for higher laser power of 4 × 10 14 W/cm 2 and comparing TDDFT calculations within the LDA with an added SIC to SFA calculations.(b) Same as (a) but for lower laser power of 10 14 W/cm 2 , and also comparing to the corresponding noninteracting electrons model.(c) Same as Fig.2(c) in the main text but from atomic Ar-comparison of the bielliptical HHG spectrograms calculated with TDDFT using either LDA or PBE XC with an added SIC.