# Three-dimensional chiral Veselago lensing

S. Tchoumakov<EMAIL_ADDRESS>Univ. Grenoble Alpes, CNRS,
Grenoble INP, Institut Néel, 38000 Grenoble, France  J. Cayssol  Univ.
Bordeaux, CNRS, LOMA, UMR 5798, F-33405 Talence, France  A. G. Grushin  Univ.
Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 38000 Grenoble, France

###### Abstract

The effect by which light focuses upon entering a medium with a negative
refractive index, known as Veselago lensing, may enable optical imaging below
the diffraction limit. Similarly, focusing electrons across a $pn$-junction
could realize a technologically promising electronic Veselago lens. However,
its scope remains limited by the lack of three-dimensional platforms and its
insensitivity to computational degrees of freedom, like spin or chirality.
Here we propose a single-material three-dimensional electronic Veselago lens
that selectively focuses electrons of a given chirality. Using the chiral
anomaly of topological semimetals it is possible to create a sharp
$pn$-junction for a single chirality, a chiral Veselago lens, and tune it with
a magnetic field to an ideal lensing condition. We estimate that chiral
Veselago lensing is observable in non-local transport and spectroscopy
experiments. In particular we show that the chiral Veselago lens leads to
giant non-local magnetoresistance.

Introduction. The similarities between the light-ray construction in optics
and the semi-classical trajectories of electrons [1, 2] host the potential for
new applications in electronics, such as electronic lenses, interferometers or
beam-splitters at the micrometer scale [3, 4, 5]. To reach this potential,
controlling the involved material interfaces is crucial. In optics, the
trajectory of light is deflected at the interface between two media, allowing
to focus, guide and disperse light controllably. In particular, between media
with opposite handedness, where the optical index changes sign, lensing can
occur even for a flat interface, a phenomenon coined Veselago lensing [6]. In
electronics, a $pn$ junction can act like an electronic Veselago lens by
focusing two-dimensional (2D) electrons with pseudo-relativistic, linear band
dispersion [3, 7, 8, 9, 10, 11, 12], an effect that has only been realized in
graphene [13, 14, 15, 16].

The prospect of creating useful devices based on electronic Veselago lenses
faces three fundamental problems. The first problem is that we lack a three-
dimensional (3D) Veselago lens. Theoretically, a $pn$ junction sharper than
the Fermi wavelength based on 3D relativistic metals, like Weyl or Dirac
metals [17], can Veselago-lens [7, 12], in analogy with graphene. While bulk
3D metals have the advantage of a larger carrier density, this property is
also a drawback as they are not easily gated. The second problem is creating a
clean interface. A $pn$ junction created via chemical doping or by interfacing
$p$-doped and $n$-doped samples will likely be disordered by impurities or
lattice mismatch. The interface could accumulate charge, affecting electronic
transport in undesirable ways [18]. The last problem is that any proposed or
realized electronic Veselago lens seems insensitive to computational degrees
of freedom, like spin or chirality. These three problems materialize in the
challenge to realize a single-material 3D Veselago lens that could select
specific electronic degrees of freedom.

Figure 1: Chiral Veselago lens based on the chiral anomaly. A uniform magnetic
field $\mathbf{B}$ is applied to a Weyl semimetal, while the electric field
$\mathbf{E}$ is only applied for $z<0$, with two metallic gates (left
orange/gold gates). This configuration generates a charge imbalance between
Weyl quasiparticles of opposite chiralities (red and blue cones) for $z<0$ due
to the chiral anomaly, while maintaining the total charge constant. This
creates an ideal $pn$ junction for a single chirality (blue cones), that
Veselago-lenses a measurable non-local current at $z>0$, tuned by changing
$\mathbf{B}$, and measured by narrow electrodes (in orange/gold at right).

In this work we take a significant step to solve these problems by proposing
how to Veselago-lens a single chirality of 3D Weyl quasiparticles
controllably. The chirality is a quantum mechanical degree of freedom of 3D
Weyl quasiparticles, labelled by $\pm$ depending on whether the spin is
aligned or anti-aligned with the electron’s momentum, which can be used for
computation [19]. By locally activating the chiral anomaly, an effect that
overpopulates one chirality with respect to the other in the presence of
colinear electric and magnetic fields, it is possible to create a sharp
$pn$-junction for a single chirality. The resulting chiral Veselago lensing
can be detected either by spectroscopic probes, like scanning-tunneling
microscopy (STM), or non-local transport.

Because of the chiral anomaly, the lensing can be tuned with electric and
magnetic fields to modulate the intensity of the image charge in STM or the
image current in non-local transport. We also show that ordinary,
quadratically dispersing electrons can Veselago-lens stronger than Weyl
electrons, albeit without the degree of control offered by the chiral anomaly
and suffering from a charging layer at the $pn$ interface. Therefore the
chiral anomaly, unique to Weyl semimetals, is the optimal tool to realize a
clean 3D chiral Veselago lens.

Figure 2: Veselago lensing of a single chirality in spectroscopy and non-local
transport. (a) Schematic set-up to measure Veselago lensing spectroscopically
(see [20] for a two-chirality set-up). (b) The upper plot shows the non-local
polarisability $\Pi(z,z^{\prime})$ for an $nn$ junction as a function of $z$
for $(x,y)=(0,0)$ for both non-relativistic (dashed) and Weyl band structures
(solid). The contour plot shows $\Pi(z,z^{\prime})$ in the $(x,z)$ plane for
$y=0$ for the Weyl case. (c) Same as (b) for a $pn$ junction. In (b-c)
$\Pi(z,z^{\prime})$ is given in units of $\mu N_{0}$, with
$N_{0}=k_{F}^{3}/(2\pi^{2}\mu_{L})$ the density of states on the left of the
junction. (d) Schematic set-up to measure Veselago lensing in non-local
transport (see [20] for a two-chirality set-up). (e) and (f) show the non-
local conductivity $\sigma_{xx}(z,z^{\prime})$ for $nn$ and $pn$ junctions,
respectively. (b,e) are calculated with $(m,\mu,\ell)_{L}=(1,0.5,50)$ and
$(v,\mu,\ell)_{R}=(2,4,50)$ for the non-relativistic case, and
$(v,\mu,\ell)_{L}=(1,1,50)$ and $(v,\mu,\ell)_{R}=(1,2,50)$ for the Weyl case.
(c,f) are calculated with $(m,\mu,\ell)_{L}=(1,0.5,50)$ and
$(v,\mu,\ell)_{R}=(-1,-0.5,50)$ for the non-relativistic case, and
$(v,\mu,\ell)_{L}=(1,1,50)$ and $(v,\mu,\ell)_{R}=(1,-1,50)$ for the Weyl
case. In (e-f) the conductivity is given in units of the total conductivity
$\sigma_{0}$ over the mean-free path $\ell$ (see [20]).

Model. Weyl semimetals host pairs of linearly dispersing bands separated in
momentum space, known as Weyl cones. Each cone can be described by the
Hamiltonian

$\displaystyle\hat{H}_{W}=(v_{F}{\bf k}-{\bf b})\cdot\hat{\bm{\sigma}}-\mu,$
(1)

where $\bm{\sigma}$ is a vector of Pauli matrices, $\mu$ is the chemical
potential, $v_{F}$ is the Fermi velocity, and ${\bf b}$ determines the
momentum space separation between Weyl cones. The eigenstates of $\hat{H}_{W}$
depend on the sign of $v_{F}$, which defines the chirality
$\mathrm{sgn}(v_{F})=\chi=\pm$. Weyl cones come in pairs of opposite $\chi$
and, in certain materials, they can be tilted [21] or anisotropic [17]. Here
we focus on materials where these effects are negligible, as discussed at the
end.

In presence of electric ($\mathbf{E}$) and magnetic (${\bf B}$) fields, the
transport of Weyl quasiparticles depends on their respective chirality, as
reflected by the continuity equations

$\displaystyle\frac{\partial n_{\chi}}{\partial t}+\bm{\nabla}\cdot{\bf
j}_{\chi}=\chi\frac{e^{2}}{2\pi^{2}\hbar^{2}}\mathbf{E}\cdot{\bf B},$ (2)

where $n_{\pm}$ and ${\bf j}_{\pm}$ are, respectively, the charge and current
distributions of carriers with $\chi=\pm$. The source term, on the right-hand
side, creates an imbalance $\Delta n$ between the chiral charges $n_{+}$ and
$n_{-}$, without altering the total charge $n=n_{+}+n_{-}$. This term, known
as the chiral anomaly, is responsible for anomalous transport responses, such
as the enhancement of the magneto-conductance of Weyl semimetals [22, 17].

The chiral anomaly competes with inter-valley scattering, leading to a charge
imbalance between Weyl nodes with opposite chiralities in the steady-state.
The inter-valley scattering time $\tau$ can range from a picosecond in TaAs to
a nanosecond in Na3Bi [23, 24]. The chiral charge imbalance $\Delta
n=n_{+}-n_{-}=\frac{\tau e^{2}}{2\pi^{2}\hbar^{2}}\mathbf{E}\cdot{\bf B}$
equals the carrier density $n$ for a magnetic field $B_{c}\equiv
2\mu^{3}/(3v^{3}e^{2}\tau E)$. For $n=10^{21}$ cm-3, $\tau=10^{-9}$ s [24, 25,
23] and $E\approx 10^{6}$ V.m-1, this crossover field is $B_{c}=1$ T. Beyond
$B_{c}$, the chiral anomaly reverses the type of carriers of one of the Weyl
nodes, e.g. from electrons to holes. Hence, if $B>B_{c}$ is applied
homogeneously and the $E$ is not homogeneous, e.g. stepwise as illustrated in
Fig. 1, the chiral anomaly creates an effective $pn$ junction for a single
Weyl chirality.

The main goal of this work is to explain how to realize and detect a single-
chirality Veselago lens, under the above conditions. To this end, we discuss
spectroscopy and non-local transport responses of a single Weyl cone to
determine the ideal conditions for single-chirality Veselago lensing.
Combining the contributions from the two chiralities we discuss how chiral
Veselago lensing can be detected using the simple set-up in Fig. 1. In the
Supplemental Material [20] we propose two related, albeit less practical,
devices to detect directly chiral Veselago lensing in spectroscopy and in non-
local transport.

We model the change in the carrier density with a step-wise chemical
potential, where
$\mu_{L}=\left(\mu^{3}\pm\frac{3}{2}v^{3}e^{2}\tau\mathbf{E}\cdot{\bf
B}\right)^{1/3}$ for $z<0$ [26] and $\mu_{R}=\mu$ for $z>0$ ($L$ and $R$
denote left and right of the interface, respectively).

Spectroscopic signatures of chiral Veselago lensing. To visualize and
understand the ideal conditions for Veselago lensing, we calculate the charge
distribution created by a local potential due to an impurity, or the tip of a
STM. For a local potential $V({\bf r})=\delta({\bf
r}_{\parallel})\delta(z-z_{0})\delta V$ there is a redistribution of charge
$\delta n(z)=\Pi(z,z_{0})\delta V$ that depends on the polarisability [27]

$\displaystyle\Pi(z,z^{\prime})=-\frac{1}{2\pi}\int d\omega~{}{\rm
Tr}\left(\hat{G}(z,z^{\prime})\hat{G}(z^{\prime},z)\right),$ (3)

written in terms of the Green’s function $\hat{G}(z,z^{\prime})$, where we
account for scattering with a mean free path $\ell\gg 1/k_{F}$ (see [20] for
details).

In Fig. 2(b,c) we plot the polarisability as a function of the $x$ and $z$
coordinates. The coordinates are given in units of $1/k_{F}$, and the
polarisability in units of $\mu_{L}N_{0}$, with
$N_{0}=k_{F}^{3}/(2\pi^{2}\mu_{L})$ being the density of states at $z<0$,
where the impurity is located. We compare the case of a single Weyl chirality
with a 3D electron gas, with a quadratic band structure [20]. In a $nn$ (or
$pp$) junction, the charge distribution shows the usual Friedel oscillations
on both sides of the junction, see Fig. 2(b). The period is different on both
sides due to the change in the Fermi wavevector when crossing the interface
[27].

The main difference between the $nn$ junction in Fig. 2(b) and the $pn$
junction in Fig. 2(c) is the appearance of an image charge as a consequence of
Veselago lensing. Veselago lensing occurs because the in-plane group velocity
changes sign across the interface, ${\bf v}_{L}=-{\bf v}_{R}$, due to the
conservation of the in-plane wavevector ${\bf k}_{\parallel}$. This condition
can be met for both a 3D electron gas and a Weyl semimetal $pn$ junction, but
the intensity of the image charge is larger for the former than the latter due
to the slower decay of Friedel oscillations in a normal electron gas [28, 29,
30, 31].

However, a drawback of using a 3D electron gas is that the electron density
accross the $pn$ junction is not constant, and one should expect a built-in
interface potential over a finite distance $d$. When $k_{F}d\gg 1$, we find
[20] that the amplitude of the Veselago image rapidly decays [14]. For typical
3D metals $d>10~{}\mu$m so $k_{F}d\gg 1$ [32] implying that the corresponding
image charge is unobservable in practice [14]. In contrast, in the chiral
Veselago lens in Fig. 1 the total charge remains constant, charges are only
transferred between cones of opposite chirality. The electric field extends
beyond the parallel plates on a distance of the order of the distance between
the two plates, so we can expect a sharp interface potential for thin film
samples.

So far we have assumed symmetric $pn$ junctions, i.e. those where $k_{F}$ is
equal on both sides of the interface. Deviations from this condition blur and
change the location of the image charge [20]. Reaching ideal lensing is
unrealistic with 3D electron gases because finding junctions with equal
$k_{F}$ but opposite carrier types is challenging in practice. Later we will
argue how the chiral anomaly of 3D Weyl semimetal aids to tune into this ideal
condition by varying the magnetic field.

Non-local transport signatures of chiral Veselago lensing. While the image
charge is convenient to understand how Veselago lensing can be enhanced (see
Supplemental Material [20] for a proposed device), surface states
contributions must be factored out [33, 34, 35, 36] to reveal Veselago
lensing. A bulk, non-local transport measurement, such as that depicted in
Fig. 2(d), is in this sense a simpler set-up. In a local electric field ${\bf
E}(z)={\bf E}_{0}\delta(z-z_{0})$ the electronic current
$j_{\mu}(z)=\sigma_{\mu\nu}(z,z_{0})E_{\nu}$ is obtained from the non-local
conductivity [37, 38]

$\displaystyle\sigma_{\mu\nu}(z,z^{\prime})=\int\dfrac{dS_{z}dS_{z^{\prime}}}{\pi\mathcal{A}}{\rm
Tr}\left(\hat{j}_{\mu}{\rm Im}\hat{G}(z,z^{\prime})\hat{j}_{\nu}{\rm
Im}\hat{G}(z,z^{\prime})\right),$ (4)

where $\hat{j}_{\mu}$ are the components of the current operator, with
$\mu=x,y,z$, and $S_{z}$ and $S_{z^{\prime}}$ are planes at $z$ and
$z^{\prime}$ with areas $\mathcal{A}$. The non-local conductivity is a complex
quantity, that accounts for the dephasing between the two probes. In Eq. (4)
we only show its real part, which can be measured by averaging the
conductivities obtained after permuting the positions $z$ and $z^{\prime}$ of
the leads (see Supplemental Material [20]).

In Figs. 2(e,f), we show the non-local conductivity for a $nn$ (or $pp$) and a
$np$ junction. For a $nn$ junction, the current is positive and decays
exponentially away from the input electric field, as shown in Fig. 2(e). In
contrast, for the $pn$ junction, the current changes sign when crossing the
interface and its magnitude peaks close to the image charge, signaling the
presence of a Veselago lens [39]. Note that the negative current does not
violate energy conservation, since the total current is positive, and is a
consequence of the change of type of the main carriers.

As with the image Veselago charge, the image Veselago current is larger for a
3D electron gas than for a Weyl semimetal, where it oscillates close to zero,
but negative on average. These oscillations are absent in the 3D electron gas
which is a single-band model, and thus we attribute them to interband
excitations that lead to hole-like and electron-like regions. Lastly, we find
that ideal lensing is achieved when both sides of the $pn$ junction are tuned
to the have the same $k_{F}$, as for the image Veselago charge.

Figure 3: Transport signatures of the chiral Veselago lens depicted in Fig. 1.
(a) The total non-local conductivity $\sigma_{xx}(z)$ (dark curve) for
$B=2B_{c}$. The contribution of each chirality is plotted as colored curves.
$\sigma_{xx}(z)$ is uniform for $z<0$ where the electric field is constant and
drops for $z>0$, with a negative conductivity for Weyl quasiparticles with
chirality $\chi=-$. (b) Giant magnetoresistance due to the chiral anomaly for
a probe located at $z=10/k_{F}$; see dashed line in (a). The magnetic field is
in units of $B_{c}=2\mu^{3}/(3v^{3}e^{2}\tau E)\approx 1T$ (see text). Chiral
Veselago lensing starts at $B=B_{c}$, which marks the change from electrons to
holes for chirality $\chi=-$ (inset red cone). For $B>B_{c}$ carriers with
$\chi=-$ contribute with a negative current that is maximal at $B=2B_{c}$,
when the Fermi wavevector is opposite of both sides of the junction. The
strong increase in magnetoresistance, marked by the blue region, is thus a
fingerprint of chiral Veselago lensing.

3D Chiral Veselago lens. We have previously considered the transport of each
chirality individually, while in the set-up of Fig. 1, the two chiralities are
present. In this set-up the electric field $E$ is finite for $z<0$, and zero
for $z>0$. In Fig. 3(a) we show the resulting total non-local conductivity and
the contribution of each chirality, calculated using the expression

$\displaystyle\sigma_{xx}(z)=\int_{-\infty}^{0}dz^{\prime}~{}\sigma_{xx}(z,z^{\prime}).$
(5)

For $z>0$, $\sigma_{xx}(z)$ decays slowly and the contribution of each
chirality is about $10\%$ of the bulk conductivity at $zk_{F}=10$. Since the
conductivity is negative for one chirality and positive for the other, the two
contributions compensate and lead to a giant magnetoresistance with respect to
the configuration without chiral anomaly. This follows from our previous
discussion: because of the chiral anomaly, carriers of one chirality (here
$\chi=-$) experience a $pn$ junction while carriers of the opposite chirality
(here $\chi=+$) experience a $nn$ junction. Since $E$ extends for all $z<0$,
the $pn$ junction leads to an extended, rather than localized, Veselago image.
Nonetheless, the ideal condition for an optimal image is still
$k_{F,L}=k_{F,R}$. Accordingly, we find that the negative current is maximal
when this condition is met.

The total $\sigma_{xx}(z)$ of the junction is a sum of the contributions of
each Weyl cone. As seen in Fig. 3(a), the total $\sigma_{xx}(z)$ is positive
throughout the junction (dark curve), since the total charge is conserved,
hiding the signature of Veselago lensing. However, Veselago lensing becomes
apparent when studying the non-local resistivity as a function of $B$, as
shown in Fig. 3(b). At $B<B_{c}$ there is no chiral Veselago lensing and the
magnetoresistance increases moderately. This is due to the increasing mismatch
between the carrier densities on each side of the junction. Once $B>B_{c}$, an
effective $pn$ junction is formed for one chirality (here $\chi=-$) that leads
to a negative contribution to current and thus to a giant increase in the
magnetoresistivity (shaded region in Fig. 3(b)). The amplitude of the negative
current is maximal at $B=2B_{c}$ where $k_{F,L}=k_{F,R}$ for $\chi=-$, and
Veselago lensing is optimal. For $B>2B_{c}$ the amplitude of the negative
current decreases, and the magnetoresistance increases moderately. The strong
increase in magnetoresistance within $[B_{c},2B_{c}]$, highlighted in Fig.
3(b), is thus a fingerprint of chiral Veselago lensing.

Discussion. We have determined how the chiral anomaly of Weyl semimetals can
be used to create a Veselago lens that distinguishes electrons by their
chirality. We have uncovered that the chiral Veselago lens based on Weyl
semimetals is advantageous for several reasons. First, a chiral Veselago lens
eliminates the charge build up at the interface, since it requires a single
material, and the total charge remains constant along the sample. Second, the
chiral anomaly allows to control the $k_{F}$ of each chirality using a
magnetic field. This allows to detect chiral Veselago lensing by measuring
non-local conductivity and tune it to the ideal condition for lensing. This
property goes beyond current graphene-based proposals, since the chiral
anomaly is unique to 3D Weyl semimetals.

We have neglected several effects in our computations. First, we neglected the
orbital localization in a magnetic field, since this occurs over the magnetic
length scale $\ell_{B}\approx 25.6$ nm$/\sqrt{B{\rm[T]}}\gg z,\lambda_{F}$.
Second, tilted Weyl cones [21] could prevent a perfect matching of $k_{F}$.
However, the chiral anomaly is often observed in systems where the Weyl cones
are generated by magnetic field from a lightly doped semimetal, e.g. GdPtBi,
ZrTe5, or KZnBi [40, 41, 17, 42] where the effect of tilting is negligible.
Lastly, Fermi arcs will be negligible in non-local transport, as it is a bulk
probe.

The estimates we provide indicate that chiral Veselago lensing is observable
under moderate electric and magnetic fields. Non-local transport experiments
have probed the chiral anomaly, albeit without $pn$-junctions and lensing [23,
43], suggesting that our proposal can be realized using current technology.
Our work opens the possibility to realize a device that can control a
computational degree of freedom using Veselago lensing in 3D materials,
enabled by the chiral anomaly.

Acknowledgements.- We thank B. Gotsman, H. Schmid, A. Molinari and J. Gooth
for discussions about experimental details. A. G. G and S. T acknowledge
financial support from the European Union Horizon 2020 research and innovation
program under grant agreement No 829044 (SCHINES). A. G. G. is also supported
by the ANR under the grant ANR-18-CE30-0001-01 (TOPODRIVE).

## References

  * Yan  _et al._ [2020] C. Yan, M. Pepper, P. See, I. Farrer, D. Ritchie, and J. Griffiths, Scientific Reports 10 (2020), 10.1038/s41598-020-59453-x.
  * Spector  _et al._ [1990] J. Spector, H. L. Stormer, K. W. Baldwin, L. N. Pfeiffer, and K. W. West, Applied Physics Letters 56, 1290 (1990).
  * Cheianov  _et al._ [2007] V. V. Cheianov, V. Fal’ko, and B. L. Altshuler, Science 315, 1252 (2007).
  * Westström and Ojanen [2017] A. Westström and T. Ojanen, Phys. Rev. X 7, 041026 (2017).
  * Karalic  _et al._ [2020] M. Karalic, A. Štrkalj, M. Masseroni, W. Chen, C. Mittag, T. Tschirky, W. Wegscheider, T. Ihn, K. Ensslin, and O. Zilberberg, Phys. Rev. X 10, 031007 (2020).
  * Veselago [1968] V. Veselago, Soviet Physics Uspekhi 10, 509 (1968).
  * Hills  _et al._ [2017] R. D. Y. Hills, A. Kusmartseva, and F. V. Kusmartsev, Phys. Rev. B 95, 214103 (2017).
  * Libisch  _et al._ [2017] F. Libisch, T. Hisch, R. Glattauer, L. A. Chizhova, and J. Burgdörfer, Journal of Physics: Condensed Matter 29, 114002 (2017).
  * Tian and Wang [2017] H. Tian and J. Wang, Journal of Physics: Condensed Matter 29, 385401 (2017).
  * Zhang  _et al._ [2017a] S.-H. Zhang, J.-J. Zhu, W. Yang, and K. Chang, 2D Materials 4, 035005 (2017a).
  * Zhang and Yang [2019] S.-H. Zhang and W. Yang, New Journal of Physics 21, 103052 (2019).
  * Yang  _et al._ [2021] Y. Yang, Y. Bi, L. Peng, B. Yang, S. Ma, H.-C. Chan, Y. Xiang, and S. Zhang, Optica 8, 249 (2021).
  * Lee  _et al._ [2015] G.-H. Lee, G.-H. Park, and H.-J. Lee, Nature Physics 11, 925 (2015).
  * Chen  _et al._ [2016] S. Chen, Z. Han, M. M. Elahi, K. M. M. Habib, L. Wang, B. Wen, Y. Gao, T. Taniguchi, K. Watanabe, J. Hone, A. W. Ghosh, and C. R. Dean, Science 353, 1522 (2016).
  * Zhang*  _et al._ [2021] X. Zhang*, W. Ren*, E. Bell, Z. Zhu, K. Watanabe, T. Taniguchi, E. Kaxiras, M. Luskin, and K. Wang, “Gate-tunable veselago interference in a bipolar graphene microcavity,” (2021), arXiv:2106.09651 [cond-mat.mes-hall] .
  * Brun  _et al._ [2019] B. Brun, N. Moreau, S. Somanchi, V.-H. Nguyen, K. Watanabe, T. Taniguchi, J.-C. Charlier, C. Stampfer, and B. Hackens, Phys. Rev. B 100, 041401 (2019).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Reviews of Modern Physics 90, 15001 (2018), arXiv:1705.01111 .
  * Zhang and Yates [2012] Z. Zhang and J. T. Yates, Chemical Reviews 112, 5520 (2012), pMID: 22783915.
  * Kharzeev and Li [2019] D. E. Kharzeev and Q. Li, “The chiral qubit: quantum computing with chiral anomaly,” (2019), arXiv:1903.07133 [quant-ph] .
  * [20] The Supplemental Material includes the detailed calculation of the Green’s functions of a junction of a 3D electron gas and of a Weyl semimetal. It also contains detailed calculations of the bulk non-local conductivity of a Weyl semimetal and figures to illustrate how Veselago lensing is affected by a smooth junction and to changes in the carrier densities. Lastly, it contains two additional device proposals to detect the chiral Veselago lensing, one using an STM, and another using non-local transport.
  * Soluyanov  _et al._ [2015] A. A. Soluyanov, D. Gresch, Z. Wang, Q. Wu, M. Troyer, X. Dai, and B. A. Bernevig, Nature 527, 495 (2015).
  * Landsteiner [2016] K. Landsteiner, Acta. Phys. Pol. B 47 (2016).
  * Parameswaran  _et al._ [2014] S. A. Parameswaran, T. Grover, D. A. Abanin, D. A. Pesin, and A. Vishwanath, Phys. Rev. X 4, 031035 (2014).
  * Behrends  _et al._ [2016a] J. Behrends, A. G. Grushin, T. Ojanen, and J. H. Bardarson, Phys. Rev. B 93, 075114 (2016a).
  * Ashby and Carbotte [2014] P. E. C. Ashby and J. P. Carbotte, Phys. Rev. B 89, 245121 (2014).
  * Behrends  _et al._ [2016b] J. Behrends, A. G. Grushin, T. Ojanen, and J. H. Bardarson, Phys. Rev. B 93, 075114 (2016b).
  * Giuliani and Vignale [2005] G. Giuliani and G. Vignale,  _Quantum Theory of the Electron Liquid_ (Cambridge University Press, 2005).
  * Chang  _et al._ [2015] H.-R. Chang, J. Zhou, S.-X. Wang, W.-Y. Shan, and D. Xiao, Phys. Rev. B 92, 241103 (2015).
  * Hosseini and Askari [2015] M. V. Hosseini and M. Askari, Phys. Rev. B 92, 224435 (2015).
  * Kaladzhyan  _et al._ [2019] V. Kaladzhyan, A. A. Zyuzin, and P. Simon, Phys. Rev. B 99, 165302 (2019).
  * Godbout  _et al._ [2020] L. J. Godbout, S. Tchoumakov, and W. Witczak-Krempa, Phys. Rev. B 102, 205202 (2020).
  * Pierret [1996] R. Pierret,  _Semiconductor Device Fundamentals_ (Addison-Wesley, 1996) pp. 200–234.
  * Inoue  _et al._ [2016] H. Inoue, A. Gyenis, Z. Wang, J. Li, S. W. Oh, S. Jiang, N. Ni, B. A. Bernevig, and A. Yazdani, Science 351, 1184 (2016).
  * Batabyal  _et al._ [2016] R. Batabyal, N. Morali, N. Avraham, Y. Sun, M. Schmidt, C. Felser, A. Stern, B. Yan, and H. Beidenkopf, Science Advances 2 (2016), 10.1126/sciadv.1600709.
  * Kourtis  _et al._ [2016] S. Kourtis, J. Li, Z. Wang, A. Yazdani, and B. A. Bernevig, Phys. Rev. B 93, 041109 (2016).
  * Yuan  _et al._ [2019] Q.-Q. Yuan, L. Zhou, Z.-C. Rao, S. Tian, W.-M. Zhao, C.-L. Xue, Y. Liu, T. Zhang, C.-Y. Tang, Z.-Q. Shi, Z.-Y. Jia, H. Weng, H. Ding, Y.-J. Sun, H. Lei, and S.-C. Li, Science Advances 5 (2019), 10.1126/sciadv.aaw9485.
  * Zhang and Butler [1995] X.-G. Zhang and W. H. Butler, Phys. Rev. B 51, 10085 (1995).
  * Baranger and Stone [1989] H. U. Baranger and A. D. Stone, Phys. Rev. B 40, 8169 (1989).
  * LaGasse and Cress [2020] S. W. LaGasse and C. D. Cress, Nano Letters 20, 6623 (2020), pMID: 32787176.
  * Cano  _et al._ [2017] J. Cano, B. Bradlyn, Z. Wang, M. Hirschberger, N. P. Ong, and B. A. Bernevig, Phys. Rev. B 95, 161306 (2017).
  * Song  _et al._ [2021] J. Song, S. Kim, Y. Kim, H. Fu, J. Koo, Z. Wang, G. Lee, J. Lee, S. H. Oh, J. Bang, T. Matsushita, N. Wada, H. Ikegami, J. D. Denlinger, Y. H. Lee, B. Yan, Y. Kim, and S. W. Kim, Phys. Rev. X 11, 021065 (2021).
  * Wang  _et al._ [2012] Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu, H. Weng, X. Dai, and Z. Fang, Phys. Rev. B 85, 195320 (2012).
  * Zhang  _et al._ [2017b] C. Zhang, E. Zhang, W. Wang, Y. Liu, Z.-G. Chen, S. Lu, S. Liang, J. Cao, X. Yuan, L. Tang, _et al._ , Nature communications 8, 1 (2017b).
  * Morrow [1987] R. A. Morrow, Phys. Rev. B 35, 8074 (1987).
  * Einevoll and Hemmer [1988] G. T. Einevoll and P. C. Hemmer, Journal of Physics C: Solid State Physics 21, L1193 (1988).
  * Pfirsch and Sudan [1991] D. Pfirsch and R. N. Sudan, Journal of Mathematical Physics 32, 1774 (1991).

##  Appendix A Green’s function of a junction between two non-relativistic
electron gas

We consider a normal electron gas without spin where the dynamics are
described by the Schrödinger Hamiltonian

$\displaystyle\hat{H}=\frac{k^{2}}{2m}-\mu,$ (6)

which depends on two parameters: the band mass $m$ and the chemical potential
$\mu$. Note that in this equation we have set units $\hbar=1$.

Since the mass term is also position dependent, we symmetrize this Hamiltonian
and use [44, 45]

$\displaystyle\hat{H}=-\frac{1}{2}\nabla\left(\frac{1}{m(z)}\nabla\right)-\mu.$
(7)

The current operator is defined as

$\displaystyle\hat{\bf J}=\frac{e\hbar}{2mi}\psi^{\dagger}({\bf
r})\left(\nabla-\nabla^{\dagger}\right)\psi({\bf r}).$ (8)

###  A.1 Bulk Green’s function

It is useful to write three forms of the Green’s function: in momentum space
$G^{\infty}(i\omega_{n},{\bf k})$, in real space $G^{\infty}(i\omega_{n},{\bf
x},z)$ and in mixed representation $G^{\infty}(i\omega_{n},{\bf
k}_{\parallel},k_{z})$. We define the mixed representation by considering
momentum space along the interface and real space transverse to it. The
$\infty$ symbol on the Green’s functions is used to remind that it correspond
to the solution for a bulk sample, where parameters are constant over all
space.

  1. (i)

in momentum space

$\displaystyle G^{\infty}(i\omega_{n},{\bf
k};m,\mu)=\frac{1}{\frac{q^{2}}{2m}-\mu-i\omega_{n}},$ (9)

  2. (ii)

in real space

$\displaystyle G^{\infty}(i\omega_{n},{\bf r};m,\mu)$
$\displaystyle=\frac{1}{(2\pi)^{3}}\int_{0}^{\infty}dqq^{2}~{}\int_{0}^{\pi}d\theta\sin(\theta)~{}d\phi\frac{e^{iqr\cos(\theta)}}{\frac{q^{2}}{2m}-\mu-i\omega_{n}}=\frac{m}{2\pi^{2}}\int_{0}^{\infty}dqq^{2}~{}\int_{-1}^{1}du~{}\frac{e^{iqru}}{q^{2}-q_{0}^{2}},$
(10)
$\displaystyle=\frac{-im}{2\pi^{2}r}\int_{0}^{\infty}dqq~{}\frac{e^{iqr}}{q^{2}-q_{0}^{2}}=\frac{m}{2\pi
r}e^{i{\rm sgn}(m\omega_{n})q_{0}r},$ (11)

where $q_{0}=\sqrt{2m(\mu+i\omega_{n})}$.

  3. (iii)

in mixed representation

$\displaystyle G^{\infty}(i\omega_{n},{\bf
k}_{\parallel},z;m,\mu)=\frac{1}{2\pi}\int_{0}^{\infty}dq\frac{e^{iqr}}{\frac{q^{2}}{2m}+\frac{k_{\parallel}^{2}}{2m}-\mu-i\omega_{n}}=\frac{m}{\pi}\int_{0}^{\infty}dq\frac{e^{iqr}}{q^{2}-q^{2}(k_{\parallel})}=\frac{i{\rm
sign}(m\omega_{n})m}{q(k_{\parallel})}e^{i{\rm
sgn}(m\omega_{n})q(k_{\parallel})r},$ (12)

where $q({\bf k}_{\parallel})=\sqrt{2m(\mu+i\omega_{n})-k_{\parallel}^{2}}$.

###  A.2 Junction Green’s function

In the situation of a junction where parameters change abruptly from
$(m_{L},\mu_{L})$ at $z<0$ to $(m_{R},\mu_{R})$ at $z>0$, one can decompose
the corresponding Green’s function over the bulk Green’s function to satisfy
the boundary conditions. The boundary conditions are obtained by integrating
the Schrödinger equation on a small segment around the interface and
necessitate the continuity of the Green’s function to be well defined

$\displaystyle G_{L}(i\omega_{n},{\bf
k}_{\parallel},z,0)=G_{R}(i\omega_{n},{\bf k}_{\parallel},z,0),$ (13)
$\displaystyle\frac{1}{m_{1}}\partial_{z^{\prime}}G_{L}(i\omega_{n},{\bf
k}_{\parallel},z,0)=\frac{1}{m_{2}}\partial_{z^{\prime}}G_{R}(i\omega_{n},{\bf
k}_{\parallel},z,0).$ (14)

where we introduce the Green’s functions $G_{L}(i\omega_{n},{\bf
k}_{\parallel},z,z^{\prime})$ and $G_{R}(i\omega_{n},{\bf
k}_{\parallel},z,z^{\prime})$ defined for $z^{\prime}<0$ and $z^{\prime}>0$
respectively, i.e. left and right side of the junction. In order to compute
the Green’s functions that satisfy the boundary condition (14), we expand the
Green’s functions $G_{L}$ and $G_{R}$ over the Green’s functions
$G_{L}^{\infty}$ and $G_{R}^{\infty}$ of the corresponding bulk material, i.e.
the Green’s functions where parameters in Eq. (6) are constants equal to
($m_{L},\mu_{L}$) and ($m_{R},\mu_{R}$) respectively. For example, for $z<0$
we write

$\displaystyle G_{L}(i\omega_{n},{\bf
k}_{\parallel},z,z^{\prime})=G^{\infty}_{L}(i\omega_{n},{\bf
k}_{\parallel},z-z^{\prime})+RG^{\infty}_{L}(i\omega_{n},{\bf
k}_{\parallel},z+z^{\prime}),$ (15) $\displaystyle G_{R}(i\omega_{n},{\bf
k}_{\parallel},z,z^{\prime})=TG^{\infty}_{R}(i\omega_{n},{\bf
k}_{\parallel},z-z^{\prime}),$ (16)

with $R$ and $T$ constants in $z^{\prime}$ but functions in $z$. The boundary
condition (14) then implies

$\displaystyle G^{\infty}_{L}(i\omega_{n},{\bf
k}_{\parallel},z)+RG^{\infty}_{L}(i\omega_{n},{\bf
k}_{\parallel},z)=TG^{\infty}_{R}(i\omega_{n},{\bf k}_{\parallel},z)$ (17)
$\displaystyle\frac{1}{m_{L}}\partial_{z}G^{\infty}_{L}(i\omega_{n},{\bf
k}_{\parallel},z)-R\frac{1}{m_{L}}\partial_{z}G^{\infty}_{L}(i\omega_{n},{\bf
k}_{\parallel},z)=T\frac{1}{m_{R}}\partial_{z}G^{\infty}_{R}(i\omega_{n},{\bf
k}_{\parallel},z),$ (18)

which leads to

$\displaystyle R$
$\displaystyle=\frac{J_{L}^{\infty}/G_{L}^{\infty}-J_{R}^{\infty}/G_{R}^{\infty}}{J_{L}^{\infty}/G_{L}^{\infty}+J_{R}^{\infty}/G_{R}^{\infty}},$
(19) $\displaystyle T$
$\displaystyle=\frac{2J_{L}^{\infty}/G_{R}^{\infty}}{J_{R}^{\infty}/G_{R}^{\infty}+J_{L}^{\infty}/G_{L}^{\infty}},$
(20)

where we introduce the current function
$J^{\infty}=\frac{1}{m}\partial_{z}G^{\infty}(i\omega_{n},{\bf
k}_{\parallel},z)$. More explicitly this leads to reflection and transmission
coefficients that remind of Fresnel coefficients in optics

$\displaystyle
R=\frac{q_{L}(k_{\parallel})/|m_{L}|-q_{R}(k_{\parallel})/|m_{R}|}{q_{L}(k_{\parallel})/|m_{L}|+q_{R}(k_{\parallel})/|m_{R}|},$
(21) $\displaystyle
T=\frac{2q_{R}(k_{\parallel})/|m_{R}|}{q_{L}(k_{\parallel})/|m_{L}|+q_{R}(k_{\parallel})/|m_{R}|}e^{i{\rm
sign}(\omega_{n})({\rm sign}(m_{L})q_{L}(k_{\parallel})-{\rm
sign}(m_{R})q_{R}(k_{\parallel}))|z|},$ (22)

with
$q_{i}(k_{\parallel})=\sqrt{2m_{i}(\mu_{i}+i\omega_{n})-k_{\parallel}^{2}}$.

We obtain the Green’s function for $z>0$ following the same procedure and
obtain similar expression with a permutation of the $L\leftrightarrow R$
indices.

We compute Friedel oscillations and non-local conductivity by performing a
numerical Fourier transform of the above Green’s function, which is in
momentum space, to have the Green’s function in real space. In order to
simplify the calculation, we perform part of the Fourier transform
analytically, using rotation invariance, for example for $z<0$ we have

$\displaystyle G_{L}(i\omega_{n},{\bf
r},z,z^{\prime})=G^{\infty,1}(i\omega_{n},{\bf
r},z-z^{\prime})+\frac{1}{2\pi}\int
dk_{\parallel}k_{\parallel}RG^{\infty,1}(i\omega_{n},{\bf
k}_{\parallel},z+z^{\prime})J_{0}(kr_{\parallel}),$ (23) $\displaystyle
G_{R}(i\omega_{n},{\bf r},z,z^{\prime})=\frac{1}{2\pi}\int
dk_{\parallel}k_{\parallel}TG^{\infty,2}(i\omega_{n},{\bf
k}_{\parallel},z-z^{\prime})J_{0}(kr_{\parallel}).$ (24)

###  A.3 Surface states

We can look for the surface states solution as poles of the Green’s function
and, more specifically, as poles of the reflection and transmission
coefficients. The denominators of $R$ and $T$ are proportional to

$\displaystyle
D(i\omega_{n},k_{\parallel})=\sqrt{2m_{L}(\mu_{L}+i\omega_{n})-k_{\parallel}^{2}}/|m_{L}|+\sqrt{2m_{R}(\mu_{R}+i\omega_{n})-k_{\parallel}^{2}}/|m_{R}|,$
(25)

and the poles are obtained as a solution to
$D(\omega+i0^{+},k_{\parallel})=0$. The solution is obtained only if the
argument of both square roots is negative and with opposite imaginary parts,
that is if we satisfy the two conditions

  1. $\bullet$

${\rm sgn}(m_{L}m_{R})<0$: we have a $pn$ junction,

  2. $\bullet$

$2m_{L}(\mu_{L}+\omega)<k_{\parallel}^{2}$ and
$2m_{R}(\mu_{R}+\omega)<k_{\parallel}^{2}$: which does not imply any condition
on $\mu_{L}$ and $\mu_{R}$.

Then one obtains

$\displaystyle\frac{k_{\parallel}^{2}-2m_{L}(\mu_{L}+\omega)}{m_{L}^{2}}=\frac{k_{\parallel}^{2}-2m_{R}(\mu_{R}+\omega)}{m_{R}^{2}}$
(26) $\displaystyle\implies$
$\displaystyle\omega_{k_{\parallel}}=\frac{1}{2}\bigg{(}\frac{1}{m_{L}}+\frac{1}{m_{R}}\bigg{)}k_{\parallel}^{2}-\frac{\mu_{R}/m_{R}-\mu_{L}/m_{L}}{1/m_{R}-1/m_{L}}.$
(27)

We observe that we have a flat band behaviour (no dispersion of the surface
states, diverging density of states) if $m_{L}=-m_{R}$, which is the case of
perfect Veselago lensing, so the surface states do not disperse. Note that the
effective chemical potential of surface states is an average of the chemical
potential on each side of the interface.

The consistency condition for a pole in $R$ and $T$ is valid only for in-plane
wavevector such that

$\displaystyle k_{\parallel}^{2}>{\rm
min}\bigg{(}m_{L}\frac{\mu_{L}/m_{L}-\mu_{R}/m_{R}}{1/m_{L}-1/m_{R}}-2\mu_{L},-m_{R}\frac{\mu_{L}/m_{L}-\mu_{R}/m_{R}}{1/m_{L}-1/m_{R}}-2\mu_{R}\bigg{)},$
(28)

so the spectrum of surface states has a hole at small wavevectors. Also the
decay of surface states in each bulk is given by
$q_{1,2}(\omega)=\sqrt{2m_{1,2}(\mu_{1,2}+\omega_{k_{\parallel}}+i0^{+})-k_{\parallel}^{2}}$,
with opposite signs in each bulk.

##  Appendix B Green’s function for a junction between two Weyl semimetals

We consider a Dirac electron gas with Hamiltonian

$\displaystyle\hat{H}_{W}=\sum_{i=1}^{3}v_{i}(k_{i}-\delta
k_{i})\hat{\sigma}_{i}-\mu,$ (29)

which depends on the three Dirac velocities $v_{x},v_{y},v_{z}$, the three
momentum shifts $\delta k_{x},\delta k_{y},\delta k_{z}$ and the chemical
potential $\mu$. Note that in this equation we relate units of time and energy
by setting $\hbar=1$. Here, we define the current operator as

$\displaystyle\hat{J}=ev_{F}\Psi^{\dagger}({\bf r})\hat{\sigma}\Psi({\bf r}).$
(30)

###  B.1 Bulk Green’s function

It is useful to write three forms of the Green’s function: in momentum space
$G^{\infty}(i\omega_{n},{\bf k})$, in real space $G^{\infty}(i\omega_{n},{\bf
x},z)$ and in mixed representation $G^{\infty}(i\omega_{n},{\bf
k}_{\parallel},k_{z})$. As for the 3D electron gas we define the mixed
representation by considering momentum space along the interface and real
space transverse to it. The $\infty$ symbol on the Green’s functions is used
to remind that it correspond to the solution for a bulk sample, where
parameters are constant over all space.

  1. (i)

in momentum space

$\displaystyle G^{\infty}(i\omega_{n},{\bf
k};m,\mu)=\frac{\mu+i\omega_{n}+\hat{H}_{0}}{\sum_{i}h_{i}^{2}-(i\omega_{n}+\mu)^{2}}.$
(31)

  2. (ii)

in real space

$\displaystyle G^{\infty}(i\omega_{n},{\bf r})$
$\displaystyle=\frac{1}{(2\pi)^{3}}\int d^{3}{\bf
k}~{}\frac{\mu+i\omega_{n}+\hat{H}_{0}}{\sum_{i}h_{i}^{2}-(i\omega_{n}+\mu)^{2}}e^{i{\bf
k}\cdot{\bf r}}$ (32) $\displaystyle=(\mu+i\omega_{n})I_{0}({\bf
r})\hat{\mathbbm{1}}+v_{x}(-i\partial_{x}-\delta k_{x})I_{0}({\bf
r})\hat{\sigma}_{x}+v_{y}(-i\partial_{y}-\delta k_{y})I_{0}({\bf
r})\hat{\sigma}_{y}-iv_{z}\partial_{z}I_{0}({\bf r})\hat{\sigma}_{z},$ (33)

with

$\displaystyle I_{0}({\bf r})$ $\displaystyle=\frac{1}{(2\pi)^{3}}\int
d^{3}{\bf k}~{}\frac{e^{i{\bf k}\cdot{\bf
r}}}{\sum_{i}h_{i}^{2}-(\mu+i\omega_{n})^{2}}=\frac{e^{i(\delta k_{x}x+\delta
k_{y}y)}}{4\pi|v_{x}v_{y}v_{z}|\tilde{r}}e^{i{\rm
sgn}(\omega_{n})(\mu+i\omega_{n})\tilde{r}},$ (34)

where $\tilde{r}=\sqrt{(x/v_{x})^{2}+(y/v_{y})^{2}+(z/v_{z})^{2}}$. From this
we obtain,

$\displaystyle G^{\infty}(i\omega_{n},{\bf
r})=\left[(\mu+i\omega_{n})\hat{\mathbbm{1}}+\frac{i}{\tilde{r}^{2}}\left(1-i{\tilde{r}}{\rm
sgn}(\omega_{n})(\mu+i\omega_{n})\right){\tilde{\bf
r}\cdot\hat{\sigma}}\right]I_{0}({\bf r}).$ (35)

  3. (iii)

in mixed representation

$\displaystyle G^{\infty}(i\omega_{n},{\bf
k}_{\parallel},z)=\frac{1}{2\pi}\int_{\mathbbm{R}}dk_{z}~{}\frac{i\omega_{n}+\hat{H}_{D}}{\sum_{i}v_{i}^{2}k_{i}^{2}-(\mu+i\omega_{n})^{2}}e^{ik_{z}z}=\frac{1}{2\pi}\int_{\mathbbm{R}}dk_{z}~{}\frac{\mu+i\omega_{n}+v_{x}k_{x}\hat{\sigma}_{x}+v_{y}k_{y}\hat{\sigma}_{y}+v_{z}k_{z}\hat{\sigma}_{z}}{v_{z}^{2}k_{z}^{2}-((\mu+i\omega_{n})^{2}-(v_{x}^{2}k_{x}^{2}+v_{y}^{2}k_{y}^{2}))}e^{ik_{z}z},$
(36)

which gives

$\displaystyle G^{\infty}(i\omega_{n},{\bf
k}_{\parallel},z)=\frac{i}{2v_{z}}\left(\frac{(\mu+i\omega_{n})+v_{x}k_{x}\hat{\sigma}_{x}+v_{y}k_{y}\hat{\sigma}_{y}}{v_{z}q({\bf
k}_{\parallel})}+{\rm sgn}(z)\hat{\sigma}_{z}\right)e^{iq({\bf
k}_{\parallel})|z|},$ (37)

where $q({\bf k}_{\parallel})={\rm
sgn}(\omega_{n}\mu)\sqrt{(\mu+i\omega_{n})^{2}-(v_{x}^{2}k_{x}^{2}+v_{y}^{2}k_{y}^{2})}\big{/}|v_{z}|$.

###  B.2 Junction Green’s function

In the case the location of cones in momentum space and the chemical potential
change abruptly at $z=0$, from $(v_{L},\delta{\bf k}_{L},\mu_{L})$ for $z<0$
to $(v_{R},\delta{\bf k}_{R},\mu_{R})$ for $z>0$, the Hamiltonian

$\displaystyle\hat{H}_{W}$ $\displaystyle=\sum_{i=1}^{2}v_{i}(k_{i}-\delta
k_{i}(z))\hat{\sigma}_{i}-\mu(z),$ (38)

can be studied separately for $z<0$ and $z>0$. For this reason, we introduce
the Green’s functions $G_{L}(i\omega_{n},{\bf k}_{\parallel},z,z^{\prime})$
and $G_{R}(i\omega_{n},{\bf k}_{\parallel},z,z^{\prime})$ defined for
$z^{\prime}<0$ and $z^{\prime}>0$ respectively, i.e. on the left and right
side of the junction.

The integration of the Weyl equation (38) on a small segment around the
interface, at $z=0$, leads to the following boundary condition between $G^{L}$
and $G^{R}$

$\displaystyle G^{L}(i\omega_{n},{\bf
k}_{\parallel},z,0)=G^{R}(i\omega_{n},{\bf k}_{\parallel},z,0).$ (39)

The Green’s function of the junction can be expanded over the bulk Green’s
functions, $G_{L}^{\infty}$ and $G_{R}^{\infty}$, that are respectively the
Green’s functions for a bulk material with parameters $(v_{L},\delta{\bf
k}_{L},\mu_{L})$ and $(v_{R},\delta{\bf k}_{R},\mu_{R})$. This is similar to
method of images in electromagnetism and we use it to satisfy the boundary
condition (39). In the following we illustrate this expansion for $z<0$

$\displaystyle G_{L}(i\omega_{n},{\bf
k}_{\parallel},z,z^{\prime})=G^{\infty}_{L}(i\omega_{n},{\bf
k}_{\parallel},z-z^{\prime})+\hat{R}_{-}G^{\infty}_{L}(i\omega_{n},{\bf
k}_{\parallel},z+z^{\prime})=G^{\infty}_{L}(i\omega_{n},{\bf
k}_{\parallel},z-z^{\prime})+G^{\infty}_{L}(i\omega_{n},{\bf
k}_{\parallel},z+z^{\prime})\hat{R}_{+},$ (40) $\displaystyle
G_{R}(i\omega_{n},{\bf
k}_{\parallel},z,z^{\prime})=\hat{T}_{-}G^{\infty}_{R}(i\omega_{n},{\bf
k}_{\parallel},z-z^{\prime})=G^{\infty}_{R}(i\omega_{n},{\bf
k}_{\parallel},z-z^{\prime})\hat{T}_{+},$ (41)

where we introduce the reflection and transmission operators $\hat{R}_{\pm}$
and $\hat{T}_{\pm}$. We introduce the notation $\pm$ in order to account for
the fact that, in general, $\hat{R}$ and $\hat{T}$ do not commute with
$G^{\infty}$.

We simplify our calculation by expanding the reflection and transmission
operators as

$\displaystyle R_{\pm}$
$\displaystyle=\alpha_{0\pm}\hat{\mathbbm{1}}+\alpha_{1\pm}{\bf
h}_{\parallel,L}\cdot\hat{\sigma},$ (42) $\displaystyle T_{\pm}$
$\displaystyle=\big{[}\beta_{0\pm}\hat{\mathbbm{1}}+\beta_{1\pm}{\bf
h}_{\parallel,R}\cdot\hat{\sigma}\big{]}\frac{v_{z,L}}{v_{z,R}}e^{i(q_{L}({\bf
k}_{\parallel})-q_{R}({\bf k}_{\parallel}))z},$ (43)

where ($\alpha_{0\pm},\alpha_{1\pm},\beta_{0\pm},\beta_{1\pm})$ are solve from
the boundary condition (39). Here we also introduce ${\bf
h}_{\parallel}=\sum_{i=1}^{2}v_{i}(k_{i}-\delta k_{i}){\bf e}_{i}$, the in-
plane component of the Hamiltonian.

The boundary condition is satisfied for the following set of coefficients

$\displaystyle\alpha_{0+}$ $\displaystyle=\frac{(\tilde{\bf h}_{L}-\tilde{\bf
h}_{R})\cdot\tilde{\bf h}_{R}-\tilde{q}_{R}i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}{\tilde{\bf h}_{L}^{2}\tilde{\bf
h}_{R}^{2}-(1+\tilde{q}_{L}\tilde{q}_{R})\tilde{\bf h}_{L}\cdot\tilde{\bf
h}_{R}+(\tilde{q}_{L}+\tilde{q}_{R})i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}\tilde{\bf h}_{L}^{2}v_{z,L},$ (44) $\displaystyle\alpha_{1+}$
$\displaystyle=\frac{-i}{q_{L}(0)}\frac{(\tilde{\bf h}_{L}-\tilde{\bf
h}_{R})\cdot\tilde{\bf h}_{R}-\tilde{q}_{R}i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}{\tilde{\bf h}_{L}^{2}\tilde{\bf
h}_{R}^{2}-(1+\tilde{q}_{L}\tilde{q}_{R})\tilde{\bf h}_{L}\cdot\tilde{\bf
h}_{R}+(\tilde{q}_{L}+\tilde{q}_{R})i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}\tilde{\bf q}_{L}v_{z,L},$ (45) $\displaystyle\beta_{0+}$
$\displaystyle=-\frac{(\tilde{q}_{L}+\tilde{q}_{R})\tilde{\bf
h}_{L}\cdot\tilde{\bf h}_{R}-(1+\tilde{q}_{L}\tilde{q}_{R})i\tilde{\bf
h}_{L}\times\tilde{\bf h}_{R}}{\tilde{\bf h}_{L}^{2}\tilde{\bf
h}_{R}^{2}-(1+\tilde{q}_{L}\tilde{q}_{R})\tilde{\bf h}_{L}\cdot\tilde{\bf
h}_{R}+(\tilde{q}_{L}+\tilde{q}_{R})i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}\tilde{\bf q}_{R}v_{z,R},$ (46) $\displaystyle\beta_{1+}$
$\displaystyle=\frac{-i}{q_{R}(0)}\frac{(\tilde{\bf h}_{L}-\tilde{\bf
h}_{R})\cdot\tilde{\bf h}_{L}+\tilde{q}_{L}i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}{\tilde{\bf h}_{L}^{2}\tilde{\bf
h}_{R}^{2}-(1+\tilde{q}_{L}\tilde{q}_{R})\tilde{\bf h}_{L}\cdot\tilde{\bf
h}_{R}+(\tilde{q}_{L}+\tilde{q}_{R})i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}\tilde{\bf q}_{R}v_{z,R},$ (47)

or for

$\displaystyle\alpha_{0-}$ $\displaystyle=\frac{(\tilde{\bf h}_{L}-\tilde{\bf
h}_{R})\cdot\tilde{\bf h}_{R}+\tilde{q}_{R}i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}{\tilde{\bf h}_{L}^{2}\tilde{\bf
h}_{R}^{2}-(1+\tilde{q}_{L}\tilde{q}_{R})\tilde{\bf h}_{L}\cdot\tilde{\bf
h}_{R}-(\tilde{q}_{L}+\tilde{q}_{R})i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}\tilde{\bf h}_{L}^{2}v_{z,L},$ (48) $\displaystyle\alpha_{1-}$
$\displaystyle=\frac{i}{q_{L}(0)}\frac{(\tilde{\bf h}_{L}-\tilde{\bf
h}_{R})\cdot\tilde{\bf h}_{R}+\tilde{q}_{R}i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}{\tilde{\bf h}_{L}^{2}\tilde{\bf
h}_{R}^{2}-(1+\tilde{q}_{L}\tilde{q}_{R})\tilde{\bf h}_{L}\cdot\tilde{\bf
h}_{R}-(\tilde{q}_{L}+\tilde{q}_{R})i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}\tilde{\bf q}_{L}v_{z,L},$ (49) $\displaystyle\beta_{0-}$
$\displaystyle=-\frac{(\tilde{q}_{L}+\tilde{q}_{R})\tilde{\bf
h}_{L}\cdot\tilde{\bf h}_{R}+(1+\tilde{q}_{L}\tilde{q}_{R})i\tilde{\bf
h}_{L}\times\tilde{\bf h}_{R}}{\tilde{\bf h}_{L}^{2}\tilde{\bf
h}_{R}^{2}-(1+\tilde{q}_{L}\tilde{q}_{R})\tilde{\bf h}_{L}\cdot\tilde{\bf
h}_{R}-(\tilde{q}_{L}+\tilde{q}_{R})i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}\tilde{\bf q}_{R}v_{z,R}$ (50) $\displaystyle\beta_{1-}$
$\displaystyle=\frac{i}{q_{R}(0)}\frac{(\tilde{\bf h}_{L}-\tilde{\bf
h}_{R})\cdot\tilde{\bf h}_{L}-\tilde{q}_{L}i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}{\tilde{\bf h}_{L}^{2}\tilde{\bf
h}_{R}^{2}-(1+\tilde{q}_{L}\tilde{q}_{R})\tilde{\bf h}_{L}\cdot\tilde{\bf
h}_{R}-(\tilde{q}_{L}+\tilde{q}_{R})i\tilde{\bf h}_{L}\times\tilde{\bf
h}_{R}}\tilde{\bf q}_{R}v_{z,R}.$ (51)

where the tilde notation implies that a quantity is normalized by
$q_{0,L/R}=\mu_{L/R}+i\omega_{n}$ (for example $\tilde{\bf h}_{L}={\bf
h}_{L}/q_{0,L}$ or $\tilde{\bf h}_{R}={\bf h}_{R}/q_{0,R}$). We also use the
simplified notation ${\bf h}_{L}\times{\bf h}_{R}\equiv({\bf h}_{L}\times{\bf
h}_{R})\cdot{\bf e}_{z}$. All these expressions are obtained for $z<0$; if
instead $z>0$ then we should switch $i\rightarrow-i$.

The Fourier transforms of $G^{\mp}({\bf k}_{\parallel})$ to position-space are
performed numerically and we use the expressions above to draw Figs. 2 and 3.

###  B.3 Smooth junction Green’s function in the WKB approximation

In the previous section we solved the Green’s function of an abrupt interface
between regions with different parameter $(v,\delta{\bf k},\mu)$ in the
Hamiltonian (29). We now consider the situation where the parameters smoothly
vary from $(v,\delta{\bf k},\mu)_{L}$ for $z\rightarrow-\infty$ to
$(v,\delta{\bf k},\mu)_{R}$ for $z\rightarrow\infty$. The parameters are
supposed to be only $z$ dependent and in the following we solve for the
Green’s function of the Weyl equation (29) in the WKB approximation following
the method in [46].

We start by introducing the inverse of the thermal Green’s function with
smoothly varying parameters

$\displaystyle\hat{H}$
$\displaystyle=h_{x}(z)\hat{\sigma}_{x}+h_{y}(z)\hat{\sigma}_{y}+\bigg{(}-iv(z)\partial_{z}-\frac{i}{2}v^{\prime}(z)-\delta
k_{z}(z)\bigg{)}\hat{\sigma}_{z}-(\mu(z)+i\omega_{n}+i{\rm
sign}(\omega_{n})/(2\tau(z)))\hat{\mathbbm{1}},$ (52)

to which we search for an eigensolution of the form $\psi(z)=\Psi(z)e^{iS(z)}$
where the amplitude $\Psi(z)$ varies slowly compared to $e^{iS(z)}$. This then
satisfies

$\displaystyle\left(h_{x}(z)\hat{\sigma}_{x}+h_{y}(z)\hat{\sigma}_{y}+v(z)\bigg{(}k-\delta
k_{z}(z)\bigg{)}\hat{\sigma}_{z}-(\tilde{\mu}(z)+\lambda)\hat{\mathbbm{1}}\right)\Psi(z)=0,$
(53)

where $\tilde{\mu}(z)=\mu(z)+i\omega_{n}+i{\rm sign}(\omega_{n})/(2\tau(z))$,
$k=\partial S/\partial z$. We can diagonalize this equation on the states
$\phi_{\sigma=\pm}$ defined by the spinors

$\displaystyle\phi_{\sigma}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\left(1+\sigma\frac{h_{z}}{\delta
h}\right)^{1/2}\\\ \sigma\left(1-\sigma\frac{h_{z}}{\delta
h}\right)^{1/2}e^{i\theta},\end{array}\right)$ (56)

with $e^{i\theta}=(h_{x}+ih_{y})/|h_{\parallel}|$ which leads to the
projectors

$\displaystyle\phi_{\sigma}^{\dagger}(z^{\prime})\phi_{\sigma}(z)=\frac{1}{2}\left(\begin{array}[]{cc}\left(1+\sigma\frac{h_{z}}{\delta
h}\right)^{1/2}\left(1+\sigma\frac{h_{z}^{\prime}}{\delta
h^{\prime}}\right)^{1/2}&\sigma\left(1+\sigma\frac{h_{z}}{\delta
h}\right)^{1/2}\left(1-\sigma\frac{h_{z}^{\prime}}{\delta
h^{\prime}}\right)^{1/2}e^{-i\theta^{\prime}}\\\
\sigma\left(1-\sigma\frac{h_{z}}{\delta
h}\right)^{1/2}\left(1+\sigma\frac{h_{z}^{\prime}}{\delta
h^{\prime}}\right)^{1/2}e^{i\theta}&\left(1-\sigma\frac{h_{z}}{\delta
h}\right)^{1/2}\left(1-\sigma\frac{h_{z}^{\prime}}{\delta
h^{\prime}}\right)^{1/2}e^{i(\theta-\theta^{\prime})}\end{array}\right),$ (59)

where $\delta h=(h_{\parallel}^{2}+h_{z}^{2})^{1/2}$ with
$h_{z}=v_{z}(\partial_{z}S-\delta k_{z})$. We then have for each eigensolution
of Eq. (53)

$\displaystyle\frac{\partial S_{\sigma}}{\partial z}=k_{\sigma}=\delta
k_{z}(z)+\frac{h(z)}{v_{z}(z)}=\delta
k_{z}+\frac{1}{v_{z}}\left((\tilde{\mu}+\lambda_{\sigma})-h_{\parallel}^{2}\right)^{1/2}$
(60) $\displaystyle\implies$ $\displaystyle
S(z)=\int^{z}_{z_{0}}du~{}\left(\delta
k_{z}(u)+\frac{1}{v_{z}(u)}((\tilde{\mu}(u)+\lambda_{\sigma})^{2}-h_{\parallel}^{2}(u))^{1/2}\right).$
(61)

The phase factor $S(z)$ is defined up to a constant of integration, associated
to the starting position $z_{0}$, and that we replace with parameter

$\displaystyle a_{\sigma}\equiv k_{\sigma}(z_{0})=\delta
k_{0}+\frac{1}{v_{z,0}}((\tilde{\mu}_{0}+\lambda_{\sigma})^{2}-h_{\parallel,0}^{2})^{1/2},$
(62)

and write the eigenvalues $\lambda_{\sigma}$ as a function of $a_{\sigma}$

$\displaystyle\lambda_{\sigma}=-\tilde{\mu}_{0}+\sigma\delta
h_{0}=-\tilde{\mu}_{0}+\sigma(h_{x,0}^{2}+h_{y,0}^{2}+v_{z,0}^{2}(a_{\sigma}^{2}-\delta
k_{0})^{2})^{1/2}.$ (63)

In this approximation we discard the trajectories that cycle, i.e. the
trajectory is uniquely defined by the parameter $a$. Also, since there is a
one-to-one correspondence between $\lambda$ and $a$, we can use either
parameter to define the trajectory.

It can be shown that the amplitude $\Psi(z)$ is given by the van Vleck
determinant, usually related to the local density of states of semiclassical
trajectories [46]. When parametrizing the trajectory with the eigenvalue
$\lambda$, we thus have

$\displaystyle\Psi_{\sigma}(z)=\left(\frac{\partial S}{\partial
z\partial\lambda_{\sigma}}\right)^{1/2}=\left(\frac{1}{v_{z}(z)}\frac{\tilde{\mu}(z)+\lambda_{\sigma}}{\left((\tilde{\mu}(z)+\lambda_{\sigma})^{2}-h_{\parallel}^{2}(z)\right)^{1/2}}\right)^{1/2}.$
(64)

Then the WKB Green’s function is obtained by integrating the previous solution
over all eigenvalues $\lambda$

$\displaystyle G(i\omega_{n},{\bf
k}_{\parallel},z,z^{\prime})=\frac{1}{2}\sum_{\sigma}\int\frac{d\lambda_{\sigma}}{2\pi}~{}\frac{1}{\lambda_{\sigma}}\phi_{\sigma}^{\dagger}(z^{\prime})\phi_{\sigma}(z)\Psi_{\sigma}(z^{\prime})\Psi_{\sigma}(z)e^{i\int^{z^{\prime}}_{z}du~{}\left(\delta
k_{z}(u)+\frac{1}{v_{z}(u)}((\tilde{\mu}(u)+\lambda_{\sigma})^{2}-h_{\parallel}^{2}(u))^{1/2}\right)},$
(65)

where the factor $1/2$ appears because we consider a bi-spinor [46]. We can
compute the above integral with the Cauchy theorem by choosing a contour such
that the integrrand is negligible at large $|\lambda|$, that is such that
${\rm Im}(S(\lambda))>0\implies\partial S/\partial z\sim{\rm
sign}((z^{\prime}-z)v_{z}(u)\mu(u)\omega_{n})$. We obtain

$\displaystyle G(i\omega_{n},{\bf
k}_{\parallel},z,z^{\prime})=\frac{i}{2}e^{i\int_{z}^{z^{\prime}}du~{}\left(\delta
k_{z}(u)+{\rm
sign}\left[(z^{\prime}-z)\omega_{n}v_{z}(u)\mu(u)\right]\frac{1}{v_{z}(u)}\left(\tilde{\mu}^{2}(u)-h_{\parallel}^{2}(u)\right)^{1/2}\right)}\times$
(66)
$\displaystyle\left(\begin{array}[]{cc}\left(a(z,z^{\prime})/v_{z}^{\prime}\right)^{1/2}\left(b(z,z^{\prime})/v_{z}\right)^{1/2}&\left(1/(a(z,z^{\prime})v_{z}^{\prime}h_{z}^{\prime
2})\right)^{1/2}\left(b(z,z^{\prime})/{v_{z}}\right)^{1/2}\left(h_{x}^{\prime}-ih_{y}^{\prime}\right)\\\
\left(a(z,z^{\prime})/v_{z}^{\prime}\right)^{1/2}\left(1/(b(z,z^{\prime})v_{z}h_{z}^{2})\right)^{1/2}\left(h_{x}+ih_{y}\right)&\left(1/(a(z,z^{\prime})v_{z}^{\prime}h_{z}^{\prime
2})\right)^{1/2}\left(1/(b(z,z^{\prime})v_{z}h_{z}^{2})\right)^{1/2}\left(h_{x}^{\prime}-ih_{y}^{\prime}\right)\left(h_{x}+ih_{y}\right)\end{array}\right),$
(69)

where

$\displaystyle h_{z}(z)={\rm
sign}(\omega_{n}\mu(z)v_{z}(z))\left(\tilde{\mu}^{2}(z)-h_{\parallel}^{2}(z)\right)^{1/2},$
(70) $\displaystyle\tilde{\mu}(z)=\mu(z)+i\omega_{n}+i{\rm
sign}(\omega_{n})/(2\tau(z)),$ (71) $\displaystyle
a(z,z^{\prime})=\frac{\tilde{\mu}_{z^{\prime}}}{h_{z}^{\prime}}+{\rm
sign}(z^{\prime}-z),$ (72) $\displaystyle
b(z,z^{\prime})=\frac{\tilde{\mu}_{z}}{h_{z}}+{\rm sign}(z^{\prime}-z).$ (73)

Note that we cannot further simplify these expressions since the square root
has a branch cut on $(-\infty,0]$ so $\sqrt{ab}=\sqrt{a}\sqrt{b}$ is only
defined if

$\displaystyle\begin{array}[]{lc}{\rm(1):}&{\rm Im}(a){\rm Im}(b)\leq 0,\\\
{\rm(2):}&{\rm Im}(a){\rm Im}(b)>0~{}{\rm and}~{}{\rm Im}(a){\rm
Im}(ab)>0.\end{array}$ (76)

##  Appendix C Current response in the bulk

In this section we compute the non-local conductivity in various situation for
a planar electric field ${\bf E}({\bf x})={\bf E}_{0}\delta(z-z_{0})$ and
measuring the current in a plane $j_{\mu}({\bf
x})=j_{\mu}(z_{1})\delta(z-z_{1})$, such that
$j_{\mu}(z_{1})=\sigma_{\mu\mu}(z_{1},z_{0})E_{0\mu}$.

###  C.1 Non-local conductivity and time-reversal symmetry

The conductivity tensor can always be decomposed over a symmetric,
$\sigma_{1}$, and an antisymmetric, $\sigma_{2}$, component with respect to
time-reversal symmetry [38]

$\displaystyle\sigma_{1,\mu\nu}$
$\displaystyle=-\sigma_{0}\int_{\mathbbm{R}}d\epsilon~{}\left(-\frac{df}{d\epsilon}\right){\rm
Tr}\left[\hat{J}_{\mu}\left(G_{+}({\bf x},{\bf x}^{\prime})-G_{-}({\bf x},{\bf
x}^{\prime})\right)\hat{J}_{\nu}\left(G_{+}({\bf x}^{\prime},{\bf
x})-G_{-}({\bf x},{\bf x})\right)\right]$ (77)
$\displaystyle\sigma_{2,\mu\nu}$
$\displaystyle=\sigma_{0}\int_{\mathbbm{R}}d\epsilon~{}\left(-\frac{df}{d\epsilon}\right)\left\\{{\rm
Tr}\left[\hat{J}_{\mu}G_{+}({\bf x},{\bf x}^{\prime})\hat{J}_{\nu}G_{-}({\bf
x}^{\prime},{\bf x})\right]-{\rm Tr}\left[\hat{J}_{\mu}G_{-}({\bf x},{\bf
x}^{\prime})\hat{J}_{\nu}G_{+}({\bf x}^{\prime},{\bf x})\right]\right\\}$
$\displaystyle-\sigma_{0}\int_{\mathbbm{R}}d\epsilon~{}f(\epsilon)\left\\{{\rm
Tr}\left[\hat{J}_{\mu}\partial_{\epsilon}G_{+}({\bf x},{\bf
x}^{\prime})\hat{J}_{\nu}G_{+}({\bf x}^{\prime},{\bf x})\right]+{\rm
Tr}\left[\hat{J}_{\mu}G_{-}({\bf x},{\bf
x}^{\prime})\hat{J}_{\nu}\partial_{\epsilon}G_{-}({\bf x}^{\prime},{\bf
x})\right]\right.~{}$
$\displaystyle\left.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-{\rm
Tr}\left[\hat{J}_{\mu}\partial_{\epsilon}G_{-}({\bf x},{\bf
x}^{\prime})\hat{J}_{\nu}G_{-}({\bf x}^{\prime},{\bf x})\right]-{\rm
Tr}\left[\hat{J}_{\mu}G_{+}({\bf x},{\bf
x}^{\prime})\hat{J}_{\nu}\partial_{\epsilon}G_{+}({\bf x}^{\prime},{\bf
x})\right]\right\\},$ (78)

with $G_{\pm}=(\hat{H}\pm i0^{+})^{-1}$. In general $\sigma_{2}$ is not zero
but one can check that, for $\mu=\nu$, there is the symmetry

$\displaystyle\sigma_{2,\mu\mu}({\bf x},{\bf
x}^{\prime})=-\sigma_{2,\mu\mu}({\bf x}^{\prime},{\bf x}),$ (79)

so

$\displaystyle\sigma_{1,\mu\mu}({\bf x},{\bf
x}^{\prime})=\frac{1}{2}\left(\sigma_{\mu\mu}({\bf x},{\bf
x}^{\prime})+\sigma_{\mu\mu}({\bf x}^{\prime},{\bf x})\right).$ (80)

That is, longitudinal components of the conductivity tensor can be obtained by
averaging the longitudinal conductivities obtained when permuting the two
probes. In the main text and below we focus on $\sigma_{1,\mu\mu}$ and drop
the reference to the index $1$.

###  C.2 Kubo formalism

The expression of the non-local conductivity is obtained in the Kubo formalism
in [37, 38]. In the situation of a planar electric field ${\bf E}({\bf
x})={\bf E}_{0}\delta(z-z_{0})$ and measuring the current in a plane
$j_{\mu}({\bf x})=j_{\mu}(z_{1})\delta(z-z_{1})$, such that
$j_{\mu}(z_{1})=\sigma_{\mu\mu}(z_{1},z_{0})E_{0\mu}$, we obtain

$\displaystyle\sigma_{\mu\nu}(z,z^{\prime})=\int\dfrac{dS_{z}dS_{z^{\prime}}}{\pi\mathcal{A}}{\rm
Tr}\left(\hat{j}_{\mu}{\rm Im}\hat{G}(z,z^{\prime})\hat{j}_{\nu}{\rm
Im}\hat{G}(z,z^{\prime})\right).$ (81)

where the current operator for the Weyl equation is $\hat{\bf
j}=v(z)\hat{\bm{\sigma}}$ and where ${\rm
Im}\hat{G}(z,z^{\prime})\equiv\frac{1}{2}\left(\hat{G}(i0^{+},{\bf r},{\bf
r}^{\prime})-\hat{G}(i0^{-},{\bf r},{\bf r}^{\prime})\right)$ with

$\displaystyle\hat{G}(i0^{+},{\bf r},{\bf r}^{\prime})-\hat{G}(i0^{-},{\bf
r},{\bf r}^{\prime})=\frac{1}{4\pi v^{2}|{\bf r}-{\bf r}^{\prime}|}$
$\displaystyle\left[\tilde{\mu}e^{i\tilde{\mu}|{\bf r}-{\bf
r}^{\prime}|/|v|}\hat{\mathbbm{1}}-\tilde{\mu}^{*}e^{-i\tilde{\mu}^{*}|{\bf
r}-{\bf r}^{\prime}|/|v|}\hat{\mathbbm{1}}+i\frac{v({\bf r}-{\bf
r}^{\prime})}{|{\bf r}-{\bf
r}^{\prime}|^{2}}\cdot\hat{\sigma}\left(e^{i\tilde{\mu}|{\bf r}-{\bf
r}^{\prime}|/|v|}-e^{-i\tilde{\mu}|{\bf r}-{\bf
r}^{\prime}|/|v|}\right)\right.$ $\displaystyle\left.+{\rm sign}(v)\frac{{\bf
r}-{\bf r}^{\prime}}{|{\bf r}-{\bf
r}^{\prime}|}\cdot\hat{\sigma}\left(\tilde{\mu}e^{i\tilde{\mu}|{\bf r}-{\bf
r}^{\prime}|/|v|}+\tilde{\mu}^{*}e^{-i\tilde{\mu}|{\bf r}-{\bf
r}^{\prime}|/|v|}\right)\right],$ (82)

where we use the notation $\tilde{\mu}=\mu+i/(2\tau)$. Since the Green’s
function is translation invariant in the $(x,y)$ plane, we can easily get rid
of one of the surface integrals in Eq. (81). We perform the second surface
integral in polar coordinates and obtain

$\displaystyle\sigma_{xx}$ $\displaystyle=-\frac{2\hbar v_{x}^{2}}{(4\pi
v^{2})^{2}\pi}\int\frac{dR}{R}d\theta~{}\left[\left(\tilde{\mu}e^{i\tilde{\mu}R/|v|}-\tilde{\mu}^{*}e^{-i\tilde{\mu}^{*}R/|v|}\right)^{2}+\frac{-X^{2}+Y^{2}+Z^{2}}{R^{2}}\left(\left(\tilde{\mu}+\frac{i|v|}{R}\right)e^{i\tilde{\mu}R/|v|}+\left(\tilde{\mu}^{*}-\frac{i|v|}{R}\right)e^{-i\tilde{\mu}^{*}R/|v|}\right)^{2}\right].$
(83)

where $(X,Y,Z)={\bf r}-{\bf r}^{\prime}$ and $X=R\cos(\theta)$,
$Y=R\cos(\theta)$. This expression shows that $\sigma_{xx}=\sigma_{yy}$ so we
can write $\sigma_{xx}=\frac{1}{2}(\sigma_{xx}+\sigma_{yy})$ and get rid of
the angular integral by removing the $X^{2}-Y^{2}$ term. We then have

$\displaystyle\sigma_{xx}$ $\displaystyle=-\frac{4\hbar v_{x}^{2}}{(4\pi
v^{2})^{2}}\int\frac{dR}{R}~{}\left[\left(\tilde{\mu}e^{i\tilde{\mu}R/|v|}-\tilde{\mu}^{*}e^{-i\tilde{\mu}^{*}R/|v|}\right)^{2}+\frac{Z^{2}}{R^{2}}\left(\left(\tilde{\mu}+\frac{i|v|}{R}\right)e^{i\tilde{\mu}R/|v|}+\left(\tilde{\mu}^{*}-\frac{i|v|}{R}\right)e^{-i\tilde{\mu}^{*}R/|v|}\right)^{2}\right]$
(84) $\displaystyle=-\frac{2\hbar v_{x}^{2}}{(4\pi
v^{2})^{2}}\int_{1}^{\infty}\frac{du}{u}~{}\left[\left(\tilde{\mu}e^{i\tilde{\mu}Zu/|v|}-\tilde{\mu}^{*}e^{-i\tilde{\mu}^{*}Zu/|v|}\right)^{2}+\frac{1}{u^{2}}\left(\left(\tilde{\mu}+\frac{i|v|}{Zu}\right)e^{i\tilde{\mu}Zu/|v|}+\left(\tilde{\mu}^{*}-\frac{i|v|}{Zu}\right)e^{-i\tilde{\mu}^{*}Zu/|v|}\right)^{2}\right]$
(85) $\displaystyle=-\frac{4\hbar v_{x}^{2}}{(4\pi
v^{2})^{2}}\int_{1}^{\infty}{du}~{}\left(\frac{\tilde{\mu}^{2}e^{i(2\tilde{\mu}Z/|v|)u}+\tilde{\mu}^{*2}e^{i(2\tilde{\mu}^{*}Z/|v|)u}-2|\tilde{\mu}|^{2}e^{i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|u}}{u}\right.$
$\displaystyle\left.+\frac{\tilde{\mu}^{2}e^{i(2\tilde{\mu}Z/|v|)u}+\tilde{\mu}^{*2}e^{i(2\tilde{\mu}^{*}Z/|v|)u}+2|\tilde{\mu}|^{2}e^{i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|u}}{u^{3}}+\frac{2i|v|}{Z}\frac{\tilde{\mu}e^{2i\tilde{\mu}Z/|v|u}-\tilde{\mu}^{*}e^{-2i\tilde{\mu}^{*}Z/|v|u}+(\tilde{\mu}^{*}-\tilde{\mu})e^{i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|u}}{u^{4}}\right.$
$\displaystyle\left.-\frac{v^{2}}{Z^{2}}\frac{e^{2i\tilde{\mu}Z/|v|u}+e^{-2i\tilde{\mu}^{*}Z/|v|u}-2e^{i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|u}}{u^{5}}\right)$
(86) $\displaystyle=-\frac{4\hbar v_{x}^{2}}{(4\pi
v^{2})^{2}}\left[\tilde{\mu}^{2}(E_{1}(-2i\tilde{\mu}Z/|v|)+E_{3}(-2i\tilde{\mu}Z/|v|))+\tilde{\mu}^{*2}(E_{1}(2i\tilde{\mu}^{*}Z/|v|)+E_{3}(2i\tilde{\mu}^{*}Z/|v|))\right.$
$\displaystyle\left.+2|\tilde{\mu}|^{2}(E_{3}(-i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|)-E_{1}(-i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|))+\frac{2i|v|}{Z}\left(\tilde{\mu}E_{4}(-2i\tilde{\mu}Z/|v|)-\tilde{\mu}^{*}E_{4}(2i\tilde{\mu}^{*}Z/|v|)\right)\right.$
$\displaystyle\left.+(\tilde{\mu}^{*}-\tilde{\mu})E_{4}(-i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|))-\frac{v^{2}}{Z^{2}}\left(E_{5}(-2i\tilde{\mu}Z/|v|)+E_{5}(2i\tilde{\mu}^{*}Z/|v|)-2E_{5}(-i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|)\right)\right].$
(87)

We can also compute $\sigma_{zz}$ as a function of $Z=|z-z^{\prime}|$ and
obtain

$\displaystyle\sigma_{zz}=-\frac{2\hbar v_{x}^{2}}{(4\pi
v^{2})^{2}}\int\frac{dR}{R}~{}\left[\left(\tilde{\mu}e^{i\tilde{\mu}R/|v|}-\tilde{\mu}^{*}e^{-i\tilde{\mu}^{*}R/|v|}\right)^{2}+\left(1-2\frac{Z^{2}}{R^{2}}\right)\left(\left(\tilde{\mu}+\frac{i|v|}{R}\right)e^{i\tilde{\mu}R/|v|}+\left(\tilde{\mu}^{*}-\frac{i|v|}{R}\right)e^{-i\tilde{\mu}^{*}R/|v|}\right)^{2}\right]$
(88) $\displaystyle=-\frac{2\hbar v_{x}^{2}}{(4\pi
v^{2})^{2}}\int_{1}^{\infty}\frac{du}{u}~{}\left[\left(\tilde{\mu}e^{i\tilde{\mu}Z/|v|u}-\tilde{\mu}^{*}e^{-i\tilde{\mu}^{*}Z/|v|u}\right)^{2}+\left(1-\frac{2}{u^{2}}\right)\left(\left(\tilde{\mu}+\frac{i|v|}{Zu}\right)e^{i\tilde{\mu}Z/|v|u}+\left(\tilde{\mu}^{*}-\frac{i|v|}{Zu}\right)e^{-i\tilde{\mu}^{*}Z/|v|u}\right)^{2}\right]$
(89) $\displaystyle=-\frac{2\hbar v_{x}^{2}}{(4\pi
v^{2})^{2}}\int_{1}^{\infty}du~{}\left(\left({\tilde{\mu}^{2}e^{2i\tilde{\mu}Z/|v|u}+\tilde{\mu}^{*2}e^{-2i\tilde{\mu}^{*}Z/|v|u}}\right)\left(\frac{2}{u}-\frac{2}{u^{3}}\right)-\frac{4\tilde{\mu}\tilde{\mu}^{*}e^{-i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|u}}{u^{3}}\right.$
$\displaystyle+\frac{2i|v|}{Z}\left({\tilde{\mu}e^{i\tilde{\mu}Z/|v|u}-\tilde{\mu}^{*}e^{-i\tilde{\mu}^{*}Z/|v|u}+(\tilde{\mu}^{*}-\tilde{\mu})e^{i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|u}}\right)\left(\frac{1}{u^{2}}-\frac{2}{u^{4}}\right)$
$\displaystyle\left.-\frac{v^{2}}{Z^{2}}\left({e^{2i\tilde{\mu}Z/|v|u}+e^{-2i\tilde{\mu}^{*}Z/|v|u}-2e^{i\delta\tilde{\mu}Z/|v|u}}\right)\left(\frac{1}{u^{3}}-\frac{2}{u^{5}}\right)\right)$
(90) $\displaystyle=-\frac{2\hbar v_{x}^{2}}{(4\pi
v^{2})^{2}}\left(2\tilde{\mu}^{2}(E_{1}(-2i\tilde{\mu}Z/|v|)-E_{3}(-2i\tilde{\mu}Z/|v|))+2\tilde{\mu}^{*2}(E_{1}(2i\tilde{\mu}^{*}Z/|v|)-E_{3}(2i\tilde{\mu}^{*}Z/|v|))-2|\tilde{\mu}|^{2}E_{3}(-i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|)\right.$
$\displaystyle\left.+\frac{2i|v|}{Z}\left[\tilde{\mu}(E_{2}(2i\tilde{\mu}Z/|v|)-2E_{4}(2i\tilde{\mu}Z/|v|))-\tilde{\mu}^{*}(E_{2}(-2i\tilde{\mu}^{*}Z/|v|)-2E_{4}(-2i\tilde{\mu}^{*}Z/|v|))\right.\right.$
$\displaystyle\left.\left.+(\tilde{\mu}^{*}-\tilde{\mu})(E_{2}(-i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|)-2E_{4}(-i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|))\right]-\frac{v^{2}}{Z^{2}}(E_{3}(-2i\tilde{\mu}Z/|v|)-2E_{5}(-2i\tilde{\mu}Z/|v|)\right.$
$\displaystyle\left.+E_{3}(2i\tilde{\mu}^{*}Z/|v|)-2E_{5}(2i\tilde{\mu}^{*}Z/|v|)-2E_{3}(-i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|))+2\times
2E_{5}(-i(\tilde{\mu}-\tilde{\mu}^{*})Z/|v|))\right).$ (91)

We show the behaviour of the non-local conductivity obtained from these two
expressions in Fig. S1. These expressions also match with our numerical
results for the non-local conductivity far from the junction.

###  C.3 Semiclassical limit.

In the semiclassical limit the wave-like nature of quasiparticles is neglected
such that one can approximate the exponentials with complex arguments with a
functions with the same decay length and same volume integral but without the
quantum oscillations

$\displaystyle e^{i\alpha Zu}\rightarrow\frac{i{\rm
Im}(\alpha)}{\alpha}e^{-{\rm Im}(\alpha)Zu}$ (92)

from which we also deduce the semiclassical limits of terms involving powers
of $1/Z=1/|z-z^{\prime}|$ by integrating this equation over $u$ to have

$\displaystyle\frac{e^{i\alpha Zu}}{Z^{n}}\rightarrow\left(\frac{\alpha}{i{\rm
Im}(\alpha)}\right)^{n-1}\frac{e^{-{\rm Im}(\alpha)Zu}}{Z^{n}}.$ (93)

This semi-classical limit is evaluated for a non-relativistic electron gas in
Ref. [37], and for a Weyl electron gas we obtain

$\displaystyle\sigma_{xx}^{\rm cl.}(Z)$ $\displaystyle=\frac{\hbar
v_{x}^{2}}{(4\pi v^{2})^{2}}\left(\left({\rm Re}(\tilde{\mu})^{2}+2{\rm
Im}(\tilde{\mu})^{2}\right)E_{1}(2{\rm Im}(\tilde{\mu})Z/|v|)-{\rm
Re}(\tilde{\mu})^{2}E_{3}(2{\rm Im}(\tilde{\mu})Z/|v|)\right),$ (94)
$\displaystyle\sigma_{zz}^{\rm cl.}(Z)$ $\displaystyle=\frac{2\hbar
v_{x}^{2}}{(4\pi v^{2})^{2}}\left({\rm Im}(\tilde{\mu})^{2}E_{1}(2{\rm
Im}(\tilde{\mu})Z/|v|)+{\rm Re}(\tilde{\mu})^{2}E_{3}(2{\rm
Im}(\tilde{\mu})Z/|v|)\right),$ (95)

that indeed match the conductivity in the Kubo formula but without the
oscillating behaviour (see Fig. S1).

Figure S1: Bulk non-local conductivity of the Weyl electron gas for a planar
electric field at $z=0$. (a) $\sigma_{xx}$ and (b) $\sigma_{zz}$ are obtained
for $(v_{F},\mu,\ell)=(1,1,50)$, and are shown in units of the total bulk
conductivity $\sigma_{0}$. We show the non-local conductivities from the Kubo
formula in blue, and from the semi-classical approximation in orange, the two
expression match relatively well up to the oscillations in $\sigma_{xx}$.

###  C.4 Total conductivity

The global conductivity tensor is obtained after performing the integral over
all the volume

$\displaystyle\langle\sigma_{\mu\mu}\rangle=\frac{1}{L_{z}}\int
dzdz^{\prime}~{}\sigma_{\mu\mu}(z,z^{\prime}).$ (96)

In the configuration with an electric field in the $z$ direction, we obtain

  1. 1.

For the non-relativistic electron gas:

$\displaystyle\sigma_{xx}=\sigma_{yy}=\frac{1}{8\pi^{2}}\frac{{\rm
Re}(q)^{2}-m\mu/3}{{\rm Im}(q)}\frac{e^{2}}{\hbar},$ (97)
$\displaystyle\sigma_{zz}=\frac{1}{6\pi^{2}}\frac{\mu}{{\rm
Im}(q)}\frac{e^{2}}{\hbar}$ (98)

so the total conductivity is

$\displaystyle\sigma_{0}=\frac{1}{3}(\sigma_{xx}+\sigma_{yy}+\sigma_{zz})=\frac{1}{12\pi^{2}}\frac{{\rm
Re}(q)^{2}}{{\rm Im}(q)}\frac{e^{2}}{\hbar}$ (99)

which is obtained in Ref. [37].

  2. 2.

For the Weyl electron gas:

$\displaystyle\sigma_{xx}=\sigma_{yy}=\sigma_{zz}=\frac{3+4\mu^{2}\tau^{2}}{6\pi^{2}|v|\tau}\frac{e^{2}}{\hbar}$
(100)

so the total conductivity is

$\displaystyle\sigma_{0}=\frac{1}{3}(\sigma_{xx}+\sigma_{yy}+\sigma_{zz})=\frac{3+4\mu^{2}\tau^{2}}{6\pi^{2}|v|\tau}\frac{e^{2}}{\hbar}.$
(101)

##  Appendix D Conditions for optimal Veselago lensing

In this section we explore how the smoothness of a $pn$ junction and how the
mismatch in carrier densities affect the amplitude of Veselago lensing. We
also show devices that realize a chiral Veselago lens and how they can be used
to explore Veselago lensing from local sources.

###  D.1 Veselago lensing for a smooth junction

In the main text we discuss the situation of a sharp $pn$ junction, based on
the absence of a built-in potential in chiral Veselago lens. In the situation
of a smooth junction, with a characteristic length $d$ we instead have the
behaviour reproduced in Fig. S2, obtained for $k_{F}d=3/2$ with other
parameters as in the main text $v_{F,L}=v_{F,R}=1$, $\mu_{L}=-\mu_{R}=1$ and
for a mean free path $\ell=50/k_{F}$. We see that the amplitude of the
Veselago image is strongly reduced compared to the image in an sharp
interface, even for such a small value of $d$.

Figure S2: (a) Friedel oscillations $\Pi(z,z^{\prime})$ and (b) non-local
conductivity $\sigma_{xx}$ of a $pn$ junction for a smooth $pn$ junction, with
extension $d=3/(2k_{F})$. The density plot in (a) corresponds to the
polarisability of the smooth $pn$ junction of the Weyl electron gas, obtained
from the WKB approximation (see Sec. B.3).

###  D.2 Veselago lensing for different densities

In Figs. S3(a-c) we show the non-local conductivity $\sigma(z,z^{\prime})$ as
a function of the location of the input potential different, at $z$, and the
location of the output current, at $z^{\prime}$, for different values of the
Fermi wavevector for $z>0$ (see Fig. 2(d) in the main text).

We find that image current through the Veselago lens is largest when the
carrier density on each side of the interface is the same (see Fig. S3(a)) and
this is the situation we choose to illustrate in the main text. A $pn$
junction with different carrier densities still shows a peaked conductivity,
due to Veselago lensing, but the amplitude of this peak decays faster than for
identical carrier densities.

Figure S3: Non-local conductivity $\sigma_{xx}$ of a $pn$ junction with
different carrier densities in the regions $z<0$ and $z>0$, given by the Fermi
wavevectors $k_{F,L}$ and $k_{F,R}$ for (a) $k_{F,R}=k_{F,L}$, (b)
$k_{F,R}=0.75k_{F,L}$ and (c) $k_{F,R}=0.5k_{F,L}$. The non-local conductivity
along the dashed line in (a) corresponds to the data in Fig. 2(f) of the main
text.

###  D.3 Devices to probe chiral Veselago lensing

In the main text we consider that the chiral anomaly is induced locally from
an extended electrode that applies a uniform electric field for $z<0$, in
presence of a uniform magnetic field. This is depicted in Fig. 1 of the main
text, where we consider that the extended electric field acts as a source for
current.

Other configurations are also possible to explore Veselago lensing as a
consequence of the chiral anomaly and that we illustrate in Fig. S4(a,b). In
both figures the principle is the same, we consider that the input voltage or
electric field is applied on the surface normal to the extended electric field
responsible for the chiral anomaly.

Figure S4: Devices to probe chiral Veselago lensing. (a) A local potential is
applied with an STM tip (or impurity) on the surface of the device and leads
to charge fluctuations that are measured by a second tip. (b) A planar
electric field is applied on the surface of the device and leads to a current
distribution that resembles that in Fig. 2(e,f) of the main text. In both (a)
and (b), the side electrodes apply an electric field at $z<0$ that is colinear
with the magnetic field, inducing a $pn$-junction, and thus a Veselago lens,
for a single chirality.