Anomalous caustics and Veselago focusing in 8-Pmmn borophene p-n junctions with arbitrary junction directions

Negative refraction usually demands complex structure engineering while it is very natural for massless Dirac fermions (MDFs) across the \textit{p-n} junction, this leads to Dirac electron optics. The emergent Dirac materials may exhibit hitherto unidentified phenomenon due to their nontrivial band structures in contrast to the isotropic MDFs in graphene. Here, as a specific example, we explore the negative refraction induced caustics and Veselago focusing of tilted MDFs across 8-\textit{Pmmn} borophene \textit{p-n} junctions. To this aim, we develop a technique to effectively construct the electronic Green's function in \textit{p-n} junctions with arbitrary junction directions. Based on analytical discussions and numerical calculations, we demonstrate the strong dependence of interference pattern on the junction direction. As the junction direction perpendicular to the tilt direction, Veselago focusing or normal caustics (similar to that in graphene) appears resting on the doping configuration of the \textit{p-n} junctions, otherwise anomalous caustics (different from that in graphene) occurs which is manipulated by the junction direction and the doping configuration. Finally, the developed Green's function technique is generally promising to uncover the unique transport of emergent MDFs, and the discovered anomalous caustics makes tilted MDFs potential applications in Dirac electron optics.

Negative refraction usually demands complex structure engineering while it is very natural for massless Dirac fermions (MDFs) across the p-n junction, this leads to Dirac electron optics. The emergent Dirac materials may exhibit hitherto unidentified phenomenon due to their nontrivial band structures in contrast to the isotropic MDFs in graphene. Here, as a specific example, we explore the negative refraction induced caustics and Veselago focusing of tilted MDFs across 8-Pmmn borophene p-n junctions. To this aim, we develop a technique to effectively construct the electronic Green's function in p-n junctions with arbitrary junction directions. Based on analytical discussions and numerical calculations, we demonstrate the strong dependence of interference pattern on the junction direction. As the junction direction perpendicular to the tilt direction, Veselago focusing or normal caustics (similar to that in graphene) appears resting on the doping configuration of the p-n junctions, otherwise anomalous caustics (different from that in graphene) occurs which is manipulated by the junction direction and the doping configuration. Finally, the developed Green's function technique is generally promising to uncover the unique transport of emergent MDFs, and the discovered anomalous caustics makes tilted MDFs potential applications in Dirac electron optics.

I. INTRODUCTION
Negative refraction is one unusual class of wave propagation 1 , which leads to novel interference and hold great potential in new applications 2 . But its realization usually demands complex structure engineering with stringent parameter conditions 3 . Graphene as the first two-dimensional material, hosts the relativistic massless Dirac fermions (MDFs) and possesses a lot of exotic physics and possible applications 4 . Especially, the two-dimensional nature of graphene is beneficial to the fabrication of the planar p-n junction (PNJ), which is the basic component of many electronic devices 5,6 . The propagation of MDFs across graphene PNJ has a close analogy to optical negative refraction at the surface of metamaterials but exhibits negative refraction in a more simple and tunable manner 7 . Recently, negative refraction in the graphene PNJ has been verified experimentally 8,9 . As a result, there is wide interest to study the negative refraction induced interference of MDFs in the PNJ structure [10][11][12][13][14][15][16][17][18][19][20][21] .
The great success of graphene attracts people to search for new two-dimensional materials in which quasiparticles can be described as MDFs 22,23 , i.e., Dirac materials. Dirac materials emerge quickly and usually host very novel MDFs, in which new physics and application potential are expected [22][23][24] . The Dirac fermions can be classified finely into four categories 25 , i.e., type-I, type-I tilted, type-III (critical tilt), and type-II ones. These four categories of Dirac fermions are expected to exhibit distinct physical properties due to different geometries of their Fermi surface. To our knowledge, type-I tilted Dirac fermions has been predicted to appear in very rare systems including quinoid-type graphene and -(BEDT-TTF) 2 I 3 26 , hydrogenated graphene 27 , and 8-Pmmn borophene [28][29][30][31] . In such Dirac materials, 8-Pmmn borophene as elemental monolayer material exhibits high mobility and anisotropic transport 32 , which has potential applications in electronics and electron optics. The well-established continuum model of 8-Pmmn borophene make it be very suit-able to the model study. This anisotropy and tilt of type-I Dirac fermions brings about unique features to various physical properties of 8-Pmmn borophene, including plasmon 33,34 , the optical conductivity 35 , Weiss oscillations 36 , oblique Klein tunneling 37 , metal-insulator transition induced by strong electromagnetic radiation 38 , and RKKY interaction 39,40 . Thus, it is expected that the anisotropy and tilt will strongly affect the negative refraction of MDFs cross 8-Pmmn borophene PNJ.
In this study, we investigate the negative refraction induced interference in 8-Pmmn borophene PNJs. Following the seminal study in graphene 7 , we focus on two typical interference pattern, i.e, caustics and Veselago focusing. Due to the anisotropy and tilt of MDFs in 8-Pmmn borophene, the interference pattern has inevitable dependence on the junction direction. Hence, we develop a construction technique of the electronic Green's function (GF), which is applicable to the PNJ with an arbitrary junction direction. We find: (I) Veselago focusing or normal caustics appears resting on the doping configuration of the PNJ when its junction direction is perpendicular to the tilt direction of MDFs. (II) To the other junction, anomalous caustics unique to anisotropic and tilted MDFs occurs, and its pattern strongly depends on the junction direction and the doping configuration of the PNJ. We expect that the developed GF technique is used to study the other transport properties of continuing emergent MDFs, and the discovered anomalous caustics can be observed in near future pointing out the potential application of novel MDFs (in 8-Pmmn borophene or other Dirac materials) in Dirac electron optics.
The rest of this paper is organized as follows. In Sec. II, we introduce the model structures and Hamiltonian of the 8-Pmmn borophene PNJ, and present the detailed construction technique of the electronic GF. In Sec. III, according to the junction direction of PNJs, we demonstrate Veselago focusing and anomalous caustics by combing numerical calculations and analytical derivation. Finally, we summarize this study in Sec. IV. We introduce the Cartesian coordinate system x ′ −y ′ accompanying the intrinsic HamiltonianĤ 0 of 8-Pmmn borophene. In 8-Pmmn borophene,Ĥ 0 describes the anisotropic and tilted MDFs, and has the valley dependence. However, we focus the new physics induced by the anisotropy and tilt, so the valley dependence ofĤ 0 is neglected in our model study. Around the Dirac point,Ĥ 0 has the form 30,33,36 wherep x ′ ,y ′ are the momentum operators, σ 1,2 and σ 0 are Pauli matrices and 2 × 2 identity matrix, respectively. The anisotropic velocities are v 1 = 0.86v F , v 2 = 0.69v F , and v t = 0.32v F with v F = 10 6 m/s. Throughout this paper, we set = v F ≡ 1 to favor our dimensionless derivation and calculations, and they can be used to define the length unit l 0 and the energy unit ε 0 through v F = l 0 ε 0 , e.g., ε 0 = 0.66 eV when l 0 = 1 nm. The eigenenergies and eigenstates ofĤ 0 are, respectively, and Here, λ = + (λ = −) denotes the conduction (valence) band, k ′ = (k x ′ , k y ′ ) is the wave vector and r ′ = (x ′ , y ′ ) is the position vector. In particular, θ S λ (k ′ ) is the azimuthal angle of the in-plane pseudospin vector S λ (k ′ ) = Ψ λ,k ′ |(σ 1 , σ 2 )|Ψ λ,k ′ relative to x ′ -axis, which has the form 37 : The anisotropy and tilt of MDFs lead to the noncollinear feature of wave vector and group velocity for a general state, which has profound modifications to physical properties of 8-Pmmn borophene comparing to those of graphene, e.g, oblique Klein tunneling 37 . And this nonlinear feature can be clearly shown by the definition of group velocity The two-dimensional nature of 8-Pmmn borophene is suitable to the fabrication of the planar PNJ. One typical 8-Pmmn borophene PNJ is shown schematically by Fig. 1(a). In Fig.  1(a), we use the normal (tangential) direction of junction interface to define x (y) axis of the Cartesian coordinate system x − y which has a rotation angle φ relative to the coordinate system x ′ − y ′ for the intrinsic HamiltonianĤ 0 , so φ can be regarded as the junction direction of φ-junction. In the coordinate systems x − y and x ′ − y ′ , an arbitrary vector A can be expressed as A = (A x , A y ) and A ′ = (A x ′ , A y ′ ) and they are related to each other A T = U(A ′ ) T through the unitary matrix The Hamiltonian of 8-Pmmn borophene PNJ in Fig. 1

(a) iŝ
where is the gate-induced scalar potential in the n (p) region as shown by Fig. 1(b) and no loss of generality we assume V 0 > 0, and Θ(x) is the step function: where a positive (negative) doping level corresponding to electron (hole) doping, so ε n > 0 and ε p < 0 for the PNJ.

B. Green's function
The propagation of anisotropic and tilted MDFs in the PNJ structure can be described by the corresponding propagator or of the HamiltonianĤ. Usually, the position and energy dependence of GF are not shown explicitly for brevity in our study, i.e., G ≡ G(E F , r 2 , r 1 ). For the quasi-one dimensional systems such as the considered PNJ in this study, the GF can be constructed through the on-shell spectral expansion 41 , i.e., it is just the spectral expansion of the states on the Fermi surface instead of all the eigenstates including on-shell and offshell ones of the system in the conventional spectral expansion method for the GF [42][43][44] . More importantly, according to the on-shell spectral expansion 41 , GF has a physical transparent formalism which favors its analytical construction through the intrinsic states and its subsequent scattering states by the scattering mechanism in the quasi-one dimensional system, e.g., the junction interface of the PNJ. , and extra GF constructed through the right-going intrinsic states and its left-going reflection states by the PNJ interface. (c)/(d) G is the sum of intrinsic GF G 0 constructed through rightgoing/left-going intrinsic states of the p region and extra GF constructed through the left-going intrinsic states and its right-going reflection states by the PNJ interface. (e)/(f) G is constructed through the right-going/left-going intrinsic states of the n/p region and its right-going/left-going transmission states of the p/n region. Fig. 2 shows the piecewise construction of the GF G through the intrinsic and scattering states in the PNJ structure. The intrinsic state are the eigenstates of n and p regions in the PNJ while the scattering states are the renormalized eigenstates by the scattering coefficients. In Fig. 2, it is the interface of the PNJ leading to the scattering of incident intrinsic states. Considering the different incident states (maybe right-going or left-going) and the relative position of r 1 and r 2 , there are total 6 pieces for the GF as shown by Fig. 2(a)-(f). To illustrate the operation of the construction technique, no loss of generality, we assume the right-going electron states incident from the left n region of the PNJ, then the resulting states can be expressed as: Here, k α,± = (k α,±,x , k y ) are the right-going (+) and left-going (−) wave vectors in the α region, and the conservation of the tangential momentum k y has been used. For brevity, we omit the subscript λ of the eigenstates since its value can be fully determined by the subscript α of the momentum through λ = sgn(ǫ α ), and we also use this simplification to the pseudospin vector in the following. Since k y is conserved, the next step is to derive k α,±,x . By using the unitary matrix (cf. Eq. 6) for the vector transfromation between two coordinate systems, we have where Here, γ 2 = γ 2 1 − γ 2 2 with γ 1 = v 2 /v 1 and γ 2 = v t /v 1 , and ǫ α = ε α /v 1 with ε n = E F + V 0 and ε p = E F − V 0 . Eq. 11 is a quadratic equation with one unknown which gives two roots for k α,±,x : So k α,±,x is derived, and then k ′ α,± is given by Eq. 10. And r(k y ) and t(k y ) are, respectively, reflection and transmission coefficients, which had been derived for the PNJ with a general junction direction 37 : t(k y ) = exp(iθ S(k n,− ) ) − exp(iθ S(k n,+ ) ) exp(iθ S(k n,− ) ) − exp(iθ S(k p,+ ) ) .
Here, v n,x ≡ v n,x (k n,+ ) is the x component of the group velocity v n (k n,+ ) in the Cartesian coordinate system x − y. When x 2 < 0, Eq. 15 is the mathematical description of Fig. 2(a)/(b), which is the sum of intrinsic GF G 0 constructed through rightgoing/left-going intrinsic states of the n region (i.e., the free GF G 0 ) and extra GF constructed through the right-going intrinsic states and its left-going reflection states by the PNJ interface. The free GF G 0 ≡ G 0 (E F , r 2 , r 1 ) is constructed completely by the intrinsic states and has the form: which is consistent with the result derived through the straightforward Fourier transformation of the GF in momentum space 40 . When x 2 > 0, Eq. 15 is the mathematical description of Fig. 2(e), which is constructed through the rightgoing intrinsic states of the n region and its right-going transmission states of the p region. Most importantly, when is the propagation phase accumulated through the negative refraction of incident right-going intrinsic states in the n region to become right-going transmission states in the p region 16 . Similarly, to assume the left-going electron states incident from the right p region of the PNJ, one can obtain the mathematical description of Fig. 2(c), (d) and (f). Noting here, the developed construction technique of GF is universal to the PNJ with a general junction direction.

III. RESULTS AND DISCUSSIONS
In this section, we present the numerical results to display the negative refraction induced interference pattern of anisotropic and tilted MDFs across the 8-Pmmn borophene PNJ and discuss the underlying physics. No loss of generality, we focus on the GF matrix element |G 11 (E F , r 2 , r 1 )|, since the interference pattern originates from the phase accumulation across the PNJ as shown in the following and has no qualitative difference among different GF matrix elements.
Previous to the detailed calculations, we firstly present the intuitive physics for the interference pattern in 8-Pmmn borophene PNJ. In contrast to the isotropic MDFs across the graphene PNJ 7 , the anisotropy and tilt of MDFs should lead to the dependence of interference pattern on the junction direction of the 8-Pmmn borophene PNJ. The junction direction dependence implies the distinct momentum matching of states across borophene PNJ, then determines the appearance and the effectiveness of interference pattern. Fig. 3 shows the dependence of momentum matching on the junction direction φ of the borophene PNJ. Six typical junction directions are plotted, i.e., φ = 0, π, ±π/6, and ±π/2 as shown in Fig.  3(a)-(f). The available states on the electron (green) and hole (blue) Fermi surfaces for the momentum matching subject to the conserved momentum k y (cf. red dashed lines in Fig. 3). In , and (f) φ = −π/2. The available states on the electron (green) and hole (blue) Fermi surfaces for the momentum matching subject to the conserved momentum k y as shown by the red dashed lines. The colored lines with arrows represent the group velocities v i,r,t for the incident, reflection and transmission states, which are nonlinear with the conserved momenta. Here, we use ε n = ε p = 0.1 as the doping configuration for p-n junctions.
n region of the PNJ, we also plot the group velocities v i,r,t for the right-going incident, left-going reflection and right-going transmission states. The noncollinear features between group velocities and conserved momenta are obvious. Comparing to the partial matching in Fig. 3(a)-(d), all states on the electron and hole Fermi surfaces may contribute to interference pattern in Fig. 3(e) and (f), then lead to the effective interference. More importantly, in Fig. 3(e) and (f), the electron and hole Fermi surfaces have the mirror symmetry about the junction interface, so Veselago focusing is expected in these two special junctions 45 . As a representative interference pattern, Veselago focusing are studied intensively and are understood deeply 7,10-17 . Therefore, the interference pattern in the special junctions with φ = ±π/2 being perpendicular to the tilt direction (namely, y ′ axis), can be as the starting point of our discussions. Then, we turn to the junctions with φ = 0 and φ = π (being parallel to the tilt direction) and the general junctions.
A. φ = ±π/2: Veselago focusing and normal caustics We have picked out two special junctions perpendicular to the tilt direction (namely, y ′ axis) of MDFs in 8-Pmmn borophene. In these two special junctions, we plot the ma- trix element of GF |G 11 (E F , r 2 , r 1 )| in Fig. 4. For two different junction directions (i.e., φ = ±π/2), three different doping configurations are considered (i.e., E F = 0, ±0.03). For both special junctions, Veselago focusing as expected 45 (normal caustics very similar to that in graphene PNJ 7 ) appears when E F = 0 (E F = ±0.03). For these two special junctions, we can analytically reveal the underlying physics of Veselago focusing and caustics. The tilt of MDFs in 8-Pmmn borophene leads to the breaking of mirror symmetry about the x ′ axis, so two cases of junction direction along jx ′ axis with j = ± are inequivalent (cf. momentum matching in Fig. 3(e) and (f)). For the junction direction along jx ′ direction, i.e., φ = − jπ/2, k j α,± = (k j α,±,x , k y ) has the form (cf. Eq. 13 for momentum components k j α,±,x ): Here, λ = sgn(ǫ α ) conforms to Eq. 2 of eigenenergies. To show the formation of interference pattern of anisotropic and tilted MDFs across the PNJ, we examine the classical trajectory determined by the phase accumulation ϕ np (cf. Eq. 17). The classical trajectory with specific k y,c is determined by ∂ k y ϕ j np (k y )| k y,c = 0 as 16 ∂ k y k j p,+,x x 2 − ∂ k y k j n,+,x x 1 + (y 2 − y 1 ) = 0.
To differentiate Eq. 18, one can obtain As a result, the classical trajectory with k y,c satisfies the equation: Noting here, Eq. 21 has no dependence on j, so it has the same form for two special junctions with φ = ±π/2. And Eq. 21 reproduces the result for isotropic MDFs when v t = 0 and v 1 = v 2 , i.e., γ = γ 1 = 1 37 .

Symmetric doping: Veselago focusing
To consider the symmetric doping, i.e., E F = 0 or ǫ n = −ǫ p = V 0 /v 1 , the classical trajectory is Obviously, when x 1 + x 2 = 0 and y 2 = y 1 , Eq. 22 has no dependence on the momentum of classical trajectory, this implies that all classical trajectories from r 1 = (x 1 , y 1 ) converge to its mirror point r 2 = (−x 1 , y 1 ) about the junction interface, i.e., the appearance of Veselago focusing. Meanwhile, k j n,+,x = −k j p,+,x and ∂ m k y k j n,+,x = −∂ m k y k j p,+,x with m being a positive integer, so the zero-order term ϕ j np (k y,c ) = −k j n,+,x (x 2 +x 1 )+ k y,c (y 2 − y 1 ) = 0 and the arbitrary order derivative of propagation phase ϕ j np on the classical trajectory also vanishes identical to the case for isotropic MDFs in graphene 37 . However, noting that the different color scales for the Veselago focusing and caustics on the upper and bottom rows of Fig. 4, implying the different focusing magnitude. This reflects the inequivalence of two special junctions as mentioned previously, which has different density of states for the focusing states.

Asymmetric doping: normal caustics
To consider the asymmetric doping configuration (i.e., E F 0), Veselago focusing disappears, but the classical trajectories can form another kind of interference pattern (i.e., caustics) if the quadratic term of ϕ np (k y ) vanishes 37 . Using Eq. 18, the vanishing quadratic term ∂ 2 k 2 y ϕ np (k y )| k y,c = 0 leads to To solve the above equation, we derive where c = n 2 (x 2 /x cump ) 2/3 with n = −ǫ p /ǫ n and x cump = −nx 1 . As a result, the classical trajectory is Eq. 25 reproduces the result for isotropic MDFs 37 when γ = γ 1 = 1, and it just indicates the normal caustics. Due to the anisotropy (v 1 v 2 ) and tilt (v t 0) of MDFs in 8-Pmmn borophene, γ 1 leads to the modification of caustics. From a different perspective, this implies the tunability of caustics by changing the anisotropy and tilt of MDFs. In addition, the analytical Eq. 25 is used to plot the red lines in Fig. 4, and they are very consistent with the numerical results.
B. φ = 0 and π: anomalous caustics Fig. 5 shows |G 11 (E F , r 2 , r 1 )| in PNJs with junction directions (i.e., φ = 0, π) parallel to tilt direction (namely y ′ axis). In Fig. 5, an outstanding feature is that the upper and bottom panels have the mirror symmetry about the x axis. This comes from one hidden symmetry: To perform inversion operations y → −y and k y → −k y for φ-junction, its momentum matching becomes the case of (π − φ)-junction (cf. Fig. 3). Using this hidden symmetry, one can also explain the mirror symmetry of interference pattern about the x axis in Fig. 4 since π/2junction or −π/2-junction is symmetrically related to itself.
Comparing to Fig. 4 for special junctions, Veselago focusing does not exist and caustics becomes anomalous in Fig. 5, i.e., anomalous caustics. The MDFs in 8-Pmmn borophene have two features: anisotropy (v 1 v 2 ) and tilt (v t 0). To assume v t = 0, the electron and hole Fermi surfaces have the mirror symmetry about the PNJ interface when E F = 0, so Veselago focusing exists 18,45 . The absence of Veselago focusing in Fig. 5 is due to the breaking of mirror symmetry by the finite tilt (cf. Fig. 3(a) and (b)). Nevertheless, the anisotropy and tilt both contribute to the caustics, even in the special junctions (cf. Eq. 25). The φ-dependent momentum matching further introduces the dependence of the anomalous caustics on the junction direction. For the junction direction along jy ′ direction, we have k j α,± = (k j α,±,x , k y ) and k j α,±,x are given by Eq. 13: which implies Using Eq. 17, the vanishing quadratic term ∂ 2 (28) The solution of the above equation is (29) Substituting Eqs. 27 and 29 into Eq. 19, the anomalous caustics can be determined and is shown by the red lines in Fig.  5. The analytical and numerical results agree very well with each other.
In addition, Fig. 5 also shows the dependence of anomalous caustics on the doping configuration. With decreasing E F from left to right on the upper (bottom) row, the anomalous caustics goes upward (downward). This implies the tunable anomalous caustics by doping configuration. It attracts us to further increase the magnitude of E F to examine the interference pattern. Fig. 6 shows the interference pattern of |G 11 (E F , r 2 , r 1 )| same as Fig. 5 except that two higher magnitude of Fermi levels (E F = ±0.06) are considered for (φ = 0)junction. On the caustics, Fig. 5 has a single branch, whereas there are two branches in Fig. 6. In the sense of the same branch number for caustics in Fig. 4 and Fig. 6, the anomalous feature is suppressed by the high value of E F . And there should be one critical E F at which the anomalous caustics changes its branch number. In all, the anomalous caustics has a strong dependence on doping configuration and exhibits rich patterns.

C. General φ: anomalous caustics
In general cases, the junction directions φ of PNJs are neither perpendicular nor parallel to the tilt direction of MDFs in 8-Pmmn borophene. It is rather difficult to examine the transport properties of these general φ-junctions, because their low symmetry requires large unit cell to perform numerical simulation, e.g., through the recursive GF based on the lattice Hamiltonian 46 . By using the developed construction technique of GF in Sec. IIB, the transport problem in general φ-junctions can be solved efficiently. Fig. 7 shows FIG. 7. Anomalous caustics in the PNJ with general junction directions shown by |G 11 (E F , r 2 , r 1 )|. Two different junction directions (i.e., φ = ±π/6) and three different Fermi levels (i.e., E F = 0, ±0.03) are considered.
|G 11 (E F , r 2 , r 1 )| by considering two different junction directions (i.e., φ = ±π/6) and three different Fermi levels (i.e., E F = 0, ±0.03). In light of the mirror symmetry about the x axis for φ-junction and (π − φ)-junction, one only needs to explicitly discuss the junction directions φ ∈ [−π/2, π/2], so φ = ±π/6 are typical enough. Fig. 7 is very similar to Fig.  5, featured by the anomalous caustics. And with decreasing E F from left to right on the upper or bottom row, the anomalous caustics goes upward. Therefore, the anomalous caustics depends strongly on the junction directions and doping configuration, favoring its tunability.

IV. CONCLUSIONS
In summary, we investigate the negative refraction induced interference of anisotropic and tilted MDFs across 8-Pmmn borophene PNJs. Because of the anisotropy and tilt, the calculation of GF of the PNJ with an arbitrary junction direction is needed, for which we develop an effective construction technique. The developed construction technique of GF can be generalized to study the transport of other novel MDFs, e.g., type-II Dirac fermions 25,47 , and semi-Dirac fermions 48 . Comparing to the seminal work 7 , here we focus negative refraction induced Veselago focusing and caustics. If the junction direction is perpendicular to the tilt direction of MDFs, Veselago focusing or normal caustics appears resting on the doping configuration of the PNJ. More importantly, to the other junction, we discover one new phenomenon, i.e., anomalous caustics, which can be manipulated by junction direction and doping configuration. The interference pattern is unique to anisotropic and tilted MDFs, which is demonstrated numerically and analytically. This model study is generally applicable to the ideal p-n junctions with perfect border, while the imperfect border (e.g., the disorder-induced losses and smooth junctions) will lead to the diffusive scattering and reduce the transmission of electrons across the electrons. To implement the Dirac fermions into electron optics, people make a great effort to improve the experimental fabrication of ideal p-n junction. For example, very ideal atomically sharp p-n junction have been created [49][50][51] . In addition, with the rapid experimental advance of borophene [52][53][54][55] and the demonstration of negation refraction 8,9,21 , we expect the anomalous caustics to be observable in the near future, then this study make novel MDFs (e.g., in 8-Pmmn borophene) promising to engineer Dirac electron optics devices.