Efficient Control of High Harmonic Generation in Carbon Nanotubes using the Aharonov-Bohm Effect

TL;DR

This paper demonstrates that applying a magnetic field to carbon nanotubes can control high-harmonic generation by inducing a bandgap via the Aharonov-Bohm effect, enabling tunable THz emission.

Contribution

It introduces a novel method to manipulate HHG in CNTs through magnetic field-induced bandgap engineering using the Aharonov-Bohm effect.

Findings

Magnetic field induces a bandgap in metallic CNTs, turning them into gapped Dirac systems.

Gapped CNTs show significantly increased HHG intensity in the THz regime.

Gapless CNTs do not produce HHG due to constant electron velocity.

Abstract

We show that high-harmonic generation (HHG) in carbon nanotubes (CNTs) can be efficiently controlled using the Aharanov-Bohm (AB) effect. When a static magnetic field (B) is applied along the tube, electronic wave functions acquire complex phases along the circumferential direction (AB effect), which modifies the band structure. In particular, when the magnetic field is applied to metallic CNTs, which can be regarded as one-dimensional massless Dirac systems, realistic values of B lead to a nonzero gap in the THz regime. We demonstrate that such change from gapless to gapped Dirac systems drastically increases the HHG intensity in the THz regime. In the gapless Dirac system, the velocity of each electron never changes under the electric field, and thus there is no HHG. On the other hand, the gap opening activates both the interband and itraband currents, which strongly contribute to HHG. Our work demonstrates a unique way to manipulate HHG in nanotubes by tuning electronic wave functions using the magnetic field and the tube structure.

Equations

$$\hat{H}(t)=-t_{\mathrm{hop}}\sum _{⟨ij⟩}{e}^{i\frac{q}{\hslash }\mathbf{A}(t)\cdot {\mathbf{r}}_{ij}}{\hat{c}}_i^{†}{\hat{c}}_j-\mu \sum _i{\hat{n}}_i.$$

$$\mathbf{h}(\mathbf{k})=\left[\begin{array}{cc}-\mu & F(\mathbf{k})\\ F^{\ast }(\mathbf{k})& -\mu \end{array}\right],$$

$${\partial }_t{\mathbf{\rho }}_{\mathbf{k}}(t)$$

$$\displaystyle+\Bigl{(}\frac{1}{T_{1}}-\frac{1}{T_{2}}\Bigl{)}{\bf U}({\bm{k}}(t))\bm{\rho}^{\rm H}_{{\rm off},{\bm{k}}}(t){\bf U}^{\dagger}({\bm{k}}(t)),$$

$${\partial }_t{J}_{||}(t)=\frac{1}{n}\frac{t_{\mathrm{hop}}}{\hslash }{\sigma }_{0}E(t)-\frac{1}{T_1}{J}_{||}(t),$$

$${\sigma }_{\mathrm{intra}}^{(3)}\simeq -F\frac{1}{A_B^2}\frac{1}{{\zeta }_3},{\sigma }_{\mathrm{inter}}^{(3)}\simeq F\frac{1}{A_B^2}\frac{1}{{\xi }_3},$$

$$\displaystyle F\equiv i\sigma_{0}t_{\rm hop}\hbar^{2}\frac{1}{n}\frac{1}{\xi_{1}}\Bigl{(}\frac{\lambda_{0}}{2\zeta_{2}}+\frac{\lambda_{2}}{\xi_{2}}\Bigl{)}$$

$$\hat{H}(t)=-t_{\mathrm{hop}}\sum _{⟨ij⟩}{e}^{i\frac{q}{\hslash }\mathbf{A}(t)\cdot {\mathbf{r}}_{ij}}{\hat{c}}_i^{†}{\hat{c}}_j-\mu \sum _i{\hat{n}}_i,$$

$${\mathbf{L}}_0\equiv {\mathbf{a}}_1+{\mathbf{a}}_2,{\mathbf{T}}_0\equiv {\mathbf{a}}_2-{\mathbf{a}}_1,$$

$$\mathbf{k}=\frac{i_L}{n}{\mathbf{l}}_0+\frac{i_T}{N_T}{\mathbf{t}}_0,$$

$$\hat{H}(t)={\sum ^{~}}_{\mathbf{k}}\left[\begin{array}{cccc}{\hat{c}}_{\mathbf{k}{A}_1}^{†}& {\hat{c}}_{\mathbf{k}{A}_2}^{†}& {\hat{c}}_{\mathbf{k}{B}_1}^{†}& {\hat{c}}_{\mathbf{k}{B}_2}^{†}\end{array}\right]\tilde{\mathbf{h}}(\mathbf{k}(t))\left[\begin{array}{c}{\hat{c}}_{\mathbf{k}{A}_1}\\ {\hat{c}}_{\mathbf{k}{A}_2}\\ {\hat{c}}_{\mathbf{k}{B}_1}\\ {\hat{c}}_{\mathbf{k}{B}_2}\end{array}\right],$$

$$\tilde{\mathbf{h}}(\mathbf{k})=\left[\begin{array}{cccc}-\mu & 0& \alpha (\mathbf{k}{)}^{\ast }& \gamma (\mathbf{k}{)}^{\ast }\\ 0& -\mu & \gamma (\mathbf{k}{)}^{\ast }& \alpha (\mathbf{k}{)}^{\ast }\\ \alpha (\mathbf{k})& \gamma (\mathbf{k})& -\mu & 0\\ \gamma (\mathbf{k})& \alpha (\mathbf{k})& 0& -\mu \end{array}\right],$$

$$\alpha (\mathbf{k})=-t_{\mathrm{hop}}[e^{-i{\mathbf{b}}_1\cdot \mathbf{k}}+e^{-i{\mathbf{b}}_3\cdot \mathbf{k}}],\gamma (\mathbf{k})=-t_{\mathrm{hop}}{e}^{-i{\mathbf{b}}_2\cdot \mathbf{k}}.$$

$$\left[\begin{array}{c}{\hat{c}}_{\mathbf{k}{A}_{+}}\\ {\hat{c}}_{\mathbf{k}{B}_{+}}\\ {\hat{c}}_{\mathbf{k}{A}_{-}}\\ {\hat{c}}_{\mathbf{k}{B}_{-}}\end{array}\right]\equiv \frac{1}{\sqrt{2}}\left[\begin{array}{cccc}1& 1& 0& 0\\ 0& 0& 1& 1\\ e^{i\frac{\pi }{6}}& -e^{i\frac{\pi }{6}}& 0& 0\\ 0& 0& e^{-i\frac{\pi }{6}}& -e^{-i\frac{\pi }{6}}\end{array}\right]\left[\begin{array}{c}{\hat{c}}_{\mathbf{k}{A}_1}\\ {\hat{c}}_{\mathbf{k}{A}_2}\\ {\hat{c}}_{\mathbf{k}{B}_1}\\ {\hat{c}}_{\mathbf{k}{B}_2}\end{array}\right].$$

$$\hat{H}(t)=$$

$${\sum ^{~}}_{\mathbf{k}}\left[\begin{array}{cccc}{\hat{c}}_{\mathbf{k}{A}_{+}}^{†}& {\hat{c}}_{\mathbf{k}{B}_{+}}^{†}& {\hat{c}}_{\mathbf{k}{A}_{-}}^{†}& {\hat{c}}_{\mathbf{k}{B}_{-}}^{†}\end{array}\right]{\mathbf{h}}^{\prime }(\mathbf{k}(t))\left[\begin{array}{c}{\hat{c}}_{\mathbf{k}{A}_{+}}\\ {\hat{c}}_{\mathbf{k}{B}_{+}}\\ {\hat{c}}_{\mathbf{k}{A}_{-}}\\ {\hat{c}}_{\mathbf{k}{B}_{-}}\end{array}\right]$$

$${\mathbf{h}}^{\prime }(\mathbf{k})=\left[\begin{array}{cccc}-\mu & F(\mathbf{k})& 0& 0\\ F(\mathbf{k}{)}^{\ast }& -\mu & 0& 0\\ 0& 0& -\mu & F(\mathbf{k}+{\mathbf{l}}_0)\\ 0& 0& F(\mathbf{k}+{\mathbf{l}}_0{)}^{\ast }& -\mu \end{array}\right].$$

$$\hat{H}(t)=\sum _{\mathbf{k}}\left[\begin{array}{cc}{\hat{c}}_{\mathbf{k}A}^{†}& {\hat{c}}_{\mathbf{k}B}^{†}\end{array}\right]\mathbf{h}(\mathbf{k}(t))\left[\begin{array}{c}{\hat{c}}_{\mathbf{k}A}\\ {\hat{c}}_{\mathbf{k}B}\end{array}\right]$$

$$\mathbf{h}(\mathbf{k})=\left[\begin{array}{cc}-\mu & F(\mathbf{k})\\ F^{\ast }(\mathbf{k})& -\mu \end{array}\right].$$

$${\hat{J}}_{\parallel }(t)$$

$$\displaystyle\sum_{{\bm{k}}}\begin{bmatrix}\hat{c}^{\dagger}_{{\bm{k}}A}&\hat{c}^{\dagger}_{{\bm{k}}B}\end{bmatrix}\begin{bmatrix}0&-e^{i{\bm{k}}(t)\cdot{\bm{b}}_{\alpha}}\\ e^{-i{\bm{k}}(t)\cdot{\bm{b}}_{\alpha}}&0\end{bmatrix}\begin{bmatrix}\hat{c}_{{\bm{k}}A}\\ \hat{c}_{{\bm{k}}B}\end{bmatrix}\Bigl{\}},$$

$${\hat{H}}_{\mathbf{k}}^{(\mathbf{K},{\mathbf{K}}^{\prime })}(t)=\left[\begin{array}{cc}{\hat{c}}_{\mathbf{k}A}^{†}& {\hat{c}}_{\mathbf{k}B}^{†}\end{array}\right]{\mathbf{h}}_{\pm }(\mathbf{k}(t))\left[\begin{array}{c}{\hat{c}}_{\mathbf{k}A}\\ {\hat{c}}_{\mathbf{k}B}\end{array}\right]$$

$${\mathbf{h}}_{\pm }(\mathbf{k})=\left[\begin{array}{cc}-\mu & -\frac{\sqrt{3}a{t}_{\mathrm{hop}}}{2}[\pm k_x+i{k}_y]\\ -\frac{\sqrt{3}a{t}_{\mathrm{hop}}}{2}[\pm k_x-i{k}_y]& -\mu \end{array}\right].$$

$${\hat{J}}_{\parallel ,\mathbf{k}}^{(\mathbf{K},{\mathbf{K}}^{\prime })}=\mp \frac{\sqrt{3}aq{t}_{\mathrm{hop}}}{2\hslash }({\mathbf{e}}_{\parallel }\cdot {\mathbf{e}}_x)\left[\begin{array}{cc}{\hat{c}}_{\mathbf{k}A}^{†}& {\hat{c}}_{\mathbf{k}B}^{†}\end{array}\right]\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\left[\begin{array}{c}{\hat{c}}_{\mathbf{k}A}\\ {\hat{c}}_{\mathbf{k}B}\end{array}\right]$$

$$-\frac{\sqrt{3}aq{t}_{\mathrm{hop}}}{2\hslash }({\mathbf{e}}_{\parallel }\cdot {\mathbf{e}}_y)\left[\begin{array}{cc}{\hat{c}}_{\mathbf{k}A}^{†}& {\hat{c}}_{\mathbf{k}B}^{†}\end{array}\right]\left[\begin{array}{cc}0& i\\ -i& 0\end{array}\right]\left[\begin{array}{c}{\hat{c}}_{\mathbf{k}A}\\ {\hat{c}}_{\mathbf{k}B}\end{array}\right].$$

$$\hslash {\partial }_t{\mathbf{\rho }}_{\mathbf{k}}(t)$$

$$[{\partial }_t{\mathbf{\rho }}_{\mathbf{k}}(t){]}_{\mathrm{corr}}=-\frac{{\mathbf{\rho }}_{\mathbf{k}}(t)-{\mathbf{\rho }}_{\mathrm{eq},\mathbf{k}(t)}}{T_1}$$

$$\displaystyle\;\;\;\;\;+\Bigl{(}\frac{1}{T_{1}}-\frac{1}{T_{2}}\Bigl{)}{\bf U}({\bm{k}}(t))\bm{\rho}^{\rm H}_{{\rm off},{\bm{k}}}(t){\bf U}^{\dagger}({\bm{k}}(t)).$$

$$\hslash {\partial }_t{\mathbf{\rho }}_{\mathbf{k}}^{\mathrm{H}}(t)$$

$$+\hslash [{\partial }_t{\mathbf{\rho }}_{\mathbf{k}}^{\mathrm{H}}(t){]}_{\mathrm{corr}},$$

$$[{\partial }_t{\mathbf{\rho }}_{\mathbf{k}}^{\mathrm{H}}(t){]}_{\mathrm{corr}}$$

$$J_{\parallel ,\mathrm{intra}}(t)$$

$$J_{\parallel ,\mathrm{inter}}(t)$$

$$\displaystyle P_{\rm inter}=\frac{q}{\hbar N}\sum_{{\bm{k}},a}e_{\parallel,a}{\rm Tr}\Bigl{[}\bm{\rho}^{\rm H}_{\bm{k}}(t){\bf d}^{\rm off}_{a}({\bm{k}}(t))\Bigl{]},$$

Full text

The recent development of strong lasers in the terahertz (THz) and mid-infrared regimes enables us to study various nonlinear and nonequilibrium phenomena in condensed matter Kruchinin et al. (2018); Giannetti et al. (2016); Basov et al. (2017); Cavalleri (2018); Oka and Kitamura (2019); de la Torre et al. (2021); Koshihara et al. (2022). Among them, high-harmonic generation (HHG), where higher harmonics of the laser injected to materials are radiated, is a fundamental and technologically important example Corkum and Krausz (2007); Krausz and Ivanov (2009); Ghimire and Reis (2019). HHG was originally studied in gas systems Ferray et al. (1988); Corkum (1993); Lewenstein et al. (1994). Recently, HHG in condensed matter has been observed, and the scope of the HHG research is widely extended to conventional semiconductors, band insulators Ghimire et al. (2011); Schubert et al. (2014); Vampa et al. (2014); Luu et al. (2015); Vampa et al. (2015); Langer et al. (2016); Hohenleutner et al. (2015); Ndabashimiye et al. (2016); Otobe (2016); Tancogne-Dejean et al. (2017); Ikemachi et al. (2017); Liu et al. (2017); You et al. (2017); Kaneshima et al. (2018); Sekiguchi et al. (2022), topological materials Giorgianni et al. (2016); Yoshikawa et al. (2017); Hafez et al. (2018); Silva et al. (2019); Chacón et al. (2020); Cheng et al. (2020); Schmid et al. (2021); Baykusheva et al. (2021a) and strongly correlated systems Silva et al. (2018); Murakami et al. (2018); Tancogne-Dejean et al. (2018); Imai et al. (2020); Murakami et al. (2021); Orthodoxou et al. (2021); Udono et al. (2022); Bionta et al. (2021); Shao et al. (2022); Hansen et al. (2022); Uchida et al. (2022); Murakami et al. (2022); Hansen and Madsen (2022); Grånäs et al. (2022).

Condensed matter is considered as an interesting playground of the HHG research, at least, for two good reasons. Firstly, the periodic potential from atoms, which is absent in gaseous systems, can lead to peculiar motion of electrons under strong electric field. An interesting example is the Dirac system, which is characterized by the non-parabolic band dispersion. Due to the nonlinear dynamics of electrons in the vicinity of the Dirac node, such systems are theoretically expected to exhibit large nonlinearity Slepyan et al. (1999); Mikhailov (2007); Wright et al. (2009); Al-Naib et al. (2014, 2015); Cheng et al. (2015); Chizhova et al. (2017); Mrudul and Dixit (2021); Sato et al. (2021); Baykusheva et al. (2021b). HHG was experimentally studied for two-dimensional (2D) Dirac systems realized in graphene Yoshikawa et al. (2017); Hafez et al. (2018) and in the surface of topological insulators Giorgianni et al. (2016); Schmid et al. (2021); Baykusheva et al. (2021a), and for the 3D Dirac system in Cd3As2 Cheng et al. (2020); Kovalev et al. (2020). It was reported that the efficiency of HHG is indeed large compared to other systems Hafez et al. (2018); Cheng et al. (2020); Schmid et al. (2021), although the origin of HHG is not necessary processes theoretically expected Paul et al. (2013); Hafez et al. (2018). The second reason is the controllability of HHG. Generally, properties of condensed matter can be tuned by system parameters such as temperature and doping level, which can be also applicable to HHG. For example, peculiar HHG dependence on crystal orientations Ghimire et al. (2011); You et al. (2017), electronic phases Silva et al. (2018); Uchida et al. (2022); Murakami et al. (2022); Bionta et al. (2021); Shao et al. (2022), physical pressure Tamaya and Kato (2021) and doping Nishidome et al. (2020); de Vega et al. (2020); Hansen and Madsen (2022) has been proposed theoretically or been observed experimentally.

In this work, we propose a unique way to efficiently control HHG in carbon nanotubes (CNTs) combining physics of the Dirac system and the Aharonov-Bohm (AB) effect, see Fig. 1(a). The proposed setup provides an ideal platform to study HHG in Dirac systems, systematically. When the static magnetic field ($B$) is applied along the tube, electrons in the CNT gain complex phases along the circumferential direction, which modifies the band structure (the AB effect) Ajiki and Ando (1993, 1994); Zaric et al. (2004); Matsunaga et al. (2008). In particular, metallic CNTs, which can be regarded as 1D massless Dirac systems, become massive Dirac systems when $B\neq 0$. This change qualitatively modifies the electron dynamics, and drastically enhances HHG.

In the dipolar gauge Li et al. (2020); Schüler et al. (2021); Murakami and Schüler (2022), the tight-binding Hamiltonian for the CNT with the magnetic field and the AC electric field ($E(t)$) along the tube is expressed as

[TABLE]

$\hat{c}^{\dagger}_{i}$ is the creation operator of an electron at the $i$th site and we neglect the spin degrees of freedom. $\langle ij\rangle$ indicates the nearest neighbor sites, $t_{\rm hop}$ is the hopping parameter, $q$ is the charge of an electron, $\hbar$ is the Dirac constant, $\mu$ is the chemical potential and $\hat{n}_{i}=\hat{c}^{\dagger}_{i}\hat{c}_{i}$. The model is defined on the graphite sheet imposing the periodic boundary condition with respect to the chiral vector ${\bm{L}}$, which is characterized by the chiral index $(n,m)$ as ${\bm{L}}=n{\bm{a}}_{1}+m{\bm{a}}_{2}$, see Fig. 1(b) Hamada et al. (1992). The effects of the magnetic and electric fields are included via the Peierls phase through the vector potential ${\bm{A}}(t)={\bm{A}}_{B}+{\bm{A}}_{E}(t)$. The contribution from the static magnetic field is ${\bm{A}}_{B}=A_{B}{\bm{e}}_{\perp}$ and that from the electric field is ${\bm{A}}_{E}(t)=A_{E}(t){\bm{e}}_{\parallel}$. $A_{B}=\frac{B\cdot|{\bm{L}}|}{4\pi}$ and $E(t)=-\partial_{t}A_{E}(t)$. ${\bm{e}}_{\parallel}$ and ${\bm{e}}_{\perp}$ are unit vectors of the tube and circumferential directions, respectively, see Fig. 1(b). In the following, we focus on the armchair CNT with the chiral index $(n,m)$ with $n=m$.

The Fourier transform of the Hamiltonian (6) is $\hat{H}(t)=\sum_{\bm{k}}\hat{\bm{\psi}}^{\dagger}_{\bm{k}}{\bf h}({\bm{k}}(t))\hat{\bm{\psi}}_{\bm{k}}$, where

[TABLE]

and ${\bm{k}}(t)\equiv{\bm{k}}-\frac{q}{\hbar}{\bm{A}}(t)$, see the Supplementary Material (SM) 111See Supplemental Material at [URL will be inserted by publisher] for the details of the model and numerical method; the derivation of Eqs. (4) and  (5); the doping effects; the dependence of the excitation frequency; the dependence on the chiral index; the effects of the relaxation times. The Supplemental Material also contains Refs. Khosravi et al. (2009); Dimitrovski et al. (2017); Dvuzhilova et al. (2021); Wilhelm et al. (2021). for details. Here, $\hat{\bm{\psi}}^{\dagger}_{\bm{k}}=[\hat{c}^{\dagger}_{{\bm{k}}A},\hat{c}^{\dagger}_{{\bm{k}}B}]$, and $\hat{c}^{\dagger}_{{\bm{k}}A}$ and $\hat{c}^{\dagger}_{{\bm{k}}B}$ are the Fourier transform of $\hat{c}^{\dagger}_{i}$ on the A and B sublattices, respectively. $F({\bm{k}})=-t_{\rm hop}[e^{i{\bm{k}}\cdot{\bm{b}}_{1}}+e^{i{\bm{k}}\cdot{\bm{b}}_{2}}+e^{i{\bm{k}}\cdot{\bm{b}}_{3}}]$, see Fig. 1(b) for ${\bm{b}}_{i}$. The Hamiltonian (16) is the same as that of graphene, and the set of ${\bm{k}}_{B}(\equiv{\bm{k}}-\frac{q}{\hbar}{\bm{A}}_{B})$ for the CNT corresponds to some lines along ${\bm{k}}_{\perp}$ on the BZ of graphene, see Fig. 1(c). Note that ${\bm{k}}(t)={\bm{k}}_{B}-\frac{q}{\hbar}{\bm{A}}_{E}(t)$, which moves along ${\bm{k}}_{\parallel}$ in time. For $B=0$, ${\bm{k}}_{B}$’s for the CNT pass the Dirac points $K$ and $K^{\prime}$ of graphene, hence the CNT bands host gapless Dirac bands, see Figs. 1(c) and (d). Meanwhile, for $B\neq 0$, ${\bm{k}}_{B}$’s are shifted by $-\frac{q}{\hbar}{\bm{A}}_{B}$ and do not passes the Dirac points. Hence the CNT bands host gapful Dirac bands. Around the Dirac points, the Hamiltonian (16) is effectively described by the Dirac model, see Fig. 1 (d) and SM Note (1) for details. In our setup, the dynamics of electrons around the Dirac points leads to HHG, which is captured by this model.

Diameters of single-wall CNTs are typically $0.7\sim 2$nm Yanagi et al. (2008). The static magnetic field applied to CNTs in optical experiments can be several tens of Tesla Zaric et al. (2004). In the following, we set $t_{\rm hop}=3$eV and $a=0.246$nm and focus on the CNT with $(n,m)=(15,15)$, whose diameter is $2.04$nm. The band gap for $B=30$T becomes $29.7$meV ($\simeq 7.5$THz), see Fig. 1 (d). We mainly focus on half filling ($\mu=0$) and apply the gaussian electric field described by $A_{E}(t)=\frac{E_{0}}{\hbar\Omega}\exp[-\frac{(t-t_{0})^{2}}{2\sigma^{2}}]\sin(\Omega(t-t_{0}))$. We set the electric-field frequency $\hbar\Omega=10$meV ($\simeq 2.5$THz), the pulse center $t_{0}=3.95$ps and the pulse width $\sigma={658}$fs, in the main text.

We simulate the time evolution of the tight-binding and Dirac models using the semiconductor Bloch equation (SBE)  Huttner et al. (2017); Yue and Gaarde (2022)

[TABLE]

where $\bm{\rho}_{\bm{k}}(t)$ is the single-particle density matrix consisting of $\rho_{\alpha\beta,{\bm{k}}}(t)=\langle\hat{c}^{\dagger}_{{\bm{k}}\beta}(t)\hat{c}_{{\bm{k}}\alpha}(t)\rangle$ with $\alpha,\beta=A,B$. The last two terms express the relaxation processes within the relaxation-time approximation. $\bm{\rho}_{{\rm eq,}{\bm{k}}}$ is the equilibrium density matrix for the initial temperature ($T$). ${\bf U}({\bm{k}})$ is a unitary matrix diagonalizing ${\bf h}({\bm{k}})$ as ${\bf U}^{\dagger}({\bm{k}}){\bf h}({\bm{k}}){\bf U}({\bm{k}})={\bm{\epsilon}}({\bm{k}})$ with ${\bm{\epsilon}}({\bm{k}})={\rm diag}[\epsilon_{0}({\bm{k}}),\epsilon_{1}({\bm{k}})]$. It offers the descrition in the band (Houston) basis. In this basis, the density matrix becomes $\bm{\rho}^{\rm H}_{{\bm{k}}}(t)\equiv{\bf U}^{\dagger}({\bm{k}}(t))\bm{\rho}_{\bm{k}}(t){\bf U}({\bm{k}}(t))$, whose off-diagonal components are $\bm{\rho}^{\rm H}_{{\rm off},{\bm{k}}}(t)$ Murakami and Schüler (2022). The diagonal components of $\bm{\rho}^{\rm H}_{\bm{k}}$ indicate the occupation of the bands, and the off-diagonal components indicate the interband hybridization. Thus, $T_{1}$ represents the relaxation time of the band occupation, while $T_{2}$ represents the dephasing time of the interband coherence. We set $T_{1}=98.8$fs and $T_{2}=19.8$fs to be consistent with previous studies for graphene and CNTs  Habenicht et al. (2007); George et al. (2008); Tani et al. (2012); Heide et al. (2021); Sato et al. (2021).

The light radiation originates from the oscillation of electric polarization (or corresponding current) Jackson (1998). Since the current along the circumferential direction does not yield polarization, the HHG spectrum (normalized by the number of sites) is evaluated from $J_{\parallel}(t)$ (the current along ${\bm{e}}_{\parallel}$) as $I_{\rm HHG}(\omega)=|\omega J_{\parallel}(\omega)|^{2}$. Here, $J_{\parallel}(\omega)=\int dte^{i\omega t}J_{\parallel}(t)$, and $J_{\parallel}(t)$ is evaluated from $\rho_{\alpha\beta,{\bm{k}}}(t)$ using the current operator $\hat{J}_{\parallel}(t)\equiv-\frac{1}{N}\frac{\delta\hat{H}(t)}{\delta A_{E}(t)}$. $2N$ is the number of atoms (sites) for the CNT. In the conditions studied here, HHG mainly originates from the bands closest to the Dirac points, i.e. solid lines in Fig. 1(d).

Figure 2(a) shows $I_{\rm HHG}$ of the armchair CNT with and without the magnetic field obtained from the tight-binding model. Only odd harmonics components appear in $I_{\rm HHG}(\omega)$, because the system (6) is symmetric under a time-space operation, $t\rightarrow t+\frac{\pi}{\Omega}$ and ${\bf e}_{\parallel}\rightarrow-{\bf e}_{\parallel}$ during the time-periodic excitation, see Ref. Neufeld et al. (2019) for general argument. When $B=0$ the intensity for $l(\geq 3)$th harmonics is very small. When $B\neq 0$ the HHG intensity for $l\geq 3$ is strongly enhanced by more than several orders of magnitude. Figure 2(b) shows the intensity at the HHG peaks $I_{\rm peak}(l)$ as a function of $B$ for $l=3,5,7$. $I_{\rm peak}(l)$ is defined as the maximum value of $I_{\rm HHG}(\omega)$ for $\omega\in[l\Omega-\Omega,l\Omega+\Omega]$. $I_{\rm peak}(l)$ is non-monotonic against $B$, since a gap too large compared to $\Omega$ suppresses the excitation and thus HHG, see the inset of Fig. 2(a). Namely, there exists the optimal value of $B$ maximizing $I_{\rm HHG}(\omega)$ for a given $E_{0}$. Furthermore, the HHG intensity in the CNT per atom for this optimal $B$ is much larger than that in graphene per atom, see Fig. 2(a). This is understood as follows. In the CNT, we can tune $B$ such that the bands closest to the Dirac points yield HHG most efficiently for a given $E_{0}$, and the contributions from the rest bands are negligible. Meanwhile, in graphene, which is identical to a CNT with $n,m\rightarrow\infty$, the number of bands along ${\bm{k}}_{\parallel}$ contributing to HHG increases. However, the contribution from each band are not optimal for a given electric-field condition, and they may interfere and cancel each other on average. We also note that the drastic enhancement of HHG with $B$ can be observed even for doped CNTs or for higher-frequency excitations, see SM Note (1).

The efficiency of the enhancement depends on the chiral index $n(=m)$. Firstly, the momentum shift due to the AB effect is $\frac{q}{\hbar}A_{B}$, which is proportional to $B$ and $L(\equiv|{\bm{L}}|\propto n)$. Therefore, in order to realize a given value of $A_{B}$ (or the gap), smaller $B$ is required for larger $n$, see SM Note (1) for examples. Thus, as far as only the bands closest to the Dirac points are involved in HHG, larger $n$ is favorable for the efficient control of HHG. Meanwhile, with increasing $n$, the number of bands along ${\bm{k}}_{\parallel}$ increases, and many bands start to contribute to HHG. Then, the impact of the AB effect on HHG becomes less prominent. In the extreme limit of $n\rightarrow\infty$ (graphene), the AB effect has no impact on physical properties.

To understand the origin of the enhancement of HHG with the magnetic field, we compare the HHG intensity evaluated from the tight-binding model ($I_{\rm HHG}$) and that from the Dirac model ($I_{\rm HHG}^{\rm Dirac}$), see Fig. 3(a)(b). We only consider the contributions from the bands closest to $K$ and $K^{\prime}$ for the Dirac model. The results from the Dirac model match well those from the tight-binding model both for $B=0$ and $B=30$T. The only exception is the 3rd harmonic generation at $B=0$, which originates from the nonlinear part of the band dispersion. Namely, the 1D gapless Dirac model shows no HHG at $B=0$. Simply speaking, this is because the transition of electrons between bands in the 1D gapless Dirac system under electric field is always perfect and electrons never change the velocity, see Fig. 3(c). Without the relaxation, this leads to $J_{||}(t)\propto A_{E}(t)$ Ishikawa (2010); Óscar Zurrón-Cifuentes et al. (2020). If we take into account the relaxation, we obtain

[TABLE]

from the SBE (21), where $\sigma_{0}=\frac{\sqrt{3}a^{2}q^{2}}{2\pi\hbar}$, see SM Note (1). This is a linear differential equation, leading to no HHG. The last term of Eq.(4), which we call $[\partial_{t}J_{\parallel}(t)]_{\rm corr}$, originates from the relaxation term of the SBE. The present form originates from the relaxation-time approximation in Eq. (21). In principle, if $[\partial_{t}J_{\parallel}(t)]_{\rm corr}$ depends on $J_{\parallel}$ nonlinearly, the 1D gapless Dirac system can show HHG. Still, in a previous experiment of HHG in CNTs Nishidome et al. (2020), metallic CNTs show very weak HHG, which is consistent with the consequence from Eq. (4).

To obtain further insight, we evaluate the contribution from the intraband current ($I_{\rm intra}=|\omega J_{\parallel,\rm intra}(\omega)|^{2}$) and that from the interband current ($I_{\rm inter}=|\omega J_{\parallel,\rm inter}(\omega)|^{2}$) using the Dirac model. We define the intraband current as $J_{\parallel,\rm intra}(t)\equiv\frac{q}{\hbar N}\sum_{\bm{k}}{\rm Tr}[\bm{\rho}^{H}_{\bm{k}}(t)\partial_{k_{\parallel}}{\bm{\epsilon}}({\bm{k}}(t))]$ and the interband current as $J_{\parallel,\rm inter}(t)\equiv J_{\parallel}(t)-J_{\parallel,\rm intra}(t)$. The former is associated with the change in the band occupation, while the latter is associated with that of the interband polarization (the off-diagonal component of $\bm{\rho}^{H}_{\bm{k}}(t)$), see SM Note (1) for details. When $B=0$, the perfect transition between the conduction and valence bands happens without inducing the interband polarization, leading to $I_{\rm HHG}=I_{\rm intra}$, see Fig. 3(c). Meanwhile, when $B\neq 0$, the gap opens and the dipole moments between the bands become finite, which allows electrons to stay in the same band and tunnel to the other band, see Fig. 3(c). Thus, both the intraband and interband currents are activated to contribute to HHG, see Fig. 3(b). In the present case, the gap and the excitation frequency are comparable, which makes the contributions from these currents comparable. The HHG contributions from $J_{\parallel,\rm intra}(t)$ and $J_{\parallel,\rm inter}(t)$ cancel each other due to their interference.

Finally, we discuss the dependence on the strength of the applied AC electric field $E_{0}$, see Fig. 1(a). In Fig. 4, we show the intensity at the HHG peaks $I_{\rm peak}(l)$ as a function of $E_{0}$ for $l=3,5,7$. In the perturbative regime with respect to $E_{0}$, we expect $I_{\rm peak}(l)\propto(E_{0})^{2l}$. (For small $B$, this perturbative behavior is recovered with $E_{0}$ below the range shown in Fig. 4.) For all harmonics, the deviation from the perturbative regime happens at larger $E_{0}$ for larger $B$, i.e. larger gap. $I_{\rm peak}(l)$ tends to be saturated due to the higher order corrections. Because of this, the value of $B$ maximizing $I_{\rm peak}(l)$ for a given $E_{0}$ depends on the value of $E_{0}$. Namely, when $E_{0}$ is small, the optimal $B$ is small, while for larger $E_{0}$, the optimal $B$ becomes large. We note that these properties are insensitive to the choice of $T_{1}$ and $T_{2}$, see SM Note (1).

One interesting observation is that the optimal value of $B$ continuously decreases with decreasing $E_{0}$ and it becomes a few Tesla for $E_{0}\sim 1$ kV/cm. This indicates that $\sigma^{(l)}$ monotonically increases with decreasing $B$. Here, $\sigma^{(l)}$ represents the $l$th order optical conductivity, i.e. $J_{\parallel}(l\Omega)\simeq\sigma^{(l)}(E_{0})^{l}$. However, this peculiar behavior turns out to be sensitive to the choice of $T_{1}$ and $T_{2}$, see SM Note (1). When $T_{1}=T_{2}$, the optimal $B$ for small $E_{0}$ remains around $20$T, where the gap is almost twice of $\Omega$. The crucial dependence on $T_{1}$ and $T_{2}$ can be attributed to the change in cancelation between the interaband and intraband currents. This point can be directly seen in the analytic expression of $\sigma^{(3)}$. If we express contributions from $J_{\rm inter}$ and $J_{\rm intra}$ as $\sigma^{(3)}_{\rm inter}$ and $\sigma^{(3)}_{\rm intra}$, the leading terms in the limit of $B\rightarrow 0$ read

[TABLE]

at $T=0$. Here $\zeta_{m}=m\hbar\Omega+i\hbar/T_{1}$ and $\xi_{m}=m\hbar\Omega+i\hbar/T_{2}$ and $\lambda_{m}=\int dx\frac{x^{m}}{\sqrt{1+x^{2}}^{7}}$. The expression tells that $\sigma^{(3)}$ diverges for $B\rightarrow 0$ when $T_{1}\neq T_{2}$, while the divergence is absent when $T_{1}=T_{2}$ due to the cancelation between $\sigma^{(3)}_{\rm inter}$ and $\sigma^{(3)}_{\rm intra}$. Our numerical simulation suggests that the same happens also for higher harmonics with $l\geq 5$. Namely, for $T_{1}=T_{2}$, the contributions from the interband and intraband currents cancel each other efficiently, while for $T_{1}\neq T_{2}$, the cancelation becomes weaker, see SM Note (1). Although we used $T_{1}=98.8$fs and $T_{2}=19.8$fs to be consistent with previous studies, to fully determine the realistic choice of $T_{1}$ and $T_{2}$ for the present energy scale, microscopic evaluation of relaxation and dephasing processes is required Malic et al. (2011); Kemper et al. (2013); Stefanucci and Leeuwen (2013); Schüler et al. (2020); Ridley et al. (2022); Rostami and Cappelluti (2021).

To summarize, we showed that HHG in metallic CNTs can be drastically manipulated in the THz regime using the AB effect. Tuning the strength of the magnetic field, the HHG intensity can be maximized, where HHG per atom is more efficient than in graphene. Such controllability of HHG may be useful to develop an efficient HHG device in the THz regime. In addition, with realistic parameters of CNTs and external fields, a metallic CNT under magnetic fields can be regarded as a tunable 1D Dirac system for THz excitations. Therefore, it serves as a useful platform to systematically study HHG in Dirac systems. Furthermore, the unique idea of manipulating HHG via the AB effect can be applied for all kinds of nanotubes including semiconducting CNTs, boron nitride nanotubes Rubio et al. (1994); Chopra et al. (1995) and transition-metal dichalcogenide nanotubes Tenne et al. (1992).

I Detail of numerical simulation

I.1 Tight Binding Model

As explained in the main text, the tight-binding (TB) model for the armchair CNT with the static magnetic field and the AC electric field along the tube can be constructed on a honeycomb (graphite) sheet;

[TABLE]

where the periodic boundary condition (PBC) with respect to the chiral vector ${\bm{L}}(=n{\bm{a}}_{1}+n{\bm{a}}_{2})$ is imposed. We set ${\bm{A}}(t)={\bm{A}}_{B}+{\bm{A}}_{E}(t)$. In this TB model for CNTs, the effects of the curvature of the tube are ignored for simplicity Hamada et al. (1992). Due to the curvature of the tube, the Dirac points may shift. Still, for the armchair CNT, which we focus on in this paper, the shift is along $k_{\parallel}$ and the system remains metallic Hamada et al. (1992). Thus, the physics of the manipulation of HHG discussed in this paper is hardly affected. The primitive cell is defined by

[TABLE]

see Fig. 5. Note that ${\bm{L}}=n{\bm{L}}_{0}\propto{\bm{e}}_{\perp}$ and ${\bf T}_{0}\propto{\bm{e}}_{\parallel}$. We have four sites (atoms) in the primitive cell. Due to the PBC with respect to ${\bm{L}}$, we have $n$ primitive cells along ${\bm{e}}_{\perp}$. We also impose the PBC along ${\bm{e}}_{\parallel}$ and take $N_{T}(\rightarrow\infty)$ primitive cells along ${\bm{e}}_{\parallel}$. Using the reciprocal vectors ${\bm{l}}_{0}\equiv 2\pi\frac{{\bm{L}}_{0}}{|{\bm{L}}_{0}|^{2}}$ and ${\bm{t}}_{0}\equiv 2\pi\frac{{\bf T}_{0}}{|{\bf T}_{0}|^{2}}$, the wave vectors can be defined as

[TABLE]

where $i_{L}=0,1\cdots n-1$ and $i_{T}=0,1\cdots N_{T}-1$. We name the four sites in the primitive cell as $A_{1},B_{1},A_{2},$ and $B_{2}$ as in Fig. 5.

Now, we consider the expression of the Hamiltonian in the momentum (Fourier) space. We introduce the electron operators in the momentum space as $\hat{c}^{\dagger}_{{\bm{k}}\alpha}=\frac{1}{\sqrt{N_{c}}}\sum_{i\in\alpha}e^{i{\bm{k}}\cdot{\bm{r}}_{i}}\hat{c}^{\dagger}_{i}$ for $\alpha=A_{1},A_{2},B_{1},B_{2}$. Here ${\bf r}_{i}$ is the position vector of the $i$th site, $N_{c}=n\cdot N_{T}$ is the total number of the cells. The inverse tranform is $\hat{c}^{\dagger}_{\bm{r}_{i}}=\frac{1}{\sqrt{N_{c}}}\sum_{\bm{k}}e^{-i{\bm{k}}\cdot{\bf r}_{i}}\hat{c}^{\dagger}_{{\bm{k}}\alpha}$ for ${\bm{r}}_{i}\in\alpha$. Then the Hamiltonian (6) is expressed as

[TABLE]

where

[TABLE]

${\bm{k}}(t)={\bm{k}}-\frac{q}{\hbar}{\bm{A}}(t)$ and

[TABLE]

Here, $\tilde{\sum}_{\bm{k}}$ means to calculate the sum over ${\bm{k}}$ expressed as Eq. (8) with $i_{L}=0,1\cdots n-1$ and $i_{T}=0,1\cdots N_{T}-1$.

Next, we show that the Hamiltonian for each ${\bm{k}}$ can be block diagonalized into $2\times 2$ matrices, each of which corresponds to the Hamiltonian of graphene at a certain momentum. We introduce the new basis for electrons at each ${\bm{k}}$ as Murakami and Schüler (2022)

[TABLE]

Then, we have

[TABLE]

with

[TABLE]

We note that each $2\times 2$ matrix in ${\bf h}^{\prime}({\bm{k}})$ has the same form as that of graphene at momentum ${\bm{k}}$ and ${\bm{k}}+{\bm{l}}_{0}$ Murakami and Schüler (2022). If we extend the region of ${\bm{k}}$ for CNTs as Eq. (8) with $i_{L}=0,1\cdots 2n-1$ and $i_{T}=0,1\cdots N_{T}-1$ and remove $\pm$ in $\hat{c}_{{\bm{k}}A_{\pm}}$ and $\hat{c}_{{\bm{k}}B_{\pm}}$, we have the expression used in the main text;

[TABLE]

with

[TABLE]

In the following, we express the Hamiltonian (15) for each ${\bm{k}}$ as $\hat{H}_{\bm{k}}(t)$. The expression tells that for CNTs one needs to choose 2n lines along ${\bm{t}}_{0}$ (${\bm{k}}_{\parallel}$ in the main text) from the momentum space of graphene and the corresponding Hamiltonian for each point is the same as graphene. The magnetic field ${\bm{A}}_{B}$ shifts the position of lines in the direction of ${\bm{l}}_{0}$ (${\bm{k}}_{\perp}$ in the main text).

The explicit expression of the current is

[TABLE]

where $N=2nN_{T}=2N_{c}$. In the following, we express the contribution from each ${\bm{k}}$ as $\hat{J}_{\parallel,{\bm{k}}}(t)$, i.e. $\hat{J}_{\parallel}(t)=\frac{1}{N}\sum_{\bm{k}}\hat{J}_{\parallel,{\bm{k}}}(t)$.

I.2 Dirac model around $K$ and $K^{\prime}$

The TB model for graphene hosts two Dirac points ${\bm{K}}$ and ${\bm{K}}^{\prime}$, where $F({\bm{K}})=0$ and $F({\bm{K}}^{\prime})=0$. To be specific, ${\bm{K}}={\bm{l}}_{0}+\frac{1}{3}{\bm{t}}_{0}$ and ${\bm{K}}^{\prime}=\frac{2}{3}{\bm{t}}_{0}$. The Hamiltonian  (15) around these points can be approximated by the Dirac model. By expanding $F({\bm{k}})=0$ around ${\bm{K}}$ or ${\bm{K}}^{\prime}$ and regarding $e^{-i{\bm{K}}\cdot{\bm{b}}_{1}}\hat{c}^{\dagger}_{{\bm{k}}B}$ ($e^{-i{\bm{K}}^{\prime}\cdot{\bm{b}}_{1}}\hat{c}^{\dagger}_{{\bm{k}}B}$ ) as new $\hat{c}^{\dagger}_{{\bm{k}}B}$, $\hat{H}_{\bm{k}}(t)$ is approximated as

[TABLE]

with

[TABLE]

Here, ”+” is for ${\bm{K}}$ and ”-” is for ${\bm{K}}^{\prime}$, $k_{x}$ and $k_{y}$ are the x and y components of ${\bm{k}}$. In this expression, we shift the origin of the momentum to ${\bm{K}}$ or ${\bm{K}}^{\prime}$. The corresponding expression of $\hat{J}_{\parallel,{\bm{k}}}$ is

[TABLE]

Here, ${\bm{e}}_{x}$ and ${\bm{e}}_{y}$ are the unit vectors for the x and y directions.

I.3 Semiconductor Bloch Equation

In the dipolar gauge Li et al. (2020); Schüler et al. (2021); Murakami and Schüler (2022), the semiconductor Bloch equation (SBE, i.e. the von Neumann equation for the single particle density matrix) is expressed as

[TABLE]

where $[\partial_{t}\bm{\rho}_{\bm{k}}(t)]_{\rm corr}$ indicates the collision integral. The collision integral expresses the scattering, relaxation and dephasing processes originating from electron-electron interactions, electron-phonon couplings and/or impurity scatterings. In principle, these terms can be evaluated microscopically Stefanucci and Leeuwen (2013); Schüler et al. (2020); Ridley et al. (2022). However, the direct microscopic evaluation is computationally expensive. Therefore, in this paper, we employ the relaxation time approximation as

[TABLE]

Here $T_{1}$ corresponds to the relaxation of the band population, while $T_{2}$ represents the interband dephasing time. We assume that initially ($t=0$) the system is in equilibrium at temperature $T$ and there is no electric field ${\bm{A}}_{E}(0)=0$. $\bm{\rho}_{{\rm eq},{\bm{k}}}$ indicates the equilibrium single particle density matrix for ${\bm{k}}$ at this initial temperature. The meaning of these terms becomes clear if one expresses the SBE using the basis that diagonalizes the Hamiltonian at each time (the Houston basis). For the expression, we introduce a unitary matrix ${\bf U}({\bm{k}})$ that diagonalizes ${\bf h}({\bm{k}})$ as ${\bf U}^{\dagger}({\bm{k}}){\bf h}({\bm{k}}){\bf U}({\bm{k}})={\bm{\epsilon}}({\bm{k}})$ with ${\bm{\epsilon}}({\bm{k}})={\rm diag}[\epsilon_{0}({\bm{k}}),\epsilon_{1}({\bm{k}})]$. The single-particle density matrix in the new basis is $\bm{\rho}^{\rm H}_{{\bm{k}}}(t)\equiv{\bf U}^{\dagger}({\bm{k}}(t))\bm{\rho}_{\bm{k}}(t){\bf U}({\bm{k}}(t))$. The diagonal component of $\bm{\rho}^{\rm H}_{\bm{k}}$ indicates the occupation of the bands, and the off-diagonal component means the interband hybridization (coherence). In equilibrium, $\bm{\rho}^{\rm H}_{{\rm eq},{\bm{k}}}$ only hosts the diagonal terms. In this basis, the SBE equivalent to Eq. (22) is

[TABLE]

with

[TABLE]

Here, $\bm{\rho}_{{\rm diag}}$ indicates a matrix consisting of diagonal components of $\bm{\rho}$, while $\bm{\rho}_{{\rm off}}$ indicates a matrix consisting of off-diagonal components of $\bm{\rho}$. ${\bf d}_{a}({\bm{k}})=i{\bf U}^{\dagger}({\bm{k}})[\partial_{a}{\bf U}({\bm{k}})]$ is the (non-abelian) Berry connection for the direction $a$, which plays the role of dipole matrix elements. $\partial_{a}=\partial/\partial k_{a}$ is the momentum derivate for the direction $a$. Note that Eqs. (21) and (23) are equivalent to the conventional expression of SBE in the length gauge Wilhelm et al. (2021); Yue and Gaarde (2022); Murakami and Schüler (2022).

As explain in the main text, the origin of the radiation is the time-dependent change in the electric polarization $P\;(\equiv q\sum_{i}{\bm{e}}_{\parallel}\cdot{\bm{r}}_{i}\langle\hat{n}_{i}\rangle)$, or the corresponding current Huttner et al. (2017); Yue and Gaarde (2022). As in previous works Ishikawa (2010); Dimitrovski et al. (2017), we introduce the intraband and interband currents based on the Houston basis as

[TABLE]

The intraband current is related to the diagonal components of $\bm{\rho}^{\rm H}_{\bm{k}}(t)$. Hence this current originates from the change in the band occupation. On the other hand, the interband current originates from the off-diagonal components of $\bm{\rho}^{\rm H}_{\bm{k}}(t)$ (the interband components). Given that the current is the time derivative of the polarization, the interband current is closely related the change in the interband polarization $P_{\rm inter}$. It is the polarization originating from the interband components as

[TABLE]

In the actual implementation, we evaluate currents ($J_{\parallel},J_{\parallel,\rm intra},J_{\parallel,\rm inter}$) using the difference of the single-particle density matrix from equilibrium $\delta{\bm{\rho}}_{\bm{k}}(t)=\bm{\rho}_{\bm{k}}(t)-\bm{\rho}_{{\rm eq},{\bm{k}}(t)}$. Note that these currents are zero in equilibrium. This strategy has two benefits. Firstly, one can save the calculation cost since away from the Dirac points $\delta{\bm{\rho}}_{\bm{k}}(t)$ is essentially zero and we do not need to evaluate them. Secondly, this can be directly applied for the Dirac model. Namely, the currents for CNTs and graphene can be evaluated from the Dirac models using $\delta\bm{\rho}_{{\bm{k}}}(t)$ at ${\bm{K}}$ and ${\bm{K}}^{\prime}$ and summing up these contributions. If we evaluate the currents from the Dirac models using $\bm{\rho}_{{\bm{k}}}(t)$ the results depend on the momentum cutoff.

I.4 Comments on previous studies

Although the manipulation of HHG in a single-wall CNT using the Aharanov-Bohm effect is our main point, there have been some studies on the effects of the static magnetic field on HHG in CNTs Khosravi et al. (2009); Dvuzhilova et al. (2021). In Ref. Khosravi et al., 2009, they consider the mangetic field applied to the perpendicular direction of the tube and numerically show the sensitivity of the HHG spectrum for the excitation with the light of the wave length 800nm. Thus, the Aharanov-Bohm effect and the Dirac physics are not relevant here. In Ref. Dvuzhilova et al., 2021, they study the propagation of the Airy beam with 1 fs duration in the photonic crystal consisting of many CNTs under magnetic field. We note that the setup and the timescale involved are very different from our case.

II Derivation of Eq. (4) in the main text

We focus on the case of $B=0$ and the contribution from the bands passing the Dirac points. The AC electric field is applied along ${\bm{e}}_{\parallel}(=\begin{bmatrix}\cos\theta\\ \sin\theta\end{bmatrix})$. The momentum on the line passing $K$ or $K^{\prime}$ can be expressed as ${\bm{k}}=k{\bm{e}}_{\parallel}$. The corresponding Hamiltonian of the Dirac model is

[TABLE]

where $k(t)\equiv k-\frac{q}{\hbar}A_{E}(t)$, $\theta_{+}=\theta+\pi$ for ${\bm{K}}$ and $\theta_{+}=-\theta$ for ${\bm{K}}^{\prime}$. We introduce a new basis (the Houston basis) as

[TABLE]

Then we can express the Hamiltonian as

[TABLE]

and the current as

[TABLE]

We introduce the single-particle density matrix in the present basis as $\rho^{R}_{k}(t)=\langle\hat{c}_{kR}^{\dagger}(t)\hat{c}_{kR}(t)\rangle$ and $\rho^{L}_{k}(t)=\langle\hat{c}_{kL}^{\dagger}(t)\hat{c}_{kL}(t)\rangle$. The corresponding SBE becomes

[TABLE]

with

[TABLE]

Here $X=R,L$ and $\rho^{X}_{\rm eq}(k(t),T(t))$ represents the equilibrium single-particle density matrix at the momentum $k(t)$ and the temperature $T(t)$. Note that this expression is slightly extended from Eq. (22). Namely, we consider the time-dependent temperature for the equilibrium state that the system is relaxed to at each time. Furthermore, the SBE tells that no interband component in the density matrix is induced by the electric field, which already indicates the absence of the interband current. Using the integrated form of the SBE, the deviation of the single-particle density matrix from equilibrium $\delta\rho_{\bm{k}}^{X}(t)=\rho_{\bm{k}}^{X}(t)-\rho^{X}_{\rm eq}(k(t),T(t))$ can be expressed as

[TABLE]

Note that, at $t=0$, the system is in equilibrium at the temperature $T(0)$ and there is no electric field.

Now, if we focus on the contributions to the current from the line passing $K$ and $K^{\prime}$ points, we have

[TABLE]

where

[TABLE]

Note that the temperature effect vanishes in the integral of the first two terms in the second line. Finally, we obtain

[TABLE]

with $[\partial_{t}J_{\parallel}(t)]_{\rm corr}\equiv\frac{1}{n}\frac{\sqrt{3}a^{2}qt_{\rm hop}}{4\pi\hbar}\int dk\{[\partial_{t}\rho_{k}^{R}(t)]_{\rm corr}-[\partial_{t}\rho_{k}^{L}(t)]_{\rm corr}\}$. Taking the derivative and substituting Eq. (33) to $[\partial_{t}\rho_{k}^{X}(t)]_{\rm corr}$, we have

[TABLE]

This is Eq. (4) in the main text. Our derivation tells that, even if we improve the relaxation time approximation considering the time-dependent temperature, HHG is absent in the gapless 1d Dirac system.

III Derivation of $\sigma^{(3)}$

We derive the third-order optical conductivity by perturbatively solving the SBE for the Dirac model as in Ref. Cheng et al., 2015. To this end, we change the gauge from the dipolar gauge, which is mainly used in this paper, to the length gauge expressed with the band basis Murakami and Schüler (2022). From the Dirac Hamiltonian Eq. (3) in the main text, we reach the SBE in the length gauge

[TABLE]

where

[TABLE]

Here, $k\in[-\infty,\infty]$ represents the momentum along ${\bf k}_{\parallel}$, $\bm{\rho}^{\rm LB}_{k}$ is the $2\times 2$ single-particle density matrix in the length gauge with the band basis, $\tilde{t}_{\rm hop}\equiv\frac{\sqrt{3}at_{\rm hop}}{2}$ and $K_{B}\equiv\frac{q}{\hbar}A_{B}$. In the following, we set $\kappa\equiv\sqrt{k^{2}+K_{B}^{2}}$, $r_{k}\equiv\frac{K_{B}}{2\kappa^{2}}$, denote the band index of the matrix as $s=0,1$, and omit ”LB” from $\bm{\rho}^{\rm LB}_{k}$.

To solve the SBE perturbatively, we expand $\bm{\rho}_{k}$ as

[TABLE]

with $\rho^{(l)}_{s_{1}s_{2},k}(t)\propto E^{l}$. Note that we have $\rho_{s_{1}s_{2},k}^{(0)}(t)=\delta_{s_{1},s_{2}}n_{s_{1},k}$, where $n_{s,k}=[e^{\epsilon_{s,k}/{k_{B}T}}+1]^{-1}$ and $\epsilon_{s,k}=(-)^{s}\tilde{t}_{\rm hop}\kappa-\mu$.

By substituting Eq. (41) into Eq. (39), we have

[TABLE]

where $\Gamma_{1}\equiv\hbar/T_{1}$, $\Gamma_{2}\equiv\hbar/T_{2}$ and $\epsilon_{g,k}\equiv\epsilon_{0,k}-\epsilon_{1,k}$. We express the SPDM as

[TABLE]

where $\omega_{0}\equiv\omega_{2}+\omega_{3}$ and $\omega\equiv\omega_{1}+\omega_{2}+\omega_{3}$. Using Eq. (42), one can evaluate $\mathcal{P}^{(l)}$ sequentially from $l=1$. For $l=1$, we have

[TABLE]

where $\nu_{3}\equiv\hbar\omega_{3}+i\Gamma_{1}$, $\vartheta_{3}\equiv\hbar\omega_{3}+i\Gamma_{2}$, $\Delta n_{k}\equiv n_{0,k}-n_{1,k}$. For $l=2$, we have

[TABLE]

where $\nu_{0}\equiv\hbar\omega_{0}+i\Gamma_{1}$ and $\vartheta_{0}\equiv\hbar\omega_{0}+i\Gamma_{2}$.

For $l=3$, we have

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

Here $\nu\equiv\hbar\omega+i\Gamma_{1}$ and $\vartheta\equiv\hbar\omega+i\Gamma_{2}$.

Now we evaluate the third order conductivity. By taking into account the contribution from the two Dirac points ($K$ and $K^{\prime}$), the third order current is expressed as

[TABLE]

where $v_{s_{2}s_{1},k}$ denotes an element of the velocity operator

[TABLE]

We define the unsymmetrized third order conductivity as

[TABLE]

where we have

[TABLE]

Note that in Eq. (52) the contribution from $s_{1}\neq s_{2}$ corresponds to the interband current ($\sigma^{(3)}_{\rm er}$), while that from $s_{1}=s_{2}$ corresponds to the intraband current ($\sigma^{(3)}_{\rm ra}$).

In the following, we focus on $T=0$ and half filling ($\mu=0$). In this case, we have $\Delta n_{k}=-1$ and $P_{1;s,k}=P_{3;s,k}=P_{5;s,k}=P_{7;s,k}=0$. Then, we have $\sigma^{(3)}_{\rm ra}=S_{2}+S_{4}$ and $\sigma^{(3)}_{\rm er}=S_{6}+S_{8}$, where

[TABLE]

The concrete expressions for these terms are

[TABLE]

The third harmonic generation corresponds to the case of $\omega_{1}=\omega_{2}=\omega_{3}=\Omega$. In the main text, we denote $\sigma^{(3)}(\Omega,\Omega,\Omega)$ simply as $\sigma^{(3)}$.

Now we consider the limit of $B\rightarrow 0$ ($K_{B}\rightarrow 0$). One can show that the leading terms are

[TABLE]

These terms diverge as $B\rightarrow 0$, and lead to the expression of $\sigma^{(3)}$ shown in the main text. The rest contributions do not diverge (regular). The expression of $\sigma^{(3)}$ tells that the intraband component strongly depends on $T_{1}$, while the interband component strongly depends on $T_{2}$. This can be understood by the following points. The intraband component is related to the occupation of electrons in the conduction and valence bands, and $T_{1}$ is the time scale of the occupation returning to its equilibrium value. On the other hand, the interband component is related to the modification of the interband polarization, i.e. the off-diagonal (interband) component of the density matrix in the Houston basis. These off-diagonal components are zero in equilibrium, and $T_{2}$ is the time scale that the components vanish.

IV Supplementary results

In this section, we show the supplementary results to support the discussion in the main text.

IV.1 Doping effects

In Fig. 6(a), we show the HHG spectra of the armchair CNT with $(n,m)=(15,15)$ away from half filling ($\mu=0.1$eV) with and without the magnetic field $B$. The excitation condition is the same as that in the main text. Although the position of the chemical potential can be much larger than the gap and the excitation frequency, the drastic enhancement of HHG caused by the magnetic field is observed. In Fig. 6(b), we show the intensity of the HHG peaks as a function of $B$. Unlike in the half-filled case, there is no clear optimal value of $B$ which maximizes the peak intensity. This may be attributed to that, in the doped system, charge carriers contributing to the intraband current exist from the beginning. This situation is different from the half-filled system, where charge carriers need to be created by the AC electric field.

IV.2 Excitation with higher frequency

In Fig. 7(a), we show the HHG spectra of the armchair CNT with $(n,m)=(15,15)$ at half filling for the excitation frequency $\hbar\Omega=40$meV. Figure 6(b) shows the corresponding intensity of the HHG peaks as a function of $B$. The drastic enhancement of HHG can be still observed even when the gap is smaller than $\hbar\Omega$. For example, the gap for $B=10$T is about $10$meV. Note that the peaks in the HHG spectrum is less sharp for $\hbar\Omega=40$meV than for $\hbar\Omega=10$meV, since the relaxation and dephasing times become relatively large with respect to the period of the excitation $\frac{2\pi}{\Omega}$. In the range of $B$ in Fig. 7(b) , we cannot observe the optimal value of $B$ unlike in the case of $\hbar\Omega=10$meV. However, it should exist in the larger-$B$ regime. Note that, for $\hbar\Omega=10$meV, the optimal $B$ is around $30\sim 40$T, where the gap is about $30\sim 40$meV.

IV.3 Dependence on chiral index $n$

As we discussed in the main text, in order to achieve the same amount of momentum shift due to the AB effect and thus to open the same size of the gap, an armchair CNT with larger $n$ requires weaker magnetic field $B$. To exemplify this, we compared the HHG spectra for different values of the chiral index $n$ and magnetic field in Fig. 8. To compare the structure of the HHG spectra for different $n$, $I_{\rm HHG}(\omega)$ is multiplied by $n^{2}$. Remember that $I_{\rm HHG}(\omega)$ is defined as the HHG intensity normalized by the number of atoms of a CNT. In Fig. 8, we take the case of $n=15$ and $B=30$T as a reference. For $n=5$, the case with the magnetic field $B=90$T shows essentially the same HHG spectrum as the reference. As we discussed in the main text, this is because i) in the present excitation condition, only the bands closest to $K$ and $K^{\prime}$ contributes to HHG and ii) the gap of the corresponding massive Dirac system is controlled by $A_{B}\propto B\cdot n$. With $B=30$T, the gap is smaller compared to the reference case. For $n=30$, the case with $B=15$T shows essentially the same HHG spectrum as the reference due to the same reason. With $B=30$T, the gap becomes larger than the reference case.

IV.4 Effects of $T_{1}$ and $T_{2}$

In Fig. 9, we show the intensity at the HHG peaks $I_{\rm peak}(l)$ as a function of the electric-field strength $E_{0}$ for $(T_{1},T_{2})=(98.8{\rm fs},98.8{\rm fs})$. This figure can be directly compared with Fig. 4 in the main text, which shows the results for $(T_{1},T_{2})=(98.8{\rm fs},19.8{\rm fs})$. The generic behavior is similar to Fig. 4 in the main text. Namely, the higher order corrections lead to the saturation of the HHG intensity. (The additional increase of $I_{\rm peak}(3)$ around $E_{0}\simeq 100$kV/cm for $B=2$T originates from the deviation from the linear band dispersion.) The deviation from the perturbative regime happens at smaller $E_{0}$ for smaller $B$, i.e. for smaller gap. The optimal value of $B$ that maximizes $I_{\rm peak}$ shifts to larger values with increasing $E_{0}$. On the other hand, unlike the results for $(T_{1},T_{2})=(98.8{\rm fs},19.8{\rm fs})$, the optimal value of $B$ does not decrease continuously with decreasing $E_{0}$. For small $E_{0}$, the HHG intensity is largest around $B=20$T, where the gap is almost twice of $\Omega$.

In order to understand the origin of the different behavior in the weak-electric-field regime between the cases of $(T_{1},T_{2})=(98.8{\rm fs},98.8{\rm fs})$ and $(T_{1},T_{2})=(98.8{\rm fs},19.8{\rm fs})$, we compare the contributions to HHG from the interaband and intraband currents, i.e. $I_{\rm inter}$ and $I_{\rm intra}$, for a weak magnetic field, see Fig. 10. When $T_{1}=T_{2}$, $J_{\rm inter}$ and $J_{\rm intra}$ cancel each other very efficiently, and thus the total intensity $I_{\rm HHG}$ becomes much smaller than $I_{\rm inter}$ and $I_{\rm intra}$, see Fig. 10(a). On the other hand, when $T_{1}>T_{2}$, the cancelation between $J_{\rm inter}$ and $J_{\rm intra}$ is less efficient and $I_{\rm HHG}$ becomes closer to $I_{\rm inter}$ and $I_{\rm intra}$, see Fig. 10(b). Thus, the total HHG intensity becomes much larger in the latter case, although both of $I_{\rm inter}$ and $I_{\rm intra}$ in the former case can be larger than those in the latter case. This result is consistent with the expression of $\sigma^{(3)}$. Namely, for $T_{1}\neq T_{2}$ it diverges with $B\rightarrow 0$ due to the absence of the cancelation between $J_{\rm inter}$ and $J_{\rm intra}$.