Terahertz driven extremely nonlinear bulk photogalvanic currents in non-resonant conditions

We report on the generation of bulk photocurrents in materials driven by non-resonant bi-chromatic fields that are circularly polarized and co-rotating. The nonlinear photocurrents have a fully controllable directionality and amplitude without requiring carrier-envelope-phase stabilization or few-cycle pulses, and are generated with photon energies much smaller than the band gap (reducing heating in the photo-conversion process). We demonstrate with ab-initio calculations that the photocurrent generation mechanism is universal and arises in gaped materials (Si, diamond, MgO, hBN), in semi-metals (graphene), and in two- and three-dimensional systems. Photocurrents are shown to rely on sub-laser-cycle asymmetries in the nonlinear response that build-up coherently from cycle-to-cycle as the conduction band is populated. Importantly, the photocurrents are always transverse to the major axis of the co-circular lasers regardless of the material's structure and orientation (analogously to a Hall current), which we find originates from a generalized time-reversal symmetry in the driven system. At high laser powers (~10^13 W/cm^2) this symmetry can be spontaneously broken by vast electronic excitations, which is accompanied by an onset of carrier-envelope-phase sensitivity and ultrafast many-body effects. Our results are directly applicable for efficient light-driven control of electronics, and for enhancing sub-band-gap bulk photovoltaic effects.

resonant or near-resonant transitions, (iii) their amplitude is insensitive to the two-color phase and is finelytunable by changing the pulse duration, (iv) their directionality is controlled by the two-color phase regardless of the material system.The photocurrents are shown to arise from excited-state occupation imbalances in the Brillouin zone (BZ) that build-up from cycle-to-cycle, and analogously to anomalous Hall currents [9] are always transverse to the main axis of the two-color field.However, in contrast to anomalous Hall currents no bias or magnetic field is required in our case.We explain the origin of this effect from intuitive physical arguments based on electron sub-cycle dynamics, and by a generalized time-reversal symmetry in the laserdressed system.Lastly, we show that for higher laser powers the generalized time-reversal symmetry is spontaneously broken by strong excitations to the conduction band (CB) that cause a breakdown of the independent-particle approximation, and also lead to CEP sensitivity.These results should be useful for alloptical ultrafast control of electron dynamics, for energy conversion with sub-band-gap light in bulk and spatially uniform materials, and for developing novel ultrafast spectroscopy techniques based on nonlinear photocurrents.
We begin by introducing the ω-2ω laser field that comprises two co-rotating circularly polarized beams: where E 0 is the electric field amplitude of the ω beam taken here in the range of 0.085-0.85V/Å (equivalent to laser powers of ~10 11 -10 13 W/cm 2 ), f(t) is a dimensionless envelope function (see SI for details), ω is the fundamental frequency, ϕ is the ω-2ω relative phase, Δ is the amplitude ratio between the beams, and  ̂+ is a right circularly-polarized unit vector.Note that we have employed the dipole approximation, which is valid in our wavelength ranges.A representative Lissajous curve for the time-dependent polarization of this ω-2ω field is presented in Fig. 1(a).Its main characteristic features are: (i) it exhibits a mirror symmetry along its x-axis, (ii) the field peak power and time-dependent polarization substantially differ along the positive and negative ends of the x-axis, i.e. at the start and middle of an optical cycle where the two fields constructively/destructively interfere (corresponding to t=0 and t=T/2, respectively, see Fig. 1(a)).Importantly, the direction of the field's mirror axis is CEP-independent (it only depends on ϕ).The interaction of E(t) in eq.(1) with matter is described from first principles by using time-dependent density functional theory (TDDFT) within the adiabatic local density approximation (LDA) in the real-space code, octopus [37][38][39], with which the time-dependent electronic current expectation value, J(t), is calculated (all technical details are delegated to the SI).
As a starting point, we analyze driven electron dynamics in diamond subject to E(t) in eq. ( 1).Diamond represents a benchmark system since it has a simple lattice structure with two atoms per unit cell, and is a wide band gap insulator with an indirect gap of 5.48 eV and a direct gap at the Γ point of 7.3 eV [40,41].These points easily allow extensive calculations with wavelengths well below the band gap (within the LDA the direct gap is 5.65 eV, slightly lower than the experimental value).Figure 1(b) presents the calculated current where the with E(t) is polarized in the diamond (111) planes and the field is comprised of 2500nm-1250nm light in the strong-field regime (I 0 =7.5×10 12 W/cm 2 ).A residual photocurrent is observed along the y-axis after the laser pulse ends (axes are denoted with respect to the field's coordinate system).This effect cannot result from CEP-sensitivity since the pulse is of long duration with a zero time integral over the vector potential.Moreover, the driving frequencies are much below the band gapit would take over 11 ω photons to induce a transition from the valence band (VB) to the CB.Thus, this effect substantially differs from previously discussed perturbative and/or resonant phenomena [24][25][26][27][28]. Importantly, the current does not decay in our simulation (see Fig. 1(b)).This is true even in the presence of electron-electron interactions incorporated in the simulation, indicating that electrons occupy a true excited state of the system (similar behavior occurs with a Perdew-Burke-Enzerhof (PBE) functional [42], see SI).The current is thus only expected to decay due to scattering channels not included in our model [43].To understand the mechanism behind this effect we explore its dependence on the various laser parameters.Figure 2(a) presents the calculated photocurrent amplitude vs. the applied laser power, I 0 , which multiplies both the ω and 2ω components.As expected, the current originates from nonlinear interactionsthere is an onset at I 0 =10 12 W/cm 2 which increases initially parabolically with I 0 , indicating that the current is proportional to four-photon processes in the perturbative regime.Indeed, four photons are the minimal number to traverse the band gap for this wavelength (1200nm) while still mixing both photon types (corresponding to annihilation of either three 2ω photons and one ω photon, or two 2ω photons and two ω photons).This mixing is crucial since without it there is no symmetry breaking as in Fig. 1(a) that allows photocurrents to be created (because the field components are individually circularly-polarized with no preferential direction in space).
For powers above I 0 ≈5×10 12 W/cm 2 in Fig. 2(a) the parabolic behavior changes to linear proportionality to I 0 , paradoxically suggesting that two-photon processes are dominant.Since two photons cannot promote electrons to the CB, higher nonlinear effects are at play.Further increasing the power changes this behavior to an unpredictable pattern that also induces nonzero x-polarized currents.We will explore this behavior in detail below, but first we will explain the origin of the photocurrent in the intermediate strong-field regime (for 10 12 <I 0 <10 13 W/cm 2 ). Figure 2(b) presents the current direction vs. the orientation of the major axis of E(t) in the (111) planes.The photocurrent remains transverse to the x-axis defined by the Lissajous in Fig. 1(a).This is despite the fact that different local lattice symmetries are probed as the field rotates (e.g.breaking local reflection or rotational symmetries).Figure 2(b) further shows that throughout this control, the current maintains its amplitude and only exhibits mild modulations.Similar results are obtained when E(t) is polarized in the (110) planes that have a different local symmetry (4-fold rather than 6-fold, see SI).These features thus seem analogous to those of anomalous Hall currents, though notably here the laser field has components along both x and y directionsit is just the main mirror axis of the Lissajous that dictates the direction of the photocurrent.Figure 3(a,b) presents the integrated and normalized electron occupations in the CB after the laser pulse ends.That is, it shows the occupations along k x ',k y ',k z ' axes in k-space after integrating over the other two axes (see SI for details).Note that the 'prime' coordinate system denotes the diamond primitive unit cell coordinates rather than those of the Lissajous, thus a current along the y-axis of the field is expected to show up as imbalances in occupations along all three k x ',k y ',k z ' axes.The occupations around the Γ point are asymmetricelectron occupations break the inversion symmetry in k-space which is generally a consequence of time-reversal symmetry in the Hamiltonian.Since E(t) breaks time-reversal symmetry, k-inversion symmetry is broken as expected, allowing photocurrents to be generated.We also note that almost all k-points in the BZ contribute to this asymmetry, even far away from the Γ point.
To understand why photocurrent components parallel to the x-axis are not generated, we analyze the dynamics of electrons driven in the strong-field regime assuming that the current is generated by two sequential steps: (i) First, the strong laser field excites an electron from the VB to the CB.This step is most likely to occur near the field peak power and near the Γ point (see highlighted points in Fig. 1(a)).(ii) Second, the excited electron is driven in the CB by the laser field, generating intraband currents and occupying different regions in the BZ.Since the Lissajous in Fig. 1(a) is symmetric about its x-axis, one expects that roughly an equivalent occupation of states with positive/negative k y is formed during the acceleration in step (ii).However, because the events that initiate tunneling are very different, a discrepancy is formed in step (i) trajectories that are initiated for positive k y are favorable, leading to a residual current.Since no such discrepancies exists along the y-axis of the field (e.g., at t=T/6 and t=5T/6, see Fig. 1(a)), no current is generated along the x-axis.
A different explanation for this effect follows from the analysis of the spatio-temporal symmetries of E(t) [44].Due to its mirror axis (see Fig. 1(a)), E(t) is invariant under a coupled operation of time-reversal accompanied by a reflection of its y-polarization: E(t)=σ x ×E(-t) (here σ x denotes a reflection E y →-E y ).Since for k-space occupations time-reversal is equivalent to inversion, this symmetry dictates a mirror plane along the k yz -plane (the product of inversion and σ x ).That is, occupations must be symmetric for k x →-k x due to a generalized time-reversal exhibited by the field, forbidding x-polarized photocurrents (note that 'x' labels the coordinate system of the field).We note that this analysis is formally valid only if neglecting other effects that may lead to symmetry breaking, e.g. the finite pulse duration or the local lattice structure.Nevertheless, the numerical results are consistent with this conclusion for a wide parameter range and different materials.This indicates that in this regime the structure of the laser field is crucial.Figure 3(c) shows the calculated photocurrent amplitudes for different pulse durations for a fixed laser power.The current increases linearly with the pulse durationthe asymmetries in k-space occupations buildup coherently from cycle-to-cycle, allowing a finely-tuned control over the current and further emphasizing the importance of sub-cycle dynamics.
Having established the main results in diamond, we explore other systems (most results are delegated to the SI for brevity).First, we show that similar effects arise in silicon, a common photovoltaic energy conversion material [45] (see Fig. S4(a) in the SI).Second, we show that nonlinear photocurrents arise also in MgO (see Fig. S4(b) in SI), which is a wide band gap oxide with primarily ionic bonding.Hence, the current generation is largely independent of the solid's characteristics.Third, we investigate a monolayer of hBN subject to E(t), showing that similar photocurrents arise (see Fig. S5 in SI).We note that this result differs from the ones obtained with a counter-rotating bi-circular field that does not generate a photocurrent [19,46].Lastly, we explore driven electron dynamics in graphene.The interactions in graphene are resonant due to its metallic nature.Nonetheless, photocurrents arise just as has been shown experimentally for the few-cycle case [33,34], but with stronger asymmetries in CB occupations (see Fig. S6 in the SI).Remarkably, the K and K' points in the BZ are still connected by a generalized time-reversal symmetry, such that their total integrated occupations are identical, but the occupation patterns allow a photocurrent to be generated.
Finally, we discuss the extreme regime for I 0 >10 13 W/cm 2 (which is experimentally accessible) in diamond (Fig. 2(a)).Here x-polarized photocurrents are also allowed.Moreover, the current increase vs.I 0 is erratic and indicative of highly nonlinear effects.We find that several phenomena occur that contribute to this behavior and result from spontaneous symmetry breaking to the generalized time-reversal symmetry (that forbids the x-polarized currents for I 0 <10 13 W/cm 2 ).In this regime the excitation from the first VB to the CB is so fast and vast that the VB maxima at the Γ point has largely depleted even before the laser pulse ends (see Fig. 4(a)).Thus, as the system evolves electrons must traverse larger gaps that are contributed from farther k-points, which breaks time-reversal symmetry.Intuitively, this leads to a sub-cycle asymmetry that allows x-polarized photocurrents, since occupations in the first half of an optical cycle no-longer exactly cancel with those in the second half.In an extreme case, the VB may be fully depleted in a single half-cycle, resulting in very large x-photocurrents.We note that it is unclear if all, or just some, of these effects are observable below the material damage threshold (which can be manipulated by changing the laser parameters or the material).
Besides the onset of x-polarized photocurrents, two other effects come into play.First, dynamical electron-electron interactions play a major role.The excitation of many electrons to the CB significantly renormalizes the effective exchange and correlation interactions from cycle-to-cycle.Figure 4(b) shows photocurrent amplitudes calculated by two different levels of theoryeither freezing or not the Hartree and exchange-correlation (HXC) potentials to their ground-state form (the frozen HXC approximation is equivalent to the independent-particle approximation (IPA)).Discrepancies between these theories only arise due to the inclusion of electron-electron interactions induced by the laser, and are usually very small [47][48][49][50][51][52].Below 10 13 W/cm 2 these are negligible (<5%).Above 10 13 W/cm 2 there are significant deviations in the obtained photocurrent (>20%, see Fig. 4(b)), as well as in the intensity of harmonic emission (>50%, see Fig. S3 in the SI).The deviations increase with the field strength, in accordance with the CB population (e.g.compare Fig. 4(a) to Fig. 4(b)).We note that it is not yet clear whether or not this effect is well-described within adiabatic TDDFT, thus results could be affected by the level of theory describing electron-electron interactions.Second, otherwise negligible envelope effects become significant.Figure 4(c) shows the photocurrent amplitudes in this regime vs. the CEP, validating that there is strong CEP sensitivity even though the pulse durations are relatively long, and no such sensitivity is observed for weaker power (see SI Fig. S2).Both of these phenomena could be probed by measuring the photocurrent and its direction in space. .Inset shows the electron occupation on the VBs during interaction with the laser pulse for the exemplary case of I0=5×10 13 W/cm 2 , showing that the VB is quickly depopulated even during the turn-on of the laser pulse.(b) Deviation in photocurrent amplitudes between the IPA and the TDDFT calculation vs. field power (for λ=1200nm, ϕ=0, Δ=√3, I0=5×10 13 W/cm 2 ).(c) Calculated photocurrent amplitudes vs. CEP, showing that CEP-sensitivity arises in the extremely strong-field regime.Calculated for the same parameters as (a), except with Δ=√2.
To summarize, we have investigated the interaction of solids with a bi-chromatic ω-2ω laser with corotating circularly-polarized components.We showed that photocurrents are generated in the bulk without external bias, and due to a coherent excitation of electrons to the CB that builds-up from cycle-to-cycle.A generalized time-reversal symmetry exhibited in the system imposes that only transverse photocurrents are generated, allowing ideal control over the current amplitude and directionality.This effect is general to threeand two-dimensional systems, gapped insulators, semi-metals, and does not strongly depend on the material properties.We also analyzed an extreme strong-field regime where time-reversal is spontaneously broken, leading to ultrafast laser-induced correlations, as well as CEP-sensitivity from long pulses.Both these phenomena pave the way to novel forms of ultrafast spectroscopy for electron dynamics that rely on nonlinear photocurrents.Importantly, our results demonstrate a method for sub-band-gap bulk photovoltaic energy conversion that is not bounded by any theoretical limit.Looking forward, these ideas could also be employed to enhance the photogalvanic effect in perovskites [22], and for improved control in petahertz electronics [53].

ASSOCIATED CONTENT
Supplemental material is available for this paper that includes details on the methodology used in calculations as well as some additional results.system to Δx=Δy=Δz=0.33 Bohr in diamond, 0.35 Bohr in Si, 0.3 bohr in MgO, 0.4 Bohr in monolayer hBN, and 0.4 Bohr in graphene.We employed a Γ-centered k-grid in each system, which converged the timedependent current expectation value: 36×36×36 k-grid in diamond, 24×24×24 k-grid in Si, 32×32×32 k-grid in MgO, 36×36×1 k-grid in monolayer hBN, and 120×120×1 k-grid in graphene.
For TDDFT calculations, we solved the time-dependent KS equations within the adiabatic approximation, represented in real-space and in the velocity gauge, given in atomic units by: , the time-dependent electron density.The KS wave functions were propagated with a time step of Δt=0.2 a.u.which converged the time-dependent current in all material systems within the LDA.For TD-PBE calculations a time step of 0.08 a.u. was used.The initial state was taken to be the system's ground state.The propagator was represented by a Lanczos expansion.
For the independent particle approximation (IPA) calculations, the same methodology was utilized but where the KS potential was kept frozen to its initial form such that   [(, )] ≡   [(,  = 0)].
For the two-dimensional systems of monolayer hBN and graphene, we added additional vacuum spacing above and below the monolayer of 40 Bohr in each direction of the non-periodic axis, and absorbing boundaries were employed along this axis with a width of 12 Bohr.
The time-dependent current expectation value was calculated directly from the time-dependent KS states as: where (, ) is the microscopic time-dependent current density: , and Ω represents the volume integral over the primitive cell.The photocurrent was calculated as  = ( → ∞) when the laser pulse has ended.For results presented in the SI, the harmonic spectrum was calculated as the Fourier transform of the first derivative of the current: , which was evaluated numerically with an 8'th order finite-difference approximation for the temporal derivative, and fast Fourier transforms (FFT).Band and k-point occupations were computed by calculating the projections of the time-dependent KS states on the field-free states at t=0.For Figures 3(a,b) in the main text, the occupations were integrated in kspace and normalized to the density of states at each k-point, g(k).
The envelope function of the employed laser pulse, f(t), was taken to be of the following 'super-sine' form [51]: where σ=0.75, T p is the duration of the laser pulse which was taken to be T p =8T (unless stated otherwise), where T is a single cycle of the fundamental carrier frequency.This form is roughly equivalent to a supergaussian pulse, but where the field starts and ends exactly at 0 which is more convenient numerically.

II. ADDITIONAL RESULTS IN DIAMOND
We present here additional results complementary to those presented in the main text.First, we show that the photocurrents arise even when utilizing a PBE XC functional, and similarly to the TD-LDA calculations presented in the main text, they do not decay over the timescale of the simulation.Figure S1(a) shows the time-dependent current expectation value calculated in diamond with PBE XC.The current arises only along the y-axis transverse to the field's mirror plane and remains constant in time even long after the pulse has ended.This supports the conclusions in the main text that the excited current-carrying state is an eigenstate of the many-body system and will only decay via radiative or scattering channels that are not included in our simulation.
Next, Fig. S1(b) shows the calculated time-dependent current expectation value driven in a diamond lattice within the (110) planes (unlike figures in the main text that present currents driven in the (111) planes).The (110) planes have a local four-fold symmetry that is different than that in the (111) planes.Still, the current arises along the field's y-axis, demonstrating similar control over the current direction.This supports the conclusion that the current generation mechanism is largely unaffected by the local lattice structure and symmetry.We next show that these currents are CEP-independent, just as one would expect for long-duration pulses.Fig. S2 shows calculated photocurrent amplitudes in diamond vs. the CEP for the same pulse durations as presented in the main text in Fig. 4(c).The resulting photocurrents are CEP-independent.This result further highlights that different physical effects come into play in higher laser powers, as discussed in the main text.
Lastly, we explore here deviations from the IPA for strongly-driven diamond in the regime of laser powers above 10 13 W/cm 2 .Results in the main text (Fig. 4(b)) showed that there are considerable deviations in the obtained photocurrents upon inclusion/exclusion of dynamical electron-electron interactions (up to 20%), as well as an onset of x-polarized photocurrents.We show here that in these conditions similar effects arise in the HHG emission.In this case deviations are even larger and can reach values of 100% for some harmonics, indicating the breakdown of the IPA. Figure S3(a) presents the calculated HHG spectra in this regime from the co-circular ω-2ω field for I 0 =5×10 13 W/cm 2 , showing these extreme deviations.Figure S3(b) shows that when averaged over the first 60 harmonic orders, these deviations increase in correspondence with the laser power, similar to the photocurrent deviations.Overall, this further validates the onset of strong laserinduced many-body effects in this intense regime.

Figure S3
. Breakdown of the IPA in strongly-driven diamond in the HHG spectrum.(a) Calculated HHG spectra within the IPA and full TDDFT calculation for laser conditions λ=1200nm, Δ=√3, I0=5×10 13 W/cm 2 .Gray lines indicate positions of integer harmonic orders for this wavelength.The plot shows that considerable deviations arise in the HHG power in this regime due to dynamical correlations.(b) Average deviation in harmonic power between the full TDDFT calculation and the IPA vs. the laser power (averaged over first 60 harmonics).

III. ADDITIONAL RESULTS IN OTHER SYSTEMS
We present additional results of photocurrent excitation in other material systems.We start by exploring several three-dimensional bulk solids.Figure S4 presents calculated time-dependent photocurrents in Si (in the (111) planes), and MgO (in the (111) planes), after interaction with intense co-circular ω-2ω pulses.These results are analogous to those presented in the main text for diamond and highlight that the photocurrent and its generation mechanism is largely independent of the material system and its chemical and physical properties.In particular, the current is always observed transverse to the field's mirror axis, and is generated even though the laser pulses are of long duration and have photon energies much smaller than the band gap in each case.We next investigate photocurrents in two-dimensional systems.Figure S5 presents investigation of photocurrent generation in monolayer hBN where the field is polarized within the monolayer (xy plane).From Fig. S5(a) we can tell that similar photocurrents arise, and that they are transverse to the field's mirror axis.We note that this occurs even though hBN monolayers lack inversion symmetry.Figure S5(b) shows the kspace occupations of the CB after the laser pulse has ended, clearly indicating that there is no k-inversion symmetry (i.e. a residual current is generated).Notably, in this case no symmetries remain in the k-space occupations, but the current is still mainly polarized along the y-axis, indicating that its directionality is largely determined by the field's properties rather than those of the material system.We also note that, as expected, the occupations of the K and K' valleys are now asymmetric, allowing for a mechanism to generate valley polarization.This is in accordance with similar results presented with counter-rotating bi-circular fields [19].We move on to graphene, which is a hexagonal Dirac semi-metal with band touching at the K and K' points.Previous theoretical and experimental work showed that one can generate strongly-driven photocurrents in graphene with monochromatic few-cycle laser pulses [33][34][35][36].Figure S6 explores current injection from the co-rotating ω-2ω field, showing that similar effects occur in the long-pulse regime (the field is polarized within the graphene sheet).This effect is seen both for shorter wavelengths (1200-600nm, see Fig. S6(a,b)), and for much longer wavelengths in the THz regime (2500-1250nm, see Fig. S6(c,d)).We note that for all tested wavelengths the k-space occupations exhibit an exact mirror symmetry, in accordance with the formal analysis presented in the main text (and because graphene is inversion-symmetric).Thus, the total populations of the K and K' valleys are identical and connected through a mirror symmetry, which is mediated by the generalized time-reversal symmetry of the ω-2ω co-rotating field.Nonetheless, photocurrents arise because the occupation patterns are extremely asymmetric between the K and K' valleys (this is also evident locally in the region near each K and the K' band touching point).

Figure 1 .
Figure 1.Bulk photocurrent generation in diamond by THz co-circular ω-2ω fields.(a) Schematic Lissajous curve for the timedependent polarization of the field in eq.(1) (for ϕ=0, Δ=√2).Small arrows indicate the direction of time, and the red arrow indicates the direction of the photocurrent (I).Different instances in time are highlighted to emphasize the mirror symmetry of the Lissajous about its x-axis, but not the y-axis.(b) Calculated current expectation value (see SI for details) for field polarization is in the (111) planes (with λ=2500nm, ϕ=0, Δ=√2, I0=7.5×10 12 W/cm 2 ).
Photocurrent angle

Figure 3 .
Figure3.Conduction band occupations in diamond after interaction with ω-2ω co-circular laser, and photocurrent proportionality to the pulse duration.(a,b) CB occupations (denoted by ρ) along kx' and kz' axes, respectively, integrated over the other two axes, presented in the coordinate system of the lattice (ky' is not shown due to similarity to kx').ρ is normalized to the density of states, g(k').Purple curves show the CB occupation, while blue curves show the deviation from inversion symmetry (i.e. the occupations subtracted from their inverted image), indicating which states contribute to the photocurrent.Calculated for λ=1200nm, ϕ=0, Δ=√3, I0=7.5×10 12 W/cm 2 .(c) Photocurrent amplitude vs. pulse duration.Inset shows calculation for a 25-cycle pulse showing the build-up of photocurrent from cycle-to-cycle.Calculated for the same parameters as (a) except Δ=√2. Photocurrent(a.u.)

Figure 4 .
Figure 4. Spontaneous symmetry breaking in extremely nonlinear photocurrents in diamond.(a) Total excitation of electrons to the CB per primitive cell vs. I0.Inset shows the electron occupation on the VBs during interaction with the laser pulse for the exemplary case of I0=5×1013 W/cm 2 , showing that the VB is quickly depopulated even during the turn-on of the laser pulse.(b) Deviation in photocurrent amplitudes between the IPA and the TDDFT calculation vs. field power (for λ=1200nm, ϕ=0, Δ=√3, I0=5×1013 W/cm 2 ).(c) Calculated photocurrent amplitudes vs. CEP, showing that CEP-sensitivity arises in the extremely strong-field regime.Calculated for the same parameters as (a), except with Δ=√2.

Figure S1 .
Figure S1.Additional calculations of photocurrents in diamond.(a) Time-dependent current expectation value calculated using PBE XC, showing results similar to those obtained within the LDA (for laser parameters λ=1200nm, Δ=√3, I0=7.5×10 12 W/cm 2 ).(b) Timedependent current expectation value calculated where E(t) is polarized in the diamond (110) planes, showing that that the photocurrent is still directed along the field's y-axis regardless of the different local symmetry compared to the (111) planes (calculated for similar laser parameters as in (a), except that Δ=1).Arrows indicate the residual photocurrent.

Figure S5 .
Figure S5.Photocurrent generation in monolayer hBN.(a) Time-dependent current expectation value calculated for the laser parameters λ=2500nm, Δ=√2, I0=10 12 W/cm 2 ).Arrow and dashed line indicate the residual photocurrent.(b) CB occupations in kspace after the pulse has ended showing an asymmetric occupation pattern around the K and K' points.

Figure S6 .
Figure S6.Photocurrent generation in graphene.(a) Time-dependent current expectation value calculated for the laser parameters λ=1200nm, Δ=√2, I0=10 11 W/cm 2 ).(b) CB occupations in k-space after the pulse has ended corresponding to the parameters in (a).(c,d) Same as in (a,b), but for λ=2500nm.Arrows and dashed lines indicate the residual photocurrents.White dashed lines mark a mirror symmetry for k-space occupations that connects the K and K' points, and forbids x-polarized photocurrents.
is the KS-Bloch state at k-point k and band n, () is the vector potential of the laser electric field within the dipole approximation, such that −  () = (), c is the speed of light in atomic units (c≈137.036),and   (, ) is the time-dependent KS potential given by:   is the charge of the I'th nuclei and   is its coordinate,   is the XC potential that is a functional