Enhanced higher harmonic generation from nodal topology

TL;DR

This paper demonstrates that the unique topological features of nodal loop semimetals enhance higher harmonic generation in the THz regime, offering a new physical property linked to their topology with potential technological applications.

Contribution

It is the first work to show that nodal linkages in semimetals protect and enhance higher harmonic generation, connecting topology with nonlinear optical responses.

Findings

Nodal linkages enforce non-monotonic intra-band responses.

Enhanced HHG observed in known nodal materials.

Potential applications in THz technology and material characterization.

Abstract

Among topological materials, nodal loop semimetals (NLSMs) are arguably the most topologically sophisticated, with their valence and conduction bands intersecting along arbitrarily intertwined nodes. But unlike the well-known topological band insulators with quantized edge conductivities, nodal loop materials possess topologically nontrivial Fermi surfaces, not bands. Hence an important question arises: Are there also directly measurable or even technologically useful physical properties characterizing nontrivial nodal loop topology? In this work, we provide an affirmative answer by showing, for the first time, that nodal linkages protect the higher harmonic generation (HHG) of electromagnetic signals. Specifically, nodal linkages enforce non-monotonicity in the intra-band semi-classical response of nodal materials, which will be robust against perturbations preserving the nodal topology. These nonlinearities distort incident radiation and produce higher frequency peaks in the teraHertz (THz) regime, as we quantitatively demonstrate for a few known nodal materials. Since THz sources are not yet ubiquitous, our new mechanism for HHG will greatly aid applications like material characterization and non-ionizing imaging of object interiors.

Figures

b). To explore whether the absence of linkages penalize HHG, we first construct its tight-binding model according to a DFT-fitted ansatz detailed in the Methods [Sup ](#bib.bib40). As evident in Figs. 4d,e, the resultant $\langle\mathbb{v}\rangle_{\mathbb{E}}$ response non-linearity and hence HHG for the three touching NLs of YH3 is even weaker than that of the single NL Ti3Al, even before spin-orbit coupling effects gap out NLs into elogated Fermi tubes. The reason for the poor HHG is that in this case, the multiple NLs “smooth” out the energy dispersion around the nodal touchings, thereby re

Equations

$$\mathbb{J}=\int d^3\mathbb{k}g(\epsilon (\mathbf{k}))⟨k|J|\mathbf{k}⟩,$$

$$\frac{dg}{dt}=\frac{\partial g}{\partial t}+\dot{\mathbb{k}}\cdot \frac{\partial g}{\partial \mathbf{k}}=\frac{f-g}{\tau },$$

$$\dot{r}$$

$$\dot{\mathbb{k}}$$

$$g(\mathbb{k},t)=F\left(\epsilon \left(\mathbb{k}-{\int }_{-\infty }^t\mathbb{E}(t^{\prime })d{t}^{\prime }\right)\right)=F\left(\epsilon \left(\mathbb{k}-\mathbb{A}(t)\right)\right),$$

$$\mathbb{J}[\mathbb{A}(t)]$$

$$H_{\text{NL}}(\mathbb{k})=\left(\cos k_x+\cos k_y+\cos k_z-r\right){\sigma }_x+\sin k_z{\sigma }_y+M{\sigma }_z,$$

$${\epsilon }_{\text{NL}}(\mathbb{k})=\pm \sqrt{M^2+{\sin }^2{k}_z+{\left(\cos k_x+\cos k_y+\cos k_z-r\right)}^2},$$

$${\mathbf{\Omega }}_{\text{NL}}(\mathbb{k})=\frac{M\cos k_z}{2|{\epsilon }_{\text{NL}}(\mathbb{k}){|}^3}\left[\begin{array}{c}\sin k_y\\ -\sin k_x\\ 0\end{array}\right]$$

$$H_{\text{Hopf}}(\mathbf{k})=Re[f(v,w)]{\sigma }_x+Im[f(v,w)]{\sigma }_y+M{\sigma }_z$$

$${\epsilon }_{\text{Hopf}}(\mathbb{k})=\pm \sqrt{M^2+{\chi }^2+4(e_1{e}_2-e_3{e}_4{)}^2},$$

$$\frac{d{J}_i}{d{A}_j}$$

$${\epsilon }_{\text{link}}(k_x,k_y)$$

$$J_j(A_i(t))$$

$$H_{\text{Co2MnGa}}=Re[f(\mathbb{k})]{\sigma }_1+Im[f(\mathbb{k})]{\sigma }_3,$$

$$u_{\text{even}}$$

$$u_{\text{odd}}$$

$$\begin{array}{rl}1& =\int \frac{{\mathrm{d}}^3{\mathbf{k}}^{\prime }}{(2\pi {)}^3}\frac{2\pi }{\hslash }|⟨{\mathbf{k}}^{\prime }|\hat{V}|\mathbf{k}⟩{|}^2\delta ({\epsilon }_{\mathbf{k}}-{\epsilon }_{{\mathbf{k}}^{\prime }})\left({\tau }_{\mathbf{k}}^{(i)}-\frac{v_{{\mathbf{k}}^{\prime }}^{(i)}}{v_{\mathbf{k}}^{(i)}}{\tau }_{{\mathbf{k}}^{\prime }}^{(i)}\right).\end{array}$$

$$J_{\hat{e}}(\mathbb{A})$$

$$\sum _{\alpha }\frac{d}{d{k}_{\parallel }^{\alpha }}\int \left[\epsilon \left({\mathbb{k}}^{\alpha }+\Delta {\mathbb{k}}_{\parallel }^{\alpha }+\mathbb{A}\right)-\epsilon \left({\mathbb{k}}^{\alpha }-\Delta {\mathbb{k}}_{\parallel }^{\alpha }+\mathbb{A}\right)\right]d{\mathbf{k}}_{\perp }$$

$$H_{\text{alt}}$$

$${\epsilon }_{1,2}=\pm \sqrt{{\sin }^2{k}_z+(\sqrt{{\sin }^2{k}_x+{\sin }^2{k}_y}-b{)}^2},$$

$${\epsilon }_{3,4}=\pm \sqrt{{\sin }^2{k}_z+(\sqrt{{\sin }^2{k}_x+{\sin }^2{k}_y}+b{)}^2}$$

$$\varepsilon_{\text{Ti${}_{3}$Al}}(\mathbb{k})=\frac{1}{2}\bigg{(}h_{1}+h_{2}+\sqrt{(h_{1}-h_{2})^{2}+h^{2}}\bigg{)}$$

$${\epsilon }_{\text{Ti3Al}}{|}_{\sqrt{k}\le 0.16}$$

$${\epsilon }_{\text{Ti3Al}}{|}_{\sqrt{k}>0.16}$$

$${\epsilon }_{\text{YH3}}(\mathbb{k})=\sqrt{(g_1^2+h_1^2)(g_2^2+h_2^2)(g_3^2+h_3^2)}$$

$$g_1=\sin (k_z),g_2=\sin (k_x),g_3=\sin (k_y)$$

$$h_1=a_1(r_1(\cos k_x){)}^3+s_1(\cos k_y){)}^3+t_1(\cos z){)}^3-m_1)$$

$$h_2=a_2(\cos k_x+\cos k_y+\cos k_z-m_2)$$

$$h_3=a_3(\cos k_x+\cos k_y+\cos k_z-m_3)$$

$$H(\mathbf{k})=\left(\begin{array}{cccccc}{\xi }_1^d& 0& 0& {\xi }_{11}^{dp}& 0& {\xi }_{13}^{dp}\\ 0& {\xi }_2^d& 0& 0& {\xi }_{22}^{dp}& {\xi }_{23}^{dp}\\ 0& 0& {\xi }_3^d& {\xi }_{31}^{dp}& {\xi }_{32}^{dp}& 0\\ {\xi }_{11}^{dp}& 0& {\xi }_{31}^{dp}& {\xi }_1^p& {\xi }_{12}^p& {\xi }_{31}^p\\ 0& {\xi }_{22}^{dp}& {\xi }_{32}^{dp}& {\xi }_{12}^p& {\xi }_2^p& {\xi }_{23}^p\\ {\xi }_{13}^{dp}& {\xi }_{23}^{dp}& 0& {\xi }_{31}^p& {\xi }_{23}^p& {\xi }_3^p\end{array}\right)(\mathbf{k}),$$

Full text

Appendix SI Validity of the ballistic approximation

To justify that the ballistic regime $\Omega\tau\gg 1$ we assumed indeed corresponds to the Terahertz regime in reasonably clean samples of known nodal materials, we discuss below how the scattering time $\tau$ can be estimated. Since experimental data on $\tau$ is rare, especially for nodal materials which are relative novel, we shall content ourselves with the following estimation approach. First, consider a representative single nodal loop (e.g. Ca3P2 Xie et al. (2015)). For small loop radius, one can work with a continuum model and hence in toroidal coordinates. Consider Coulomb impurities as the dominant scatterers, and take the scatterer concentration to be $\sim 10^{24}\,$m*-3*, i.e. of the same order as that of Weyl semimetals such as Cd3As2 Ramakrishnan et al. (2015). For screening, we take the long-wavelength limit to obtain a Thomas-Fermi dielectric function $\varepsilon(\bm{q})=1+\frac{q_{s}^{2}}{\bm{q}^{2}}$, where $q_{s}^{2}=e^{2}\nu(E_{F})/\varepsilon_{0}$ is the square of the inverse screening length, with $\nu$ the density of states, $\varepsilon_{0}$ the vacuum permittivity, and $e$ the elementary charge. To obtain a rough estimate for the relaxation time, one can proceed by assuming that the relaxation time depends only on the Fermi energy. Denoting it as $\tau_{\bm{k}}$, it can can then be calculated Schliemann and Loss (2003) numerically by solving

[TABLE]

for anisotropic band dispersion $\varepsilon_{\bm{k}}$. In the above, $\hat{V}$ is the screened impurity potential, and the superscript $(i)$ indicates the direction of the applied field. The resulting relaxation time depends on the Fermi energy, but its values can be placed within the range $\tau\,\sim\,5\times 10^{-13}\,s$, which is similar to that of Graphene as obtained via RF admittance measurements Pallecchi et al. (2011). Then, with a terahertz driving field of frequency $\sim 100\,$THz, one has, roughly, $\Omega\tau\,\sim\,50$, validating the ballistic approximation.

Appendix SII Proof of the non-monotonicity of the response of a topological linkage

To show that singularities of the $\mathbb{v}=\nabla_{\mathbb{k}}\varepsilon$ vector field can indeed ensure significant non-linearity in the current response $\mathbb{J}(\mathbb{A})$, we approximate the Fermi-Dirac distribution $F(\varepsilon)$ as a step function777This introduces negligible errors at room temperature., and rewrite the velocity current contribution in a chosen $\hat{e}$ direction as

[TABLE]

where $k_{\parallel}=\mathbb{k}_{\parallel}\cdot\hat{e}=\mathbb{k}\cdot\hat{e}$ and ${\bf k}_{\perp}=\mathbb{k}-\mathbb{k}_{\parallel}$ are the components of ${\bf k}$ parallel and perpendicular to $\hat{e}$. The BZ is sliced into strips parallel to $\hat{e}$ and indexed by $\mathbb{k}_{\perp}$, such that the thin nodal Fermi tube is mathematically broken into thin cross sections $\alpha$ in each $\mathbb{k}_{\perp}$ strip. Each $\alpha$ has a width of $2|\Delta\mathbb{k}_{\parallel}^{\alpha}|$ that is determined via $\varepsilon(\mathbb{k}^{\alpha}\pm\Delta\mathbb{k}_{\parallel}^{\alpha})=\mu$, where $\mathbb{k}^{\alpha}$ is the central position of $\alpha$. The current is thus proportional to the sum of all energy differences $\Delta\varepsilon_{\alpha}({\bf k}_{\perp}^{\alpha})=\varepsilon\left(\mathbb{k}^{\alpha}+\Delta\mathbb{k}_{\parallel}^{\alpha}+\mathbb{A}\right)-\varepsilon\left(\mathbb{k}^{\alpha}-\Delta\mathbb{k}_{\parallel}^{\alpha}+\mathbb{A}\right)$ between the two “sides” of the Fermi tube. This rewriting of current response as shifted boundary terms actually remains valid in arbitrarily high dimensions Lee et al. (2018b), and can be elegantly expressed in terms of generalized Bianchi identities Petrides et al. (2018).

Nonlinear response in $J_{\hat{e}}$ can be investigated through the gradient $\frac{dJ_{\hat{e}}}{dA_{\hat{e}}}=$

[TABLE]

where $\hat{e}\cdot\mathbb{v}_{F}$ is the component of the Fermi velocity of NL $\alpha$ along the applied field. More generically, the response tensor can be obtained by replacing the second derivative by the Hessian, and projecting the Fermi velocity along the applied field. To summarize, Eq. S3 expresses the current response gradient as the sum of the dispersion curvature over all occupied states, inversely weighted by the Fermi velocity.

Appendix SIII Alternative model for the single nodal loop

In the main text, we have mentioned that a single NL has a response nonlinearity that is “unprotected”. This is because its dispersion $\varepsilon(\mathbb{k})$ can be deformed, without affecting its ring of nodes, such that its velocity vector field $\mathbb{v}=\nabla_{\mathbb{k}}\varepsilon$ exhibits very different cancellation extents. To illustrate this, we introduce an alternative toy model for the single NL:

[TABLE]

where $0<b<1$. This 4-band model possesses energy bands

[TABLE]

[TABLE]

The intersection of the bands $\varepsilon_{1,2}$ of Eq. S5 forms a single NL in the $k_{x}$-$k_{y}$ plane. But due do their proximity with bands $\varepsilon_{3,4}$, band intersections also exist at the origin, imitating the effect of another node at the center of the NL in the $k_{x}$-$k_{y}$ plane despite the absence of another linked NL. A comparison with the dispersion of the 2-band single NL model is shown in Fig. S2. That said, the $\mathbb{v}$ singularity at the center can be easily deformed away from the plane of the NL, thereby reducing the non-linear response, unlike in the Hopf NL where the central singularity is an unavoidable feature of the nodal linkage.

Appendix SIV Further details on the nodal materials and their lattice models

SIV.1 Ti3Al

The energy dispersion of Ti3Al is modeled by the expression Zhang et al. (2018)

[TABLE]

where $h_{i}=A_{i}(k_{x}^{2}+k_{y}^{2})+B_{i}k_{z}^{2}+M_{i}$ and $h=2Ck_{z}$ with parameter values $A_{1}=-9.66\text{eV}$Å2, $A_{2}=11.37\text{eV}$Å2, $B_{1}=36.22\text{eV}$Å2, $B_{2}=-25.71\text{eV}$Å2, $M_{1}=0.12\text{eV}$, $M_{2}=-0.52\text{eV}$ and $C=22.34\text{eV}$Å2.

To conveniently represent its cross section in the $k_{z}=0$ plane, it is also useful to use its piecewise polynomial approximation

[TABLE]

where $k=k_{x}^{2}+k_{y}^{2}$ and $c_{1}=0.0911$, $c_{2}=-4.0615$, $c_{3}=-0.2339$, $c_{4}=9.715$, $c_{5}=-41.8543$ $c_{6}=81.1788$ and $c_{7}=-59.1387$. Units of energy and momentum $k$ are in $eV$ and $\AA^{-1}$ respectively.

SIV.2 YH3

YH3 consists of three intersecting NLs along each orthogonal plane in the 3D BZ, and can be modeled by the dispersion

[TABLE]

where the factors $g_{i}^{2}+h_{i}^{2}$ each yield an individual NL Zhou et al. (2018). They are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where DFT-fitted parameters Shao et al. (2018) are given by $m_{1}=2.99$, $a_{1}=2$, $r_{1}=s_{1}=t_{1}=1.032$, $n_{1}=3$, $m_{2,3}=2.96$ and $a_{2,3}=4$. Units of energy and momentum $k$ are $eV$ and Å*-1* respectively. This ansatz is sufficient in taking into account of the slight anisotropy in the NLs (Fig. S4).

SIV.3 Co2MnGa

SIV.3.1 Tight-binding model

We take the tight-binding model of Co2MnGa from Ref. Chang et al. (2017), which we copy here for ease of reading. It consists of six bands: Three $d$-orbitals from Mn and three $p$-orbitals from Ga. In reciprocal space with the basis $(d_{xz},d_{yz},d_{xy},p_{x},p_{y},p_{z})$, it is given by:

[TABLE]

with:

[TABLE]

SIV.3.2 Velocity field in the nodal Co2MnGa network

The velocity field of Co2MnGa  as illustrated in Fig. S5 here, is the basis of detailed analysis of the shape of its response current. We consider the current contribution of the green and yellow nodal structures in the top panel of Fig.5(a) of the main text. In Fig.5(b)(A) we mark several points for the drift $p_{0}$, and study the resulting drifted Fermi surface, as well as the change of the velocity occupied by it. The sharp increase from (i) to (ii) is a consequence of the rapidly varying velocity field, as soon as the Fermi surface is driven even slightly away from equilibrium. This sharp increase soon saturates and begins to fall as the flower petal is shifted away from the intense red patches situated near the touching points of the inner and outer chains [(ii) to (iv)]. Going from (iv) to (v), the current eventually changes sign because the intense blue patches are now being sampled, The current slowly becomes positive again [(v) to (vi)] when the flower petal is shifted to the outer chain [Fig. 6(b)].

In contrast to the current due to the green nodal structure [Fig.5(b)(A)], the gradient of the current is not as steep for the yellow structure [Fig.5(c)(A)]. Furthermore, the yellow circular ring also contribute exclusively negative velocity [Fig.5(c)(i-vi)] within the drift considered here, unlike the green one where the regions occupied by the ring vary in sign [Fig.5(b)(i-vi)]. These observations, together with the evolution of the velocity field similar to that discussed in the previous paragraph, can account for the current pattern of the yellow nodal structure.

The above analysis demonstrates how the strong inhomogeneity of the velocity field, which results from the knotted links and touching loops, gives rise to its characteristic response. Indeed, on the two high symmetry planes displayed in Fig. S5(d-e), the flower petals are clearly discernible as contours of sign-changing points. As shown in Fig. S5(b), the Fermi surface of Co2MnGa  fits matchingly the void of the velocity field. By visualizing the velocity fields of simpler artificial nodal links (not shown), we observe that the more complicated a nodal network, the more often the velocity field switches direction. Our analysis holds best for small chemical potentials that results in thin NLs, even though the response curve usually remains qualitatively preserved as the chemical potential is increased (Fig. S6), unless the NLs become so thick that they fuse and generate different nodal topologies.

SIV.4 Individual plots of the Fourier harmonics from HHG

In Fig. S7, we present more detailed data of the individual higher harmonics generated of our selected nodal materials in increasing level of complexity.

SIV.5 Comments on our theoretical approach

Here, we comment on our approach and connect it to a few experimental breakthroughs on HHG in solids. Throughout this work, we compute the electrical current by solving the semi-classical equations of motion (EOMs) to obtain the electronic velocity, which is then integrated over the phase space, weighted by the distribution function dictated by the Boltzmann equation.

From explaining the Wiedemann-Franz law to predicting the Bloch oscillation, such an approach has been hugely successful since the dawn of the quantum theory of solids Ziman (2001). The semi-classical EOMs with “anomalous velocity” corrections are not only conceptually appealing but also predictively powerful in explaining diverse phenomena such as the various Hall effects and quantum oscillations Xiao et al. (2010). As for the Boltzmann equation, it is expected to work best for weakly interacting systems which can be well described by the Landau theory of Fermi liquid Kubo (1973). While the importance of correlation effects in nodal systems is certainly material-dependent, the fact that DFT+ARPES calculations of Co2MnGa match so nicely the experimental results Belopolski et al. (2017) may hint that a tight binding description of quasiparticles, and hence the application of Boltzmann equation, should be a good starting point.

In general, solving the Boltzmann equation is highly nontrivial, complicated by the collision integral that assumes different forms and requires various approximation schemes depending on the scattering mechanisms. Confronted with this difficulty, progress can be made by working within the relaxation time approximation (RTA), where the collision integral is assumed to take the form $(F-f)/\tau$. This approximation is motivated by the consideration that the various types of scattering (electron-phonon, electron-impurity, electron-electron) will relax the system to a distribution $F$ (often chosen to be the Fermi-Dirac distribution) within a timescale $\tau$ which can be computed if a microscopic model for the scattering processes is given Ho et al. (2018).

In our work, the collision integral is completely neglected. One way to justify this is as follows. Invoking the constant RTA, one has $(\partial_{t}+\dot{\bm{p}}(t)\cdot\partial_{\bm{p}})f=(F-f)/\tau$, where $f$ is the distribution function, $\bm{p}$ the momentum. In presence of a periodic electric field of frequency $\Omega$ [which enters the equation via $\dot{\bm{p}}(t)$], it is reasonable to assume a periodic steady state, provided the driving force is not too weak and not too slow compared to the dissipative relaxing term. Then, a periodic average of the equation implies that the LHS is of the order $O(\Omega\tau)$, and the RHS of order $O(1)$, so that for $\Omega\tau\gg 1$, the collision term may be ignored altogether. This conforms to the intuition that scattering is not so relevant in presence of a fast driving field Liu et al. (2017).

As with most physical sciences, theoretical descriptions of HHG in solids admit a varying degree of sophistication. For example, one may wish to examine the spatial distribution of the harmonics by taking a Wannier representation for the valence bands Osika et al. (2017), or to study the contributions of interband current by adding an interband polarization term Vampa et al. (2015), or to account for the Coulomb interaction via a dephasing term Golde et al. (2006). For the purpose of illustrating the HHG due to topologically enforced non-linearities in linked NLSMs, we considered only the intraband current within a non-interacting picture. This approach closely mirrors the theoretical model developed in the groundbreaking observation of HHG in bulk ZnO Ghimire et al. (2011). There, a Peierl’s substitution of a continuous wave (cw) electric field is performed on a single band of ZnO. The current is computed by multiplying the velocity with the electron charge and density. Simply by taking into account longer-range hoppings, qualitative features including the linear energy cutoff of the experimental HHG can already be captured. In our current expression, $J[\mathbb{A}(t)]=e\int_{\mathrm{BZ}}\frac{\mathrm{d}^{3}\bm{p}}{(2\pi\hbar)^{3}}f\left[\varepsilon(\bm{p}-\bm{p}(t))\right]\bm{v}(\bm{p})$, upon a coordination transformation, the electric field $\bm{p}(t)$ dependence turns up in the velocity, and the expression becomes similar to that in Ghimire et al. (2011), with the exception that we integrate over the BZ instead of evaluating it at a certain momentum. We note that similar semi-classical approaches have subsequently been employed in Schubert et al. (2014) to discuss Bloch oscillations in HHG, as well as in Liu et al. (2017) to study the role of Berry curvature in HHG of monolayer MoS2. With regards to scattering and many-body effects, rather reassuringly, they should play only minor roles in HHG experiments Liu et al. (2017); Ndabashimiye et al. (2016), in accordance with the intuition that most scattering cannot keep up with the fast driving field. Together with the fact that interactions are mostly insignificant in symmetry-protected topological materials, our theoretical approach may prove sufficient to describe HHG in nodal-line semimetals.