## Abstract

By using the effective-index method (EIM) and the finite element methods (FEM), a surface plasmon polariton (SPP) waveguide structured by a dielectric ridge placed on monolayer black phosphorus is proposed and analyzed in the infrared spectral region. It is found that the strong anisotropic dispersion of the black phosphorus (BP) gives rise to direction-dependent waveguide modes in a dielectric-loaded black phosphorus waveguide (DLBPW). The effective mode index, propagation loss and cutoff wavelength of higher order modes are investigated along the armchair (AC) and the zigzag (ZZ) directions of the black phosphorus. Moreover, the propagation characteristics of single-mode are investigated for different widths of the dielectric ridge and different polarization directions of the black phosphorus. Via tuning the carrier density, the electromagnetic field propagation features can be effectively modified. Also, the coupling effect by adding more dielectric bridges can tune the propagation properties. These results will lead to great applications in black phosphorus-based optical integrated devices.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

As the photonic components and circuits become more integrated, they start to suffer from the diffraction limit [1]. Surface plasmon polaritons (SPPs) [2–4], which can confine and control electromagnetic waves at the subwavelength scale, propagate along metal-dielectric interfaces and can overcome the diffraction limit, which enables highly integrated devices in nanophotonics. Noble metals are usually used to support SPPs, and various metal SPPs waveguides have been proposed, such as dielectric-loaded waveguides [5–7], gap plasmonic waveguides [2,8], plasmonic metallic gratings [3] and channel plasmonic waveguides [9,10], which promotes the development of photonic applications. Recently, graphene, a two-dimensional (2D) material, has been demonstrated to propagate the surface plasmon waves which attract enormous attention in the optoelectronic field for infrared and terahertz wavelengths. Compared with metal SPPs, graphene surface plasmons (GSPs) can strongly confine infrared radiation in an atomic monolayer at subwavelength scales due to strong light-matter interaction [11]. Moreover, GSPs could be tuned by chemical doping, and/or the gate voltage [12], which becomes a promising candidate for novel photonic devices. A variety of graphene-based GSP structures have been designed and analyzed, such as graphene nanoribbons [13,14], dielectric loaded graphene waveguides [15–18], and graphene groove [19,20]. Among these structures, dielectric loaded graphene waveguides were investigated intensively owing to their simple fabrication. Black phosphorus (BP), a newly emerging 2D material, exhibits outstanding features in electronics and optics since it has an intrinsic band gap which is different from the zero-gaped graphene [21–23]. BP with a particular atomic structure, a hexagonal lattice with a puckered structure, results in in-plane anisotropic electronic and optical properties [23]. Therefore, BP has been intensely investigated for potential applications in light-matter interaction [24], photovoltaic devices [25] and field effect transistors [26,27]. Recently, SPPs in a BP layer have been studied and analyzed including collective plasmonic excitations [28], localized surface plasmon resonance [29], plasmonic modes in the phosphorene pairs [30] and coupling between graphene and anisotropic black phosphorus surface plasmons [31]. These results reveal that SPPs in the BP layer exhibit stronger field enhancement and higher light confinement, which provide a theoretical basis for BP based optoelectronic devices.

In this paper, we theoretically investigate the SP modes of the dielectric-loaded black phosphorus waveguide (DLBPW) in the infrared spectral region. Compared with dielectric loaded graphene plasmon waveguide, the operation frequency region in DLBPW is broader. Another more intriguing feature is that anisotropy can make the waveguide structure with BP more tunable and applicable. The structure of the DLBPW is easy to fabricate and can be combined with engineered substrate structures which can provide more freedom to manipulate the SPPs. Therefore, it is significant to characterize the properties of SPPs in the DLBPW. The characteristics of the SP modes, including the effective mode index, propagation loss and cutoff wavelength of higher order modes are investigated by using the effective-index method (EIM) and the finite element method (FEM). The anisotropic dispersion of the BP gives rise to orientation dependent modes, which are considerably different from the dielectric loaded metal waveguide and the dielectric loaded graphene waveguide.

## 2. Structure and theoretical analysis

The proposed dielectric-loaded black phosphorus waveguide (DLBPW) structure is shown in Figs. 1(e) and 1(f). A monolayer of BP is placed in between a dielectric ridge and a dielectric substrate. The anisotropic structure of the BP leads to the fact that the waveguide shows different propagation characteristics when varying the angle under which the BP layer is placed under the ridge. The electromagnetic wave propagates along the z-direction (strip direction), and the normal direction of the BP is the y-direction. The angle $\theta $ between the BP ZZ direction and the z-axis can be arbitrary. However, there are two special cases as shown in Fig. 1 where ZZ is parallel to the z-axis $(\theta ={0}^{\circ})$ as shown in Fig. 1(e) and AC parallel to the z-axis $(\theta =\text{9}{0}^{\circ})$ as shown in Fig. 1(f). The width and height of the dielectric ridge are $W$and $h,$ respectively. The operating wavelength is $\lambda .$ For simplicity, in our numerical calculation, constant and non-dispersive relative permittivities of ${\epsilon}_{r1}={\epsilon}_{r2}=2.89$ for the dielectric ridge and substrate are taken [29]. The cladding is air, and the corresponding relative permittivity is ${\epsilon}_{\text{air}}=1.$ The relative permeability in DLBPW are ${\mu}_{r1}={\mu}_{r2}={\mu}_{\text{air}}=1.$

Originating from the puckered honeycomb lattice, the effective mass along the AC direction is smaller than the value in the ZZ-direction. Therefore, the optical conductivity in BP in the AC-direction is larger than that in the ZZ-direction, which results in anisotropic plasmonic properties along different directions. The direction-dependent conductivities of the BP can be described using the semi-classical Drude model [28]:

where $jj$represents the AC- or ZZ-direction.$\hslash $is the reduced Planck's constant, ${D}_{jj}=\pi {e}^{2}n/{m}_{jj}$ is the Drude weight, and $n={10}^{13}{\text{cm}}^{-2}$ is the electron doping density which can be tuned by the gate voltage. $\omega $is the incident optical frequency.$\eta =10\text{meV}$ is the relaxation rate, and ${m}_{jj}$ is the electron effective mass along the AC-direction and ZZ-directions which can be described as ${m}_{AC}={\hslash}^{2}/\left(2{\gamma}^{2}/\Delta +{\eta}_{c}\right)$ and ${m}_{ZZ}={\hslash}^{2}/\left(2{\nu}_{c}\right)$ [28]. For the monolayer BP we have $\gamma =4a/\pi \text{eVm,}$$\Delta \text{=}2\text{eV}$${\eta}_{c}={\hslash}^{2}/0.4{m}_{0},$ ${v}_{c}={\hslash}^{2}/1.4{m}_{0},$ where ${m}_{0}$ is the free electron mass and $a=0.223\text{nm}$ is the scale length of BP.By employing the equivalent relative permittivity with effective width $d,$ the BP’s dielectric function can be described as ${\epsilon}_{jj}=\text{5}\text{.76}+i{\sigma}_{jj}/\left({\epsilon}_{0}\omega d\right)$ which depends on the optical conductivity [29].

The conductivities of the BP with the electric field polarized in the propagating direction (z-axis) can be described as [32]

_{$\theta $}is the angle between the ZZ-direction of BP and the z-axis (direction of electric field). ${\sigma}_{AC}$ and${\sigma}_{ZZ}$ represent the conductivity of the electric field along the AC and ZZ polarization direction. When $\theta $ is at any angle, the conductivity is the total contribution in both directions.

In the presence of the complex optical conductivity, the BP can support and propagate transverse magnetic (TM) plasmonic waves which show similar properties to the ones in graphene [30]. Assuming the fields of the incident light are propagating along the planar BP sheet sandwiched between two dielectrics with the relative permittivity ${\epsilon}_{r1}$ and ${\epsilon}_{r2}$ possess the form of $\mathrm{exp}\left(i\beta z-i\omega t\right)$ where the electromagnetic wave propagates along the z-direction. Here, $\beta $ is the propagating wavevector in the three-layer waveguide structure. For the TM mode, the magnetic field in the x-axis direction can be described as:

Different dielectric coating layer (ridge) can be used to tune the SPP propagating characteristics. Similar to dielectric-loaded graphene waveguides [15], the dielectric strip results in a higher refractive index for the SP modes on the BP-dielectric interface compared to the BP-air interface, gives rise to SPPs bound by the dielectric strip. In a DLBPW structure, the dielectric strip serves as the core of the sandwich, and there are two special cases as shown in Figs. 1(e) and 1(f). The effective-index method (EIM) can be used to calculate the propagation parameters.

In the first step of the EIM, the refractive index ${N}_{\text{core}}$ of the SP mode in the BP sandwiched between two dielectrics can be obtained from Eq. (6). Likewise, the effective mode index ${N}_{\text{clad}}$ of the SP mode can be obtained similarly.

In the second step of the EIM, the eigenequation of the equivalent dielectric planar waveguide for the m-th order guided transverse electric (TE) mode is given as:

where${\xi}_{1}\text{=}\sqrt{{N}_{\text{core}}^{2}{k}_{0}^{2}-{\beta}^{2}}$,${\xi}_{2}\text{=}\sqrt{{\beta}^{2}-{N}_{\text{clad}}^{\text{2}}{k}_{0}^{2}}$.In the DLBPW, the real part of the effective refractive index is${N}_{\text{eff}}=\mathrm{Re}\left(\beta /{k}_{0}\right).$ The imaginary part is $\mathrm{Im}\left(\beta /{k}_{0}\right)$.

The cutoff condition of the guided modes is ${\xi}_{2}=0,$ and the cutoff wavelength of m-th order guided mode is:

The influence of the height of the dielectric strip has not been taken into account in EIM, because of the SP mode is highly confined at the interface of the BP. In our numerical calculation, the height of the dielectric strip is 180 nm, which indicates that the electromagnetic wave field surpasses the top of the dielectric strip could be neglected in an analytical calculation.

## 3. Results and discussion

The BP optical conductivity depends on the electric field polarized in a different direction. Figures 1(a) and 1(c) show the effective mode indices of the SP modes as a function of the wavelength for the angles $\theta ={0}^{\circ}$(i.e., ${\sigma}_{ZZ})$ and $\theta =\text{9}{0}^{\circ}$ (i.e., ${\sigma}_{AC})$ by the EIM and the FEM. The insets of Figs. 1(a) and 1(c) show the corresponding amplitudes of ${E}_{y}$ at different wavelengths. It can be seen that, at the wavelength of $40\text{\mu m,}$ the fundamental mode along the ZZ-direction $(\theta ={0}^{\circ})$ is ${N}_{\text{eff}}=165,$ while along the AC-direction $(\theta =\text{9}{0}^{\circ})$ ${N}_{\text{eff}}=65,$ which indicates the plasmonic waves with the electric field polarized in the ZZ-direction possesses larger effective refractive indices than that in the AC-direction. Likewise, the propagation loss along the ZZ-direction is more significant than that along the AC-direction because the conductivity along the ZZ-direction is smaller than the AC-direction.

The fundamental mode (m = 0) has a higher real part of the effective mode index than the higher order modes, which indicates a shorter wavelength. With the increase of the wavelength, the effective mode index decreases monotonically, and the electric field becomes more substantial (see inset of Fig. 1(a)). The propagation loss of higher order mode decreases sharply with increasing wavelength. Together with the real part of the effective mode index, it indicates the miserable confinement of the electromagnetic field.

The EIM results show a good agreement with the FEM simulation results for the real part. However, when the wavelength is close to the cutoff wavelength, the analytical results show slight deviations from the FEM results. For that, there are two main reasons. First, in the EIM analysis, the BP is characterized by interface conductivity with zero thickness. In the FEM simulation, the BP is treated as a dielectric with finite thickness. Second, the EIM makes an idealized assumption for the continuity condition of the electromagnetic wave at the interface and neglects the change of electromagnetic field in the corner regions. When the electromagnetic field is confined in a core, the EIM results will perfectly match with the FEM simulation results at high frequencies. In the low-frequency region, the electromagnetic field cannot be perfectly constrained in the strip. Therefore, the difference between the two methods gradually shows up. This means, in areas where electromagnetic waves cannot be constrained entirely, the accuracy of the EIM method is reduced. In contrast, the results of the FEM method are more accurate.

It can be seen from Figs. 1(b) and 1(d) that the FEM results deviate from the EIM results more in the shorter wavelength which can be understood from the scattering factor $(\eta )$ in the BP optical conductivity. When ignoring the scattering factor $(\eta =0),$ only the imaginary part of the BP intraband optical conductivity is obtained which leads ${N}_{\text{core}}$ and ${N}_{\text{clad}}$ being real. The present structure is similar as the traditional planar waveguide. The waveguide modes can be obtained with real effective refractive index $({N}_{\text{eff}}).$ However, when $\eta \ne 0,$ the refractive indices and the propagating wave vector become complex. And from our numerical calculation, $\mathrm{Im}{N}_{\text{eff}}$ isn’t confined between $\mathrm{Im}{N}_{\text{core}}$ and $\mathrm{Im}{N}_{\text{clad}},$ and its much smaller than the real value $\mathrm{Re}{N}_{\text{eff}}$ which is consistent with the case in Ref [33]. $\mathrm{Im}{N}_{\text{eff}}$ exhibits complex dependence on the total BP optical contributions both from the real and imaginary parts.

By solving Eq. (8), the fundamental mode and the multi-modes operation region as a function of the width of the dielectric strip can be calculated. The results are shown in Fig. 2. At a fixed wavelength, the numbers of guided modes increase with increasing the width of the dielectric strip. The increase indicates that the higher order modes can propagate in the DLBPW with a large width. The number of guided modes decrease with increasing the wavelength at a fixed width $W$. Especially, when the width is small and the wavelength is large enough, the DLBPW can only propagate the fundamental mode. Also, the DLBPW shows a larger fundamental mode operation region for $\theta =\text{9}{0}^{\circ}$ in comparison to the case of $\theta ={0}^{\circ}$, which demonstrates that different placement of the BP on dielectric media can be used to tune the mode operation regions.

As already mentioned, on the BP layer surface, the electric field can be polarized along the ZZ-, AC- or any direction which leads to the modified surface conductivity. In Fig. 3, at a fixed wavelength, the effective fundamental mode indices as a function of this angle $\theta $ decrease monotonically with increasing this angle. It is attributed to the BP possessing a larger conductivity along the AC-direction $(\theta =\text{9}{0}^{\circ})$ than that along the ZZ-direction $(\theta ={0}^{\circ}).$ The strong anisotropic dispersion of the BP gives rise to the direction-dependent modes of the DLBPW, which is expected to realize in designing novel devices. At a fixed angle, effective mode indices increase as the wavelength decreases at the expense of higher propagation loss. Therefore, it is essential to consider this in practical applications.

Figure 4 shows the effective refractive indices of the fundamental mode in the DLBPW as a function of the height and various widths of the dielectric ridge. ${N}_{\text{eff}}$ and $\mathrm{Im}\left(\beta /{k}_{0}\right)$ increase with increasing width. As shown in Fig. 4(a), it has a significant impact on ${N}_{\text{eff}}$when the height is small. When the height is significantly more than $80\text{nm}$, the effect of the height on ${N}_{\text{eff}}$ is no longer visible. When the height is small, the electromagnetic field overflows, and there is a significant edge effect. The edge effect limits the equivalent refractive index. When the height reaches $80\text{nm}$, the electromagnetic field is minimal, and therefore the increase of the height has no longer a significant effect. When the height is 0, surface plasmons occur, and when the height is gradually increased, the waveguide appears, the intermediate process is mixed. As shown in Fig. 4(b), the electromagnetic wave loss shows a peak with the height.

In our opinion that behavior can be useful for the DLBPW and can provide valuable applications for integrated optoelectronic devices.

Figure 5 shows the real part of the effective refractive indices of the fundamental mode in the DLBPW with a width of $200\text{nm}\text{.}$ By varying the electron doping, the effective refractive indices of the fundamental mode can be effectively tuned. Figure 5 shows that ${N}_{\text{eff}}$ decrease with the increase of the electron doping $n.$ For the wavelength of $30\text{\mu m}$ at the electron doping being $n={10}^{13}{\text{cm}}^{-2},$ the real part of the effective refractive index is more extensive and larger than ${N}_{\text{eff}}=\text{200,}$ which indicates that the SP mode in the DLBPW is severely confined and is smaller than the GSPs mode in the dielectric-loaded graphene plasmon waveguide (DLGPW) [15]. It can be beneficial to highly integrated optoelectronic devices. The real part of the effective refractive indices is much more sensitive to the electron doping at $\theta ={0}^{\circ}$ (Fig. 5(a)). The main reason is that the dependence of the conductivity along ZZ-/AC-direction on the electron doping is different, see Eq. (2).

When two waveguides are integrated on a chip in a parallel alignment, a coupling occurs due to the mode field overlap. Figure 6 shows the effective refractive indices of the fundamental mode in the DLBPW coupling. ${N}_{\text{eff}}$ decreases with the increase of $s$ and approaches the value of a single DLBPW as the distance increases. The coupling can lead to stronger constraints. The imaginary part increases with the increase of $s,$ and gradually approaches the imaginary part of a single DLBPW, as can see in Figs. 6(b) and 6(d). Such behavior indicates that the loss increases as the coupling decreases. This finding means the coupling enhances the constraint and reduces losses. In summary, the coupling can also be used as a means of adjusting the performance of the waveguide.

## 4. Conclusion

Within this paper, we propose the concept of the dielectric-loaded black phosphorus waveguide (DLBPW) in the infrared spectral region. The SP modes of the DLBPW are analyzed using the effective-index method (EIM) and the finite element method (FEM). The effective refractive index, propagation loss, the cutoff wavelength of higher order modes is calculated. The inserted monolayer BP (black phosphorene) is polarized in the armchair (AC) and/or zigzag (ZZ) directions. The characteristic of the single-mode is investigated for different width of the dielectric ridge and different polarized directions of black phosphorene. Via tuning the electron doping, the propagation properties of the fundamental mode in DLBPW can be tuned flexibly. From one dielectric ridge to two dielectric ridges, the coupling between them can also effectively modify the propagation features. Furthermore, the DLBPW provides direction-dependent propagation properties, which will lead to novel optical functional devices.

## Funding

National Natural Science Foundation of China (11547030 and 11574319), the Center of Science and Technology of Hefei Academy of Science (2016FXZY002), and the Department of Science and Technology of Yunnan Province (2016FC001).

## Acknowledgments

We are particularly grateful to Dr. Wieser Robert for his helpful discussion in the process of our writing this article.

## References

**1. **D. K. Gramotnev and S. I. Bozhevolnyi, “Plasmonics beyond the diffraction limit,” Nat. Photonics **4**(2), 83–91 (2010). [CrossRef]

**2. **W. L. Barnes, A. Dereux, and T. W. Ebbesen, “Surface plasmon subwavelength optics,” Nature **424**(6950), 824–830 (2003). [CrossRef] [PubMed]

**3. **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). [CrossRef] [PubMed]

**4. **M. L. Brongersma and V. M. Shalaev, “Applied physics. The case for plasmonics,” Science **328**(5977), 440–441 (2010). [CrossRef] [PubMed]

**5. **T. Holmgaard and S. I. Bozhevolnyi, “Theoretical analysis of dielectric-loaded surface plasmon-polariton waveguides,” Phys. Rev. B Condens. Matter Mater. Phys. **75**(42), 245405 (2007). [CrossRef]

**6. **K. Hassan, J. C. Weeber, L. Markey, A. Dereux, A. Pitilakis, O. Tsilipakos, and E. Kriezis, “Thermo-optic plasmo-photonic mode interference switches based on dielectric loaded waveguides,” Appl. Phys. Lett. **99**(24), 241110 (2011). [CrossRef]

**7. **H. S. Chu, E. P. Li, P. Bai, and R. Hegde, “Optical performance of single-mode hybrid dielectric-loaded plasmonic waveguide-based components,” Appl. Phys. Lett. **96**(22), 221103 (2010). [CrossRef]

**8. **A. Ishikawa, S. Zhang, D. A. Genov, G. Bartal, and X. Zhang, “Deep subwavelength terahertz waveguides using gap magnetic plasmon,” Phys. Rev. Lett. **102**(4), 043904 (2009). [CrossRef] [PubMed]

**9. **S. I. Bozhevolnyi, V. S. Volkov, E. Devaux, J. Y. Laluet, and T. W. Ebbesen, “Channel plasmon subwavelength waveguide components including interferometers and ring resonators,” Nature **440**(7083), 508–511 (2006). [CrossRef] [PubMed]

**10. **J. Dintinger and O. J. F. Martin, “Channel and wedge plasmon modes of metallic V-grooves with finite metal thickness,” Opt. Express **17**(4), 2364–2374 (2009). [CrossRef] [PubMed]

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

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

**13. **J. Christensen, A. Manjavacas, S. Thongrattanasiri, F. H. 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]

**14. **E. Forati and G. W. Hanson, “Surface plasmon polaritons on soft-boundary graphene nanoribbons and their application in switching/demultiplexing,” Appl. Phys. Lett. **103**(13), 133104 (2013). [CrossRef]

**15. **W. Xu, Z. H. Zhu, K. Liu, J. F. Zhang, X. D. Yuan, Q. S. Lu, and S. Q. Qin, “Dielectric loaded graphene plasmon waveguide,” Opt. Express **23**(4), 5147–5153 (2015). [CrossRef] [PubMed]

**16. **X. Zhou, T. Zhang, L. Chen, W. Hong, and X. Li, “A graphene-based hybrid plasmonic waveguide with ultra-deep subwavelength confinement,” J. Lightwave Technol. **32**(21), 3597–3601 (2014).

**17. **J. Gosciniak and D. T. Tan, “Graphene-based waveguide integrated dielectric-loaded plasmonic electro-absorption modulators,” Nanotechnology **24**(18), 185202 (2013). [CrossRef] [PubMed]

**18. **W. Xu, Z. H. Zhu, K. Liu, J. F. Zhang, X. D. Yuan, Q. S. Lu, and S. Q. Qin, “Toward integrated electrically controllable directional coupling based on dielectric loaded graphene plasmonic waveguide,” Opt. Lett. **40**(7), 1603–1606 (2015). [CrossRef] [PubMed]

**19. **P. Liu, X. Zhang, Z. Ma, W. Cai, L. Wang, and J. Xu, “Surface plasmon modes in graphene wedge and groove waveguides,” Opt. Express **21**(26), 32432–32440 (2013). [CrossRef] [PubMed]

**20. **P. Wan, C. Yang, and Z. Liu, “Channel hybrid plasmonic modes in dielectric-loaded graphene groove waveguides,” Opt. Commun. **420**, 72–77 (2018). [CrossRef]

**21. **S. Balendhran, S. Walia, H. Nili, S. Sriram, and M. Bhaskaran, “Elemental analogues of graphene: silicene, germanene, stanene, and phosphorene,” Small **11**(6), 640–652 (2015). [CrossRef] [PubMed]

**22. **T. Low, A. Rodin, A. Carvalho, Y. Jiang, H. Wang, F. Xia, and A. C. Neto, “Tunable optical properties of multilayer black phosphorus thin films,” Phys. Rev. B Condens. Matter Mater. Phys. **90**(7), 075434 (2014). [CrossRef]

**23. **V. Tran, R. Soklaski, Y. Liang, and L. Yang, “Layer-controlled band gap and anisotropic excitons in few-layer black phosphorus,” Phys. Rev. B Condens. Matter Mater. Phys. **89**(23), 235319 (2014). [CrossRef]

**24. **J. Lu, J. Yang, A. Carvalho, H. Liu, Y. Lu, and C. H. Sow, “Light–matter interactions in phosphorene,” Acc. Chem. Res. **49**(9), 1806–1815 (2016). [CrossRef] [PubMed]

**25. **M. Buscema, D. J. Groenendijk, G. A. Steele, H. S. van der Zant, and A. Castellanos-Gomez, “Photovoltaic effect in few-layer black phosphorus PN junctions defined by local electrostatic gating,” Nat. Commun. **5**(1), 4651 (2014). [CrossRef] [PubMed]

**26. **M. Buscema, D. J. Groenendijk, S. I. Blanter, G. A. Steele, H. S. van der Zant, and A. Castellanos-Gomez, “Fast and broadband photoresponse of few-layer black phosphorus field-effect transistors,” Nano Lett. **14**(6), 3347–3352 (2014). [CrossRef] [PubMed]

**27. **Y. Du, H. Liu, Y. Deng, and P. D. Ye, “Device perspective for black phosphorus field-effect transistors: contact resistance, ambipolar behavior, and scaling,” ACS Nano **8**(10), 10035–10042 (2014). [CrossRef] [PubMed]

**28. **T. Low, R. Roldán, H. Wang, F. Xia, P. Avouris, L. M. Moreno, and F. Guinea, “Plasmons and screening in monolayer and multilayer black phosphorus,” Phys. Rev. Lett. **113**(10), 106802 (2014). [CrossRef] [PubMed]

**29. **Z. Liu and K. Aydin, “Localized surface plasmons in nanostructured monolayer black phosphorus,” Nano Lett. **16**(6), 3457–3462 (2016). [CrossRef] [PubMed]

**30. **H. Lu, Y. Gong, D. Mao, X. Gan, and J. Zhao, “Strong plasmonic confinement and optical force in phosphorene pairs,” Opt. Express **25**(5), 5255–5263 (2017). [CrossRef] [PubMed]

**31. **J. Nong, W. Wei, W. Wang, G. Lan, Z. Shang, J. Yi, and L. Tang, “Strong coherent coupling between graphene surface plasmons and anisotropic black phosphorus localized surface plasmons,” Opt. Express **26**(2), 1633–1644 (2018). [CrossRef] [PubMed]

**32. **F. Xia, H. Wang, and Y. Jia, “Rediscovering black phosphorus as an anisotropic layered material for optoelectronics and electronics,” Nat. Commun. **5**(1), 4458 (2014). [CrossRef] [PubMed]

**33. **Y. Q. Liu and P. K. Liu, “Excitation of surface plasmon polaritons by electron beam with graphene ribbon arrays,” J. Appl. Phys. **121**(11), 113104 (2017). [CrossRef]