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 χ⁽²⁾ and χ⁽³⁾, 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|²)E [20], where J_L = σ_gE is the conventional linear surface current with σ_g modeled by Kubo formula [21–23]. J_NL = σ_g,NL|E|²E is the Kerr nonlinear surface current with σ_g,NL being the third-order nonlinear surface conductivity which is given by [24,25]

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]

 

Fig. 1 Schematic of the nonlinear graphene multilayer. The position of graphene is denoted by n with n ∈ [1, N] and N is the total number of graphene sheets. The interlayer space is denoted by d, σ_g is the surface conductivity of graphene, and ε_d is the relative permittivity of quartz medium between graphene.

Download Full Size | PDF

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⁻⁴, 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|² = |E_x|² + |E_z|². The peak value of the intensity reads I₀ = |E|²_max, which is utilized to feature the intensity of the supermodes.

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.

 

Fig. 2 Propagation constants of the supermodes for different number of graphene sheets. (a) real and (b) imaginary parts of the propagation constants as the peak intensity I₀ = 2 V²/μm². (c) real and (d) imaginary parts of the propagation constants as the peak intensity I₀ = 40 V²/μm². (e) real and (f) imaginary parts of the propagation constants of the asymmetric supermodes. The arrows represent the increasing order of the label s or s′.

Download Full Size | PDF

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₀ = 2 V²/μm², 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₀ = 40 V²/μm². 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₀ = 40 V²/μm² 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₀ = 40 V²/μm². 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.

 

Fig. 3 Transverse electric field (E_x) distributions of the supermodes for the peak intensity I₀ = 40 V²/μm². (a) Mode profiles of four supermodes as N = 3. (b) Mode profiles of six supermodes as N = 4. The red dashed lines represent the position of graphene.

Download Full Size | PDF

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 3N/2 nonlinear supermodes as N is even and (3N−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.

 

Fig. 4 Electric field intensity (|E|²) distribution of the supermodes in the nonlinear graphene multilayer as N = 4. The field intensity is given by I₀ = 40 V²/μm².

Download Full Size | PDF

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₀ ≈34 V²/μm², 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₀ = 16 and 36 V²/μm², 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].

 

Fig. 5 (a) and (b) The nonlinear propagation constant β versus the peak intensity I₀ in the nonlinear graphene multilayers of N = 3 and 4, respectively.

Download Full Size | PDF

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.03 V²/μm² at μ_c = 0.1 eV and λ = 10 μm for the asymmetric supermode of s′ = 2. It increases to I₀ = 325 V²/μm² 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₀ = 873 V²/μm² 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.

 

Fig. 6 (a) and (b) The peak intensity threshold of the asymmetric modes versus different wavelength λ and chemical potential μ_c as N = 4. (a) The threshold I_c for mode of s′ = 1. (b) The threshold I_c for mode of s′ = 2.

Download Full Size | PDF

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.