Infra-Red Active Dirac Plasmon Serie in Potassium Doped-Graphene (KC8) Nanoribbons Array on Al2O3 Substrate

A theoretical formulation of the electromagnetic response in graphene ribbons on dielectric substrate is derived in the framework of the ab initio method. The formulation is applied to calculate the electromagnetic energy absorption in an array of potassium-doped graphene nanoribbons (KC8-NR) deposited on a dielectric Al2O3 substrate. It is demonstrated that the replacement of the flat KC8 by an array of KC8-NR transforms the Drude tail in the absorption spectra into a series of infrared-active Dirac plasmon resonances. It is also shown that the series of Dirac plasmon resonances, when unfolded across the extended Brillouin zones, resembles the Dirac plasmon. The Dirac plasmon resonances’ band structure, within the first Brillouin zone, is calculated. Finally, an excellent agreement between the theoretical absorption and recent experimental results for differential transmission through graphene on an SiO2/Si surface is presented. The theoretically predicted micrometer graphene nanoribbons intercalation compound (GNRIC) in a stage-I-like KC8 is confirmed to be synthesized for Dirac plasmon resonances.

One of the biggest challenges in applied plasmonics or photonics is finding a way to excite and manipulate the 2D plasmons directly by the incident electromagnetic field. Even though 2D plasmons modes produce a strong localized electric field, that field is evanescent and therefore cannot be excited directly by light. In single-layer graphene, the Dirac plasmon can be excited only indirectly, e.g., by exciting the localized plasmons on the AFM tip, which then excites the Dirac plasmon in the graphene [14]. However, subwavelength nanostructures such as graphene nanoribbons (GNRs) support 'plasmon resonances' with very localized electric field that can be radiated into the surrounding area; thus, it can also be pumped directly by an external radiation.
The measurements of the electromagnetic field transmission through the GNR arrays on SiO 2 /Si substrate clearly show the existence of strong plasmon resonances in the THz and mid-infrared frequency range (depending on nanoribbon thickness) [15,16]. Moreover, the infrared near-field imaging of GNRs on an Al 2 O 3 substrate shows that, in addition to the conventional plasmonic resonances, GNRs support the edge plasmons (distributed along the graphene edges [17]). The plasmon resonances in submicrometer multilayer graphene ribbons on an Si/SiO 2 substrate caused by different doping concentrations (graphene Fermi energy) enable the tuning of the IR reflectivity [18]. In the above experiments, the plasmon resonances appear in the THz and IR frequency ranges, thus hybridization with IR active Fuchs-Kliewer surface optical (SO) phonons [19] in polar substrates was also considered.
Thus far, the theoretical description of the optically active plasmon resonances in submicrometer graphene ribbons, despite providing useful information, such as plasmon wave functions or analytical dispersion relations, is limited to semi-analytical modeling, mostly based on the simple Drude model conductivity, which does not take into account the substrate polarization [20,21]. The ab initio calculations of the energy loss function in semimetallic (zigzag) and semiconducting (armchair) nanoribbons [22] provide information about the interesting interplay between the intraband and interband plasmons. The ab initio calculations of the dielectric response in different GNRs, taking into account the electron scattering with SO phonons in various polar and nonpolar substrates, provided very useful information about plasmon propagation length and plasmon-phonon hybridization [23]. However, these ab initio studies have only focused on a few nanometers thick GNRs, while the radiative plasmon resonances were not studied. Interesting experimental and theoretical studies show that the THz absorbance of the graphene monolayer can be considerably enhanced by depositing the graphene on a dielectric substrate of specific dielectric permittivity and thickness [24,25]. Another very interesting phenomenon of the modulation of the THz graphene absorption is achieved by applying an optical pump signal, which modifies the conductivity of the graphene sheet [25]. However, the sharp plasmon resonances, which occur only in graphene ribbons, are not observed in these investigations.
A way of approaching the creation of a Dirac plasmon (DP) in GNRs has been explored via the synthesis of alkali-metal-doped graphene on metallic substrates, being extensively studied experimentally and theoretically [26][27][28][29]. These experiments have shown that, by doping graphene with electron donors, the Dirac plasmon resonances can be excited and extensively studied. Furthermore, the use of advanced multilayer graphene nanoribbons will help control the plasmonic resonances derived from the perpendicular electric field in those nanostructures [18].
In this paper, we explore the electromagnetic response in an array of potassium-doped graphene nanoribbons (KC 8 -NR) deposited on a (Al 2 O 3 ) dielectric substrate. The singlelayer KC 8 (KC 8 -SL) optical conductivity tensor σ 0 µν (ω) and the bulk Al 2 O 3 macroscopic dielectric function s (ω) are calculated from first principles. Special attention is paid to the series of Dirac plasmon resonances (DPR) in d = 50, 100 and 200 nm thick KC 8 -NR arrays. We show that the series of DPR consists of a series of dipolar or infrared-active DPR and a series of non-dipolar or dark DPR. We demonstrate that the DPR in the first Brillouin zone (1stBZ) when unfolded in the extended Brillouin zone (exBZ) resembles the Dirac plasmon (DP) in the KC 8 -SL. For smaller separations between the nanoribbons, thus with a stronger interaction between, which causes dispersion, we calculate the resulting DPR band structure within the 1stBZ. Finally, we apply the proposed formulation to calculate the electromagnetic absorption in the doped graphene microribbons on the SiO 2 /Si surface. The results are then compared with experimental measurements of differential transmission through the same sample [15], and KC 8 micrometer length GNRs intercalation compounds were synthesized.
The rest of the paper is organized as follows. In Section 2, we present the theoretical model used to calculate the electromagnetic energy absorption A in the array of KC 8 -NR. In Section 3, we present the ab initio computational details and the results for the KC 8 -SL optical conductivity σ 0 yy (ω) and the bulk Al 2 O 3 macroscopic dielectric function s (ω). In Section 4, we present the results for the absorption spectra A in different arrays of KC 8 -NR and DPR band structure. In Section 5, we show the comparison with available experimental results. Section 6 contains some conclusive remarks. Unless stated otherwise, atomic units are used, i.e., e =h = m = 1, where e is the electronic charge,h is the reduced Plank constant, m is the electron mass and c is the speed of light in vacuum.

Calculation of the Electromagnetic Energy Absorption A
In this section, we briefly describe the method of calculation of the electromagnetic energy absorption in the system consisting of an array of KC 8 -NR deposited on a dielectric Al 2 O 3 substrate (according to the formulation developed by [30]). The dielectric substrate occupies the region z < 0, while the KC 8 -NR of width d and period l are arranged so that their graphene layers are placed at z = z 0 above the dielectric surface, as illustrated in Figure 1. The substrate polarization is described by the dielectric function s , while the polarization of the dielectric media occupying the region z > 0 is described by the dielectric function 0 . If we assume that the sample is driven by an external electromagnetic field of unit amplitude, frequency ω and wave vector k, with incidence perpendicular to the surface, the polarization e is parallel to the surface and ω = kc, where c is the speed of light. The electromagnetic energy absorption in the array of KC 8 -NR can be obtained by using the following expression [30]: Considering that the thickness of the KC 8 -SL is significantly smaller than the IR or visible light wavelength (λ = 2πc ω > 100 nm), and taking into account only the bare polarizations (µ = x or y), Equation (2) can be simplified as Here, the tensor σ µν represents the screened conductivity of the KC 8 -NR array, which is the solution of the Dyson equation [31]: where σ 0 µν is the nonlocal irreducible conductivity tensor of the KC 8 -NR array and Γ αβ is the propagator of the bare electric field corrected by the presence of the Al 2 O 3 substrate. The part of the electromagnetic energy absorbed by the Al 2 O 3 substrate is neglected, which is a reasonable approximation considering that the Al 2 O 3 is mostly transparent in the frequency interval of interest. Both tensors, σ 0 and Γ, are described explicitly below.
If the nanoribbon width is much larger than the unit cell constant (d >> a), the nonlocal effects in the x − y plane are negligible and the optical response of the KC 8 -NR can be approximated by the local optical conductivity. Moreover, since the thickness of the KC 8 -SL is significantly smaller than the wavelength λ > 100 nm, it can be treated as a 2D crystal, localized, e.g., in the graphene, z = z 0 plane. In this approximation, the conductivity tensor becomes: where ρ = (x, y) is the 2D position vector, and σ 0 µν (ω) is the 2D optical conductivity of the KC 8 -SL. All this enables the Fourier expansion of σ 0 µν : Here, the 2D reciprocal vectors are G = (0, g), with g = 2πn l ; n = 0 ± 1, ±2, . . ., and where Q = (Q x , Q y ) is the 2D transfer wave vector. The propagator Γ remains translationally invariant in the x − y direction so it can be Fourier transformed as where Γ µν,gg (Q, ω, z, z ) = Γ µν (Q + G, ω, z, z )δ gg .
Using the expansions (7) and (9) and assuming that the screened conductivity σ can be transformed the same way as the conductivity σ 0 (expansion Equation (7)), the Dysons Equation (4) transforms into matrix equation for the screened conductivity After inserting the Fourier expansion of the screened conductivity (Equation (7), where σ 0 → σ) into Equation (3), and using the identity σ µµ,g−g 0 g −g 0 (Q + G 0 , ω) = σ µµ,gg (Q, ω), we obtain the final expression for the electromagnetic energy absorption rate rate per unit area We use the expression (12) to determine the intensity of the electromagnetic modes beyond the optical limit as well, e.g., in the nonradiative or evanescent region (Q > ω/c), simply by using the conductivity σ calculated for a finite wavevector (Q = 0).

Electric Field Propagator
The propagator of the electric field can be written aŝ where the propagator of the 'free' electric field (or free photon propagator) is [31,32] 2πω The propagator of the scattered electric field in the region z > 0 is [32] Here, the unit vectors of the s(TE) polarized electromagnetic field are e 0,± and Q 0 and z are the unit vectors in the Q and z directions, respectively. More specifically, Γ sc represents the electric field produced by an external point dipole which is reflected at the dielectric surface. Therefore, '−' represents the incident electric field, while '+' represents the reflected electric field. Γ 0 represents the 'direct' electrical field produced by point the dipole so that the superscript '0' represents the spherical forward propagating field. The reflection coefficients of the s(TE) and p(TM) polarized electromagnetic waves at the media/substrate interface are r s = (β 0 − β s )/(β 0 + β s ) and r p = (β 0 s − β s 0 )/(β 0 s + β s 0 ), respectively. The complex wave vectors in the perpendicular (z) direction are β 0,s = ω 2 c 2 0,s (ω) − Q 2 .

Calculation of RPA Optical Conductivity σ µν (ω)
Since we study a 2D crystal which consists of just two atomic layers, its electromagnetic response is strongly dispersive in the perpendicular z direction. For this reason, we define the spatially dependent conductivity where the conductivity matrix is defined as [30] (16) Here, the current vertices are and the current produced by the transitions between the Bloch states φ * nK → φ mK+Q is defined as where G z = 2πn/L; n ∈ Z represents the reciprocal vector in the z direction, K = (K x , K y ) is the 2D wave vector and φ nK and E nK are the Bloch wave functions and energies obtained by the DFT calculations. The spin quantum number 's' is merged with the band quantum number, i.e., n ≡ (n, s), Ω = S × L is the normalization volume, S is the normalization surface, L is the superlattice unit cell in the z direction and f nK = [e (E nK −E F )/kT + 1] −1 is the Fermi-Dirac distribution at temperature T. η intra and η inter represent the phenomenological intraband and interband damping parameters, respectively. The two-dimensional conductivity used in this study is defined as It should be noted that this defined conductivity do not depend on lattice parameter L. The KC 8 -SL is a conductive 2D crystal so it is appropriate to divide its RPA optical conductivity into intraband and interband contributions which are both determined from the optical limit of the nonlocal conductivities . According to (16), the nonlocal intraband (n = m) conductivity is [30] σ intra where the effective number of the charge carriers is Here, K ∈ 1.SBZ indicates that summation is performed within the first surface Brillouin zone. The nonlocal interband (n = m) conductivity is [30] An alternative modeling of the KC 8 -SL conductivity could be done in analogy with the graphene conductivity modeling in the Dirac cone approximation [33], but taking into account the parabolic K(σ) band crossing the Fermi level.

Calculation of Substrate Macroscopic Dielectric Function s (ω)
We assume that the dielectric media is vacuum (i.e., 0 = 1) and that the substrate is aluminium-oxide (Al 2 O 3 ) described by the macroscopic dielectric function s (ω). To calculate s (ω), we start from the 3D Fourier transform of the independent electron response function where k ∈ 1.BZ indicates that summation is provided within first Brillouin zone. The charge vertices are defined as Here, k = (k x , k y , k z ), q = (q x , q y , q z ) and G = (G x , G y , G z ) are the 3D wave vector, the transfer wave vector and the reciprocal lattice vector, respectively, and the integration is performed over the normalization volume Ω. We use the response matrix (22) to determine the dielectric matrix as where the bare Coulomb interaction is v GG (q) = 4π |q+G| 2 δ GG . Finally, the macroscopic dielectric function is determined by inverting the dielectric matrix

Computational Details
The KS wave functions φ nK and energies E nK used to calculate the RPA conductivities σ µν and the substrate macroscopic dielectric function (ω) are determined using the planewave self-consistent field DFT code (PWSCF) within the QUANTUM ESPRESSO (QE) package [34]. For all crystal structures (KC 8 -SL, doped graphene and bulk Al 2 O 3 ), the coreelectrons interaction is approximated by the norm-conserving pseudopotentials [35,36]. The exchange correlation (XC) potentials in the KC 8 -SL and Al 2 O 3 are approximated by the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) functional [37] and in the graphene by the Perdew-Zunger local density approximation (LDA) functional [38]. The ground state electronic density in KC 8 -SL is calculated using the 8 × 8 × 1 Monkhorst-Pack K-mesh [39], the plane-wave cut-off energy is 60Ry and we use the hexagonal Bravais lattice, where a = 4.922 Å and the separation between the KC 8 layers is L = 2.5a. Since the graphene unit cell is doped by holes, with the doping concentration 1.5 × 10 13 cm −2 , the graphene ground state electronic density is calculated using a dense 101 × 101 × 1 K-mesh and the plane-wave cut-off energy 60Ry. The Bravais lattices is hexagonal, where a = 2.461 Å and the separation between the graphene layers is L = 5a. The ground state electronic density of the bulk Al 2 O 3 is calculated using 9 × 9 × 3 K-mesh, the plane-wave cut-off energy is 50Ry and the Bravais lattices is hexagonal (12 Al and 18 O atoms in the unit cell) with the lattice constants a = 4.76 Å and c = 12.99 Å.
The optical conductivity (18)- (21) in the KC 8 -SL is calculated using a 201 × 201 × 1 K-mash and the band summations (n, m) are performed over 100 bands. The damping parameters are η intra = 10 meV and η inter = 40 meV and the temperature is T = 25 meV. The graphene optical conductivity is calculated using a 601 × 601 × 1 K-mash and the band summations are performed over 20 bands. The damping parameters are η intra = 1 and 15 meV, η inter = 25 meV and the temperature is T = 25 meV. The response function (22) of the Al 2 O 3 is calculated using a 21 × 21 × 7 k-point mesh and the band summations (n, m) are performed over 120 bands. The damping parameter is η = 100 meV and the temperature is T = 10 meV. For the optically small wave vectors q ≈ 0 used in this modeling, the crystal local field effects are negligible, so the crystal local field effects cut-off energy is set to zero. Figure 2 shows the ab initio optical conductivity σ 0 yy in the KC 8 -SL. The intraband contribution σ intra is turquoise shaded, while the interband contribution σ inter is orange shaded. The two pronounced peaks at ω ≈ 4 and 14 eV correspond with the interband transitions between the graphene C(π) and C(σ) bands. The insert in Figure 2 shows the KC 8 band structure, and we can see that the KC 8 band structure does not differ much from the graphene band structure. The only influence of the K adatoms is the appearance of the potassium parabolic K(σ) band (turquoise dashed lines show the parabolic fit of the K(σ) band for the effective mass m * = 0.92) which abundantly donates electrons to the graphene C(π) band (denoted by magenta dashed lines) but in the way that it still remains partially filled. This causes the Fermi level shift by 1 eV above the Dirac point so that the onset for the interband transitions between the graphene C(π) bands appear at 2 eV (also denoted by brown dashed line in Figure 2). At the same time, this causes the appearance of two intraband excitations channels, K(σ) and C(π), which appear as strong Drude peak (shaded by turquoise color at ω ≈ 0). Accordingly, this provides large effective number of charge carriers, n yy = 0.021a −2 0 [7.58 × 10 14 cm −2 ], resulting in a very intensive Dirac plasmon with zero direct interband damping.  Figure 3a shows the ab initio macroscopic dielectric function (25) of the bulk Al 2 O 3 crystal. We can see that 1 is almost constant ( 1 ≈ 3) for low frequencies (ω < 3 eV), i.e., in the IR and even in the visible range, while 2 is zero up to the band gap energy (E g ∼6 eV). This suggests that Al 2 O 3 is a good choice for the substrate for the IR plasmonics, since its electronic excitations are far above the IR plasmons, and its IR active SO phonons (at ω TO < 100 meV) [40] are still below the IR plasmons considered here. Therefore, in the frequency range of interest (red frame in Figure 3a), there is no dissipation of the electromagnetic energy in the substrate (it is transparent) and the dielectric function is constant.
Finally, in Figure 3b, we demonstrate the influence of the KC 8 -SL interband transitions and the influence of the dielectric substrate polarization on the Dirac plasmon dispersion relations. The dispersion relations are derived from the maxima of the real part of the screened conductivity for σ 0 = σ intra and Γ = Γ 0 (brown dashed line), σ 0 = σ intra + σ inter and Γ = Γ 0 (orange dashed-dotted line) as well as σ 0 = σ intra + σ inter and Γ = Γ 0 + Γ sc (solid turquoise line).

Results
In this section, we first demonstrate how replacing the KC 8 -SL by the KC 8 -NR of various widths d influences the absorption spectra A. Then, we demonstrate that the series of the Dirac plasmon resonances in the KC 8 -NR when unfold in ExBZ resembles the Dirac plasmon in the KC 8 -NR. Finally, we present the DPR band structure within the 1stBZ. We focus on the e =ŷ (perpendicular to the NR) polarized electromagnetic field, and the separation between the graphene planes and the dielectric surfaces is fixed to z 0 = 3.0 Å.
The blue lines in Figure 4 show the absorption spectra (12)  We can see that, after cutting the KC 8 -SL into nanoribbons, the Drude asymmetric tail transforms into a series of IR-active DPR n = 1, 3, 5, 7, . . . with the energy depending on the nanoribbon width d. As expected, as d increases, the energy of the DPR decreases and the energy difference between them becomes smaller. Considering that the sample is driven by the electromagnetic field, homogeneous in the y direction, the peaks appearing in the absorption spectra obviously represent optically active dipolar modes. According to the continuity equationρ ind = − ∂ ∂y j ind y and j ind = σ ⊗ E, where E is the external field given by Equation (1), the induced density can be calculated from ρ ind ∼ ∂σ yy /∂y. Indeed, Figure 5a, which shows the induced electronic densities at frequencies ω corresponding to the absorption peaks n = 1, 3, 5, 7 in Figure 4a, clearly demonstrates their dipolar character.   On the other hand, the induced densities of the other non-dipolar modes n = 2, 4, 6, . . . are zero, so they represents the dark modes, not visible in the absorption spectrum. For the period chosen here (l = 2d), the separation between the nanoribbons d is quite large, and the interaction between the dipolar modes in the different nanoribbons is negligible; thus, they can be considered as almost decoupled resonances of the individual nanoribbons. However, we show below that these modes are still weakly dispersive as we increase the wave vector Q y , suggesting their small but finite interaction.   Figure 6d. The nanoribbon period is again l = 2d. The red doted lines denote the energies of the IR active modes (n = 1, 3, 5, . . .), as they appear in Figure 4, while the turquoise doted lines denote the energies of the dark modes (n = 2, 4, 6, . . .). Green dotted lines denote the dispersion relation of the Dirac plasmon in the KC 8 -SL. We can see that the principal mode n = 1 is the most dispersive one (within each BZ) and the most intensive in the first two BZ 0 < Q y < 4π/l. The rest of the modes, n = 2, 3, 4, 5, . . ., are less dispersive (flat patterns within some of the exBZs), and they are the most intensive through the few extended BZ, but precisely in the way that resembles the DP dispersion relation. This is particularly noticeable for the larger widths d, as can be seen by comparing Figure 6c,d. Moreover, we can see that the modes n = 3, 5, . . . in the extended BZ are 'folded' (although with much lower intensity) into the radiative region (Q y ≈ 0) where they become IR active. This is clearly noticeable, e.g., for the n = 3 mode in Figure 6a-c, where just a small fraction of that mode, in exBZ, 'projects' into the radiative region (Q y = 0). This is in accordance with the results in Figure 4, where just the principal mode n = 1 is the most intensive, while the higher modes n = 3, 5, 7 . . . are significantly suppressed. These results show us that we can fold the fragments of the Dirac plasmon (at Q y ≈ 2π(2k + 1)/l; k = 0, 1, 2, . . .) into the radiative region and make them IR-active resonances, in a controlled way, by cutting the KC 8 -SL into an array of nanoribbons, as sketched in Figure 5b,c. For example, by changing the nanoribbon width, we can choose which part of the Dirac plasmon in the KC 8 -SL we want to fold into the radiative region. It should be noted that this procedure is also valid in the opposite direction, starting from a single nanoribbon antenna. A single nanoribbon supports nondispersive and localized (both IR active and dark) plasmon resonances, but, after the nanoribbons are arranged in a lattice, the plasmon resonances unfold over the extended BZ resembling the DP, regardless of the period l. All these manipulations are experimentally feasible, which could have a significant impact on the applied plasmonics.

DPR Band Structure
We now analyze the dispersion relations of the DPR within the 1stBZ, Q y ∈ [−π/l, π/l], i.e., the DPR band structure. In the previous examples, the separation between the nanoribbons is quite large, so the coupling between the DPR in the different ribbons is weak and, consequently, the DPR are weakly dispersive within the 1stBZ. However, we show above that the dispersion across the exBZs is strong so that it resembles the Dirac plasmon in the KC 8 -SL. The origin of this dispersivity is simple: for larger wave vector Q y , the spatial variation of the external field partially or fully fits the spatial variation (or symmetry) of higher excited modes n = 2, 3, 4 . . ., regardless of whether they are dipolar or non-dipolar modes, and finally it efficiently excites these modes. On the other hand, the homogenous external field (Q y = 0) selectively and quite inefficiently excites the higher dipolar modes n = 3, 5, 7, . . . Therefore, the inter-zonal dispersivity always exists, even though the interaction between the ribbons is negligible. However, for smaller separations, the interaction between the nanoribbons is getting stronger and DPRs become dispersive within the 1stBZ (or within some individual exBZ). Figure 7 shows the absorption intensities A in the KC 8 -NR array deposited on the Al 2 O 3 surface for different wave vectors Q y ∈1stBZ. The nanoribbon width is d = 200 nm and the period is l = 220 nm. The photon dispersion ω = Qc is also shown (green dotted lines) in order to denote the radiative region ω > Qc. It can be noticed that this small separation (l − d = 20 nm) causes substantial dispersivity of the principal n = 1 mode, resulting with the band width of about W = 80 meV and the band gap opening of about E g = 60 meV. The higher DPR (n = 2, 3, 4 . . .) are less dispersive, probably because they produce short-ranged electric field so the inter-ribbon interaction is weaker. It is interesting that the dark mode n = 2 seems to split into two branches as Q y decreases. This may be the evidence of the surface or 'edge' plasmons localized at the KC 8 -NR boundaries [17]. The edge plasmons are the counterparts of the standard surface plasmons appearing on metal surfaces. They are the extra solutions of the Maxwell's equation, due to the symmetry breaking caused by the edge, and have an evanescent character, in contrast to the DPR which oscillate across the nanoribbon (see Figure 5a). The DPR band-gap can be manipulated by changing the KC 8 -NR parameters which opens the possibility for trapping the photons in the principal n = 1 band and achieving the Dirac plasmon Bose-Einstein condensate. Therefore, the doped graphene nanoribbons enable direct light-plasmon interaction, which can be exploited in many plasmonic or photonics applications, while at the same time it can serve as a polygon for exploring fundamental physical phenomena such as strong light-matter interactions, which has been intensively explored recently [41].

Comparison with Experiment
In order to verify the accuracy and experimental feasibility of the above results, we compare them with some recent experimental results. Since optical absorption experiments on the KC 8 -NR arrays still do not exist, we compare our results with the experimental results for the differential transmission ∆T = T − T CNP through the array of doped graphene micro-ribbons on Si/SiO 2 substrate [15], where T CNP is the transmission coefficient through the device at the charge neutral point (CNP). ∆T is directly related to our infrared absorption spectrum A. In our calculations, the graphene is doped by holes, where the hole concentration is chosen to be 1.5 × 10 13 cm −2 [15] (E F = −0.374 eV with respect to the Dirac point). The Si/SiO 2 substrate is described by the dielectric constant s = 6.5, which is between 3.8 in SiO 2 and 11.7 in Si. The separation between the graphene and the SiO 2 surface is taken to be z 0 = 4 Å [42].
In order to gain insight into the measured data for wider energy range, in Figure 8a, we compare the experimental result for ∆T (blue circles) with our results for A calculated for two intraband damping parameters η intra = 1 meV (brown line) and η intra = 15 meV (turquoise line) and on extended frequency scale. The graphene ribbon width is d = 1 µm and the period is l = 2 µm. We can see that for the smaller damping parameters absorption spectra shows DPR n = 1, 3, 5, . . ., which for the larger damping smooth out into an asymmetrical lineshape which is in excellently agreement with the experimental data. Therefore, we can conclude that the experimental lineshape mainly consists of the principal dipolar mode n = 1, and its asymmetry is a consequence of the excitations of higherorder dipolar modes n = 3, 5, . . . Figure 8b shows the theoretical absorption spectra A in an array of graphene ribbons of widths d = 1 µm (blue line), 2 µm (red line) and 4 µm (green line) and compares them with the experimental results for ∆T (blue, red and green circles). The period is l = 2d. These results undoubtedly confirm that the broad experimental peaks corresponds to the n = 1 DPR, while the higher-order dipolar DPR n = 3, 5, 7, . . . give the spectrum an asymmetric shape. Moreover, these results determine the natural intraband damping parameter η intra = 15 meV, which is a consequence of the electron-phonon (instrinsic and SO phonons) interactions, scattering on impurities and other crystal imperfections.  All this suggests that the electron-phonon interaction (or maybe some other scattering mechanisms) is likely to play a significant role in profiling the higher-order plasmon resonances n = 3, 5, . . . . Below, we show that the synthesis of the potassium intercalated graphene (KC 8 ) ribbons is indeed possible and explore how the potassium adatoms influence the strength of the electron-phonon coupling. The latter is very important because alkali metals can sometimes increase and sometimes decrease the strength of the electron coupling to the graphene E 2g phonon [43], which is, as already mentioned, very important for the damping of the plasmon resonances.
Our scanning electron microscopy (SEM) analysis of GNRs placed on a carbon tape and analyzed with 5 kV Helios NanoLab DualBeam scanning electron microscope showed multilayer GNR structures with lengths of several microns and widths ranging around 100 nm. Figure 9a confirms a flat-rippled multilayer nanoribbon morphology. These GNRs were synthesized via CVD following the procedure described in [44] and further intercalated by conducting a two-zone vapor transport method, as described in [45]. The intercalation compounds was obtained by the combination of an alkali metal (K) placed in a glass vial with GNRs sealed under high vacuum conditions at 10 −6 mbar in a proportion of 3.2 mg of GNRs per 1.3 mg of potassium (∼KC 8 GNRIC). The characteristic Raman spectrum from the GNRs (Figure 9b-i) revealed the characteristic D-band and G-band located at ∼1338 cm −1 and ∼1574 cm −1 , respectively. The D-band exhibited a larger intensity caused by the edge ripple proportion, while the D/G ratio was found to be 1.25 characteristic of graphene nanoribbons [44]. At ∼2674 cm −1 , we observed the 2D-line characteristic of graphitic GNRs [44,46]. An intercalation process using the discussed pristine sample was performed obtaining a KC 8 GNRIC. This sample was kept under vacuum conditions to avoid oxidation during the Raman measurements. The Raman spectrum obtained from the KC 8 GNRIC in Figure 9c shows the characteristic broad Fano-line-shape composed by a G-line at ∼1505 cm −1 caused by the intercalation of potassium layers in between the graphene ribbons. This characteristic G-line in intercalation compounds originates from strong electron-phonon coupling (EPC) interactions existing between the potassium atoms and the graphene layers, as we reported previously for graphite intercalation compounds [47]. It is proven that by doping graphene with electron donors the Dirac plasmon resonances can be excited [18]. Thus, here we could introduce the fact that, by obtaining a highly e − -doped intercalation compound (i.e., KC 8 ), we must obtain: (i) a strong EPC that will be responsible for superconductivity in stage I GICs according to the BCS theory directly related to the G-line phonon frequency and to the adiabatic (ω A ) and non-adiabatic (ω N A ) phonon frequencies of the GIC; (ii) a strong EPC in GNRIC will serve to excite plasmonic resonances derived from the perpendicular electric field in those nanostructures. To estimate the renormalized electron-phonon scattering line width (Γ EPC ) in GNRIC, we consider the G-line phonon frequency (ω G from E 2g ) of the Raman spectrum in Figure 9c in the following equation [47]: where ω G is the measured G-line frequency from the Fano function (1505 cm −1 ), ω A is the adiabatic phonon frequency (1223 cm −1 ) and ω N A corresponds to the non-adiabatic phonon frequency (1534 cm −1 ) [48]. From this equation we obtained, Γ EPC =243 cm −1 for KC8 GNRIC is indicative of a potential superconducting behavior of the material as it behaves linearly with the measured FWHM (Γ ph =262 cm −1 ). An individualized GNRIC can be evinced in Figure 9d to confirm no further damage to the structure of the graphitic ribbon. Finally, our theoretical results are in excellent agreement with the experimental results, confirming the credibility of the presented method and the above-stated conclusions. The Fano-like line shape derived from the strong coupling between K atoms and the graphene is a fine characteristic of a stage I intercalation compound; (d) Individualized GNRs after potassium intercalation. The structure and shape of the GNRs was not affected after the intercalation as observed in the micrograph. GNRs were confirmed to be longer than 10-20 microns with a width ∼100-200 nm.

Conclusions
We developed the ab initio theoretical formulation of the electromagnetic response in doped graphene nanoribbons and used it to calculate the optical absorption in an array of potassium-doped graphene nanoribbons deposited on a dielectric Al 2 O 3 surface. We demonstrated that the replacement of the single-layer doped graphene by the graphene ribbons of different period l causes the 'projection' of the Dirac plasmon into the radiative region, turning it into a series of IR-active Dirac plasmon resonances. This encourages the fabrication of graphene nanoribbons with the desired electromagnetic response in the IR or THz frequency range, which could be used in plasmonic, photonic or optoelectronic applications. We showed that the DPR band structure (band gap E g and band width W) can be tuned by changing the period l. By creating a large band gap E g , one can enable trapping of the photons in the principal n = 1 band and achieve the Dirac plasmon Bose-Einstein condensate. Therefore, the graphene ribbons can be exploited in many plasmonic or photonic applications, but at the same time they can serve as a polygon for exploring fundamental physical phenomena such as strong light-matter interactions.