Generation of Anisotropic Massless Dirac Fermions and Asymmetric

Feb 23, 2017 - ... than the interatomic distances would create a superlattice structure, mixing the ... by a generalized Weyl equation(35) (in general...
0 downloads 0 Views 2MB Size
Letter pubs.acs.org/NanoLett

Generation of Anisotropic Massless Dirac Fermions and Asymmetric Klein Tunneling in Few-Layer Black Phosphorus Superlattices Zhenglu Li, Ting Cao, Meng Wu, and Steven G. Louie* Department of Physics, University of California at Berkeley, Berkeley, California 94720, United States Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, United States S Supporting Information *

ABSTRACT: Artificial lattices have been employed in a broad range of two-dimensional systems, including those with electrons, atoms, and photons, in the quest for massless Dirac fermions with high flexibility and controllability. Establishing triangular or hexagonal symmetry, from periodically patterned molecule assembly or electrostatic gating as well as from moiré pattern induced by substrate, has produced electronic states with linear dispersions from two-dimensional electron gas (2DEG) residing in semiconductors, metals, and graphene. Different from the commonly studied isotropic host systems, here we demonstrate that massless Dirac fermions with tunable anisotropic characteristics can, in general, be generated in highly anisotropic 2DEG under slowly varying external periodic potentials. In the case of patterned few-layer black phosphorus superlattices, the new chiral quasiparticles exist exclusively in certain isolated energy window and inherit the strong anisotropic properties of pristine black phosphorus. These states exhibit asymmetric Klein tunneling, in which the transmission probability of the wave packets with normal incidence is no longer unity and can be tuned and controlled. In general, the direction of wave packet incidence for perfect transmission and that of the normal incidence are different, and the difference can reach more than 50° under an appropriate barrier orientation in black phosphorus superlattices. Our findings provide insight into the understanding and possible utilization of these novel emergent chiral quasiparticles. KEYWORDS: Black phosphorus, anisotropic 2D electron gas, two-dimensional superlattice, anisotropic Dirac fermions, asymmetric Klein tunneling

T

questions: (1) can these novel quasiparticles arise from highly anisotropic host states, and (2) what new physical insights and phenomena one can obtain from such systems? In this work, we show that massless Dirac fermions with a chiral character can in general be generated from anisotropic 2DEG systems. We moreover propose that semiconducting few-layer black phosphorus23,24 is an ideal system to achieve this aim upon electron or hole doping.25,26 Electron and hole carriers in black phosphorus exhibit strong anisotropy with their effective masses along the two orthogonal crystal axes (m*x and my*) differing by an order of magnitude,27−29 forming highly anisotropic two-dimensional systems. With properly designed periodic potential (that can be realized with patterned molecule assembly or electrostatic gating under laboratory conditions), we predict that highly anisotropic massless Dirac fermions are generated and that the ratio of the group velocities along the two crystalline directions reaches (mx*/my*)1/2 ∼ 1/3. These anisotropic quasiparticles moreover exist within an

he unusual relativistic-like quasiparticles in graphene, known as massless Dirac fermions, intrinsically come from the symmetry of its honeycomb crystal structure1,2 and therefore can be designed and manipulated in other systems as first theoretically predicted.3,4 Many experimental efforts have been successfully performed to generate Dirac fermions in common isotropic two-dimensional electron gas (2DEG) on metal or in semiconductor quantum wells under slowly varying periodic potentials formed by molecule assembly or patterned gate,5−10 as well as with ultracold atoms trapped in honeycomb optical crystal structures11,12 and photons confined in photonic crystals.13,14 Similarly, lattice mismatch between graphene and substrate can introduce a long-range moiré superlattice potential, giving rise to the new generation of Dirac points in graphene.15−22 These platforms combine advantages of an extended degree of control and various foreign modifications,8 making much physics accessible such as the interplay between Dirac fermions with strong correlation6 and spin−orbit coupling,9 and the Hofstadter butterfly effect.17−21 A simple symmetry argument prescribes a trigonally modulating potential to generate massless Dirac fermions from an isotropic system.3,8 However, previous studies do not address the © XXXX American Chemical Society

Received: November 28, 2016 Revised: February 11, 2017 Published: February 23, 2017 A

DOI: 10.1021/acs.nanolett.6b04942 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters H2DEG(p) =

py2 px2 + 2mx* 2m*y

(1)

where px and py are the crystal momenta along the x and y directions, respectively. Equation 1 is valid around both the conduction band minimum and valence band maximum, where the band disperses much faster along the x direction (armchair direction) than along the y direction (zigzag direction) (Figure 1a and b). An external periodic potential varying much slower than the interatomic distances would create a superlattice structure, mixing the states in the Γ valley. For an isotropic 2DEG, with m*x = m*y , a triangular potential would strongly mix states near three Kj (j = 1, 2, 3) points (with Kj = pj/ℏ, determined by the periodicity and orientation of the external potential) with degenerate energy forming a triangle in reciprocal space, leading to Dirac fermions.3,9 In an anisotropic 2DEG, mx* ≠ my* (for electron-doped few-layer black phosphorus, we take mx* = 0.15me and my* = 1.18me, where me is the free electron mass27); however, in general, one can still choose three Kj points (on which the states have the same energy) and mix these states with some scattering strength from a designed potential, that is, solving an inverse problem. As a demonstration, we pick up the following three Kj points, K1 = (K 0 , 0)

Figure 1. Few-layer black phosphorus. (a) Crystal structure and orientation. (b) Schematic band structure near the band gap around the Γ point. Both the lowest conduction band (blue curve) and the highest valence band (red curve) disperse much faster along the kx direction (armchair direction) than the ky direction (zigzag direction). (c) Geometry of a triangle defined by three points with degenerate energy in the reciprocal space, fitting into the band dispersion. (d) Elliptical isoenergetic contour. Three Kj points of degenerate energies on the energy surface define three Gj vectors. The reciprocal space of the superlattice is spanned by G1 and G2.

⎛ K 3 K 2 = ⎜− 0 , ⎜ 2 2 ⎝

⎛ K 3 K3 = ⎜ − 0 , − ⎜ 2 2 ⎝

⎞ m*y K 0⎟ mx* ⎟⎠

⎞ m*y K 0⎟ mx* ⎟⎠

(2)

where K0 is a free parameter scaling the periodicity of the external potential U(r). (For concrete illustration, we shall set the periodicity of U(r) as 2.35 nm in plotting the figures in the main text.) In generating eq 2, we demand the Kj points satisfy the condition (Figure 1c and d): E(ℏK1) = E(ℏK2) = E(ℏK3) where E(ℏKj) is the energy eigenvalue of H2DEG(ℏKj). With the degeneracy condition satisfied, the choice of the three Kj points still has many possibilities and may lead to rich phases described by a generalized Weyl equation35 (in general tilted Dirac cones or open Fermi surfaces), and even higher pseudospin fermions can be generated with more K j included.36,37 Here, we consider a special case with the three Kj points satisfying K1 + K2 + K3 = 0 (eq 2) which guarantees electron−hole symmetry in the Dirac cones generated later. The corresponding external potential which would lead to mixing of states near those three Kj points then has reciprocal lattice vectors (and superlattice Brillouin zone (BZ)) given by G1 = K2 − K1 and G2 = K3 − K1, while G3 = G1 − G2 is a dependent vector (Figure 1d). As a first example, we consider a simple external potential of the following form,3,9

isolated energy window separating from other states and hence are expected to be quite measurable and controllable. Unlike the Klein tunneling process in graphene with carriers normal incident upon a potential barrier always experiencing perfect transmission,30,31 we find that the anisotropic massless Dirac fermions generated in black phosphorus superlattices allow for an asymmetric Klein tunneling, in which the normal incident wave packets are no more unimpeded and can be generally tuned and controlled. Specifically, in the case of black phosphorus superlattices, the directions of normal incidence and perfect transmission can be differed by more than 50°. The anisotropic massless Dirac fermions and the asymmetric Klein tunneling can be manifested in various experiments including scanning tunneling microscope, electrical transport, quantum oscillations, quantum Hall, magnetoresistance, and so forth. Moreover, these novel tunable features could provide new ingredients in applications32,33 such as electron optics, where the system serves as anisotropic crystals and negative refraction and anomalous reflection of the chiral Dirac fermions are expected. Monolayer and few-layer black phosphorus is a direct gap semiconductor with a band gap of the order of 1 eV at the Brillouin zone center. We shall use this case as a generic example of an anisotropic massive 2DEG. Let us begin with the effective Hamiltonian for the Γ valley of a few-layer black phosphorus (or monolayer phosphorene),29,34

U (r) = 2W [cos(G1·r) + cos(G2 ·r) + cos(G3·r)]

(3)

where W defines the strength of the external potential. The real-space distribution of this potential is plotted in Figure 2a, with W = 0.1 eV. We define k as a small wavevector, i.e., |k| ≪ |Kj|, and expand the superlattice wave function near the three Kj points as a linear combination of |j⟩ = ei(Kj+k)·r. In this basis, the Hamiltonian up to the first order in k is given by H′2DEG = H0 + H1 with ⟨j|H0|j′⟩ = W (1 − δjj′) and ⟨j|H1|j′⟩ = ℏvj·kδjj′, where δjj′ is the Kronecker delta function and vj = (vjx,vjy) = B

DOI: 10.1021/acs.nanolett.6b04942 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters

Figure 2. Real-space distribution of a superlattice of sinusoidal potentials defined in eq 3 with W = 0.1 eV. a1 and a2 are the superlattice lattice vectors of the external periodic potential. |a1| = |a2| = 2.35 nm, and the angle between the two vectors is 63.4°. The two panels on the right are the line profiles of the external potential. (b) Band structures of the states under the sinusoidal periodic potential from nonperturbative numerical calculations. Anisotropic massless Dirac fermion are generated. (c) Band structures along two directions passing through one Dirac point, showing linearly dispersing features with strongly anisotropic group velocities. (d) DOS of the new Dirac fermion system generated from the external periodic potential, with the real electron spin degree of freedom included. Dashed line represents the DOS from the lowest-energy conduction band in pristine few-layer black phosphorus before applying the external potential.

⎛ ℏK jx ℏK jy ⎞ ⎜ ⎟. The eigenvalues of the states for k = 0 , ⎝ mx* m*y ⎠ (Hamiltonian H0) are −W, −W, 2W. The resulting two 1 1 degenerate states are 2 (0, 1, − 1)T and 6 (2, −1, −1)T , and

dispersion varies with the direction of k, resulting in an anisotropic Dirac cone band structure. The anisotropy in the |vy|

group velocity, |v | = γ0 , originates from the anisotropic effective x

masses intrinsic to the original anisotropic 2DEG. We thus have demonstrated the generation of anisotropic Dirac fermions in general, and in few-layer black phosphorus in particular, for a sinusoidal potential defined in eq 3 for small k using a perturbative analysis. We have also numerically solved in the whole superlattice BZ for the full band structure (Figure 2b and c) and find that the Dirac cone indeed disperses for a large part of the BZ much faster along the kx direction than along the ky direction. We now examine more complex shapes for the applied external potential within a unit cell of the superlattice. Considering possible experimental fabrications, we perform calculations on three likely cases: a conical potential of radius 0.3 nm with 2.35 nm periodicity mimicking molecule assembly,7,10 a rectangular potential of 2.67 nm long and 1.65 nm wide with 7.84 nm periodicity representing patterned electrostatic gating,3,5,9 and a harmonic potential of radius 0.4 nm with 2.35 nm periodicity. We expand all quantities with plane waves to construct the Hamiltonian. The diagonal matrix elements are the unperturbed band energies from eq 1, and the off-diagonal matrix elements are from the external potential, which is given in a Fourier expansion as U(r) = ∑GU(G)eiG·r. Numerical diagonalization of the Hamiltonian gives the energy and wave function of states in the whole BZ. Generally, if the potential has small deviations in the corresponding Fourier components given in eq 3, the generated Dirac points will be slightly shifted away from the three Kj,12 and the above analysis remains valid approximately. Numerical simulations (i.e.,

we may construct a 2 × 2 subspace with them.3,9 The Hamiltonian at finite k in the new basis (setting the zero of ℏv energy at −W) reads H′ = 20 ( −σzkx − γ0σxk y) where v0 = ℏK*0 , mx

γ0 =

mx* m*y

, and σj (j = x,y,z) are Pauli matrices. With a unitary

transform,3,9 we arrive at H(k) =

ℏv0 (σxkx + γ0σyk y) 2

(4)

The Hamiltonian H resembles the low-energy effective Hamiltonian in graphene,1,2 except that an extra factor γ0 renormalizes the second term. The energy eigenvalues and the corresponding eigenvectors are therefore given by E(k) =

λℏv0 2 kx + γ02k y2 2

⎛ 1 ⎞ ψλ(k) = ⎜ iϕ ⎟ ⎝ λe s ⎠

(5)

(6)

where λ = ±1 denoting the upper and lower bands, and ϕs is the polar angle of the pseudospin vector associated with ψλ(k) with respect to kx. Here we use the notation “s” for pseudospin. The energy dispersion is linear, and these quasiparticles behave as massless Dirac fermions under the modulation of the periodic external potential. The group velocity of the linear C

DOI: 10.1021/acs.nanolett.6b04942 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters

Figure 3. (a) Isoenergetic contour of the anisotropic Dirac cone. Wavevector k, pseudospinor s, and group velocity v are noncollinear in general. (b) With the elliptical Fermi surface plotted in the x−y coordinate system, a potential barrier is created with barrier normal along the x′-direction in the x′−y′ coordinate system. Vectors and angles are denoted by the prime symbol if in the x′−y′ coordinate system, and by r if representing reflected waves. (c) An n−p−n junction (relative to the Dirac point) created by a potential barrier of height V0 in region II. (d) Tunneling process through the potential barrier. u, ur, q, qr, θ, θr, and λ′ in region II correspond to v, vr, k, kr, ϕ, ϕr, and λ in region I.

Figure 4. Transmission probability (T) versus the incident angle of group velocity (ϕv′) with respect to the potential barrier normal, showing asymmetric Klein tunneling. E0 = 10 meV is used. (a) Two 1 different potential barrier orientations. At α = αm = arctan γ , the

diagonalization of the Hamiltonian matrix) show that all three external potential models produce sizable massless Dirac cones in the superlattice BZ with quite achievable parameters in experiments (see Supporting Information for more results and discussions, including superlattice periodicity dependence). We have therefore given a general scheme for generating massless Dirac fermions in arbitrarily anisotropic 2DEG from the perspective of inverse design. One important prerequisite for ease in probing and utilizing the emergent Dirac fermion states is to have these states wellseparated from other states. In contrast, tunable and anisotropic Dirac fermions can also be generated in graphene using onedimensional periodic potential4,38 due to the chiral nature of states in the original Dirac cone. However, those new states are typically obscured by other states in the same energy window.4,38 In the approach proposed in this work, the resulting density of states (DOS) has a clear V-shape feature and vanishes at the energy where the two new Dirac cones meet, well-separated from other states (Figure 2d, see Supporting Information). Thus, they can be more easily characterized in many experiments involving properties of Fermi surface, such as scanning tunneling microscope, electrical transport, quantum oscillations, quantum Hall, magneto-

0

normal-incidence direction and the perfect-transmission direction are maximally differed by ϕ′v,m = −50.8°. (b) Symmetric and asymmetric Klein tunneling profiles corresponding to the two geometries in panel a. (c, d) Asymmetric Klein tunneling with various parameters for α = αm .

resistance, and electron optics. In the model potential described by eq 3, the condition W > ℏ2K20/16mx should at least be satisfied to avoid obscuring states stemming from original 2DEG dispersions. Comparison with the DOS of the pristine few-layer black phosphorus in Figure 2d indicates that the system is only slightly modified by the superlattice potential, and the intrinsic band gap of few-layer black phosphorus (∼1 eV) is still well-preserved. Therefore, in such a black phosphorus superlattice, the semiconducting phase and the massless Dirac fermion phase can be reversibly switched by tuning the Fermi level within an achievable carrier density range in the order of 1012−1013 cm−2 (see Supporting Information). Considering the sensitivity of black phosphorus to air, we note that the encapsulation of the material by using hexagonal boron nitride could be helpful, as it virtually does not D

DOI: 10.1021/acs.nanolett.6b04942 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters

Fourier components of a potential deviate from eq 3 slightly, eq 7 is then an approximation. Now we consider transmission through a potential barrier with a height of V0 and width of D0 (D0 ≫ |a1|,|a2|). Let us align the minor and major axes of the elliptical Fermi surface along the kx and ky directions, respectively (see Figure 3). We place a potential barrier in the x′−y′ coordinate system (rotated at some arbitrary angle α from the x−y coordinate system) and make the barrier axis infinitely long along the y′ direction with the width D0 lying along the x′ direction. In the continuum limit, translational symmetry along the y′ direction will conserve the ky′ wavevector component of a propagating wave impinging on the potential barrier, based on which we look for the solutions in the three regionsbefore (I), inside (II), and after (III) the barrier (see Supporting Information). The various notations are defined in the caption of Figure 3. The transmission probability in the ballistic transport limit is obtained by matching the carrier wave function at the boundaries.2,30 The solutions of the three regions in the x′− y′ coordinate system take the form

affect the electronic structure around the band edges of black phosphorus.39 Next, we explore one of the most counterintuitive phenomena in graphene, Klein tunneling. That is, normalincident carriers in graphene experience a perfect transmission through a potential barrier independent of the potential width and height, made possible by the chiral nature of the linearly dispersing Dirac fermions and charge-conjugation symmetry.30 In an anisotropic Dirac cone, for states with wavevector k (measured from the Dirac point), the wavevector k, the pseudospinor s, and the group velocity v are all generally noncollinear at a given k point.4 In fact, they are given (for the potential in eq 3) by k = (kx, ky), s = λℏ(kx , γ0k y)/ kx2 + γ02k y2 , λv

and v = 20 (kx , γ02k y)/ kx2 + γ02k y2 (λ = ±1) and related by (Figure 3a) tan ϕv = γ0 tan ϕs = γ02 tan ϕk = γ02

ky kx

(7)

where ϕk, ϕs, and ϕv are referenced to kx axis (note that s depends on λ but ϕs does not, by definition). If the relevant

⎫ ⎧ ⎛ 1 ⎞ i(k x ′+ k y ′) r ⎛ 1 ⎞ i(kxr′x ′+ ky′y ′) x′ y′ ⎪ ⎪ ψ (x′, y′) = 1 ⎜ + ′ < e x e , 0 ⎜ ⎟ ⎟ r 2 ⎝ λeiϕs ⎠ 2 ⎝ λeiϕs ⎠ ⎪ ⎪ I ⎪ ⎪ ⎪ ⎪ a ⎛ 1 ⎞ i(qx′x ′+ ky′y ′) b ⎛ 1 ⎞ i(qxr′x ′+ ky′y ′) ⎨ ψII(x′, y′) = + , 0 < x′ < D0 ⎬ ⎜ iθ r ⎟e ⎜ iθ ⎟e s s 2 ⎝ λ′e ⎠ 2 ⎝ λ′e ⎠ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ t ⎛ 1 ⎞ i(kx′x ′+ ky′y ′) ψIII(x′, y′) = x′ > D0 ⎪ , ⎜ iϕ ⎟e ⎪ s 2 ⎝ λe ⎠ ⎭ ⎩

The transmission amplitude is t(ϕv′) =

λλ′e

−ikx ′D

(e

iθsr

iθs

tunneling behavior compared with graphene. For a wave packet, it is the group velocity, not wavevector, that describes the direction of center-of-mass motion and energy flow. In graphene, normal incident wave packet (or energy flow) is unimpeded; however, in the anisotropic case, the transmission probability of a normal incident wave packet is not unity and can be tuned and controlled. Furthermore, the normal incident direction and the perfect transmission direction are in general different. The difference between the two directions depends 1 on α and is maximized at αm = arctan γ , taking on a value of

iϕsr

iϕs

− e )(e − e ) A

(9)

where r

r

r

r

A = e−iqx′D(eiθs+ iθs + eiϕs + iϕs − λλ′eiθs + iϕs − λλ′eiθs+ iϕs) r

r

r

r

(8)

r

− e−iqx′D(eiθs+ iθs + eiϕs + iϕs − λλ′eiθs + iϕs − λλ′eiθs+ iϕs ) (10)

The transmission probability through the potential barrier is T(ϕv′) = tt*. In the case of graphene, besides the resonant unit transmission at specific angles owing to a Fabry-Pérot like effect, normal-incident electron wave packet in graphene always shows perfect transmission because of suppression of backscattering30 owing to the chiral character of the Dirac fermions. In the anisotropic case, the pseudospin-momentum locking (given by eqs 6 and 7) gives rise to a complete suppression of the scattering only under the condition k′ → − k′. Considering the fact that ky′ is conserved, this process will happen only when ky′ = 0. Consequently, ϕrs = π + ϕs, θrs = π + θs, ϕs = θs, and qrx′ = −qx′, resulting in T = 1, i.e., perfect transmission independent of the potential barrier height and width. With an anisotropic Dirac cone, perfect transmission therefore occurs when the incident wavevector k is along the normal to the potential barrier. However, the fact that the group velocity and the wavevector are generally noncollinear (eq 7 and Figure 3a) leads to a remarkably distinct Klein

0

ϕ′v,m = −50.8° with the effective masses considered here for black phosphorus (see Supporting Information). The directional dependence of the transmission probability T(ϕv′) is plotted in Figure 4, with varying incident angle ϕv′ of the group velocity relative to the potential barrier normal, and in principle is measurable from a directionally dependent nonlocal resistivity experiment.40 Because of the asymmetric Klein tunneling behavior, the transmission probability of normal incident wave packets can be controlled with different barrier alignment, height, and width, whereas as discussed above the perfect transmission direction depends on α only. In the cases where the Fermi surface is smaller inside the potential barrier region than the outside regions, then at certain angles, no propagating wave solutions exist. But evanescent wave solutions are allowed by letting qx′ = iκ′ and qrx′ = iκr′ (κ′, κ r′ ∈ ). Therefore, a beam of normalincident wave packets can be completely turned off by tuning E

DOI: 10.1021/acs.nanolett.6b04942 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters

(9) Sushkov, O. P.; Castro Neto, A. H. Phys. Rev. Lett. 2013, 110, 186601. (10) Wang, S.; Tan, L. Z.; Wang, W.; Louie, S. G.; Lin, N. Phys. Rev. Lett. 2014, 113, 196803. (11) Wunsch, B.; Guinea, F.; Sols, F. New J. Phys. 2008, 10, 103027. (12) Tarruell, L.; Greif, D.; Uehlinger, T.; Jotzu, G.; Esslinger, T. Nature 2012, 483, 302−305. (13) Rechtsman, M. C.; Zeuner, J. M.; Tünnermann, A.; Nolte, S.; Segev, M.; Szameit, A. Nat. Photonics 2012, 7, 153−158. (14) Rechtsman, M. C.; Zeuner, J. M.; Plotnik, Y.; Lumer, Y.; Podolsky, D.; Dreisow, F.; Nolte, S.; Segev, M.; Szameit, A. Nature 2013, 496, 196−200. (15) Pletikosić, I.; Kralj, M.; Pervan, P.; Brako, R.; Coraux, J.; N’Diaye, A. T.; Busse, C.; Michely, T. Phys. Rev. Lett. 2009, 102, 056808. (16) Yankowitz, M.; Xue, J.; Cormode, D.; Sanchez-Yamagishi, J. D.; Watanabe, K.; Taniguchi, T.; Jarillo-Herrero, P.; Jacquod, P.; LeRoy, B. J. Nat. Phys. 2012, 8, 382−386. (17) Ponomarenko, L. A.; Gorbachev, R. V.; Yu, G. L.; Elias, D. C.; Jalil, R.; Patel, A. A.; Mishchenko, A.; Mayorov, A. S.; Woods, C. R.; Wallbank, J. R.; Mucha-Kruczynski, M.; Piot, B. A.; Potemski, M.; Grigorieva, I. V.; Novoselov, K. S.; Guinea, F.; Fal’ko, V. I.; Geim, A. K. Nature 2013, 497, 594−597. (18) Dean, C. R.; Wang, L.; Maher, P.; Forsythe, C.; Ghahari, F.; Gao, Y.; Katoch, J.; Ishigami, M.; Moon, P.; Koshino, M.; Taniguchi, T.; Watanabe, K.; Shepard, K. L.; Hone, J.; Kim, P. Nature 2013, 497, 598−602. (19) Hunt, B.; Sanchez-Yamagishi, J. D.; Young, A. F.; Yankowitz, M.; LeRoy, B. J.; Watanabe, K.; Taniguchi, T.; Moon, P.; Koshino, M.; Jarillo-Herrero, P.; Ashoori, R. C. Science 2013, 340, 1427−1430. (20) Yang, W.; Chen, G.; Shi, Z.; Liu, C.-C.; Zhang, L.; Xie, G.; Cheng, M.; Wang, D.; Yang, R.; Shi, D.; Watanabe, K.; Taniguchi, T.; Yao, Y.; Zhang, Y.; Zhang, G. Nat. Mater. 2013, 12, 792−797. (21) Yu, G. L.; Gorbachev, R. V.; Tu, J. S.; Kretinin, A. V.; Cao, Y.; Jalil, R.; Withers, F.; Ponomarenko, L. A.; Piot, B. A.; Potemski, M.; Elias, D. C.; Chen, X.; Watanabe, K.; Taniguchi, T.; Grigorieva, I. V.; Novoselov, K. S.; Fal’ko, V. I.; Geim, A. K.; Mishchenko, A. Nat. Phys. 2014, 10, 525−529. (22) Wang, E.; Lu, X.; Ding, S.; Yao, W.; Yan, M.; Wan, G.; Deng, K.; Wang, S.; Chen, G.; Ma, L.; Jung, J.; Fedorov, A. V.; Zhang, Y.; Zhang, G.; Zhou, S. Nat. Phys. 2016, 12, 1111. (23) Li, L.; Yu, Y.; Ye, G. J.; Ge, Q.; Ou, X.; Wu, H.; Feng, D.; Chen, X. H.; Zhang, Y. Nat. Nanotechnol. 2014, 9, 372−377. (24) Liu, H.; Neal, A. T.; Zhu, Z.; Luo, Z.; Xu, X.; Tománek, D.; Ye, P. D. ACS Nano 2014, 8, 4033−4041. (25) Li, L.; Ye, G. J.; Tran, V.; Fei, R.; Chen, G.; Wang, H.; Wang, J.; Watanabe, K.; Taniguchi, T.; Yang, L.; Chen, X. H.; Zhang, Y. Nat. Nanotechnol. 2015, 10, 608−613. (26) Li, L.; Yang, F.; Ye, G. J.; Zhang, Z.; Zhu, Z.; Lou, W.; Zhou, X.; Li, L.; Watanabe, K.; Taniguchi, T.; Chang, K.; Wang, Y.; Chen, X. H.; Zhang, Y. Nat. Nanotechnol. 2016, 11, 593−597. (27) Qiao, J.; Kong, X.; Hu, Z.-X.; Yang, F.; Ji, W. Nat. Commun. 2014, 5, 4475. (28) Fei, R.; Yang, L. Nano Lett. 2014, 14, 2884−2889. (29) Rodin, A. S.; Carvalho, A.; Castro Neto, A. H. Phys. Rev. Lett. 2014, 112, 176801. (30) Katsnelson, M. I.; Novoselov, K. S.; Geim, A. K. Nat. Phys. 2006, 2, 620−625. (31) Young, A. F.; Kim, P. Nat. Phys. 2009, 5, 222−226. (32) Rickhaus, P.; Maurand, R.; Liu, M.-H.; Weiss, M.; Richter, K.; Schönenberger, C. Nat. Commun. 2013, 4, 2342. (33) Wilmart, Q.; Berrada, S.; Torrin, D.; Nguyen, V. H.; Fève, G.; Berroir, J.-M.; Dollfus, P.; Plaçais, B. 2D Mater. 2014, 1, 011006. (34) Pereira, J. M., Jr.; Katsnelson, M. I. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 92, 075437. (35) Goerbig, M. O.; Fuchs, J.-N.; Montambaux, G.; Piéchon, F. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 045415. (36) Liu, Z.; Wang, J.; Li, J. Phys. Chem. Chem. Phys. 2013, 15, 18855.

V0, which is also a missing feature in pristine graphene or any isotropic Dirac fermion system. The phenomena found in this work are general and applicable to any anisotropic host systems. Moreover, the emergence of highly tunable and easily accessible anisotropic massless Dirac fermions in few-layer black phosphorus superlattices should provide a range of interesting experimental investigations and a new direction for possible device applications.



ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.nanolett.6b04942. discussions in three realistic potential profiles, dependence of superlattice periodicity, isolation of the new Dirac fermions within an energy window, and details of solving for asymmetric Klein tunneling (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Zhenglu Li: 0000-0002-3851-9241 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the Theory of Materials Program at the Lawrence Berkeley National Laboratory through the Office of Basic Energy Sciences, U.S. Department of Energy under Contract No. DE-AC02-05CH11231 which provided theoretical analyses and calculations of generated anisotropic Dirac fermions; and by the National Science Foundation under Grant No. DMR-1508412 which provided numerical simulations of Klein tunneling. Computational resources have been provided by the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy, and the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI1053575. The authors thank L. Li, A. Rubio, F. H. da Jornada, G. Antonius, and Y.-W. Son for fruitful discussions.



REFERENCES

(1) Geim, A. K.; Novoselov, K. S. Nat. Mater. 2007, 6, 183−191. (2) Castro Neto, A. H.; Guinea, F.; Peres, N. M. R.; Novoselov, K. S.; Geim, A. K. Rev. Mod. Phys. 2009, 81, 109. (3) Park, C.-H.; Louie, S. G. Nano Lett. 2009, 9, 1793−1797. (4) Park, C.-H.; Yang, L.; Son, Y.-W.; Cohen, M. L.; Louie, S. G. Phys. Rev. Lett. 2008, 101, 126804. (5) Gibertini, M.; Singha, A.; Pellegrini, V.; Polini, M.; Vignale, G.; Pinczuk, A.; Pfeiffer, L. N.; West, K. W. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 241406. (6) Singha, A.; Gibertini, M.; Karmakar, B.; Yuan, S.; Polini, M.; Vignale, G.; Katsnelson, M. I.; Pinczuk, A.; Pfeiffer, L. N.; West, K. W.; Pellegrini, V. Science 2011, 332, 1176−1179. (7) Gomes, K. K.; Mar, W.; Ko, W.; Guinea, F.; Manoharan, H. C. Nature 2012, 483, 306−310. (8) Polini, M.; Guinea, F.; Lewenstein, M.; Manoharan, H. C.; Pellegrini, V. Nat. Nanotechnol. 2013, 8, 625−633. F

DOI: 10.1021/acs.nanolett.6b04942 Nano Lett. XXXX, XXX, XXX−XXX

Letter

Nano Letters (37) Bradlyn, B.; Cano, J.; Wang, Z.; Vergniory, M. G.; Felser, C.; Cava, R. J.; Bernevig, B. A. Science 2016, 353, aaf5037. (38) Park, C.-H.; Yang, L.; Son, Y.-W.; Cohen, M. L.; Louie, S. G. Nat. Phys. 2008, 4, 213−217. (39) Cai, Y.; Zhang, G.; Zhang, Y.-W. J. Phys. Chem. C 2015, 119, 13929−13936. (40) Lee, G.-H.; Park, G.-H.; Lee, H.-J. Nat. Phys. 2015, 11, 925− 929.

G

DOI: 10.1021/acs.nanolett.6b04942 Nano Lett. XXXX, XXX, XXX−XXX