Nonlinear spin current generation in noncentrosymmetric spin-orbit coupled systems

TL;DR

This paper introduces a theoretical mechanism for generating nonlinear spin currents in noncentrosymmetric spin-orbit coupled systems using electric fields, with potential for enhanced efficiency and control in spintronics applications.

Contribution

It proposes a new nonlinear spin current generation mechanism based on electric field application, supported by Boltzmann transport theory, applicable to topological insulators and Rashba/Dresselhaus systems.

Findings

Nonlinear spin current proportional to E^2 induced by electric fields.

Enhanced spin current magnitude compared to previous methods.

Rectification of spin current with AC electric fields like terahertz light.

Abstract

Spin current plays a central role in spintronics. In particular, finding more efficient ways to generate spin current has been an important issue and studied actively. For example, representative methods of spin current generation include spin polarized current injections from ferromagnetic metals, spin Hall effect, and spin battery. Here we theoretically propose a new mechanism of spin current generation based on nonlinear phenomena. By using Boltzmann transport theory, we show that a simple application of the electric field $\bf{E}$ induces spin current proportional to $\bf{E^2}$ in noncentrosymmetric spin-orbit coupled systems. We demonstrate that the nonlinear spin current of the proposed mechanism is supported in the surface state of three-dimensional topological insulators and two-dimensional semiconductors with the Rashba and/or Dresselhaus interaction. In the latter case, the angular dependence of the nonlinear spin current can be manipulated by the direction of the electric field and by the ratio of the Rashba and Dresselhaus interactions. We find that the magnitude of the spin current largely exceeds those in the previous methods for a reasonable magnitude of the electric field. Furthermore, we show that application of AC electric fields (e.g. terahertz light) leads to the rectifying effect of the spin current where DC spin current is generated. These findings will pave a new route to manipulate the spin current in noncentrosymmetric crystals.

Figures

Equations

$$-eE\frac{\partial f}{\partial k_x}=-\frac{f-f_0}{\tau },$$

$$f_1=e\tau E\frac{\partial f_0}{\partial k_x}=e\tau E\frac{\partial \epsilon }{\partial k_x}\frac{\partial f_0}{\partial \epsilon },$$

$$f_2=e\tau E\frac{\partial f_1}{\partial k_x}=(e\tau E{)}^2\frac{{\partial }^2{f}_0}{\partial {k_x}^2}.$$

$$j_{\mu ,s_{\nu }}^{(n)}=\sum _I\int \frac{{\mathrm{d}}^2\mathbf{k}}{(2\pi {)}^2}⟨I,\mathbf{k}\left|{\hat{ȷ}}_{\mu s_{\nu }}\right|I,\mathbf{k}⟩f_n^I,$$

$$j_{\mu ,s_{\nu }}^{\text{odd}}$$

$$j_{\mu ,s_0}^{\text{even}}$$

$$j_{\mu ,s_{\nu }}^{\text{even}}$$

$$j_{\mu ,s_0}^{\text{even}}$$

$${\hat{\mathcal{H}}}_{\text{TI}}=\frac{k^2}{2m}+v(k_x{\sigma }_y-k_y{\sigma }_x).$$

$${\epsilon }^{\pm }(\mathbf{k})=\frac{k^2}{2m}\pm vk$$

$$\frac{\partial \hat{\mathcal{H}}}{\partial k_x}=\frac{k_x}{m}+v{\sigma }_y,\frac{\partial \hat{\mathcal{H}}}{\partial k_y}=\frac{k_y}{m}-v{\sigma }_x.$$

$${\hat{ȷ}}_{x{s}_x}$$

$${\hat{ȷ}}_{y{s}_x}$$

$$⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{x{s}_x}\right|\pm ,\mathbf{k}⟩$$

$$⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{x{s}_y}\right|\pm ,\mathbf{k}⟩$$

$$⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{y{s}_x}\right|\pm ,\mathbf{k}⟩$$

$$⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{y{s}_y}\right|\pm ,\mathbf{k}⟩$$

$$j_{x,s_y}^{(2)}$$

$$\displaystyle\vec{j}^{(2)}_{\theta s}\equiv\left(\begin{array}[]{l}j^{(2)}_{\theta s_{x}}\\ j^{(2)}_{\theta s_{y}}\\ \end{array}\right)\equiv\left(\begin{array}[]{l}j^{(2)}_{xs_{x}}\cos\theta+j^{(2)}_{ys_{x}}\sin\theta\\ j^{(2)}_{xs_{y}}\cos\theta+j^{(2)}_{ys_{y}}\sin\theta\\ \end{array}\right).$$

$${\hat{\mathcal{H}}}_{\text{RD}}=\frac{k^2}{2m}+\alpha (k_x{\sigma }_y-k_y{\sigma }_x)+\beta \left(k_x{\sigma }_x-k_y{\sigma }_y\right),$$

$${\epsilon }^{\pm }\left(\mathbf{k}\right)=\frac{k^2}{2m}\pm k\sqrt{{\alpha }^2+{\beta }^2-2\alpha \beta \sin 2\varphi }.$$

$$\frac{\partial \hat{\mathcal{H}}}{\partial k_x}=\frac{k_x}{m}+\alpha {\sigma }_y+\beta {\sigma }_x,\frac{\partial \hat{\mathcal{H}}}{\partial k_y}=\frac{k_y}{m}-\alpha {\sigma }_x-\beta {\sigma }_y.$$

$${\hat{j}}_{x{s}_x}$$

$${\hat{j}}_{y{s}_x}$$

$$⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{x{s}_x}\right|\pm ,\mathbf{k}⟩$$

$$⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{x{s}_y}\right|\pm ,\mathbf{k}⟩$$

$$⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{y{s}_x}\right|\pm ,\mathbf{k}⟩$$

$$⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{y{s}_y}\right|\pm ,\mathbf{k}⟩$$

$$j_{x,s_y}^{(2)R}$$

$$=$$

$${{⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{x{s}_x}^D\right|\pm ,\mathbf{k}⟩}_D|}_{\beta \to \alpha }$$

$${{⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{x{s}_y}^D\right|\pm ,\mathbf{k}⟩}_D|}_{\beta \to \alpha }$$

$${{⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{y{s}_x}^D\right|\pm ,\mathbf{k}⟩}_D|}_{\beta \to \alpha }$$

$${{⟨\pm ,\mathbf{k}\left|{\hat{ȷ}}_{y{s}_y}^D\right|\pm ,\mathbf{k}⟩}_D|}_{\beta \to \alpha }$$

$$\frac{j_s^{(2)}\times 2/\hslash }{j_c^{(1)}/e}=\frac{3\gamma e\tau E}{\hslash }$$

$$j_{x,s_x}^{(2)}$$

$$j_{y,s_x}^{(2)}$$

Full text

I Introduction

Spins and their flow in solids have attracted recent intensive attentions from the viewpoints of both fundamental physics and spintronics applications. The conventional and direct way to generate spins or spin current in solids is to inject the spin polarized current from metallic ferromagnets Datta1990a ; Gardelis1999 ; Schmidt1999 ; Hu2001 . Meanwhile, recent researches have been focusing on the electric manipulation of spin and spin current without using the magnets, where the relativistic spin-orbit interaction (SOI) plays an essential role. For such an example, the spin Hall effect supports the conversion of the charge current to the spin current Hirsch1999 ; Zhang2000 ; Murakami2003 ; Kato2004 ; Murakami2004 ; Murakami2004a ; Sinova2004 ; Wunderlich2005 ; Engel2005 ; Sugimoto2006 ; Valenzuela2006 ; Sinova2015 . In the presence of the SOI, the spin Hall conductivity $\sigma_{H}^{s}$ becomes nonzero due to the extrinsic mechanism such as the skew scattering Hirsch1999 ; Zhang2000 ; Kato2004 or the intrinsic mechanism by the Berry phase of the Bloch wave functions Murakami2003 ; Sinova2004 ; Murakami2004 ; Murakami2004a ; Wunderlich2005 . These two mechanisms induce the $\sigma_{H}^{s}$ proportional to $O(\tau)$ and $O(1)$, respectively, in terms of the transport lifetime $\tau$. Spin battery is another method to produce the spin current, where the precession of the ferromagnetic moment is excited by the magnetic resonance absorption, and the damping of this collective mode results in the flow of the spin current to the neighboring system through the interface Saitoh2006a ; Ando2011 ; Dushenko2016 ; Lesne2016 ; Kondou2016 . Interband spin selective optical transition under the irradiation of the circularly polarized light also induces the spin polarized current which is known as the circular photogalvanic effect Ganichev2001 ; Ganichev2002a . These methods have been successfully applied to study the variety of phenomena, but the experimental signals associated with the spin current are quite small and the device structure to detect them is limited. A more efficient way to create the spin current based on another physical origin has been desired for the purpose of spintronics application.

In this paper, we theoretically propose that a simple application of the electric field produces the nonlinear spin current proportional to the square of the electric field ($E^{2}$) and also the square of the transport lifetime ($\tau^{2}$), due to an interplay of the SOI and broken inversion symmetry. Therefore, it can produce larger spin current compared with previous methods. This effect is supported by nontrivial spin texture in energy bands that appears in inversion broken systems with the SOI, e.g., the surface Weyl state of three-dimensional (3D) topological insulators (TI) and two-dimensional (2D) semiconductors with the Rashba and/or Dresselhaus SOI. This new mechanism also offers the rectification of the spin current, i.e., the generation of the DC spin current from AC electric fields. These proposed mechanisms are based on the nonlinear current responses in noncentrosymmetric systems which is captured in the semiclassical treatment using Boltzmann equation as follows.

Noncentrosymmetric systems support nonlinear charge current proportional to $E^{2}$. The canonical example is a p-n junction, where the difference of $I-V$ characteristics between the right and left directions leads to the charge current proportional to $E^{2}$. However, for the periodic systems with conserved crystal momentum $\bm{k}$, the situation is less trivial. This is because the time-reversal symmetry $\mathcal{T}$ imposes the condition on the energy dispersion, i.e., $\varepsilon_{\alpha}(\bm{k})=\varepsilon_{\bar{\alpha}}(-\bm{k})$ with $\bar{\alpha}$ being the opposite spin to $\alpha$. Therefore, even with the broken inversion symmetry $\mathcal{I}$, there remains a certain symmetry between $\bm{k}$ and $-\bm{k}$ as long as one is concerned about the charge degrees of freedom. Thus, in Boltzmann transport phenomena where the charge current is determined by the energy dispersion only, it is necessary to further break the time reversal symmetry in addition to $\mathcal{I}$, e.g., by the external magnetic field $\bm{B}$ or the spontaneous magnetization $\bm{M}$, in order to to realize the nonreciprocal charge responses Rikken2001 ; Krstic2002 ; Rikken2005 ; Pop2014 ; Morimoto2016a ; Yasuda2016 . Exceptions necessarily require that the information of the wave functions enters into the transport properties through e.g. the Berry phase Sodemann2015 ; Morimoto2016b . However, it should be noted that these Berry phase contributions are not the leading order effect in semiclassics. Namely, the dominant one, which is proportional to $(\tau E)^{2}$ in the clean limit, is the contribution captured by the Boltzmann equation.

On the other hand, the situation is dramatically different for the spin current. In this case, one needs to distinguish the spin components of the energy bands. The spin split bands in noncentrosymmetric systems with the SOI could produce the spin current proportional to $(\tau E)^{2}$ even without breaking the $\mathcal{T}$ symmetry. The difference of the required symmetry for the charge current and the spin current is discussed in detail in the section III. Since this effect arises from the Boltzmann transport, the generated nonlinear spin current becomes very large (with $\propto\tau^{2}$) compared with previous methods mentioned above.

We note that the nonlinear spin current in transition metal dichalcogenides (TMDs) was also studied theoretically Yu2014 . While ref. Yu2014 is focused on the band structure with the Ising-type spin splitting along the fixed ($z$-) direction, our theory is applicable to cases with general SOIs that lacks the $S_{z}$ conservation. Especially, Rashba system, being intensively studied in the context of the spintronics, is a typical example that breaks $S_{z}$ conservation. Considering the ubiquitousness of the Rashba system which emerges universally at interfaces and even in the bulkIshizaka2011 ; Sakano2013 , the applicability to such system is a great advantage of the present study for future spintronics studies. Furthermore, the nonlinear spin current in the present study is $2$ or $3$ orders of magnitude larger compared with ref. Yu2014 since the latter is proportional to a small higher order coefficient, namely, the trigonal warping. The detailed comparison to ref. Yu2014 is discussed in the section VI. The present nonlinear spin current also ensures controllability of the spin polarization of the flowing spin current through the direction of the electric field and/or the Rashba-Dresselhaus ratio.

II Theoretical methods

II.1 Boltzmann equation

First we derive the general formula for nonlinear spin current in the semiclassical regime by using Boltzmann equation. We consider a system with the electric field $E$ applied in the $x$ direction. The Boltzmann equation for the distribution function $f$ is given by

[TABLE]

in the relaxation time approximation ($\tau$ being the relaxation time of electron), where $f_{0}$ is the original distribution function in the absence of $E$. (We have set $\hbar=1$ and adopt the convention $e>0$ throughout this paper.) In order to study the (nonlinear) current response in each order in $E$, we expand the distribution function as $f=f_{0}+f_{1}+f_{2}+...$, where $f_{n}\propto E^{n}$ . The iterative substitution in the Boltzmann equation yields $f_{n}=\left(e\tau E\frac{\partial}{\partial{k_{x}}}\right)^{n}f_{0}$ Sodemann2015 ; Yasuda2016 ; Morimoto2016b ; Yu2014 . In particular, the distribution function of the first order in $E$ is given by

[TABLE]

and that of the second order in $E$ is Sodemann2015 ; Yu2014

[TABLE]

The second order term $f_{2}$ typically shows modulation of electron occupation having the quadrupole structure as illustrated in Fig. 1(a).

II.2 definition of spin current

The conventional definition of the spin current operator is given by the anticommutator of the velocity ($\propto\frac{\partial\mathcal{\hat{H}}}{\partial k_{\mu}}$) and the spin ($\propto\sigma_{\nu})$, $\hat{\jmath}_{\mu s_{\nu}}\equiv\frac{1}{4}\left\{\frac{\partial\mathcal{\hat{H}}}{\partial k_{\mu}},\sigma_{\nu}\right\}$ Murakami2004 ; Murakami2004a ; Sinova2004 . Hence, the spin current of the $n$th order in $E$ is given by

[TABLE]

where $\mu$ is the direction of flow, $\nu$ is the direction of the spin polarization, $I$ is the band index, $f_{n}^{I}$ is the $n$th order distribution function for $I$th band. In the following, we focus on the second order nonlinear spin current $j_{\mu,s_{\nu}}^{(2)}$ that appears in noncentrosymmetric systems. Intuitively, an interplay of quadrupole modulation of $f_{2}$ and nontrivial spin texture due to the SOI [as illustrated in Fig. 1(a)] leads to the nonlinear spin current $j_{\mu,s_{\nu}}^{(2)}$ such as shown in Fig. 1(b) as we will see in detail in the section IV.

III symmetry argument

The nonlinear charge and spin current ($j_{\mu,s_{0}}$ and $j_{\mu,s_{\nu}}$, respectively) are constrained by the time reversal symmetry $\mathcal{T}$. To see this, we suppose that the Hamiltonian satisfies $\mathcal{H}(\vec{k},\vec{\sigma})=\mathcal{H}(-\vec{k},-\vec{\sigma})$, and hence, every eigenstate has its time-reversal symmetry partner that carries the opposite momentum and opposite spin. First, the charge current $v_{\mu}=\frac{\partial\mathcal{\hat{H}}}{\partial k_{\mu}}$ is odd under $\mathcal{T}$ ($\mathcal{T}:v_{\mu}\rightarrow-v_{\mu}$) while the spin current is even ($\mathcal{T}:\hat{\jmath}_{\mu s_{\nu}}\rightarrow\hat{\jmath}_{\mu s_{\nu}}$). Next the distribution functions $f_{n}$ is even for even $n$ and odd for odd $n$, because $f_{n}=\left(e\tau E\right)^{n}\frac{\partial^{n}f_{0}}{\partial k_{\mu}^{n}}\sim(v_{\mu})^{n}$. Therefore, it follows that all odd orders of the spin current are zero and that all even orders of the charge-current are zero in the presence of the time-reversal symmetry:

[TABLE]

In particular, we find that the second order charge current vanishes while the second order spin current can be nonvanishing. Finally, a similar argument applies when a system has the inversion symmetry $\mathcal{I}$ with $\mathcal{H}(\vec{k},\vec{\sigma})=\mathcal{H}(-\vec{k},\vec{\sigma})$. Since the spin direction is not flipped by the inversion operator (and hence, $\mathcal{I}:\hat{\jmath}_{\mu s_{\nu}}\rightarrow-\hat{\jmath}_{\mu s_{\nu}}$), all charge and spin nonlinear current in the even order are zero:

[TABLE]

These symmetry analyses indicate that the nonlinear spin current $\propto E^{2}$ in the Boltzmann transport requires broken inversion, but it does not require broken time-reversal symmetric systems. In the following sections, we study a few examples of noncentrosymmetric systems with the SOI that support the nonlinear spin current.

IV Surface state of the 3D TI

We start with the surface of a 3D TI. It is described by the Hamiltonian $\hat{\mathcal{H}}_{\text{TI}}=v(k_{x}\sigma_{y}-k_{y}\sigma_{x}),$ where $v$ is the velocity of the Weyl cone. The energy dispersion is $\varepsilon^{I}=Ivk$ with $I=\pm$, and the spin polarization for each branch in the $k$ space is $\left\langle\pm,\bm{k}\left|\vec{\sigma}\right|\pm,\bm{k}\right\rangle=\pm(-\sin\phi,\cos\phi,0)$, where $k_{x}=k\cos\phi,k_{y}=k\sin\phi$. We show the Fermi surface (FS) and the spin direction for the upper branch together with the second order distribution function in Fig. 1(a). By using spin current operators, $\hat{\jmath}_{xs_{x}}=\hat{\jmath}_{ys_{y}}=0$ and $\hat{\jmath}_{xs_{y}}=-\hat{\jmath}_{ys_{x}}=\frac{1}{2}v,$ we can show that $j_{\mu,s_{\nu}}^{(n)}\propto\int\frac{\mathrm{d}^{2}\bm{k}}{(2\pi)^{2}}\frac{\partial^{n}f_{0}^{\pm}}{\partial k_{x}^{n}}=0.$ Namely, all the spin currents are zero.

However, nonzero spin currents are generated in the presence of the parabolic term $k^{2}/(2m)$ in the Hamiltonian $\hat{\mathcal{H}}_{\text{TI}}$;

[TABLE]

The emergence of the parabolic dispersion is expected in general when the system has a band asymmetry between the electron and hole bands.

The energy dispersion is given by

[TABLE]

and the Fermi surface is formed by one of these two branches depending on the sign of the chemical potential $\mu$. The Fermi momentum is determined as $k_{F}^{\pm}=\mp mv\pm\sqrt{2m\mu+m^{2}v^{2}}$ with $\pm$ corresponding to the sign of $\mu$. The velocity operators in this case are given as

[TABLE]

The spin current operators are given by

[TABLE]

which are summarized as $\hat{\jmath}_{\mu s_{\nu}}=\frac{k_{\mu}}{2m}\sigma_{\nu}$ up to irrelevant constant terms. And their expectation values for each branch of Weyl cone are

[TABLE]

As expected from the symmetry argument, all the linear spin currents vanish after the $\phi$ integration; $j_{x,s_{x}}^{(1)}=j_{x,s_{y}}^{(1)}=j_{y,s_{x}}^{(1)}=j_{y,s_{y}}^{(1)}=0.$ This result can be shown explicitly as follows. All the expectation values of the spin currents are the zeroth or the second order in $\cos\phi$ or $\sin\phi$ while $f_{1}\propto\cos\phi$. The product of these two terms are first or third order in $\cos\phi$ or $\sin\phi$ which vanishes by the $\phi$ integration.

Second order spin currents are calculated by the integration by part at zero temperature as

[TABLE]

By similar calculation shown in Appendix, we have $j_{y,s_{x}}^{(2)}=\frac{1}{5}j_{x,s_{y}}^{(2)}$ and $j_{x,s_{x}}^{(2)}=j_{y,s_{y}}^{(2)}=0$. Note that the signs of the spin currents depend on the sign of chemical potential $\mu$.

Nonzero spin current generation is naturally understood in terms of the spin direction at the FS and the second-order distribution function possessing a quadrupole structure: See Fig. 1(a). Namely, the distribution function is positive toward $\pm x$ direction and hence both the $+y$ spin flowing in the $+x$ direction and the $-y$ spin flowing in the $-x$ direction are accelerated by the application of $E$ parallel to the $x$ direction. In total, $j_{x,s_{y}}^{(2)}$ becomes positive. Similarly, since $f_{2}$ is negative toward $\pm y$ direction, both the $-x$ spin flowing in the $+y$ direction and the $+x$ spin flowing in the $-y$ direction are negatively accelerated, thus resulting in the positive $j_{y,s_{x}}^{(2)}$.

In order to clarify the real space texture of the generated spin current, we define the spin current toward the $\theta$ direction as

[TABLE]

We show the polar plot of $\vec{j}^{(2)}_{\theta s}$ in Fig. 1(b), where the blue line shows the amplitude of the spin current $|\vec{j}^{(2)}_{\theta s}|$ while the black arrows show the direction of the spin polarization $\vec{j}^{(2)}_{\theta s}/|\vec{j}^{(2)}_{\theta s}|$. Using the fact $j^{(2)}_{xs_{y}}=5j^{(2)}_{ys_{x}}$ and $j^{(2)}_{xs_{x}}=j^{(2)}_{ys_{y}}=0$, the magnitude of the spin current is given by $|\vec{j}^{(2)}_{\theta s}|\equiv\sqrt{(j^{(2)}_{\theta s_{x}})^{2}+(j^{(2)}_{\theta s_{y}})^{2}}=|j^{(2)}_{ys_{x}}|\sqrt{25\cos^{2}\theta+\sin^{2}\theta}$, which well describes the blue curve in Fig. 1(b).

V Rashba-Dresselhaus system

Rashba and Dresselhaus type SOIs are present in wide classes of materials without inversion symmetry. The Hamiltonian including the Rashba type and the linear Dresselhaus type SOIs is

[TABLE]

where $m$ is the electron effective mass, $\alpha$ is the Rashba SOI strength and $\beta$ is the Dresselhaus SOI strength. There are two bands indexed by $I=\pm$,

[TABLE]

The spin polarization in the $k$ space is $\left\langle\pm,\bm{k}\left|\vec{\sigma}\right|\pm,\bm{k}\right\rangle=\pm(\cos\varphi,-\sin\varphi,0)$ , where $\varphi\equiv\mathrm{arg}\left[(\beta k_{x}-\alpha k_{y})+i(\beta k_{y}-\alpha k_{x})\right]$. We show FSs and the spin textures for various values of $\tan^{-1}(\beta/\alpha)$ in Fig. 2 while keeping $\alpha^{2}+\beta^{2}=1$. FSs are anisotropic for the general Rashba-Dresselhaus system. In this case there are two FSs in contrast to the case of the surface state of TI.

The anisotropic Fermi momentum for the upper band is $k_{F+}^{+}=-mA+\sqrt{m^{2}A^{2}+2m\mu},$ while those for lower bands are $k_{F-}^{\pm}=+mA\pm\sqrt{m^{2}A^{2}+2m\mu}$, with $A(\phi)=\sqrt{\alpha^{2}+\beta^{2}-2\alpha\beta\sin 2\phi}$. For $\mu>0$, $k_{F+}^{+}$ and $k_{F-}^{+}$ form Fermi surfaces, while $k_{F-}^{+}$ and $k_{F-}^{-}$ do for $\mu<0$. Note that the Fermi surface for $\mu<0$ vanishes for $\phi$ such that $m^{2}A(\phi)^{2}+2m\mu<0$.

Velocity operators are given as

[TABLE]

From these, we have spin current operators as

[TABLE]

which are again summarized as $\hat{\jmath}_{\mu s_{\nu}}=\frac{k_{\mu}}{2m}\sigma_{\nu}+const.$. And their expectation values for each band are

[TABLE]

Using these results, we numerically calculated the second-order spin current for some values of $\tan^{-1}(\beta/\alpha)$ and the direction of the applied electric field $\theta_{E}$, where $\bm{E}=E(\cos\theta_{E},\sin\theta_{E})$. The polar plot of the spin current is summarized in Fig. 3. Note that the distribution function under the application of the electric field in general direction is obtained by a simple substitution $E\frac{\partial}{\partial k_{x}}\rightarrow\bm{E}\cdot\frac{\partial}{\partial\bm{k}}$ in eqs. (2) and (3). We have numerically confirmed that all the first-order spin currents are zero, which is consistent with the symmetry requirement. We have also confirmed the chemical potential dependence is negligible when $\mu>0$. Detailed arguments for the Rashba system ($\tan^{-1}(\beta/\alpha)=0$, the leftmost column in Fig. 3) and the Dresselhaus system ($\tan^{-1}(\beta/\alpha)=\pi/2$ , the rightmost column in Fig. 3) are given below.

V.1 Rashba system

We first investigate the pure Rashba system, for which $\alpha\neq 0$, $\beta=0$ and $\theta_{E}=0$. The eigenstates and the spin polarization are the same as those in the surface state of 3D TI; $\left\langle\pm,\bm{k}\left|\vec{\sigma}\right|\pm,\bm{k}\right\rangle=\pm(\cos\varphi,-\sin\varphi,0)=\pm(-\sin\phi,\cos\phi,0).$ We show the spin textures of the FSs in the pure Rashba system in Fig. 2(a). The spin texture forms vortex structures, whose directions are opposite between the inner and outer FSs. All the first-order spin currents are analytically shown to vanish by the $\phi$ integration; $j_{x,s_{x}}^{(1)R}=j_{x,s_{y}}^{(1)R}=j_{y,s_{x}}^{(1)R}=j_{y,s_{y}}^{(1)R}=0$ where the superscript $R$ indicates Rashba system. This is consistent with the symmetry argument. Furthermore, second-order spin currents are calculated as

[TABLE]

Similarly, we have $j_{y,s_{x}}^{(2)R}=\frac{1}{5}j_{x,s_{y}}^{(2)R}$ and $j_{x,s_{x}}^{(2)R}=j_{y,s_{y}}^{(2)R}=0$ This relation is the same as that in the case of the TI. We note that the sign of the spin currents is opposite compared to that in the TI.

The polar plot of the second-order spin current in the Rashba system is shown in Fig. 3 in the panel corresponding to $\tan^{-1}(\beta/\alpha)=0$ and $\theta_{E}=0$. The shape of the pattern is completely the same as that in 3D TI (Fig. 1(b)) but the spin polarization is opposite.

V.2 Dresselhaus system

We next investigate the Dresselhaus system, for which $\beta\neq 0$, $\alpha=0$ and $\theta_{E}=0$. We show the spin direction of the FSs in the pure Dresselhaus system in Fig. 2(d). The spin texture forms hedgehog structures, whose directions are opposite between the inner and outer FSs; $\left\langle\pm,\bm{k}\left|\vec{\sigma}\right|\pm,\bm{k}\right\rangle=\pm(\cos\varphi,-\sin\varphi,0)=\pm(\cos\phi,-\sin\phi,0)$. The eigenenergy and the distribution functions between the Rashba Hamiltonian and the Dresselhaus Hamiltonian are the same. Only the difference is the expectation value of the spin current operators.

We find the relation between expectation values of spin current operators of the Rashba and the Dresselhaus systems as

[TABLE]

where super and subscripts $R/D$ denote the Rashba/Dresselhaus systems. These relations and the equivalence of the band dispersion guarantee all the linear spin currents to be zero as expected. Furthermore, the second-order spin currents are given by $j_{x,s_{x}}^{(2)\text{D}}=5j_{y,s_{y}}^{(2)\text{D}}=j_{x,s_{y}}^{(2)\text{R}}=5j_{y,s_{x}}^{(2)\text{R}},$ and $j_{x,s_{y}}^{(2)\text{D}}=j_{y,s_{x}}^{(2)\text{D}}=0$ where superscripts $D$ is for Dresselhaus system. The polar plot of the second-order spin current in the Dresselhaus system is shown in Fig. 3. See the panel corresponding to $\tan^{-1}(\beta/\alpha)=\pi/2$ and $\theta_{E}=0$ therein. The peanuts-like shape is completely the same as those in 3D TI and the Rashba system, but the spin polarization reflects the hedgehog structure at FSs.

V.3 Carrier density and temperature dependences

Now we consider the dependence of the spin current on the carrier density $n$ and temperature $T$. We show the carrier density and temperature dependence of $j_{x,s_{y}}^{(2)\text{R}}(=j_{x,s_{x}}^{(2)\text{D}}=5j_{y,s_{x}}^{(2)\text{R}}=5j_{y,s_{y}}^{(2)\text{D}})$ in Fig. 4. We take the Rashba system for example here, but the generic features are common for other cases also. Equation (30) and Fig. 4(a) indicate that the magnitude of the spin current at the zero temperature increases as the increase of carrier density $n$ and becomes constant for $n>n_{D}$ with $n_{D}=m^{2}\alpha^{2}/\pi$ being the carrier density corresponding to the Dirac point. According to eq.(30), the magnitude of the spin current is proportional to the difference of the Fermi momentum defined for each FS, which is constant above the Dirac point. The constant spin current above the Dirac point indicates that the effect of finite temperature is tiny as shown in Fig. 4(b).

VI Discussion

We have demonstrated that the spin current of the second order in $E$ is generated in noncentrosymmeric systems with nontrivial spin textures in the momentum space. We also note that the amplitude of the spin current is $2$ or $3$ orders of magnitude larger than the previous proposal on TMDsYu2014 , indicating that our mechanism can generate nonlinear spin current more efficiently. In TMDs, the anisotropic Fermi surface due to the trigonal warping plays the crucial role in the spin current generation. Authors of ref. Yu2014 claimed that the generated nonlinear spin current normalized by the linear charge current is

[TABLE]

where $\gamma$ is the coefficient of the trigonal warping which has the dimension of the length. (See eqs. (3) and (5) in ref. Yu2014 . We have replaced the coefficient $\beta$ in ref. Yu2014 by $\gamma$ to avoid the confusion. ) The $\gamma$ values are summarized in the Table 1 in ref. Yu2014 , which is of the order of $0.1\sim 1\mathrm{\AA}$ for MoS2 and GaSe. To show that our proposed method is more efficient mechanism to generate the spin current, we calculated the same ratio for our system and define $\gamma$ parameter by eq.(35). For Rashba system with $\mu>0$ for example, the linear charge current $j_{c}^{(1)}=\frac{e^{2}\tau E}{2\pi}(2\mu+m\alpha^{2})$ and the second order spin current is given in eq.(30). The $\gamma$ value is calculated as $\gamma=\frac{5\alpha}{12}\frac{1}{2\mu+m\alpha^{2}}<\frac{5}{12m\alpha}.$ The maximum value of this ratio is achieved by setting $\mu=0$. In this case the $\gamma$ value is $630\mathrm{\AA}$ for GaAs by substituting $\alpha\simeq 0.1$ eVÅSimmons2015 and $m\simeq 0.3m_{e}$ Glover1973 . For the bulk Rashba semiconductor BiTeI, $\gamma$ parameter is $8.1\mathrm{\AA}$ by assuming $\alpha\simeq 3.9\times$ eVÅIshizaka2011 and $m\simeq 0.15m_{e}$Sakano2013 . Here, $m_{e}$ is the electron mass in the vacuum. For the surface of the TI, the linear charge current is $j_{c}^{(1)}=\frac{e^{2}\tau E}{4\pi m}\left[-mv+\sqrt{2m\mu+m^{2}v^{2}}\right]\sqrt{2m\mu+m^{2}v^{2}}$ and second order spin current is given in eq.(14). The $\gamma$ parameter is $\gamma=\frac{5}{12}\frac{1}{\sqrt{2m\mu+m^{2}v^{2}}}\sim\frac{5}{12mv}$ by assuming $2m\mu\ll m^{2}v^{2}$. This value is about $17\mathrm{\AA}$ for 3D TI $\mathrm{Bi_{2}Se_{3}}$ by using $v\simeq 6.2\times 10^{5}$ m/s Zhang2009 and $m=0.53m_{e}$ Kim2016 . These three values are much larger than that discussed in ref Yu2014 . Thus, we can conclude that our proposed method is more efficient mechanism to generate the nonlinear spin current.

Generation of the spin current proportional to $E^{2}$ indicates that the DC spin current is induced by the AC electric field $E(t)=Ee^{i\omega t}$. The time-dependent Boltzmann equation yields the second order distribution function which is composed of two terms; the time-independent term and the one with $2\omega$ frequency Sodemann2015 . The latter one vanishes in time-average, while the former gives us a finite rectified spin current, which is calculable by the equivalent procedure as in the present study. This rectified spin current can be induced for example by shining the terahertz light.

Under the irradiation of the light on systems with spin-splitted bands, the circular photogalvanic effect also contributes to the spin current associated with the charge current. The interband transition with optical selection rule gives us an unbalanced distribution of the positive and negative momenta on the spin splitted band resulting in the spin polarized photocurrent Ganichev2001 ; Ganichev2002a . Similarly, photocurrent is also generated by spin galvanic effect. The optical spin accumulation by the absorption of the circularly polarized light results in the photocurrent induction in the assymetric spin flip scattering processes Ganichev2002b ; Ganichev2004 . However, these phenomena can be excluded by using the linearly polarized terahertz light which does not selectively excite electrons with lifted spin degeneracy.

We next estimate the magnitude of the spin current for various systems. We define the 3D spin conductivity as $j_{x,s_{y}}^{(2)}/(E\hbar/(2ec))$, where $c$ is the lattice constant for thickness direction and we assume the reasonable value of the magnitude of electric field $E\simeq 10^{2\sim 5}$ V/m. The spin conductivity is the order of $2\times 10^{2\sim 5}\ \mathrm{\Omega^{-1}m^{-1}}$ for GaAs by substituting $\tau\simeq 2.5$ ps and $c=5.7\mathrm{\AA}$. It is also the order of $7\times 10^{0\sim 3}\ \mathrm{\Omega^{-1}m^{-1}}$ for BiTeI with $\tau\simeq 0.072$ ps Wang2013a and $c=6.9\mathrm{\AA}$. As for 3D TI, it is the order of $1.3\times 10^{2\sim 5}\ \mathrm{\Omega^{-1}m^{-1}}$ for Bi2Se3 by substituting $\tau\simeq 3.1$ ps Glinka2013 , $c=29\mathrm{\AA}$ and $\mu=0.1$ eV. These values are larger than the typical value of the spin Hall conductivity $10^{0\sim 4}\ \mathrm{\Omega^{-1}m^{-1}}$ Sinova2015 . The effect of finite temperature summarized in Fig. 4 has a peculiar feature. For the typical sheet carrier density of the order of $\sim 10^{13}\mathrm{cm^{-2}}$, the carrier density $n/n_{D}$ is the order of $10^{4}$ for GaAs and $1$ for BiTeI. The room temperature in Fig. 4, $300\mathrm{K}/(m\alpha^{2}/k_{B})$, is about $400$ for GaAs and $0.1$ for BiTeI. As seen in Fig. 4, we may conclude that the spin current never reduces even at room temperature.

Finally, we discuss the validity of the present work. Our derivation of the second order distribution function is based on the expansion with respect to $\tau eE\frac{\partial}{\partial\hbar k_{x}}\sim\tau eE/(\hbar k_{0})$, where $k_{0}\sim m\alpha/\hbar^{2}$ is the typical momentum of, for instance, the Rashba system. For the convergence of the expansion, the electric field $E$ must satisfy $E\ll m\alpha/(e\hbar\tau)$. This condition has two physical interpretations. One is that the energy due to the electric field $eE/k_{0}$ must be much smaller than the disorder broadening $\hbar/\tau$, and the other is the distance between two FSs $m\alpha/\hbar^{2}$ must be much larger than the shift of the distribution function in the momentum space $\tau eE/\hbar$ to avoid the level mixing by the applied electric field. The upper limit of the electric field is of the order of $10^{3}$ V/m for GaAs, $10^{6\sim 7}$ V/m for BiTeI and $10^{5\sim 6}$ V/m for Bi2Se3; the latter two values are sufficiently large for usual terahertz experiments ($\sim 10^{5}$ V/m). Note that two FSs come very close when $\alpha~{}\sim\beta$. In this case, the distance between two FSs in the momentum space becomes $m|\alpha-\beta|/\hbar^{2}\sim 0$, and hence our results are not valid near the persistent helix phase, $\alpha=\beta$.

Acknowledgment — We thank M. Kawasaki and Y. Tokura for fruitful discussions. This work was supported by the Grants-in-Aid for Scientic Research from MEXT KAKENHI (Grant Nos.JP25400317, JP15H05854 and JP17K05490) (ME), the Gordon and Betty Moore Foundation’s EPiQS Initiative Theory Center Grant (TM), and JSPS Grant-in-Aid for Scientic Research (No. 24224009, and No. 26103006) from MEXT, Japan, and ImPACT Program of Council for Science, Technology and Innovation (Cabinet office, Government of Japan) (NN). KWK acknowledges support from “Overseas Research Program for Young Scientists” through Korea Institute for Advanced Study (KIAS). This work is also supported by CREST, JST (JPMJCR16F1).

Appendix A Derivation of second order spin current

In this section, we show the derivation of the spin currents; $j_{x,s_{x}}^{(2)},j_{y,s_{x}}^{(2)},j_{y,s_{y}}^{(2)}$ for the surface of 3D TI and the Rashba system which are skipped in the main text.

A.1 Surface of 3D TI

For surface states of 3D TIs, the spin currents are given by

[TABLE]

[TABLE]

and

[TABLE]

A.2 Rashba system

For Rashba systems, the spin currents are given by

[TABLE]

[TABLE]

and

[TABLE]

A.3 $1/m$ expansion in the surface of TI

In this subsection, we investigate the effect of the parabola term in the Hamiltonian of the surface of 3D TI. When we expand the second order spin current with respect to $1/m$, we obtain

[TABLE]

The leading term is the expression for $j^{(2)}_{xs_{y}}$, which is obtained only by considering the correction to the current operator due to the $k^{2}$ dispersion. Namely, when we evaluate the expectation value $\left\langle I,\bm{k}\left|\hat{\jmath}_{\mu s_{\nu}}\right|I,\bm{k}\right\rangle$, we use $\hat{\mathcal{H}}=k^{2}/(2m)+v(k_{x}\sigma_{y}-k_{y}\sigma_{x})$ for the Hamiltonian but $\varepsilon^{\pm}=\pm vk$ for the distribution function. In this situation, for example, $j_{x,s_{y}}^{(2)}$ is given by

[TABLE]

This indicates that the spin current at the TI surface arises from the interplay between the surface Weyl state exhibiting a nontrivial spin texture and the effect of the $k^{2}$ dispersion introducing the $k$ linear term in the current operator.

A.4 Numerical calculation in Rashba-Dresselhaus system

The second order spin current in the coexistence of the Rashba and the Dresselhaus terms, for example, $j_{x,s_{x}}^{(2)}$ is calculated as

[TABLE]

The analytical integration over $k$ is possible for given values of $\phi$ with the use of the relations

[TABLE]

[TABLE]

Then, the integral over $\phi$ is evaluated numerically. The similar calculations are carried out for the other components of the second order spin current. This explassion is used in the Fig. 3 in the main text.