Linear and nonlinear spin current response of anisotropic spin-orbit coupled systems

We consider a 2DEG with an arbitrary crystal orientation defined by the unit normal ${\hat{\textbf{n}}} = n_x{\hat{\textbf{x}}}+n_y{\hat{\textbf{y}}}+ n_x{\hat{\textbf{z}}}$, with underlying basis vectors ${\hat{\textbf{x}}}, {\hat{\textbf{y}}}, {\hat{\textbf{z}}},$ pointing along the crystal axes $[100], [010], [001]$, respectively. The system is in the presence of Rashba and linear-in-k Dresselhaus SO couplings. The corresponding Hamiltonian is determined by the function $\varepsilon_0(k_x,k_y) = \hbar^2 k^2/2m$ and the SO vector field $d_p(k_x,k_y) = \mu_{p\nu}k_{\nu}\,(p = x,y,z; \nu = x,y)$. For a given normal ${\hat{\textbf{n}}}$, the xyz coordinate system can be rotated in order to obtain a new $x^{\prime} y^{\prime} z^{\prime}$ coordinate system, with the zʹ-axis pointing along the direction of ${\hat{\textbf{n}}}$, and orthonormal basis $\hat{{\textbf{l}}},\hat{{\textbf{m}}},\hat{\textbf{n}}$, $\hat{{\textbf{l}}} = {\hat{\textbf{m}}}\times{\hat{\textbf{n}}}$. Thus, for each orientation the matrix of SO material parameters µ_ij must be understood as referred to the new coordinates, as well as the condition ${\textbf{k}}\cdot{\hat{\textbf{n}}} = 0$, which is transformed to the condition $k^{\prime}_z = 0$. We shall use the symbol R+D$[hkl]$ to indicate the Rashba (R) and Dresselhaus (D) SO couplings when the sample is grown along the crystallographic direction $[hkl]$.

The k-space available for vertical transitions is determined by the condition $\varepsilon_{-}(\vec{k}) \leqslant \varepsilon_\mathrm{F} \leqslant \varepsilon_{+}(\vec{k})$ ('Pauli blocking'), and the conservation of energy $\varepsilon_+(k_x,k_y)-\varepsilon_-(k_x,k_y) = \hbar\omega$, for a given Fermi energy $\varepsilon_\mathrm{F}$ and exciting frequency ω. This means that only those points $(k_x,k_y)$ lying on the resonance curve $C_\mathrm{r}(\omega) = \{(k_x,k_y)| 2d(k_x,k_y) = \hbar\omega\}$ which satisfy the inequality $k_\mathrm{F}^+(\theta)\lt k\lt k_\mathrm{F}^-(\theta)$ will contribute to the JDOS (see figure 1(a)). Here $k_\mathrm{F}^{\lambda}(\theta)$ are the Fermi contours, defined by the equation $\varepsilon_{\lambda}(k_x,k_y) = \varepsilon_\mathrm{F}$ which, written in polar coordinates, are given by $k^{\lambda}_\mathrm{F}(\theta) = \sqrt{2m \varepsilon_\mathrm{F}/\hbar^2+ k^{2}_\mathrm{so}(\theta)} - \lambda k_\mathrm{so}(\theta)$, where $k_\mathrm{so}(\theta) = mg_{[hkl]}(\theta)/\hbar^2$ is a characteristic wave number of the SO interaction. The function $g_{[hkl]}(\theta)$ accounts for the anisotropy of the energy splitting $2d(k_x,k_y) = 2kg_{[hkl]}(\theta)$, with $g_{[hkl]}(\theta) = |\boldsymbol{\mu}_x\cos\theta+\boldsymbol{\mu}_y\sin\theta|$, where the vectors $\boldsymbol{\mu}_\nu$, with components $(\boldsymbol{\mu}_{\nu})_{p} = \mu_{p\nu}$, are the columns of the matrix $\mu_{i\nu}$. The above restrictions imply that the allowed transitions are possible only in the energy window $\hbar\omega_+(\theta)\leqslant\hbar\omega\leqslant\hbar\omega_-(\theta)$, where $\hbar\omega_{\lambda}(\theta) = 2d(k_\mathrm{F}^{\lambda}(\theta)) = 2k_\mathrm{F}^{\lambda}(\theta)g_{[hkl]}(\theta)$ is the minimum (maximum) photon energy $\hbar\omega_+(\theta)$ ($\hbar\omega_-(\theta)$) required to induce a vertical transition between states lying along the direction θ in the k-space (inset in figure 1(a)). As a consequence, the JDOS for the system with SO interaction R+D$[hkl]$ reads as

$$\begin{equation} J_{+-}(\omega) = \frac{\hbar\omega}{16\pi^{2}} \int \!\! \mathrm{d}\theta\,\,\frac{\Theta[1-|\xi(\omega,\theta)|]}{g^{2}_{[hkl]}(\theta)}\ , \end{equation} \tag{ 7 }$$

where $\Theta(x)$ is the Heaviside unit step function, and $\xi(\omega,\theta) = [\omega-\frac{1}{2}(\omega_-(\theta)+\omega_+(\theta))]/[\frac{1}{2}(\omega_-(\theta)-\omega_+(\theta))]$. According to this expression the edges of the absorption spectrum will be at the photon energies $\hbar\omega = \textrm{min}_{\theta}\{\hbar\omega_+(\theta)\}$ and $\hbar\omega = \textrm{max}_{\theta}\{\hbar\omega_-(\theta)\}$, which will occur along some directions $\theta_\lt$ and $\theta_\gt$ where the energy splitting is minimum and maximum, respectively. In general, the set of van Hove singularities can be identified geometrically by analyzing how the resonance curve $C_\mathrm{r}(\omega)$ enters, intersects, and leaves the region of allowed transitions $k_\mathrm{F}^+(\theta)\lt k\lt k_\mathrm{F}^-(\theta)$ (shaded area in figure 1(a)) as frequency varies. The resonance curve is a rotated ellipse with equation $|k_x\boldsymbol{\mu}_x+k_y\boldsymbol{\mu}_y| = \hbar\omega/2$. After a rotation by an angle ζ defined by $\tan 2\zeta = 2\boldsymbol{\mu}_x\cdot\boldsymbol{\mu}_y/[|\boldsymbol{\mu}_x|^2-|\boldsymbol{\mu}_y|^2]$, the ellipse becomes of the standard form with principal axes given by $Q_x(\omega) = \hbar\omega/2g_{[hkl]}(\theta_\lt)$ and $Q_y(\omega) = \hbar\omega/2g_{[hkl]}(\theta_\gt)$, where $\theta_\gt = \theta_\lt+\pi/2$. In order to maintain the angle $\theta_\lt$ as the direction of global minimum energy separation (as defined above) we choose $\theta_\lt = \zeta$ if $g_{[hkl]}(\zeta) \lt g_{[hkl]}(\zeta + \pi/2)$, and $\theta_\lt = \zeta + \pi/2$ otherwise. The critical van Hove energies are given by those values of $\hbar\omega$ at which the axes $Q_x(\omega)$ and $Q_y(\omega)$ of the resonance ellipse touch tangentially the Fermi contours $k_\mathrm{F}^{\pm}(\theta)$ along $\theta_\lt$ and $\theta_\gt$. Figure 1(a) illustrates the situation for the R+D[123] case. The critical frequencies are then defined by the matchings $Q_x(\omega^+) = k_\mathrm{F}^+(\theta_\lt)$, $Q_x(\omega_a) = k_\mathrm{F}^-(\theta_\lt)$, $Q_y(\omega_b) = k_\mathrm{F}^+(\theta_\gt)$, and $Q_y(\omega^-) = k_\mathrm{F}^-(\theta_\lt)$:

$$\begin{align} \hbar\omega^+ & = 2 k^{+}_\mathrm{F}(\theta_{\lt}) g_{[hkl]}(\theta_{\lt}), \end{align} \tag{ 8 }$$

$$\begin{align} \hbar\omega_a & = 2 k^{-}_\mathrm{F}(\theta_{\lt}) g_{[hkl]}(\theta_{\lt}), \end{align} \tag{ 9 }$$

$$\begin{align} \hbar\omega_b & = 2 k^{+}_\mathrm{F}(\theta_{\gt}) g_{[hkl]}(\theta_{\gt}), \end{align} \tag{ 10 }$$

$$\begin{align} \hbar\omega^- & = 2 k^{-}_\mathrm{F}(\theta_{\gt}) g_{[hkl]}(\theta_{\gt}). \end{align} \tag{ 11 }$$

Note that while $\omega^+$ is always smaller than ω⁻, the order relation between ω_a and ω_b can change for different orientations.

Zoom In Zoom Out Reset image size

In figure 1(b) we show the JDOS for different crystal orientations as a function of frequency, keeping the same SO parameters values. The overall shape and size of the spectra reveals a strong dependence on the direction of sample growth. For the case R+D[111], the Hamiltonian is formally identical to that of a system with Rashba coupling only, which presents an isotropic splitting of the states. Thus, the resonance curve and Fermi contours are concentric circles and the JDOS displays the well known box-like shape with only two spectral features ($\hbar\omega^{\pm}$). However, for other orientations the k-space for allowed optical transitions becomes no longer isotropic, and two types of shapes may appear, depending on the relative values of ω_a and ω_b . When $\omega_a \lt \omega_{b}$ the spectra can develop a convex shape between these critical energies, as shown in figure 1(b) for the R+D[001] and R+[123] systems. On the other hand, when $\omega_a \gt \omega_b$ the JDOS presents a linear dependence instead, as is illustrated by the R+D[110] case. This can be explained by observing how the resonance curve $C_\mathrm{r}(\omega)$ overlaps the region of allowed transitions bounded by the Fermi contours in each case. When $\omega_{a}\lt\omega_{b}$, the semi-axis $Q_{x}(\omega)$ of the curve $C_\mathrm{r}(\omega)$ will reach the line $k^{-}_\mathrm{F}(\theta)$ first before the semi-axis $Q_{y}(\omega)$ intersects the line $k^{+}_\mathrm{F}(\theta)$, which means that for $\omega_{a} \lt \omega \lt \omega_b$ there is a portion of the curve which does not contribute to the JDOS. In contrast, when $\omega_a \gt \omega_b$ we have the opposite situation, the semi-axis $Q_{x}(\omega)$ touch the Fermi contour $k^{-}_\mathrm{F}(\theta)$ after the axis $Q_{y}(\omega)$ contacts the Fermi contour $k^{+}_\mathrm{F}(\theta)$. This imply that there is a range of frequencies, $\omega_b \lt \omega \lt \omega_a$, for which the ellipses $C_\mathrm{r}(\omega)$ lie entirely within the allowed zone (shaded area in figure 1(a)), causing a linear increase in the JDOS. The dependence of the critical frequencies (8)–(11) on the Miller indices hints the growth direction as an additional element to have spectra with specific width, shape or frequency region.

The spin conductivity tensor at the fundamental frequency for the SO coupled 2DEG, as obtained from Kubo formula (3), is (no sum over repeated indices is implied)

$$\begin{align} \operatorname{Re}\sigma^{\ell}_{ij}(\omega) & = \sigma^{\ell}_{ij}(0) - \frac{e}{8\pi} \left( \boldsymbol{\mu_{x}} \times \boldsymbol{\mu_{y}} \right)_{\ell} \frac{\hbar\omega}{8m/\hbar^2} \frac{1}{2\pi} \nonumber\\ & \quad \times \int^{2\pi}_{0} \frac{\mathrm{d}\theta}{g^4(\theta)} \ln \left| \frac{\left[ \omega + \omega_{+}(\theta) \right]\left[ \omega - \omega_{-}(\theta) \right]}{\left[ \omega - \omega_{+}(\theta) \right]\left[ \omega + \omega_{-}(\theta) \right]}\right| \nonumber \\ & \quad \times \left[ \epsilon_{ijz} + \sin2\theta \left( \delta_{iy} - \delta_{ix} \right)\delta_{ij} + \cos2\theta \left( 1-\delta_{ij} \right) \right] \end{align} \tag{ 12 }$$

$$\begin{align} \operatorname{Im}\sigma^{\ell}_{ij}(\omega) & = -\frac{e}{8\pi} \left( \boldsymbol{\mu_{x}} \times \boldsymbol{\mu_{y}} \right)_{\ell} \frac{\hbar\omega}{16m/\hbar^2}\nonumber\\ & \quad \times \int^{2\pi}_{0} \! \frac{\mathrm{d}\theta}{g^4(\theta)}\left[ \epsilon_{ijz} + \sin2\theta \left( \delta_{iy} - \delta_{ix} \right) \delta_{ij} \right.\nonumber\\ & \quad \left. + \cos2\theta \left( 1-\delta_{ij} \right) \right] \Theta [ 1-|\xi(\omega,\theta)|], \end{align} \tag{ 13 }$$

where

$$\begin{align} \sigma^{\ell}_{ij}(0) & = - \frac{e}{8\pi} \left( \boldsymbol{\mu_{x}} \times \boldsymbol{\mu_{y}} \right)_{\ell} \frac{1}{2\pi} \int^{2\pi}_{0} \! \frac{\mathrm{d}\theta}{g^2(\theta)}\nonumber\\ & \quad \times \left[ \varepsilon_{ijz} + \sin2\theta \left( \delta_{iy} - \delta_{ix} \right)\delta_{ij} + \cos2\theta \left( 1-\delta_{ij} \right) \right], \end{align} \tag{ 14 }$$

is the dc spin conductivity. Note that $\sigma^{\ell}_{yy}(\omega) = -\sigma^{\ell}_{xx}(\omega)$. Complex integration gives the result

$$\begin{equation} \sigma^{\ell}_{ij}(0) = -\frac{e}{8\pi}\frac{(\boldsymbol{\mu}_x\times \boldsymbol{\mu}_y)_{\ell}}{|\boldsymbol{\mu}_x\times\boldsymbol{\mu}_y|} \left(\!\!\begin{array}{cc} -\mbox{Im}\,(z_+) & 1+\mbox{Re}\,(z_+) \\ \\ -1+\mbox{Re}\,(z_+) & \mbox{Im}\,(z_+) \end{array}\right) \ , \end{equation} \tag{ 15 }$$

where

$$\begin{equation} z_+ = -\frac{[A-\sqrt{A^2-(B^2+C^2)}]}{B^2+C^2}(B+\mathrm{i}C) \ , \end{equation} \tag{ 16 }$$

with $A = (|\boldsymbol{\mu}_x|^2+|\boldsymbol{\mu}_y|^2)/2$, $B = (|\boldsymbol{\mu}_x|^2-|\boldsymbol{\mu}_y|^2)/2$, and $C = \boldsymbol{\mu}_x\cdot\boldsymbol{\mu}_y$. The general expression (15) extends the result reported in [54], which is valid for SO vector fields with $d_z({\textbf{k}}) = 0$ ($\mu_{zx} = \mu_{zy} = 0$) only.

Since $\left( \boldsymbol{\mu_{x}} \times \boldsymbol{\mu_{y}} \right)_{x}$ and $\left( \boldsymbol{\mu_{x}} \times \boldsymbol{\mu_{y}} \right)_{y}$ vanish for $\mu_{z\nu} = 0$ we have that the only systems that can support a linear spin current with perpendicular-to-plane spin polarization strictly, are those grown along the [001] and [111] directions; any other crystal orientation will have in-plane spin polarized current components. Moreover, the vanishing of the common factor $\boldsymbol{\mu_{x}} \times \boldsymbol{\mu_{y}}$ implies the absence of an induced spin current via electric-dipole interaction in the 2DEG with R+D$[hkl]$ SO coupling. This reminds the well known effects due to the recovery of the SU(2) symmetry predicted in systems with R+D[001], like the infinite spin lifetime due to fixed spin precession axis [20] or the formation of a persistent spin helix state [19, 21].

As we mentioned before, Kammermeier et al [45] found that for an arbitrary crystal orientation is still possible to have conditions for spin-preserving symmetries due to the interplay of Rashba and Dresselhaus SO interaction. The requirement for that is to have samples with two Miller indices equal in modulus and a particular relation between the Rashba and Dresselhaus parameters (in our language, a proper combination of the elements of the SO matrix µ_ij ) [45]. It can be verified that for these special conditions, the factor $\boldsymbol{\mu_{x}} \times \boldsymbol{\mu_{y}}$ is zero. Without loss of generality, lets choose, after Kammermeier, the orientation $\vec{\hat{n}} = (\eta,\eta,n_z)$, with $\vec{\hat{m}} = (-1,1,0)/\sqrt{2}$ and $\vec{\hat{l}} = (n_z, n_z, -2\eta)/\sqrt{2}$, where $n_z^2 = (1-2\eta^2)$. The corresponding vectors of SO parameters become $\boldsymbol{\mu_x} = (0, -\alpha + \gamma k^{2}_{n}(1-9\eta^2)n_z, 0)$ and $\boldsymbol{\mu_y} = (\alpha + \gamma k^{2}_{n}(1+3\eta^2)n_z, 0, -\gamma k^{2}_{n} \sqrt{2}\eta (1-3\eta^2))$, where α and $\gamma k_n^2$ are the Rashba and linear Dresselhaus coupling strengths, respectively; here γ is the bulk Dresselhaus parameter characterizing the cubic Dresselhaus interaction and $k_n^2 = \langle({\textbf{k}}\cdot{\hat{\textbf{n}}})^2\rangle$ is an expectation value appearing in the linearization for symmetric narrow quantum wells [55]. When these satisfy the relation $\alpha/\gamma k_n^2 = (1-9\eta^2)n_z$ ($\mu_{yx} = 0$), the SO vector field becomes collinear, ${\textbf{d}}({\textbf{k}}) = \gamma k_n^2(3\eta^2-1)(-2n_z,0,\sqrt{2}\,\eta)k_y$. Under these conditions, $\boldsymbol{\mu_{x}} \times \boldsymbol{\mu_{y}} = \mu_{yx}(\mu_{zy},0,-\mu_{xy}) = 0$, implying the vanishing of the spin current. Note that when $\mu_{zy} = 0$ (which is true for the [001] and [111] growth directions, with η = 0 and $\eta^2 = 1/3$, respectively) with $\mu_{yx}\neq 0$ and $\mu_{xy}\neq 0$ (that is $\alpha/\gamma k_n^2\neq \pm 1$ for [001] and $\alpha/\gamma k_n^2\neq -2/\sqrt{3}$ for [111]), the polarization of the spin current is along the direction ${\hat{\textbf{n}}}$ (zʹ). On the other hand, if $\mu_{xy} = 0$ with $\mu_{yx}\neq 0$ and $\mu_{zy}\neq 0$ the polarization of the spin current is along ${\hat{\textbf{l}}}$ (x^ʹ) direction. In contrast, the spin current polarized parallel to the ${\hat{\textbf{m}}}$ direction is null regardless the magnitude of the SO strength parameters, $\sigma^y_{ij}(\omega) = 0$. Interestingly, although the condition $\mu_{xy} = 0$, which rewrites as $\alpha/\gamma k^{2}_{n} = -(1+3\eta^2)n_z$, does not corresponds to a collinear SO vector field, still produces an absence of out-of-plane-polarized spin current.

In figure 2, a typical spectrum of the linear spin current conductivity $\sigma^z_{ij}(\omega)$ is displayed, in this case for the R+D[123] system. As anticipated by the JDOS, we can identify the presence of van Hove features at the critical energies (8)–(11), defined by the Pauli blocking and the energy conservation condition for vertical transitions (figure 1(b)). Given that the critical frequencies (8)–(11) depend not only on the magnitude and relative value of the SO material parameters but also on the crystal orientation, a comparative study of the spectra for different orientations would allow to choose the linear spin current response with particular spectral characteristics like width, shape, size, or frequency region.

Zoom In Zoom Out Reset image size

The second-order Kubo formula (6) leads to the SH spin conductivity of the R+D$[hkl]$ system,

$$\begin{align} \begin{array}{l} \sigma _{ijl}^{\ell ,(2\omega )}(\omega ) = \frac{{{e^2}}}{{{{(\hbar \omega )}^2}}}\frac{{({\hbar ^2}/m)}}{{{{(2\pi )}^2}}}({{\boldsymbol{\mu }}_{\textbf{x}}} \times {{\boldsymbol{\mu }}_{\textbf{y}}})\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{ \times _p}\int_0^{2\pi } \frac{{{\textrm{d}}\theta }}{{{g^6}(\theta )}}{{\hat k}_i}{{\hat k}_\nu }\left\{ {{{({{\boldsymbol{\mu }}_{\textbf{x}}} \times {{\boldsymbol{\mu }}_{\textbf{y}}})}_p}{\mu _{\ell \nu }}({\delta _{jl}} - {{\hat k}_j}{{\hat k}_l})C(\omega ,\theta )} \right.\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, - {\mkern 1mu} {_{\ell pq}}{\mu _{q\rho }}{{\hat k}_\rho }({{\boldsymbol{\mu }}_{\textbf{j}}} \cdot {{\boldsymbol{\mu }}_{\boldsymbol{\nu }}})({{\hat k}_x}{\delta _{ly}} - {{\hat k}_y}{\delta _{lx}})[C(\omega ,\theta )\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { - 4C(2\omega ,\theta )] + (j \leftrightarrow l)} \right\} \end{array} \end{align} \tag{ 17 }$$

where $\hat{k}_{i}(\theta) = \cos\theta \delta_{ix} + \sin\theta \delta_{iy}$ and

$$\begin{align} C(x,\theta) & = -\frac{1}{4} \left( \frac{2mg^2(\theta)}{\hbar^2} + \frac{\hbar x}{4} \operatorname{ln} \bigg\vert \frac{[x+\omega_{+}(\theta)][x-\omega_{-}(\theta)]}{[x-\omega_{+}(\theta)][x+\omega_{-}(\theta)]} \bigg\vert \right)\nonumber\\ &\quad -\mathrm{i}\pi\frac{\hbar x}{16}\,\Theta[1-|\xi(x,\theta)|]. \end{align} \tag{ 18 }$$

As in the linear response, the system grown along $\vec{\hat{n}} = (\eta,\eta,n_z)$ will not support a spin current at 2ω when $\mu_{yx} = 0$, or equivalently when the SO field $d_{p}({\textbf{k}}) = \mu_{p\nu}k_{\nu}$ becomes collinear. However, for the same special class of orientations, to analyze the spin conductivity for each spin index $\ell$, we have to focus in the following factors

$$\begin{align} (\boldsymbol{\mu_{x}})_{\ell} & = \mu_{yx}\delta_{\ell y} , \end{align} \tag{ 19 }$$

$$\begin{align} (\boldsymbol{\mu_{y}})_{\ell} & = \left( \mu_{xy}\delta_{\ell x} + \mu_{zy} \delta_{\ell z} \right), \end{align} \tag{ 20 }$$

$$\begin{align} \big[ \left( \boldsymbol{\mu_{x}} \times \boldsymbol{\mu_{y}} \right) \times \boldsymbol{\mu_{x}} \big]_{\ell} & = \mu^{2}_{yx} \big( \mu_{xy}\ \delta_{\ell x} + \mu_{zy}\ \delta_{\ell z} \big) , \end{align} \tag{ 21 }$$

$$\begin{align} \big[ \left( \boldsymbol{\mu_{x}} \times \boldsymbol{\mu_{y}} \right) \times \boldsymbol{\mu_{y}} \big]_{\ell} & = -\mu_{yx} \left( \mu^{2}_{xy} + \mu^{2}_{zy} \right) \delta_{\ell y}. \end{align} \tag{ 22 }$$

We note that there are some possibilities to control the polarization of the nonlinear spin current. For $\mu_{yx}\neq 0$, the condition $\mu_{xy} = 0$ ($\mu_{zy} = 0$) corresponds to have a spin current flowing with the spin lying in the plane defined by ${\hat{\textbf{m}}}$ and ${\hat{\textbf{n}}}$ (${\hat{\textbf{l}}}$ and ${\hat{\textbf{m}}}$), that is $\sigma^{x,(2\omega)}_{ijl}(\omega) = 0$ ($\sigma^{z,(2\omega)}_{ijl}(\omega) = 0$). On the other hand, to have $\sigma^{y,(2\omega)}_{ijl}(\omega) = 0$ is possible only for $\mu_{yx} = 0$. Therefore, we have that for a system with crystal orientation $\vec{\hat{n}} = (\eta,\eta,n_z)$, and given that the linear $\sigma^y_{ij}(\omega) = 0$, the generation of a spin current polarized along the $\vec{\hat{m}}$ (yʹ) direction will depend quadratically on the electric field, because it is induced as a second-order response strictly.

Figure 3 shows the SH spin current conductivity $\sigma^{z,(2\omega)}_{ijl}(\omega)$ for the R+D[123] system. As expected, beside the peaks related with critical points at the energies $\hbar\omega^{\pm},\hbar\omega_a,\hbar\omega_b$, van Hove singularities appear also at their subharmonics. Similarly to the first-order spin current conductivity, the magnitude, direction, and spin polarization of the nonlinear spin current could be modified through frequency variation or by the relative value of the SO strengths. Moreover, as for the linear response, a comparative study between several orientations, would allow to decide in favor of particular spectral characteristics.

Zoom In Zoom Out Reset image size

Another aspect worth noting is that the R+D$[hkl]$ systems will present a nonlinear spin Hall effect, through $\sigma^{z,(2\omega)}_{yxx}(\omega) \neq 0$ and $\sigma^{z,(2\omega)}_{xyy}(\omega) \neq 0$, except for R+D[001] and R+D[111] cases. The necessary condition for this phenomenon is to have $d_z(\vec{k}) \neq 0$ $(\mu_{z\nu} \neq 0)$, which is characteristic of [110] samples, one of the low Miller indices usually studied.

An estimation of the relative size of the linear and SH spin currents (taken from the example shown in figures 2 and 3), out of vanishing conditions, gives ${\cal{J}}^{(2)}/{\cal{J}}^{(1)}\sim (\sigma_0^{(2)}E/\sigma_0^{(1)})N(\omega) = [(eE/k_\mathrm{R})/\varepsilon_\mathrm{R}]N(\omega)$, where $N(\omega)$ is a frequency dependent number, and $k_\mathrm{R} = m\alpha/\hbar^2,\,\varepsilon_\mathrm{R} = \alpha k_\mathrm{R}$. The factor $(eE/k_\mathrm{R})/\varepsilon_\mathrm{R}$ is a measure of the coupling between the SO dipole $e/k_\mathrm{R}$ and the electric field, compared to the characteristic SO energy. At $\hbar\omega\approx 4$–5 meV (within THz regime), $N\sim 10^{-3}$–10⁻², and for an electric field $E\sim 10^3$–10⁵ V m⁻¹ [14, 15], we estimate ${\cal{J}}^{(2)}/{\cal{J}}^{(1)}\sim 10^{-3}$–1. For a comparison to the first-order charge current we estimate ${\cal{J}}^{(2)}(2/\hbar)/(J^{(1)}/e)\sim [(eE/k_\mathrm{R})/\varepsilon_\mathrm{R}]\times 10^{-4} \sim 10^{-2}$ at $\hbar\omega = 6\,$meV for $E\sim 10^5$ V m⁻¹. These estimations suggest that the induce SH spin current can reach not negligible size at reasonable values of the electric field, for some frequencies in the range of few THz.

To end this section we comment on the consequences of having cubic-in-k terms in the spin–orbit coupling. The inclusion of k³-Dresselhaus interaction breaks the SU(2) symmetry achieved under Kammermeier conditions for a linear-in-k Hamiltonian. However, under these same conditions, the total SO field (R+D with linear and cubic terms) for the [111] and [110] orientations is still collinear, as was reported by Kammermeier et al [45]. This leads to vanishing first and second order spin currents, as can be shown from the general expressions (3) and (6) in section 2, by noting that every term in them involves the vector ${\textbf{D}}_i\equiv{\textbf{d}}({\textbf{k}}) \times \partial {\textbf{d}}({\textbf{k}})/\partial k_i$. When the SO field ${\textbf{d}}({\textbf{k}})$ is collinear, ${\textbf{d}}\parallel \partial {\textbf{d}}/\partial k_i$ and then ${\textbf{D}}_i = {\textbf{0}}$. For other growth directions one may guess that such a vanishing will not occur, but it is hard to predict any result without explicit calculations. The same is true for the analysis of the specific polarization of the spin currents, even for the special [111], [110] cases.

In this section we considered another anisotropic model, used to study the effect of an anisotropic Rashba splitting on the longitudinal optical conductivity of 2D puckered structures [51] like black phosphorus and group IV monochalcogenides [56–58]. Here the reduced symmetry of the splitting of the states is due to mass anisotropy. The kinetic contribution to the low energy Hamiltonian is that of an anisotropic free electron gas [49, 59], $\varepsilon_{0}(\vec{k}) = \hbar^2 k^{2}_{x}/2m_x + \hbar^2 k^{2}_{y}/2m_y$, while the Rashba SO field is taken as $\vec{d}(\vec{k}) = \alpha ( \sqrt{m_d/m_y} k_y \vec{\hat{x}} - \sqrt{m_d/m_x} k_x \vec{\hat{y}} ) + \Delta \vec{\hat{z}}$, where α is the SO strength and $m_d = \sqrt{m_xm_y}$ is the geometric mean of the masses m_x and m_y along the x- and y- directions [50]. The model has been extended to include an energy parameter, $\Delta\geqslant 0$, which gives rise to a gap in the energy dispersion; this gapped anisotropic Rashba low-energy model could describe a gapped conduction band of the phosphorene monolayer [52]. We note that there is no evidence of the presence of a gap in this kind of systems (see for example [60, 61]), however being ours a theoretical approach, we artificially extend the model used in [51] to explore the effects of a $\Delta\hat{\sigma}_z$ term, given the analytical solvability of the related optical response problem. A gap could be caused by some interaction-induced phenomena, associated to spontaneous symmetry breaking mechanism or to an extrinsic symmetry breaking effect, like the presence of magnetic impurities [62]. Dyrdal et al considered a 2DEG with Rashba spin splitting coupled to a magnetic substrate via exchange interaction, which add a gap term to the Hamiltonian [63]. In a Dirac system, a finite nonlinear Hall conductivity was obtained whenever a gap is open in the system via electron-interaction-driven inversion-symmetry breaking [64]. Another invoked mechanism to open a gap is the photoinduced mass via intense circularly polarized light [65]. We observe that by taking $\Delta = 0$ or $\Delta\neq 0$ we can move from a model with TRS to a model with broken TRS and that the results for the first are readily derived from the analytical expressions.

Written in polar coordinates, the conduction ($\lambda = +$) and valence ($\lambda = -$) bands are $\varepsilon_{\lambda}(k,\theta) = \hbar^2 k^2 g^2(\theta)/2m_d + \lambda \sqrt{\alpha^2 k^2 g^2(\theta) + \Delta^2}$. The function $g(\theta) = [(m_d/m_x) \cos^2\theta + (m_d/m_y) \sin^2\theta]^{1/2}$ measures the separation of energy-constant curves along the direction θ in the k-space; note that $g(\theta) = 1$ when $m_x = m_y$, which corresponds to the well known case of a magnetized 2DEG with Rashba coupling in quantum wells of semiconductor heterostructures [63]. The energy difference between the bands is $\varepsilon_{+}(\vec{k}) - \varepsilon_{-}(\vec{k}) = 2d(\vec{k}) = 2 \sqrt{\alpha^2 k^2 g^2(\theta)+\Delta^2}$ and the constant energy-difference curve, $C_\mathrm{r} (\omega) = \{(k_x,k_y) \vert \varepsilon_{+}(\vec{k}) - \varepsilon_{-}(\vec{k}) = \hbar\omega \}$ is the ellipse with equation $\varepsilon_{0}(k_x,k_y) = [(\hbar\omega/2)^2 - \Delta^2] / 2\varepsilon_\mathrm{A}$, where $\varepsilon_\mathrm{A} = m_d\alpha^2 / \hbar^2$ is a characteristic energy associated to the Rashba interaction.

The shape of the band $\varepsilon_{-}(\vec{k})$ depends on the ratio $p = \Delta/\varepsilon_\mathrm{A}$. When p < 1, the band acquires a mexican hat shape, with a local maximum $-\Delta$ at the origin ${\textbf{k}} = {\textbf{0}}$ and two local minimums of value $\varepsilon_\mathrm{min} = -(\varepsilon^{2}_\mathrm{A}+\Delta^2)/2\varepsilon_\mathrm{A}$ at k-points lying on the ellipse $\varepsilon_{0}(k_x,k_y) = (\varepsilon^{2}_\mathrm{A} - \Delta^2)/2\varepsilon_\mathrm{A}$. Otherwise, the valence band develops only a minimum at the origin. This means that there are several distinct positions for the Fermi level: (i)  $\varepsilon_\mathrm{F}\gt\Delta$ (figure 4(a)), and then two Fermi contours are generated

$$\begin{align} k^{\pm}_\mathrm{F}(\theta) = \frac{1}{\alpha g(\theta)} \left[ \left( \sqrt{\varepsilon^{2}_\mathrm{A} +\Delta^2 + 2\varepsilon_\mathrm{A}\varepsilon_\mathrm{F} } \mp \varepsilon_\mathrm{A} \right)^2 - \Delta^2 \right]^{1/2}, \end{align} \tag{ 23 }$$

from equations $\varepsilon_+(k,\theta) = \varepsilon_-(k,\theta) = \varepsilon_\mathrm{F}$; (ii)  $|\varepsilon_\mathrm{F}|\lt\Delta$ (figure 4(b)), where there is only one Fermi contour $k^{-}_\mathrm{F}(\theta)$, lying on the valence band; and (iii) if p < 1, $\varepsilon_\mathrm{min} \lt \varepsilon_\mathrm{F} \lt -\Delta$ (figure 4(c)), where the Fermi lines arise only from the valence band through the equation $\varepsilon_-(k,\theta) = \varepsilon_\mathrm{F}$, with roots $k^{-}_\mathrm{F}(\theta)$ and

$$\begin{align} q^{-}_\mathrm{F}(\theta) = \frac{1}{\alpha g(\theta)} \left[\left( \varepsilon_\mathrm{A} - \sqrt{\varepsilon^{2}_\mathrm{A} +\Delta^2 + 2\varepsilon_\mathrm{A}\varepsilon_\mathrm{F} }\right)^2 - \Delta^2 \right]^{1/2}. \end{align} \tag{ 24 }$$

Note that the situation (i) or (iii) includes the gapless case $\Delta = 0$, the Fermi level being then positive or negative, respectively.

Zoom In Zoom Out Reset image size

Interestingly, the Fermi lines described by (23) and (24) are concentric ellipses vertically (horizontally) oriented if $m_y \gt m_x$ ($m_x \gt m_y$), see insets in figure 4. The energy separation of the bands at these lines is independent of the direction θ in k-space, taking the values $2d(k^{\pm}_\mathrm{F}(\theta)) = 2\left( \sqrt{\varepsilon^{2}_\mathrm{A} + \Delta^2 + 2\varepsilon_\mathrm{A}\varepsilon_\mathrm{F}} \mp \varepsilon_\mathrm{A} \right)$ or $2d(q^{-}_\mathrm{F}(\theta)) = 2\left(\varepsilon_\mathrm{A} - \sqrt{\varepsilon^{2}_\mathrm{A} + \Delta^2 + 2\varepsilon_\mathrm{A}\varepsilon_\mathrm{F}} \right)$, according to the position of the Fermi level.

The lowest (highest) energy of the spectrum of allowed interband transition will be denoted by $\hbar\omega_+$ ($\hbar\omega_-$), see figure 4. When $\vert\varepsilon_\mathrm{F}\vert \lt \Delta$ we have that the lowest possible transition occurs at the energy gap $\hbar\omega_+ = 2\Delta$. When $\varepsilon_\mathrm{F} \gt \Delta$ or $\varepsilon_\mathrm{min} \lt \varepsilon_\mathrm{F} \lt -\Delta$, we have $\hbar\omega_+ = 2d(k^{+}_\mathrm{F})$ or $\hbar\omega_+ = 2d(q^{-}_\mathrm{F})$, respectively. The highest possible energy transition is always given by $\hbar\omega_{-} = 2d(k^{-}_\mathrm{F})$. Moreover, the resonance curve $C_\mathrm{r}(\omega)$ is an ellipse with the same shape and orientation than the Fermi lines, but differing only in size given its frequency dependence. As a consequence, there will be only two critical points in the JDOS, given by the frequencies $\omega_+$ and $\omega_-$ at which the ellipse $C_\mathrm{r}(\omega)$ enters and leaves, respectively, the region of allowed transitions (shaded areas in the insets of figure 4), and which defines an absorption window $\hbar\omega_+(\varepsilon_\mathrm{F})\lt\hbar\omega\lt\hbar\omega_-(\varepsilon_\mathrm{F})$. All this is apparent in the JDOS

$$\begin{align} J_{+-}(\omega;\varepsilon_\mathrm{F})& = \frac{\hbar\omega}{8\pi\alpha^2} \Theta(1-\vert \eta(\omega,\varepsilon_\mathrm{F}) \vert), \ \ \ \eta(\omega,\varepsilon_\mathrm{F}) \nonumber\\ & = \frac{\omega - (\omega_{-} + \omega_{+})/2}{(\omega_{-} - \omega_{+})/2}\,, \end{align} \tag{ 25 }$$

which is displayed as a color map in figure 5, for a system having p < 1. For a given value of $\varepsilon_\mathrm{F}$ the color gradation shows the linear dependence on the exciting frequency.

Zoom In Zoom Out Reset image size

According to the formula (3), the induced spin current response function of the anisotropic Rashba model becomes

$$\begin{align} \sigma^{\ell}_{ij}(\omega)& = - \delta_{\ell z} \frac{e}{8\pi} \frac{1}{2\varepsilon_{A}} \left[ \varepsilon_{ijz} - i \frac{m_d}{m_i}\! \left(\frac{2\Delta}{\hbar\tilde{\omega}}\right) \delta_{ij} \right] \nonumber\\ &\quad \times \left[ A(\tilde{\omega}) + \frac{1}{2}\hbar(\omega_{-}-\omega_{+}) \right], \end{align} \tag{ 26 }$$

where

$$\begin{equation} A(x) = \frac{\hbar x}{4}\left[1-\left(\frac{2\Delta}{\hbar x}\right)^2\right] \operatorname{ln}\left[\frac{(x+\omega_{+})(x-\omega_{-})} {(x-\omega_{+})(x+\omega_{-})}\right]. \end{equation} \tag{ 27 }$$

Remarkably, the linear spin currents generated in this system will have electrons with out-of-plane spin orientations only. This a consequence of the breaking of TRS by the term $d_z({\textbf{k}})$, which is a non null constant in the present model. This makes the product of matrix elements $\langle \lambda | \mathcal{\hat{J}}^{\ell}_{i} | -\lambda \rangle \langle -\lambda | \hat{v}_j | \lambda \rangle$ an odd (even) function in k-space when $\ell = x,y$ ($\ell = z$). Thus, the Kubo expression (3) integrates to a non zero value when $\ell = z$ only. The longitudinal components of the spin current conductivity (figure 6(a)) are proportional to the gap parameter, $\sigma^{z}_{ii}(\omega)\propto\Delta$, and therefore they vanish for the gapless case. Moreover, these diagonal components are inversely proportional to $\sqrt{m_i}$, such that $m_x\sigma^z_{xx}(\omega) = m_y\sigma^z_{yy}(\omega)$, as a consequence of the mass anisotropy. On the other hand, the Hall components are nonzero, regardless of the value of Δ, indicating the generation of a spin Hall effect in the gapped or ungapped system. In addition, $\sigma^{\ell}_{xy}(\omega) = -\sigma^{\ell}_{yx}(\omega)$, as can be seen in figure 6(b). The effect of the position of the Fermi level with respect to the gap, manifests through the critical energies $\hbar\omega_{\pm}(\varepsilon_\mathrm{F})$. The variation of Fermi energy leads mainly to a change of the window $\hbar (\omega \_{\textrm{ - }}{\omega _ + })$, just like in JDOS.

Zoom In Zoom Out Reset image size

If ${\textbf{E}}^{\omega}_0 = E_0^{\omega}(\cos\varphi{\hat{\textbf{x}}}+\sin\varphi{\hat{\textbf{y}}})$ is the amplitude of the external field, the spin current can be written as the sum of a component along ${\textbf{E}}^{\omega}_0$ and a component perpendicular to it,

$$\begin{align} {\boldsymbol{\cal{J}}^z}(\omega) & = [\cos^2\varphi\,\sigma^z_{xx}(\omega)+\sin^2\varphi\, \sigma^z_{yy}(\omega)]{\textbf{E}}^{\omega}_0+[\sin\varphi\cos\varphi (\sigma^z_{xx}(\omega)\nonumber\\ & \quad -\sigma^z_{yy}(\omega))+\sigma^z_{xy}(\omega)]({\textbf{E}}^{\omega}_0\times {\hat{\textbf{z}}}), \end{align} \tag{ 28 }$$

which reduces to ${\boldsymbol{\cal{J}}^z}(\omega) = \sigma^z_{xy}(\omega)({\textbf{E}}^{\omega}_0\times {\hat{\textbf{z}}})$ when $\Delta = 0$. Expression (28) suggests how the induced spin current could be manipulated through the frequency dependence of the tensor $\sigma^z_{ij}$ and the direction of the applied in-plane electric field.

For the SH spin current conductivity we obtain the following expressions from (6) (no sum over repeated indices is implied),

$$\begin{align} \sigma^{\ell,(2\omega)}_{ijj}(\omega) = \delta_{\ell j} \varepsilon_{ijz} \frac{e^2 / 8\pi}{(\hbar\tilde{\omega})^2} \frac{1}{k_\mathrm{A}} \left( \frac{m_d}{m_j} \right)^{1/2} F(\tilde{\omega}), \end{align} \tag{ 29 }$$

$$\begin{align} \sigma^{\ell,(2\omega)}_{iii}(\omega)& = \frac{e^2/8\pi}{(\hbar\tilde{\omega})^2} \frac{1}{k_\mathrm{A}} \left(\frac{m_{d}}{m_i}\right)^{3/2} \bigg\{\varepsilon_{\ell iz}\ G(\tilde{\omega}) + 2i \left(1-\delta_{\ell z}\right)\nonumber\\ & \quad \times \left(\frac{2\Delta}{\hbar\tilde{\omega}}\right) \Big[A(\tilde{\omega}) - A(2\tilde{\omega})\Big] \bigg\}, \end{align} \tag{ 30 }$$

$$\begin{align} \sigma^{\ell,(2\omega)}_{iji}(\omega) & = \sigma^{\ell,(2\omega)}_{iij}(\omega) = \frac{e^2/16\pi}{(\hbar\tilde{\omega})^2} \frac{1}{k_\mathrm{A}} \left(\frac{m_d}{m_i}\right)^{1/2}\nonumber\\ & \quad \times\bigg\{\delta_{\ell i}\varepsilon_{ijz}\ H(\tilde{\omega}) + 2i (1-\delta_{\ell z})\left(\frac{2\Delta}{\hbar\tilde{\omega}}\right)\nonumber\\ & \quad \times \Big[A(\tilde{\omega}) - A(2\tilde{\omega})\Big]\bigg\}, \end{align} \tag{ 31 }$$

where $k_\mathrm{A} = \varepsilon_\mathrm{A}/\alpha = m_d\alpha / \hbar^2$, A(x) is given by (27), and

$$\begin{align} F(\tilde{\omega}) & = \Bigg\{1 + \frac{1}{4}\left[1-\left(\frac{2\Delta}{\hbar\tilde{\omega}}\right)^2 \right] \Bigg\} A(\tilde{\omega}) \nonumber\\ & \quad - \left[1-\left(\frac{2\Delta}{2\hbar\tilde{\omega}}\right)^2\right]A(2\tilde{\omega}) +\frac{1}{8}\hbar(\omega_{-}-\omega_{+}), \end{align} \tag{ 32 }$$

$$\begin{align} G(\tilde{\omega})& = \Bigg\{1-\frac{3}{4} \left[1- \left(\frac{2\Delta}{\hbar\tilde{\omega}} \right)^2 \right] \Bigg\} A(\tilde{\omega})\nonumber\\ & \quad - \Bigg\{4 - 3\left[ 1-\left(\frac{2\Delta}{2\hbar\tilde{\omega}}\right)^2 \right] \Bigg\} A(2\tilde{\omega}) - \frac{3}{8}\hbar (\omega_{-} - \omega_{+}), \end{align} \tag{ 33 }$$

$$\begin{align} H(\tilde{\omega})& = \Bigg\{2 - \frac{1}{2}\left[ 1 - \left(\frac{2\Delta}{\hbar\tilde{\omega}}\right)^2 \right] \Bigg\} A(\tilde{\omega})\nonumber\\ & \quad - 2 \Bigg\{2 - \left[ 1 - \left(\frac{2\Delta}{2\hbar\tilde{\omega}}\right)^2 \right] \Bigg\} A(2\tilde{\omega}) - \frac{1}{4}\hbar(\omega_{-} - \omega_{+}). \end{align} \tag{ 34 }$$

A number of conclusions can be derived from these expressions. In contrast to the linear response, the above expressions imply that there is no z-polarized SH spin current induced in the system, $\sigma^{z,(2\omega)}_{ijl}(\omega) = 0$. On the other hand, equation (30) shows that if $\Delta = 0$ the longitudinal response $\sigma^{\ell,(2\omega)}_{iii}$ describes a spin current with the spin oriented perpendicularly to the electric field (and to the spin current), while one with parallel spin orientation and current is possible when $\Delta\neq 0$. The Hall components $\sigma^{\ell,(2\omega)}_{ijj}$ (29) generate spin currents with spin always oriented in the direction of the field (normal to the spin current), regardless of the presence of an energy gap. As for the components in (31), $\sigma^{\ell,(2\omega)}_{iij}(\omega)$ and $\sigma^{\ell,(2\omega)}_{iji}(\omega)$, they are associated to a nonlinear current with spin orientation parallel to it when $\Delta = 0$, otherwise such an orientation is not fixed for $\Delta\neq 0$.

The terms in braces in (29)–(31) involve the masses through the combination $m_d = \sqrt{m_xm_y}$ only. Thus, the in-plane anisotropy becomes apparent in every non-zero tensor component through a factor of the type $(m_x / m_y)^{\pm \nu/4}$, with $\,\nu = 1$ or 3. Figure 7 shows the frequency dependence of the in-plane polarized SH components $\sigma^{x,(2\omega)}_{ijl}(\omega)$ and $\sigma^{y,(2\omega)}_{ijl}(\omega)$. As expected, the spectral structure around $\hbar\omega_{\pm}$ is now accompanied by new features around the subharmonics $\hbar\omega_{\pm}/2$. The overall structure can be modified through Fermi energy variation, given the behavior of the functions $\hbar\omega_{\pm}(\varepsilon_\mathrm{F})$ observed in figure 5.

Zoom In Zoom Out Reset image size

The SH spin current response (29)–(31) presents characteristic differences with respect to the linear response, as is the case in the model of section 3. Diagonal components are present even in absence of an energy gap, subharmonic structure appears in the spectrum, and the spin polarization is no longer out-of-plane, so that a nonlinear spin Hall effect with in-plane spin orientation is generated. These features may be useful in nonlinear optical spintronic devices.