Designer curved-space geometry for relativistic fermions in Weyl metamaterials

Weyl semimetals are recently discovered materials supporting emergent relativistic fermions in the vicinity of band-crossing points known as Weyl nodes. The positions of the nodes and the low-energy spectrum depend sensitively on the time-reversal (TR) and inversion (I) symmetry breaking in the system. We introduce the concept of Weyl metamaterials where the particles experience a 3d curved geometry and gauge fields emerging from smooth spatially varying TR and I breaking fields. The Weyl metamaterials can be fabricated from semimetal or insulator parent states where the geometry can be tuned, for example, through inhomogeneous magnetization. We derive an explicit connection between the effective geometry and the local symmetry-breaking configuration. This result opens the door for a systematic study of 3d designer geometries and gauge fields for relativistic carriers. The Weyl metamaterials provide a route to novel electronic devices as highlighted by a remarkable 3d electron lens effect.

In this work, we study a class of structures which we call Weyl metamaterials. In these systems, the chiral Weyl fermions familiar from Weyl semimetals, depicted in Fig. 1a, are moving in an artificial 3d curved-space geometry. The great advantage of these systems is that, first of all, the effective geometry experienced by the particles can be widely tuned by external conditions. Since the band crossings in Weyl semimetals result from accidental degeneracy, the positions of the nodes and the low-energy spectrum can be tuned by a large amount by varying time-reversal (TR) and inversion (I) symmetry breaking perturbations, as depicted in Fig. 1b. In addition, a 3d curved geometry can be realized by atomically smooth TR or I breakings whose variation is significant only on the scale of the sample size, as illustrated in Fig. 1c. These two factors combined with the obvious richness offered by 3d space give Weyl metamaterials flexibility and enormous advantage over 2d materials [5,19,20] as a platform for relativistic curved-space physics. Weyl metamaterials can be fabricated from 3d semimetal or insulator parent states by generic inhomogeneous TR and I breaking couplings. Hence the con- Figure 1: a: Weyl semimetals exhibit one or more bandcrossing pairs which act as sources and sinks of the Berry curvature. The low-energy spectrum near the nodes is described by the Weyl Hamiltonian H = ±vk · σ resulting in a conical relativistic dispersion. b: Locations of the Weyl nodes and the the conical envelopes depend on the magnitude of the TR and I breaking in the system and are not fixed by symmetry. For example, by varying the magnetization Mi, one can widely tune the spectrum. c: A smoothly varying magnetization in real space (bottom) give rise to local Weyl cones that are deformed depending on the texture. This translates to an effective curved-space geometry experienced by relativistic carriers (top). d: Semiclassical trajectories of carriers in the sample can be engineered by designing suitable TR and I breaking textures. cept reaches substantially beyond the strain-engineering through elastic deformations [21][22][23][24][25].
Starting from a generic four-band model, we derive a curved-space Weyl Hamiltonian H = e µ i (k µ − a µ )σ i describing the low-energy dynamics of Weyl metamaterials. The emergent gauge field a µ and the frame fields e µ i which encode the metric tensor g µν are solved as functions of the TR and I breaking fields. These expressions provide the fundamental tool for systematic reverse-engineering of synthetic curved-space geometries and gauge fields in Weyl metamaterials. Our results are applicable to a wide class of different lattice and continuum models as illustrated by examples. The general theory is illustrated by proposing simple magnetic textures giving rise to remarkable electron lensing effects. More generally, Weyl metamaterials pave the way for novel 3d electronic devices through curvature engineering as sketched in Fig. 1d.
Engineering Weyl metamaterials from general parent states. The central idea of our work is to start from a generic insulating or semimetal phase and introduce smooth spatially varying TR and I breaking (TRIB) fields that push the parent system into an inhomogeneous Weyl semimetal phase. The TR and I breaking terms are not assumed to be small in any sense, however their spatial dependence is assumed to be smooth. Besides an intrinsic Weyl semimetal, the parent state for a Weyl metamaterial could be a Dirac semimetal or TR invariant topological insulator [2,30] as it is well known that these systems generically undergo a transition to Weyl semimetals when the I [33,34] or TR symmetry is broken. A possible mechanism to break TR in topological insulators is magnetic doping, which has been demonstrated in thin films [35][36][37]. As the minimum models of Weyl semimetals have four energy bands, we write down the Clifford representation of the most general TR and I invariant four-band models and identify all TR and I breaking terms as was also done in Ref. [32]. Without any TRIB fields, the most general four-band Hamiltonian is where the parameters n, and m are symmetric in the momentum m(k) = m(−k), while the kinetic term is antisymmetric κ(k) = −κ(−k). We have denoted the unit matrix by I. Furthermore, γ µ denotes the five 4 × 4 gamma matrices satisfying the anticommutation relations {γ µ , γ ν } = 2δ µν . Since the Hamiltonian (1) respects TR and I symmetry, it follows that γ 1,2,3 are odd under both I and TR, while γ 4 must be even. From this, it is then evident that γ 5 = γ 1 γ 2 γ 3 γ 4 must also be odd under both operations. To write down the most general Hermitian 4 × 4 Hamiltonian, we need to include ten additional matrices constructed using commutators of the gamma matrices, i.e. γ ij ≡ −i[γ i , γ j ]/2. The symmetries of γ ij can be easily deduced from the symmetries of their constituent gamma matrices. As it turns out, we can group all ten matrices into four distinct sets: three sets consisting of three mutually anticommuting matrices each, and a fourth set consisting only of the matrix γ 45 . The sets, as well as the transformation properties of each set under TR or I, are summarized in table I -in fact, we can see that each group breaks either TR or I symmetry.
The general 4 × 4 Hamiltonian is expressed as where the functions u, w, u and f characterize the TRIB fields and fix the position of the Weyl nodes and the low- energy spectrum. To extract the low-energy Weyl Hamiltonian, we must block diagonalize (2). To make analytical progress, we concentrate on the case u = 0, f = 0. This restriction leaves room for substantial generality, allowing for a full three-component TR breaking field u and I breaking field w. To illustrate the wide applicability of the general formulation, we give in the Supplementary Information (SI) three concrete examples of Weyl metamaterials: one based on a topological insulator-ferromagnet layer structure [31], a popular theoretical model describing magnetically-doped topological insulators [29], and a 3d Dirac semimetal [9,10] with broken inversion symmetry .
Weyl block reduction. Here we derive the local Weyl Hamiltonian from the original four-band model (2). The block reduction is achieved by the unitary transformation theory developed in SI and is given by where Weyl nodes are now contained within either the upper or lower block, while the other block is gapped. The low-energy physics in the vicinity of the Weyl nodes is contained in the local Weyl approximation H W = d 0 (k)σ 0 + d(k) · σ. Here, we have collected the Pauli matrices σ 1,2,3 into a vector σ, and written the 2 × 2 identity matrix as σ 0 . For concreteness, let us relegate the general case u = 0, w = 0 to Sec. I of the SI and for now consider the specific case H = κ(k) · γ + mγ 4 + u · b, which could describe a TR and I invariant insulator or semimetal in the presence of spatially varying magnetization M = u. The local Weyl approximation yields d 0 = 0 and where u = |u| andû = u/u. This system has Weyl nodes k W satisfying κ(k W ) = ± √ u 2 − m 2û whenever m < |u|. The gapped block and the results for the general case u = 0 and w = 0 are given in the SI. In the local Weyl approximation H W , the d-vector depends not only on k but has an implicit position dependence through u which is slowly varying in space. The Weyl reduction (3), which is exact for position-independent TRIB fields, is an approximation when u and w depend on position. The unitary transformations that block diagonalizes the four-band model become spatially dependent through u(r), w(r). Regarding k in H 0 as a canonical conjugate operator of r, the transformation also produces terms proportional to gradients (as well as higher-order spatial derivatives) of u(r), w(r). The gradient terms could, for example, give rise to off-diagonal elements in Eq. (3). These terms would give corrections to H W of the order of O (∂ i u) 2 , (∂ i w) 2 which, along with other gradient corrections, are very small for slowly varying fields. This is confirmed by direct numerical comparisons between the local density of states (LDOS) of the exact four-band model and the local Weyl approximation H W , and typically suggest an excellent agreement between the two. The LDOS, defined by ν(r, E) = En |Ψ n (r)| 2 δ(E−E n ), probes the local properties of the system and therefore is an ideal tool to investigate the local Weyl approximation. We have calculated the LDOS for representative Weylmetamaterial models with H = κ(k) · γ + mγ 4 + u · b, where κ i (k) = t sin k i and m and t are constants and different TR-breaking textures u. This is illustrated in Fig. 2 which establishes an excellent agreement between the exact result and the local Weyl approximation in the energy window E g determined by the energy gap of the gapped block. Away from the Weyl node but within E g , the smoothed LDOS exhibits quadratic envelope typical for Weyl semimetals. At the Weyl node, the LDOS typ-ically becomes modified due to the emergent gauge field effect, as clarified in the next section.
Designer curved-space geometry for Weyl particles. To gain better insight into the low-energy physics, we expand H W around the local Weyl points. For each Weyl point, we get a low-energy Hamiltonian of the form Here we have employed the Einstein summation convention -Greek indices represent integers from 0 to 3, while Latin indices take integer values from 1 to 3, unless otherwise stated. We define the frame fields (or tetrads) as where k W (u, w) is the momentum corresponding to the local Weyl point around which we have expanded. Introducing a two-component spinor Ψ, the low-energy Schrödinger equation then takes the form H W Ψ = EΨ which when squared can be formally written in the form g µν p µ p ν Ψ = 0 with appropriately defined canonical momenta p µ . When the term proportional to the unit matrix vanishes d 0 = 0 in H W , this equation describes a massless relativistic particle with dispersion E 2 = g ij p i p j moving in a curved space. Here p i = k i − k W i and the contravariant metric g µν = η ab e µ a e ν b is defined by the standard relation [38] using tetrad fields e µ ν which coincide with those defined in Eq. (6) with the addition e 0 0 = 1. The covariant metric g µν is defined as the inverse of g µν and comes in to play in the semiclassical dynamics through Christoffel symbols.
The frame fields e i ν for the model (4) can be straightforwardly evaluated and are given in the SI. We can also easily find the inverse (co-frame) elementsẽ that satisfy e α νẽ µ α = δ µ ν and e µ αẽ α ν = δ µ ν . The contravariant metric tensor is given by where the elements of the 3 × 3 matrix [g] ij = g ij are evaluated at the local Weyl node k = k W . In the form (8) we have assumed that m is constant in the vicinity of k W , the more general case where m depends on k is given in the SI. The relations (4), (6), and (8), along with their generalizations given in the SI, provide the fundamental relation between the effective geometry and the physical symmetry-breaking configuration and serve as a starting point in designing synthetic curved-space geometries for the Weyl quasiparticles. The formula Eq. (8) has a similar role in Weyl metamaterials as the relation between the metric and electromagnetic permeability and permitivity tensors in transformation optics [26][27][28]. It shows concretely how a Weyl metamaterial acts as a transformation media where the geometry is tuned by u.
We postulate an effective long-wavelength quantum Hamiltonian by making the replacement k → −i∇ in Eq. (5). A direct substitution would leave H W non-Hermitian. This can be remedied by the standard prescription of symmetrizing the non-commuting product of gradient and the frame fields. Assuming for simplicity that d 0 = 0, we obtain an effective curved-space Weyl Hamiltonian Here we have defined an effective gauge field a = k W (r) which gets a small correction due to the symmetrization. Though the correction is parametrically small since it depends on the derivatives of the TRIB fields, we include it to preserve Hermiticity exactly. Numerical comparisons again show that the effective Hamiltonian (9) produces the low-energy LDOS of the exact four-band model accurately, as seen in Figs. 2c and d. The linear approximation also illuminates the behavior of the exact model. The flat part of the LDOS in the vicinity of the Weyl node appears due to a chiral mode appearing in the gap of the bulk Weyl cones due to the Landau quantization induced by the local effective magnetic fieldB = ∇×k W . The inset of Fig. 2c shows how a constant field of magnitude |B| give rise to a chiral mode representing the zeroth Landau level dispersing in the direction the field. The chiral mode has a 1d dispersion and a constant density of states. This interpretation is confirmed by the fact that the calculated LDOS is approximately constant in the energy window 2EB = 2vB |B| which correspond to the gap between −1 and 1 Landau levels where vB is the Fermi velocity along the magnetic field. In addition, as shown in Fig. 2d, for textures which the effective magnetic field vanishes, the LDOS displays a quadratic trend down to the Weyl node. Semiclassical dynamics. To study the quasiparticle dynamics, we employ the semiclassical theory [39]. We assume that the system is doped above the Weyl node and that the Landau quantization due to the real and effective gauge fields can be neglected. Then we can ignore interband dynamics and concentrate on a single positive energy band. The semiclassical wave-packet motion for inhomogeneous systems [39] is described by the equations of motioṅ where E = (g ij k i k j ) 1 2 is the relativistic curved-space dispersion and E and B are electric and magnetic fields including the emergent fieldB. Quantum particles have orbital degrees of freedom encoded in the frame fields, leading to quantum-geometric effects given by the Berry forces on top of the classical effects. The novelty compared to the previously extensively studied flat case is that, in addition to the curved-space dispersion, we must consider a 6 × 6 phase-space Berry-curvature tensor Ω instead of its 3 × 3 momentum block. In the studied model, the Berry tensor has components where d i are the coefficients of the Pauli matrices in Eq. (3). In the linear approximation, the Berry curvatures can be expressed in terms of the frame fields and their derivatives by substituting d j = e i j k i . In the lowest non-trivial order, we can drop the Ω r i r l term since it is the only term that involves two spatial derivatives and is parametrically small. Away from the Weyl nodes, the Berry curvatures vanish as Ω kik l ∼ 1 |k| 2 , Ω r i k l ∼ 1 |k| , so by increasing the doping, one recovers the classical geometric effects together with the emergent magnetic field. However, in general, both classical and quantum geometric effects need to be considered together. On the other hand, we have proved in the SI that in the absence of fields E = 0, B = 0, the ballistic motion is governed by the geodesic equationr l + Γ l ijṙ iṙj = 0 in the lowest nontrivial order in spatial derivatives. Here Γ l ij is the Christoffel symbol calculated from the metric g ij . The geodesic motion is a direct manifestation of the curved geometry and indicates how the carrier motion can be tailored through geometry. The carrier dynamics in disordered solids also include scattering processes, which could be implemented by incorporating the semiclassical dynamics into a Boltzmann equation. However, this goes beyond the scope of the present work where we illustrate the designer geometries by solving Eq. (10).
3D Weyl electron lens. The negative index of refraction in optical metamaterials enable fabrication of exotic devices [28] such as the Veselago lens [40] which can also be realized for electrons in graphene [41]. Inspired by these ideas, we propose a 3d Weyl electron lens as an application of Weyl metamaterials.
This phenomenon can be realized, for example, in a two-node model H = v i k i γ i + mγ 4 + u · b, which could represent a topological insulatorferromagnet sandwich structure or magnetically-doped topological insulator with slowly rotating magnetization u = u [cos ϕ sin θ(r), sin ϕ sin θ(r), cos θ(r)] expressed in spherical coordinates. The texture is parametrized by the function θ(r) = ωr and depicted in Fig. 3 a which illustrates the carrier trajectories in the setup. The carriers of one chirality entering the sample in the vicinity of r = 0 will converge periodically in the vicinity of focal points. This motion results from a combined effect of the curved geometry and the effective magnetic field that can be straightforwardly calculated in our general framework. In Fig. 3 b we have solved the trajectories from the semiclassical equations corresponding to initial conditions where parallel beams enter the sample at z = 0. Thus, the image of the object at infinity is reproduced at focal points. The situation in Fig. 3 c corresponds to the case where a point-like object is located at z = 0 and imagined at the focal points. The carriers with the opposite chirality experience an opposite magnetic field and diverge away from x = y = 0 axis as shown in Fig. 3 c. Thus the considered geometry will realize a chiralityselective 3d electron lens. The qualitative features of lens motion can be reproduced by analytical calculations neglecting the Berry curvature effects presented in the SI. However, the exact numerical solution shows that the motion is not exactly planar due to quantum geometric effects. Interestingly, a similar 3d lens effect can also be realized in unidirectional magnetic textures as discussed in the SI. We also discovered a 2d lens effect for a simpler magnetic texture that only rotates in plane. The 3d Weyl lens illustrates that Weyl metamaterials posses richness beyond the 2d materials considered for curvature and gauge field engineering. Interesting future challenges include, for example, finding realizations for 3d electronic invisibility devices.
Conclusion. In this work, we proposed Weyl metamaterials as a highly tuneable platform for relativistic fermions in curved space. We provided a theoretical description of these systems and established an explicit relation between the local TR and I breaking and the effective geometry and gauge fields experienced by the carriers. Our results can be employed in systematic reverseengineering of 3d curved geometries and studying quantum effects [42] in designer geometries. One interesting direction is to study the curved-space modifications to the exotic response properties such as the chiral anomaly and its manifestations in Weyl semimetals [43][44][45][46][47][48]. From the point of view of applications, Weyl metamaterials offer a novel route to 3d electronic devices through geometry engineering.
Acknowledgements. The authors acknowledge the Academy of Finland for support.
Author contributions. Both authors contributed to all aspects of this work.  In this supplement, we derive the unitary transformations that reduce the 4 × 4 Hamiltonians for different models into blocks of (2 × 2) Weyl Hamiltonians.
It turns out that it is beneficial to first consider the Hamiltonian in Eq. (2) with w = 0, i.e.
As can be easily shown by inspection, the matrices in b have identical properties to those of the regular 2 × 2 Pauli matrices, i.e.
These matrices, like the Pauli matrices, thus exponentiates to the special unitary group SU(2), which is two-to-one homomorphic to the orthogonal group O(3). In practice, this means that the transformation corresponds to rotating the vector u by an angle ϕ around the axis defined by the unit vectorn. Furthermore, by observing that the components of the other vectors in table I in the main text along with γ can be written as p = −γ 5 b, b = γ 4 b, and γ = −γ 45 b, respectively, we can conclude that the exponential map has the same effect (it commutes with γ 4 and γ 5 ) on terms containing these vectors as well. The first step is to rotate u in Eq. (I.1) to point along the x axis. This can for example be achieved by first rotating u into the xz plane and then rotating it around the y axis, aligning it with the x axis. These two transformations are effected by which is expressed in terms of the polar (θ) and azimuthal (ϕ) angles of u. Applying this transformation to the whole Hamiltonian yields To block-diagonalize this, we need to eliminate either γ 1 or γ 4 . This is achieved using R = e −i ρ 2 γ14 with tan ρ = (κ(k) ·û)/m. The final result is then which splits into blocks in the basis given by for example With this choice, the d vectors for the upper and lower block becomes Depending on the sign of m(k), the nodes will either be located in the upper or the lower block. The d vector corresponding to the Weyl block will regardless have the form where the potential minus sign on the first and second component have been taken care of using a unitary transformation σ z . The above vector vanishes whenever all three components are zero, i.e.
Furthermore, note that if d W (k W ) = 0, it follows immediately that d W (−k W ) = 0, which is just a manifestation of the fact that Weyl nodes come in pairs. To arrive at an effective curved-space Weyl Hamiltonian, we linearize d W around the Weyl nodes. The Weyl Hamiltonian can then be written in terms of frame fields as where the frame fields e j i can be found to be where the sign in the last component corresponds to either k W (upper sign) or −k W (lower sign). The co-frame fields e are easy to evaluate whenever e j 0 = e 0 j = 0, and are given bỹ Finally, the metric tensor can be straightforwardly calculated to be where we have have omitted writing out the arguments. The whole expression must then be evaluated at the relevant Weyl node.

I breaking with preserved TR
In this section we show how to obtain similar results for a system with only inversion-symmetry breaking. The starting Hamiltonian is (I.17) The last term can be rotated using U(w) to give U(w)HU † (w) = K · γ + m(k)γ 4 + wγ 14 , (I. 18) where K can be directly extracted from Eq. (I.5). On this we can then simply apply a rotation exp(iνγ 23 /2) with tan ν = K 3 /K 2 to get Employing, for example, the representation we see that this Hamiltonian is block diagonal with blocks given by It follows that the Weyl block is where, if sgn(K 2 (k W )) < 0, we apply a transformation σ z to get rid of the minuses in front of K 1 =ŵ · κ(k) and m.
The condition for the existence of the Weyl nodes k W become In this case, the minimum number of nodes is four. We can now linearize this around the Weyl points, yielding We can now immediately read off the frame fields to be e j 1 =ŵ · ∂ qj κ(q)| q=k W (I.25) and the metric tensor

Block reduction with TR and I breaking
The problem becomes significantly more cumbersome when both u and w are non-vanishing. In this case, it is useful to write w = w ⊥ + w , where w ⊥ · u = 0 and w u. The goal is to find a transformation that removes all matrices belonging to p from the Hamiltonian, such that we are left with the special case discussed in the first part of this appendix.
Before tackling the full problem, first consider the case where w ⊥ = 0. Then it follows that U(u) also aligns w with the x axis. The Hamiltonian is then U(u)HU † (u) = K · γ + mγ 4 + uγ 23 + wγ 14 , (I. 29) where K can be directly extracted from Eq. (I.5). We can now remove γ 14 by a simple unitary transformation F = e iδγ1/2 , provided that tan δ = w/m: (I. 30) Evidently, F also introduces new b terms, but at this point, we can simply refer back to the first half of this appendix, as the problem is mathematically identical to the case w = 0.
Solving the full problem can thus be done if we can find a transformation that aligns the vector associated with b and the vector associated with p. For the sake of clarity, we will not concern ourselves with that of H 0 until at the very end (as it turns out, all transformations merely rotates κ), so our primary focus will now lie on (I.31) As usual, we begin by applying U(u), which gives us We can then make a rotation around b 1 with the unitary transformation V = e iµb1/2 which for tan µ =w ⊥,3 /w ⊥,2 eliminates the p 3 term: At this point, we make two transformations in succession; we first apply V b = e iαb3/2 and then V p = e iβp3/2 which gives a seemingly unwieldy where w = sgn(w · u)w . We can now determine the angles α and β by requiring that the terms proportional to b 2 and p 2 vanish. One possible solution is then (I.34) From this then follows that to which we can apply our earlier methods. The effects of these transformations on H 0 are then straightforwardly calculated. As usual, U(u) gives us which upon application of V turns into VUH 0 U † V † = K 1 γ 1 + (K 2 cos µ + K 3 sin µ)γ 2 With the identifications q 1 = K 1 cos α + (K 2 cos µ + K 3 sin µ) sin α (I.39) we can write our full Hamiltonian as and refer back to the case with w u, which combined with our earlier discussion on parallel w and u gives us on which we can apply the machinery developed for the case of w = 0 (Eq. (I.1)). As a final comment on the general case, we note that the equations for the Weyl nodes now take the form κ(k ± W ) = ±rû + s(ŵ ×û) + t(û × (ŵ ×û)), Both of these cases reduce to the correct equation in the w = 0 limit.

II. Realizations of Weyl metamaterials
To illustrate the utility of our theoretical framework and the physics of Weyl metamaterials, we consider three much-studied models as platforms for Weyl metamaterials. These include the topological insulator-ferromagnet layer structure [31], the popular toy model introduced by Vazifeh and Franz [29], and a Dirac semimetal model of Cd 3 As 2 [9] with broken inversion symmetry. As detailed in Ref. [31], the topological insulator-ferromagnet heterostructure can be modelled by the Hamiltonian where v F is the Fermi velocity, k the momentum,∆(k z ) = ∆ S τ x + ∆ D (cos k z τ x − sin k z τ y ), and m is a time-reversal breaking magnetization. The parameters ∆ S and ∆ D characterize the tunnelling between different surfaces within and between neighboring layers of the heterostructure, respectively. While not unique, a suitable choice for the representation of the Dirac gamma matrices allows us to write Eq. (II.1) as where we have promoted the model from having magnetization only in the z direction mb 3 → m · b.
Here we see that m corresponds to the TR-breaking field u in Table I and letting it vary as a function of position, we obtain an example of a TR-breaking metamaterial. By applying our general framework, it is straightforward to study the curved-space low-energy physics of the model. We next consider the toy model introduced in [29]. The model represents a 3d (topological) insulator with a finite magnetization and has a Hamiltonian where we have slightly deviated from the notation in [29] for the TR and I breaking terms. Here, the mass term is M k = µ − 2t 3 α=1 cos k α , where µ and t are the chemical potential and hopping amplitude, respectively. The parameters λ, λ z characterize the strength of the spin-orbit hopping and u can be identified with magnetization. Expressing this Hamiltonian via Dirac gamma matrices gives us H = 2λ(γ 2 sin k y + γ 1 sin k x ) + 2λ z γ 3 sin k z + γ 4 M k + u 0 ε + u · b. (II.4) Within the scope of our present formalism, we cannot account for the I-breaking ε term and must hence restrict ourselves to u 0 = 0, but otherwise this is now in a form which we can apply our formalism to, starting from Eq. (I.1). For example, we immediately get that the equations for the Weyl nodes are given by where λ 1,2 = λ and λ 3 = λ z . The frame fields and the metric can be now extracted from our general formulas. Our third example is a 3D Dirac semimetal [9] with broken I symmetry. The Hamiltonian describes the spectrum of Cd 3 As 2 (or Na 3 Bi) near the Γ point and is given by . The above Hamiltonian can be directly translated into our gamma-matrix representation via the substitution τ z σ j = γ j for j = 1, 2 and τ z σ z = γ 4 . We can add an inversion-breaking field to this model by introducing a coupling w · p, such that we get H = ε 0 (k) + Ak x γ 1 + Ak y γ 2 + M (k)γ 4 + w · p. (II.7) While the above form with the I-breaking vector parameter w is a bit abstract, it is possible to implement it by, for example, elastic deformations due to strain [22][23][24][25]. Starting from Eq. (I.17), we can then then derive the desired quantities, viz. the frame fields and the induced metric.

III. Semiclassical motion and geodesic equation
The geodesic equation can be derived from the semiclassical equations of motion in Eq. (10). In the absence of electromagnetic fields and Berry curvatures, the equations reduce tȯ which -given that E is a constant of motion -can by a rescaling of the time t → εt be brought into the forṁ We recognize these as the Hamilton equations of motion, which follow from extremizing the action with the Lagrangian given by This problem is equivalent to extremizing the length of a path in a curved space given by the metric g ij , and it immediately follows that r obeys the geodesic equation with the connection Γ given by (III.6) Reintroducing the Berry curvatures changes this picture; in the presence of Berry curvatures, the energy is no longer a constant of motion, and would also seemingly make the semiclassical paths deviate from the geodesic equations. However, it turns out that to linear order in spatial derivatives, the Berry curvatures do not alter the picture meaningfully. We have thatṙ (III.7) If we neglect terms that are second order or higher in spatial derivatives, we can remove Ω rr and write (I + Ω kr ) −1 ≈ I − Ω kr . The equations can then be written aṡ By rewriting the first equation in component form and differentiating once more with respect to time, we havë where the last step follows from using the fact that E = g ij k i k j , and dropping terms with more than one spatial derivative in the last line. To show that this is equivalent to a geodesic equation to linear order, we will solve for k in Eq. (III.8):ṙ l = ∂ k l E − Ω k l r s ∂ ks E + Ω k l ks ∂ r s E = 1 E g lj k j − Ω k l r s g sj k j + 1 2 Ω k l ks (∂ r s g 1 ij)k i k j , (III. 11) or in matrix formṙ From this, we get which can then be inserted into Eq. (III.11), resulting in r l = (∂ r t g lj )g js + 1 2 g ln (∂ r n g ts ) ṙ sṙt −→r l + Γ l stṙ sṙt = 0, (III.14) where is the Christoffel symbol.

IV. Analytical solution for 3D lens geometry
In this appendix we present an analytical solution for the oscillatory behavior of the trajectories in the limit of small deviations from x 2 + y 2 = 0 and k x /k z , k y /k z 1. The spatial part of the metric for the 3D lens is given by where the omitted elements follow from the symmetry of g. For small θ, the metric reduces to From Eq. (IV.1), it is straightforward to calculate its spatial derivatives which in our approximation reduce to and ∂ z g = 0.
We can now write down the equations of motions using Eqs. (10): where we have ignored the contributions from the the Berry curvatures, as their effects are small on short enough time scales. We have additionally taken into account the effective magnetic field arising from a position-dependent Weyl node,B = ∇ × k W /e. Since this magnetic field acts with different signs on the opposite nodes, we have introduced the parameter c = ±1 to distinguish the Weyl nodes (c = 0 would correspond to pure geodesic motion). In the derivation of this linear differential equation, we have only kept terms up to linear order in k x and k y and used the fact thaṫ k z = 0. The solution to this differential equation can be formally written as Ψ(t) = exp(At)Ψ(t = 0), where A is the matrix on the right-hand side. From this, it is evident that if A has imaginary eigenvalues, we will obtain oscillatory solutions. To simplify, we now assume that v x = v y ≡ v. The eigenvalues are then from which it follows that we have sinusoidal solutions typically only for one chirality c = 1. Furthermore, sincė z ≈ v 2 z (1 − m 2 /u 2 )k z /E, we have that the length of one oscillation with period T = 2π/|λ| is (IV.7) The distance between two focal points is ∆z/2 which depends on the energy of particles weakly through k z in the square root. The analytical predictions are confirmed by exact numerics. The Berry-curvature effects that are neglected here affect the long time dynamics, most clearly by resulting in a non-planar motion that is not evident for short times. The analytical results are applicable also for different textures when ω is replaced by a characteristic measure of the variation of u z component near the lens axis.

V. Alternative lens realizations
In the main text, we considered an electron lens in a two-node model H = v i k i γ i + mγ 4 + u · b with a skyrmion-like magnetic texture. This is not necessary for the lens effect to occur as we can see in Fig. 1, which illustrates how a similar effect can be realized by a unidirectional texture u = [0, 0, u(r)]. If u(r) decreases radially away from the lens axis as in Fig. 1 a, the carrier trajectories near the axis are qualitatively similar to the rotating texture. However, even simpler textures result in nontrivial lensing; the case where u(r) varies only in the x direction, as shown in Fig. 1 b, still leads to a 2d lens effect. It is encouraging from a practical point of view that no elaborate 3d textures are needed in realizing remarkable effects. In Figs. 1 c and d, we have plotted the particle trajectories corresponding to a texture like in b. We have used a Gaussian envelope although the precise form is not crucial for the qualitative features. The motion in the x − z plane looks qualitatively similar to the 3d-lens motion but only in that plane -not for any arbitrary plane parallel to the z-axis as in a 3d lens.