Orientation-dependent crossover from retro to specular Andreev reflections in semi-Dirac materials

In the framework of Bogoliubov–de Gennes equation, we theoretically study the transport properties in normal-superconducting junctions based on semi-Dirac materials (SDMs). Owing to the intrinsic anisotropy of SDMs, the configuration of Andreev reflection (AR) and the differential conductance are strongly orientation-dependent. For the transport along the linear dispersion direction, the differential conductance exhibits a clear crossover from retro AR to specular AR with increasing the bias-voltage, and the differential conductance oscillates with the interfacial barrier strength without a decaying profile. Conversely, for the transport along the quadratic dispersion direction, the boundary between the retro AR and specular AR becomes ambiguous when the orientation angle increases, and the differential conductance decays with increasing the momentum mismatch or the interfacial barrier strength. We illustrate the pseudo-spin textures to reveal the underling physics behind the anisotropic coherent transport properties. These results enrich the understanding of the superconducting coherent transport in SDMs.


Introduction
In a normal-superconducting (NS) junction with ideal contacts, the transport properties are dominated by the Andreev reflection (AR) in the subgap energy regime of E Δ 0 , with E the incident energy and Δ 0 the superconducting gap [1,2]. In accordance with the origination of the AR-involved electron and hole, the AR can be classified into two typical categories, i.e., the intra-band and inter-band phase-coherent processes. In most conventional-metal-based NS junctions, the chemical potential in the N region satisfies μ N Δ 0 , thus in the subgap regime the AR is a intra-band phase-coherent scattering process, during which an incident electron from the N region is approximately reto-reflected as a hole [1][2][3][4]. While in the NS junctions based on Dirac materials, μ N can be continuously tuned to the regime of μ N < E [5][6][7]. Consequently, a conduction-band electron incident from the N region is specularly reflected back as a valence-band hole, leading to a inter-band phase-coherent scattering process known as specular-AR [7][8][9][10][11][12][13][14][15][16][17]. Remarkably, when the N region is weakly doped, a crossover from reto-AR to specular-AR occurs with increasing E, manifesting itself as a dip at E = μ N in the E-dependent conductance spectrum. This signature has been experimentally observed in bilayer-graphene-based NS junctions [7]. In addition, in Dirac-material-based NS junctions, due to the novel momentum-spin/pseudo-spin textures of Dirac fermions, the subgap differential conductance oscillates with the interfacial barrier strength without a decaying profile [12], as confirmed by recent experimental measurements [18].

Model and method
In the leading order of k, the single-particle effective Hamiltonian of SDM reads h(k) = vk xσx + ηk 2 yσ y [32][33][34], whereσ x,y label the Pauli matrices operating on the pseudo-spin space, v is the Fermi velocity, and η ≡ 2 /(2m) with m the effective mass. The low-energy dispersion of the normal state is E = 2 v 2 k 2 x + η 2 k 4 y , which is anisotropic in the momentum space. To address the effects of the intrinsic anisotropy on the subgap transport properties, we consider a SDM-based NS junction with the normal direction of NS interface defining as n = (cos θ, sin θ), where θ indicates the interface orientation angle between directions of the x-axis and n, as that carried out in the superconducting hybrids based on monolayer phosphorus [50] and type-II Weyl semimetals [16,17]. In this way, the NS interface is located at x = −y tan θ, as schematically shown in figure 1. For the sake of convenience, we take a coordinate transformation of (x,ỹ) T =R(x, y) T , where the superscript T denotes the transposition of the row vector and the rotation operatorR ≡ e iσ y θ . In the new coordinate system, as shown in figure 1, the N and S regions are placed within the ranges ofx < 0 andx > 0, respectively, and the chemical potentials in the N and S regions are assumed to be tuned independently. We consider the transport properties along thex-direction, and assume that the translational symmetry in theỹ-direction is preserved, so that the transverse momentumk y can be treated as a good quantum number [8][9][10][11][12][13][14][15]. Moreover, under this assumption the influences of boundary effects on the transport can be rationally neglected. In the S region, we take the intra-sublattice/orbit s-wave pairing, as proposed by recent theoretical work [48,49]. In practice, the superconductivity in the S region can be induced by an s-wave superconductor via the proximity effect, as implemented in similar NS junctions based on graphene [7,18,19] and Weyl semimetals [22,23].
Under these lines, the Bogoliubov-de Gennes (BdG) Hamiltonian is given by [48,49] acting on the pseudo-spin ⊗ Nambu space. In the new coordinate system, the single-particle effective Hamiltonian takes the form ofĥ(k) = v(k x cos θ −k y sin θ)σ x + η(k x sin θ +k y cos θ) 2σ y . The chemical potentialμ(x) =μ N Θ(−x) +μ S Θ(x), where Θ(x) is the Heaviside step function andμ S (μ N ) denotes the chemical potential of the S(N) region. In this paper, we assume that the relation ofμ S μ N is satisfied, so that the leakage of Cooper pairs from the S to N regions can be safely neglected [8][9][10][11][12][13][14][15]. In doing so, the pair potential can be effectively modeled by a step function, i.e.,Δ(x) =Δ 0σ0 e iφ Θ(x), with φ the superconducting phase andσ 0 a 2 × 2 identity matrix operating on the pseudo-spin space. In addition, we assume the relation ofΔ 0 μ S holds throughout this work, to ensure the validity of mean-field approximation. We note that the BdG Hamiltonian shown in equation (1) can also be derived from the lattice model of a graphene-like system in proximity to an s-wave superconductor, and the details are given in appendix A.
To rescale related parameters as dimensionless ones, it is convenient to introduce three natural scales, i.e., the momentum scale p 0 ≡ 2mv, the energy scale E 0 ≡ p 2 0 /(2m), and the length scale l 0 ≡ /p 0 [34,35]. In doing so, one can define the dimensionless quantities as Hereafter, we take the proposed dimensionless quantities to express corresponding formulism.
In the N region, by diagonalizing the BdG Hamiltonian given in equation (1), the rescaled dispersion can be written as where ζ e(h) = (k e(h) cos θ − q sin θ) 2 + (k e(h) sin θ + q cos θ) 4 , and the quantity with a subscript e(h) denotes the electron-like (hole-like) one. The group velocity in the N region can be parameterized as As can be seen, the group velocity is θ-dependent and the two components exhibit distinct behaviors with respect to the momentum, reflecting the anisotropic aspect of SDMs. According to equation (2), when the interface orientation angle θ = 0 • , the low-energy excitations of the N region disperse quadratically in the transport direction. While in the case of θ = 0 • , the low-energy excitations exhibit linear dispersions with respect to the longitudinal momentum k e(h) . Since the situations of θ = 0 • and θ = 0 • possess distinct scattering configurations, it is necessary to consider them separately. Without loss of generality, in both situations of θ = 0 • and θ = 0 • we consider the scattering configurations resulting from electron-like propagating incident modes. The scenarios stemming from hole-like incident modes can be analyzed in the same way.

Scattering configuration with θ = 0 •
In the case of θ = 0 • , equation (2) gives four roots of k e(h) for a set of fixed E, μ N , and q, implying that there are four possible electron-like (hole-like) scattering modes in the N region. We note that to get a propagating incident mode, the condition of q ∈ [q C e,1 (θ), q C e,2 (θ)] should be satisfied, where q C e,1 (θ) and q C e,2 (θ) are, respectively, the upper and lower critical transverse momenta of an electron-like mode. Under the condition of q ∈ [q C e,1 (θ), q C e,2 (θ)], the four electron-like modes can be classified into two propagating modes k e,1 with vx e,1 > 0 and k e,2 with vx e,2 < 0, and two evanescent ones κ e,1 with Im(κ e,1 ) < 0 and κ e,2 with Im(κ e,2 ) > 0. The situation for the hole-like mode is more complicated. Defining the range of q to get a propagating hole-like mode as [q C h,1 (θ), q C h,2 (θ)], we have [q C h,1 (θ), q C h,2 (θ)] ⊆ [q C e,1 (θ), q C e,2 (θ)] for E 0 and μ N 0. As a consequence, equation (2) determines two types of k h depending on the value of q. In the case of q ∈ [q C h,1 (θ), q C h,2 (θ)], there are two propagating hole-like modes k h,1 with vx h,1 > 0 and k h,2 with vx h,2 < 0, and two complex ones κ h,1 with Im(κ h,1 ) < 0 and κ h,2 with Im(κ h,2 ) > 0. While the four roots of k h are all complex when q ∈ [q C e,1 (θ), q C e,2 (θ)] and q / ∈ [q C h,1 (θ), q C h,2 (θ)], we denote two of which by k h,2 and κ h,1 , satisfying Im(k h,2 ) < 0 and Im(κ h,1 ) < 0, respectively. In the S region, there are eight evanescent modes in the subgap regime, we choose four of them indicated by k eq,1 (2) and k hq,1 (2) , with Im(k eq,1 (2) ) > 0 and Im(k hq,1 (2) ) > 0, respectively. Although the evanescent modes do not contribute to the current, they need to be included in the wave functions to get correct boundary conditions. With this in mind, the wave function in the N and S regions are, respectively, parameterized as where r ee(he),1 and r ee(he),2 are, respectively, the reflection amplitudes of propagating and evanescent modes in the NR(AR) processes, and the coefficients t eq(hq),1 and t eq(hq),2 represent the transmission amplitudes for the propagating and evanescent electron-like (hole-like) quasiparticles, respectively. The details of related basis scattering states are given in appendix B.1.
In practice, there may exist defects and lattice mismatch at the boundary between the N and S regions, which may profoundly influence the transport properties. To model the effects originating from the interfacial defects and lattice mismatch, we introduce a interfacial barrier modeled byŨ(x) ≡Ũ 0 δ(x), with δ(x) being the Delta function. In doing so, the boundary conditions can be formulated as where M(θ) = ∂ X − Z sin 2 θσ yτ 0 and Z =Ũ 0 E 0 l 0 indicating the normalized interfacial barrier strength. The probabilities for the NR and AR processes are, respectively, defined as with τ z the Pauli matrix operating on the Nambu space. We note that the probabilities of evanescent modes R ee,2 and R he,2 are vanishing and do not contribute to the current. Therefore, in accordance with the Blonder-Tinkham-Klapwijk (BTK) formula [4], the zero-temperature differential conductance along the X-direction can be written as where W is the width along the Y-direction of the junction. It is convenient to normalize the conductance by G 0 (θ) = e 2 W πh [q C e,2 (θ) − q C e,1 (θ)], indicating the conductance of a related SDM-based NN junction in the ballistic limit. In the special case of θ = 90 • , the cut-off of transverse momenta for an electron-like mode are given by q C e,2 (90

Scattering configuration with θ = 0 •
We now turn to the scattering problem in the NS junction with θ = 0 • , where the concerned transport is along the linear dispersion direction. We note that in the current direction, the low-energy excitations in SDMs are similar as that in Dirac materials, thus the scattering amplitudes can be calculated in term of the standard procedure employed in Dirac-material-based NS junctions [8][9][10][11][12][13][14][15], see appendix B.2 for details. For the scattering problem resulting from an electron-like incident mode, in the limit regime of μ S max(μ N , E, Δ 0 ), the reflection amplitudes for the processes of NR and AR are, respectively, given by . As can be seen, both r ee and r he are independent of Z, implying that in the limit of μ S max(μ N , E, Δ 0 ), the transport properties along the linear dispersion direction are insensitive to the interfacial barrier, similar as that in a graphene-based NS junction [12]. However, beyond the limit regime of μ S max(μ N , E, Δ 0 ), the scattering probabilities and hence the differential conductance become sensitive to Z, as shown in figure 5(a).
The zero-temperature differential conductance along the linear dispersion direction can be obtained from equation (7), by substituting R ee,1 and R he,1 , respectively, with In addition, the cut-off of transverse momenta can be reduced into

Results and discussion
In this section, we proceed to analyze the numerical results and reveal manifestations of the anisotropic band structure of SDM in the subgap transport properties. According to the typical values of v and m taken in a series of references [33][34][35][36], we set Δ 0 = E 0 /1000 in the numerical calculation. Furthermore, we assume that the S region is heavily doped to the regime of μ S Δ 0 , and set μ S = 500Δ 0 in the numerical calculation for definiteness.

Orientation-dependent crossover from the retro-AR to specular-AR
As a starting point, we focus on the subgap differential conductance of two special situations with θ = 0 • and θ = 90 • . The most prominent scenario is that the subgap differential conductance in the two special cases exhibits distinct crossover behaviors from the retro-AR to specular-AR, as shown in figures 2(a) and (b), reflecting the intrinsic anisotropy of SDMs. In general, for μ N < Δ 0 , the subgap differential conductance vanishes at eV = μ N due to the absence of AR for any angle of incidence. Since the subgap differential conductance is dominated by the retro-AR and specular-AR, respectively, in the regimes of eV < μ N and eV > μ N , the line of eV = μ N serves as a boundary between the retro-AR and specular-AR in the subgap regime [7][8][9][10][11][12][13][14][15]. In the case of θ = 0 • , when eV = μ N the subgap differential conductance is non-vanishing even in the presence of large momentum-mismatch between the N and S regions (i.e., the large difference between μ N and μ S ), as depicted in figure 2(c). Since the subgap differential conductance only vanishes at eV = μ N , the line of eV = μ N serves as an unambiguous boundary between the retro-AR and specular-AR, as shown in figure 2(a). However, for θ = 90 • , in the regime of μ N < Δ 0 the subgap differential conductance is strongly suppressed by the large momentum-mismatch, so that G(90 • )| eV<Δ 0 G(90 • )| eV=μ N = 0, as illustrated in figure 2(d). Consequently, in the case of θ = 90 • , the boundary between the retro-AR and specular-AR disappears, as shown in figure 2(b).
The anisotropic aspect of the subgap transport also manifests itself in the distinct behaviors of the zero-bias differential conductance for θ = 0 • and θ = 90 • . As shown in figure 2(c), when θ = 0 • the differential conductance always exhibits a zero-bias peak in the regime of μ N Δ 0 . Remarkably, in the case of 0 < μ N Δ 0 , G(0 • ) G 0 (0 • ) | eV=0 takes a universal value of 8 5 (e.g., for μ N /Δ 0 = 0.1, 0.5, 1), implying that the zero-bias differential conductance is insensitive to the momentum-mismatch. On the contrary, as can be seen in figure 3(b), the zero-bias differential conductance strongly depends on the momentum-mismatch for θ = 90 • and G(90 • )| eV=0 almost vanishes in the regime of μ N Δ 0 . Therefore, the behavior of zero-bias differential conductance is strongly orientation-dependent.
To address the unique anisotropy of the subgap differential conductance in SDM-based NS junctions, we compare our results with those in similar NS junctions based on the Dirac materials with isotropic dispersions. Since the excitations in SDMs disperse linearly for θ = 0 • , we compare G(0 • ) with that in NS junctions based on graphene and topological insulators which possess linear low-energy excitations in normal states [11][12][13]. As shown in figure 2(c), although the configuration of G(0 • ) is similar as that in NS junctions based on graphene [11,12] and topological insulators [13], there is a significant difference on the value of G(0 • ) G 0 (0 • ) | eV=0 . For eV = 0, the reflection probabilities in the regime of 0 < μ N μ S can be analytically expressed as R ee | eV=0 = q 4 In terms of the BTK formula [4], we have 8 5 , differing from the value of 4 3 in graphene-and topological-insulator-based NS junctions [11][12][13]. This consequence can be ascribed to the unique band structures of SDMs intermediate between the linear and quadratic spectra. In the case of θ = 90 • , the transport is along the quadratic dispersion direction. Because bilayer graphene harbors quadratic low-energy excitations, we compare G(90 • ) with the differential conductance in bilayer-graphene-based NS junctions. As proposed in references [7][8][9], in bilayer-graphene-based NS junctions with weakly doped N regions, the zero-bias differential conductance keeps finite and the boundary between the retro-AR and specular-AR is unambiguous. In contrast, as charted out in figures 2(b) and (d), in the SDM-based NS junction with μ N Δ 0 , the zero-bias differential conductance G(90 • )| eV=0 almost vanishes and the  crossover behavior from the retro-AR to specular-AR disappears. The novel subgap transport properties in SDM-based NS junctions can be ascribed to the unique pseudo-spin textures of SDMs.
To understand the underlying physics behind the anisotropic subgap transport properties, it is instructive to reveal the pseudo-spin textures of SDMs. To do so, we define the pseudo-spin vector as P ≡ σ ⊗ τ 0 , which is a unit vector consisting of the expectation values of operator σ ⊗ τ 0 in normalized basis scattering states, where σ ≡ (σ x , σ y ) [51][52][53][54][55]. Specifically, the pseudo-spin carried by a propagating electron-like mode can be formulated as The pseudo-spin textures of electron-like conduction band for θ = 0 • and θ = 90 • are schematically shown in figure 3. As can be seen, the pseudo-spin configuration is highly anisotropic, in stark contrast to that of Dirac materials [56][57][58]. In the NS junction with θ = 0 • , when q is small enough, the pseudo-spin carried by an electron-like incident mode is nearly opposite to that of the normally reflected one. Consequently, the conservation of pseudo-spin suppresses the NR, giving rise to an enhanced AR even in the presence of strong momentum-mismatch. Moreover, since the conductance is mainly contributed by the modes with small q, the differential conductance keeps finite for eV = μ N , leading to a clear boundary between the retro-AR and specular-AR. Whereas for θ = 90 • , as depicted in figure 3(b), the incident and normally reflected modes carry the same pseudo-spin, so that the NR is favorable in the presence of large momentum-mismatch. Therefore, in the case of θ = 90 • , the AR and subgap differential conductance are strongly suppressed by the large momentum-mismatch, i.e., G(90 • )| eV<Δ 0 G(90 • )| eV=μ N = 0, thus blurring the boundary between the retro-AR and specular-AR. We now turn to the transport properties in a SDM-based NS junction extending along an arbitrary direction. When the N region is weakly doped satisfying μ N < Δ 0 , as shown in figure 4(a), for θ 75 • , the subgap differential conductance exhibits a clear crossover from retro-AR to specular-AR with increasing the bias-voltage. While for θ 75 • , the boundary between retro-AR and specular-AR is ambiguous. Therefore, the crossover from the retro-AR to specular-AR is blurred with increasing the orientation angle θ. For the NS junction with a heavily doped N region, we find that the subgap differential conductance non-monotonically varies with θ, as depicted in figure 4(b). This signature is quite different from that in similar NS junctions based on type-II Weyl semimetals [16] and monolayer phosphorus [50], where the differential conductance decays monotonically when the orientation angle increases.

Effects of the interfacial barrier
Owing to the intrinsic anisotropy of SDMs, the effects of the interfacial barrier on the differential conductance are orientation-dependent. For θ = 0 • , the pseudo-spin carried by the incident and normally reflected modes are nearly opposite to each other when the incident angle is small enough, resulting in the prohibition of NR and the enhancement of AR. As a consequence, the zero-bias differential conductance periodically oscillates with Z without a decaying profile, as charted out by figure 5(a). In addition, by inspecting figure 5(a) one can find that the oscillation period of Z-dependent G(0 • )| eV=0 increases by dropping the chemical potential μ N . In the limit of μ S max(μ N , eV, Δ 0 ), the oscillation period goes to infinity, so that the differential conductance becomes insensitive to the interfacial barrier. On the contrary, when θ = 0 • the differential conductance exponentially decays with increasing the interfacial barrier strength, as shown in figures 5(b)-(d). In this regard, the influences of the interfacial barrier on the subgap transport are orientation-dependent, in stark contrast to that in superconducting hybrids based on type-II Weyl semimetals. In type-II-Weyl-semimetal-based superconducting hybrids, the differential conductance oscillates with the interfacial barrier strength without a decaying profile, although the low-energy excitations in type-II Weyl semimetals are also anisotropic due to the tilting of Weyl cones [17]. The anisotropic aspect of the Z-dependent differential conductance in SDM-based NS junctions can be ascribed

Conclusion
In conclusion, in the framework of BdG equation, we have figured out the manifestations of the intrinsic anisotropy of SDMs in the AR configurations and subgap differential conductance. For the transport along the linear dispersion direction (θ = 0 • ), the subgap differential conductance exhibits a clear crossover from the retro-AR to specular-AR with increasing the bias-voltage, and the zero-bias differential conductance is insensitive to the momentum mismatch when the N region is weakly doped. Moreover, when the interfacial barrier strength increases, the zero-bias conductance exhibits an oscillating configuration without a decaying profile. However, for the transport along the quadratic dispersion direction (θ = 0 • ), the boundary between the retro-AR and specular-AR fades out with increasing the orientation angle, and the zero-bias differential conductance rapidly drops when the momentum mismatch or the interfacial barrier strength increases. We have elucidated the anisotropic coherent transport properties resorting to the unique pseudo-spin textures of SDMs. These results would provide intriguing insights for the coherent transport in SDMs-based superconducting hybrid structures, and we anticipate more interesting results for the Andreev bound states and supercurrents in SDMs-based Josephson junctions.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.

Appendix A. Derivation of the BdG Hamiltonian
We derive the BdG Hamiltonian from the lattice model of a graphene-like system in proximity to an s-wave superconductor. Physically, the semi-Dirac dispersion can be realized in a graphene-like systems with breaking the hexagonal symmetry [26,32,33]. Thus we take an artificial honeycomb lattice model with anisotropic nearest-neighbor (NN) hopping as a prototype. The NN tight-binding model containing an intra-sublattice/orbit Bardeen-Cooper-Schrieffer superconducting pairing term is given by where the primitive lattice vectors m 1 = a 2 (3, where the parameter γ k is given by with the NN lattice vectors being defined as δ 1 = a 2 (1, 3), and δ 3 = a(−1, 0). Accordingly, the dispersion for the normal state is = ±|γ k |.
For the critical case of t = 2t, vanishes at M = ( 2π 3a , 0), and the normal state reaches the semi-Dirac phase. Now we linearize the Hamiltonian H for a small k around the M point by setting k x → k x + 2π 3a and k y → k y , with k x a 1 and k y a 1 being satisfied. Up to the leading order in k, we arrive at γ M+k 3iak x e 2iπ/3 + 3 4 a 2 k 2 y e 2iπ/3 , where the relation of k = −k has been used. RewriteH q in the form ofH q = q ψ † q H q ψ q , with the basis and By defining v ≡ 3at 2 , η ≡ 3a 2 t 8 , and μ/2 → μ, the upper block can be formulated as where the anticommutation relation has been employed. Taking a unitary transformation H BdG ≡ U † HU with The BdG Hamiltonian can be compactly written as where we rewrite Δ ≡ Δ 0 e iφ with Δ 0 the amplitude of pairing potential and φ the superconducting phase, σ 0 is a 2 × 2 identity matrix, and the Pauli matricesσ x,y andτ x,y,z act on the pseudo-spin and Nambu spaces, respectively.

Appendix B. Calculation of the scattering states in NS junctions based on semi-Dirac materials
In this appendix we present necessary calculation details regarding the wave functions and related quantities in SDMs-based NS junctions.

B.1. Scattering states for θ = 0
We assume the translational symmetry in the Y-direction is preserved, so that the transverse momentum q can be treated as a good quantum number. By solving the rescaled BdG equation H BdG (−i∂ X , q)ϕ = Eϕ, the basis states in the N region are given by 1(2) cos θ − q sin θ) + i(k h,1 (2) where the related quantities k e(h),1 (2) and κ e(h),1 need to be numerically determined in accordance with equation (2) shown in the main text.