Semiempirical modeling of the effects of the intrinsic and extrinsic optical phonons on the performance of the graphene-based devices

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 correctly predicting the behavior of graphene-based plasmonic devices.

www.nature.com/scientificreports/ dielectrics 13,19 . A substrate like SiC, hBN, or SiO 2 , with one or more (n) dominant vibrational S-POPs modes, 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. Notice, the semi-classical spatial dispersion model tolerable results for frequencies far from vibrational excitations in which the phonon effects are unimportant. 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,19 (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. One of the main limitations of the built-in graphene model in CST is that it is written based on the intraband and interband equations and does not include the effects of vibrational modes. However, we modify the model to the impacts of intrinsic and extrinsic phonons. 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 two-dimensional (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 19 . Here, we take in a 1D periodic array of infinitely long (in the y-direction) graphene ribbons of width w x and periodicity d x (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 19 . 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 SiO 2 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. Figure 1 shows the schematic of a typical structure containing a 1D periodic array of graphene ribbons of width w x and periodicity d x , placed on a nonpolar (polar) dielectric grown on an n ++ -Si together with the aluminum layer, acting as the ohmic contact to the gate material (the dielectric layer). As observed in the figure, by applying an appropriate external bias (V G ) between the gate terminal and the aluminum contact on graphene, interconnecting the graphene ribbons, one can control the graphene chemical potential, μ C . Letters m label the different orders of the transmitted and scattered waves in arbitrary directions depicted by red arrows.  www.nature.com/scientificreports/ k xm > k 0 , k zm 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 k xm > k 0 is valid only for m ≠ 0, meaning that only the zeroth-order scattered light can propagate 20 . 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 20 to reach the extinction expression for the 1D array,

Theory
where η 0 = 376.73 is the free space impedance and Z G,s(p) represents the impedance of a 1D periodic graphene array illuminated by an s(p)-polarized wave 21 and where and R G is the graphene impedance that we will describe in the following sub-section. (3) and (4), the extinction and impedance magnitudes (α 1D and R G ) 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 the 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 22 , www.nature.com/scientificreports/ in which N phj 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 k B 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. The value of the fit parameter strongly depends on the underlying dielectric type -i.e., whether it is nonpolar or polar, and the number of POPs oscillations for the polar material. For the fit to be practicable, we need some reliable experimental data to compare. Hence, we introduce the fit parameter (β) solely for DLC and SiO 2 , for which experimental data are available 13 . In other words, the model is a semiempirical one that applies to any nonpolar/polar dielectric material for which experimental data are available. 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, other groups 14,15 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.2 μm.

Results and discussion
Graphene on a non-polar dielectric. First, we considered arrays of graphene ribbons of widths (50 nm ≤ w x ≤ 130 nm) and pitches d x = 2w x placed on a non-polar (DLC) substrate of relative permittivity ε DLC = 5.7-i.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 Fig. 2a 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 Fig. 2a. 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 frequencies-i.e., less than a second-whereas 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 w x and d x , exhibits a dominant resonance peak, decreasing as w x (d x ) 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 w x 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 decreases-decreasing 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.
In 1-D periodic arrays, the plasmons' wavenumber q pl ≈ π / (w x− w 0 ), in which w 0 ≈ 28 nm is the electrically inactive width of structured graphene 13 . The open circles and crosses ( ×), in Fig. 2b, show the dependence of the extinction peak's frequency versus the surface plasmons wavenumber, q pl , obtained from (3) and the experimental data 13 , as shown in 2(a). The dashes represent the dispersion relation obtained for the similar graphene ribbons arrays on nonpolar substrates 5 , www.nature.com/scientificreports/ where e is the elementary charge, ε 0 is the free-space permittivity and ε env is the average dielectric constant of the upper and the lower media surrounding the graphene array, and with q′ = 2πε 0 ε env (ℏω op ) 2 /e 2 μ C as the plasmons wavenumber at the frequency ω = ω op , vanishing at q pl ≪ q′ -i.e., γ(q pl ≪ q′) → 0. Figure 2b demonstrates how accurately the plasmonic peaks obey the dispersion relation. Another critical factor in determining the extinction peak frequency is the graphene's chemical potential for a given array. According to (7a), ω pl ∝ √ µ C , conveying an increase in μ C causes a blueshift in peak frequency.
So, we have considered an array of widths w x = 100 nm and calculated the extinction spectra while varying the chemical potential in the range of 0.3 eV ≤ μ C ≤ 0.5 eV. Figure 3a illustrates the calculated spectra, using (3). Then, we plotted the extinction peak frequency versus the chemical potential for all seven arrays of Figs. 2 in 3b, demonstrating a more transparent presentation of the blueshifts caused by the increase in μ C . As seen from this figure, one can conclude that the narrower the width of the ribbons, the larger the blueshift caused by an increment in μ C . This figure gives us a piece of valuable information about how to choose a design parameter. Considering a CO 2 laser of operating frequency ~ 28.5 THz (0.117 eV) 19 as the light source, the best choice for w x of the graphene ribbons with chemical potential in the range of μ C ~ 0.3-0.4 eV would be 80-100 nm.
Graphene on the polar dielectric. Electrons in graphene placed on polar substrates are exposed to both intrinsic and extrinsic phonons. Besides the effects of intrinsic phonons mentioned in the previous subsection, extrinsic phonons cause mode splitting around their vibrational frequencies and turn single plasmonic mode into n + 1 coupled plasmon-phonon modes, where n is the number of the extrinsic vibrational frequencies 5 .
Here, we consider two arrays of graphene ribbons of w x = 240 and 160 nm and μ C = 0.6 and 0.45 eV, placed on SiO 2 slabs -i.e., a polar dielectric with n = 3 and ħω 1, 2, 3 = 100 meV (24.63 THz), 134 meV (32.68 THz), and 144.8 meV (35.32 THz), just the same as those used in two independent experimental studies 13,19 . The red and green dashes in Fig. 4 represent the numerical results obtained from (3), fitted to the experimental results depicted by crosses (blue 13 and magenta 19 , 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, N ph in (6) for this case is the total number of phonons provided for all vibrational frequencies.
(7b) www.nature.com/scientificreports/ The labels m 1 − m 4 , 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 m 3 . 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, m 3 , peaks significantly. In other words, the mode m 3 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.  www.nature.com/scientificreports/ Next, we further compared the model with the experimental results of 13 , varying the widths of the graphene ribbons in the range of 60 ≤ w x ≤ 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 (0.041-0.41 eV) (Fig. 5a) and the dependencies of the corresponding peak's frequency on the plasmons wavenumber, q, as depicted in Fig. 5b. The flat parts of the extinction spectra in Fig. 5a 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 w x , in agreement with the experimental results of 13 . Those features are (i) as w x decreases, 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 m 1 to m 2 and then to m 4 , before the resonance peak for m 4 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 m 4 increases as ω m4 → ω op and beyond, due to the damping through intrinsic optical phonon emission. Figure 5b shows the dependencies of the extinction peaks frequencies on the plasmons wavenumber that inversely varies with the ribbons widths-i.e., q pl ∝ (w x− w 0 ) −1 . The open magenta circles and blue crosses represent the peak 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. A 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 frequencies of hybrid modes and phonons 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 Fig. 6). 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 w x decreases, for a given chemical potential, varying from one mode to another. Moreover, a common www.nature.com/scientificreports/ feature observed for the resonant peaks m 2 , m 3 , and m 4 is not true for m 1 , 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 m 4 in Fig. 5a and further observation into the inherent phonon damping, the impact of optical phonons at mode damping becomes clearer. Figure 7 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 the 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, SiO 2 , 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, which the experimental data and CST simulation confirm.

Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the corresponding author upon reasonable request.