Enhanced plasmonic nanofocusing of terahertz waves in tapered graphene multilayers

We investigate the plasmonic nanofocusing of terahertz waves in tapered graphene multilayers separated by dielectrics. The nanofocusing effect is significantly enhanced in the graphene multilayer taper compared with that in a single layer graphene taper due to interlayer coupling between surface plasmon polaritons. The results are optimized by choosing an appropriate layer number of graphene and the field amplitude has been enhanced by 620 folds at λ = 50 μm. Additionally, the structure can slow light to a group velocity ~1/2815 of the light speed in vacuum. Our study provides a unique approach to compress terahertz waves into deep subwavelength scale and may find great applications in terahertz nanodevices for imaging, detecting and spectroscopy. ©2016 Optical Society of America OCIS codes: (250.5403) Plasmonics; (040.2235) Far infrared or terahertz; (310.6628) Subwavelength structures, nanostructures; (230.4170) Multilayers; (230.7370) Waveguides. References and links 1. T. L. Cocker, V. Jelic, M. Gupta, S. J. Molesky, J. A. J. Burgess, G. De Los Reyes, L. V. Titova, Y. Y. Tsui, M. R. Freeman, and F. A. Hegmann, “An ultrafast terahertz scanning tunnelling microscope,” Nat. Photonics 7(8), 620–625 (2013). 2. K. Moon, H. Park, J. Kim, Y. Do, S. Lee, G. Lee, H. Kang, and H. Han, “Subsurface nanoimaging by broadband terahertz pulse near-field microscopy,” Nano Lett. 15(1), 549–552 (2015). 3. F. D’Angelo, Z. Mics, M. Bonn, and D. Turchinovich, “Ultra-broadband THz time-domain spectroscopy of common polymers using THz air photonics,” Opt. Express 22(10), 12475–12485 (2014). 4. P. K. Mandal and V. Chikan, “Plasmon-phonon coupling in charged n-type CdSe quantum dots: A THz time-domain spectroscopic study,” Nano Lett. 7(8), 2521–2528 (2007). 5. S. Koenig, D. Lopez-Diaz, J. Antes, F. Boes, R. Henneberger, A. Leuther, A. Tessmann, R. Schmogrow, D. Hillerkuss, R. Palmer, T. Zwick, C. Koos, W. Freude, O. Ambacher, J. Leuthold, and I. Kallfass, “Wireless sub-THz communication system with high data rate,” Nat. Photonics 7(12), 977–981 (2013). 6. Y. H. Li, A. Rashidinejad, J. M. Wun, D. E. Leaird, J. W. Shi, and A. M. Weiner, “Photonic generation of W-band arbitrary waveforms with high time-bandwidth products enabling 3.9 mm range resolution,” Optica 1(6), 446–454 (2014). 7. J. H. Jeong, B. J. Kang, J. S. Kim, M. Jazbinsek, S. H. Lee, S. C. Lee, I. H. Baek, H. Yun, J. Kim, Y. S. Lee, J. H. Lee, J. H. Kim, F. Rotermund, and O. P. Kwon, “High-power Broadband Organic THz Generator,” Sci. Rep. 3, 3200 (2013). 8. S. D. Jackson, “Towards high-power mid-infrared emission from a fibre laser,” Nat. Photonics 6(7), 423–431 (2012). 9. J. Dai and X. C. Zhang, “Terahertz wave generation from thin metal films excited by asymmetrical optical fields,” Opt. Lett. 39(4), 777–780 (2014). 10. S. A. Maier, S. R. Andrews, L. Martín-Moreno, and F. J. García-Vidal, “Terahertz surface plasmon-polariton propagation and focusing on periodically corrugated metal wires,” Phys. Rev. Lett. 97(17), 176805 (2006). 11. A. Rusina, M. Durach, K. A. Nelson, and M. I. Stockman, “Nanoconcentration of terahertz radiation in plasmonic waveguides,” Opt. Express 16(23), 18576–18589 (2008). 12. W. B. Qiu, X. H. Liu, J. Zhao, S. H. He, Y. H. Ma, J. X. Wang, and J. Q. Pan, “Nanofocusing of mid-infrared electromagnetic waves on graphene monolayer,” Appl. Phys. Lett. 104(4), 041109 (2014). 13. C. H. Gan, “Analysis of surface plasmon excitation at terahertz frequencies with highly doped graphene sheets via attenuated total reflection,” Appl. Phys. Lett. 101(11), 111609 (2012). #261950 Received 28 Mar 2016; revised 15 Jun 2016; accepted 17 Jun 2016; published 22 Jun 2016 © 2016 OSA 27 Jun 2016 | Vol. 24, No. 13 | DOI:10.1364/OE.24.014765 | OPTICS EXPRESS 14765 14. A. Toma, S. Tuccio, M. Prato, F. De Donato, A. Perucchi, P. Di Pietro, S. Marras, C. Liberale, R. Proietti Zaccaria, F. De Angelis, L. Manna, S. Lupi, E. Di Fabrizio, and L. Razzari, “Squeezing terahertz light into nanovolumes: nanoantenna enhanced terahertz spectroscopy (nets) of semiconductor quantum dots,” Nano Lett. 15(1), 386–391 (2015). 15. M. A. Seo, H. R. Park, S. M. Koo, D. J. Park, J. H. Kang, O. K. Suwal, S. S. Choi, P. C. M. Planken, G. S. Park, N. K. Park, Q. H. Park, and D. S. Kim, “Terahertz field enhancement by a metallic nano slit operating beyond the skin-depth limit,” Nat. Photonics 3(3), 152–156 (2009). 16. M. Schnell, P. Alonso-Gonzalez, L. Arzubiaga, F. Casanova, L. E. Hueso, A. Chuvilin, and R. Hillenbrand, “Nanofocusing of mid-infrared energy with tapered transmission lines,” Nat. Photonics 5(5), 283–287 (2011). 17. A. A. Govyadinov and V. A. Podolskiy, “Metamaterial photonic funnels for subdiffraction light compression and propagation,” Phys. Rev. B 73(15), 155108 (2006). 18. F. Bonaccorso, Z. Sun, T. Hasan, and A. C. Ferrari, “Graphene photonics and optoelectronics,” Nat. Photonics 4(9), 611–622 (2010). 19. L. Ju, B. Geng, J. Horng, C. Girit, M. Martin, Z. Hao, H. A. Bechtel, X. Liang, A. Zettl, Y. R. Shen, and F. Wang, “Graphene plasmonics for tunable terahertz metamaterials,” Nat. Nanotechnol. 6(10), 630–634 (2011). 20. Q. L. Bao, H. Zhang, B. Wang, Z. H. Ni, C. H. Y. X. Lim, Y. Wang, D. Y. Tang, and K. P. Loh, “Broadband graphene polarizer,” Nat. Photonics 5(7), 411–415 (2011). 21. A. Vakil and N. Engheta, “Transformation optics using graphene,” Science 332(6035), 1291–1294 (2011). 22. S. Thongrattanasiri and F. J. García de Abajo, “Optical field enhancement by strong plasmon interaction in graphene nanostructures,” Phys. Rev. Lett. 110(18), 187401 (2013). 23. A. Woessner, M. B. Lundeberg, Y. Gao, A. Principi, P. Alonso-González, M. Carrega, K. Watanabe, T. Taniguchi, G. Vignale, M. Polini, J. Hone, R. Hillenbrand, and F. H. L. Koppens, “Highly confined low-loss plasmons in graphene-boron nitride heterostructures,” Nat. Mater. 14(4), 421–425 (2014). 24. A. K. Geim and K. S. Novoselov, “The rise of graphene,” Nat. Mater. 6(3), 183–191 (2007). 25. V. W. Brar, M. S. Jang, M. Sherrott, J. J. Lopez, and H. A. Atwater, “Highly confined tunable mid-infrared plasmonics in graphene nanoresonators,” Nano Lett. 13(6), 2541–2547 (2013). 26. P. Li, T. Wang, H. Böckmann, and T. Taubner, “Graphene-enhanced infrared near-field microscopy,” Nano Lett. 14(8), 4400–4405 (2014). 27. B. Zhu, G. Ren, Y. Gao, Y. Yang, Y. Lian, and S. Jian, “Graphene-coated tapered nanowire infrared probe: a comparison with metal-coated probes,” Opt. Express 22(20), 24096–24103 (2014). 28. 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). 29. J. Christensen, A. Manjavacas, S. Thongrattanasiri, F. H. L. Koppens, and F. J. G. de Abajo, “Graphene plasmon waveguiding and hybridization in individual and paired nanoribbons,” ACS Nano 6(1), 431–440 (2012). 30. 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). 31. S. Ke, B. Wang, H. Huang, H. Long, K. Wang, and P. Lu, “Plasmonic absorption enhancement in periodic cross-shaped graphene arrays,” Opt. Express 23(7), 8888–8900 (2015). 32. P. Y. Chen and A. Alù, “Atomically thin surface cloak using graphene monolayers,” ACS Nano 5(7), 5855–5863 (2011). 33. D. K. Gramotnev and S. I. Bozhevolnyi, “Nanofocusing of electromagnetic radiation,” Nat. Photonics 8(1), 14–23 (2014). 34. B. F. Zhu, G. B. Ren, Y. X. Gao, Y. Yang, B. L. Wu, Y. D. Lian, and S. S. Jian, “Nanofocusing in the graphene-coated tapered nanowire infrared probe,” J. Opt. Soc. Am. B 32(5), 955–960 (2015). 35. A. Wiener, A. I. Fernández-Domínguez, A. P. Horsfield, J. B. Pendry, and S. A. Maier, “Nonlocal effects in the nanofocusing performance of plasmonic tips,” Nano Lett. 12(6), 3308–3314 (2012). 36. H. Choo, M. K. Kim, M. Staffaroni, T. J. Seok, J. Bokor, S. Cabrini, P. J. Schuck, M. C. Wu, and E. Yablonovitch, “Nanofocusing in a metal-insulator-metal gap plasmon waveguide with a three-dimensional linear taper,” Nat. Photonics 6(12), 838–843 (2012). 37. E. Verhagen, A. Polman, and L. K. Kuipers, “Nanofocusing in laterally tapered plasmonic waveguides,” Opt. Express 16(1), 45–57 (2008). 38. J. A. Gerber, S. Berweger, B. T. O’Callahan, and M. B. Raschke, “Phase-resolved surface plasmon interferometry of graphene,” Phys. Rev. Lett. 113(5), 055502 (2014). 39. E. D. Palik, Handbook of Optical Constants of Solids I (Academic, 1985). 40. R. W. Yu, R. Alaee, F. Lederer, and C. Rockstuhl, “Manipulating the interaction between localized and delocalized surface plasmon-polaritons in graphene,” Phys. Rev. B 90(8), 085409 (2014). 41. V. Atanasov and A. Saxena, “Electronic properties of corrugated graphene: the Heisenberg principle and wormhole geometry in the solid state,” J. Phys-Condens. Mat. 23(17), 175301 (2011). 42. V. Atanasov and A. Saxena, “Tuning the electronic properties of corrugated graphene: confinement, curvature, and band-gap opening,” Phys. Rev. B 81(20), 205409 (2010). 43. P. Tassin, T. Koschny, M. Kafesaki, and C. M. Soukoulis, “A comparison of graphene, superconductors and metals as conductors for metamaterials and plasmonics,” Nat. Photonics 6(4), 259–264 (2012). 44. M. I. Stockman, “Nanofocusing of optical energy in tapered plasmonic waveguides,” Phys. Rev. Lett. 93(13), 137404 (2004). #261950 Received 28 Mar 2016; revised 15 Jun 2016; accepted 17 Jun 2016; published 22 Jun 2016 © 2016 OSA 27 Jun 2016 | Vol. 24, No. 13 | DOI:10.1364/OE.24.014765 | OPTICS EXPRESS 14766 45. F. De Angelis, G. Das, P. Candeloro, M. Patrini, M. Galli, A. Bek, M. Lazzarino, I. Maksymov, C. Liberale, L. C. Andreani, and E. Di Fabrizio, “Nanoscale chemical mapping using three-dimensional adiabatic compression of surface plasmon polaritons,” Nat. Nanotechnol. 5(1), 67–72 (2010). 46. S. He, X. Zhang, and Y. He, “Graphene nano-ribbon waveguides of record-small mode area and ultra-high effective refractive ind


Introduction
Terahertz (THz) waves categorize between the far infrared and millimeter radio waves, which have shown great applications in multiple fields such as THz imaging [1,2], spectroscopy analysis [3,4] and wireless communications [5,6].However, there exist considerable challenges to develop THz sources with high power density and nanoscale size [7][8][9].This issue seriously hinders the applications of THz waves in nonlinear optics and nanoscale devices.Recently, nanofocusing of THz waves based on surface plasmon polaritons (SPPs) has attracted great attentions due to the possibility of field enhancement and light manipulation below the diffraction limit [10][11][12].SPPs on metals such as gold and silver manifest weak field confinement at THz frequencies because the frequencies are far away from that of the bulk plasmon oscillation of metals [12,13].To bridge the mismatch, metal and metal/dielectric metastructures have been proposed to generate spoof SPPs, such as nanoantenna [14], corrugated metal wires [10], nanoslit [15,16] and metamaterial photonic funnels [17].However, the large loss and weak confinement of metal plasmons have prevented further enhancement and compression of the field [16].
Recently, a two-dimensional semi-metal material with zero bandgap, known as graphene, has become a fascinating candidate for supporting infrared SPPs due to its unique electronic and photonic properties [18][19][20].The TM-polarized SPPs can be tightly confined on graphene surface with the SPP wavelength being much smaller than the incident one, which has great potential for localizing THz wave on subwavelength scale and achieving strong field enhancement [21][22][23].Another unique feature of graphene lies in that the surface conductivity can be conveniently tuned by gate voltage, electric field and chemical doping [24,25].Therefore, it has attracted great attentions for plasmonic nanofocusing of THz wave in monolayer graphene structures.For example, the monolayer graphene sheets can offer a 7-fold enhancement of the mid-infrared evanescent field [26] and the THz waves can experience a field enhancement of ten times stronger in the graphene coated taper probe than in the metal coated ones [27].For graphene multilayer nanostructures, there have been demonstrated to be some novel properties due to strong field coupling, such as low loss plasmonic supermode and negative coupling effect [28][29][30].
In this work, we study the plasmonic nanofocusing of THz waves in tapered graphene multilayers.We firstly derive a theoretical model to describe the SPPs propagation in the tapered graphene waveguide with adiabatic approximation, in which the validity of the adiabatic approximation is improved by taking the end reflection of the graphene waveguide into account.Combining with the numerical calculated results, it is found that the field confinement can be largely increased at the taper tip and the amplitude of THz wave can be enhanced by a factor of 620 at λ = 50 μm, both have been greatly enhanced compared with those in single layer graphene taper.As illustrated in Fig. 1(a), a graphene multilayer scroll (multiple graphene monolayers separated by dielectric) is tapered to concentrate the incident THz waves propagating along z-axis.Figure 1(b) shows the xz and xy cross sections of the graphene taper.The taper is geometrically characterized by the following parameters (labelled in the xz cross section): the input aperture diameter (2R a ), the output aperture diameter (2R b1 for the inner layer and 2R b for the outer layer), taper length (L), and graphene layer number (N).For a graphene multilayer taper as N > 1, the graphene interlayer space at the input and output port can be calculated as

Theoretical model
, respectively.Graphene is treated as an infinitesimally thin layer and the graphene conductivity σ g , is governed by the Kubo formula [31,32].As shown in Fig. 1(b), for a taper composed of N graphene monolayers (Black labels), the space can be divided into N + 1 areas (Red labels).The relative permittivity of the n-th area is assumed to be ε n (0≤ n ≤N).The dispersion relation of the graphene multilayer scroll can be obtained by the transfer-matrix (see Appendix for details): where n represents the graphene layer number, 1≤n≤N, κ i = (β 2 -ε i k 0 2 ) 1/2 with β being the wave vector of SPPs in the z direction, k 0 = 2π/λ, and n eff = β/k 0 being the effective mode index.I m (κr) and K m (κr) are the modified Bessel functions of the first and second kind (m = 0, 1).A m and B m are arbitrary constants determined by boundary conditions ).The adiabatic approximation (WKB method) is commonly employed to deal with SPPs propagating along the tapered structures, and the theory assumes no reflection from the end of the waveguide [33][34][35].Actually, the end reflections could become noticeable because there exists strong impedance mismatch at the waveguide/air interfaces, which have been demonstrated in previous experimental reports [36][37][38].By solving the Helmholtz Equation with adiabatic approximation and taking the end reflections of the taper into our theoretical model, we can obtain the formula for electric/magnetic field distributions when SPPs propagate in the tapered graphene multilayer (see Appendix for details): where Φ(x,y,z) and Φ 0 (x,y,z) denote the corresponding fields (E r , E z or H φ ) of the graphene taper and graphene scroll, respectively.Φ 0 (x,y,z) can be expressed by Eqs. ( 6)-( 8) in the Appendix.It is "+" for E r and "-" for E z and H φ according to the coordinate system.r a and r b are the reflection coefficients of the input and output ports respectively.

Plasmonic nanofocusing in single layer graphene taper
Now we study the nanofocusing properties of the single layer graphene taper.In the following study, we assume that the taper is filled with Si and surrounded by air (ε 1 = ε 2 = … = ε N-1 = ε Si = 11.7 [39], ε N = 1).The incident THz wavelength is set as λ = 50 μm.In order to efficiently couple the light into the graphene taper, the input aperture diameter is set as 2R a = 2 μm, comparable to half of the SPPs modal wavelength (λ SPP /2 ≈1.5 μm).The output diameter is set as 2R b = 10 nm [35,36].The chemical potential of graphene is set as μ c = 0.5 eV and the corresponding relaxation time τ = 0.5 ps at 300 K [40].It should be noted that the graphene bending induced geometric potential is considered in the calculation [41,42], as analyzed in the Appendix.Figure 2(a) shows the variation of the effective index along the propagation direction.We can see that the real part (Re[n eff ]) decreases gradually at first and then rises rapidly near the taper tip, approaching to a maximum of 30 at the tip.The mode lateral extent at the tip, which is used to characterize the field confinement, is given by [43] we arrive at l SPP = 265 nm by using above parameters, demonstrating a subwavelength localization of the THz wave by plasmonic nanofocusing.It also appears that the imaginary part (Im[n eff )] manifests similar behavior along the propagation direction.To see the SPP field more clearly near the taper tip, we calculate the field distributions by numerical method based on the commercial finite-element software (COMSOL Multiphysics).It has been a great challenge to couple the light into graphene plasmons [38].In order to efficiently excite the graphene SPPs, eigenmodes of the input port (solved by Eqs. ( 6)-( 8) in Appendix) are used as the incident field in the calculation.The maximum mesh element size on the taper tip was set to be smaller than ΔR b /10 to obtain a satisfactory numerical convergence.The field distributions of the dominating transverse component E r and the normalized field ) over a small area (labelled by dashed rectangle) at the tip are shown in Fig. 2(c).|E 0 | is defined as the amplitude maximum of the incident field and the incident power is normalized to be 1 W [27].It is observed that the SPP field dominantly distributes on the exterior surface of the graphene single layer.There is clearly a concentration of the field as a hotspot is created when the wave is guided to the taper tip, indicating the occurrence of plasmonic nanofocusing.The lateral extension of the hotspot is observed to be more than two hundred nanometers, which agrees well with the theoretical calculation.To quantitatively determine the field enhancement, the field distribution on the graphene exterior surface is calculated with adiabatic approximation and compared with numerical simulations.The reflection coefficients of the taper ports are calculated by solving S-parameter in Comsol Multiphysics (see Appendix).The results are shown in Fig. 2(d).As expected, the field amplitude at the tip has been enhanced for more than two orders of magnitude due to plasmonic nanofocusing.The field variation profiles of the theoretical and numerical results are in accordance, validating our theoretical model in characterizing the propagation and nanofocusing of SPPs.In fact, the adiabatic approximation is valid under the condition |d(Q 1 −1 )/dz| << 1, where |d(Q 1 −1 )/dz| is the well-known adiabatic parameter with Q 1 = Re(β) [33,44].Figure 2(b) plots the variation of |d(Q 1 −1 )/dz| along the propagation direction.The adiabatic parameter becomes much larger than 1 near the taper tip, indicating that there exists reflection in the adiabatic region, but it can still be negligible compared with the complete reflection at the graphene/air boundaries [21,38].As a result, the validity of the adiabatic approximation has been expanded, which can be confirmed by the great accordance between the theoretical and numerical results.In addition, the numerical calculated field amplitude is larger than the theoretical result at the taper tip, because there exists an evanescent field enhancement due to dielectric mismatch at the end of taper, which is not considered in the theoretical model.To get more accurate results for the field enhancement, the following studies are mainly performed based on numerical simulations.
The nanofocusing depending on the structure parameters of the single layer graphene taper are investigated here.Figure 2(e) shows the field amplitude enhancement at the taper tip (|E tip |/|E 0 |) as a function of the tip radius R b .For practical applications, the taper length is selected to be 3 μm and the smallest tip radius is set to be 5 nm [36,45].It can be observed that the field is strongly enhanced, especially for R b < 20 nm.The inset shows the mode lateral extent l spp as a function of the tip radius R b .As the tip radius increases, the field extent at the taper tip increases gradually, resulting in the decrease of field nanofocusing.It indicates that the strongly enhanced field is originated from the high spatial localization of SPPs at the taper tip.For a graphene scroll without taper, the plasmonic nanofocusing will completely vanish.Figure 2(f) shows the field amplitude enhancement as a function of the taper length L. The field amplitude is oscillated with the taper length.The oscillation is due to the interference of SPPs in the tapered waveguide and the peak value can be obtained when the taper length is on resonance.The resonant lengths are calculated to be 2.3 μm, 5.1 μm and 7.2 μm by the resonant condition 2Re(k spp (z))dz + φ = 2mπ (0≤z≤L).φ is the phase shift caused by the reflection at the input and output ports and the value can be obtained from the reflection coefficients.We can see the theoretical calculated resonant lengths agree well with the numerical data.The overall profile of the result shows that the largest field amplitude enhancement is about 240 by nanofocusing in the single layer graphene taper.In multilayer graphene structures, there exist multiple eigenmodes due to field coupling between the SPPs from individual layers [28].Figure 3(a)-3(c) show the mode profiles of E r (Eq.( 7) in Appendix) for the 1-, 2-, and 3-layer scroll respectively.The outer and inner layer radii are set as R b = 5 nm and R b1 = 4 nm, which are used for the output aperture radii of the graphene multilayer tapers studied below.For a specific graphene scroll, the number of the calculated effective indices and modes is equal to the layer number N, similar with that in a graphene multilayer sheet array [28].Particularly, we can see that the mode with the largest Re(n eff ) can be best confined inside the scroll, while the mode with the smallest Re(n eff ) cannot be confined well, indicating a great distinction between the different modes for plasmonic nanofocusing.To summarize, the modal eigenvalues (effective indices) for the N-layer scrolls (N from 1 to 10) are plotted in Fig. 3(d) (real part) and 3(e) (imaginary part).The modes are labelled as s = 1, 2, ...N, in terms of the increasing of Re(n eff ).We focus on the modes with the largest Re(n eff ) (s = N), because they have the best confinement as mentioned above.As shown in Fig. 3(d), Re(n eff ) for mode s = N (red circles) significantly increases from 30 to more than 2800, which is due to the interlayer coupling between SPPs so that the THz field can be strongly confined in the taper [46].Correspondingly, the mode lateral extent reduces from 265 nm to 2.8 nm according to Eq. ( 5), suggesting the field compression is significantly enhanced in the graphene multilayer taper.Compared with the plasmonic nanofocusing in metal-dielectric and graphene monolayer structures, in which the field can be compressed to about several tens to several hundreds of nanometers, the field compression in the graphene multilayer structure has been greatly enhanced [11,16,25,47].This result is of particular significance since the graphene multilayer taper could achieve much better field confinement and thus having superior potentials in obtaining an enhanced plasmonic nanofocusing.To demonstrate the predicted enhancement, the plasmonic nanofocusing of THz wave is investigated in the graphene multilayer taper.The two-layer structure is firstly studied, in which there exist two modes due to different mode couplings of SPPs. Figure 4(a) and 4(b) show the variation of Re(n eff ) and Im(n eff ) for the two modes (green for s = 1 and red for s = 2), respectively.In particular, Re(n eff ) of mode s = 2 increases to a maximum of 680 at the tip, obviously much larger than that of mode s = 1.Correspondingly, the compressed mode lateral extension at the taper tip for s = 2 is calculated to be 11.7 nm. Figure 4(c) shows the adiabatic parameter |d(Q 1 −1 )/dz| varies along the propagation direction for s = 2.The variation becomes much slower than that in the single layer structure, implying that the SPP wavelength changes more slowly in the graphene multilayer taper.The normalized field distributions of E r component and |E| for the graphene two-layer taper are illustrated in Fig. 4(d).The most notable feature is that the SPPs transport to the taper end with optical field dominantly confined between the graphene layers.While for the monolayer graphene or metal tapers, the field cannot be confined well in such a nanoscale area [28,45].It further confirms that the field localization in graphene multilayer taper has been strongly enhanced, making it possible to obtain a largely enhanced nanofocusing of THz waves.The oscillation pattern in the field distribution of |E|/|E 0 | is due to constructive and destructive interference between the forward and backward propagating SPPs, as mentioned above.Figure 4(e) shows the theoretical and numerical calculated field amplitude (|E|/|E 0 |) profiles on the surface of the inner layer graphene.The theoretical and numerical results are in accordance that there is a large field enhancement at the taper tip as well.The field amplitude of the numerical result is larger at the tip, the same with that in the single layer graphene taper.In addition, the field enhancement for mode s = 2 is calculated to be remarkably larger than that of s = 1, as shown in the Appendix.Figure 4(f) shows the field amplitude enhancement as a function of the taper length L for s = 2 as ΔR b = 1 nm and 2 nm.The overall enhancement for ΔR b = 1 nm is larger than that of ΔR b = 2 nm, because the SPP field can be localized in a much smaller area at smaller ΔR b .For each ΔR b , more than one resonant peaks present over the range of L calculated, and the field enhancement decreases as L increases.To compare the nanofocusing effect between the N-layer (N = 1, 2, 3, …) graphene taper reasonably, the lengths of all the tapers are chosen to fulfill the resonance condition and to be those nearest to 3 μm.For example, the length is about 2.85 μm as N = 2.The corresponding field enhancement is about 470, which is about 200 larger than that of the single layer graphene taper.

Enhanced plasmonic nanofocusing in graphene multilayer taper
The field amplitude enhancement for N = 3, 4 and 5 is also calculated as functions of L, respectively.The results are shown in Fig. 5(a)-5(c) and the field distributions at the corresponding resonant peaks are shown in Fig. 5(d)-5(f).We can see that the field distributions in the graphene multilayer tapers (N = 3, 4, 5) are similar with that of N = 2.As layer number N increases, less optical energy is observed to leak out from the taper tips.The field lateral extents are calculated to be 6.8 nm, 5.2 nm and 4.4 nm for N = 3, 4, 5 respectively, which indicates that the field confinement becomes stronger as the graphene layer number increases.Moreover, there are four peaks as N = 5, while there are only three as N = 3, over the same area labelled by black dashed rectangle in Fig. 5(d).The reason is that, as the number of graphene layer increases, the real part of the effective index increases, and results in the decrease of spatial oscillation period.It suggests that the SPP field becomes more sensitive to the taper length as the graphene layer increases.
The field enhancement for the mode s = N as a function of the layer number N is shown in Fig. 5(g).It reaches the maximum at N = 3 and then decreases.The enhanced nanofocusing of SPPs in the tapered graphene multilayer originates from the balance between field concentration and field damping.Figure 5(h) shows inverse of the normalized effective mode area A eff (1/A eff , red curve) at the taper tip and the propagation loss η (green curve), characterizing field energy concentration and damping respectively, as functions of the graphene layer N. The normalized effective mode area is defined as A eff = W(x,y)dxdy/(λ 2 W max ), where W(x,y) is the energy density of the waveguide cross section [27].η is defined as the field energy attenuation as the SPPs propagate from the input to the output port, i.e., η = exp(−k 0 Im(n eff )dz).We can see that as N increases, both the field energy concentration and attenuation keep increasing and tend to be saturated.Compared with the single layer graphene taper (A eff = 2.6 × 10 −7 ), the effective mode area of the multilayer taper experiences a great decrease (A eff = 1.1 × 10 −8 for N = 2), indicating a largely enhanced field nanofocusing.For the graphene few-layer tapers, the field concentration increases faster than that of field damping, resulting in the increase of field enhancement.As N further increases, the increase of field damping becomes more and more noticeable and goes beyond that of field concentration, thus the field enhancement presents to be decreased with N. Generally, the largest field amplitude enhancement can be obtained to be 620 as N = 3, which has been increased for more than 380 compared with that in single layer graphene taper.
We also pay attention to the greatly increased effective indices in the tapered graphene multilayers.Re(n eff ) has a large increase from 30 to about 2800 for N = 10 at the taper tip.Correspondingly, the SPP wavelength λ SPP can be compressed to 1/2800 of the vacuum wavelength.The group velocity of THz wave in the waveguide, calculated by ν g = c/[d(nω)/dω], is slowed down to ~1/2815 of the light speed in vacuum [48].Such strong confinement and slow velocity of SPPs in the graphene multilayer taper is of great promise for applications in THz nanophotonics and slow-wave optics.It is well-known that the optical properties of graphene can be conveniently tuned by varying the chemical potential with chemical doping and gate voltage.In order to expand the investigation to a wider THz frequency range, the field amplitude enhancement |E tip |/|E 0 | is mapped as functions of THz wavelength and graphene chemical potential for N = 1, 2, and 3 respectively, as shown in Fig. 6.For a fixed graphene layer N, the field enhancement on resonance increases with the wavelength and chemical potential, because the SPP propagation loss in graphene waveguide is inversely proportional to both of them.The overall results indicate that the best nanofocusing can be obtained at N = 2 for longer THz wavelength and larger chemical potential, and the largest field amplitude enhancement is obtained to be about 1800.It means that the field intensity, that scales as the square of field amplitude, can be enhanced by 5 to 6 orders of magnitude.If these fields are used for nonlinear optics, the field intensity would be enhanced by more than 10 orders of magnitude by plasmonic nanofocusing in graphene multilayer taper.

Conclusion
In conclusion, the plasmonic nanofocusing of THz wave in tapered graphene multilayer is studied.Compared with the single layer graphene taper, the field localization and field intensity of THz wave have been significantly enhanced due to interlayer coupling between SPPs in the graphene multilayer taper.The largest field amplitude enhancement is obtained to be 620 at λ = 50 μm, which increases further to 1800 at longer THz wavelength and larger graphene chemical potential.In addition, the plasmonic wavelength can be reduced to 1/2800 of the THz wavelength and the group velocity can be slowed down to 1/2815 of the light speed in vacuum at the tip of graphene multilayer taper.
The key practical concerns lie in the fabrications and applications of such a graphene multilayers coated taper.It is noticed that there have been several techniques to coat the graphene monolayers on nanofibers [49], AFM probes [50] and some other nanostructures [51].Based on these techniques, a bottom-up fabrication method that coating the graphene monolayers layer by layer will be practicable.For applications like efficient nonlinear detection and imaging [52,53], it is desirable to perform with a strong and nanometer-sized field.As concluded above, plasmonic nanofocusing of THz wave in the graphene multilayer taper provides a unique approach to largely enhance and compress the field, which will have considerable potential for THz nonlinear optics.
Graphene is treated as an infinitesimally thin layer with surface conductivity σ g characterized by Kubo formula.The analysis is performed on the fundamental transverse magnetic (TM) mode because it is always supported in graphene without cutting off.The non-zero components of the TM mode are E z , E r , H φ .E z component in the region between the (n -1)th and (n + 1)th graphene scroll can be written as the superposition of electromagnetic (EM) waves in opposite directions (as shown in Fig. 7(a)): where ) 1/2 being the transversal propagation, k 0 = 2π/λ being the wave number in free space, and β = n eff k 0 being the complex propagation constant of the SPP along z axis.n eff is the effective mode index and ε i is the dielectric constant in the i-th area.According to the Maxwell Equation, E r and H φ can be respectively expressed by: Since the tangential components of the EM fields at graphene layer is continuous, the boundary conditions at the n-th graphene layer (r = R n ) can be expressed as ). Substituting Eqs. ( 6) and ( 8) into the boundary conditions, the field amplitude coefficients in the neighboring areas can be written as the transfer matrix: , where 0 0 n is the layer number of the graphene, 1≤n≤N.Specially, n = 1 represents the single layer waveguide.According to our definition that A 0 = 0, B N = 0.With Eqs. ( 6)-( 11), the dispersion relationship and the transversal mode field in graphene scroll waveguide can be obtained.The adiabatic approximation (WKB method) is usually employed to deal with the optical field propagating along a tapered waveguide [33][34][35].To derive the final formula, we start with the Helmholtz Equation 22 2 0 ( , , ) ( , , ) ( , , ) 0.
For variables separation we assume E z (x, y, z) = Ψ(x, y, z)Θ(z), where Ψ(x, y, z) denotes the mode profile varying slowly with the z axis.Considering that steady.The field at arbitrary position z can be equivalent to the superposition of two parts (shown in Fig. 7(b)): the part from the input port (E a ), which propagates forward along z axis; the other part from the output port (E b ), which is originated from the input part and propagates backward along z axis.According to eq. ( 18), the relationship between E a and E b could be expressed by (0) 0 0 0 ( ) where L is the waveguide length, Φ 0 is the modal field of the incident light.The field with a distance z from the input port can be expressed by With Eqs. ( 19) and ( 20), we arrive at the final expression for the SPP field propagating in the graphene tapered waveguide ] ( , , ), 1 exp( 2) where Φ(x, y, z) and Φ 0 (x, y, z) denote the corresponding field (E r , E z or H φ ) of the graphene multilayer taper and scroll, respectively.It is "+" for E r and "-" for E z and H φ according to the coordinate system.

Influence of graphene bending induced geometric potential on the plasmonic nanofocusing
Monolayer graphene sheet is a two-dimensional material that the electrons can be perfectly confined in the surface.When the graphene layer is curved into scroll, the electronic band structure will be affected and result in the geometric potential.According to the previous study [41,42], the geometric potential scales linearly with the bending curvature and shifts the Fermi energy of the graphene.The shifted chemical potential is expressed as where μ c0 is the chemical potential of the graphene sheet, v F ≈c/300 is the Fermi velocity of the carriers in graphene, and R is the curvature radius of the graphene scroll.c is the velocity of light in vacuum.Figure 8(a) plots the chemical potential of the graphene scroll varies as function of the curvature radius.The initial chemical potential is set to be 0.5 eV.We can see that when the curvature radius is large enough, the influence of the bending induced geometric potential is negligible (e.g.Δμ c = 6 meV for R = 100 nm).When R is less than 100 nm, the geometric potential decreases rapidly as the radius decreases and the influence on the graphene chemical becomes distinctive.To investigate the influence of geometric potential, the nanofocusing of the single layer graphene taper with and without considering the geometric potential are calculated for comparison.Fig. 8(b) and Fig. 8(c) show the plots of effective indices and nanofocusing for a single layer graphene taper with and without considering the influence of geometric potential respectively.We can see that the propagation loss is a little larger when the bending effect is considered, due to the decreased chemical potential.As a result, the field enhancement is a little smaller when the geometrical potential is considered.Obviously, it will be more reasonable to take the influence of the geometric potential into account for practical applications.

Reflection coefficients of the graphene scroll
As the air does not support SPPs propagating, the wave will be strongly reflected with a little radiation loss at the graphene/air interface [21,38].The reflection coefficients for the graphene single layer and multilayer scroll are calculated by solving the S-parameter in COMSOL Multiphysics., where R 1 represents the inner layer radius and R 2 represents the interlayer space.The port reflectance is nearly to be 1 for the both scrolls, which confirms that the SPP is almost completely reflected at the graphene/air interface, in accordance with that in monolayer graphene sheet [38].

Comparison of nanofocusing in two-layer graphene taper for modes s = 1 and s = 2
As mentioned, there are two eigenmodes for a graphene two-layer taper, labeled as s = 1 and s = 2. Figures 10(a)-10(d) show the normalized field distributions of E r component and |E|/|E 0 | for the two modes, respectively.Obviously, the field of s = 1 dominantly distributes on the exterior surface of the taper, which can not be confined well.On contrast, the field of s = 2 is strongly confined at the taper tip. Figure 10(e) shows the effective indices varying long the tapered direction for the two modes.The mode lateral extents at the taper tip, are calculated to be 361.7 nm and 11.7 nm for s = 1 and s = 2 respectively, indicating a much better field localization of s = 2.To accurately compare the plasmonic nanofocusing effect, the normalized field distributions (|E|/|E 0 |) of the two modes on the inner-layer graphene surface are calculated by simulation, as shown in Figure 10(f).The field amplitude enhancement of mode s = 2 is more than 100 larger than that of s = 1, which further confirms the better nanofocusing of mode s = 2.

Fig. 1 .
Fig. 1.(a) Schematic for plasmonic nanofocusing of THz wave in tapered graphene multilayers.(b) xz and xy cross sections of the graphene multilayer taper.The geometric parameters of the taper are labelled in the xz cross section.The graphene monolayers and dielectric areas are respectively labelled with black and red numbers in the xy cross section.

Fig. 2 .
Fig. 2. Plasmonic nanofocusing of THz wave in a single layer graphene taper.(a) Plot of calculated effective index along the taper.The red and green curves represent the real part (Re[n eff ]) and imaginary part (Im[n eff ]), respectively.(b) Plot of adiabatic parameter variation along the taper.(c) Normalized modal profiles of E r (left) and |E|/|E 0 | over the area labelled by red dashed line (right) calculated by numerical simulation.(d) Normalized field amplitude |E|/|E 0 | on the exterior surface of graphene taper calculated with adiabatic approximation (red dashed curve) and numerical simulation (green curve).(e) Field amplitude enhancement (|E tip |/|E 0 |) as a function of the tip radius R b .Inset: Plot of SPPs mode lateral extent l spp as a function of tip radius R b .(f) Field amplitude enhancement (|E tip |/|E 0 |) as a function of the taper length L. Without specifications, the taper structure parameters are set as R a = 1 μm, R b = 5 nm, L = 3 μm.

Fig. 3 .
Fig. 3. (a)-(c) Theoretical calculated mode profiles of E r for graphene scrolls, N = 1, 2, 3 respectively.(d) Real part and (e) imaginary part of the effective indices.The red circles denote the modes with the largest Re(n eff ) (s = N).For all the scrolls, R b = 5 nm, R b1 = 4 nm.

Fig. 4 .
Fig. 4. Plasmonic nanofocusing in graphene two-layer taper.(a) Real part and (b) imaginary part of the effective indices for modes s = 1 and s = 2 along the taper.(c) The variation of adiabatic parameter for mode s = 2 along the taper.(d) Normalized filed distributions of E r component (left) and |E|/|E 0 | over the dashed line labelled area (right)calculated by simulation.(e) Normalized field amplitude |E|/|E 0 | on the surface of the inner layer graphene calculated by adiabatic approximation (red dashed curve) and numerical simulation (green curve).(f) Field amplitude enhancement |E tip |/|E 0 | as a function of taper length as ΔR b = 1 nm and 2 nm, respectively.Without specifications, the taper structure parameters are set as R a = 1 μm, R b = 5 nm, R b1 = 4 nm, L = 3 μm.

Fig. 5 .
Fig. 5. Plasmonic nanofocusing in graphene multilayer taper.(a)-(c) Plots of calculated field enhancement |E tip |/|E 0 | as functions of taper length for the 3, 4, 5-layer graphene taper, respectively.(d)-(f) The corresponding field distributions of E r component and |E|/|E 0 | for the 3, 4, 5-layer graphene taper.The taper lengths are labelled with arrows in (a)-(c).(g) Plot of field amplitude enhancement |E tip |/|E 0 | dependence on the graphene layer number N. (h) Plot of inverse of normalized effective mode area (1/ A eff ) at the taper tip (red dots) and propagation loss η (green dots) as functions of the graphene layer number N. R a = 1 μm, R b = 5 nm, R b1 = 4 nm.

Fig. 7 .
Fig. 7. (a) Schematic for the SPPs field superposition between interlayer graphene.(b)Schematic for the field of SPPs propagating and reflecting in the tapered waveguide.

Fig. 8 .
Fig. 8. Influence of bending induced geometric potential on the plasmonic nanofocusing.(a) Plot of graphene chemical potential as a function of curvature radius R. (b) Plot of effective index along the SPPs propagation direction for single layer graphene taper with (Red line) and without (Blue line) considering the influence of geometric potential.(c) Plot of normalized field distribution |E|/|E 0 | for single graphene taper with (Red line) and without (Blue line) considering the influence of geometric potential.

Fig. 9 .
Fig. 9. Calculated reflection coefficients for graphene scroll port, N = 1 and 2 respectively.(a) Plot of reflection coefficients for a single layer graphene scroll as function of scroll radius R. (b) Real part and (c) imaginary part of the reflection coefficients for graphene two-layer scroll as function of scroll sizes.R 1 and R 2 are the radius of the inner layer scroll and the interlayer space between of the two layers, respectively.

Figure 9 (
Figure 9(a) shows the calculated reflection coefficients for single layer graphene scroll as function of scroll radius.Figures 9(b) and 9(c) show the reflection coefficients for mode s = 2 of graphene two-layer scroll mapped as function of R 1 and R 2, where R 1 represents the inner layer radius and R 2 represents the interlayer space.The port reflectance is nearly to be 1 for the both scrolls, which confirms that the SPP is almost completely reflected at the graphene/air interface, in accordance with that in monolayer graphene sheet[38].