Horizon physics of quasi-one-dimensional tilted Weyl cones on a lattice

To simulate the dynamics of massless Dirac fermions in curved spacetimes with one, two, and three spatial dimensions we construct tight-binding Hamiltonians with spatially varying hoppings. These models represent tilted Weyl semimetals where the tilting varies with position, in a manner similar to the light cones near the horizon of a black hole. We illustrate the gravitational analogies in these models by numerically evaluating the propagation of wave packets on the lattice and then comparing them to the geodesics of the corresponding curved spacetime. We also show that the motion of electrons in these spatially varying systems can be understood through the conservation of energy and the quasi-conservation of quasimomentum. This picture is confirmed by calculations of the scattering matrix, which indicate an exponential suppression of any noncontinuous change in the quasimomentum. Finally, we show that horizons in the lattice models can be constructed also at finite energies using specially designed tilting profiles.


I. INTRODUCTION
Analogies in physics, between seemingly unrelated systems, can not only lead to advances in the understanding of these systems, but also contribute to the development of new applications. Among the known analogies, a particularly fruitful connection has been made between gravitational physics and condensed matter systems. Such an analogy was first proposed to study the extraordinary consequences of the curvature of spacetime on quantum fields, namely the Hawking and Unruh radiations, which are in practice too weak to be measured [1][2][3][4][5]. This has attracted considerable amount of interest in recent years from various directions and gravitational analogies have been considered in many contexts including electronic, acoustic, optical and even magnetic and superconducting settings . Particularly, implementations using Bose-Einstein condensates have been used to mimic Hawking and Unruh radiations [30][31][32].
Recent attempts have been made to exploit developments in the prediction and synthesis of Weyl semimetals (WSMs), whose low-energy electronic states can be described by the Weyl equation [33]. In particular, it has been shown that a position-dependent tilting of the Weyl cone in a way to create neighboring regions with type-I and type-II WSM can lead to a black hole analogue [34][35][36][37][38][39]. This analogy can be understood by thinking of the Weyl cones in the material as the light cones of a curved spacetime, thus the boundary of the WSMs of different types as the horizon. Exactly at the transition point between type-I and type-II, the Weyl cones are tan-gent to the zero-energy surface, a situation called type-III WSM. One study has shown that the band structure of Zn 2 In 2 S 5 has such a property that makes it promising for the realization of a black hole horizon [38]. Proposals for tilting the cone as a function of real space include the use of structural distortions, spin textures, and external position-dependent driving [36,[40][41][42][43][44]. Most works until recently assumed that, given that the tilting changes smoothly, one can define band dispersions varying as a function of space. Since WSMs are defined on a lattice, and the tilting variation implies the lack of translational symmetry, and therefore of a reciprocal space, this kind of assumption is imprecise. This was addressed subsequently, by studying a class of tight-binding models with position-dependent nearest-neighbor hopping concentrating on the single-band one-dimensional (1D) cases [45,46].
Here we go beyond previous works by introducing twoband tilted Weyl cones defined on lattices, in one, two and three spatial dimensions and investigate their properties. Taking into account the lattice explicitly is motivated by physical Weyl semimetals having band structures that are ultimately defined on discrete lattices, which play an important role near horizons when wavepackets start slowing down exponentially [45]. Moreover, lattice formulations allow for numerically tracking wavepacket dynamics with the very efficient Chebyshev expansion method [47,48]. We discuss the consequences of different ways of creating Weyl cones on a lattice as such models are not unique. We connect the low momentum limit of these models to gravitational physics and show they are equivalent to Dirac equations in a curved spacetime background. We emphasize that, due to the two-band nature of our models, the connection to the relativistic Dirac equation is more explicit than previously studied incarnations. Then, we explore the propagation of wave packets in these models and compare them to geodesics of gravitational systems. Finally, we calculate the scattering matrix in the 1D models and discuss the large quantitative and qualitative differences arising between the different types of models.

II. TILTED WEYL CONE MODELS
In order to simulate horizons in lattice models we construct inhomogeneous tight-binding models hosting tilted Weyl cones, whose tilting depends on the position in real space. The horizon appears at the points where the tilted cone goes through the Fermi level, i.e. where the WSM goes from type-I to type-II. We want the local low energy effective Hamiltonian to correspond to the following tilted Weyl continuum Hamiltonian where k is the quasi-momentum, σ i are the Pauli matrices acting on a pseudospin representing different orbitals, and t is the tilting vector. The energy dispersion of this Hamiltonian is For t = 0 we get an untilted Weyl cone. In the onedimensional case (see Fig. 1), for t > 0 the spectrum is tilted clockwise and at t = 1 one of the branches of the cone becomes completely flat and the group velocity in this branch becomes zero. In this section, first we discuss how to construct simple lattice models in one-(1D), two-(2D) and threedimensions (3D) that host one or several tilted Weyl cones. Then, we show how these Hamiltonians are connected to the Dirac equation in a curved space-time when the tilting changes as a function of position. With this we get lattice models that describe the motion of electrons in a curved background. as in Eq. (5b) for three different tilting parameters t. The three tiltings correspond to the untilted type I node, the node at the horizon and the overtilted type II node.

A. 1D lattice models
We start with presenting 1D models that have tilted Weyl cones at low energies. In principle, there are many ways to put the Eq. (1) Hamiltonian on a lattice with the same low energy effective Hamiltonian around k = 0, but with different behavior at |k| 0. Here, we study the following two models where c x = (c x↑ , c x↓ ) is the annihilation operator of an electron at site x ∈ Z with pseudo-spin ↑ or ↓. For a constant position-independent tilt, the Bloch Hamiltonians of these after Fourier transformation are given as which give the following dispersion relations (see Fig. 2) Both these models converge to Eq. (1) for k → 0 and describe the same kind of tilted Weyl cone, but they differ away from k = 0. In particular, H 1D 1 has a second Weyl cone at the edge of the Brillouin zone, unlike H 1D 2 which only has one cone. In both cases if the tilting parameter is less (more) than one we get type I (type II) Weyl nodes.
In the overtilted type II case there is an electron (hole) pocket forming to the right (left) of the node. The main difference between the two models, that will be relevant in later results, is that the red band in the overtilted region does not cross zero between k = 0 and k = π in H 1D 1 , but it does in H 1D 2 . In the first model the electron and hole pockets extend throughout the whole Brillouin zone and they connect the two Weyl nodes, while in H 1D 2 the pockets are finite and located next to the single node, which is the situation encountered in the band structure of type II Weyl semimetals.

B. 2D and 3D lattice models
Similar to the 1D case there are infinitely many ways to create 2D and 3D lattice models with a tilted Weyl cone around k = 0. Generalizing the two systems studied in 1D we get the following Bloch Hamiltonians in 2D Their dispersion relations are shown in Fig. 3. H 2D 1 has 4 Weyl nodes in the entire Brillouin zone, while H 2D 2 again only has a single Weyl node.
In 3D the Hamiltonian corresponding to H 1D 2 does not exist, because the σ z Pauli matrix was already used. This is consistent with the Nielsen-Ninomiya theorem [49] which does not allow a single Weyl node in a 3D lattice. The 3D equivalent of H 1D which has 8 Weyl nodes. At k z = 0 the dispersion relation is identical to that of H 2D 1 , represented in Fig. 3(a).

C. Dirac equation in curved spacetime
In the long-wavelength limit, the dynamics of the lattice models can be described by a low-energy continuum model. In particular, the effective Hamiltonians for the lattice Hamiltonians in Eqs. (4), (6) and (7) can be written (using Einstein notation) as (see Appendix A for details) where i runs on the spatial dimensions of the system (d To explore this gravitational analogy in a more precise way, we introduce the metric proposed in Ref. [50] where τ denotes the temporal coordinate, t 2 = t i t i and dx 2 = dx i dx i . The massless Dirac equation for this background metric can we written in the form [51] where γ 0 and γ i are the gamma matrices in flat spacetime and by definition γ ab = [γ a , γ b ]/2 for a, b ∈ (0, 1, · · · , d). To link this equation to the Hamiltonian (8), we take the following representations for the gamma matrices: in (1+1)D: in (1+2)D: while in (1+3)D one may choose the Weyl representation The Dirac equation on the gravitational background (9) coincides with the low-energy dynamical equation for the lattice models provided the last term in Eq. (10) vanishes. In (1+1)D, this is always the case. In higher dimensions, we constrain the vector t i (x) to have a vanishing curl which cancels the last term in (10). Keeping this term would produce a spatially-varying on-site potential in the tight-binding model. Since the term is proportional to the derivative of the tilting, its effect is negligible in systems with slowly varying tilting. The region in space described by t 2 − 1 = 0 can be thought of as an event horizon. Notice that Eq. (9) represents a wide class of metrics, some of which are of particular importance in the context of gravitational physics (see, e.g., a recent discussion in [52]).
To simulate different metrics we introduce spatial inhomogeneity in the tilting parameter. In the real space Hamiltonians (3) we make the tilting parameter position dependent. This is different from previous models where the position dependence is in the Fermi velocity [45,46,53]. In this paper we will only consider effectively 1D horizons, so t x will be a function of x and t y/z = 0. In all cases the horizon will be defined by the points where t x = 1.

III. WAVE PACKET DYNAMICS
Now that we defined the systems of interest, we move on to study the propagation of wave packets in lattice in Eq. (6b) for three different tilting parameters in the x direction tx, and ty = 0. The three tiltings correspond to the untilted type I node, the node at the horizon and the overtilted type II node. models with horizons. We study the different behaviors for our models with different horizons. We will show how the world lines of wave packets in our models match the geodesics of specific metrics, and we will simulate the propagation of Gaussian wave packets based on the Chebyshev expansion method. The details of this method are explained in Appendix B. The videos for all simulations discussed in the paper are included in the Supplemental Material (SM).

A. Linear horizon in 1D
First, we discuss 1D models with linearly changing tiltings t x = 2x/L, where L is the length of the system. With this we linearly increase the tilting from t = 0 to t = 2 with the t = 1 horizon exactly at L/2. The systems we study are L = 1000 long. The starting wave packet is localized at x 0 = 100 in the red band of Fig. 2, and propagates to the right towards the horizon.
The red band becomes more and more flat approaching the horizon in both H 1D 1 and H 1D 2 . As a result, the group velocity of the wave packet gets smaller and smaller. While the wave packet is slowing down it becomes narrower. The time dependence of the position of the wave packet is shown in Fig. 4. For the figures we rescaled the position and time using ξ = 1 − 2x/L and τ = 2t/L. With this the wave packet moves from the initial position ξ 0 = 0.8 to the horizon ξ = 0.
The wave packet dynamics results can be understood with a philosophy similar to that of the WKB approximation. Since the tilting parameter varies continuously and slowly with position we can think of the system as being locally the infinite H 1D 1/2 (k). Since the tilting is proportional to σ 0 the eigenvectors at different tiltings are all the same, from which it follows that states in the red or blue bands will stay eigenstates at whatever point of the chain. This means, that a wave packet that is initially in the red (blue) band will stay in the red (blue) band. The system is not strictly homogeneous, thus the momentum is not a conserved quantity, but since the inhomogeneity comes from a slow change in space, large jumps of the momentum are exponentially suppressed (this will be further clarified in Sec. IV). The quantity that is strictly conserved is the energy of the wave packet. Using the above mentioned rules together with the dispersion relations in Fig. 2, the results of the 1D simulations can be explained.
First, we discuss the H 1D 1 system. Focusing on the red band we see that it becomes completely flat at the horizon. If we take a packet with a small energy E > 0, the horizon is invisible for it since after a certain tilting there is no point in the red band that is at this energy. With a continuous slow change of the momentum, the state with energy E will go from k = 0 to k = π and the group velocity will go from positive to negative while the wave packet bounces back from the horizon (see Fig. 5(a)). The closer the energy is to 0, the closer the wave packet reaches the horizon and in the limiting case of E = 0 the wave packet gets pinned to the horizon. On longer time scales because of the lattice length-scale that unavoidably comes into play [45], the wave packet slowly disintegrates and the center of mass will go towards the left. The behavior of the wave packets at different energies can be seen in Fig. 4 with the green curves.
The H 1D 2 system is qualitatively very different from the H 1D 1 system. Here, at every position there is an E > 0 state in the red band. This means that a finite energy wave packet will go through the horizon with a continuous change in momentum (see Fig. 5(b)). The zero energy wave packet will approach slowly the horizon because of the zero group velocity at the horizon but eventually it will go through it. The behavior of the wave packets at different energies can be seen in Fig. 4 with the purple curves.
The time dependence of the wave packet approaching the horizon is consistent with the results in Ref. [54]. There, a one-band model with nearest-neighbour hopping was studied, where the hopping varies with the position. This can be understood, by coarse-graining, as a model with a linear dispersion close to zero energy, whose slope varies with position. The slope was set to zero at the origin of the lattice, and finite away from it, thus mimicking the evolution of a light cone when moving away from a horizon. Since it has only one band, the hopping is exactly zero where the horizon is meant to be, meaning that the right and left side of the lattice are completely disconnected. In our two-band model this is not the case and, as we saw in the H 1D 2 system, the wave packet can propagate through the horizon.
It was shown in Ref. [54], both using numerical calculations and an analytical derivation based on a semiclassical approximation, that zero-energy wave packets propagating in the one-band model precisely follow the geodesics of a (1+1)D dilaton gravity. In the case where the hopping evolves linearly with position, these geodesics can be expressed as where in our case α = 2/L. Using the previously defined scaled position and time this simply reads as As we can see in Fig. 4 this dependence is present for our models too, with deviations only at larger timescales due to finite size effects, when the size of the wave packet becomes comparable to the lattice constant.

B. Linear horizon in 2D and 3D
In 2D and 3D, wave packets propagate qualitatively differently from the 1D case. To understand the difference, we first consider a 2D system where t x = 0 and t y = 0 everywhere. A Gaussian wave packet centered Figure 6. Propagation of wave packets in the 2D systems without tilt. The figure shows snapshots of the wave packet for different starting momenta in the system H 2D 1 . The length and width of the system are L = 500 and W = 500, the initial wave packet has a width of 40 sites. The anisotropy in the top row is due to the specific pseudospin configuration of the initial wave packet.
at a large enough k will propagate based on the group velocity at k and will maintain its Gaussian shape. But getting closer to k = 0 this is no longer true. A wave packet with finite size in real space will also have a finite size in momentum space. This means that at low enough momenta the wave packet encapsulates the Weyl node and thus it will have components propagating in all directions radially. This behavior of the wave packet propagation is shown in Fig. 6 with snapshots (in the SM the full animations are available).
In all simulations we choose the pseudospin components corresponding to the eigenvector of the infinite Hamiltonian with the average momentum of the wave packet. Because the eigenvectors are momentum dependent and because the wave packet includes multiple momenta a zero momentum wave packet includes a superposition of positive and negative energy states which leads to Zitterbewegung of the electrons [55][56][57][58].
The world lines of the wave packets in 2D and 3D are shown in Figs. 7 and 8. Here we only focus on the center of mass of the wave packets on the x axis. As we can see the results are very similar to the 1D results in Fig. 4. The main difference are the oscillations due to the Zitterbewegung at small times.

C. Power law horizons in 1D
So far the tilting was always linearly dependent on the position, now we consider more generic cases and study how the different choices affect the wave packet propagation. We only discuss the H 1D 1 system, for the H 1D 2 model similar results can be obtained.
We consider models that follow a power law close to the horizon where t = 1. For the position dependence of the tilting parameter we use where γ is the exponent. At γ = 1 we recover the model studied in the previous sections. This tilting dependence is shown in Fig. 9(a). The world lines of the wave packet for different γ pa- rameters are shown in Fig. 9(b). For linear and supralinear dependencies (γ ≥ 1) we get wave packets that infinitely approach the horizon. The higher the exponent the slower the approach is. For sublinear dependencies (γ < 1) the wave packet goes very close to the horizon but ultimately bounces back from it. These results are consistent with the results obtained for the single band model in Ref. [54].

D. Horizons for finite energy in 1D
In Sections III A we saw that in the H 1D 1 model, for the linear profile a wave packet at finite energies always deviates from the exponential time dependence and bounces back from the horizon. This happens because the effective horizon at finite energies is no longer at t = 1 but at lower values, thus the wave packet can not reach the center of the chain. Focusing on the red band in Fig. 2(a) the dispersion relation and group velocity are which implies v 2 = (1 − t) 2 − E 2 . Hence, for a wave packet of energy E this means that there is an effective horizon at Knowing this we can construct a tilting profile that has a horizon at finite energy similar to the zero energy one. In order to reproduce the zero energy world line we want a linearly decreasing group velocity For 2x ≤ L this can be satisfied by choosing where k 0 is the initial momentum and x 0 is the initial position of the wave packet. For 2x > L we can choose an arbitrary tilting profile since that part of the chain will not be accessible by the wave packet. For simplicity we use the same dependency but flipped. The tilting parameter as a function of position is shown in Fig. 9(c). Using these horizons the wave packet world lines are shown in Fig. 9(d). As we can see the world lines that were bouncing back with the linear profile at finite momenta now are approaching the horizon similarly to the k = 0 wave packet. Close to the horizon because of lattice effects we get the disintegration of the wave packet similarly to the k = 0 case.

IV. SCATTERING MATRIX
To characterize the different possible scenarios in the models discussed above we turn to scattering theory. Taking the same systems studied in the previous section we attach two semi-infinite leads on both ends. We restrict the discussion to 1D systems, because with the change in tilting always being in one direction only, we can always choose periodic boundary conditions in the other directions to get effectively 1D systems. The leads are made using the same model as the system with fixed tilting t = 0 (t = 2) for the left (right) lead. We then calculate the scattering matrix between the two leads numerically using the Kwant package [59]. The scattering matrix (S) encodes the probability amplitudes of scattering from incoming modes (i) into outgoing modes (o).
where ε is the energy of the considered modes. First, we discuss the H 1D 1 system. The propagating modes at a given ε > 0 energy in this system are shown in Fig. 10a. There are two incoming and two outgoing modes in each side of the scattering region resulting in a scattering matrix that is 4 × 4. The scattering matrix at ε = 0 is shown in Fig. 10c. We can see that each incoming mode is scattered into an outgoing mode with probability one. Incoming modes from the red bands are totally reflected into outgoing modes in the red band of the same lead. Incoming modes of the blue bands are completely transmitted into outgoing modes in the blue band of the opposite lead. This behavior is consistent with the reasoning of the previous section: the red band becomes completely flat so states in this band cannot cross the horizon, while the blue band stays qualitatively the same and the states in this band can go through the horizon. The same is valid for all energies as long as there are propagating modes in both leads.
Then, we study the H 1D 2 system. The propagating modes at a given ε > 0 energy in this system are shown in Fig. 10b. The left lead has one incoming and one outgoing mode while the right lead has two and two resulting in a scattering matrix that is 3 × 3. The scattering matrix at ε = 0 is shown in Fig. 10d. The incoming mode in the blue band is totally transmitted into the outgoing mode in the blue band. At zero energy the incoming modes in the red band are equally split between the two outgoing modes in the red band. In this system there is a big difference between the red incoming mode on the left and right sides. Starting from the left side we get total transmission, while starting from the right side we get total reflection.
In the H 1D 2 system the energy dependence of the scattering amplitudes is more complicated than the H 1D 1 system. The zeros and ones are unaffected but the equal splitting between the red bands is only valid at ε = 0. Let us consider the i l incoming mode. At a large enough positive energy if we follow the same energy modes throughout the scattering region we see that they appear at continuously increasing momenta until we reach a large enough tilting at the other side of the horizon where a second same energy mode appears at a very different negative momentum. Since large momentum changes are not allowed because of the slow varying of the tilting the i l mode will be totally transmitted into the o r1 mode at large enough positive energies. The closer we go to zero energy we see that at the point of the appearance of the second mode with the same energy, the difference between the two momenta is decreased, and at ε = 0 it is equal to zero. This means that the i l mode can scatter into both the o r1 and o r2 with continuous momentum change, thus they have equal scattering probabilities. In between the zero and the high enough energies the scattering probabilities show an exponential behavior (see Fig. 11) The same was observed for the H 1D 2 Hamiltonian for a similar tilting dependency in Ref. [27,28]. The energy dependence of the scattering rate can be expressed as We see that the momentum jumps are exponentially suppressed. Increasing the system size makes the change of tilting smoother, which in turn makes the quasi conservation of momentum stricter, and causes the probability amplitudes to reach their asymptotic behavior faster. The Hawking fragmentation coined in Ref. [27] can thus be understood as intermode scattering between multiple states at the same energy with different momenta. Increasing the energy of the incoming mode increases the momentum difference between the scattered modes, which leads to the energy dependence shown in Fig. 11.

V. DISCUSSION AND CONCLUSIONS
We introduced lattice models for tilted WSMs where the tilting varies smoothly with position. These models at small momenta effectively correspond to the massless Dirac equation in a curved spacetime background. To understand the role of high momenta and inter-Weylnode effects, we studied two different types of systems: the first type (H 1 ) has 2 d nodes where d is the dimension of the system, while the second type (H 2 ) only has a single Weyl node.
Using these models, we simulated wave packet propagations for three types of tilting profiles in 1D, 2D and 3D, using the Chebyshev expansion method. When the tilting varies linearly with position, and in the zeroenergy limit of both H 1 and H 2 , the wave packets follow the geodesics of (1+1)D de Sitter spacetime, as found previously in a single-band model [45]. At finite energies (E > 0) though, H 1 and H 2 give very different behaviours. H 1 causes wave packets to bounce back, similarly to the single-band model, while H 2 causes wave packets to be transmitted. This is explained by the differences in the band structure at large momenta, which allow for scattering to same-energy states in H 2 but not in H 1 . When the tilting varies as a power law with position, we found that the transition between models with and without a horizon found in the single-band model is also present in two-band models. Finally, we designed a specific tilting profile that moves the horizon away from zero energy, and can be tuned to obtain eternal slowing down for any initial wave packet.
We showed that these results can be understood using the local dispersion relations, the conservation of energy, and only allowing continuous change for the quasimomentum. To further show the validity of this approach, we computed the scattering matrix for 1D systems, which corresponds well to this analysis. H 1 only has fully transmitted or reflected modes, and no scattering takes place. While in H 2 , due to the presence of multiple states at a single energy, intermode scattering does occur, in agreement with scattering amplitudes obtained in a related model [27], which we now relate directly to wave packet trajectories.
In summary, our results bridge the gap between the single-band description and concrete proposals for the experimental realisation of Weyl gravitational analogues. Indeed, they point to the rich physics arising from the presence of a second band, and the presence or absence of other Weyl nodes at large momenta. These can in some cases dominate the observed phenomena and therefore indicate limitations of the analogy between tilted WSMs and gravitational systems which should be taken into account in experimental setups. Finally, we identified the possibility of tuning the tilting profile to provide horizons for finite-energy wave packets, thus expanding the possibilities offered by this type of gravitational analogue. This algorithm does not require the eigenvalue problem of the Hamiltonian to be solved. Using sparse matrices the recurrence relation can be evaluated efficiently, allowing us to study much larger systems than those accessible with methods that require the diagonalization of the Hamiltonian.