## Abstract

We investigate the surface vector plasmonic lattice solitons (PLSs) in semi-infinite graphene-pair arrays (GPAs). The surface vector PLSs are composed of two components which are associated with different band gaps. Both components undergo mutual self-trapping at the boundary of the semi-infinite structure when the self-focusing nonlinearity of graphene and the light diffraction reach a balance. Thanks to the strong confinement of SPPs, the surface vector PLSs can be squeezed into a deep-subwavelength width of ~0.003*λ*. By comparing with bulk solitons, the surface PLSs are more readily to excite by external waves and more sensitive to the surrounding environment. The study may develop promising applications in all-optical switching devices and optical sensors on deep-subwavelength scale.

© 2017 Optical Society of America

## 1. Introduction

Graphene, a single monolayer of carbon atoms arranged in 2D honeycomb lattice, has attracted much attention with its unique features in electronics and optics [1–5]. The optical response of graphene is related to the surface conductivity and sensitive to the chemical potential. At relatively low frequencies, graphene behaves like metal with the intraband transition of electrons dominating [6]. In comparison with that on metal, surface plasmon polaritons (SPPs) in graphene exhibit superior features such as stronger field confinement, lower loss, and flexible tunability by external fields or gate voltage [6–9]. These features make graphene a promising material for novel optical nanodevices [10–12].

Recently, the nonlinear properties of graphene have also attracted extensive attention such as the third-harmonic generation, four-wave mixing, and saturable absorption [13–19]. In particular, the third-order nonlinear electrodynamic effects of graphene have been comprehensively studied by Mikhailov *et al*. [16]. The Kerr-type nonlinearity is one typical class of the third-order nonlinear effects which presents in the dependence of graphene surface conductivity on the light intensity [20–23]. Compared with a graphene monolayer, the graphene multilayer structures can exhibit nonlinear effects more significantly and diversely. For example, the graphene sheet arrays (GSAs) could support plasmonic lattice solitons (PLSs). They are spatial solitons and can remain stable distribution of transverse field during propagation.

Graphene sheet arrays are usually composed of periodic-arranged graphene sheets embedding in dielectric. The structure can provide a fertile platform for manipulating light propagation [24–27]. The PLSs supported by GSAs are formed as the graphene nonlinearity and the SPP tunneling between adjacent graphene sheets reach a balance. Plasmonic lattice solitons in the infinite GSAs have been investigated, including scalar PLSs in monolayer graphene sheet arrays [23,28] and vector PLSs in graphene-pair arrays [29]. Following the studies, surface scalar PLSs localized at and near the boundary of semi-infinite monolayer graphene sheets arrays have also been deliberated [23,30]. The nonlinear surface state is associated with the properties of the corresponding Bloch waves in the infinite counterpart. The surface vector solitons are composed of two or more mutually trapped scalar components, which have been investigated in periodic dielectric structure [31–33]. However, the plasmonic analogy has not been studied yet. According to the previous works [23,28–30], we can predict that the semi-infinite graphene-pair arrays could support surface vector PLSs with the components associated with different band gaps. The surface vector PLSs will combine advantages from both vector solitons and surface solitons, including light-control-light property, relative strong localization, and the convenience of excitation [29,30].

In this work, we shall investigate the surface vector PLSs in semi-infinite graphene-pair arrays. There are two propagation bands of SPPs in the structure. The surface vector PLSs are composed of two components belonging to different band gaps with different frequencies. When both light components excite the self-focusing nonlinearity of graphene at the interface to balance their respective diffraction into the arrays, the surface vector PLS is formed with both components propagating stably. In the structure, the surface vecor PLSs can exhibit all-light control on deep-subwavelength scale (~0.003*λ*) due to the strong confinement of SPPs on graphene and the boundary effect. The threshold powers for two components of the surface vector PLSs are also studied. By changing the chemical potential of graphene, the threshold power can be modulated flexibly.

## 2. Transverse distribution and propagation of surface vector PLSs

The structure of the semi-infinite graphene-pair arrays is shown in Fig. 1(a), where the graphene pairs (*x* > 0) embedding in the dielectric (*ε _{d}* = 2.13). The distance between two sheets of every graphene pair is

*d*

_{1}= 30 nm. The period of structure is

*d*=

*d*

_{1}+

*d*

_{2}where

*d*is the spacing between adjacent graphene pairs (

_{2}*d*

_{2}= 50 nm). Considering the transverse magnetic polarization (TM), the Maxwell’s equations can be rewritten in the matrix form [34]

*ε*stands for the relative permittivity along the

_{r}(x)*x*axis. Equation (1) is aimed at solving the eigen problem. The transverse fields,

*H*and

_{y}*E*, are two elements of the eigenvector, while the propagation constant

_{x}*k*is the eigenvalue. In the linear case, the solutions of Eq. (1) form two bands of eigenvalues, which represent the linear propagation constants of SPPs in the structure. The relation between propagation constant

_{z}*β*and Bloch momentum

*φ*=

*k*with

_{x}d*k*the transverse Bloch wave vector, namely the dispersion relation of SPPs in GPAs, is shown in Fig. 1. There are two propagation bands for the structure in linear case as the wavelength is either

_{x}*λ*

_{1}= 9.8 μm or

*λ*

_{2}= 10 μm. For

*λ*

_{1}= 9.8 μm, the first and second bands occupy the range of propagation constants

*β*

_{1}= 95.4 ~96.6 μm

^{−1}and 82.1 ~85.4 μm

^{−1}, respectively. Accordingly, the finite gap is from 85.4 μm

^{−1}to 95.4 μm

^{−1}. For

*λ*= 10 μm, the first and second bands have the propagation constant rang of

_{2}*β*

_{2}= 91.3 ~92.8 μm

^{−1}and

*β*

_{2}= 76.5 ~80.7 μm

^{−1}, respectively. The finite gap is from 80.7 μm

^{−1}to 91.3 μm

^{−1}. Two components of the surface vector PLSs respectively originate from the Brillouin zone center

*φ*= 0 for

*λ*

_{1}= 9.8 μm and the Brillouin zone edge

*φ*=

*π*for

*λ*

_{2}= 10 μm, as presented with solid dots in Figs. 1(b) and 1(c). The first component stays in the semi-infinite gap (

*β*

_{1}> 96.6 μm

^{−1}), while the second component corresponds to the finite band gap (80.7 μm

^{−1}<

*β*

_{2}< 91.3 μm

^{−1}). Both components exhibit normal diffraction which will be inhibited by the self-focusing nonlinearity of graphene.

In the numerical calculation, graphene is treated as a thin film with an equivalent thickness Δ ≈1 nm. The relative equivalent permittivity of graphene could be given by *ε _{g}* = 1 + i

*σ*/(

_{g}η_{0}*k*Δ) [14], where

_{0}*η*and

_{0}*k*are the impendence and propagation constant in vacuum. The nonlinearity of graphene originates from its surface conductivity

_{0}*σ*. It is proportional to the tangential electric field which can be written as

_{g}*σ*=

_{g}*σ*+

_{g,linear}*σ*|

^{NL}*E*|

_{z}^{2}. The nonlinear conductivity

*σ*is obtained from the general quantum theory of the graphene third-order nonlinear conductivity [16],

^{NL}*σ*

_{0}^{(3)}=

*e*

^{4}

*ћV*

_{F}^{2}/(4

*πμ*

_{c}^{4}) with

*μ*is the chemical potential of graphene and the Fermi velocity

_{c}*V*≈

_{F}*c*/300.

*S*

^{(}

^{m}^{/}

^{n}^{)}is dimensionless and represents the combination of intraband and interband contributions. The terms

*n*= 0 and

*m*= 0 refer to the net intraband contribution and net interband contribution, respectively. The other two terms are due to the combination of both intraband and interband contributions. The details are presented in the reference [16]. The linear part of the surface conductivity

*σ*(

_{g,linear}*λ*,

*μ*,

_{c}*τ*) is also obtained from the model [16], which is consistent with the result of Kubo formula [35,36], where

*λ*is the incident wavelength in air and

*τ*is the momentum relaxation time. Here, these parameters are respectively set as

*λ*

_{1}= 9.8 μm,

*λ*

_{2}= 10 μm,

*μ*= 0.16 eV,

_{c}*τ*= 0.5 ps [37–40]. In the condition, the graphene permittivity increases as the light intensity increases, indicating that the graphene nonlinearity belongs to the self-focusing nonlinear type.

Taking the graphene nonlinearity into consideration, Eq. (1) turns into a nonlinear eigen problem, which could be solved by using the self-consistent method [41]. We set two Gaussian field distributions around the array boundary as initial values for the numerical iterative calculation, where the peak intensities are set as 180 V^{2}/μm^{2} and 220 V^{2}/μm^{2} for the first and second components respectively. In the calculation, the transverse magnetic field *H _{y,m}* and transverse electric field

*E*are obtained by solving the nonlinear eigen problem of Eq. (1), at meanwhile, the tangential electric field

_{x,m}*E*is derived by the relation which is also from Maxwell’s Equations

_{z,m}*m*= 1, 2 corresponding the first and second field components. The relative permittivity

*ε*in the zone occupied by graphene is equal to the value of

_{r}(x)*ε*= 1 + i

_{g}*σ*/(

_{g}η_{0}*k*Δ). While in the dielectric zone, the relative permittivity

_{0}*ε*equals to the constant value of

_{r}(x)*ε*= 2.13. As the relative permittivity of graphene depends on the field intensity, the relative permittivity

_{d}*ε*corresponding to graphene location is updated with the solved field distributions (|

_{r}(x)*E*|

_{z}^{2}= |

*E*

_{z,}_{1}|

^{2}+ |

*E*

_{z,}_{2}|

^{2}) of the previous step during every iterative process. Repeating the process of solving the field distributions of two components and updating the permittivity, we can obtain the transverse distribution of the surface vector PLSs until the eigenvalues are stable.

The transverse field distributions of the two components are obtained as shown in Fig. 2. Here, the loss of graphene is not taken into account at first. It has also been verified that the distributions do not change much even though the loss is included. For the first component, the tangential electric field (*E _{z}*) and the intensity distribution |

*E*|

^{2}= |

*E*|

_{x}^{2}+ |

*E*|

_{z}^{2}are presented in Figs. 2(a) and 2(b). Its propagation constant is obtained as

*β*

_{1}= 98.2 μm

^{−1}which stays in the semi-infinite band gap (

*β*

_{1}> 96.6 μm

^{−1}). The power and the effective width of surface PLS are respectively given by [42,43]

*w*

_{1}= 0.027 μm, which is equivalent to 0.0027

*λ*, while the input power

*P*

_{1}= 39.8 W/m. For the second component, the propagation constant

*β*

_{2}is about 82.4 μm

^{−1}, belonging to the finite gap (80.7 μm

^{−1}<

*β*

_{2}< 91.3 μm

^{−1}). Its normalized tangential electric field (

*E*) and intensity distribution (|

_{z}*E*|

^{2}) are shown in Figs. 2(c) and 2(d). The effective width

*w*

_{2}= 0.04 μm, which is equivalent to 0.004

*λ*. The power of the second component

*P*

_{2}= 48.7 W/m. However, in the infinite plasmonic lattice, the effective widths of two components are

*w*

_{1}= 0.03 μm and

*w*

_{2}= 0.046 μm under the same input powers. In Fig. 2, the width of the bulk PLS represented by the red dash line is larger than that based on the semi-infinite lattice. It indicates that surface vector PLSs are easier to excite and exhibit stronger power confinement in comparison with the bulk counterparts in the infinite GPAs. Besides, as the surface vector PLSs distribute at the interface of graphene arrays and an open half-space, they can interact with other materials and be modulated by the external environment more readily. As a result, the surface vector PLSs may find applications in optical sensors and environment detection.

When we modulate the intensity of arbitrary component, the propagation constants of both components are changed simultaneously, as shown in Figs. 3(a) and 3(b). As the light intensity is enhanced, the values of both propagation constants increase, and thus stay deeper in their corresponding band gaps. We also consider the influence of light intensity on the widths and propagation distances of two components. As the intensity of the first component increases, both the widths of the two components decrease, as shown in Fig. 3(c). As the intensity increases, the nonlinear effect is strengthened, leading to more localized soliton distribution. When the loss of graphene is considered, the soliton propagation length is obtained by *L*_{1,2} = 1/(2Im(*β*_{1,2})). As shown Fig. 3(d), the real propagation lengths of two components increase as the intensity of the first component increases.

To demonstrate the propagation of the surface vector PLSs, we substitute the nonlinear eigenmode solutions (see Fig. 2) into the full Maxwell’s equations. The soliton propagation is simulated by using the modified split-step Fourier beam propagation method [34]. In the lossless case, two components propagate stably as shown in Figs. 4(a) and 4(b). As any one of the two components is removed, the other one will diffuse into GPAs, as shown in Figs. 4(c) and 4(d). It means that two components of the surface vector PLSs should exist simultaneously to remain the stable propagation. Their mutual influence during the light propagation will benefit to developing all-optical switches. We also present the propagation of the surface vector PLS in the real lossy case, as shown in Figs. 4(e) and 4(f). In the nonlinear condition, the propagation constants of two components are obtained as *β*_{1} = 98.2 + 1.1*i* μm^{−1} and *β*_{2} = 82.4 + 1.3*i* μm^{−1}. Accordingly, the propagation lengths of both components are about 0.4 μm in Figs. 4(e) and 4(f).

## 3. Threshold power of surface vector PLSs

We note that the first component of surface vector PLS originate from the Brillouin zone center *φ* = 0 with *λ*_{1} = 9.8 μm and stays in the semi-infinite gap (*β*_{1} > 96.6 μm^{−1}), while the second component (*λ _{2}* = 10 μm) corresponds to the Brillouin zone edge

*φ*=

*π*and emerges from the finite gap (80.7 μm

^{−1}<

*β*

_{2}< 91.3 μm

^{−1}).

When only considering the first component, the relation between power and propagation constant of the formed scalar soliton is shown in Fig. 5(a). From the relation curve, the threshold power is *P _{c}* = 46.5 W/m which is the minimum power to excite scalar surface soliton of the first component. From Fig. 5(a), all the nonlinear surface localized modes staying in the semi-infinite gap are separated from the first propagation band. The localized field of the scalar soliton redistributes the structure permittivity due to the intrinsic nonlinear effect. The structure with updated permittivity distribution also guides eigen modes for

*λ*

_{2}= 10 μm. In the condition only injecting the light of the first component, the eigen mode for

*λ*

_{2}= 10 μm and locating at the Brillouin zone edge

*φ*=

*π*is regarded as a linear guided mode [32], which corresponds to the second component. Accordingly, the guided mode is modulated by the scalar solitons of the first component. Figure 5(b) shows the propagation constant of the linear guided mode, which corresponds to the second component, varies with that of the scalar solitons of the first component. Their propagation constants staying in the gap are also separated from the second propagation band. Focusing on the scalar soliton of the second component, the power varies with the propagation constant as shown in Fig. 5(c). The threshold power is about 73.5 W/m. Similarly, the guided modes for the first component are also influenced by the scalar light field of the second component. The relation between the propagation constant of linear guided modes for the first component and that of the input scalar field for the second component is shown in Fig. 5(d). It should be noted that the composite state of the scalar soliton and the corresponding guided mode could also be regarded as surface vector solitons [32].

To further investigate the interaction between the two components of the surface vector PLS, we plot the transverse distribution of the scalar soliton for the first component in Fig. 6(a). The nonlinear mode corresponds to the dot a_{1} in Fig. 5(a), which has the peak intensity *I*_{1} = 180 V^{2}/μm^{2}. The transverse profile of the mode guided by the scalar soliton is shown in Fig. 6(b), which corresponds to the dot a_{2} in Fig. 5(b). As the intensity amplitude of the guided mode increases, it acts with the surface scalar soliton, leading to the formation of mutual-trapping stable light beams, namely the surface vector solitons. As the peak intensity increases to *I*_{2} = 220 V^{2}/μm^{2}, the surface vector PLS is formed, as shown in Fig. 2. The effective width of the scalar soliton presented in Fig. 6(a) is about 0.05 μm, which is almost twice larger than that of 0.027 μm in Fig. 2(b), although both of them have the same peak intensity *I*_{1} = 180 V^{2}/μm^{2}. The added second component can enhance the localization of the first component. For the surface scalar soliton shown in Fig. 6(a), the power *P*_{1} = 47.4 W/m, while the power of the vector soliton component in Fig. 2(b) is *P*_{2} = 39.8 W/m, which is even smaller than the threshold power of the single first component (*P _{c}* = 46.5 W/m). We can see that the added light beam can also decrease the threshold power for the first light component and make the surface PLS inspired more readily. We also consider the influence of the chemical potential on the threshold power for each single component. As shown in Fig. 7, the threshold power for arbitrary component increases as the chemical potential increases. The enhancement of the chemical potential leads to stronger curvature of the dispersion-relation curve [28,29], which means that the diffraction becomes stronger and needs more nonlinear effect to balance.

## 4. Conclusions

In conclusion, we have investigated the surface vector PLSs in semi-infinite graphene-pair arrays. There are two propagation bands of SPPs in the structure. The components of surface vector PLSs have different frequencies and correspond to distinct band gaps. They undergo mutual self-trapping through the balance between SPPs diffraction and the nonlinearity effect of graphene. Due to the strong confinement of SPPs in graphene, the width of surface vector PLSs can reach deep-subwavelength scale. By varying the chemical potential of graphene, the threshold power can be tuned flexibly. The study may find applications in all-optical switches, optical sensors, and soliton-based navigation on deep-subwavelength scale.

## Funding

Program 973 (No. 2014CB921301); National Natural Science Foundation of China (NSFC) (Nos. 11674117 and 11304108); 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. **K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, “Electric field effect in atomically thin carbon films,” Science **306**(5696), 666–669 (2004). [CrossRef] [PubMed]

**2. **Y. Zhang, Y. W. Tan, H. L. Stormer, and P. Kim, “Experimental observation of the quantum Hall effect and Berry’s phase in graphene,” Nature **438**(7065), 201–204 (2005). [CrossRef] [PubMed]

**3. **A. K. Geim and K. S. Novoselov, “The rise of graphene,” Nat. Mater. **6**(3), 183–191 (2007). [CrossRef] [PubMed]

**4. **A. K. Geim, “Graphene: status and prospects,” Science **324**(5934), 1530–1534 (2009). [CrossRef] [PubMed]

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

**6. **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]

**7. **J. Christensen, A. Manjavacas, S. Thongrattanasiri, F. H. L. Koppens, and F. J. de Abajo, “Graphene Plasmon Waveguiding and Hybridization in Individual and Paired Nanoribbons,” ACS Nano **6**(1), 431–440 (2012). [CrossRef] [PubMed]

**8. **Y. Fan, B. Wang, H. Huang, K. Wang, H. Long, and P. Lu, “Plasmonic Zitterbewegung in binary graphene sheet arrays,” Opt. Lett. **40**(13), 2945–2948 (2015). [CrossRef] [PubMed]

**9. **J. S. Gomez-Diaz, M. Tymchenko, and A. Alù, “Hyperbolic Plasmons and Topological Transitions Over Uniaxial Metasurfaces,” Phys. Rev. Lett. **114**(23), 233901 (2015). [CrossRef] [PubMed]

**10. **H. Yan, X. Li, B. Chandra, G. Tulevski, Y. Wu, M. Freitag, W. Zhu, P. Avouris, and F. Xia, “Tunable infrared plasmonic devices using graphene/insulator stacks,” Nat. Nanotechnol. **7**(5), 330–334 (2012). [CrossRef] [PubMed]

**11. **A. Andryieuski, A. V. Lavrinenko, and D. N. Chigrin, “Graphene hyperlens for terahertz radiation,” Phys. Rev. B **86**(12), 121108 (2012). [CrossRef]

**12. **Z. Li, K. Yao, F. Xia, S. Shen, J. Tian, and Y. Liu, “Graphene Plasmonic Metasurfaces to Steer Infrared Light,” Sci. Rep. **5**(1), 12423 (2015). [CrossRef] [PubMed]

**13. **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]

**14. **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]

**15. **S. Hong, J. I. Dadap, N. Petrone, P. Yeh, J. Hone, and R. M. Osgood Jr., “Optical third-harmonic generation in graphene,” Phys. Rev. X **3**(2), 021014 (2013). [CrossRef]

**16. **S. A. Mikhailov, “Quantum theory of the third-order nonlinear electrodynamic effects of graphene,” Phys. Rev. B **93**(8), 085403 (2016). [CrossRef]

**17. **S. A. Mikhailov, “Nonperturbative quasiclassical theory of the nonlinear electrodynamic response of graphene,” Phys. Rev. B **95**(8), 085432 (2017). [CrossRef]

**18. **A. Marini, J. D. Cox, and F. J. Garcia de Abajo, “Theory of graphene saturable absorption,” Phys. Rev. B **95**(12), 125408 (2017). [CrossRef]

**19. **P. Lan, M. Ruhmann, L. He, C. Zhai, F. Wang, X. Zhu, Q. Zhang, Y. Zhou, M. Li, M. Lein, and P. Lu, “Attosecond Probing of Nuclear Dynamics with Trajectory-Resolved High-Harmonic Spectroscopy,” Phys. Rev. Lett. **119**(3), 033201 (2017). [CrossRef] [PubMed]

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

**21. **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]

**22. **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]

**23. **Y. V. Bludov, D. A. Smirnova, Y. S. Kivshar, N. M. R. Peres, and M. I. Vasilevskiy, “Discrete solitons in graphene metamaterials,” Phys. Rev. B **91**(4), 045424 (2015). [CrossRef]

**24. **H. Huang, B. Wang, H. Long, K. Wang, and P. Lu, “Plasmon-negative refraction at the heterointerface of graphene sheet arrays,” Opt. Lett. **39**(20), 5957–5960 (2014). [CrossRef] [PubMed]

**25. **Y. Fan, B. Wang, H. Huang, K. Wang, H. Long, and P. Lu, “Plasmonic Bloch oscillations in monolayer graphene sheet arrays,” Opt. Lett. **39**(24), 6827–6830 (2014). [CrossRef] [PubMed]

**26. **C. Qin, B. Wang, H. Long, K. Wang, and P. Lu, “Bloch mode engineering in graphene modulated periodic waveguides and cavities,” J. Opt. Soc. Am. B **32**(8), 1748–1753 (2015). [CrossRef]

**27. **Y. Fan, B. Wang, K. Wang, H. Long, and P. Lu, “Plasmonic Zener tunneling in binary graphene sheet arrays,” Opt. Lett. **41**(13), 2978–2981 (2016). [CrossRef] [PubMed]

**28. **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]

**29. **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]

**30. **Z. Wang, B. Wang, H. Long, K. Wang, and P. Lu, “Surface plasmonic lattice solitons in semi-infinite graphene sheet arrays,” J. Lightwave Technol. **35**(14), 2960–2965 (2017). [CrossRef]

**31. **J. Hudock, S. Suntsov, D. Christodoulides, and G. Stegeman, “Vector discrete nonlinear surface waves,” Opt. Express **13**(20), 7720–7725 (2005). [CrossRef] [PubMed]

**32. **I. Garanovich, A. A. Sukhorukov, Y. S. Kivshar, and M. Molina, “Surface multi-gap vector solitons,” Opt. Express **14**(11), 4780–4785 (2006). [CrossRef] [PubMed]

**33. **Y. V. Kartashov, F. Ye, and L. Torner, “Vector mixed-gap surface solitons,” Opt. Express **14**(11), 4808–4814 (2006). [CrossRef] [PubMed]

**34. **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]

**35. **P. Y. Chen and A. Alù, “Atomically thin surface cloak using graphene monolayers,” ACS Nano **5**(7), 5855–5863 (2011). [CrossRef] [PubMed]

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

**37. **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]

**38. **F. Wang, C. Qin, B. Wang, H. Long, K. Wang, and P. Lu, “Rabi oscillations of Plasmonic Supermodes in Graphene Multilayer Arrays,” IEEE J. Sel. Top. Quantum Electron. **23**(1), 4600105 (2017). [CrossRef]

**39. **C. Qin, B. Wang, H. Long, K. Wang, and P. Lu, “Non-reciprocal Phase Shift and Mode Modulation in Dynamic Graphene Waveguides,” J. Lightwave Technol. **34**(16), 3877–3883 (2016).

**40. **H. Huang, S. Ke, B. Wang, H. Long, K. Wang, and P. Lu, “Numerical study on plasmonic absorption enhancement by a rippled grapheme sheet,” J. Lightwave Technol. **35**(2), 320–324 (2017). [CrossRef]

**41. **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]

**42. **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]

**43. **Y. Kou, F. Ye, and X. Chen, “Multiband vector plasmonic lattice solitons,” Opt. Lett. **38**(8), 1271–1273 (2013). [CrossRef] [PubMed]