Effects of Intrinsic and Extrinsic Surface Optical Phonons on the Performance of the Graphene based Devices: an Accurate and Straightforward Modeling

Surface plasmons in graphene have mainly been affected by intrinsic optical phonons due to the vibrations of the carbon atoms and surface polar optical phonons (S-POPs) of the underlying dielectric surface. This plasmon hybridization dramatically changes the features of the plasmonic devices. However, a complete theoretical model for the graphene impedance to consider the optical phonons effects is yet remained to be developed. Here, we show how to derive a model for graphene impedance to include such impacts on graphene surface plasmons. Earlier models suffer from two limitations — i.e., the inability to show (i) the transformation of a single pure plasmonic mode into multiple hybrid plasmon-phonon excitations and (ii) the damping effect for energies beyond that of the intrinsic optical phonons due to the phonon emission. Our new model overcomes these two limitations. Then, we calculate the extinction spectra for a one-dimensional periodic array of graphene ribbons obtained through the impedance boundary condition method, addressing these obstacles. These spectra are directly related to graphene impedance, modeled using the dielectric function we developed in our earlier work. The extinction spectra show the presented model overcoming the limitations, firmly fitting the experimental data reported by others. Furthermore, we introduce our developed model for graphene to the CST Studio software to verify the accuracy of our extinction relation and impedance model. This study can be a step forward to correctly predicting the behavior of graphene-based plasmonic devices. Electron-phonon interaction is an inevitable source of scattering, strongly affecting carrier transport. In particular, plasmon dispersion in graphene has been altered significantly by intrinsic and extrinsic optical phonons exposure. Plasmons with energies higher than intrinsic optical phonons (≈ 0.2 eV) damp through inherent phonon emission into the intraband single-particle excitation region. Furthermore, the extrinsic polar optical phonons of the surface (i.e., S-POPs) of the underlying polar-dielectric substrate significantly affect the dispersion of the surface plasmons of graphene. This interaction turns the primary plasmonic mode into several hybrid plasmon-phonon modes, depending on the number of S-POP excitation modes on the substrate surface, as we have demonstrated earlier. Electron-phonon interactions enhance the graphene impedance, altering its extinction spectra. In other words, electron scatterings result in loss of the energy and momentum of the carrier and hence, hinder the electronic response to an applied electric field due to the impedance enhancement. The scattering rate is directly associated with the number of optical phonons. Hence, to examine the impedance properly, the impact of the atomic vibrations must be considered, particularly at room temperature when the number of phonons is high enough. The resistance formulas for graphene proposed earlier, obtained through interband and intraband transitions, using electrical conductivity, do not include phonons impacts. Furthermore, to theoretically predict the measured light transmittance through a sample including the resistive loss, one requires an accurate impedance model that can lead to the knowledge of the exact extinction quantities, foreseeing the empirical characteristics of plasmonic devices. The existing models for graphene resistance suffer from two main limitations. They cannot predict the plasmon damping phenomenon caused by the intrinsic optical phonon emission, demonstrated by. Moreover, they cannot foresee the impacts of the extraneous S-POPs on the plasmon excitations for graphene loaded on a polar dielectric— like SiC, hBN, or SiO2, having one or more S-POPs excitation modes. A substrate with a few dominant vibrational modes (n) splits a single graphene plasmonic mode into n+1 hybrid plasmon-phonon modes. This effect is prevalent near the splitting energies, becoming negligible far from those values. So long as the working frequency is near the phonon’s frequency, we should worry about the correctness of the impedance model in predicting the exact output extinction ratio. This work aims to overcome these discrepancies by introducing a new model for graphene impedance, taking the presence of intrinsic optical phonons and extrinsic S-POPs into account. In general, the impedance is directly related to two values: (i) The loss function, obtained by the imaginary part of the dielectric function, can be developed by random phase approximation (RPA) — i.e., Im{εRPA}−1 ; (ii) The number of intrinsic optical phonons and extrinsic S-POPs. Our accurate and straightforward model includes these parameters. We also acquire an exact relation for evaluating extinction values, benefiting from the impedance boundary condition method, showing an excellent agreement between theory and experiment by utilizing the modeled impedance. Our new model indeed predicts the damping of hybrid modes for energies higher than that of the intrinsic optical phonons and the mode splitting caused by the extrinsic S-POPs. Moreover, we define our impedance model for graphene in CST Microwave Studio software and obtain the desired extinction values. The results show that the presented formulas work rightly, paving the way for accurate prediction of the experimental results. Different approaches are employed to excite plasmons in graphene, such as patterning into one or twodimensional (1D or 2D) periodic arrays, like ribbons and discs, light scattering from structures adjacent to graphene, or coupling to a grating placed on top or below the graphene. The first approach is the most straightforward based on today’s technological progress. In this approach, by adjusting parameters like periodicity, width, or the spaces between nearby arrays, one can quickly obtain desired plasmonic peak frequency. Here, we take in a 1D periodic array of infinitely long (in the y-direction) graphene ribbons of width wx and periodicity dx (along the x-direction). Moreover, we obtain an extinction relation when a normal s-polarized (i.e., polarization perpendicular to the ribbons) lightwave illuminates periodic plasmonic patterns, exciting surface plasmons on the graphene surface. A p-polarized light (i.e., polarization parallel to the stripes) cannot recognize the periodicity and hence does not excite surface plasmons. We use the impedance boundary condition method to reach the desired extinction relation and apply the RPA dielectric function to introduce desired impedance. The organization of the rest of the paper is as follows. In the next section, after a brief introduction of Green’s function for an array of ribbons (i.e., the potential response produced by line sources generating the stripes), we obtain the extinction ratio formula. Then, we model the graphene impedance in the presence of the optical phonons. Next, we present the numerical results, including extinction spectra of two periodic plasmonic arrays, one placed on a diamond-like-carbon (DLC) substrate — i.e., non-polar substrate —and the other on SiO2 as a polar substrate with three dominant vibrating modes. Then, we compare the results of our theoretical model with the existing experimental results, showing perfect agreement among them. Moreover, we present the outcomes of software simulations to demonstrate our model’s correctness. The material introduced in CST software comes in Supplementary Information. Finally, we complete the paper with conclusive remarks. Theory Extinction relation. Green’s function can represent the solution of Maxwell’s equations in 1D periodic structures, such as periodically arrayed graphene ribbons, considering a line source shifted by an appropriate phase,     1D , , , 2 zm

can lead to the knowledge of the exact extinction quantities, foreseeing the empirical characteristics of plasmonic devices.
The existing models for graphene resistance suffer from two main limitations. They cannot predict the plasmon damping phenomenon caused by the intrinsic optical phonon emission, demonstrated by [12][13][14][15][16][17][18] . Moreover, they cannot foresee the impacts of the extraneous S-POPs on the plasmon excitations for graphene loaded on a polar dielectric-like SiC, hBN, or SiO2, having one or more S-POPs excitation modes. A substrate with a few dominant vibrational modes (n) splits a single graphene plasmonic mode into n+1 hybrid plasmon-phonon modes. This effect is prevalent near the splitting energies, becoming negligible far from those values 5,13 . So long as the working frequency is near the phonon's frequency, we should worry about the correctness of the impedance model in predicting the exact output extinction ratio.
This work aims to overcome these discrepancies by introducing a new model for graphene impedance, taking the presence of intrinsic optical phonons and extrinsic S-POPs into account. In general, the impedance is directly related to two values: (i) The loss function, obtained by the imaginary part of the dielectric function, can be developed by random phase approximation (RPA)i.e., Im{εRPA} −1 5,13,18 ; (ii) The number of intrinsic optical phonons and extrinsic S-POPs. Our accurate and straightforward model includes these parameters. We also acquire an exact relation for evaluating extinction values, benefiting from the impedance boundary condition method, showing an excellent agreement between theory and experiment by utilizing the modeled impedance. Our new model indeed predicts the damping of hybrid modes for energies higher than that of the intrinsic optical phonons and the mode splitting caused by the extrinsic S-POPs. Moreover, we define our impedance model for graphene in CST Microwave Studio software and obtain the desired extinction values. The results show that the presented formulas work rightly, paving the way for accurate prediction of the experimental results.
Different approaches are employed to excite plasmons in graphene, such as patterning into one or twodimensional (1D or 2D) periodic arrays, like ribbons and discs, light scattering from structures adjacent to graphene, or coupling to a grating placed on top or below the graphene. The first approach is the most straightforward based on today's technological progress. In this approach, by adjusting parameters like periodicity, width, or the spaces between nearby arrays, one can quickly obtain desired plasmonic peak frequency 18 . Here, we take in a 1D periodic array of infinitely long (in the y-direction) graphene ribbons of width wx and periodicity dx (along the x-direction). Moreover, we obtain an extinction relation when a normal s-polarized (i.e., polarization perpendicular to the ribbons) lightwave illuminates periodic plasmonic patterns, exciting surface plasmons on the graphene surface. A p-polarized light (i.e., polarization parallel to the stripes) cannot recognize the periodicity and hence does not excite surface plasmons 18 . We use the impedance boundary condition method to reach the desired extinction relation and apply the RPA dielectric function to introduce desired impedance.
The organization of the rest of the paper is as follows. In the next section, after a brief introduction of Green's function for an array of ribbons (i.e., the potential response produced by line sources generating the stripes), we obtain the extinction ratio formula. Then, we model the graphene impedance in the presence of the optical phonons. Next, we present the numerical results, including extinction spectra of two periodic plasmonic arrays, one placed on a diamond-like-carbon (DLC) substratei.e., non-polar substrate -and the other on SiO2 as a polar substrate with three dominant vibrating modes. Then, we compare the results of our theoretical model with the existing experimental results, showing perfect agreement among them. Moreover, we present the outcomes of software simulations to demonstrate our model's correctness. The material introduced in CST software comes in Supplementary Information. Finally, we complete the paper with conclusive remarks.

Theory
Extinction relation. Green In (2a) and (2b), kx(y) and k0 = 2π/λ represent the x(y)-component of the incident wavevector and the related wavenumber, with wavelength λ. Under a normal s-polarized illumination kx(y) = 0. So, according to 2(b), for kxm > k0, kzm becomes imaginary, degrading the Green's function rapidly in the far-field region, z ≫ z′, making the reflected and transmitted lights invisible. For the sub-wavelength unit cell, the condition of kxm > k0 is valid only for m ≠ 0, meaning that only the zeroth-order scattered light can propagate 19 . Assume an array is positioned in an x-y plane at z′ = 0, ignore the terms m ≠ 0, and use the impedance boundary condition method 19  1 cos and RG is the graphene impedance that we will describe in the following sub-section. (3) and (4), the extinction and impedance magnitudes (α1D and RG) are interrelated. So, to estimate the extinction coefficient accurately, we must genuinely predict the impedance magnitude. As we mentioned earlier, we have recently developed 5 exact plasmon-phonon dispersion relations for graphene placed on polar substrates, using dielectric function expanded by RPA (See Supplementary Information for RPA). Loss function, related to the Im{εRPA} −1 , directly relates to the graphene impedance. The behavior of loss function depends on the intrinsic optical phonons and extrinsic S-POPs 5 , which means the impedance also depends on the numbers of the inherent optical phonons and extraneous vibrational modes. So the impedance takes its value based on the type of the underlying dielectric substrate. In other words, the graphene impedance relation linearly depends on the total number of optical phonons obtained by the Bose-Einstein distribution 21 ,

Impedance Relation. Based on
in which Nphj and ћ phj stand for the number and energy of phonons of the j-th excited mode (j = op for the intrinsic phonons and 1, 2, 3 for the extrinsic S-POPs), ωphj is the related radian frequency, ћ is the reduced Planck's constant, and kB and T are the Boltzmann constant and the ambient temperature, respectively. At low temperatures, the number of optical phonons is small enough not to participate in plasmon hybridization and damping. The phonon's association starts to be significant as the temperature increases and, then, the impedance enhances due to the increased electron-phonon scattering. Therefore, the impedance directly links to the total number of phonons too. Hereafter, we work at room temperature, T = 300 K. Knowing all these, we introduce a new formula for the graphene impedance that is directly proportional to the loss factor and the total number of phonons, in which RPA follows (S9) and (S10) in Supplementary Information for non-polar and polar dielectric substrates, and β (m −2 ) is a fit parameter, providing a perfect match between the experimental and theoretical extinction data.
To arrive at (6), we carefully examined all parameters that the impedance depends on, using experimental results reported earlier 6,13 , whereas previous studies solely considered the dependency of the impedance on either the number of phonons 6 or the loss parameter 13 . Yet, another group 14 modeled the impedance using the Drude model, ignoring the phonons' effect that is particularly crucial for operating near or above the phonon frequencies. Hence, to the best of our knowledge, (6) is the most accurate and straightforward equation describing the impedance that can be used to predict the behavior of the graphene-based plasmonic devices operating below the mid-infrared wavelength of ~ 6.9 μm.

Results and Discussion
Graphene on a non-polar dielectric. First, we considered arrays of graphene ribbons of widths (50 nm ≤ wx ≤130 nm) and pitches dx=2wx placed on a non-polar (DLC) substrate of relative permittivity DLC = 5.7i.e., similar structures used in an experimental study by [13] and calculated the extinction spectra, using (3) - (6) and (S9) with β = 5.29×10 14 m −2 , to have the best fit to the experimental data reported by [13]. The red dashes and blue cross signs (x) in Figure 1(a) depict our numerical results and the experimental data 13 for graphene ribbons of the given widths and chemical potentials, μC = 0.3 eV. Then, introducing graphene as a material with impedance (6) in the CST Microwave Studio software, we calculate the reflection, R, and absorption, A, coefficients for the array of graphene ribbons on the dielectric, which results in the extinction spectrum through (A+R), as shown by the green dots in Figure 1(a). For more information about using CST built-in BASIC interpreter (VBA Macro Language) for graphene as a material, see Supplementary Information. The same approach is also applicable to other commercial software. This figure demonstrates an excellent agreement between all three sets of data, confirming the correctness of our developed theory. Besides being a compact formula, Eq. (3) requires a much shorter time for calculating the extinction spectrum over a wide range of frequenciesi.e., less than a secondwhereas CST simulation for the desired frequency range consumes a much longer time. As shown in the figure, the extinction spectrum for each array, with given wx and dx, exhibits a dominant resonance peak, decreasing as wx (dx) becomes narrower. Moreover, the corresponding spectral width becomes thicker, experiencing a blueshift. One may attribute these to the increase in the damping of surface plasmons. Various scattering mechanisms like electron scattering by background ionized impurities, optical phonons, the edge of the arrays, and electron-hole generation within the single-particle excitation region may dampen the plasmons, becoming less probable as wx increases. Background scattering highly depends on the fabrication process, and for a well-made graphene layer, the electron lifetime due to impurities takes the value of 85 fs 13 . The surface plasmons induced on a superstructure are affected by intrinsic phonons. They damp through emitting inherent phonons and enter the phonon emission (PE) continuum if their frequency becomes more significant than op (i.e., the phonon's frequency) 12,13 . For < op, where the phonon emission is not considerable, the edge scattering, resulting in high-energy high-momentum intervalley phonon emission 2 , is the dominant effect among damping processes. So, plasmon peaks lessen as the width decreasesdecreasing width equals heightening edge effects. In addition, plasmons fall within the interband single-particle excitation continuum and decay into electron-hole pairs at frequencies much higher than that of inherent optical phonons'. Our desired frequency range does not include electron-hole excitations.
Another critical factor in determining the extinction peak frequency is the graphene's chemical potential for a given array. According to (7a), pl , pl q   conveying an increase in μC causes a blueshift in peak frequency. So, we have considered an array of widths wx= 100 nm and calculated the extinction spectra while varying the chemical potential in the range of 0.3 eV ≤ μC ≤ 0.5 eV. Figure 2(a) illustrates the calculated spectra, using (3). Then, we plotted the extinction peak frequency versus the chemical potential for all seven arrays of Figure 1 in Figure 2 [13,18]. The red and green dashes in Figure 3 represent the numerical results obtained from (3), fitted to the experimental results depicted by crosses (blue 13 and magenta 18 , respectively), using the fit parameter β = 3.063×10 12 m −2 . The green and red dots in this figure depict the corresponding numerical results obtained by CST simulations, using Eq. (6) with the dispersion relations for the intrinsic and the three relevant extrinsic phonons modes we have already developed in [5]. Moreover, Nph in (6) for this case is the total number of phonons provided all vibrational frequencies, Figure. 3. Extinction spectra for two sets of graphene ribbons (i.e., w = 240 (160) nm and μC = 0.6 (0.45) eV) placed on SiO2 with surface optical phonons of energies ħω1, 2, 3 = 100, 134, and 144.8 meV. The red (green) dashes show the corresponding results obtained from (3), the blue (magenta) crosses represent the experimental results reported by [13] ( [18]), and the green (red) dots show the CST outcome, using (6). The labels m1−m4 denote the peaks of the four coupled phonons-plasmons modes.
The labels m1−m4, seen in this figure, denote the peaks of the four coupled plasmons-phonons modes for both arrays, with center frequencies ωm1 < ω1, ω1 < ωm2 < ω2, ω2 < ωm3 < ω3 and ω3 < ωm4. A notable feature shown in this figure is the damped peak of the third coupled mode denoted by m3. The third mode starts to rise at a frequency greater than ω2. However, the proximity of the S-POPs frequencies ω2 and ω3 prevent the third hybrid mode, m3, peaks significantly. In other words, the mode m3 is too weak to be practically identified as an exciting hybrid mode. Here again, the closeness of the results from our theoretical model to the experimental results confirms the accuracy of the developed model. Nonetheless, the minor deviation between the fit with two different experimental results could be due to a slight dependence of the fit parameter on μC.
Next, we further compared the model with the experimental results of [13], varying the widths of the graphene ribbons in the range of 60 ≤ wx ≤ 240 nm and keeping the chemical potentials constant, μC = 0.6 eV. Then, we calculated the extinction spectra in the frequency range of 10-100 THz (Figure 4(a)) and the dependencies of the corresponding peak's frequency on the plasmons wavenumber, q, as depicted in Figure 4(b). The flat parts of the extinction spectra in Figure 4(a) almost coincided, making them indistinguishable. Hence, to make them distinguishable, we shifted the vertical axis of each curve arbitrarily. The vertical dots in this figure denote the position of the inherent optical phonons frequency, ωop. As seen in this figure, four apparent futures are accompanying the reduction in wx, in agreement with the experimental results of [13]. Those features are (i) as μC increases, each resonance peak for a given array exhibits a blueshift that differs from one to another; (ii) unlike the other three modes, the blueshift exhibited by ωm3 is insignificant because of the closeness of ω2 and ω3, as mentioned earlier; (iii) the weight of the mode intensity transfers from m1 to m2 and then to m4, before the resonance peak for m4 reaches the onset of the PE continuum (ωm4→ωop), beyond which it starts to damp through the emission inherent phonons (a similar phenomenon observed by [5]); (iv) the linewidth of m4 increases as ωm4→ωop and beyond, due to the damping through intrinsic optical phonon emission.  (3), crosses depict the experimental results extracted from [13], and the four dashed lines represent the modes dispersions calculated using (23) of [5]. The red dots depict the dispersion of the SPRs with no optical phonons present. Figure 4(b) shows the dependencies of the extinction peaks frequencies on the plasmons wavenumber that inversely varies with the ribbons widthsi.e., qpl ∝ (wx−w0) −1 . The open magenta circles and blue crosses represent the peaks frequencies obtained from (3) and the experimental results of [13], respectively. The four dashed curves show the dispersion curves obtained from (23) of [5]. Moreover, the red dots show the dispersion of the surface plasmons resonance (SPRs) on the graphene array ignoring the effects of the intrinsic and extrinsic optical phonons. Comparison of these data reveals the coupled modes peaks and SPR frequencies coincide far from splitting energies, where S-POPs have negligible effects, providing only electrons participate in the excitations. Hence, the hybridization of the plasmon-phonon prevails as the hybrid modes' and phonons' frequencies are in proximity.
Finally, we investigated the effects of the graphene chemical potential, μC, on the center frequencies of the extinction peaks (ωm1-4), varying the ribbons widths (see Figure 5). The Horizontal dashes represent the vibrational frequencies of the corresponding S-POPs, as indicated in each part, showing the relations ωm1 < ω1, ω1 < ωm2 < ω2, ω2 < ωm3 < ω3 and ω3 < ωm4 remain independent of the ribbons widths and chemical potentials. The proximity of ω2 and ω3 is why ωm3 is the least affected resonant frequency even by the extreme change in μC for the narrowest ribbons. The plots in this figure help a designer easily find the appropriate design parameters (d and μC) to achieve the desired hybrid mode frequency. These plots also confirm the blueshift exhibited by each extinction peak as wx decreases, for a given chemical potential, varying from one mode to another. Moreover, a common feature observed for the resonant peaks m2, m3, and m4 is not true for m1, that the narrower the ribbons in the array, the larger the blueshift exhibited by the resonant peak as μC increases. Furthermore, a variation in the chemical potential has influenced ωm4 the most for the narrowest ribbons. By extracting full width at half maximum (FWHM) of the hybrid mode m4 in Figure 4a and further observation into the inherent phonon damping, the impact of optical phonons at mode damping becomes clearer. Figure 6 depicts the extracted FWHM versus the resonance peak frequency. The vertical denotes the locus of ωop. It is evident that for peaks' frequencies higher than ωop, the resonance peak damps more rapidly due to the inherent phonon emission.

Conclusion
Utilizing dielectric function obtained by random phase approximation and the number of optical phonons, we have modeled the resistance of graphene placed on a non-polar, DLC, and polar, SiO2, dielectric and, then, calculated extinction spectra, which match with the experiment perfectly. Our results demonstrate that the S-POPs convert a single pure plasmonic mode into multiple hybrid plasmon-phonon resonances, which have never been shown theoretically. Moreover, our presented model predicts the damping caused by optical phonon emissions. Our article provided a new formula of extinction, based on the impedance boundary condition method, in which the experimental data and CST simulation confirm it.