Quantum dynamics in 1D lattice models with synthetic horizons

We investigate the wave packet dynamics and eigenstate localization in recently proposed generalized lattice models whose low-energy dynamics mimics a quantum field theory in (1+1)D curved spacetime with the aim of creating systems analogous to black holes. We identify a critical slowdown of zero-energy wave packets in a family of 1D tight-binding models with power-law variation of the hopping parameter, indicating the presence of a horizon. Remarkably, wave packets with non-zero energies bounce back and reverse direction before reaching the horizon. We additionally observe a power-law localization of all eigenstates, each bordering a region of exponential suppression. These forbidden regions dictate the closest possible approach to the horizon of states with any given energy. These numerical findings are supported by a semiclassical description of the wave packet trajectories, which are shown to coincide with the geodesics expected for the effective metric emerging from the considered lattice models in the continuum limit.


I. INTRODUCTION
Connections between different fields of physics have proven fruitful by opening entirely new research avenues in recent years. Dualities between gravitational and many-body theories have become important tools in the study of quantum critical systems [1], while the search for electronic states following the Dirac equation, prevalent in high-energy physics, led to the discovery of topological insulators, one of the most major fields in condensed matter physics in decades. In this context, replicating some of the physics of curved spacetime in condensed matter settings has been proposed as a promising way of on the one hand understanding gravitational problems in a simplified setting, and on the other of searching for novel, gravity-inspired, physical effects within condensed matter theory [2]. In particular, in a seminal work Unruh proposed the construction of an analog black hole horizon and its radiation using a fluid flowing with a spatially varying speed profile that is partly subsonic and partly supersonic [3,4]. Similarly, many proposals for analog gravity setups emulating a broad range of emergent curved spacetimes have been put forward in a variety of electronic, acoustic, optical and even magnetic and superconducting settings [2,. Some of these proposals have already been implemented in experiments, mainly using Bose-Einstein condensates [27][28][29]. In all of these proposed and realised black hole analogues however, the role played by the atomic lattice (periodic or otherwise) remains largely unexplored, even though it is an essential component of any condensed matter system.
Recently, it was shown that a black hole analogue may be realized in Weyl semimetals (WSMs) by tilting the Weyl cone as a function of real space, transitioning from a type-I to a type-II WSMs [30][31][32][33][34][35]. The tilt causes part of the band structure close to the Weyl node to become progressively flatter as the type I-type II transition point is approached. This is a direct analogy for the tilting of a light cone close to a black hole, in which case the surface on which the light cone tips across the time axis defines an event horizon, beyond which all light is trapped [30]. It has been shown that Zn 2 In 2 S 5 sits precisely at the transition, and was coined to be a 'type-III' Weyl semimetal [34]. Tuning the tilt of the cone across real space could be achieved using structural distortions, spin textures, or external position-dependent driving [32,[36][37][38][39][40].
To circumvent the difficulty of defining spatial variations in the tilt of Weyl cones, which themselves require reciprocal space and translational symmetry to be defined, previous studies generally assume that the control parameter responsible for tilting changes smoothly, and that a band structure varying as a function of real space can be defined, despite the absence of translation symmetry. Here, we take a more rigorous approach, and define a Hamiltonian in real space with a hopping that varies as a function of position. This has the benefit of allowing us to explicitly study new effects that arise from the presence of the lattice and the potential limitations it imposes on the dynamics in the emergent analogue gravitational system.
We consider generalized nearest-neighbor tightbinding (TB) models with position-dependent hopping and show that their low-energy dynamics is similar to that of a Dirac field with a position-dependent velocity, which mimics the presence of a background curved spacetime [41,42]. Rather than the Dirac equation and Weyl nodes in 3D we consider single band systems in 1D with progressive band-flattening in real space which are fully tractable numerically, allowing direct comparison with analytical semi-classical treatments of the problem. Then, concentrating on the case of power-law position-dependent TB models, we demonstrate how and when the model serves as an analogue for horizon physics, as witnessed by the critical slowing down of wave packet dynamics. We support the numerical results with a semiclassical picture valid at large wave lengths, but for arbitrary spatial variations of the hopping integrals. We find the formal solution for semiclassical trajectories in the most general case, accompanied by explicit solutions for the power-law dependencies to compare with the numerical results. Finally, we show that all the eigenstates in the models with power-law variation of the hopping are also localized in a power-law manner. Consistent with the observed wave packet dynamics, the low-lying states localize on the horizon, whereas high-energy states are exponentially suppressed in a region near the horizon.
In the following, we first introduce the general model and its low-energy sector with gravitational analogies (Sec. II). Then, we calculate numerically the time evolution of wave packets for power-law hopping models in Sec. III. We derive semiclassical equations of motion along with their solutions and compare them to numerical results in Sec. IV. This is followed by a discussion of the eigenstates in Sec. V, before we conclude in Sec. VI.

II. TB MODELS AND THEIR GRAVITATIONAL ANALOGY
Consider electrons on a one-dimensional lattice of N sites, with nearest-neighbour hopping only: Here, t n is a position-dependent hopping parameter whose amplitude increases with n. This Hamiltonian has a particle-hole symmetry (PHS) represented by the unitary operatorÛ The PHS can be broken by adding a potential energy, but we concentrate on models with PHS throughout this work.

A. Low-energy limit
The low-energy properties of these lattice models mimic those of a Dirac field in curved spacetime. This is made plausible by the observation that approximating the full lattice by disconnected (periodic) sections, results in local band structures ε(n, k) ∼ −2t n cos k [43]. Each of these has two Fermi points k F = ± π 2 at halffilling. Accordingly, we can define the local Fermi velocity v F (n) = ∂ k ε(n, k = ±k F ) ∼ ±2t n . Motivated by this observation, we construct a precise correspondence of the low-energy properties of the lattice model to a Dirac field with position-dependent velocity by introducing the transformation a n = ν=±ψ ν (x n ) e iνk F n . (3) We additionally take the continuum limit, in which The second term includes fast oscillations e 2iνk F x , and can be neglected in the limit of slowly-varying fieldŝ ψ ν (x). Therefore, the lattice model is equivalent in the continuum limit to a model for the 1D massless Dirac fieldΨ = (ψ + ,ψ − ) T governed bŷ with a space-dependent velocity such that v(x) = 2t(x), and Pauli matrix σ z acting on the spinor. By variation we obtain the equation of motion which is identical to that describing the dynamics of a 1D massless Dirac field in the presence of the background (1 + 1)D metric This metric can possess a horizon at positions where v(x) = 0. Based on this property, we mainly consider power-law variations of the hopping integrals, which yields effective local velocities of the form v(x) = v 0 x γ . In the following section we will study the wave packet dynamics on lattices with power-law hopping variation and discuss the results in light of their gravitational resemblances.

III. WAVE PACKET DYNAMICS
Motivated by the equivalence between the low-energy dynamics of the position-dependent lattice model and Dirac particles in a curved space, in this part we explore numerically the wave packet dynamics on these lattices. We focus on the cases where the hopping grows as a power-law with position: The initial momenta are p0 = −π/2 (top) and p0 = −0.9 × π/2 (bottom). In the first case, the wave packet slows down and localizes at the origin of the lattice, where it disintegrates. In the second case, it is reflected before reaching the origin.
such that the maximum hopping in the system is equal to one. A key aspect of black hole physics is that an observer at infinity will see a wave packet falling towards a black hole become sharper and slower as it approaches the horizon, before asymptotically reaching the horizon shaped as a Dirac distribution. As we expect the velocity of a wave packet in the lattice model to be associated with the local strength of the hopping, the equivalent of a horizon in the lattice model may occur where the Fermi velocity approaches zero, i.e. in the vicinity of the bond of the lattice between site n = 1 and a 'virtual' site n = 0. To compute the time evolution of an initial wave packet in the lattice model we work in the basis of the diagonalized Hamiltonian and use the expression where | signifies the th eigenvector of H, and τ denotes time. We define a Gaussian wave packet by: where w is the width, n 0 the initial position and p 0 the initial velocity of the wave packet. Note that, for the sake of simplicity, we introduce the rescaled parameters n = n/(N − 1),τ = τ /(N − 1), andw = w/(N − 1). Let us first consider a linearly increasing hopping parameter (γ = 1). In that case, we find two different pos- sible types of behavior for the wave packet, depending on whether p 0 is equal to or different from −π/2. Example time evolutions of both cases are presented in Fig. 1. In each case, the wave packet starts by sharpening and slowing down as it moves towards n = 0. Wave packets with p 0 = −π/2 never reach the origin of the lattice, and instead come to a standstill at non-zero n before moving away from the origin and broadening again. In contrast, the peak position of wave packets with p 0 = −π/2 continues to approach the origin of the lattice indefinitely. As their peak comes close to n = 0, these wave packets start to form ripples in their tails, which move away from the origin. Eventually, the wave packet consists almost entirely of these ripples, but conserves a maximum amplitude at the origin of the lattice.
The observed asymptotic localization of wave packets at the origin coincides with the key feature expected for wave packet dynamics in the presence of a horizon. One key difference with what is expected close to a black hole horizon, however, is the formation of ripples. This feature of the model can be understood as a consequence of the discreteness of the lattice and the unitarity of time evolution. Indeed, consider two different wave packets both with p 0 = −π/2. If they could asymptotically localise at the origin of the lattice, they would become indistinguishable from each other, and it would then be impossible to propagate them back to their distinct original configurations by reversing time. Since this cannot be the case in our system, which has unitary time evolution, the two wave packets have to develop specific features close to the origin of the lattice, which act as signatures of unitary time evolution. These features here come in the shape of ripples, propagating away from the origin of the lattice. Wave packets with p 0 = −π/2 do not exhibit eternal slowdown, pointing to the fact that the horizon physics can only be probed at a critical initial momentum.
We now turn to the case of γ = 1/2, describing a square-root position-dependence of the hopping. The time evolution of the wave packet amplitude for two values of p 0 is displayed in Fig. 2. In this case neither of the wave packets localizes at the zero-velocity point, and both turn around and move out to infinity at late times. Additionally, we don't observe any formation of ripples, unlike in the case γ = 1. One similarity with this case, however, is that wave packets with p 0 = −π/2 reach the origin of the lattice, while wave packets with p 0 = −π/2 reverse direction at a nonzero distance from the origin. All these features suggest that while the model with γ = 1 gives rise to a horizon, the model with γ = 1/2 does not. In the next section we will combine numerical results with a semiclassical analysis and show that indeed γ = 1 is a critical value below which the model does not contain any horizon physics. In contrast, for γ ≥ 1, the point n = 0 resembles the horizon of a black hole with low-energy particles eternally slowing down upon approaching it.

IV. SEMICLASSICAL DYNAMICS
To further analyze the time evolution of wave packets, it is constructive to compare the numerical results to a semiclassical description for the trajectories of the wave packet center of mass. We introduce a continuous functionψ(x) which coincides with the wave function of the lattice model at the discrete lattice points, so that ψ(x n ) = ψ n with x n = n/(N − 1).
For a general position-dependent hopping model, the equation for the energy eigenvalue can be written as the recursive relation Assuming that we can expandψ(x) in a Taylor series, we can relate ψ n±1 toψ(x n ) exactly, by where δx = 1/(N − 1), and we replaced derivatives usingp m = (−id/dx) m . We also summed over odd and even powers separately, which can be done formally by expressing them sum in terms of sines and cosines of the momentum operator. The result is nothing but the well-known expression of the translation operator T ∆x = e −i∆xp . Now assuming t n ≡ t(x n ) in Eq. (11) with a well-behaved function t(x), the eigenvalue equation can be re-written as The right-hand side in Eq. (13) can be interpreted as the continuum Hamiltoniañ where we used the fact with t (x) = dt/dx. Now, neglecting the commutation relations betweenx andp, we obtain semiclassical EOM for the expectation values x and p: Here, we rescaled time and momentum as τ → τ /δx and p → p/δx. Differentiating Eq. (17) and replacing the derivative of p by the right-hand side of Eq. (18) yields which is a straightforward second order differential equation for the dynamics of the position.

A. General solutions for the trajectories
In this part, we present the formal solution to the semiclassical Eqs. (17) and (18) for general form of the hopping. We first note that defining an auxiliary function This allows us to write Eq. (19) in the simplified form Therefore, F 2 − 2 t(x) 2 is just a constant A and replacing F with its original definition, we end up with the first order differential equation with a formal solution for the most general case. The integration constant A can be fixed using Eqs. (17) and (22) at τ = 0. This yields A = − 2 t(x 0 ) cos p 0 2 , with x 0 and p 0 indicating the position and momentum at τ = 0. This also fixes the signs in Eqs. (22) and (23) to be −sgn[t(x) sin p].
Notice that at a turning point p = 0 when momentum undergoes a sign change, and at points where the sign of the hopping parameter switches, the sign in Eqs. (22) and (23) (17) and (18), we obtain the new equation which directly relates the position and the momentum. It has the solution This relation shows that t(x) cos p is a constant of motion and, in fact, we can assign E wp = −2t(x) cos p as the conserved average energy of the wave packet in a semiclassical sense. In particular, we see that for initial value p 0 = ±π/2, the momentum of the wave packet remains constant throughout the time evolution. Since it also implies E wp = 0 for all times, this can also be thought of as a consequence of energy conservation. Finally, Eq. (25) determines the position of the turning point of the wave packet (when p = 0) and in particular the minimum distance from the horizon, as t(x min ) = t(x 0 ) cos p 0 .

B. Trajectories for power-law hopping
Although Eqs. (23) and (25) give a general solution for the semiclassical equations, the former is just a formal expression in terms of an integral. Here, we therefore focus on the specific case of power-law hopping, defined as t(x) = x γ , for which the trajectories read using the hypergeometric function 2 F 1 (a, b, c; z) with three real parameters a, b, c, and the variable z. This expression gives a real value only when −4x 2γ /A > 1 or equivalently x > x 0 (cos p 0 ) 1/γ , in agreement with the turning point being given by x min = x 0 (cos p 0 ) 1/γ . For the special cases of γ = 1 or γ = 1/2, corresponding respectively to linear and square-root forms of position-dependence, the solution of Eq. (26) simplifies to Inverting this result and writing the position x in terms of the time τ , while also choosing values for B by matching different parts of the solution, yields for linear and square-root position-dependence respectively. Substituting these back into Eq. (25), the evolutions of the corresponding momenta are found to be cos p γ=1/2 = √ x 0 cos p 0 C. Zero-energy wave packets and equivalence to geodesics In the limit of p 0 = −π/2, the constant A vanishes, and the general spatial trajectory of Eq. (22) becomes Not surprisingly, the semiclassical trajectories in this limit coincide with the lightlike geodesics (ds 2 = 0) of the general metric in Eq. (7). In the case of power-law hopping variations, the integral equation (23) simplifies to which in turn leads to Notice that these expressions agree with Eqs. (27) and (30) after substituting p 0 = −π/2. two distinct types of behavior, depending on the value of the exponent γ. If γ ≥ 1, the wave packet faces an eternal deceleration and only asymptotically reaches the horizon. In contrast, if γ < 1, the wave packet reaches the point x = 0 in a finite time. Although its velocity vanishes there, it does not become stuck. Instead, we observe back-scattering from that point similar to a classical particle bouncing from a hard wall. The value γ = 1 thus separates a region of parameter space where a wave packet with momentum p 0 = − π 2 reflects off the origin of the lattice from one where it localizes at the origin. In general relativity, the eternal slowdown of particles is a key feature of particles approaching a black hole horizon as seen by a distant observer. Therefore, the transition at γ = 1 found here separates lattice models with and without a synthetic horizon. The special case of γ being precisely one has previously been shown to mimic a (1 + 1)D anti-de Sitter spacetime [42].

D. Comparing semiclassical and numerical results
In order to see the effect on the trajectories of changing the initial momentum p 0 , we show in Fig. 4 the time evolution of the wave packet maximum obtained numerically for γ = 1 and γ = 1/2. The numerical results are in good agreement with the semiclassical trajectories given by Eqs. (27) and (28). For the special case of p 0 = −π/2 and γ = 1, the wave packet asymptotically approaches the horizon at x = 0 at large times τ . For all other momenta, the wave packet bounces back at a nonzero distance x min = x 0 cos p 0 from the horizon (see also Fig. 5). For γ = 1/2, the wave packet never local- We now turn to the deviations from the semiclassical picture to quantify its breakdown, as signalled by the disintegration of the wave packets stuck to the horizon in Fig. 1. We define a wave packet ψ G (τ ) whose position x(τ ) and momentum p(τ ) are given by Eqs. (27) and (18), and whose width follows the same time dependence as x(t). Notice that this wave packet is not a solution to the dynamics, but serves as a reference or idealized case to compare the actual dynamics of a wave packet on the discrete lattice to, for the case of linear hopping. The numerical overlap between these two wave packets is shown as a function of time in Fig. 6(a). It is practically constant and equal to one up to the point where ripples start to appear in the time-evolved wave function, at which point the overlap starts to decrease. The large oscillations seen at late times can be explained by the evolution of ψ G alone: at large t, the width of ψ G is smaller than the lattice size, and therefore the norm of ψ G evaluated on the discrete lattice oscillates, depending on whether its peak position is on a lattice point or between sites.
To further analyse the relation between the non-zero lattice spacing and the decreasing overlap and formation of ripples, Fig. 6 shows the overlap between ψ G (t) and ψ(t) for multiple values of the width. The onset time for the decrease of the overlap goes up with increasing width of the initial wave packet. For initial widths of 20, 30, 40, and 50 lattice spacings, the width of ψ G when the overlap reaches 1/2 is 1.31, 0.96, 0.86, and 0.71 lattice spacings respectively. These values are close to the lattice spacing which therefore acts as an effective critical value of the width. We therefore argue that the lattice plays a key role in the formation of the ripples, which are the main observable difference between the dynamics of wave packets with p 0 = −π/2 in the lattice model and those in relativistic continuum theories.

V. EIGENSTATES AND THEIR LOCALIZATION
So far, we studied the dynamics of wave packets in position-dependent lattice models and compared them to a semiclassical picture to highlight possible gravitational analogies. It has also been previously shown that the eternal slowdown of zero-energy wave packets upon approaching the horizon is always associated with the presence of a divergent density of states (DOS) at zero energy, in the N → ∞ limit [42,44]. In particular, it was pointed out that the transition from back-scattering to slowdown at γ = 1 coincides with a transition in the spectral properties of the lattice Hamiltonian at zero energy. In this part, we therefore turn to the eigenstates of the position-dependent lattice models whose properties can shed further light on the features found in the wave packet dynamics as well as the density of states.
The Hamiltonian of Eq. (1) for power law hopping can be written as an N × N matrix with non-vanishing elements Then the eigenvalue problem H|Ψ ε = ε|Ψ ε with |Ψ ε = (ψ 1 , · · · , ψ N ) T can be written as a set of N coupled equations in the recursive form, (34) with the boundary conditions ψ 0 = ψ N = 0.

A. Exact form for zero-energy states
For the special case of ε = 0 (zero modes), we can find the exact form of the discrete wavefunction amplitudes ψ n . Eq. (34) for zero-energy states becomes (n − 1) γ ψ n−1 + n γ ψ n+1 = 0, which yields Repeating the sequence above, we obtain the wave function amplitudes on even and odd sites Equation (35) for n = 1 readily shows that ψ 2 = 0 for zero-energy states, and, subsequently, from Eq. (38) we find that the wave function amplitude identically vanishes on all even sites.
For an even number of lattice points N e = 2N , the boundary condition at the second end will read ψ Ne+1 = ψ 2N +1 = 0 which cannot be fulfilled unless all odd lattice points have vanishing amplitudes. As a result, for an even number of lattice points there is no zero mode at all. For an odd number of lattice points on the other hand, the boundary conditions on both ends of the chain force the amplitude to vanish on even lattice points. Therefore for odd N o = 2N + 1 we find a single zero mode with the wave function given by Eq. (37).
The qualitative behavior of this wave function in the limit of large n, can be found using the Stirling's approximation formula n! ≈ √ 2πn(n/e) n , with Euler's constant e, to be This form applies to the tail of the zero-energy wave function in a large lattice (N 1).
Finally approximating ω n − ω n−1 ≈ ω n+1 − ω n ≈ dω/dx, with ω(x) a continuous function such that ω(x = n) ≡ ω n , we obtain the differential equatioñ or equivalently The general solution of this equation (for γ = 1) reads For the special cases γ = 1 and γ = 1/2, this gives the respective approximate eigenstates respectively. It should be noted that although they approximate the exact eigenstates, the approximate wave functions ψ n are not necessarily normalized or even orthogonal to each other. In addition, there is no limitation on the energy ε, except for the fact that we only get mathematically well-behaved results for an energy range consistent with the exact bandwidth. Barring these unavoidable shortcomings and the appearance of some phase differences between ψ n and the exact eigenstates for low n, we find good agreement between the eigenstates obtained numerically by exact diagonalization and the real part of the approximate analytical results of Eqs. (45) and (46), as shown in Figs. 7 and 8. Both for γ = 1 and γ = 1/2, two qualitative features of the wave functions stand out: (i) a power-law localization and (ii) regions of suppressed amplitude. States with energies close to zero have an envelope approaching its maximum at the vicinity of n = 0, where the hopping becomes vanishing small. In contrast, the envelope of eigenstates with nonzero energy is suppressed and becomes vanishing small for a finite range around n = 0. The extent of the suppressed regions grows with the absolute value of energy, |ε |, and is bordered by a region with power-law behavior for the envelope of the wave function.
The appearance of forbidden regions is a universal feature for all models with power-law variation of the hopping studied here. This can be understood by noticing that in both Eqs. (42) and (43) the function ω(x) (or its discrete counterpart ω n ) acquires an imaginary part for n < n c,ε ∼ ε/2 1/γ . As a result, the real part of the trial wave functions in Eq. (40) show an exponentially decaying position-dependence for These regions appear shaded in the right panels of Figs. 7 and 8. A more intuitive picture for the existence of forbidden regions can be found by recalling the local band structure picture, with ε(n, k) ∼ −2 n/(N − 1) γ cos k. This shows that for a state with non-zero energy ε, the region with 2 n/(N − 1) γ < |ε| becomes classically forbidden, as there is no available locally extended states there. The only way to penetrate that region is then by quantum tunneling, with its associated exponential decay of wave function amplitudes. Moreover, the existence of forbidden regions for states with non-zero energy provides an alternative explanation for the back-scattering of wave packets with p 0 = −π/2, which have non-zero average energy.
C. Exact solution of the lattice model with γ = 1/2 The lattice model with t n = √ n, i.e. γ = 1/2, is of special interest as one can establish an exact expression for the eigenstates in terms of Hermite polynomials, as we will show in the following. In this case, the eigenvalue equation reads with the boundary conditions ψ 0 = ψ N = 0. Below, we show that this recurrence relation can be related to the relation for the well-known Hermite polynomials To show this, let us change to variablesψ n , defined as ψ n+1 =ψ n / √ 2 n n!. The recurrence relation in terms of the new variables then becomes −εψ n−1 which simplifies to Comparison of Eqs. (49) and (52) then yields with the normalization constant as derived in detail in appendix A. The boundary condition ψ N +1 = H N (ε √ 2 ) = 0 implies that the eigenvalues ε m are directly related to the roots of N -th Hermite polynomial. Since H N is guaranteed to have N real roots, this always yields all N eigenvalues for the problem. In appendix A we provide the asymptotic behavior (n 1) of the exact solution found here, and compare it with the corresponding limit of Eq. (46). Unsurprisingly, we find good agreement between the asymptotic forms of the exact and approximate analytical solutions, which can serve as further justification for the method used in Sec. V B.

VI. SUMMARY
We considered a family of nearest-neighbor tightbinding models with position-dependent hopping and showed that the dynamics of zero-energy wave packets in them coincide with that of Dirac fields in a static curved spacetime. Using both numerical simulation of the wave packet dynamics on the lattice and a semiclassical picture valid at long wave lengths, we further elucidate the analogies between the position-dependent lattice models and a Dirac particle subjected to a gravitational background.
For power-law variation of the hopping with an exponent γ ≥ 1, we showed that zero-energy wave packets with momentum p 0 = −π/2 eternally slow down while approaching the point of vanishing hopping. This is reminiscent of the dynamics of a particle approaching a black hole horizon as seen by a distant observer. In contrast, for lower values of γ, the model does not produce horizon physics and wave packets back-scatter from the point of zero hopping. We also showed that wave packets with non-zero average energy, or p 0 = −π/2, never reach the horizon and instead reflect back at a non-zero distance which increases with the absolute value of energy.
To understand the observed wave packet dynamics beyond the semiclassical picture, we studied the eigenstates both numerically and in an approximate analytical way. We found that the low-energy eigenstates have a powerlaw localization of their wave function envelopes. This behavior of the zero-energy eigenstates results in the formation of ripples and the localization of wave packets observed in the simulated dynamics. The states with non-zero energy, in contrast, show two types of behavior at long and short distances from the horizon. While the long distance behavior is qualitatively similar to lowenergy states, at short distances we see exponential lo-calization. The latter comes from the fact that regions with low values for the hopping become classically forbidden for states non-zero energy. These results explain the back-scattering of wave packets with non-zero energy, as they cannot penetrate the forbidden region.
In this appendix, we provide supplementary details about the exact solution for the model with square-root position dependence of the hopping.
First, to normalize the wave functions given by Eq. (53), the Christoffel-Darboux formula is used. It applies to sequences of orthogonal polynomials in general, and for Hermite polynomials reads [45], For the particular case of the zero energy state, we obtain which uses the observation that H N −1 (0) = 2 (N −1)/2 (N − 2)!! for odd N . For even N , H N −1 (0) vanishes, and there is no zero mode. To reach the final approximate form, the Stirling formula at large N has been employed.
We can also find an approximate form for the normalization constant using the asymptotic form of Hermite polynomials at large N derived below, which yields This relation is valid for small and intermediate values of ε as long as ε √ N . It shows that for large N and finite energies the normalization factor exponentially decreases with ε.
We next consider the asymptotic (n 1) properties of the eigenstates. We invoke the Hermite differential equation which shows that for large n but small and intermediate z, the Hermite polynomials behave like sin( √ 2nz) and cos( √ 2nz). Using a more careful analysis, it has been found that [46] e − z 2 2 H n (z) ∼ (A8) This provides a good approximation for large n, provided that n n c,ε =ε 2 /4. The oscillatory part of this expressions can alternatively be derived by taking the real part of the large n limit of Eq. (46).