Diffraction catastrophes and semiclassical quantum mechanics for Veselago lensing in graphene

TL;DR

This paper investigates how trigonal warping affects electron focusing in graphene n-p junctions, revealing valley-dependent caustics and polarization effects, supported by semiclassical theory and tight-binding simulations.

Contribution

It introduces a semiclassical framework to analyze trigonal warping effects on Veselago lensing, predicting valley-specific focusing and transmission features in graphene.

Findings

Trigonal warping causes cusp caustics and valley-dependent focusing.

Transmission peaks split into two, indicating valley polarization.

Semiclassical predictions match tight-binding simulation results.

Abstract

We study the effect of trigonal warping on the focussing of electrons by n-p junctions in graphene. We find that perfect focussing, which was predicted for massless Dirac fermions, is only preserved for one specific sample orientation. In the general case, trigonal warping leads to the formation of cusp caustics, with a different position of the focus for graphene's two valleys. We develop a semiclassical theory to compute these positions and find very good agreement with tight-binding simulations. Considering the transmission as a function of potential strength, we find that trigonal warping splits the single Dirac peak into two distinct peaks, leading to valley polarization. We obtain the transmission curves from tight-binding simulations and find that they are in very good agreement with the results of a billiard model that incorporates trigonal warping. Furthermore, the positions of the transmission maxima and the scaling of the peak width are accurately predicted by our semiclassical theory. Our semiclassical analysis can easily be carried over to other Dirac materials, which generally have different Fermi surface distortions.

Figures

Equations

$$H_{\alpha}^{\pm}=\pm\Big{(}|\mathbf{p}|+\alpha\mu|\mathbf{p}|^{2}\cos\big{[}3(\phi_{\mathbf{p}}+\theta)\big{]}\Big{)}+U(\mathbf{x}),$$

$$S_{np}(x,y,p_y)=L_1{p}_{x,e,\alpha }(p_y)+x{p}_{x,h,\alpha }(p_y)+y{p}_y,$$

$$x_{\text{cusp},\alpha }=L_1\frac{U_0-E}{E}\left(1-\alpha \frac{4U_0}{3t}\right)\stackrel{U_0=2E}{====}{L}_1-\frac{8\alpha E}{3t}{L}_1.$$

$$U_{0,\text{cusp},\alpha }=E\frac{L_1+L_2}{L_1}\left(1+\alpha \frac{4E}{3t}\frac{L_2}{L_1}\right)\stackrel{L_1=L_2}{====}2E+\frac{8\alpha E^2}{3t}.$$

$$H_{\text{full}}=\left(\begin{array}{cc}0& -tu(\mathbf{k})\\ -tu(\mathbf{k}{)}^{\ast }& 0\end{array}\right),\text{ where }u(\mathbf{k})=\sum _{{\mathbf{\delta }}_i}{e}^{i\mathbf{k}\cdot {\mathbf{\delta }}_i},$$

$${\mathbf{\delta }}_1=a_{CC}\left(0,1\right),{\mathbf{\delta }}_2=\frac{a_{CC}}{2}\left(\sqrt{3},-1\right),{\mathbf{\delta }}_3=\frac{a_{CC}}{2}\left(-\sqrt{3},-1\right).$$

$${\mathbf{v}}_1=a_{CC}\left(\sqrt{3},0\right),{\mathbf{v}}_2=\frac{a_{CC}}{2}\left(\sqrt{3},3\right).$$

$$\mathbf{K}=\frac{2\pi }{3\sqrt{3}{a}_{CC}}\left(-1,\sqrt{3}\right),{\mathbf{K}}^{\prime }=\frac{2\pi }{3\sqrt{3}{a}_{CC}}\left(1,\sqrt{3}\right).$$

$$H_K=\left(\begin{array}{cc}0& \frac{3}{2}t{a}_{CC}(k_x-i{k}_y)-\frac{3}{8}t{a}_{CC}^2(k_x+i{k}_y{)}^2\\ \frac{3}{2}t{a}_{CC}(k_x+i{k}_y)-\frac{3}{8}t{a}_{CC}^2(k_x-i{k}_y{)}^2& 0\end{array}\right).$$

$$H_{K^{\prime }}=\left(\begin{array}{cc}0& \frac{3}{2}t{a}_{CC}(k_x+i{k}_y)+\frac{3}{8}t{a}_{CC}^2(k_x-i{k}_y{)}^2\\ \frac{3}{2}t{a}_{CC}(k_x-i{k}_y)+\frac{3}{8}t{a}_{CC}^2(k_x+i{k}_y{)}^2& 0\end{array}\right).$$

$$E_{\alpha }^{\pm }=\pm \left(\frac{3}{2}t{a}_{CC}|\mathbf{k}|+\alpha \frac{3}{8}t{a}_{CC}^2|\mathbf{k}{|}^2\cos (3{\varphi }_{\mathbf{k}})\right),$$

$$E_{\alpha }^{\pm }=\pm \left(|\mathbf{p}|+\alpha \mu |\mathbf{p}{|}^2\cos (3{\varphi }_{\mathbf{p}})\right).$$

$$E_{\alpha}^{\pm}=\pm\Big{(}|\mathbf{p}|+\alpha\mu|\mathbf{p}|^{2}\cos\big{[}3(\phi_{\mathbf{p}}+\theta)\big{]}\Big{)}.$$

$$H_{\alpha}^{\pm}=\pm\Big{(}|\mathbf{p}|+\alpha\mu|\mathbf{p}|^{2}\cos\big{[}3(\phi_{\mathbf{p}}+\theta)\big{]}\Big{)}+U(\mathbf{x}),$$

$$U(\mathbf{x})=U_0\Theta (x),$$

$$S_{np}(x,y,p_y)=L_1{p}_{x,e,\alpha }(p_y)+x{p}_{x,h,\alpha }(p_y)+y{p}_y.$$

$$y=-L_1\frac{\partial p_{x,e,\alpha }}{\partial p_y}-x\frac{\partial p_{x,h,\alpha }}{\partial p_y}.$$

$$x_{\text{cst}}=-L_1\frac{{\partial }^2{p}_{x,e,\alpha }/\partial p_y^2}{{\partial }^2{p}_{x,h,\alpha }/\partial p_y^2},$$

$$\frac{\partial E_{\alpha }^{\pm }}{\partial p_x}\frac{\partial p_x}{\partial p_y}+\frac{\partial E_{\alpha }^{\pm }}{\partial p_y}=0.$$

$$\frac{{\partial }^2{p}_x}{\partial p_y^2}=-\frac{{\partial }^2{E}_{\alpha }^{\pm }/\partial p_y^2}{\partial E_{\alpha }^{\pm }/\partial p_x},$$

$$\frac{\partial E_{\alpha }^{\pm }}{\partial p_x}=\pm \frac{p_x+2\alpha \mu p_x^2}{|p_x|},\frac{{\partial }^2{E}_{\alpha }^{\pm }}{\partial p_y^2}=\pm \frac{1-7\alpha \mu p_x}{|p_x|},\frac{{\partial }^2{p}_x}{\partial p_y^2}=-\frac{1-7\alpha \mu p_x}{p_x+2\alpha \mu p_x^2}.$$

$$p_{x,e,K}=\frac{1}{2\mu }\left(1-\sqrt{1-4\mu E}\right),p_{x,e,K^{\prime }}=\frac{1}{2\mu }\left(-1+\sqrt{1+4\mu E}\right),$$

$$p_{x,h,K}=\frac{1}{2\mu }\left(1-\sqrt{1+4\mu (U_0-E)}\right),p_{x,h,K^{\prime }}=\frac{1}{2\mu }\left(-1+\sqrt{1-4\mu (U_0-E)}\right)$$

$$x_{\text{cusp},\alpha }=L_1\frac{U_0-E}{E}\left(1-8\alpha \mu U_0\right)+\mathcal{O}({\mu }^2).$$

$$\Delta x_{\text{cusp}}=x_{\text{cusp},K}-x_{\text{cusp},K^{\prime }}=16\mu L_1{U}_0\frac{U_0-E}{E}+\mathcal{O}({\mu }^2).$$

$$x_{\text{cusp},\alpha }=L_1\frac{U_0-E}{E}\left(1-\alpha \frac{4U_0}{3t}\right).$$

$$\Delta x_{\text{cusp}}=x_{\text{cusp},K}-x_{\text{cusp},K^{\prime }}=L_1\frac{8U_0}{3t}\frac{U_0-E}{E}.$$

$$\Psi (x,y)={\int }_{-\infty }^{\infty }\text{d}{p}_{y}f(x,y,p_y)e^{i{S}_{np}(x,y,p_y)/h}.$$

$$\frac{\partial S_{np}}{\partial p_y}=\frac{{\partial }^2{S}_{np}}{\partial p_y^2}=\frac{{\partial }^3{S}_{np}}{\partial p_y^3}=0.$$

$$S_{np}(x,y,p_y)=S^{(4)}(x,y,p_y)+\mathcal{O}({\beta }^5)=q_0(x,y)+q_1(x,y)\beta +\frac{q_2(x,y)}{2}{\beta }^2+\frac{q_3(x,y)}{6}{\beta }^3+\frac{q_4(x,y)}{24}{\beta }^4+\mathcal{O}({\beta }^5),$$

Full text

Veselago lenses Veselago (1968) are special types of lenses, which are made of materials with a negative refractive index. Such lenses can overcome the diffraction limit Pendry (2000) and can nowadays be realized in metamaterials Smith et al. (2000); Houck et al. (2003); Grbic and Eleftheriades (2004), chiral metamaterials Tretyakov et al. (2003); Pendry (2004); Zhang et al. (2009); Xiong et al. (2010) and photonic crystals Parimi et al. (2004); Cubukcu et al. (2003). An electronic analog of a Veselago lens can be created using n-p junctions, with the classical trajectories of the charge carriers playing the role of the rays in geometrical optics. Such junctions focus electrons, because the group velocity for holes is in the direction opposite to their phase velocity, whereas the two velocities are in the same direction for electrons. However, in conventional semiconductors, such interfaces are unsuitable because of their high reflectivity, owing to the presence of a depletion region.

Cheianov et al. Cheianov et al. (2007) realized that graphene does not have this drawback. It has zero bandgap, as the valence and conduction bands touch at two non-equivalent corners of the Brillouin zone, known as $K$ and $K^{\prime}$. The low-energy charge carriers are ballistic over large distances and follow the Dirac equation Wallace (1947); McClure (1956); Slonczewski and Weiss (1958); Semenoff (1984); Castro Neto et al. (2009); Katsnelson (2013). This gives rise to Klein tunneling: normally incident electrons are transmitted with unit probability Klein (1929); Katsnelson et al. (2006); Cheianov and Fal’ko (2006); Shytov et al. (2008); Tudorovskiy et al. (2012); Reijnders et al. (2013); Young and Kim (2009); Stander et al. (2009); Chen et al. (2016), which makes the interface exceptionally transparent. Recently, Veselago lensing in graphene was experimentally observed Lee et al. (2015); Chen et al. (2016).

Theoretical studies have used the Dirac equation to investigate focussing by flat Cheianov et al. (2007); Reijnders and Katsnelson (2017) and circular Cserti et al. (2007); Péterfalvi et al. (2010); Wu and Fogler (2014) junctions and in zigzag nanoribbons Choi et al. (2014). It was shown that when the electron and hole charge carrier densities are equal, a flat interface is able to focus all trajectories into a single point Cheianov et al. (2007). According to catastrophe theory Berry and Upstill (1980); Poston and Stewart (1978); Arnold et al. (1982); Arnold (1975), such a situation is exceptional: any perturbation of the Hamiltonian will ruin this ideal focus. Indeed, in tight-binding simulations of an n-p junction Milovanović et al. (2015), where transmission was studied as a function of potential strength, a direction-dependent broadening of the transmission peak was observed. However, the precise origin of this broadening was not clarified.

An important correction to the Dirac Hamiltonian is the second-order term in the expansion of graphene’s tight-binding Hamiltonian, known as the trigonal warping term Ajiki and Ando (1996); Ando et al. (1998); Katsnelson (2013). When adding this term, the Hamiltonian becomes dependent on the crystallographic direction of the sample and becomes different for the two valleys, leading to different classical trajectories. Using this principle, a valley beam splitter based on an n-p-n’ junction was devised Garcia-Pomar et al. (2008), in which the trajectories in the $K$- and $K^{\prime}$-valleys are deflected in different angular directions. Creating valley polarization Rycerz et al. (2007); Xiao et al. (2007); Gorbachev et al. (2014) is important for graphene valleytronics applications, where the valley index is used to encode information in a way similar to spintronics.

In this paper, we develop a complete semiclassical theory of Veselago lensing by an n-p junction in the presence of trigonal warping. We show that the ideal focus generally disappears and analyze the different diffraction catastrophes, known as caustics Berry and Upstill (1980); Poston and Stewart (1978); Arnold et al. (1982); Arnold (1975), that are formed. We show that electrons from the $K$- and $K^{\prime}$-valleys are generally focussed at different positions, even at moderate energies. We obtain these positions using semiclassical methods Connor and Farrelly (1981a); Reijnders and Katsnelson (2017); S. Yu. Dobrokhotov et al. (2014) based on the Pearcey function Pearcey (1946); Connor and Farrelly (1981b); Connor and Curtis (1982) and find very good agreement with tight-binding simulations using the Kwant code Groth et al. (2014). We also show that an initial sublattice polarization leads to tilting of the focus. Furthermore, we obtain the transmission as a function of the hole carrier density from both tight-binding simulations Groth et al. (2014) and a semiclassical billiard model Beenakker and van Houten (1989); Milovanović et al. (2013) that incorporates trigonal warping. We find that trigonal warping explains the previously observed peak broadening Milovanović et al. (2015) and that our semiclassical theory accurately predicts the positions of the transmission maxima and the scaling of the peak widths. Since the transmission maximum occurs at a different potential for the two valleys, we find that a graphene n-p junction can act as a valley filter and that one can manipulate the polarization by changing the potential strength. This is similar to the situation in chiral metamaterials Tretyakov et al. (2003); Pendry (2004); Zhang et al. (2009); Xiong et al. (2010), where one can select a certain polarization by tuning the frequency of the incident light.

Trigonal warping adds terms that are quadratic in momentum to the linear term in the massless Dirac Hamiltonian Ajiki and Ando (1996); Ando et al. (1998); Katsnelson (2013). We make all terms in the Hamiltonian dimensionless Reijnders et al. (2013) by scaling energies with the electron energy $E$ and lengths with $L_{1}$, the distance from the source to the junction, see Fig. 1(a). Using first-order perturbation theory, we obtain the classical Hamiltonian in dimensionless variables Ajiki and Ando (1996); Ando et al. (1998); Katsnelson (2013); Berlyand and S. Yu. Dobrokhotov (1987); Belov et al. (2006):

[TABLE]

where $\mu=E/6t$, with $t$ the nearest-neighbor hopping, indicates the relative importance of the quadratic term and $\alpha$ equals +$1$ ($-1$) for valley $K^{\prime}$ ($K$). The angle $\phi_{\mathbf{p}}=\arctan(p_{y}/p_{x})$ and $\theta$ determines the orientation of the sample, with $\theta=0$ corresponding to zigzag edges along the $x$-direction. The dimensionless semiclassical parameter $h=3ta_{CC}/2EL_{1}$, with $a_{CC}$ the carbon-carbon distance.

We consider a sharp, one-dimensional junction: $U(\mathbf{x})=U_{0}\Theta(x)$, with $\Theta(x)$ the Heaviside step function. Then $p_{y}$ is conserved and the classical action for an electron with energy $E$ emitted from $(-L_{1},0)$ equals Cheianov et al. (2007); Reijnders and Katsnelson (2017)

[TABLE]

with $p_{x,e,\alpha}$ ($p_{x,h,\alpha}$) the longitudinal momentum in electron (hole) region, i.e. $H_{\alpha}^{+}(p_{x,e,\alpha}(p_{y}),p_{y})=E$. The classical trajectories, i.e. the points where $\partial S_{np}/\partial p_{y}=0$, and the caustics, where $\partial^{2}S_{np}/\partial p_{y}^{2}$ vanishes as well, are plotted in Fig. 1 for various parameters. The Dirac Hamiltonian, $\mu=0$ in Eq. (1), is symmetric in $p_{x}$ and $p_{y}$. When $U_{0}=2E$, this leads to an ideal focus Cheianov et al. (2007), at which all derivatives of $S_{np}$ with respect to $p_{y}$ vanish Reijnders and Katsnelson (2017). This generally changes when we include trigonal warping. For $\theta=0$, the symmetry in $p_{y}$ is preserved, leading to reflection symmetry in the $x$-axis. However, the symmetry in $p_{x}$ is broken and we obtain cusp caustics, shown in Fig. 1(c) and (d). Because we also break the symmetry between the valleys $K$ and $K^{\prime}$, their cusp points, at which $\partial^{4}S_{np}/\partial p_{y}^{4}\neq 0$, are at different positions on the $x$-axis. Tuning the potential, we obtain sections of the butterfly caustic, see Fig. 1(e), and eventually pass through the butterfly singularity, at which $\partial^{4}S_{np}/\partial p_{y}^{4}=0$, but $\partial^{6}S_{np}/\partial p_{y}^{6}\neq 0$ Reijnders and Katsnelson (2017); Poston and Stewart (1978). For generic $\theta$, Fig. 1(f), both symmetries are broken and the cusp point is no longer on the $x$-axis. Only when $\theta=\pi/6$ (armchair edges), the symmetry in $p_{x}$ is restored and we recover an ideal focus at $U_{0}=2E$.

From here on, we consider zigzag edges along the $x$-direction ($\theta=0$), unless otherwise indicated, as they illustrate the generic situation. We obtain an expression for $x_{\text{cusp},\alpha}$ by solving $\partial^{2}S_{np}/\partial p_{y}^{2}=0$ for $x$, and setting $p_{y}=0$ Reijnders and Katsnelson (2017). Expanding the result up to second order in $\alpha\mu$, we find, in units with dimensions

[TABLE]

Hence, the cusp point for the valley $K$ ($K^{\prime}$) is always to the right (left) of the focus for the Dirac Hamiltonian, given by the first term. Although the above expansion is only sufficient for low energies, it clearly indicates that the effect is sizeable.

We study these effects by performing tight-binding simulations with Kwant Groth et al. (2014). We use the setup of Ref. Milovanović et al. (2015), shown in Fig. 1(a), where current enters the sample through an injector lead and is able to exit through a collector lead and through drain leads on each side Sup . Considering large $L_{2}$ and $W_{c}=W$, we compute the sample wavefunction. Figure 2(a)–(c) shows the resulting density, averaged over sublattices and summed over lead modes Reijnders and Katsnelson (2017); Brey and Fertig (2006). For narrow leads, i.e. large $h_{\text{lead}}=\hbar v_{F}/EW_{i}$, we observe cusp caustics in the valleys $K$ and $K^{\prime}$, see Fig. 2(a). By adjusting $E$ and $U_{0}$, we also observe sections of the butterfly caustic, see Fig. 2(b). Figure 2(c) shows that we obtain a bright focussing spot for wide leads (small $h_{\text{lead}}$), as predicted in Ref. Reijnders and Katsnelson (2017). Subsequently, we fit a Gaussian to a subset of the points on the $x$-axis. Averaging over subsets, we obtain the position of the maximum and an error estimate.

Within the semiclassical approximation, the wavefunction near a cusp can be constructed in terms of the Pearcey function Connor and Farrelly (1981a); Reijnders and Katsnelson (2017); S. Yu. Dobrokhotov et al. (2014); Sup . We find the parameters for this function by expanding the action (2) up to fourth order in $p_{y}$. Figure 2(d)–(f) shows the predictions $x_{\text{max},\alpha}$, obtained from maximizing this wavefunction, as a function of electron energy, sample orientation and $L_{1}$. They are in very good agreement with the positions extracted from the simulations. Note that the approximation fails for low energies, or near $\theta=\pi/6$, as we are too close to the ideal focus, at which $\partial^{4}S_{np}/\partial p_{y}^{4}=0$. However, in these cases, $x_{\text{cusp},\alpha}$ provides a good estimate.

Experimentally, one typically measures the transmission $T$ as a function of the potential strength. We therefore fix $L_{1}$ and $L_{2}$, set $W_{c}=W_{i}$ and compute the intravalley scattering for varying $U_{0}$ using Kwant. We compare these results with a semiclassical billiard model Beenakker and van Houten (1989); Milovanović et al. (2013) that incorporates trigonal warping in the trajectories and has uniformly distributed initial positions and transversal momenta. Figure 3 shows that the agreement between the two methods is very good. Note that the tight-binding transmission is slightly higher than the billiard result on the right (left) of the peak for valley $K$ ($K^{\prime}$), most likely because of interference. Higher order effects and intervalley scattering only mildly influence the results. Figure 3(e) shows that, for small $h_{\text{lead}}$, passing through the butterfly singularity does not significantly enhance the transmission.

Valley splitting is observed at all simulated energies and the total transmission peak is broadened with respect to the result of a Dirac billiard model, which explains earlier observations Milovanović et al. (2015). The transmission maximum for the $K$-valley is always at a lower potential. Setting $x_{\text{cusp},\alpha}=L_{2}$ in Eq. (3) and solving for $U_{0}$, we obtain, up to $\mathcal{O}(\mu^{2})$, in units with dimensions,

[TABLE]

We also obtain predictions $U_{0,\text{max},\alpha}$ for the positions of the maxima from our semiclassical theory. Solving $x_{\text{max},\alpha}=L_{2}$, we obtain values that are in excellent agreement with the simulated peak positions, see Fig. 3(g). As $\partial^{4}S_{np}/\partial p_{y}^{4}$ is lower at the $K$-peak, our theory also correctly predicts that it has lower maximal transmission.

The valley polarization can be increased by going to higher energies, which increases the peak splitting. Alternatively, one can increase $L_{2}$, see Fig. 3(d). Finally, one can also increase both $L_{1}$ and $L_{2}$, see Fig. 3(f), which decreases the peak width $\Delta U_{0}$. We obtain a semiclassical prediction for the scaling behavior of $\Delta U_{0}$ by performing a Taylor expansion of the equation $x_{\text{max},\alpha}=L_{2}$ around $U_{0,\text{cusp},\alpha}$. For $L_{1}=L_{2}=L$, we find that $\Delta U_{0}\propto L^{-1/2}$ in leading order, which is in perfect agreement with the billiard results.

Finally, we briefly discuss symmetry breaking due to an initial sublattice polarization Reijnders and Katsnelson (2017). When we add a constant mass $m_{\text{lead}}$ in the lead, we observe tilting of the main focus, see Fig. 4. Inspecting the Fourier transform, we observe that it is not symmetric in $p_{y}$, verifying earlier predictions Reijnders and Katsnelson (2017). To certify that this effect can be attributed to the Dirac character of the electrons, and not to trigonal warping, we extract the wavefunction at $x=-L_{1}$ from Kwant and compute its evolution using the continuum Dirac equations Reijnders and Katsnelson (2017). This reproduces the tilting, but changes the position of the focus, confirming our hypothesis.

The symmetry breaking can be quantified by setting $W_{c}=W$ and subsequently splitting the collector lead at $y=0$. We then extract the relative transmission $\eta_{b,\text{out},\alpha}$ through the bottom ($y<0$) collector lead and subtract the relative intensity $\eta_{b,\text{in},\alpha}$ located on the sites with negative $y$ for the incoming modes. Figure 4(c) shows that the amount of symmetry breaking strongly depends on $L_{2}$, as expected for a tilted focus. In agreement with predictions made for the Green’s function Reijnders and Katsnelson (2017), it increases with increasing mass, decreases with increasing energy, and changes sign when we change the sign of the mass. Simulating transmission as a function of potential strength, initial polarization seems to slightly decrease the valley polarization.

In short, we have developed a complete semiclassical theory of Veselago lensing in the presence of trigonal warping. We have observed both cusp and butterfly singularities, which provide an interesting relation between the physics of graphene and catastrophe theory Berry and Upstill (1980); Poston and Stewart (1978); Arnold et al. (1982); Arnold (1975). Our theory shows excellent agreement with tight-binding simulations and explains previously observed peak broadening Milovanović et al. (2015). We believe that our predictions could be experimentally verified, using, for instance, a transverse magnetic focussing setup Chen et al. (2016). Although this would require rather high energies, these have been experimentally realized Wehling et al. (2008); Nair et al. (2013); Lee et al. (2016); Khademi et al. (2016). We remark that a smooth junction will lead to additional peak broadening Phong and Kong (2016); Milovanović et al. (2015); Chen et al. (2016), but this makes a theoretical study much more complicated Reijnders et al. (2013).

We emphasize that our analysis of the disintegration of the ideal focus does not fundamentally depend on the fact that the Fermi surface distortion is trigonal: any distortion that breaks the symmetry in $p_{x}$ will lead to the formation of cusp caustics. Therefore, our semiclassical analysis can easily be carried over to other Dirac materials, such as topological insulators Hasan and Kane (2010); Qi and Zhang (2011). Since these generally exhibit stronger band bending, the deviations from the Dirac behavior will be much stronger and therefore visible at lower energies.

S1 Supplemental material

S2 Detailed derivations

For the convenience of the reader, this section contains the details of the derivations that have been discussed in the main text. The semiclassical theory that we develop is based on methods that have been well-established in the literature. For the Dirac Hamiltonian, a similar theory was presented in Ref. Reijnders and Katsnelson (2017).

S2.1 Hamiltonian with trigonal warping

In reciprocal space, the full tight-binding Hamiltonian for graphene is given by (see e.g. Ref. Katsnelson (2013))

[TABLE]

where the nearest-neighbor vectors $\bm{\delta}_{i}$ are given by

[TABLE]

Furthermore, we choose our lattice vectors $\mathbf{v}_{i}$ as

[TABLE]

The electron and hole bands touch at two non-equivalent Dirac points in the Brillouin zone, named $K$ and $K^{\prime}$:

[TABLE]

Expanding the Hamiltonian (S1) up to second order in $\mathbf{k}$ around the $K$-point, we obtain Ajiki and Ando (1996); Ando et al. (1998); Katsnelson (2013)

[TABLE]

Keeping only the linear term, we obtain the Dirac Hamiltonian. The quadratic term is the trigonal warping correction. Expanding the Hamiltonian (S1) around the $K^{\prime}$-point, we find

[TABLE]

Considering the trigonal warping term as a perturbation to the Dirac Hamiltonian, we compute the first order correction to the dispersion relation using perturbation theory and obtain Ajiki and Ando (1996); Ando et al. (1998); Katsnelson (2013)

[TABLE]

where the valley index $\alpha$ equals $-1$ for the $K$-valley and $+1$ for the $K^{\prime}$-valley. Furthermore, $\phi_{\mathbf{k}}$ is the angle in $k$-space, defined by $\tan\phi_{\mathbf{k}}=k_{y}/k_{x}$. Let us now consider an electron with energy $E$ and let $L_{1}$ be the typical length scale of the problem, in our case the distance from the injector lead to the junction interface. We then define the dimensionless variables $\tilde{E}_{\alpha}^{\pm}=E_{\alpha}^{\pm}/E$, $\tilde{\mathbf{x}}=\mathbf{x}/L_{1}$, $\tilde{\mathbf{k}}=\mathbf{k}L_{1}$ and the dimensionless semiclassical parameter $h=3ta_{CC}/(2EL_{1})$. The dimensionless coefficient $\mu=E/(6t)$ controls the relative importance of the second term and we require $\mu^{2}$ to be small in order for our perturbative expansion to be valid. Introducing $\tilde{\mathbf{p}}=h\tilde{\mathbf{k}}$ and omitting the tildes, we obtain

[TABLE]

We remark that $\phi_{\mathbf{k}}=\phi_{\mathbf{\tilde{k}}}=\phi_{\mathbf{\tilde{p}}}$. We will, from here on, omit the tildes and exclusively work with these dimensionless variables, unless otherwise indicated.

Let us now consider a passive rotation by an angle $\theta$, which means that we rotate the coordinate system over an angle $\theta$. Then the nearest-neighbor vectors $\bm{\delta}_{i}$ should be multiplied (from the left) by the matrix of a rotation over an angle $-\theta$. Because of this, the expressions for the Dirac points should be multiplied by the same factor. In turn, this leads to an additional factor $e^{3i\theta}$ in front of the quadratic term in the expansion of $u(\mathbf{k})$ around the $K$-point, and to an additional factor $e^{-3i\theta}$ in its expansion around the $K^{\prime}$-point. Because of these additional factors, our dimensionless expression for the energy then becomes

[TABLE]

Using standard semiclassical arguments Berlyand and S. Yu. Dobrokhotov (1987); Belov et al. (2006), one can then show that the classical Hamiltonian functions for electrons and holes are given by

[TABLE]

It is important to note that this Hamiltonian is no longer invariant with respect to the orientation of the sample, as indicated by the presence of the orientation angle $\theta$. In our choice of coordinates, $\theta=0$ corresponds to zigzag edges along the $x$-axis and $\theta=\pi/6$ corresponds to armchair edges along the $x$-axis. Furthermore, we remark that the lattice is invariant under a $\pi/3$ rotation, but that such a rotation interchanges the sublattices $A$ and $B$ and the valleys $K$ and $K^{\prime}$.

S2.2 Caustics

We consider an n-p junction with the one-dimensional potential

[TABLE]

where $\Theta(x)$ is the Heaviside step function. As the potential is one-dimensional, the transversal momentum $p_{y}$ is a conserved quantity. For $x<0$, the charge carriers are electrons with longitudinal momentum $p_{x,e,\alpha}$ that satisfies $E_{\alpha}^{+}(p_{x,e,\alpha},p_{y})=E$. For $x>0$, we have holes with longitudinal momentum $p_{x,h,\alpha}$ that satisfies $E_{\alpha}^{-}(p_{x,h,\alpha},p_{y})+U_{0}=E$.

Let us study electrons that are emitted by a point source at $(-L_{1},0)$ and are incident on this potential. It easily follows that the classical action is given by Cheianov et al. (2007); Reijnders and Katsnelson (2017)

[TABLE]

The trajectories, which are the stationary points of the action, i.e. the points where $\partial S_{np}/\partial p_{y}=0$, are given by

[TABLE]

The caustic is the set of points $(x_{\text{cst}},y_{\text{cst}})$ where the second derivative $\partial^{2}S_{np}/\partial p_{y}^{2}$ vanishes as well. Differentiating Eq. (S13), we obtain

[TABLE]

and we can obtain $y_{\text{cst}}$ by entering this value back into Eq. (S13). The cusp point $(x_{\text{cusp}},y_{\text{cusp}})$ is the point on the caustic where the third derivative $\partial^{3}S_{np}/\partial p_{y}^{3}$ is identically zero.

Let us first look at the Dirac case for a moment. Setting $\mu=0$ in Eq. (S9), we see that we have rotation symmetry in momentum space. In particular, we have reflection symmetry in $p_{x}$ and in $p_{y}$. As $p_{y}$ is conserved at the barrier interface and the group velocity for holes is antiparallel to their momentum, symmetry in $p_{x}$ means that $p_{x,h,\alpha}=-p_{x,e,\alpha}$ when $U_{0}=2E$. This leads to an ideal focus, as the electrons are focussed in the mirror image of the source Cheianov et al. (2007); Reijnders and Katsnelson (2017). Furthermore, the symmetry $p_{y}\to-p_{y}$ ensures that the classical trajectories are also symmetric in the $x$-axis.

When we include trigonal warping, the dispersion relation becomes dependent on the sample orientation $\theta$. For $\theta=0$, which means that we have zigzag edges along the $x$-axis, the symmetry in $p_{x}$ is broken, but the symmetry in $p_{y}$ is preserved. Therefore, we no longer have an ideal focus at $U_{0}=2E$, but the trajectories are still symmetric in the $x$-axis. On the other hand, when $\theta=\pi/6$, we do have symmetry in $p_{x}$, but the symmetry in $p_{y}$ is broken. Indeed, we immediately see that $\cos(3[\phi_{p}+\pi/6])=-\sin(3\phi_{p})$ and setting $\bar{\mathbf{p}}=(-p_{x},p_{y})$, we observe that $\sin(3\phi_{\bar{\mathbf{p}}})=\sin(3[\pi-\phi_{\mathbf{p}}])=\sin(3\phi_{\mathbf{p}})$. Therefore, we have an ideal focus at $U_{0}=2E$ when we have armchair edges along the $x$-axis. However, the trajectories are no longer symmetric in the $x$-axis. For generic $\theta$, both symmetries are broken and there will be no ideal focus at $U_{0}=2E$.

Let us now take a closer look at $\theta=0$, as the symmetry in $p_{y}$ considerably simplifies the calculations. Since $E_{\alpha}^{\pm}(p_{x},-p_{y})=E_{\alpha}^{\pm}(p_{x},p_{y})$, we also have $p_{x,e,\alpha}(p_{y})=p_{x,e,\alpha}(-p_{y})$, and the same identity holds for $p_{x,h,\alpha}$. Therefore, their second derivatives are symmetric with respect to $p_{y}$, which implies that $x_{\text{cst}}$ is. Hence, we infer from Eq. (S13) that the caustic is symmetric in the $x$-axis, which, in particular, means that the cusp points are. As the third derivative $\partial^{3}S_{np}/\partial p_{y}^{3}$ is antisymmetric in $p_{y}$, it vanishes at $p_{y}=0$ and there is always a cusp point on the $x$-axis. In what follows, we investigate the position $x_{\text{cusp}}$ of this particular cusp point. We first note that when we expand $S_{np}$ at $(x_{\text{cusp}},0)$ in $p_{y}$ around $p_{y}=0$, the terms of odd order vanish. Let us now determine the various deratives of $p_{x}$ with respect to $p_{y}$ by differentiating the energy relation (S9). Fixing the energy $E$ and considering $p_{x}$ to be a function of $p_{y}$, we obtain

[TABLE]

Since we consider $\theta=0$, the second term vanishes at $p_{y}=0$, and we find that $\partial p_{x}/\partial p_{y}$ equals zero at that point. Of course, we could already have concluded this from its antisymmetry in $p_{y}$. Taking the derivative of the above result, and using that the first derivative of $p_{x}$ vanishes, we obtain

[TABLE]

where all derivatives are to be evaluated at zero $p_{y}$. Evaluating the energy derivatives using Eq. (S10), we find

[TABLE]

We note that the last expression does not contain a $\pm$-sign and is therefore the same for both electrons and holes.

To find the position of the cusp, we need to know the momenta $p_{x,e,\alpha}$ and $p_{x,h,\alpha}$ of the right-moving charge carriers. For electrons, for which me need to take the plus sign in Eq. (S10), we set $U_{0}=0$, require that $p_{x}>0$ and set $\phi_{\mathbf{p}}=0$. We then obtain

[TABLE]

where we have used that for low energies the leading term in the Taylor expansion should give $p_{x,e,\alpha}=E$. For holes, we take the minus sign in Eq. (S10) and set $U=U_{0}>E$, $p_{x}<0$ and $\phi_{\mathbf{p}}=\pi$. Requiring that $p_{x,h,\alpha}=U_{0}-E$ in the zeroth order expansion in $\mu$, we have

[TABLE]

Inserting the results (S18) and (S19) into Eq. (S14), we obtain an expression for $x_{\text{cusp},\alpha}$. This expression is rather cumbersome, but can be easily implemented into a computer algebra system. It is used to generate the cusp lines in Fig. 2(d)–(f).

To gain more insight into the splitting, we perform a first order Taylor expansion in $\mu$ and obtain

[TABLE]

For small $\mu$, we therefore find the distance between the two cusp points as

[TABLE]

Going back to units with dimensions, we have, up to second order in $\mu$,

[TABLE]

and the distance between the two cusp points equals

[TABLE]

The implications of these equations have been extensively discussed in the main text.

S2.3 Position of the maximum

Within the Dirac approximation, one can construct an analytical solution for scattering by a sharp n-p junction Cheianov et al. (2007); Reijnders and Katsnelson (2017). The wavefunction in the hole region can be expressed as an integral over the transversal momentum $p_{y}$, which labels the trajectories:

[TABLE]

The action $S_{np}$ is given by Eq. (S12), with $\mu=0$ for the Dirac case, and the amplitude function $f(x,y,p_{y})$ depends on the details of the source Reijnders and Katsnelson (2017). This expression is not only valid for the Green’s function, but also for current injection through a lead.

Based on general, partially semiclassical, principles, expression (S24) should still be valid when we include trigonal warping. However, in this case, we have $\mu\neq 0$ in the action $S_{np}$ and the amplitude $f(x,y,p_{y})$ is likely a more complicated function. Within the semiclassical approximation, which is valid for small $h$, we can compute the maximum of the wavefunction (S24) near a cusp point using the Pearcey approximation Connor and Farrelly (1981a); Reijnders and Katsnelson (2017); S. Yu. Dobrokhotov et al. (2014). Although this approximation typically predicts too large values for the wavefunction, it accurately predicts the position of the maximum, even for relatively large values of $h$ Reijnders and Katsnelson (2017). A key observation is that, within this approximation, the position of the maximum does not depend on the function $f(x,y,p_{y})$, but only on the action $S_{np}$. When we include trigonal warping, we can therefore obtain a reliable result for the position of the main focus without solving the full problem. However, this approximation is only valid when we consider sufficiently large values of $a_{4}$. Otherwise, we are dealing with a higher order singularity, and the wavefunction should be approximated using a different special function.

Before applying the Pearcey approximation to our problem, let us briefly review it and introduce the relevant notations. The main contribution to the integral (S24) is given by the critical points, where $\partial S_{np}/\partial p_{y}$ vanishes. This gives rise to the classical trajectories, which we discussed in the previous subsection. At the caustic, the second derivative also vanishes, and at the cusp the third derivative is also identically zero:

[TABLE]

We obtain the simplest approximation to the wavefunction (S24) near a cusp point by making a Taylor expansion of the action $S_{np}$ around the critical point $(x_{\text{cusp}},y_{\text{cusp}},p_{y,0})$ up to the first nonzero term. This gives

[TABLE]

where $\beta=p_{y}-p_{y,0}$. Because of Eq. (S25), the expansion coefficients satisfy:

[TABLE]

where $\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}=(x,y)-(x_{\text{cusp}},y_{\text{cusp}})$.

The leading order approximation to the wavefunction can then be expressed Connor and Farrelly (1981a); Reijnders and Katsnelson (2017); S. Yu. Dobrokhotov et al. (2014) in terms of the Pearcey function $\text{P}^{\pm}(x)$, which is defined by the integral

[TABLE]

where the superscript plus or minus corresponds to the sign in front of the $t^{4}$ term. This function possesses two important symmetries Connor and Farrelly (1981b); Connor and Curtis (1982), which simplify computations:

[TABLE]

These can be verified directly using definition (S28). Furthermore,

[TABLE]

where $\text{P}_{v}^{\pm}(u,v)$ is the partial derivative with respect to the second argument, and $\text{P}_{u}^{\pm}(u,v)$ the partial derivative with respect to the first argument. Numerical values of the Pearcey function can be computed rather efficiently by deforming the integration contour Connor and Curtis (1982) or by numerically solving a differential equation Connor and Farrelly (1981b).

Within the Pearcey approximation, one then obtains the following approximation to the integral (S24) Connor and Farrelly (1981a); Reijnders and Katsnelson (2017); S. Yu. Dobrokhotov et al. (2014):

[TABLE]

where the coefficients are defined by Eq. (S27). Taking the norm of this approximation, we indeed see that the position of its maximum does not depend on the function $f$, but only on the expansion coefficients of the action, i.e. on purely classical quantities. We see from Eq. (S29) that the Pearcey function is symmetric, and a numerical study shows that its maximum lies on the $x$-axis. One can show explicitly that Connor and Farrelly (1981b)

[TABLE]

where $J_{\alpha}(x)$ is the Bessel function. The maximum of $|P^{+}(x,0)|$ can be determined numerically and lies at $x_{0,+}=-2.19863$. By the symmetry relations (S29), we therefore find that the maximum of $|P^{\pm}(x,0)|$ lies at $x_{0,\pm}=\mp 2.19863$. Therefore, when $a_{4}$ is sufficiently large, we can find a good approximation for the position of the maximum by solving the set of equations

[TABLE]

We now return to the case $\theta=0$, for which we explicitly computed the position of the caustic in the previous section. The coefficient $a_{4,\alpha}$ is given by

[TABLE]

where all quantities are to be evaluated at $p_{y}=0$. We can obtain an expression for the fourth derivative of $p_{x}$ in the same way as we obtained an expression for the second derivative in Eqs. (S15)–(S17), using the fact that the odd order derivatives of $E_{\alpha}^{\pm}$ with respect to $p_{y}$ vanish. We find that

[TABLE]

Inserting this into Eq. (S34), we obtain a complex expression for $a_{4,\alpha}$. Performing a Taylor expansion in $\mu$, we find

[TABLE]

At $U_{0}=2E$, the term that is zeroth order in $\mu$ vanishes. We then have $a_{4,K}>0$, which means that the maximum lies left of the cusp point for the $K$-valleys, and $a_{4,K^{\prime}}<0$, i.e. the maximum lies right of the cusp point for the $K^{\prime}$-valley. In both cases, the maximum lies closer to the point $x=L_{1}$ than the cusp point.

The coefficients $\mathbf{b}_{n,\alpha}$ in Eq. (S27) can be computed as

[TABLE]

Specializing to $\theta=0$, we obtain

[TABLE]

where we have used that $\partial p_{x}/\partial p_{y}=0$. We can make an expansion of $b_{2,\alpha,x}$ and find

[TABLE]

Combining Eqs. (S33) and (S38), we observe that the maximum is on the $x$-axis, and that its position can be found by solving

[TABLE]

Hence, the position of the maximum equals

[TABLE]

We emphasize that we cannot use this approximation for low energies (small $\mu$) at $U_{0}=2E$, since $a_{4,\alpha}$ is small in that case, which implies that we are considering a section of a higher order caustic. In fact, we are too close to the Dirac case, where an ideal focus is observed at $U_{0}=2E$.

Finally, let us briefly return to units with dimensions and investigate the dependence of this expression on $L_{1}$ at constant $E$ and $U_{0}$. When restoring dimensions, the dimensionless variables $\tilde{x}_{\text{max},\alpha}$ and $\tilde{x}_{\text{cusp},\alpha}$, which are used in Eq. (S41), are replaced by $x_{\text{max},\alpha}/L_{1}$ and $x_{\text{cusp},\alpha}/L_{1}$. Furthermore, we note that $x_{\text{cusp},\alpha}$ scales linearly with length. The coefficient $a_{4,\alpha}$ is therefore independent of the length scale of the system, as is $b_{2,\alpha,x}$. Hence, the only part of the second term that scales with length is the semiclassical parameter $h$, which is proportional to $1/L_{1}$. Therefore, we conclude that the dependence of $x_{\text{max},\alpha}$ on $L_{1}$ is given by

[TABLE]

where $c_{1}$ and $c_{2}$ are constants that do not depend on $L_{1}$. This behavior is indeed observed in Fig. 2(e).

S2.4 Higher order singularities

For $\theta=0$, let us also consider the points where $a_{4,\alpha}$ vanishes and at which we pass through a higher order singularity. In the simplest approximation, we can analyze these points by looking at the first order Taylor expansion in $\mu$, Eq. (S36). For the $K$-valley, the term that is linear in $\mu$ is always positive, whereas it is always negative for the $K^{\prime}$-valley. The first term does not depend on the valley index, but is negative for $U_{0}>2E$, and positive for $U_{0}<2E$. Therefore, the coefficient $a_{4,\alpha}$ vanishes at $U_{0}>2E$ for the $K$-valley and at $U_{0}<2E$ for the $K^{\prime}$-valley. Neglecting terms of $\mathcal{O}(\mu^{2})$ in the expansion of $a_{4,\alpha}$ and solving for the potential $U_{0,\text{hs},\alpha}$ where it vanishes, we obtain

[TABLE]

When we consider the Dirac Hamiltonian, i.e. $\mu=0$, this reduces to $U_{0,\text{hs},\alpha}=2E$, as we already extensively discussed. Let us also compute the sixth order derivative of the action at the cusp, to see whether it is nonzero at the points where $a_{4,\alpha}$ vanishes. We have

[TABLE]

where the derivatives are to be evaluated at zero $p_{y}$. Using the same procedure as before to compute the sixth derivative of $p_{x}$ with respect to $p_{y}$ at $p_{y}=0$, we obtain

[TABLE]

Inserting the momenta (S18) and (S19) into this result, we obtain an expression for $a_{6,\alpha}$. We can then specialize to the point $U_{0,\text{hs},\alpha}$, where $a_{4,\alpha}$ vanishes. Inserting the result (S43) and expanding the result up to first order in $\mu$, we obtain

[TABLE]

We see that this coefficient is indeed nonzero when $\mu$ is nonzero, indicating that the butterfly singularity is the highest type of singularity that can be expected in the system when trigonal warping is included. This is confirmed by a numerical analysis, where we numerically solve for $U_{0,\text{hs},\alpha}$ without neglecting higher order terms in $\mu$ and insert the result into $a_{6,\alpha}$. For the $K^{\prime}$-valley, one then has $a_{6,\text{hs},\alpha}>0$ and for the $K$-valley, one has $a_{6,\text{hs},\alpha}<0$. The caustics that are observed in the set of classical trajectories, see Figs. 1(e) and S7, agree with this conclusion Poston and Stewart (1978).

Determining the position of the maximum for this higher order singularity requires a special function involving $\exp(iu^{6})$ in the integrand. We will not go into this here, since the generic singularity in our system is the cusp. Instead, we determine the position of the butterfly singularity. Entering $U_{0,\text{hs},\alpha}$ into Eq. (S14) and expanding the result up to first order in $\mu$, we obtain

[TABLE]

which reduces to the Dirac result $x_{\text{hs},\alpha}=L_{1}$ when we set $\mu=0$.

S2.5 Position of the maximum of the transmission peak when varying $U_{0}$

Let us now consider the case $\theta=0$ and have a look at the transmission as a function of potential strength. Fixing $L_{1}$ and $L_{2}$ and varying the potential $U_{0}$, one finds a value $U_{0,\text{max},\alpha}$ at which the transmission is maximal. We expect this maximum to be reached when the position $x_{\text{max},\alpha}$ of the maximum equals $L_{2}$. Therefore, we can find an approximation to $U_{0,\text{max},\alpha}$ by replacing $x_{\text{max},\alpha}$ by $L_{2}$ in Eq. (S40) and by subsequently solving for $U_{0}$. Since it is not easy to solve this equation analytically, we use numerical methods to find a solution.

Alternatively, we can consider the potential $U_{0,\text{cusp},\alpha}$ for which $x_{\text{cusp},\alpha}$ equals $L_{2}$. This will give us an estimate of the potential at which the maximum is reached. This estimate will generally be too low (high) for the valley $K$ ($K^{\prime}$), since the maximum in real space lies to the left (right) of the cusp point. We remark that when $L_{1}=L_{2}$ and we are considering low energies, this is the only method with which we can say something about the position of the maximum, as Eq. (S40) is not valid for this case due to the smallness of $a_{4,\alpha}$. Looking at Eq. (S14), we see that the only dependence on the length scale of the system is through $L_{2}/L_{1}$. Since these are dimensionless quantities, and their ratio is independent of the length scale of the system, our first conclusion is that $U_{0,\text{cusp},\alpha}$ does not depend on this length scale. Let us now compute $U_{0,\text{cusp},\alpha}$ for low energies, which correspond to small values of $\mu$. We can then use the expansion (S20) for $x_{\text{cusp},\alpha}$ and solve the equation $x_{\text{cusp},\alpha}=L_{2}$ for $U_{0}$. Expanding the final result for $U_{0,\text{cusp},\alpha}$ in $\mu$, we obtain

[TABLE]

For the Dirac case, where $\mu=0$, we obtain $U_{0,\text{cusp},\alpha}=E(1+L_{2}/L_{1})$, which is the solution of the equation $L_{1}(U_{0}-E)/E=L_{2}$ that one easily establishes for the Dirac equation. From Eq. (S48), we then immediately see that the maximum for the $K$-valley occurs at a lower potential than for the Dirac equation, whereas the maximum for the $K^{\prime}$-valley occurs at a higher potential than for the Dirac equation.

S2.6 Height of the transmission peak when varying $U_{0}$

We also investigate the maximal height of the transmission peak that is obtained by varying $U_{0}$, in order to see whether it differs for the valleys $K$ and $K^{\prime}$. Although the Pearcey approximation, Eq. (S31), generally predicts too high values Reijnders and Katsnelson (2017) for the absolute value of the wavefunction when we consider small $h$, we may be able to obtain some qualitative information about this difference.

Looking at Eq. (S31), we see that the peak height is proportional to the inverse of $|a_{4,\alpha}|$ and to the amplitude function $f(x_{\text{cusp},\alpha},y_{\text{cusp},\alpha},p_{y,0})$, which should be evaluated at $U_{0,\text{max},\alpha}$. However, we do not expect that the latter function will significantly contribute to a difference in peak height, considering its functional form for the Dirac equation, see Ref. Reijnders and Katsnelson (2017). For the Dirac equation, the function $f$ essentially consists of two parts. The first part comes from the electron source and depends on the angle of incidence and the sublattice polarization. The second part is the transmission coefficient through the barrier, which is to be evaluated at $p_{y,0}$. Because of the symmetry of the system, $p_{y,0}=0$ for both valleys. In the Dirac picture, this corresponds to normal incidence, i.e. the incidence angle $\phi=0$, and hence all phase factors that depend on the angle equal unity. Furthermore, at normal incidence we have Klein tunneling, i.e. the transmission is unity regardless of the height of the potential. Although trigonal warping certainly introduces corrections in the amplitude function, we do not expect that these will greatly influence either of the two factors at the cusp point. Therefore, we believe that the amplitude function $f$ at the cusp point depends only weakly on both the potential $U_{0}$ and the valley index $\alpha$. We thus believe that the factor that has the largest influence on a possible difference in peak height between the two valleys is the coefficient $|a_{4,\alpha}|$, on which we therefore focus our attention.

As with the position of the maximum, one generally needs to resort to numerical methods to obtain the coefficient $|a_{4,\alpha}|$ at $U_{0,\text{max},\alpha}$. However, to gain qualitative understanding of the effect, and to be able to obtain an analytical expression, we can also evaluate the coefficient $a_{4,\alpha}$ at $U_{0,\text{cusp}\alpha}$, given by Eq. (S48). Expanding the result up to first order in $\mu$, we obtain

[TABLE]

For $L_{1}=L_{2}$, we see that the first term vanishes. The absolute value of the second term is equal for both values of $\alpha$, implying that, for small values of $\mu$ and hence for low energies, the transmission peaks for both valleys will be equal in height when we consider $L_{1}=L_{2}$.

If we want to see whether this changes when we go to higher energies, we need to consider the second term in the expansion of $a_{4,\alpha}$. Unfortunately, when we expand $x_{\text{cusp},\alpha}$ to second order in $\mu$, the equation $x_{\text{cusp},\alpha}=L_{2}$ becomes third order in $U_{0}$, which makes the expression for $U_{0,\text{cusp},\alpha}$ much more complicated. This also makes the second order expansion of $a_{4,\alpha}$ in $\mu$ at $U_{0,\text{cusp},\alpha}$ very cumbersome. Nonetheless, we can say something about it by taking another look at our previous equations. The key observation is that in all our equations $\mu$ appears together with $\alpha$ in the combination $\alpha\mu$. This holds for the Hamiltonian (S10) and for the momenta (S18) and (S19) and therefore also for all derived quantities, such as $x_{\text{cusp},\alpha}$ and $a_{4,\alpha}$. This implies that we are actually dealing with an expansion in $\alpha\mu$ rather than in just $\mu$, which means that the second term in the expansion is proportional to $\alpha^{2}\mu^{2}$. Since $\alpha^{2}=1$, this term breaks the symmetry between the two valleys that is present in the linear term. Therefore, when we include the quadratic term in the expansion, $|a_{4,\alpha}|$ is different for the two values of $\alpha$.

Unfortunately, this analysis does not show whether the sign in front of the second term is positive or negative. A numerical study of the expansion shows that it is generally positive, i.e. that $|a_{4,K}|$ is larger than $|a_{4,K^{\prime}}|$. We find the same effect when we numerically evaluate $a_{4,\alpha}$ at $U_{0,\text{max},\alpha}$ for typical energies and length scales. This predicts that, for high energies, the height of the transmission peak for the $K$-valley is lower than the height for the $K^{\prime}$-valley, which is indeed seen in our numerical experiments.

We remark that the situation changes when we consider incident electrons with a uniform angular distribution, rather than a distribution that is uniform in momenta. For a uniform angular distribution, we need to integrate over the incidence angle $\phi$ rather than over the momentum $p_{y}$. Changing variables, we obtain an additional factor $\partial\phi/\partial p_{y}$ in the amplitude. It is clear that this does not change the classical trajectories, and that the position and width of the maximum will be the same for both distributions. However, the height of the transmission peak is affected by this factor. The angle $\phi$ is determined by the two components of the group velocity of the electrons, i.e.

[TABLE]

In the Pearcey approximation, we need to evaluate $\partial\phi/\partial p_{y}$ at $p_{y,0}=0$. After some calculus, using the fact that $\partial H_{\alpha}/\partial p_{y}$ vanishes at $p_{y}=0$, we find

[TABLE]

We therefore see that this factor is larger for the $K$-valley than for the $K^{\prime}$-valley, which could counter the effect of the coefficient $a_{4}$. Since the values given by the Pearcey approximation are generally too large for small values of $h$, these equations cannot tell us what the combined effect of these two factors will be. Our billiard simulations show that for a uniform angular distribution both peaks are approximately equal in height.

S2.7 Scaling of the width of the transmission peak

We have already stated that the maximal transmission will be obtained for the potential $U_{0,\text{max},\alpha}$ for which $x_{\text{max},\alpha}=L_{2}$. Changing the potential, we obtain a peak in the transmission, which has a full width at half maximum (FWHM) $\Delta U_{0,\alpha}$. Let us now consider the special case $L_{2}=L_{1}=L$ and see how this width $\Delta U_{0,\alpha}$ changes when we vary $L$. To this end, we first need to find the potential $U_{0,\text{left},\alpha}$, where the transmission equals a certain amount, for instance half, of its maximal value. Let us therefore denote by $x_{1,\pm}$ ($x_{2,\pm}$) the point on the $x$-axis to the left (right) of $x_{0,\pm}$ where the absolute value of the Pearcey function $P^{\pm}(x,0)$ equals $1/\sqrt{2}$ of its maximal value. By symmetry, these points are related by $x_{1,\pm}=-x_{2,\mp}$. Before, we obtained an estimate for $x_{\text{max},\alpha}$ by solving Eq. (S40). To obtain $U_{0,\text{left},\alpha}$, we replace $x_{0,\pm}$ in Eq. (S40) by $x_{2,\pm}$ and change $x_{\text{max},\alpha}$ to $L$. Then $U_{0,\text{left},\alpha}$ is given by the potential $U_{0}$ that solves this equation for a given energy. Similarly, replacing $x_{0,\pm}$ by $x_{1,\pm}$, we find $U_{0,\text{right},\alpha}$.

By the above procedure, we obtain two equations, namely

[TABLE]

where

[TABLE]

We remark that, since we are using dimensionless units, $L=1$ and $x_{\text{cusp},\alpha}$, given by Eq. (S14), does not depend on the length scale of the system. Since both $U_{0,\text{left},\alpha}$ and $U_{0,\text{right},\alpha}$ are close to $U_{0,\text{cusp},\alpha}$, we can perform a first order Taylor expansion of $x_{\text{cusp},\alpha}(E,U_{0})$ in $U_{0}$ around $U_{0,\text{cusp},\alpha}$. When we also expand $g(E,U_{0})$ to zeroth order around $U_{0,\text{cusp},\alpha}$, we obtain the set of equations

[TABLE]

We note that $U_{0,\text{cusp},\alpha}$ solves the equation $x_{\text{cusp},\alpha}(E,U_{0})=L$, whence the first term on the right-hand side cancels the term on the left-hand side in both equations. Subtracting the two equations above, and performing some elementary operations, we find that

[TABLE]

Since the predictions given by the Pearcey approximation are usually too large, we do not expect the estimate (S55) to be very precise. However, it does give us information on the scaling behavior of $\Delta U_{0,\alpha}$ with $L$. Since $U_{0,\text{cusp},\alpha}$ does not depend on the length scale of the system, as we established above, the only dependence on length is in the semiclassical parameter $h$. Therefore, we conclude that, as we vary the potential $U_{0}$, the width of the transmission peak scales as $L^{-1/2}$.

One may object that the scaling behavior expressed in Eq. (S55) was obtained using a first order Taylor expansion of $x_{\text{cusp},\alpha}$ and a zeroth order expansion of $g(E,U_{0})$. When we also include the first-order term in the Taylor expansion of $g(E,U_{0})$ around $U_{0,\text{cusp},\alpha}$, the expression for $\Delta U_{0,\alpha}$ becomes more complicated. However, when we expand the result, we still find that the leading-order term is proportional to $L^{-1/2}$, although we also obtain corrections of order $L^{-1}$ and higher. This statement remains true when we add the second-order terms in both Taylor expansions. Therefore, the leading-order term of the width of the transmission peak scales as $L^{-1/2}$.

Finally, let us briefly say something about the difference in peak width between the $K$ and $K^{\prime}$ valleys. It is hard to perform an analytical study of the Eqs. (S52) to find an analytical expression for the width of the peak. However, we can study their solutions numerically for both valleys. Because of the limitations of the Pearcey approximation, the width that is predicted by such a study is probably not very accurate. However, a numerical analysis does predict that the peak width will generally be larger for the $K^{\prime}$-valley than for the $K$-valley, in agreement with our numerical experiments.

S3 Details of the numerical simulations and data processing

This section contains additional information on the numerical implementation and the way the data were processed.

S3.1 Setup and data processing for the simulations using the Kwant code

For the numerical simulations using the Kwant code Groth et al. (2014), we create a sample using the setup shown in Fig. 1(a). It consists of an injector lead of width $W_{i}$ through which the particles enter and a collector lead of width $W_{c}$ through which they exit. Particles can also exit through the two drain leads on the side, which diminish scattering by the sides of the sample and thereby prevent infinite internal reflections. The injector lead is modelled as a semi-infinite ribbon, with translation vector $m\mathbf{v}_{1}+n\mathbf{v}_{2}$, where $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ are the lattice vectors defined in Eq. (S3). In our standard setup, $m=1$ and $n=0$, which means that the leads have zigzag edges along the translationally invariant direction. When $n\neq 0$, the translation vector makes an angle $\theta$ with the vector $\mathbf{v}_{1}$. This angle corresponds to the angle $\theta$ of the passive rotation that we considered in the detailed derivations. When $\theta=0$, which means that we are in the standard setup, we can easily separate the modes in the $K$ and the $K^{\prime}$-valley by their momentum Brey and Fertig (2006). When we consider rotated samples, we only consider those angles $\theta$ for which we can clearly separate the modes in different valleys. For $\theta=0$, we have checked that the lead wavefunctions obtained from the Kwant code are in agreement with the solutions of the continuum Dirac model Brey and Fertig (2006).

For the computations of the wavefunction in real space, $W_{c}$ equals $W$, the width of the sample. We choose the length $L_{2}$ from the junction interface to the collector lead sufficiently long, and have checked that it does not significantly influence the outcomes. The sample width $W$ was typically chosen as $2L_{1}$, with $L_{1}$ the distance from the injector lead to the junction interface. We consider this width realistic and at the same time sufficiently large to not substantially influence the outcome. From a computation of the $S$-matrix, we conclude that intervalley scattering is present, but is not the dominant scattering process. The amount of intervalley scattering, measured as a fraction of the total transmission into the collector lead, is larger for lead modes with higher transversal momentum, most likely due to interaction with the drain leads at the edges. It ranges from approximately 1% for the lowest lead mode, which is almost completely transmitted into the collector lead, up to 60% for higher lead modes, which are largely transmitted into the drain leads and therefore have a much smaller effect on the outcome. We come back to the size of intervalley scattering and provide a more quantitative picture later on, when we discuss scattering as a function of potential strength. When drain leads are absent, the position of the focus is somewhat changed, but valley splitting is still observed. We ascribe the change of position to additional internal reflections in the electron region, which make their way into the hole region and cause additional interference.

In order to be able to compare our computed wavefunction with semiclassical results, we do not plot the square of the absolute value of the wavefunction on each site, but instead average it over the two sublattices to obtain the envelope function. We do this, for each mode separately, by averaging the wavefunction over a site and its three neighboring sites. Numerically, we replace $|\Psi_{i,\alpha,n}|^{2}$, with $\Psi_{i,\alpha,n}$ the wavefunction on site $i$ resulting from lead mode $n$ in valley $\alpha$, by $|\Psi_{i,\text{av},\alpha,n}|^{2}=\tfrac{1}{2}|\Psi_{i,\alpha,n}|^{2}+\tfrac{1}{6}\sum_{\langle ij\rangle}|\Psi_{j,\alpha,n}|^{2}$.

We obtain the total wavefunction in valley $\alpha$ by adding the averaged wavefunctions of all of the modes in that valley, i.e. $|\Psi_{\text{av},\alpha}|^{2}=\sum_{n}|\Psi_{\text{av},\alpha,n}|^{2}$, where we have suppressed the site index $i$. Subsequently, we determine the position of the focus on the line $y=0$. Since the maximum of the wavefunction is sensitive to the details of the sample and shows rather strong fluctuations, we do not consider it a suitable measure for the position of the focus. Instead, we consider a collection of data sets, and fit a Gaussian to each of them. To obtain this collection, we first compute the minimum and maximum of $|\Psi_{\text{av},\alpha}|^{2}$ and denote it by $I_{\text{min}}$ and $I_{\text{max}}$, respectively. Subsequently, we draw a line at $I_{\text{min}}+x(I_{\text{max}}-I_{\text{min}})$, where $x$ varies, and extract all points between its left-most and right-most intersections with $|\Psi_{\text{av},\alpha}|^{2}$. We form a collection of data sets by letting $x$ vary in steps of 0.01, typically between 0.3 and 0.7. After the fitting, we extract the position of the maximum of each of the Gaussians. By computing the average and standard deviation, we obtain a reliable value and an error for the position of the focus. Hence, the error bar that is shown in the figures does not reflect the width of the focus, but rather how well its position is defined. We refer to Fig. S2 for two examples of the fitting procedure.

For the rotated samples, we also extract the total wavefunction on the $x$-axis, and determine the position of the peak for this data. This is probably not the position of the true focus, since we know from a study of the trajectories that the cusp point does not lie at $y=0$, see Fig. 1. Another indication for this is given by the tilting of the focus in the plots of the wavefunction for the rotated samples, see Fig. S1. However, the differences are seen to be rather small, and we obtain a good indication of the position of the focus by considering the data on the $x$-axis.

In most simulations, we consider an atomically sharp potential barrier. We have also done a few simulations for a smooth potential, see Fig. S2, which was implemented as

[TABLE]

We consider $L_{\text{NP}}$ to be a good measure for the length scale of the potential. Alternatively, one can compute the tangent to the potential $U(x)$ at $x=0$ and take the length scale of the potential to be the difference between the points $x_{-}$ and $x_{+}$, which are defined as the points where the tangent crosses zero and $U_{0}$, respectively. This leads to a width of $x_{+}-x_{-}=2L_{\text{NP}}/3$.

We obtain data for a varying potential by considering a collection of samples with different potentials. We typically divide the energy interval into 200 steps, and numerically compute the $S$-matrix for each value of the potential. In this case $W_{c}=W_{i}$, and $W=3L_{1}+W_{i}$ for $L_{1}=L_{2}=100$ nm and $W=2L_{1}+W_{i}$ when $L_{1}$ is larger. We obtain the tunneling coefficients for intravalley and intervalley scattering by adding the squared norms of the appropriate $S$-matrix elements and by subsequently dividing by the total number of modes in both valleys. We observe that intervalley scattering only mildly affects the results, see Fig. S3.

S3.2 Details for the billiard model and probability distributions

In the billiard simulations, electrons are modelled as billiard balls Beenakker and van Houten (1989); Milovanović et al. (2013, 2015). We consider $N$ electrons, with initial positions between $-W_{i}/2$ and $W_{i}/2$ that are drawn from a uniform distribution. Their initial transversal momenta are either drawn from a uniform distribution, or from a distribution that is uniform in the emission angle (see below). Subsequently, we compute their classical trajectories, and check whether the particles reach the collector lead, which has a width $W_{c}=W_{i}$. Furthermore, for each electron, we draw a random number between zero and one for each junction interface on its trajectory. When this number is smaller than the tunneling probability, the electron is transmitted by the interface; otherwise, it is reflected. We count the number $n_{T}$ of electrons that reach the collector lead, and define the tunneling probability as $T=n_{T}/N$. We observe a standard $1/\sqrt{N}$ convergence of $T$, roughly independent of the dispersion relation and the distribution of the transversal momenta. We use $N=250000$ in our simulations, for which the standard deviation of $T$ is approximately equal to 0.001.

We first consider a uniform distribution of the initial transversal momenta. For a nonrelativistic (Schrödinger) particle in a box, the momenta are equally spaced, which gives rise to a uniform momentum distribution in the continuum limit. This model is therefore valid when $W_{i}\gg\lambda_{\text{el}}$, with $\lambda_{\text{el}}$ the (de Broglie) wavelength of the electrons, see also Ref. Milovanović et al. (2013). Although the momenta in a graphene ribbon are not precisely equally spaced Brey and Fertig (2006), it was shown that a uniform momentum distribution gives correct predictions for graphene Hall barr experiments Milovanović et al. (2013, 2015). In the Dirac approximation, a uniform momentum distribution $f_{p_{y}}(p_{y})$ between $-p_{y,\text{max}}/2$ and $p_{y,\text{max}}/2$ corresponds to a distribution $f_{\phi}(\phi)=\cos\phi/(2\sin\phi_{\text{max}})$ of emission angles, where $p_{y,\text{max}}=E\sin\phi_{\text{max}}$. This can be easily verified using the formula

[TABLE]

and the relation $p_{y}=E\sin\phi$. We have verified that, within the Dirac approximation, both samplings indeed give the same result for the transmission.

We also consider a distribution that is uniform in the emission angle. This distribution, with a maximal emission angle of 45 degrees, was used in Ref. Lee et al. (2015), where the authors used billiard simulations to model their experiments. Following their approach, we also set the maximal emission angle to 45 degrees. For the Dirac equation, the emission angles can be converted into momenta using the relation given above. When we include trigonal warping, this conversion is more complicated, since the emission angle is determined by the group velocity. In this case, we obtain the transversal momenta by numerically finding the root of Eq. (S50).

The classical trajectories of the electrons are determined from Eq. (S13). The longitudinal momenta $p_{x}$ are computed by solving the equation $H_{\alpha}^{\pm}=E$ for $p_{x}$, given $p_{y}$. In the Dirac case this can be done analytically, but when we include trigonal warping we need to use numerical methods. We also perform a calculation using the full nearest-neighbor tight-binding Hamiltonian, given by Eq. (S1), and find that the tunneling probabilities only differ for high energies, see Fig. S3. Furthermore, we study the influence of a finite-size sample by allowing for a cutoff in the transversal direction. Taking out the electrons for which $|y(x=0)|>L_{1}+W_{i}/2$ in the simulations with uniformly distributed transversal momenta, we effectively limit the emission angle to 45 degrees, but find no significant influence on the results. We ascribe this to Klein tunneling, which collimates the electron beam, whence the largest contribution to the transmission comes from electrons that are emitted with small angles. We therefore plot the results without cutoff. Likewise, choosing maximal emission angles larger than 45 degrees in the model with uniformly distributed emission angles has a minimal effect on the position of the transmission maximum, although it does influence its value. Taking out electrons for which $|y(x=0)|>L_{1}+W_{i}/2$ at a maximal emission angle of 45 degrees also has very little influence.

The tunneling coefficient through an n-p junction is well-known for the Dirac equation Katsnelson et al. (2006); Reijnders and Katsnelson (2017). To obtain the tunneling coefficient when we include trigonal warping, we compute the eigenvectors of the matrix Hamiltonian, see Eqs. (S5) and (S6), and normalize them to unit current. Subsequently, we match them at the barrier interface. We observe that replacing the trigonal warping tunneling coefficient by the Dirac tunneling coefficient does not essentially change the observed transmission, see Fig. S3. Since the latter can be computed much more efficiently than the former, we use the Dirac tunneling coefficient in our simulations. This means that trigonal warping is only incorporated into the trajectories.

S3.3 Fitting the data for a varying potential

The data that we obtain when we vary the potential are rather smooth, both for the billiard simulations and for the simulations using Kwant. To further smoothen the data, we apply a Gaussian filter with radius $r$. For the Kwant data, this radius varies in integer steps between four and fourteen and for the billiard data it varies between two and seven. Subsequently, we compute the maximum for each of these smoothened data sets. We then average over the obtained maxima, and compute the standard deviation. For all cases, the standard deviation is smaller than the plot marker. To check our results, we also fit a Gaussian to a reasonable subset of the raw data and extract the position of its maximum. We find good agreement between the results of the two methods. The graphs show the results of the first procedure.

To determine the width of the peak as a function of $L_{1}=L_{2}=L$, we once again smoothen the data by applying a Gaussian filter with radius $r$. For the Kwant data, this radius varies between four and fourteen. For the billiard data, it varies between two and four. We then compute the potentials $U_{0,\text{left}}$ and $U_{0,\text{right}}$, left and right of the transmission maximum, respectively, for which the transmission equals 80 percent of its maximal value. Averaging over the values obtained for different Gaussian filters and subtracting $U_{0,\text{left}}$ from $U_{0,\text{right}}$, we obtain our final result.

S3.4 Data processing for simulations with a mass in the lead

In our simulations using the Kwant code, we include a mass term $m_{\text{lead}}$ in the lead by adding an on-site potential $m_{\text{lead}}$ for sites that belong to sublattice $A$, and $-m_{\text{lead}}$ for sites that belong to sublattice $B$. The obtained lead wavefunctions are in agreement with the results of the continuum model.

The computations of the wavefunction in real space are performed in the same way as before. To quantify the asymmetry in the averaged wavefunction, we consider slices of the sample with a given $x$-coordinate, and subsequently compute the fraction $\eta_{b,x,\alpha}$ of the total intensity on the slice that is located on the sites with $y\leq 0$, that is,

[TABLE]

Since the wavefunction at $x=-L_{1}$ is not symmetric in the line $y=0$, because the lead wavefunction of the incoming mode is not, we subtract the fraction $\eta_{b,-L_{1},\alpha}$ at $x=-L_{1}$ to see the change in the intensity distribution.

Furthermore, we consider the continuum evolution of the initial wavefunction outputted by the Kwant code. To this end, we extract the wavefunction from the Kwant simulation at the left-most point of the sample ($x=-L_{1}$). We split the sites into those belonging to sublattice $A$ and those belonging to sublattice $B$. Subsequently, we compute the evolution according to the continuum (Dirac) model, using Eq. (42) from Ref. Reijnders and Katsnelson (2017). We also compute the Fourier transform of the initial wavefunction using the equations from the same paper. Comparing these simulations with the simulations of the sample wavefunction using the Kwant code, we can establish which symmetry breaking effects can be explained by the continuum model, and which effects are due to trigonal warping.

We also consider a second type of simulations to probe the asymmetry induced by a mass in the lead. Here, we consider two collector leads instead of one, with the first one, called the top collector lead, between $y=0$ and $y=W/2$, and the second one, called the bottom collector lead, between $y=-W/2$ and $y=0$. We then compute the $S$-matrix for scattering into both collector leads and extract the transmission probabilities $T_{\text{bc},n,\alpha}$, which shows which part of incoming mode $n$ in valley $\alpha$ is scattered into the bottom collector lead, and $T_{\text{tc},n,\alpha}$, which indicates the fraction scattered into the top collector lead. Summing over all modes in a given valley, we obtain the total transmission probabilities $T_{\text{bc},\alpha}=\sum_{n}T_{\text{bc},n,\alpha}$ and $T_{\text{tc},\alpha}=\sum_{n}T_{\text{tc},n,\alpha}$. We then define the fraction that is transmitted into the lower lead as

[TABLE]

In order to obtain a measure for the asymmetry that is induced by the presence of the mass term, one should subtract the initial fraction of the intensity that is at negative $y$. Unfortunately, this is not a simple endeavor, since the total wavefunction in the lead is a sum of forward-moving modes and reflected modes and is not easily computed. Therefore, we confine ourselves to the forward-moving modes. For each mode, we extract the total intensity located on the sites with $y\leq 0$. Summing over all incoming modes and dividing by the total intensity, we obtain the fraction $\eta_{b,\text{in},\alpha}$ of the total intensity that is located on the sites with negative $y$. The procedure is the same as in Eq. (S58), the only difference being that we do not average over sublattices this time. We obtain the dependence of the asymmetry $\eta_{b,\text{out},\alpha}-\eta_{b,\text{in},\alpha}$ on $L_{2}$ by computing the $S$-matrix for a collection of samples with varying $L_{2}$ and fixed energy, mass and potential. We have checked that the measure $\eta_{b,\text{out},\alpha}-\eta_{b,\text{in},\alpha}$ is roughly independent of the width of the sample and is roughly equal to zero for zero mass (see Fig. 4), making it a suitable measure for the asymmetry. Furthermore, it evolves smoothly with changing width of the lead.

Finally, we compute the dependence of the transmission on the potential strength for samples with $W_{c}=W_{i}$ in the presence of a mass $m_{\text{lead}}$ in the lead. The computations using the Kwant code are done in the same way as before, the only difference being the mass in the lead. Since the billiard model does not include information on the initial wavefunction, the results for this model are not affected by the presence of a mass term in the lead.

S4 Additional simulation results

S4.1 Kwant simulations in real space

S4.2 Determination of the maximum

S4.3 Higher order effects in the numerical simulations

S4.4 Transmission as a function of potential strength

We note that in the simulations of the billiard model with uniformly distributed transversal momenta, reported in Fig. S4(b) and (c) and in Fig. S5, the transmission probabilities for the valleys $K$ and $K^{\prime}$ are equal at $U_{0}=2E$. To explain this, we first note that the sample that we are considering is symmetric in the junction interface ($x=0$), as $L_{1}=L_{2}=L$ and $W_{c}=W_{i}$. Second, we note that for $\theta=0$ and at $U_{0}=2E$ there is a special relation between electron and hole modes in the two valleys. We have $p_{x,e,K}(p_{y})=-p_{x,h,K^{\prime}}(p_{y})$ and $p_{x,h,K}(p_{y})=-p_{x,e,K^{\prime}}(p_{y})$. To convince oneself of these relations, one can visualize the energy contours of the dispersion relation (S8) and invoke electron-hole symmetry.

Alternatively, one can also show the above relations explicitly, though the derivation is slightly cumbersome. Consider an electron in valley $K$, with energy $E$ and transversal momentum $p_{y}$. Let $p_{x,e,K}(p_{y})$ be the associated longitudinal momentum, for which the relation $E=|\mathbf{p}|-\mu|\mathbf{p}|^{2}\cos(3\phi_{\mathbf{p}})$ holds, see Eq. (S10). We now note that $-\cos(3\phi_{\mathbf{p}})=\cos(3[\pi-\phi_{\mathbf{p}}])$. Multiplying both sides by minus one and adding $2E$, we obtain the relation $E=-|\mathbf{p}|-\mu|\mathbf{p}|^{2}\cos(3[\pi-\phi_{\mathbf{p}}])+2E$. Looking at Eq. (S10), we see that this is exactly the dispersion relation for a hole with momentum $-p_{x,e,K}(p_{y})$ in valley $K^{\prime}$ at $U_{0}=2E$. Therefore, we arrive at the relation $p_{x,e,K}(p_{y})=-p_{x,h,K^{\prime}}(p_{y})$. The second relation can be establish analogously.

Let us now consider an electron in the $K$-valley, which is emitted from the point $y_{i}$ in the injector lead with transversal momentum $p_{y,s}$ and longitudinal momentum $p_{x,e,K}(p_{y,s})$. This electron hits the junction interface at the point $y_{e}$, after which its longitudinal momentum becomes $p_{x,h,K}(p_{y,s})$. We first assume that it is transmitted to the point $y_{c}$ in the collector lead. Let us now also consider an electron in the $K^{\prime}$-valley, which is emitted from the point $y_{c}$ in the injector lead with transversal momentum $p_{y,s}$. Then, because of the above relations, its longitudinal momentum equals $p_{x,e,K^{\prime}}(p_{y,s})=-p_{x,h,K}(p_{y,s})$. Since $L_{1}=L_{2}$, this electron hits the junction interface at the same coordinate $y_{e}$, after which its longitudinal momentum becomes $p_{x,h,K^{\prime}}(p_{y,s})=-p_{x,e,K}(p_{y,s})$. Because $L_{1}=L_{2}$, this particle subsequently reaches the collector lead at the point $y_{i}$ and is therefore transmitted. This means that we have established a one-to-one relationship between trajectories in the $K$-valley that are transmitted and trajectories in the $K^{\prime}$-valley that are transmitted. When the electron in the $K$-valley that we considered does not reach the collector lead, the above procedure gives us an invalid trajectory in the $K^{\prime}$-valley, since this new trajectory does not start from a point in the lead.

We also note that, since we have used the Dirac tunneling probability in our simulations, the tunneling probability at the barrier interface does not depend on the valley. Then, because of the one-to-one relationship between transmitted trajectories in the two valleys, the amount of electrons transmitted in valley $K$ equals the amount of electrons transmitted in valley $K^{\prime}$. Of course, this supposes that our simulations have converged, i.e. that we have sampled a sufficiently large number of initial positions and initial transversal momenta.

With the above argument, we have shown that, at $U_{0}=2E$, the transmission probabilities for the valleys $K$ and $K^{\prime}$ are equal within the framework of a billiard model with uniformly distributed transversal momenta. In principle, this argument should also hold for a billiard model with uniformly distributed emission angles. However, in our simulations, we have limited the maximal emission angle to $\pi/4$. This breaks the above argument, as it is not certain that the two related trajectories in the valleys $K$ and $K^{\prime}$ both have emission angles smaller than $\pi/4$. Indeed, in Fig. S4(a) and (d), we observe that for a billiard model with uniformly distributed emission angles the two transmission probabilities are not equal at $U_{0}=2E$.

S4.5 Numerical results for the higher order singularity

S4.6 Breaking of sublattice symmetry in the lead