Isolated zeros destroy Fermi surface in holographic models with a lattice

We study the fermionic spectral density in a strongly correlated quantum system described by a gravity dual. In the presence of periodically modulated chemical potential, which models the effect of the ionic lattice, we explore the shapes of the corresponding Fermi surfaces, defined by the location of peaks in the spectral density at the Fermi level. We find that at strong lattice potentials sectors of the Fermi surface are unexpectedly destroyed and the Fermi surface becomes an arc-like disconnected manifold. We explain this phenomenon in terms of a collision of the Fermi surface pole with zeros of the fermionic Green's function, which are explicitly computable in the holographic dual.


Introduction
The physics of strongly correlated electron systems remains a major puzzle in modern condensed matter theory. The possible deviations from conventional Fermi liquid behavior are simultaneously extremely interesting and extremely hard to study. Nonetheless, evidence coming from experiments in high temperature superconductors and other strange metallic systems points out that such non-Fermi liquid systems do exist in nature and display many unconventional phenomena. The defining feature of non-Fermi liquid behavior is the absence of long lived quasiparticles anchored on a well-defined Fermi surface, which could be used as building blocks for Fermi liquid perturbation theory. Such signatures of a destruction of the quasiparticles are seen in the angle resolved photoemission (ARPES) studies of experimentally realized strange metals [1]. In high T c superconductors in the normal phase, the spectral width, or the inverse lifetime, of the fermionic excitation at the Fermi level becomes unnaturally broad in the anti-nodal directions, whereas one still observes well defined quasiparticles at the nodes [2,3] -the nodal-antinodal dichotomy. In the pseudogap phase the phenomenon of so-called Fermi arcs is even more striking: the sharp Fermi surface simply ends at a point in the Brillouin zone, which is topologically forbidden for a Fermi liquid-like system [4]. A careful theoretical understanding of these phenomena has been hampered by an absence of conceptually new approaches that do not rely on a stable quasiparticle description.
The holographic duality provides such a conceptually novel way to treat the strongly correlated systems without the need to postulate the quasiparticle description to begin with; for a review see [5]. It has been shown in the earlier works [6][7][8] that one can obtain a spectrum with or without the long lived fermionic quasiparticles depending on the parameters of the holographic model. The natural question arises of whether it is possible to achieve the transition between these regimes as a function of direction in the Brillouin zone, as it is observed in real materials. This is the question we address in the present work. Most of the earlier works in holography on single fermion spectral functions are restricted to isotropic setups, with the rare exceptions including [9][10][11][12][13][14]. In order to study anisotropy and the effects of the Brillouin zone boundary, we introduce a periodic modulation of the chemical potential, which mimics the ionic lattice, breaks the rotational symmetry down to a discrete group and introduces Bloch momenta. 1 A similar study has been performed recently in [10,14] which also includes a spontaneous breaking of translational invariance. In our case we restrict ourselves to a simpler setup that includes only an explicit periodic lattice, without spontaneous striped order. This allows us to perform the analysis of our results in the cleanest possible way. Similarly, in order to isolate the effects of the periodic lattice, we don't consider any non-minimal interaction terms in the Lagrangian of fermions coupled to gravity in holographic dual description.
As expected, we observe the generic effects of a periodic potential: the fermionic dispersion relation becomes multivalued in the first Brillouin zone due to the presence of lattice copies from the neighbouring zones. Umklapp gaps appear from the interaction between these copies and results in the formation of Fermi pockets. These basic effects had already been observed before in various holographic models with periodic potentials [9][10][11]. It simply shows the universality of Umklapp at the boundary of the Brillouin zone in fermionic responses in periodic potentials. However, we report here that, for substantially strong lattice potentials (much stronger then considered previously in i.e. [9,11]), a novel physical effect appears: the partial destruction of the Fermi surface due to the interaction of poles and zeros in the Green's function. A noticable spectral weight suppression in strong holographic mixed spontaneous and explicit lattices was observed earlier in [10,14]. However the effect which we observe is different: in our purposely simple tractable model we can completely identify that there is not just suppression, but that the Fermi surface is actually destroyed due to a collision of the defining pole with a zero in the Green's function. 2 .
Zeros of the fermionic Green's function have been observed in holographic models in several contexts. One kind of zeros, the "alternative quantization zeros", has been pointed out in the early works using bottom-up models [6][7][8] as well as in top-down constructions [16]. The existence of these zeros is understood straightforwardly within the holographic approach: they originate from the fact that, for a range of the parameters, a particular holographic theory can be treated as a dual to two distinct quantum theories on the boundary [17]. These two treatments are the "direct" and "alternate" quantization of the boundary operators and the simultaneous existence of both leads to the appearance of zeros in the Green's function of one theory, precisely at the point where the other one has poles. As we will see below, in the presence of a strong lattice potential these "alternative quantization zeros" approach the pole corresponding to the putative Fermi surface. This proximity kills the peak in the spectral response. The origin of this type of zeros is quite clear from the holographic point of view and universal in that context, but their interpretation in terms of conventional condensed matter theory remains elusive. We shall comment on possible interpretations in the discussion section. The phenomenon we observe is somewhat similar to the "pole-skipping" in holographic correlators for hydrodynamic energy-density modes, discussed recently in [18-21, 21, 22] and for fermionic modes in [23]. There one also observes the line of poles in the spectral density being cut by the line of zeros.
This holographic understanding of controlled zeros in the Green's function has already been exploited in attempts to describe the zeros in the spectrum arising from Mott physics [24][25][26][27][28][29]. One mechanism relies on an extra Pauli or dipole coupling present in the fermionic Lagrangian. It has been shown that in the particular case of massless bulk fermions this can partially convert the Fermi surface into the line of zeros. We intentionally do not include the extra coupling in our model, which allows us to distinguish the phenomenon we observe here from the one mentioned in those works. On the other hand, zeros in spectral functions can have several origins but true Mott physics is intrinsically linked to the presence of a lattice and translational symmetry breaking as studied here.
The spectral signature of zeros colliding with poles/peaks is a very identifiable characteristic and for that reason of high interest. More recently, a new zero-pole collision has been found [30] unrelated to holography. It was shown that in the presence of a quantum critical continuum coupled with two systems with discrete spectra, the spectrum of one such system has a characteristic zero at the resonance of the other. And this zero may collide with a pole. This effect is in the same class as the Fano resonance, where the spectrum of a continuum theory interacting with a discrete system has a zero at the resonance frequency of the latter. Alongside with our main finding we observe this new class of "resonant zeros" in our holographic model. The reason is that the near horizon geometry generically encodes a certain type of the quantum critical continuum and the periodic potential gives rise to the many copies of the discrete particle dispersion spectra within the first Brillouin zone. The fermion spectral function precisely probes a discrete sector coupled via continuum to other discrete spectra. The effect of these zeros is also interesting, but not so spectacular as the one from "alternative quantization" ones. These discrete-continuum-discrete resonant zeros cut through the Fermi surface and destroy the quasiparticle peak at a particular point, but they do not remove extended intervals from the Fermi surface.
Our finding that "alternative quantization" zeros interfere with Fermi-surface poles in holographic models with strong lattice potentials is theoretical and it is our thorough understanding of a peculiar aspect of holographic theories that allows us to unambiguously identify this mechanism. The result strikingly resembles the phenomenon of Fermi arcs, seen in the pseudogap. The creation of Fermi arcs was already the motivation for the holographic studies [24][25][26][27][28][29], but in our work they are directly tied to the presence of the lattice and anchored to the directional pattern of the Brillouin zone. Undoped Mott insulators are known to have zero responses in the single particle fermionic spectral function. In conventional condensed matter theory there are attempts to explain the formation of Fermi arcs in the pseudogap phase of the cuprates originating from these Mott zeros in the spectral density [31][32][33][34][35][36]. In these approaches, it is argued that, due to strong interactions, the self-energy of the quasiparticle diverges forcing the "dressed" Green's function to vanish at certain points in the phase space. These zeros do not violate the Luttinger count [37], but render the Fermi surface disconnected. It would be interesting to determine whether there is any connection between our results and this other approach, but as we shall discuss in the conclusion, there are a number of fundamental open questions that will require significant further study.
The paper is organized as follows, in the first two sections we give an overview of basic features of ordinary Fermi surfaces in a periodic potential (Sec. 2) as well as the holographic fermionic response in the case of a simple isotropic background (Sec. 3). In Section 4, we describe the model and the method. The main result is presented in Section 5 and we draw our conclusions in Section 6. The appendices are devoted to the details of our treatment in the homogeneous black hole background (App. A), the construction of periodic backgrounds (App. B) and the analysis of the fermionic response (App. C). We also explain the Brillouin zone representation of the Green's function in Appendix D.

Umklapp scattering and Fermi pockets in unidirectional potential
We start by recalling the basic features of the fermionic spectral function in a periodic potential. We wish to illustrate that at strong lattice potentials there are several non-linear responses, that are normally not considered in linear analysis of Bloch wave physics and Umklapp. Nevertheless these non-linear responses follow from a straightforward calculation. We consider a non-interacting Dirac fermionψ in 2+1 dimensions in the presence of a periodically modulated chemical potential µ(x) = µ[1 + λ sin(px)] along the x-direction only. In order to facilitate the analysis of our main results, we keep the y-direction homogeneous. The Dirac equation reads (here the bars denote a 2 + 1 dimensional toy-model in order to avoid confusion with our main treatment below) where σ i are the Pauli matrices. We introduce frequency, momentum in y-direction as well as the Bloch momentum in x-direction ψ(x, y, t) ≡ e i(kxx+kyy−ωt)ψ k (x),ψ k (x) ∼ψ k (x + 2nπ/p), n ∈ N, (2.3) to get the equation on the Bloch wave functionψ k (x), which is by construction periodic with the same period as the potential µ(x): Here, the Bloch wave functionψ k = {ψ ↑ k ,ψ ↓ k } is a 2-component spinor and the Dirac operator is a 2 × 2 matrix differential operator. Since the Bloch wave function has the same period as the potential, it is convenient to expand it in the Fourier series: where δω ≡ ω − µ for brevity. The Green's function in the Fourier mode representation is simply the inverse of the Dirac operator: G lm = (M −1 ) lm . In ARPES experiments, the most important element of the Green's function is the G 00 component, which characterizes the projection of the fermionic linear response function in a crystal on the plane-wave states of the incident photon and photo-electron. Therefore, in what follows, we will be focusing on the spectral function associated with the G 00 component of the fermionic Green's function A(ω, k) = Im TrG 00 (ω, k), where the trace is taken over the spin states. Throughout the paper we assume a practical definition of the Fermi surface as the locus of maxima of the spectral density in the momentum plane at the Fermi level. In Fig. 1, we show various examples of the spectral density in the toy model (2.4) for various parameters of the background potential. 3 In the left column we consider a weak periodic potential (λ = 0.08) for three values of the size of the Brillouin zone p. For a weak periodic potential all the results can be easily understood from linear analysis. When the BZ is large enough, the Fermi surface does not reach the BZ boundary and has the same isotropic shape as in the case without modulation with Fermi momentum k f . For a smaller BZ (smaller p, middle panel), lattice copies of the Fermi surface from the neighbouring Left column: At weak modulation amplitude λ = 0.08, the shape of the Fermi surface appears circular as in the absence of a lattice. However, lattice copies appear in the neighboring Brillouin zones (BZ). When the BZ (red dashed gridlines) becomes smaller than the Fermi momentum, these copies overlap and an Umklapp gap is opened at the intersection point, giving rise to Fermi pockets. Right column: At strong modulation amplitude λ = 1 nonlinear effects are seen. The shape of the Fermi surface is now affected even far from the boundary of the BZ. When BZ gets smaller, the FS is first squeezed and then strong Umklapp gaps are opened, leading to the small dumbbell-like Fermi pockets near k x axis and the nearly k x independent flat "band" at finite k y .
Brillouin zones are seen. These copies are visible in the G 00 component due to the offdiagonal mixing terms iλµ/2 in Eq. (2.6). Therefore, the spectral density in these lattice copies is suppressed by a factor of λ and, in case of the weak potential, they are hardly visible. When the Brillouin zone size is further reduced (p/2 < k f ), the neighbouring Fermi surfaces start to overlap and Umklapp gaps are opened at the Brillouin zone boundary. From the point of view of the inversion of matrix (2.6), this is simply a linear eigenvalue repulsion effect due to non-zero off-diagonal terms iλµ/2. This occurs around the point when the eigenvalues of the top-left and bottom-right blocks of the matrix (2.6) become identical. In this case, the series of circular Fermi surfaces turn into a series of Fermi pockets in addition to a an outer band parallel to x-axis. This "band" is an artefact of the unidirectional modulation of the potential, which we introduced for simplicity. In the case of a realistic crystal lattice in both x-and y-directions, Fermi pockets would form in the y-direction as well.
For strong modulating potential λ = 1, right column of Fig. 1, the situation is more involved. Even though qualitatively the response is the same, quantitatively non-linear effects will now also affect the shape of the Fermi surface. Firstly, of course, the lattice copies of the Fermi surface acquire a larger spectral weight due to the stronger mixing and are visible already for large Brillouin zones (top row). As the BZ decreases, the strong interaction effects that change the shape of the Fermi surface are seen. It becomes "squeezed" as in the middle row of Fig. 1. Finally, once the Fermi surfaces overlap, the umklapp gap that opens is so large that the outer "band" gets pushed far away and is almost flat (independent of k x ), while the heavily deformed Fermi pockets are stretched along the BZ boundary, deforming in a dumbbell-like shape, see the bottom row of Fig. 1.
As we will see below, these different types of the Fermi surface geometries that we find in the toy model (2.4) with non-interacting electrons on top of the periodic potential will also appear in the fermionic response of the strongly coupled holographic model. This will aid us in distinguishing effects that are due to strong self-interactions from effects that are due to strong lattice potential.

Holographic Fermi surfaces and zeros
Next we recall the universal presence of zeros in the fermionic spectral response in a finitedensity holographic model of a strongly interacting system of fermions. This is also so in the absence of a lattice and the simplest example of such a system -the homogeneous model at finite temperature and chemical potential -is described by a Reissner-Norström black hole in the 3+1 dimensional curved AdS space [5]. The corresponding metric reads and z h being the radius of the black hole horizon which is related to the temperature and chemical potential in the dual theory: The gauge field potential of the charged black hole is simply In what follows we will set z h = 1 by choosing the appropriate measuring units.
In order to study the fermionic response in the holographic framework, one considers the Dirac equation on the curved space-time (3.1) [5][6][7][8][38][39][40][41]: where Ψ = (ψ ↑ , ψ ↓ ) T is a 4-component 3+1 dimensional Dirac spinor, which we break into two 2-component subparts corresponding to the different spin states on the boundary. Here e f µ are the tetrad vectors; ω a bµ is the spin connection 4 ; η ac is the Minkowski (−, +, +, +) metric; A µ is the gauge field, and σ ab ≡ 1 2 [Γ a , Γ b ] is the generator of Lorentz transformations, where we choose the Γ-matrices to be The tunable parameters here are q and m -the charge and mass of the bulk fermion field. In the homogeneous background (3.1) one can immediately expand in plane waves along {x, y} directions. It is also convenient to rescale the spinor and introduce the new fields ζ a : where g µν is the metric and g its determinant. After this redefinition, the equation for ζ ↑ reads Taking advantage of the isotropy of the RN black hole, we choose the coordinates in such a way that k y = 0 and the Weyl spinors ζ ↑ and ζ ↓ decouple. For ζ ↓ , the k x momentum term has the opposite sign. Therefore, in the subsector k y = 0 the two spinors describe leftand right-moving modes in the x-direction. This simplification generically doesn't happen in the presence of a lattice, however it still arises when the fermion propagates along the unidirectional potential, as we will see in the next section 5 .
4 Spin connection is defined as via ∂µe a ν + ω a bµ e b ν − Γ τ µν e a τ = 0, where Γ τ µν is Christoffel symbol. 5 This is a consequence of our choice of the Γ-matrices (3.6). Their reduction on the boundary, defined by the action of 2+1 dimensional Lorentz generator σ µν , µ, ν ∈ (t, x, y) on a positive subspace of Γ z coincides withγ µ in (2.4). The different spin states which we consider are the eigenstates ofγ y and they decouple in case when the fermions propagate along x-axis, since in the corresponding 1+1 dimensional problem they possess different chiralities.
The analysis near the AdS boundary z → 0 shows the spinor components behave as where the constants a 0 r , b 0 s can be fixed by solving the equations of motion, see Appendix A, and can in principle depend on parameters (ω, k). Given the expansion (3.9) of the bulk fermion profile one identifies the independent coefficients a and b with respectively the source and the expectation value of the corresponding fermionic boundary operator Ψ of conformal dimension ∆ = 3 2 + m [6,7,38,39,[42][43][44]. This is called "direct quantization" (Direct Q). It is important to note however, that this identification is not unique: for m ∈ [0, 1 2 ) one can consider an "alternative quantization" (Alt. Q, which we denote with tilde:Ψ) and consider b as the source and a as the expectation value of the operator Ψ with conformal dimension∆ = 3 2 − m. In this regime of m ∈ [0, 1 2 ), a single bulk model can correspond to the two distinct boundary theories, depending on the type of quantization chosen. A useful way of studying alternative quantization is to consider the direct framework but to extend the range of m to negative values of m ∈ (− 1 2 , 0]. In this way the roles of a and b are switched in (3.9).
Within the linear response approximation, the holographic identification of source and response allows one to obtain the two-point function of the fermionic operator under consideration: the fermionic Green's function G = Ψ † Ψ = (response)/(source). In general the Green's function is a matrix, and to divide out the sources requires a few steps. In the notation of [40], given the relation between the sources and responses 6 the (retarded) fermionic Green's function reads where γ t = iσ 1 with our choice of the Γ-matrices (3.6) [6,7,38,40,45,46]. Therefore, in direct quantization G ↑↑ = b ↑ /a ↑ , while in alternative quantizationG ↑↑ = a ↑ /b ↑ . In other words, the two point functions in both quantizations are related by In particular and importantly, this means that the poles of the Green's function in the alternatively-quantized boundary theory correspond to the zeros of the Green's function in case of the direct quantization. The final point is to determine which type of the Green's function we are looking at. This is fixed by the boundary conditions at the black hole horizon, see Appendix A. The retarded Green's function corresponds to a purely infalling solution at the horizon: (3.13) 6 Note that in the boundary fermionic theory the a ↑ is a source to b ↓ and vice versa, see [45,46] In short, in order to evaluate the Green's function, one has to solve the Dirac equation (3.8) in the full bulk geometry (3.1) and find the ratio between the a and b branches of the AdS boundary expansion (3.9) of the solution. The first-order ordinary differential equation (3.8) together with the boundary condition (3.13) is solved with the numerical shooting method, as we explain in Appendix A, and for given values of the background parameters {T, µ}, fermionic parameters {q, m}, frequency and momentum {ω, k x = k}, the asymptotic form (3.9) is read off from the solution. This gives us the Green's function G ↑↑ (ω, k) from Eq. (3.12). In complete analogy we can evaluate G ↓↓ by solving the appropriate equation, or simply use the symmetry mentioned above G ↓↓ (k) = G ↑↑ (−k). Finally, we evaluate the fermionic spectral density as which is the central object of our study. In what follows we mostly focus on the features of the Fermi surface, defined as the locus of the spectral density peaks at zero frequency: The fermionic spectral density (3.14) at zero frequency in the RN black hole (3.1) computed this way is shown on Fig. 2. We have chosen the fermionic bulk mass m = 1/4, the charge q = 1 and T /µ ≈ 0.005. The Fermi surface with direct quantization is shown on the left panel. Moreover, since m < 1/2 alternative quantization is also possible and the result is shown on the right. The Fermi surfaces are circular, as expected for fermionic excitations at finite chemical potential in an isotropic background. We also observe the specific holographic feature of the appearance of multiple nested Fermi surfaces [7]. This is a generic feature in holographic models: the number of the observed Fermi surfaces depends on the mass and the charge of the bulk fermion as well as on the background chemical potential, and in principle it can be arbitrary [5].
Here we wish to call attention to another interesting phenomenon which is evident when comparing the results of direct and alternative quantizations: each plot in Fig. 2 includes the position of the Fermi surface in the other quantization in dashed lines. We see that the Fermi surfaces in the direct and alternative quantization alternate. This effect is more clearly visible on the left panel of Fig. 3, where the k y = 0 cut of the spectral density is shown. Indeed, the solid and dashed vertical grid lines, indicating the Fermi surfaces in the direct and alternative quantization, respectively, alternate: there is always an alternative quantization FS in between two direct ones.
At this point it is important to remember that, due to the inverse relation between the Green's function in Direct Q. and Alt. Q. (3.12), a pole in the alternative Green's functioñ G always corresponds to a zero in the Direct Green's function G. While a pole near the real axis produces a discernible peak in the spectral density A, a zero is only reflected in a depletion of spectral density with a minimum set by the imaginary part of the position of zero in the complex plane, which is proportional to temperature. This depletion, unlike the peaks, is harder to spot in density plots like Fig.2. The reason is simply that in order for the depletion in spectral density to be visible, the latter must have a finite background value. However, the spectral density at zero frequency is set by temperature itself and therefore the depletion is not seen. Therefore, the most direct way of detecting the zero in the Green's function is indeed to study the peaks in its alternatively quantized counterpart. This is how we will identify the zeroes of the Green's function in the remainder of this paper.
Another reason for the absence of detectable depletions of the spectral density is that the trace of the Green's function contains several additive terms, only one of which is suppressed. Indeed, in a particularly simple example with k y = 0 the Green's function is diagonal in the spin representation Therefore, the peaks of a single component, shown in dashed in Fig. 3, are clearly seen in the total Im TrG, while the depletion, if any, would be shaded with the finite value coming from the opposite-spin component. One could alternatively capture the position of zeros by looking at the real part of the Green's function at real ω, which changes sign at exactly this point. However this method suffers from the same problem: the contribution of the two spin components add together in the trace of the Green's function and one has to diagonalize it in the spin space in order to distill the position of zeros. The Green's function zeros become particularly important in the case where the peaks in both direct and alternative quantization of the same spin component (shown with arrows in Fig. 3) come closer to each other. In this situation, the pole and the zero of the Green's function in the complex ω plane would recombine and the residue of the pole would vanish In this way, a peak in the spectral density would be "eaten" by the approaching depletion producing a distinct observable phenomenon. This pattern is similar to what happens in a Fano resonance in a continuum coupled to a discrete system, but here it is manifested as a destruction of the Fermi surface at ω = 0. In the simple isotropic model described above this does not happen since the poles and zeros are always separated. In what follows we will show that the situation changes when a strong periodic background chemical potential is considered. Before moving forward, another comment regarding the existence of the zeros in the Green's function is in order. On the right panel of Fig. 3 we show the results obtained for a different mass m = 3/4 for which only a single quantization in the dual CFT is possible [6,43]. In other words, one cannot prescribe a physical meaning to the Alt. Q Green's function. Formally, however, we can still evaluate the Alt. Q expression by inverting the Direct Q result. The artificially computed Alt. Q spectral density clearly exhibits the peaks even in this case, which therefore correspond to the zeros of the Direct Q Green's function. Therefore, we conclude that the general mechanism of appearance of the zeros, is independent of whether or not the alternatively quantized picture is well defined.

Fermionic spectral function in a holographic lattice
We now turn to the study on how the fermionic spectral function discussed in Sec. 3 is affected by the periodic chemical potential, and in particular when the strength of the lattice potential is strong. The bulk model with periodically modulated chemical potential -the holographic ionic lattice -has been introduced in [47] and studied in great detail in [48,49]. In what follows we will adhere to the framework used previously by some of us in [50][51][52]. We consider the holographic model with Einstein-Maxwell action where F = dA is the gauge filed strength tensor, R is the Ricci scalar and Λ -the cosmological constant. We introduce the spatially modulated chemical potential (c.f. (2.2)) The non-isotropic black hole solution can now no longer be constructed analytically, but must be found numerically. In order to find the gravitational background solution it is sufficient to consider the following metric ansatz The blackening factor f (z) is that of the RN black hole (3.1) and all the other ansatz functions depend on both z and x coordinates. Given that, at the horizon (Z 2 − T 2 ) z=1 = 0, the temperature is still given by (3.3). Using the DeTurck trick [53][54][55] and the numerical methods for solving partial differential equations (PDEs) developed in [50,51,56] we obtain the background gravitational solutions for a given temperature T at a fixed chemical potential. The lattice wave-vector is p and λ is the amplitude of the chemical potential modulation in units of µ 0 . The details of the numerical procedure and precision control are discussed in Appendix B. With the gravitational background solution at hand we proceed to solve the Dirac equation (3.5). The co-frame is now less trivial than for the isotropic case: The equation for the ζ ↓ component is obtained by a parity transformation: Given a solution of (4.6) as a function of both z and x, we extract the near boundary coefficients a(x) and b(x) defined in Eq. (3.9). In order to obtain the linear response S-matrix of Eq. (3.10), we expand these periodic expansion coefficients in Fourier series similarly to the method shown in Sec. 2: Therefore, the S-matrix is a infinite matrix with both spin and Fourier indices (c.f. (2.6)).
The Green's function is evaluated in the same way as in Eq. (3.11), except that it is now a matrix. It is worth mentioning that, since the Fourier basis exponents differ from each other by a shift with exactly one unit of the lattice wave vector p, the indices m, l can be interpreted as the Brillouin zone index. Therefore, the Green's function G ml is simply a matrix in the Brillouin zone representation. We comment on some features of its structure in Appendix D.
As it was discussed earlier in Sec.2, the most interesting component of the Green's function for us is G 00 , which measures the overlap between the response function in the material and the plane waves of the ARPES probe (2.7). In order to measure G 00 and the associated spectral density, we consider a plane wave source a(x) = 1 as a boundary condition when solving the bulk equations of motion (4.6) and read off the homogeneous component of the response b(x), see also Appendix C.
In order to evaluate the Green's function in the alternative quantization picture, one has to invert the full infinite matrix and the simple formula (3.12) turns intõ In order to perform this inversion one has to evaluate all the components of the Direct Q Green's function, which is hard in practice even when one truncates the Fourier series.
Therefore, instead of using this method we evaluate the alternatively quantized Green's function by solving the Dirac equation with the negative fermionic mass parameter (see Appendix A). As explained above, this switches the roles of a(x) and b(x) in the linear response calculation. In Appendix D we check that the two approaches give identical results.

Destruction of the Fermi surface by the zeros in the Green's funciton
Let us now analyze the results which we get for the fermionic spectral function in the holographic model described in Sec. 4. In what follows we study a series of the background gravitational solutions for various Brillouin zone sizes p, but with fixed charge density ρ. We also fix the temperature T / √ ρ, 7 and the amplitude of the chemical potential modulation λ (4.2). We consider two cases with weak and strong potential modulation, in direct analogy with the toy model study performed in Sec. 2. We use the same temperature and parameters for the bulk fermion as in the homogeneous RN-black hole case addressed in Sec. 3: As seen in Fig. 2, the size of the Fermi surface in the absence of the periodic potential for these parameters is k f ≈ 0.94 √ ρ.

Weak lattice potential
We start with a weakly modulating lattice with strength λ = 0.1. On Fig. 4, we show the momentum distribution of the spectral density, which in this case displays the circular Fermi surface with exactly the same size as the homogeneous one (c.f. Fig. 2). The weak lattice potential does not affect the shape of the Fermi surface. We choose the size of the Brillouin zone to be smaller then the Fermi momentum p < 2k f , therefore the umklapp copies of the FS overlap and, in perfect agreement with the observations in toy-model of Sec. 2, we see umklapp gaps opening at the BZ boundaries. Another feature which is expected is the suppression of the spectral density in the neighboring Fermi surfaces. Indeed, these are almost invisible in the linear scale plot on top of Fig. 4. However, the bottom panel shows the logarithm of the spectral density, which makes the λ-suppressed lattice copies of the Fermi surface clearly visible. In a nut shell, the results obtained in the holographic model with a weak periodic potential are in perfect agreement with both the standard logic of fermion physics in a periodic potential, and the shape and size of the holographic Fermi surface in the homogeneous background. This serves as a consistency check of our approach and numerical techniques. Similarly, a familiar pattern is also seen on the energy-momentum resolved spectral function shown on Fig. 5. As expected, the usual holographic dispersion relation for fermionic excitations is recovered. In this case, it consists of the two nested Dirac cones in the vicinity of the Fermi level. The spectral lines quickly broaden due to the quantum critical self energy Σ(ω) ∼ ω 2ν k F [5][6][7][8]. However, there is one distinctive new feature due to the presence of the holographic lattice. There is a localized depletion of the spectral weight along the lines which correspond to the dispersion bands of the neighbouring Brillouin zones. This feature is a consequence of a general effect pointed out earlier in [30]. The density plot shows a cut through the same data as Fig.4 in the k x , ω plane, at k y = 0.001 √ ρ. Around (k x , ω) ≈ (±0.4, 0.02), the dark lines crossing through the inner Dirac cone arise from the "resonant" type zeroes [30] discussed in Sec. 1.
Namely, when multiple systems with discrete spectra interact with each other by means of the quantum critical continuum, the spectral function of one of the systems develops isolated zeros, or depletions, at the positions of the energy levels of the other systems. In the case of the holographic model with a lattice, the discrete systems are the umklapp copies of the fermionic dispersion while the quantum critical continuum originates from the near horizon geometry [8]. This interesting effect may, in principle, also lead to a destruction of the Fermi surface at an isolated point, however it will not be relevant in the present study.

Strong lattice potential
Let us now turn to the analysis of the case of strong potential modulation λ = 8, depicted in Fig. 6. Similarly to the toy model study of Sec. 2, we consider a series of backgrounds with different potential wave-vector p which sets the sizes of the Brillouin zone. The parameters considered allow to access the alternative quantization, which has an inverse Green's functionG = G −1 . In order to analyze the spectral density in Alt. Q, we solve the Dirac equations (4.6) on top of the same gravitational background lattice, but taking the negative sign for the bulk fermion mass m = −1/4 (see Appendix D for the details of obtaining the Alt Q fermionic response). We show the results of the direct and alternative quantizations side by side on Fig. 6. Note, however, that we choose logarithmic scale for the Alt. Q results in order to resolve all the features of the FS. On the first row of Fig. 6, the BZ is much larger then the Fermi surface in direct quantization (top left, p = 5 √ ρ ≈ 5k f ). Therefore, since the umklapp surfaces are far away, the FS is not deformed even in this regime of strong lattice potential. This is expected; when momentum is much smaller than the BZ, i.e. in the long wave length limit, the periodic structure of the potential becomes irrelevant and only the mean value of the chemical potential µ 0 plays a role in the formation of the Fermi surface. 8 On the other hand, the outer FS in the alternate quantization is larger (k f ≈ 1.5 √ ρ, see Fig. 2).
Therefore, on the top right panel of Fig. 6, we see that, in this case, the FS is already deformed by the BZ boundaries, similar to the non-interacting toy model in Sec. 2. When the BZ becomes smaller, p ≈ 3k f , second row of Fig. 6, the neighbouring Fermi surfaces come close to each other and get deformed due to the strong lattice potential in direct quantization. This is exactly the same situation observed in the toy model of Sec. 2. For Alt. Q (right panel), we readily observe the formation of Fermi pockets and the flat outer band. This is again similar to the toy model at strong lattice potential.  figure 6, it is clearly seen that the Fermi surface loses its sharpness and even disappears around k y = 0.
The novel interesting phenomenon arises when the Brillouin zone is squeezed even further. On the bottom row of Fig. 6 we show the situation with p = 1.4 √ ρ. Quite strikingly, we see that as the FS is squeezed further more, the sharp spectral density peaks indicating the shape of the Fermi surface disappear along the boundary of the BZ in direct quantization. A more detailed view in Fig. 7 shows this explicitly. On the other hand, the alternative quantization plot displays heavily squeezed dumbbell-like Fermi pockets, which are centered at the BZ boundaries. As discussed in Sec. 3, we have a thorough understanding of some of the peculiar features of the holographic fermionic spectral response. This allows us to figure out the origin of the depletion of spectral function. The zeros of the direct quantization are the corresponding poles of the alternative quantization Green's function. In Fig. 8, we show high resolution results of the FS for Direct and Alternative quantizations. In each plot, we show with the dashed green lines the positions of the poles in the other quantization scheme. Most strikingly, the position of the secondary Fermi surface in the Alt. Q overlaps with the primary Fermi surface in the Direct Q, as seen on the left panel of Fig. 8. Therefore, there are zeros in Direct Q that get pushed towards the poles defining the FS. This effect destroys the Fermi surface in an extended region near the BZ boundary. This cancellation of poles and zeros is the fundamental reason for the novel phenomenon we observe.
We can expose the "zero-eats-pole" effect in more detail by analysing cuts along the k x -axis of the spectral density. These Momentum Distribution Curves (MDC) are shown in Fig. 9, where we plot the Direct and Alternative quantizations, in exactly the same fashion as we did in Fig. 3. A clear depletion is observed in the Direct Q MDC corresponding to the secondary FS peak in the Alt. Q scheme. As the size of the Brillouin zone (red dashed grid line) is decreased as shown in bottom left plot of Fig. 9, this depletion approaches the peak shown by filled red pointer and absorbs it completely at p = 1.4 √ ρ as shown in the bottom right plot of Fig. 9. Therefore, we conclude that for strong lattices, the Fermi surfaces in both direct and alternative quantizations are deformed in such a way that the poles and zeros overlap and cancel each other. This destruction has a dramatic effect on the very existence of the quasiparticle excitations in part of the kinematic region. This is best seen on Energy Distribution Curves (EDCs), which we show on Fig. 10. Here, we plot the frequency dependence of the spectral weight at two perpendicular directions in the Fermi surface: f . When the BZ boundary is far away (top panel of Fig. 10), the EDCs are practically identical in both directions, confirming that the FS is almost isotropic. In the bottom left panel, the anisotropy is now manifest, since the FS is deformed by the lattice. Finally, the most drastic effect is seen on the bottom right panel of Fig. 10. On the one hand, in the k y -direction, there is still a sharp peak corresponding to a quasiparticle excitation. However, the spectral density in the k x direction (red line) is totally incoherent -there is no excitation with definite energy which propagates along the k x -direction.

Discussion
In this work we have studied the fermionic spectral function in a holographic model with periodic ionic lattices. Our most important finding is a novel phenomenon which appears as a destruction of coherent spectral weight peaks of the Fermi surface along the directions of the lattice vector. We have shown that the origin of this phenomenon is the interaction between poles and zeros of the fermionic Green's function. More precisely, due to nonlinear effects from the strong background lattice potential, these are pushed together and suppress each other. What is tantalizing, is that the patterns of the spectral density observed here look very similar to the results of ARPES experiments in strongly correlated materials in which nodal-antinodal dichotomy and Fermi arcs are observed. There one also observes the destruction of the coherent spectral weight in certain directions. It is therefore very interesting to further investigate whether our results can clarify these unconventional phenomena.
A warning is in order however; one should handle the results of holographic models with great care. As we discussed in Sec. 3, the presence of zeros in our holographic treatment follows from the simple logic relying on the existence of two possible choices for the dual CFT operators corresponding to each quantization, and the existence of a Fermi surface in the alternate quantization. However, as we noted on Fig. 3, when only one quantization is At large lattice momentum, when the BZ boundary is far away from the FS, the the peaks in direct and alternate quantization alternate along the MDCs (c.f. Fig. 3). The positions of Direct Q peaks and Alt Q dips are correlated. When the BZ boundary is brought closer, the peaks and zeroes seem to absorb each other. For p = 1.4 √ ρ, the quasiparticle peak is eaten by the zero completely.
allowed, zeros are still present and perhaps their existence is a more fundamental feature, which is yet to be explored. The fact that the zeros and (the multiple) poles of the fermionic Green's function appear in simple models as alternating series suggests that, as we drive some parameter (chemical potential or temperature) to zero, they might coalesce and form a branch cut in the complex plane. This signals the truly unparticle shape of the Green's When the BZ boundaries are far away, as in the top picture, the EDCs look very similar: they both show a well-defined, sharp peak. The dip on the right shoulder of the peak corresponds to the alternate quantization zero. When the BZ boundary is brought closer, first some asymmetry in the sharpness of the peaks appears. When the BZ boundaries are brought very close, the peak for the EDC near the BZ boundary gets destroyed by the zero completely.
function characteristic of the ultraviolet CFT of the holographic duality. This idea is rather at odds with the physics of experimentally realized strange metals, because the UV CFT of a holographic model is clearly not the UV theory governing the behaviour of electrons in real materials. Moreover, the multiple pole features of the fermionic response in the holographic model arguably rely on the large-N approximation, which may or may not be efficient in reality. However, one can regard our results with a larger perspective; we observe that the destruction between zeros and poles requires the two fundamental ingredients: the very existence of zeros in the fermionic Green's function and strong lattice potentials, which influences the position of poles and zeros by bringing them closer together so that eventually they annihilate each other. Even though the understood origin of the zeros in our particular holographic model cannot be directly mapped to known condensed matter systems, there are models in condensed matter theory, which display similar features, in particular Mott insulators [31,37]. Our finding suggests that, if upon doping the insulating Mott zeros of the Green's function remain present, their proximity to the quasi-particle poles may indeed be the origin of this spectacular experimental phenomenon. Doped Mott insulators are exactly the systems where the real Fermi arcs are observed. In this regard holographic studies certainly confirm its use as a convenient theoretical laboratory for studying the phenomenology of the strongly correlated systems which goes far beyond the applicability of the Fermi liquid theory.
In order to render the solutions regular on both ends of the integration interval it is therefore convenient to rescale the coordinate as well. More specifically, we define With these modifications, the problem reduces to solving two coupled first-order linear ODEs (3.8) on the interval r ∈ [0, 1] subject to the boundary conditions (A.2), (A.4). There are 3 integration constants a, b and h, one of which can be set to unity due to the linearity. The remaining two are fixed by the two boundary conditions giving a unique solution for fixed ω, k x , k y . We use the numerical shooting method to obtain this solution. We shoot from both ends of the interval, using the expansion series (A.2) and (A.4) as the initial conditions, and look for the values of free constants b, h for which the solutions match in the arbitrary point in the interior of the interval. The advantage of this method is that we have a direct control over the response constant b and do not need to extract it from subleading terms in the near boundary expansions of the solution as it is usually done when shooting from horizon only.

B Numerical calculus and precision control for gravity background
In order to obtain the non-homogeneous gravitational background with a periodic boundary condition for the gauge field (4.2), we have to solve the full set of the Einstein equations. These follow from Eq. (4.1), which gives a set of 6 coupled non-linear partial differential equations (PDEs) in coordinates (x, z), which are not elliptic. In practice we solve for functionsQ µν ,Â µ defined as such that the RN black hole solution corresponds to the choiceQ µν (x, z) = 0,Â t (x, z) = 1. We use the DeTurck trick, outlined in [53][54][55] and used for this type of holographic models in [47,49]. It allows to recast the Einstein's equations as a boundary value problem for a set of nonlinear elliptic equations. We choose the RN black hole as a reference metric, which allows us to use for the non-homogeneous background the same expression of temperature as a function of chemical potential as that of the RN black hole (3.3).
The boundary conditions in x-direction are periodic, as dictated by the symmetry of the background lattice potential. It is convenient to rescale the spatial coordinate x → (2π/p)x in order to fix the integration interval to unity:x ∈ [0, 1). At the z → 0 boundary we require the metric to be asymptotically AdS with a nontrivial inhomogeneous source for the gauge field (c.f. (4.2)) z → 0 :Q µν (x, 0) = 0,Â t (x, 0) = 1 + λ cos(2πx). (B. 3) The horizon boundary conditions follow from the asymptotic expansion of the Einstein equations near the horizon f (1) = 0. We expand the unknown functions in the Taylor series down to the first derivative order and substitute these expansions into the Einstein equations. This gives us 5 equations in the subleading order which relate the derivatives of the fields to their horizon values, i.e. generalized Robin boundary conditions. On top of that we get one algebraic equationQ tt (x, 1) =Q zz (x, 1) at the leading order, which guarantees the solution to be static [47][48][49].
The numeric algorithm follows closely the treatment of our earlier works [50,51,56] with several technical improvements. We used a finite difference discretization of the equations on a homogeneously spaced grid and realized a Newton-Raphson procedure to solve the nonlinear boundary value problem. While supplementary procedures were implemented, including a relaxation scheme with Orszag regularization, no other methods matched the speed and accuracy of the direct Newton-Raphson method [59,60], see also [56].
We ran a comprehensive series of precision and accuracy tests to make a judicious choice of computational parameters. We employed the norm of the DeTurck vector and the value of the trace of the Einstein tensor as measures of convergence. The thermodynamic potential was used as a physical observable, from which we could estimate the relative precision of our calculations. From this analysis, we were able to conclude that the configuration which gives us the best balance between speed and accuracy of computations for both the backgrounds and the fermions is a grid of n x × n y = 34 × 80 points, where we use 8th-order accuracy central finite difference derivatives in our finite difference scheme. This achieved relative errors of less than 10 −5 in the thermodynamic potential, while still allowing for the fermionic equations to be computed quickly and accurately enough. In order to make sensible comparisons, we work in the canonical ensemble with a constant charge density. Since this is not a parameter in our equations, but rather an observable we can only extract after solving the background, we employed a root-finding algorithm with a relative tolerance of 10 −3 to fix the charge density.
We implement our numerical routines in Python 3.6 using a set of packages that can be found in the standard SciPy stack [61]. Under the hood, these use the SuperLU [62] package for solving the sparse linear equations. We have performed several cross checks with a similar code in Mathematica [63], which was used earlier in [50][51][52], in order to prove the reliability of the package. The code was designed such that it could run a large number of small instances in parallel on any number of machines. This technique was suitable for running on the Lorentz Institute Maris cluster

C Numerical calculus for Dirac equation
There are few subtleties which one encounters when solving the Dirac equation on the non-homogeneous background (4.6) numerically. First, in order to make the fermionic solutions regular on both horizon and asymptotic boundary we perform the same set of redefinitions of the wave-functions (A.1) and holographic coordinate (A.5) as it was done in the homogeneous case discussed in Appendix A.
Moreover, there are several strategies which one can use in order to solve the set of first order differential equations numerically. One of them is to solve a Cauchy problem, integrate the equations starting from horizon and read off the boundary asymptotes. This is somewhat similar to the shooting method we discussed above and it was used in the early work [11]. However, when dealing with PDEs, one has to consider solving for all possible Fourier modes on the horizon and tracing out all the components of the response matrix S nl (3.10) before one gets access to the single S 00 component, which we are mostly interested in. Therefore we find this approach not suitable. Instead, we consider the boundary value problem, imposing the (position dependent) infalling boundary conditions (A.4) at the horizon as well as setting the source a(x) on the asymptotic boundary z → 0 to a desired Fourier mode. In this way, setting for instance a(x) = e ipl with fixed l and measuring the full profile of b(x) = n e ipn , we get the information on all components of the response matrix in a row S ln , n ∈ [−N/2, N/2], where N is the size of the grid in the x-direction. This boundary value approach is more useful since we can obtain the desired value S 00 in a single iteration by considering simply a(x) = 1 as boundary condition.
However, a problem arises when one tries to implement the boundary value PDE solving code in the case of a Dirac equation. The boundary value problem requires setting the boundary conditions for each field on both ends of the integration interval. For a first order differential equation, like the Dirac equation, this over determines the problem. We avoid this obstacle by formally setting a trivial boundary conditions of the form '0 = 0' for half of the fields on the boundary. Indeed, after performing the rescaling of the fields (A.1) we guarantee that the sub-leading spinor component will behave as ξ r ∼ z 2m at the boundary. Therefore, setting ξ r (0) = 0 does not impose any extra constraint in the problem. This precise trick does not work however in case of alternative quantization, when we formally take m < 0 and the "response" branch ∼ z 2m diverges on the boundary. In the particular case of m = −1/4 studied here, this divergence is mild enough and we simply get rid of it by further redefining the "response" fermion component with an extra factor of z. More precisely, in case of m = −1/4 instead of (A.1) we use ζ r → z −m−1 (1 − z) −i ω 4πT ξ r . This changes the near boundary behavior of the ξ r component to ξ r ∼ z 1/2 . As before, the Dirichlet boundary condition ξ r = 0 is trivial and does not over constrain the problem.
We use a similar logic on the horizon; in addition to the leading behaviour ξ r (x) = iξ s (x), Eq. (A.4), following from the ingoing boundary condition, we include the subleading terms in the expansion of the equations of motion,. These, relate the derivatives of the functions to their boundary values. Since these relations are obtained from the equations of motion themselves, they do not introduce extra constraints and we are left with the correct amount of boundary conditions.
Another problem with the first order differential equations is that the matrix which represents the discretized problem on a lattice does not have a positively defined spectrum of eigenvalues typical of elliptic problems. Therefore, one can not rely on the iterative methods, because these are not guaranteed to converge. In our case, the Dirac equation is linear and one does not require to use iterations of any kind: the problem is solved "in one shot" by inverting the master matrix once. Even though we managed to make use of this approach, it may not always be applicable, especially in cases where the pseudospectral collocation is used and hence the matrix is dense, or simply if the grid is dense and inversion of the huge matrix is not feasible. In this case one can improve the situation by substituting the first order equations (4.6) with the second order elliptic ones. Indeed, if we represent the equations (4.6) for ξ ↑ and their P-transformed counterparts for ξ ↓ as we can apply to one of the equations the differential operator of the other equation, and vice versa, to get the following second order system In this form, the second order differential operator for each spinor component is the Laplacian in curved space. This representation allows us to use the full power of the iterative numerical techniques designed for elliptic equations. We have observed that the direct inversion of the master matrix, similar to the one we used in the first order case, becomes numerically less demanding due to improved features of the differential operators. When using the second order equations (C.2) one has to take care that no ghost solutions are obtained which do not solve the original problem. We can guarantee this by using the expanded first order Dirac equations as the boundary condition on the horizon. Unlike the first order case, these boundary conditions are not trivial, but they rather impose the constraint that the first order equations are satisfied at the horizon, and this constraint is further propagated in the bulk by the second order system. We checked in our numerical calculations that the two approaches, with the first and the second order differential equations, give identical results and this cross-check serves as a good confirmation of the validity of our numerical treatment.

D Green's function in the Bloch momentum representation
Here, we study spatial features of the Green's function in the Bloch momentum representation and the way the poles in the alternative quantization manifest themselves as zeroes in direct quantization. We focus on the response matrix (4.8). In practice, by representing the Dirac equation as a boundary value problem, as outlined in Appendix C, we have a direct control over a(x) and can, for instance, source any given harmonic mode. In most cases we just switch on a 0 = 1, corresponding to the constant source and obtain b(x) as the series b (0) (x) = n e inpx S n0 . (D.1) Then, we extract the b 0 homogeneous component, which gives S 00 = b 0 /1 (note that we normalized a 0 = 1 here). In principle, we can go further and study a full rectangular sector of a (formally infinite) S-matrix. In order to do so, we consider several harmonic sources a (l) (x) = e ilpx and evaluate the responses for these cases b (l) (x) = n e inpx S nl . (D. 2) It is now clear that by extracting N Fourier modes of the response profiles b (l) (x) for a set of L harmonic sources we get access to the full N × L patch of the S-matrix. Formally, in order to treat the alternative quantization as exchanging the of roles between b(x) and a(x) and obtain theS 00 component of the Alt. Q response matrix, we have to guess what kind of the boundary "response" condition a(x) would lead to a constant "source" b(x) = b 0 . In other words, our goal in this case is to find a set of coefficients a l , such that And after taking the a 0 component we arrive at the expression:S 00 ≡ b 0 /a 0 = (S −1 ) 00 , from which the equation (4.9) of the main text follows. As we see here, if S 00 has a diverging value, it will enter the determinant of S and therefore force the (S −1 ) 00 component to vanish in complete analogy to the simpler homogeneous case we studied in Sec. 3. Given that using the harmonic sources we can evaluate a large enough sector of the S-matrix, we can invert it approximately and obtain a required value of (S −1 ) 00 . However this treatment is not feasible in practice. The other way of obtaining the alternative quantization picture is by directly setting the boundary condition for b(x) = 1 and read out the profile of a(x). In this way we directly measure theS 00 component and no matrix inversion is needed. This approach is equivalent to setting the bulk fermion mass to its negative value (3.9), as we discussed above and it is much less demanding, therefore we predominantly use it in this work. Nonetheless, the rescaling of the fermionic wave function discussed in Appendix C is different in this case. Consequently, the equations differ from the direct quantization ones. We find it important to check whether the two approaches do indeed lead to the equivalent results. On Fig. 11 we show that indeed the result obtained from inverting the S-matrix for a set of 34 harmonic sources (left panel) coincides with the one which we get by changing the sign of the bulk fermion mass (right panel). One extra comment on the structure of the S-matrix is in order. This clarifies the relation between the shifts in parameter k and the values of the spectral function in the different Brillouin zones. It is useful to recall that the k-parameter is a part of the definition of the Bloch wave-function (4.  Figure 11: Comparison of the two treatments of the alternative quantization picture The comparison between computation in alternate quantization and doing the inversion procedure from the direct quantization data. The main features (sharpness, shape, and the reduction of spectral weight near (k x , k y ) = (0.65, 0.40) √ ρ) are demonstrably reproduced. The inversion procedure does suffer from limited numerical precision, as the contributions from all the different Green's function components range over many orders of magnitude. The non-smoothness of the left-hand side arises because the width of the FS is smaller than the density of data points, causing the computation to not quite hit the peak. This can be improved in principle, however due to the prohibitive computational cost of this procedure we don't show it here.
Theã(x) andb(x) functions are still periodic in the unit cell, therefore (D.7) is a perfectly valid Bloch wave representation as well. However, the Bloch momentum is now different and the relation betweenb This identity serves as another useful check of our numerical procedures and our results are in agreement with it.