## Abstract

We investigate the plasmonic lattice solitons (PLSs) in nonlinear graphene sheet arrays (GSAs) composed of spatially separated graphene sheets embedded in dielectric. Both the nonlinearities of graphene and dielectric are considered. The self-focusing PLSs at the Brillouin zone edges can be yielded by balancing the normal diffraction of surface plasmon polaritons (SPPs) via either the nonlinear effect of graphene or self-focusing dielectric. The self-defocusing PLSs corresponding to anomalous diffraction of SPPs at the Brillouin zone center could be yielded by the nonlinearity of self-defocusing dielectric alone. The width and propagation distance of the PLSs are dependent on the period of the GSAs and the chemical potential of graphene. Thanks to the strong confinement of SPPs, the PLSs in GSAs can be squeezed into an effective width as small as *λ*/250. The study may find applications in optical circuits and switches on deep-subwavelength scale.

© 2015 Optical Society of America

## 1. Introduction

Optical spatial solitons refer to stationary propagation of light beam with a constant transverse profile as the diffraction and nonlinear effects arrive at a dynamic balance in nonlinear materials [1,2 ]. They have continuously attracted attention due to the significant applications in optical interconnect, frequency conversion, and soliton-based navigation [3,4 ]. In recent years, the discrete spatial solitons [5–8 ], supported by discrete optical systems such as periodic waveguide arrays, have been intensively investigated. As the spatial solitons closely associate with diffraction, the waveguide arrays allow for diffraction engineering [9] and thus provide a flexible platform for yielding and manipulating solitons [7,9 ]. So far, the discrete solitons have been proposed and demonstrated in metallic waveguide arrays and are termed as plasmonic lattice solitons (PLSs), where the nonlinearity is often provided by dielectric between metallic films [10–14 ]. Due to the concentration of surface plasmon polaritons (SPPs) on metals, the PLSs could possess a subwavelength effective width.

Since the practical preparation of graphene, it has spurred widespread interest due to its unique features in electronics and optics [15–21 ]. The SPPs on graphene have stronger confinement and relatively less propagation loss in comparison with that on metals. Following the dielectric and metallic waveguide arrays, the graphene sheet arrays (GSAs) composed of alternatively stacked graphene and dielectric layers are expected to manifest analogical features. It has been reported that a variety of linear optical effects could be realized in GSAs, such as negative refraction, Talbot effect, and Bloch oscillations [22–24 ]. Since graphene is also a nonlinear material [25–29 ], the nonlinear GSAs may provide a new platform to support many nonlinear effects that have been demonstrated in dielectric and metallic waveguide arrays. Moreover, the nonlinear effects should take on a new look due to the unique properties of graphene.

In this work, we shall study the PLSs based on the SPP Bloch modes in GSAs. In the linear GSAs, the SPPs undergo anomalous diffraction at the Brillouin zone center, where the lateral wave vector of SPPs approaches to zero. The anomalous diffraction makes the SPPs tend to be focusing in the GSAs. On the other hand, the normal diffraction at the Brillouin zone edges with large lateral wave vector of SPPs leads to defocusing propagation. When considering the nonlinearity in the GSAs, the SPP diffraction could be balanced by the nonlinear self-focusing and self-defocusing effects. In practice, the nonlinearity comes from either the graphene or the dielectrics between graphene sheets in the array. As graphene is a self-focusing material, the normal diffraction can be restrained to form the PLSs. Because the dielectrics could be self-focusing or self-defocusing material [30,31 ], they are able to balance normal and anomalous diffraction, respectively, resulting in formation of self-focusing or self-defocusing PLSs [9]. Furthermore, the PLSs experience smaller width and larger propagation distance when both the nonlinearities of graphene and self-focusing dielectric are involved. We also study the PLSs shaping by varying the structure period and the chemical potential of graphene.

The contents are arranged as follows. In Sec. 2, we give a brief review on the general properties of linear GSAs, aiming at showing the diffraction properties of SPPs in different regions of the Brillouin zone. Section 3 discusses the formation of PLSs in the GSAs when only the nonlinearity of graphene takes effect. Next, we consider the influence of the dielectric nonlinearity in Sec. 4. We firstly employ self-focusing and self-defocusing dielectrics without considering the graphene nonlinearity to realize self-focusing and self-defocusing PLSs at the edge and center of the Brillouin zone, respectively. Then we utilize both the graphene nonlinearity and self-focusing dielectric to obtain highly compact PLSs in the GSAs. The influences of the period of GSAs and the chemical potential of graphene on the PLSs are studied in Sec. 5.

## 2. SPPs in linear GSAs

We start from briefly reviewing the linear GSAs. The diagram of the GSAs is shown in Fig. 1(a)
, where the graphene sheets are periodically embedded in the host dielectric with a period of *d*. The relative permittivity of the dielectric is denoted by *ε _{d}*. The surface conductivity

*σ*(

_{g}*λ, μ*) of graphene is governed by the Kubo formula [32,33 ], where

_{c}, τ*λ*is the incident wavelength in air,

*μ*is the chemical potential of graphene, and

_{c}*τ*is the momentum relaxation time. Concerning the transverse magnetic polarization, the dispersion relation of SPPs in GSAs is given by [34,35 ]

*φ*=

*k*is the Bloch momentum [7,36 ] with

_{x}d*k*being the Bloch wave vector along the

_{x}*x*direction,

*κ*= (

*k*

_{z}^{2}−

*ε*

_{d}k_{0}

^{2})

^{1/2},

*ξ*=

*η*

_{0}

*σ*/(

_{g}*iε*),

_{d}k_{0}*k*is the wave vector of SPPs,

_{z}*k*

_{0}= 2

*π*/

*λ*with

*λ*being the wavelength in air, and

*η*

_{0}is the impedance of air. For a fixed wavelength

*λ*= 10 μm, we can plot the wave vector of SPPs as a function of the Bloch momentum in stem of the dispersion relation, namely, the diffraction relation of SPPs. The diffraction relation can be regarded as the band structure in the waveguide arrays [37], which is depicted in Fig. 1(b) with the solid curve. The parameters are given by

*d*= 30 nm,

*ε*= 2.25,

_{d}*μ*= 0.15 eV, and

_{c}*τ*= 0.5 ps. As the SPP Bloch modes extend over the entire GSA, the localized incident waves injected from several graphene sheets have to suffer diffraction during propagation. If the incident wave packet possesses a lateral wave vector located at the Brillouin zone center (

*φ*= 0), the diffraction is anomalous due to the hyperbolic shape of the diffraction curve [9]. The diffraction at the Brillouin zone edges (

*φ*= ±

*π*) is normal since the shape of the diffraction curve is parabolic, analogous to diffraction in free space. Figure 1(b) also plots the propagation loss of SPPs in the GSA (black solid curve). As

*k*=

_{x}*π/d*, Re(

*k*) is about 106.1 μm

_{z}^{−1}. The corresponding wavelength of the SPP Bloch mode is about 59.2 nm and the propagation length in respect of field intensity is 350 nm, which is about six times as large as the SPP wavelength. The figure of merit of the ratio is even larger than that of SPPs in monolayer graphene [38].

The diffraction curve can be manipulated by changing the structure period and the chemical potential of graphene. As shown in Fig. 1(c), the relation is rather flat when the period is large and the coupling between adjacent graphene sheets is weak. Once the period decreases approaching to the mode width of SPPs in a single-layer graphene, which is equal to *ξ*, the diffraction curve appears a sharp corner at the Brillouin zone center, forming the so- called Dirac point in the band structure [39,40
]. In the vicinity of the Dirac point, SPPs experience strong coupling [34]. The influence of the chemical potential on the diffraction curve is shown in Fig. 1(d). The Dirac point emerges as the chemical potential of graphene becomes larger. The crucial value of the chemical potential is about 0.2 eV as *d* = 30 nm.

We also perform first-principles simulations to achieve the eigen wave vectors of SPPs in the GSAs, which is more applicable for nonlinear problems. In the simulation, the Maxwell’s equations are transformed to the matrix form [11]

*ε*(

_{r}*x*) denotes the local relative permittivity along the lateral dimension of the structure,

*H*and

_{y}*E*are the transverse magnetic field and electric field. In the calculation, graphene is treated as a thin film with an equivalent thickness of Δ ≈1 nm. Then the relative equivalent permittivity of graphene could be given by

_{x}*ε*= 1 +

_{g}*iσ*

_{g}η_{0}/(

*k*

_{0}Δ) [26]. The equation can be solved by using the finite-difference frequency-domain method [41]. We employ the Yee’s grids to discretize the distributions of the field and relative permittivity. Thus solving Eq. (2) becomes the problem of calculating the eigen values of a large sparse matrix. For the linear case, the diffraction curve calculated by Eq. (2) is also shown in Fig. 1(b) with circles, which agrees fairly with the rigorous data. For nonlinear problems, this method is readily utilized in the self-consistent iteration algorithm [42].

## 3. PLSs yielded by graphene nonlinearity

We firstly consider the PLSs resulted from the nonlinearity of graphene. The surface conductivity of graphene can be expressed by *σ _{g}* =

*σ*+

_{g,linear}*σ*|

^{NL}*E*|

_{z}^{2}, and

*σ*is the nonlinear conductivity, which is given by [25,28 ]

^{NL}*ω*denotes the frequency and the Fermi velocity

*V*≈

_{F}*c*/300. Taking into account the nonlinear surface conductivity, the equivalent permittivity is also intensity-dependent. Thus Eq. (2) turns into a nonlinear problem, which can be solved by using the self-consistent method [42, 43 ]. As the relative permittivity of graphene depends on the field distribution, we initially solve the eigen values for the linear case according to Eq. (2). The relative permittivity is then updated in terms of the field distribution (eigen vectors). Repeating the process of solving the field and updating the permittivity, we can obtain the field distribution of PLSs until the eigen values are stable.

Normally, the nonlinear change of the graphene permittivity is positive as the SPP intensity increases, indicating the self-focusing nonlinearity of graphene. In order to realize PLSs in GSAs, the excited mode should experience normal diffraction in the GSA. In the diffraction curve shown in Fig. 1(b), the normal diffraction locates at the edge of the Brillouin zone where the curve is nearly parabolic. When the mode at *k _{x}* =

*π*/

*d*is excited, the self-focusing nonlinearity of graphene will balance the normal diffraction. The lateral electric field distribution of the PLSs is shown in Fig. 2(a) . In the calculation, the maximum value of intensity |

**E**|

^{2}is set as 200 V

^{2}/μm

^{2}, where |

**E**|

^{2}= |

*E*|

_{x}^{2}+ |

*E*|

_{z}^{2}. The other parameters are the same with that used in Fig. 1(b). The profile generally follows odd coupling of SPPs in adjacent graphene and manifests a good confinement. The corresponding intensity distribution is shown in Fig. 2(b). The full width at half maximum (FWHM) of the intensity profile approximately equals to 40 nm, which reaches a deep subwavelength of 0.004

*λ*.

As the nonlinear eigen mode of the PLS is injected into the GSA, the evaluation of the profile in the GSA can be demonstrated by using the modified split-step Fourier beam propagation method [11]. In the calculation, we consider the full nonlinear Maxwell’s equations, where the linear part of the standard beam propagation method is altered by projecting the injected wave packet on the linear modes of GSAs, and the nonlinear part is approximately the sum of the front several terms of its series expansion in matrix form (see Appendix). To illustrate clearly the propagation of the PLSs, we ignore the loss of graphene at first. The PLS propagation in a lossless GSA is shown in Fig. 3(a)
. One sees evidently that the profile remains unchanged during propagation. As the nonlinear effect is not activated, the SPPs diffract into a wide direction as shown in Fig. 3(b). In case the loss is accounted (*τ* = 0.5 ps), the diffraction-free propagation of the PLSs can still be observed except for a decay of the energy, as shown in Fig. 3(c). The propagation distance is about 0.35 μm. One can introduce gain into dielectric materials by using quantum dots [44], or implement graphene doping [45] to compensate for the loss.

## 4. PLSs yielded by the nonlinearity of both graphene and dielectric

Since graphene is a self-focusing material, only the SPP modes experiencing normal diffraction at the edges of the Brillouin zone can be utilized to yield PLSs. That demands to excite the SPP modes with a large lateral wave vector. In order to generate PLSs more conveniently, we also consider the nonlinearity of dielectrics between graphene, which may help to yield PLSs based on the SPP mode at the center of the Brillouin zone. To reveal the role of the dielectric nonlinearity, we investigate the PLSs formed by only the dielectric nonlinearity at first and then both nonlinearities together with graphene.

The relative permittivity of dielectric can be expressed as *ε _{d}* = [(

*ε*

_{d}_{,linear})

^{1/2}+

*n*

_{2}|

**E**|

^{2}]

^{2}, where

*n*

_{2}is the self-focusing or self-defocusing nonlinear coefficient, depending on whether it is positive or negative. Here we choose the nonlinear coefficient

*n*

_{2}= 4.8 × 10

^{−4}μm

^{2}/V

^{2}. The maximum value of intensity |

**E**|

^{2}is set as 100 V

^{2}/μm

^{2}, corresponding to the variation of the refractive index Δ

*n*= 0.048 [30]. The transverse electric field and intensity distributions of the PLSs at

_{max}*k*=

_{x}*π*/

*d*are shown in Figs. 4(a) and 4(b) , respectively. The FWHM of the intensity profile equals to 37 nm, as shown in Fig. 4(b). The value is very close to that of the PLSs by using graphene nonlinearity. It means that the nonlinear effects derived from self-focusing dielectrics and graphene could contribute independently to the formation of PLSs and their contributions are equal. When excited at

*k*= 0, the SPPs undergo anomalous diffraction. In order to restrain the diffraction and yield PLSs, the self-defocusing dielectrics have to be utilized. Here we choose

_{x}*n*

_{2}= −4.8 × 10

^{−4}μm

^{2}/V

^{2}

_{.}The electric field and intensity distributions of the PLSs are depicted in Figs. 4(c) and 4(d), respectively. The FWHM of the PLSs is about 92 nm, which is much larger than the FWHM of the self-focusing solitons (

*k*=

_{x}*π*/

*d*). This is because of the stronger diffraction at

*k*= 0 since the curvature of the diffraction relation at the Brillouin zone center is much larger than that at the edge, as shown in Fig. 1(b). It should be mentioned that graphene nonlinearity alone cannot balance the anomalous diffraction of SPPs at

_{x}*k*= 0.

_{x}Now we investigate the PLSs yielded by the nonlinearity of both graphene and dielectric. For the SPP mode excited at *k _{x}* =

*π*/

*d*, the normal diffraction could occur. Thus both the nonlinear effects of graphene and self-focusing dielectrics between graphene can restrain the diffraction, leading to the formation of self-focusing PLSs. The transverse electric field (

*E*) and the intensity distribution of the PLSs are shown in Figs. 5(a) and 5(b) . The FWHM of the PLSs is about 36 nm, which is slightly smaller than that for only self-focusing dielectric or graphene nonlinearity. It indicates that the self-focusing nonlinearity from graphene partly enhances the nonlinearity leading to more concentrative localized wave packet. For

_{x}*k*= 0, the SPPs experience anomalous diffraction. To yield the self-defocusing PLSs, the self-defocusing dielectrics have to be used since the nonlinearity of graphene plays a negative role. The transverse electric field and the intensity distribution are shown in Figs. 5(c) and 5(d). The FWHM of the PLSs is about 112 nm, which is remarkably larger than that for only the self-defocusing dielectric nonlinearity. The self-focusing effect of graphene weakens the self-defocusing effect of dielectric, leading to a broadened beam width of the PLSs. In addition, we obtain the values of propagation constant

_{x}*k*= 114.75 + 1.59

_{z}*i*(μm

^{−1}) for self-focusing PLSs, and

*k*= 81.71 + 2.36

_{z}*i*(μm

^{−1}) for self-defocusing PLSs. By comparing the imaginary parts of the propagation constants, the loss of self-focusing solitons is less than that of self-defocusing solitons, which accords with the propagation losses of the modes at different regions of the Brillouin zone, as shown in Fig. 1(b). The case is contrary to the result in metal-dielectric arrays [11]. Figure 6(a) illustrates the propagation of the self-focusing PLSs in the GSAs. The propagation distance is about 0.31 μm, which is larger than that of 0.21 μm of the self-defocusing solitons as shown in Fig. 6(b). In other words, the self-focusing PLSs in GSAs undergo both a smaller width and a lager propagation distance in comparison with the self-defocusing PLSs.

## 5. Tunability of the propagation of PLSs

In this section, we investigate the influence of the array period, chemical potential of graphene, and incident intensity on the self-focusing PLSs generated at *k _{x}* =

*π*/

*d*. Both the nonlinearity of graphene and self-focusing dielectrics take effect in forming the PLSs. Figure 7(a) shows the intensity distribution of PLSs as the period varies. The chemical potential of graphene is given by

*μ*= 0.15 eV. One sees from the figure that the intensity distribution tends to gradually spread out as the period increases. The field is mostly confined in the center of two graphene sheets like that shown in Fig. 5(b). Note that the structure is expanding as the period increases. At the same time, the propagation constant of self-focusing PLSs decreases as the period increases, as shown in Fig. 7(b). The real part of the propagation constant decreases dramatically as the period increases to

_{c}*d*= 30 nm. The propagation constant drops slowly when the period increases further. This is similar to the result in the linear case as shown in Fig. 1(c). The nonlinearity has no effect on the trend of propagation constant variation but make it increase overall to higher values. It could also be seen that the propagation loss has a minimum value as

*d*≈16 nm.

The dependence of the lateral profile of PLSs on the chemical potential of graphene is shown in Fig. 7(c). The period is given by *d* = 30 nm. The width of solitons keeps increasing as the chemical potential increases until it approaches to 0.35 eV. When the chemical potential increases further, the width almost keeps at a constant. It demonstrates that there is a finite small-value range of the chemical potential for PLSs flexible tunability when the structure period is fixed. At the same time, both the real part and the imaginary part of the propagation constant decrease as the chemical potential increases, as depicted in Fig. 7(d). By comparing with the linear result in Fig. 1(d), the self-focusing nonlinearity also improves overall the values of the propagation constant but the variation trend is the same.

We also consider the influence of the field intensity on the self-focusing PLSs. The effective width of the PLSs is given by *w* = (∫*x*
^{2}|**E**|^{2}
*dx*/∫|**E**|^{2}
*dx*)^{1/2}. The dependence of the effective width on the intensity is plotted in Fig. 8(a)
. One sees that the PLSs manifest a better confinement as the input light intensity increases. At the same time, the propagation constant including both the real part and the imaginary part increases as the intensity increases, as shown in Fig. 8(b). In fact, the change of the propagation constant is not remarkable since the intensity undergoes a five-time increase.

## 6. Conclusions

In conclusion, we have investigated the PLSs in GSAs by introducing the nonlinearities from both graphene and nonlinear dielectric. Due to the strong confinement of SPPs on graphene, the PLSs with deep-subwavelength effective width are obtained as the nonlinearity effect and SPPs diffraction reach an appropriate balance. The self-defocusing PLSs based on the SPPs at the Brillouin zone center can be yielded with the aid of the nonlinear effect of self-defocusing dielectrics between graphene. The self-focusing PLSs for SPPs at the Brillouin edges can be yield by the nonlinear effects of either graphene or self-focusing dielectrics. As the nonlinearities of both graphene and self-focusing dielectrics take effect, the PLSs possess smaller width and larger propagation distance. Moreover, the confinement and propagation loss of the PLSs can be optimized by changing the period of the GSAs and chemical potential of graphene. As for the self-focusing PLSs, the effective width can be less than *λ*/250. The study provides a new platform to control SPPs and will benefit the design of nonlinear GSAs and selection of materials for yielding PLSs. It will also find great applications in optical circuits, switches, and imagers on deep-subwavelength scale.

## Appendix

In this appendix, we briefly describe the simulation of field propagation by using the modified split-step Fourier beam propagation method [11]. The Maxwell’s equations can be written in the matrix form

*N*= 5 which is sufficient to include the contribution of nonlinearity.

## Acknowledgments

The authors thank Dr. Yongmin Liu for the discussion of the numerical calculations. This work is supported by the 973 Program (No. 2014CB921301), the National Natural Science Foundation of China (NSFC) (Nos. 11304108 and 11104095), and the Specialized Research Fund for the Doctoral Program of Higher Education of China (No. 20130142120091).

## References and links

**1. **V. E. Zakharov and A. B. Shabat, “Exact theory of two-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media,” Sov. Phys. JETP **34**(1), 62–69 (1972).

**2. **H. Eisenberg, Y. Silberberg, R. Morandotti, A. R. Boyd, and J. S. Aitchison, “Discrete spatial optical solitons in waveguide arrays,” Phys. Rev. Lett. **81**(16), 3383–3386 (1998). [CrossRef]

**3. **Z. Chen, M. Segev, and D. N. Christodoulides, “Optical spatial solitons: historical overview and recent advances,” Rep. Prog. Phys. **75**(8), 086401 (2012). [CrossRef] [PubMed]

**4. **S. Blair and K. Wagner, “Spatial soliton angular deflection logic gates,” Appl. Opt. **38**(32), 6749–6772 (1999). [CrossRef] [PubMed]

**5. **D. N. Christodoulides and R. I. Joseph, “Discrete self-focusing in nonlinear arrays of coupled waveguides,” Opt. Lett. **13**(9), 794–796 (1988). [CrossRef] [PubMed]

**6. **Y. S. Kivshar, “Self-localization in arrays of defocusing waveguides,” Opt. Lett. **18**(14), 1147–1149 (1993). [CrossRef] [PubMed]

**7. **D. N. Christodoulides, F. Lederer, and Y. Silberberg, “Discretizing light behaviour in linear and nonlinear waveguide lattices,” Nature **424**(6950), 817–823 (2003). [CrossRef] [PubMed]

**8. **J. W. Fleischer, T. Carmon, M. Segev, N. K. Efremidis, and D. N. Christodoulides, “Observation of discrete solitons in optically induced real time waveguide arrays,” Phys. Rev. Lett. **90**(2), 023902 (2003). [CrossRef] [PubMed]

**9. **H. S. Eisenberg, Y. Silberberg, R. Morandotti, and J. S. Aitchison, “Diffraction management,” Phys. Rev. Lett. **85**(9), 1863–1866 (2000). [CrossRef] [PubMed]

**10. **E. Feigenbaum and M. Orenstein, “Plasmon-soliton,” Opt. Lett. **32**(6), 674–676 (2007). [CrossRef] [PubMed]

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

**12. **F. Ye, D. Mihalache, B. Hu, and N. C. Panoiu, “Subwavelength plasmonic lattice solitons in arrays of metallic nanowires,” Phys. Rev. Lett. **104**(10), 106802 (2010). [CrossRef] [PubMed]

**13. **I. V. Irosh, P. A. Belov, A. A. Zharov, I. V. Shadrivov, and Y. S. Kivshar, “Nonlinear Tamm states in nanostructured plasmonic metamaterials,” Phys. Rev. A **86**(2), 023819 (2012). [CrossRef]

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

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

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

**17. **C. Berger, Z. Song, X. Li, X. Wu, N. Brown, C. Naud, D. Mayou, T. Li, J. Hass, A. N. Marchenkov, E. H. Conrad, P. N. First, and W. A. de Heer, “Electronic confinement and coherence in patterned epitaxial graphene,” Science **312**(5777), 1191–1196 (2006). [CrossRef] [PubMed]

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

**19. **A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, “The electronic properties of graphene,” Rev. Mod. Phys. **81**(1), 109–162 (2009). [CrossRef]

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

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

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

**23. **Y. Fan, B. Wang, K. Wang, H. Long, and P. Lu, “Talbot effect in weakly coupled monolayer graphene sheet arrays,” Opt. Lett. **39**(12), 3371–3373 (2014). [CrossRef] [PubMed]

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

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

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

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

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

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

**30. **R. W. Boyd and G. L. Fischer, *Nonlinear Optical Materials* (Elsevier Science Ltd, 2001).

**31. **H. Hu, K. Wang, H. Long, W. Liu, B. Wang, and P. Lu, “Precise determination of the crystallographic orientations in single ZnS nanowires by second-harmonic generation microscopy,” Nano Lett. **15**(5), 3351–3357 (2015). [CrossRef] [PubMed]

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

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

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

**35. **L. A. Falkovsky and S. S. Pershoguba, “Optical far-infrared properties of a graphene monolayer and multilayer,” Phys. Rev. B **76**(15), 153410 (2007). [CrossRef]

**36. **T. Pertsch, T. Zentgraf, U. Peschel, A. Bräuer, and F. Lederer, “Anomalous refraction and diffraction in discrete optical systems,” Phys. Rev. Lett. **88**(9), 093901 (2002). [CrossRef] [PubMed]

**37. **J. Schilling, “Uniaxial metallo-dielectric metamaterials with scalar positive permeability,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **74**(4), 046618 (2006). [CrossRef] [PubMed]

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

**39. **B. Wang, X. Zhang, K. P. Loh, and J. Teng, “Tunable broadband transmission and phase modulation of light through graphene multilayers,” J. Appl. Phys. **115**(21), 213102 (2014). [CrossRef]

**40. **S. H. Nam, A. J. Taylor, and A. Efimov, “Diabolical point and conical-like diffraction in periodic plasmonic nanostructures,” Opt. Express **18**(10), 10120–10126 (2010). [CrossRef] [PubMed]

**41. **M. L. Brongersma and P. G. Kik, *Surface Plasmon Nanophotonics* (Springer, 2007).

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

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

**44. **R. D. Schaller, M. A. Petruska, and V. I. Klimov, “Tunable near-infrared optical gain and amplified spontaneous emission using PbSe nanocrystals,” J. Phys. Chem. B **107**(50), 13765–13768 (2003). [CrossRef]

**45. **M. Jablan, H. Buljan, and M. Soljačić, “Plasmonics in graphene at infrared frequencies,” Phys. Rev. B **80**(24), 245435 (2009). [CrossRef]