Strong-field physics in three-dimensional topological insulators

We investigate theoretically the strong-field regime of light-matter interactions in the topological-insulator class of quantum materials. In particular, we focus on the process of nonperturbative high-order harmonic generation from the paradigmatic three-dimensional topological insulator bismuth selenide (${\mathrm{Bi}}_2{\mathrm{Se}}_3$) subjected to intense midinfrared laser fields. We analyze the contributions from the spin-orbit-coupled bulk states and the topological surface bands separately and reveal a major difference in how their harmonic yields depend on the ellipticity of the laser field. Bulk harmonics show a monotonic decrease in their yield as the ellipticity increases, in a manner reminiscent of high harmonic generation in gaseous media. However, the surface contribution exhibits a highly nontrivial dependence, culminating with a maximum for circularly polarized fields. We attribute the observed anomalous behavior to (i) the enhanced amplitude and the circular pattern of the interband dipole and the Berry connections in the vicinity of the Dirac point and (ii) the influence of the higher-order, hexagonal warping terms in the Hamiltonian, which are responsible for the hexagonal deformation of the energy surface at higher momenta. The latter are associated directly with spin-orbit-coupling parameters. Our results thus establish the sensitivity of strong-field-driven high harmonic emission to the topology of the band structure as well as to the manifestations of spin-orbit interaction.

Figure 1

(a) Crystal structure of ${\mathrm{Bi}}_2{\mathrm{Se}}_3$, comprising alternating Bi and Se layers, stacked along the $z$ direction. Five consecutive layers form a quintuple layer (QL, cp. red rectangle), the building block of the lattice. Each QL comprises five atoms: two equivalent Bi sites, two equivalent Se sites [Se(1)], and a third Se atom, Se(2), which assumes the role of an inversion center. The hexagonal lattice constants are $a=4.14$ Å and $c=28.70$ Å. The space-dependent wave function of the surface states (squared magnitude), ${\mathrm{\Psi }}_{\mathrm{Surf}}({\mathit{k}}_{\parallel };z)$, is sketched as a red-shaded surface and illustrates the employed boundary conditions. (b) Schematic representation of the $C_{3v}$-symmetric ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ (111) surface (rhombohedral convention), exposing a top Se layer and underlying Bi and Se' layers. (c) Sketch of the 3D Brillouin zone of bulk ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ (black) with the four time-reversal-invariant points indicated ($\mathrm{\Gamma },L,F,Z$). The projected 2D BZ of the (111) surface is shown as a red hexagon, with labeled high-symmetry points $\overline{\mathrm{\Gamma }},\overline{K},\overline{M}$.

Figure 2

(a) Energy dispersion of the bulk states $E_{\mathrm{B}}^{\pm }(k_x,k_y,k_z)$ along the high-symmetry directions $\overline{\mathrm{\Gamma }M}$ and $\overline{\mathrm{\Gamma }K}$ in inverse space, shown for different $k_z$ values. The black lines correspond to the plane $k_z=0$, projected band dispersion curves pertaining to increasing $k_z$ are given in progressively lighter blue colors, whereby the $k_z$ is varied by $\mathrm{\Delta }{k}_z=1.5\times {10}^{-3}{\AA }^{-1}$. The abscissa covers the range from $k_x=0$ to $k_x=0.36{\AA }^{-1}$ ($\mathrm{\Gamma }K$ direction) and from $k_y=0$ to $k_y=0.31{\AA }^{-1}$ ($\mathrm{\Gamma }M$ direction). (b) Energy dispersion of the surface modes $E_{\text{2D}}^{\mathrm{S};(\pm )}(k_x,k_y)$ resulting from the employed TBM model, given in the disk defined by $k_{\parallel }\le 0.4{\AA }^{-1}$. Near the $\overline{\mathrm{\Gamma }}$ point, the dispersion is nearly linear. At higher momenta, the hexagonal warping effect is seen as a consequence of the higher-order contributions. (c), (d) Spin polarization of the surface lower (c) and upper (d) Dirac cones over a selected portion of the BZ. The white arrows indicate the magnitude and direction of the in-plane polarization, whereas the color coding corresponds to the magnitude of the spin polarization in $\widehat{z}$ direction (out of plane). We note that the Berry curvature in momentum space follows a similar pattern as the out-of-plane spin polarization.

Figure 3

HHG spectra of ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ driven by circularly polarized fields, for both bulk (a) and surface states (b). The relative orientation of the (111) surface (real space) and the MIR polarization vector is sketched in the inset, whereby the propagation direction points toward the reader. The MIR pulse is left-circularly polarized with $I_0=0.0025\mathrm{TW}/{\mathrm{cm}}^2$ and a FWHM duration of 12 cycles. The helicity of the emitted harmonics is encoded in color: left, corotating orders are shown in magenta (full lines), right, counter-rotating in cyan (dashed lines). Bulk and surface states obey different selection rules: $\omega =(6n\pm 1){\omega }_0$ for the bulk (a) versus $\omega =(3n\pm 1){\omega }_0$ for the surface (b). In both calculations, the dephasing time is set at $T_2=1.25$ fs.

Figure 4

(a) Sketch of the excitation geometry, where the major axis of the MIR ellipse (orange) remains perpendicular to the mirror axis (dashed line) of the crystal throughout the measurement. (b)–(d) show the calculated ellipticity response for harmonic orders 5, 7, and 9. The calculations (12-cycle Gaussian pulse, $I_0=0.0085\mathrm{TW}/{\mathrm{cm}}^2$, ${\lambda }_{\mathrm{MIR}}$=7.5 $\mu \mathrm{m}$, $T_2=1.25$ fs) pertain to the contributions from the bulk [(B), dashed lines, diamonds] and from the surface states [(S), solid lines, circles].

Figure 5

(a), (b) Real (a) and imaginary (b) parts of the interband transition matrix element ${\mathit{d}}_{vc}\left({\mathit{k}}_{\parallel }\right)$ between the surface valence and conduction bands. The streamlines (white) indicate the local direction of the vector fields in momentum space. (c) HHG spectra of the TSSs illuminated by a LCP MIR pulse with $I_0=0.0075\mathrm{TW}/{\mathrm{cm}}^2$, ${\lambda }_{\mathrm{MIR}}\sim 7.5\mu \mathrm{m}$, $T_2=1.25$ fs, and a duration of 12 cycles. The blue (dashed) curve corresponds to a calculation which considers only the real part of the interband dipole ${\mathit{d}}_{cv}\left(\mathit{k}\right)$, i.e., leading term in Eqs. (25) and (26). The red (full) curve pertains to the full expression for the dipole.

Figure 6

Upper row: Momentum-resolved temporal snapshot of the population distribution in the upper Dirac cone at different time points during the interaction with a LCP MIR driving field ($I_0=0.0025\mathrm{TW}/{\mathrm{cm}}^2$, 12 cycles), shown at three different time points (a)–(c) of the pulse envelope. The top panels show the $x$ component of the electric-field amplitude of the MIR driving pulse (in atomic units). Bottom row: (d)–(f) show the corresponding population distribution for one of the degenerate components of the bulk conduction states (${\psi }_{\mathrm{B},\nu =1}^{+}$) for the same conditions as for the TSSs. Note that ${\rho }_{cc}^{{\mathit{k}}_{\parallel }}\in [0,1]$ for the TSSs and ${\rho }_{c(\nu =1)c(\nu =1)}^{{\mathit{k}}_{\parallel }}={\rho }_{c(\nu =2)c(\nu =2)}^{{\mathit{k}}_{\parallel }}\in [0,1/2]$ for the BSs.

Figure 7

HHG spectra emitted from the surface states driven by a 12-cycle left-circularly polarized pulse with a peak intensity of $I_0=0.004\mathrm{TW}{\mathrm{cm}}^{-2}$, whereby one of the TBM parameters $A_{12}$ (a) or $A_{14}$ (b) is varied (s. legend). The spectra corresponding to the parameters listed in Table 1 are shown in black (full lines).

Figure 8

Ellipticity dependence of the HHG yields of HOs 5 (full lines), 7 (dashed lines), and 9 (dash-dotted lines) of a MIR pulse with ${\lambda }_{\mathrm{MIR}}=7.5\mu \mathrm{m}$ and different driving peak intensities: $I_0=0.0025\mathrm{TW}/{\mathrm{cm}}^2$ in (a) and (d), $I_0=0.0045\mathrm{TW}/{\mathrm{cm}}^2$ in (b) and (e), and $I_0=0.01\mathrm{TW}/{\mathrm{cm}}^2$ in (c) and (f). The upper row corresponds to the surface states, bulk states are plotted in the bottom row. As illustrated in the inset of panel (e), the main axis of the polarization ellipse $\widehat{\mathbf{e}}$ (orange ellipse) is set perpendicular to the mirror plane (${\widehat{\sigma }}_{\mathrm{refl}}^{\left(y\right)}$, dashed line). The pulse has a Gaussian profile with a FWHM duration of 12 cycles.

Figure 9

Simplified representation of the lattice structure illustrating the nearest-neighbor vectors $\pm \mathit{a}$ and $\pm \mathit{b}$.

Figure 10

HHG spectra of ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ driven by linearly polarized fields, for both bulk and surface states. The relative orientation of the (111) surface (real space) and the MIR polarization vector is sketched in the insets. (a) and (c) show spectra of the bulk states for a linearly polarized MIR laser pulse oriented along the $\overline{\mathrm{\Gamma }M}$ and $\overline{\mathrm{\Gamma }K}$ directions, respectively. All emitted harmonics are parallel with respect to ${\mathit{E}}_{\mathrm{MIR}}\left(t\right)$ (cyan, full). The laser field has a peak intensity of $I_0=0.002\mathrm{TW}/{\mathrm{cm}}^2$ and a FWHM duration of 15 cycles. (b) and (d) show the corresponding spectra for the surface states. The polarization of the even harmonics flips from parallel (cyan, full) when ${\mathit{E}}_{\mathrm{MIR}}\left(t\right)\parallel \overline{\mathrm{\Gamma }M}$ (i.e., ${\mathit{E}}_{\mathrm{MIR}}\left(t\right)\parallel {\widehat{\sigma }}_{\mathrm{refl}}^{\left(y\right)}$) to perpendicular (magenta, dashed) with respect to the driving field when ${\mathit{E}}_{\mathrm{MIR}}\left(t\right)\parallel \overline{\mathrm{\Gamma }K}$ (i.e., ${\mathit{E}}_{\mathrm{MIR}}\left(t\right)\perp {\widehat{\sigma }}_{\mathrm{refl}}^{\left(y\right)}$). In all calculations, the dephasing time is set at $T_2=1.25$ fs.

Figure 11

(a) Vector field plot of the difference of the Berry connections of upper and lower Dirac cones of the TSS. The color quantifies the absolute magnitude of $\mathrm{\Delta }{\mathit{\xi }}_{cv}\left({\mathit{k}}_{\parallel }\right)$. The streamlines (white) indicate the local direction of the vector fields in momentum space. (b), (c) Stream plots of the derivatives of the $x$ (b) and the $y$ components (c) of the phase of the interband dipole matrix element ${\mathit{d}}_{cv}\left({\mathit{k}}_{\parallel }\right)$.

Figure 12

HHG spectra emitted from the bulk states driven by a 12-cycle left-circularly polarized pulse with a peak intensity of $I_0=0.004\mathrm{TW}{\mathrm{cm}}^{-2}$, whereby one of the TBM parameters $A_{12}$ (a) or $A_{14}$ (b) is varied (s. legend). The spectra corresponding to the parameters listed in Table 1 are shown in black (full line).