Anisotropic Photonics Topological Transition in Hyperbolic Metamaterials Based on Black Phosphorus

Based on in-plane anisotropy of black phosphorus (BP), anisotropic photonics topological transition (PTT) can be achieved by the proposed hyperbolic metamaterials structure, which is composed of alternating BP/SiO2 multilayer. Through effective medium theory and calculated iso-frequency contour, PTT can be found by carefully choosing the incident plane and other parameters. With the finite element method and transfer matrix method, a narrow angular optical transparency window with angular full width at half maximum of 1.32° exists at PTT. By changing the working wavelength, thickness of SiO2, or electron doping of black phosphorus, the incident plane of realizing PTT can be modulated, and anisotropic PTT is achieved.

Recently, atomically thin two-dimensional (2D) materials, such as graphene, black phosphorus (BP), and hexagonal boron nitride (BN), have shown many extraordinary electronic and photonic properties [29]. As 2D materials have developed into a research hotspot, HMMs based on 2D materials has gradually attracted wide attention. As one of the most popular 2D materials, graphene can replace the role of metal in HMMs of the multilayer structure due to negative permittivity in infrared region [30]. Graphene-based HMMs have been extensively investigated in radiative heat transfer [31], negative refraction [32], slow light effect [33], and perfect absorber [34]. Compared with graphene, BP also presents metallic behavior with negative permittivity in mid-infrared region. Differently, Nanomaterials 2020, 10 BP have strongly inherent in-plane anisotropy, and many interesting applications have been widely explored based on this property, including anisotropic acoustic plasmons [35], polarization sensitive resonators [36], photonic spin hall effect [37], and chirality [38]. In addition, researches on BP-based HMMs has also been reported in anisotropic absorber and biaxial HMMs [39,40], but this has yet to be fully investigated. Particularly, BP-based HMMs have not been used to realize anisotropic photonic topological transition (PTT).
In this work, we theoretically propose to construct a BP-based HMMs structure consisting of alternating BP/SiO 2 multilayer, which can realize anisotropic PTT thanks to inherent in-plane anisotropy of BP. The results obtained by finite element method (FEM) simulations and transfer matrix method (TMM) both demonstrate that a narrow angular optical transparency window can be achieved at PTT. And the structure's IFC can transit from open hyperboloid to closed ellipsoid by changing the angle ϕ of the incidence plane. When the working wavelength λ = 5 µm, the ϕ of incident plane where PTT appears can be modulated from 0 • to 61 • by changing electron doping of BP and can be modulated from 79 • to 0 • by changing the thickness of SiO 2 . In addition, the angle ϕ for PTT can be modulated from 30 • to 56 • by changing the working wavelength. Our findings pave a new way in anisotropic angle-dependent optical applications.

Design and Theories
As shown in Figure 1a, the proposed BP-based HMMs structure consists of periodic multilayer of BP and SiO 2 layers, which is surrounded by air. The top view of Figure 1a is shown in Figure 1b. As depicted in Figure 1, a p-polarized light at wavelength λis incident on the side of the proposed BP-based HMMs structure. Here, incident plane is the p-z plane, and the angle between the p-z and y-z plane is ϕ. The thickness of SiO 2 is t d = 300 nm. The thickness of BP is given by t bp = n × a z /2, where n (= 3) is the number of layers of BP and a z (= 10.7 Å) is lattice constant in the out-of-plane direction [41]. The permittivity of SiO 2 is 1.82. d (= t d + t bp ) is the thickness of periodic unit of the BP-based HMMs structure. As shown in the inset of Figure 1a, the atoms in monolayer BP are covalently bonded to form a puckered honeycomb structure which leads to unique in-plane anisotropic optical property. In our proposed BP-based HMMs structure, the x and y directions correspond to the zigzag (ZZ) and armchair (AC) directions of BP, respectively. For BP, translational symmetry is broken in the z direction, and it has a direct energy gap at the Г point [41]. Thus, a low-energy in-plane Hamiltonian is used to describe the systematic behavior around the Г point based where γ and β describe the effective couplings between the conduction and valence bands. E c and E υ are the first conduction and valence band edge energies in BP. Near the Г point, in-plane effective electrons masses along the AC and ZZ directions can be obtained by [42]: For 3-layer BP, the layer-dependent bandgap ∆ is 1.1 eV [43]. The other parameters used here are η c =h 2 /0.4m 0 , υ c =h 2 /0.4m 0 , and γ = 4a/π eVm. a (= 0.223 nm) is the scale length of the BP and π/a is the width of the Brillouin zone. m 0 = 9.10938 × 10 −31 kg is the standard electron rest mass. The in-plane anisotropic conductivity of BP can be described by a simple semiclassical Drude model. The conductivity of in-plane BP along the AC and ZZ crystalline directions are given as [44] σ AC,ZZ = iD AC,ZZ π(ω + iη/ ) where D AC, ZZ = πe 2 ρ/m AC, ZZ , is the Drude weight, which is dependent on the electron charge. η (= 10 meV) is used to define the BP relaxation rate. ρ (= 5 × 10 13 cm −2 ) is the electron doping. With the angle ϕ, the conductance matrix connecting the surface current and electric field can be expressed as σ = [σ pp , σ ps ; σ sp , σ ss ], where σ pp = σ AC cos2ϕ + σ ZZ sin2ϕ, σ ss = σ AC sin2ϕ + σ ZZ cos2ϕ, and σ sp = σ ps = (σ ZZ − σ AC ) sinϕcosϕ [45]. The cross conductivity σ sp vanishes for isotropic 2D materials such as graphene. Hence, the effective permittivity of BP in p-axis and z-axis directions can be derived by where ε r ( = 5.76) is the relative permittivity of BP [39] and ε 0 is the vacuum permittivity.
Here, the working wavelength is λ = 5 µm. Obviously, the length d of periodic unit of the BP-based HMMs structure is sufficiently small compared to the working wavelength, so the proposed multilayer structure can be modeled as an anisotropic effective medium by the effective medium theory (EMT) [46]. In p-z plane, we define the p-axis component of effective permittivity as ε p and the z-axis component of effective permittivity as ε z . Based on the EMT, ε p and ε z can be expressed as follows, In k space, the IFC of BP-based HMMs for p-polarized light is given by where k p and k z are the components of the wavevector along pand z-directions respectively, k 0 (= ω/c) is the wavevector of light in the vacuum.
Nanomaterials 2020, 10, x FOR PEER REVIEW 2 of 11 Differently, BP have strongly inherent in-plane anisotropy, and many interesting applications have been widely explored based on this property, including anisotropic acoustic plasmons [35], polarization sensitive resonators [36], photonic spin hall effect [37], and chirality [38]. In addition, researches on BP-based HMMs has also been reported in anisotropic absorber and biaxial HMMs [39,40], but this has yet to be fully investigated. Particularly, BP-based HMMs have not been used to realize anisotropic photonic topological transition (PTT).
In this work, we theoretically propose to construct a BP-based HMMs structure consisting of alternating BP/SiO2 multilayer, which can realize anisotropic PTT thanks to inherent in-plane anisotropy of BP. The results obtained by finite element method (FEM) simulations and transfer matrix method (TMM) both demonstrate that a narrow angular optical transparency window can be achieved at PTT. And the structure's IFC can transit from open hyperboloid to closed ellipsoid by changing the angle φ of the incidence plane. When the working wavelength λ = 5 μm, the φ of incident plane where PTT appears can be modulated from 0° to 61° by changing electron doping of BP and can be modulated from 79° to 0° by changing the thickness of SiO2. In addition, the angle φ for PTT can be modulated from 30° to 56° by changing the working wavelength. Our findings pave a new way in anisotropic angle-dependent optical applications.

Design and Theories
As shown in Figure 1a, the proposed BP-based HMMs structure consists of periodic multilayer of BP and SiO2 layers, which is surrounded by air. The top view of Figure 1a is shown in Figure 1b. As depicted in Figure 1, a p-polarized light at wavelength λ is incident on the side of the proposed BP-based HMMs structure. Here, incident plane is the p-z plane, and the angle between the p-z and y-z plane is φ. The thickness of SiO2 is td = 300 nm. The thickness of BP is given by tbp = n * az/2, where n (= 3) is the number of layers of BP and az (= 10.7 Å ) is lattice constant in the out-of-plane direction [41]. The permittivity of SiO2 is 1.82. d (= td + tbp) is the thickness of periodic unit of the BP-based HMMs structure. As shown in the inset of Figure 1a, the atoms in monolayer BP are covalently bonded to form a puckered honeycomb structure which leads to unique in-plane anisotropic optical property. In our proposed BP-based HMMs structure, the x and y directions correspond to the zigzag (ZZ) and armchair (AC) directions of BP, respectively. For BP, translational symmetry is broken in the z direction, and it has a direct energy gap at the Г point [41]. Thus, a low-energy in-plane Hamiltonian is used to describe the systematic behavior around the Г point based The top view of (a). A p-polarized light in p-z plane is incident on the side of the proposed BP-based HMMs structure from air. ϕ is the angle between the incident p-z plane and y-z plane.
According to Equation (5), we calculate the p-axis and z-axis components of effective permittivity under different ϕ when the working wavelength λ is 5 µm. We find that ε z (= 1.885) is a constant positive value and ε p is a variable with ϕ. As shown in Figure 2a, Re (ε p ) changes from negative to positive values as ϕ increases from 0 • to 90 • , and Im (ε p ) is always near 0. In addition, when Re (ε p ) = 0, we can regard it as PTT point. Here, the PTT refers to the optical topological transition of HMMs' IFC instead of the topological phase transition. When the IFC of structure transitions between open hyperboloid and closed ellipsoid, PTT point will exist. The regime, which is very close to PTT point, is also known as the epsilon-near-zero (ENZ) regime [21]. At PTT point, HMMs can significantly suppress the diffraction and scattering of incident light, which can provide a new way for efficiently manipulating light-matter interactions at nanoscales. For our proposed BP-based HMMs structure, it will change from Type II HMM to elliptic dispersion near PTT point. Based on Equation (6), we calculate the complex wavevector k p /k 0 of BP-based HMMs' IFC at ϕ = 46 • , as shown in Figure 2b,c. For an ideal lossless medium, HMMs' IFC will degenerate to two points at PTT point [23], which means that only the light with pure wavevectors along p-axis is through metamaterials and an angular transparency window is achieved along the p-axis direction. Here, for our BP-based HMMs, the topology of IFC at PTT point maintains a narrow hyperboloid as shown Figure 2b. Thus, the light with very small k z wavevectors are allowed to propagate inside the metamaterials. As shown in Figure 2c, Im (k p /k 0 ) as function of k z /k 0 exhibits a conical dispersion. The inset of Figure 2c shows that the conical dispersion achieves a degenerate state at the origin and the light with wavevector along p-axis has Im(k p /k 0 ) = 0. This indicates that the light with wavevector along p-axis will not be affected by absorption losses. However, for the light with small k z wavevectors, the intrinsic loss of materials will continue to exist due to Im(k p /k 0 ) is not close to zero in this moment. In general, even though the light with small k z wavevectors can propagate inside the metamaterial, the existence of intrinsic loss makes it possible to suppress the light propagation away from p-axis direction, which is helpful for maintaining a narrow angular optical transparency window. Based on the above theoretical analysis, our proposed HMMs provides a possible way to realize a narrow angular optical transparency window when the incident p-polarized light in p-z plane due to PTT. is very close to PTT point, is also known as the epsilon-near-zero (ENZ) regime [21]. At PTT point, HMMs can significantly suppress the diffraction and scattering of incident light, which can provide a new way for efficiently manipulating light-matter interactions at nanoscales. For our proposed BP-based HMMs structure, it will change from Type Ⅱ HMM to elliptic dispersion near PTT point. Based on Equation (6), we calculate the complex wavevector kp/k0 of BP-based HMMs' IFC at φ = 46°, as shown in Figure 2b,c. For an ideal lossless medium, HMMs' IFC will degenerate to two points at PTT point [23], which means that only the light with pure wavevectors along p-axis is through metamaterials and an angular transparency window is achieved along the p-axis direction. Here, for our BP-based HMMs, the topology of IFC at PTT point maintains a narrow hyperboloid as shown Figure 2b. Thus, the light with very small kz wavevectors are allowed to propagate inside the metamaterials. As shown in Figure 2c, Im (kp/k0) as function of kz/k0 exhibits a conical dispersion. The inset of Figure 2c shows that the conical dispersion achieves a degenerate state at the origin and the light with wavevector along p-axis has Im(kp/k0) = 0. This indicates that the light with wavevector along p-axis will not be affected by absorption losses. However, for the light with small kz wavevectors, the intrinsic loss of materials will continue to exist due to Im(kp/k0) is not close to zero in this moment. In general, even though the light with small kz wavevectors can propagate inside the metamaterial, the existence of intrinsic loss makes it possible to suppress the light propagation away from p-axis direction, which is helpful for maintaining a narrow angular optical transparency window. Based on the above theoretical analysis, our proposed HMMs provides a possible way to realize a narrow angular optical transparency window when the incident p-polarized light in p-z plane due to PTT.

Results and Discussion
As shown in Figure 3a, a p-polarized light is incident on the BP-based HMMs structure at an incident angle θ in the p-z plane. Here, the angle between the p-z and y-z plane is φ = 46° and working wavelength λ = 5 μm. We simulate the transmission, reflection, and absorption of the

Results and Discussion
As shown in Figure 3a, a p-polarized light is incident on the BP-based HMMs structure at an incident angle θ in the p-z plane. Here, the angle between the p-z and y-z plane is ϕ = 46 • and working wavelength λ = 5 µm. We simulate the transmission, reflection, and absorption of the multilayer structure with 138 µm width along the p-axis. In this work, all numerical simulation results are obtained by the commercial software COMSOL Multiphysics (COMSOL Multiphysics 5.4, Stockholm, Sweden) based on FEM. In our work, periodic boundary conditions are adopted in the z-axis direction. We verify the position of PTT point by EMT. In addition, we theoretically analyze the propagation features of the proposed BP-based HMMs structure by using transfer matrix method (TMM) under p-polarized light. Based on Maxwell's equations and boundary conditions, the magnetic field between adjacent layers can be related via a transfer matrix, which can be obtained as follows [47] be obtained as follows [47] cos( ) Here, we divide BP-based HMMs structure into k layers along the p-axis direction. The subscript i corresponds to the propagation of light in the i-th layer with a thickness of di and qp = kp/(k0)εz.
According to Equation (6) where q0 = cosθ, the reflection (R(θ)) and transmission (T(θ)) of the structure can be obtained by |r(θ)| 2 and |t(θ)| 2 , respectively. Further, the absorption can be written as A(θ) = 1−R(θ)−T(θ).  Here, we divide BP-based HMMs structure into k layers along the p-axis direction. The subscript i corresponds to the propagation of light in the i-th layer with a thickness of d i and q p = k p /(k 0 )ε z .
According to Equation (6), we can get k p = ε z k 2 0 − (ε z /ε p )k 2 z . Here, k z (= k 0 sinθ) is the z component of incident light wavevector. The total transfer matrix (M(θ)) connecting the fields at the incident end and the exit end can be expressed as By means of the TMM, the reflection and transmission coefficients can be calculated as t(θ) = 2 M 11 (θ) + M 22 (θ) + q 0 M 12 (θ) + (1/q 0 )M 21 (θ) (10) Nanomaterials 2020, 10, 1694 6 of 10 where q 0 = cosθ, the reflection (R(θ)) and transmission (T(θ)) of the structure can be obtained by |r(θ)| 2 and |t(θ)| 2 , respectively. Further, the absorption can be written as A(θ) = 1−R(θ)−T(θ). Figure 3b depicts the FEM-simulated and TMM-calculated optical transmission (T), reflection (R) and absorption (A) of p-polarized light as a function of incident angle (θ) for the proposed BP-based HMMs structure when ϕ = 46 • . Obviously, the results obtained by numerical simulation (FEM) are in good agreement with the theoretical calculation (TMM). The results indicate that a narrow angular optical transparency window with angular full width at half maximum (FWHM) of 1.32 • appears. The transmittance can reach 99.7% at θ = 0 • , but the transmittance is only 0.1% at θ = 2 • . The reason is that the incident light will be affected by the material's inherent losses and the energy will be attenuated with further propagation in the medium except for the wavevector of light along p-axis direction. Moreover, the absorption in the Figure 3b drops to almost zero when θ = 0 • , which can also verify the results obtained in Figure 2c. So, the PTT point of proposed BP-based HMMs structure happens when ϕ = 46 • , working wavelength λ = 5 µm, electron doping ρ = 5 × 10 13 cm −2 , and the thickness of SiO 2 t d = 300 nm.
Then, the influence of electron doping of BP on the position of PTT point for the proposed BP-based HMMs structure will be discussed. Figure 4a shows numerically simulated optical transmission as a function of the angle ϕ between the incidence plane (p-z plane) and y-z plane, and incident angle θ when the electron doping of BP is 2.42 × 10 13 , 3 × 10 13 , 4 × 10 13 , and 10 × 10 13 cm −2 , respectively. PPT happens when the narrowest angular optical transparency window is achieved. Obviously, the angle ϕ, which is corresponding to the PTT, is increasing from 0 • to 61 • as the electron doping of BP increases from 2.42 × 10 13 to 10 × 10 13 cm −2 when the other parameters unchanged. The FWHMs for narrowest angular optical transparency window are all smaller than 1.328 • in Figure 4a. Figure 4b,c show Re(ε p ) and Im(ε p ) as a function of electron doping ρ of BP and angle ϕ based on Equation (5), respectively. Here, the purple dots represent the position of the PTT point obtained by simulation calculation. The cyan solid line in the Figure 4b represents Re (ε p ) = 0, and the cyan dashed line in the Figure 4c represents Im (ε p ) = 0.1. We find that these purple dots which are obtained by simulation satisfy Re (ε p ) ≈ 0 and Im (ε p ) < 0.1. It means that the PTT point obtained by simulation calculation is consistent with that obtained by theoretical calculation. Besides, as shown in the inset of Figure 4a, we can also achieve a wide range of narrow angular optical transparency window from ϕ = 0 • to ϕ = 10 • due to the existence of ENZ regime when ρ = 2.42 × 10 13 cm −2 , and the FWHM of the angular optical transparency window stays around 1.3 • .
The influence of the thickness t d of SiO 2 on the position of PTT point is also discussed in detail. Figure 5a shows numerically simulated optical transmission as a function of the angle ϕ between the incidence plane (p-z plane) and y-z plane, and incident angle θ when the thickness of SiO 2 is 30, 170, 480, and 620 nm, respectively. The FWHMs for narrowest angular optical transparency window are all smaller than 1.334 • in Figure 5a. The angle ϕ, which is corresponding to the PTT point, is decreasing from 79 • to 0 • as the thickness of SiO 2 increases from 30 to 620 nm when the other parameters unchanged. Figure 5b,c show Re(ε p ) and Im(ε p ) as a function of t d and angle ϕ based on Equation (5), respectively. Here, the purple dots represent the position of the PTT point obtained by simulation calculation. The cyan solid line in the Figure 5b represents Re (ε p ) = 0, and the cyan dashed line in the Figure 5c represents Im (ε p ) = 0.1. It is found that these purple dots obtained by simulation satisfy Re (ε p ) ≈ 0 and Im (ε p ) < 0.1. It means that the PTT point obtained by simulation calculation agrees with that obtained by theoretical calculation. Besides, as shown in the inset of Figure 5a, a wide range of narrow angular optical transparency window is also achieved from ϕ = 0 • to ϕ = 10 • due to the existence of ENZ regime when t d = 620 nm. and the cyan dashed line in the Figure 4c represents Im (εp) = 0.1. We find that these purple dots which are obtained by simulation satisfy Re (εp) ≈ 0 and Im (εp) < 0.1. It means that the PTT point obtained by simulation calculation is consistent with that obtained by theoretical calculation. Besides, as shown in the inset of Figure 4a, we can also achieve a wide range of narrow angular optical transparency window from φ = 0° to φ = 10° due to the existence of ENZ regime when ρ = 2.42 × 10 13 cm −2 , and the FWHM of the angular optical transparency window stays around 1.3°. The influence of the thickness td of SiO2 on the position of PTT point is also discussed in detail. Figure 5a shows numerically simulated optical transmission as a function of the angle φ between the incidence plane (p-z plane) and y-z plane, and incident angle θ when the thickness of SiO2 is 30, 170, 480, and 620 nm, respectively. The FWHMs for narrowest angular optical transparency window are all smaller than 1.334° in Figure 5a. The angle φ, which is corresponding to the PTT point, is decreasing from 79° to 0° as the thickness of SiO2 increases from 30 to 620 nm when the other parameters unchanged. Figure 5b Figure 5a, a wide range of narrow angular optical transparency window is also achieved from φ = 0° to φ = 10° due to the existence of ENZ regime when td = 620 nm. Finally, we discuss the influence of working wavelength λ on the position of PTT point for the proposed BP-based HMMs structure. Figure 6a shows numerically simulated optical transmission as a function of the angle φ between the incidence plane (p-z plane) and y-z plane, and incident angle θ when working wavelength λ is 4, 4.5, 5, and 6 μm, respectively. The FWHM for narrowest angular optical transparency window is approximately 1.3° in Figure 6a. Obviously, the angle φ, which is corresponding to the PTT point, is increasing from 30° to 56° as λ increases from 4 to 6 μm when the other parameters unchanged. Figure 6b,c show Re(εp) and Im(εp) as a function of λ and angle φ based on Equation (5), respectively. Here, the purple dots represent the position of the PTT point obtained by simulation calculation. The cyan solid line in the Figure 6b represents Re (εp) = 0, and the cyan dashed line in the Figure 6c represents Im (εp) = 0.1. These purple dots obtained by simulation satisfy Re (εp) ≈ 0 and Im (εp) < 0.1. Thus, the PTT point obtained by simulation calculation is consistent with that obtained by theoretical calculation. Finally, we discuss the influence of working wavelength λ on the position of PTT point for the proposed BP-based HMMs structure. Figure 6a shows numerically simulated optical transmission as a function of the angle ϕ between the incidence plane (p-z plane) and y-z plane, and incident angle θ when working wavelength λ is 4, 4.5, 5, and 6 µm, respectively. The FWHM for narrowest angular optical transparency window is approximately 1.3 • in Figure 6a. Obviously, the angle ϕ, which is corresponding to the PTT point, is increasing from 30 • to 56 • as λ increases from 4 to 6 µm when the other parameters unchanged. Figure 6b,c show Re(ε p ) and Im(ε p ) as a function of λ and angle ϕ based on Equation (5) Re (ε p ) ≈ 0 and Im (ε p ) < 0.1. Thus, the PTT point obtained by simulation calculation is consistent with that obtained by theoretical calculation. Nanomaterials 2020, 10, x FOR PEER REVIEW 8 of 11 Figure 6. (a) Numerically simulated optical transmission as a function of the angle φ between the incidence plane (p-z plane) and y-z plane, and incident angle θ when λ = 4, 4.5, 5, and 6 μm, respectively. (b) Re(εp) and (c) Im(εp) as function of λ and φ.

Conclusions
In summary, anisotropic PTT is theoretically and numerically investigated based on the proposed BP-based HMMs structure consisting of alternating BP/SiO2 multilayer. By tailoring the IFC of BP-based HMMs from open hyperboloid to closed ellipsoid, both the theoretical calculations and numerical simulations show that a narrow angular transparency window appears at PTT point for p-polarized light. Moreover, we find that the angle φ, at which PTT appears, can be affected by working wavelength λ, thickness td of SiO2, or electron doping ρ of BP. It is believed that our work provides a new way in various angle-dependent optical applications, such as privacy protection and detectors with ultra-high signal-to-noise ratios.

Conclusions
In summary, anisotropic PTT is theoretically and numerically investigated based on the proposed BP-based HMMs structure consisting of alternating BP/SiO 2 multilayer. By tailoring the IFC of BP-based HMMs from open hyperboloid to closed ellipsoid, both the theoretical calculations and numerical simulations show that a narrow angular transparency window appears at PTT point for p-polarized light. Moreover, we find that the angle ϕ, at which PTT appears, can be affected by working wavelength λ, thickness t d of SiO 2 , or electron doping ρ of BP. It is believed that our work provides a new way in various angle-dependent optical applications, such as privacy protection and detectors with ultra-high signal-to-noise ratios.