Non-equilibrium nature of non-linear optical response: Application to the bulk photo voltaic effect

TL;DR

This paper investigates the non-linear optical bulk photovoltaic effect, emphasizing its non-equilibrium nature and the role of electron-phonon interactions, revealing non-local responses and the limitations of equilibrium theories.

Contribution

It introduces a semiclassical approach and quantum master equation analysis to understand the non-linear response, highlighting the importance of out-of-time-ordered correlators.

Findings

Non-linear optical response involves out-of-time-ordered correlators.

Scattering rate significantly influences non-linear DC response.

Non-local effects are crucial for spatially inhomogeneous excitation profiles.

Abstract

The bulk photovoltaic effect is an example of a non-linear optical response that leads to a DC current that is relevant for photo-voltaic applications. In this work, we theoretically study this effect in the presence of electron-phonon interactions. Using the response function formalism we find that the non-linear optical response, in general, contains three operator correlation functions, one of which is not ordered in time. This latter correlator cannot be computed from equilibrium field theory. Using a semiclassical approach instead, we show that the bulk photovoltaic effect can be attributed to the dipole moment of the generated excitons. We then confirm the validity of the semiclassical result (which agrees with the non-interacting result) for non-linear DC response from a quantum master equation approach. From this formalism we find that, in contrast to usual linear response, the scattering rate has a strong implicit effect on the non-linear DC response. Most interestingly, the semiclassical treatment shows that the non-linear DC response for spatially inhomogeneous excitation profiles is strongly non-local and must involve the aforementioned out-of-time-ordered correlators that cannot be computed by equilibrium field theory.

Equations

$$\hat{H}(t)={\hat{H}}_0+\int dxA(r)\hat{J}(r)$$

$$\begin{array}{rl}& ⟨\hat{J}(r)⟩\simeq ⟨\hat{J}(r){⟩}_0+{\int }_{-\infty }^{\infty }d{x}^{\prime }{A}_0(x^{\prime }){\chi }_{lin}(r,x^{\prime })\\ & +{\int }_{-\infty }^{\infty }d{x}^{\prime }{A}_0^2(x^{\prime }){\chi }_{n,2}(r,x^{\prime })\\ & +{\int }_{-\infty }^{\infty }d{x}_{1}d{x}_2{A}_0(x_1)A_0(x_2){\chi }_{n,3}(r,x_1,x_2)\\ & +{\int }_{-\infty }^{\infty }d{x}_{1}d{x}_2{A}_0(x_1)A_0(x_2){\chi }_{otoc}(r,x_1,x_2)\end{array}$$

$${\chi }_{lin}(r,x^{\prime })=\iota {\int }_{-\infty }^{t}d{t}^{\prime }{e}^{\eta t^{\prime }}\alpha (t^{\prime }){⟨\text{commutator}{\mathcal{J}}_0(r^{\prime }){\mathcal{J}}_0(r)⟩}_0$$

$${\chi }_{n,2}(r,x^{\prime })=\iota {\int }_{-\infty }^{t}d{t}^{\prime }{e}^{2\eta t^{\prime }}{\alpha }^2(t^{\prime }){⟨\text{commutator}{\mathcal{J}}_1(r^{\prime }){\mathcal{J}}_0(r)⟩}_0$$

$$\begin{array}{rl}& {\chi }_{n,3}(r,x_1,x_2)\\ & =-\frac{1}{2}{\int }_{-\infty }^{t}d{t}_{1}d{t}_2{e}^{\eta (t_1+t_2)}\alpha (t_1)\alpha (t_2){⟨\text{anticommutator}\mathcal{T}\{{\mathcal{J}}_0(r_1){\mathcal{J}}_0(r_2)\}{\mathcal{J}}_0(r)⟩}_0\end{array}$$

$$\begin{array}{rl}& {\chi }_{otoc}(r,x_1,x_2)={\int }_{-\infty }^{t}d{t}_{1}d{t}_2{e}^{\eta (t_1+t_2)}\alpha (t_1)\alpha (t_2){⟨{\tilde{\mathcal{J}}}_0(r_1){\tilde{\mathcal{J}}}_0(r){\tilde{\mathcal{J}}}_0(r_2)⟩}_0\end{array}$$

$${\hat{H}}_0=\sum _k{\Delta }_k\left(c_k^{†}{c}_k-v_k^{†}{v}_k\right)$$

$$\begin{array}{rl}& |\delta {\Psi }_k(t){⟩}_I\simeq A_0(J_{0,k}^x-\iota J_{0,k}^y)\frac{\sin \left(\frac{\Omega }{2}-{\Delta }_k\right)t}{\left(\frac{\Omega }{2}-{\Delta }_k\right)}|{\overline{eh}}_k{⟩}_I\\ & \equiv A_0|J_{0,k}^{\perp }|e^{\iota {\varphi }_k}\frac{\sin \left(\frac{\Omega }{2}-{\Delta }_k\right)t}{\left(\frac{\Omega }{2}-{\Delta }_k\right)}|{\overline{eh}}_k{⟩}_I\end{array}$$

$$\begin{array}{rl}& P(t)\equiv \sum _k⟨\delta {\Psi }_k(t)|\left(r_e-r_h\right)|\delta {\Psi }_k(t)⟩\\ & =\frac{A_0^2}{2\pi }{\int }_{BZ}dk\frac{{\sin }^2\left(\frac{\Omega }{2}-{\Delta }_k\right)t}{{\left(\frac{\Omega }{2}-{\Delta }_k\right)}^2}{\left({\vec{J}}_{0,k}\text{crossproduct}{\vec{J}}_{1,k}\right)}_z\end{array}$$

$$f(t)\equiv \eta e^{-\eta t}$$

$$\begin{array}{rl}& ⟨P⟩\equiv {\int }_0^{\infty }dtP(t)f(t)=\frac{A_0^2}{\pi }{\int }_{BZ}dk\frac{{\left({\vec{J}}_{0,k}\text{crossproduct}{\vec{J}}_{1,k}\right)}_z}{{\eta }^2+{\left(\Omega -2{\Delta }_k\right)}^2}.\end{array}$$

$$\begin{array}{rl}& J_{shift}^{DC}=\frac{⟨P⟩}{\tau }=\frac{A_0^2}{\pi }{\int }_{BZ}dk\frac{\eta }{(\Omega -2{\Delta }_k{)}^2+{\eta }^2}({\vec{J}}_{0,k}\text{crossproduct}{\vec{J}}_{1,k}{)}_z.\end{array}$$

$$\underset{\eta \to 0}{\lim }{J}_{shift}^{DC}=A_0^2{\int }_{BZ}dk({\vec{J}}_{0,k}\text{crossproduct}{\vec{J}}_{1,k}{)}_z\delta (\Omega -2{\Delta }_k)$$

$${\partial }_t\rho (x,t)=-{\partial }_x(J^d+J_{shift}^{DC})$$

$${\partial }_t\rho =D{\partial }_x^2\rho +f(x)$$

$$f(x)=-{\partial }_x{J}_{shift}^{DC}=2\alpha A_0(x){\partial }_x{A}_0(x)$$

$$\begin{array}{rl}& \rho (x,t)={\int }_{-\infty }^{\infty }d{x}^{\prime }f(x^{\prime }){\int }_0^{t}d\tau \frac{e^{-\frac{(x-x^{\prime }{)}^2}{4D\tau }}}{\sqrt{4\pi D\tau }}\end{array}$$

$$J^d(x)={\int }_{-\infty }^{\infty }d{x}^{\prime }f(x^{\prime })\frac{(x-x^{\prime })}{4\sqrt{\pi D}}{\int }_0^{t}d\tau {\tau }^{-\frac{3}{2}}{e}^{-\frac{(x-x^{\prime }{)}^2}{4D\tau }}$$

$$J^d(x)=\frac{1}{2}{\int }_{-\infty }^{\infty }d{x}^{\prime }f(x^{\prime })sgn(x-x^{\prime })$$

$$J_{shift}^{DC}=-J^d=D{\partial }_x\rho$$

$$J^d(x)=\alpha {\int }_{-\infty }^{\infty }d{x}^{\prime }d{x}^{\prime \prime }{A}_0(x^{\prime })A_0(x^{\prime \prime })sgn(x-x^{\prime }){\partial }_{x^{\prime }}\delta (x^{\prime }-x^{\prime \prime })$$

$${\chi }_{otoc}(r,x_1,x_2)\simeq {\partial }_{x_1}\delta (x_1-x_2)\frac{(x-x_1)}{4\sqrt{\pi D}}{\int }_0^{t}d\tau {\tau }^{-\frac{3}{2}}{e}^{-\frac{(x-x_1{)}^2}{4D\tau }}$$

$$J^d=\frac{n{e}^2{\tau }_d}{m}\frac{V}{L}\propto \frac{A_{0}V}{L}$$

$$P_{max,BPVE}=(J_{shift}^{DC}{)}^{2}R\propto \frac{A_0^{3}L}{{\tau }_d}.$$

$${\hat{H}}_F=\sum _k\left({ϵ}_k^u{f}_k^{u†}{f}_k^u+{ϵ}_k^d{f}_k^{d†}{f}_k^d\right)$$

$${\hat{H}}_p=\sum _{q,s}{\omega }_q^s{a}_q^{s†}{a}_q^s$$

$$\begin{array}{rl}& {\hat{U}}_0\equiv {\hat{U}}_F\otimes {\hat{U}}_p=e^{-\iota {\hat{H}}_{F}t}\otimes e^{-\iota {\hat{H}}_{p}t}.\end{array}$$

$$\hat{V}=\sum _{q=k-g}^{k,g,s}\left(V_1{c}_k^{†}{c}_g+V_2{v}_k^{†}{v}_g\right)\left(a_q^s+a_{-q}^{s†}\right)$$

$${\partial }_t{\hat{\rho }}_I=-\iota \text{commutator}{\hat{V}}_I{\hat{\rho }}_I.$$

$$\begin{array}{r}{\partial }_t{\hat{\rho }}_I^e=-Tr{\left\{\text{commutator}{\hat{V}}_I(t){\int }_{-\infty }^{t}d{t}^{\prime }\text{commutator}{\hat{V}}_I(t^{\prime }){\hat{\rho }}_I(t)\right\}}_p.\end{array}$$

Full text

I Introduction

The fluctuation-dissipation theorem that relates the perturbative response of observables to equilibrium correlators of observables is one of the corner stones of quantum many-body physics. Recently, a number of non-linear phenomena such as non-linear optical response, photogalvanic effect, pump-probe spectroscopy have been used to characterize solid state materials. For example, certain second harmonic generation coefficients have been used to detect a hidden inversion symmetry breaking nematic phase Harter2017 . Similar correspondence has been shown between second harmonic generation and inversion symmetry breaking in certain topological materials with hexagonal warping term and a gap for quasiparticle excitations Li2019 ; Nagaosa2017 . The non-linear response of the DC current in a Weyl material to circularly polarized light, the so-called circular photogalvanic effect has been shown to be quantized Moore2017 . In fact, the quantized circular photogalvanic response has been directly related to the topological monopole charge of Weyl materials Moore2017 .

The circular photogalvanic effect is one example of a broader class of bulk photovoltaic effects where the non-linear response to optical excitation results in a DC current. The bulk or anomalous photovoltaic effect, which was discovered experimentally nature1946 as a large voltage in optically illuminated semiconductors has been attributed to a large number potential sources such as asymmetric scattering from impurities, phonons or electron-electron interactions Struman1977 ; Belinicher1978 ; Strumanrev1980 ; Struman'sbook ; Fridkin'sbook ; Spivak2009 . The bulk photovoltaic effect (BPVE) for linearly polarized light has been argued to be a result of shift currents Sipe2000 in bulk materials, where the current is generated by the sequential tunneling of electrons to one direction. Theoretically, using a non-linear response formalism, the BPVE was argued to arise intrinsically even in single crystalline materials that break inversion symmetry Kraut1981 . It was found that the magnitude of this response is directly related to the Berry curvature over parts of the Brillouin zone where the optical perturbation is resonant Morimoto2016 ; Kraut1981 . Aside from fundamental interest as a probe of Berry curvature, the BPVE also has potential applications to solar cell devices provided materials with high efficiency can be found Cook2017 ; Tan2016 ; Rangel2017 ; Rappe2012 ; Zheng2016 ; Nakamura2017 ; Wang2017 ; Gong2018 ; Ogawa2017 .

Despite the use of characterizing symmetries and topology of materials, the understanding of non-linear optical responses such as the BPVE is quite limited outside the purview of non-interacting systems. As elaborated in Sec. II, the response function formalism for non-linear effects such as the BPVE contains terms that are neither time-ordered nor anti-time ordered. Such out of time ordered correlated terms that have also been shown to appear in non-linear response of superconductors Larkin1968 cannot be computed using standard equilibrium Feynman diagram formalism and in this sense are truly non-equilibrium properties.

The utility of the BPVE for solar cell devices ultimately relies on the ability to annihilate photons and convert them to electrical power. However, in the equilibrium diagram formalism, each power of the vector-potential $A(t)$ is associated with either the absorption or emission of a photon. The creation of a DC current at second order in a vector potential would then have to be interpreted as the sequential absorption and emission of a photon. Within this picture, the generation of electrical energy appears quite paradoxical. An electrical current can only generate electrical energy in the presence of a finite voltage that can arise in the presence of a scattering process, which is ignored in the non-linear response theory for shift currents Sipe2000 ; Kraut1981 ; Morimoto2016 . While scattering processes play a crucial role in certain extrinsic mechanisms for the BPVE Belinicher1978 ; Struman1977 , the non-linear response mechanism is typically considered to ignore such scattering processes Sipe2000 ; Kraut1981 ; Morimoto2016 . Thus, it is desirable to extend the non-linear response formalism for the BVPE to include such scattering processes in a way that is completely consistent with energy conservation.

In this work, we systematically study the non-linear response approach to the BVPE in a minimally interacting system, where we include the interaction of the electrons with the phonons in the material. We assume that the electron-phonon interaction is symmetric so as to not directly contribute to the shift current. In Sec. II, we start by developing a response function formalism for the BPVE and demonstrating the role of out of time-ordered correlators (OTOC) Larkin1968 . However, the OTOC relevant to this work involves three operators and is not connected to Lyapunov exponents as is the case with the OTOC in Ref. Larkin1968 . The presence of such correlators invalidates the standard computation using equilibrium Feynman diagrams. In order to gain intuition, we first develop a semi-classical picture of the BPVE (Sec. III) based on the processes shown in Fig. 1. This picture is shown to lead to the conventional expression for the DC current in the case of homogeneous optical excitation. However, the application of inhomogeneous excitation profiles are found to lead to results that cannot be described by equilibrium response functions. Such a semiclassical treatment is only justified in the strong electron-phonon scattering limit relative to the optical excitation. We go beyond this limit in Sec. IV using a Floquet Master equation approach. We find that this approach, in addition to reproducing the known small vector potential limit of the DC response also predicts a linearly scaling DC response in the limit of small electron-phonon scattering. Finally, in Sec. V, we study the cross-over between the small and large vector potential limit by computing the response numerically within the same formalism for the Rice-Mele model Rice1982 .

II Second order optical response

We consider a generic Hamiltonian for a two band semiconductor coupled to electromagnetic (EM) field, given by

[TABLE]

where $\hat{H}_{0}$ is the two band system Hamiltonian and $r=(x,t)$ includes both space and time co-ordinates. $\hat{J}(r)$ is the current operator which couples to the EM vector potential, $A(r)=A_{0}(x)\,e^{\eta t}\left(e^{\iota\Omega t}+e^{-\iota\Omega t}\right)\equiv A_{0}(x)\,e^{\eta t}\alpha(t)$ where $\eta=0_{+}$ is the rate at which the EM field is turned on.

We can expand $\hat{J}(r)$ in powers of $A$ as: $\hat{J}(r)\simeq\hat{J}_{0}(r)+A(r)\hat{J}_{1}(r)$. The response $\langle\hat{J}(r)\rangle$, up to second order in $A_{0}(x)$, can be written as (see Appendix A)

[TABLE]

where $\chi_{lin}$ is linear response given by

[TABLE]

Here the current operators are in interaction picture defined as $\mathscr{J}_{0,1}(r)\equiv e^{\iota\hat{H}_{0}t}\hat{J}_{0,1}(r)e^{-\iota\hat{H}_{0}t}$. Among the non-linear response functions, $\chi_{n,2}$ is the part of non-linear response that can be viewed as a linear response of the current to a perturbation $A^{2}_{0}(r)\,\hat{J}_{1}(r)$ and is given by:

[TABLE]

The term proportional to $\chi_{n,3}$, while non-linear and involving three current operators, is given by the time-ordered product as

[TABLE]

$\chi_{otoc}$ involves out of time ordered correlator of cubic current operators given by

[TABLE]

In these expressions, $\langle..\rangle_{0}$ means expectation value w.r.t. the ground state of the system. Since $\chi_{otoc}(t)$ contains out of time ordered correlator (OTOC), Feynman diagrams cannot be used to evaluate it. One requires Keldysh formalism Larkin1968 to compute such correlators in presence of any interaction.

Another way to see that the zero frequency part of $\langle\hat{J}(r)\rangle$ is likely not computable as a response around equilibrium is to quantize the electromagnetic field in terms of photons (as mentioned in the introduction). Expanding $A(t)=a^{\dagger}e^{i\omega t}+h.c$ and substituting, A(t) into Eq. 2, we see that $\langle\hat{J}(r)\rangle_{DC}\propto\langle aa^{\dagger}\rangle$ , which represents absorption and remission of a photon. This suggests that no energy is absorbed, even though the DC current in the presence of a voltage would do work. This would likely represent a non-equilibrium situation suggesting the need for a non-equilibrium theory. Consequences of this somewhat imprecise argument will become clear from the semiclassical treatment of the DC current in the next section.

III SemiClassical description

In order to understand the absorption of energy more clearly, we include an explicit dissipation process through phonon scattering with a semiclassical approach. In this section we treat electron-phonon scatterings as probabilistic events along with the quantum mechanical evolution of the electronic state due to external EM field. For simplicity we consider two band Hamiltonian for a semiconductor given by

[TABLE]

where $c_{k}$($v_{k}$) is conduction (valence) band state operator and $2\Delta_{k}$ is the band gap.

III.1 Shift current in homogeneous excitation

The incident EM radiation leads to the creation of electron-hole (e-h) pairs. The time dependent amplitude of such a pair can be extracted from the EM coupling in Eq. 1 to be

[TABLE]

where $|0\rangle$ is the ground state and $|\bar{eh}_{k}\rangle\equiv c^{\dagger}_{k}v_{k}|0\rangle$ is e-h pair state. $\vec{J}_{0}\cdot\vec{\sigma}\equiv\hat{J}_{0}=\partial_{A}\hat{H(t)}|_{A=0}$ such that $\vec{\sigma}$ is the vector of Pauli matrices in the band basis (i.e. $\sigma^{z}|c_{k}\rangle=|c_{k}\rangle$, $\sigma^{z}|v_{k}\rangle=-|v_{k}\rangle$) and $|J^{\perp}_{0,k}|\equiv\sqrt{(J^{x}_{0,k})^{2}+(J^{y}_{0,k})^{2}}$.

In non-centrosymmetric materials, the resulting e-h pair can carry a dipole moment. So the electric field generated dipole moment evolves in time accordingly as:

[TABLE]

where $\vec{J}_{1,k}\equiv-\partial_{k}\vec{J}_{0,k}$. In this semiclassical approach the coherent evolution under the electric field is interrupted by phonon scatterings modeled as Poisson processes. The coherent evolution time $t$ corresponding to the Poisson process is described by an exponential distribution function as

[TABLE]

where $\eta$ is the total rate for electron and hole scatterings as detailed in the later part of this subsection. In the weak electric field limit, once the scattering happens the e-h pair disperses equally in either direction and no further dipole moment can be generated. Thus integrating the quantum dipole moment in Eq. 9 with $f(t)$ leads to an average dipole moment generated due to e-h pair, given by

[TABLE]

The DC current can be determined by the average rate of dipole moment generation over the characteristic scattering time, $\tau\equiv\frac{1}{\eta}$. Using this relation we obtain the expression for the shift current to be

[TABLE]

Here $\eta$ is a physical energy scale corresponding to electron phonon scattering strength. Considering the limit of $\lim_{\eta\to 0}\frac{\eta}{(\Omega-2\Delta_{k})^{2}+\eta^{2}}=\pi\delta(\Omega-2\Delta_{k})$ in Eq. 12 we obtain

[TABLE]

which matches with the known expression for shift current given in Moore2017 (see Appendix B).

The semiclassical approach presented above, is valid in the limit of strong e-ph coupling or weak electric field ($A_{0}|J^{\perp}_{0,k}|\ll\eta$) where the scattering processes are instantaneous events with time scale $\tau\ll\frac{1}{A_{0}|J^{\perp}_{0,k}|}$. In Sec. IV we will use a more quantum mechanical approach to study the case of longer scattering times $\tau$.

Let us now consider in detail the scattering of electrons by phonons following the coherent exciton generation discussed above. Assuming the phonon scattering matrix elements are momentum independent, the momenta and the corresponding velocities of the electrons and holes are uncorrelated. The resulting motion of the electrons and holes are diffusive starting from the average position. The electrons and holes can recombine and emit a phonon. A current is generated if the electron and hole annihilate with a pair that is different from the way they were created. This is the notion of a shift current. However, if the electron and hole of the generated pair recombine with each other before they can diffuse away, for local electron-hole recombination, the recombined excitons do not have any dipole moment. Therefore, in calculating the shift current, the average dipole moment of the exciton, remains unchanged right until the exciton merges into the electron-hole gas.

III.2 Diffusion current response for inhomogeneous excitations

Let us now consider the response of current to an inhomogeneous optical excitation profile applied to the system which is often the case in practical setups Struman'sbook . We assume that the intensity of the optical excitation is non-zero over a finite region $R$ in Fig. 2 and vanishes outside. We note that at finite temperature, there is finite (although exponentially small as temperature vanishes) carrier density due to thermal excitations on top of EM excited electron and holes. Recombination with these thermally excited carriers leads to a finite lifetime for both electrons and holes, even though the total charge density is conserved. The shift current $J_{shift}^{DC}=\alpha A^{2}_{0}(x)$, in the case of an inhomogeneous intensity $A^{2}_{0}(x)$, varies in space. This combined with the conservation of charge leads to an inhomogeneous charge density. The gradient in charge leads to a diffusion current $J^{d}=-D\partial_{x}\rho$, where the diffusivity of charge $D=\sigma/\chi$ and $\sigma$ and $\chi$ are the conductivity and compressibility of the unperturbed carriers. The density profile is determined from the equation of charge conservation as

[TABLE]

Substituting the form for the diffusion current transforms this equation into a diffusion equation given by

[TABLE]

where the source term is given by the gradient of shift current and has the form:

[TABLE]

Since there is no charge density initially (i.e. $\rho(x,0)=0$), the solution of Eq. 15 is given by

[TABLE]

The diffusion current can be obtained by differentiating the solution in Eq. 17 w.r.t. $x$ to be

[TABLE]

Performing a change of variable, $p=\frac{|x-x^{\prime}|}{2\sqrt{D\tau}}$, the limits for $p$ becomes $\infty$ and $\frac{|x-x^{\prime}|}{2\sqrt{Dt}}$ respectively. In the long time limit set by the illumination length and diffusion coefficient($t\gg\frac{|x-x^{\prime}|^{2}}{D}$), we can approximate the lower limit to be zero $\frac{|x-x^{\prime}|}{2\sqrt{Dt}}\to 0$. In this long time limit, performing the Gaussian integral we obtain the expression for the diffusion current to be

[TABLE]

Substituting Eq. 16 in Eq. 19 and integrating by parts we obtain $J^{d}(x)=-J^{DC}_{shift}(x)$, i.e. the diffusion current cancels the shift current in steady state. A similar situation is discussed in Ref. Strumanrev1980 where the photocurrent is compensated by the ohmic current in a circuit with open boundary condition. The cancellation of the shift current by the diffusion current relates the shift current to the density profile as:

[TABLE]

The charge density $\rho(x)$ can be measured capacitively using the gate shown in Fig. 2. This provides a way to measure the shift current without contacts Yakovenko .

Using $A_{0}(x^{\prime})=\int dx^{\prime\prime}\,A_{0}(x^{\prime\prime})\delta(x^{\prime}-x^{\prime\prime})$ we can rewrite the expression for the diffusion current (Eq. 19) as

[TABLE]

From this expression, it is clear that the diffusion current in case of inhomogeneous excitation profile (i.e. $A_{0}(x)$ not constant in space) is non-local. This implies that the resultant DC current from the response, given by Eq. 2, must contain a spatially non-local contribution.

The non-local response is quite unexpected in the light of Sec. II. This is because the time-ordered response functions $\chi_{n,2}$ and $\chi_{n,3}$ can be related by analytic continuation to imaginary time correlation functions, which in turn can be written as derivatives of an effective action Sachdev . This is similar to finite temperature correlation functions in classical statistical mechanics, where correlation functions decay on the scale of a correlation length Goldenfeld . One subtlety to keep in mind here is that $\chi_{n,2}$ and $\chi_{n,3}$ in Eqs. 4-5 also involve time. However, analytic continuation from imaginary time suggests that the time-ordered correlations at frequencies larger than temperature are also expected to be local in space and time on a scale determined by correlation lengths and times Sachdev . Thus one does not expect $\chi_{n,2}$ or $\chi_{n,3}$ to contribute to the non-local response in Eq. 21.

Assuming that $\chi_{n,2}$ and $\chi_{n,3}$ are local for the reasons outlined in the previous paragraph, the only term in Eq. 2 that can incorporate non-locality of $J^{d}$ is the term involving OTOC, as appears in the expression of $\chi_{otoc}$ in Eq. 6. The total current is written as $\langle J\rangle_{DC}=J^{DC}_{shift}+J^{d}$. Since the shift current response $J^{DC}_{shift}$ is local, we can compare the non-local contribution contained in $\chi_{otoc}$ with the non-local diffusion response (Eq. 18) to write the non-local part of $\chi_{otoc}$ as

[TABLE]

for $|x-x_{1}|,|x-x_{2}|\gg\xi_{th}$, where $\xi_{th}$ is a thermal correlation length. Thus, the semiclassical approach can be used to constrain the non-linear response correlator (Eq. 6), which cannot be computed from equilibrium field theory methods. In principle, this response can be computed from a Keldysh field theory approach Larkin1968 , but it is clear that these correlators will be qualitatively more non-local compared to equilibrium correlators.

III.3 Finite voltage and size effect on BPVE

Finite voltage that may be generated in the presence of spatial variations can have an effect on the electron-hole gas as mentioned at the end of Sec. III.1. The carrier density of this electron hole plasma, $n$, is limited since the rate of exciton generation, for zero temperature, is balanced by recombination rate. We emphasise here that the time scale for recombination processes to lead to a steady state is inversely proportional to some power of EM amplitude ($\propto\frac{1}{A_{0}^{\alpha}},\;\alpha>0$) till which the response function formalism is not applicable. We circumvented this problem in the last sub-section by working at finite temperature. The analysis of the balancing of the recombination process at zero temperature considered in this sub-section, assuming $\alpha=1$ for free diffusion of carriers, leads to a steady state carrier density with $n^{2}\propto A_{0}^{2}$. The shift current in the absence of a voltage can only be stable for perfect translation invariance with periodic boundary conditions. Any variations of parameters in the system that would change the shift current, leads to a finite voltage. Let us now consider boundary conditions different from periodic so that a finite voltage $V$ can be generated over length $L$ of the material. The voltage $V$ will lead to a combination of drift and diffusion current among the carriers which is given by

[TABLE]

where $\tau_{d}$ is the scattering time for the charge carriers.

The combined system with a BPVE together with the drift current can be thought of as a current source $J^{d}$ in parallel with a resistance $R\propto L/(\tau_{d}A_{0})$. The maximum voltage generated by this process is $V_{max}=J^{DC}_{shift}R\propto LA_{0}/\tau_{d}$. This voltage can become anomalously large in the limit of large L which is consistent with experimental measurements reported in Ref. Fridkin2001 and the discussions in Ref. Strumanrev1980 . The maximum power that can be extracted from this system is

[TABLE]

Note that total incident power scales as $P_{inc}\propto A_{0}^{2}L$. Thus, despite the high open circuit voltage $V_{max}$, the efficiency scales as $P_{max,BVPE}/P_{inc}\propto A_{0}$ which is quite small and in agreement to the experimental results as discussed in Ref. Fridkin2001 .

IV Quantum approach

In this section we go beyond the limit of weak optical excitation strength by using a more quantum mechanical master equation based approach to the photocurrent. The problem of the semiconductor in EM field can be described by Eq. 1 where $\hat{H}_{0}$ is the Bloch Hamiltonian as in Eq. 7. The EM coupled Hamiltonian (Eq. 1) is time periodic with period $T\equiv\frac{1}{\Omega}$. The time periodic Hamiltonian in Eq. 1 can be expressed in Floquet basis as

[TABLE]

where $\epsilon^{u,d}_{k}\in[0,\Omega]$ are the two quasi energies and $f^{u,d}_{k}$ correspond to Floquet state operators (see Appendix C for derivation of Floquet spectrum).

In the quantum mechanical treatment, we include the phonons explicitly with a Hamiltonian that is written as

[TABLE]

where $a^{s}_{q}$ is annihilation operator for phonon of mode $s$ and momentum $q$. In absence of any e-ph coupling the evolution operator for combined e-ph system is given by

[TABLE]

The phonons couple to the electrons through the e-ph interaction term in the Hamiltonian given by

[TABLE]

where $V_{1}$($V_{2}$) is the amplitude of intra band scatterings within conduction (valence) band.

The equation of motion for e-ph combined density operator $\hat{\rho_{I}}(t)$ , in the interaction picture basis that is rotating according to Eq. 27, is given by

[TABLE]

Note that $\hat{V}_{I}$ is obtained from Eq. 28 by transforming the field operators $c_{k},v_{k}$ and $a^{s}_{q}$ to the interaction picture basis. In a standard procedure (outlined in Appendix D) to obtain the dynamical equation for the reduced electronic density operator, $\hat{\rho}^{e}_{I}$, Eq. 29 is solved up to second order in $|V_{1,2}|$ following which the phonon operators are traced out under Markovian approximation. The resultant equation of motion for $\hat{\rho}^{e}_{I}$ becomes

[TABLE]

$\hat{\rho}^{e}_{I}$ can be represented more concretely in terms of its matrix elements in Floquet basis:

[TABLE]

where in the last line we have used the cyclic property of trace.

IV.1 Equation for equilibrium distribution

To obtain an equation of motion for $\Pi^{ab}_{k}$ we differentiate Eq. 31 and use Eq. 30 in the similar way described in Tom2015 . In the process we assume zero temperature for the phonon bath and translational invariance for the electronic system (see Appendix D) to obtain the final form of the equation of motion for $\Pi^{ab}_{k}$ to be

[TABLE]

where $F_{k}$ can be thought of re-normalized Hamiltonian (local in $k$) and is given by,

[TABLE]

Here the Fourier components of the matrix elements of $V$ and $\tilde{V}$ are given by,

[TABLE]

where $|\psi^{\alpha}_{k}(n)\rangle$ is the nth frequency component of the Floquet state $|\psi^{\alpha}_{k}(t)\rangle$ and $\hat{V}^{S}=\sum_{p,p^{\prime}}\left(V_{1}\,c^{\dagger}_{p}c_{p^{\prime}}+V_{2}\,v^{\dagger}_{p}v_{p^{\prime}}\right)$. $G$ is the constant density of states for the phonon bath and $\Theta(x)$ is the Heaviside step function.

The structure of the tensor $D_{k}$,

[TABLE]

in Eq. 32, can be understood by considering the diagonal limit of $\Pi_{k}$. In this limit the product $\Pi^{\beta\beta}_{g}(1-\Pi^{\alpha\alpha}_{k})$ along with the corresponding matrix elements represent the probability of scattering from $|\psi^{\beta}_{g}\rangle$ to $|\psi^{\alpha}_{k}\rangle$ state. The $\Theta$ function ensures the initial state has higher quasi energy than the final state thus only allowing processes that dissipate energy to the phonon bath. Similarly the tensor $A_{k}$ is involved in reverse scattering processes (for instance from $|\psi^{\alpha}_{k}\rangle$ to $|\psi^{\beta}_{g}\rangle$ scatterings) and is given by,

[TABLE]

The steady state version of Eq. 32 can be obtained by setting $\partial_{t}\Pi^{ab}_{k}=0$. The result is a set of non linear algebraic equations involving $\Pi^{\alpha\beta}_{k}\Pi^{\gamma\delta}_{g}$ terms and thus finding exact, simultaneous solutions for $\left\{\Pi^{\alpha\beta}_{k}\right\}$ is numerically intensive.

We note since the EM amplitude is much smaller than the bandwidth of the electronic energy spectrum ($A_{0}|J^{\perp}_{0,k}|\ll|\Delta(0)-\Delta(\pi)|$), $\Pi^{\alpha\beta}_{k}$ differs from its ground state value only in a small region (degeneracy region) of $O[A_{0}]$ in Brillouin zone. Hence scattering contributions between two degeneracy regions (being smaller by order of $A_{0}$) can be ignored.

Since $\Pi^{\alpha\beta}_{k}$ is close to its ground state value for almost all of the Brillouin zone, we can use the ground state values for $\Pi^{\alpha\beta}_{g}$ in the tensors $A^{\alpha\beta}_{k},D^{\alpha\beta}_{k},F^{\alpha\beta}_{k}$ to linearize the steady state equation of motion which can be readily solved analytically. The linearized equation of motion, for the $k$ points inside the degeneracy region, is obtained to be (see Appendix E for details)

[TABLE]

where $\tilde{\epsilon}^{ab}_{k}=\frac{1}{2}(\epsilon^{a}_{k}-\epsilon^{b}_{k})\delta^{ab}$ is the Hamiltonian in Floquet basis. $M_{k}$ and $\Lambda_{k}$ can be deduced from $A_{k}$ and $D_{k}$ as defined in Eq. 35-36, to be

[TABLE]

Here $p$ and $q$ are effective scattering rates for electrons and holes defined as:

[TABLE]

where $f_{>}$ is the fraction of the Brillouin zone such that $\Omega>2\Delta_{k}$. $\hat{S}_{v}$ and $\hat{S}_{c}$ are the projector operators on the valence and conduction band subspaces respectively defined as

[TABLE]

The solutions of Eq. 37 can be used to calculate current responses.

IV.2 Expression for current response

The matrix elements of current operator in the basis of Floquet states are

[TABLE]

Solving the linear equations in Eq. 37 and calculating $C^{\alpha\beta}_{k}$ according to Eq. 41 (Appendix F contains expressions for the steady state matrix elements, $\Pi^{\alpha\beta}_{k}$ and $C^{\alpha\beta}_{k}$), we obtain expression for the current response. The AC current can be obtained by computing the $\pm\Omega$ frequency components of $C^{\alpha\beta}_{k}$ and combining with the steady state solution for $\Pi^{\alpha\beta}_{k}$ (see Appendix F). The resultant response for the AC current up to leading order in electric field amplitude is

[TABLE]

where $t_{k}$ is given by

[TABLE]

Similarly, computing the zero frequency component of $C^{\alpha\beta}_{k}$, we find the DC current up to second order in electric field magnitude to be

[TABLE]

The total DC current is obtained integrating Eq. 44 over Brillouin zone as

[TABLE]

If the EM perturbation is the smallest parameter (i.e. $A_{0}|J^{\perp}_{0,k}|\ll p+q$), from Eq. 43 $t_{k}\simeq 1$ which implies

[TABLE]

Recognizing the total scattering rate as $\eta=p+q$ and using the limit, $\lim_{A_{0}\to 0}\epsilon_{k}=\frac{\Omega}{2}-\Delta_{k}$, we recover the semi classical result obtained in Eq. 13.

$J^{DC}_{shift}(k)$ can also be evaluated in the opposite limit of weak scattering (i.e. $p,q\ll A_{0}|J^{\perp}_{0,k}|$) to obtain

[TABLE]

Putting this expression into Eq. 45 it is easy to see that the resultant DC current scales linearly with the electric field in this limit, i.e. $J^{DC}_{shift}\propto A_{0}$. This result contradicts Kubo formula by producing a linear scaling with the external oscillating EM field for the DC current.

V Numerical evaluation of the weak damping limit

According to the analytic results in Eq. 46 and Eq. 47, the DC current response scales quadratically at small $A_{0}$ and linearly at large $A_{0}$ (smaller and larger than e-ph scattering strengths). To study how the crossover occurs between the two regimes, we perform the integral in Eq. 45 using the expression for $J^{DC}_{shift}(k)$ given by Eq. 44 for the Rice-Mele model Rice1982 . We also numerically study the scaling property of AC current magnitude for the same model. The Rice-Mele model is described by the following tight binding Hamiltonian Rice1982

[TABLE]

AAAABBBtttt’t’t’

The EM vector potential $A(t)$ introduces phases to the hoping terms according to Peierls substitution as

[TABLE]

Transforming to Fourier space, the perturbed Hamiltonian, in the atomic ($a_{k},b_{k}$) basis, becomes

[TABLE]

In terms of system parameters, $\left(\vec{J}_{0,k}\crossproduct\vec{J}_{1,k}\right)_{z}=-\frac{d}{8\Delta_{k}}\left(h^{\prime 2}-h^{2}\right)$ where $\vec{J}_{0,k}\equiv\partial_{A}\vec{H}_{k}(A)|_{A=0}$ and $\vec{J}_{1,k}\equiv\partial^{2}_{A}\vec{H}_{k}(A)|_{A=0}$.

The AC current for this Rice-Mele model is obtained by numerical integration of the expression given by Eq. 42 over the Brillouin zone with fixed values for $p,q$. The resulting values for $|J^{AC}|/|J^{\perp}_{0,k_{D}}|$ are plotted as a function of $2A_{0}|J^{\perp}_{0,k_{D}}|/(p+q)$ on a log scale in Fig. 3 where $k_{D}$ is the momentum of the degeneracy point (with the property, $2\Delta_{k_{D}}=\Omega$). The points on the plot fit to a straight line with slope 1. Thus the AC current response scales linearly in a range of values from $A_{0}|J^{\perp}_{0,k_{D}}|\ll(p+q)$ to $A_{0}|J^{\perp}_{0,k_{D}}|\gg p+q$. This scaling is as expected from conventional linear response theory. On the other hand same plot for the shift current, $|J^{DC}_{shift}|/|J^{\perp}_{0,k_{D}}|$, shows a crossover of scaling with $A_{0}$ as shown in Fig. 4. In the limit of weak electric field, $\left(ln\left|2A_{0}|J^{\perp}_{0,k_{D}}|/(p+q)\right|<-1\right)$, the points fit to a straight line of slope 2 whereas in the limit of weak scattering $\left(ln\left|2A_{0}|J^{\perp}_{0,k_{D}}|/(p+q)\right|>1\right)$, the corresponding fit line has slope 1. Therefore, the DC current response indeed shifts scaling from quadratic to linear through a region where the electric field and the scattering strengths are comparable $\left(ln\left|2A_{0}|J^{\perp}_{0,k_{D}}|/(p+q)\right|\sim 0\right)$. Thus the e-ph coupling strength sets an energy scale for the EM field strength for a non-linear to linear crossover of the shift current response.

VI Conclusion

We have studied an example of a non-linear response i.e. the BPVE in a system with electron-phonon interactions. Since the response is dominated by resonant transitions, we have limited our treatment to a two band model. As we found in Sec. II, the BPVE contains OTOC terms that cannot be computed by simple diagrams. However, in Sec. III, we find that the BPVE can be understood semiclassically from the dipole moment of generated excitons. The semiclassical limit assumes that the scattering rate is large compared to the driving force. We remedy this in Sec. IV using a master equation approach. From this we find that in the small scattering rate limit the DC current scales linearly with the vector potential instead of quadratically. This signals that this result is beyond the response function formalism in Sec. II. The AC current doesn’t show any such anomalous scaling behavior.

The DC current we obtain in the small $A_{0}$ limit is consistent with the non-interacting result for BPVE obtained previously. The breakdown seen in Sec. IV at larger $A_{0}$ shows that electron-phonon interaction affects the result in an implicit way. However, what is interesting from comparing to the general response in Sec. II is that only the $\chi_{n,2}$ term contribute to the shift current in configurations with homogeneous excitation profiles. This term is really conventional linear response of $\hat{J}_{0}$ to $\hat{J}_{1}$ contributions in the BPVE and can be computed from equilibrium field theory. However, as shown in Sec. III.2, the situation is quite different for inhomogeneous excitation profiles. In this case, owing to inhomogeneous carrier density generated by the incident light an additional diffusion current response is generated. This modulation of the density profile provides an alternative way to measure the shift current capacitively without the use of contact leads. The diffusion current contribution is not captured by the simple shift current that can be computed from $\chi_{n,2}$ and requires the computation of $\chi_{otoc}$. As discussed in Sec. II, since $\chi_{otoc}$ is not a time-ordered correlator, it cannot be computed from conventional equilibrium field theory but rather requires a Keldysh technique. In this work, we have circumvented this using a semi-classical approach that is of limited validity in the context of weak interactions to place a constraint on the long ranged part of $\chi_{otoc}$ (Eq. 22). The full quantum mechanical computation of $\chi_{otoc}$ using Keldysh field theory for strongly interacting systems where the semiclassical estimate would break down, would be an interesting future direction.

We acknowledge interesting comments from Victor Yakovenko at the inception of this project. This work has been supported by the National Science Foundation NSF DMR1555135 (CAREER) and JQI-NSF-PFC (supported by NSF grant PHY-1607611). Jay Sau is grateful to the Aspen Center for Physics where part of the work was completed.

Appendix A Derivation of second order response

We will derive the current response, $\langle\hat{J}(r)\rangle$, up to second order in the excitation amplitude, $A_{0}(x)$, due to the EM perturbation, $\int dx\,A(r)\hat{J}(r)$, as given by Hamiltonian in Eq. 1. The space-time dependence of $A(r)$ is defined in Sec. II. The response can be written as the trace of current operator with the density matrix, as

[TABLE]

where $\mathscr{J}(r),\;\tilde{\rho}(t)$ both are in interaction picture, i.e. evolved by the unperturbed Hamiltonian for the system $\hat{H}_{0}$ as in Eq. 1. The evolution of $\tilde{\rho}(t)$ is thus given by:

[TABLE]

where $\tilde{\rho}_{-\infty}$ is initial equilibrium density matrix for the system, before the EM field was turned on. $\tilde{U}(t)$ is the interaction picture time evolution operator which can be expanded up to second order in $A_{0}(r)$ from the time-ordered exponential as follows:

[TABLE]

Using this expansion for $\tilde{U}(t)$ in Eq. 51 and using cyclic property of trace, we obtain for $\langle\hat{J}(r)\rangle$, up to second order in $A_{0}(x)$:

[TABLE]

where the linear response function, $\chi_{lin}$ is defined as

[TABLE]

The non-linear response function $\chi_{n,2}$ is given by

[TABLE]

$\chi_{n,3}$ contains equilibrium time-ordered correlators of three current operators and is given by

[TABLE]

$\chi_{otoc}$ contains out of time ordered correlator as

[TABLE]

$\langle..\rangle_{0}=Tr\left\{\tilde{\rho}_{-\infty}..\right\}$ refers to the ground state expectation value. Since there’s no current in equilibrium, the first term vanishes (i.e. no current response in absence of any excitation) and one obtains Eq. 2.

Appendix B Derivation of shift vector

In this appendix we establish the agreement of our result for the DC current as obtained in Eq. 13 to the conventional form of the same [12, 13, 14] as given by,

[TABLE]

where $u^{cv}_{k}=\langle c|\hat{u}_{k}|v\rangle$ and $\hat{u}_{k}$ is the velocity operator defined as, $\hat{u}_{k}=\partial_{k}\hat{H}_{0,k}$. $\alpha^{cc}_{k}$ ($\alpha^{vv}_{k}$) is the Berry connection for the $|c_{k}\rangle$ ($|v_{k}\rangle$) state and $\phi_{k}=-\text{Im}[\log u^{cv}_{k}]$ . Recalling our result, Eq. 13:

[TABLE]

we recognize that we need to show

[TABLE]

First we observe that

[TABLE]

where from the definition of current operators, we have used $\hat{J}_{0,k}=-\hat{u}_{k}$ and $\hat{J}_{1,k}=\partial_{k}\hat{u}_{k}$. Using chain rule for derivatives, we expand the off-diagonal matrix element of $\partial_{k}\hat{u}_{k}$ as:

[TABLE]

$\langle\partial_{k}c|\hat{u}_{k}|v\rangle$ can be expanded as:

[TABLE]

Now inserting Identity operator, $\hat{I}=|c\rangle\langle c|+|v\rangle\langle v|$ in between for the second and third terms in Eq. 63 we obtain

[TABLE]

where $a^{\alpha\beta}_{k}=\langle\alpha|\partial_{k}\beta\rangle$, $u^{vv}_{k}=\partial_{k}\epsilon^{v}_{k}$. Interchanging $c$ and $v$ in Eq. 64 we obtain

[TABLE]

We observe that the Berry phases are related to the $a^{\alpha\beta}_{k}$ as $\alpha^{\alpha\beta}_{k}=-\iota a^{\alpha\beta}_{k}$. The off-diagonal term $a^{cv}_{k}$ can be related to the velocity matrix element as

[TABLE]

Using complex conjugate of Eq. 65in Eq. 62 we obtain:

[TABLE]

Writing $u^{cv}_{k}=|u^{cv}_{k}|e^{-\iota\phi_{k}}$ we obtain for the imaginary part, $\text{Im}\left[u^{vc}_{k}\left(\partial_{k}\hat{u}_{k}\right)^{cv}\right]$:

[TABLE]

where the last equality is the connection of the shift vector to the berry connection as defined in previous literature [12, 13, 14].

Appendix C Floquet Theory for two band system

For a system with time periodic Hamiltonian of frequency $\Omega$, the solutions for the time dependent Schrodinger equation (TDSE) have the following form

[TABLE]

where the quasi energy $\epsilon^{\alpha}_{k}$ can be chosen to belong in a $[0,\Omega]$ interval and $|\phi^{\alpha}_{k}(t)\rangle$ is time periodic with period $\frac{2\pi}{\Omega}$.

Substituting the form of Eq. 69 in TDSE, one obtains the Floquet equations in frequency space as

[TABLE]

where $\hat{H}^{k}_{m-n}$ corresponds to the $(m-n)\Omega$ frequency component of the Hamiltonian, $\hat{H}_{k}$. The eigenvalue solutions for Eq. 70 produce $\epsilon^{\alpha}_{k}$ and $|\psi^{\alpha}_{k}(n)\rangle$.

For the two band case discussed in the main text, the EM perturbation in frequency space can be written as

[TABLE]

where the subscript $n$ correspond to $n\Omega$ frequency component for the $c_{k}$, $v_{k}$ operators .

For the unperturbed Hamiltonian (Eq. 7) the quasi energies are simply

[TABLE]

with the corresponding states as

[TABLE]

If $\Omega=2\Delta_{k}$, the quasi energies are degenerate and hence to obtain the effect of perturbation (Eq. 71) an exact diagonalisation needs to be performed in the vicinity of those degeneracy points.

Up to leading order in $A_{0}$, the shifts can be obtained by diagonalizing the following matrix given by

[TABLE]

The corresponding eigenvectors can be used to obtain the Floquet state solutions in time domain to be

[TABLE]

Appendix D Derivation of equation of motion for $\Pi^{ab}_{k}$

Expanding the solution of Eq. 29 up to second order in $|V_{1,2}|$ we obtain:

[TABLE]

Since, the first term contains single power of bath operators, its trace can be dropped. In that case we note, $|\hat{\rho}_{I}(t^{\prime\prime})-\hat{\rho}_{I}(t^{\prime})|\approx O[|V_{1,2}|^{2}]$ and so up to second order in $|V_{1,2}|$, we can replace $\hat{\rho}_{I}(t^{\prime\prime})$ by $\hat{\rho}_{I}(t^{\prime})$ in Eq. 76 to obtain a second order perturbative result for $\hat{\rho}^{e}_{I}(t)$ as

[TABLE]

Differentiating Eq. 77 w.r.t. time one obtains Eq. 30.

Using cyclic property of trace, from Eq. 31, $\partial_{t}\Pi^{ab}_{k}$ produces

[TABLE]

Noting that $U_{F}(t)\;f^{\alpha\dagger}_{k}f^{\beta}_{g}\;U^{\dagger}_{F}(t)=e^{\iota(\epsilon^{\alpha}_{k}-\epsilon^{\beta}_{g})t}\;f^{\alpha\dagger}_{k}f^{\beta}_{g}$ and $U_{p}(t)\;a^{s}_{q}\;U^{\dagger}_{p}(t)=a^{s}_{q}\;e^{\iota\omega^{s}_{q}t}$, we obtain

[TABLE]

where the interaction matrix elements $V^{\alpha\beta}_{kg}$ between the Floquet states have the form:

[TABLE]

And the integral over time for the interaction gives

[TABLE]

where $\eta$ is an arbitrarily small imaginary number introduced to regularize the integrals. Considering only the Cauchy principal value one obtains energy conserving delta functions as follows

[TABLE]

Going back to Eq. 78, the first commutator in the second term gives

[TABLE]

In order to calculate the second commutator we need to trace out the bath operators from the product of $\commutator{\Theta^{ab}_{k}}{\hat{V}_{I}(t)}\int_{-\infty}^{t}dt^{\prime}\,\hat{V}_{I}(t^{\prime})$. We use Markovian approximation which assumes that the bath’s response time ($\tau_{p}$) is much shorter than the dynamical time scale due to interaction Hamiltonian, $\tau_{p}\ll\frac{1}{|V|}$. Consequently, the phonon density matrix remains effectively constant and can be decoupled from the combined density operator as

[TABLE]

The decoupling implies $\langle aa^{\dagger}f^{\dagger}ff^{\dagger}f\rangle\simeq\langle aa^{\dagger}\rangle\langle f^{\dagger}ff^{\dagger}f\rangle$.

Now for the bath using zero temperature expectation values, i.e. $\langle a^{s}_{q}a^{s^{\prime}\dagger}_{q^{\prime}}\rangle=\delta^{ss^{\prime}}_{qq^{\prime}}$ and assuming continuum energy spectrum with a constant density of states, $G$, for each momentum, i.e. $\sum^{s}=G\int_{0}^{\infty}d\omega_{q}$, we obtain for the trace of $\commutator{\Theta^{ab}_{k}}{\hat{V}_{I}(t)}\int_{-\infty}^{t}dt^{\prime}\,\hat{V}_{I}(t^{\prime})$,

[TABLE]

Similar expression can be obtained for the term, $Tr\left\{U_{0}\left(\int_{-\infty}^{t}dt^{\prime}\,\hat{V}_{I}(t^{\prime})\right)\commutator{\Theta^{ab}_{k}}{U_{0}\hat{V}_{I}(t)U^{\dagger}_{0}}\right\}_{p}$.

While tracing out the electrons we assume the electronic system is non interacting and it has translational invariance so that

[TABLE]

Using these expressions in Eq. 85 and putting back to Eq. 78 we obtain the rate equation as in Eq. 32.

Appendix E Linearization procedure of master equation

As argued in the main text, away from the degeneracy region only the valence bands are completely occupied which implies

[TABLE]

where the states are written in frequency representation.

Putting $\Pi^{\alpha\beta}_{g}$ values as input, Eq. 32 becomes linear in $\Pi^{\alpha\beta}_{k}$ since we can ignore scatterings in between two degeneracy regions as its contribution is order of magnitude smaller.

Furthermore, using the state expressions we obtain for the interaction matrix elements,

[TABLE]

Similarly all other matrix element combinations ($V^{\alpha\beta}_{kg}(n)V^{\gamma\delta}_{gk}(-n)$) for different $g$ points ($\Omega>2\Delta_{g}\;\;\Omega<2\Delta_{g}$) can be obtained.

The additional $\Theta$ function puts constraint on the values of $n$ that goes into the summation formula in the steady state rate equation. Consequently the contributions from $g:\Omega<2\Delta_{g}$ become $O[A_{0}^{2}]$ and thus can be dropped. Incorporating all these, one obtains the $M_{k}$ and $\Lambda_{k}$ matrix elements as in Eq. 37.

Appendix F Steady state solutions and Current matrix elements

Using the expressions for Floquet states (Eq. 75) one can easily calculate the $M_{k}$ and $\Lambda_{k}$ matrix elements in Eq. 38 to find the solution for Eq. 37 in the form, $\Pi_{k}\equiv t_{k}\hat{I}+\vec{R}_{k}\cdot\vec{\sigma}$, given by

[TABLE]

The current operator, up to first order in $A_{0}$, is given by

[TABLE]

Using the Floquet state expressions (Eq. 75) along with the current operator definition in frequency space (Eq. 90), one can obtain the zero frequency component for the matrix elements of current operator, $\tilde{C}^{\alpha\beta}_{k}(0)$, in $\{\hat{I},\vec{\sigma}\}$ basis to be

[TABLE]

The $\pm\Omega$ frequency components of the matrix elements of current operator, $\tilde{C}^{\alpha\beta}_{k}(\pm\Omega)$, are given by

[TABLE]