Energy loss spectrum and surface modes of two-dimensional black phosphorus

The structural features and the electron energy loss spectrum of black phosphorus (BP) have been experimentally analyzed and they are discussed based on a theoretical calculation. The low-energy loss spectra of typical samples reveal that the emerging high-mobility two-dimensional material BP often exhibits both bulk and surface plasmon modes. The surface modes of BP are strongly thickness dependent. Electrodynamic analysis indicates that the Fuchs–Kliewer-like surface plasmon modes consist of two branches with different charge symmetry: the upper side and lower side have the same charge polarity as the lower branch and the opposite charge polarity to the upper branch. This study provides fundamental insight into the characteristic nature of BP plasmonics.


Introduction
Surface plasmon is a type of surface electromagnetic mode propagation along the interface between the metal like medium and the dielectrics. The amplitude of the electric field decays exponentially away from the interface. Surface plasmonics of particular interest because they offer an efficient method for light manipulation [1]. In the past decade, as an important subject of surface plasmonic investigations, graphene plasmonics has provided a variety of gate-controlled optoelectronic applications, such as light harvesting and optical sensing [2,3]. In addition to the unique electronic structure of graphene, graphene plasmonic devices can confine light in atomic thickness, resulting in dramatically enhanced local fields owing to graphene's two-dimensional (2D) nature. However, replaceable 2D plasmonic materials are still lacking. Exploration of new 2D plasmonic building blocks in novel material systems is highly desirable to advance the field of plasmonics. As a result, a number of 2D materials in addition to graphene have been successfully prepared, including hexagonal boron nitride, transition metal dichalcogenides, and black phosphorus (BP) [4][5][6]. Among these emerging 2D systems, BP is one of the most promising candidates because of its remarkable electronic and photonic properties, such as its high carrier mobility (5200 cm 2 V −1 s −1 at room temperature), thickness-dependent direct bandgap, and gate-tunable optoelectronic properties [7][8][9][10]. A variety of experimental studies have been performed to investigate the structural and physical properties of BP by scanning tunneling microscopy [11], atomic force microscopy, Raman spectroscopy [12], and the electron loss spectroscopy [13,14]. Although theoretical and experimental studies have demonstrated the potential of BP as a promising plasmonic platform [15][16][17][18][19][20][21], systematic analysis of the thickness-dependent surface modes in the low energy loss spectrum is still lacking. In the present work, we reveal that few-layer BP exhibits both bulk and surface plasmon modes in the visible and ultraviolet regions by electron energy loss spectroscopy (EELS) in transmission electron microscopy (TEM). The surface modes of BP are strongly thickness dependent. We performed an ab initio calculation of the dielectric function of BP. Based on the anisotropic dielectric function and electrostatic approximation, we also calculated the contribution of the Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. thickness-dependent anisotropic surface modes to the low energy loss spectrum. The theoretical results are in good qualitative agreement with the experimental results.

Experimental
Crystalline BP is a typical sensitive material when exposed to ambient O 2 and H 2 O, and it shows visible structural degradation [22]. Although the mechanism is still not completely understood, the presence of electrons/light has been shown to either initiate or accelerate this degradation [14]. Thus, the pristine BP samples in this study were handled with minimal ambient light exposure, and prior structural/compositional analyses were performed to ensure that the BP flakes selected for study were not significantly altered. In the experimental measurements, BP crystalline plates with thicknesses of 1-10 layers in distilled water solvent were dropped onto a 400 mesh copper grid on which carbon nanotubes were grown. After the distilled water vaporized, the BP nanoplates were attached to the nanotubes. The EELS measurements were performed with a JEOL-2100F electron microscope equipped with a Gatan spectrometer operating at a voltage of 200 kV. A collection semi-angle of more than 100 mrad was used. The TEM images were recorded with a JEM-ARM200 TEM (JEOL Inc.). High-resolution real space microscopy and EELS were performed for free-standing BP supported by carbon nanotubes on the edges. The advantage of using the carbon nanotubes as a scaffold for the flakes is that it can avoid introducing the background from the substrate in the EELS. Also, the surface electromagnetic modes are sensitive to the dielectric functions of the mediums above and below the interface along which it propagates. So, the spectroscopic feature and plasmonic behavior of the BP supported on the carbon nanotubes are different from the ones of the BP deposited on a flat substrate.
A TEM image of the sample at low magnification is shown in figure 1(a). It clearly shows that the BP crystalline plates are supported by the carbon nanotubes. A high-resolution real space image and the EELS spectrum of the BP plates are shown in figures 1(c) and (d), respectively. The crystal structure and lattice symmetry were investigated by TEM observations and electron diffraction along the main axis directions. All of the obtained experimental data can be well indexed to the Cmca orthorhombic unit cell with lattice parameters of a=3.313 Å, b=4.376 Å, and c=10.478 Å. The structural image in figure 1(c) is for an exfoliated BP flake viewed along the [001] crystallographic direction, where P atomic columns can be directly observed. The schematic structural model is shown in figure 1(b), which is consistent with previously reported structures as AA stacking [13,23]. The core energy loss spectrum in figure 1(d) clearly shows that the phosphorus L 23 peak is located at 132 eV. Its shape is consistent with what has been observed with the pristine BP [14].
The low-energy loss spectrum recorded at 200k magnification with an exposure time of 1 s in the image mode is shown in figure 2. The thickness of the crystalline samples was analyzed as follows. First, the zero loss peaks were isolated and the elastic scattering counts I 0 were recorded. The total counts I t were also recorded. The mean free path was estimated by the equation [24] F E E E E 106 ln 2 , where F is the relativistic factor, E 0 is the incident electron energy, b is the semi-collection angle, and E m is the mean energy loss estimated by the empirical formula E Z 7.6 ,

Data analysis
The thickness-dependent energy loss spectrum is dominated by the bulk loss of BP at around 19 eV, and the energy loss associated with excitation of the surface modes on the upper and lower surface of the crystal plate appears as a broad plateau. With increasing thickness, both features in the EELS spectrum become stronger. The plasmons around 20 eV in BP has been reported recently by Nicotra et al [16] for both the zigzag and armchair directions. The weak shoulder observed in the EELS spectrum around 10 eV is identified as a single-particle transition from the 2 G + band, derived from p z orbitals, whose wave function has alternating signs between nearest neighbors. To extract the contributions of the surface modes, we removed the zero loss peaks by the reflected tail fit function. The bulk contributions in the experimental EELS spectra were then calculated using the jellium model. The bulk plasmon frequencies E p and full widths at half-maximum of the plasmon peaks E p D were fitted to the experimental data using the following equation: where E is the energy. The fitted bulk losses were then removed and the surface loss contributions remain, which are shown in figure 3.
We performed an ab initio calculation of the electronic structure of BP with the Vienna ab initio simulation package. Two layers of phosphorus atoms stack as the AA type. A 4.37 Å×3.31 Å×10.47 Å unit cell with eight atoms was constructed. The projector augmented-wave approach and the generalized gradient approximation exchange-correlation function were used in the calculation. We used 500 eV as the energy cutoff of the wave functions and a 21×21×7 K-point mesh. The obtained anisotropic dielectric functions of BP in the X, Y, and Z directions shown in figure 4 were used as the input parameters in the following calculation of the EELS spectrum.
We will now discuss the theoretical calculation and analysis of the EELS spectrum of BP. Because of the difficulty of an accurate calculation of the EELS spectrum for the biaxial material, the analytical analysis in the present study is confined to the calculation of the EELS spectrum of the uniaxial-type thin film. The EELS data of the uniaxial crystals have the same in-plane dielectric functions as the dielectric functions along the x-and y-axis directions, whereas the out-of-plane dielectric function (i.e. the dielectric function of the z-axis direction) is different. The biaxial nature of BP was considered and EELS calculation was separately performed for the XZ and YZ planes. The analytical analysis was performed in the manner of nonretarded approximation. Similar analysis of the electron energy loss probability of uniaxial thin films has been performed for graphite [25,26] and hexagonal boron nitride [27], in which the theoretical results are in good agreement with the experimental measurements. We started with separate calculations of the guided modes, bulk loss, and Begrenzungseffekt. The total electron energy loss probability is the sum of the three contributions.
The typical results for a 5 nm thick hypothetical uniaxial thin film with the in-plane dielectric function equal to the calculated dielectric functions for the x-and y-axis directions and the out-of-plane dielectric function equal to the calculated z direction dielectric function is shown in figure 5. A bulk plasmon peak at 19 eV is the prominent feature in both situations. The plateaus below 19 eV in the loss functions are contributed by the surface guided modes confined by the upper and lower boundaries of the thin flakes. The energy transfer from the bulk mode to the guided modes is described by the Begrenzungseffekt. The total loss functions are the sum of the loss functions obtained from the bulk modes, guided modes, and Begrenzungseffekt, which have similar profiles as the experimental data.
To further understand the structural features of the EELS spectrum, we plot the dispersion of the guided modes and bulk modes in log scale in figure 6. The dispersions of the relevant modes calculated with the X and Z direction dielectric functions are shown in figures 6(a) and (c), and the counterparts calculated with the Y and Z direction dielectric functions are shown in figures 6(b) and (d). Figures 6(a) and (b) are the bulk mode dispersion at 19 eV. Figures 6(c) and (d) are the guided mode dispersion below 19 eV. The guided modes split into two branches in both sets of plots. The lower branches below 14 eV in both sets of plots correspond to the lower branch of the surface plasmon pointed out by Kliewer and Fuchs [28,29], while the upper branches between 14 and 19 eV are the upper branch of the surface plasmon. The calculated dielectric functions of the X, Y, and Z  Another important feature is that the surface loss function has a local maximum at about 0.7 eV. This is consistent with the fact that small but obvious peaks appear at a frequency below 2 eV in the experimental data.  We will now discuss the fundamental features of the total loss functions and guided mode loss functions for samples with different thicknesses ( figure 7). The bulk loss monotonically increases with increasing thickness, as expected. However, the guided mode behavior is more complicated. The main peak of the surface loss moves to higher energy with increasing thickness, which is consistent with the experimental data for the surface loss part of the spectrum (figure 3). The thickness dependence of the data shown in figure 3 mainly comes from two aspects. First, the shape change of the curves with thickness of 7.76, 10.58, 15.79 nm mainly comes from the thickness dependence of the surface mode dispersion as point out by [27]. Second, what we measured are actually the stacks of BP flakes, especially for the thicker regions. So the surface mode contributions from the upper and lower sides of the slabs add together. It also induces thickness dependence to our data.

Conclusions
BP has been investigated by high-resolution electron microscopy and EELS. The orthorhombic crystal structure and crystal lattice constants were determined by real space measurement, as well as diffraction analysis. We also calculated the anisotropic dielectric function of BP by an ab initio method. Quasistatic analysis was performed and the EELS spectrum was analyzed within the frame of classical electrodynamics. The biaxial nature of BP was considered and EELS calculation was separately performed with the dielectric functions of BP in X, Z direction and Y, Z direction. The contributions of the bulk and surface plasmons were separately measured and calculated. The Fuchs-Kliewer-like surface plasmon modes within the BP slab are observed and the data are well supported by analytical analysis. The fast electron can act as a pulsed white light source with evanescent components and efficiently stimulate the surface plasmons along the surface of BP. The lower branch of Fuch Kliewer modes with smaller energies than 14 eV has the same charge polarity on the upper surface and lower surface of the slab while the upper branch of Fuch Kliewer modes with larger energies than 14 eV has the opposite charge polarity on the upper surface and lower surface of the slab. BP shows potential as a future 2D optoelectronic material and it has possible applications in plasmonic devices.