Recent progress in the theory of bulk photovoltaic effect

The bulk photovoltaic effect (BPVE), sometimes also called the photogalvanic effect (PGE), refers to the electric current generation in a homogeneous material under light illumination, in contrast to the traditional photovoltaics where a heterojunction, such as a p–n junction, is needed to separate the photo-generated carriers (Fig. 1).^1–4 It has attracted increasing interest among the communities of material science and material engineering due to its potential to surpass the Shockley–Queisser limit governing traditional solar cells,^5,6 and the simplified device geometry due to the homogeneity is also promising for fabricating light detectors. Meanwhile, intense theoretical investigations have been made to understand the physical origin of BPVE,^7,8 which seems to be a rather peculiar phenomenon, as one would naively expect that in the absence of built-in field (as exists in a p–n junction), the oscillating light field would only drive the charge carriers to oscillate periodically without inducing a net current. Thus, the breaking of such intuition implies that we must go beyond the linear response regime to fully account for the BPVE.

To our knowledge, the first investigations of BPVE were conducted in the late 1960s and early 1970s. Experiments were carried out to measure the BPVE photocurrent for several materials that do not have an inversion center, and theories based on time-dependent perturbation theory were formulated for different mechanisms.^9–11 A review of the research work at that time can be found in the book by Sturman and Fridkin.¹ These early works contribute significantly to the understanding of BPVE, where important concepts such as shift current and ballistic current were introduced as possible mechanisms for the DC photocurrent, but the simplified models employed by those works hindered the further interpretation and prediction of this effect for real materials in a broader spectral range. Progress was made toward calculating BPVE in real materials by adopting different electronic structure theories,^12–15 but a truly first-principles calculation with direct comparison with experiments had been lacking. In 2012, Young and Rappe revisited the theories of BPVE and adapted the shift current theory into a formula that is amenable to first-principles calculations.⁴ Good (though not perfect) agreement between the first-principles simulations and the experimental results for BaTiO₃ reinvigorated the theoretical study of BPVE due to the demonstrated predictive power of first-principles calculations.^16,17 Later, first-principles models for other BPVE mechanisms have been reported,^18,19 which, in general, improve the accuracy of the theoretical prediction, and different design routes for enhancing the BPVE in materials are suggested based on the developed first-principles calculations.^20–22 In addition to the exciting progress in numerically calculating BPVE, there is also a renewed theory development that is mainly based on Floquet theory and non-equilibrium Green's function (NEGF) formalism instead of the well-developed perturbation theory.^8,23,24 These newly developed theories allow for the study of BPVE for finite systems and their temporal and spatial behavior. Also, the topological nature of BPVE is explored by several works, which propose that BPVE can serve as a way to probe the topological phase transition.^8,25

In this review, we survey the recent progress in theories and numerical calculations in the field of the bulk photovoltaic effect, aiming to introduce the basic concepts as well as the latest understanding of this effect. The rest of this review is organized as follows: In Sec. II, we will have a brief review of some early and recent experiments measuring BPVE and the characteristics of the observed photocurrent, which are pertinent to the development of the BPVE theory. Then, in Sec. III, we start to discuss the theories proposed for BPVE, including the theory of shift current, ballistic current, and injection current. The recently developed Floquet and NEGF formalisms will also be introduced there. We will present the numerical implementations for each mechanism and talk about the technical details. Following the development of BPVE theory, in Sec. IV, we will go over the strategies of improving BPVE in real materials that are guided by the theory framework introduced in Sec. III. In Sec. V, we will summarize the content introduced in this review and give our perspective on the future development of the bulk photovoltaic effect.

Although there is no consensus on which experiment observed the BPVE for the first time, it is clear that in the 1970s numerous experiments demonstrated the existence of a steady DC in homogeneous materials lacking inversion symmetry under uniform illumination. These observations inspired a plethora of theoretical studies to understand this phenomenon.¹ More recently, there have been experiments trying to clarify the nature of the observed photocurrent, with a focus on trying to separate the contributions from different mechanisms.^26,27 Meanwhile, it has been shown experimentally by Spanier et al.⁶ that the BPVE can indeed surpass the Shockley–Queisser limit. Therefore, to understand the signature and significance of BPVE, we devote this section to reviewing some of the important experiments.

The signatures of BPVE perhaps were detected as early as 1930s, where the photoelectret effect was observed in ferroelectric materials.^29,30 In photoelectrets, the light illumination could induce a long-lasting change in the polarization, and a built-in electric field can be observed. This polarization difference between the ground state and excited states may share the same origin with the BPVE.³¹ Later, a large open-circuit voltage (⁠$Voc$⁠) induced by light was found in SbSI_0.35Br_0.65 and BaTiO₃ between 1970 and 1972,^32,33 which strongly indicated the existence of BPVE in these ferroelectric materials, as traditional photovoltaic effects cannot have $Voc$ larger than the bandgap. In 1974, Glass³⁴ provided the concrete evidence of BPVE by showing a steady-state DC in iron-doped LiNbO₃ and linear scaling with light intensity. Although these early works are undoubtedly crucial for BPVE, the most important experiments in terms of recent theory development can be argued to be the ones conducted by Koch et al. in 1975 and 1976 for BaTiO₃, where the dependence of the photocurrent on light polarization, frequency, and intensity were clearly reported.^28,35 It is this set of experiments that most of recent first-principles simulations compare with and motivate further development in theoretical and numerical methodology.^4,18

To measure the bulk photovoltaic properties of BaTiO₃, several single-crystal samples were fabricated and poled with an electric field about 5 kV/cm while cooling down through the Curie temperature to ensure uniformity. The samples were usually 0.02–0.05 cm thick, 0.1–0.2 cm wide, and 0.1–0.3 cm long. Then, a 488 nm laser focused to 0.03 cm diameter was scanned through the sample from one electrode to the other. What is significant is that the open-circuit voltage was non-zero, larger than the bandgap of BaTiO₃, and almost uniform across the sample, except when quite near to the electrode. This is different from the traditional photovoltaic effect, where the charge carrier separation can only happen at an interface of distinct materials, and the open-circuit voltage is usually smaller than the bandgap. This position-independent behavior was largely unexplored until a recent study²³ (see Sec. III D) in which the real-space distribution of BPVE was simulated explicitly.

More importantly, the spectral dependence of the open-circuit voltage (or equivalently the short-circuit current³⁵) was measured in a range of photon energies above the bandgap of BaTiO₃, as shown in Fig. 2(a). Two orthogonal linear polarizations of light were used, and distinct spectral behaviors were obtained. In particular, the zzZ component (lower-case letters representing the light polarization and upper-case letter representing the current direction) exhibited a sign change around 390 nm. When changing the intensity of light, a linear scaling was found, showing that the observed photocurrent is a second-order response to the electric field. This combination of behaviors is unique to BPVE and cannot be understood by the traditional photovoltaic effect, so any theories trying to explain BPVE are built on these observations.

Since the early experiments were conducted showing the existence of BPVE in non-centrosymmetric (breaking P-symmetry) materials, including the ones by Koch,^28,35 two major complementary theories were proposed to explain the observed photocurrent, namely, the shift current and ballistic current. A detailed introduction and discussion of these can be found in Sec. III, but for the purpose of illustrating the ideas underlying the experiments discussed here, it suffices to know that it is believed that the shift current will be less susceptible to the magnetic field whereas the ballistic current can give rise to a Hall current as any classical charge current does.³⁶ Therefore, to validate the shift current and ballistic current theory, people have designed experiments trying to separate the two types of current with the help of a magnetic electric field.

An outstanding work was conducted by Burger et al.²⁶ in which Bi₁₂GeO₁₂ was chosen as the target material. It belongs to the T space group (#197), which dictates that only the yzX component of the photocurrent (and all permutations of the indices) could have non-vanishing value for linearly polarized light. For circularly polarized light, a similar symmetry argument shows that only xX, yY, and zZ will be non-zero, with a lower-case letter showing the propagation direction of the light. These symmetry properties are particularly useful, since in the absence of magnetic field, a linear light whose polarization lies within the yz-plane or circular light that is propagating along the x-axis could only generate photocurrent flowing along the x-direction, so that a Hall current after turning on a magnetic field along the y-axis can be uniquely identified along the z-direction. If, however, the system already had non-magnetic BPVE current along directions other than the x-axis, then extra effort and caution would have to be taken to separate the Hall current from the “intrinsic” response.

Their procedure for the current separation is as follows (Fig. 3): under circularly polarized light, there is only photocurrent of ballistic type (see Sec. III), so one can extract the mobility of the carriers μ_nth of ballistic current via Hall effect. For linearly polarized light, however, there could exist both shift current and ballistic current, and only the ballistic current is believed to respond to the magnetic field. Thus, one can apply a magnetic field again, obtaining the Hall current under linear light, and then calculate the non-magnetic ballistic current with the help of μ_nth. After subtracting the non-magnetic ballistic current from the total non-magnetic current, the non-magnetic shift current can be acquired, and the contributions from these two mechanisms can, thus, be separated.

Though such a procedure looks reasonable, there are two caveats. For one, this experiment is designed under the assumption that shift current does not respond to the static and uniform magnetic field. This was first argued by Ivchenko³⁶ in which he stated that as long as the cyclotron frequency of the magnetic field is much smaller than the difference between the light frequency and bandgap, then the shift current will be barely impacted by the magnetic field, without further proof. Thus, the validity of this assumption remains to be examined. For another, it is worth mentioning that in the same paper by Ivchenko, the authors demonstrated that the magnetic field can break the time-reversal symmetry and induce a new current, which will be proportional to the non-magnetic shift current for a two-band model. Thus, in another work by Burger et al.,²⁷ they took the new current into account and discussed several different scenarios for the relative magnitude of this current, leaving the exact separation of ballistic current and shift current still an open question. Despite these caveats, this work constitutes an important step forward toward experimentally verifying various BPVE mechanisms.

In addition to the experiments designed to understand the fundamental physics of the BPVE, there is also great research interest toward exploiting the BPVE in real-world applications.^37–39 Since the BPVE is not governed by the rules of traditional photovoltaics, it is, in principle, not limited by the Shockley–Queisser limit that is imposed on traditional solar cells.

The Shockley–Queisser limit explains and quantifies that for high efficiency, the ideal bandgap for traditional solar should not be too large or too small.⁵ If the bandgap is too large, then a large portion of the sunlight spectrum is unable to be absorbed. If the bandgap is too small, then even though more sunlight can be absorbed, the photoexcited electrons will initially occupy higher-energy states but will rapidly thermalize and relax to the conduction band bottom before they can be harvested by the electrodes. Thus, for the solar spectrum, there exists a perfect bandgap that maximizes the power conversion efficiency, and for any specific bandgap value, there exists a maximum power conversion efficiency. However, for the bulk photovoltaic effect, both shift current and ballistic current mechanisms involve non-thermalized carriers giving rise to a current (see Sec. III), so the Shockley–Queisser limit no longer applies.

To demonstrate the capability of the BPVE to achieve high power conversion efficiency, Spanier et al.⁶ employed a tip-enhanced geometry⁴⁰ that can effectively harvest the non-thermalized electrons in BaTiO₃ (Fig. 4). Such geometry makes use of the fact that the electrode tip can screen the polarization bound charge in a very confined region so that it can create a very large electric field around the tip, which will let the non-thermalized electrons ionize more electrons from valence bands and effectively generate an incident photon-to-collected electron (IPCE) efficiency larger than unity. This process should be distinguished from charge separation caused by the band bending at the ferroelectric–electrode interface that falls into the category of traditional photovoltaics. Rather, it manifests an efficient usage of the hot carriers with already asymmetrically distributed momenta that are caused by the BPVE (more specifically, the ballistic current mechanism). As a result, the power conversion efficiency of BaTiO₃ for this geometry is 4.8%, around 50% higher than the Shockley–Queisser limit for materials with a 3.2 eV bandgap, which shows the potential of using BPVE to design next-generation high-efficiency photovoltaic devices.

Over the past 50 years, different theories have been proposed to understand the nature of BPVE, and most of them are based on the time-dependent perturbation theory, either in density matrix form^2,41,42 or in second-quantized form.^3,43 These theories are especially successful in explaining the DC photocurrent for bulk materials under linearly or circularly polarized light, but they do not address the temporal response or account for spatially inhomogeneous light illumination. Thus, in recent years, Floquet theory combined with non-equilibrium Greens function (NEGF) methods have been developed aiming to address these issues.^8,23,24,44 In this section, we will first review the original theories of BPVE based on perturbation theory, and then the Floquet theory and NEGF methods will be discussed at some length to provide perspective on their advantages and shortcomings.

Experimentally, the scaling of BPVE photocurrent with the light intensity is linear,²⁸ so for linearly polarized light, a phenomenological description^1,45 of BPVE can be written as

where j^q is the photocurrent density along the Cartesian direction q, E_r and E_s are the components of the electric field of light, and $σrsq$ is the response tensor or susceptibility tensor that characterizes the BPVE for a certain system. There are some other conventions of labeling the response tensor in the literature, such as $σrsq$ or $σqrs$⁠, where q represents the current propagation direction. In this review, we will stick with the conventions that q always appears in the superscript or in the last place. As $ErEs$ is proportional to the light intensity, this expression can correctly describe the scaling behavior of BPVE with light intensity. It is noted that such second-order response has already imposed the symmetry constraint that only non-centrosymmetric (no inversion center) structures can have BPVE. To see this, imagine applying an inversion operation to the system. Polar vectors, such as j^q, E_r, and E_s, will acquire a minus sign,

If the system possesses inversion symmetry, then the response tensor $σrsq$⁠, an intrinsic property of the material, will return to itself after the inversion operation. As a result, inversion symmetry leads to

which indicates that j^q will always be zero for a centrosymmetric structure. Therefore, breaking inversion symmetry is a prerequisite for BPVE, and accordingly BPVE can be used in detecting phase transitions involving inversion-symmetry breaking.⁴⁶ 

Furthermore, Eq. (1) can help predict the form of $σrsq$ if there exist some crystal symmetries. The analysis is essentially an extension of what we have done above, i.e., the response tensor should not only be invariant under a symmetry operation belonging to the symmetry group of the crystal, but it should also be a linear combination of other tensor components after the symmetry operation.⁴⁷ One can use this equality to determine nonzero and independent tensor components. As an example, we can analyze the form of $σrsq$ for tetragonal PbTiO₃. Tetragonal PbTiO₃ belongs to the $C4v$ group, so it contains the identity, two fourfold rotations, one twofold rotation, and four mirror operations. If we denote the symmetry operation as $a↔$⁠, then under the symmetry operation, the transformed operator will be $(σrsq)′=|a↔|∑ijkariasjaqkσijk=σrsq$⁠. For the C₂ rotation about the z axis, the $σxxy$ will transform as

Thus, by symmetry, we already know that $σxxy$ would always be zero. By repeating the same process for all symmetry operations in $C4v$ and all tensor components, we will find that the $σrsq$ should take the following form (in Voigt notation):

We can see that symmetry analysis is very powerful to give us rich information about the BPVE response tensor even before any numerical calculation is done. More interestingly, if any tensor component that should be zero by symmetry is found to be nonzero in experiments or numerical simulations, then it reveals that some symmetry breaking processes must have occurred.⁴⁸ 

To develop a microscopic theory for $σrsq$⁠, one can start with a non-interacting many-body system where the two-body interaction is effectively treated in a mean-field fashion, a strategy that is widely used in modern electronic structure calculations such as density functional theory (DFT) and the Hartree–Fock approximation.⁴⁹ Then, one is interested in how the equilibrium density matrix of this system will evolve under the perturbation from light. We would especially like to know the form of the resulting non-equilibrium steady-state density matrix. More concretely, under the dipole approximation, the electron–light interaction in the velocity gauge^50,51 can be expressed as the following minimal coupling form:

and the full Hamiltonian can be written as

Here, $v̂$ is the velocity operator and $Â$ is the vector potential of light, which can be rewritten as $Â(t)=−iωÊ(t)$ in the velocity gauge.⁵¹ $Ĥ0$ describes the non-interacting Hamiltonian with known energy spectrum $εn$ and eigenstates $|ψn⟩$⁠. Then, the equilibrium (non-perturbed) density matrix (operator) can be constructed as

where p_i is the probability of being in the many-body state i. We would like to know the steady-state density matrix $ρ̂(t)$ under continuous illumination because the steady-state current can be computed as

Note that the form of Eq. (6) is written in the interaction picture, and in this picture, $ρ̂I(t)$ can be calculated perturbatively as

where $Ĥe−lightI(t)$ is the perturbation in interaction picture.⁵¹ As in BPVE theory, the response is second-order, we only retain the third term in Eq. (7), and then compute the current j^q via Eq. (6). After a certain amount of algebra, the steady-state current can be explicitly written as

where $flk$ is the Fermi–Dirac distribution function, $vnlr(k)≡⟨nk|v̂r|lk⟩$ is the velocity matrix, and η is an infinitesimally small value (⁠$0+$⁠) appearing in the adiabatic turning-on $eηt$⁠.⁴⁹ Note that we are considering a perfect crystal in the thermodynamic limit, so the dependence of the eigenstates on the crystal momentum k has been made explicit here.

Equation (8) is one central result for the BPVE theory as it expresses the steady-state current response tensor with quantities that can be obtained from numerical models such as quadratic band structure models, tight binding models or, from first-principles calculations. It is, therefore, tempting to conduct numerical calculations of BPVE based on this expression. However, such calculations would be cumbersome, due to the summation over band index m. A closer inspection of Eq. (8) will reveal that there is no selection rule for m, meaning that, in principle, one should include an infinite number of bands when summing over m. In practice, even though the number of bands in the summation would always be truncated, the long tail due to the function of form $1/(εnk−εmk−2iη)$ will still require a very large number of bands for converged results, which would cause formidable computational cost. Thus, most numerical calculations will not directly use Eq. (8), but instead employ some further simplified forms.

To simplify Eq. (8), we will split it into two contributions: the “three-band” contribution where $n≠m$ corresponding to the off diagonal part of $ρ̂(t)$⁠, and “two-band” contribution where n = m, corresponding to the diagonal part of $ρ̂(t)$⁠. It turns out that these two contributions will appear under different conditions and thus carry distinct physical meanings.

We will first focus on the three-band contribution, which has a more well-known name, shift current. The reason why it is called “shift” current will become clear later. After imposing the condition that $n≠m$⁠, the summation over m can be carried out analytically so as to avoid the necessity of including a large number of bands in numerical calculations. The general procedure is to make use of the identity,

and then replace $vmnq(k)$ in Eq. (8). The rationale for why it can enable the analytical summation of m is that the Hamiltonian $Ĥ0$ in the commutator will give a term $(εnk−εmk)$⁠, which can exactly cancel the principal part in $1/(εnk−εmk−iη)$⁠. Care must be taken, however, when summing over m after making this substitution, because now this summation indeed includes the terms n = m. Thus, we need to manually exclude the terms involving $⟨nk|[Ĥ0,x̂]|nk⟩$⁠. Now, with the help of the expression for position operator $x̂$ in a periodic system by Blount,⁵² 

where $unk$ is the lattice periodic part of the eigenfunction of $Ĥ0, ψnk(x)=eik·xunk(x)$⁠, we can finally rewrite the three-band contribution as

Here, $χnq(k)≡⟨unk|i∇kunk⟩$ is the Berry connection, and $ϕnlr(k,k)$ is the phase of $⟨nk|v̂r|lk⟩$⁠. This expression is composed of two parts,

with

which can regarded as the k-resolved transition rate, and

which has a unit of length and can be regarded as the coordinate shift of carriers in real-space during the transition. It is for this reason that Eq. (14) is named as shift vector, and the three-band contribution equation (11) is usually called shift current. Note that Eq. (11) is now in a two-band form having only n and l, but in essence it is still a three-band expression as the summation of the m is encapsulated in the k-derivative terms. In other words, the first-order expansion of the k-derivative terms involves another summation over all states.² 

Shift current has many interesting properties that are distinguished from classical charge currents. For one, it is independent of carrier lifetime and is robust against the scattering by disorder.^2,8 It is not a current carried by classical moving particles as it is exclusively from the coherence of the density matrix, which has no interpretation in the classical picture. Instead, it is a manifestation of wave-packet evolution when transitions between different electronic states are happening. For another, it contains quantum information (the so-called geometrical information) of the electronic structure, as the phases of the wave functions are considered explicitly in the shift vector $Rrq(n,l;k)$⁠, whereas for classical charge carriers, only group velocities (diagonal elements of the velocity matrix) and occupations (diagonal elements of the density matrix) are relevant. Thus, its quantum nature has attracted a great deal of attention, and its connection to the modern theory of polarization and topological materials have been explored due to their common relation to Berry connection.⁵³ 

It is now feasible to compute shift current for real materials reliably via first-principles calculations, providing a possible route to quantify its contribution to the experimentally observed photocurrent.^{4,16,17,54,55} Natos and Sipe demonstrated that the shift current can be calculated from first-principles theories,¹³ and Young and Rappe revolutionized this field by showing that the first-principles prediction of shift current from density functional theory (DFT) can be directly compared to experiments.⁴ The latter formalism bears the caveat that the numerical differentiation of wave functions with respect to k might break the gauge invariance (global phase of wave functions) of the shift vector equation (14). Inspired by the strategy employed in the modern theory of polarization,^56,57 the gauge invariance is preserved by transforming the direct derivative into a logarithmic derivative, and the shift currents of tetragonal BaTiO₃ and PbTiO₃ were thus computed using DFT.

Figure 5(a) shows the $σxxz$ (gray) and $σzzz$ (black) components of the current response tensor. It is clear that for both PbTiO₃ and BaTiO₃ the largest response is at frequencies much larger than the band edge, which cannot be captured by simple model calculations which usually only consider energy regions around the band edges.^1,58,59 Also, the shift current calculated from first principles demonstrates large polarization dependence, where $σxxz$ and $σzzz$ differ not only in magnitude, but also in sign, which is consistent with experimental observations. To make quantitative comparison with experiments by Koch et al.,²⁸ Young and Rappe⁴ also calculated the total shift current flowing through the system by taking into account the absorption of light and the sample dimensions,² 

where J^q is the total current, $αrs$ is the absorption coefficient characterizing how much light can be absorbed and how deep the electric field of light can penetrate, and w is the sample width. In the experiments,^28,35 the irradiation intensity is $0.35–0.6 mW/cm2$⁠, from which the electric field can be deduced, and the sample width is $0.1–0.2 cm$⁠. Combined with the theoretical absorption coefficient $αrs$⁠, a quantity that is readily evaluated from first principles in the form of Fermi's golden rule,⁶⁰ the total shift current can be computed and compared against the experimental photocurrent, as shown in Fig. 5(b). What is remarkable is that despite a small mismatch of the frequencies, the calculated shift current of BaTiO₃ can reproduce all the salient features at the band edge, including the overall magnitude, line shape as well as sign reversal. Thus, Young and Rappe inferred that the main contribution to BPVE is shift current, at least for BaTiO₃. In addition, shift current is predicted as a mechanism to generate pure spin current (PSC), the first proposal to apply BPVE in spintronic devices. This will be further discussed in Sec. III A 3.

However, a follow-up study by Fei et al.⁶¹ shows that the contribution of shift current to BPVE might be exaggerated since the absorption coefficient computed from single-particle approximation (assumed by DFT) will be underestimated. After improving the band structure with the GW approximation and introducing the exciton correction to the absorption coefficient,⁶² they found that the calculated shift current will be scaled down such that there is a larger discrepancy between the experimental and theoretical spectra (Fig. 6). This shows that in addition to the shift current, other mechanisms could also participate in the generation of photocurrent. Indeed, as stated earlier, the shift current only originates from the off diagonal elements of the density operator $ρ̂(t)$⁠, but the contribution from the diagonal part, the two-band contribution, has not been considered. Therefore, it remains to examine the contribution from the diagonal part and whether it will improve the BPVE theory.

Under linearly polarized light, the two-band contribution can be obtained from Eq. (8) by imposing the condition that n = m,

If the band structure possesses time-reversal symmetry, which is the case for nonmagnetic materials, then it can show that the three-velocity term $vnlrvlnsvnnq$ will undergo a sign reversal for $−k$⁠: $vnlr(−k)vlns(−k)vnnq(−k)=−vlnr(k)vnls(k)vnnq(k)$⁠, In the meantime, the Fermi–Dirac function and the delta function will be even for k and $−k$⁠. As a result, when considering the response to a linearly polarized light where the rs and sr component cannot be distinguished, the integration of k over the Brillouin zone in Eq. (16) will be exactly zero, meaning that no contribution will exist for the diagonal part of the density operator. (For magnetic systems, $jq,diag$ will no longer vanish and is referred to as injection current, which will be discussed in more detail in Subsection III A 3.)

However, $jq,diag$ is nonzero when additional scattering processes are present. To see this, we formally rewrite Eq. (16) into a form whose physical meaning is manifest,

where we have explicitly considered the transition from the valence band v to conduction band c in a semiconductor due to light excitation. The minus sign of $vvkq$ comes from $fck−fvk=0−1=−1$⁠. $Γcv,krs(ω)$ is the carrier generation rate that contains the transition intensity $vnlr(k)vlns(k)$ and the energy selection rule $δ(εnk−εlk−ℏΩ)$⁠. Note that we have discretized the integration in the Brillouin zone by a summation over k points. This is simply the expression for current in the framework of the Boltzmann transport equation, in which the current is equal to the carrier velocity multiplied by its distribution function,¹⁸ and it is expected that without any other interaction, $Γcv,−krs(ω)=Γcv,krs(ω)$⁠. If we include additional interactions when computing the carrier generation rate, that is, we extend the Fermi's golden rule to higher orders, then it is likely that $Γcv,−krs$ is no longer equal to $Γcv,krs$ in a non-centrosymmetric system, and we call the current from the asymmetric carrier generation ballistic current. One should not confuse ballistic current with ballistic transport⁶³ as they carry distinct but related meanings. For ballistic transport, it means that the carriers can flow for a certain length without any scatterings, whereas the ballistic current can only exist in the presence of coherent scatterings during the optical excitation which induce the population asymmetry (the flow of carriers after the optical excitation will encounter no scatterings for a period τ₀, which is similar to ballistic transport).

To systematically investigate the effect of interaction on the carrier generation rate, we can express the overall generation rate $Γrs$ in terms of the velocity–velocity (current–current) correlation function χ^rs, which essentially counts the number of excited electron/holes by assuming that each absorbed photon will generate an electron–hole pair,

where the real-time retarded correlation function χ^rs can be obtained from the corresponding imaginary-time correlation function $χTrs$ via the analytical continuation $χrs(ω)=χTrs(iωn→ω+i0+)$⁠, where $0+$ is a infinitesimal positive number.^51,64 In this approach, the interaction effect can be included in the correlation function $χTrs$⁠,

from which we can see that the overall carrier generation rate can be decomposed into $k–$resolved carrier generation rate $Γrs(ω)=∑cvkΓcv,krs(ω)$⁠, and it can be evaluated perturbatively with respect to different interactions.

Various processes could lead to asymmetric $Γcv,krs$⁠, such as electron–phonon interaction,¹⁸ electron–hole interaction,¹⁹ and scattering from defects. Among these, the electron–phonon interaction and electron–hole interaction are of most interest because they are intrinsic to any semiconductor, regardless of the quality of the crystal. Thus, most work investigating ballistic current will focus on these two interactions. We would like to have a few more words about these scattering processes being intrinsic as some people would instead regard ballistic current as an extrinsic mechanism for BPVE. Such claim comes from the comparison with shift current where only a perfect and static lattice is considered, so shift current is considered as an intrinsic mechanism, and ballistic current is classified as extrinsic due to the participation of additional processes. However, this classification will be misleading since any realistic perfect material would have lattice vibration and Coulomb interaction, so we think that ballistic current should also be intrinsic if intrinsic scattering processes are considered. Nevertheless, different authors could have different philosophies for this classification, and readers should be careful about what they mean by intrinsic and extrinsic.

The ab initio calculations of ballistic current were realized by Dai et al., where the electron–phonon mechanism (termed as phonon-ballistic current) was taken into account.¹⁸ By treating the electron–phonon interaction via the Frölich e–ph Hamiltonian,

where $Φ̂q̃μ=âq̃μ+â−q̃μ†$ is the phonon field operator, $âq̃μ$(⁠$âq̃μ†$⁠) are the phonon annihilation (creation) operators, and $gμkk′nn′$ is the electron–phonon coupling matrix, the perturbative expansion of Eq. (19) can be performed with the help of Feynman diagrammatic technique where each perturbative term $ΔΓcv,krs(ω)$ can be represented by a connected diagram. For the lowest order non-zero terms, there are three diagrams involved in the optical transition as shown in Figs. 7(a)–7(c), and it can be proved that only Fig. 7(a) will lead to an asymmetric carrier generation. Evaluating this term with Feynman rules, performing the analytical continuation, and using Eq. (19), the asymmetric part of the carrier generation rate $Γcv,krs,asym(ω)=12(ΔΓcv,krs(ω)−ΔΓcv,−krs(ω))$ can be expressed in terms of velocity matrices $vnlr(k)$⁠, electron–phonon coupling matrices $gμkk′nn′$⁠, and band energies $εnk$⁠. The complete form of $Γcv,krs,asym(ω)$ can be found in the work,¹⁸ and in combination with Eq. (17), the phonon ballistic current can be computed from first-principles calculations.

Similarly, in another work by us,¹⁹ the electron–hole interaction (named as exciton ballistic current) is considered on the same footing as the electron–phonon interaction. However, unlike the electron–phonon interaction where it suffices to keep only lowest order terms, the long-range character of the Coulomb interaction will require in principle infinite orders of terms in the perturbative expansion. Luckily, most diagrams can be shown not to contribute to the asymmetric scattering, and the sum of infinite orders of ladder diagrams, the ones involved in asymmetric scattering, can be done exactly. A certain amount of algebra will lead to rather simple results for the sum of ladder diagrams,

where

Here, $Vq̃cc′,vv′$ is the screened Coulomb interaction in the basis of eigenstates of H₀, and $q̃$ is the Fourier component of the Coulomb interaction.⁶⁵ Equation (22) can be solved numerically to yield $ṽcv,k(ω)$⁠, with which one can calculate the ballistic current from electron–hole interaction from Eqs. (16), (18), and (21). A different approach to computing exciton ballistic current is also presented in the same work¹⁹ where the Bethe–Salpeter equation is solved, from which the carrier generation rate can be computed from the exciton wave functions.

According to the procedures prescribed above, the phonon and exciton ballistic currents are calculated for BaTiO₃ and can be found in Figs. 8 and 9. Note that in Fig. 8 we are plotting the total current [as in Eq. (15)] whereas in Fig. 9 we are plotting the response tensor $σxxz$ and $σzzz$⁠. Clearly, the discrepancy between the experimental photocurrent and theoretical shift current for the xxZ component can be partially filled by the phonon ballistic current, but for the zzZ component where the shift current is already in good agreement with experiments, the phonon ballistic current barely changes the theoretical photocurrent. This shows that in addition to shift current, phonon ballistic current is also a major mechanism in BPVE, and it remains to check whether the exciton ballistic current can further improve the BPVE theory. Unfortunately, Fig. 9 shows that exciton ballistic current can be two orders smaller than the phonon counterpart, and the smallness makes it hard to connect the features found in the diagrammatic approach with those in the many-body approach. Thus, even though we have included infinite orders of Coulomb interaction when computing the asymmetric generation rate [Eq. (18)], the canceling among the diagrams makes its overall contribution much smaller than that from the electron–phonon interaction, where only the lowest-order diagram is taken into account. A similar calculation for monolayer MoS₂ shows similar insights.¹⁹ To summarize, when evaluating the intrinsic contribution from ballistic current, it is usually safe to only consider the electron–phonon interaction. The agreement between theory and experiment at the present level of theory is incomplete, so to further improve BPVE theory, scatterings from other sources, such as defects, should be included when computing the asymmetric carrier generation rate.

As alluded above, for nonmagnetic (time-reversal-symmetric, T-symmetry) systems, the injection current will vanish under linearly polarized light.^66–70 As a result, most theoretical and experimental work regarding injection current is centered around the circularly polarized light, which has a slightly different expression from Eq. (16).^46,71,72 We will discuss circular photogalvanic effect (CPGE) further in Sec. III B. Nonetheless, in recent years, more 2D magnetic materials have been discovered that inspired a renewed interest in their electronic and optical properties, especially in their photovoltaic effect. As T-symmetry is usually broken in magnetic systems, the symmetry argument about the parity of $vnlr(k)$ and $εnk$ no longer applies, which brings about the possibility of observing the injection current even under the linearly polarized light.

One example of the 2D magnetic materials attracting attention is CrI₃, a ferromagnetic insulator.^67,68 In the bilayer case, it can exhibit two phases: ferromagnetic (FM) and antiferromagnetic (AFM) (Fig. 10). The latter will break both the inversion symmetry and T-symmetry, causing an asymmetric band structure at k and $−k$⁠. Furthermore, the band velocities at k and $−k$ will not cancel so that there is no requirement that the carrier generation rates at opposite k points will be equal. Therefore, by using Eq. (16), the injection current of CrI₃ can be calculated [Fig. 10(c)]. A similar investigation of MnBi₂Te₄ has also been done,⁷⁰ demonstrating a giant injection current (two orders higher than the shift current of PbTiO₃ and BaTiO₃.)

Two things to note about these calculations: (1) The AFM phase of CrI₃ is special in that its centrosymmetry is broken by spins. So, an inversion operation about the interlayer inversion center will keep the lattice the same but reverse the spin directions. Thus, the spin-orbital coupling (SOC) is required to make sure that the AFM and reverse-AFM will correspond to different energies. Otherwise, neglecting the SOC will make the band structure still symmetric for k and $−k$⁠.^67,68 (2) The relaxation times τ₀ used in these works^68,70 are obtained from experimental values of materials belonging to the same family and are somewhat arbitrary. So, the large injection current observed in these calculations is partly due to the large relaxation time. A more consistent treatment for relaxation time would be calculating τ₀ from first principles as is done by Dai et al.¹⁸ when computing the ballistic current. Their calculations show that the constant relaxation time approximation is reasonable, showing weak dependence on band indices and crystal momenta, but the value would differ from material to material. For example, the computed momentum relaxation time of BaTiO₃ is 2 fs, compared with 100 and 600 fs used in MnBi₂Te₄ and CrI₃, respectively. Thus, when interpreting the magnitude of injection current calculated with the constant relaxation time approximation, one has to be careful about the choice of the τ₀, and this consideration also holds for circular injection current (Sec. III B).

Another scenario where the linear injection current could be important is the photo-generation of pure spin current (PSC).^70,73 In this case, T-symmetry-breaking is not needed because there is no charge current. What is required for linear injection PSC, however, is a sizable SOC which will make the band structure at k and $−k$ have different spin characters, as shown in Fig. 11. If we use Eq. (16) to compute the spin-up and spin-down injection currents separately, the band structure of each spin is asymmetric even though the overall band structure is symmetric. Then, spin-up carriers will have a net current whose direction is opposite to spin-down carriers, generating a pure spin current. One can immediately see why strong SOC is a prerequisite for injection PSC as for systems with no or weak SOC, the spin-polarized band structure will become symmetric again. Therefore, it is expected to observe large PSC in materials having heavy elements, such as CdS, SnTe, and transition metal dichalcogenides (TMD). On the hand, in terms of generating PSC, shift current is also predicted to be a viable mechanism.⁷⁴ In contrast to the linear injection current, the shift current mechanism does not necessarily require SOC. Instead, it can exist in antiferromagnetic systems. Take the PT-symmetric (breaking P- and T-symmetry individually but preserving the PT as a whole) hematite Fe₂O₃ as an example, the PT-symmetry will make shift vectors at the same k point have opposite directions for spin-up and spin-down components despite the same transition rate, thus giving rise to a pure spin current. The general symmetry requirements for using shift current to generate PSC has been identified in the work,⁷⁴ and its existence is demonstrated by performing first-principles calculations on several antiferromagnetic materials, such as Fe₂O₃ and BiFeO₃.

Having discussed several major mechanisms for BPVE under linearly polarized light, it could be inspiring to unify them in the perspective of a kinetic model. A kinetic model takes into account various processes (excitation, recombination, etc.) contributing to the time-evolution of the density matrix, including the diagonal (occupation) and off diagonal (coherence) matrix elements. In principle, it is able to describe temporal, steady-state, and equilibrium time evolution, while most experiments measure observables in steady-state or equilibrium states, in which the density matrix elements should possess stable values and can, thus, be used to compute these observables. Therefore, to study the steady-state DC due to BPVE, it is desirable to find the influence of all relevant photoinduced and thermal processes and connect them in a quantum Liouville equation in order to establish the steady-state values in a kinetic model.

This kinetic model was originally conceived by Belinicher et al., where several important processes were considered, including light excitation, electron–phonon coupling, and defect scattering.⁹ As pointed out by Belinicher et al., the foundational idea of the kinetic model is that the time-evolution of the density operator can be described by the quantum Liouville equation (written in the Schrödinger picture),

where $He−light$ is quantized (the vector potential is expressed as photon creation and annihilation operators) and drives the excitation (including stimulated recombination if T > 0) and spontaneous recombination, while $Ĥe−phonon$ and $Ĥe−defect$ are responsible for the thermalization in the full kinetic cycle, and they can also participate in the excitation and recombination processes as well in the full kinetic cycle. The quantum treatment of light can enable the spontaneous recombination because it will allow the electronic system to be coupled to the light that is at every possible frequency, so there is a driving force for electron and hole to recombine even though their energy difference is different from that of the incident light. The electrical current calculated from $j=e/mTr[ρ̂ĵ]$ can, thus, be categorized as excitation, thermalization, and recombination current according to the processes participating in the current generation.

The progress reported in this review can be described as developing first-principles computational approaches that provide some terms in the kinetic model. The rest of these terms can be approximated, with sensible functional forms. For example, in deriving the shift current, we are essentially computing the current generated in the excitation process, so for this purpose, we take $Ĥe−light$ as composed of only monochromatic light (which is equivalent to taking it as a classical monochromatic electromagnetic field), and we approximate the thermalization and spontaneous recombination process related to $Ĥe−phonon$ and $Ĥe−defect$ by the constant-relaxation-time approximation,

with $ℏ/τ0=η$⁠. The last term concerns the dissipation that takes into account the thermalization, which would otherwise be taken care of by $Ĥe−phonon$ and $Ĥe−defect$⁠, and the recombination (the carriers no longer recombine spontaneously through $Ĥe−light$ since it now represents a classical field and cannot absorbed the emitted photons). To compute shift current, we consider how the off diagonal elements of $ρ̂(t)$ evolve according to Eq. (24), which will lead to Eqs. (8) and (11). On the other hand, if we are interested in computing the ballistic current, which is from the diagonal part of $ρ̂(t)$ during the excitation process, then we need to include $Ĥe−phonon$ and $ĤCoulomb$ in Eq. (24) and consider their contributions to the diagonal elements of $ρ̂(t)$only at the excitation process. To be more specific, at the moment of the optical excitation, the scatterings from $Ĥe−phonon$ and $ĤCoulomb$ will interfere with the scatterings from $Ĥe−light$ and give rise to the phonon and exciton ballistic current, while the thermalization and recombination processes⁷⁵ are still approximated by $iℏ[ρ̂(t)−ρ̂0]/τ0$⁠.

Some additional types of bulk photovoltaic current originating from the more general expression equation (23) have been formulated already,⁹ while their reformulation into first-principles calculations are still ongoing. It has been shown by Belinicher et al. that in addition to the excitation shift current, there also exist real-space shift currents associated with the thermalization and recombination. The recombination shift current is easily evaluated, since its form should be exactly the same as the excitation shift current equation (11) except that the distribution functions will be replaced by the non-equilibrium steady-state distribution. Hence, it is mostly concentrated at the band edge states with sign opposite to the excitation shift current. The thermalization shift current, or phonon shift current, accompanies the electron–phonon scattering processes after the optical excitation. The related shift vector is similar to Eq. (14) with the phase being changed to the phase of the electron–phonon coupling matrix elements.⁹ This contribution could be important, since the phonon ballistic current plays an important role in BaTiO_3.^18,19 Other possibilities also exist, for example, in the spontaneous recombination process if electron–phonon and electron–hole interactions are taken into account, so there are still numerous opportunities in this regard. Being able to evaluate the kinetic model from first principles will greatly enhance the accuracy and the predictive power of the BPVE theory, and this kinetic model devised for steady-state could be potentially extended to compute the temporal evolution of the density matrix in order to study ultrafast experiments.

At this point, readers who are familiar with the anomalous Hall effect (AHE) may have noticed that the shift current and ballistic current in BPVE have direct parallels with the side jump and skew scattering mechanisms in AHE.⁷⁶ It is interesting however that no comparison has been made for the two phenomena; in this section, we make the connection explicit. AHE in itself is a vast topic and includes many different aspects for its theory development, so this section is by no means comprehensive. Readers are encouraged to read the review by Nagaosa et al.⁷⁶ to learn more about AHE.

AHE is the phenomenon that when measuring the Hall effect in a ferromagnetic metal, the Hall current deviates from the Lorentz law and is usually very large. From this description, it is clear that AHE usually refers to the linear transverse response to static electric field, and for it to happen, breaking T-symmetry is required. For BPVE, however, it is the second-order response to the oscillating electric field (light), and it is not restricted to the current generation in a direction transverse to the electric field of light. Rather than breaking T-symmetry, breaking inversion symmetry is required for BPVE, which is due to the characteristics of second-order response as discussed in the beginning of Sec. III A. So, it can certainly be seen that AHE and BPVE describe two distinct phenomena.

However, when formulating the theories for AHE and BPVE, one can often find that the ideas developed in AHE can be borrowed to understand BPVE. There are two extrinsic mechanisms induced by defects in AHE, namely, skew-scattering and side-jump. It is generally accepted that skew-scattering means the breaking of detailed balance (transition rate of $k→−k$ no longer equals to that of $−k→k$⁠) in the presence of scattering processes (mostly scattering from magnetic impurities or disorders) with strong SOC and breaking of T-symmetry, which makes that the carrier generation rate have a preferred direction and thus asymmetric. This is very similar to the idea rooted in the ballistic current for BPVE, but in BPVE the strong SOC and breaking of T-symmetry are not required, and scattering processes usually refer to electron–phonon scattering and electron–hole scattering.

The side-jump mechanism describes the coordinate shift of the electron wave packet when scattered by a magnetic impurity with SOC, which resembles the shift current mechanism for BPVE. Indeed, the expression for coordinate shift in side-jump mechanism is exactly same as the shift vector in shift current mechanism except that the transition rate is now governed by the impurity scattering in AHE instead of the optical excitation in BPVE.⁷⁷ One can try to replace the transition rate in side-jump with the optical transition rate (transition dipole matrix) and then recover the shift current expression in a semiclassical treatment of AHE. Thus, it is expected that the further development of BPVE theory could continue to be inspired by the better-understood AHE phenomenon. Conversely, a novel nonlinear AHE effect due to the so-called Berry curvature dipole is predicted in non-magnetic but non-centrosymmetric systems, which clearly has its inspiration from the theory of BPVE.⁷⁸ 

Gauge invariance is an important principle in the context of nonlinear optics. Although it is physically imposed that any observable must be gauge-invariant, it is sometimes not so easy to identify the equivalence between formulas derived from two different gauges. In the case of shift current, for instance, we started with the velocity gauge [$Ĥ(r̂,p̂)→Ĥ(r̂,p̂−eA)$] and arrived at Eq. (8) (with the constraint that $m≠n$ to describe the shift current⁴¹) while Eq. (11) is the expression that can be directly obtained from the length gauge³ [$Ĥ(r̂,p̂)→Ĥ(r̂,p̂)+eE·r̂$]. The process going from Eq. (8) to Eq. (11) via Eq. (9) and Eq. (10) can be regarded as proof that these two gauges will give equivalent results,^2,43 but it is not obvious at all why Eq. (8) and Eq. (11) are equivalent solely by comparing the two expressions. Therefore, when obtaining analytical expressions, it is always advisable to select a gauge at convenience in the beginning and stick to it throughout the derivation.

Despite the formal equivalence between velocity-gauge and length-gauge expressions [Eqs. (8) and (11)], in practical calculations they do not yield identical results. This is because in Eq. (8), the summation of intermediate states m runs over an infinite set of states, while in practice some kind of truncation has to be made. This truncation breaks the equivalence between two gauges, and the influence can be drastic. For example, for a two-band model, using Eq. (8) will lead to zero shift current regardless of the model details, whereas Eq. (11) can give nonzero shift current for inversion-breaking systems. In addition, velocity gauge expressions within the truncated space will give artificial divergence as the frequency ω goes to zero even for a wide-gap system, while eliminating this artificial divergence using some sum-rule is again to convert the velocity-gauge expression to length-gauge expression.⁴³ The length gauge equation (11) does not suffer from the truncation problem since the Dirac delta function indicates that only bands relevant to the optical excitation will contribute to the shift current. Therefore, generally speaking, length-gauge will give more accurate results when it comes to numerical calculations.

That being said, velocity gauge also has several advantages over the length gauge. Conceptually, velocity-gauge expressions have clearer, which may render clearer physical interpretation.^43,79 Numerically, as the length-gauge usually requires the derivatives of some wave functions with respect to crystal momentum k, computing them is more challenging as a continuous gauge in k (wave function phase) is required, but not always guaranteed or provided by numerical calculations. So, clever manipulation of Eq. (11) is needed to make it insensitive to the smoothness of wave functions with k (e.g., Wilson loop methods, see Refs. 4 and 80). In the velocity-gauge expression equation (8), however, the required quantities are band energies and velocity matrices, which are easily obtained from standard ab initio calculations or tight-binding models. It has been demonstrated that including more bands in the truncated space can both efficiently suppress the artificial divergence as $ω→0$ and approach length-gauge results.⁸¹ Thus, the velocity gauge (when well converged) is simpler for numerical implementation.

To summarize, it needs to be emphasized that gauge-invariance is ensured in nonlinear optical effects, but it could be a difficult task to prove the gauge-invariance, and the invariance could be broken numerically. For more accurate results, we recommend using the formulas obtained from length-gauge or velocity-gauge with some sum-rules.

After finishing the discussion of the basic aspects of the BPVE under linearly polarized light, it is natural to examine the photo-response under the circularly polarized light. This phenomenon has a more well-known name in the area of spectroscopy, the circular photogalvanic effect (CPGE).^46,71,72 So, in this section, we refer to the circular BPVE as CPGE to be consistent with the existing literature. Different from the linear BPVE equation (1), the phenomenological description of CPGE can be written as

Here, γ_ql is the response tensor for CPGE, q is again the Cartesian direction of the photocurrent, and l is the propagation direction of the circularly polarized light. Note that E is vector form of the electric field, so for circularly polarized light, it has the form $E=(Ex, e±iπ2Ey, 0)$ with $Ex=Ey=E0$ [assuming the light is propagating along z, and +(−) representing left (right) helicity]. Then, $i[E×E*]l=±(ErEs−EsEr)$⁠, so the second-order response must differentiate between $ErEs$ and $EsEr$ in order to have non-zero response.

Now looking at the ballistic current equation (17) and the more explicit expressions for the carrier generation rate,¹⁸ one can find that there is no requirement that the rs component has to be equal to sr, so the ballistic current from the electron–phonon interaction and electron–hole interaction can also appear for circularly polarized light. More interestingly, the intrinsic diagonal contribution (in the absence of extra scattering processes) will be non-zero as well for the same reason, and this non-vanishing injection current is actually what people usually refer to as the CPGE. For shift current, however, the derivation of Eq. (11) has already symmetrized the rs and sr components due to the consideration that linearly polarized light cannot distinguish rs and sr. Therefore, Eq. (11) cannot be used directly to analyze the shift current under circularly polarized light, and different expressions need to be worked out in this case. Below, we examine the injection current and shift current more closely for left-handed circularly polarized light, and the results for the right helicity can be obtained with a sign reversal.

The derivation of shift current under circularly polarized light is largely the same as that under linearly polarized light, except that now different components of the vector field will have different phases by $π/2$⁠: $ÂCP(t)=Ârer cos ωt+Âses sin ωt$⁠, where $er$ and $es$ are the unit vectors representing the polarization directions. Then, the perturbation will become

and the density operator equation (7) will be expanded under the new perturbation, and the off diagonal contribution (responsible for shift current) can be separated by requiring $n≠m$ in Eq. (6). Following exactly the same procedure as for linear shift current, we can arrive at an important expression for the off diagonal contribution,⁸² 

where $p$ represents the principal part integration and all the other symbols carry the same meaning as those in Eq. (11).

With T-symmetry, the terms multiplied by the delta function $δ(εnk−εlk+ℏΩ)$ can be shown to vanish by considering the parities at k and $−k$ as well as by switching the dummy indices n, l, and Ω. However, the terms multiplied by the principal part $p1εnk−εlk+ℏΩ$ will survive, giving rise to a mysterious non-resonant (sub-bandgap) response, which we will denote as $jppq,CP−sh$⁠. One is therefore attempted to conclude that there exists sub-bandgap shift current for circularly polarized light. However, there is no concrete experimental result showing the existence of such a non-resonant shift current, so this term is apparently unphysical and has to be canceled by some arguments or by other unknown terms. Fortunately, we can indeed identify such a term in the diagonal contribution of density operator $ρ̂I$⁠, which makes the shift current vanish under circularly polarized light for systems possessing T-symmetry, and it is the combination of $jppq,CP−sh$ with this term from $ρnnI$ that is responsible for a distinct additional contribution in metals (see Sec. III B 3).

The diagonal contribution $ρnnI$ in Eq. (6) can be derived in a similar fashion as the off diagonal contribution, but under circularly polarized light, it is worth pausing at an intermediate step and examining the derivation,

Apparently, when the infinitesimal $η→0$⁠, Eq. (28) will diverge (ballistic current will also diverge if there is no relaxation, $τ0→∞$⁠). However, in the constant relaxation-time approximation $η=ℏ/τ0$ where τ₀ is the relaxation time, as long as the excited carriers are scattered (which is always the case in real materials), the τ₀ will be a finite value and thus the injection current will not diverge. On the other hand, even in the clean limit where $τ0→∞$⁠, the injection rate$djq,CP−diag/dτ0$ will remain finite and constant; this is why people name this current “injection current” in the first place because it seems to indicate that the light is injecting carriers into the system with a constant rate regardless of the relaxation time (this is of course only valid when η is still reasonably small). So, the usual interpretation or definition for injection current is

which has an intrinsic term, injection rate, that does not depend on the scattering mechanisms of materials, and an extrinsic term, the relaxation time, that varies with relaxation mechanisms.^3,43,82

Proceeding with the definition equation (29), we can get a well-known expression for circular injection current or CPGE for semiconductors,

and a slightly modified expression for metal can be easily obtained by considering the partial occupation of the Fermi–Dirac function. Like shift current, the circular injection current is also related to the phase of wave functions from the $Imvcvr(k)vvcs(k)$ and is reminiscent of the Kubo formula for Hall conductivity. Indeed, it has been shown that the trace of the circular injection current response tensor can be quantized for a two-band model of Weyl semimetal,²⁵ as shown in Fig. 12 (we always assume perfect linear dispersion, i.e., within the Lifshitz energy). The exact quantization makes use of an idealized Weyl node (which is perfectly linear with k) and the fact that only two bands are involved in the transition, which would allow us to rewrite the trace of $dγql/dτ0$ extracted from Eq. (30) [comparing with Eq. (25)] in a form that is dependent of the monopole charge C_L of the Weyl node,

Later, an experiment measuring the CPGE on a Weyl semimetal CoSi found a rough quantization.⁷¹ To be specific, at certain photon energies, only some regions adjacent to Weyl points in the Brillouin zone are responsible for the observed CPGE, and the optical transitions occur between the two bands composing the Weyl nodes. Accordingly, the circular injection current (γ, not $dγ/dτ0$⁠, as τ₀ is not accessible from this experiment) is seen to have dips and plateaus [Fig. 12(d)] that agree with the predicted quantization. However, it should be stressed that this quantization is only limited to two-band models,²⁵ and unlike quantum Hall conductivities whose quantization is robust, the exact quantization equation (31) will be removed in the presence of ubiquitous disorder and interactions.⁸³ As a result, the imperfect quantization in CoSi might not be due to experimental resolution but has a deeper physical origin. Readers who are interested in the experimental aspect of the relation between topology and BPVE are encouraged to resort to Ref. 84.

There are other variations of Eq. (30) that are used to explain or predict novel phenomena. In the work by Ji el al.,⁴⁶ the spatial inhomogeneity of the light spot is taken into account to rationalize the sign flip of photocurrent when measured at different positions of the sample. In another work by Fei et al.,⁶⁷ Eq. (30) is adapted so that the spin freedom is shown explicitly, and then they demonstrate that the circular injection current can be used as a mechanism to generate spin current in PT-symmetric antiferromagnetic insulators (breaking P- and T-symmetry individually but preserving the PT-symmetry). The basic idea is that the circularly polarized light will make the $Imvcvr(k)vvcs(k)$ have different signs for spin-up and spin-down electrons under PT-symmetry, causing their respective photocurrents to flow along opposite directions. In contrast to the PSC generated for the nonmagnetic materials using linear injection current, which requires a large SOC (Fig. 11),⁷³ the spin current generated via circular injection current mechanism is insensitive to SOC and therefore fundamentally different from the linear scenario, though it is quite similar to the linear shift current mechanism in generating PSC. A summary of using BPVE for generating charge and spin current can be found in the work.⁶⁹ 

To close the discussion on circular injection current, we want to emphasize that some important information will be lost when using the interpretation of injection rate $jq,CP−diag/dτ0$⁠, as pointed out by Gao.⁸² To see this, let us rewrite the energy term in Eq. (28) as

The second term at the RHS of Eq. (32) is weakly dependent on τ₀ and will, thus, be dropped when taking the derivative $jq,CP−diag/dτ0$⁠, but it turns out that the current associated with this term, which we denote as $jppq,CP−diag$⁠, is totally physical and can cancel the sub-bandgap contribution from $jppq,CP−sh$⁠.

To see how $jppq,CP−diag$ and $jppq,CP−sh$ combine, the first step is to realize that by symmetrizing n and l in Eqs. (27) and (28), the terms involving the derivative of velocity matrices in $jppq,CP−sh$ can be rewritten in the form of

where the k-dependence has become implicit for conciseness, and the Berry connection terms can be shown to vanish by switching n and l. In the limit $η→0$⁠, the terms involving the diagonal velocity matrices and the energy denominators in $jq,ppCP−diag$ can be recast into

Thus, combining $jq,ppCP−diag$ and $jppq,CP−sh$⁠, we can get a new term, which we call Fermi surface response:

where integration by part has been used and the boundary terms will be zero since the Brillouin zone is a closed manifold.

For a semiconductor, the Fermi–Dirac distribution function (⁠$fnk$ and $flk$⁠) will be either 0 or 1 and constant throughout the Brillouin zone, so their derivatives over k will be zero. As a result, $jq,Fermi$ will vanish in semiconductors, meaning that no sub-bandgap shift current will occur.⁹ On the other hand, it can be seen that $jq,Fermi$ will survive for a metal where the bands crossing the Fermi surface will be partially occupied, so that $fnk$ and $flk$ can indeed have non-zero derivative with k. This is another contribution to BPVE, in addition to the circular injection current, and $ωjq,Fermi$ can be shown to be quantized as well for a single Weyl node if ω is much smaller than the separation between the crossing point and the Fermi level. Yet in contrast to the circular injection current which has an energy selection rule, $jq,Fermi$ originates from the electronic excitation at Fermi surface and will always contribute to the current regardless of the frequency of light. Thus, it can be regarded as a non-resonant contribution. Moreover, the Fermi surface contribution is an intrinsic mechanism similar to shift current in the sense that it is independent of the scattering time and thus insensitive to impurities.

Our derivation thus far is based on the belief that photocurrent will appear only when the photon energy $ℏω$ is equal or larger than the bandgap E_g for gapped systems. This is physically reasonable because the optical absorption will only occur when $ℏω>Eg$⁠, which will allow the photo energy to be converted to electrical energy. Except the obvious scenarios where in-gap defect states could induce sub-bandgap photocurrent,^26,27 there have been reports demonstrating that sub-bandgap photocurrent can indeed exist more intrinsically, and we are going to show why they do not conflict with our derivation.

The work by Kaplan et al. shows that for linearly polarized light, the small imaginary numbers in the denominators of Eq. (17) do not necessarily all equal to η, i.e., they could be $1/(εnk−εmk−iη1)$ and $1/(εnk−εlk−ℏΩ−iη2)$ with $η2≠2η1$⁠. The justification for this modification is that there is no guarantee that different bands have the same relaxation time⁸⁵ (or the adiabatic switching-on rate is different for different transition processes), and therefore, for systems without time-reversal symmetry, the exact cancelation similar to Eqs. (32)–(34) no longer works since the second term of the RHS of Eq. (32) will have a different weight, and one may get a non-vanishing sub-bandgap photocurrent.

The second type of sub-bandgap current deals with metallic systems with a Fermi surface and an optical bandgap.^86,87 This means that the even though some bands are partially occupied, there is no optical transition below the optical bandgap because the optical transition is vertical in k-space. Then, the arguments and derivations in Sec. III B 3 can be applied here and demonstrate that a sub-bandgap signal is predicted due to the Fermi-surface contribution. It should be stressed that optical bandgap for metallic systems is different from our previous considerations on perfect semiconductors and insulators with the Fermi level lying within the bandgap. The latter case has no Fermi surface, and as a result there is no sub-bandgap photocurrent.

The third scenario where sub-bandgap photocurrent can exist is in strongly correlated systems. In the first place, wherever there exists strong correlation, it is unclear whether the BPVE theory developed for independent-particle systems can be applied. Nonetheless, we can still employ some of the intuition we have developed so far. For example, it is reported that excitonic states can contribute to shift current, which can be understood as the change of coordinate from a delocalized state to a localized excitonic state.⁸⁸ Since the eigen-energies of bound exciton states are usually below the bandgap (defined as the difference between the conduction band minimum and the valence band maximum), this will contribute to the sub-bandgap photocurrent. Another work has examined the photocurrent of ferroelectric excitonic insulators where the strong correlation will even break the ground-state symmetry. As a result of the phase transition, photocurrent will appear below the bandgap corresponding to the non-interacting system.⁸⁹ However, one should notice that for such scenario, if one defines the bandgap as the minimum energy required for optical transition to happen for a given system (ground-state to excitonic state or ground-state to excited state for the symmetry-breaking system), then these photocurrents are still essentially above bandgap.

Until now, we have covered all the major contributions to BPVE for linearly and circularly polarized light, from the perspective of time-dependent perturbation theory. All the mechanisms are second-order in the electric field of light, so they can explain the linear scaling of the experimental photocurrent with light intensity. The more exotic photon-drag effect by considering the non-vertical optical transitions (non-zero momentum carried by light) has also been discussed by several authors.^1,90,91 Moreover, there could be higher order contributions, such as the jerk current originating from the third-order response to electric field (though it is essentially discussing second-order response to the oscillating electric field from light and first-order response to a co-existing static electric field).^92,93 Note that the co-existence of static and oscillating electric field is also implicitly considered in the works^22,94 where the atomic displacements can be driven by a static electric field, which would further influence the BPVE. Readers who are interested in these mechanisms are encouraged to refer to the original literature cited therein. We also note here that there have been many discussions about the relation between BPVE and quantum geometry or topology, but this is largely beyond the scope of this review. So, we will just provide several representative works here for interested readers.^{25,83,95–99}

The theories presented in Secs. III A and III B are all based on time-dependent perturbation theory by treating the optical field as a small perturbation. However, recently a different BPVE formalism has been formulated from the Floquet theory, where the optical field is not necessarily weak. The benefit of the Floquet theory is not immediately obvious, but when combined with non-equilibrium transport theory, it can be easily adapted to investigate the BPVE in a finite system with explicit attachment of electrodes as well as randomly distributed disorder.^8,23,44 Here, we outline the framework of the Floquet theory of BPVE in this section and demonstrate that it can lead to the same shift current expression for linearly polarized light. In Sec. III D, we will review its application in finite systems such as Anderson insulators.

We start by giving a brief review of Floquet theory. For a general quantum system whose state is described by a state vector $|Ψ(t)⟩$⁠, its dynamics is determined by the time-dependent Schrödinger equation,^100,101

The Floquet theorem states that for a periodic Hamiltonian

where T is the periodicity, there exists a solution in the following form:

Here, $εα$ is some quasi-energy that must be solved for, $|uα(t)⟩$ is a function periodic in time: $|uα(t+T)⟩=|uα(t)⟩$⁠, and ${α}$ is a set of parameters that label different solutions. This is reminiscent of the Bloch theorem which states that for a potential periodic in real space, $V̂(r+R)=V̂(r)$⁠, the solutions of the time-independent Schrödinger equation must have the form $|ψ(r)⟩=eik·r|u(r)⟩$⁠, where $|u(r)⟩$ is a lattice-periodic function.

Since $|uα(t)⟩$ has the periodicity T, it can be Fourier transformed into: $|uα(t)⟩=∑ne−inωt|uαn⟩$⁠, with $ω=2π/T$⁠. Now, with the definition of Floquet mode $|uαn⟩$⁠, Eq. (36) can be Fourier transformed,

One can see that the dimension of the original Hamiltonian $Ĥ(t)$ has been augmented by the Floquet indices n and m, and in fact, they span all integers from $−∞$ to $+∞$⁠. However, since each element in the matrix $Ĥmn$ corresponds to the transition probability from nth Floquet mode to mth Floquet mode, in practice if one only considers the low-energy excitations from the ground state, then the Floquet indices will be truncated to a small value to reflect the consideration that the higher-order excitations are neglected. As another note, the Floquet modes are usually unknown, but the matrix $Ĥmn$ has explicit forms provided that $Ĥ(t)$ can be expressed in a known basis, so one must diagonalize $Ĥmn$ to find the Floquet modes $|uαn⟩$ and the quasi-energies $εα$⁠.