Pure spin photocurrent in non-centrosymmetric crystals: bulk spin photovoltaic effect

Spin current generators are critical components for spintronics-based information processing. In this work, we theoretically and computationally investigate the bulk spin photovoltaic (BSPV) effect for creating DC spin current under light illumination. The only requirement for BSPV is inversion symmetry breaking, thus it applies to a broad range of materials and can be readily integrated with existing semiconductor technologies. The BSPV effect is a cousin of the bulk photovoltaic (BPV) effect, whereby a DC charge current is generated under light. Thanks to the different selection rules on spin and charge currents, a pure spin current can be realized if the system possesses mirror symmetry or inversion-mirror symmetry. The mechanism of BSPV and the role of the electronic relaxation time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document}τ are also elucidated. We apply our theory to several distinct materials, including monolayer transition metal dichalcogenides, anti-ferromagnetic bilayer MnBi2Te4, and the surface of topological crystalline insulator cubic SnTe.


Derivations of NLO spin current conductivity
In this section we derive the NLO spin current conductivity from quadratic response theory, in a similar fashion to Refs. [1,2]. The Hamiltonian of the system can be written as where H 0 is the unperturbed Hamiltonian, while V is a perturbation, from e.g. the interaction between light and the electrons. Without the interaction term V , the equilibrium density matrix should be Note that [ρ 0 , H 0 ] = 0. When V is turned on, the equation of motion of the density matrix ρ(t) should be (von Neumann equation) The last term − ρ−ρ0 τ is a dissipation that describes the interaction between the system and the heat bath: the system always has the tendency to return to ρ 0 . Here we adopted the constant relaxation time approximation and use a uniform τ for all states. In practice, τ should be different for different states (e.g., electrons and holes may have different τ ). More rigorously, one should use ∂ρ ∂t | col , which incorporates all relaxation processes, such the scattering of the electrons/holes with phonons, impurities, etc, and the spin relaxation process. But we made simplifications by using ∂ρ ∂t | col ≈ − ρ−ρ0 τ . In order to solve Eq. (S3), we use a trick similar to the transformation between the Schrodinger picture and the interaction picture. Letρ(t) = e t τ e i H 0 t ρ(t)e −i H 0 t , then we have = · · · (S5) By iteratively puttingρ(t) into the bracket on the rightmost of the equation above, we can obtaiñ ρ(t) =ρ (0) (t) +ρ (1) (t) +ρ (2) (t) + · · · (S6) Then going back fromρ(t) to ρ(t) with ρ(t) = e − t τ e −i H 0 tρ (t)e i H 0 t , we have . Then ρ (1) can be calculated as Obviously, Then, the second order ρ (2) We have For an arbitrary operator θ, the thermal expectation value of θ should be θ = Tr(θρ) (S15) The first order response is The the second order response is The last equity can be obtained with an interchange of dummy variables as (n → l, l → m, m → n).
For the nonlinear spin current effect, the interaction is where A is the vector potential while E is the electric field. After a second quantization, we have V nm (ω) in the basis of Bloch waves For the spin current in the a direction with spin component i, we should have θ nm = j a,i nm = 1 2 {v a , s i } nm , where s i = 2 σ i is the spin operator.
Finally, the nonlinear direct spin current is And the spin conductivity can be written as where we have replaced E mn with ω mn . Eq. (S21) is exactly the same as Eq. (2) in the main text.
2 Reproduce the shift and injection charge current formulae 2.1 NLO charge conductivities in the length gauge In Sec. 1, we adopted the so-called velocity gauge in deriving the NLO spin/charge conductivity, where the interaction with the oscillating electric field is described by the vector potential and velocity as in Eq. (S18). A charge current version of Eq. (S21), where one sets j a, has been derived in Refs. [1,2] and has been used in e.g. Refs. [3,4]. On the other hand, a different version of the formulae [5] for calculating the NLO charge conductivity, which uses the length gauge, seems more popular in literatures. In the length gauge, the interaction with light is described by where r is the position (length) operator. We first give a brief description on the formulae in length gauge, and then show that they are equivalent with our Eq. (S22).
In a time reversal symmetric system, the NLO charge current is divided into two parts, namely the shift current and the injection current, The shift current conductivity is is the generalized derivative of position operator r with respect to k and ξ a nn = i u nk |∂ k a |u nk , where |u nk is the cell-periodic part of the wavefunction |nk . Note that for Bloch waves, it is not straightforward to define the full matrix of the position operator r mn = m|r|n , since the Bloch waves are infinite in space. Usually r mn is defined with only the interband elements But the definition of ξ is valid for intraband elements with m = n as well.
Under a linearly polarized light b = c, Eq. (S24) can be further simplified as The injection current response function is where ∆ a nm = v a nn − v a mm and [r c mn , r b nm ] = r c mn r b nm − r b mn r c nm . The shift current can be induced by a linearly polarized light (b = c), and is a static current. On the other hand, the injection current cannot be induced by a linearly polarized light, and it should grow with time, thus it is somewhat like injected into the system. If the carrier lifetime is τ , then there should be a dissipation term −j injection /τ . And in the static limit (t → ∞), so the effective conductivity should be defined as τ η a bc (0; ω, −ω). Eqs. (S24, S28) are derived in Ref. [5] and is widely used in many papers.
Here we would like the remark that the physical mechanism of the shift current and injection current are more evident in the length gauge format. Eqs. (S27, S28) are reminiscent of the Fermi's golden rule. The Dirac functions indicate that light induces resonant interband transitions, and the transition rate are proportional to |r b mn | 2 for LPL, and [r c mn , r b nm ] for CPL. Then the shift current comes from the fact that the wavefunction centers of electrons and holes differ by a factor of R a nm , which leads to a electric dipole eR a nm . The time derivative of the electric dipole is the current. On the other hand, the injection current comes from the fact that the electrons and holes have different velocities, and the net current is determined by the velocity difference ∆ a nm = v a nn − v a mm . In order to show the equivalence, the first step is to factorize the denominator of Eq. (S22) with Sokhotski-Plemelj Formula. That is, in the limit of τ → ∞ (δ = 1/τ → 0), one has 1

Equivalence between the velocity gauge and the length gauge
here δ(x) is the Dirac delta function and P stands for the Cauchy principal value in k integration.
The next step is to rearrange the numerator of Eq.
according to their behavior under T operation. As discussed in the main text, one has T N 0abd (k) = . Consequently, after a summation over ±k, only the imaginary part of N 0abd (k) would contribute to the final result, and thus we can ignore the real part of N 0abd (k) and treat it as a purely imaginary quantity.
Reproduce the shift current. First, we note that the a dc current should be a real quantity. 1 Here we focus on the first term in the bracket of Eq. (S22), the second term can be analyzed in the same fashion.
For a LPL, E b and E c has no phase difference and only the real part of the conductivity σ a bc contributes to the NLO current. As discussed above, the numerator N 0abd (k) can be treated as a purely imaginary quantity, thus one needs to consider the imaginary part of the denominator, which is The second term in symmetric with respect to the permutation of m, n (after summation over m, n), while the imaginary part of v a nl v b lm v c mn is asymmetric with respect to m, n. Therefore only the first term in Eq. (S31) would not vanish in the final result. Now it's straightforward to check that that the integrand in Eq. (S22) is terms with format where we have replaced ω in the prefactor of Eq. (S22) with ω ml because of the the delta function.
On the other hand, the integrand of Eq. (S24) is The equivalence between Eq. (S32) and Eq. (S33) can be proved by using Eq. (S26) and plugging in the Note that the first term in Eq. (S34), F b lm;a = r a lm ∆ b ml +r b lm ∆ a ml ω lm , does not contribute to the final result because f lm r c ml F b lm;a is asymmetric in m, l. Reproduce the injection current. For a circularly polarized light, E b and E c should have a phase difference of i. In order to get a real current, we need to consider the imaginary part of σ a bc . In this case we need to examine the real part of the denominator of Eq. (S22), which is One can see that the first and second term in Eq. (S35) corresponds to m = n and m = n, respectively.
In case that τ → 0, the contribution from the first term is much smaller than the second term and thus we only consider the second term.
Putting m = n (n = l) in the first (second) term of Eq. (S22), one has the integrand as In the last step we have taken the asymmetry part (bc ↔ −cb) of the integrand. One can see that Eq.
(S36) is the same as the integrand of Eq. (S28) plus a τ factor, as expected.

Shift-and injection-like mechanisms for the spin current
As discussed in Sec. 2, the well-known shift and injection formulae in a T -conserved system can be reproduced from the charge current part (j a,s 0 = ev a ) of Eq. (S21). In this section, we would show that for the spin current with j a,s i = 1 2 {v a , s i }, i = 0, Eq. (S21) can also be broken down into shift-like and injection-like parts. By shift-(injection-) like, we mean that the conductivity scale with τ 0 (τ 1 ), where τ is the carrier relaxation time.

P-broken, T -conserved systems
We first study a T -conserved system. As in Sec. 2, we still need to factorize the denominator of Eq.
(S21) with Eq. (S30), and the arrange the numerator according to its behavior under T operation. Unlike N 0abd (k), for the spin current, one has T N iabd (k) = N iabd * (−k). Consequently, after a summation over ±k, only the real part of N iabd (k) contributes to the total conductivity, in contrast to the charge current, where the imaginary part of N 0abd (k) contributes.
where we have taken the symmetry part (bc ↔ cb) as r b lm , r c ml = r b lm r c ml + r c lm r b ml because we are dealing with LPL. One can see that for a T -conserved system under LPL, the mechanism for spin current generation in actually injection-like. A prominent feature is that the NLO spin conductivity is (approximately) proportional to the lifetime τ .
Circularly polarized light. Under LPL, one needs to keep the imaginary part of the denominator as in Eq. (S31). This time the second term in Eq. (S31) would also contribute to the final result because the real part of j a,s i nl v b lm v c mn is also symmetric in m, n. As a result, under CPL, the contribution to the spin current has a shift-like part, which depends linearly on τ .

P-broken, T -broken, PT -conserved systems
For a PT -conserved system, one needs to study the imaginary part of N iabd (k) because PT N iabd (k) = −Ñ iabd * (−k). With similar analysis as in the previous section for T -conserved system, one can find that the generation of spin current under LPL (CPL) has shift (injection)-like mechanism, respectively.
The same analysis also applies to charge current in a PT -conserved system, and the results are listed in Table 2 of the main text.

P-broken, T -broken, PT -broken system
In this kind of systems, both the real and imaginary parts of N iabd (k) can contribute to the conductivity for spin and charge currents. As a result, under either LPL or CPL, the shift-and injection-like mechanisms would both contribute.

The definition of spin current
There are lots of debates on the definition of the spin curren. In the main text we adopted the conventional definition j 1 = 1 2 (vs + sv). But this definition is not entirely correct, although it is convenient, physically appealing, and extensively employed in many works until today. For example, this spin current is not conserved. Also, as suggested in Rashba's work [9], this definition would lead to a non-zero spin current even if an inversion asymmetric insulator is in equilibrium (under no electric field, light, etc.). There are lots of debates, and there are also works claiming that we do not need to modify this definition. See e.g., Ref. [10].
In Ref. [11], the authors proposed another definition, which is dt . This definition can be conserved in some systems. However, it requires that "spin generation in the bulk is absent due to symmetry reasons". In other words, it requires the bulk integration Here we roughly estimate the difference between the spin current defined with j 1 = 1 2 (vs + sv) and j 2 = d(rs) dt . Compared with j 1 , j 2 has an additional term that comes from the torque on the spins [11,12]. This term is proportional to [H, s]. Therefore, the relative difference j2−j1 j1 can be roughly H · s i , where · indicates matrix norm. This can be naively understood in the following way. We have where {a, b} = ab + ba ensures hermiticity. The first term is just j 1 , while the second term comes from the torque on the spins.
The ratio between these two terms is (very roughly) H · s i . Rigorously speaking, the position operator r needs extra care in infinite solid-state systems.
We have calculated and plotted α z in the first Brillouin zone for MoS 2 ( Figure S1a) and MnBi 2 Te 4 ( Figure S1b). One can see that α is on the order of 0.1 ∼ 0.2. From this point of view, one may deduce that the difference between j 1 and j 1 is indeed not negligible, but in general cases, it would not qualitatively change the main results.  For spin current generation, the temperature rise due to the energy dissipation is generally not significant. Here we take monolayer MoS 2 as an example. For the bulk spin photovoltaic effect studied in this work, the main energy consumption is the photon absorption due to interband transitions (electron-hole pair generation), and the absorbance is c 0 d, where 0 is the vacuum permittivity, i is the imaginary part of the dielectric function, σ r is the real part of the optical conductivity. d is the thickness of the material, which is taken as 0.6 nm for MoS 2 . The energy consumption rate per unit area is where I = 0 c 2 E 2 is the intensity of the light. From our ab initio calculations, at ω = 3 eV one has σ r (ω) = 4 × 10 5 Ω/m. In the main text, we showed that light with electric field strength on the order of E = 100 V/m would be able to generate a detectable spin current. With this field strength the energy absorption power is P = 1.2 × 10 −4 W/cm 2 , which is rather small. Here we calculate the temperature rise under an electric field of E = 1 MV/m, which is much larger, but readily available with laser technology. At this field strength, one has P = 1.2 × 10 4 W/cm 2 . Assume that MoS 2 is put on a substrate with thermal conductivity κ and thickness l subs . If a continuous wave (CW) light is used, then the steady steady-state temperature rise can be roughly estimated from Assuming that l subs = 1 µm and κ = 100 W · m −1 · K −1 , then ∆T CW ≈ 1.2 K, which is not significant.
On the other hand, if a pulse laser is used, then the temperature rise can be estimated from Where S is the area of a unit-cell, k B is the Boltzmann constant, while τ pulse is the duration of the pulse, and is taken as 1 ps here. One can find that ∆T pulsed ≈ 0.7 K, which is also not significant as well.
We would like to note that, not all energies absorbed by the materials goes to the phonon (lattice) system. They may be re-emitted as photons when e.g., electrons and holes recombine. Therefore, the temperature rise in the ion system may be even lower than ∆T CW and ∆T pulsed estimated above.

Maximum spin current density
With the MoS 2 parameters, when the external field strength is E = 1 MV/m, the spin current density is 10 8 A m 2 2e , and the temperature rise in the sample is estimated to be on the order of 1 K, which is not high. If we further increase the electric field strength, then two main issues will arise. 1) On the experimental side, a strong laser may destroy the sample. This might happen for electric field strength above 10 MV/m when one uses a continuous wave laser. In this case, the spin current density is around 10 10 A m 2 2e , which is quite large, while the temperature rise can be as high as 100 K. Note that the temperature rise can be mitigated with better thermal management. On the other hand, the sample may survive in an even stronger electric field if the pulsed laser is used. For example, with a femtosecond laser, the electric field can be as high as 100 MV/m, with a temperature rise of In summary, with a pulsed laser, our theoretical picture may work up to an electric field strength

Detection of the spin current
In the main text we discussed how the spin current can be detected. The schemes are illustrated in Fig. S2.

Computational benchmarks
In this section we provide computational benchmarks on our NLO spin/charge conductivity formula Eq. (S21). b. We reproduce the NLO charge conductivity for monolayer GeS, which was studied in Ref. [13]. Our results are shown in Fig. S4 . One can see that it is in great agreement with Fig. 2 in Ref. [13].

NLO conductivity as a function of lifetime τ
The analysis in Sec. 3 demonstrates that under LPL, the charge (spin) currents of MoS 2 are shift (injection) like, while the charge (spin) currents of MBT are injection (shift) like (see also Table 2 in the main text). For shift (injection) like mechanism, the conductivity scales as τ 0 (τ 1 ). We have numerically tested the τ -dependence for MoS 2 and MBT. The results can be found in Figs. S5 and S6, respectively. The expected dependence on τ can be clearly observed.  The spin conductivity is approximately linearly dependent on τ , while the charge conductivity remains a constant. Such behavior manifests that in a T -symmetric system, the charge current comes from shift-like mechanism, while the spin current comes from injection-like mechanism under LPL. The conductivities are at ω = 3.8 eV. Figure S6: NLO spin (a) and charge (b) conductivity of MBT under LPL as a function of carrier relaxation time τ . The spin conductivity is approximately independent of τ , while the charge conductivity is approximately linearly dependent on τ . Such behavior manifests that in a PT -symmetric system, the spin current comes from shift-like mechanism, while the charge current comes from the injection-like mechanism under LPL. The conductivities are at ω = 0.6 eV.

NLO conductivity as a function of SOC strength
In the main text, we argued that in MBT, the NLO charge conductivity should increase with the SOC strength λ. This is computationally verified in Fig. S7. Figure S7: Peak values of NLO spin (a) and charge (b) conductivity of MBT as a function of SOC strength. As discussed in the main text, the spin conductivity is nonzero even when SOC is completely turned off (λ = 0), while the charge conductivity is zero when λ = 0 due to inversion-spin rotation symmetry and grows with λ. Note that SOC strongly modifies the band structure of MBT, thus the spin and charge conductivity have relatively more complicated dependence on λ, as compared that in the case of MoS2.

NLO spin current under circularly polarized light
In the main text, the NLO spin currents under LPL are presented. They correspond to the symmetric real part of Eq. (S21). In this section we show the spin current under CPL as in Figs.