## Abstract

In-plane photonic spin splitting effect is investigated in tunneling terahertz waves through an epsilon-near-zero metamaterial sandwiched between monolayer black phosphorus (BP). The strong in-plane anisotropy of BP layers will induce in-plane asymmetric spin splitting. The asymmetric spin splitting can be flexibly tuned by the angles between the incident plane and the armchair crystalline directions of the top and bottom BP layers, i.e., *ϕ*_{1} and *ϕ*_{2}. Based on this, an angle-resolved barcode-encryption scheme is discussed. For the special case of *ϕ*_{1} = *ϕ*_{2} = 0 or 90°, the transmitted beam undergoes Goos-Hänchen shift, which varies with the carrier density of BP. We believe these findings can facilitate the development of novel optoelectronic devices in the Terahertz region.

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

## 1. Introduction

Photonic spin splitting, refer to the spatial separation of two opposite spin components of the reflected/transmitted beam [1–4], has attracted significant attention owing to its application in precise metrology and quantum information [5–7]. The spin splitting can occur in directions both parallel and perpendicular to the plane of incidence, i.e., so-called the in-plane and out-of-plane spin splitting (IPSS and OPSS) [4–6]. Both IPSS and OPSS can be considered as results of spin-orbit coupling [1,3,4,8]. It has been demonstrated by Götte and Dennis in 2012 that the in-plane and out-of-plane spin splitting can be considered as analogous but reverse effects [9]. Different from Goos-Hänchen (GH) shift occurring in total reflection, the IPSS can appear in both cases of partial and total reflections [3,4]. In 2017, the upper limit of IPSS was derived, equal to the incident beam waist [4]. This upper limit can be obtained by optimizing the incident linear polarization state at Brewster angle [4]. Very recently, large IPSS near the critical angle of total reflection was theoretically proposed and experimentally verified [10].

Epsilon-near-zero (ENZ) metamaterial has significant interest due to its novel light-matter interactions [11]. According to the electromagnetic boundary condition, the vanishing epsilon in ENZ material will lead to strong discontinuity in the normal electric field, thus the strong field localization and enhancement [12]. The nonlinear and Purcell effects can therefore be boosted by ENZ metamaterials [13–15]. The vanishing epsilon of ENZ material also results in nearly constant phase advance when an electromagnetic field travels through the whole material [16]. Based on this phenomena, phase-mismatch-free harmonic generation, wavefront reshaping, and ultrafast phase transitions have been developed [17]. It has been demonstrated recently that ENZ metamaterials can enhance the GH shift [18] and the spin splitting [19,20].

To tune the spin splitting, we sandwich the ENZ metamaterial by monolayer black phosphorus. The light-matter interaction in this structure can be tuned by modulating the conductivity of BP, which is proportional to the electron carrier density [21]. BP possesses striking in-plane anisotropy property, enabling novel polarization-dependent and angle-resolved optoelectronic devices [22–26]. The asymmetric spin splitting induced by the in-plane anisotropy of BP was predicted very recently [27]. The reflected beam from a BP layer sitting on a silicon substrate undergoes both in-plane and out-of-plane asymmetric spin splitting near Brewster angle. However, the in-plane displacements of two opposite spin components are always positive, cannot be tuned flexibly [27]. The spin splitting of the transmitted beam from the structure are tiny small. The large asymmetric spin splitting in the transmitted beam is still missing [26,27].

Here, a novel structure combining the ENZ metamaterial and BP layers is proposed to realized large and tunable in-plane asymmetric spin splitting in the transmitted terahertz beam. The ENZ metamaterial acts as an optical cavity, resulting in an enhancement of light-matter interaction in BP layers. The BP layers can tune the asymmetric spin splitting via electron carrier density or angles between the incident plane and the armchair crystalline directions of BP layers.

## 2. Theory

Figure 1(a) shows the schematic of in-plane asymmetric spin splitting. An ENZ metamaterial is sandwiched by BP layers on both its top and bottom sides. The photonic property of BP layers can be described by surface conductivities. Under the Drude model, the conductivity along the armchair and zigzag crystalline directions are respectively [27,28]

Here,*D*

_{arm,zig}= πe^{2}

*ρ*/

*m*is the Drude weight with

_{arm,zig}*m*=

_{arm}*ћ*

^{2}/(2

*γ*

^{2}/Δ

*+ η*),

_{c}*m*=

_{zig}*ћ*

^{2}/2

*v*being the electron mass along the armchair or zigzag directions respectively. Parameters

_{c}*η =*10 meV,

*γ =*4

*a*/

*π*eVm, Δ

*=*2eV,

*η*

_{c}= ћ^{2}/0.4

*m*,

_{0}*v*

_{c}= ћ^{2}/1.4

*m*

_{0}. The scale length of BP

*a*= 0.223 nm. The electron carrier density

*ρ*can be changed by electric doping via a bias voltage [21–23].

Owing to the in-plane anisotropy of BP layer, the transmission varies when the incident plane changes with respect to the BP structure. A rotation angle *ϕ* is introduced, which is the angle between the incident plane and the armchair axis of BP crystal. The rotation angle can be changed practically by rotating the device with fixed incident beams. As shown by [24], the crystalline axes can be determined by analyzing the reflection of normal incident visible light with a polarizer and a CCD. The armchair and zigzag axes of BP are associated with maximum and minimum brightness of R channel, respectively [24].

For a rotation angle of *ϕ*, the conductance matrix connecting the surface current and electric light field can be given by $\overleftrightarrow{\sigma}$ = [*σ _{xx}*,

*σ*;

_{xy}*σ*,

_{yx}*σ*], where

_{yy}*σ*=

_{xx}*σ*cos

_{arm}^{2}

*ϕ + σ*sin

_{zig}^{2}

*ϕ*,

*σ*=

_{yy}*σ*sin

_{arm}^{2}

*ϕ + σ*cos

_{zig}^{2}

*ϕ*, and

*σ*=

_{yx}*σ*= (

_{xy}*σ*-

_{zig}*σ*)sin

_{arm}*ϕ*cos

*ϕ*[27]. The cross conductivity

*σ*results from the in-plane anisotropic of BP, which vanishes for isotropic 2D materials such as graphene.

_{yx}Transfer Matrix method is a powerful tool in the analysis of light propagation through layered media [29,30]. For dielectric media, the propagation of *s* and *p* waves are described by two 2 × 2 matrices independently. However, the cross conductivity of BP will induce a coupling between *p* and *s* waves. Therefore, a 4 × 4 matrix should be employed to describe the coupling.

Consider the propagation of light across an interface formed by a BP layer that separates two dielectrics with refractive indexes *n _{j}* and

*n*

_{j}_{+2}, as shown by Fig. 1(b). The electric field in media

*j*and

*j*+ 1 are respectively

*A*and

_{l}*C*are the amplitudes of rightward waves, while

_{l}*B*and

_{l}*D*are of leftward waves.

_{l}*k*is the

_{x}*x*component of wavevector, which is the same in different media. ${k}_{lz}={\left[{k}_{l}^{2}-{k}_{x}^{2}\right]}^{1/2}$ with

*k*being the wavenumber in is the

_{l}*l*medium. The electromagnetic boundary conditions at the interface between media

*j*and

*j*+ 1 are where

**H**

_{j}_{,}

_{j}_{+1}are the magnetic field in media

*j*and

*j*+ 1. By substituting Eq. (2) to Eq. (3), we have

*ω*and

*μ*are circular frequency and permeability, respectively. Therefore, the transmission matrix for a dielectric interface containing BP layer is

*σ*and

_{xy}*σ*vanish, the upper-right and lower-left parts of the transmission matrix vanish, thus the 4 × 4 transmission matrix can be rewritten into two 2 × 2 matrices, which connect

_{yx}*A*,

_{j}*B*with

_{j}*A*

_{j}_{+1},

*B*

_{j}_{+1}and

*C*,

_{j}*D*with

_{j}*C*

_{j}_{+1},

*D*

_{j}_{+1}, respectively.

As the photonic property of BP layers have been described by surface conductivities, the thicknesses of BP layers are neglected. The propagation matrix in a homogenous medium *j* with a thickness of *d _{j}* is [29]

*N*BP layers shown in Fig. 1(c), the transfer matrix can be obtained by transmission and propagation matrices. Denote the electric field coefficients on left side of the first interface by

*A*

_{1},

*B*

_{1},

*C*

_{1},

*D*

_{1}and those on the right side of the end interface by

*A*

_{N+}_{1},

*B*

_{N+}_{1},

*C*

_{N+}_{1},

*D*

_{N+}_{1}. These two sets of field coefficients are then related by a 4 × 4 transfer matrix

*M*, namely,

*M*=

*T*

_{01}

*P*

_{1}

*T*

_{12}

*P*

_{2}…

*T*

_{N}_{-1,}

_{N}P_{N}T_{N}_{,}

_{N}_{+1}. Equation (7) is the transfer matrix to describe the light propagation through layered media with BP. For a multilayer structure without BP, both the transmission and propagation matrices can be rewritten respectively into two 2 × 2 matrices, the final matrix

*M*can be therefore reduced into two 2 × 2 independent matrices [30].

According to Eq. (7), the Fresnel transmission coefficients can be obtained. By setting *C*_{1} = *B*_{N + 1} = *D*_{N + 1} = 0, we have

*A*

_{1}=

*B*

_{N + 1}=

*D*

_{N + 1}= 0, we have

In our case, two BP layers are considered. The electron carrier densities and rotation angles (Angle between the armchair axis of BP crystal and the *x _{g}*-axis) of top and bottom BP layers are respectively

*ρ*

_{1},

*ϕ*

_{1}and

*ρ*

_{2},

*ϕ*

_{2}, as shown in Fig. 1(a). The carrier densities are set to be the same, i.e.,

*ρ*

_{1}=

*ρ*

_{2}=

*ρ*, through the whole paper.

Considering a horizontally (*H*) polarized Gaussian incident beam tunneling through the BP multi-layer structure. The spectrum of incident beam is${\tilde{E}}_{i}=\mathrm{exp}\left[-({k}_{x}^{2}+{k}_{y}^{2}){w}_{0}^{2}/4\right]|H\u3009$ where *w*_{0} being the beam waist. According to [31,32], the transmitted spectrum can be given by ${\tilde{E}}_{t}=Q{\tilde{E}}_{i}$. The transformation matrix **Q** can be expressed as

*t*+

_{sp}*t*)cot

_{ps}*θ*, K = (

*t*-

_{pp}*t*)cot

_{ss}*θ*, ${t}_{ij}^{\text{'}}$are the first derivate of

*t*.

_{ij}*κ*

_{x}_{,}

*=*

_{y}*k*,

_{x}_{y}/

*k*

_{0}with

*k*

_{0}being the wavenumber in vacuum. In the circular polarization basis, the right- and left-handed circular polarization (RCP and LCP) components of the transmitted beams for

*H*incident polarization are:

*x*-axis. With respect to the geometric prediction, the centroid displacements can be defined as ${X}_{\pm}={\displaystyle \iint {\tilde{E}}_{t}^{\pm}{\partial}_{{k}_{rx}}{\tilde{E}}_{t}^{\pm}|}d{k}_{x}d{k}_{y}/{\displaystyle \iint |{\tilde{E}}_{t}^{\pm}{|}^{2}}d{k}_{x}d{k}_{y}$ [32]. After some straightforward calculation, we obtain

*t*. Generally, the cooperation effect of the first and second terms will induce an asymmetric spin splitting.

_{sp}## 3. Results and discussion

Firstly, we consider the case with the rotation angles of top and bottom BP layers being *ϕ*_{1} = *ϕ*_{2} = 0. Thus, the cross conductivities are zero, resulting in the vanishing of the spin-dependent term in Eq. (11). Therefore, the transmitted beam undergoes pure GH shift without spin splitting. Figure 2(a) shows the GH shifts changing with the incident angle *θ* for different carrier density *ρ* when the refractive index and thickness of ENZ metamaterial are *n* = 0.1 and *d* = 0.1 mm, respectively. Here, we focus our attention on the terahertz region by setting the frequency incident beam being 1 THz, since this region is underdeveloped but ripe for exploitation [33]. For each carrier density, the GH shift reaches a peak at *θ* = 6.12°. The peak GH shift increases with the carrier density, which is clearly shown by the blue line in Fig. 2(c). Thus, the GH shift can be enhanced by the BP layers, since *ρ* = 0 corresponds to the case without BP. At other incident angle such as *θ* = 3.5°, the GH shift may decreases with the carrier density.

The GH shift depends on the refractive index *n* and thickness *d* of the ENZ metamaterial. As shown by Fig. 2(b), with the increase of *n*, the peak value of GH shift decreases and the peak position moves, indicating that a larger shift can be obtain by a smaller refractive index of the metamaterial. Figure 2(d) gives the GH shift increases monotonously with thickness *d*, when *ρ* = 5 × 10^{17}m^{−2}, *θ* = 6.12°.

By rotating the BP layer, the transmitted beam undergoes asymmetric spin splitting. The displacements of two opposite spin components of the transmitted beam are given in Fig. 3(a) as function of the incident angle when *ϕ*_{1} = *ϕ*_{2} = 45°, *d* = 4.2 μm, *n* = 0.1. In this special case, the displacements of RCP component *X*_{+} are almost vanish for all carrier densities. However, the displacements of LCP component *X*_{-} can take large values. When the carrier density increases from 2 × 10^{17} to 5 × 10^{17}m^{−2}, the peak displacement decreases gradually, and change its sign at *ρ* = 4.7 × 10^{17}m^{−2}. Meanwhile, the peak position moves from 71.5° to 50.5° gradually. At the peak positions, owing to the strong light-interaction in BP layers, the transmission coefficients |*t _{ps}*| are almost identical with coefficients |

*t*| (see Fig. 3(b)). And they are with nearly −90° phase difference (

_{pp}*φ*= arg[

*t*/|

_{pp}*t*]), as shown by Fig. 3(b). Therefore, the first three terms of

_{ps}*W*

_{-}in Eq. (12) almost vanishes, leading to large displacement in

*X*

_{-}. The phenomenon was also found at air-chiral metamaterial interface [34].

The refractive index and the thickness of the ENZ metamaterial affect the asymmetric spin splitting in the transmitted beam. Figures 4(a) and 4(b) show the displacements *X* _{±} changing with the incident angle *θ* for the refractive index *n* = 0.05, 0.1. 0.2 (a), and for *n* = 0.1, 0.1 + 0.01*i*, 0.1 + 0.1*i*(b), respectively. One can find from Fig. 4(a) that, the maximum value of displacements |*X* _{±} | as well as the maximum displacement difference |*X*_{+}-*X*_{-}| decreases with the increase of the real part of the refractive index Re[*n*]. Therefore, a near-zero permittivity can enhance the asymmetric spin splitting. The near-zero permittivity can be achieved by both artificial and natural materials. Polar dielectrics like LiF, NaCl, and GaAs are the natural terahertz ENZ materials with the real part of permittivity close to zero near the polaritonic resonance frequency [35]. Effective terahertz ENZ materials have been realized by rectangular waveguide [36], graphene-dielectric composite structure [37], etc.

The loss of the metamaterial can be described by the image part of the refractive index Im[*n*] [18]. As shown by Fig. 4(a), Im[*n*] will reduce the displacements *X* _{±} . However, for the Im[*n*] smaller than 0.01, the displacements *X* _{±} are almost unchanged. Figure 4(c) shows the displacements *X* _{±} changing with the metamaterial thickness *d* when *n* = 0.1. The displacements *X* _{±} will approach to asymptotic values when the thickness *d* increases. For a larger incident angle, both *X*_{+} and *X*_{-} possess a sharper peak.

The asymmetric spin splitting in the transmitted beam depends strongly on the rotation angle of BP layers, as shown in Fig. 5. In Fig. 5, the rotation angles of the top and bottom layers are equal, i.e., *ϕ* = *ϕ*_{1} = *ϕ*_{2}. The displacements *X*_{+} and *X*_{-} are mirror symmetric about *ϕ* = 90°. The displacements of the RCP component *X*_{+} takes large values *ϕ*>90°, while *X*_{-} for *ϕ*<90° . When *θ* = 44°, *X*_{+} has a positive peak, and become two positive peaks for *θ* = 47°. When *θ*>50°, a positive and a negative peaks can be found in the pattern of *X*_{+}. And distance between this two peaks enlarges with the increase of the incident angle *θ*.

In the above case, one spin component has large displacement while the other almost vanishes. Figure 6(a) gives the displacements *X* _{±} changing with the incident angle for the thickness *d* = 0.15 mm the carrier density *ρ* = 0 and 5 × 10^{17}m^{−2}, respectively. Without BP layers (*ρ* = 0), *X*_{+} and *X*_{-} coincide each other. When *ρ* = 5 × 10^{17}m^{−2}, however, *X*_{+} and *X*_{-} are different excepting for *ϕ* = 0, 90°, or 180°. At these three angles, the two opposite spin components of the transmitted beam do not split. However, the two spin components may split if the incident beam is elliptically polarized [4]. At *θ* = 6.21°, both *X*_{+} and *X*_{-} are positive. However, *X*_{+} and *X*_{-} are opposite in sign at *θ* = 11.5°. Figures 6(b) and 6(c) show the displacements *X* _{±} changing with the rotation angle for *θ* = 6.21° and 11.5°, respectively. The displacements vary with rotation angle flexibly.

The asymmetric displacements for LCP and RCP components of the transmitted beam along *x*-axis can be regarded as two independent channels for information processing [27]. The centroid displacements larger than a certain threshold value can be considered as “1”, while smaller than that value as “0”. The encoding rule varies with the choice of the threshold values. For the case of *θ* = 11.5°, three different codes: “11”, “10”, “01” can be obtained by tuning the rotation angle with a threshold of 0. All four different codes: “11”, “10”, “01”, “00” are found for the case of *θ* = 6.21° with a threshold of 1.8 mm. This is comparing to best three different codes in [27] and [38].

Twisted bilayer graphene have attracted much attention recently owing to its unconventional superconductivity at the magic-angle [39]. When the rotation angles of the top and bottom BP layers are different, the transmitted beam shows interesting spin splitting phenomena. Figure 7 shows the displacements *X* _{±} as functions of *ϕ*_{1} and *ϕ*_{2} when *θ* = 14.14° and *d* = 0.15 mm. *X*_{+} and *X*_{-} are almost symmetric about the centroid point of (*ϕ*_{1},*ϕ*_{2}) = (90°,90°). Both *X*_{+} and *X*_{-} will change sign when *ϕ*_{2} crosses 90°. With the change of *ϕ*_{1}, only the magnitude of *X* _{±} will change.

Figure 8 shows Displacements *X* _{±} changing with the rotation angle *ϕ*_{2} for different *ϕ*_{1} when *θ* = 10° (a), 20° (b),30° (c), respectively. When *θ* = 10°, there is one peaks for *X* _{±} in the *ϕ*_{2} region of 0-90°. When *θ* = 20° and 30°, the *X*_{-} become two peaks in the same region, while the *X*_{+} are nearly flat.

The ENZ metamaterials sandwiched with BP layers show unique photonic spin splitting effect. Although the CVD growth of 2D black phosphorus film has been demonstrated in 2016, the large-scale of monolayer BP is still challenge [40]. Hopeful methods have been proposed recently to address this challenge [41]. The unique photonic spin splitting effect based on BP can be readily extended to other anisotropic 2D materials such as ReSe_{2} and ReS_{2}, for which high quality large-scale atomic films are well reachable [42].

## 4. Conclusions

In conclusion, we have extended the transformation matrix to describe the propagation in a multi-layer dielectric stack containing BP layers. A novel air-BP-metamaterial-BP-air structure is proposed to enhance the in-plane photonic spin splitting effect. The transmitted beam through the structure undergoes asymmetric spin splitting due to the strong in-plane anisotropy of BP layers. The asymmetric splitting depends strongly on the rotation angles of BP layers, namely, *ϕ*_{1} and *ϕ*_{2.} However, only the GH shift occurs without spin splitting for *ϕ*_{1} = *ϕ*_{2} = 0 or 90°. Based on the tunable asymmetric spin splitting, a two-channel angle-resolved barcode encoding is achieved. The photonic spin splitting effect in BP layers promises novel angle-resolved optoelectronic devices in the Terahertz region.

## Funding

National Natural Science Foundation of China (61705086, 61675092, 61475066); Natural Science Foundation of Guangdong Province (2017A030313375, 2016TQ03X962, 2017A010102006, 2016A030311019, 2016A030313079, 2017A030313359); Science & Technology Project of Guangzhou (201803020023, 201707010396, 201704030105, 201605030002, 201604040005).

## References

**1. **O. Hosten and P. Kwiat, “Observation of the spin hall effect of light via weak measurements,” Science **319**(5864), 787–790 (2008). [CrossRef] [PubMed]

**2. **K. Y. Bliokh, C. T. Samlan, C. Prajapati, G. Puentes, N. K. Viswanathan, and F. Nori, “Spin - Hall effect and circular birefringence of a uniaxial crystal plate,” optia **3**, 1039 (2016). [CrossRef]

**3. **X. Qiu, Z. Zhang, L. Xie, J. Qiu, F. Gao, and J. Du, “Incident-polarization-sensitive and large in-plane-photonic-spin-splitting at the Brewster angle,” Opt. Lett. **40**(6), 1018–1021 (2015). [CrossRef] [PubMed]

**4. **W. Zhu, J. Yu, H. Guan, H. Lu, J. Tang, J. Zhang, Y. Luo, and Z. Chen, “The upper limit of the in-plane spin splitting of Gaussian beam reflected from a glass-air interface,” Sci. Rep. **7**(1), 1150 (2017). [CrossRef] [PubMed]

**5. **X. Zhou, L. Sheng, and X. Ling, “Photonic spin Hall effect enabled refractive index sensor using weak measurements,” Sci. Rep. **8**(1), 1221 (2018). [CrossRef] [PubMed]

**6. **X. Zhou, Z. Xiao, H. Luo, S. Wen, Z. Xiao, Z. H. Luo, and S. Wen, “Experimental observation of the spin Hall effect of light on a nanometal film via weak measurements,” Phys. Rev. A **85**(4), 043809 (2012). [CrossRef]

**7. **X. Zhou, X. Ling, H. Luo, and S. Wen, “Identifying graphene layers via spin Hall effect of light,” Appl. Phys. Lett. **101**(25), 251602 (2012). [CrossRef]

**8. **K. Y. Bliokh and Y. P. Bliokh, “Conservation of angular momentum, transverse shift, and spin Hall effect in reflection and refraction of an electromagnetic wave packet,” Phys. Rev. Lett. **96**(7), 073903 (2006). [CrossRef] [PubMed]

**9. **J. B. Götte and M. R. Dennis, “Generalized shifts and weak values for polarization components of reflected light beams,” New J. Phys. **14**(7), 073016 (2012). [CrossRef]

**10. **X. Zhou, L. Xie, X. Ling, S. Cheng, Z. Zhang, H. Luo, and H. Sun, “Large in-plane asymmetric spin angular shifts of a light beam near the critical angle,” Opt. Lett. **44**(2), 207–210 (2019). [CrossRef] [PubMed]

**11. **A. Alù, M. G. Silveirinha, A. Salandrino, and N. Engheta, “Epsilon-near-zero metamaterials and lectromagnetic sources: Tailoring the radiation phase pattern,” Phys. Rev. B Condens. Matter Mater. Phys. **75**(15), 155410 (2007). [CrossRef]

**12. **I. Liberal and N. Engheta, “Near-zero refractive index photonics,” Nat. Photonics **11**(3), 149–158 (2017). [CrossRef]

**13. **L. H. Nicholls, F. J. Rodríguez-Fortuño, M. E. Nasir, R. M. Córdova-Castro, N. Olivier, G. A. Wurtz, and A. V. Zayats, “Ultrafast synthesis and switching of light polarization in nonlinear anisotropic metamaterials,” Nat. Photonics **11**(10), 628–633 (2017). [CrossRef]

**14. **M. Z. Alam, I. De Leon, and R. W. Boyd, “Large optical nonlinearity of indium tin oxide in its epsilon-near-zero region,” Science **352**(6287), 795–797 (2016). [CrossRef] [PubMed]

**15. **A. V. Chebykin, A. A. Orlov, A. S. Shalin, A. N. Poddubny, and P. A. Belov, “Strong Purcell effect in anisotropic *ε*-near-zero metamaterials,” Phys. Rev. B Condens. Matter Mater. Phys. **91**(20), 205126 (2015). [CrossRef]

**16. **I. Liberal, A. M. Mahmoud, Y. Li, B. Edwards, and N. Engheta, “Photonic doping of epsilon-near-zero media,” Science **355**(6329), 1058–1062 (2017). [CrossRef] [PubMed]

**17. **X. Niu, X. Hu, S. Chu, and Q. Gong, “Epsilon-Near-Zero Photonics: A New Platform for Integrated Devices,” Adv. Opt. Mater. **6**(10), 1701292 (2018). [CrossRef]

**18. **Y. Xu, C. T. Chan, and H. Chen, “Goos-Hänchen effect in epsilon-near-zero metamaterials,” Sci. Rep. **5**(1), 8681 (2015). [CrossRef] [PubMed]

**19. **W. Zhu and W. She, “Enhanced spin Hall effect of transmitted light through a thin epsilon-near-zero slab,” Opt. Lett. **40**(13), 2961–2964 (2015). [CrossRef] [PubMed]

**20. **M. Jiang, W. Zhu, H. Guan, J. Yu, H. Lu, J. Tan, J. Zhang, and Z. Chen, “Giant spin splitting induced by orbital angular momentum in an epsilon-near-zero metamaterial slab,” Opt. Lett. **42**(17), 3259–3262 (2017). [CrossRef] [PubMed]

**21. **X. Wang and S. Lan, “Optical properties of black phosphorus,” Adv. Opt. Photonics **8**(4), 618 (2016). [CrossRef]

**22. **R. Peng, K. Khaliji, N. Youngblood, R. Grassi, T. Low, and M. Li, “Midinfrared Electro-optic Modulation in Few-Layer Black Phosphorus,” Nano Lett. **17**(10), 6315–6320 (2017). [CrossRef] [PubMed]

**23. **B. Deng, R. Frisenda, C. Li, X. Chen, A. Castellanos-Gomez, and F. Xia, “Progress on Black Phosphorus Photonics,” Adv. Opt. Mater. **6**(19), 1800365 (2018). [CrossRef]

**24. **N. Mao, J. Tang, L. Xie, J. Wu, B. Han, J. Lin, S. Deng, W. Ji, H. Xu, K. Liu, L. Tong, and J. Zhang, “Optical Anisotropy of Black Phosphorus in the Visible Regime,” J. Am. Chem. Soc. **138**(1), 300–305 (2016). [CrossRef] [PubMed]

**25. **Y. Wang, F. Zhang, X. Tang, X. Chen, Y. Chen, W. Huang, Z. Liang, L. Wu, Y. Ge, Y. Song, J. Liu, D. Zhang, J. Li, and H. Zhang, “All-Optical Phosphorene Phase Modulator with Enhanced Stability Under Ambient Conditions,” Laser Photonics Rev. **12**(6), 1800016 (2018). [CrossRef]

**26. **W. Shu, J. Zhang, H. Luo, S. Chen, W. Zhang, W. Wu, S. Wen, and X. Ling, “Photonic spin Hall effect on the surface of anisotropic two-dimensional atomic crystals,” Photon. Res. **6**(6), 511 (2018). [CrossRef]

**27. **H. Lin, B. Chen, S. Yang, W. Zhu, J. Yu, H. Guan, H. Lu, Y. Luo, and Z. Chen, “Photonic spin Hall effect of monolayer black phosphorus in the Terahertz region,” Nanophotonics **7**(12), 1929–1937 (2018). [CrossRef]

**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. **A. Yariv and P. Yeh, *Photonics: Optical Electronics in Modern Communications*, 6th ed. New York, (Oxford University, 2007).

**30. **T. Zhan, X. Shi, Y. Dai, X. Liu, and J. Zi, “Transfer matrix method for optics in graphene layers,” J. Phys. Condens. Matter **25**(21), 215301 (2013). [CrossRef] [PubMed]

**31. **M. Liu, L. Cai, S. Chen, Y. Liu, H. Luo, and S. Wen, “Strong spin-orbit interaction of light on the surface of atomically thin crystals,” Phys. Rev. A (Coll. Park) **95**(6), 063827 (2017). [CrossRef]

**32. **K. Y. Bliokh and A. Aiello, “Goos–Hänchen and Imbert–Fedorov beam shifts: an overview,” J. Opt. **15**(1), 014001 (2013). [CrossRef]

**33. **X. Zhang, A. Shkurinov, and Y. Zhang, “Extreme terahertz science,” Nat. Photonics **11**(1), 16–18 (2017). [CrossRef]

**34. **M. Jiang, H. Lin, L. Zhuo, W. Zhu, H. Guan, J. Yu, H. Lu, J. Tan, and Z. Chen, “Chirality induced asymmetric spin splitting of light beams reflected from an air-chiral interface,” Opt. Express **26**(6), 6593–6601 (2018). [CrossRef] [PubMed]

**35. **M. Kafesaki, A. A. Basharin, E. N. Economou, and C. M. Soukoulis, “THz metamaterials made of phonon-polariton materials,” Photon. Nanostructures **12**(4), 376–386 (2014). [CrossRef]

**36. **V. Pacheco-Peña, N. Engheta, S. Kuznetsov, A. Gentselev, and M. Beruete, “Experimental Realization of an Epsilon-Near-Zero Graded-Index Metalens at Terahertz Frequencies,” Phys. Rev. Appl. **8**(3), 034036 (2017). [CrossRef]

**37. **M. A. K. Othman, C. Guclu, and F. Capolino, “Graphene-Dielectric Composite Metamaterials: Evolution from Elliptic to Hyperbolic Wavevector Dispersion and The Transverse Epsilon-Near-Zero Condition,” Nanophotonics **7**(1), 073089 (2013). [CrossRef]

**38. **T. Tang, J. Li, L. Luo, P. Sun, and J. Yao, “Magneto-Optical Modulation of Photonic Spin Hall Effect of Graphene in Terahertz Region,” Adv. Opt. Mater. **6**(7), 201701212 (2018). [CrossRef]

**39. **Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, “Correlated Insulator Behaviour at Half-Filling in Magic-Angle Graphene Superlattices,” Nature **556**(7699), 80–84 (2018). [CrossRef] [PubMed]

**40. **J. B. Smith, D. Hagaman, and H. F. Ji, “Growth of 2D black phosphorus film from chemical vapor deposition,” Nanotechnology **27**(21), 215602 (2016). [CrossRef] [PubMed]

**41. **T. Niu, “New properties with old materials: Layered black phosphorous,” Nano Today **12**, 7–9 (2017). [CrossRef]

**42. **A. Dathbun, Y. Kim, S. Kim, Y. Yoo, M. S. Kang, C. Lee, and J. H. Cho, “Large-Area CVD-Grown Sub-2 V ReS2 Transistors and Logic Gates,” Nano Lett. **17**(5), 2999–3005 (2017). [CrossRef] [PubMed]