Plasmon modes of graphene nanoribbons with periodic planar arrangements

Among their amazing properties, graphene and related low-dimensional materials show quantized charge-density fluctuations--known as plasmons--when exposed to photons or electrons of suitable energies. Graphene nanoribbons offer an enhanced tunability of these resonant modes, due to their geometrically controllable band gaps. The formidable effort made over recent years in developing graphene-based technologies is however weakened by a lack of predictive modeling approaches that draw upon available {\it ab initio} methods. An example of such a framework is presented here, focusing on narrow-width graphene nanoribbons organized in periodic planar arrays. Time-dependent density-functional calculations reveal unprecedented plasmon modes of different nature at visible to infrared energies. Specifically, semimetallic~(zigzag) nanoribbons display an intraband plasmon following the energy-momentum dispersion of a two-dimensional electron gas. Semiconducting~(armchair) nanoribbons are instead characterized by two distinct intraband and interband plasmons, whose fascinating interplay is extremely responsive to either injection of charge carriers or increase in electronic temperature. These oscillations share some common trends with recent nanoinfrared imaging of confined edge and surface plasmon modes detected in graphene nanoribbons of $100$-$500$~nm width.

Among their amazing properties, graphene and related low-dimensional materials show quantized charge-density fluctuations-known as plasmons-when exposed to photons or electrons of suitable energies. Graphene nanoribbons offer an enhanced tunability of these resonant modes, due to their geometrically controllable band gaps. The formidable effort made over recent years in developing graphene-based technologies is however weakened by a lack of predictive modeling approaches that draw upon available ab initio methods. An example of such a framework is presented here, focusing on narrow-width graphene nanoribbons organized in periodic planar arrays. Time-dependent density-functional calculations reveal unprecedented plasmon modes of different nature at visible to infrared energies. Specifically, semimetallic (zigzag) nanoribbons display an intraband plasmon following the energy-momentum dispersion of a two-dimensional electron gas. Semiconducting (armchair) nanoribbons are instead characterized by two distinct intraband and interband plasmons, whose fascinating interplay is extremely responsive to either injection of charge carriers or increase in electronic temperature. These oscillations share some common trends with recent nanoinfrared imaging of confined edge and surface plasmon modes detected in graphene nanoribbons of 100-500 nm width.
A comprehensive characterization of the dielectric properties of such systems is, however, lacking.
Here, we provide an ab initio study of plasmon excitations in regular planar arrays of GNRs, sorting a wide range of frequencies, from the lower THz to extreme ultraviolet (UV). We use time-dependent (TD) density functional theory (DFT) in random-phase approximation (RPA), emphasizing the 4ZGNR and 5AGNR geometries, which are respectively characterized by four zigzag chains [ Fig. 1(a)] and five dimer lines [ Fig. 1(c)] across the GNR width [26].
The dangling bonds of each GNR array are passivated by hydrogen atoms on both sides, with the C-C and C-H bond lengths being fixed to their nominal values, which differ by less than 1% from the corresponding geometrically optimized values (Figs. S1 and S2 in [34]). The GNR arrays are separated by an in-plane vacuum width of 15Å [Figs. 1(a) and 1(c)], while periodic boundary conditions are used for the direction parallel to the GNR axis, which mimics a situation of long suspended ribbons with fixed edges on the far ends. The equilibrium electronic structures of the systems (Sec. I in [34]) are computed using the local-density approximation (LDA [35]) with norm-conserving pseudopoten-arXiv:1602.00185v4 [cond-mat.str-el] 9 Sep 2016 tials [36] and the plane-wave (PW) basis [37,38]. The three-dimensional periodicity required by PW-DFT is generated by replicating the GNR arrays over an out-ofplane distance L of 15Å, which ensures negligible overlap (but not negligible interaction) of charge density between the replicas.
Accordingly, the different geometry of the assemblies [Figs. 1(a) and 1(c)] produces electronically distinct band dispersions and densities of states (DOS). 4ZGNR [ Fig. 1(b)] appears as a semimetal with the valence and conduction bands overlapping close to the X point. The quasiflat dispersions near the intrinsic Fermi-level E F give rise to strong peaks in the DOS, as opposite to MG where the linear dispersing valence and conduction levels yield a vanishing DOS at E F . 5AGNR [ Fig. 1(d)] is a semiconductor with the valence and conduction electrons having parabolic-like dispersions around a small gap of ∼ 0.36 eV at the Γ point that result in two peaks in the DOS. It should be noted that local spin density calculations suggest the opening of a band gap larger than 0.1 eV in ZGNRs [26][27][28]. Additionally, GW approaches predict larger band gaps in both ZGNRs and AGNRs by roughly 1 eV with respect to local density calculations [27]. Nonetheless, band-gap values of the same order of the LDA band gap of 5AGNR have been measured for some GNRs as wide as about 20 nm grown on Au(111) [30]. Thus, the application of an RPA scheme to the LDA band structure of 5AGNR can be of help in interpreting plasmon measurements on currently synthesized GNR structures [25]. Complementarily, the LDA analysis of a virtually gapless GNR, i.e., 4ZGNR, is particularly instructive to emphasizing the different role played by doping.  The starting point of our TDDFT approach is the density-density response function of noninteracting elec-trons in the GNRs, as given by the Adler-Wiser formula [39,40] which is a corollary of the Kubo formula for a periodic system with associated reciprocal lattice and Brilloiun Zone (BZ) [41]. Hartree atomic units are used throughout this work, unless otherwise stated. In Eq. (1) the electron energies ε νk and states |νk are indexed by the band number ν and the wave vector k in the first BZ. These are taken to be the KSeigensystems of our PW-DFT approach, leading to the electronic structure of Figs. 1(b) and 1(d). The KS wave functions, normalized to unity in the crystal volume Ω, are expressed as linear combinations of PWs that depend on the reciprocal-lattice vectors G associated to the replicated GNR lattices (Sec. I and II in [34]). The correlation matrix elements are given by ρ kq νν (G) = νk|e −i(q+G)·r |ν k + q . The population of single-particle levels is established by the Fermi-Dirac distribution f νk , which we evaluate by sampling temperatures from 300 to 900 K. The factor of 2 accounts for the spin degeneracy, while η is a small (positive) lifetime broadening parameter [42].
Polarization effects are activated by a test electron or photon with incident energy ω and in-plane momentum q that weakly perturbs our systems. These are described by the density-density response function of interacting electrons, which can be obtained in the TDDFT framework as χ GG = χ 0 GG + (χ 0 vχ) GG [43,44]. In the RPA, one neglects short-range exchangecorrelation effects by simply replacing the unknown v by the bare Coulomb terms v 0 GG = 4πδ GG /|q + G| 2 . A serious drawback stems from the long-range character of the Coulomb potential, which allows nonnegligible interactions between repeated planar arrays even at large distances. To cut off this unwanted phenomenon, we replace v 0 GG by the truncated Fourier integral [18,19,[45][46][47] v GG = wherev 0 gg is the Fourier transform of v 0 GG along the out-of-plane axis, while g and G denote the in-plane and out-of-plane components of G.
Within linear response theory, the inelastic cross section corresponding to a process where the external perturbation creates an excitation of energy ω and wave vector q+G is related to the diagonal elements of the inverse dielectric matrix: ( −1 ) GG = δ GG +(vχ) GG . Collective excitations (plasmons) are dictated by the zeros in the real part of the macroscopic dielectric function (permittivity): M = 1/( −1 ) 00 . The so-called energy-loss (EL) function is proportional to the imaginary part of the inverse permittivity: E loss = −Im[( −1 ) 00 ]. Nonlocal field effects are included in E loss through the off-diagonal elements of χ GG [48].
The peculiar electronic structure of 4ZGNR [ Fig. 1(b)] and 5AGNR [ Fig. 1(d)], as compared to the well-known band dispersion of MG, is reflected in the EL spectra of the intrinsic systems shown in Fig. 2. Undoped 4ZGNR and 5AGNR have two high-energy excitations (for ω > 3 eV) that follow one-electron transitions connecting the k points with high DOS in the π-π * , σ-π * and π-σ * bands. These are counterparts to the π and σ-π interband plasmons observed in intrinsic MG [Figs. 2(a) and 2(b)], few-layer graphene and graphite [49][50][51][52]. Specifically, the π and π-σ structures of the GNRs exhibit a discontinuous dispersion vs q and ω, as they are split into more branches [Figs. 2(c) and 2(d), plus Figs. S3(c) and S3(d) in [34]]. This is due to the finite width of the GNRs in the periodic array, generating several one-dimensional bands of π and σ character.
In tr a b a n d In te rb an d The number and dispersion of these bands is also strongly influenced by the GNR width, chirality and inplane distance: the wider the ribbon the more regular the high-energy interband peaks, which approach the π and σ-π peaks of MG as the ribbon width tends to infinity [Figs. 2(a) and 2(b), plus Figs. S3(a) and S3(b) in [34]]. Then, the main designing "ingredients" of the GNR arrays may be finely tuned to reach a specific energy for the π and σ-π excitations, which in turn may be used to change the response of a GNR based device working in the visible (VIS) to UV regime. The low-energy ends of the spectra (for ω < 3 eV) exhibit an extra peak in both metallic and semiconducting GNRs, which is strictly absent in MG at the absolute zero. These structures are more detailed in Fig. 3(a), 3(e), 4(a) and 4(e) below. The large DOS value close to E F in 4ZGNR [ Fig. 1(b)] yields a concentration of n * = 3.96×10 12 cm −2 conduction electrons, which allows the appearance of an intraband plasmon where the charge carriers located on each ribbon of the array oscillate as a single 2D gas [Figs. 2(c), 3(a) and Figs. S3(c), S4(a) in [34]]. This observation is confirmed by the typ-ical square-root-like dispersion of two-dimensional plasmons [53] for low q [ Fig. S4(d) in [34]] that has been mostly observed in extrinsic MG [17], which even in the intrinsic case allows for a weak intraband mode at room temperature assisted by the conduction-electron concentration n * = 1.15×10 11 cm −2 .
On the other hand, the energy gap at E F in 5AGNR [ Fig. 1(d)] yields a negligibly small intraband mode due to the tiny concentration of conduction electrons at room temperature (n * = 8.70×10 8 cm −2 ). The latter is not detectable in Figs. 2(d) and 4(a), but can be, in principle, characterized (Sec. III in [34]). In contrast, an interband two-dimensional plasmon is clearly recorded in the low-energy spectrum of 5AGNR, as testified by the intense signal in Figs. 2(d), 4(a) and 4(d); this corresponds to a collective mode that is triggered by transitions between the valence and conduction DOS peaks at Γ [ Fig. 1(d)]. The collective nature of the newly detected modes in 4ZGNR and 5AGNR is proved in Figs. 3(e) and 4(e), respectively, where we see that each excitation peak in the EL spectrum corresponds to a zero in the real-permittivity Re( m ), at a frequency where the imaginary-permittivity Im( m ) is small. Moreover, as is the case for the high-energy π excitations, these low-energy modes arise from transitions involving the πbands, which means that their intensities and energymomentum dispersions can be modulated according to external factors that change the band levels, such as the already mentioned ribbon width, in-plane distance and chirality. ΔEF=0.00 ΔEF=0.25eV n*=3.96×10 12 cm -2 n*=1.58×10 13 cm -2 n*=1.98×10 13 cm -2 n*=2.73×10 13 cm -2 3. EL spectrum and complex permittivity of intrinsic and extrinsic 4ZGNRs at room temperature. (a-d) Eloss vs ω < 1.5 eV and q < 0.11Å −1 . (e-h) Re( m ), Im( m ) and Eloss vs ω < 1.5 eV at q = 0.039Å −1 . In (a-d) the same color code or intensity scale as in Fig. 2 is used, with the green dots denoting the (ω, q) dispersion of the intraband plasmon.
Let us now see how the dielectric properties of the GNR arrays behave with injecting or ejecting electrons by doping or gating. Extrinsic systems are simulated here by slightly changing the level populations in Eq. (1), in such a way that band dispersions and single-particle KS orbitals are negligibly altered by the applied variations of the f νk factors. For doping levels ∆E F not larger than ∼ 1 eV the high-energy end (ω > 3 eV) of our EL spectra is practically unaffected. On the contrary, unprecedented new features are recorded at the low-energy end (ω < 3 eV).
In Fig. 3(a)-3(d) (and Fig. S4(a)-S4(c) in [34]) we show the low-ω and low-q region of the EL spectrum of 4ZGNR arrays, zooming on the undoped case and analyzing three positive doping levels, whose charge-carrier concentrations are consistent with the measurements of Ref. [54]. We observe a single dispersive structure, the intraband plasmon, which is a genuine collective mode, with the EL peak corresponding to a zero in Re( m ) and a small value of Im( m ) [ Fig. 3(e-h)]. We also notice minor differences in the four EL spectra, with the plasmon energy slightly increasing with increasing ∆E F up to 1.0 eV [ Fig. S4(d) in [34]].
For these low wave vectors (q < 0.02Å −1 ), the intraband mode is the most intense contribution, while the interband plasmon is depressed because the doping partially fills the conduction band near Γ thus inhibiting quasivertical (q → 0, ω) interband transitions. In the 0.02 < q < 0.06Å −1 region, both the intraband and interband plasmons coexist. At larger values of q, the interband plasmon becomes the most intense peak while the intraband plasmon is strongly damped.
A slightly larger value of the doping (∆E F = 0.4 eV) leads to an even more intriguing situation: the single dispersive peak visible in Figs. 4(d) and 4(h) has a double nature, as testified by the kink in peak dispersion and the abrupt decrease in intensity (increase in width) found at q ∼ 0.05Å −1 [ Fig. S5(f) in [34]]. Indeed, interband transitions between the high-DOS points of Fig. 1(d) for q < 0.04Å −1 are strongly quenched by electron population of conduction levels; thus, the intense peak showing the √ q dispersion is mostly originated by the intraband plasmon. Conversely, for q > 0.04Å −1 the intraband plasmon enters a region where it is damped by interband transitions; as a result, most of the spectral weight is concentrated on the interband plasmon, while the overdamped intraband plasmon only appears as a faint peak. Another remarkable effect is the high sensitivity of the I n t r a b a n d I n t e r b a n d arb un arb un arb un arb un intraband plasmon to the type of doping; opposite doping levels, such as the ∆E F = ±0.2 and ∆E F = ±0.3 values of Fig. S5(h) in [34], produce significantly different charge-carrier concentrations and plasmon dispersion curves for energies larger than 0.1 eV and transferred momenta above 0.02Å −1 . Such a sensitivity is ascribed to the slight asymmetry of the valence and conduction electron levels of 5AGNR close to the band gap [inset in Fig. 1(d)].
On the other hand, the interband plasmon is much less influenced by the doping type, because the valence and conduction DOS peaks have similar intensities [ Fig. S5(h) in [34]]. The interplay of the two modes is also modulated by changes in the incident momentum direction, relative to the ribbon axis, due to the tensor character of the GNR dielectric response (Fig. S6 in [34]). In the lower THz region, interband transitions are forbidden by the band gap. The intraband plasmon follows an ω-vs-q dispersion that is consistent with the semiphenomenological relation of Ref. [31], and appears at the same scale as the low-q dispersions of Ref. [33], computed from the TB approximation where 5AGNR is virtually gapless [ Fig. S7(a) in [34]]. Then, the TDDFT response of narrow-width GNRs has the correct q→0, ω→0 limiting behavior predicted by (non-ab initio) approaches on larger GNR structures, currently available for experiments.
As a final remark we observe that in semiconducting GNRs, due to their small band gap, the temperature plays a role in dictating the populations of the levels close to E F (Sec. IV in [34]). Accordingly, the intraband plasmon mode can be triggered by working at temperatures larger than ∼ 500 K, which may have a crucial role in relation with power consumption of nanodevices. Chargecarrier concentrations generated by temperature increase are nevertheless much smaller than those obtained with doping or gating. For this reason, no particular interference is recorded between intraband and interband plasmon modes (Fig. S8in [34]).
In summary, we have discussed the dielectric properties and plasmon dispersions in planar GNR arrays scrutinizing the excitation energy regime going from the THz to the UV scale, by an ab initio strategy based on TDDFT+RPA.
On the THz regime, we have detected new collective modes of different nature. Semimetallic GNRs display an intraband 2D plasmon with large intensity relative to the high-energy plasmons even in the intrinsic case. Semiconducting GNRs experience a fascinating interplay of intraband and interband collective modes, whose relative intensities and dispersions are strongly influenced by the actual occupation of single-particle levels near the Fermi energy.
Some recent calculations have reported the existence of two extrinsic plasmons in MG and bilayer graphene (BLG) [17][18][19][20]. In particular, the plasmon coupling in BLG [20] shares some common features with 5AGNR at similar doping conditions; i.e., one of the two modes disperses like q 1/2 and the other is quasivertical.
Indeed, a first experimental evidence of an edge (interband) plasmon superimposed to a conventional (intraband) plasmon has been given in patterned GNRs grown on Al 2 O 3 [25]. The two modes are well resolved in space on GNR samples of 480 nm width, at a working frequency of ∼ 0.15 eV and a doping level of ∼ 0.3 eV. In our narrow-width GNR, the interband and intraband features are resolved in momentum space, only.
Our calculations demonstrate that it is possible to construct new materials with plasmonic resonances that are tunable to suit a specific demand in both the VIS-UV and THz regimes, by means of chemical doping, electronic gating, and also a careful choice of the geometry. These findings if confirmed by further experiments will widen the perspectives on applications of GNR arrays for the engineering of nanophotonic and nanoelectronic devices.
C.V.G. acknowledges the financial support of "Secretaria de Educación Superior, Ciencia, Tecnología e Innovación" (SENESCYT-ECUADOR). All authors thank Dr V.M. Silkin from the University of the Basque Country for his invaluable support in developing the TDDFT code. * corresponding author: antonello.sindona@fis.unical.it i

Supplemental Material
The main text has been concerned with the dielectric response of two narrow-size graphene nanoribbons (GNRs), organized in periodic planar arrays with the zigzag (Z) and armchair (A) geometries. Time-dependent densityfunctional theory has been carried out for both intrinsic (undoped) and extrinsic (doped, gated) systems within the random phase approximation (RPA). The energy loss spectra of the GNR arrays have been reported and discussed for a range of incident energies below 20 eV and transferred momenta smaller than 0.8Å −1 . These supplementary notes are organized as follows: Sec. I provides the details on ground state computations and geometry optimizations for the GNR arrays, denoted in the main text as 4ZGNR and 5AGNR [26]; Sec. II deals with the calculations of the plasmon structure of the intrinsic systems in comparison with that of pristine graphene; Sec. III gives some additional arguments on the effect of doping in 4ZGNR and 5AGNR, focusing on the energy-momentum dispersion of the intraband and interband plasmons; the semiclassical limit of the approach is validated by zooming on the lower THz region; Sec. IV is devoted specifically to the role played by increasing the electronic temperature on 5AGNR.

LOCAL DENSITY CALCULATIONS AND GEOMETRY OPTIMIZATION
Density functional calculations were performed using the plane-wave (PW) basis set [37], represented by the space functions PW k+G (r) = Ω −1/2 0 e i(k+G)·r , where k is a wave vector in the first Brillouin Zone (BZ), G a reciprocal vector, and Ω 0 denotes the unit-cell volume of the real-space lattice. The number of PWs was limited by the cut off condition: |k + G| 2 /2 < 25 Hartree.
The local density approximation (LDA) was used, as parametrized by the Perdew-Zunger form of the uniformgas correlation energy [35]. Norm-conserving pseudopotentials of the Troullier-Martins type [36] were adopted.  Fig. 1(a), which are plotted in Fig. 1(b) of the main text.
The basic parameters of the unit cells of 4ZGNR and 5AGNR, i.e., the C-C and C-H bond lengths and angles, were determined by geometry optimization using the Broyden-Fletcher-Goldfard-Shanno method [37].The resulting average C-C and C-H distances turned out to be slightly different from the nominal values of 1.42 and 1.09Å, respectively, used in the main text to produce the non optimized structures of Figs. 1(a) and 1(c).
Periodic planar arrays of 4ZGNR and 5AGNR were then generated with both the optimized and nonoptimized unit cells assuming an in-plane vacuum distance of 15Å. The three-dimensional periodicity required by PW-DFT was set by replicating the GNR arrays with an out of plane lattice constant of 15Å, corresponding to the unit-cell volumes Ω 0 ∼ 950Å 3 for 4ZGNR and Ω 0 ∼ 1470Å 3 for 5AGNR.  Fig. 1(c), which are plotted in Fig. 1(d) of the main text.
Both geometry optimization and ground state calculations were carried out using an unshifted (Γ-centered) Monkhorst-Pack (MP) grid [38] made of 60 × 10 × 1 wave vectors, which results in a uniform sampling of the irreducible area of the first BZ (see Fig. S6).
The band energies obtained with the optimized ii cells [denoted OPT in Figs. S1(b) and S2(b)] were found to differ negligibly from those of the non optimized structures [denoted NOM in Figs. S1(b) and S2(b)] that are also reported in Figs. 1(b) and 1(d) of the main text. Regular planar arrays made of wider ribbons were tested as well (though not reported here) for more systematic future studies.
Electronic structure calculations were also run on monolayer graphene (MG) for comparison purposes (see Figs. 2 and S3). In this case, the honeycomb lattice structure of C atoms (separated by a nearest neighbor distance of 1.42Å) was replicated with an out-of-plane vacuum distance of 15Å, i.e., a unit-cell volume Ω 0 ∼ 110Å 3 . The first BZ of the two-dimensional material was represented on an MP grid of 60×60×1 k points. In the ΓX domain, we calculated the KS energies ε νk and wave functions being normalized to unity in the crystal volume Ω = N Ω 0 , with N expressing the number of sampled k points (in the MP mesh) and the G-sum being limited by the cut off condition set forth above. The full plasmon structure of the intrinsic GNRs was explored at frequencies ω ≤ 20 eV [Figs. 2(c), 2(d), S3(c) and S3(d)]. Momentum values q parallel to the GNR axis were considered in the range of ∼0.02 to ∼0.8Å −1 .
The converged ground state density was determined with ∼ 15000 PW coefficients {c νk+G } per occupied wave function, leading to the electronic structure of Figs. 1(b) and 1(d). n(r) was subsequently used in a non self-consistent run to determine the KS energies and wave functions on an MP mesh of 180 × 1 × 1 k points. 120 bands (indexed by ν) were included to have a reliable dielectric response in the considered frequency and momentum ranges.
EL calculations for MG were run on an MP grid of 180×180×1 wave vectors in the irreducible first BZ of the two-dimensional material, including 80 bands and ∼ 5000 PW coefficients per wave function [Figs. 2(a), 2(b), S3(a) and S3(b)].
The non interacting density-density response (or unperturbed susceptibility) χ 0 GG was computed for MG, 4ZGNR, and 5AGNR at room temperature. As reported in Eq. (1) of the main text, χ 0 GG includes all possible one-electron processes between occupied and empty band levels, separated in energy by ε νk − ε ν k+q . The dynamical screening features are provided by the retarded Green's functions (ω+ε νk −ε ν k+q +iη) −1 , which contains the positive infinitesimal η. Transition rates between the band levels are made of products of correlation matrix elements weighted by the Fermi-Dirac factors f νk − f νk+q . The unperturbed susceptibility of 4ZGNR and 5AGNR was represented with 121 reciprocal lattice vectors, sorted in length and cut off by the condition |G| = |g| 2 + G 2 < 7.2Å −1 , with g and G denoting the inplane and out-of-plane components of G, respectively. 60 q-values were sampled along ΓX, and 3001 frequencies below 20 eV. η was replaced with a (finite) lifetime broadening parameter of 0.02 eV. The unperturbed susceptibility of MG was represented with 50 reciprocal lattice vectors, being such that |G| < 3.2Å −1 [17]. The ΓK and ΓM high-symmetry paths were sampled with 60 and 90 q-values, respectively. 3001 ω-values below 30 eV were selected with η = 0.02 eV.
Accordingly, unperturbed susceptibility matrices were calculated for all sampled frequencies and momenta. With the χ 0 GG matrices at hand, the interacting densitydensity response function (susceptibility) χ GG was obtained from the central equation of TDDFT (see main text). The latter is a Dyson-like equation, which can be rewritten as by making explicit the summation indices. The effective potential matrix elements in Eq. (S4) may be calculated from the truncated Fourier integrals given in Eq. (??), which yields The v GG terms have been proved to efficiently cut off the spurious interaction between replicas of planar graphenebased systems [18,[45][46][47]. Finally, the EL spectra were computed by where is the solution of Eq. (S4). Fig. S3 provides a complementary view to Fig. 2, with the EL spectra of MG, 4ZGNR and 5AGNR being reported for selected q values in the corresponding momentum domains. As pointed out in the main text, the one-peak structure of the π and σ-π plasmons in MG is turned into well resolved sub-peak structures in the GNRs. These originate from the several π and σ bands of the 1D materials [ Figs. 1(b), 1(d), S1(a) and S2(b)].
The low-energy intraband and interband modes for 4ZGNR and 5AGNR, respectively, appear to be well resolved and have higher intensities with respect to the π and σ-π modes. In particular the intraband mode of 4ZGNR is generated by a conduction-electron concentration n * of 3.96×10 12 cm −2 [ Fig. 2(c)].
Such a mode is absent in MG at the absolute zero because of the vanishing behavior of the density of states (DOS) at the inequivalent K points of the first BZ. The change in the Fermi-Dirac statistics at room temperature leads to n * = 1.15×10 11 cm −2 [Figs. 2(a) and 2(b)], which in turns is responsible for a weak structure in the optical region of MG. This is barely visible at the lowest q's in the plots of Figs. 2(a), 2(b), S3(a) and S3(b), due to both the on-set of the π structure and the energy resolution (0.01 eV) used to sample the 0 − 20 eV range.

DOPED SYSTEMS
To explore the dielectric properties of the GNRs at infrared frequencies, we used an MP mesh of 1000 × 1 × 1 wave vectors with ∼ 30 bands, covering the frequency region ω ≤ 2 eV. Besides the dielectric properties of the pristine GNR arrays, we also studied how these properties change with doping. Indeed doping can happen as an unwanted effect due to chemical contamination or it    can be introduced in a controlled way, i.e., after functionalization of the GNR or through a gating potential.
Doping or gating was simulated by changing the occupation factors f νk and f νk+q in χ 0 GG (see Eq. (1)), according to the concentration of injected or ejected electrons. It can be safely assumed that the doping values sampled in this Letter have a negligible effect on the KS electronic structure [17,18].
All other settings were the same as for the intrinsic systems. Figs In semimetals like mono-and bilayer graphene, even a small doping yields the formation of a quasi twodimensional gas of charge carriers (electrons or holes). The fact that the conduction or valence band becomes partially occupied or unoccupied causes the appearance of low-energy intraband single-particle transitions and collective excitations.
In 4ZGNR, the non vanishing DOS at the Fermi level [ Fig. 1(b)] causes the appearance of a well resolved intraband mode being weakly affected by the doping (gating). As shown in Fig. S4(d), the plasmon resonance en-ergies slightly increase with increasing the injected electron concentration, following the typical squareroot-like dispersion of the two-dimensional plasmon in MG as q → 0. These features are separately visible in Figs.

3(a)-3(d).
It is worth mentioning that the effective charge-carrier density contributing to the intraband mode in 4ZGNR is of the same order of magnitude for all sampled doping levels, i.e., 10 13 cm −2 . Therefore, moderate energy shifts, up to 0.5 eV, are recorded in the intraband plasmon dispersions [ Fig. S4(d)]. In contrast, doping levels of 0.5 and 1 eV in MG lead to a more pronounced charge-carrier density increase, of 10 13 and 10 14 cm −2 , respectively. Accordingly, a substantial variation in the plasmon energy of the extrinsic two-dimensional plasmon of MG has been reported [17].
In 5AGNR, the presence of a small gap with high DOS [ Fig. 1(d)] allows the appearance of an interband plasmon that seems to be weakly influenced by electron injection or ejection.
On the other hand, positive or negative Fermi energy shifts generate an intraband plasmon that is extremely sensitive to the doping level and may interfere with the interband plasmon. These other features are highlighted v in Fig. S5(h), where the intraband and interband plasmon dispersions of Figs. S5(a)-S5(g) are plotted together vs q and ω. Some of them are also separately visible in Figs. 4(a)-4(d). Notice that a tiny intraband contribution should be present in undoped 5AGNR because of the Fermi-Dirac population of low lying conduction states at room temperature.
Indeed, 5ANGR is characterized by a conductionelectron concentration of ∼10 9 cm −2 at T = 300 K, which is negligibly small with respect to that of 4ZGNR and MG. Nonetheless, local maxima in E loss are recorded for q < 0.04Å −1 , yielding the (ω, q) dispersion in Fig. S5(h). Also detectable in Fig. S5(h) is the effect of the slight asymmetry in the valence and conduction bands of 5AGNR close to band gap [inset in Fig. 1(d)], which appears to influence the plasmon structure especially in the intra-band components. This is due to the different charge-carrier concentrations originated from positive and negative doping values.
Thus, the intraband plasmon of 5AGNR depends of the doping type (sign) even at frequencies of the order of ∼ 50 THz.
We should also mention that in a realistic situation of suspended GNR arrays on top of a substrate, Fermi level shifts of 0.5 − 1 eV, obtained for example by applying a gate voltage in air, can lead to a collapsing of the sample. Nevertheless, the carrier density concentration predicted in our 4ZGNR-with a doping level of 1.0 eV-is similar to the experimental concentration of 2.5×10 13 cm −2 achieved by depositing alternating wafer-scale graphene sheets and thin insulating layers [54]. Furthermore, electrical conductance measurements at small and room temperatures, performed on a set of doped parallel GNRs, with widths ranging from 14 to 63 nm, has given a chargecarrier concentration of 3.6×10 12 cm −2 [55]. This value is in close agreement with the predicted charge carrier concentration of our doped 5AGNR arrays (Figs. 4 and S5).
The interplay between the intraband and interband modes can be also detected by probing 5ANGR with oblique incident momenta of in-plane modulus q forming different angles θ relative to the ribbon axis. This is shown in Fig. S6 where the EL spectra have been computed by re-sampling the first BZ of the system with an MP mesh of 270×40×1 points. We see that specific combinations of modulus and angle in the incident momentum lead the two modes to be superimposed, whereas other choices leave the two-peak structure separated in momentum space.
The role of positive doping is also spotted in Fig. S7(a), where a comparison of the intraband plasmon dispersion of 5AGNR is attempted with the calculations of Ref. [33]. The latter were based on a tight-binding (TB) gapless two-band model at the absolute zero. A two-dimensional Coulomb potential was adopted, depending on a localization parameter of the p z orbitals. The correlation matrix elements (S3) were approximated by the low-q limit form ρ kq νν (0) ≈ νk|ν k+q . Of course, the detailed properties of the plasmon modes are critically dependent on the valence and conduction band dispersions, which are significantly different in DFT-LDA and TB calculations. The TB electronic features of 5AGNR allow only for intraband modes that, nevertheless, appear at the same scale as the TDDFT intraband plasmons of doped 5AGNR.
To have a closer look at the THz region, we performed TDDFT calculations on an MP mesh of 12500 × 1 × 1 points, including only the contributions of valence and conduction states of 5AGNR in Eq. (1). Well converged results were obtained with a broadening lifetime parameter η of ∼ 0.5 THz.
The EL spectra for 5AGNR at T = 300 K and ω < 25 THz are shown in Fig. S7(b). The plasmon dispersion curves are reported in Fig. S7(c) and compared to the semiclassical approach of Ref. [31], where plasmon resonances in doped GNR arrays of 10 − 100 nm width were described by the following analytical expression: Eq. (S8) (given in Hartree au) depends on the concentration n * of conduction electrons, whose effective  Then, our TDDFT calculations on a narrow-width GNR have the correct q→0,ω→0-limit that fits with (non-ab initio) predictions on larger GNR structures currently available for experiments. Such TB approaches, due to their simplified description of the band energies and wave functions, are only able to predict the existence of one plasmon mode, the extrinsic intraband plasmon.

TEMPERATURE EFFECTS IN 5AGNR
As shown in the main text and in the previous section, the large tunability in the intraband plasmon mode of 5AGNR is critically dependent on the electron or hole occupancy of conduction or valence states, lying within an energy window of 0.5 − 1 eV around the Fermi level. Such a population is given by the Fermi-Dirac statistical factors f νk and f νk+q in Eq. (1) that are significantly influenced, within the considered energy range, by even moderate temperature changes below ∼ 1000 K.
This effect inevitably plays a major role in any nanodevice design approach. To quantify and characterize it, we ran EL calculations on undoped 5AGNR with electronic temperature values larger than 500 K. All other settings were the same as for room temperature calculations. The resulting EL spectra are reported in Fig. S8 for energies below ∼ 1 eV and momenta smaller than 0.03Å −1 .
At T = 300 K [ Fig. S8(a)], the interband plasmon is clearly visible. Nevertheless a faint intraband peak may be spotted at energies below ∼ 0.1 eV, and the intraband plasmon dispersion can be computed. As the temperature is increased to 500 K [ Fig. S8(b)], the intraband plasmon peak begins to appear in the same intensity scale as the interband plasmon. At higher temperatures, say, T = 700, 900 K [Figs. S8(c) and S8(d)], the intraband plasmon is well resolved and also well separated from the interband plasmon.
The increase of the intraband plasmon intensity is readily understood considering how the electronic temperature affects the population of the KS states near the Fermi level. At room temperature the number of electrons that are capable to overcome the 0.36 eV gap is small, with a concentration of the order of 10 9 cm −2 . Thus, interband excitations are dominant. At T = 500 K the electron population of the conduction band becomes appreciable, with a concentration of roughly 2×10 10 cm −2 , generating the small peak in the EL spectrum. As the temperature further increases (T ≥ 700 K), the smearing width of the Fermi-Dirac distribution function increases, the conduction electron concentrations become larger than 10 11 cm −2 , and the intraband plasmon fully appears in the EL spectrum. Charge-carrier concentrations triggered by temperature increase are nevertheless much smaller than those obtained with doping or gating. For this reason, no particular interference is recorded in Fig. S8 between intraband and interband plasmon modes. It should be also mentioned that in the realistic case of GNR arrays suspended on top of a substrate, the substrate phonons become relevant and significantly affect the plasmon damping at temperatures above 500 K. vii T=300 K n*= 8.70×10 8 cm -2 T=500 K n*= 2.23×10 10 cm -2

FINAL REMARKS
We have discussed the response of planar GNR arrays with armchair and zigzag shaped edges to a probe particle (i.e, an electron or a photon) of energy below 20 eV and in-plane momentum smaller than 0.8Å −1 . The use of an ab initio strategy based on TDDFT in the RPA has let us characterize both the intrinsic and extrinsic plasmon structures of the systems, which could not have possibly been predicted by less sophisticated TB approaches. At high-energies (> 3 eV), we have found the two standard interband excitations of carbon-based materials, namely the π and σ-π plasmons, whose energymomentum dispersions are strongly influenced by the GNR geometry. At lower energies (< 1.5 eV), we have detected one intraband plasmon mode for the semimetallic ZGNR. In the same energy range, we have identified two modes in the semiconducting AGNR, i.e., one interband and one intraband plasmons. The latter shows a large sensitivity to doping level as well as electronic temperature. The interplay between the two collective excitation in AGNR, which co-exist and interact with each other, represent the main result of our study. A similar interference effect, between a two-dimensional plasmon and acoustic-like plasmon, has been reported in a recent study on bilayer graphene [20]. We have checked that the intraband and interband plasmons appear also in larger AGNRs. These findings require further confirmation by measurements, although the experiments of Ref. [25] seem to suggest that interband and intraband plasmon in low-gap GNRs exist, can be controlled and exploited for future nanodevice technology.