Low-frequency divergence and quantum geometry of the bulk photovoltaic effect in topological semimetals

We study the low-frequency properties of the bulk photovoltaic effect in topological semimetals. The bulk photovoltaic effect is a nonlinear optical effect that generates DC photocurrents under uniform irradiation, allowed by noncentrosymmetry. It is a promising mechanism for a terahertz photodetection based on topological semimetals. Here, we systematically investigate the low-frequency behavior of the second-order optical conductivity in point-node semimetals. Through symmetry and power-counting analysis, we show that Dirac and Weyl points with tilted cones show the leading low-frequency divergence. In particular, we find new divergent behaviors of the conductivity of Dirac and Weyl points under circularly polarized light, where the conductivity scales as $\omega^{-2}$ and $\omega^{-1}$ near the gap-closing point in two and three dimensions, respectively. We provide a further perspective on the low-frequency bulk photovoltaic effect by revealing the complete quantum geometric meaning of the second-order optical conductivity tensor. The bulk photovoltaic effect has two origins, which are the transition of electron position and the transition of electron velocity during the optical excitation, and the resulting photocurrents are respectively called the shift current and the injection current. Based on an analysis of two-band models, we show that the injection current is controlled by the quantum metric and Berry curvature, whereas the shift current is governed by the Christoffel symbols near the gap-closing points in semimetals. Finally, for further demonstrations of our theory beyond simple two-band models, we perform first-principles calculations on magnetic Dirac semimetal MnGeO$_3$and Weyl semimetal PrGeAl. Our work brings out new insights into the structure of nonlinear optical responses as well as for the design of semimetal-based terahertz photodetectors.

We study the low-frequency properties of the bulk photovoltaic effect in topological semimetals. The bulk photovoltaic effect is a nonlinear optical effect that generates DC photocurrents under uniform irradiation, allowed by noncentrosymmetry. It is a promising mechanism for a terahertz photodetection based on topological semimetals. Here, we systematically investigate the low-frequency behavior of the second-order optical conductivity in point-node semimetals. Through symmetry and power-counting analysis, we show that Dirac and Weyl points with tilted cones show the leading low-frequency divergence. In particular, we find new divergent behaviors of the conductivity of Dirac and Weyl points under circularly polarized light, where the conductivity scales as ω −2 and ω −1 near the gap-closing point in two and three dimensions, respectively. We provide a further perspective on the low-frequency bulk photovoltaic effect by revealing the complete quantum geometric meaning of the second-order optical conductivity tensor. The bulk photovoltaic effect has two origins, which are the transition of electron position and the transition of electron velocity during the optical excitation, and the resulting photocurrents are respectively called the shift current and the injection current. Based on an analysis of two-band models, we show that the injection current is controlled by the quantum metric and Berry curvature, whereas the shift current is governed by the Christoffel symbols near the gap-closing points in semimetals. Finally, for further demonstrations of our theory beyond simple two-band models, we perform first-principles calculations on the shift and injection photocurrent conductivities as well as geometric quantities of antiferromagnetic MnGeO3 and ferromagnetic PrGeAl, respectively, as representatives of real magnetic Dirac and Weyl semimetals. Our calculations reveal gigantic peaks in many nonvanishing elements of photoconductivity tensors below photon energy ∼0.2 eV in both MnGeO3 and PrGeAl. In particular, we show the ω −1 enhancement of the shift conductivity tensors due to the divergent behavior of the geometric quantities near the Dirac and Weyl points as well as slightly gapped topological nodes. Moreover, the low-frequency bulk photovoltaic effect is tunable by carrier doping and magnetization orientation rotation. Our work brings out new insights into the structure of nonlinear optical responses as well as for the design of semimetal-based terahertz photodetectors.

I. INTRODUCTION
Topological semimetals are emerging as efficient infrared and terahertz photodetectors [1]. In contrast to semiconductors whose absorption spectum is bounded below by the band gap, semimetals can detect radiations down to the terahertz range because of their gapless spectrum. A promising mechanism for the generation of photocurrents in semimetals is the bulk photovoltaic effect. It refers to the generation of photocurrent under uniform irradiation of light due to the intrinsic inversion asymmetry of the system. Since the bulk photovoltaic effect does not require a bias voltage for breaking inversion symmetry, dark current noise can be suppressed [1].
To achieve high photosensitivity, we need to understand how to obtain large photoconductivity. It is believed that band topology plays an important role [2][3][4][5][6]. The bulk photovoltaic effect occurs due to the inversion-asymmetric transition of electron position or velocity during the optical excitation, and the resulting photocurrents are respectively called the shift current and the injection current [7]. In nonmagnetic systems, linearly polarized light induces shift currents while circularly polarized light induces injection currents. Remarkably, both the linear shift [4] and circular injection [5] currents were found to be intimately related to the topological quantities, the Berry connection and the Berry curvature, respectively. These discoveries have led to various theoretical and experimetnal studies searching for topological enhancement near the gap-closing points [8][9][10][11][12][13][14][15][16][17][18][19].
However, while there is a concrete proportional quantitative relationship between the injection current and the Berry curvature, no such a quantitative relation exists between the shift current and the Berry connection. For example, a Dirac point in two dimensions has a quantized π Berry phase, thus having a nontrivial Berry connection. Nevertheless, such a Dirac point does not generate a shift current because of its inversion symmetry. Furthermore, no simple quantitative relation was found between the shift current and the shift vector [20], a gauge-invariant quantity related to the Berry connection, without some special requirements like the momentum independence of the dipole matrix elements [21].
The bulk photovoltaic effect in magnetic topological semimetals is more poorly understood, although a recent work has revealed some general aspects [22]. Due to time-reversal symmetry breaking in magnetic systems, linearly (circularly) polarized light can generate injection (shift) currents as well as shift (injection) currents. There are some works highlighting the generation of large linear injection currents in magnetic systems [22][23][24][25][26][27], but the relationship between the response and band topology has not been understood. Moreover, there has been very little attention to the circular shift current, while the first concept of the shift current appeared as a response to circularly polarized light [28][29][30]. One reason for this is that injection currents are typically stronger than the shift currents. Since there is a rapid progress in the experimental observation [31][32][33] and theoretical proposal [34][35][36][37][38][39][40][41][42] of various magnetic topological semimetals, addressing optical properties in magnetic topological semimetals is now a timely subject.
In this work, we reveal general low-frequency properties of the shift and injection currents in magnetic and nonmagnetic point-node semimetals. By employing symmetry and powercounting analysis, we determine the leading low-frequency divergence near the gap-closing point, as summarized in Table I. We show that tilted Dirac and Weyl points can generate the leading divergence. While the bulk photovoltaic response in tilted Dirac and Weyl points have been studied [13,[23][24][25], here we approach them in a unified view and discover new aspects. Our theory covers type-I and type-II spectrum of Dirac and Weyl points in any dimensions. In particular, our analysis include type-II Dirac points in two dimensions and type-I and type-II Dirac points in three dimensions, which are relevant to magnetic Dirac semimetals, whose shift and injection currents have not been previously studied. It is widely known that the protection of Dirac points against opening the gap requires inversion symmetry, which forbids the bulk photovoltaic response. That is true in nonmagnetic systems. However, in magnetic systems, symmetry under the combination of spatial inversion P and time reversal T , which is P T symmetry, is enough for the protection [39][40][41][42], so inversion symmetry can be broken. Also, we study the largely overlooked circular shift current in Weyl and Dirac systems. The circular shift current grows fast as the photon frequency gets smaller, and it scales as ω −1 in three dimensions. This indicates that Dirac semimetals in three dimensions can show divergent photovoltaic responses like Weyl semimetals, although the Berry curvature is identically zero due to P T symmetry. In two dimensions, it grows faster as ω −2 .
We find that tilted Dirac and Weyl points show what we call the separation of responses, where photocurrents of different origins manifest through different current directions. For example, circular shift and circular injection currents flow in different directions. This can be useful in the detection of shift currents in the coexistence of the stronger injection currents. The separation of responses can occur from the symmetry transformation property under magnetic operations M T or C 2 T , where M and C 2 are mirror and twofold rotation, so it remains robust in the system beyond the k-linear approximation as long as those symmetries are preserved.
Our symmetry and power-counting analysis are enough for understanding the general pattern of the response for systems with a linear spectrum. However, for a deeper understanding of the response, we propose a new perspective on the lowfrequency bulk photovoltaic effect. We uncover the full geometric nature of shift and injection currents. Here, as well as the Berry curvature, another geometric quantity called the quantum metric has a crucial role. Since geometry has more information than topology about quantum states, we observe the consequence of band topology in a more broad perspective through geometric quantities. We show that the linear injection conductivity is determined only by the quantum metric near the gap closing, in the same way as the circular injection conductivity is determined by the Berry curvature. This completes the geometric understanding of the injection currents at the low-frequency regime. Furthermore, we show that the shift conductivities are related to a more interesting quantity, which is the Christoffel symbols. Unlike the Berry connection, which has a complicated relationship to the magnitude of the shift current, the Christoffel symbols directly control the magnitude of the response. In this viewpoint, the enhancement of the shift and injection current responses near the gap-closing point can be attributed to the divergence of the geometric quantities at the geometric singularity, i.e., the gapclosing point. Furthermore, our unique perspective allows to view the bulk photovoltaic effect as novel tools for experimental measurements of quantum geometry in materials.
Our geometric interpretation of the bulk photovoltaic effect is clearly distinguished from the one in Ref. [22], which also discuss the role of the quantum metric. In Ref. [22], the shift current is decomposed into four parts, and one of them was interpreted as the geometric contribution. The injection current is identified to have no geometric origin. These are in contrast to our interpretation in which both shift and injection currents are fully geometric responses. Only through our direct quantitative relationships the bulk photovoltaic effect can be identified as a useful measure of quantum geometry in experiments as well as a way to theoretically understand the low-frequency divergent behavior.
Finally, for further demonstrations of our theory beyond simple two-band models, we perform first-principles relativistic band theoretical calculations of the shift and injection photocurrent conductivities as well as geometric quantities in antiferromagnetic (AF) MnGeO 3 and ferromagnetic (FM) PrGeAl (as will be reported in Sec. VI below), respectively, as representatives of real magnetic Dirac and Weyl semimetals. In AF MnGeO 3 , although both T and P symmetries are broken, the combined P T symmetry is preserved [42]. Thus, AF MnGeO 3 was recently predicted to a rare magnetic Dirac semimetal. [42] As our theory predicts (Table I), we find nonzero elements of circular shift and linear injection photoconductivity tensors in AF MnGeO 3 . In fact, several nonzero elements exhibit huge peaks below photon energy of ∼0.2 eV. Our calculations reveal that at least there are three Dirac points in the vicinity of the Fermi level and two of the Dirac points are accompanied by a slightly gapped Dirac point each. The calculated quantum metric and Christoffel symbol of first kind exhibit divergent behaviors near both the gapless and gapped Dirac points, thus leading to the geometric enhancement in linear injection and circular shift currents, respectively. In contrast, both T and P T symmetries are broken in FM PrGeAl [43], and hence all four types of the bulk photovoltaic effect may emerge in FM PrGeAl, as our theory predicts (Table I). Indeed, we find that many nonzero elements   (9) and (10).
In T -(P T -) symmetric systems, only the responses with the positive T -(P T -) parity appears. All of the four responses can appear when T and P T symmetries are both broken. Based on their symmetry properties, we call the linear injection and circular shift current responses as P T -symmetric responses and the linear shift and circular injection current responses as T -symmetric responses. Symplectic Christoffel symbols indicate the symplectic analog of the Christoffel symbol of the first kind (See Eq. (33)). Here we consider only interband-transitive processes in the clean limit where the relaxation rate Γ = τ −1 is smaller than the photon frequency ω. In this case, the injection current is typically larger than the shift current, and they become comparable when ω approaches Γ.
of all shift and injection conductivity tensors show gigantic peaks in the low-frequency range up to 0.1 eV. There are at least 160 type I and type II Weyl points within ±0.1 eV of the Fermi level in FM PrGeAl [43]. Our calculations indicate that the low-frequency peaks can be further increased by up to a factor of 5 by raising the chemical potential from the Fermi level to the energy of certain Weyl points. Our calculations also reveal the divergent behaviors of both quantum geometric tensor and quantum geometric connection (symplectic Christoffel symbol and Christoffel symbol of first kind) near the Weyl points and also anticrossing topological nodes, thus leading to the gigantic low-frequency shift and injection currents in FM PrGeAl. Furthermore, we notice that FM PrGeAl is a soft ferromagnet [44]. Thus, the magnetization direction can be easily rotated away from the easy c-axis to, e.g., aaxis, and this may cause a topology change of the Weyl point distribution in the Brilloin zone [43], thereby resulting in significant changes in the shift and injection photoconductivities. Also, since the circular shift and injection conductivity tensors are antisymmetric with respect to the magnetization direction, they would change sign when the magnetization direction is reversed. All these observations show that the gigantic lowfrequency shift and injection photocurrents in FM PrGeAl can be tuned by either carrier doping or magnetization direction rotation.
The outline of this paper is as follows. We explain the shift and injection currents as a second-order optical response in Sec. II. Then, we study the symmetry and low-frequency divergence of the shift and injection currents in Sec. III. Section IV enriches our analysis by revealing the quantum geometrical nature of the low-frequency responses. We elaborate more on the symmetry and divergence with concrete models and numerical calculations in Sec. V. In Sec. VI, we present our first-principles calculations on the shift and injection photocurrent conductivities as well as geometric quantities of AF MnGeO 3 and FM PrGeAl. In particular, we analyze the gigantic peaks in the calculated low-frequency photoconductivity spectra in MnGeO 3 and PrGeAl in terms of the divergent behaviors of the geometric quantities near the gapless and slightly gapped topological nodal points. Finally, we discuss several issues about the low-frequency divergence in Sec. VII.

II. SHIFT AND INJECTION CURRENTS
Let us expand current density j in increasing power of the electric field E as The first term is the familiar linear response, and the other terms are nonlinear responses. Since the current density oscillates in-phase with the electric field in the linear response, the DC photocurrent generation is inherently a nonlinear optical effect. In our work, we assume that the electric field is sufficiently small such that perturbation theory works (we discuss in Sec. VII how small it should be). While even-order responses to E vanishes in centrosymmetric systems, they are allowed in noncentrosymmetric systems. The bulk photovoltaic effect studied in the present paper is thus primarily a second-order response. The second-order optical response under the uniform illumination of light has the form j c (ω 1 + ω 2 ) = σ c;ab (ω 1 + ω 2 ; ω 1 , ω 2 )E a (ω 1 )E b (ω 2 ) in general. Let us focus on the DC (direct current) generation where j c DC = j c (0), and σ c;ab DC (ω) = σ c;ab (0; ω, −ω). In the clean limit -where the interband relaxation rate is smaller than the photon frequency and the band gap, interband transitions are described by two processes: shift and injection 1 .
The shift and injection currents correspond to the current generated by the change of the electron position and velocity, respectively, during the interband transition of electrons [7]. One can see this by noting that the shift and injection conductivities have the form of the Fermi Golden rule [7] (See Appendix C for explicit calculations). Explicitly, they have the form [46,47] σ c;ab is the energy difference between bands m and n, H|n = ω n |n , r a mn = m|i∂ a |n and v c mn = −1 m|∂ c H|n , and we use the notation ∂ a = ∂/∂k a . R c;a mn = r c mm − r c nn + i∂ c log r a mn is called the shift vector -characterizing the interbandtransition of the displacement, and ∆ c mn = v c mm − v c nn is the interband-transition of the velocity. τ is the relaxation time that saturates the injection current: withtout it, the injection of moving electrons and holes leads to a constant growth in time. We take the electron charge as −e (i.e., e > 0). Eq. (4) is valid with and without time-reversal symmetry.

III. SYMMETRY AND POWER-COUNTING ANALYSIS
To get a general perspective on the low-energy properties of the shift and injection currents, we review their symmetry properties and then study the pattern of low-frequency divergences. The key properties presented in this section are summarized in Table I, and they serve as basic ingredients for the analysis in Sec. V.

A. Symmetry of the shift and injection conductivities
Let us first decompose the second-order DC conductivity into its real and imaginary parts: Using E * (ω) = E(−ω) and Eq. (2), one can see that the conductivity can be symmetrized such that Therefore, we always consider conductivity tensors satisfying Eq. (6). We note that the expressions in Eq. (4) are already symmetrized. Assuming the form E(t) = |E|e −iωt (cos φ, sin φ, 0) + c.c. for linearly polarized light and E(t) = |E|e −iωt (1, i, 0) + c.c. for circularly polarized light, the real part of the second-order optical conductivity is The real part of the conductivity is responsible for the current generation regardless of the polarization,while the imaginary part of the conductivity is responsible for the current generated by the circularly polarized light. If one measures the current difference between the ones generated by leftcircularly polarized light and right-circularly polarized light ∝ σ c;xy − σ c;yx , only the imaginary part contributes. In the following, we call the real part σ c;ab L as linear conductivity and the imaginary part σ c;ab C as circular conductivity. From the definition in Eq. (2) and the transformation properties of the current and electric field, it is clear that the second-order optical conductivity transforms like a third-rank tensor under spatial transformations. That is, σ c1;a1b1 2 . However, one should be careful when taking time reversal for relaxational processes as is well-known from the Onsager reciprocity relations in linear response theory [48][49][50]. For example, it seems like Eq. (2) implies that time reversal reverses the sign of the DC conductivity for linearly polarized light. However, we need to additionally reverse the sign of phenomenological relaxation rate Γ, in order to make the decay in time to the growth in time. The correct time reversal for the second-order DC conductivity is See Appendix D. When applying this equality, the delta function and τ should be interpreted as and Γ −1 , so they reverses sign under Γ → −Γ. Thus, the real part of the shift and injection conductivity tensor transforms as and the imaginary part transforms as under the spacetime symmetry transformation (t, Alternatively, one can verify these transformation rule by examining transformations of R c,a mn , r a nm , ∆ c mn from the form in Eq. (4). See Appendix E for a derivation.
Knowing these symmetry transformation properties, one can use the MTENSOR [51] in the Bilbao Crystallographic Server to see which tensor components are required to vanish by symmetry for any of 80 magnetic point groups. For our purpose, the most important symmetries are simply time reversal T and spacetime inversion P T symmetries (note again 2 We define σ c;ab DC by j c DC = σ c;ab DC E a E b , where j c 1 = Mc 1 cj c and E a 1 = Ma 1 aEa 2 . that P is always broken in the present paper). We summarize their role in Table. I. While only linear shift and circular injection currents can be generated in time-reversal-symmetric systems, they vanish in P T -symmetric systems, so linear injection and circular shift currents can be generated. timereversal symmetry and spacetime inversion symmetry are thus complementary to each other, as pointed out in [22]. This behavior is manifested in the geometric quantities related to the responses. As we show below in Sec. IV, T -symmetric responses are related to both of the Berry curvature and the quantum metric, while P T -symmetric responses are related to the quantum metric only. In general magnetic systems without T and P T symmetries, all of the four types of currents can be generated.
When there is M T or C 2 T symmetry instead, where M and C 2 are mirror and twofold rotation, a phenomenon of "separation of responses" occurs, meaning that different directions manifest different types of responses. It is because those symmetries act like time reversal in some directions and act like spacetime inversion in the other directions. For example, M x T acts like time reversal in the y and z directions, whereas it acts as spacetime inversion in the x direction. In this case, x-polarized light generates a shift current along y and z while generating an injection current along x. We demonstrate the separation of responses through model calculations in Sec. V and also through first-principles calculations in Sec. VI.

B. Power-counting analysis of the low-energy divergence
Let us now examine the low-energy divergence of the second-order responses in semimetals. We can estimate the divergence by counting the power of photon frequency in Eq. (4). Since the delta function has dimension ω −1 , R c and r a mn has dimension k −1 , and ∆ c has dimension ω/k, the shift and injection conductivity scales as where k is the characteristic wave vector. When the system has dispersion E ∝ k α , k ∼ E 1/α ∼ ω 1/α . Thus, smaller α is preferred to get large optical responses for small ω in 2D, while it is independent on the dispersion in 3D. In lattice systems, stable Weyl points always (and also Dirac points protected by symmorphic symmetries [52]) appear pairwise according to the Nielsen-Ninomiya theorem [53,54]. Therefore, a response in semimetals should be a sum of responses from different Weyl or Dirac points. However, different gap-closing points are located at different energy levels in general, so it is expected that only a particular point contributes to the low-energy response significantly [5]. An exact cancellation or reinforcement among different gapclosing points can occur due to symmetries, but it can be considered straightforwardly from the symmetry transformation properties of the conductivity tensor. In this regard, the re-sponse from a single gap-closing point can be associated with the low-energy response of a whole system.
Thus, we consider a gap-closing point with a linear dispersion, described by a Dirac Hamiltonian where Γ a are mutually anticommuting matrices, such that the spectrum has the form In this case, we have through a dimensional analysis. This divergence is expected to occur in the absence of a symmetry cancellation. However, Eq. (12) has too many symmetries, so we need to break them in general. A Dirac point has inversion symmetry (by definition it is nonchiral), so the second-order optical response is forbidden. Also, as noted in Ref. [13,22,25], a Weyl point described by Eq. (12) does not show a second-order optical response by linearly polarized light. It is because Eq. (12) has SO(d) rotational symmetry in d spatial dimensions. In 3D, the absence of mirror symmetry (chirality) in Weyl semimetals allows one unique nonvanishing independent component σ 3;12 inj,C under circularly polarized light. Because a Weyl point in 3D has time-reversal symmetry around the gap-closing point, the circular shift current is also forbidden, and the generated DC current is the circular injection current. Even for a more general linear dispersion described by H = d a,i=1 v ai k a Γ i , which has apparently less symmetry, only circular injection currents for a Weyl point can be nonvanishing. It is because the conductivity for this Hamiltonian is given by σ c;ab DC,0 is the conductivity for v ai = vδ ai , which is the case for Eq. (12) (See Appendix G for a derivation).
The only way to generate the leading divergence for shift currents and linear injection currents is to tilt the Dirac or Weyl cone, as shown in Fig. 1. Since it allows anisotropic optical excitations around the gap-closing point, photocurrents can flow whose direction depends on the direction of the tilting. To see that other symmetry breaking gives subleading power in ω −1 , let us add symmetry-breaking perturbations to the Dirac Hamiltonian.
Here, λ p,a is a constant parameter, Γ a =0 are mutually anticommuting matrices, and we also include Γ 0 as the identity matrix. Since the dimensionless parameter is λ p,a k p ω −1 ∼ λ p,a ω p−1 , responses due to perturbations in λ p,a have weaker Here, v is the velocity of a Dirac or Weyl fermion at zero tilting, and v is the overall velocity shift by µ → µ + v kx, giving rise to a tilting of the cone. The gray planes in (b,c) show the Fermi level. (d,e) Optically active region in momentum space for (d) type-I and (e) type-II cases. Small black dots at the center indicates the location of the gap-closing point. Red circles around the point shows the surface satisfying ω = ωc − ωv = 2 vk, where k = 0 at the gap closing. Both conduction and valence bands are unoccupied in the white region, only the valence band is occupied in the light gray region, and both bands are occupied in the gray region. Electrons can be optically excited only on the solid red arcs (i.e., θ− < θ < θ+, where θ is the absolute value of the polar angle in either 2D or 3D), which is in the light gray region. low-frequency divergences for p > 1. For example, let us consider a quadratic correction, resulting in the dispersion relation E = ±( vk + λ 2 k 2 ). Optical excitations occur in the region satisfying . By inserting this to Eq. (11), we have corrections to the leading divergence by a fraction (d − 3)λ 2 ω/v 2 1 for small ω. At the p = 0 order, λ 0,a comes as either the chemical potential, the shift of the location of the gap-closing point, or a mass term 3 , i.e., we can write the Hamiltonian as None of them generate nonvanishing shift and linear injection currents. µ breaks no symmetry, k−k 0 can be redefined as k such that gap closes at k = 0, and the mass term can serve as an inversion symmetry operator so that it forbids second-order optical responses. Let us now consider p = 1. As we show above, the linear spectrum without tilting has zero shift and linear injection current responses. Thus, the remaining possibility is tilting the cone by adding λ 1,0 = −v cos θ: this gives an overall tilting of energy levels by µ → µ + v k x .
Our analysis shows that the shift and injection conductivity tensors of tilted massless Dirac and Weyl points have the form Let us explain how Fs depend on v /v and 2µ/ ω in general. When µ = 0, i.e., when the Fermi level is away from the gapclosing point, the chemical potential sets the lowerbound for frequency, so Fs do not diverge as ω → 0.
When the Fermi level is exactly at the gap-closing point, i.e., µ = 0, Fs show significantly different behaviors for |v /v| < 1 and |v /v| > 1, which are called type-I and type-II [55], respectively [See Fig. 1(b,c)]. In the type-I case, tilting cannot generate shift and injection currents when µ = 0. It is because, in this case, the Fermi surface is a point, so anisotropic excitation cannot occur. Thus, only circular injection currents can be generated, which do not need tilting for its generation.
In contrast, in the type-II case, the Fermi surface has a finite size at µ = 0, so anisotropic excitation can occur in principle. However, it depends on whether the response is T -symmetric or P T -symmetric. While the T -symmetric responses (linear shift and circular injection) has a nontrivial response at µ = 0 [13], the P T -symmetric responses (linear injection and circular shift) has a vanishing response at µ = 0. This is related to the fact that a Hamiltonian with only k-linear terms has an emergent CP T symmetry, where C is the particle-hole operator, and CP T = 1: At µ = 0, P T -symmetric responses should be zero since they respect an effective C symmetry also (where (P T ) −1 takes the role of an effective C operator) -which reverses the direction of the current, while T -symmetric responses lack C symmetry such that they can be nontrivial.
This shows that magnetic and nonmagnetic systems have the same low-frequency divergent behavior when µ = 0 exactly. Nevertheless, as far as µ is small but not exactly zero, we can still expect enhanced P T -symmetric responses at small frequency ω ∼ 2µ/ by a factor in front of Fs in Eq. (16). Therefore, P T -symmetric responses in magnetic systems also can show a divergent behavior associated with the ω → 0 limit with a fixed ratio of 2µ/ ω.
We investigate the symmetry properties of the tilted Dirac and Weyl points more closely in Sec. V by explicitly calculating the conductivity tensors. Before that, we derive the general formula for the shift and injection conductivities for arbitrary Dirac Hamiltonians and provide their geometric aspects in the following sections. It adds more perspective on the transformation rule Eqs. (9) and (10) and the divergent behavior near the gap-closing point.

IV. QUANTUM GEOMETRIC ASPECTS
In the previous section, we analyze the overall trend of shift and injection currents using symmetry and power-counting analysis. Here we show that every detail of the conductivity profile for the shift and injection currents is determined by quantum geometric quantities in the low-frequency regime. It was pointed out in Refs. [4,6] that the shift current is related to quantum geometry because the shift vector includes the Berry connection -a geometric quantity. Following these works, the geometric nature of the shift and injection currents was previously discussed in several works [5,22]. However, no simple quantitative relationship between the response and the geometric quantities has been found except for the circular injection conductivity [5]. In this section, we show that the shift and injection conductivities are proportional to geometric quantities that have natural geometric meaning on the Bloch sphere. These relationships are not limited to massless Dirac and Weyl points and are exact for any two-band system or P T -symmetric four-band system. It implies that interbandtransitive photovoltaic responses at the low-frequency regime probe the quantum geometry of materials. In this perspective, the low-frequency divergence of the shift and injection current responses of gap-closing points can be attributed to their geometrically singular nature. Also, time reversal symmetry transformation of the conductivity tensors, which are quite confusing, can be simply understood from the transformation properties of the geometric quantities. In the following, we first derive the formula that relates injection and shift currents with the Bloch vector of general Dirac Hamiltonians with arbitrary matrix size in any spatial dimensions. Then, we provide a quantum geometric interpretation of our formula.

A. Shift and injection conductivity for Dirac Hamiltonians
We consider the low-energy effective model systems described by the following d M × d M Dirac Hamiltonian where Γ i are mutually anticommuting matrices. This Hamiltonian described a single Dirac or Weyl point when f i (k) = k i , but here we do not need to assume linear dispersion and consider general form of f i (k)s. In particular, the above Hamiltonian describes general two-band Hamiltonians when d M = 2, where Γ i=1,2,3 are three Pauli matrices, and it describes general four-band P T -symmetric Hamiltonians (with ..,5 are five Gamma matrices. Let us express the injection and shift conductivities in terms of f i s. This makes theoretical analysis and numerical calculations convenient. As for the injection current, one can integrate the delta function easily by using ∆ c mn = ∂ c ω mn such that ∆ c mn δ(ω mn −ω) = ∂ c Θ(ω mn −ω). After that, we obtain for ω > 0, where ω cv is the energy gap between the conduc-tion and valence bands,n is the surface normal vector, and is the so-called quantum geometric tensor [56] (See Appendix H for a derivation of the second equality). Here, The vanishing of J ij for d M = 4 is due to the presence of P T symmetry. Q ba is called the quantum geometric tensor because its real and imaginary parts are related to the quantum metric g ba and the Berry curvature F ba by The relationship between the circular injection current and the Berry curvature was found in Ref. [5]. On the other hand, the role of the quantum metric in determining the linear injection current was not discussed in the literature. We explain more on the geometric meaning of the quantum geometric tensor in Sec. IV B.
The shift current has more complicated form, and it can be related to the matrix elements of the derivatives of the Hamiltonian as [57] R c,a mn r b nm r a where w ac mn = −1 m|∂ a ∂ c H|n is the diamagnetic term. The second line involves virtual transitions among three different bands, so it vanishes in our Dirac system that has only two energy levels with energy ω c = f 0 +f and 4 In fact, the integrand should be (C bca − C * acb )/2 since the integrand has the form (R c,a mn − R c,b nm )r b nm r a mn , but C bca gives the same value of the real part of the current. For linearly polarized light, the conductivity is symmetric under a ↔ b. In this case, C bca + C * acb + a ↔ b is real, such that it contributes to the imaginary part of the conductivity, and thus to the imaginary part of the current. For circularly polarized light, since C bca + C * acb − a ↔ b is imaginary, the conductivity takes a real value, so it contributes to the imaginary part of the current for circularly polarized light.
for ω > 0, where A special case of this formula was derived in Ref. [57] for two-band models with time-reversal symmetry. Our formula in Eq. (23) extends the existing formula to describe arbitrary systems described by the Dirac Hamiltonian in Eq. (17). We show below that C bca has geometric meaning as a connection.
B. Geometry on the generalized Bloch sphere Let us explain the geometric meaning of the quantum geometric tensor Q ba as geometric quantities defined on the generalized Bloch sphere. This point of view is helpful for understanding the geometric meaning of C bca as well as that of Q ba .
We consider the following general Dirac Hamiltonian in Eq. (17). Then, the generalized Bloch vector f (k) is a map where d Γ is the number of Gamma matrices. This map defines a pull-back of the quantum geometric tensor from the f -space to the Brillouin zone. Let us recall that the quantum geometric tensor has the following form where and f = |f |. This is a pull-back of the quantum geometric tensor q ij defined in the f -space to the momentum space by a transformation ∂ a f i of tangent vectors where ∂ i = ∂/∂f i . The quantum metric and the Berry curvature are given by where η ij and ij are the real and imaginary parts of q ij , i.e., q ij = η ij − i ij . In this viewpoint, the quantum metric and the Berry curvature are pull-backs of the metric η ij and the symplectic form ij defined on the Bloch sphere. The metric η ij measures length ds through ds 2 = η ij df i df j , and the symplectic form ij measures the oriented area dA through dA = ij dx i dx j in the f -space. While the geometric quantity η ij is defined on the whole f space, it can be regarded to measure the length on the unit sphere with f = 1, and it is irrelevant for the radial direction f = |f | 5 . For example, ds 2 = (1/4)(dθ 2 + sin 2 θdφ 2 ) in polar coordinates when d M = 2. To see this without a coordinate transformation, first note that P ij = δ ij −f ifj is the projection to the plane perpendicular tof . The metric measures the length only along the angular directions on a sphere with a fixed f . Also, the f −2 factor normalizes the length such that only the angle between two points on a sphere is measured. Similarly, the symplectic form also measures the area on the unit sphere. In this sense, the quantum geometric tensor is a geometric quantity defined on the generalized Bloch sphere (f = 1).
Another geometric quantity called the Levi-Civita connection also can be constructed on the generalized Bloch sphere. Its components are called the Christoffel symbols, and they can be written in two ways -the first and the second kind. The Christoffel symbols of the second kind γ k ij are defined by where e i is the projection of the unit vector along the f i direction to the tangent space by P ij = δ ij −f ifj . It measures how vectors and tensors are changed as we move them parallel to the direction of the curved surface (which is the generalized Bloch sphere in our case) 6 . We have It is related to the metric tensor η ij by where Using the Christoffel symbols of the second kind and the quantum geometric tensor, we define the Christoffel symbols of the first kind as Here we distinguish the first and second Christoffel symbols by the uppercase and lowercase letters for the first component, while we do not distinguish the uppercase and lowercase for other quantities. We can also define a similar quantity using the symplectic form rather than the metric tensor by [59] To write γ kij andγ kij in a unified way, we introduce We call c kij as the quantum geometric connection in analogy with the quantum metric tensor. The Levi-Civita connection does not transform like a tensor under coordinate transformations [60], which is due to the derivative acting on tensorial quantities in the definition of the Christoffel symbols of the second kind [See Eq. (29)].
The Christoffel symbols of the second kind defined on the generalized Bloch sphere γ k ij are related to those defined in the Brillouin zone where the second term shows the non-tensorial transformation property. It follows that the quantum geometric connection in the Brillouin zone has the form It is identical to the quantity defined in Eq. (23), as one can see by using ∂ a f = k f −1 f k ∂ a f k . This quantity, the quantum geometric connection, reveals the geometric nature of the low-frequency shift current in the most transparent way. Let us note that, in general relativity, the equivalence principle requires that the Levi-Civita connection does not appear directly as an observable quantity, because it is not invariant under a coordinate transformation. However, here we do not have such an equivalence principle for the Bloch vector f , so it is allowed to observe the Levi-Civita connection (or quantum geometric connection).
C. More on the geometric aspect of the shift current Equation (22) shows that the linear (circular) shift current corresponds to the imaginary (real) part of the quantum geometric connection. Thus, the circular shift current reveals the Christoffel symbol of the first kind Γ cab = On the other hand, the linear shift current is related with the Berry curvature as well as the quantum metric (through the Christoffel symbols of the second kind). Combined with the geometric property of the injection current in Eq. (18), it shows that P T -symmetric responses originate from the quantum metric only and T -symmetric responses are controlled by both the Berry curvature and the quantum metric. When the diamagnetic term w ac mn in Eq. (F7) vanishes, the relation between the shift conductivity and the quantum metric and Berry curvature can be made more direct from for ω > 0 when the diamagnetic term vanishes. This formula can be applied, e.g., to Dirac and Weyl Hamiltonians that are at most linear in momentum. Note that the real and imaginary part of the conductivity in Eq. (36) has the form of the Berry curvature dipole [61] and the quantum metric dipole [62], respectively.

D. Generalization to Multibands
Let us discuss generalizing our geometric interpretation to include three or more bands (when bands are doubly degenerate due to P T symmetry, this means that we consider six or more bands). The shift and injection conductivity takes the form n∈occ m∈unocc k I c;ab nm δ(ω − ω mn ) for ω > 0. Because the energy conservation imposed by the delta function chooses a particular set of an unoccupied band m for an occupied band n, the interband-transitive optical response is, in general, not associated with a property of the occupied band alone. On the other hand, for example, the quantum geometric tensor Q ba is defined by suming over all occupied n and unoccupied m indices of the matrix elements by n∈occ m∈unocc r b nm r a mn , and so it becomes a property of the ground state n∈occ m∈all r b nm r a mn − n∈occ m∈occ r b nm r a mn , depending only on the occupied states.
In our analysis, though, we focus on Dirac and Weyl points where the quantum geometric tensor diverges at the gapclosing points. Thus, the geometric quantities of the occupied bands are dominated by the property of the two crossing bands n = 1 and m = 2, through a large value of r b 12 r a 21 and their derivatives, and they manifest through the shift and injection currents for small ω. Similarly, when bands are Kramers degenerate due to P T symmetry, the matrix elements involving the indices for the four crossing bands are dominant contributions. It means that, at low frequencies, we have a good geometric approximation for the conductivity tensors by where Q ba and C bca are quantum geometric tensor and quantum geometric connection, respectively, defined by  where g ba = Re [Q ba ] as above. In general, the injection conductivity tensors Eq. (37) differs from the exact expression in Eq. (4) because the former has the information of all band indices rather than the specific bands n and m involved in the optical transition. Moreover, additional differences come in the shift conductivity tensors due to the virtual transition terms: C bca = −i n∈occ m∈unocc R c,a mn r b nm r a mn +virtual transitions (See Appendix J). When the optical excitation occurs very close to a gap-closing point, however, only the band indices near the gap closing significantly contribute to the geometric quantities, effectively selecting specific band indices. Also, virtual transition terms are suppressed by a factor (ω/∆E) 2 [5], where ∆E is the characteristic energy difference between the crossing bands and the other bands. Thus, Eq. (37) becomes a good approximation near the gap closing. Let us note that, in two-band or P T -symmetric four-band models, Eq. (37) becomes exact and corresponds to expressions Eqs. (18) and (22) above.
On the other hand, insulators or ordinary metals do not have geometric singularities in general, and the geometric approximation Eq. (37) does not apply to them so well. Nevertheless, recalling that the Berry curvature of each band (rather than that of the whole occupied bands) is well-defined, we can hope for a possibility of defining a well-defined geometric quantity associated with a pair of bands also. Let us see whether it makes sense to give a geometric meaning to the matrix element r a nm r b mn ≡ g ab;nm − iF ab;nm /2 by focusing on the real part g ab;nm (we note that it is different from the nonabelian quantum metric [63,64] of the occupied bands, defined by (g ab ) n1n2 = m∈unocc (r a n1m r b mn2 ). Since g ab;nm is a positive-semidefinite symmetric rank-2 tensor. i.e., g aa;nm ≥ 0 for all a for given n and m, this quantity is meaningful as a metric tensor although its interpretation is not clear. This point of view can help to understand the structure of the circular shift current. One can see that the matrix element of the circular shift current R c,a nm r b nm r a mn can be written as Γ bca;nm + virtual transitions, where Γ bca;nm is the Christoffel symbol of the first kind defined from the metric g ba;nm (See Appendix J for a derivation). Therefore, one may still regard the shift current as originating from a Christoffel-symbol-like quantity, when the virtual transitions, terms involving virtual transitions, is negligible.

V. MODEL CALCULATIONS
Our theoretical analysis reveals the circular shift current as an interesting new component of the interband bulk photovoltaic response in magnetic systems. Also, the full general-ity of our theory allows us to understand the shift and injection currents in Dirac and Weyl semimetals in arbitrary spatial dimensions in a unified way. In this section, we investigate the shift and injection current responses of tilted Dirac and Weyl points more closely with explicit calculations of the conductivity tensors. We first deal with massless Dirac and Weyl Hamiltonians up to linear order in momentum, which cover both type-I and type-II spectra in arbitrary spatial dimensions. In addition to the symmetry and divergence properties investigated above, we show the phenomenon of separation of responses, meaning that nonvanishing P T -symmetric responses and the T -symmetric responses do not coexist in the same component. It can occur generically in tilted Dirac and Weyl systems having C 2x or M x T symmetry, where x is the direction of tilting. Next, we study a model of Dirac surface states of magnetic topological insulators, which includes k 2 and k 3 terms in the Hamiltonian.

A. Tilted Weyl and Dirac semimetals: k-linear order
Let us first revisit the model of a tilted Weyl point [13,25] to understand its general second-optical response in more detail. The Hamiltonian has the form The Weyl point described by this Hamiltonian is called type-I when |v /v| < 1 and type-II when |v /v| > 1 [55]. The spectrum for the two cases are shown in Fig. 1(b,c). When light with frequency ω > 0 is illuminated, optical excitation occur when two bands have energy difference ω cv = 2 vk = ω due to the energy conservation and only the lower band is occupied, i.e., ω v = −µ + v k cos θ − vk < 0 and ω c = −µ + v k cos θ + vk > 0 (Fig. 1), where k = |k| and k x = k cos θ. When v = 0, the range of θ does not cover the whole sphere and is confined to a subspace θ − < θ < θ + in general [ Fig. 1], where θ ± are functions of dimensionless parameters 2µ/ ω and v /v. The minimal angle θ − is always zero for a type-I Weyl point, but it is typically nonzero for a type-II Weyl point [ Fig. 1(d,e)]. This asymmetric excitation leads to a nonzero optical conductivity given by The form of F c;ab (θ) for nonvanishing components are summarized in Table. II. Because of the SO (2) for the imaginary component of the conductivity. F c;ab (θ) follows the same symmetry. Remarkably, P T -symmetric responses (linear injection and circular shift) and T -symmetric responses (linear shift and circular injection) do not coexist in the same component. To understand this, let us note that P T = C 2y M y T = C 2z M z T . Since our model has C 2y T and C 2z T symmetries, nonvanishing P T -symmetric responses appear in the components that are invariant under M y and M z , whereas T -symmetric responses appear in the components that reverses sign under M y and M z . Figure 2 shows some representative components calculated from quantum geometric quantities by Eqs. (18), (19), and (36). There are some features that need to be discussed. First, as we explained in Sec. III only the circular injection current is nonvanishing at the neutral filling µ = 0 in the type-I case where |v /v| < 1. It is because, in this case, the Fermi surface is a point, so anisotropic excitation cannot occur. The circular injection response is quantized because of the quantized Berry curavature from a Weyl point, as found in Ref. [5] [See Eq. (18)]. Other responses are significant near ω = 2µ, where the anisotropic excitation occurs.
However, there are significant differences between the P Tsymmetric responses and the linear shift response. P Tsymmetric responses have peaks at ω = |2µ|, i.e., when the excitation occurs on a full hemisphere, while T -symmetric responses vanish at ω = |2µ| and changes sign there [13]. For the linear injection current, the peak at ω = |2µ| is natural because the transition of the electron velocity during the excitation from the valence band v to the conduction band c, is all positive or all negative on the hemisphere. There is no simple analogous way to understand the trend of shift current response based on the shift vector, but the vanishing of the linear shift current response at ω = |2µ| can be attributed to the T -symmetric nature. Since T symmetry requires that the current generated from one hemisphere equal to the current generated from the other hemisphere, both hemisphere should generate zero currents because linear shift currents are not generated when excitations occur on the full sphere. There are also interesting differences between the P T -symmetric responses and T -symmetric responses at µ = 0 in the type-II case where |v /v| > 1. As explained in Sec. III, emergent CP T symmetry at µ = 0 requires that the former has a vanishing response while the latter can have a nontrivial response. In other words, P T -symmetric responses do not distinguish type-I and type-II cones at µ = 0 whereas T -symmetric responses distinguish them.
As we understand a single Weyl point, it is straightforward to extend our knowledge to Dirac points in two and three dimensions. In both 2D and 3D, the protection of a Dirac point requires P T symmetry. In 2D, the Dirac Hamiltonian has the form where τ i=x,y,z are Pauli matrices for the orbital degrees of freedom, and spinless P T = σ z K symmetry forbids the mass term mσ z . In 3D, the Dirac Hamiltonian has the form H 3D = −µ − v k x + v(k x τ x + k y τ y σ x + k z τ z ), and the twofold degeneracy (Kramers degeneracy) of bands at every momentum require P T = iσ y K symmetry. Here, τ i=x,y,z and σ i=x,y,z are Pauli matrices for the orbital and spin degrees of freedom, respectively. In 3D, even in the presence of P T symmetry, two mass terms are allowed, which are m 1 τ y σ y and m 2 τ y σ z . We need threefold or fourfold rotational symmetry to protect the Dirac point in 3D by disallowing mass terms. Our massless Dirac model have continuous θ rotational symmetry around the x axis under C θ = e iθ(σx+τxσx/2) , satisfying C θ H(k)C −1 θ = H(C θ k). Keeping either threefold C 3x or fourfold C 4x rotational symmetry in crystals preserves the gap closing [52,65,66], which we assume here. Because gapless Dirac points have P T symmetry, they can only have linear injection or circular shift current responses. These responses have the same pattern shown for a Weyl point.
Here we emphasize again that while multiple responses coexist in magnetic Weyl and Dirac semimetals, each response occurs through different components of the conductivity. This helps measure each response separately. In particular, it facilitates the measurement of the circular shift current in magnetic Dirac and Weyl semimetals. Table. II and Eq. (7) shows that the current generated along the y direction under the illumination of circularly polarized light propagating in the z direction is only the circular shift currents. The photocurrent along the y direction should thus be identified with the circular shift current.

B. Dirac surface state
As an application to a more realistic model with knonlinear terms, we study the single Dirac surface state of magnetic topological insulators. Let us begin with the following effective Hamiltonian studied in Refs. [23,24].
Here, ∆ = 0 is due to spin ordering along the y direction, and it breaks M x = iσ x , rotation C 2z = −iσ z , and time reversal T = iσ y K symmetries. Since this in-plane ordering preserves C 2z T symmetry, it does not open the gap, and it just shifts the location of Dirac points by −∆/ v from the time-reversal-invariant momentum. The shifting tilts the Dirac cone because of the quadratic term: if we write (k x , k y ) = (−∆/ v + q x , q y ), the Hamiltonian has the form H = −µ−( ∆/mv)q x + v (q x σ y − q y σ x ) up to linear order in q, which is studied above. Assuming C 3z symmetry of the nonmagnetic state, we add a hexagonal warping term in order to account for the crystalline symmetry of the real system.
This term breaks M y = iσ y symmetry and C 2z T symmetry that are preserved by the spin ordering, so it opens small band gap (about 0.8 meV for parameters given below). Since T and C 2z T symmetries are both broken, all four types of shift and injection currents can be generated in this system. However, the residual M x T symmetry imposes that the separation of responses remains exact: nonvanishing components of the conductivity are σ y;xx shift,L , σ y;yy shift,L , σ x;xy shift,L for linear shift current, σ x;xx inj,L , σ x;yy inj,L , σ y;xy inj,L for linear injection current, σ y;xy shift,C for circular shift current, and σ x;xy inj,C for circular injection current.
For a numerical calculation, we take µ = 50 meV, m = 0.13m e where m e is the free electron mass, v = 2.5 eVÅ −1 , λ = 250 eVÅ −3 , and τ −1 = 1 meV and use Eqs. (18) and (22). Figure 3 show the calculated photoresponsivity κ c;ab = 2σ c;ab / 0 c, which has the dimension of the photocurrent density per unit intensity of light [57]. The peak value (occurring at ω ∼ 2µ = 100 meV) of the linear injection part is the strongest, and the others are smaller by two orders of magnitude. However, since the circular shift current grows as ω −2 while the linear injection current grows as ω −1 , the circular shift current can be comparable to or larger than the linear injection current when the peak is located below 10 meV. On the other hand, the small-frequency divergence of the linear shift current is weaker because it is due to λ = 0 [10] and is not from the tilting, and thus the peak value scales like O(ω 0 ) as µ is lowered. Therefore, the y component photocurrent generated by a circularly polarized light, j y = (κ yxx shift,L + κ yyy shift,L − 2κ y;xy shift,C )I, is dominated by the circular shift (κ y;xy shift,C ) current when ω ∼ 2µ < 100 meV. The magnitude of the linear shift current and the circular shift current can be compared in experiments since the circular parts can be measured from the current difference between the left-circularly polarized light and the right-circularly polarized light.
In experiments, the value of the observed photocurrents can be smaller than the value predicted here. For example, the photocurrents observed in Ref. [24] shows photoresponsivity of about 5 nAcm −1 W −1 at ω ∼ 250 meV, which is two orders of magnitude smaller than the calculated value here. While several factors can contribute to this reduction, one is from the cancellation between the top and bottom surfaces. This cancellation can be reduced by increasing the thickness of the sample because light attenuates more while propagating within the bulk such that the photocurrent is generated significantly on only one surface that is directly illuminated.

VI. FIRST-PRINCIPLES CALCULATIONS FOR REAL TOPOLOGICAL SEMIMETALS
For further demonstrations of our theory beyond simple two-band models, we perform first-principles calculations on the shift and injection photocurrent conductivities as well as geometric quantities of antiferromagnetic MnGeO 3 and ferromagnetic PrGeAl, respectively, as representatives of real magnetic Dirac and Weyl semimetals. We notice that such calculations on the bulk photovoltaic effects in real magnetic topological semimetal have not been reported yet despite the fact that topological semimetals are expected to be efficient infrared and terahertz photodetectors [1].

A. Antiferromagnetic Dirac semimetal MnGeO3
MnGeO 3 forms a centrosymmetric rhombohedral structure [see Fig. 9(a)] with space group R3 [67], and consequently, it would not exhibit any bulk photovoltaic effects. Interestingly, it becomes antiferromagnetic below 38 K [67] and the AF structure (magnetic space group −3 ) [see Fig. 9(a)] breaks both T and P symmetries while preserving the combined P T symmetry [42], thus leading to AF-induced bulk photovoltaic effects with linear injection and circular shift currents (see Table I). Furthermore, it was recently predicted to be a Dirac semimetal with P T symmetry-protected Dirac points (DPs) [42].
In AF MnGeO 3 , because of its P T symmetry, there are only nonvanishing circular shift photocurrent and linear injection photocurrent, as mentioned before (see Table I). Furthermore, the −3 magnetic space group admits only three nonvanishing independent matrix elements (i.e., xxz = −xzx = yzy, xyz = −xzy = −yxz and zxy = −zyx) of the circular shift conductivity tensor and six nonvanishing independent matrix elements (i.e., xxx = −xyy = −yxy, xyz = −yxz, xxz = yyz, xxy = yxx = −yyy, zxx = zyy and zzz) of the linear injection conductivity tensor [51]. Hereafter we use the shorthand notation cab for σ c;ab . We display these nonvanishing conductivity elements in the low photon energy range in Fig. 4. For simplicity, we plot only the four pronounced xxx, xxz, zxx and zzz elements of the linear injection conductivity tensor in Fig. 4. τ −1 = 10 meV is assumed. We notice that the magnitudes of the linear injection conductivity elements (σ c;ab ) are gigantic in the photon energies below 0.25 eV (Fig. 4). The magnitudes are order of τ e 3 /(2π 2 ) = 500 µA/V 2 , which are one order of magnitude larger than those in architypical polar semiconductors CdS and CdSe [68]. Circular shift photocurrents (Fig. 4) are also 10 times larger than linear shift currents in semiconductors CdS and CdSe [68]. This is remarkable because it demonstrates that the AF magnetisminduced linear injection and circular shift photocurrents, respectively, can be as large as circular injection and linear shift currents in nonmagnetic noncentrosymmetric materials. Fur- thermore, this means that the photocurrents in AF semimetals can be controlled via manipulating the magnetism with, e.g., an applied magnetic field [26]. AF MnGeO 3 hosts at least three DPs near the Fermi level E F along the k z axis, as shown in Fig. 5. In particular, there is a DP just above E F (at 1.7 meV) and being located close to the Γ point [at k DP 1 = (0, 0, −0.00364)2π/a]. This could explain the large values of the calculated photocurrents, as shown in Figs. 4(a) and 4(c). To further examine the important contributions of the DPs to the photocurrents, we also calculate the conductivity spectra with the Fermi level set to the DP2 Dirac point energy (i.e., E = 46.5 meV). DP2 is located at k DP 1 = (0, 0, 0.11112)2π/a above the k z = 0 plane in the k-space [see Fig. 5(a)]. We notice that both the shapes and magnitudes of all the conductivity spectra except conductivity element σ xxz , roughly remain the same. For ex-ample, the gigantic peak of about 1000 µA/V 2 at ∼100 meV in the Re(σ zzz ) linear injection current spectrum appears in both cases [see green curves in Figs. 4(a) and 4(b)]. Nonetheless, its sharp negative peak at 50 meV disappears in the case where the chemical potential µ is tuned to µ = E DP 2 and a sharp positive peak of larger than 700 µA/V 2 occurs at 20 meV instead. Interestingly, there is a sharp positive peak at photon energy of 5 meV in the Im(σ xxz ) circular shift current calculated by setting µ = E DP 2 . This is due to the ω −1 behavior of the shift conductivity near the Dirac point. Comparing red curves in Figs.1 (c) and (d) for Im(σ xxz ), one can also see that the large negative peak moves from photon energy 65 meV to 25 meV [see Fig. 4(d)].
While an enhanced joint density of states (JDOS), ρ(ω) = k m,n f nm δ(ω mn − ω), is a possible origin of large conductivity tensors in insulators [57], it cannot explain the peaks shown here. Figure 4(e) shows that JDOS is suppressed at low-frequencies rather than being enhanced. As we show above, linear injection and circular shift currents are, respectively, related to geometric quantities quantum metric (g ab ) and Christoffel symbols of the first kind (Γ c;ab ) at low-energies through Eq. (37). Thus, the large enhancement of conductivity tensors has geometric origin. To demonstrate this, we display g ab and Γ c;ab at µ = E F and also at µ = E DP 2 along the symmetry lines in Fig. 5. For most of the DPs, a DP is associated with a gapped DP (gDP) located approximately at the inverted position in the k-space. For example, the associated gPD for DP2 is gDP1 at k DP 1 = (0, 0, −0.11051)2π/a, and the DP2 energy level falls within the gDP1 band gap [see Fig. 5(a)]. On the other hand, there is no gDP associated with DP1. Figures 5(b) and 5(d) clearly show that for µ = E F , g ab and Γ c;ab have sharp peaks near DP1 along the k x and k y directions. g xx and g yy also peak sharply at the positions of the DP2 and gDP1 along the k z axis even though µ = E F = E DP 2 . This indicates that the gigantic linear injection and circular shift currents stem, respectively, from the large values of the quantum metric and Christoffel symbol in the vicinity of the DP1 Dirac point. Furthermore, g zz has prominent peaks in the vicinity of (but slightly away from) the DP2 and gDP1 points along the k x and k y directions, which are mainly caused by the interband transitions from the lower (occupied) Dirac cone to higher (empty) Dirac cone with transition energies of ∼0.1 eV [see Fig. 5(b)]. These g zz peaks thus give rise to the gigantic peak in linear injection current Re(σ zzz ) at around 0.1 eV [see green curves in Fig. 4(a)].
For µ = E DP 2 , as expected, g ab and Γ c;ab exhibit sharp peaks close to both the DP2 and gDP1 points [see Figs. 5(c) and 5(e)]. In particular, g zz has huge positive peaks near the DP2 and gDP1 points along all three cartesian coordinate directions, thus resulting in the prominent peak at 0.1 eV in linear injection conductivity Re(σ zzz ) [ Fig. 4(b)], as in the µ = E F case [ Fig. 4(a)]. Nevertheless, Fig. 5(e) shows that the Γ c;ab peaks near the DP2 and gDP1 have the opposite signs. This could explain why the circular shift conductivity elements Im(σ xyz ) and Im(σ zxy ) become smaller when the Fermi level is raised from E F to E DP 2 [Figs. 4(c) and 4(d)] because the contributions from the DP2 and gDP1 cancel each other to some extent. In contrast, the linear injection current elements remain almost unchanged mainly because the g ab peaks have the same signs. Nonetheless, there is a huge peak in the Γ xxz at gDP1, which is absent at DP2, because now the Fermi level falls within the gDP1 gap [see Fig. 5(e)]. This results in the large low energy peaks in the Im(σ xxz ) conductivity element [see red curves in Fig. 4(d)].

B. Ferromagnetic Weyl semimetal PrGeAl
PrGeAl forms a body-centered tetragonal structure (see Fig. 9b) with noncentrosymmetric space group I4 1 md and point group 4mm. [69] It becomes ferromagnetic at T C = 16 K [44] and the FM structure (Fig. 9b) has no T symmetry nor P T symmetry. Therefore, all four types of photocurrents may emerge in FM PrGeAl (Table I). Interestingly, it was recently predicted to be a rare ferromagnetic noncentrosymmetric Weyl semimetal [43]. Furthermore, the Weyl nodes and surface Fermi arcs in PrGeAl were observed in very recent photoemission spectroscopy experiments [44]. Thus, FM PrGeAl provides a valuable platform for studying all types of bulk photovoltaic effects in Weyl semimetals.
The crystalline point group of PrGeAl is 4mm. Thus, there are three inequivalent nonvanishing matrix elements (i.e., xxz = yyz, zxx = zyy and zzz) of linear shift conductivity and one nonvanishing element (xxz = yyz = −xzx = −yzy) of circular injection conductivity in PrGeAl above T C = 16 K. [51]. The magnetic point group of FM PrGeAl is 4m m . Consequently, in addition, there are two nonvanishing elements (i.e., xyz = −xzy = −yxz and zxy = −zyx) of circular shift conductivity and one nonvanishing element (xyz = −yxz) of linear injection conductivity in PrGeAl below T C . [51] Here, the presence of M x T and M y T symmetries forbids linear shift and injection conductivities to be simultaneously nonvanishing in the same component. The same is true for the circular polarization. The calculated conductivity spectra of these nonvanishing matrix elements are displayed in Fig. 6. We notice that all the shift current elements have large peaks below photon energy 0.25 eV [see Figs. 6(a) and 6(b)]. This is due to the ω −1 enhancement of the shift conductivity tensors at low frequencies. Indeed, the magnitudes of these peaks below 0.1 eV are comparable to those in architypical nonmagnetic Weyl semimetal TaAs [16], which also has the I4 1 md space group. Linear shift current element Re(σ zzz ) is also large in the higher energy range between 0.4 eV and 1.0 eV. This peak is not due to magnetic order because it is observed in T -symmetric responses. It is hardly related with the responses of Weyl points. The peak is far from the Fermi level, and also the conductivity components showing the large peaks is not generated by linearly dispersing Weyl points [See Table. II for example]. We do not aim to explain its origin here. The calculated circular shift conductivity elements below 0.2 eV are comparable to that of the linear shift current (Fig. 6), and remarkably, are one order of magnitude larger than those in AF Dirac semimetal  Similar to the linear shift conductivity, circular injection conductivity Im(σ xxz ) has a gigantic broad peak between 0.5 eV and 1.0 eV, and the peak value is almost one order of magnitude larger than that of linear injection conductivity in AF MnGeO 3 (see Fig. 4). It also has large peaks below 0.15 eV, although the magnitudes of these peaks are a few times smaller than the gigantic peak in the higher energy region. Linear injection conductivity Re(σ xyz ) in FM PrGeAl is also large and is about ten times larger than that of AF MnGeO 3 in the very low energy region. This one order of magnitude enhancement can be related to the presence of many Weyl points, as we analyze further below.
As in MnGeO 3 , JDOS is suppressed at low energies [ Fig. 6(d)] Thus, the origin of the large low-frequency photocurrents should be attributed to the geometric enhancement. We have calculated all four geometric quantities along the symmetry lines in the Brillouin zone [ Fig. 9(d)], as displayed in Fig. 6. We find that all four quantities have sharp peaks in the vicinity of the Σ and Σ 1 points, where there are several anticrossing nodal points with the Fermi level falling within their gaps [see Fig. 6(e)]. Since these gaps near the anticrossing points are mostly within 0.1 eV, the large peaks in the geometric quantities thus give rise to the sharp photoconductivity peaks around 0.1 eV and below. In particular, Christoffel symbols Γ xyz and Γ zxy have gigantic peaks near the Σ and Σ 1 points with the same signs, thus resulting in sharp peaks in circular shift conductivities Im(σ xyz ) and (σ zxy ), respectively, below photon energy 0.1 eV [ Fig. 6(a)].
FM PrGeAl has been reported to be a rare noncentrosymmetric FM Weyl semimetal with at least 160 Weyl points (WPs) within ±0.1 eV of the Fermi level [43]. Furthermore, there are both type I and II WPs among them. Here we consider two WPs, one type I and the other type II, and study their influences on the photocurrents. The type I WP (WP1) sits at k = (0.160, 0.204, −0.003)2π/a in the BZ and is located at 53.9 meV above E F (E WP1 ). The type II WP (WP2) is at k = (0.015, 0.255, 0.217)2π/a and its energy is  67.4 meV above E F (E WP2 ). They correspond, respectively, to Weyl points W 1 3 and W 1 2 reported in Ref. [43]. To study the influences of the WPs on the photocurrents, we calculate all the nonvanishing conductivity elements and all the geometric quantities with the Fermi level set to E WP1 and also to E WP2 .
The conductivity and geometric quantity spectra obtained by setting the Fermi level to the WP1 Weyl point energy (E WP1 ) are plotted in Fig. 7. Remarkably, compared with the µ = E F case in Fig. 6, all the photoconductivity elements within the photon energy of 50 meV increase by a factor of 2 or more. In particular, linear shift conductivity element Re(σ zxx ) is enhanced by a factor of 5 [ Fig. 7(b)]. Other changes include that the peak at ∼70 meV in circular shift conductivity elements Im(σ xyz ) and Im(σ zxy ) get significantly reduced [ Fig. 7(a)] and that the positive peak at ∼40 meV in linear shift conductivity element Re(σ zzz ) changes to the negative peak [ Fig. 7(a)]. Nonetheless, the features of all the conductivity elements above 100 meV remain essentially unchanged.
The changes caused by moving the Fermi level to the WP1 WP energy mentioned above, can be largely explained by Here q = k − kWP1 denotes the momentum displacement from the WP1 point.
to the five-fold increased linear shift conductivity Re(σ zxx ). Note that although in the µ = E F case Γ zxx also has sharp peaks near the Σ and Σ 1 points, these peaks have both positive and negative signs [see Figs. 6(i)] and thus their contributions to the linear shift conductivity cancel each other to some extent.
The results obtained by setting the Fermi level to the energy of the type II WP2 Weyl point are plotted in Fig. 8. We notice that as for the WP1 WP case, all the photoconductivity spectra remain more or less unchanged for photon energies larger than 100 meV. As for the WP1 WP case, the peak height of linear shift conductivity Re(σ zxx ) at very low energy of ∼10 meV gets doubled, although the sign of the peak changes from positive to negative [see Fig. 8(b)]. However, in contrast to the WP1 Weyl point case, many low energy peaks especially of shift current conductivity elements, become smaller [Figs. 8(a) and 8(b)]. In particular, the peak at 10 meV in both circular shift elements Im(σ xyz ) and Im(σ zxy ) decreases by a factor of nearly 2, and also the peak sign changes to the opposite. These different changes in the low energy photoconductivity spectra betwen the type I and type II Weyl point cases may be attributed to their different energy dispersions and hence the different distributions of the four geometric quantities near the Weyl points (see Figs. 7 and 8). Figure 7(d) shows that in the type I WP1 case, the upper and lower Weyl cone bands lie, respectively, above and below the WP1 energy. As a re-   8(d)]. This hinders the low energy inter-Weylband transitions along these two directions when E F is set to E WP2 , and thus reduces the peaks in the geometric quantities near the WP2 point, thereby leading to several smaller peaks in the photocurrents below ∼20 meV. Of course, in a real topological semimetal such as the present system, situation could be more complicated. Figure 8(

VII. DISCUSSION
Let us discuss some issues related to the smallness of the frequency. We note that our theory is reliable for the photon frequency above 1 THz because we assume ω Γ. Since the typical relaxation time in semimetals is one picosecond, ω/2π > 1 THz should be taken for our theory to apply. Thus, the divergence of the responses in our model should not be interpreted as a physical divergence, and it is cut off at ω ∼ Γ. Moreover, intraband (i.e., non-transitive) secondorder optical responses exist in magnetic systems, and they become comparable to the interband responses when ω < ∼ Γ. To estimate the magnitude of the intraband responses, let us consider the semiclassical second-order optical response σ c;ab It has the same symmetry property as the linear injection current (so it appears in magnetic systems where time-reversal symmetry is broken) and scales as ω −2 µ 2−d near tilted Dirac or Weyl points. Since the ω −2 factor enhances the semiclassical response at low-frequences, it becomes comparable to the peak value of the linear injection response at ω < ∼ τ −1 = Γ if we take µ ∼ ω.
Another thing to care about at low frequencies is the validity of the perturbation theory. We suppose that the photovoltaic responses are dominated by the second-order responses since the electric field is generally weak enough. However, at small frequencies, A = E/ω becomes large, so higher-order responses can become significant. According to the Floquet theory analysis in Ref. [4], perturbation theory works well for ω −1 eEv Γ, where we assume that the interband velocity matrix element has the same order of magnitude as the intraband velocity. It sets a lower bound for ω.
For v ∼ 1 eVÅ −1 and Γ ∼ 1 THz, the lower bound is ω min ≈ I/(W/cm 2 ) THz. While this bound is very small for typical lasers with intensity I < 1 W/cm 2 , it should be taken into account for an analysis of the experimental results obtained by using high-intensity pulse lasers.
In the non-perturbative regime, light absorption occurs faster than the interband relaxation, so electrons are excited until the system reaches a new equilibrium where the absorption and emission of light finds a balance. The relaxation plays an important role here because a mismatch between the absorption and emission, which generate opposite photocurrents, requires a finite relaxation [4]. Even in a very clean system with an extremely small Γ, the low-frequency divergent behavior we discuss above is cut off by this reason.
At finite temperature, thermal fluctuation can significantly reduce low-frequency photovoltaic response in tilted-Dirac or Weyl semimetals [13]. When the peak value is considered, the relevant frequency scale where the peak value is reduced to half is 2|v/v |k B T / for type I and 2k B T / for type II. It sets a quite high lower-bound, since 2k B T / =12 THz at room temperature (cf. 0.17 THz at 4.2K and 3.2 THz at 77 K). When the tilting is small such that |v /v| < 0.4, the cutoff scale is larger than 30 THz at room temperature, which is the highest edge of the terahertz radiation. Therefore, cooling will be needed in order for tilted-Dirac or Weyl semimetals to work as an efficient terahertz photodetector.
Throughout this work, we focus on the DC generation. However, let us note that our results can also be applied to the second harmonic generation, where an alternating current of frequency 2ω is generated by a uniform illumination of light with frequency ω. Since the second harmonic gen-eration is associated with the shift vector [4,13], its lowenergy divergence has the same form as that for the shift current [13], and thus it can be related to the quantum geometric connection. We believe that similar geometric interpretations of other quantities are also possible. It will be an interesting topic to explore quantum geometric properties of the third-order optical conductivity. Since it has four components, we expect the existence of a third-order optical response controlled by the Riemann curvature tensor R a bcd = ∂ c Γ a db − ∂ d Γ a cb + e (Γ a ce Γ e db − Γ a de Γ e cb ). Note added.-Recently, Hikaru Watanabe and Yoichi Yanase studied the bulk photovoltaic responses in magnetic systems independently of ours. They also figure out the circular shift current and discuss its enhancement near gap-closing points. Our results are consistent with the results in their preprint when there is an overlap.
• −e is the electron charge. e > 0. • Geometric quantities in momentum space (indices a = 1, 2, 3 are for the (k 1 , k 2 , k 3 ) basis) g ab metric tensor F ab /2 symplectic form Q ab = g ab − iF ab /2 quantum geometric tensor Γ a bc Christoffel symbol of the second kind Γ abc Christoffel symbol of the first kind Γ abc symplectic Christoffel symbol C abc = Γ abc − iΓ abc quantum geometric connection (A2) • Typical quantities in second-order optical conductivity.
− v a nn , R c;a mn = r c mm − r c nn + i∂ c log r a mn , r a mn = i m|∂ a |n , r a mn;c = ∂ c r a mn − i(r c mm − r c nn )r a mn . (A3)

Appendix B: Frequently used identities
In the appendix we use the following two identities.
where we integrate by parts in the first line. The following is a specific example of the above identity with O = H/ , and it is used very often in the following appendices.
where M m←n is the probability of the transition from n to m. We note that this Fermi Golden rule calculation gives an overall (−1) sign that is absent in the result of in [19,22] while they also take −e as the electron charge. Also, our expressions differ from the result in [19,22] by the sign of ω in the electric fields.
Appendix E: Derivaion of Eqs. (9) and (10) Let M be a point-group symmetry operation: Then, r a mn (k) → r a mn (k) = g ab r b mn (M −1 k).
We can see this from where we use that where H (k) =MH(k)M −1 , and H(k) = e −ik·r He ik·r is the Bloch Hamiltonian. From these properties, one can find that where σ c;ab DC is defined from the M-transformed state |u nk and Hamiltonian H (k).
We can repeat this process for time reversal |u nk → |ψ nk =T |u n,−k . (E7) For the dipole matrix elements, where we use the antiunitary property of time reversal in the third line, and we integrate by parts in the fourth line.
We have To see this, let us take occupied bands n 1 and n 2 . Then, m∈unocc r a n1m r b mn2 = m∈unocc n 1 |i∂ a |m m|i∂ b |n 2 = m∈unocc ∂ a n 1 |m m|∂ b n 2 = ∂ a n 1 |∂ b n 2 − p∈occ ∂ a n 1 |p p|∂ b n 2 = 1 2 ∂ a n 1 |∂ b n 2 − p∈occ ∂ a n 1 |p p|∂ b n 2 + (a ↔ b) where the symmetric part is the nonabelian quantum metric of the occupied states, and the antisymmetric part is the nonabelian Berry curvature of the occupied states, where r a n1n2 = i n 1 |∂ a |n 2 is the nonabelian Berry connection of the occupied states. It follows that

Dirac Hamiltonian
Here, we consider the following d M × d M Dirac Hamiltonian Let us note that this kind of Hamiltonian describe genetric two-band systems when d M = 2 and generic P T -symmetric four-band systems with (P T ) 2 = −1 when d M = 4. Let us suppose that half of the bands are occupied. Then, since ω unocc − ω occ = 2 i f 2 i = 2f , we have Also, the velocity matrix elements are given by It follows that The first term is where we use the property of Gamma matrices {Γ i , Γ j } = 2δ ij and define As for the second term, we note that for m, n ∈ occ, One can see this as follows. First, if we dot product m|Γ j |n with f , we have m|Γ j |n f j = m|f j Γ j |n = −f δ mn . On the other hand, for a vector n perpendicular to f , m|Γ j |n n j = m|n j Γ j |n = 0 because n j Γ j anticommutes with f i Γ i . Since m|Γ j |n is parallel to f i , we have Eq. (H13). To sum up, we have In two-band systems, we have J ij = − ijkfk because [Γ i , Γ j ] = 2i ijk Γ k , where Γ i = σ i is a Pauli matrix. On the other hand, in four-band systems described by a Dirac Hamiltonian, which are P T -symmetric four-band systems, we have J ij = 0.
Appendix I: Compatibility between metric, symplectic form, and connection

Metric connection
A connection γ k ij is called metric compatible if the inner product done by the metric is invariant under parallel transport [60]: In components, it reads This can be seen like this. 0 = (∇ k η) ij = e i · ∇ k η · e j = ∇ k (e i · η · e j ) − ∇ k e i · η · e j − e i · η · ∇ k e j = ∂ k η ij − γ l ki η lj − η il γ l kj , where we use the notation η ij = e i · η · e j and the definition ∇ k e i = γ l ki e l . A metric compatible connection is called a metric connection. If a metric connection γ k ij satisfies the torsion-free condition γ k ij = γ k ji (torsion is the antisymmetric part of γ k ij ), it is uniquely determined to be and it is called the Levi-Civita connection.

Symplectic connection
Similarly, a connection is called a symplectic connection when it satisfies In components, it reads ∂ k ij − il γ l kj + jl γ l ki = 0.
Unlike the metric compatibility, this condition alone cannot uniquely determine the connection even after imposing the torsion-free condition γ k ij = γ k ji . However, adding metric compatibility makes the connection unique. On the generalized Bloch sphere, the Levi-Civita connection is the unique connection that is compatible with both the metric and the symplectic form. nm r a mn,c + r b nm,c r a mn = r b nm ∂ c r a mn − ir b nm r a mn (r c mm − r c nn ) + (∂ c r b nm )r a mn − ir b nm r a mn (r c nn − r c mm ) = ∂ c (r b nm r a mn ), we obtain Γ bca;nm = 1 2 (∂ c g ba;nm + ∂ a g bc;nm − ∂ b g ca;nm ) = 1 2 Re ∂ c (r b nm r a mn ) + ∂ a (r b nm r c mn ) − ∂ b (r c nm r a mn ) = 1 2 Re r b nm r a mn,c + r b nm,c r a mn + r b nm r c mn,a + r b nm,a r c mn − r c nm r a mn,b − r c nm,b r a mn = 1 2 Re 2r b nm r a mn,c + r b nm (r c mn,a − r a mn,c ) − (r c nm,b − r b nm,c )r a mn − r c nm (r a mn,b − r b mn,a ) , where we add 0 = r b nm r a mn,c − r b nm,c r a mn − (r b nm r a mn,c − r b nm,c r a mn ) in the last line and use Re(r b nm,a r c mn ) = Re(r c nm r b mn,a ) for the last term. Let us note that r c mn,a − r a mn,c terms in the parenthesis are virtural transitions as can be seen from Eq. (F6). Thus, we have Γ bca,nm = Re r b nm r a mn,c + virtual transitions.
Equivalently, we find that the matrix element of the circular shift current is given by Re r b nm r a mn,c − r b nm,c r a mn = Γ bca;nm − Γ acb;nm + virtual transitions.
An analogous proof for the complex part will be more challenging because it requires the calculation of the inverse quantum metric becauseΓ bca = 1 2 F bd (g −1 ) de Γ eca . Manganese germanate MnGeO 3 forms a rhombohedral ilmenite structure with a P symmetric space group R3 [67] (see Fig. 9a). The unit cell contains two formulae units (f.u.). It becomes antiferromagnetic (AF) below 38 K with the same unit cell. [67]. Although the AF structure (Fig. 9a) breaks both T and P symmetries, AF MnGeO 3 with magnetic space group −3 has the combined P T symmetry. [42]. PrGeAl crystallizes in the body-centered tetragonal structure (see Fig. 9b) with a broken P symmetry space group I4 1 md. [69] The unit cell also contains two f.u. It becomes ferromagnetic (FM) at 16 K [44] and the FM structure (Fig. 9b) has no T symmetry nor P T symmetry.
The electronic band structure and magnetic properties of MnGeO 3 and PrGeAl are calculated based on first-principles density functional theory with the generalized gradient approximation (GGA) [70]. The experimental structural parameters for MnGeO 3 [67] and PrGeAl [69] are used in the present calculations. To better describe the Coulomb correlation among Mn 3d electrons and also among Pr 4f electrons, we adopt the GGA+U scheme [71]. Following the recent studies [42,43], we use the effective U value of 4.0 eV for both Mn 3d and Pr 4f electrons. The calculations are performed using the accurate projector-augmented wave (PAW) method [72], as implemented in the Vienna ab-initio simulation package (VASP) [73,74]. The fully relativistic PAW potentials are adopted in order to include the spin-orbit coupling (SOC) effect. Large plane-wave cutoff energies of 450 eV and 500 eV are used for MnGeO 3 and PrGeAl, respectively. For the Brillouin zone integration, k-point meshes of 12 × 12 × 12 and 16 × 16 × 16 are used for MnGeO 3 and PrGeAl, respectively. All the calculations are performed with an energy convergence within 10 −6 eV between the successive iterations.
Here we consider MnGeO 3 and PrGeAl, respectively, in the AF and FM states with the magnetic moments being par- allel to the c-axis. The calculated relativistic band structures of MnGeO 3 and PrGeAl are displayed in Fig. 10(a) and Fig. 10(b), respectively. They are both semimetals with low density of states at the Fermi level (E F ) of 0.315 states/eV/f.u. and 0.133 states/eV/f.u., respectively. The calculated Mn spin magnetic moment in MnGeO 3 is 4.29 µ B and that of Pr in PrGeAl is 1.89 µ B . The calculated band structures agree well with that reported previously [42,43]. Nonlinear optical photocurrents are calculated based on the linear response formalism with the independent particle approximation, as described above. Specifically, the DC shift and injection photocurrent conductivity tensors are calculated using Eq. (4) given in Sec. II. Since a large number of kpoints are needed to get accurate NLO responses [75,76], we use the efficient Wannier function interpolation method based on maximally localized Wannier functions (MLWFs) [77][78][79]. For MnGeO 3 , 68 MLWFs per unit cell of Mn d, Ge p and O p orbitals are constructed by fitting to the GGA+U+SOC band structure. In PrGeAl, 76 MLWFs per unit cell of Pr d and f orbitals as well as Ge p and Al sp orbitals are adopted. The band structures obtained by the Wannier interpolation are nearly identical to that from the GGA+U+SOC calculations (see Fig. 10). The shift current conductivity tensors are then evaluated by taking very dense k-point meshes of 200 ×200 × 200 for MnGeO 3 and of 160 ×160 × 220 for PrGeAl. We find that the conductivity tensors obtained using such dense k-point meshes converge within a few percent. Here we consider the "cold" semimetals, i.e., the Fermi-Dirac function in Eq. (4) is taken to be a step function. Furthermore, the Dirac δ function is replaced by a Lorentzian function with broadenning width of τ −1 = 10 meV. Accordingly, in the injection conductivity calculations, we use the relaxation time τ given by τ −1 = 10 meV.