The Berry phase rectification tensor and the solar rectification vector

Several remarkable relations between the quantum geometry and topology of electronic states in crystalline solids and their non-linear optical and transport properties have been unraveled over the last few decades. One of the second order effects intimately intertwined with the geometry of electronic wave-functions is the shift current [1–10], which is a rectified current arising in response to AC electric fields that drive optical transitions in crystals without inversion symmetry. The shift current arises from band-off-diagonal components of the velocity operator, in contrast to the injection current, which arises from the band-diagonal difference of group velocities between empty and occupied states [2, 3, 5–10]. In materials with time reversal symmetry (TRS), the shift current differs from the injection current in that only the shift current can lead to rectified currents in response to linearly polarized electric fields, making it attractive for photovoltaic applications [6, 11–14] since it can lead to a DC current for electric fields with no net degree of circular polarization such as those arriving from the Sun.

Another second order effect with an intimate connection to the geometry of electronic wave functions is the non-linear Hall effect driven by the Berry curvature dipole (BCD) [10, 15–17], that has been recently observed experimentally [18–23]. Like the shift current, the BCD non-linear Hall effect can also give rise to a rectified current in response to a linearly polarized electric field, but with the crucial difference that this effect occurs only in metals at low frequencies in a Drude-like fashion [10]. These two distinct effects, however, appear to have a connection, as suggested by the quantum-rectification sum rule (QRSR) derived in [10]. The QRSR states that the frequency integral of the total rectification conductivity for any time-reversal-invariant material is completely independent of energy dispersions and entirely controlled by the quantum geometric Berry connections of the bands. The contributions to the QRSR arise from the BCD effect, which accounts for all the intra-band weight in a Drude-like peak, and the shift current, which accounts for all the inter-band contributions. In contrast, injection currents do not contribute to the QRSR under time reversal invariant conditions.

In this work, we will show that the QRSR defines a three-index Berry phase rectification tensor (BPRT), which depends only on the Berry phase geometry of a time-reversal-invariant material. In d-dimensions the BPRT has units of $\mathrm{(length)}^{d-3}$, and therefore, it is dimensionless for 3D materials. We will introduce an operational definition of the BPRT by considering the change of polarization in response to ideal short pulses of electric fields (see figure 1) and compute it from first principle calculations for the Weyl semi-metal TaAs [24–27] and demonstrate that, remarkably, the net contribution from the BCD [28] is comparable to the inter-band contribution from shift currents. We will also compute the BPRT for the insulator LiAsSe₂, which is a promising photo-voltaic material [6]. Additionally we will introduce the notion of the solar rectification vector (SRV), whose magnitude serves as a technologically useful figure of merit for bulk photo-current generation, and is defined from the average DC current produced under ideal black-body solar radiation, analogous to the model employed by Shockley and Quessier (SQ) in their seminal study of p-n junction solar cells [29]. We will demonstrate that LiAsSe₂ could generate photo-currents as large as 5 mA cm⁻², which is about 5% of the ideal maximum from the SQ model [29]. The direction of the SRV for these materials is shown as a yellow vector in figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As we will illustrate, one of the appealing features of our formulas for computing the BPRT and the SRV, is that they can be readily combined with realistic first principle calculations of material band structures. Therefore, these formulas have the potential to be used to implement a more comprehensive and systematic search of materials with technologically desirable bulk photovoltaic properties.

Consider an electronic system which is in equilibrium for $t\lt0$, and at t = 0 is subjected to a quench of the vector potential, leading to a delta-function pulse in the electric field of the form:

$$\begin{equation} \mathbf{A}(t) = \mathbf{A}_0\Theta(t), \qquad \mathbf{E}(t) = -\mathbf{A}_0\delta(t), \end{equation} \tag{ 1 }$$

where Θ(t) is the heavy-side function. After the pulse, the net change of electric polarization can be obtained from:

$$\begin{equation} \Delta\mathbf{P} = \int_0^\infty \mathbf{j}(t) dt. \end{equation} \tag{ 2 }$$

This protocol is illustrated in figure 1. To second order in electric fields the current is related to the linear and second order conductivity as:

$$\begin{align} j^\gamma(t) &= \int_{-\infty}^{\infty}dt_1\sigma_{(1)}^{\gamma\alpha}(t,t_1)E^\alpha(t_1)\nonumber \\ & \quad+\int_{-\infty}^{\infty}dt_1\int_{-\infty}^{\infty}dt_2\sigma_{(2)}^{\gamma\alpha\beta}(t,t_1,t_2)E^\beta(t_1)E^\alpha(t_2), \end{align} \tag{ 3 }$$

where Greek letters denote space indices and here and throughout the paper sum over repeated indices is understood unless otherwise stated. Therefore, the expansion of the polarization in powers of A₀ is:

$$\begin{equation} \Delta P^\gamma = -\sigma^{\gamma\alpha}A_0^\alpha+\frac{e^3}{\hbar^2}R^{\gamma\alpha\beta}A_0^{\beta}A_0^{\alpha} + \ldots, \end{equation} \tag{ 4 }$$

where $\sigma^{\gamma\alpha}$ is the zero frequency linear conductivity, and $R^{\gamma\alpha\beta}$ defines the BPRT of the system, which is given by (see appendix A):

$$\begin{equation} \frac{e^3}{\hbar^2}R^{\gamma\alpha\beta} = \int_{-\infty}^{\infty}dt \sigma_{(2)}^{\gamma\alpha\beta}(t,0,0) = \int_{-\infty}^{\infty}\frac{d\omega}{2\pi} \sigma_{(2)}^{\gamma\alpha\beta}(-\omega,\omega).\end{equation} \tag{ 5 }$$

Equation (4) expresses the notion that the BRPT measures the net displacement of the average position of an electronic system to second order in the strength of delta function pulse in electric field. As demonstrated in [10] for the ground state of any non-interacting system with a time-reversal-invariant band structure, $R^{\gamma\alpha\beta}$ only contains information about the Berry connections of the bands. Specifically, in the clean limit of a time reversal invariant band structure, it is given by:

$$\begin{align} R^{\gamma\alpha\beta} = & \frac14\Big\{\left\langle\partial^\alpha\Omega^{\beta\gamma}\right\rangle+ \left\langle\left[A^{\alpha},i\partial^{\beta}A^{\gamma}\right]\right\rangle \nonumber\\ & \qquad +\left\langle\left[A^\alpha,\left[A^\beta,\bar{A}^\gamma\right]\right]\right\rangle+\left(\alpha\leftrightarrow \beta\right)\Big\}, \end{align} \tag{ 6 }$$

where $\partial^\alpha \equiv \partial/\partial k^\alpha$, ${\Omega}_{n}^{\alpha\beta} = \partial^\alpha {A}_{nn}^{\beta}-\partial^\beta {A}_{nn}^{\alpha}$ is the abelian Berry curvature of the nth band, and ${\bar{A}}^{\alpha}_{nm} = {{A}}^{\alpha}_{nm}(1-\delta_{nm})$ is the off-diagonal non-Abelian Berry connection and average is defined as:

$$\begin{equation} \langle \cdots \rangle = \sum_n \int \frac{d^d\mathbf{k}}{(2\pi)^d}f_n\langle n \mid(\cdots)\mid n\rangle, \end{equation} \tag{ 7 }$$

where f_n is the Fermi–Dirac occupation function of the nth Bloch band. The first term in equation (6) arises from the intra-band BCD while the two other terms arise from the inter-band shift current. Therefore the BPRT connects these two effects and allows for a concrete measure of the shift current effect, as a collective displacement of average position of the electron system in response to an electric field pulse. Such an interpretation of the BPRT suggests that time domain spectroscopic techniques [30, 31] are promising tools to measure the BPRT in materials. By construction the BPRT is symmetric in its last two indices, and the BCD does not contribute to fully diagonal components of the form $R^{\alpha\alpha\alpha}$. For a related, but more restricted, sum rule see [32].

Notice that the intra-band contribution to the BPRT, namely the Berry dipole contribution proportional to $\left\langle\partial^\alpha\Omega^{\beta\gamma}\right\rangle$ in equation (6), is gauge invariant as it only depends on the gauge invariant Berry curvature. Since the full BPRT tensor is also gauge invariant [10], it follows that the remainder of the tensor is separately gauge invariant. This remainder contains the contributions of inter-band transitions arising from shift-current effects [10]. Therefore, there is no gauge ambiguity in separating the intra- and inter-band contributions to the BPRT, and we will compute them separately for the semimetal TaAs in section 5.

In this section we introduce a figure of merit for bulk photovoltaic materials that captures their amount of DC current generated from Sun-light under ideal conditions, and therefore it is a more relevant quantity for technological applications. Consider an statistical ensemble electric fields that is time-translationally invariant, whose average and two-time autocorrelation functions are given by:

$$\begin{equation} \overline{E^\alpha(t_1)} = 0,\qquad \overline{E^\beta(t_1)E^\alpha(t_2)} = I^{\beta\alpha}(t_1-t_2).\end{equation} \tag{ 8 }$$

Now, from equation (3), it follows that the ensemble averaged current is given by:

$$\begin{equation} \overline{j^\gamma(t)} = \iint_{-\infty}^{\infty}dt_1dt_2\sigma_{(2)}^{\gamma\alpha\beta}(t,t_1,t_2)I^{\beta\alpha}(t_1-t_2). \end{equation} \tag{ 9 }$$

Then the average current will be time independent and given by:

$$\begin{equation} \overline{j^\gamma} = \int_{-\infty}^{\infty}d\omega\sigma_{(2)}^{\gamma\alpha\beta}(-\omega,\omega)I^{\beta\alpha}(\omega). \end{equation} \tag{ 10 }$$

Where $I^{\beta\alpha}(\omega)$ is the Fourier transform of the electric field auto-correlation function. In the next subsection, we will compute such electric field auto-correlation function for the ensemble of solar black-body radiation.

3.1. Autocorrelations for black-body radiation

The photon eigenmode representation of the electric field operator in Heisenberg representation and residing in a periodic box of volume V is:

$$\begin{align} \mathbf{E}(\mathbf{r},t) ={\,} & i \sum_{\mu,k}\sqrt{\frac{\hbar\omega_\mathbf{k}}{2\varepsilon_0V}}\bigg[\mathbf{e}_\mu(\mathbf{k})a_{\mathbf{k}}^\mu e^{-i(\omega_{\mathbf{k}}t-\mathbf{k}\cdot\mathbf{r})}\nonumber \\ &\qquad\quad\qquad\qquad -\mathbf{e}^*_\mu(\mathbf{k})a^{\dagger\mu}_{\mathbf{k}} e^{i(\omega_{\mathbf{k}}t-\mathbf{k}\cdot\mathbf{r})}\bigg], \end{align} \tag{ 11 }$$

where $\omega_{\mathbf{k}} = c|\mathbf{k}|$, $a_\mathbf{k}^\mu$, $a_\mathbf{k}^{\dagger\mu}$ are photon destruction and creation operators satisfying standard commutation relations, µ = x, y labels the directions of light polarization and the unit vectors satisfy:

$$\begin{align} \mathbf{e}_x\cdot\mathbf{k} = \mathbf{e}_y\cdot\mathbf{k} = \mathbf{e}_x\cdot\mathbf{e}_y = 0, \end{align} \tag{ 12 }$$

$$\begin{align} \sum_{\mu}\mathbf{e}^\alpha_\mu(\mathbf{k})\mathbf{e}^{\beta*}_{\mu}(\mathbf{k}) = \delta_{\alpha\beta}-\hat{\mathbf{k}}_{\alpha}\hat{\mathbf{k}}_{\beta}.\end{align} \tag{ 13 }$$

The quantum ensemble average of the photon number is given the by the Planck distribution:

$$\begin{equation} \langle a_\mathbf{k}^{\mu\dagger}a_{\mathbf{k}^{\prime}}^{\mu^{\prime}}\rangle = \delta_{\mathbf{k}\mathbf{k}^{\prime}}\delta_{\mu\mu^{\prime}} n_{\mathbf{k}} ,\qquad n_{\mathbf{k}} = \frac{1}{e^{\beta\hbar\omega_{\mathbf{k}}}-1}, \end{equation} \tag{ 14 }$$

where β = 1/k_B T is the inverse temperature of the black body. To compute the auto-correlation function of the electric field we use the following prescription for quantum-to-classical conversion of the averages:

$$\begin{equation} \overline{E^\alpha(\mathbf{r},t)E^\beta(\mathbf{r}^{\prime},t^{\prime})} \leftrightarrow \left\langle :E^\alpha(\mathbf{r},t)E^\beta(\mathbf{r}^{\prime},t^{\prime}):\right\rangle.\end{equation} \tag{ 15 }$$

Where in the left-hand side we have the correlations denoted by an over-bar, in which the electric field is viewed as a classical field, and in the right-hand side we have quantum mechanical averages in which the electric field is viewed as the operator appearing in equation (11). We have added a normal ordering of the quantum electric fields which removes the vacuum Casimir-type fluctuations which should be negligible in the thermodynamic limit $V \rightarrow \infty$. With this we obtain:

$$\begin{align} & \overline{{E}^\alpha(\mathbf{r},t){E}^{\beta}(\mathbf{r},t^{\prime}) } \nonumber \\ & \quad = \frac{\hbar}{\varepsilon_0V}\sum_{\mathbf{k}}(\delta_{\alpha\beta}-\hat{\mathbf{k}}_{\alpha}\hat{\mathbf{k}}_{\beta})\omega_\mathbf{k}n_\mathbf{k}\cos[\omega_\mathbf{k}(t-t^{\prime})].\end{align} \tag{ 16 }$$

In the above formula all modes of light contribute to produce a fully isotropic correlator of electric fields at a given point r. This would be the situation if the point of interest was fully immersed inside the isotropic black-body radiation ensemble. However, when the radiation that arrives at r is highly directional we can modify the formula above by a one with an extra weighing function inside the momentum integral $f(\hat{\textbf{k}})$:

$$\begin{align} & \overline{{E}^\alpha(\mathbf{r},t){E}^{\beta}(\mathbf{r},t^{\prime}) } \nonumber \\ & \quad = \frac{\hbar}{\varepsilon_0V}\sum_{\mathbf{k}}(\delta_{\alpha\beta}-\hat{\mathbf{k}}_{\alpha}\hat{\mathbf{k}}_{\beta})\omega_\mathbf{k}n_\mathbf{k}\cos[\omega_\mathbf{k}(t-t^{\prime})]{\,}f(\hat{\mathbf{k}}),\end{align} \tag{ 17 }$$

which phenomenologically takes into account the preferred directionality of the incoming light at a given point r. In the special case in which the function $f(\hat{\textbf{k}}) = 1$ only for a very small region of solid angles, ΔΩ, around a specific direction $\hat{\textbf{k}} = \hat{\textbf{k}}_0$, we have that:

$$\begin{align} & \overline{{E}^\alpha(\omega_1){E}^{\beta}(\omega_2) } \nonumber \\ & \quad = \frac{\Delta\Omega}{4\pi}\frac{\hbar(\delta_{\alpha\beta}-\hat{\mathbf{k}}_0^\alpha\hat{\mathbf{k}}_0^\beta)}{(2\pi)^2\varepsilon_0(\beta\hbar c)^3}\frac{(\hbar\beta|\omega_1|)^3}{e^{\hbar\beta|\omega_1|}-1}\delta(\omega_1+\omega_2).\end{align} \tag{ 18 }$$

For the light emitted from the Sun and arriving on earth, $\Delta\Omega = 4\pi (R_\mathbf{E}/R_\mathbf{S})^2$ and $\hat{\textbf{k}}_0$ would correspond to the unit vector defining the direction from Sun to earth, $R_\mathbf{S}$ is the radius of the Sun, and $R_\mathbf{E}$ is the distance of the Sun to the earth.

The formula above is expected to describe the ideal solar radiation before entering the earth atmosphere, but the solar irradiance spectrum on the surface of earth can still be reasonably well approximated by it (see e.g. [33]). Upon entering the atmosphere the spectral intensity gets reduced at specific frequencies corresponding to absorption of molecules in the atmosphere. There will also be a randomization of the transversality of the electric field due to elastic scattering. This could be accounted for by averaging the projector over some distribution of $\hat{\textbf{k}}_0$, but for simplicity we will simply replace the transverse projector in equation (18) by a delta function with all components of the electric field, namely we take the electric field correlator to be:

$$\begin{equation} \overline{{E}^\alpha(\omega_1){E}^{\beta}(\omega_2) } = \frac{\hbar\Delta\Omega}{2(2\pi)^3}\frac{\delta_{\alpha\beta}\delta(\omega_1+\omega_2)}{\varepsilon_0(\beta\hbar c)^3}\frac{(\hbar\beta|\omega_1|)^3}{e^{\hbar\beta|\omega_1|}-1}.\end{equation} \tag{ 19 }$$

As we will see in the section 3.2, as along as the material is oriented so that its 'solar rectification vector' (to be defined in next subsection) is orthogonal to $\hat{\textbf{k}}_0$, one obtains the same current from the equation (19) as the one one would obtain from the equation (18) that takes into account the transverse nature of incident light.

Therefore, in conclusion, the electric field auto-correlation function associated with solar radiation, and which is needed in order to compute the average DC current in equation (10), can be written as:

$$\begin{align} I^{\beta\alpha}(\omega) = \frac{\hbar\delta^{\beta\alpha}}{(2\pi)^2\varepsilon_0(\beta_S\hbar c)^3}\left(\frac{R_{\mathbf{S}}}{R_\mathbf{E}}\right)^2\mathcal{I}(\hbar\beta_S|\omega|),\end{align} \tag{ 20 }$$

$$\begin{align} \mathcal{I}(x) = \frac{x^3}{e^x-1},\end{align} \tag{ 21 }$$

where $\mathcal{I}(x)$ is the dimensionless Planck spectral distribution function, $\beta_S = 1/(k_BT_S)$ is the inverse temperature of the Sun T_S ≈ 5778K, $R_\mathbf{S}$ is the radius of the Sun, and $R_\mathbf{E}$ is the Sun-earth distance. As we will see in the next subsection, the geometric factor $(R_{\mathbf{S}}/R_\mathbf{E})^2$ will disappear from the normalized figure of merit that measures the electric current normalized to the incoming photon flux. They are needed however if one would like to compute the actual magnitude of the induced photo-current.

3.2. Definition of solar rectification vector

Equations (10) and (20) provide the total current density induced by exposing a material to solar radiation. To define a figure of merit for current generation, it is convenient to normalize such total current density (number of electrons per unit area per unit time) by the photon flux F_s (number of photons per unit area per unit time $[F_s] = m^{-2} s^{-1}$), so as to count how many current-carrying electrons are induced per incoming photon. Let us therefore also compute the photon flux on earth. The photon density, n is given by:

$$\begin{align} n^\mu(\mathbf{r}) = a^{\mu\dagger}(\mathbf{r})a^\mu{(\mathbf{r})},\end{align} \tag{ 22 }$$

$$\begin{align} a^\mu(\mathbf{r},t) = \frac{1}{\sqrt{V}}\sum_{\mathbf{k}}a^\mu_{\mathbf{k}}e^{-i\omega_{\mathbf{k}}t+i\mathbf{k}\cdot\mathbf{r}}.\end{align} \tag{ 23 }$$

The photon current can be obtained from the continuity equation for photon density:

$$\begin{equation} \partial_tn^\mu = -\boldsymbol{\nabla}\cdot\mathbf{F}^{\mu}.\end{equation} \tag{ 24 }$$

Which leads to the following expression for the photon current per unit of area per unit time:

$$\begin{align} F_x = \sum_\mu \mathbf{F}^\mu \cdot \hat{\textbf{x}} = \frac{c}{V}\sum_{\mu,\mathbf{k},k_x\gt0}\left(\hat{\mathbf{k}}\cdot\mathbf{x}\right)a^{\dagger\mu}_\mathbf{k}a^\mu_{\mathbf{k}}.\end{align} \tag{ 25 }$$

The expression above would give us the number of photons moving on a single direction (e.g. toward the right) on a given surface of a black-body, which can be pictured to be the surface of the Sun. On the surface of the Earth the photon flux is given by:

$$\begin{align} F_s & = \left(\frac{R_\mathbf{S}}{R_{\mathbf{E}}}\right)^2\frac{c}{V}\sum_{\mu,\mathbf{k},k_x\gt0}\left(\hat{\mathbf{k}}\cdot\mathbf{x}\right)\left\langle a^{\dagger\mu}_\mathbf{k}a^\mu_{\mathbf{k}}\right\rangle \nonumber \\ & = \left(\frac{R_\mathbf{S}}{R_{\mathbf{E}}}\right)^2\frac{2\zeta(3)}{(2\pi)^2}\frac{c}{(\beta\hbar c)^3}, \end{align} \tag{ 26 }$$

here ζ is the Riemann zeta-function.

Therefore, we define a three-component vectorial figure of merit called the SRV, and denoted by $\eta^\gamma$ (γ = x, y, z), as the ratio of the electric current density and the photon flux, given by:

$$\begin{equation} \eta^\gamma \equiv \frac{\overline{j^\gamma}}{eF_s}.\end{equation} \tag{ 27 }$$

Notice that the SRV defined above depends only on the temperature of the Sun and the intrinsic rectification properties of the material. By combining equations (18) and (26) with the definition from equation (27), one arrives at the following explicit form of the SRV:

$$\begin{equation} \eta^\gamma = \frac{\hbar}{2\zeta(3)ec\varepsilon_0}\int_{-\infty}^{\infty}d\omega~\sigma_{(2)}^{\gamma\alpha\alpha}(-\omega,\omega)\mathcal{I}(\hbar\beta_S|\omega|).\end{equation} \tag{ 28 }$$

Where $\mathcal{I}(x)$ is defined in equation (21). We see that $\eta^\gamma$ is a dimensionless vector for 3D materials and its magnitude, $\eta = \sqrt{\sum_i\eta_i^2}$, counts the number of electrons contributing to the rectified current per incoming photon, therefore serving as a dimensionless figure of merit for the ideal bulk photo-current generated by a material under short-circuit configurations. Because the Planck distribution has vanishing intensity at small frequency, in the clean limit, the BCD Drude-like peak does not contribute to the rectified current and all the contributions arise from shift currents in time-reversal invariant materials, while injection current contributions vanish (see appendix B).

It is instructive to contrast this figure of merit with a comparable figure of merit for current generation for the ideal model of a p-n junction solar-cell proposed by SQ [29]. In the most ideal version of the SQ model one assumes that every photon with energy above the band-gap is absorbed and generates an electron and a hole that are successfully collected at opposite electrodes contributing to the current generation of the material. Therefore in the SQ model the total current density is:

$$\begin{equation} j_{\mathrm{SQ}} = e F_s \eta_{\mathrm{SQ}}.\end{equation} \tag{ 29 }$$

Here $\eta_{\mathrm{SQ}}$ is the analogue of the magnitude of the SRV, but defined for a conventional solar cell within the SQ model, and simply measures fraction of incident photons that have energy above the band-gap and thus can be absorbed by the material. The flux of such photons with energies above the band-gap, Δ₀, can be computed in an analogous fashion to the total photon flux but inserting an extra factor to enforce that their energy is larger than Δ₀, as follows:

$$\begin{equation} F(\Delta_0) = \left(\frac{R_\mathbf{S}}{R_{\mathbf{E}}}\right)^2\frac{c}{V}\sum_{\mu,\mathbf{k}\atop k_x\gt0}\Theta(ck-\Delta_0)\left(\hat{\mathbf{k}}\cdot\mathbf{x}\right)\left\langle a^{\dagger\mu}_\mathbf{k}a^\mu_{\mathbf{k}}\right\rangle,\end{equation} \tag{ 30 }$$

where Θ is the Heaviside function. The ratio of absorbed photons, $\eta_{\mathrm{SQ}}$, in equation (29), can thus be obtained as $\eta_{\mathrm{SQ}} = F(\Delta_0)/F_s$. SQ were primarily interested in optimizing the total power delivered by the solar cell and not the total current. It is however clear that the maximum current density is attained for $\eta_{\mathrm{SQ}} = 1$ which is achieved in the limit of vanishing gap. As we will see LiAsSe₂ has a value of η = 0.05, and thus can produce a photocurrent that is 5% of the ideal SQ limit, in which every incoming photon produces an electron and a hole that do not recombine and are successfully captured at opposite electrodes contributing to the net current.

Our discussion of photo-current generation so far has assumed that light propagates inside the material as it would do in vacuum. This is a reasonable approximation only if the incoming radiation is not substantially modified inside the material. This can be justified for example if the material is much thinner than the decay length for light absorbtion in the material, $L \sim \lambda_{in}/\mathrm{Im}[n]$, where λ_in is the typical wavelength of the incoming radiation and n is the complex index of refraction of the material. However, a more realistic calculation of the photo-current would include the frequency dependence of light absorbtion and propagation in the material, which is often accounted for by the Glass coefficient [11, 34].

In this section we take into account that the intensity of electric field decays inside the material, and derive the expression for the total current. For simplicity, we assume that incident radiation is normal to the solar cell surface so that that the electric field is parallel to the surface, namely along the XY plane in figure 4, and that the induced photo-current also flows along this XY plane. To estimate the electric field penetration, we use linear response theory. The geometry of the problem is shown in figure 4, and by the standard procedure of boundary condition matching from Maxwell equations (see e.g. [35]) leads to the following equations:

$$\begin{align} \mathbf{E}_\mathrm{I}(\mathbf{z},t) = \mathbf{E}_\mathrm{I}e^{i(k_1{z}-\omega t)},\quad \mathbf{E}_\mathrm{R}(\mathbf{z},t) = \mathbf{E}_\mathrm{R}e^{-i(k_1{z}+\omega t)}, \end{align} \tag{ 31 }$$

$$\begin{align} \mathbf{H}_\mathrm{I}(\mathbf{z},t) = \mathbf{H}_\mathrm{I}e^{i({k}_1{z}-\omega t)},\quad \mathbf{H}_\mathrm{R}(\mathbf{z},t) = \mathbf{H}_\mathrm{R}e^{-i({k}_1{z}+\omega t)},\end{align} \tag{ 32 }$$

$$\begin{align} \mathbf{E}_{\mathrm{I,R}} = \left(E^x_{\mathrm{I,R}},E^y_{\mathrm{I,R}},0\right),\quad \mathbf{H}_{\mathrm{I,R}} = \left(H^x_{\mathrm{I,R}},H^y_{\mathrm{I,R}},0\right),\end{align} \tag{ 33 }$$

$$\begin{align} \mathbf{E}_\mathrm{T}(\mathbf{z},t) = \left({E}_{\mathrm{T}_{x}}e^{i(k_{2x}{z}-\omega t)},{E}_{\mathrm{T}_{y}}e^{i(k_{2y}{z}-\omega t)},0\right),\end{align} \tag{ 34 }$$

$$\begin{align} \boldsymbol{\nabla}\times\mathbf{H} = -i\mathbf{J}-\varepsilon_0\omega\mathbf{E},\qquad \boldsymbol{\nabla}\times\mathbf{E} = \omega\mu_0\mathbf{H},\end{align} \tag{ 35 }$$

$$\begin{align} J^x(\omega) = \sigma^{xx}_{(1)}(\omega)E^{x}(\omega),\quad J^y(\omega) = \sigma^{yy}_{(1)}(\omega)E^{y}(\omega).\end{align} \tag{ 36 }$$

For simplicity we assume the magnetic permitivities of the vacuum and the bulk photovoltaic material to be frequency independent constants. We also assume that the material has a point group so that its linear conductivity is diagonal in its set of principal axes, and we take one these axes to be the z-axis that coincides with the direction of incoming light. The frequency dependence of these diagonal components of the linear conductivity of the material can be obtained from standard linear response theory [10, 32] and are given by:

$$\begin{align} \sigma^{\beta\beta}_{(1)}(\omega) &= -\frac{e^2}{\hbar}\sum_{nm}\int \frac{d^3\mathbf{k}}{(2\pi)^3}\Bigg\{\delta_{nm}\frac{\partial \epsilon_{n}}{\partial k ^\beta}\frac{i\frac{\partial f_{n}}{\partial k ^\beta}}{\omega+i\Gamma}+i(f_{n}-f_{m})\nonumber \\ & \qquad\qquad\qquad\qquad\quad(\epsilon_{n}-\epsilon_{m})\frac{A^{\beta}_{nm}A^\beta_{mn}}{\omega-\epsilon_{n}+\epsilon_m+i\Gamma}\Bigg\},\end{align} \tag{ 37 }$$

where Γ is relaxation rate and _n is the nth Bloch band energy dispersion. In the expression above β index is not summed over. From the above one can obtain the dispersion relations of light outside and inside the material to be:

$$\begin{align} k_1^2 = \mu_0\omega^2\varepsilon_0,\quad k_{2x}^2 = \mu_0\omega^2\varepsilon^x_2(\omega),\quad k_{2y}^2 = \mu_0\omega^2\varepsilon^y_2(\omega), \end{align} \tag{ 38 }$$

$$\begin{align} \varepsilon^x_2 (\omega) = \varepsilon_0+i\frac{\sigma^{xx}_{(1)}(\omega)}{\omega},\quad \varepsilon_2^y (\omega) = \varepsilon_0+i\frac{\sigma^{yy}_{(1)}(\omega)}{\omega}.\end{align} \tag{ 39 }$$

By taking into account standard electromagnetic boundary conditions [35], we can obtain the following relation describing the transmitted, $\mathrm{T}$, component of the electric field inside the material for a given incident component, $\mathrm{I}$:

$$\begin{align} \mathbf{E}_\mathrm{T}(\omega,z) \nonumber\\ & = \frac{2~n_1}{n_1+n^x_2(\omega)}(E^x_\mathrm{I}(\omega),0,0) e^{i(\mathrm{Re}\left[{k}_{2x}\right]{z}-\omega t)}e^{-\mathrm{Im}\left[{k}_{2x}\right]{z}} \nonumber \\ & \quad+\frac{2~n_1}{n_1+n^y_2(\omega)}(0,E^y_\mathrm{I}(\omega),0) e^{i(\mathrm{Re}\left[{k}_{2y}\right]{z}-\omega t)}e^{-\mathrm{Im}\left[{k}_{2y}\right]{z}},\end{align} \tag{ 40 }$$

$$\begin{equation} n_1 = c\sqrt{\mu_0\varepsilon_0},~ n^\beta_2 = c\sqrt{\mu_0\varepsilon^\beta_2(\omega)},~ k_{2\beta} = \omega n^\beta_2(\omega)/c.\end{equation} \tag{ 41 }$$

Since the above relation between transmitted and incoming light is linear, we can compute ensemble averages of the transmitted electric field from the known ensemble averages of the incoming electric field. We take the ensemble average of the incoming electric field to be that of the ideal black-body radiation from equation (18), where we keep track explicitly of the transverse character of incoming light. From this it follows that the electric field auto-correlation function inside the bulk of the solar cell is:

$$\begin{align} & \overline{E^\alpha_\mathrm{T}(\omega,z)E^\beta_\mathrm{T}(-\omega,z) }\nonumber\\ & = \delta_{\alpha\beta}\left|\frac{2n_1}{n_1+n^\beta_2(\omega)}\right|^2 e^{-2\frac{\omega}{c}\mathrm{Im}[n^\beta_2(\omega)]z}\overline{E^\beta_\mathrm{I}(\omega)E^\beta_\mathrm{I}(-\omega) }, \end{align} \tag{ 42 }$$

where in the above expression the repeated indices are not being summed over and indices are understood to take only (x, y) components. Thus the total current is:

$$\begin{align} \frac{\overline{I^\gamma}}{W} = \sum_{\beta = x,y}\int_{-\infty}^{\infty}d\omega |T^\beta(\omega)|^2\frac{\sigma_{(2)}^{\gamma\beta\beta}(-\omega,\omega)}{\alpha^\beta(\omega)}I^{\beta\beta}(\omega), \end{align} \tag{ 43 }$$

$$\begin{align} T^\beta(\omega) = \frac{2n_1}{n_1+n^\beta_2(\omega)},\quad \alpha^\beta(\omega) = 2\frac{\omega}{c}\mathrm{Im}[n^\beta_2(\omega)], \end{align} \tag{ 44 }$$

where W is the width of the sample, $n^\beta_2(\omega)$ is the complex index of refraction from equation (41), and $I^{\beta\beta}(\omega)$ is given in equation (20). The formula above generalizes that from [11, 34] by adding the frequency dependence of the reflection coefficient at the material surface, the spectral distribution of incident light as well as the directional dependence of the conductivity of the material. We expect it to provide accurate estimates for realistic bulk-photovoltaic solar cell devices.

In this section we combine our general considerations with first principle calculations of the non-linear opto-electronic properties of the 3D Weyl semi-metal TaAs and the insulator LiAsSe₂. Because it is semi-metallic, TaAs has BCD and shift current contributions whereas LiAsSe₂ only has shift current contributions. Our primary interest in TaAs is to investigate the interplay between the intraband BCD and interband shift-currents to its BPRT, although we also compute its SRV, while in the case of LiAsSe₂ our primary interest is the estimation of its SRV due to its promise as a photovoltaic material, but we also compute its BPRT. Throughout this section we will be using the ideal formalism that neglects the modification of the incoming radiation inside the material, namely we will compute the BPRT and SRV from equations (6) and (27), respectively.

The band structures of LiAsSe₂ and TaAs were calculated by the full-potential local-orbital [36, 37] minimum-basis code with the experimental crystal structure. Details of the calculations of band structure, Berry phases and non-linear rectification conductivity in the clean limit are presented in appendix C. The frequency dependent non-linear conductivity of LiAsSe₂ is shown in figure 3 (for the one of TaAs see appendix C). Our results are in agreement with previous reports [6, 38]. In particular, the σ^zxx and σ^zzz components for TaAs can exceed 300 µA V⁻² at infrared frequencies, while the σ^zzz component of LiAsSe₂ can exceed 100 µA V⁻² in a wide range frequency from ∼1.25 to ∼1.65 eV, which as we will see, leads to a sizable SRV. The space symmetry I4₁md of TaAs leads to the following non-zero components of the BPRT:

$$\begin{equation} \begin{array}{c} R^{xxz} = R^{yyz} = R^{xzx} = R^{yzy} = 0.13, \\ R^{zxx} = R^{zyy} = 0.05,\quad R^{zzz} = -0.04. \end{array} \end{equation} \tag{ 45 }$$

The values above contain the contribution from BCD and shift currents. The BCD part of the contribution to BPRT is:

$$\begin{equation} \begin{array}{c} R^{xxz}_{BCD} = R^{yyz}_{BCD} = R^{xzx}_{BCD} = R^{yzy}_{BCD} = 0.09, \\ \\ R^{zxx}_{BCD} = R^{zyy}_{BCD} = -0.18. \end{array} \end{equation} \tag{ 46 }$$

Therefore, remarkably, we see that in spite of having all its spectral weight localized in a delta function Drude-like peak, in TaAs the BCD contributes to the net rectification tensor with a magnitude that is comparable with the net contribution from shift currents, and sometimes they have opposite signs as it is the case for the $R^{zxx} = R^{zyy}$. The SRV for TaAs has a single non-trivial component $\boldsymbol{\eta}_{\mathrm{TaAs}} = (0,0,-0.001)$, which is oriented along its polar axes as depicted in figure 2(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the other hand, the lower Cc space group symmetry of LiAsSe₂ allows for a larger set of independent components:

$$\begin{equation} \begin{array}{c} R^{xxx} = 0.06,\qquad R^{xyy} = -0.05, \qquad R^{xzz} = -0.09,\\ R^{zxx} = 0.01,\qquad R^{zyy} = -0.05, \qquad R^{zzz} = 0.18,\\ R^{xxz} = R^{xzx} = -0.002, \qquad R^{yyx} = R^{yxy} = -0.04,\\ R^{yyz} = R^{yzy} = -0.02,\qquad R^{zzx} = R^{zxz} = -0.8. \end{array} \end{equation} \tag{ 47 }$$

The SRV has two components $\boldsymbol{\eta}_{\mathrm{LiAsSe_2}} = (-0.02,0,0.04)$, with the higher SRV magnitude η = 0.05, or in other words, it produces a photocurrent that is is 5% of the ideal limit for p-n junction solar cells.

We have introduced an operational definition of the BPRT as the second order polarization response of a material to an ideal delta function pulse of applied electric field. In the case of a time reversal invariant band structure this tensor is an intrinsic property of the material which only depends on the Berry connections of the Bloch bands and not on their energies. This tensor contains intra-band contributions which arise from the BCD and inter-band contributions which arise from the shift currents. This operational definition of the tensor is in principle applicable in the presence of electron–electron interactions, which is an interesting direction for future research. We have computed these tensors for TaAs, where we find that remarkably the Berry Dipole contributions are comparable to the shift current contributions.

We have also introduced a technologically relevant figure of merit for bulk photo-current generation termed the SRV, which counts the number of electrons contributing to the rectified current per incoming photon, under ideal solar black-body radiation, in analogy to the model of solar cells in the classic work of SQ [29]. Using this, we have found that LiAsSe₂ can produce a large bulk current of about 0.05 electrons per incident photon. We have also generalized the formula of photo-current production defining the Glass coefficient [11, 34] to include the black-body spectral distribution of intensity, the reflection coefficient at the surface of the material and the directional dependence of the linear conductivity of the material, which should provide accurate estimates of realistic photo-current production in solar cells based on bulk rectification mechanisms.

We would like to mention in passing that the formulae for rectification conductivities of [1, 2] are incomplete and are missing terms even for the inter-band shift currents. Details of the discrepancy and comparison with other various formulae employed in the literature [1, 2, 4, 5, 9, 32, 39–43] are presented in (see appendix D). For the complete correct formulae for the general second order conductivity of any band structure see [10].

All data that support the findings of this study are included within the article (and any supplementary files).

We thank Qian Niu for pointing to us the operational interpretation of the QRSR in terms of electric field pulses. We also thank Qian Niu, Cong Xiao, Jorge Facio, Jeroen van-den-Brink, Dennis Wawrzik, and Jhih-Shih You for valuable discussions.

In time domain, the current that is second order in electric fields can be written as:

$$\begin{equation} j_{2} ^\gamma(t) = \iint_{-\infty}^{\infty}dt_1 dt_2 E^{\beta}(t_2)E^{\alpha}(t_1)\sigma^{\gamma\beta\alpha}_{(2)}(t,t_1,t_2).\end{equation} \tag{ 48 }$$

The second order conductivity appearing above was derived in [10], and is given by:

$$\begin{equation} \sigma^{\gamma\beta\alpha}_{(2)}(t,t_1,t_2) = -\frac{e^3}{\hbar^2}e^{\Gamma(t_2-t)}\langle[r^\beta(t_2),[r^\alpha(t_1),v^{\gamma}(t)]]\rangle, \end{equation} \tag{ 49 }$$

for $t\geq t_1\geq t_2$ and 0 otherwise. This correlation function has global time translation invariance, namely:

$$\begin{equation} \sigma^{\gamma\beta\alpha}_{(2)}(t,t_1,t_2) = \sigma^{\gamma\beta\alpha}_{(2)}(t-\tau,t_1-\tau,t_2-\tau), \end{equation} \tag{ 50 }$$

which implies it has only two independent arguments when we shift the arguments by t (τ = t). Thus using the following convention for Fourier transforms:

$$\begin{align} f(t) = \int_{-\infty}^{\infty}d\omega f(\omega)e^{-i\omega t}, \end{align} \tag{ 51 }$$

$$\begin{align} f(\omega) = \int_{-\infty}^{\infty}\frac{dt}{2\pi}f(t)e^{i\omega t}, \end{align} \tag{ 52 }$$

we introduce the two argument conductivity in a frequency domain as:

$$\begin{align} & \sigma^{\gamma\beta\alpha}_{(2)}(0,\tau_1,\tau_2) \nonumber \\ & = \iint \frac{d\omega_1d\omega_2}{(2\pi)^2} e^{-i\omega_1\tau_1-i\omega_2\tau_2}\sigma^{\gamma\beta\alpha}_{(2)}(\omega_1,\omega_2).\end{align} \tag{ 53 }$$

The equation (53) is the definition of the $\sigma^{\gamma\beta\alpha}_{(2)}(\omega_1,\omega_2)$. Now, the expression above leads to the following expression for the current:

$$\begin{align} j_{2} ^\gamma(\omega) = & \int_{-\infty}^{\infty}d\omega_1 E^{\beta}(\omega-\omega_1)E^{\alpha}(\omega_1)\nonumber\\ & \times\sigma^{\gamma\beta\alpha}_{(2)}(-\omega_1,\omega_1-\omega).\end{align} \tag{ 54 }$$

Therefore the net rectified or DC current can be obtained from the above as the current at frequency ω = 0, and is given by:

$$\begin{equation} j_{2,\mathrm{DC}} ^\gamma = \int_{-\infty}^{\infty}d\omega E^{\beta}(-\omega)E^{\alpha}(\omega)\sigma^{\gamma\beta\alpha}_{(2)}(-\omega,\omega).\end{equation} \tag{ 55 }$$

By using the above definitions and conventions one can arrive at the equation (4) from the main text, namely, we have that:

$$\begin{align} \int_{-\infty}^\infty dt \sigma^{\gamma\beta\alpha}_{(2)}(t,0,0)& = \int_{-\infty}^\infty dt \sigma^{\gamma\beta\alpha}_{(2)}(0,-t,-t)\nonumber \\ & = \iiint_{-\infty}^\infty dt\frac{d\omega_1d\omega_2}{(2\pi)^2} \sigma^{\gamma\beta\alpha}_{(2)}(\omega_1,\omega_2)e^{i(\omega_1+\omega_2)t}\nonumber \\ & = \int_{-\infty}^{\infty}\frac{d\omega}{2\pi} \sigma_{(2)}^{\gamma\alpha\beta}(-\omega,\omega).\end{align} \tag{ 56 }$$

Now, the total current of the system is given by:

$$\begin{align} j^\gamma(t) = & \int_{-\infty}^{\infty}dt_1\sigma_{(1)}^{\gamma\alpha}(t,t_1)E^\alpha(t_1) \nonumber \\ & +\int_{-\infty}^{\infty}dt_1\int_{-\infty}^{\infty}dt_2\sigma_{(2)}^{\gamma\alpha\beta}(t,t_1,t_2)E^\beta(t_1)E^\alpha(t_2), \end{align} \tag{ 57 }$$

after averaging over ensemble, namely using relations:

$$\begin{equation} \overline{E^\alpha(t_1)} = 0,\qquad \overline{E^\beta(t_1)E^\alpha(t_2)} = I^{\beta\alpha}(t_1-t_2), \end{equation} \tag{ 58 }$$

we obtain the following expressions:

$$\begin{align} \overline{j^\gamma(t)} & = \iint_{-\infty}^{\infty}dt_1dt_2~\sigma_{(2)}^{\gamma\alpha\beta}(0,t_1-t,t_2-t)I^{\beta\alpha}(t_1-t_2)\nonumber\\ & = \int_{-\infty}^{\infty}d\omega\sigma_{(2)}^{\gamma\alpha\beta}(-\omega,\omega)I^{\beta\alpha}(\omega).\end{align} \tag{ 59 }$$

The second order rectification conductivity in general is given by:

$$\begin{align} \sigma^{\gamma\beta\alpha}_{(2)}(-\omega,\omega) & = \sigma^{\gamma\beta\alpha}_{\mathrm{J}}(-\omega,\omega)+\sigma^{\gamma\beta\alpha}_{\mathrm{BCD}}(-\omega,\omega)\nonumber\\ &\quad +\sigma^{\gamma\beta\alpha}_{\mathrm{I}}(-\omega,\omega)+\sigma^{\gamma\beta\alpha}_{\mathrm{S}}(-\omega,\omega), \end{align} \tag{ 60 }$$

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 ($\forall n,m,n\neq m : \Gamma \ll \epsilon_{n}-\epsilon_m$). Each contribution is given by:

$$\begin{equation} \sigma^{\gamma\beta\alpha}_{\mathrm{J}}(-\omega,\omega) = \frac{e^3}{\hbar^2}\int\frac{d\mathbf{k}}{(2\pi)^3}\sum_{nm}\frac{\frac{\partial \epsilon_n}{\partial k^\gamma}\frac{\partial^2 }{\partial k^\alpha\partial k^\beta}f_n\delta_{nm}}{\omega^2+\Gamma^2},\end{equation} \tag{ 61 }$$

$$\begin{align} \sigma^{\gamma\beta\alpha}_{\mathrm{BCD}}(-\omega,\omega) & = -\frac12\frac{e^3}{\hbar^2}\frac{1}{\omega+i\Gamma}\int\frac{d\mathbf{k}}{(2\pi)^3}\nonumber\\ & \quad\times\sum_{nm}\hat{A}^{\gamma}_{mn}\hat{A}^{\alpha}_{nm}\frac{\partial}{\partial k^\beta}(f_{m}-f_n)+\left(\begin{array}{c}\alpha\leftrightarrow \beta \\ \omega\leftrightarrow-\omega\end{array}\right). \end{align} \tag{ 62 }$$

$$\begin{align} \sigma^{\gamma\beta\alpha}_{\mathrm{I}}(-\omega,\omega)& = \frac{e^3}{\hbar^2}\int\frac{d\mathbf{k}}{(2\pi)^3}\nonumber \\ & \quad\times\sum_{nm}\frac{(f_{m}-f_n)\hat{A}^{\beta}_{nm}\hat{A}^{\alpha}_{mn}(\frac{\partial}{\partial k^\gamma}\epsilon_n-\frac{\partial}{\partial k^\gamma}\epsilon_m)}{(\omega-\epsilon_n+\epsilon_m)^2+\Gamma^2}\nonumber\\ & \quad+\left(\begin{array}{c}\alpha\leftrightarrow \beta\\ \omega\leftrightarrow-\omega\end{array}\right),\end{align} \tag{ 63 }$$

$$\begin{align} \sigma^{\gamma\beta\alpha}_{\mathrm{S}}(-\omega,\omega) &= \frac12\frac{e^3}{\hbar^2}\int\frac{d\mathbf{k}}{(2\pi)^3}\nonumber\\ &\quad\times\sum_{nm}\Bigg\{\hat{A}^{\gamma}_{mn}\frac{\partial}{\partial k^\alpha}\frac{(f_{n}-f_m)\hat{A}^{\beta}_{nm}}{\omega-\epsilon_{n}+\epsilon_{m}+i\Gamma} \nonumber\\& \quad + i\frac{(f_{n}-f_m)\hat{A}^{\beta}_{nm}}{\omega -\epsilon_{n}+\epsilon_{m}+i\Gamma}\sum_c \bigg[\hat{A}^\alpha_{mc}\hat{\bar{A}}^{\gamma}_{cn}-\hat{\bar{A}}^{\gamma}_{mc}\hat{A}^\alpha_{cn}\bigg]\Bigg\}\nonumber\\ &\quad+\left(\begin{array}{c}\alpha\leftrightarrow \beta\\ \omega\leftrightarrow-\omega\end{array}\right).\end{align} \tag{ 64 }$$

For purposes of ensemble averaging over incoming random light, we need to consider only diagonal components of the nonlinear rectification conductivity, namely $\sigma_{(2)}^{\gamma\alpha\alpha}$.

To simplify the notation for further analysis we write:

$$\begin{equation} Q^{\gamma\beta\alpha}_{nm} = \hat{A}^{\beta}_{nm} \frac{\partial}{\partial k^\alpha}\hat{A}^{\gamma}_{mn}-i\hat{A}^{\beta}_{nm}\sum_c \left[\hat{A}^\alpha_{mc}\hat{\bar{A}}^{\gamma}_{cn}-\hat{\bar{A}}^{\gamma}_{mc}\hat{A}^\alpha_{cn}\right].\end{equation} \tag{ 65 }$$

Thus, the conductivity which is the one that enters into the SRV, is:

$$\begin{equation} \sigma^{\gamma\alpha\alpha}_{\mathrm{J}}(-\omega,\omega) = \frac{e^3}{\hbar^2}\int\frac{d\mathbf{k}}{(2\pi)^3}\sum_{nm}\frac{\frac{\partial \epsilon_n}{\partial k^\gamma}\frac{\partial^2 }{\partial k^\alpha\partial k^\alpha}f_n\delta_{nm}}{\omega^2+\Gamma^2},\end{equation} \tag{ 66 }$$

$$\begin{align} \sigma^{\gamma\alpha\alpha}_{\mathrm{BCD}}(-\omega,\omega) &= \frac{1}{2}\frac{e^3}{\hbar^2}\frac{\Gamma}{\omega^2+\Gamma^2}\nonumber \\ & \quad\times\int\frac{d\mathbf{k}}{(2\pi)^3}\sum_{m}f_{m}\frac{\partial}{\partial k^\beta}\Omega_m^{\alpha\gamma},\end{align} \tag{ 67 }$$

$$\begin{align} \sigma^{\gamma\alpha\alpha}_{\mathrm{I}}(-\omega,\omega) = & \frac{e^3}{\hbar^2}\int\frac{d\mathbf{k}}{(2\pi)^3}\nonumber\\ & \qquad\times\sum_{nm}\frac{(f_{m}-f_n)\hat{A}^{\alpha}_{nm}\hat{A}^{\alpha}_{mn}(\frac{\partial}{\partial k^\gamma}\epsilon_n-\frac{\partial}{\partial k^k}\epsilon_m)}{(\omega-\epsilon_n+\epsilon_m)^2+\Gamma^2}, \end{align} \tag{ 68 }$$

$$\begin{align} \sigma^{\gamma\alpha\alpha}_{\mathrm{S}}(-\omega,\omega) = & \frac{e^3}{\hbar^2}\int\frac{d\mathbf{k}}{(2\pi)^3}\sum_{nm}(f_{m}-f_n)\nonumber\\ & \qquad\times\frac{\mathrm{Re}\left[Q^{\gamma\alpha\alpha}_{nm}\right](\omega-\epsilon_{n}+\epsilon_{m})+\Gamma\mathrm{Im}\left[Q^{\gamma\alpha\alpha}_{nm}\right]}{(\omega-\epsilon_{n}+\epsilon_{m})^2+\Gamma^2}.\end{align} \tag{ 69 }$$

TRS implies:

$$\begin{align} \mathcal{T}\mathbf{A}_{mn}(k)\mathcal{T}^{-1} = \mathbf{A}_{nm}(-k),\end{align} \tag{ 70 }$$

$$\begin{align} \mathcal{T}\mathbf{r}_{mn}(k)\mathcal{T}^{-1} = \mathbf{r}_{nm}(-k),\end{align} \tag{ 71 }$$

$$\begin{align} \mathcal{T}\mathbf{v}_{mn}(k)\mathcal{T}^{-1} = -\mathbf{v}_{nm}(-k).\end{align} \tag{ 72 }$$

TRS makes the Jerk, Injection and the part proportional to $\mathrm{Re}\left[Q^{\gamma\alpha\alpha}\right]$ of the Shift current vanish after momentum integration. Thus the contribution to the total rectification conductivity is given by BCD and resonant part of the shift current:

$$\begin{align} \sigma^{\gamma\alpha\alpha}_{\mathrm{BCD}}(-\omega,\omega)&+ \sigma^{\gamma\alpha\alpha}_{\mathrm{S}}(-\omega,\omega)\nonumber\\ & = \frac{e^3}{\hbar^2}\int\frac{d\mathbf{k}}{(2\pi)^3}\Bigg\{\frac{\Gamma}{\omega^2+\Gamma^2}\sum_{m}f_{m}\frac{\partial}{\partial k^\alpha}\Omega_m^{\alpha\gamma} \nonumber\\ & \quad+\sum_{nm}(f_{m}-f_n)\frac{\Gamma\mathrm{Im}\left[Q^{\gamma\alpha\alpha}_{nm}\right]}{(\omega-\epsilon_{n}+\epsilon_{m})^2+\Gamma^2}\Bigg\}.\end{align} \tag{ 73 }$$

In $\Gamma\rightarrow 0$ limit, we can simplify the expression using the following relation:

$$\begin{equation} \lim_{a\rightarrow0^+}\frac{a}{d_z^2+a^2} = \pi \delta(d_z).\end{equation} \tag{ 74 }$$

Now we can write a general version of the sum rule by including the spectral weight, which we assume to be a smooth function of ω, to obtain:

$$\begin{align} \int d\omega F(\omega)\sigma^{\gamma\alpha\alpha}(-\omega,\omega) &= \pi\frac{e^3}{\hbar^2}\sum_{\alpha,nm}f_{nm}F(|\epsilon_{nm}|)\nonumber\\ &\quad\times\Bigg\{A^{\alpha}_{nm}i\partial^\alpha A^\gamma_{mn}+\sum_l A^{\alpha}_{nm}\nonumber\\ &\quad\left[A^{\alpha}_{ml}A^{\mu}_{ln}-A^{\gamma}_{ml}A^{\alpha}_{ln}\right]\Bigg\},\end{align} \tag{ 75 }$$

where $f_{nm} = f_n - f_m$ and all zero-frequency contributions vanish under the assumption of absence of spectral weight ($F(\omega\rightarrow 0) = 0$).

Now in the case of the operational definition of the BPRT, we consider the net change of the polarization coursed by pulse-like electric field:

$$\begin{align} \mathbf{E}(t) = \mathbf{E}_0\tau_0\delta(t),\qquad \mathbf{E}(\omega) = \frac{\mathbf{E}_0\tau_0}{2\pi},\end{align} \tag{ 76 }$$

$$\begin{align} \Delta \mathbf{P} = \int_0^\infty \mathbf{j}(t)dt.\end{align} \tag{ 77 }$$

In this set-up, without assuming TRS the resulting change of the polarization is:

$$\begin{align} \Delta P^\gamma = & \frac{E_0^\alpha E_0^\beta}{4}\tau_0^2\frac{e^3}{\hbar^2}\nonumber \\ & \times\Bigg\{ \frac{1}{\Gamma}\left(\left\langle\partial^\alpha\partial^\beta v^{\gamma}\right\rangle+\sum_{ab}f_{ab}A^\beta_{ab}A^{\alpha}_{ba}\Delta^\gamma_{ba}\right) \nonumber \\ & \quad +\left\langle\left[A^{\beta},i\partial^{\alpha}A^{\gamma}\right]\right\rangle+\left\langle\partial^\beta\Omega^{\alpha\mu}\right\rangle \nonumber\\ & \quad\quad +\left\langle\left[A^\beta,\left[A^\alpha,\bar{A}^\gamma\right]\right]\right\rangle+\left(\alpha\leftrightarrow \beta\right)\Bigg\}, \end{align} \tag{ 78 }$$

where $\sum_{ab} = \sum_{ab}\int d \mathbf{k}/ (2\pi)^d$. As we see this result contains two extra terms proportional to 1/Γ which arise from the Jerk and injection current terms, but which vanish under the assumption of time-reversal symmetry.

The exchange and correlation energies were considered in the generalized gradient approximation with Perdew–Burke–Ernzerhof functional [36]. We projected the Bloch wave functions into high-symmetry atomic-orbital-like Wannier functions with diagonal position operator and considered all the electrons, which guarantees the full band overlap from ab-initio and Wannier functions in the energy window from −20 to 20 eV. The calculated band structures of TaAs and LiAsSe2 are shown in figure 5. Based on the highly symmetric Wannier functions, we constructed tight-binding (TB) model Hamiltonians and calculated the BCD and second order optical conductivity tensor. The γ (with α, β, γ = x, y, z) component of Berry curvature for the nth band at k point is computed by following:

$$\begin{equation} \Omega_{n}^\gamma(\mathbf{k}) = \frac{\varepsilon^{\gamma\alpha\beta}}{2}\left(\frac{\partial \hat{A}_{n}^{\beta}}{\partial k^\alpha}-\frac{\partial \hat{A}_{n}^{\alpha}}{\partial k^\beta}\right), \end{equation} \tag{ 79 }$$

where, $\hat{A}_{n}^{\beta} = \langle{n},\mathbf{k}|i\partial ^\beta|n,\mathbf{k}\rangle$ is the Berry connection of nth band of Hamiltonian $H(\mathbf{k})$ and $\varepsilon^{\gamma\alpha\beta}$ is the Levi-Civita tensor. The BCD was calculated from the Berry curvature in a differential manner and with integral in the whole $\mathbf{k}$-space [17]:

$$\begin{equation} D_{\alpha\beta} = \int \frac{d\mathbf{k}}{(2\pi)^3}\sum_n f_n\frac{\partial\Omega^\beta(\mathbf{k})}{\partial k^\alpha} = \left\langle\frac{\partial\Omega^\beta(\mathbf{k})}{\partial k^\alpha}\right\rangle, \end{equation} \tag{ 80 }$$

where $\Omega^\delta_n = \varepsilon^{\delta\alpha\beta}\partial A^\beta/\partial k^\alpha$ is the Berry curvature vector and f_n is a Fermi–Dirac distribution. The relation between two and three index BCD tensors is:

$$\begin{equation} \varepsilon^{\alpha\gamma\delta}D^{\beta\delta} = \partial^\beta \Omega^{\alpha\gamma}.\end{equation} \tag{ 81 }$$

Zoom In Zoom Out Reset image size

The second order optical conductivity tensor is calculated by the following expression:

$$\begin{align} \sigma^{\gamma\alpha\beta}_{(2)}(-\omega,\omega) &= \frac{\pi e^3}{2\hbar^2}\int\frac{d\mathbf{k}}{(2\pi)^3}\sum_{nm}(f_n-f_m)\nonumber \\ &\quad\qquad\quad\times\Big\{A^\beta_{nm}i\partial^\alpha A^\gamma_{mn}+A^\beta_{nm}\left[A^\alpha,\bar{A}^\gamma\right]_{mn} \nonumber\\ &\qquad\quad\quad+(\alpha\leftrightarrow\beta)\Big\}\delta\left(\omega-\frac{\epsilon_{n}-\epsilon_m}{\hbar}\right), \end{align} \tag{ 82 }$$

here _n is nth band energy. The expression above is equivalent to expressions discussed in [5, 10] (see appendix D). The numerical results are consistent with the symmetry analysis. The plots showing the frequency dependence of the second order rectification conductivities of TaAs and LiAsSe₂ are shown in figure 5. These plots show only the interband shift currents, since the BCD is simply a delta function peak at zero frequency.

After a k-grid convergency test, it was found that the change of BCD for TaAs is less than 5% with k-grid increasing from 360 × 360 × 360 to 480 × 480 × 480. A similar k-grid convergency was obtained with a k-grid of increasing from 240 × 240 × 240 to 360 × 360 × 360 for the second order optical conductivity in TaAs and LiAsSe₂.

Here we discuss the equivalence of our formulae for the rectification conductivity from equation (60) and others employed in the literature [1, 2, 4, 5, 9, 32, 39, 40]. Importantly, we show that shift currents derived in classic papers [1, 2] are only valid for Galilean systems, where spin orbit coupling is neglected, but not in general.

D.1. Comparison with [5, 39]

Here we will demonstrate that the interband rectification current in time reversal invariant materials arising from $\sigma_I + \sigma_S$ from equations (61)–(64), is identical to that from [5, 32]. This is the part of the rectification conductivity that controls the inter-band contribution to the QRSR and BPRT, but it is missing the intraband BCD contribution. This interband part however is the only one contributing to the SRV.

We define the following notation:

$$\begin{align} r^{a}_{nm} = \frac{v^a_{nm}}{i\epsilon_{nm}} = \bar{A}^a_{nm}, \qquad(n\neq m), \end{align} \tag{ 83 }$$

$$\begin{align} r^{a;b}_{nm} = \partial^b r^a_{nm}-ir^a_{nm}(A^b_{nn}-A^b_{mm}), \end{align} \tag{ 84 }$$

where $\epsilon_{nm} = \epsilon_n-\epsilon_m$ and (a, b, c) are space indices. Alternatively, this can be written as:

$$\begin{align} r^{a;b}_{nm} = & \frac{i}{\epsilon_{nm}}\Bigg\{\frac{v^a_{nm}\Delta^b_{nm}}{\epsilon_{nm}}+\frac{v^b_{nm}\Delta^a_{nm}}{\epsilon_{nm}}-\vartheta^{ab}_{nm}\nonumber\\ & \qquad\qquad+\sum_{p\neq n,m}\left(\frac{v^a_{np}v^b_{pm}}{\epsilon_{pm}}-\frac{v^a_{pm}v^b_{np}}{\epsilon_{np}}\right)\Bigg\}, \end{align} \tag{ 85 }$$

where $\Delta^a_{nm} = v_{nn}^a-v_{mm}^a$, and

$$\begin{align} \vartheta^{ab}_{nm}& = \langle{n}|\partial^a\partial^b H|m\rangle\nonumber\\ & = \partial^bv^a_{nm}+i\sum_p(v^a_{np}r^b_{pm}-r^b_{np}v^a_{pm}).\end{align} \tag{ 86 }$$

By using the identity identity $[r^a,r^b] = 0$, one obtains:

$$\begin{align} \partial^a A^b_{nm}-i[A^a,\bar{A}^b]_{nm} = \partial^bA^a_{nm}-iA^a_{nm}(A^b_{nn}-A^b_{mm}), \end{align} \tag{ 87 }$$

which leads to an alternative form to denote the generalized derivative:

$$\begin{align} r^{a;b}_{nm}& = \partial^b A^a_{nm}-iA^a_{nm}(A^b_{nn}-A^b_{mm})\nonumber\\ & = \partial^a A^b_{nm}-i\sum_p\left\{A^c_{np}\bar{A}^a_{pm}-\bar{A}^a_{np}A^c_{pm}\right\}, \end{align} \tag{ 88 }$$

for n ≠ m. By using equations (83) and (88) one can obtain that the conductivity from [5, 39] can be written as:

$$\begin{align} \sigma^{abc}(0,\omega,-\omega) = -\frac{i\pi e^3}{2\hbar^2}&\int\frac{d\mathbf{k}}{(2\pi)^3}\sum_{nm}f_{nm} \left(r^b_{mn}r^{c;a}_{nm}+(b \leftrightarrow c)\right)\nonumber\\ &\quad\times\delta(\frac{\epsilon_{mn}}{\hbar}-\omega), \end{align} \tag{ 89 }$$

this is the same as the sum of the real part of injection and shift currents $\sigma_{S} + \sigma_I$ assuming time-reversal invariant conditions from equation (60), discussed in section (see appendix B), which were the ones employed in the main text. In the main text, we are also omitting the '0' in the conductivity arguments, but it is implicit every time we refer to a 'rectification' conductivity.

D.2. Comparison with [1, 2, 9]

We will here compare our formula from the current paper and [10] to the formulae from the classic [1, 2] (see e.g. [38, 40] for an application of this) and those of more recent work of [9] (see also [41–43] for closely related formulae to [9]). We will show that [1, 2] are missing crucial corrections that can only be neglected if one assumes the system to have an underlying parabolic dispersion (namely ignores spin–orbit coupling). On the other hand we will show that, upon proper regularization of infinitesimal imaginary parts in the frequency denominators, the formulae of [9] are equivalent to those we employ in the current paper and in our previous work [10].

Following [9] we write the second order perturbation to the Hamiltonian, as:

$$\begin{align} \hat{H} = \hat{H}_0 + \hat{V}(t), \end{align} \tag{ 90 }$$

$$\begin{align} \hat{V}(t) = \frac{e}{\hbar}A^\alpha \hat{h}^\alpha+\frac12\left(\frac{e}{\hbar}\right)^2A^{\alpha}(t)A^{\beta}(t)\hat{h}^{\alpha\beta}+\cdots, \end{align} \tag{ 91 }$$

where $\hat{H}$ is total Hamiltonian, $\hat{H}_0$ is unperturbed band Hamiltonian, $\mathbf{E}(\omega) = i\omega\mathbf{A}(\omega)$ is electric field and

$$\begin{equation} h^{\mu}_{ab} = v^{\mu}_{ab} = \partial^{\mu}\epsilon_a \delta_{ab}+i\epsilon_{ab}A^{\mu}_{ab}, \end{equation} \tag{ 92 }$$

$$\begin{align} h^{\mu\alpha}_{ab} ={\,}& \partial^{\alpha}v^{\mu}_{ab}+iv^{\mu}_{ab}(A^{\alpha}_{bb}-A^{\alpha}_{aa}) \nonumber\\ & \quad+{{\sum_{c}}}^{\prime}\left\{\frac{v^{\alpha}_{ac}v^{\mu}_{cb}}{\epsilon_{ca}}-\frac{v^{\mu}_{ac}v^{\alpha}_{cb}}{\epsilon_{bc}}\right\}, \end{align} \tag{ 93 }$$

$$\begin{align} h^{\mu\alpha\beta}_{ab} ={\,}& \partial^{\beta}h^{\mu\alpha}_{ab}+ih^{\mu\alpha}_{ab}(A^{\beta}_{bb}-A^{\beta}_{aa})\nonumber \\ & +{{\sum_{c}}}^{\prime}\left\{\frac{v^{\beta}_{ac}h^{\mu\alpha}_{cb}}{\epsilon_{ca}}-\frac{h^{\mu\alpha}_{ac}v^{\beta}_{cb}}{\epsilon_{bc}}\right\}, \end{align} \tag{ 94 }$$

where ${\sum}^{\prime}$ means that one omits terms when denominator vanishes and latin sub-indices, (a, b, c) denote band labels, while upper greek indices (α, β, γ) denote space components.

From the above, the second order rectification conductivity is found to be:

$$\begin{align} \sigma^{\mu\alpha\beta}(0,\omega,-\omega)& = \frac{e^3}{\omega^2}\sum_{abc}\Bigg\{f_a h^{\mu\alpha\beta}_{aa}+f_{ab}\frac{h^{\alpha}_{ab}h^{\mu\beta}_{ba}}{\omega-\epsilon_{ab}}\nonumber\\ &\quad +f_{ab}\frac{h^{\beta}_{ab}h^{\mu\alpha}_{ba}}{-\omega-\epsilon_{ab}}+f_{ab}\frac{h^{\alpha\beta}_{ab}h^{\mu}_{ba}}{\epsilon_{ba}}+\frac{h^{\alpha}_{ab}h^{\beta}_{bc}h^{\mu}_{ca}}{\epsilon_{ac}}\nonumber\\ &\quad\times\left(\frac{f_{ab}}{\omega-\epsilon_{ba}}-\frac{f_{cb}}{\omega-\epsilon_{bc}}\right)+\frac{h^{\beta}_{ab}h^{\alpha}_{bc}h^{\mu}_{ca}}{\epsilon_{ac}}\nonumber\\ &\quad\times\left(\frac{f_{ab}}{-\omega-\epsilon_{ba}}-\frac{f_{cb}}{-\omega-\epsilon_{bc}}\right)\Bigg\}.\end{align} \tag{ 95 }$$

Here $f_{ab} = f_a-f_b$, with f_a the Fermi–Dirac occupation of band a. We will compare the expression above with the expression from [1, 2]. [1, 2] employed the same gauge of [9], but assumed that the electrons have an underlying parabolic kinetic energy in addition to the periodic potential, and thus entirely neglects spin–orbit coupling effects. Therefore, the second order perturbed Hamiltonian for a single electron in [1, 2] has the form:

$$\begin{align} \hat{H} = \hat{H}_0+\hat{V}(t),\qquad A^\alpha = \frac{E^\alpha_0}{\omega}\cos{\omega t} \end{align} \tag{ 96 }$$

$$\begin{align} \hat{V}(t) = -\frac{e}{m_0} \hat{p}^\alpha A^\alpha+\frac{e^2}{2m_0}A^\alpha A^\alpha\hat{\unicode{x1D7D9}}, \end{align} \tag{ 97 }$$

where m₀ is the bare electron mass. The second order rectification conductivity of [1, 2] is therefore found to be:

$$\begin{align} \sigma^{\mu}_{\alpha\beta}(\omega) = & \frac{|e|^3}{8\pi^2\omega^2}\nonumber\\ & \times\mathrm{Re}\left\{\sum_{\Omega = \pm\omega}\sum_{abc}f_{ba}\frac{v^{\alpha}_{ab}v^{\beta}_{bc}v^\mu_{ca}}{(\epsilon_{ac}-i\delta)(\epsilon_{ab}+\Omega-i\delta))}\right\}.\end{align} \tag{ 98 }$$

For purposes of comparing both expressions it is convenient to split them into pieces. The conductivity from [1, 2] we split as follows:

$$\begin{align} \mathrm{KB}_1 = \frac{e^3}{\omega^2}\mathrm{Re}\left[\sum_{abc}f_{ab}\frac{v^{\alpha}_{ab}v^{\beta}_{bc}v^\mu_{ca}}{(\epsilon_{ac}-i\delta)(\epsilon_{ab}+\omega-i\delta))}\right], \end{align} \tag{ 99 }$$

$$\begin{align} \mathrm{KB}_2 = \frac{e^3}{\omega^2}\mathrm{Re}\left[\sum_{abc}f_{ab}\frac{v^{\alpha}_{ab}v^{\beta}_{bc}v^\mu_{ca}}{(\epsilon_{ac}-i\delta)(\epsilon_{ab}-\omega-i\delta))}\right], \end{align} \tag{ 100 }$$

where KB stands for Kraut-Baltz from [1, 2]. We also split the conductivity from equation (95) [9] into the following terms:

$$\begin{align} \mathrm{PM}_0 = \frac{e^3}{\omega^2}\sum_{ab}\Bigg\{f_a h^{\mu\alpha\beta}_{aa} +f_{ab}\frac{h^{\alpha\beta}_{ab}h^{\mu}_{ba}}{\epsilon_{ba}}\Bigg\}, \end{align} \tag{ 101 }$$

$$\begin{align} \mathrm{PM}_1 = \frac{e^3}{\omega^2}\sum_{ab}\Bigg\{f_{ab}\frac{h^{\alpha}_{ab}h^{\mu\beta}_{ba}}{\omega-\epsilon_{ab}}+f_{ab}\frac{h^{\beta}_{ab}h^{\mu\alpha}_{ba}}{-\omega-\epsilon_{ab}}\Bigg\}, \end{align} \tag{ 102 }$$

$$\begin{align} \mathrm{PM}_2 = \frac{e^3}{\omega^2}\sum_{abc}\Bigg\{\frac{h^{\alpha}_{ab}h^{\beta}_{bc}h^{\mu}_{ca}}{\epsilon_{ac}}\left(\frac{f_{ab}}{\omega-\epsilon_{ba}}-\frac{f_{cb}}{\omega-\epsilon_{bc}}\right)\Bigg\}, \end{align} \tag{ 103 }$$

$$\begin{align} \mathrm{PM}_3 = \frac{e^3}{\omega^2}\sum_{abc}\Bigg\{\frac{h^{\beta}_{ab}h^{\alpha}_{bc}h^{\mu}_{ca}}{\epsilon_{ac}}\left(\frac{f_{ab}}{-\omega-\epsilon_{ba}}-\frac{f_{cb}}{-\omega-\epsilon_{bc}}\right)\Bigg\}, \end{align} \tag{ 104 }$$

where PM stands for Parker–Morimoto–Orenstein–Moore. Also we want to use one more index for each sub-term in these expressions. Namely, we will denote by $\mathrm{PM}_{1,1}$ the firs term from expression $\mathrm{PM}_{1}$, and analogously for other terms.

By comparing equations (90), (91) and equations (96) and (97), we see that for systems considered in [1, 2] certain additional constrains are satisfied for the matrix elements as a consequence of the underlying parabolic dispersion. More specifically, the interband rectification conductivity from [1, 2], can be obtained from the formula of [9], upon enforcing the following relations:

$$\begin{align} \hat{h}^\alpha = -\frac{\hbar}{m_0} \hat{p}^\alpha, \end{align} \tag{ 105 }$$

$$\begin{align} \hat{h}^{\alpha\beta} = \frac{\hbar^2}{e^2}\frac{e^2}{2m_0}\delta^{\alpha\beta}\hat{\unicode{x1D7D9}}, \end{align} \tag{ 106 }$$

$$\begin{align} \hat{h}^{\alpha\beta\gamma} = \hat{h}^{\alpha\beta\gamma\delta} = \cdots = 0.\end{align} \tag{ 107 }$$

From the expressions above one can see that the terms we have labeled PM₀ and PM₁ would vanish, because the tensor $\hat{h}^{\alpha\beta}_{ab}\sim \delta_{ab}\delta^{\alpha\beta}$, is diagonal in the (a, b) band indices, and therefore its contribution to the rectification conductivity vanishes after multiplying by the difference of Fermi occupation functions, $f_{ab} = f_a-f_b$, which are only non-zero for a ≠ b. Now one can verify that the remaining terms agree between the two theories, specifically as follows:

$$\begin{align} \mathrm{Re}\left[\mathrm{PM}_{21}\right] & = \frac{e^3}{\omega^2}\mathrm{Re}\left[\sum_{abc}\frac{h^{\alpha}_{ab}h^{\beta}_{bc}h^{\mu}_{ca}}{\epsilon_{ac}}\frac{f_{ab}}{\omega-\epsilon_{ba}}\right] \nonumber \\ & = \frac{e^3}{\omega^2}\mathrm{Re}\left[\sum_{abc}\frac{v^{\alpha}_{ab}v^{\beta}_{bc}v^{\mu}_{ca}}{\epsilon_{ac}}\frac{f_{ab}}{\epsilon_{ab}+\omega}\right] = \mathrm{KB}_1, \end{align} \tag{ 108 }$$

$$\begin{align} &\mathrm{Re}\left[\mathrm{PM}_{32}\right] = -\frac{e^3}{\omega^2}\mathrm{Re}\left[\sum_{abc}\frac{h^{\beta}_{ab}h^{\alpha}_{bc}h^{\mu}_{ca}}{\epsilon_{ac}}\frac{f_{cb}}{-\omega-\epsilon_{bc}}\right]\nonumber\\ & = \left|\begin{array}{c}\text{complex}\\ \text{conjugation}\end{array}\right| = -\frac{e^3}{\omega^2}\mathrm{Re}\left[\sum_{abc}\frac{h^{\beta}_{ba}h^{\alpha}_{cb}h^{\mu}_{ac}}{\epsilon_{ac}}\frac{f_{cb}}{-\omega-\epsilon_{bc}}\right]\nonumber\\ & = \left|a\leftrightarrow c\right| = -\frac{e^3}{\omega^2}\mathrm{Re}\left[\sum_{abc}\frac{h^{\beta}_{bc}h^{\alpha}_{ab}h^{\mu}_{ca}}{\epsilon_{ca}}\frac{f_{ab}}{-\omega-\epsilon_{ba}}\right]\nonumber\\ & = \frac{e^3}{\omega^2}\mathrm{Re}\left[\sum_{abc}\frac{h^{\beta}_{bc}h^{\alpha}_{ab}h^{\mu}_{ca}}{\epsilon_{ac}}\frac{f_{ab}}{\epsilon_{ab}-\omega}\right]\nonumber\\ & \qquad = \frac{e^3}{\omega^2}\mathrm{Re}\left[\sum_{abc}\frac{v^{\beta}_{bc}v^{\alpha}_{ab}v^{\mu}_{ca}}{\epsilon_{ac}}\frac{f_{ab}}{\epsilon_{ab}-\omega}\right] = \mathrm{KB}_2.\end{align} \tag{ 109 }$$

Expressions for $\mathrm{Re}[\textrm{PM}]_{22}$ and $\mathrm{Re}[\textrm{PM}]_{31}$ can be obtained from $\mathrm{Re}[\textrm{PM}]_{21}$ and $\mathrm{Re}[\textrm{PM}]_{32}$ respectively by swapping indices $\alpha\leftrightarrow\beta$. Such a symmetrization automatically implied by conductivity once it is contracted with electric field indices, as is used for current expression in [2]. Thus we conclude the conductivity from [2] is equivalent to the sum of terms $\mathrm{Re}[\textrm{PM}]_{2}$ and $\mathrm{Re}[\textrm{PM}]_{3}$.

Note that in general (for example in the presence of spin–orbit coupling) $\mathrm{Re}[\textrm{PM}]_{0}$ and $\mathrm{Re}[\textrm{PM}]_{1}$ do contribute to the conductivity. These terms include not only intra-band terms, but also some contribution to interband shift currents terms that are missing in [1, 2], and are given by:

$$\begin{align} &\mathrm{Re}\left[\frac{e^3}{\omega^2}\sum_{ab} f_{ab}\frac{iv^{\alpha}_{ab}v^{\mu}_{ba}(A^{\beta}_{aa}-A^{\beta}_{bb})}{\omega-\epsilon_{ab}}\right]\nonumber\\ & \quad+\mathrm{Re}\left[\frac{e^3}{\omega^2}\sum_{ab} f_{ab}\frac{v^{\alpha}_{ab}}{\omega-\epsilon_{ab}}{{\sum_{c}}}^^{\prime}\left\{\frac{v^{\beta}_{bc}v^{\mu}_{ca}}{\epsilon_{cb}}-\frac{v^{\mu}_{bc}v^{\beta}_{ca}}{\epsilon_{ac}}\right\}\right]\nonumber\\ & \qquad\qquad+\left(\begin{array}{c}\alpha\leftrightarrow\beta \\ \omega\leftrightarrow-\omega\end{array}\right).\end{align} \tag{ 110 }$$

Or in notations used in [1, 2]:

$$\begin{align} &\mathrm{Re}\left[\frac{e^3}{2\omega^2}\sum_{\Omega = \pm\omega}\sum_{ab} f_{ab}\frac{iv^{\alpha}_{ab}v^{\mu}_{ba}(A^{\beta}_{aa}-A^{\beta}_{bb})}{\Omega+\epsilon_{ba}-i\delta}\right]\nonumber\\ &\quad+\mathrm{Re}\left[\frac{e^3}{\omega^2}\sum_{\Omega = \pm\omega}\sum_{ab}{{\sum_{c}}}^{\prime} f_{ab}\frac{v^{\alpha}_{ab}}{\Omega+\epsilon_{ba}-i\delta} \frac{v^{\beta}_{bc}v^{\mu}_{ca}}{\epsilon_{cb}}\right]\nonumber\\ & \qquad\qquad+\left(\begin{array}{c}\alpha\leftrightarrow\beta \end{array}\right).\end{align} \tag{ 111 }$$

To show that the extra shift current contribution, which we demonstrate above, is crucial—we will compute our Quantum Rectification Sum Rule for shift currents described in [1, 2] and [9]. With the proper regularization of the infinitesimal imaginary parts of the frequency denominators of both of these references, the frequency integration of interband conductivities are:

$$\begin{align} & 2\frac{\hbar^2}{e^3}\frac{1}{\pi}\int_0^{\infty} \mathrm{Re}\left[\mathrm{PM}^{\mu\beta \alpha}(0,-\omega,\omega)\right] d\omega \nonumber\\ &\quad = \sum_af_{a}\left(i[A^{\beta}\partial^\alpha A^{\mu}]_{aa}+[A^{\beta} ,[A^\alpha,\bar{A}^{\mu}]_{aa}\right)\nonumber\\ &\qquad\qquad+\left(\begin{array}{c}\alpha\leftrightarrow\beta \end{array}\right),\end{align} \tag{ 112 }$$

$$\begin{align} & 2\frac{\hbar^2}{e^3}\frac{1}{\pi}\int_0^{\infty} \mathrm{Re}\left[\mathrm{KB}^{\gamma\beta \alpha}(0,-\omega,\omega)\right]d\omega \nonumber\\ &\quad = \sum_{abc}f_{ab}\frac{\epsilon_{bc}}{\epsilon_{ab}}\mathrm{Re}\left[A^{\alpha}_{ab}A^{\beta}_{bc}A^{\mu}_{ca}\right].\end{align} \tag{ 113 }$$

Therefore we see that while the expression for interband contribution to the sum rule obtained from [9] (equation (112)) is the same as the one derived in our previous paper [10], the one that would be derived from [1, 2] (equation (113)) differs from them. As previously discussed, this is because [1, 2] are missing some terms which only vanish if one assumes a parabolic dispersion which completely neglects the spin–orbit coupling.