Optical selection rules of zigzag graphene nanoribbons

TL;DR

This paper develops an analytical tight-binding model to understand the optical selection rules in zigzag graphene nanoribbons, revealing how wave function parity influences optical transitions and their relation to nanotube spectra.

Contribution

It provides the first analytical expressions for wave functions and optical transition rules in zigzag graphene nanoribbons, highlighting differences and correlations with armchair nanotubes.

Findings

Optical selection rules depend on wave function parity, with ΔJ odd for valence-conduction transitions.

The absorption spectra of ribbons and nanotubes are correlated through van Hove singularities.

Matrix elements show smooth behavior at Dirac points and edge state intersections.

Abstract

We present an analytical tight-binding theory of the optical properties of graphene nanoribbons with zigzag edges. Applying the transfer matrix technique to the nearest-neighbor tight-binding Hamiltonian, we derive analytical expressions for electron wave functions and optical transition matrix elements for incident light polarized along the structure axis. It follows from the obtained results that optical selection rules result from the wave function parity factor $(-1)^J$, where $J$ is the band number. These selection rules are that $\Delta J$ is odd for transitions between valence and conduction subbands and that $\Delta J$ is even for transitions between only valence (conduction) subbands. Although these selection rules are different from those in armchair carbon nanotubes, there is a hidden correlation between absorption spectra of the two structures that should allow one to use them interchangeably in some applications. The correlation originates from the fact that van Hove singularities in the tubes are centered between those in the ribbons if the ribbon width is about a half of the tube circumference. The analysis of the matrix elements dependence on the electron wave vector for narrow ribbons shows a smooth non-singular behavior at the Dirac points and the points where the bulk states meet the edge states.

Figures

Tables

Table 1. Table 1: The absorption peak positions of ZGNR ( 10 ) 10 (10) in the region ω > γ 𝜔 𝛾 \omega>\gamma compared to the estimate ω ~ j subscript ~ 𝜔 𝑗 \widetilde{\omega}_{j} given by Eq. ( 88 ) and transition energies between the states j ​ ( v ) → ( j + 1 ) ​ ( c ) → 𝑗 𝑣 𝑗 1 𝑐 j(v)\rightarrow(j+1)(c) denoted by □ □ \Box ’s and △ △ \triangle ’s in Fig.  8 (e). The index i 𝑖 i numbers the peaks in Fig.  8 (a). The last column presents the energy differences between the numbered subband extrema in Fig.  8 (a). All quantities are measured in the hopping integral γ 𝛾 \gamma .

$i$	$j$	${\omega }_i$	${\tilde{\omega }}_j$	${\square }_v\to {\square }_c$	${△}_v\to {△}_c$	${\square }_i-{△}_{i+1}$
1	2	$1.074$	$1.081$	$1.089$	$1.076$	$1.058$
2	3	$1.509$	$1.511$	$1.527$	$1.518$	$1.491$
3	4	$1.821$	$1.819$	$1.839$	$1.833$	$1.799$
4	5	$1.983$	$1.980$	$2.000$	$1.998$	$1.959$

Equations

$$H=\left(\begin{array}[]{cccccc}0&\gamma q&0&0\\ \gamma q&0&\gamma&0\\ 0&\gamma&0&\gamma q\\ 0&0&\gamma q&0\end{array}\right)$$

$$c_{j-1}\gamma -c_{j}E+c_{j+1}\gamma q$$

$$c_{j-1}\gamma q-c_{j}E+c_{j+1}\gamma$$

$$\displaystyle\left(\begin{array}[]{c}c_{j}\vphantom{-\dfrac{1}{q}}\\ c_{j+1}\vphantom{-\dfrac{1}{q}}\end{array}\right)$$

$$\displaystyle\left(\begin{array}[]{c}c_{j}\\ c_{j+1}\end{array}\right)$$

$$T_1$$

$$\displaystyle\left(\begin{array}[]{c}c_{2p-1}\\ c_{2p}\end{array}\right)$$

$$\displaystyle\left(\begin{array}[]{c}c_{2p}\\ c_{2p+1}\end{array}\right)$$

$$C_{2p+1}$$

$$T=T_{2}\,T_{1}=\left(\begin{array}[]{cc}-\dfrac{1}{q}&\dfrac{\alpha}{q}\\ -\dfrac{\alpha}{q}&\dfrac{\alpha^{2}-q^{2}}{q}\end{array}\right)\,.$$

$${\lambda }^2+\left(\frac{1}{q}+q-\frac{{\alpha }^2}{q}\right)\lambda +1=0.$$

$${\lambda }_{1,2}=A\pm \sqrt{A^2-1},$$

$$A=\frac{{\alpha }^2-q^2-1}{2q}=-\cos \theta .$$

$${\lambda }_{1,2}=-e^{\mp i\theta },$$

$$\alpha =\frac{E}{\gamma }=\pm \sqrt{q^2-2q\cos \theta +1}.$$

$$E=\pm \gamma \sqrt{4{\cos }^2\frac{k}{2}-4\cos \frac{k}{2}\cos \theta +1},$$

$$V_1$$

$${\xi }_{1,2}=\frac{1+q{\lambda }_{1,2}}{\alpha }.$$

$$S$$

$$T^{w}=\dfrac{1}{\xi_{2}-\xi_{1}}\left(\begin{array}[]{cc}\xi_{2}\lambda_{1}^{w}-\xi_{1}\lambda_{2}^{w}&\lambda_{2}^{w}-\lambda_{1}^{w}\\ \xi_{1}\xi_{2}\left(\lambda_{1}^{w}-\lambda_{2}^{w}\right)&\xi_{2}\lambda_{2}^{w}-\xi_{1}\lambda_{1}^{w}\end{array}\right)\,.$$

$$\sin w\theta -2\cos \frac{k}{2}\sin \left[(w+1)\theta \right]=0.$$

$$C_{2p+1}=T^p{C}_1$$

$$\displaystyle=\dfrac{1}{2i}\left(\begin{array}[]{c}\lambda_{1}^{p}-\lambda_{2}^{p}\\ \lambda_{1}^{p}\xi_{1}-\lambda_{2}^{p}\xi_{2}\end{array}\right)$$

$$c_{2p}$$

$$c_{2p+1}$$

$$c_{2p}$$

$$c_{2p+1}$$

$$c_{2p}^{(j)}$$

$$c_{2p+1}^{(j)}$$

$$c_{2n-1}^{(j)}$$

$$c_{2p}^{(j)}$$

$$c_{2p-1}^{(j)}$$

$$c_{2p}^{(j)}$$

$$\{u_{2p-1,2p-1},u_{2p,2p}\}=\{(-1{)}^p,(-1{)}^p\}{|}_{p=1,\dots ,w}.$$

Full text

graphene nanoribbon, selection rules, absorption

I Introduction

Graphene nanoribbons with zigzag edges are quasi-one-dimensional nanostructures based on graphene Novoselov et al. (2004) that are famous for their edges states. These states were theoretically predicted for ribbons with the zigzag edge geometry by Fujita Fujita et al. (1996) and for a slightly modified zigzag geometry by Klein Klein (1994), although the history could be dated back to the pioneering works on polymers Kivelson and Chapman (1983); Tanaka et al. (1987). Since then, edge states in zigzag ribbons have been attracting much attention from the scientific community Nakada et al. (1996); Miyamoto et al. (1999); Wakabayashi et al. (1999); Lin and Shyu (2000); Wakabayashi (2001); Brey and Fertig (2006); Sasaki et al. (2006); Chang et al. (2006); Son et al. (2006a, b); Kohmoto and Hasegawa (2007); Yang et al. (2007a); Malysheva and Onipko (2008); Onipko (2008); Wakabayashi et al. (2010); Dutta and Pati (2010); Chung et al. (2011); Wakabayashi and Dutta (2012); Deng et al. (2013); Deng and Wakabayashi (2014); Chung et al. (2016), because such peculiar localization of the states at the edge of the ribbon should result in the edge magnetization due to the electron-electron interaction. Although the effect was proved to be sound against an edge disorder Nakada et al. (1996), such an edge magnetization had not been experimentally confirmed until quite recently Magda et al. (2014). A fresh surge of interest to physics of zigzag nanoribbons is expected due to the recent synthesis of zigzag ribbons with atomically smooth edges Ruffieux et al. (2016) and a rapid development of the self-assembling technique Schulz et al. (2017).

The edge states in zigzag ribbons have been predicted to be important in transport Deng et al. (2013); Deng and Wakabayashi (2014); Luck and Avishai (2015), electromagnetic Brey and Fertig (2007), and optical properties Lin and Shyu (2000); Chang et al. (2006); Hsu and Reichl (2007). Although considerable attention has been given to zigzag ribbons’ optical properties Lin and Shyu (2000); Chang et al. (2006); Hsu and Reichl (2007); Sasaki et al. (2009); Duan et al. (2010); Sasaki et al. (2011); Gundra and Shukla (2011a); Chung et al. (2011); Sanders et al. (2012, 2013); Zhu et al. (2013); Chung et al. (2016), including many-body effects Yang et al. (2007b); Prezzi et al. (2008); Gundra and Shukla (2011a), the effect of external fields Chang et al. (2006); Gundra and Shukla (2011a), curvature Chung et al. (2016), wave function overlapping integrals Sanders et al. (2012, 2013), the finite length effect Berahman and Sheikhi (2011), and the role of unit cell symmetry Gundra and Shukla (2011b), a number of problems have not been covered yet. In particular, it is known that the optical matrix element of graphene is anisotropic at the Dirac point Grüneis et al. (2003); Saito et al. (2004) due to the topological singularity inherited from the wave functions Ando et al. (2002); Ando (2005). However, the fate of this singularity in the presence of the edge states, i.e., in zigzag nanoribbons, has not been investigated. This requires analysis of the optical transition matrix element dependence on the electron wave vector, in contrast to the usual analysis limited solely to the selection rules.

It was obtained numerically by Hsu and Reichl that the optical selection rules for zigzag ribbons are different from those in armchair carbon nanotubes Hsu and Reichl (2007). By matching the number of atoms in the unit cell of a zigzag ribbon and an armchair tube, it was demonstrated that the optical absorption spectra of both structures are qualitatively different Hsu and Reichl (2007). However, a comparison of these structures based on the matching of their boundary conditions, similar to what has been accomplished for the band structures White et al. (2007) and optical matrix elements Portnoi et al. (2015) of armchair graphene nanoribbons and zigzag carbon nanotubes, has not been reported yet.

The distinctive selection rules of zigzag graphene nanoribbons were noticed as early as 2000 by Lin and Shyu Lin and Shyu (2000). This remarkable and counter-intuitive result, especially when compared to the optical selection rules of carbon nanotubes Ajiki and Ando (1994); Grüneis et al. (2003); Jiang et al. (2004); Goupalov (2005); Malić et al. (2006); Zarifi and Pedersen (2009), was obtained numerically and followed by a few attempts to provide an analytical explanation Chung et al. (2011); Sasaki et al. (2011).

Within the nearest-neighbor approximation of the $\pi$-orbital tight-binding model the optical selection rules for graphene nanoribbons with zigzag edges is a result of the wave function parity factor $(-1)^{J}$, where $J$ numbers conduction (valence) subbands. This factor has been obtained numerically as a connector of wave function components without explicit expressions for the wave functions being presented Chung et al. (2011). Concurrently, the factor $(-1)^{J}$, responsible for the optical selection rules, is missing in some papers providing explicit expressions for the electron wave functions (see Appendix of Ref. Wakabayashi et al., 2010). Although it emerged occasionally in later works dealing with the transport and magnetic properties of the ribbons Wakabayashi and Dutta (2012); Deng and Wakabayashi (2014), its important role was not emphasized and its origin remains somewhat obscure. At the same time, Sasaki and co-workers obtained the optical matrix elements which, although providing the same selections rules, are very different from those in Ref. Chung et al., 2011. Moreover, despite being reduced to the low-energy limit around the Dirac point, the matrix elements in Ref. Sasaki et al., 2011 remain strikingly cumbersome.

It is the purpose of the present paper to demonstrate a simple way of obtaining analytical expressions for optical transition matrix elements in the orthogonal tight-binding model. The essence of this work is an analytical refinement of the paper by Chung et al. Chung et al. (2011), which provides an alternative explanation of the selection rules to that given in terms of pseudospin Sasaki et al. (2011). However, we do not simply derive analytically the results of the study Chung et al. (2011) showing their relation to the zigzag ribbon boundary condition and secular equation, but extend the approach to the transitions between conduction (valence) subbands considered by Sasaki et al. Sasaki et al. (2011). Unlike both mentioned studies, we go beyond a “single point” consideration of the optical matrix elements and analyze the matrix elements as functions of the electron wave vector. The presence of possible singularities in these dependencies at $k=2\pi/3$, corresponding to the Dirac point, and at the transition point $k_{t}$, where the edge states meet bulk states, is in the scope of our study. It is also the purpose of this paper to investigate relations between zigzag ribbons’ and armchair nanotubes’ optical properties by matching their boundary conditions in lieu of matching the number of atoms in the unit cells as was done by Hsu and Reichl Hsu and Reichl (2007).

This paper is organized as follows. In Sec. II we present the tight-binding Hamiltonian and solve its eigenproblem by the transfer matrix method, following the original paper by Klein Klein (1994), in this section, many analogies can be drawn with the treatment of finite length zigzag carbon nanotubes Compernolle et al. (2003); optical transition matrix elements are derived within the so-called gradient (effective mass) approximation, and optical selection rules are obtained. The analytical results are discussed and supplemented by a numerical study in Sec. III. Finally, the summary is provided in Section IV. We relegate to the appendixes some technical details on ribbon wave functions and supplementary results on matching periodic and “hard wall” boundary conditions.

II Analytical tight-binding model

II.1 Hamiltonian eigenproblem

Let us consider a zigzag ribbon within the tight-binding model, which is the orthogonal $\pi$-orbital model taking into account only nearest-neighbor hopping integrals. The atomic structure of a graphene nanoribbon with zigzag edges is presented in Fig.  1. A ribbon with a particular width can be addressed by index $w$, numbering trans-polyacetylene chains — so-called “zigzag” chains.

For such a ribbon, the tight-binding Hamiltonian can be constructed in the usual way by putting $k_{x}\rightarrow 0$, where $k_{x}$ is the transverse component of the electron wave vector. We avoid the procedure described by Klein Klein (1994), since it results in a Hamiltonian for which concerns were raised by Gundra and Shukla Gundra and Shukla (2011b). Thus, for the ribbon with $w=2$, it reads

[TABLE]

where $\gamma$ is the hopping integral and $q=2\cos(k/2)$ with $k=k_{y}a$ being the dimensionless electron wave vector and $a=|\mathbf{\bm{a}}_{1}|=|\mathbf{\bm{a}}_{2}|=2.46$ Å being the graphene lattice constant. The Hamiltonian $H$ has a tridiagonal structure, therefore its eigenproblem can be solved by the transfer matrix method, which is a general mathematical approach for analytical treatment of tridiagonal and triblock diagonal matrix eigenproblems Molinari (1997). This approach was developed and widely used for investigation of one-dimensional systems Kerner (1954); Schmidt (1957); Hori and Asahi (1957); Matsuda (1962). An alternative approach may be based on continuants, which also have been using for the investigation of conjugated $\pi$-carbons such as polyenes and aromatic molecules Lennard-Jones (1937); Coulson (1938) and carbon nanotubes Mestechkin (2005) (see also Refs. Rutherford, 1948, 1952; Muir, 1960).

We use $H$ to derive the relations between the eigenvector components presented in the paper by Chung et al. Chung et al. (2011). In particular, we pay special attention to the origin of the $(-1)^{J}$ factor and its relation to the eigenstate parity. In the rest of this section, we solve the eigenproblem for $H$.

II.1.1 Eigenvalues: proper energy

In this part of the section, we find eigenvalues by the transfer matrix method Kerner (1954); Schmidt (1957); Hori and Asahi (1957); Matsuda (1962). The eigenproblem for the Hamiltonian given by Eq. (1) can be written as follows:

[TABLE]

where $p=1,\dots,w$, $w=N/2$, and $N$ is the number of atoms in the ribbon unit cell. Each of the equations above can be rewritten in the transfer matrix form Schmidt (1957):

[TABLE]

where $\alpha=E/\gamma$. Introducing

[TABLE]

and substituting $j$ into (9) yield

[TABLE]

whence the following recursive relation can be readily noticed:

[TABLE]

and the following transfer matrix equation can be obtained:

[TABLE]

Thus, the transfer matrix in question is

[TABLE]

The characteristic equation for finding the eigenvalues of $T$, $\mbox{det}\left(T-\lambda I\right)=0$, is a quadratic one:

[TABLE]

This equation has the following solution:

[TABLE]

where

[TABLE]

A new variable $\theta$ has been introduced above to reduce the eigenvalues $\lambda_{1,2}$ to the complex exponent form, which is favourable for further calculations:

[TABLE]

where the upper (lower) sign is used for $\lambda_{1}$ ($\lambda_{2}$). We must note that another choice of variable $\theta$, i.e., $A=\cos\theta$, is also possible, but it results in the inverse numbering of the proper energy branches. The minus sign is a better choice because it allows one to avoid a change of the lowest (highest) conduction (valence) subband index when the ribbon width increases.

Equation (43) allows one to express the proper energy in terms of $\theta$ and $q$:

[TABLE]

Taking into account that $q=2\cos(k/2)$, for the proper energy, we obtain

[TABLE]

where $\theta$ is to be found from the secular equation for the fixed ends boundary condition as in the case of a finite atomic chain Schmidt (1957); Hori and Asahi (1957); Matsuda (1962). The physical interpretation of the parameter $\theta$ is to be given further. We note that Eq. (46) has similar form not only to the graphene energy band structure Wallace (1947); Saito et al. (1992, 1998) but also to the eigenenergies of the finite length zigzag carbon nanotubes Compernolle et al. (2003) [cf. with Eq. (32) therein].

II.1.2 Secular equation

For the fixed end boundary condition, which, in the context of the electronic properties being considered, is better referred to as the “hard wall” boundary condition, the general form of the secular equation is $(T^{w})_{22}=0$ Matsuda (1962). This equation can be obtained by imposing the constraint $c_{0}=c_{N+1}=0$ on Eq. (39), where $p=w$, which physically means the vanishing of the tight-binding electron wave functions on sites [math] and $N+1$, or equivalently on zigzag chains [math] and $w+1$ as illustrated in Fig. 1. Hence, for the secular equation, the $w$-th power of the transfer matrix $T$ is needed. The simplest way of calculating $T^{w}$ is $T^{w}=S\Lambda^{w}S^{-1}$, where $\Lambda$ is the diagonal form of $T$ and $S$ is the matrix making the transformation to a new basis in which $T$ is diagonal. The eigenvalues of $T$ are given by Eq. (44), therefore, $\Lambda$ can be easily written down. Concurrently, the $S$ matrix can be constructed from eigenvectors of $T$ written in columns. By setting the first components of the vectors to be equal to unity, one can reduce them to

[TABLE]

where the following notation is used:

[TABLE]

Then the matrix $S$ and its inverse matrix $S^{-1}$ can be written as follows:

[TABLE]

Expressions (57) are of the same form as in the atomic ring problem Hori and Asahi (1957). Using (57), the $T^{w}$ calculation yields

[TABLE]

Now by the aid of (52) and (44) from (58), we can find the explicit form of the secular equation for $\theta$:

[TABLE]

The equation above is very much like that analyzed by Klein Klein (1994) for so-called “bearded” zigzag ribbons, therefore, the same basic analysis can be carried out.

As can be seen from Fig. 2, all nonequivalent solutions of Eq. (59) reside in the interval $\theta\in(0,\pi)$. When the slope of $q\sin[(w+1)\theta]$ at $\theta=0$ is greater than that of $\sin w\theta$ , i.e., $(q\sin[(w+1)\theta])^{\prime}_{\theta=0}>(\sin w\theta)^{\prime}_{\theta=0}\Rightarrow 2\cos(k/2)>w/(w+1)$, there are $w$ different solutions in the interval, which give $2w$ branches of the proper energy (46). This is indicated in Figs. 2 (a) and 2 (b). However, as seen from Figs. 2 (c) and 2 (d), when $2\cos(k/2)\leq w/(w+1)$, one solution is missing and Eq. (46) defines only $2w-2$ branches. The missing solution can be restored by analytical continuation $\theta=i\beta$, where $\beta$ is a parameter to be found. In this case, the secular equation (59) and the proper energy (46) must be modified accordingly by changing trigonometric functions to hyperbolic ones.

The above introduced parameter $\theta$ ($\beta$) can be interpreted as a transverse component of the electron wave vector and the secular equation (59) can be referred to as its quantization condition.

II.1.3 Eigenvectors: wave functions

Let us now find eigenvectors of the Hamiltonian given by Eq. (1). To obtain the eigenvector components, we choose the initial vector $C_{1}=(c_{0},c_{1})$ as a linear combination of the transfer matrix eigenvectors that satisfies the “hard wall” boundary condition $c_{0}=0$ : $C_{1}=(V_{1}-V_{2})/(2i)$. It is to be mentioned here that the opposite end boundary condition, $c_{N+1}=0$, is ensured by Eq. (59). The chosen $C_{1}$ yields

[TABLE]

or, equivalently,

[TABLE]

Substituting (44) and (52) into (63) and keeping in mind the definition of $\alpha$, one readily obtains

[TABLE]

It is worth pointing out that for the starting $p=1$ from the equations above one gets components $c_{2}$ and $c_{3}$. Although it may seem strange because of the missing $c_{1}$, this is how it should be for $c_{1}$ has already been specified by the proper choice of the initial vector $C_{1}$.

Equation (65) can be further simplified (see Appendix A) so that for the eigenvector components, one has

[TABLE]

where we have introduced the index $j$ to number various values of $\theta$, which are solutions of Eq. (59). As one may have noticed the above expressions still have one drawback: $p=1$ defines components $c_{2}$ and $c_{3}$, while it would be much more convenient if $p=1$ would instead specify $c_{1}$ and $c_{2}$. To obtain desired dependence of the eigenvector components on the variable index, one needs to redefine in Eq. (67) the index $p\rightarrow n-1$:

[TABLE]

and then put $n\rightarrow p$. The latter is permissible since $n$ is a dummy index that can be denoted by any letter. Note that due to the change of the terms order in the sine function one $(-1)$ factor in the coefficient $c^{(j)}_{2p+1}$ above cancels, therefore, $j-1$ in the exponent has been replaced by $j$. Thus, for the Hamiltonian (1), we end up with the following eigenvectors:

[TABLE]

where we have got rid of $(-1)$ in $c^{(j)}_{2p}$. Since the whole eigenvector $\left|c^{(j)}\right\rangle=\left(c^{(j)}_{1},c^{(j)}_{2},\ldots,c^{(j)}_{N}\right)$ can be multiplied by any number, one can choose this number to be $(-1)$. Having multiplied $\left|c^{(j)}\right\rangle$ by $(-1)$, one has to change $\pm$ to $\mp$ in the coefficient $c^{(j)}_{2p+1}$, therefore in Eq. (68), the upper “$-$” stands now for the conduction band, while the lower “$+$” for the valence band. The $(-1)^{p}$ factor, however, cannot be eliminated in a similar way because it determines the signs of various components differently. Nevertheless, this factor is of no significance, too, for it can be eliminated by a unitary transform $U$, which is a diagonal matrix with the main diagonal defined as

[TABLE]

For $w=2$, it reads

[TABLE]

As follows from (69), $U$ is both a unitary and an involutory matrix. It can be straightforwardly checked that applying the unitary transform (69) to the eigenvector of $H$ given by Eq. (1), i.e., $\left|\tilde{c}^{(j)}\right\rangle=U\left|c^{(j)}\right\rangle$, we obtain eigenvectors of the Hamiltonian $\tilde{H}=UHU^{\dagger}$. For $w=2$, the explicit form of the new Hamiltonian is

[TABLE]

The general form of the eigenvectors of $\tilde{H}$ is the same as (68) but without $(-1)^{p}$ factor:

[TABLE]

Equations (72) and (68) present components of non-normalized eigenvectors $\left|c^{(j)}\right\rangle$. Normalization constant $N_{j}$ for these vectors can be found from the normalization condition $N^{2}_{j}\left\langle c^{(j)}\right.\left|c^{(j)}\right\rangle=N^{2}_{j}\sum_{p=1}^{w}c^{(j)\ast}_{2p-1}c^{(j)}_{2p-1}+c^{(j)\ast}_{2p}c^{(j)}_{2p}=1$, which yields

[TABLE]

We do not use “ $\widetilde{\phantom{c}}$ ” two distinguish the two types of eigenvectors mentioned above because, by definition, unitary transform preserves the dot product, therefore the normalization constant is the same in both cases.

As in the case of the secular equation, eigenvectors and normalization constants for the missing solution $\theta$ are obtained by the substitution $\theta\rightarrow i\beta$, which results in wave functions being exponentially decaying from the ribbon edges to its interior. These wave functions describe the so-called edge states Klein (1994); Fujita et al. (1996); Nakada et al. (1996). In contrast to them, the wave functions given by normal solutions $\theta_{j}$ extend over the whole ribbon width, therefore they describe the so-called extended or bulk states. It can be shown that normalized eigenvectors’ components for extended and edge states seamlessly match in the transition point $k_{t}$ defined by $2\cos(k/2)=w/(w+1)$ (see Appendix B).

The matching of the bulk and edge state wave functions is shown in Fig. 3, where the wave functions of the zigzag ribbon with $w=15$ are plotted as functions of the atomic site positions $x_{2p-1}=(\sqrt{3}a/2)(p-1)$ and $x_{2p}=(a/2\sqrt{3})+x_{2p-1}$ normalized by the ribbon width $W=x_{2w}$.

Figure 3 presents wave functions for several energy branches $J(s)$, where $J$ is the energy branch number and $s=c$ or $v$ refers to the conduction or valence branch, respectively. As one can see, a bulk state wave function $\left|1(v)\right\rangle$, Fig. 3 (a), transforms into a wave function $\left|1(v)\right\rangle$ predominantly concentrated at the ribbon edges and decaying towards the ribbon center, Fig. 3 (c), by becoming a linear function of $x_{i}/W$ at $k=k_{t}$ as shown in Fig. 3 (b). One can also see that the parity factor can be associated with the mirror or inversion symmetry of the electron wave function. For conduction subbands, if the parity factor $(-1)^{J}$ is positive, then the wave function is symmetric with respect to the inversion center denoted by the large black point as seen for $\left|2(c)\right\rangle$ and $\left|4(c)\right\rangle$ in Figs. 3 (a) and 3 (c). This means the wave function is odd. However, if $(-1)^{J}$ is negative, then the wave function is even, i.e., it is symmetric with respect to the reflection in the dashed dotted line signifying the ribbon center. This happens for $\left|3(c)\right\rangle$ in Fig. 3 (b). For the valence subbands, the behavior is opposite: if $(-1)^{J}$ is negative then the state wave function is odd, as can be seen from Fig. 3 for the subband $1(v)$, but it is even for positive parity factor $(-1)^{J}$. Such behavior is in agreement with the general properties of motion in one dimension Landau and Lifshitz (1977). The parity factor attributed to the mirror symmetry with respect to the line bisecting the ribbon longitudinally (see Fig. 1) has been discussed in the literature Kivelson and Chapman (1983); Hsu and Reichl (2007); Sanders et al. (2012, 2013). In this view, it should be noted that the unit cells of ribbons with odd $w$ do not have such a reflection symmetry (see Fig. 1 for $w=3$), nevertheless as we see from Fig. 3, for such ribbons, the wave functions can still be classified as even or odd in aforementioned sense. This suggests that the symmetry argument developed in Ref. Gundra and Shukla, 2011b as a criterion for the usage of the gradient approximation, which is to be discussed in the next section, is not complete, since in that form it applies only to ribbons with even $w$. Finally, we notice that the state wave functions can be classified by a number of twists of the envelope functions presented in Fig. 3 by dashed and dashed dotted curves. The number of such twists (nodes) is equal to $J(c)$ and $J(v)-1$ for the conduction and valence subbands $J(s)$, respectively. This behavior is similar to what is expected from the oscillation theorem Landau and Lifshitz (1977).

II.2 Optical transition matrix elements

In this section, we study the optical properties of graphene nanoribbons with zigzag edges. Optical transition matrix elements are worked out in the gradient (effective mass) approximation Blount (1962); Johnson and Dresselhaus (1973); Lew Yan Voon and Ram-Mohan (1993); Lin and Shung (1994) and optical selection rules are obtained. However, before moving to the matrix elements of the ribbons, we shall introduce details of optical absorption spectra calculations where these matrix elements are to be used.

Within the first-order time-dependent perturbation theory the transition probability rate between two states, say $\left|\Psi_{f}\right.\rangle$ and $\left|\Psi_{i}\right.\rangle$ having energy $E_{f}$ and $E_{i}$, respectively, is given by the golden rule Anselm (1981):

[TABLE]

where $\delta(\dots)$ is the Dirac deltafunction, and $\hat{H}_{\text{int}}\left(t\right)$ is a time-dependent interaction Hamiltonian coupling a system in question to that causing a perturbation, which is periodic in time with frequency $\omega$. Considering an incident plane electromagnetic wave as a perturbation, one can show in the dipole approximation, $e^{i\mathbf{\bm{k}}\cdot\mathbf{\bm{r}}}\approx 1$, that

[TABLE]

where $\hat{\mathbf{\bm{v}}}$ is the velocity operator, $E_{0}$ is the electric field strength amplitude and $\mathbf{\bm{e}}_{p}$ is the vector of electromagnetic wave polarization. Thus, optical transition matrix elements can be reduced to the velocity operator matrix elements (VMEs).

The total number of transitions per unit time in solids irradiated by electromagnetic wave at zero temperature is a sum of $A_{i\rightarrow f}$ over all initial (occupied) states in the valence band and final (unoccupied) states in the conduction band. To account for losses such as impurity and electron-phonon scattering, the deltafunction in Eq. (74) is replaced by a Lorentzian. The difference in occupation numbers of the initial and final states due to the finite temperature is introduced by the Fermi-Dirac distribution. Then, for the absorption coefficient due to the interband transitions, one has

[TABLE]

where $E_{m,s}(k)$ is the dispersion of the electron in the $m$-th conduction ($s=c$) or valence ($s=v$) subband, $f(E_{m,s}(k))$ is the Fermi-Dirac distribution function, $M_{n(s),m(s^{\prime})}(k)$ is the optical transition matrix element being a function of the electron wave vector, $\Gamma$ is the phenomenological broadening parameter ($0.004\,\gamma$) Lin and Shyu (2000). Note that for nonzero temperature summation over initial states should also include states in the conduction band, therefore indices $s,s^{\prime}$ have been introduced above. The frequency of an incident wave, $\omega$, as well as the electron energy, is measured in the hoping integral $\gamma$.

Similar to Ref. Chung et al. (2011), we follow the prescription of the gradient approximation Blount (1962); Lew Yan Voon and Ram-Mohan (1993) to obtain the velocity operator right from the system Hamiltonian:

[TABLE]

whence for a one-dimensional case,

[TABLE]

with $H$ being the Hamiltonian of the unperturbed system. Note that the derivative $\partial H/\partial k$ is different from $\partial H/\partial\mathbf{\bm{A}}$ mentioned in Ref. Sasaki et al. (2011), where $\mathbf{\bm{A}}$ is the vector potential. The former has a clear relation to the minimal coupling $\mathbf{\bm{k}}\rightarrow\mathbf{\bm{k}}+(e/\hbar)\mathbf{\bm{A}}$ via the expansion $H(\mathbf{\bm{k}}+(e/\hbar)\mathbf{\bm{A}})=H(\mathbf{\bm{k}})+(e/\hbar)\nabla_{\mathbf{\bm{k}}}H\cdot\mathbf{\bm{A}}+\ldots$, where higher-order terms can be neglected for small $\mathbf{\bm{A}}$. Such an approach is equivalent to the effective mass treatment since the commutator $\left[\ldots\right]$ in Eq. (77) implies that the crystal momentum $\mathbf{\bm{k}}$ is an operator:

[TABLE]

which commutes with the position operator in the same way as real momentum $\mathbf{\bm{p}}$, i.e. $\left[x,k_{x}\right]=i$. Note, however, that there is no formal restriction to low energies around the Dirac point, $k=2\pi/3$, as in the $\mathbf{\bm{k}}\cdot\mathbf{\bm{p}}$ theory with the effective mass approximation for graphene DiVincenzo and Mele (1984); Ando (2006), carbon nanotubes Ajiki and Ando (1993, 1996), or graphene nanoribbons Brey and Fertig (2006); Palacios et al. (2010).

In what follows, we proceed with the calculation and analysis of the velocity operator matrix elements (VMEs) in the gradient (effective mass) approximation. Introducing the following vector:

[TABLE]

the VME is evaluated as

[TABLE]

where indices “$c$” and “$v$” denote the conduction and valence band, respectively, and the eigenvectors $\left|c^{(n,m)}\right\rangle$ are meant to be normalized. In Eq. (80), the graphene lattice constant $a$ emerged because, in contrast to the general expression (78), the electron wave vector $k$ is now treated again as a dimensionless quantity.

Let us calculate VMEs for the Hamiltonian $\tilde{H}$ of the form presented by Eq. (71). Similar calculations for $H$ results in the same final expression. Due to the nature of unitary transforms it is not essential which of the Hamiltonians and corresponding eigenvectors one uses. The components of vectors $\left|\tilde{\zeta}^{(j)}\right\rangle$ are

[TABLE]

with upper “$+$” ( lower “$-$”) being used for conduction (valence) subbands. Substituting Eqs. (72) and (II.2) into (II.2), one obtains

[TABLE]

where $S_{n,m}$ is a sum. A similar form of the matrix element was obtained in Ref. Chung et al., 2011 but explicit expressions for the sum $S_{n,m}$ and normalization constants $N_{n}(N_{m})$ were not provided and potential singularities in VME due to $N_{j}$ and $S_{n,m}$ dependence on $k$ were not analysed. Such an analysis has not been carried out elsewhere including Ref. Sasaki et al., 2011.

It is known that the topological singularity in the graphene wave functions Ando et al. (2002); Ando (2005) leads to anisotropic optical matrix element and absorption in the vicinity of the Dirac point Grüneis et al. (2003); Saito et al. (2004). This anisotropy is eliminated in the matrix element of carbon nanotubes Grüneis et al. (2003); Jiang et al. (2004), but the matrix element can exhibit singular behavior at the Dirac point of the tube’s Brillouin zone if a perturbation such as strain, curvature Portnoi et al. (2015) or external magnetic field Portnoi et al. (2009); Moradian et al. (2010); Chegel and Behzad (2014) is applied. The sharp dependence of the zigzag ribbon VME on the electron wave vector around $k=\pm 2\pi/3$ could be triggered by the presence of the edge states. This possibility, however, has not been analysed yet. The VME behavior at the transition point $k_{t}$ has not been investigated either. Being of practical interest Portnoi et al. (2015) this requires a thorough analysis of possible singularities in the VME dependence on $k$. The $S_{n,m}$ sum is given by

[TABLE]

In Eq. (83), normalization constants have been added since the vectors given by Eq. (72) and used for obtaining Eq. (II.2) are not normalized. It is important to allow for normalization constants in the VMEs because otherwise due to their $\theta_{j}$ and therefore $k$ dependency the VME curve’s behavior in the vicinity of the transition point $k_{t}$ is incorrect. It is also worth noting that for $\theta_{n}=\theta_{m}$, or equivalently for $S_{n,n}$, there is an indeterminacy of $\tfrac{0}{0}$ type in the summation result of Eq. (84). This indeterminacy can be easily resolved by L’Hospital’s rule, which yields

[TABLE]

In a similar fashion, one can check that for $\theta_{n}\rightarrow 0$, $S_{n,n}\rightarrow 0$. Note, however, that if $\theta_{n}\rightarrow 0$, then the normalization constant $N_{n}$ given by Eq. (73) becomes infinitely large, thereby introducing indeterminacy into the VME. For transitions between the valence and conduction subbands this indeterminacy is not essential for it is multiplied by an exact zero, originating from the square brackets in Eq. (83), which ensures a zero final result.

As can be seen from Eq. (83), $M_{n(c),n(v)}$ is zero, whereas $M_{n(c),n+1(v)}\sim N_{n}N_{m}S_{n,m}\sin(k/2)$. Thus, optical selection rules are: if $\Delta J=n-m$ is an even integer, then transitions are forbidden, whereas if $\Delta J=n-m$ is an odd integer, then transitions are allowed. The influence of the factor $S_{n,m}$ together with the normalization constants $N_{n}$ and $N_{m}$ on the transition probability, omitted in Ref. Chung et al., 2011, will be discussed in detail in Sec. III. In the remainder of this section, we consider transitions between only conduction (valence) subbands, which are considered in Ref. Sasaki et al., 2011 but are beyond the scope of Ref. Chung et al., 2011.

If the temperature is not zero, then there is a nonzero probability to find an electron in the conduction subband states. Therefore, an incident photon can be absorbed due to transitions between conduction subbands. The same is true for valence subbands, which are not fully occupied. That is why, as has been pointed out above, the summation in Eq. (76) is to be carried out over transitions between conduction (valence) subbands too. Thus, for the absorption coefficient calculation, one also needs VMEs for such transitions. Making use of Eqs. (72) and (II.2), we obtain

[TABLE]

where “$+$” and “$-$” are used for VME of transitions between conduction, $s=c$, and valence, $s=v$, subbands. For the specified transitions, the optical selection rules are the following: transitions are allowed if $\Delta J$ is an even number and they are forbidden otherwise. These matrix elements and corresponding selection rules should be important in spontaneous emission (photoluminescence) calculations Oyama et al. (2006).

In the case of $n=m$, VME given by Eq. (86) is nothing else but the group velocity of an electron in the $n$-th band. If $n=m=1$, then $\theta_{n}=\theta_{m}\rightarrow 0$ as $k$ approaches the transition point $k_{t}$. As a result, in Eq. (86), the indeterminacy arises in precisely the same manner as discussed above for Eq. (83). In the present case, however, it is essential since the expression in square brackets of Eq. (86) is not an exact zero. The indeterminacy can be resolved by the application of L’Hospital’s rule twice. This burden, however, can be bypassed by calculating the VME by the aid of simplified expressions for eigenvectors at the $k_{t}$ provided in Appendix B. Such a calculation yields

[TABLE]

where the upper (lower) sign is used for the conduction (valence) subband. It is easily seen from the expression above that in the limit of a wide ribbon the electron group velocity at $k_{t}\approx 2\pi/3$, i.e., approaching to the Dirac point, is $\mp v_{F}/2$.

Velocity matrix elements for transitions involving edge states can be easily obtained from Eqs. (83) and (86) with $S_{n,m}$ given by Eq. (84) after $\theta\rightarrow i\beta$ replacement being applied. It should be noticed that the Eqs. (83) and (86) obtained here are incomparably simpler than their analogues in Ref. Sasaki et al. (2011) [cf. with Eqs. (18) and (19) therein]. In the next section, we discuss and investigate numerically the obtained results.

III Numerical Results and Discussion

III.1 Electronic properties

The physical properties of graphene nanoribbons are often related to those of carbon nanotubes (CNTs). In particular, one usually compares the electronic properties of graphene nanoribbons with those of carbon nanotubes Hsu and Reichl (2007); Sasaki et al. (2009). In most cases, such a comparison is based merely on the fact that an unrolled carbon tube transforms into a graphene ribbon. However, this approach is a crude one. Firstly, because only zigzag (armchair) ribbons with even number of carbon atom pairs can be rolled up into armchair (zigzag) tubes. Secondly, because a more relevant and subtle comparison requires the matching of boundary conditions. It has been shown by White et al. White et al. (2007) that periodic and “hard wall” boundary conditions can be matched for armchair ribbons and zigzag carbon nanotubes if the width of the ribbons is approximately equal to half of the circumference of the tubes. In Fig. 4, we demonstrate that a similar correspondence of the electronic properties takes place for zigzag graphene nanoribbons with $w$ zigzag chains, ZGNR($w$), and armchair carbon nanotubes, ACNT$(w+1,w+1)$ and ACNT$(w,w)$ depending on which parts of the Brillouin zones are matched (see Appendix C). The impossibility of matching a zigzag ribbon with just one of the tubes arises from the secular equation (59) linking transverse wave vector $\theta$ with the longitudinal wave vector $k$.

For sure, due to the presence of the edge states, one should not expect the transport properties of undoped ribbons to be the same as those of tubes, but the equivalence of the optical properties seems to be quite natural thing. However, this is not the case. As was shown numerically Lin and Shyu (2000); Hsu and Reichl (2007); Chung et al. (2011, 2016) and has been demonstrated above analytically, the optical selection rules of zigzag ribbons are different from those of armchair tubes Ajiki and Ando (1994); Vuković et al. (2002); Grüneis et al. (2003); Jiang et al. (2004); Malić et al. (2006) (see also Appendix D). This leads to transitions between the edge states being forbidden, which should also have important implications for zigzag ribbon based superlattices Saroka and Batrakov (2016); Saroka et al. (2015, 2014). A somewhat similar picture is observed in the bilayer graphene quantum dots of triangular shape, where the edge states are dispersed in energy around the Fermi level Abdelsalam et al. (2016).

III.2 Optical properties

III.2.1 Optical transition matrix elements

To scrutinize the velocity operator matrix elements (VMEs) for allowed transitions we focus on the zigzag ribbon with $w=10$. In Figs. 5 and 6, we plotted the VMEs given by Eqs. (83) and (86) as functions of the electron wave vector in the first Brillouin zone (BZ). Figure 5 includes results for an armchair tube for the sake of comparison. All plots are normalized by the graphene Fermi velocity $v_{F}=\sqrt{3}a\gamma/(2\hbar)$. The arbitrary phase factor of the VMEs, which does not affect their absolute values, was chosen such that it favours plots’ clarity. As in previous sections, we follow the adopted two index notation for the ribbon bands: $J(s)$, where $J=1,\ldots,w$ is the band number and $s=c,v$ is the band type with “$c$” and “$v$” standing for conduction and valence bands, respectively. With this notation in mind one can see that the VME curves for transitions $j(v)\rightarrow(j+1)(c)$ [$(j+1)(v)\rightarrow j(c)$], where $j=1,\dots,w-1$ are shown in Fig. 5 (a). The VME curves for transitions $1(v)\rightarrow 2n(c)$ [$2n(v)\rightarrow 1(c)$], where $n=1,\dots,w/2\text{ or }(w-1)/2$, and for transitions between conduction (valence) subbands only, i.e., $1(s)\rightarrow(2n-1)(s)$, where $n=1,\dots,w/2\text{ or }(w-1)/2$ are presented in Figs. 6 (a) and 6 (b), respectively. As one can see, the VME curves deviate significantly from the previously reported $\sin(k/2)$ behavior Chung et al. (2011); Sasaki et al. (2011), according to which extrema are to be at $k=\pi$, i.e., at the edge of the BZ. The deviation is due to the $S_{n,m}$ and $N_{j}$ given by Eqs. (84) and (73) [see also Eq. (98)], respectively. The shift of the VME curve extrema from the BZ edge is larger for low-energy transitions. Interestingly enough, the positions of these extrema in BZ do not coincide with those of the energy band extrema resulting in the van Hove singularities in the density of states. The curves labeled by $1$ in Figs. 5 (a) and 6 (a) represent direct transitions from the edge states to the closest in energy bulk states. These curves have the largest magnitudes among the ribbons VMEs. However, even for them the maximum absolute values are well below $v_{F}$, in sharp contrast to what is seen in Fig. 5 (b) for ACNT$(11,11)$ (cf. with Refs. Jiang et al. (2004); Portnoi et al. (2009)). Though it is difficult to ignore the fact that the shapes of the VME curves $2$ to $9$ in Fig. 5 (a) are very similar to those obtained for ACNT VMEs in Fig. 5 (b). The most profound curves in Fig. 6 (b) are also labeled by $1$, but they do not have corresponding transitions depicted in the panel to the right. This is because these curves are, in fact, the electron group velocities in $1(v)$ and $1(c)$ subbands given by Eq. (86). As can be seen, at the transition points $k_{t}$ and $k_{t}^{\prime}$ marked by vertical lines the group velocity curves have magnitudes about $v_{F}/2$. This is in accordance with Eq. (87). Ignoring the group velocity curve, one finds that the most prominent magnitudes of VME have transition $1(c)\rightarrow 3(c)$ [$3(v)\rightarrow 1(v)$]. The probability rate described by VMEs of $1(s)\rightarrow(2n-1)(s)$, where $n=2,\dots$, transitions is comparable to that of transitions $1(v)\rightarrow 2n(c)$ [$1(c)\rightarrow 2n(v)$], where $n=2,\dots$, labeled by $2$, $3$, $4$, etc., in Fig. 6 (a) and 6 (b). However, these transitions are less intense compared to $1(v)\rightarrow 2(c)$ [$2(v)\rightarrow 1(c)$], or majority of the $j(v)\rightarrow(j+1)(c)$ [$(j+1)(v)\rightarrow j(c)$], where $j=1,\dots,w-1$, transitions presented in Fig. 5 (a). A regular smooth behavior of all matrix elements at the K(K′) and $k_{t}$ ($k_{t}^{\prime}$) points is worth highlighting, especially for those including $1(s)$ subbands. We noticed, however, that for increasing ribbon width (up to $w=25$) the VME curve peaks for transitions involving $1(s)$ subbands gain a sharper form, therefore a singular VME behavior may still be expected for $1(v)\rightarrow 2(c)$ [$2(v)\rightarrow 1(c)$] transitions in ribbons with $w>25$.

III.2.2 Absorption

It follows from Figs. 5 and 6 (see also Appendix E) that the absorption spectra of zigzag ribbons are mostly shaped by $v\rightarrow c$ transitions with $\Delta J=1$ presented in Fig. 5 (a). However, other transitions may play an important role at certain conditions created by interplay of the doping (or temperature) and ribbon width. To check this, we investigated optical absorption spectra given by Eq. (76) for narrow ribbons with $w=2\ldots 10$. In what follows we discuss ZGNR($6$) for it has the most prominent features and additionally it has been recently synthesized with atomically smooth edges Ruffieux et al. (2016).

Figure 7 compares the absorption spectra of ZGNR($6$) for various positions of the Fermi level, $E_{F}$. As one can see, depending on $E_{F}$ the absorption spectrum has $4$ or $5$ pronounced peaks, which we label in ascending order of their frequency as $A$, $B$, $C$, $D$, and $E$. Peaks $D$ and $E$ are not sensitive to the doping, whereas peaks $A$, $B$, and $C$ are. In contrast to peaks $A$ and $C$ undergoing suppression with increasing $E_{F}$, peak $B$ significantly strengthens. Such different behavior of the three peaks is explained by their different nature.

Let us start with the most interesting case of the peak $B$ at $\omega=0.9\gamma$, which corresponds to the wavelength of about $400$ nm if $\gamma\approx 3$ eV. This peak stems from transitions $1(c)\rightarrow 3(c)$. At $T=0$ K, the valence subbands are fully occupied, therefore, we can safely exclude from the consideration transition $3(v)\rightarrow 1(v)$, which must be blocked due to the exclusion principle. The steep doping dependence of the peak $B$ observed in Fig. 7 (a) has two causes. Firstly, dispersion of subbands $1(c)$ and $3(c)$ and resulting density of states $\sim\left(\partial E_{j,s}(k)/\partial k\right)^{-1}$ presented in Fig. 7 (d). Secondly, the nonzero VMEs for transition $1(c)\rightarrow 3(c)$ in the $k$-interval $\left(2\pi/3,\pi\right)$, as shown in Fig. 7 (b).

Without doping the peak $B$ is absent in the absorption spectrum because both subbands $1(c)$ and $3(c)$ are empty. The introduction of doping results in large number of edge states in the almost flat subband $1(c)$ being occupied with electrons. If the point of the Fermi level intersection with the subband $1(c)$ is denoted as $k_{F}$, then one can say that $k_{F}$ rapidly shifts towards the K point upon ribbon doping. In Figs. 7 (b) and 7 (c), the values of $k_{F}$ for $E_{F}=0.001\gamma$, $0.004\gamma$, and $0.02\gamma$ are marked by vertical lines labeled as \raisebox{-1.2pt}{2}⃝, \raisebox{-1.2pt}{3}⃝, and \raisebox{-1.2pt}{4}⃝, correspondingly. As seen in Fig. 7 (b) at $E_{F}=0.001\gamma$, i.e., $k_{F}=$\raisebox{-1.2pt}{2}⃝, VME of $1(c)\rightarrow 3(c)$ transition represented by curve ‘$2$’ is close to the maximum magnitude, nevertheless, the intensity of the peak $B$ in Fig. 7 (a) presented by curve \raisebox{-1.2pt}{2}⃝ is not that large. The low intensity at such a level of doping is related to the fact that the subband $3(c)$ has a dispersion to the right of the vertical line \raisebox{-1.2pt}{2}⃝, which leads to transitions although being strong contribute into absorption at different frequencies. Upon further increase of the $E_{F}$ up to $0.02\gamma$, i.e., $k_{F}=$\raisebox{-1.2pt}{4}⃝, the VME for $1(c)\rightarrow 3(c)$ transition decreases in magnitude to about $v_{F}/2$. However, due to the flatness of subband $3(c)$ in the vicinity of the band minimum [thick black point in Fig. 7 (c)], all the transitions between lines \raisebox{-1.2pt}{2}⃝ and \raisebox{-1.2pt}{4}⃝ contribute into absorption nearly at the same frequency, which corresponds to the van Hove singularity in the density of states shown in Fig. 7 (d). This results in the sharp enhancement of the peak $B$.

The filling of the subband $1(c)$ with electrons affects all the transitions: $1(c)\rightarrow 3(c)$, $5(c)$, etc. However, in ZGNR($6$) the higher order transition $1(c)\rightarrow 5(c)$ is buried in the peak $C$ for it has lower density of states compared to the subband $4(c)$. To observe higher order transitions one has to take a wider ribbon. Any of the ribbons $w=8,9,10$ can be chosen but ribbon with $w=9$ is the best choice for there transitions $1(c)\rightarrow 5(c)$ results in a clear peak at $\omega\approx\gamma$.

According to our calculations, ZGNR($6$) and ZGNR($7$) are the best choices for a detection of the tunable peak due to $1(c)\rightarrow 3(c)$ transitions. The latter is in agreement with the results of Sanders et al. Sanders et al. (2012, 2013) based on the matrix elements of the momentum and with the wave function overlapping taken into account. For wider ribbons, the peak broadens and loses intensity due to combined effect of the VME and density of states reduction.

As for peaks $A$ and $C$ at $\omega=0.65\gamma$ and $\gamma$ in Fig. 7 (a), they arise from interband transitions $1(v)\rightarrow 2(c)$ [$2(v)\rightarrow 1(c)$] and $1(v)\rightarrow 4(c)$ [$4(v)\rightarrow 1(c)$], respectively. Strictly speaking, many subbands converge into $E=\pm\gamma$ at $k=\pi$, therefore some other transitions also contribute into the peak $C$. By mentioning only one type of transition, we mean the dominant contribution in terms of density of states as indicated in Fig. 7 (d). The intensity of the peak $C$ decreases with doping for it results in the subband $1(c)$ being filled with the electrons whereby transitions $4(v)\rightarrow 1(c)$ are blocked due to the exclusion principle. The same Pauli blocking also takes place for transitions $2(v)\rightarrow 1(c)$, therefore, intensity of the peak $A$ decreases too. A more gentle decrease of peak $A$ intensity compared to that of peak $C$ is due to low doping. As one can see in Fig. 7 (c), for the chosen values of the Fermi level the point $k_{F}$ does not reach position of the subband $2(c)$ minimum. For larger doping, $A$-peak intensity decreases as it happens for peak $C$, and it totally disappears if the doping is high enough to attain the $2(c)$ subband.

The effect of the finite temperature is similar to that of doping discussed above (see Appendix E).

Finally, let us compare the zigzag nanoribbon absorption spectra with those of armchair nanotubes. In Fig. 8 (a), the absorption spectra of ZGNR$(10)$ and ACNT$(11,11)$ are presented together with that of ACNT$(10,10)$. For the sake of comparison, each spectrum is not normalized by the number of atoms in the unit cell. The first peculiarity, which one can notice, is that in the ribbon all but the lowest in energy absorption peaks lose approximately half of their intensity compared to the peaks in the tubes. The second peculiarity is that ZGNR$(10)$ and ACNT$(11,11)$ have the same pattern of absorption peaks in the high frequency range $\omega>\gamma$, which is highlighted in the light blue. Both features are not accidental, as follows from the plots presented in Fig. 8 (b)-8(d) for ribbons and tubes of larger transverse size.

In order to explain the noticed difference and similarity, we focus on ZGNR$(10)$ and ACNT$(11,11)$. Obviously, a large difference in peak intensities between the tube and ribbon cannot be explained only by the velocity matrix elements being higher in the tube than in the ribbon, as follows from Fig. 5, therefore the density of states should be accounted for. Here we do not appeal to the suppression due to the momentum conservation as in Ref. Sasaki et al., 2011 for we regard all transitions, even between subbands with different indices, as direct ones. At the same time, the correlation of the absorption peaks’ positions is to be related to the van Hove singularities in the density of states too. Thus we need to have a closer look at the band structures and density of states of ZGNR$(10)$ and ACNT$(11,11)$. In Fig. 8 (e), the ZGNR$(10)$ band structure (solid curve) is compared with that of ACNT$(11,11)$ (dashed curve). Similar comparison is presented for the density of states in Fig. 8(f). The peaks numbered as $1$, $2$, $3$, and $4$ in Fig. 8 (a) result from the transitions between ACNT$(11,11)$ subband extrema marked by numbered circles in Fig. 8 (e). The same peaks in ZGNR$(10)$ originate from the transitions involving the subband extrema marked in Fig. 8 (e) by the numbered squares (triangles) for the conduction (valence) subbands. Selection rules in both structures allow transitions between the markers of the same shape. Let us be more specific and focus on the peak ‘$1$’. In ACNT$(11,11)$, this peak is due to transition between two van Hove singularities in the density of states. Although the density of states in the tube is nearly twice as high as than that in the ribbon due to the double degeneracy of the tube’s subbands, this cannot explain the difference in the intensities of the tube and ribbons absorption peaks, since, according to the selection rules, two type of transitions with the same frequency are allowed in the ribbon: $3(v)\rightarrow 2(c)[2(v)\rightarrow 3(c)]$. The difference in intensities arises due the fact that positions of the band extrema for adjacent bands in the ribbon are shifted in the $k$-space, thereby each of the specified in Fig. 8 (e) ribbon transitions happens either ‘from’ or ‘to’ the band extrema and not ‘between’ them as happens in the tube. In other words, each of these transitions involves only one van Hove singularity. In this view, the extremely high intensity of the lowest in energy absorption peak in ZGNR$(10)$ arises due to the high density of states originating from the flatness of the $1(c)$ band dispersion at $E=0$.

As one can notice from Fig. 8 (e), the tube subband extrema take middle positions in energy between extrema of adjacent ribbons subbands. This leads to the tube’s and ribbon’s transition energies being very close as illustrated by a parallelogram in Fig. 8 (e). As a result, a correlation between the absorption peak positions arises. To understand the origin of this correlation we need to analyze the positions of the van Hove singularities, which can be derived from the analytical expression for the band structure. However, for a zigzag ribbon, such an expression cannot be obtained in a closed-form from Eq. (46), since the secular equation (59) does not allow expressing of its solution in such a form; though closed-form solutions for two specific cases, $k=2\pi/3$ and $\pi$, have been reported for this type of equation Klein (1994); Compernolle et al. (2003); Wakabayashi and Dutta (2012). On the other hand, since the armchair nanotube band structure has a closed-form given by Eq. (108), the positions of the van Hove singularities and, therefore, the absorption peak positions can be easily obtained for them. Then, a simple analytical expression for ACNT$(w+1,w+1)$ peak positions,

[TABLE]

can be used as an estimation of the absorption peak positions in ZGNR($w$), when $5/6>j/(w+1)>1/6$. In Fig. 8 (a), the vertical dashed lines denote the peak positions given by Eq. (88). As one can see, outside the light blue regions peak positions do not necessarily coincide; the ribbon spectra also have additional peaks outside these regions resulting from transitions involving $1(s)$ subbands and the selection rule $v\rightarrow c$ $\Delta J=1,3,\ldots$, etc. In contrast to this, within the regions $\gamma<\omega<2\gamma$ the above-mentioned correlation takes place for all ribbons with $w>5$. To estimate the reliability of Eq. (88) in Table 1 we compared the numerically calculated peaks positions in the ZGNR$(10)$ with those resulting from Eq. (88). We supplemented these results with numerically evaluated energies of $j(v)\rightarrow(j+1)(c)$, where $j=2,3,4,5$, transitions involving one band extremum state, i.e., those which occur between the states denoted by square ($\Box$) and triangle ($\triangle$) markers in Fig. 8 (e). As seen from Table 1, a deviation of $\widetilde{\omega}_{j}$ from $\omega_{i}$ does not exceed $1\%$ of the hopping integral, i.e., $30$ meV for $\gamma\approx 3$ eV. It also follows from the Table 1 that the above presented picture is a simplified one. In reality, the absorption peaks are averages of all transitions taking place in between the two subband extrema shifted in the $k$-space so that peak positions $\omega_{i}$ and their estimates $\widetilde{\omega}_{j}$ are squeezed between the $j(v)\rightarrow(j+1)(c);(\Box,\triangle)$ transition energies and the energy differences between the corresponding van Hove singularities, $\Box_{i}-\triangle_{i+1}$.

The panels (b)-(d) in Fig. 8 show that the aforementioned correlation may extend to the low-energy region $\omega<\gamma$. This region of a ribbon’s spectrum is dominated by the transitions originating from the edge states. It is evident that the absorption peaks originating from these transitions cannot correlate with the peaks in armchair tubes. In fact, they can only hide this feature. In order to verify our assumption, in Fig. 9 we split the ZGNR$(20)$ absorption spectrum into two parts: ‘part I’ containing only transitions involving the $1(s)$ subband, i.e., edge states, and ‘part II’ containing the rest of the transitions. As can be seen from Fig. 9, it is the latter that correlates with the tubes’ absorption spectrum. Only the first absorption peak in ACNT$(21,21)$ does not have a counterpart in the ribbon spectrum. Thus, Eq. (88) has a broader applicability and with its help the hidden correlation could be verified even by absorption measurements in the optical range. Equation (88) describes zigzag ribbon peak positions when $j=2\ldots w/2$ (even $w$) or $(w-1)/2$ (odd $w$).

The revealed correlation of the absorption peak positions in armchair tubes and zigzag ribbons may be affected by excitonic effects. Excitons are known to be important in one dimensional systems due to the enhanced binding energy Loudon (1959). However, such effects rarely were a subject of investigation in the metallic families of graphene nanoribbons Yang et al. (2008); Lu et al. (2013) and carbon nanotubes Ando (1997); Jiang et al. (2007); Hartmann et al. (2011). Moreover, it seems that attention has never been paid to the high energy transitions, therefore this problem requires a thorough study. Yet, a general qualitative picture says that the positions of the presented peaks should be red-shifted by the amount of the binding energies. These energies can be linked to the system’s transverse size by an analytical phenomenological quasi-one dimensional exciton model, which has been successfully applied to semiconducting quantum wires Bányai et al. (1987); Ogawa and Takagahara (1991a, b) and carbon nanotubes Wang et al. (2005). Then, since the tubes and ribbons in question have comparable widths and diameters, the binding energies and, therefore, shifts are expected to be close for both structures (neglecting the different shapes of their cross-sections), thereby preserving the unveiled correlation in the absorption spectra. Some excitonic states may require a magnetic field for their brightening if they happen to be dark ones Shaver et al. (2007). We should also mention that the correlation reported here can be additionally hidden by a landscape of absorption peaks originating from $\sigma$-orbitals.

IV Conclusions

In summary, we considered the optical properties of zigzag graphene nanoribbons within the orthogonal $\pi$-orbital tight-binding model and effective mass approximation for polarization of the incident radiation parallel to the ribbon axis. It was analytically confirmed that the selection rules between valence and conduction subbands, $\Delta J=n-m$ is odd, and between conduction (valence) subbands only, $\Delta J=n-m$ is even, stem from the wave function parity factor, $(-1)^{J}$, where $J$ is an integer numbering the energy bands. It was also shown that this parity factor originates from the ribbon’s secular equation.

A comprehensive comparison of optical properties between carbon nanotubes and zigzag nanoribbons shows significant differences. Most importantly, the concept of cutting lines Samsonidze et al. (2003); Izumida et al. (2016), or even its generalization to ‘cutting curves’ Wakabayashi and Dutta (2012); Deng and Wakabayashi (2015), being unable to explain selection rules fails with respect to optical properties of zigzag graphene nanoribbons, while it works well for armchair carbon nanotubes. Nevertheless, a proper comparison reveals the absorption spectra of a zigzag nanoribbon and an armchair carbon nanotube have a correlation between the positions of the peaks originating from the $v\rightarrow c$ transitions between the bulk states , if $N_{t}=2N_{r}+4$, where $N_{t,r}$ is the number of atoms in the tube’s (ribbon’s) unit cell, i.e. when the ribbon width is about half of the tube circumference. Putting it differently, this correlation takes place for ZGNR$(w)$ and ACNT$(w+1,w+1)$ if $w>5$.

The analysis of the velocity operator matrix element dependencies on the electron wave vector shows that they have a smooth regular behavior at least up to $w=25$ in the whole Brillouin zone, including the Dirac ($k=\pm 2\pi/3$) and transition ($k=k_{t}$) points. However, the matrix element behavior deviates significantly from the previous estimation $\sim\sin(k/2)$. For all types of transitions the magnitude of the velocity operator matrix elements attain a maximum value for $k\in\pm(\pi/2,\pi)$.

A close examination of the absorption spectra of zigzag ribbons shows they should have temperature and doping dependent absorption peaks originating from transitions between only conduction (valence) subbands, $\Delta J=2,4,\ldots$, etc., which could be tuned, for instance, by a gate voltage. In particular, narrow zigzag ribbons with $w=6,7$ should have such prominent temperature and doping dependent absorption peaks. Although beyond the single electron tight-binding model the energy bands of zigzag ribbons are known to be modified by electron-electron interaction Magda et al. (2014) and the effect of the substrate, we believe that experimental observation of the tunable absorption should be possible as the latter effect, for instance, can be eliminated by system suspension.

Finally, we point out that the obtained velocity matrix elements of single electron transitions can be utilized in further study of excitonic effects via Elliot’s formula for absorption Elliott (1957); Elliott and Loudon (1960).

Acknowledgements.

This work was supported by the EU FP7 ITN NOTEDEV (Grant No. FP7-607521); EU H2020 RISE project CoExAN (Grant No. H2020-644076); FP7 IRSES projects CANTOR (Grant No. FP7-612285), QOCaN (Grant No. FP7-316432), InterNoM (Grant No. FP7-612624); Graphene Flagship (Grant No. 604391). The authors are very thankful to R. Keens and C. A. Downing for a careful reading of the manuscript and to A. Shytov and K. G. Batrakov for useful advice and fruitful discussions.

Appendix A Wave-function parity factor

In order to clarify the origin of the wave-function parity factor, we present in detail the simplification of the eigenvector component (65).

Equation (65) can be further simplified if one expresses $2\cos(k/2)$ in terms of the quantized momentum $\theta$ from the quantization condition (59) as

[TABLE]

and then substitutes the result into the square brackets of Eq. (65):

[TABLE]

Note that the proper energy $E$ entering Eq. (65) can also be re-casted only in terms of $\theta$ by substituting (89) into Eq. (46):

[TABLE]

Now making use of Eqs. (90) and (91), one readily obtains that

[TABLE]

where the upper (lower) sign is applied for the conduction (valence) band state. The first ratio in the expression above is a trivial one, $\tfrac{\sin\theta}{\left|\sin\theta\right|}=1$ for $\theta\in(0,\pi)$. However, the second ratio deserves special attention because, as we will see next, it is a clue to the optical properties of zigzag ribbons.

The magnitude of the second ratio is, of course, unity, but its sign depends upon $\theta$. To determine the sign of the ratio $\tfrac{\left|\sin\left[(w+1)\theta\right]\right|}{\sin\left[(w+1)\theta\right]}$ one needs to analyze it along with the quantization condition (59). Since absolute value is always positive the sign of the ratio is determined by the sign of its denominator defined by the secular equation solutions.

Let us investigate how secular equation solutions, $\theta_{j}$, are spread in the range $(0,\pi)$. For this purpose, one can continuously change the parameter $q$ from [math] to $\infty$ similar to what is presented in Fig. 2. Varying $q$ between the above mentioned limits, one finds that the two values of $q$ determine the left and right ends of the intervals in each of which one $\theta_{j}$ is confined. By putting the parameter $q=0$ into Eq. (59), we get $\sin w\theta_{j}=0$ with $\theta_{j,\text{min}}=\pi(j-1)/w$ being solutions, while $q=\infty$ yields $\sin\left[(w+1)\theta\right]=0$ with $\theta_{j,\text{max}}=\pi j/(w+1)$ as solutions; in both cases $j=1\ldots w$ enumerates solutions. It is worth noting that although the upper value of $q=2\cos(k/2)$ is limited to $2$, we can take a greater value for an estimation because an increase of $q$ above $2$ shifts the initial interval right boundaries so that the original intervals are contained within the new $\theta$-intervals depicted in Fig. 2. The left boundaries of the intervals can also be pushed further left to put all the new intervals within even wider ones:

[TABLE]

With inequalities (93) at hand it is easy to analyze the argument of $\sin\left[(w+1)\theta_{j}\right]$ for it is evident that for all $\theta_{j}$ satisfying inequalities (93) the sine function argument $(w+1)\theta_{j}$ is squeezed between $\pi(j-1)$ and $\pi j$. This leads to positive and negative signs of $\sin\left[(w+1)\theta_{j}\right]$ for odd and even $j$, respectively. Therefore, the second ratio in Eq. (92) can be written as

[TABLE]

where $j$ is an integer being interpreted as the band number.

Appendix B Edge and bulk state eigenvectors at the transition point

Let us obtain the wave functions of the edge states in the explicit form and show how it reduces at the transition point $k_{t}$ defined as a solution of the equation $2\cos(k/2)=w/(w+1)$. As has been mentioned above, to do this one needs to use substitution $\theta\rightarrow i\beta$, which upon application to (72) yields

[TABLE]

with $p=1,\ldots,w$. Note that $j=1$ for bands containing edge states, therefore the parity factor has been ruled out and $\mp$ in (72) has been replaced with $\pm$ in (95). The same substitution applied to the normalization constant (73) leads to

[TABLE]

As one can notice, the expression under the square root of (96) is negative, therefore the imaginary unit resulting form it must cancel with that in (95). Hence, for normalized eigenvector components it can be written

[TABLE]

where $p=1,\dots,w$ and

[TABLE]

Note that the eigenvector (97) does not contain $(-1)^{J}$ factor like Eq. (34) in work Wakabayashi and Dutta (2012). Even for inverse band enumeration, this factor would be $(-1)^{w}$ not $(-1)^{J}$. At the transition point, $\beta_{j}\rightarrow 0$, which results in divergence in (98) if all hyperbolic functions are expanded to the first order. However, using the original definition of the constant:

[TABLE]

where the factor of $2$ is due to the fact that $\sum_{p=1}^{w}\sinh p\beta_{j}=\sum_{p=1}^{w}\sinh\left[(w+1-p)\beta_{j}\right]$, the same first order expansion results in

[TABLE]

Thus, for normalized eigenvectors in the vicinity of the transition point, one has

[TABLE]

where

[TABLE]

The same result can be obtained starting from the eigenvectors (72) and their normalization constant specified as $N_{j}=1/\sqrt{2\sum_{p=1}^{w}\sin^{2}p\theta_{j}}$, therefore wave functions approaching $k_{t}$ from the left and from the right attain the same value. As a result of this seamless transition of one type of functions into another, the VMEs can be obtained as smooth functions of electron wave vector $k$ for the lowest conduction (higherst valence) subbands, i.e., for $j=1$.

It is to be mentioned here that the edge states can be also obtained in zigzag carbon nanotubes with finite length Compernolle et al. (2003). Unlike the case of the infinite ribbon the number of such states is finite in tubes. Recently, it has been shown that this number is related to the winding number Izumida et al. (2015, 2016). However, the state at the transition point, the charge density of which decays quadratically towards the structure center, seems to be less likely in the finite tubes.

Appendix C Periodic boundary conditions

In this part of the appendix, we demonstrate how the fixed end (‘hard wall’) boundary condition employed in this paper for zigzag ribbon investigation is related to the periodic boundary condition that is used for carbon nanotubes. A carbon nanotube of the armchair type (see Ref. Saito et al. (1998) for tubes classification) is unrolled into a graphene nanoribbon with zigzag edges. The tight-binding Hamiltonian of the armchair nanotube differs from that of the zigzag ribbon by the upper right and lower left nonzero elements. For instance, for the ribbon Hamiltonian given by Eq. (1) an equivalent tube Hamiltonian is

[TABLE]

Despite these differences the eigenproblem of such a Hamiltonian reduces to the same transfer matrix equation as Eq. (39). The periodic boundary condition, however, requires $C_{N+1}=C_{1}$, whence it follows that the secular equation is $\det(T^{w}-I)=0$. To obtain the explicit form of the secular equation, one can use (58), but there is a faster way if one uses the following relation Kerner (1954):

[TABLE]

Using the above relation and taking into account that $\det T=1$, the secular equation can be recasted as

[TABLE]

The cyclic property of the trace operation allows further simplification of the secular equation:

[TABLE]

where $\Lambda$ is a diagonal form of the transfer matrix $T$ with the diagonal elements given by $\lambda_{1,2}=e^{\pm i\theta}$, i.e. a new variable $\theta$ is defined as $A=\cos\theta$ [cf. with Eq. (43)], $S,S^{-1}$ are given by Eqs. (57). Such treatment is equivalent to that with $\lambda_{1,2}$ given by Eq. (43), the difference is in subband enumeration similar to that mentioned for the hard wall boundary condition. In Fig. 4, the tube’s band enumeration, we refer to as direct one, corresponds to $A=-\cos\theta$. The above chosen inverse enumeration, $A=\cos\theta$, is shown in the right panel of Fig. 5. It was chosen to obtain the tube’s energy bands in a form close to graphene energy bands Wallace (1947); Saito et al. (1992, 1998). Thus, for an armchair tube secular equation, we end up with

[TABLE]

whence it is evident that $\theta_{j}=2\pi j/w$ with $j$ being an integer numbering solutions and $w=N/2$ with $N$ being the number of carbon atoms in the tube’s unit cell. To obtain the tube energy bands $\theta_{j}$ should be substituted into $\pm\gamma\sqrt{q^{2}+2q\cos\theta+1}$, which yields

[TABLE]

where we use $j$ for the band numbering.

In the case of the hard wall boundary condition and variable $\theta$ introduced as above, i.e. with the reverse enumeration of the ribbon bands, the secular equation has the form:

[TABLE]

The proper energy is obtained by substituting solutions of this equation into $\pm\gamma\sqrt{q^{2}+2q\cos\theta+1}$. Solutions of (109) can be found in the zero approximation by setting $k=0$; ideally, one should set $q=2\cos(k/2)\rightarrow\infty$. This leads to $\sin\left[(w+1)\theta\right]=0$ with $\theta_{j}=\pi j/(w+1)$ being solution. Equating $\theta_{j}$ obtained for a tube and ribbon, one gets:

[TABLE]

where $N_{t,r}$ is the number of atoms in the unit cell of the tube and ribbon, respectively. As follows from (110) if

[TABLE]

then the proper energies are approximately equal at $k=0$. It is also possible to consider the opposite limit when $k=\pi$, which leads to $\theta_{j}=\pi j/w$ in the case of the ribbon. The usage of this $\theta_{j}$ results in a better match of the ribbon and tube energies close to the edge of the Brillouin zone, i.e., at $k=\pi$, if the following relation holds between the number of atoms in the structures: $N_{t}=2N_{r}$.

Appendix D Armchair nanotube selection rules

In this section, we derive selection rules for transitions in armchair carbon nanotubes (ACNTs). In spite of being known for a long time Ajiki and Ando (1994); Grüneis et al. (2003); Jiang et al. (2004); Goupalov (2005); Malić et al. (2006); Zarifi and Pedersen (2009), they have not been derived from the full tight-binding Hamiltonian. The purpose of this exercise is to provide deeper understanding of the difference in the optical properties of zigzag graphene nanoribbons and ACNTs and also to show their relation to the graphene single layer sheet.

To calculate velocity operator matrix elements, one needs the wave functions. Substitution of Eq. (107) solution $\theta_{j}=\frac{2\pi j}{w}$ into $T^{w}-I$ gives a zero matrix. Hence the boundary condition $C_{N+1}=C_{1};\rightarrow(T^{w}-I)C_{1}=0$ is fulfilled for any components of the initial vector $C_{1}$. We see that for the periodic boundary condition the initial vector $C_{1}$ can be an arbitrary one. The most reasonable choice of $C_{1}$ is one of the eigenvectors (51). Let it be $V_{2}$. Then, with $\lambda_{1,2}=e^{\pm i\theta}$ the wave-function components can be found from Eq. (39) as follows:

[TABLE]

where $p=1,\ldots,w$, $f_{j}=1+qe^{-i\theta_{j}}$, and we have changed the order of the components as it was done for Eq. (67). Introducing new function $\tilde{f}_{j}=e^{i\theta_{j}/3}f_{j}$ into Eq. (112) and applying the unitary transform $U_{j}=\{u_{2p-1,2p-1},u_{2p,2p}\}=\{e^{i\theta_{j}(p-2/3)},e^{i\theta_{j}p}\}|_{p=1\ldots w}$ to the vector $\left|c^{(j)}\right\rangle$, we obtain

[TABLE]

where $p=1,\ldots,w$. The normalization constant $N_{j}=1/\sqrt{2w}$ for $\left|\tilde{c}^{(j)}\right\rangle$ and it is independent of $\theta_{j}$.

As one can see, the unitary matrix $U_{j}$ depends on the band index $j$, therefore the new Hamiltonian that preserves the matrix element upon the transfromation of the $\left|c^{(n,m)}\right\rangle$ vectors is $\tilde{H}=U_{n}HU_{m}^{\dagger}$. However, such a Hamiltonian satisfies the time independent Schrodinger equation only if $n=m$. This is the selection rule for ACNT optical transitions, which also means all transitions $c\rightarrow c$ and $v\rightarrow v$ are forbidden.

For $\tilde{H}=U_{j}HU^{\dagger}_{j}$ the components of the vectors $\left|\tilde{\zeta}^{(j)}\right\rangle$ are

[TABLE]

with the upper “$-$” ( lower “$+$”) being used for the conduction (valence) subbands. By putting Eqs. (113) and  (D) into Eq. (II.2), and accounting for the normalization constant $N_{j}$, for allowed transitions we have

[TABLE]

Similarly, calculations for the group velocity yields

[TABLE]

where “$+$” (“$-$”) refers to the conduction (valence) subbands.

The same result is obtained from the graphene Hamiltonian and eigenvectors: $(1/\hbar)\left\langle c_{c}\right|\partial H/\partial k_{y}\left|c_{v}\right\rangle$ with $H_{11}=H_{22}=0$, $H_{12}=H^{\ast}_{21}=\gamma\left(e^{ik_{x}a/\sqrt{3}}+2e^{-ik_{x}a/2\sqrt{3}}\cos(k_{y}a/2)\right)$ and $k_{x}=2\pi j/C_{h}$, where $C_{h}$ is the tube circumference and $a=2.46$ Å is the graphene lattice constant. If $\theta_{j}=\sqrt{3}k_{x}a/2$, $k=k_{y}a$, and the tube chiral index is $w/2$, then $k_{x}=4\pi j/(\sqrt{3}aw)=2\pi j/C_{h}$. Hence, Eq. (D) can be restored by cutting graphene’s optical transition matrix elements along the lines specified by the quantization of $k_{x}$. Finally, we note that a calculation of the matrix elements with the eigenvectors (112) and the Hamiltonian (103) also provides straightforward justification of the selection rules for it results in zero matrix elements when $n\neq m$.

Appendix E Supplementary results

For the sake of completeness, in Fig. 10, we present VME curves obtained for transitions between the lower (higher) energy valence (conduction) subbands.

These transitions can be referred to as $j(s)\rightarrow(j+2)(s)$, where $j=1,\ldots,w-2$. Noticing that the curve labeled by $1$ in Fig. 10 is the same as the curve labeled by $2$ in Fig. 6 (b), one easily sees that the transitions labeled from $2$ to $7$ are much weaker compared to the transitions in Fig. 6. Unlike the VME curves in Figs. 5 (a) and 6, all curves of $j(s)\rightarrow(j+2)(s)$ transitions converge to zero at the edge of the BZ and have extrema decreasing in magnitude and shifting from the K(K′**) point towards the BZ edge for greater $j$’s.

Figure 11 shows that temperature has a similar influence on the absorption spectra to doping. The observed changes are explained in the same way as presented for Fig. 7. The peak due to the transitions $1(c)\rightarrow 3(c)$ is weaker and broader for ZGNR($9$) compared to that in ZGNR($6$). At the same time, the peak at $\omega=\gamma$ due to transitions $1(c)\rightarrow 5(c)$ is quite intense.