## Abstract

We investigate the nonlinear supermodes of surface plasmon polaritons in graphene multilayers with arbitrary number of graphene layers. Apart from the symmetric and anti-symmetric supermodes which exist in linear multilayer graphene waveguides, more asymmetric supermodes emerge in the nonlinear counterparts as the field symmetry is broken. The number of asymmetric supermodes relies largely on the layer number of graphene. There is a certain threshold of field intensity for the emergence of each individual asymmetric supermode. The threshold increases as the incident wavelength or chemical potential of graphene increases. The study may find applications in building all-optical mode converters and switches.

© 2017 Optical Society of America

## 1. Introduction

Graphene is a two-dimensional material composed of carbon atoms, which exhibits unique electronic and optical characteristics [1–3]. In terahertz and far infrared regime, graphene behaves like a metal and can support transverse magnetic (TM) polarized surface plasmon polaritons (SPPs) due to the coupling between electrons and electromagnetic waves. Graphene SPPs possess many unique features such as huge mode localization, low propagation loss, and flexible tunability [4–7]. The previous researches of graphene SPPs mainly focus on the linear optical properties of graphene [8, 9]. Recently, the nonlinear optical properties of graphene have also attracted much attention [10–13], ranging from self-action of SPPs in graphene [14], four-wave mixing in graphene flakes [15] to the generation of subwavelength spatial solitons in graphene arrays [16]. As a fundamental issue, the supermodes considering the nonlinear of graphene have been demonstrated in double-layer graphene waveguide, which can find applications in nonlinear optical couplers [17]. The supermode refers to the guided modes in composite-array waveguides initially introduced to describe the modes in phase-locked arrays of semiconductors lasers [18]. The supermodes in graphene multilayers are yielded by the coupling of individual SPP mode in each single graphene sheet. It has been demonstrated that the number of linear supermodes in multilayer graphene is equal to that of graphene sheets and the profiles of the supermodes are either symmetric or antisymmetric [19].

In this work, we investigate the nonlinear supermodes in graphene multilayers with finite number of layers. In the presence of Kerr nonlinearity of graphene excited by high light intensity, the multilayer structure can support asymmetric supermodes except for the symmetric and anti-symmetric ones. The generation condition of the asymmetric supermodes is explored by using the first-principle numerical simulations. The total number of nonlinear supermodes is related to the layer number of graphene. The dependence of the intensity threshold to yield asymmetric supermodes on the wavelength and chemical potential is also discussed in detail.

## 2. Nonlinear supermodes in graphene multilayers

We start by discussing the nonlinear optical response of graphene. Compared to conventional bulk dielectrics, the optical nonlinearities of which are described by the second- and third-order nonlinear susceptibilities χ^{(2)} and χ^{(3)}, graphene is a two-dimensional conducting material. The nonlinear property of graphene can be described by the total surface current ** J** =

*J**+*

_{L}

*J**= (σ*

_{NL}*+ σ*

_{g,L}*|*

_{g,NL}

*E|*^{2})

**[20], where**

*E*

*J**= σ*

_{L}

_{g}**is the conventional linear surface current with σ**

*E**modeled by Kubo formula [21–23].*

_{g}

*J**= σ*

_{NL}*|*

_{g,NL}**|**

*E*^{2}

**is the Kerr nonlinear surface current with σ**

*E**being the third-order nonlinear surface conductivity which is given by [24,25]*

_{g,NL}_{c}the chemical potential.

*V*=

_{F}*c*/300 is the Fermi velocity and

*c*is the speed of light. Here we choose λ = 10 μm (

*ħ*ω = 0.124 eV) and μ

_{c}= 0.15 eV. Since 3

*ħ*ω/2 > |μ

_{c}|, the third-harmonic light can induce interband transition, which can be neglected here when considering graphene SPPs [26,27]. It should also be mentioned that no second-order nonlinear response exists in graphene due to the centrosymmetric nature of the structure.

Taking the graphene optical nonlinearity into account, we stack monolayer graphene sheet and dielectric medium alternatively, forming a multilayer structure. The embedded dielectric medium is quartz with a relative permittivity of ε* _{d}* = 4 [17]. As shown in Fig. 1, the position of graphene is denoted by

*n*with

*n*∈ [1,

*N*] and

*N*being the total layer number of graphene. The interlayer space between graphene is

*d*= 30 nm. Consider the graphene SPPs propagate along

*z*direction, the nonlinear supermodes can be described by the following eigen equation [28]

*k*

_{0}and η

_{0}are the vacuum wave number and impedance, respectively. ε

*(*

_{r}*x*) denotes the local permittivity distribution along

*x*direction. The equivalent permittivity of graphene is given by ε

*= 1 +*

_{g}*i*σ

*η*

_{g}_{0}/(

*k*

_{0}Δ) [29], with σ

*being the total surface conductivity. Here graphene is treated as an equivalent thin film with a thickness of Δ ≈1 nm. β is the propagation constant and (*

_{g}*H*,

_{y}*E*)

_{x}*is the transverse mode profiles of the nonlinear supermodes.*

^{T}The nonlinear supermodes can be obtained by numerically solving the eigen equation with self-consistent iteration method [30,31]. In order to obtain the nonlinear solutions, we choose the mode profiles of the linear supermodes as the initial values in the iteration. The relative permittivity is then updated in terms of the solved electromagnetic field distribution of the previous step. Repeating the process of solving the electromagnetic field and updating the permittivity, once the difference between two adjacent calculation results of field distribution exceeds the accuracy of 10^{−4}, the iteration will be ceased and the stable eigenvalues can be obtained [32]. The stable eigenvalues are exactly the propagation constants of the nonlinear supermodes and the eigenvectors are the profiles of the supermodes. The electric field intensity of SPPs is given by |*E*|^{2} = |*E _{x}*|

^{2}+ |

*E*|

_{z}^{2}. The peak value of the intensity reads

*I*

_{0}= |

*E*|

^{2}

*, which is utilized to feature the intensity of the supermodes.*

_{max}The propagation constants of the nonlinear supermodes are shown in Fig. 2 as the total number of graphene layers *N* varies. As *N* increases, more nonlinear supermodes could be supported. As *N* is fixed, one sees that the propagation constants reveal remarkable differences for different nonlinear supermodes. Analogous to the linear cases, the mode wavelength and propagation length of the nonlinear supermodes are given by λ_{p} = 2π/Re(β) and *L*_{p} = 1/[2Im(β)], respectively [33]. Generally, the mode wavelength of the supermode decreases, the propagation length increases accordingly.

Due to the nonlinearity of graphene, the propagation constants of the supermodes depend largely on the field intensity. In Figs. 2(a) and 2(b), we choose a low peak intensity of *I*_{0} = 2 V^{2}/μm^{2}, there are *N* supermodes in a graphene multilayer with *N* layers of graphene. The index of the supermodes is labeled by *s* in an ascending order of the real part of the propagation constant. The situation is similar to that of linear case [19]. In Figs. 2(c) and 2(d), we choose the peak intensity *I*_{0} = 40 V^{2}/μm^{2}. It shows that more supermodes emerge compared to the case of lower intensity in Figs. 2(a) and 2(b). It should be mentioned that the intensity of *I*_{0} = 40 V^{2}/μm^{2} is far lower than the intensity damage threshold of graphene [34,35]. As denoted with red crosses in Figs. 2(c) and 2(d), there also exists supermodes which are asymmetric. The index of the asymmetric supermodes is labeled by *s′* in the order of the emergence as the intensity increases. To see clearly the asymmetric supermodes, the propagation constants are separated from all modes and shown in Figs. 2(e) and 2(f). Generally, in the nonlinear case, the number of asymmetric modes becomes larger as the number of graphene layers increases. The propagation constant increases as the layer number increases. At the same time, the propagation loss decreases, following the same property of the symmetric and anti-symmetric supermodes in the linear case.

The transverse electric field (*E _{x}*) distributions of the supermodes in the nonlinear graphene multilayer are depicted in Fig. 3 as

*N*= 3 and 4, respectively. The peak intensity of the supermode is given by

*I*

_{0}= 40 V

^{2}/μm

^{2}. There are four and six supermodes in the two graphene multilayers, respectively. As

*N*= 3, the first three supermodes (

*s*= 1, 2, 3) have explicit symmetric profiles with respect to the geometric center of the graphene multilayer. The last supermode (

*s′*= 1) is asymmetric and is similar to the mode of

*s*= 3, except for symmetry breaking. The same features can be found for

*N*= 4. The last two asymmetric modes (

*s′*= 1, 2) are similar to the mode of

*s*= 4. The difference lies in that the maximum intensity locates in different layers of graphene.

The reason for the emergence of asymmetric supermodes can be explained as follows. Although the graphene multilayer structure is symmetric, the surface conductivities of individual graphene sheets are different when taking into account the nonlinearity of graphene. Considering that the light is concentrated in the region with higher permittivity, the SPPs in the multilayer structure experience highest intensity at the graphene with largest surface conductivity. Note that larger surface conductivity suggests larger equivalent permittivity of graphene [29]. As long as the incident field is asymmetric, the asymmetric supermodes can thus be excited. The number of asymmetric supermodes can be determined with a phenomenological analysis. For a specific asymmetric mode, the peak field intensity appears in a specific layer of graphene. Note that the geometric structure is symmetric, there exists *N*/2 asymmetric supermodes as *N* is even and the value becomes (*N* − 1)/2 as *N* is odd. For example, as *N* = 3 shown in Fig. 3(a), the peak intensity can locate at any graphene sheets (*n* = 1, 2, 3). However, the asymmetric modes with the peak intensity at *n* = 1 and *n* = 3 refer to the same mode. So there is only one asymmetric supermode. As *N* = 4 shown in Fig. 3(b), the asymmetric modes with the peak intensity at *n* = 1 and *n* = 4 (or *n* = 2 and *n* = 3) refer to the same mode. Thus there exists two asymmetric supermodes. So the multilayer structure with *N* layers of graphene can support in total 3*N*/2 nonlinear supermodes as *N* is even and (3*N*−1)/2 supermodes as *N* is odd.

## 3. General properties of supermodes

We now focus on the propagation of the nonlinear supermodes in the graphene multilayer. As a nonlinear eigenmode is injected into the graphene multilayer, the evaluation of the mode profile in the waveguide can be demonstrated by taking advantage of the modified split-step Fourier beam propagation method [36]. There are six supermodes in the graphene multilayer as *N* = 4. All modes can propagate steadily in the structure with the mode profile remains fairly constant except for the propagation loss, as shown in Fig. 4. The linear-like supermodes are depicted in Figs. 4(a)-4(d). Most of power is carried in two symmetric layers of graphene. In contrast, the power is located at a single layer of graphene for the asymmetric supermodes, as illustrated in Figs. 4(e) and 4(f). Exciting the asymmetric supermodes may benefit to realizing waveguide routing and all-optical switches.

The propagation constant as a function of the intensity is plotted in Fig. 5. As *N* = 3, there are three supermodes as the intensity is weak. For all supermodes, the propagation constant increases as the mode intensity increases. As the intensity approaches to *I*_{0} ≈34 V^{2}/μm^{2}, there appears the asymmetric supermode of *s′* = 1, as shown in Fig. 5(a). It suggests that the excitation of the asymmetric supermode requires a threshold of intensity. The figure also shows clearly that the asymmetric supermode branches out from the highest linear-like supermode (*s* = 3). The similar process exactly happens for *N* = 4, as shown in Fig. 5(b). There are two asymmetric supermodes (*s′* = 1, 2) except for the four linear-like modes and they both branch from the mode *s* = 4. The thresholds of intensity for the two supermodes are about *I*_{0} = 16 and 36 V^{2}/μm^{2}, respectively. In addition, as long as the intensity reaches the required excitation threshold, the number of supermodes will not changes and no new asymmetric modes emerge as the intensity further increases. It should be mentioned that the intensity considered here is far lower than the damage threshold of graphene [34,35].

Figures 6(a) and 6(b) illustrate the influence of incident wavelength and chemical potential of graphene on the intensity threshold of the asymmetric supermodes for *N* = 4. Generally, the intensity thresholds of the two asymmetric supermodes (*s′* = 1, 2) experience the same dependence on the chemical potential and wavelength. The higher asymmetric supermode (*s′* = 2) shown in Fig. 5(b) has larger threshold and is harder to excite. As the chemical potential increases, the threshold increases fast. For instance, the threshold of intensity reads *I*_{0} = 0.03 V^{2}/μm^{2} at μ_{c} = 0.1 eV and λ = 10 μm for the asymmetric supermode of *s′* = 2. It increases to *I*_{0} = 325 V^{2}/μm^{2} as the chemical potential is doubled, as illustrated in Fig. 6(b). On the other hand, when the chemical potential is fixed, the threshold also increases quickly as the incident wavelength increases. For example, the threshold reaches *I*_{0} = 873 V^{2}/μm^{2} at λ = 16 μm and μ* _{c}* = 0.2 eV for the asymmetric supermode of

*s′*= 2, which is more than twice the threshold at λ = 10 μm. Consequently, the threshold of intensity could also be sensitively tuned by varying the incident wavelength and chemical potential of graphene.

Now we discuss the experimental proposals of exciting the nonlinear supermodes. The key element here is to fulfill the mode matching condition between the excitation field and nonlinear mode profiles. The graphene SPPs can be excited by nanoemitters such as quantum dots, molecules lying on graphene illuminated by the pumper laser [37,38]. The required mode profiles for the nonlinear supermodes can be obtained by individually controlling the amplitudes and relative phases of SPPs in each graphene by adjusting the size and orientation of nanoemitters. The required excitation thresholds of the asymmetric supermodes are realized by tuning the pump power or the laser spot size. In practice, the laser pump power is controlled by applying tunable optical attenuators, and the average laser spot size can be tuned by employing objectives with different numerical apertures.

## 4. Conclusion

In conclusion, we have investigated the intensity-dependent propagation constants and mode profiles of the supermodes in the nonlinear graphene multilayers. The influence of the number of graphene layers is especially concentrated. There are three types of nonlinear supermodes which are symmetric, anti-symmetric, and asymmetric. Apart from the *N* symmetric and anti-symmetric supermodes in the graphene multilayer with *N* graphene layers, the numbers of the asymmetric supermodes are *S* = *N*/2 and (*N*−1)/2 as *N* is even and odd, respectively. The intensity threshold of asymmetric modes increases sensitively as the incident wavelength or chemical potential of graphene increase. The study may find applications in developing all-optical converters and switches on the deep-subwavelength scale.

## Funding

973 Program (No. 2014CB921301); National Natural Science Foundation of China (NSFC) (Nos. 11304108 and 11674117); Natural Science Foundation of Hubei Province (2015CFA040); Specialized Research Fund for the Doctoral Program of Higher Education of China (No. 20130142120091).

## References and links

**1. **F. Bonaccorso, Z. Sun, T. Hasan, and A. C. Ferrari, “Graphene photonics and optoelectronics,” Nat. Photonics **4**(9), 611–622 (2010). [CrossRef]

**2. **R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, “Fine structure constant defines visual transparency of graphene,” Science **320**(5881), 1308 (2008). [CrossRef] [PubMed]

**3. **S. Klaiman, U. Günther, and N. Moiseyev, “Visualization of branch points in PT-symmetric waveguides,” Phys. Rev. Lett. **101**(8), 080402 (2008). [CrossRef] [PubMed]

**4. **A. Vakil and N. Engheta, “Transformation optics using graphene,” Science **332**(6035), 1291–1294 (2011). [CrossRef] [PubMed]

**5. **F. H. L. Koppens, D. E. Chang, and F. J. García de Abajo, “Graphene plasmonics: a platform for strong light-matter interactions,” Nano Lett. **11**(8), 3370–3377 (2011). [CrossRef] [PubMed]

**6. **C. Qin, B. Wang, H. Long, K. Wang, and P. Lu, “Nonreciprocal phase shift and mode modulation in dynamic graphene waveguides,” J. Lightwave Technol. **34**(16), 3877–3883 (2016).

**7. **F. Wang, C. Qin, B. Wang, H. Long, K. Wang, and P. Lu, “Rabi oscillations of plasmonic supermodes in graphene multilayer arrays,” IEEE J. Quantum Electron. **23**(1), 4600105 (2017).

**8. **H. Da and C. W. Qiu, “Graphene-based photonic crystal to steer giant Faraday rotation,” Appl. Phys. Lett. **100**(24), 241106 (2012). [CrossRef]

**9. **H. Da, Q. Bao, R. Sanaei, J. Teng, K. P. Loh, F. J. Garcia-Vidal, and C. W. Qiu, “Monolayer graphene photonic metastructures: Giant Faraday rotation and nearly perfect transmission,” Phys. Rev. B **88**(20), 205405 (2013). [CrossRef]

**10. **S. A. Mikhailov, “Non-linear electromagnetic response of graphene,” Europhys. Lett. **79**(2), 27002 (2007). [CrossRef]

**11. **K. L. Ishikawa, “Nonlinear optical response of graphene in time domain,” Phys. Rev. B **82**(20), 201402 (2010). [CrossRef]

**12. **E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, and S. A. Mikhailov, “Coherent nonlinear optical response of graphene,” Phys. Rev. Lett. **105**(9), 097401 (2010). [CrossRef] [PubMed]

**13. **A. Auditore, C. De Angelis, A. Locatelli, S. Boscolo, M. Midrio, M. Romagnoli, A.-D. Capobianco, and G. Nalesso, “Graphene sustained nonlinear modes in dielectric waveguides,” Opt. Lett. **38**(5), 631–633 (2013). [CrossRef] [PubMed]

**14. **A. V. Gorbach, “Nonlinear graphene plasmonics: amplitude equation for surface plasmons,” Phys. Rev. A **87**(1), 013830 (2013). [CrossRef]

**15. **J. Tao, Z. Dong, J. K. W. Yang, and Q. J. Wang, “Plasmon excitation on flat graphene by s-polarized beams using four-wave mixing,” Opt. Express **23**(6), 7809–7819 (2015). [CrossRef] [PubMed]

**16. **Z. Wang, B. Wang, H. Long, K. Wang, and P. Lu, “Plasmonic lattice solitons in nonlinear graphene sheet arrays,” Opt. Express **23**(25), 32679–32689 (2015). [CrossRef] [PubMed]

**17. **D. A. Smirnova, A. V. Gorbach, I. V. Iorsh, I. V. Shadrivov, and Y. S. Kivshar, “Nonlinear switching with a graphene coupler,” Phys. Rev. B **88**(4), 045443 (2013). [CrossRef]

**18. **E. Kapon, J. Katz, and A. Yariv, “Supermode analysis of phase-locked arrays of semiconductor lasers,” Opt. Lett. **9**(4), 125–127 (1984). [CrossRef] [PubMed]

**19. **C. Qin, B. Wang, H. Huang, H. Long, K. Wang, and P. Lu, “Low-loss plasmonic supermodes in graphene multilayers,” Opt. Express **22**(21), 25324–25332 (2014). [CrossRef] [PubMed]

**20. **Q. Bao and K. P. Loh, “Graphene photonics, plasmonics, and broadband optoelectronic devices,” ACS Nano **6**(5), 3677–3694 (2012). [CrossRef] [PubMed]

**21. **N. M. R. Peres, “Colloquium: the transport properties of graphene: an introduction,” Rev. Mod. Phys. **82**(3), 2673–2700 (2010). [CrossRef]

**22. **S. Ke, B. Wang, C. Qin, H. Long, K. Wang, and P. Lu, “Exceptional points and asymmetric mode switching in plasmonic waveguides,” J. Lightwave Technol. **34**(22), 5258–5262 (2016). [CrossRef]

**23. **H. Huang, S. Ke, B. Wang, H. Long, K. Wang, and P. Lu, “Numerical study on plasmonic absorption enhancement by a rippled graphene sheet,” J. Lightwave Technol. **99**(1), 1 (2017). doi:. [CrossRef]

**24. **S. A. Mikhailov and K. Ziegler, “Nonlinear electromagnetic response of graphene: frequency multiplication and the self-consistent-field effects,” J. Phys. Condens. Matter **20**(38), 384204 (2008). [CrossRef] [PubMed]

**25. **D. A. Smirnova, I. V. Shadrivov, A. I. Smirnov, and Y. S. Kivshar, “Dissipative plasmon-solitons in multilayer graphene,” Laser Photonics Rev. **8**(2), 291–296 (2014). [CrossRef]

**26. **S. A. Mikhailov and K. Ziegler, “New electromagnetic mode in graphene,” Phys. Rev. Lett. **99**(1), 016803 (2007). [CrossRef] [PubMed]

**27. **B. Wang, X. Zhang, F. J. García-Vidal, X. Yuan, and J. Teng, “Strong coupling of surface plasmon polaritons in monolayer graphene sheet arrays,” Phys. Rev. Lett. **109**(7), 073901 (2012). [CrossRef] [PubMed]

**28. **Y. Liu, G. Bartal, D. A. Genov, and X. Zhang, “Subwavelength discrete solitons in nonlinear metamaterials,” Phys. Rev. Lett. **99**(15), 153901 (2007). [CrossRef] [PubMed]

**29. **M. L. Nesterov, J. Bravo-Abad, A. Y. Nikitin, F. J. Garcia-Vidal, and L. Martin-Moreno, “Graphene supports the propagation of subwavelength optical solitons,” Laser Photonics Rev. **7**(2), L7–L11 (2013). [CrossRef]

**30. **M. Mitchell, M. Segev, T. H. Coskun, and D. N. Christodoulides, “Theory of self-trapped spatially incoherent light beams,” Phys. Rev. Lett. **79**(25), 4990–4993 (1997). [CrossRef]

**31. **O. Cohen, T. Schwartz, J. W. Fleischer, M. Segev, and D. N. Christodoulides, “Multiband vector lattice solitons,” Phys. Rev. Lett. **91**(11), 113901 (2003). [CrossRef] [PubMed]

**32. **Z. Wang, B. Wang, K. Wang, H. Long, and P. Lu, “Vector plasmonic lattice solitons in nonlinear graphene-pair arrays,” Opt. Lett. **41**(15), 3619–3622 (2016). [CrossRef] [PubMed]

**33. **F. Wang, C. Qin, B. Wang, S. Ke, H. Long, K. Wang, and P. Lu, “Rabi oscillations of surface plasmon polaritons in graphene-pair arrays,” Opt. Express **23**(24), 31136–31143 (2015). [CrossRef] [PubMed]

**34. **G. Xing, H. Guo, X. Zhang, T. C. Sum, and C. H. A. Huan, “The Physics of ultrafast saturable absorption in graphene,” Opt. Express **18**(5), 4564–4573 (2010). [CrossRef] [PubMed]

**35. **M. Currie, J. D. Caldwell, F. J. Bezares, J. Robinson, T. Anderson, H. Chun, and M. Tadjer, “Quantifying pulsed laser induced damage to graphene,” Appl. Phys. Lett. **99**(21), 211909 (2011). [CrossRef]

**36. **M. D. Feit and J. A. Fleck Jr., “Light propagation in graded-index optical fibers,” Appl. Opt. **17**(24), 3990–3998 (1978). [CrossRef] [PubMed]

**37. **A. Yu. Nikitin, F. Guinea, F. J. Garcia-Vidal, and L. Martin-Moreno, “Fields radiated by a nanoemitter in a graphene sheet,” Phys. Rev. B **84**(19), 195446 (2011). [CrossRef]

**38. **M. Fujiwara, K. Toubaru, T. Noda, H. Q. Zhao, and S. Takeuchi, “Highly efficient coupling of photons from nanoemitters into single-mode optical fibers,” Nano Lett. **11**(10), 4362–4365 (2011). [CrossRef] [PubMed]