Rabi Regime of Current Rectification in Solids

We investigate rectified currents in response to oscillating electric fields in systems lacking inversion and time-reversal symmetries. These currents, in second-order perturbation theory, are inversely proportional to the relaxation rate, and, therefore, naively diverge in the ideal clean limit. Employing a combination of the non-equilibrium Green function technique and Floquet theory, we show that this is an artifact of perturbation theory, and that there is a well-defined periodic steady-state akin to Rabi oscillations leading to finite rectified currents in the limit of weak coupling to a thermal bath. In this Rabi regime the rectified current scales as the square root of the radiation intensity, in contrast with the linear scaling of the perturbative regime, allowing to readily diagnose it in experiments. More generally, our description provides a smooth interpolation from the ideal Periodic Gibbs Ensemble describing the Rabi oscillations of a closed system to the perturbative regime of rapid relaxation due to strong coupling to a thermal bath.

Our study is motivated by the following question: what is the ultimate fate of current rectification in Bloch bands in the ideal limit where relaxation times become very large?As we will demonstrate, there is in fact a well defined periodic steady state in such a limit, that we will refer to as the "Rabi regime", in which the system sustains a finite DC rectified current.
A useful starting point to appreciate the nontrivialities of such a clean limit is to analyze the problem perturbatively in the amplitude of electric field, as commonly done in most studies (see however Ref. [12,[28][29][30][31]).Perturbation theory predicts a rectified current j, that grows as the square of the amplitude of the field, j ∝ E 2 .For frequencies above the threshold for interband transitions, such perturbative BPVE are often separated into two mechanisms known as the shift and the injection current effects [3-10, 12-17, 23-26, 29-31].The injection current originates from difference of the banddiagonal velocity of the empty and occupied bands at a given crystal momentum k.The shift current, on the other hand, originates from the difference of positions of Bloch wave-functions between the empty and occupied bands at a given k, and can be computed as the contribution arising from the band-off-diagonal velocity operator.
A crucial distinction between the shift and injection currents is that, within perturbation theory, the shift current appears to have a finite value in the "clean" limit of vanishing relaxation rate, Γ → 0, while the in- jection current appears to diverge in such limit as 1 Γ, which ultimately arises from the vanishing quasi-energy denominators appearing at higher orders of the perturbation theory for the rectified current (see Refs. [3-10, 12-17, 23-26, 29-31] and S.I.I A).Such divergence is often handled in an ad-hoc manner by computing the response of the rate of change of the current, dj dt, and assuming that such growth leads to a current saturation to a value proportional to the relaxation time τ ∼ ̵ h Γ.However, recently an interesting non-perturbative study of the CPGE in Weyl semimetals [30] demonstrated that the rectified current saturates to a finite value even in the clean limit of vanishing relaxations (Γ → 0) within a semiclassical kinetic framework.The underlying mechanism for such saturation is the Rabi dynamic broadening of absorption [32][33][34][35], which occurs when the energy scale controlling the transitions between conduction (c) and valence (v) bands exceeds the relaxation rate, eE ⋅ ⟨c r v⟩ ≫ Γ, which we refer to as the Rabi regime.
In the present work we develop a microscopic description of the currents for arbitrary values of the nonlinearity parameter eE ⋅ ⟨c r v⟩ Γ, that captures the perturbative and the Rabi regimes on an equal footing.To do so, we employ a Keldysh-Floquet formalism [36][37][38] in a generic two-band system coupled to an ideal fermionic bath, following the pioneering approach of Ref. [12] (see also Ref. [28]).As we will see, and contrary to the expectations of perturbation theory, in the Rabi regime, the traditional resonant shift current contributions vanish, whereas the injection currents approach a finite limit that scales as the square root of the radiation intensity in a sharp contrast to the perturbation theory expectations.We will also demonstrate that the Rabi regime can be viewed as an example of thermalizing synchronization of a system under an external periodic drive that can be described by the periodic Gibbs ensemble [39,40].
We derive the nonperturbative expression for currents within a two-band model (see Fig. 1(a)).The electric current operator is defined as: ĵ = ev ̵ h = ∂ Ĥ0 (k + eA(t) ̵ h) ∂A(t), where Ĥ0 (k) is the 2x2 matrix Bloch Hamiltonian, and A(t) is the vector potential from spatially uniform but time dependent electric field.Since the crystal momentum k is conserved, the problem is equivalent to a collection of independent driven two-level systems.We restrict our analysis to a monochromatic electric field with frequency ω: Here E is a vector with complex entries, allowing us to capture light of arbitrary degree of polarization, including the case of linear polarization, when all components can be chosen to be real, to fully circularly polarized light, when two orthogonal components differ by a phase of π 2. The periodicity in time allows to employ the Floquet picture (for details see S.I.I B) where multiple Floquet bands appear with a quasi-energy that is boosted by multiples of the driving frequency (see FIG. 1(a)).We simplify the problem by truncating the Floquet Hamiltonian to two bands that are in resonance, in the spirit of a rotating wave-approximation [13].This approximation is well justified when the off-diagonal terms in the Floquet Hamiltonian are smaller in comparison to the Floquet quasi-energy difference to other remote Floquet bands, namely when e E ⋅ ⟨c r v⟩ ≪ ̵ hω (see e.g.Ref. [41] and S.I.I B).Thus the approximate Floquet Hamiltonian is: where 1 stands for valence, 2 for conduction and 1,2 are effective valence and conduction band energies respectively (which could be dressed by higher order perturbative corrections with respect to bare band energies, as further discussed in S.I.I B).The subscript F stands for the representation of the operator in the Floquet picture, which is related to the ordinary Schrödinger picture as follows: In order to capture relaxation processes, we couple the system to a bath and apply the non-equilibrium Green function technique on the Keldysh contour (see S.I.I C and Refs.[13,36,[42][43][44][45][46]).We choose a simple model in which each fermionic site in the system of interest is coupled to its own fermionic bath, with a common hopping amplitude V mix (see Fig. 1(b)).The temperature of the bath is T bath = 1 (k B β) and the chemical potential is µ.The effective density matrix of the system is given by the lesser equal time Green Function G < (t, t), and can be shown to be (see S.I.I C): where ) are valence and conduction Fermi-Dirac occupation factors respectively, h ± = h x ± ih y and Γ = 2π V mix 2 is the relaxation rate.
The DC current of the system, J α = −ieTr Ĝ< F vα F ̵ h, can be decomposed into three contributions: Here α denotes the real space indices and the velocity operator in Floquet representation vα F is decomposed in the Pauli basis, namely vα F = ∑ i=x,y,z v α i σ i .In the supplementary (see S.I.I E) we compare the currents Eq.(5-7) with perturbation theory and show that Eq.( 5) and Eq.( 7) recover the resonant behaviour of the shift and injection currents respectively, whereas Eq.( 6) becomes the non-resonant component of the shift current in the limit eE ⟨c r v⟩ ≪ Γ.Now, to analyse the clean limit behaviour of the injection current Eq.( 7), it is useful to take the approximation in which both the diagonal h z matrix elements in Eq.( 2) are greater than Γ and the off-diagonal elements h x,y .This is typically well satisfied in most solids except in special situations such as resonant absorption on extremely flat bands, and we will demonstrate that this is a good approximation by explicit calculations later on.Therefore we replace 1 (h 2 + Γ 2 4) ≈ πδ(h z ) h 2 x + h 2 y + Γ 2 4, which leads to the following expression: where v 1,2 are conduction and valence band velocities respectively.We again see that in the limit of fast relaxation (eE ⋅ ⟨c r v⟩ ≪ Γ) Eq.( 8) reproduces the behaviour predicted by perturbation theory (see S.I.I A).Remarkably, however, in the clean limit (Γ → 0), the above formula predicts a finite current, in sharp constrast to the naive extrapolation of perturbative result.In other words, the relaxation rate in the denominator of the perturbative expressions acquires a non-perturbative modification by the driving electric field of the form: Therefore in the clean limit, the injection current scales as the absolute value of the electric field, J 3 ∝ E , and, accordingly, it is proportional to the square root of the radiation intensity.On the other hand, the term J 1 from Eq.( 5), which reduces to the usual resonant shift current from perturbation theory (see S.I.I A for details), can be seen to vanish in the clean limit Γ → 0 from Eq.( 5).This is noteworthy because in the perturbative regime (eE ⟨c r v⟩ ≪ Γ ) the shift current naively approaches a finite value in the Γ → 0 limit.
Synchronization and Rabi Limit of Rectification.While the Keldysh formalism allows for a description with arbitrary strength of coupling to the bath, there is a simpler way to understand the ideal behavior in the limit of vanishing coupling to the bath (Γ → 0).In fact, this limit can be understood simply as a form of Rabi oscillations associated with the inter-band transitions driven by the oscillating field.We will describe how to understand this limit within the picture of the Periodic Gibbs Ensemble (PGE) [39,[47][48][49] that captures the steady state synchronization of the system with the driving field.
Consider an initial state described by a density matrix ρ 0 .This density matrix can be decomposed in the eigenstates of the time dependent Hamiltonian, ψ α (t), and therefore the state at any later time t, is given by: where ρ αβ = Tr ρ 0 ψ α (t 0 )ψ † β (t 0 ) .Now, the Floquet theorem implies that, barring accidental degeneracies, the operators ψ α (t)ψ β (t) † are only periodic when α = β.The late-time synchronization associated with the PGE can be understood as a process in which the memory of these off-diagonal amplitudes of the density matrix in the Floquet basis disappears in a kind of thermalization process, leading to a steady state that is exactly periodic and synchronized with the drive (see S.I.I G): Remarkably the above ensemble is identical to the one that we have obtained within the Keldysh formalism in the limit of Γ → 0, when one chooses the initial state ρ 0 to be the equilibrium Fermi-Dirac density matrix in the absence of the periodic perturbation, with the chemical potential and the temperature of the bath (see S.I.I E. for more details).In fact, within the same rotating-wave approximation used to solve the Floquet problem, this density matrix is explicitly given by (see S.I.I G): This density matrix encodes the physics of Rabi oscillations (see S.I.I F for details).The above reduces exactly to the density matrix in Eq.(4) in the clean limit Γ → 0, once it is expressed in the Floquet picture (see Eq.( 3)), and therefore predicts the same rectification currents that we have previously described in the clean limit.
We would like to note that most studies of PGE to date have focused on what might be called "internal" synchronization, which considers a closed system acting as its own bath.In this context, the initial condition, ρ 0 , is freely chosen and it is not unique.In our context, however, the emergence of the PGE follows from different principles.Coupled to the bath, the system loses memory of its initial state at late times.It does so by flowing towards a unique stable periodic solution.Remarkably, in the limit of weak coupling to the bath, this steady state coincides exactly with one specifically chosen PGE, whose initial condition is the one associated with the thermal equilibrium system with an infinitesimal coupling to the bath in the absence of the periodic drive.
Therefore, although we have performed our calculations in a rather specific microscopic setting, we have been able to recover the universality of the PGE in the limit of weak coupling to the bath that we are using.Since the PGE can be justified under generalized entropy maximixation principles [39,[47][48][49], this is a compelling indication that our results describe the behavior of a large class of systems coupled to ideal heat baths.
We would like to illustrate this behavior for representative nodal fermions with Hamiltonians that are linear in momentum.These linear in k Hamiltonians have negligible shift currents (see details in S.I.I J) and therefore allow us to focus on the behavior of injection currents, which we will consider from here on in this section.We will consider two types of model that are relevant to a large class of materials.The first is an ideal 3D Weyl fermion, and our focus will be on the non-perturbative modifications to the CPGE.As we will see, our results are in perfect agreement with those obtained recently in Ref. [30].The second will be a 2D tilted Dirac massive fermion, and our focus will be to investigate the nonperturbative regime of rectification for linearly polarized light in a time reversal breaking system.
The ideal 3D Weyl Hamiltonian is: Here v 0 is a Fermi velocity and σα are Pauli matrices.This model respects TR but breaks inversion symmetry.
When the system has a finite chemical potential, light absorption occurs above a threshold frequency ̵ hω > 2 F (see Fig. 2(a)).By using the formula from Eq.( 8), one obtains the following non-perturbative approximate expression of the injection current above such a threshold (see S.I.I I for details): where E 2 = E * ⋅ E, with E understood as the complex vector defined in Eq.( 1) (see S.I.I I for comparison of this approximate formula against direct evaluation from the integral in Eq.( 8)).Eq.( 14) in the perturbative regime ) interestingly, the injection current approaches a value that is independent of the relaxation rate and it is given by: Here ζ ≈ 0.3 is a numerical pre-factor with a weak dependence on the degree of light polarization.Its value for perfectly circularly-polarized light can be computed exactly from Eq.( 8) to be ζ = 1 (2 √ 2), in agreement with Ref. [30] (see S.I.I I for details).The behavior of the rectified current in these two regimes and their crossovers are shown in Fig. 2(a,b).
We will now consider a 2D Dirac Hamiltonian given by: where m is the mass which breaks time-reversal symmetry, v x , v y are anisotropic Fermi velocities, and u x > 0 is the tilt term that breaks inversion.The above model features absorbtion within a window of frequency given by In this window the maximum current occurs when ̵ hω ≈ 2 F (see Fig. 2(c)), and for the electric field along the tilt direction.From Eq.( 8), the corresponding component of the injection current can be approximated as (see S.I.I I): Within perturbation theory this current would be 3,p to normalise the numerical non-pertubative results shown in Figs.2(c)-(d), so that deviations from 1 signal deviations from the perturbative regime.
Summary and experimental outlook.We have developed a formalism which captures on equal footing the perturbative regime of fast relaxation (eE ⟨c r v⟩ ≪ Γ) and the non-perturbative regime of strong light intensity (eE ⟨c r v⟩ ≫ Γ) of current rectification for interband transitions.In the perturbative regime, we recover the well-known behavior according to which shift currents approach a value that is independent of the relaxation rate Γ, while injection currents scale as 1 Γ.Interestingly in the opposite non-perturbative clean limit of slow relaxation (Γ → 0) the shift current vanishes, while the injection current approaches a finite value independent of Γ, but with a net current that scales as the square root of the radiation intensity, which can guide its identification in experiments.We have shown that this non-perturbative clean limit can be understood as optical Rabi oscillations synchronized with the incident radiation that realizes a time dependent generalized periodic Gibbs ensemble in a setting very different from its initial proposal.
Nodal Weyl semi-metals are promising platforms to realize the Rabi regime because their inter-band dipole matrix element diverges when approaching the Weyl node as ⟨c r v⟩ ∝ 1 k.As a consequence, they can access this non-perturbative regime above a light intensity that decreases with frequency, namely when ev 0 E > ̵ hΓω.For RhSi [65,66], using ̵ hΓ −1 ≈ 10ps [58] we estimate that the non-perturbative Rabi regime will be accessed at light intensities above 4 × 10 5 W/cm 2 for a photon energy of ̵ hω ≈ 0.5eV, but this required light intensity can be decreased as ω 2 at lower photon energies. Acknowledgements.
We wish to thank Achilleas Lazarides, Kin Fai Mak and Andrea Cavalleri for valuable discussions and correspondence, and to Nikita Leppenen and Leonid Golub for sharing their calculations that helped us understand the precise connection to their work in Ref. [30].We acknowledge financial support from the Deutsche Forschungsgemeinschaft through SFB 1143 (project-id 247310070) and cluster of excellence ct.qmat (EXC 2147, project-id 39085490).

A. Perturbation theory at small relaxations
As discussed in Refs.[14,17] the second order rectification conductivity in general can be separated into the following contributions: (18) where J stands for "Jerk", BCD for "Berry curvature dipole", I for "injection" and S for "shift current" respectively.We consider band structure with small relaxations (∀n, m, n ≠ m ∶ Γ ≪ n − m ).Each contribution is given by [14,17]: where Âα nm = Âα nm (1 − δ nm ) is off-diagonal Berry connection.
We split the conductivity above into resonant and non-resonant parts by separating the contributions into those that require a resonant condition that matches the light frequency with an energy difference in the limit Γ → 0 (namely those containing delta functions enforcing a Fermi's Golden rule), from those that are non-resonant and contain the principal parts where the frequency is not forced to match an energy difference.The "BCD" and "Shift" conductivity have both resonant and nonresonant parts, which we label by "R" and "NR" additional subscripts and are given by: where Now, their non-resonant part is: "Jerk" and "injection" components can be viewed as purely resonant, and are given by the following expressions: For purposes of comparing with the current nonperturbative formalism, we write the injection and shift currents of two-band systems system predicted by perturbation theory for the monochromatic electric field of the following form: Where DC current is defined as: the injection and shift currents are: BCD and Jerk rectification conductivities are low frequency Fermi surface terms, namely they vanish in the absence of Fermi surface at zero temperature and are associated with pole singularities at ω = 0, and therefore they are not described by the non-perturbative formalism of the main text that focuses on interband transitions.On the other hand, shift and injection can be nonzero for inter-band transitions.Notice that, interestingly, taken at face value, perturbation theory appears to predict that in the limit of small Γ the shift terms remain finite while injection diverges.

B. Floquet Formalism
We use Floquet theory to determine the nonperturbative effect of the electric field and couple the system to a bath that allows to sensibly describe steady state in the presence of relaxation processes.The microscopic Hamiltonian has the form ( ̵ h = e = 1, unless otherwise is stated): Expansion to 3rd order in vector potentials is necessary to compute electric currents correctly to order E 2 0 , since the current is the expectation value of the velocity operator: We assume that the electric field is a general periodic function with frequency ω: Since the Hamiltonian is periodic (ω = 2π T ), we apply discrete Fourier transform and now the Hamiltonian in Floquet picture has the following components: Thus the Hamiltonian structure in Floquet picture is: Additionally, one can show that in Floquet representation the velocity operator has the following form: We follow the approach of Ref. [12] and take a simplified 2 band model.We focus on inter-band resonant processes between the conduction and valence bands.The matrix elements highlighted by the rectangular box in Eq.( 46) are the ones that we keep within the "rotating wave" two band truncation, and ignore contributions of order E 3 and higher.Therefore, the Floquet Hamiltonian and velocity operator can be truncated to an effective 2 band model, which reads as: where 1 stands for "valence", 2 for "conduction ", (i, j) denotes the band index, E 1,2 are unperturbed band energies, 1,2 are dressed effective band energies and A = E ω.
The above formalism is a slight generalization of that in Ref. [12], that adds some dressing of the band energies by the drive, although this ingredient is not crucial for the key quantitative predictions made in the main text.Within our current notation, the following is the convention to convert a matrix from Floquet picture into the usual Schrödinger picture: And for the states (eigenvectors of the Hamiltonian from Eq.( 52)), the Floquet representation and the Schrödinger picture are related as:

C. General Keldysh Formalism
To model relaxation we couple our system to a simple Fermion bath at temperature T and chemical potential F .The bath couples uniformly to each site of the lattice sites where the fermions hop, as depicted in Fig. 1(b).The partition function and Lagrangian of the system are given by: where L TB describes the tight binding model of the system without the bath, L mix the tunneling from the system sites to the states in the bath, and L bath describes the bath.We will take ρ 0 to be the initial density matrix of the whole system plus bath, (a k , a † k ) are creation and annihilation operators of the bath states labeled by k, (c i , c † i ) is creation and annihilation operator of the fermion on the site "i", H ij is the tight-binding Hamiltonian of the system, bath k is the energy of the k-th state of the bath, β = 1 (k B T ), V k mix is the coupling of the system site (all sites of the physical system are coupled to identical and independent baths, as depicted in Fig. 1(b)) to the k-th state of the bath and f (ω) is the Fermi-Dirac distribution.C is the closed contour for the time integration of the Keldysh approach [42].
We can split the closed contour of time integration into two parts: the sub-script "+" or "-" denote the forward and backward parts of the contour and plays the role of an additional effective internal degree of freedom [38].From this one can write the action describing the system-bath coupling as follows: Now we apply the standard Keldysh fermionic rotation: After Fourier transforming the above we get: (71) Applying the same procedure to the rest of the action brings us to the following result: where G A R are advanced and retarded Green functions respectively and G K stands for the Keldysh Green function.They are given by: where After integrating out the fermionic bath, we obtain the following effective action for the system: In order to obtain the equal time lesser GF: we can focus on the imaginary parts of the Green functions: which leads us to the following form of effective interaction: where Γ(ω) is an effective imaginary self-energy or relaxation rate, that captures the relaxation processes and it is explicitly given by:

D. Keldysh Floquet Formalism
In the case of a periodically driven system the action satisfies: Therefore we can expand it with a the discrete Fourier transform: and arrive to the following action: Additionally, after splitting the frequency integration into segments: where ω = 2π T , and relabeling the summation index (n + m ′ = m), one obtains the result: And analogous procedure can be performed with the tight binding part of the action, yielding: where c i,a,n = c i,a (ω + nω).Now the system's full Green function, after integrating out the bath, is given by: where δ ij nm = δ ij δ nm are Kronecker deltas.The lesser Green Function G < is the one we need in order to obtain the density matrix of the system, which can related to advanced G A and retarded G R Greens functions, as follows: We have: Therefore the technical task is reduced to finding the advanced and retarded green functions of the system.

E. Truncated Green Functions
For the special choice of bath in Fig. 1(b) in which each site couples to an identical bath, the self energy from Eq.( 99) is independent of the system site indices "ij".To obtain the same truncation, in the spirit of the rotating wave approximation, as we had before for Eq.( 52), here we need to restrict to the indices (n = 0, n = −1).Thus effectively, after truncation we have: and all the notations are defined in Eq. (52)(53)(54)(55)(56)(57).
After inverting we have: The lesser GFs associated with "valence" and "conduction" bands are given by: Or explicitly: In general, the time average density matrix is given by: ρDC To proceed analytically, we assume a bath with a broad and flat density of states over the system's energy states, which allows us to neglect the frequency dependence of relaxation rate in Eq.( 84), and therefore we take Γ as a constant.Additionally, we assume that the dressed band energies are above or below Fermi energy, which allows us to take the occupation numbers as frequency independent (this is rigorously justified only in insulators, but in the case of metals, it can be viewed as a reasonable approximation that will captures the essence of the Pauli blocking effect of optical transitions).After frequency integration we have: The DC current is: where the velocity operator is: Thus the explicit result is given by: where ∑ k = ∫ d d k (2π) d and the lower (x, y, z) index denotes the decomposition of the matrix in σ-basis.The last term in the expression above is effectively an average of the velocity operator in the quasi-equilibrium state with dressed energies, and can be seen to identically vanish.
Now, we split the current into on three contributions: where which means that in order to compute the currents (J 1 , J 2 , J 3 ) with accuracy O(E 2 0 ), we can safely assume that vα F = ∂ Ĥ ∂k α and ignore the v α E 2 contribution.Finally, the velocity operator (with restored units) can be approximated as: Or in the Schrödinger picture: For future comparison let us write the density matrix derived from the Keldysh-Floquet approach (−iG < ) in the Schrödinger picture: where h ± = h x ± ih y .
F. The DC current in a clean limit In this section we will focus on the clean limit of the current Eqs.(117-119).If the electric field and Γ are small, we can use the following approximation: where Λ = E ⋅ A 12 2 + Γ 2 4 is the effective width of the Lorentzian and 12 = 1 − 2 .
By noticing that: we can write this as: where: Now we can write the components of the current (J 1 , J 2 , J 3 ) as: On resonance (ω ≈ 12 ) we can use the following approximation: The expression above combined with the identity: recovers the perturbation results Eqs.(31)(32)(33).The identification is the following: J 1 is the resonant shift current, J 2 is the non-resonant shift current, J 3 is the injection current, that in the clean limit (Γ = 0) are: . (135)

G. Periodic Gibbs Ensemble
In the spirit of the rotating-wave approximation, and in order to consider the same level of approximation in which Keldysh-Floquet formalism is developed, we take the evolution is defined by the following truncated Hamiltonian: We can parametrize Hamiltonian vector in spherical coordinates, namely h = h(cos ϕ sin θ, sin ϕ sin θ, cos θ).Using the Hamiltonian above one can find the evolution operator in the Schrödinger picture: where the first argument of the evolution operator is the final time and the second -the initial one: Now we follow the PGE procedure [39,[47][48][49].The system of interest has a conserved quantity that is the total number of particles.We construct the PGE density matrix in a way that it conserves the given quantity during the evolution.Initially occupations are given by: The time evolution of occupation numbers above is: Now we can can construct the PGE density matrix that satisfies the initial condition: Tr ρ S (0 The initial state of the system is chosen to be the thermal state (see the discussion in the main text) with the temperature and chemical potential of the bath: Solving Eq.( 146) one can obtain: Using Eq.(144) and Eq.( 148) one can rewrite the PGE density matrix from Eq.(145) as: Using the following relations: one can show that the PGE density matrix is: Which is the same matrix introduced in Eq.( 12) of the main text.

H. Rabi oscillations
We consider the 2 band system with the evolution operator from Eq.(137).We assume the system initially to be thermal with the temperature and chemical potential of the bath, namely: One can show that the evolution of the conduction and valence bands are given by: −i sin θ sin(ht)e iϕ (cos(ht) − i cos θ sin(ht))e iωt , (154) Now we can construct the Rabi density matrix as follows: To phenomenologically capture the synchronization of the system with drive, we perform a time average of the terms that have a frequencies different from the drive frequency ω, in the above the density matrix, as follows: Which leads to the following synchronized Rabi density matrix: which also can be simplified to: Which is the same matrix as that obtained from the PGE in Eq.( 152) and introduced in Eq.( 12) of the main text.

I. Injection current for 3D Weyl and 2D Dirac Fermions
Here we derive the approximate analytic expression of the injection current.The resonant injection current is given by: In the T → 0 limit, the valence/conduction band occupation difference is where Θ is the Heaviside theta function and F is the Fermi energy.

3D Weyl fermions
The Hamiltonian of 3D Weyl fermions is given by: Ĥ0 = α=x,y,z here we have re-scaled the momentum v 0 k → k, to simplify the final expression.The components of the momentum vector are k = kn = k(cos φ sin θ, sin φ sin θ, cos θ).
Consequently, the off-diagonal Berry connections and band energy difference of the Hamiltonian Eq.( 162) are given by: Note that A * 12 = A 21 .The direction of the injection current behaves as J 3 ∼ [E × E * ].We consider the frequency of the drive to be ω > 2 F , which sets f 1 − f 2 = 1.The injection current can be then approximated as: (165) The term containing the square root inside the integral is of the form (1 + X) −1 2 , with X < 1.We therefore expand this term as (1 + X) −1 2 ≈ 1 − X 2, and obtain: which after units restoring is: FIG. 3: Comparison of exact evaluation of the integral from Eq.( 165) and the approximate expression from Eq.( 167) for the injection current for 3D Weyl fermion as a function of a polarisation s, such that s → 0 for linearly polarizes light and s → 1 for perfectly circularly polarized light (see Eq.171) in the Rabi regime ev0 E ≫ Γ ̵ hω.
The above expansion can only be justified parametrically either when Γ ̵ hω ≫ ev 0 E or when the light is almost linearly polarized at arbitrary Γ, however, as we will show by explicit numerical evaluation of the integral, the approximation still works to about 17% even for perfectly circularly polarized light.
If the electric field is perfectly circularly polarised, the injection current of ideal Weyl model can be evaluated from the full integral in Eq.( 161) and yields the following expression: where hωΓ .This expression is in fact exactly identical to that derived in Ref. [30] in the limit of τ τ p → 0, where τ ,p are the energy and momentum relaxation times introduced in Ref. [30].The expression above in the limit of a large electric field E → ∞ can be approximated as: On the other hand the approximation from Eq.(167) gives: Despite the fact that we see a small discrepancy between approximation and exact calculation (J exact 3 J approx 3 ≈ 1.17), formula Eq.( 167) displays good agreement with the exact evaluation of the integral in Eq.( 165) for light that is not perfectly circularly polarized.To illustrate this we consider light with elliptical polarization that interpolates from perfectly linearly polarized to perfectly circularly polarized as follows: Parameters are chosen outside of the regime in which the approximation of Eq.( 167) is expected to be justified, namely . The result of comparison can be seen on FIG.I I 1.We see that Eq.( 167) matches the exact integral from Eq.( 165) for light that is almost linearly polarized (s ≈ 0) and deviates from it the most in the case of perfect circularly polarized light (s = 1), but only by about 17%.Therefore, we conclude that Eq.( 167) produces a good approximation of the photocurrent current over different regimes.

2D Dirac fermions
After the analogous momentum rescaling the Hamiltonian of 2D tilted Dirac fermions is given by: where α = u x v x .The momentum vector is two dimentional k = kn = k(cos φ, sin φ).The off-diagonal berry connections, band energies and velocity differences are: Assuming ω = 2 F , we see that . For simplicity we assume that the electric field is along the tilt so that the corresponding component of the injection current can be approximated as: which after units restoring is:

J. The shift current for linear in momentum models
This section focuses on contributions of shift current in models that are linear in momentum k.We demonstrate the contribution of shift currents is neglible in such models and the rectification response is dominated by the the injection current in the limit where the two band approximation is valid.

Two-band Keldysh-Floquet formalism
In the Keldysh-Floquet formalism employed in the main text, the velocity operator from Eq.(120) for the linear in momentum model simplifies to the following expression: which is a time-independent operator in Schrödinger's picture.Therefore, we see that there are no associated off-diagonal components of the velocity operator, which are the ones that would give rise to shift currents.After averaging of the density matrix from Eq.( 4) of the main text with the velocity operator written above, the resultant expression of the current picks only the contribution from the injection current: which means that in the two-band approximation of the linear in momentum model k the shift current is absent.
2. 3D Weyl model's shift current from the perturbation theory Now, let's consider the prediction of the perturbation theory for the 3D Weyl fermion model, which Hamiltonian is: According to the perturbation theory, the following expression gives the shift current part of the rectification conductivity: As one can see, the shift current contribution of 3D ideal Weyl fermion is of order Γ and can be neglected in a clean limit.
Alternatively, one could use the time-reversal and rotational symmetry argument to show that the perturbative shift current vanishes for the Hamiltonian of Eq.(180) in the clean limit (Γ → 0).First, the shift conductivity from Eq.(181) of Hamiltonian Eq.(180) has to be a three index tensor, that is symmetric under SO(3) transformations.π rotational symmetry around x, y, z axes limits only those components of the tensor Eq.( 181) to be finite, which all three indices are distinct.Additionally, the tensor must be symmetric under cyclic permutations of (x, y, z) labels, as these can be implemented as a subgroup of SO(3).These properties force the tensor to be proportional to the Levi-Civita tensor, and we conclude that the shift conductivity from Eq.(181) of Hamiltonian Eq.(180) has to be proportional to it.On the other hand, it can be shown that in the clean limit the resonant part of the shift current is even under time reversal symmetry (see table 1 of Ref. [50]).Nevertheless the part of bilinear of electric fields that contracts with the Levi-Civita tensor is its circularly polarized component, which is a timereversal-odd, and thus such components must be absent from the perturbative expressions of the shift current in the clean limit [50].In summary, the combination of three ingredients: time-reversal symmetry, SO(3) symmetry and the clean limit (Γ → 0), force the shift current of the ideal Weyl models to vanish.

2D Dirac model's shift current from the perturbation theory
Next, we focus on another linear momentum modelthe 2D Dirac model.The Hamiltonian of this model is: According to the perturbation theory, the resonant part of the shift conductivity is: ) We see that the resonant shift currents do not vanish exactly for linear in k models of 2D Dirac fermions from perturbation theory.However, the injection conductivity is typically parametrically larger than the shift conductivity in the limit in which one is justified to focus only on the two bands (Γ 2 ≪ 2 F − m 2 ≪ 2 F ), which is the same regime in which our rotating-wave two-band truncation of the Floquet-Keldysh formalism is justified.In this limit the dominant contribution of the shift conductivity typically scales as σ S,R ∼ 2 F − m 2 3  F , while the injection conductivity typically scales as σ I,R ∼ 2 F − m 2 ( 2 F Γ), and the ratio of them is σ I,R σ S,R ∼ F Γ. Consequently, the shift current contribution can be neglected.

FIG. 1
FIG. 1: a) Energy crossing between boosted valence and conduction bands in Floquet representation.b) Depiction of underlying tight binding model with physical sites (red balls) which are tunnel coupled (solid lines) among themselves and with their own identical fermionic bath (blue balls).

2 F −m 2 4 F im 2 F 2 F −m 2 4 F − 2 F −m 2 3 F 2 F 2 F −m 2 (2m 2 + 2 F ) 3 4 F im 2 F −m 2 3 F − im 2 F −m 2 3 F 2 F −m 2 (
ω ↔ −ω .(189)It turns out that in this model the shift current is finite.The value of the tilt, u x v x , determines the frequency window of non-zero resonant shift current, and the maximum value of the conductivity for small tilts occurs near the middle of this window at a frequency ω ≈ 2 F (see Fig.2(c) of the main text).The peak values of the conductivity can be analytically computed and are given by (here ̵ h = e = v x = v y = 1):σ xβα S,R (−ω, ω)