Coherent control of current injection in zigzag graphene nanoribbons

TL;DR

This paper uses Fermi's golden rule to analyze how coherent optical processes can control current injection in zigzag graphene nanoribbons, highlighting their sensitivity to edge state occupation and the effects of doping.

Contribution

It introduces a detailed theoretical framework for coherent control of current injection in zigzag graphene nanoribbons, emphasizing the role of localized edge states and doping effects.

Findings

Rich low-energy optical structure due to edge states

Doping significantly alters absorption and current injection

Resonant processes enhance control at specific energies

Abstract

We present Fermi's golden rule calculations of the optical carrier injection and the coherent control of current injection in graphene nanoribbons with zigzag geometry, using an envelope function approach. This system possesses strongly localized states (flat bands) with a large joint density of states at low photon energies; for ribbons with widths above a few tens of nanometers, this system also posses large number of (non-flat) states with maxima and minima close to the Fermi level. Consequently, even with small dopings the occupation of these localized states can be significantly altered. In this work, we calculate the relevant quantities for coherent control at different chemical potentials, showing the sensitivity of this system to the occupation of the edge states. We consider coherent control scenarios arising from the interference of one-photon absorption at $2\hbar\omega$ with two-photon absorption at $\hbar\omega$, and those arising from the interference of one-photon absorption at $\hbar\omega$ with stimulated electronic Raman scattering (virtual absorption at $2\hbar\omega$ followed by emission at $\hbar\omega$). Although at large photon energies these processes follow an energy-dependence similar to that of 2D graphene, the zigzag nanoribbons exhibit a richer structure at low photon energies, arising from divergences of the joint density of states and from resonant absorption processes, which can be strongly modified by doping. As a figure of merit for the injected carrier currents, we calculate the resulting swarm velocities. Finally, we provide estimates for the limits of validity of our model.

Tables

Table 1. Table 1: Onset energies for the lowest energy transitions for an undoped cold ZGNR. Tuples ( n , m ) 𝑛 𝑚 (n,m) indicate a transition from band m 𝑚 m to band n 𝑛 n and every onset energy indicates the position of a JDOS divergence. The peak number is as indicated in Fig.  4 .

Peak	${\xi }^{xx}$		${\xi }^{yy}$
number	E (eV)	Transition	E (eV)	Transition
1	0.149	$(2,-1),(1,-2)$	0.000	$(1,-1)$
2	0.323	$(4,-1),(1,-4)$	0.236	$(3,-1),(1,-3)$
3	0.350	$(3,-2),(2,-3)$	0.410	$(5,-1),(1,-5)$
4	0.498	$(6,-1),(1,-6)$	0.439	$(4,-2),(2,-4)$
⋮	⋮	⋮	⋮	⋮

Table 2. Table 2: Velocity matrix elements at the Dirac point 𝐊 𝐊 \bm{\mathrm{K}} . At a given κ x subscript 𝜅 𝑥 \kappa_{x} , any of these matrix elements are purely real or purely imaginary (which is explicitly indicated by the presence (absence) of the imaginary unit i 𝑖 i ). The corresponding expressions at the other Dirac point 𝐊 ′ superscript 𝐊 bold-′ \bm{\mathrm{K}^{\prime}} are identical, except that the 𝒙 ^ − limit-from bold-^ 𝒙 \bm{\hat{x}}- components of the matrix elements flip sign; the 𝒚 ^ bold-^ 𝒚 \bm{\hat{y}} -components of the matrix elements remain unchanged. The range of validity for this expressions is given in the third column.

Type	Expression .	Conditions
nConf	$v_{nm}^x\left({\kappa }_x\right)=-4v_F\left({\zeta }_m^{\text{conf}}+{\zeta }_n^{\text{conf}}\right)A_n^{\text{conf}}{A}_m^{\text{conf}}\left[\frac{{𝒦}_m^{\mathrm{conf}}\sin \left({𝒦}_n^{\mathrm{conf}}W\right)-{𝒦}_n^{\mathrm{conf}}\sin \left({𝒦}_m^{\mathrm{conf}}W\right)}{{\left({𝒦}_m^{\mathrm{conf}}\right)}^2-{\left({𝒦}_n^{\mathrm{conf}}\right)}^2}\right]$	$\left|n\right|\ge 2,\left|m\right|\ge 2,\forall {\kappa }_x$, or
$\updownarrow$		, or
mConf	$v_{nm}^y\left({\kappa }_x\right)=-i\mathrm{\hspace{0.33em}4}{v}_F\left({\zeta }_m^{\text{conf}}-{\zeta }_n^{\text{conf}}\right)A_n^{\text{conf}}{A}_m^{\text{conf}}\left[\frac{{𝒦}_m^{\mathrm{conf}}\sin \left({𝒦}_n^{\mathrm{conf}}W\right)-{𝒦}_n^{\mathrm{conf}}\sin \left({𝒦}_m^{\mathrm{conf}}W\right)}{{\left({𝒦}_m^{\mathrm{conf}}\right)}^2-{\left({𝒦}_n^{\mathrm{conf}}\right)}^2}\right]$
nEdge	$v_{nm}^x\left({\kappa }_x\right)=-4v_F\left({\zeta }_m^{\text{edge}}+{\zeta }_n^{\text{edge}}\right)A_n^{\text{edge}}{A}_m^{\text{edge}}\left[\frac{{𝒦}_n^{\mathrm{edge}}\sinh \left({𝒦}_m^{\mathrm{edge}}W\right)-{𝒦}_m^{\mathrm{edge}}\sinh \left({𝒦}_n^{\mathrm{edge}}W\right)}{{\left({𝒦}_m^{\mathrm{edge}}\right)}^2-{\left({𝒦}_n^{\mathrm{edge}}\right)}^2}\right]$	$\left|n\right|\ge 1,\left|m\right|\ge 1,{\kappa }_x>W^{-1}$
$\updownarrow$
mEdge	$v_{nm}^y\left({\kappa }_x\right)=-i\mathrm{\hspace{0.33em}4}{v}_F\left({\zeta }_m^{\text{edge}}-{\zeta }_n^{\text{edge}}\right)A_n^{\text{edge}}{A}_m^{\text{edge}}\left[\frac{{𝒦}_n^{\mathrm{edge}}\sinh \left({𝒦}_m^{\mathrm{edge}}W\right)-z_m\sinh \left({𝒦}_n^{\mathrm{edge}}W\right)}{{\left({𝒦}_m^{\mathrm{edge}}\right)}^2-{\left({𝒦}_n^{\mathrm{edge}}\right)}^2}\right]$
nConf	$v_{nm}^x\left({\kappa }_x\right)=i4v_F\left({\zeta }_m^{\text{edge}}+{\zeta }_n^{\text{conf}}\right)A_n^{\text{conf}}{A}_m^{\text{edge}}\left[\frac{{𝒦}_n^{\mathrm{conf}}\sinh \left({𝒦}_m^{\mathrm{edge}}W\right)-{𝒦}_m^{\mathrm{edge}}\sin \left({𝒦}_n^{\mathrm{conf}}W\right)}{{\left({𝒦}_m^{\mathrm{edge}}\right)}^2+{\left({𝒦}_n^{\mathrm{conf}}\right)}^2}\right]$	$\left|n\right|\ge 2,\left|m\right|=1,{\kappa }_x>W^{-1}$
$\updownarrow$
mEdge	$v_{nm}^y\left({\kappa }_x\right)=-4v_F\left({\zeta }_m^{\text{edge}}-{\zeta }_n^{\text{conf}}\right)A_n^{\text{conf}}{A}_m^{\text{edge}}\left[\frac{{𝒦}_n^{\mathrm{conf}}\sinh \left({𝒦}_m^{\mathrm{edge}}W\right)-{𝒦}_m^{\mathrm{edge}}\sin \left({𝒦}_n^{\mathrm{conf}}W\right)}{{\left({𝒦}_m^{\mathrm{edge}}\right)}^2+{\left({𝒦}_n^{\mathrm{conf}}\right)}^2}\right]$
nEdge
$\updownarrow$	${(\text{Conf}\leftrightarrow \text{Edge})}^{†}$	$\left|n\right|=1,\left|m\right|\ge 2,{\kappa }_x>W^{-1}$
mConf

Equations

$$\Psi (\mathbf{r})$$

$$+\sum _{{\mathbf{R}}_B}{\psi }_B({\mathbf{R}}_B)\phi (\mathbf{r}-{\mathbf{R}}_B).$$

$${\psi }_A(\mathbf{r})$$

$${\psi }_B(\mathbf{r})$$

$$\left[\begin{array}{cccc}0& -i{\partial }_x-{\partial }_y& 0& 0\\ -i{\partial }_x+{\partial }_y& 0& 0& 0\\ 0& 0& 0& -i{\partial }_x+{\partial }_y\\ 0& 0& -i{\partial }_x-{\partial }_y& 0\end{array}\right]$$

$$\times \left[\begin{array}{c}{F}_A^{\mathbf{K}}(\mathbf{r})\\ F_B^{\mathbf{K}}(\mathbf{r})\\ F_A^{{\mathbf{K}}^{\prime }}(\mathbf{r})\\ F_B^{{\mathbf{K}}^{\prime }}(\mathbf{r})\end{array}\right]=\frac{E}{\gamma }\left[\begin{array}{c}{F}_A^{\mathbf{K}}(\mathbf{r})\\ F_B^{\mathbf{K}}(\mathbf{r})\\ F_A^{{\mathbf{K}}^{\prime }}(\mathbf{r})\\ F_B^{{\mathbf{K}}^{\prime }}(\mathbf{r})\end{array}\right],$$

$${\mathbf{F}}^{\mathbf{K}}(\mathbf{r})$$

$${\mathbf{F}}^{{\mathbf{K}}^{\prime }}(\mathbf{r})$$

$${\Phi }_B^{\mathbf{K}}(y=0)$$

$${\Phi }_B^{{\mathbf{K}}^{\prime }}(y=0)$$

$$\gamma \left[\begin{array}{cc}0& {\kappa }_x-{\partial }_y\\ {\kappa }_x+{\partial }_y& 0\end{array}\right]\left[\begin{array}{c}{\Phi }_A^{\mathbf{K}}(y)\\ {\Phi }_B^{\mathbf{K}}(y)\end{array}\right]=E\left[\begin{array}{c}{\Phi }_A^{\mathbf{K}}(y)\\ {\Phi }_B^{\mathbf{K}}(y)\end{array}\right].$$

$${\Phi }_A^{\mathbf{K}}(y)$$

$${\Phi }_B^{\mathbf{K}}(y)$$

$$e^{-2\mathcal{K}W}=\frac{{\kappa }_x-\mathcal{K}}{{\kappa }_x+\mathcal{K}},$$

$${\kappa }_x=\mathcal{K}\coth \left(W\mathcal{K}\right),$$

$${\Phi }_A^{\mathbf{K}}(y)$$

$${\Phi }_B^{\mathbf{K}}(y)$$

$${\zeta }_n^{\text{edge}}$$

$$A^{\text{edge}}$$

$$E_n^{\mathrm{edge}}=n\gamma \sqrt{{\kappa }_x^2-({\mathcal{K}}^{\mathrm{edge}}{)}^2}.$$

$${\kappa }_x={\mathcal{K}}_n\cot \left(W{\mathcal{K}}_n\right),$$

$${\Phi }_A^{\mathbf{K}}(y)$$

$${\Phi }_B^{\mathbf{K}}(y)$$

$${\zeta }_n^{\text{conf}}$$

$$A_n^{\text{conf}}$$

$$E_n^{\mathrm{conf}}$$

$$E_{\pm 1}(W^{-1})$$

$$E_{\pm n}(W^{-1})$$

$${\mathbf{v}}_{nm}({\kappa }_x)$$

$$H=\gamma \left[\begin{array}{cc}0& -i{\partial }_x-{\partial }_y\\ -i{\partial }_x+{\partial }_y& 0\end{array}\right],$$

$$v_{nm}^x({\kappa }_x)=0\mathrm{if}{\zeta }_n\ne {\zeta }_m,$$

$$v_{nm}^y({\kappa }_x)=0\mathrm{if}{\zeta }_n={\zeta }_m.$$

$$\mathbf{E}(t)=\mathbf{E}(\omega )e^{-i\omega t}+\mathbf{E}(2\omega )e^{-2i\omega t}+c.c.,$$

$$H_{\mathrm{int}}(t)=-e\mathbf{v}\cdot \mathbf{A}(t),$$

$$U(t)$$

$$U_{\mathrm{int}}(t)=$$

$$+(i\hslash {)}^{-2}{\int }_{-\infty }^t{V}_{\mathrm{int}}(t_1)d{t}_1{\int }_{-\infty }^{t_1}{V}_{\mathrm{int}}(t_2)d{t}_2+\dots$$

Full text

nanoribbons, optical coherent control, graphene

††preprint: APS/123-QED

I Introduction

The electronic properties of low-dimensional materials depend strongly on the size and geometry of the system Ogawa_optics_low_dimensions ; Nakada_edges_early . For instance, the bandstructure of a monolayer and a stripe of graphene are significantly different. A stripe of graphene is usually referred as a graphene nanoribbon, where the boundaries impose novel conditions on the wavefunctions; for a zigzag graphene nanoribbon (ZGNR), the wavefunction vanishes on a single sublattice, A or B, at each edge. As shown earlier Nakada_edges_early ; BreyFertig ; marconciniKDP , in ZGNR, there are confined states that extend across the width of the system, incorporating states from both sublattices. There is also another class of states strongly localized at each edge, which incorporate states from either one or the other sublattice; these states are known as edge states. Although confined states are also found in other types of ribbons, such as armchair, the edge states are present only for zigzag ribbons. These edge and confined states provide many of the novel characteristics seen in ZGNR. Moreover, the energy of these states can be easily tuned by changing the ribbon width, applying external fields, and functionalizing the system Lee_functionalization ; daRocha_functionalization . Since for an undoped ZGNR the Fermi level coincides with the flat part of the edge states, then tuning the doping level allows to easily control the contribution of the edge states. Given that a 2D graphene sheet lacks of these localized states, a ZGNR offers the advantage of having optical responses that are easily tuneable. Over the last years, a number of studies have reported the special properties of these localized states Nakada_edges_early ; BreyFertig ; marconciniKDP ; Yang_Cohen_Louie_2008 ; Kunstmann_Stability ; Bellec_edges_artificial_graphene ; Delplace_edges_zak_phase and recent investigations have described more novel properties and applications  Yazyev_applications ; luck2015 ; Fujita_EELS_2015 ; pelc-brey ; gonsalbez-spin-filtered ; Stegmann . At zero energy they have an important role in the electronic transport properties of both clean and disordered ZGNR, as Luck et al. luck2015 (and references therein) have recently shown using a tight-binding formalism with a transfer-matrix approach. A detailed review of these localized states in graphene-like systems can be found in Lado et al. lado2015 . The optical properties of ZGNR and graphene nano-flakes have been studied from a number of perspectives Hsu-Reich_selection_rules ; Sasaki_optical_transitions ; Yamamoto2006_TB_optics_flakes ; Berahman_optics_chiral ; zhu_optical_conduc ; Yang_Cohen_Louie_2008 ; Prezzi_2008_optics ; Fujita_EELS_2015 , always showing the strong influence of the edge states in the dielectric function. First-principles studies of functionalization in graphene ribbons have shown Lee_functionalization that the low-energy $\pi$ electrons at the edges of the ZGNR lead to higher binding energies as compared with ribbons of different shape edges. Similar studies indicate daRocha_functionalization that the optical response of functionalized ZGNR depends strongly on the size, shape and location of the deposited molecule, suggesting functionalization as an effective way of fine-tuning the electronic and optical properties of ZGNR.

In this work, we investigate the optical injection of carriers and currents in graphene nanoribbons by means of coherent light fields at $\omega$ and $2\omega$. In general, for arbitrary beams, this technique is referred as coherent current control. It is based on the fundamental feature that if the quantum evolution of a system can proceed via several pathways, then the interference between such pathways can play a determining role in the final state of the system ManykinAfanas ; Manykin . In a semiconductor, it is possible to control the injection of carriers ccontrolDrielSipe2001 ; sun-norris2010 ; Rioux2011 ; fregoso , spins, electrical current kiran2012 , spin current Muniz_SpinTopos , and even valley current Muniz_ValleyCurrent15 , using phase-dependent perturbations, usually involving coherent beams or pulses of light. In a one-color scheme, the interference is between transition amplitudes associated with different polarizations ccontrolDrielSipe2001 . Although carrier injection can be achieved with one-color excitation, current injection cannot. This is due to symmetry considerations, since one-color current injection is characterized by a third-rank tensor, hence it is only allowed in systems that lack inversion symmetry ccontrolDrielSipe2001 . Due to the inversion symmetry in zigzag graphene ribbons, the one-color coherent control process is forbidden. In a two-color scheme, the interference is between pathways related to photon absorption processes arising from different phase related beams, one at $\omega$ and the other at $2\omega$. In this case, current injection is characterized by a fourth-rank tensor, hence it is nonzero for a ZGNR. In both schemes, the different pathways connect the same initial and final states. Here our focus is on two-color current injection, and we consider two classes of processes: the first class arises from the interference of one-photon absorption at $2\hbar\omega$ with two-photon absorption at $\hbar\omega$, and the second class arises from the interference of one-photon absorption at $\hbar\omega$ with stimulated electronic Raman scattering at $\hbar\omega$. In general, coherent control injection allows for the placement of electrons and holes in different bands and portions of the Brillouin Zone as $\omega$ is varied. Thus, as we will show, the current injection is very sensitive to the presence of both confined and edge states. In line with plausible experiments, we consider nanoribbons with a width on the order of 20 nanometers, which leads to unit cells containing a few hundreds of atoms. For this reason, we employ an envelope function strategy to calculate the relevant energies and velocity matrix elements; the rest of the calculation follows a conventional Fermi’s golden rule approach to calculate the absorption coefficients.

The article is organized as follows. In Sec. II, we describe the model Hamiltonian employed to describe the wavefunctions, the resulting bandstructure, and the selection rules for the velocity matrix elements. In Sec. III, we describe the different carrier injection and current injection coefficients, including the conventional and Raman contributions. In Sec. IV, we revisit these calculations, but for a $p$-doped system. This allows us to show the significant change in the signals that can be accomplished by altering the occupation of the edge states. In Sec. V, we provide an estimate of the limits of validity of the model employed in this work. Finally, in Sec. VI, we present our final discussions and conclusions.

II Theoretical Model

II.1 Model Hamiltonian

A zigzag graphene nanoribbon (ZGNR) is a strip of monolayer graphene CastroNetoReview2009 ; DasSarmaReview2011 that has been cut such that the edges along its length have a zigzag shape, as shown in Fig. 1. We take the ribbon to lie in the ($xy$) plane, with $\hat{\bm{x}}$ as the longitudinal direction along which the ribbon extends over all space; $\hat{\bm{y}}$ then identifies the direction across the ribbon, along which the electron states are confined.

We assume passivated carbon atoms at the longitudinal boundaries, as if hydrogen atoms were adsorbed marconciniKDP ; Fujita_EELS_2015 ; this allows the passivation of any dangling edge states and the neutralization of the spin moments at the edges Fujita_EELS_2015 . We take $W=a\sqrt{3}\,(2N+2)/6$ as the effective width, where $N$ is the total number of atoms in the unit cell, $a=a_{\mathrm{cc}}\sqrt{3}=0.246~{}\mathrm{nm}$ is the graphene lattice constant, and $a_{\mathrm{cc}}$ is the carbon-carbon distance (see Fig. 1). The edge at $y=a/\sqrt{3}$ is formed by A-atoms, while the edge at $y=W-a/\sqrt{3}$ is formed by B-atoms. The lattice vector is $\bm{a}=a\hat{\bm{x}}$ and the atomic sites are set in terms of the graphene lattice vectors, $\bm{a}_{1}=(\hat{\bm{x}}-\sqrt{3}\hat{\bm{y}})\,a/2$ and $\bm{a}_{2}=(-\hat{\bm{x}}-\sqrt{3}\hat{\bm{y}})\,a/2$. The Dirac points of monolayer graphene are projected marconciniKDP into the one-dimensional Brillouin zone of the ZGNR, $[-\frac{\pi}{a},\frac{\pi}{a})$, as $\mathbf{K}=\frac{2\pi}{3a}$ and $\mathbf{K}^{\prime}=-\frac{2\pi}{3a}$.

We express the total wavefunctions as linear combinations of atomic orbitals $\varphi$ that are centered at atomic sites A and B,

[TABLE]

Then, following Marconcini and Macucci marconciniKDP , we employ the semi-empirical $\bm{k}\cdot\bm{p}$ method to describe $\Psi(\bm{r})$ with a smooth envelope function approach. The coefficients $\psi_{A}$ and $\psi_{B}$ in Eq. (II.1) can be written as

[TABLE]

where the $F_{A(B)}^{\bm{K}(\bm{K}^{\prime})}(\mathbf{r})$ are the envelope function components associated with the $\bm{K}(\bm{K}^{\prime})$ Dirac point and the orbital at atom A(B)111The graphene’s honeycomb lattice is composed by two distinct triangular lattices, A and B. On each sub-lattice all atoms are equivalent.. In writing Eq. (2) we have replaced $\psi_{i}(\bm{R}_{i})\to\psi_{i}(\bm{r})$ for $i=\{A,B\}$, on the basis of two assumptions. First, we assume that atomic orbitals are strongly localized at their corresponding atom, and second, we assume that the envelope functions are slow-varying functions of $\bm{r}$ near the $\mathbf{K}$ ($\mathbf{K}^{\prime}$) Dirac point. These envelope functions satisfy the Dirac equation,

[TABLE]

where $\gamma=(\sqrt{3}/2)\,ta$, $t=2.70~{}\mathrm{eV}$ is the nearest-neighbor hopping parameter and $v_{F}=\gamma\hbar^{-1}$ is the graphene Fermi velocity. Because of the translational symmetry along $\hat{\bm{x}}$, each envelope function can be factorized as the product of a propagating plane wave along the length direction ($\hat{\bm{x}}$), and a function confined along the width direction ($\hat{\bm{y}}$),

[TABLE]

where $\kappa_{x}$ ($\kappa_{x}^{\prime}$) is the wavevector along the length of the ribbon, measured from the Dirac point $\bm{K}$ ($\bm{K}^{\prime}$). The passivation of the carbon atoms at the edges terminates the $\pi$ orbitals thereat, thus it is reasonable to assume that the full wavefunction vanishes at the lattice sites located at the effective edges. This leads to the following boundary conditions for the confined part of the envelope functions marconciniKDP ,

[TABLE]

These boundary conditions and the block diagonal form of the matrix in Eq. (II.1) cause the envelope functions at $\bm{K}$ to be uncoupled from their counterparts at $\bm{K}^{\prime}$; therefore they can be studied separately. With the use of Eq. (4), the Dirac equation for the $\bm{K}$ valley is

[TABLE]

The solutions of Eq. (7) are of the form marconciniKDP ,

[TABLE]

where $\mathcal{K}=\sqrt{\kappa_{x}^{2}-(E/\gamma)^{2}}$. Under the boundary conditions (Eq. (6a)), this leads to a relation between the transverse ($\mathcal{K}$) and the longitudinal ($\kappa_{x}$) wavenumbers,

[TABLE]

which shows that they are coupled for ZGNR. If $\mathcal{K}$ is taken to be real, then Eq. (10) reduces to

[TABLE]

and without loss of generality we assume $\mathcal{K}$ to be positive. Equation (11) supports two eigensolutions for $\kappa_{x}>W^{-1}$, which we label as $n=1$ for positive energies and $n=-1$ for negative energies; both correspond to states strongly confined at the edges, henceforth referred as edge states marconciniKDP ,

[TABLE]

where $L$ is a normalization length along the $\hat{\bm{x}}$ direction. We have also set $\mathcal{K}\to\mathcal{K}^{\text{edge}}$, and $A^{\text{edge}}$ is the usual wavefunction normalization coefficient,

[TABLE]

and the eigenenergy is

[TABLE]

Conversely, if we consider solutions of Eq. (10) with $\mathcal{K}$ purely imaginary, of the form $i\mathcal{K}_{n}$ with $\mathcal{K}_{n}$ real, then Eq. (10) reduces to

[TABLE]

where, without loss of generality, we take $\mathcal{K}_{n}$ to be positive. These solutions give states that extend over the full width of the ribbon, and are known simply as confined states; for these we set $\mathcal{K}_{n}\to\mathcal{K}_{n}^{\mathrm{conf}}$ and label them by $n=\pm 1,\pm 2,\pm 3,\ldots$, starting with $\pm 1$ for those with energies closest to zero. These confined states exist for any real $\kappa_{x}$, except those with band index $n=\pm 1$, which exist only for $\kappa_{x}{\leq}W^{-1}$. The dispersion relations of the confined states with band index $n=\pm 1$ connect with that of the edge states; both share the band index $n=\pm 1$ (transition from the red to the blue traces in Fig. 3). The confined states have the form

[TABLE]

where

[TABLE]

We can indicate any of the edge or confined states simply by $|n\kappa_{x}\rangle$, where if $|n|\geq 2$ the state is confined, while if $|n|=1$ then the state is confined for $\kappa_{x}{\leq}W^{-1}$, but it is an edge state if $\kappa_{x}>W^{-1}$.

Equations (16) and (22) describe the bandstructure of ZGNR, shown in Figs. 2 and 3. The edge states are flattened towards the zero energy level for $\kappa_{x}>W^{-1}$ (Fig. 3), whereas the confined states have a parabolic structure around the Dirac points, with an axis of symmetry at $\kappa_{x}=W^{-1}$, except for the two confined states nearest to zero energy, with band index $n=\pm 1$ and $\kappa_{x}\leq W^{-1}$ (Fig. 3). These confined states are associated with the Dirac cones of 2D graphene. Since the extrema of the confined states occur at $\kappa_{x}=W^{-1}$ , we can express the band energies at such value of $\kappa_{x}$ as

[TABLE]

for the edge and confined states, respectively. This indicates that the band gap scales as $W^{-1}$ and provides an estimate of the photon energy at which the absorption edge occurs with respect to the ribbon width $W$. It turns out that the sign functions appearing in the expressions for $\Phi_{A}^{\mathbf{K}}(y)$ [Eq. (12) for edge states and Eq. (18) for confined states] alternate for consecutive states, being $+1$ for the first state above zero energy, $-1$ for the next up, and so on; the situation is reversed for negative energies. This sign factor plays an important role in the selection rules of the quantities we calculate. Therefore we indicate these sign factors on the bandstructure diagram [Figs. 2 and 3): a solid line indicates that the confined part of an A-site component of the envelope function has $\zeta_{n}=+1$, whereas a dashed trace means it has $\zeta_{n}=-1$.

II.2 Velocity matrix elements

We employ the envelope functions given by Eq. (4) in order to calculate the velocity matrix elements (VME) that describe the coupling between two states $|n,\kappa_{x}\rangle$ and $|m,\kappa_{x}\rangle$ as,

[TABLE]

where $\kappa_{x}$ is a wavenumber and $n$, $m$ are band indices. The velocity operator is given by $\bm{v}=\left[\bm{r},H\right]/(i\hbar)$, which, together with the Hamiltonian in Eq. (7) for the $\mathbf{K}$ valley,

[TABLE]

leads to $\bm{v}=v_{F}(\sigma_{x},\sigma_{y})$, where $\sigma_{x}$ and $\sigma_{y}$ are the Pauli matrices and $v_{F}=\gamma/\hbar$ is graphene’s Fermi velocity. The resulting expressions are given in Appendix A, Table 2, and obey the following selection rules:

[TABLE]

We close this section by mentioning that the solutions corresponding to the Dirac point $\bm{K}^{\prime}$ are analogous to those presented here for $\bm{K}$. As shown by Marconcini et al. marconciniKDP , the wavefunctions for the A sites, Eqs. (12) and (18), at the $\bm{K}^{\prime}$ differ by a sign factor from those at $\bm{K}$. Moreover, the velocity operator at the $\bm{K}^{\prime}$ has the form $\bm{v}=v_{F}(\sigma_{x},-\sigma_{y})$. This, together with the properties of the envelope functions at both valleys, causes the $x$ component of the VME at $\mathbf{K}^{\prime}$ to have opposite sign of those at $\mathbf{K}$; the $y$ components of the VME are the same near $\mathbf{K}$ as near $\mathbf{K}^{\prime}$.

III Coherent injection and control

III.1 Framework

In this section, we describe the general framework of the two-color coherent control scheme. As mentioned in the Introduction, the quantum interference is between pathways associated with photon absorption processes arising from different phase related beams. These pathways connect the same initial and final states, as shown for the processes in Fig. 3, where we consider the two-color scheme with beams at $\omega$ and at $2\omega$. This figure depicts the two classes of processes we study in this paper.

The first, conventional processes, are those where current injection arises due to the interference of one-photon absorption (OPA) at $2\hbar\omega$ and two-photon absorption (TPA) of (two) photons with energy $\hbar\omega$ ccontrolDrielSipe2001 ; this is depicted with the set of arrows on the right of Fig. 3, under the label “CONV”. In the remaining of the discussion, we label variables associated with conventional processes with a subindex ‘C’.

The second class of processes arise in experiments on narrow band gap or gapless materials, with $\hbar\omega>E_{g}$, where $E_{g}$ is the energy band gap. Under this condition, current injection can arise due to the interference of OPA at $\hbar\omega$ and stimulated electronic Raman scattering (ERS) at $\hbar\omega$ RiouxERS . This ERS is indicated by the set of arrows at $2\hbar\omega$ and $\hbar\omega$ in the left of Fig. 3, under the label “ERS”. We refer to variables associated with this Raman processes with a subindex ‘R’. We mention that in coherent control experiments on typical semiconductors, the beam frequencies employed are such that $\hbar\omega<E_{g}<2\hbar\omega$, and, consequently, the ERS current is absent because OPA at $\hbar\omega$ is impossible.

Following van Driel and Sipe ccontrolDrielSipe2001 ; ccontrolDrielSipe2005 , we calculate the one- and two-photon carrier injection and current injection rates due to the interaction with a classical electromagnetic field

[TABLE]

in the long wavelength limit, where $\omega$ is the fundamental frequency. The interaction between the electric field and the electron system is accounted by the minimal coupling prescription in the Hamiltonian of Eq. (25); we do the usual replacement $p_{j}\to p_{j}-eA_{j}(t)$, for $j=(x,y)$, with $p_{j}=-i\hbar\partial_{j}$, and obtain the interaction Hamiltonian that acts as the perturbation,

[TABLE]

where $e=-|e|$ is the electron charge and $\bm{A}(t)$ is the vector potential associated with the electric field, $\bm{E}(t)=-\partial\bm{A}(t)/\partial t$. We treat this problem using standard time-dependent perturbation theory and Fermi’s golden rule. Since we are interested in OPA, TPA and ERS processes, the unitary evolution operator $U(t)$ is expanded perturbatively up to second order,

[TABLE]

where

[TABLE]

and

[TABLE]

Under the perturbation of Eq. (28), the evolution of the system’s state \big{|}\Upsilon\rangle is not just the ground state \big{|}0\rangle, but it also contains an amplitude of the excited state $|nm\bm{\kappa_{x}}\rangle$ (this ket corresponds to a state with an electron-hole pair),

[TABLE]

where \big{|}c_{nm\bm{\kappa}_{x}}(t)\big{|}^{2} is the probability that the system is at \big{|}nm\bm{\kappa_{x}}\rangle; the missing terms in Eq. (32) correspond to higher order excitations, which we neglect in this work. The carrier injection and the current injection rates are given by

[TABLE]

respectively, where $L$ is the normalization length introduced below Eq. (14). To describe the optical processes we are interested, we compute $c_{nm\bm{\kappa}_{x}}(t)$ up to second order (a tutorial derivation can be found in Ref. ccontrolDrielSipe2001 ; see also Ref. RiouxERS ). Then, the macroscopic expressions for these injection rates get the form,

[TABLE]

where repeated indexes indicate summation, $\omega$ is the fundamental frequency, $\dot{n}^{(1)}$ and $\dot{n}_{\mathrm{C}(\mathrm{R})}^{(2)}$ account for the first- and second-order absorption processes, respectively; overall $\dot{n}$ refers to the rate of injected carriers per unit length along the ribbon (carriers per unit length per unit time). The OPA coefficient is described by a second-order tensor, $\xi^{ab}$, while the TPA and the ERS absorption coefficients are described by fourth-order tensors, $\xi^{abcd}_{\mathrm{C}}$ and $\xi^{abcd}_{\mathrm{R}}$, respectively. Here, $\dot{J}^{a}$ includes the electron and hole contributions to the current (charge per unit time), injected per unit time along the ribbon. The current injection coefficient $\eta(\omega)$ in Eq. (38) includes the conventional and the ERS contributions, i.e., $\eta(\omega)=\eta_{C}(\omega)+\eta_{R}(\omega)$. In the following sections, we give the full expressions for these coefficients. Note that the coefficients can be chosen such that $\xi^{abcd}_{C}=\xi^{bacd}_{C}=\xi^{badc}_{C}$ and $\eta^{abcd}=\eta^{acbd}$.

III.2 First-order absorption process

We calculate the expressions for the coefficients $\xi$ and $\eta$ appearing in Eq. (35)–(38) using Fermi’s golden rule. For the one-photon absorption coefficient, we obtain

[TABLE]

where we have gone from a sum over states to an integral over reciprocal space by $L^{-1}\sum_{\kappa_{x}}\to(2\pi)^{-1}\int{d\kappa_{x}}$. In this expression the sum $\sum_{nm}$ runs over all bands, filled and empty (similarly for the other response functions considered here); $\omega_{nm}(\kappa_{x})=\hbar^{-1}E_{nm}(\kappa_{x})$ and $E_{nm}(\kappa_{x})=E_{n}(\kappa_{x})-E_{m}(\kappa_{x})$ is the energy difference between two states at a given $\kappa_{x}$. A factor of two has been included to account for spin degeneracy, which we do throughout the paper. The $x-$components of the VME at the $\bm{K}$ and $\bm{K}^{\prime}$ valleys differ just by a sign while the $y-$components of the VME are the same. Consequently, since all integrals over reciprocal space include pairs of VME, the integration over $\kappa_{x}$ can be restricted to a single valley, $\bm{K}$, and another factor of two included to account for the contribution of the $\bm{K}^{\prime}$ valley.

The occupation of the states is described by the Fermi-Dirac distribution. In all of our integrals over reciprocal space $f_{mn}(\kappa_{x})=f_{m}(\kappa_{x})-f_{n}(\kappa_{x})$, with $f_{n}(\kappa_{x})=[1+e^{(E_{n}(\kappa_{x})-\mu)/(k_{B}T)}]^{-1}$ at temperature $T$ and chemical potential $\mu$. Until the end of Sec. IV, we confine ourselves to zero temperature, hence $f_{n}(\kappa_{x})=\theta(E_{n}(\kappa_{x})-\mu)$, where $\theta$ is the Heaviside step function. Because of the selection rules for the VME, Eq. (26), the only nonzero components of the one-photon coefficient are $\xi^{xx}$ and $\xi^{yy}$, which we plot in Fig. 4 for a system at zero chemical potential. As a comparison dropoff , we include plots of $W\xi^{xx}_{2D}$, where $W$ is the effective width of the ribbon,

[TABLE]

and $\xi_{\mathrm{2D}}^{xx}=\xi_{\mathrm{2D}}^{yy}$ is the OPA coefficient for a 2D monolayer of grapheneRioux2011 ; here $\sigma_{0}=g_{s}g_{v}e^{2}/(16\hbar)$ is the universal optical conductivity of graphene, and $g_{s}=2$, $g_{v}=2$ are the spin and valley degeneracies, respectively. For ZGNR, the main difference between the two OPA coefficients is that $\xi^{yy}$ diverges at zero photon energy, due to a divergence in the joint density of states (JDOS) between bands $n=+1$ and $n=-1$. In contrast, for such a pair of bands $\xi^{xx}$ is identically zero, due to the VME selection rules. For an undoped ZGNR, $\xi^{xx}$ displays its first divergence at about 0.15 eV, which is the value of the band gap at zero Fermi level, and corresponds to the onset of the transitions $(2,-1)$ and $(1,-2)$ at $\kappa_{x}=W^{-1}$; these four states give the initiation energy for $\xi^{xx}$. In the following we indicate a transition from band $m$ to band $n$ by $(n,m)$; hence, for zero chemical potential and zero temperature, the possible transitions have $m\leq-1$ and $n\geq 1$. In general, the $\xi^{xx}$ and $\xi^{yy}$ OPA coefficients possess an infinite number of divergences that arise due to the infinite number of parabolic bands in the bandstructure. Indeed, the JDOS between states with band index $n$ and $m$,

[TABLE]

can be shown to diverge as $(E-E_{nm}^{\mathrm{gap}})^{-1/2}$ for the confined states, and as $E^{-1}$ for edge states, where $E$ is the photon energy and $E_{nm}^{\mathrm{gap}}$ is the energy band gap between bands $n$ and $m$. In frequency space, these divergences occur at photon energies $E$ such that $E=E_{nm}^{\mathrm{gap}}$; in reciprocal space, they occur at $\kappa_{x}$ points where argument of the delta function has a zero derivative. The absorption coefficients inherit these JDOS divergences if the associated velocity matrix elements are nonzero at the $\kappa_{x}$ at which $dE_{nm}/d\kappa_{x}=0$. The sensitivity of an experiment to these divergences would depend on the resolution of the photon energy and on the magnitude of the velocity matrix elements, as well as on the presence of scattering effects that are not included in this simple treatment. In every pertaining Figure, we signal the location of these JDOS divergences by small green ticks. An interesting characteristic of $\xi^{xx}$ and $\xi^{yy}$ is that the divergence at the initiation energy always involves an edge state (see Table 1); this is reasonable, as these states are involved in the minimum band gap for an undoped system.

As mentioned above, the sum over states runs over all bands, filled and empty, but for a given photon energy range (e.g. $0-0.5~{}e$V, as in Fig. 4) the sum requires a finite number of bands. We refer to this as the “full” response. In order to highlight the contribution of the edge states, we also compute the response coefficients with a restricted sum over states $\sum_{nm}$, such that $n$ or $m$ are $\pm 1$, e.g., $(n,m)=\{(1,-1),(1,-2),(2,-1),\ldots\}$; we refer to this as the “edge” contribution and in the appropriate figures we plot it with black-dashed lines. This allows us to easily identify the contribution to OPA from states at bands $\pm 1$. At low photon energies such contribution is dominant: for $\xi^{xx}$, all transitions at photon energies $\hbar\omega<0.350~{}$eV are from or to edge states; for $\xi^{yy}$, all transitions at photon energies $\hbar\omega<0.439~{}$eV are from or to edge states. Consequently, at low-photon energies the “full” and “edge” contributions are indistinguishable. This is shown in Fig. 4 (see also Table 1), where for comparison we also plot $W\xi_{2D}^{xx}$, where $\xi_{2D}^{xx}$ is the OPA coefficient of graphene calculated Rioux2011 at the same level of approximation adopted here; it is clear how the presence of the edge states in ZGNR significantly modifies the OPA. Finally, we mention that the Dirac delta functions appearing in all our expressions are treated with an interpolation scheme interpNote .

III.3 Second-order absorption processes

III.3.1 Conventional process

In this section, we start by considering the second order process related to the absorption of two photons of energy $\hbar\omega$, indicated by the rightmost arrows in Fig. 3. Carrying the perturbation calculation up the second order, we obtain the two-photon absorption (TPA) coefficient,

[TABLE]

where

[TABLE]

which we regard as the effective velocity matrix element (effective VME) for the second order conventional process (C) process. Here $\beta_{C}$ is a small constant introduced to broaden resonant processes (discussed below) and the sum over $\ell$ corresponds to the virtual electron and virtual hole contributions ccontrolDrielSipe2001 . Although this sum runs over all bands (filled and empty), a converged value is obtained for $\ell=20$ bands for a photon energy range of 0–1 $e$V. From the selection rules for the regular VME, Eq. ((26)), we obtain the selection rules for $\mathcal{V}$,

[TABLE]

and from this we identify eight nonzero $\xi^{abcd}_{C}$ components, four of them independent, namely $\xi^{xxxx}_{C},\xi^{xxyy}_{C}=\left(\xi^{yyxx}_{C}\right)^{*},\xi^{xyxy}_{C}=\xi^{xyyx}_{C}=\xi^{yxxy}_{C}=\xi^{yxyx}_{C},$ and $\xi^{yyyy}_{C}$, which we show in Fig. 5. A feature of these coefficients is that the onset of the two-photon absorption signal is at the minimum band gap between bands $(2,-1)$, except for $\xi^{xyxy}_{C}$, which has its onset at 0 eV; this follows from the selection rules for the effective VME, which are inherited from the usual VME, and indicate that the transition $(1,-1)$ is allowed.

As we found for the OPA coefficients $\xi^{ab}$, the TPA coefficients $\xi^{abcd}_{C}$ suffer from divergences, but for the TPA coefficients they are of two types: JDOS divergences and effective-VME-divergences. The latter results when the nominal virtual state lies at the average of the energies between two transition states, $|n\kappa_{x}\rangle$ and $|m\kappa_{x}\rangle$, i.e., when (see Eq. (43))

[TABLE]

Such condition corresponds to a resonant TPA and an instance where this occurs is indicated on Fig. 3 by the three dots along the vertical line at $\kappa_{x}=0$. In Fig. 5 we distinguish these two types of divergences by small vertical lines of different color; a green tick indicates the presence of a JDOS-divergence, while a red tick indicates the presence of an effective-VME-divergence. In order to broaden the latter resonances, a small damping constant $\beta_{C}$ of 20 meV was introduced in the denominator of Eq. (43). This value, which is close to the thermal energy $k_{B}T$ associated with room temperature, was chosen arbitrarily. A more detailed theory would be necessary to indicate how these resonances are really broadened; the choice we make here simply allows us to identify easily where these resonances occur in our calculations. We mention that the onset of $\xi^{xxxx}_{C}$ is due to the transitions $(2,-1)$ and $(1,-2)$, which are free from resonances because the matrix elements to the intermediate states (one of the edge bands $\pm 1$ that would lead to a divergent condition) are forbidden by the selection rules. Therefore, in the photon energy range 0 to 0.15 eV, the coefficient $\xi^{xxxx}_{C}$ is free of resonances.

We present the $\xi^{abcd}_{C}$ coefficients in Fig. 5, and identify the edge contributions to them (black-dashed lines). As we found for $\xi^{ab}$, for $\xi^{abcd}_{C}$ the edge states make a dominant contribution at low photon energies, and are involved at the onset of TPA. As a comparison dropoff , in Fig. 5, we include plots of $W\xi^{abcd}_{2D}$, where $W$ is the effective width of the ribbon,

[TABLE]

and $\xi_{\mathrm{2D}}^{xxxx}=\xi_{\mathrm{2D}}^{yyyy}=\xi_{\mathrm{2D}}^{xyxy}=\xi_{\mathrm{2D}}^{xyyx}=-\xi_{\mathrm{2D}}^{xxyy}$ are the TPA coefficients for a 2D monolayer of graphene Rioux2011 ; as before, $g_{s}=2$ and $g_{v}=2$ are the spin and valley degeneracies, respectively.

III.3.2 ERS process

Now we consider another second order process involving light at $2\omega$ and light at $\omega$, stimulated electronic Raman scattering, which can be characterized as virtual absorption at $2\hbar\omega$ followed by emission at $\hbar\omega$; see the left diagram in Fig. 3. This process exists in semiconductors when the fundamental photon energy is larger than the band gap, which is always the case for an undoped ZGNR, because the edge states provide a zero-gap system. Following an earlier treatment of graphene RiouxERS , we find the ERS carrier injection to be

[TABLE]

where the effective VME for the ERS process are

[TABLE]

As in Eq. (43), $\beta_{\mathrm{R}}$ is a small constant introduced to broaden resonant processes and the sum over $\ell$ runs over all bands (filled and empty), but a converged value is obtained for $\ell=30$ bands for a photon energy range of 0–1 $e$V. The first term in the sum of Eq. (III.3.2) corresponds to photo-emission by an electron, and the second to photo-emission by a hole RiouxERS . Note that due to the different frequencies involved in Eq. (37), symmetrization of $\mathcal{V}_{\mathrm{R}}^{\mathrm{ij}}$ is unnecessary. The selection rules for $\mathcal{V}_{\mathrm{R}}^{\mathrm{ij}}$ are the same as those for $\mathcal{V}_{\mathrm{C}}^{\mathrm{ij}}$ (Eq. (44)); note, however, that $\mathcal{V}_{\mathrm{R}}^{\mathrm{ij}}\not=\mathcal{V}_{\mathrm{R}}^{\mathrm{ji}}$, although $\mathcal{V}_{\mathrm{R}}^{\mathrm{ij}}$ and $\mathcal{V}_{\mathrm{R}}^{\mathrm{ji}}$ satisfy the same selection rule. From this we identify six nonzero terms for the ERS carrier injection coefficient, $\xi_{\mathrm{R}}^{xxxx}$, $\xi_{\mathrm{R}}^{xyyx}=\left(\xi_{\mathrm{R}}^{yxxy}\right)^{*}$, $\xi_{\mathrm{R}}^{xxyy}=\left(\xi_{\mathrm{R}}^{yyxx}\right)^{*}$, $\xi_{\mathrm{R}}^{xyxy}$, $\xi_{\mathrm{R}}^{yxyx}$, and $\xi_{\mathrm{R}}^{yyyy}$. As do the conventional coefficients, the ERS coefficients suffer from JDOS and effective-VME divergences, the later arising whenever

[TABLE]

is satisfied. These conditions correspond to resonant processes, when a state is located at an energy $|E_{nm}(\kappa_{x})|$ above (below) the final (initial) state $n$ ($m$). As in Eq. (43), a small damping constant $\beta_{R}$ of 20 meV was introduced in the denominators of Eq. (III.3.2). All of these ERS coefficients present a large number of these resonances, causing $\xi^{abcd}_{\mathrm{R}}$ to be highly sensitive to the value of the $\beta_{\mathrm{R}}$ parameter. However, these resonances are of small magnitude for the energy range chosen for Fig. 6, hence they are not apparent. As shown, three of these components have their onset at zero photon energy, because the symmetry properties of the involved matrix elements allow for transitions between the two edge states.

III.4 Current Injection

III.4.1 Injection coefficients

We begin with the expression for $\eta_{C}$, the current injection coefficient characterizing the conventional process. Here the interference between the TPA at $\hbar\omega$ with OPA at $2\hbar\omega$ (see the right diagram in Fig. 3) leads to net current injection coefficients (including electron and hole contributions) given by ccontrolDrielSipe2001

[TABLE]

From the selection rules for the regular and the effective VME, Eq. (26) and Eq. (44), we identify three nonzero current injection coefficients, $\eta_{\mathrm{C}}^{xxxx},\eta_{\mathrm{C}}^{xyyx},$ and $\eta_{\mathrm{C}}^{xxyy}=\eta_{\mathrm{C}}^{xyxy}$. Notice that for all these tensors the first Cartesian component is $x$: Due to the confinement of the ribbons along the $\hat{y}$ direction (see Fig. 1), the current injection can only flow along the $\hat{x}$ direction, and all tensor components $\eta_{\mathrm{C}}^{yabc}$ are zero.

Turning to the expression for $\eta_{R}$, the current injection coefficient characterizing the interference between the ERS discussed above and the OPA at $\omega$ (see the left diagram in Fig. 3), including both electron and hole contributions we find

[TABLE]

where $\mathcal{V}_{\mathrm{R}}$ is given by Eq. (III.3.2). On the basis of the matrix elements selection rules, we identify three nonzero ERS current injection coefficients, $\eta_{\mathrm{R}}^{xxxx}$, $\eta_{\mathrm{R}}^{xyyx}$, and $\eta_{\mathrm{R}}^{xxyy}=\eta_{\mathrm{R}}^{xyxy}$.

Over the frequency range shown in Fig. 7, the conventional and the ERS current injection coefficients are of the same order, dropping off as the inverse of the third power of the photon energy, as do the coefficients for graphene RiouxERS . Thus we only plot the total injection coefficients $\eta^{abcd}=\eta_{\mathrm{C}}(\omega)+\eta_{\mathrm{R}}(\omega)$. For comparison, we include plots of $W\eta^{abcd}_{2D}$ (with the respective values of the Cartesian indices), where

[TABLE]

and $\eta_{\mathrm{2D}}^{xxxx}=3\eta_{\mathrm{2D}}^{xxyy}=3\eta_{\mathrm{2D}}^{xyyx}$ are the net current injection coefficients for a 2D monolayer of grapheneRiouxERS ; as before, $g_{s}=2$ and $g_{v}=2$ are the spin and valley degeneracies, respectively. As we saw for carrier injection, the edge states provide the strongest contribution at the onset of current injection. Another characteristic of these coefficients is that $\eta^{xxxx}$ has its onset at the band gap between bands $(2,-1)$, while ${\eta}^{xxyy}$ and ${\eta}^{xyyx}$ have their onset at 0 eV. This is due to the selection rules that the matrix elements involved in both the conventional and ERS process satisfy, allowing transitions between bands $(1,-1)$. An important characteristic of the current injection coefficients is that they are free of JDOS divergences, because the diagonal matrix elements in their respective expressions, Eqs. (50) and (51), are identically zero at the $\kappa_{x}$ at which the minimum gap occurs. However, a number of effective VME resonances do exist at photon energies indicated by the small red ticks in Fig.7, such that Eq. (45) is satisfied. As explained before, the magnitude of these resonances is broadened by a small damping constant. These coefficients are shown in Fig. 7, where we present the net current injection arising from the addition of the conventional and ERS contributions, i.e., $\eta(\omega)=\eta_{\mathrm{C}}(\omega)+\eta_{\mathrm{R}}(\omega)$.

III.4.2 Swarm velocities

The numerical values of the coefficients $\xi^{ab}$, $\xi^{abcd}_{\mathrm{C}(\mathrm{R})}$, and $\eta^{abcd}_{\mathrm{C}(\mathrm{R})}$ do not immediately give a sense of the average velocities with which the electrons and holes are injected. Sometimes an average, or swarm velocity is introduced to indicate this ccontrolDrielSipe2001 . In the system considered here, we could introduce a swarm velocity for both the conventional and ERS processes, according to

[TABLE]

where $\Omega=2\omega$ for $\bm{\mathscr{V}}_{\mathrm{C}}$ because $\bm{\dot{J}}_{\mathrm{C}}$ arises from the interference of OPA at $2\omega$ with TPA at $\omega$, while $\Omega=\omega$ for $\bm{\mathscr{V}}_{\mathrm{R}}$ because $\bm{\dot{J}}_{\mathrm{R}}$ arises from the interference of OPA at $\omega$ with the ERS described above. Besides describing an average speed that characterizes the injected carriers, one can consider maximizing Eq. (53) by using appropriate phases in the optical beams, and adjusting the relative amplitudes of the light at $\omega$ and $2\omega$. Considering just the swarm velocity of the conventional process, such optimization leads to equal OPA and TPA, and it follows that the intensity of the fundamental beam at $\omega$ should be about half an order of magnitude larger that of the beam at $2\omega$, for a fundamental photon energy of about $0.4~{}e$V. In contrast, the swarm velocity of the ERS process depends only on the intensity of the beam at $2\omega$. Further, in trying to optimize the net swarm velocity, determined by the total current injected divided by the total carrier density injected, one finds that the beam at $2\omega$ should have an intensity about an order of magnitude larger than the beam at $\omega$. Since in typical experiments the beam at $2\omega$ is obtained by second harmonic generation of part of the beam at $\omega$, this would be impractical. Thus we calculate the conventional and Raman swarm velocities for typical sun-norris2010 beam intensities of the fundamental and second harmonic fields, shown in Fig. 8. We complement these carrier velocities with the total average velocity of the injected carriers

[TABLE]

also evaluated at typical sun-norris2010 beam intensities. These carrier velocities are shown in Fig. 8. As a reference, at the photon energy of 0.25 eV, the maximum swarm velocity of the conventional process for a monolayer of graphene is $2.9\times 10^{5}$ ms*-1*. Hence the carrier velocities in ZGNR are comparable to those on a monolayer of graphene, as might be expected.

IV Doping

In the previous sections, we investigated the carrier and current injection at zero chemical potential. Since the dispersion relations of the edge states in ZGNR have a zero band gap and are flattened for $\kappa_{x}>W^{-1}$ (Fig. 3), those states are always involved at the onset energy of all of the optical response coefficients studied here. This suggests that doping is an effective method to alter the population of these two bands and the current that can be injected by the optical transitions between them. In this section, we revisit the calculations of $\xi^{ab}$, $\xi^{abcd}_{C(R)}$ and $\eta^{abcd}$ for a negative chemical potential, corresponding to a $p$-doped system. Besides the modified contribution from the edge states, we will also see significant modification in the contributions from other bands, particularly in the region near the $\mathbf{K}$ and $\mathbf{K}^{\prime}$ points, where doping leads to either a “valley” of filled states ($n$-doped), or a “hill” of unfilled states ($p$-doped); see Fig. 3.

We consider two negative Fermi levels, $\mu_{1}=-0.1$ eV and $\mu_{2}=-0.2$ eV, which in Fig. 3 we indicate by the upper boundaries of the grey areas. The value of $-0.1$ eV is interesting because, at this chemical potential, the flat part of band $-1$ (i.e., the region where $\kappa_{x}>W^{-1}$, cf. Fig. 3) contains empty states; this condition allows transitions from lower energy bands with final states in band $-1$, but also disables transitions from band $-1$ to upper bands. The second value, $\mu=-0.2$ eV, is interesting because at this potential a “hill” of unfilled states arises in the first parabolic band (band $-2$ in Fig. 3) at energies below our nominal value of zero.

We present the results of the calculations of OPA coefficients for those values of the chemical potential in Fig. 9. In an undoped sample, the JDOS divergences in $\xi^{xx}$ at low photon energies are due to the onset of the transitions $(2,-1),(1,-2),(4,-1),$ and $(1,-4)$ (see Table 1 and Fig. 4). Since all of these transitions involve bands $\pm 1$, any nonzero chemical potential has the capacity to significantly alter the OPA at these photon energies. For instance, if the Fermi level is at $-0.1$ eV, then the flat part of band $-1$ contains empty states, and the low photon energy divergences are removed. In addition, at this chemical potential transitions of the type $(-1,n)$, for $n$ odd and $<-1$ are permitted. However, the contributions to the OPA from these new transitions are of smaller magnitude than the contribution from the $(1,-2)$ transition, which is unaffected by the $-0.1$ eV doping. For this reason, the $(1,-2)$ transition remains as the main contribution to the $\xi^{xx}$ coefficient at low photon energies at this chemical potential (see Fig. 9).

At the Fermi level $-0.2$ eV, the edge states are completely empty, as are the states at the higher points of band $-2$ near the $\mathbf{K}$ and $\mathbf{K}^{\prime}$ points. This condition allows transitions of the type $(-2,n)$, for $n$ even and $<-2$, and also forbids transitions of the type $(n,-2)$, for $n$ odd and $\geq 1$, and $\kappa_{x}$ near the $\mathbf{K}$ and $\mathbf{K}^{\prime}$ points. It is this latter restriction which significantly changes the $\xi^{xx}$ coefficient near its onset. A further decrease in the Fermi level would consistently remove the divergences in $\xi^{xx}$ at low photon energies. All these observations were confirmed with a band-by-band calculation of $\xi^{xx}$.

The effect of doping the system has a larger influence on the onset energy of $\xi^{yy}$ that on that of $\xi^{xx}$. This is because the JDOS divergences at low photon energies relevant for $\xi^{yy}$ are due to the transitions $(1,-1)$, $(3,-1),$ and $(1,-3)$ (cf. Table 1). Therefore, even for small doping, the large contribution coming from the transitions between the two edge states (bands $\pm 1$, $\kappa_{x}>W^{-1}$) is significantly decreased, and leads to a greater change of the magnitude of $\xi^{yy}$ than of the magnitude of $\xi^{xx}$. A special signature of $\xi^{yy}$ for $\mu=-0.2~{}e$V (dark-violet signal, Fig. 9 b)) is the presence of two narrow peaks at 0.045 and 0.075 eV; the first of these peaks is due to the $(-1,-2)$ transition, while the second is from the $(-2,-3)$ transition. These two transitions are active only for those $\kappa_{x}$ states at which the “hill” of band $-2$ is empty (see Fig. 3). Notably, the transition $(-2,-3)$ brings a new JDOS divergence because it is active over a range of reciprocal space that includes $\kappa_{x}=W^{-1}$, where both bands have their maximum and their energy difference $E_{nm}(\kappa_{x})$ has a zero derivative (see the discussion below Eq. (41)).

In general, all these new transitions involve more JDOS divergences if the range of $\kappa_{x}$ over which they are active includes the $\kappa_{x}$ at which the band pairs have their maxima or minima. For instance, the divergences 1–4 in Fig. 9 are the same as those in Fig. 4 and Table 1, but the divergences 5–6 arise due to the new transitions allowed at nonzero chemical potentials: in Fig. 9 a), at the chemical potential $-0.20~{}e$V, the divergence 5 at $0.179$ eV is due to the transition $(-2,-4)$, which is active over a range of $\kappa_{x}$ that includes the $\kappa_{x}$ at which bands $(-2,-4)$ have their maxima, hence a new JDOS divergence appears. Likewise for $\xi^{yy}$ in Fig. 9 b) at $\mu=-0.20~{}e$V: divergences 5 at $0.089$ eV and 6 at $0.268$ eV exist because the transitions $(-2,-3)$ and $(-2,-5)$ are active over regions of reciprocal space that include the $\kappa_{x}$ at which such bands have their maxima.

In Fig. 10, 11, and 12 we present the nonzero $\xi^{abcd}_{\mathrm{C}}$, $\xi^{abcd}_{\mathrm{R}}$, and $\eta^{abcd}$ coefficients for selected nonzero Fermi levels. As was seen for $\xi^{ab}$, doping the ZGNR has the effect of modifying the responses around their onset energy, either due to the removal of some transitions, or due to the appearance of new ones, which in the undoped system were forbidden because the initial and final states were filled [e.g. ($-1,-2$) or ($-1,-3$)]. This shows that doping is an effective way of modifying the carrier and current injection in ZGNR, where the most significant changes are due to the removal of density of states at the edge bands.

We close this section by mentioning that we performed finite temperature calculations at room temperature; this was achieved by implementing a temperature dependence of the Fermi factors through the Fermi-Dirac distribution. We found that the only significant change is in that the onset energy of the coefficients $\xi^{ab}$, $\xi^{abcd}_{C(R)}$, and $\eta^{abcd}$ are smaller. However, the magnitudes of the coefficients at energies near the lower onsets are several orders of magnitude smaller that the magnitudes of the corresponding coefficients at zero temperature near their energy onsets.

V Limits of the model

The model employed in this work inherits the limits of applicability of time-dependent perturbation theory, which is restricted to situations of low electron-hole pair densities haug04 (for high injection densities a density matrix formalism could be employed to study the dynamics). The regime of validity of the perturbation treatment used here can be estimated: we require the populated fraction of excited states accessible to a typical Gaussian pulse to be small.

V.1 Graphene sheet

As a reference, we first consider monolayer graphene. When the electric fields of the optical beams are all aligned along $\hat{\bm{x}}$, the one- and two-photon injection coefficients for a 2D graphene sheet are Rioux2011 given by Eqs. (40) and  (46). For each of $\xi^{xx}_{2D}$ and $\xi^{xxxx}_{2D}$, we set the number of carriers injected per unit area to be less than the number of states per unit area accessible to the optical beam. Then taking the beam intensity as $I(\omega)=2\epsilon_{0}c{\left\lvert E(\omega)\right\rvert}^{2}$, we arrive to

[TABLE]

where $\alpha$ is the time-bandwidth product for the optical beam (which we take as 0.44, typical for a Gaussian beam), $\Delta t$ is the pulse-duration, and $v_{F}\approx 10^{6}$ m/s is graphene’s Fermi velocity.

V.2 Zigzag nanoribbons

The estimate for the nanoribbon case is similar to the graphene sheet, aside from the fact that the areal ratios become length ratios, i.e. for each one of OPA and TPA coefficients we set the number of carriers injected per unit length to be less than the number of states per unit length accessible to the optical beam, giving us

[TABLE]

where $\alpha$ and $\Delta t$ where defined previously, $v_{e}$ is the velocity of the injected electrons in the conduction band, given by the matrix element $v_{nn}$, and $v_{h}$ is the velocity of the holes injected in the valence band, given by $v_{mm}$. Equation (58) provides the expression for the conventional ($\mathrm{C}$) and ERS processes ($\mathrm{R}$).

In order to compare the limiting intensities of our model for a graphene sheet and for ZGNR, we assume a typical pulse duration of $220$ fs and beam wavelengths of $3.2\mu m$ and $1.6\mu m$ for the $\omega$ and $2\omega$ beamssun-norris2010 . Then we identify the states that contribute at these two wavelengths, and find that, on average, $|v_{e}|+|v_{h}|\approx v_{F}$. From Eqs. (55) and (57), at $\lambda=1.6~{}\mu$m,

[TABLE]

and from Eqs. (56) and (58), at $\lambda=3.2~{}\mu$m,

[TABLE]

Equations (59) and (60) indicate that the limiting intensities of our model are similar for a graphene sheet and for a ZGNR, within an order of magnitude.

We find that, under the assumptions made in this section, the estimated limit for the beam intensities at $\omega$ in the ZGNR and the 2D graphene are about two orders of magnitude below the intensities used in some experiments sun-norris2010 on 2D graphene, where coherent current injection was observed. Due to relaxation processes, of course, the number of allowed carrier excitations below saturation is expected to be higher than our estimates, leading to larger values of the beam intensities for which a perturbation approach would be valid. Based on the estimates in Eqs. (59) and (60), if relaxation processes affect the ribbon samples as effectively as they do for 2D samples, we can expect coherent control in ZGNR to be observable at the higher intensities used in 2D graphene experiments.

VI Summary and Discussion

We have calculated the response coefficients for one- and two-photon charge injection and the two-color current injection in a graphene zigzag nanoribbon; we use the semi-empirical $\bm{k}\cdot\bm{p}$ method to describe the electron wavefunctions by smooth envelope functions.

The only nonzero one-photon injection coefficients correspond to the case of all-$x$ or all-$y$ aligned fields, i.e., $\xi^{xx}$ and $\xi^{yy}$. These two coefficients possess a rich structure of divergences, caused by divergences of the joint-density-of-states originating from the infinite set of parabolic bands present in the zigzag nanoribbon. These two coefficients have distinct selection rules for the allowed transitions.

The two-photon carrier injection coefficients drop off as the fifth power of the photon energy at large photon energies, as they do for monolayer graphene. Moreover, these coefficients possess two classes of divergencies. One corresponds to the joint-density-of-states divergences associated with the parabolic bands. The second class corresponds to divergences arising from resonant conditions, when the two-photon absorption processes arise from sequential one-photon absorption processes between real states. In our calculation here we broadened these resonances phenomenologically, but a more sophisticated treatment of these resonantly enhanced transitions is an outstanding problem on which we hope this work will encourage further study. The onset of the signals is determined by the minimum energy band gap and the selection rules for these coefficients.

We calculated the electron and hole contributions to the conventional and the stimulated electronic Raman scattering (ERS) current injection processes, finding that the only nonzero components are associated with current injected along the length of the nanoribbon, as expected. The behavior of these coefficients as a function of the photon energy follows the behavior of 2D graphene [$\sim(\hbar\omega)^{-3}$] at large photon energies, aside of the resonances present in the ribbons. We have also calculated the so-called swarm velocity of the injected electrons, which inherits a rich structure as a function of the photon energy due to the details of the structure of the injection coefficients. All these calculations were presented for a system at zero Fermi level and zero temperature. However, we also carried finite temperature calculations and found that, within this model, finite temperatures only account for changes at the onset of the signals, which are several orders of magnitude smaller than the nominal values at zero temperature.

Lower bound estimates on the permissible incident intensities for which the calculations here can be valid were presented. They are similar to those of monolayer graphene, where coherent current injection has been observed at much higher intensities than these simple estimates, which do not take into account the relaxation effects in the excited populations. Thus experiments to demonstrate coherent current injection in ZGNR seem to us to be in order.

For experiments contemplated for ribbons of different width than those studied here, it is important to note that simple scaling arguments show that the wider the ribbon, the stronger the confinement of the energy bands. As shown in this work, at low photon energies, the band gap follows a linear relation with respect to the inverse of the ribbon width. Consequently, increasing the width of the ribbon decreases the energy band gap between any pair of bands. This in turn shifts the onset energy of the response coefficients towards zero energy and increases the number of JDOS divergences per photon energy. For instance, the onset of the response coefficients when light is polarized along the length of the ribbon is determined by the bangap between bands $(1,-2)$ (see Fig. 3). For such pair of bands, a linear fit shows that the band gap depends on the ribbon width $W$ as $E^{\mathrm{gap}}_{1,-2}\approx aW^{-1}$ with $a=2.98~{}e\mathrm{V}\cdot\mathrm{nm}$. Besides altering the onset energy of the responses, a larger width also leads to a larger magnitude of the injection coefficients, larger than would be expected simply on the basis of the increase in material; e.g., a width increase of about $15\%$ doubles the size of $\eta^{xxxx}$.

As the outstanding signature of the zigzag nanoribbons are the strongly localized edge states, we have identified their contribution to the carrier- and current-injection processes. In all cases the edge states always participate in the onset of the signals. This lead us to consider a second scenario to study these localized states: given that the dispersion relations of these states are flattened towards zero energy for certain regions in $k$-space, we re-visited our calculations considering doped scenarios. We found that that even small doping levels allow for significant changes around the onset energy of the signals. This is because the large joint-density of states present between the edge states is diminished with nonzero chemical potentials. Due to the relative ease of doping graphene systems, the present work shows that zigzag nanoribbons offer an excellent opportunity to investigate scenarios in which electrical currents can be generated and controlled optically. While more sophisticated treatments of the electron states and the inclusion of electron-electron interaction Yang_Cohen_Louie_2008 ; lado2015 will undoubtedly add to the richness of the injection processes, we hope that the description given here will motivate all-optical current injection experiments. Although coherent control has been studied and observed on graphene sheets, zigzag graphene nanoribbons have the advantage of having optical responses that depend strongly on the geometry and width of the ribbon. Moreover, as shown in the literature, the localized states present in these ribbons are highly sensitive to external fields, doping and functionalization. All these characteristics endow graphene zigzag ribbons with a richness absent in simpler graphene sheets.

VII Acknowledgments

C.S. acknowledges partial support from CONACYT (Mexico) and Rodrigo A. Muniz for useful discussions. J.L.C. acknowledges the support from EU-FET grant GRAPHENICS (618086), the ERC-FP7/2007-2013 grant 336940, and the FWO-Vlaanderen project G.A002.13N. J.E.S. and C.S. acknowledge support from the Natural Sciences and Engineering Research Council of Canada (NSERC). All the authors acknowledge A. Paramekanti and M. Killi for drawing our attention to this problem.

Appendix A Velocity matrix elements