Lifetime and nonlinearity of modulated surface plasmon for black phosphorus sensing application

Black phosphorus surface plasmon (BPSP) is a new promising candidate material for electromagnetic field confinement at the subwavelength scale. Here, we theoretically investigated the light confinement, second-order nonlinearity and lifetimes of tunable surface plasmons in nanostructured black phosphorus nanoflakes with superstrates. The grating structure can enhance the local optical field of the fundamental wave (FW) and second harmonic wave (SHW) due to the surface plasmon resonance. Based on the coupled mode theory (CMT), a theoretical model for the nanostructured black phosphorus was established to study the spectrum features of FW. The lifetimes of the plasmonic resonant modes were investigated with the finite difference time domain (FDTD) simulations and CMT. Since the permittivity of black phosphorus depends on its Fermi energy and electron scattering rate, the lifetimes of plasmonic absorption modes are tunable with both the Fermi energy and scattering rate. The intensity, wavelengths and spectral width of BPSP resonance modes and their lifetimes can be precisely controlled with the Fermi energy, scattering rate, side length and refractive index of the superstrate. The sensitivity is described by varying the refractive index of the superstrate such as an aqueous solution. We have introduced a secondorder nonlinear source to investigate the SHW of nanostructured black phosphorus. This paper presents the corner/edge energy distribution and the tunable lifetime of BPSP as well as their unprecedented capability of photon manipulation for second-order nonlinearity within the deep subwavelength scale. The configuration and method are useful for research of the absorption, lifetime of light and nonlinear optical processes in black phosphorus-based optoelectronic devices, especially the modulation and sensing applications.


Introduction
6][7] Black phosphorus-based saturable absorbers are used to achieve pulse emission of lasers, Q-switching, and mode locking. 8The 2D black phosphorus, with puckered hexagonal honeycomb structure having ridges of atoms, has strong in-plane anisotropic electrical and optical properties. 9,10lack phosphorus has a thickness-dependent direct bandgap ranging up to ∼2 eV, which could possess broadband photo-response from the mid-IR to far-IR and even the visible band. 11,12The direct bandgap and plasmon of black phosphorus can be tuned with electric voltage bias, even for more than 20 layers.The black phosphorus plasmonic properties and unprecedented capability of wide-band photon manipulation have been presented within the deep subwavelength scale. 13Compared with traditional bulk plasmonic materials, it is possible to integrate devices on the atomic scale using few-layer black phosphorus due to highly localized electric fields. 14Highly confined and electrically tunable graphene plasmons, such as electromagnetically induced transparency, and Fano resonance, have been studied below the classical diffraction limit. 15The dipole-like interactions have been observed in 2D graphene materials due to the localized surface plasmon (LSP) resonance effect (the collective oscillations of charged carriers). 16,175][26][27][28][29][30][31][32][33][34] Experiments and theories have shown that linear transmission through arrays is sensitive to many factors such as array periodicity, hole shape and size, incidence angle, and film thickness.For instance, the experimental observations of the SHG signal have also been reported. 25Methods of enhancing the SHG were widely studied by double plasmonic resonance, Fano resonance, and even nanocup structure. 31They measured the SHG from a periodic array using femtosecond laser pulses and observed the SHG signal intensity.Particularly they measured the propagation time duration of the pulses through the array to clarify the slow light effect.Modulated black phosphorus surface plasmon offers a promising way to realize chip-scale integrated photonic circuits because of their ability to confine and control electromagnetic waves at the subwavelength scale with strong confinement of surface plasmon.It was found that some challenges such as serious scattering loss, which is caused by the substrate, boundary and absorbents, limit the mode lifetime, plasmon propagation length and spectral range.For the black phosphorus, it is imperative to investigate the localized surface plasma properties and their lifetime.The findings provide a possibility to realize black phosphorus-based optoelectronic devices that are highly photoresponsive, wavelength selectable and lifetime tunable.The fabrication of layered BP, especially phosphorene, may pave the way towards the realization of high-performance nano-devices, photodetectors and other photothermal properties. 35The peak position and bandwidth of the phosphorene absorption spectra are tunable via the surrounding dielectric environment of the periodic nanoribbons. 36n this paper, we numerically analyze the confined surface plasmon of the fundamental wave and the second harmonic wave in the nanostructured black phosphorus nanoflake, which forms a grating to enhance the light-matter interaction.Coupling of light from the normal direction with nanostructured black phosphorus requires an arbitrary grating coupler providing lateral photon momentum.We present a CMT analysis based on the cavity theory that describes the linear resonant modes of the cavity systems. 37,47The linear absorption spectra for various Fermi energies, the electron scattering rate, side length of black phosphorus and the refractive indexes of superstrates with the FDTD method agree well with the theoretical description.Lifetime is important for the research of intrinsic loss for the m th mode and interaction between each mode and light field.The lifetimes of plasmonic absorption modes are tunable with the Fermi energy, the scattering rate, side length of black phosphorus and refractive index of the superstrate.The intensity, wavelengths and spectral width of BPSP resonance modes can be precisely controlled.Moreover, the sensitivity is demonstrated by varying the refractive index n 1 of the superstrate such as an aqueous solution.The second-order nonlinearity in nanostructured black phosphorus is achieved due to the symmetry breaking.
With the introduction of the second-order nonlinear source, we investigate the second-order nonlinearity such as second harmonic generation (SHG), sum frequency generation (SFG) and difference frequency generation (DFG).[40][41][42] 2. The coupled BPSP system and theoretical analysis We consider a nanostructured black phosphorus nanoflake, formed by the black phosphorus square array, which is composed of the elementary cell with period P = 200 nm in the xyplane.The black phosphorus square array is placed on the SiO 2 substrate and covered with a superstrate such as an aqueous solution in Fig. 1(a).The black phosphorus nanoflake has the cuboid shape (length × width × thickness) with L × L × d.The parameters of the black phosphorus flake are set as follows: side length L = 150 nm and thickness d = 2 nm.The refractive index of the superstrate is n 1 in the model.The permittivity of a 2D black phosphorus layer can be characterized as a diagonal tensor.
All simulations are carried out for photon energies of up to 0.7 eV, 12,17,18,36 the bandgap of BP.Within this energy window, the Drude model is sufficient to model the conductivity of the multilayer BP system.The conductivity of few-layer black phosphorus is given as 7,[12][13][14][15]36 where j denotes the x-or y-direction and D j = πe 2 n sj /m j , is the Drude weight along the j-axis; e is the electron charge; m 1 and m 2 denote the in-plane electron effective masses near the Γ point within the Hamiltonian model, which can be written as m 1 = ħ 2 /(2γ 2 /Δ + η c ) and m 2 = ħ 2 /v c .Here, ħ is the reduced Planck constant, η e (eV) is the parameter describing the scattering rate in black phosphorus; n s1 and n s2 denote the elements of electronic density in black phosphorus along the x-and y-directions, respectively; ω is the frequency of the incident light.
We can introduce a phosphorene layer with volumetric anisotropic permittivity converting the 2D surface conductivity into a 3D conductivity using the relation σ 2D = d σ 3D . 36The equivalent relative permittivity diagonal tensor ε D (ε 11 , ε 22 , ε 33 ) is defined to describe the black phosphorus with thickness d.The equivalent relative permittivity tensor for the in-plane component of black phosphorus is given by the following: where the z component is ε 33 = ε r .ε 0 = 8.854 × 10 −12 F m −1 is the vacuum permittivity.ε r is the relative permittivity.For 2D black phosphorus, ε r = 5.76. 7From the classical Drude model, the value of the collision frequency γ bp of the Drude model for black phosphorus in the mid-IR and far-IR wavelength regimes can be approximately given as γ bp = η e /ħ.For nanostructured black phosphorus, the parameters of the conduction band can be set as follows: the thicknessdependent direct band gap Δ is about 0.7 eV for thickness d = 2 nm, γ = 4a/π eV m (a = 0.223 nm is the scale length), η c = ħ 2 /(0.4m 0 ), v c = ħ 2 /(0.7m 0 ), η e = 10 −3 eV, and E f − E c = 0.6 eV; electron mass m 0 = 9.10938 × 10 −31 kg.
The real part and imaginary part of the relative permittivity component ε 11 and ε 22 are shown in Fig. 1(d) and (e), respectively.5]36 The characteristics of the BPSP in nanostructured black phosphorus can be analyzed with the CMT method as shown in Fig. 1(c).When the light wave S +,in passes through the black phosphorus grating, the energy can be coupled into the black phosphorus grating due to the SPR effect.The resonance modes of BPSP can be inspired by the energy amplitude a m (m = 1,2) in the whole region of the black phosphorus grating, where da m /dt = −iωa m .a m is the m th mode of resonant modes with resonance frequency ω m .k m stands for the coupling between cavity modes.The incoming and outgoing waves are depicted by S ±,in(out) .The energy amplitude a 1 and a 2 of the BPSP resonators can be expressed as follows: 37,47 |a〉 and |K〉 represent the field amplitude of the resonant mode and the coupling coefficient between resonators, respectively.〈K| is dependent on |K〉.They can be written as In eqn (3), Ω, Γ i , Γ e , and M matrixes are given as where γ ipq = 1/τ ipq and γ wpq = 1/τ wpq ( p,q = 1,2).In Ω, Γ i , Γ w , and M matrixes, if p ≠ q ( p,q = 1,2), ω pq = 0, τ ipq = 0, and τ wpq = 0; if p = q, ω pq = ω p , τ ipq = τ ip , τ wpq = τ wp , and μ pq = 0.The 1/τ im (m = 1,2) is the decay rate due to the intrinsic loss of the m th mode in the black phosphorus grating and τ im is the lifetime in the decay process due to the intrinsic loss for the m th mode.The 1/τ wm is the decay rate due to the energy coupling from each mode in the light field and τ wm is the lifetime in the process of energy coupling between each mode and the light field.μ pq are the coupling coefficients between two resonant modes.Using boundary conditions of S −in = 0 and eqn (3)-( 5), the transmission function t (ω, E f ) = S +,out /S +,in and reflection function r(ω, E f ) = S −,out /S +,in of this black phosphorus plasmonic grating can be calculated using the following formula: where The absorption efficiency of the system can be expressed as The derived expression A(ω) in eqn (8) with CMT includes the parameters ω m , τ im , τ wm , μ 12 , and μ 21 , which are used to describe the absorbance of black phosphorus.The resonant wavelength λ m (m = 1,2) can be obtained from the FDTD simulation.The parameter ω m (m = 1,2) in eqn (8) can be obtained using ω m = 2πc/λ m .After choosing the appropriate parameters of 1/τ im , 1/τ wm , μ 12 and μ 21 , we can precisely compare the intensity, resonant wavelength and spectral width of the resonance absorption peaks A(ω) obtained by CMT with that of the absorption spectrum simulated by the FDTD method; we can then derive the value of these parameters such as ω m , 1/τ im and 1/τ wm .
The second-order nonlinearity in nanostructured black phosphorus is described by the nonlinear Lorentz force acting on the black phosphorus electrons.The nonlinear theory has been formulated based on the Vlasov-Maxwell equations.The first thing is the source of SHW.SHW comes from a secondorder nonlinear process originating from the surface atoms of the black phosphorus grating, where the broken spatial inversion symmetry will lead to non-zero second-order nonlinear susceptibility.The illumination of the high peak power femtosecond pulse laser upon the black phosphorus will induce second-order nonlinear polarization, which will radiate SHW. 31 The second-order nonlinearity in nanostructured black phosphorus is described by the nonlinear Lorentz force acting on the black phosphorus electrons.The nonlinear theory has been formulated based on the Vlasov-Maxwell equations.The equation of the second-order nonlinear source in the black phosphorus grating can be described as follows: 24 Here, j ð2Þ z = 0, k represents the x, y, and z coordinates; n 0 and m e stand for the black phosphorus-ion density and effective electron mass in the black phosphorus grating, respectively.The S (2) is the second-order nonlinear source.There are two second-order contributions to nonlinear optics, one is the electric part of the Lorentz force ρ (1) E (1) .The surface charge density ρ (1) appears, which diminishes the normal component of the electric field E (1) inside the black phosphorus.Another contribution is the magnetic part of the Lorentz force j (1) × B (1) .The moving electrons are mainly located inside the skin depth of the black phosphorus.Due to the penetrating magnetic field B (1) , the moving electrons feel the magnetic part of the Lorentz force.

Fundamental wave and second harmonic wave
Light can be coupled into the nanostructured black phosphorus flake from the normal direction.There is a propagating wave inside the black phosphorus flake with an effective wavelength matching the grating period.The electric field of the incident light is oriented across the black phosphorus flake.The surface plasmon resonance can be excited if the frequency of the incident light matches the grating period.To study the optical properties in the square array of the black phosphorus flake in Fig. 1(a), the normalized linear-optical absorption (black circles) of FW through the nanostructured black phosphorus flake is investigated with FDTD simulation as shown in Fig. 2(a).Here, the parameters are set as P = 200 nm, L = 150 nm, and d = 2 nm, η e = 10 −3 eV, and E f − E c = 0.6 eV.The red line represents the theoretical CMT result in Fig. 2

(a).
There are two BPSP resonance modes with the first resonant mode at wavelength λ 1 = 7.7 μm and the second resonant mode at wavelength λ 2 = 9.7 μm.The two BPSP resonance modes originate from the dipole resonance, which exists due to the collective oscillations of charged carriers.The parameter ω m (m = 1,2) in eqn ( 8) can be obtained with ω m = 2πc/λ m .The absorption spectra A(ω) in eqn ( 8) can be obtained with appropriate parameters 1/τ w1 , 1/τ w2 , 1/τ i1 , 1/τ i2 , μ 12 and μ 21 .We can precisely compare the intensity, resonant wavelength and spectral width of the resonance absorption peaks A(ω) obtained by CMT with that of the absorption spectrum simulated by the FDTD method in Fig. 2(a).Then, we can obtain the parameters in eqn (8) with 1/τ w1 = 8.6 × 10 10 rad s −1 , 1/τ w2 = 7.3 × 10 10 rad s −1 , 1/τ i1 = 1.41 × 10 12 rad s −1 , 1/τ i2 = 5.3 × 10 12 rad s −1 , μ 12 = μ 21 = 10 11 rad s −1 .We suppose that the coupling coefficients μ 12 and μ 12 are unchanged in the model.The theoretical CMT results (red line) are in good agreement with the results using FDTD simulation (black circles), which can ensure the precise adjustment of the BPSP.The larger spectral width is directly related to the shorter pulse length.The electric field distribution of the individual resonances overlaps spatially, which leads to a coupling of the two resonances.The key feature of the black phosphorus grating is the coupling of a photonic resonance with an electric resonance.
An important aspect of a black phosphorus plasmon is that it generally comprises an electric dipole moment.An incident electromagnetic wave can couple to the dipole and excite the corresponding plasmon oscillation.The oscillating dipole radiates and the associated energy loss significantly broadens the linewidth of the emission spectrum.The electric field inside the black phosphorus is a superposition of a wave decaying on the order of the skin depth and depolarization field.In the case of black phosphorus plasmons, which also comprise a magnetic dipole moment, the linear-optical coupling strength of the electric field to the plasmon is usually stronger than that of the magnetic field.
For varying ε(r) (i.e., ∇ε ≠ 0) in the nanostructured black phosphorus, there is ρ (1,2) = ∇•D (1,2) ≠ 0 in the region of black phosphorus, and ρ (1,2) is the bulk charge density for FW and SHW.The spectra stem from the fact that waves satisfy ∇•D (1,2) ≠ 0 and may induce electric charge oscillation ρ (1) (r)e iωt or ρ (2) (r)e i2ωt .On the surface of black phosphorus, D ð1;2Þ z has abrupt discontinuity across the boundary.The boundary condition is n•D ð1;2Þ z z = σ (1,2) where z is the unit vector along the z-direction and n = nz.σ (1,2) is the surface charge density for FW and SHW, respectively.The dipole-like interactions have been observed in nanostructured black phosphorus nanoflake due to the LSP resonance effect.
The distributions of charge density σ (1) and corresponding electric field E ð1Þ x for the FW at the resonant wavelength λ 1 = 7.7 μm are shown in Fig. 2(b) and (c), respectively.Inside the black phosphorus, one dipole is arranged at the left and right edges corresponding to the polarization of the incident light in Fig. 2(b).The patterns of charge density σ (1) at the left and right edges of black phosphorus are nearly the same, but they have opposite signs that are labeled as ±.The distribution of surface charge density σ (1) in black phosphorus grating is consistent with the local charge dipole oscillations.The electric field E ð1Þ x at wavelength 7.7 μm is seen to be entirely localized inside the black phosphorus region and the edge along the x-direction as shown in Fig. 2(c).Seemingly, the dipole in Fig. 2(b) is flanged by the opposite surface charge and the charges establish a standing BPSP resonance inside the black phosphorus region and edge along the x-direction in Fig. 2(c).Strong localization in the C 1 region indicates that the C 1 region is the principal channel of the energy transportation.After the absorption of the incident wave, a local field locates at the surface of the black phosphorus.This field component is associated with a non-zero charge density near the interface, and its magnitude decays very fast outside the black phosphorus, approximately on the scale of a Thomas-Fermi screening length, which is about 6 nm for black phosphorus. 36he charge density σ (1) and corresponding electric field E ð1Þ x for the FW at resonant wavelength λ 2 = 9.7 μm are seen to be strongly located at the four corners of the black phosphorus as shown in Fig. 2(e) and (f ), which is a typical corner effect.Inside one black phosphorus, two dipoles are arranged at the corners due to the corner effect.Seemingly, the two dipoles are flanged by opposite surface charges that are labeled as ± in Fig. 2(e).These charges establish a standing wave resonance along the edge region of the black phosphorus in Fig. 2

(f ).
Strong localization in the C 2 region indicates that the C 2 region is the principal channel of the energy transportation.The electric field E ð1Þ x has the exponential decay out of the black phosphorus flake.
The black phosphorus flake has mirror symmetry with respect to the x-axis for the fundamental wave.After the con-sideration with the polarization along the x-axis, we have the following two relations: anti-symmetrical charge density σ (1) (−x,y,z) = −σ (1) (x,y,z) and the corresponding symmetrical electric field E ð1Þ x (−x,y,z) = E ð1Þ x (x,y,z).Total electric field intensity (|E (1) | 2 ) profile shown in Fig. 2(d) implies that the total electric field intensity (|E (1) | 2 ) is mostly localized at the left and right edges corresponding to the polarization of the light and four corners for resonant wavelengths λ 1 = 7.7 μm. 14The electric field intensity for resonant wavelengths λ 2 = 9.7 μm is mostly localized around the four corners of the square due to the sharpness in shape in Fig. 2(g).
Beyond these fascinating linear-optical effects, the black phosphorus grating promises unprecedented phenomena in nonlinear optics.The excitation of BPSP produces a strong enhancement of the electric field in the surface.The generation of second harmonics in these enhanced fields is also expected.If the square array of black phosphorus flakes is irradiated with monochromatic light, the SHW can be excited in resonance.
In black phosphorus, the linear polarization P (1) = (−e)n s r (1) can be related to electrons with charge (−e) and density (n s ) being displaced by r (1) from their equilibrium position.Illumination with an electromagnetic plane wave oscillating with a single frequency ω 0 thus results in an electron velocity v (1) = ∂r (1) /∂t.In the simplest case, the velocity v (1) , the polarization P (1) and the incident local electric field E (1) have the same direction, which in turn is perpendicular to the Bloch wave vector k.Furthermore, if the electrons can be viewed individually, the local magnetic field B (1) at the position of an electron corresponds to the value and orientation of the local electric field, hence B (1) ∝ k × E (1) .
It is obvious that such an electron feels the magnetic part of the Lorentz force F L,m = (−e) v (1) × B (1) , which is of second order in perturbation.The second-order oscillation r (2) at frequency 2ω 0 resulting from this force can lead to SHW radiation.With the electric dipole approximation, the electric field of SHW can be related to the electric field of the FW such that 24 where χ (2) stands for a generalized second-order nonlinear susceptibility.
In the case of the incident wavelengths, λ 1 = 7.7 μm and λ 2 = 9.7 μm in Fig. 3(a), we can obtain the second-order nonli-nearity including SHG as well as the SFG and DFG in Fig. 3(b).For the FW wave propagating away from the substrate, the amplitude of the Fourier spectrum for E ð1Þ x with two different wavelengths, λ 1 = 7.7 μm and λ 2 = 9.7 μm, is shown in Fig. 3(a).For the SHW wave propagating away from the substrate, the amplitude of the Fourier spectrum for E We find that the second-order nonlinearity modes at wavelengths λ 3 = 3.85 μm and λ 4 = 4.85 μm are obtained from the incident monochromatic wave with wavelengths λ 1 = 7.7 μm and λ 2 = 9.7 μm due to the SHG effect, respectively.The second-order nonlinearity mode at wavelength λ 5 = 4.3 μm is the sum-frequency field signal for the FW with wavelengths λ 1 and λ 2 .The second-order nonlinearity mode at wavelength λ 6 = 38.1 μm is the difference-frequency field for the FW with two resonant wavelengths λ 1 and λ 2 .The SHW conversion efficiency η is defined in the nonlinear optical process as the expression η = |E (2) (2ω 0 )/E (1) (ω 0 )|.The second-order nonlinear conversion efficiency for the FW with two resonant wavelengths is about 10 −11 .The strongly localized field and noncentrosymmetry can induce an increase in second harmonic nonlinearity signals.
The distributions of charge density σ (2) and corresponding electric field E ð2Þ x for the SHW at the resonant wavelength λ 3 are shown in Fig. 3(g) and (h), respectively.Charge density σ (2) is mainly located at the left and right edges corresponding to the polarization of the light.The patterns of charge density σ (2) at the left and right edges of black phosphorus are nearly the same, and they have the same signs which are labeled as ±.The distribution of surface charge density σ (2) in the black phosphorus grating is consistent with the local charge dipole oscillations.The corresponding electric field E ð2Þ x at wavelength λ 3 is shown in Fig. 3(h).Seemingly, the charges density σ (2) in Fig. 3(g) is flanged by the opposite surface charge and the charges establish a standing BPSP resonance inside the black phosphorus region in Fig. 3(h).Moreover, the standing BPSP resonance has the exponential decay out of the edge along the x-direction in Fig. 3(h).The distributions of charge density σ (2) and corresponding electric field E ð2Þ x at wavelength λ 4 are shown in Fig. 3(i) and ( j), respectively.The electric field E ð2Þ x is seen to be entirely localized inside the black phosphorus region in Fig. 3( j).The distributions of charge density σ (2) and corresponding electric field E ð2Þ x at the wavelength λ 5 are depicted in Fig. 3(k) and (l), respectively.Charge density σ (2) is mainly located at the four edges.The electric field E ð2Þ x is seen to be entirely localized inside the two edge regions along the y-direction, and the electric field E ð2Þ x has the exponential decay out of the edge along the x-direction in Fig. 3(l).The distributions of charge density σ (2) and corresponding electric field E ð2Þ x at the wavelength λ 6 are shown in Fig. 3(m) and (n), respectively.Charge density σ (2) is mainly located at the four corners, four edges and the center region.The dipole-like interactions have been observed in Fig. 3(n) due to the LSP resonance effect.
The black phosphorus flake has mirror symmetry with respect to the x-axis for the second harmonic wave.After the consideration with the polarization along the x-axis, we have the following two relations: symmetrical charge density σ (2) (−x, y,z) = σ (2) (x,y,z) and the corresponding anti-symmetrical electric field E x (x,y,z).The symmetry of charge density σ (2) and corresponding electric field E ð2Þ x is different from that x with the four wavelengths λ 3 , λ 4 , λ 5 , and λ 6 .The distribution of total electric field intensity |E (2) | 2 at resonant wavelength (c) λ 3 , (d) λ 4 , (e) λ 5 , and (f ) λ 6 , respectively.The plots of charge density σ (2) and (e) corresponding electric field E of charge density σ (1) and corresponding electric field E ð1Þ x .The second harmonic wave results from the dipole-like interactions due to the LSP resonance effect (the collective oscillations of anti-symmetrical charge density σ (1) ).

Modulation and sensing application
The intensity, wavelengths and spectral width of resonance peaks can be precisely modulated with the Fermi energy, scattering rate η e and side length L of BP.The black phosphorus square array is covered with the superstrate of an aqueous solution, and the sensitivity is demonstrated by varying the refractive index n 1 of the superstrate.A significant enhancement of the chemical or gas sensing performance of noble metal incorporated BP systems has been reported. 38,39The BPSP resonance performance can also be enhanced in a black phosphorus nanoflake grating compared to that in a single black phosphorus nanoflake. 48irst, we investigate the change in BPSP in black phosphorus by varying the relative Fermi energy E f − E c , which can be adjusted with applied voltage bias or doping concentration. 40,41The parameter η e is fixed with 0. respectively.With the increment of relative Fermi energy E f − E c , the absorption of the first resonant wavelength will get higher while the absorption of the second resonant wavelength will become lower as seen in Fig. 4(a)-(c); the two resonance modes shift to shorter wavelength with increasing Fermi energy.The shift of resonance wavelength and the change in absorption amplitude can be elaborated by CMT theory.The intensity, wavelengths and spectral width of resonance peaks with the FDTD method are in good agreement with the theoretical CMT results.The evolution of simulated linear-optical absorption spectra for different relative Fermi energies E f − E c is presented in Fig. 4(d).With the relative Fermi energy of black phosphorus tuned in the range of 0.15-0.8eV, the first resonance mode can be modulated in the range of 7.1-15 μm while the second resonance mode can be modulated in the range of 8.5-15 μm.This response characteristic in Fig. 4(d) between the relative Fermi energy E f − E c and resonant modes is especially valuable for the modulated sensing application in the optoelectronic devices of black phosphorus.
The parameters 1/τ wm and 1/τ im (m = 1,2) with different  2 .In the process of energy coupling between each mode and the light field, the decay rate for 1/τ w1 and 1/τ w2 with various relative Fermi energies E f − E c are shown in Fig. 4(e).In the process of intrinsic loss for the m th mode, the decay rate for 1/τ i1 and 1/τ i2 with various relative Fermi energies E f − E c are shown in Fig. 4(f ).The lifetime τ wm (m = 1,2) can be modulated with the relative Fermi energy E f − E c .As previously mentioned, adjusting the applied voltage or doping concentration corresponds to a change in the Fermi energy E f .The lifetimes τ wm (m = 1,2) strongly depend on the relative Fermi energy E f − E c , which can be electrostatically tuned, thus resulting in the sensing application with applied voltage bias.To get more insight into the loss in BPSP structure, the absorption spectra with different modulated parameters η e (a) η e = 0.04 eV, (b) η e = 0.02 eV, and (c) η e = 0.002 eV are depicted in Fig. 5 by using the CMT theory (red line) and the FDTD method (black circles), respectively.The simulated absorption spectra with the FDTD method show good agreement with the theoretical CMT results.The absorption spectra become broader and shallower when the scattering rate η e increases from η e = 0.002 eV to η e = 0.04 eV.As the scattering rate η e increases, the absorption for the resonant wavelength 7.7 μm drops from 48% to 14%, while the absorption for the resonant wavelength 9.7 μm drops from 46% to 24% as shown in Fig. 5(a)-(c).The change in the absorption amplitude and width of the spectra with FDTD simulation can also be well elaborated by CMT theory.
The evolution of the linear-optical absorption spectra for different η e is shown in Fig. 5(d).With decreasing the scattering rate η e , the absorption intensities of the two resonance modes decay exponentially.This exponential decay response characteristic is especially valuable for the application of conversion efficiency in the optoelectronic devices of black phosphorus.
The relative Fermi energy E f − E c = 0.6 eV is fixed here.Similar process in Table 1 can be applied for the cases η e = 0.001 eV, 0.002 eV, 0.004 eV, 0.006 eV, 0.008 eV, 0.01 eV, 0.02 eV, 0.03 eV, 0.04 eV, 0.05 eV, 0.06 eV, 0.07 eV and 0.08 eV.The absorption spectra A(ω) can be obtained with appropriate parameters 1/τ w1 , 1/τ w2 , 1/τ i1 and 1/τ i2 for one case, η e .We can then precisely compare the intensity, resonant wavelength and spectral width of the resonance absorption peaks A(ω) with that of the absorption spectrum simulated by the FDTD method.The parameters λ m , 1/τ wm , and 1/τ im from CMT theory and FDTD simulation for the different scattering rate η e are tabulated in Table 2.
The parameters 1/τ wm and 1/τ im (m = 1,2) with different η e in Table 2 have  In the process of the energy coupling between each mode and the light field, the decay rate for 1/τ w1 and 1/τ w2 with various η e are shown in Fig. 5(e).In the process of intrinsic loss for the m th mode, the decay rate for 1/τ i1 and 1/τ i2 with various η e are plotted in Fig. 5(f ).Therefore, the lifetimes τ w1 , τ w2 , τ i1 and τ i2 can be modulated with the scattering rate η e .There is a good broad tunability of the BPSP in black phosphorus by changing the side length L. To get more insight into the BPSP structure, the absorption spectra with different modulated side lengths (a) L = 160 nm, (b) L = 120 nm, and (c) L = 80 nm are shown in Fig. 6 using the CMT theory (red line) and FDTD method (black circles), respectively.The change in resonance wavelength, absorption amplitude, and width of the spectra can be well elaborated by CMT theory.The simulated absorption spectra with the FDTD method are in good agreement with the theoretical CMT results.
Additionally, the optical properties of the BP grating can be modulated by their side length.The evolution of the simulated  linear-optical absorption spectra for different side lengths L was investigated with the FDTD method as shown in Fig. 6(d).
The two resonance modes are red-shifted with increasing edge length L. This quasi-linear response characteristic in Fig. 6(d) between the edge length L and resonant wavelength is especially valuable for the application in the optoelectronic devices of black phosphorus.
The parameter E f − E c and the scattering rate η e are fixed with 0.6 eV and 0.001 eV for various side lengths L, respectively.A similar process in Table 1 can be applied for the cases L = 70 nm, 80 nm, 90 nm, 100 nm, 110 nm, 120 nm, 130 nm, 140 nm, 150 nm, 160 nm, and 170 nm using CMT theory and FDTD simulation.The appropriate parameters 1/τ w1 , 1/τ w2 , 1/τ i1 and 1/τ i2 for one case with η e can be obtained after precisely comparing resonance absorption peaks A(ω) with the absorption spectrum simulated by the FDTD method.We can then obtain the parameters λ m , 1/τ wm , and 1/τ im (m = 1,2) for all the cases with side length L in Table 3.
The decay rates 1/τ wm and 1/τ im (m = 1,2) with different side lengths L in Table 3 have the fitting expressions as follows: 1/τ w1 = −1.5 × 10 7 L 2 + 3.7 × 10 9 L − 1.5 × 10 11 , 1/τ w2 = −9.3× 10 6 L 2 + 1.9 × 10 9 L + 1.2 × 10 10 , 1/τ i1 = −1.2× 10 8 L 2 + 3.1 × 10 10 L − 5.6 × 10 11 , 1/τ i2 = 4.3 × 10 8 L 2 − 4.7 × 10 10 L 2 + 2.5 × 10 12 .In the process of energy coupling between each mode and the light field, the decay rates for 1/τ w1 and 1/τ w2 with various side lengths L are shown in Fig. 6(e).In the process of intrinsic loss for the m th mode, the decay rates for 1/τ i1 and 1/τ i2 with various side lengths L are shown in Fig. 6(f).The lifetimes τ w1 , τ w2 , τ i1 and τ i2 can be changed with the modulated side length L. The expressions of the lifetimes for resonant modes in CMT are obtained from theoretical fitting of exact values in the FDTD simulation.The theoretical descriptions of the tunable lifetimes will be useful in applying the theory for future black phosphorus applications by changing the side length of BP.Since the intensities, wavelengths and spectral widths of resonance peaks can be precisely controlled with the Fermi energy E f , the scattering rate η e and side length L, we now discuss the sensing performance of the proposed nanosensor.
Here, the parameters are those used in Fig. 2(a).In the following discussion, we only modify the refractive index n 1 of superstrates such as an aqueous solution.
The absorption spectra with different refractive indexes of superstrates (a) n 1 = 1.8,(b) n 1 = 1.4,and (c) n 1 = 1.0, are shown in Fig. 7 with the CMT (red line) and FDTD method (black circles), respectively.The change in resonance wavelength, absorption amplitude and broad of spectra with FDTD simulation can be well elaborated by CMT theory.Additionally, the evolution of the simulated linear-optical absorption spectra for different refractive indexes of superstrates, n 1 , was investigated with the FDTD method as shown in Fig. 7(d).The two resonance modes are red-shifted with increasing refractive index n 1 .This quasi-linear response characteristic in Fig. 7(d) between the refractive index n 1 and the resonant wavelength is especially valuable for the sensing application of black phosphorus.
Similar processes to that in Table 1 can be applied for the cases of the refractive indexes of superstrates, where n 1 = 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, and 1.8, using CMT theory and FDTD simulation.The appropriate parameters 1/τ w1 , 1/τ w2 , 1/τ i1 and 1/τ i2 for one case, n 1 , can be obtained after precisely comparing the intensity, resonant wavelength and spectral width of resonance absorption peaks A(ω) in eqn (8) with those of absorption spectra simulated by FDTD method.The parameters λ m , 1/τ wm , and 1/τ im (m = 1,2) can then be obtained for all the cases of superstrate refractive indexes n 1 in Table 4.
The parameters 1/τ wm and 1/τ im (m = 1,2) with different refractive indexes n 1 of the superstrates in Table 4 have the fitting expressions as follows: 1/τ w1 = −3.3× 10 9 n 1 + 8.9 × 10 10 , 1/τ w2 = −2.2× 10 9 n 1 + 7.7 × 10 10 , 1/τ i1 = −4 × 10 10 n 1 + 1.4 × 10 12 , 1/τ i2 = 1.3 × 10 11 n 1 + 5.1 × 10 12 .In the process of energy coupling between each mode and light field, the decay rate for 1/τ w1 and 1/τ w2 with various refractive indexes n 1 are shown in Fig. 7(e).In the process of intrinsic loss for the m th mode, the decay rate for 1/τ i1 and 1/τ i2 with various n 1 of the superstrates are shown in Fig. 7(e).The lifetimes τ w1 , τ w2 , τ i1 and τ i2 can be modulated with the modulated refractive index n 1 of the superstrate.The expressions of the lifetimes for resonant modes in CMT were obtained from the theoretical fitting of the exact values in the FDTD simulation.The theoretical descriptions of the tunable lifetime will make it useful in applying the theory for future black phosphorus sensing applications by changing the refractive index n 1 of superstrates such as aqueous solutions.
It is known that there are two common ways to describe the figure of merit (FOM) of a sensor. 42One is related to the resonance wavelengths shift at a certain refractive index n 1 , which is defined as FOM.The other is defined as FOM i , which is related to the intensity variation at a certain resonance wavelength λ.In order to quantitatively describe the sensitivity of the purposed sensor, we investigated the change in two resonance wavelengths by varying the refractive index n 1 of the superstrate.We define the sensitivity as follows: where the resonance wavelengths λ(n 1 + Δn), λ(n 1 ), intensity T (λ, n 1 + Δn) and T (λ,n 1 ) can be obtained by the FDTD simulation.ΔT and Δλ are the intensity variation and resonant wavelength variation caused by the refractive index change, Δn, of the superstrate environment, respectively.The sensing response FOM and FOM i can be obtained when the superstrate refractive index n 1 increases from 1.0 to 1.8 as shown in Fig. 7(f).By using eqn (11), we get a FOM of 0.293 for resonance peak 1 at λ 1 and 0.303 for resonance peak 2 at λ 2 at refractive index n 1 = 1.It can be seen that the FOM for peak 1 and peak 2 increase accordingly as the refractive index n 1 increases from 1.0 to 1.4.However, FOM for peak 1 and peak 2 decrease as the refractive index n 1 continues to increase from 1.4 to 1.8.By using eqn (12), we get a FOM i of −24 for peak 1 and −3.0 for peak 2 at refractive index n 1 = 1.It can be seen that FOM i for peak 1 and peak 2 increase as the refractive index n 1 increases from 1.0 to 1.8.

Conclusion
We have investigated the modulated surface plasmons and their modulated lifetimes for nanostructured black phosphorus.The expressions of the lifetimes for linear-optical resonant modes have been investigated with the FDTD simulations and CMT.The lifetimes of plasmonic absorption modes are tunable with the Fermi energy, collision frequency, side length and superstrate.The theoretical descriptions and data fitting are useful in applying the methods for future black phosphorus modulation.The strongly confined black phosphorus plasmonic waves result from an enhancement of the local field.The strongly localized field induces a desired increase in the second-order nonlinearity.It may stimulate new insights to visualize general nonlinear nanophotonic processes in black phosphorus.The proposed configuration and results could provide guidance for the research of the lifetimes of modes in black phosphorus, and the design of black phosphorus-based optoelectronic modulation devices and sensing applications.
where m 1 and m 2 denote the elements along the x-direction and y-direction, respectively.14][15]36 In response to a local electric field E (1) , the equation of motion of an electron with inplane electron effective masses and charge (-e) is where r (1) is the electron displacement from their equilibrium position and γ bp is the collision frequency of the Drude model for black phosphorus, which represents Ohmic damping.In linear optics, for ξ = 0 in eqn (14), this leads to an oscillator due to the resonance of the collective oscillations of electrons.
In the frequency domain, driving the electron with a harmonic electric field and frequency results in the oscillation amplitude: The electronic density n s of black phosphorus with thickness d is given by a diagonal tensor: 43,[12][13][14][15] where n s1 and n s2 denote the elements of electronic density n s in black phosphorus along the x-and y-directions, respectively.ħ is the reduced Planck's constant.k B is the Boltzmann constant and T is the temperature.E f is the Fermi energy, E c is the conduction band edge energy, and E f − E c is defined as the relative Fermi level herein.The incident electromagnetic wave can couple with the dipole and excite the corresponding plasmon oscillation in black phosphorus.The electric field inside the black phosphorus is a superposition of a wave decaying on the order of the skin depth and of the depolarization field.In the case of plasmons, which also comprise a magnetic dipole moment, the linear-optical coupling strength of the electric field to the plasmon is usually stronger than that of the magnetic field.The linear polarization P (1) = (−e) n s r (1) can be related to electrons with charge (−e) and density n s being displaced by r (1) from their equilibrium position.In linear optics, the linear polarization P (1) is simply proportional to the electric field: This allows us to define an equivalent relative permittivity diagonal tensor ε D (ε 11 , ε 22 , ε 33 ) to describe the black phosphorus with thickness d.
The conductivity of few-layer black phosphorus is given as the equivalent permittivity tensor component of black phosphorus is given by the following: Fundamental wave and second harmonic wave Herein, the black phosphorus flake is treated as a thin layer with the anisotropic Drude model.The linear optical response of the black phosphorus in Maxwell equations is given by the following: 24 @B ð1Þ @t ¼ À∇ Â E ð1Þ ð21Þ @E ð1Þ @t ¼ c 2 ∇ Â B ð1Þ À 1 ε 0 j ð1Þ ð22Þ j ð1Þ ¼ e 2 n s m bp 1 ðγ bp À iωÞ E ð1Þ ð23Þ where j ð1Þ z = 0, and B (1) , E (1) and J (1) represent the magnetic flux intensity vector, electric field vector and current density vector of FW, which are marked with the superscript (1).The magnetic flux intensity and electric field of FW can be obtained by discretizing eqn ( 21) and ( 22) in space and time using the FDTD method.From eqn (23), we have the equation for the x component −iωj For the second order of the electric field, ξ ≠ 0 in eqn (14).ξ [r (1) ] 2 is the driving term for the second-order displacement r (2) .The second-order polarization is P (2) ∝ r (2) and the secondorder displacement is r (2) ∝ [r (1) ] 2 .The second-order nonlinearity in nanostructured black phosphorus is described by the nonlinear Lorentz force acting on the black phosphorus electrons.The nonlinear theory has been formulated based on the Vlasov-Maxwell equations.The equations of second-order nonlinearity in black phosphorus grating can be described as follows: 24 @B ð2Þ @t ¼ À∇ Â E ð2Þ ð26Þ @E ð2Þ @t ¼ c 2 ∇ Â B ð2Þ À 1 ε 0 j ð2Þ ð27Þ @j ð2Þ @t ¼ Àγ bp j ð2Þ þ e 2 n s1 m bp E ð2Þ þ S ð2Þ : ð28Þ Here The j (2) , B (2) , and E (2) represent the current density vector, magnetic flux intensity vector, and the electric field vector of SHW, which are marked with the superscript (2).

Numerical simulation
The 3D-FDTD method was used for the calculations.There are two computational loops for the calculations of the above firstorder eqn ( 21)-( 23), and the second-order eqn ( 26)-( 28) in the program.Yee's discretization scheme was utilized and all electric and magnetic components in a cubic grid can be defined.The perfectly matched absorbing boundary conditions were used below and on top of the computational space along the z-direction, and the periodic boundary conditions were employed on the boundaries along the x-and y-directions.We only considered one unit cell in the computational space.The incident plane wave, polarized along the x-direction, propagates along the z-direction.The grid size of 2 nm was used to mesh the BP thin film, SiO 2 substrate and the superstrate of an aqueous solution, respectively.The refractive index of the SiO 2 was set as 1.4.

Measurement nonlinear signals of FW and SHW
It is quite necessary to discuss the feasibility of the experimental demonstration.Similar to graphene and monolayer MoS 2 , the patterned black phosphorus can be fabricated using the standard procedure of nanofabrication.The processes are described as follows.First, 300 nm photoresist (i.e.ZEP520 or PMMA) is coated on the top of the black phosphor film.Then, the patterned photoresist was defined using electron beam lithography.Finally, the black phosphorus pattern was generated with reactive ion etching (i.e.O 2 plasma or CHF 3 ).The reflection and transmission of the structure for the fundamental wave can be measured by Fourier Transform Infrared Spectroscopy (FTIR).For the nonlinear response of the patterned black phosphorus, ultrafast laser pumps at two fundamental wavelengths were incident on the structure.The signal of nonlinear response was collected by the single mode fiber, which was integrated with the spectrometer. 10,22,41,44easurements of SHG signals were also reported. 25,31Other measurements of third-order nonlinear signals were also reported. 45,46

Fig. 1
Fig. 1 (a) The nanostructured black phosphorus nanoflake on SiO 2 substrate.The input light wave is polarized along the x-direction and propagates along the z-direction.(b) Side view of the black phosphorus flake.(c) The CMT model with two resonant modes of the nanostructured black phosphorus system.(d) The real part and (e) imaginary part of the relative permittivity tensor component ε 11 and ε 22 , respectively.

Fig. 3
Fig. 3 (a) The amplitude of the Fourier spectrum for FW E ð1Þ x with two wavelengths λ 1 = 7.7 μm and λ 2 = 9.7 μm, (b) the amplitude of the Fourier spectrum for SHW E
001 eV for various relative Fermi levels E f − E c .The absorption spectra with different relative Fermi energies (a) E f − E c = 0.7 eV, (b) E f − E c = 0.5 eV, and (c) E f − E c = 0.3 eV are shown in Fig. 4 by the CMT theory (red line) and FDTD method (black circles),

Fig. 4
Fig. 4 The absorption with different relative Fermi energies (a) E f − E c = 0.7 eV, (b) E f − E c = 0.5 eV, and (c) E f − E c = 0.3 eV using the CMT (red line) and FDTD method (black circles), respectively.(d) The evolution of linear-optical absorption for different relative Fermi energies with using the FDTD method.(e) 1/τ w1 and 1/τ w2 with various relative Fermi energies E f − E c .(f) 1/τ i1 and 1/τ i2 with various relative Fermi energies E f − E c .

Table 1 λ
m , 1/τ wm , and 1/τ im from CMT and FDTD simulation for the different E f − E c

Table 2 λ
m , 1/τ wm , and 1/τ im from CMT and FDTD simulation for the different η e

Table 3 λ
m , 1/τ wm , and 1/τ im from CMT and FDTD simulation for length L