Dissipative analogue of four-dimensional quantum Hall physics

The four-dimensional quantum Hall effect (QHE), which features a non-trivial configuration of three-dimensional Weyl cones on its boundaries, relies on synthetic dimensions for its simulation in experiment. Here, we propose a three-dimensional analogue of this model in the form of a dissipative Weyl semimetal (WSM) described by a non-Hermitian (NH) Hamiltonian, which in the long-time limit manifests the anomalous boundary physics of the four-dimensional QHE in the bulk spectrum. The topology of the NH WSM is captured by a three-dimensional winding number whose value is directly related to the total chirality of the surviving Weyl nodes. Upon taking open boundary conditions, instead of Fermi arcs, we find exceptional points with an order that scales with system size.


I. INTRODUCTION
Topological phases of matter have been at the forefront of research in condensed matter physics over the last decades [1][2][3][4][5][6][7][8]. The communality between models describing such phases is the presence of robust states, which appear in an anomalous configuration on the boundaries of the system and whose existence is described by a topological invariant as dictated by the bulk-boundary correspondence [9]. An early example of a topological model is the four-dimensional quantum Hall effect (QHE) [10,11], which features three-dimensional Weyl cones on its boundaries captured by the second Chern number. While many proposals exist for the experimental simulation of the four-dimensional QHE ranging from photonic setups [12][13][14] to ultracold atoms [15,16] and electrical circuits [17], only three experiments have been reported [18][19][20], which make use of topological pumping in twodimensional setups [18,19] as well as quantum simulation [20]. Due to the physical obstruction of having four spatial dimensions, all of the proposals [12][13][14][15][16][17] as well as the experiments [18][19][20] necessarily rely on synthetic dimensions rendering the physical realization of this system highly non-trivial. Here, we propose an alternative approach to probe the topology of the four-dimensional QHE, namely, through studying a topologically equivalent, three-dimensional analogue with dissipation in the form of a non-Hermitian Weyl semimetal.
Weyl semimetals (WSM) are three-dimensional topological models, which feature Fermi arcs on the surfaces that connect Weyl cones in the bulk [21][22][23][24][25][26][27]. As such, a WSM can be seen as the surface realization of the fourdimensional QHE with one important caveat: The total chirality of the Weyl cones appearing in the bulk spectrum has to disappear meaning that an anomalous configuration of Weyl cones cannot exist [28]. However, if we allow for the system to be dissipative, this obstruction can be circumvented. Indeed, we find that introducing modulated gain and loss in the WSM, which renders the Hamiltonian of the model non-Hermitian (NH) (see Ref. 29 for a recent review), results in the appearance of an anomalous configuration of Weyl cones in the bulk spectrum in the long-time limit: The spectrum of NH Hamiltonians is generally complex, where the imaginary part is associated with the inverse lifetime, such that only states with Im(E) > 0 survive when t → ∞ [30].
The three-dimensional NH WSM studied in this paper is ensured to be equivalent to the four-dimensional QHE by making use of arguments that are grounded in topology. Specifically, this means that there exists a connection between the topological invariants of the two models: Whereas the chiralities of the Weyl cones that appear in the bulk spectrum of the NH WSM are found to be directly related to the values of the second Chern number of the QHE, we also find that the total chirality of the Weyl cones in the positive imaginary energy plane is given by a three-dimensional spectral winding number [31], which captures the point-gap topology of the NH WSM, where point gaps are a purely NH phenomenon with the complex energy bands not crossing a reference point [31] [cf. Fig. 2(a)]. Our three-dimensional NH WSM thus indeed mimics the four-dimensional QHE.
In addition to studying the bulk properties of the NH WSM, we also investigate its boundary behavior. Interestingly, instead of finding Fermi arcs, we find that the boundaries feature exceptional points (EP) with an order scaling with system size, where EPs are degeneracies in the complex spectrum at which both the eigenvalues and the eigenvectors coalesce, thus rendering the Hamiltonian matrix defective at these points, where the order of an EP is set by the number of eigenvectors that collapse [32]. At the high-order EPs in the model studied here, all bulk states are found to collapse and localize to either one of the boundaries. The appearance of these EPs is in line with the recent proposal that a non-trivial point-gap topology is linked to the NH skin effect [33,34], which refers to the piling up of bulk states at the boundaries, and the NH skin effect is in turn linked to the appearance of or vicinity to EPs of this kind [35,36]. We nevertheless do not observe a piling up of bulk states away from the EPs, which may be attributed to the fact that the EPs are well isolated from the rest of the spectrum. This paper is organized as follows: In Section II we briefly summarize the four-dimensional QHE, which is followed by the derivation of the NH WSM model and a discussions of its properties in Section III. We conclude with a discussion in Section IV.
The topological invariant associated with this model is the second Chern number C 2 where the integral is taken over the Brillouin zone (BZ), abcde is the antisymmetric Levi-Civita symbol, a, b, c, d, e ∈ {1, 2, 3, 4, 5},d i ≡ d i /|d|, and leads to [11]  In agreement with the Bloch spectrum, this spectrum is doubly degenerate such that each band-gap crossing corresponds to a doubly degenerate Weyl cone. The color of the bands with the color bar shown in (b) corresponds to the inverse participation ratio (IPR), In = 4N j=1 |Ψn,j| 4 /( 4N j=1 |Ψn,j| 2 ) 2 , where Ψn,j is the eigenfunction associated with energy band En, j is the site index and there are a total of 4N sites. The IPR measures the localization of the eigenstates: It goes to zero for extended states, whereas it acquires a finite value for localized states.
The value of the second Chern number C 2 is related to the presence of three-dimensional Weyl cones on the boundaries of the four-dimensional lattice model, where the number of boundary Weyl cones corresponds to |C 2 | [11]. Indeed, when we take open boundary conditions (OBC), we find Weyl cones at the high-symmetry points of the three-dimensional boundary BZ as shown in Fig. 1. Unlike the two-dimensional QHE, the four-dimensional QHE may be time-reversal symmetric [37], and indeed we find that T H * QHE (−k)T −1 = H QHE (k) with T = iΓ 3 and T T * = −1. The model in Eq. (1) thus belongs to symmetry class AII [38][39][40][41].

III. NON-HERMITIAN WEYL SEMIMETAL
Here we introduce a three-dimensional NH model, which is topologically equivalent to the four-dimensional QHE and realizes Weyl cones with a nonzero total chirality in the positive imaginary plane of the spectrum. We start by reexpressing the Hamiltonian in Eq. (1) as where d i (k) are the components of the vector d(k) defined in Eq. (2).H(k) has a chiral structure, i.e., , and anticommutes with Γ 5 , i.e., {H(k), Γ 5 } = 0. The term proportional to Γ 5 can thus be interpreted as a mass term. Therefore, the gap in the energy spectrum of H QHE (k) closes when there is a gap closing in the spectra ofH(k) and sin k w Γ 5 simultaneously. This means that the topology of H QHE (k) is captured byH(k, k w = 0) andH(k, k w = π). Due to the chiral symmetry (given by Γ 5 ) ofH(k, k w = 0, π), we can compute a winding number to characterize the topological properties of these Hamiltonians [40,41]. This means that we need to compute the winding number for one of the NH off-diagonal components ofH(k, . We find that equivalent results can be derived for all these four NH choices: The term h(k) · σ, which determines the gap closings of these NH Hamiltonians, is Hermitian and independent of k w , which means that all four possible choices of NH Hamiltonians are topologically equivalent. Indeed, the only difference between the four Hamiltonians is a constant and/or overall minus sign that goes with the identity term. Therefore, we may focus on the topological characterization of the following NH Hamiltonian This model is time-reversal symmetric according to H NH (k) = T H T NH (−k)T −1 with T = iσ y and T T * = −1, and as such belongs to class AII † [33,42]. This is in line with the statement in Ref. 30 that a d + 1-dimensional Hermitian model in class s is characterized by a ddimensional NH model in class s † . The topological properties of the four-dimensional QHE in Eqs. (1) and (2) are thus captured by the three-dimensional NH model in Eq. (6).
The eigenvalues of H NH (k) are E ± (k) = i(h + i cos k i ) ± |h(k)|, such that gap closings, E + (k) = E − (k), appear when |h(k)| = 0 is satisfied, i.e., when k = (0, 0, 0), P [(π, 0, 0)], P [(π, π, 0)] and (π, π, π). Despite the abundance of second-order exceptional lines in the complex spectra of three-dimensional NH models [43,44], these eight gap closings correspond to "ordinary" three-dimensional Weyl cones, and are degenerate in the usual sense, i.e., the geometric and algebraic multiplicity of the gap closings are equivalent. Indeed, the vector h(k), which determines the occurrence as well as the nature of the degeneracies, is entirely real, i.e., h(k) ∈ R 3 , and as such prohibits the appearance of exceptional structures in the complex spectrum of H NH (k). The Hamiltonian H NH (k) thus describes a three-dimensional NH Weyl semimetal. Due to the conventional nature of the gap closings, we can ascribe a chirality χ to each Weyl cone in a similar fashion as in the Hermitian case: By performing a low-energy approximation around the Weyl points, we obtain an effective Hamiltonian of the form M k · σ with k = (k x , k y , k z ), where the determinant of the matrix M determines the chirality, detM = χ. This leads to the following energies and chiralities at the eight degeneracies We observe that the three Weyl cones at k = P [(π, 0, 0)] and at k = P [(π, π, 0)] have the same energy, such that we find four band-gap closings in the spectrum with chiralities 1, −1, 3 and −3 as shown in Fig. 2 where the integral is over the BZ, µνσ is the Levita- The value of the winding number can be understood from the energy eigenvalues: When |h| > 3, the energy bands do not wind around E = 0 such that w 3D = 0. When |h| < 3, on the other hand, the bands wind around E = 0 and w 3D acquires a nonzero value. In fact, we notice that the value of the winding number corresponds to the difference between the total negative chirality, χ ImE>0 − , and the total positive chirality, χ ImE>0 + , of the cones with positive imaginary energy, i.e., w 3D = χ ImE>0 − − χ ImE>0 + , and as such is indirectly related to the values of C 2 . For example, when −1 < h < 1, there are three cones with negative chirality, χ ImE>0 − = 3, and one cone with positive chirality, χ ImE>0 + = 1, with Im(E) > 0, such that w 3D = 3 − 1 = 2 in this energy region. This result is in analogy with the result in Ref. 30, where the one-dimensional winding number is found to correspond to the difference between the number of left-and rightmoving modes with positive imaginary energy appearing in the one-dimensional Hatano-Nelson model [45]. We thus find that the value of the three-dimensional winding number w 3D corresponds to the total chirality of the Weyl cones surviving in the long-time limit, such that a nonzero w 3D signals the realization of anomalous Hermitian boundary physics in the bulk on the NH model.
As the three-dimensional NH model realizes Weyl cones in the bulk of its spectrum in the same fashion as in the Hermitian case, we naively expect the appearance of Fermi arcs on the surfaces connecting the cones. In Fig. 2(b), we plot the spectrum for OBC in x, and see that instead of Fermi arcs EPs (in red) appear, whose orders scale with system size. These EPs appear when k y , k z ∈ {0, π}, and have an algebraic multiplicity of 2N , with N the total number of unit cells, while the geometric multiplicity remains 2 regardless of the system size (see Appendix B for an analytical derivation of these results). The appearance of these EPs, where all eigenstates collapse onto two eigenstates of which one is localized entirely on one boundary and the other entirely on the other (see Appendix B), is consistent with the state-ment in Refs. 33, 34, where it is shown that having a topologically non-trivial point gap results in the appearance of skin states, i.e., bulk states that are localized to the boundary. However, as the EPs are well isolated from the rest of the bands in the spectrum, the states away from the EPs do not pile up and behave as ordinary bulk states (see Appendix B for the derivation of the localization of the bulk states), as is also found in Ref. 33, where a two-dimensional version of our model is studied. This is consistent with the logic presented in Refs. 35, 36, where the presence of higher-order EPs is linked to the NH skin effect. Indeed, only when the other states can "feel" the presence of such EPs, will they also start to pile up at the boundaries. As such, this model provides an interesting insight in the interplay of the presence of EPs with high order and the piling up of bulk states.

IV. DISCUSSION
In this work, we have studied a three-dimensional analogue of the four-dimensional QHE in the form of an NH WSM. The NH model in Eq. (6) is special in the sense that its gap closings can never be exceptional as h(k) ∈ R 3 as discussed above. By introducing NH terms into the vector h(k), however, it is straightforward to find exceptional structures as shown in Fig. 3.
NH Hamiltonians are well suited to describe classical systems with dissipation, e.g., [46][47][48], but also find applications in the quantum realm, e.g., [49][50][51][52], and it may thus be possible to access the NH model in Eq. (6) in experiment. Our model could for example be realized in a three-dimensional photonic crystals, where gain and loss can be introduced through optical pumping [53][54][55]. A connection between d + 1-dimensional Hermitian topological insulators and d-dimensional NH models with non-trivial point gap topology was recently also made on an abstract level by Lee et al. in Ref. 30. As an illustration, they present a one-and two-dimensional NH examples, which emulate the edge and surface physics of a two-dimensional Chern insulator and three-dimensional chiral topological insulator, respectively, which is complemented by the higher-dimensional model treated in this work.
We note that the exotic feature displayed by the NH WSM under OBC, namely, the appearance of EPs with an order scaling with system size, has previously been associated with the breaking of bulk-boundary correspondence [35,36], as well as an extreme spectral instability against boundary conditions [35,36,[56][57][58][59][60]. It may thus be tempting to attribute the same properties to the model in this paper. The models studied in these papers, however, all have so-called line gaps, where line gaps are a straightforward generalization of the gap in Hermitian spectra with the bands not crossing a line [31], whereas the model in Eq. (6) only features point gaps. Indeed, comparing Figs. 2(a) and (b), we see that there is no spectral instability, and we may thus conclude that these phenomena-a breaking of bulk-boundary correspondence, the NH skin effect and associated higherorder EPs, and the spectral instability-do not appear as a triad for models with only point gaps in the spectrum.
Lastly, we point out that due to the pivotal role played by time in obtaining an anomalous Weyl-cone configuration in the bulk spectrum of our model, we strictly speaking need four dimensions, namely three spatial and one time dimension, to realize our NH analogue of the four-dimensional QHE. Therefore, while it is thus possible to mimic the boundary physics of d + 1-dimensional Hermitian models in d-dimensional NH models, where the dimensions refer to spatial dimensions, one crucially needs d + 1 spacetime dimensions to obtain the desired bulk behavior.

Acknowledgments
We thank Emil J. Bergholtz, Jan Carl Budich and Loïc Herviou for useful conversations. We also thank Patrick Emonts for assisting in obtaining analytical results for the OBC spectrum for larger system sizes. F.K.

Appendix A: Numerical instability at exceptional points
Here we compare the analytically and numerically computed spectra of the three-dimensional NH WSM Hamiltonian with OBC in x. In Fig. 4(a), we plot the OBC spectrum with N = 14 unit cells with the analytical results in orange and the numerical results in blue, such that the orange (analytical) spectrum is equivalent to the spectrum in Fig. 2(b). Whereas the results for the bulk spectrum are in good agreement with each other up to machine precision-the difference of the real part of the numerically and analytically obtained eigenvalues is of order ∼ 10 −15 whereas the difference of the imaginary part of the numerical and analytical results ranges from ∼ 10 −15 for most of the eigenvalues to ∼ 10 −4 -, the results for the in-gap states, i.e., the EPs, deviate significantly hinting at a numerical instability at these points. This numerical instability becomes even larger upon increasing the system size as shown in Fig. 4(b) for N = 60 unit cells, where the deviation (in blue) away from the EPs (orange) is significant. We point out that the numerical eigenvalues are obtained using the standard eigensolver in Mathematica, and that changing the eigensolver, decreasing the tolerance and increasing the maximum iteration does not alter the numerical results. Indeed, numerical computations performed in Python ( geev LAPACK routines) yield the same results. This failure of standard numerical techniques to correctly find the spectrum is due to the extreme defectiveness of the matrix Hamiltonian at the EPs, i.e., there are only two eigenvectors instead of 2N , which prevents the convergence of numerical values. Due to the apparent correspondence of the analytically and numerically computed bulk bands [cf. Fig. 4(a)], we plot the bulk bands numerically (in blue) and the EPs analytically (orange) in Fig. 4(c) for a large system (N = 60), and see that the bulk gap remains open with the increase in system size. We stress that the results in Fig. 4 tell a cautionary tale, and one needs to be extremely careful when performing numerical computations on NH matrices that have a large defectiveness. Indeed, for such matrices, analytical checks should be made.

Appendix B: Analytical computation of eigensystem
In this Appendix, we compute eigenvalues and eigenvectors for the NH WSM with OBC in x.
Eigenvalues. The Hamiltonian with OBC in x and a total of N unit cells reads (c 1,A , c 1,B , c 2,A , c 2,B , . . . , c N,A , c N,B ) T with A, B the degrees of freedom and c † x,j (c x,j ) creating (annihilating) a state with degree of freedom j in unit cell x, such that we find where with h 0 = h + cos k y + cos k z , and h 2 and h 3 defined in Eq. (6), i.e., h 2 = sin k y and h 3 = sin k z . The eigenvalues λ for this Hamiltonian are found by solving the characteristic equation to reduce this eigenvalue equation to where we identify A =H N −1,x , D =H 1,x , and Making use of the specific form of the matrices, we find and We note that the matrixH N −1 ⊥,x has nonzero entries only in its lower right 2 × 2 block, such that subtracting it fromH N −1,x only changes the lower right block, i.e., is again only non-zero in the bottom right corner. Repeating this step N times, and defining we finally arrive at the expression We thus find that the determinant of the 2N -dimensional OBC Hamiltonian can be expressed in terms of a product of N determinants of two-dimensional matrices. We note that each determinant in this product can ultimately be expressed in terms of the (inverse of the) determinant of H 1,x -each a i depends on a −1 1 = adj a 1 det −1 a 1 -, such that it is not possible to find eigenvalue solutions by simply setting each determinant, det a i , to zero separately. Instead, one should explicitly reduce this equation by making use of the explicit form of a i to obtain a polynomial equation of order 2N . Interestingly, as we will see in the following, an elegant and straightforward solution can be derived for some special points in the spectrum.
Eigenvectors at EPs. Next, we turn to the eigenvectors at these k points, i.e., where Φ = (φ A,1 , φ B,1 , . . . , φ A,N , φ B,N ) T and X = (χ A,1 , χ B,1 , . . . , χ A,N , χ B,N ) are the right and left eigenvectors, respectively, with the index A, B referring to the degrees of freedom inside the unit cell as before.
Using that H 1,x is an identity matrix with ih 0 = λ on the identity for k y , k z ∈ {0, π}, we find for the right eigenvectors and for the two boundaries, and for the bulk. From these equalities, we obtain the following for n ∈ 2, . . . , N − 1, such that we immediately find The Hamiltonian matrix H N,x thus has two linearly independent eigenvector solutions when k y , k z ∈ {0, π}, namely, Φ 1 = (1, 1, 0, . . . , 0) T , and Φ 2 = (0, . . . , 0, −1, 1) T . When repeating the same exercise for the left eigenvectors, we find for the two boundaries, and for the bulk. Therefore, following the same steps as before, we arrive at χ A,n = χ B,n = 0, ∀n ∈ {2, . . . , N − 1}, such that we again obtain two linearly independent solutions, i.e., X 1 = (1, −1, 0, . . . , 0) T , and X 2 = (0, . . . , 0, 1, 1) T . These left eigenvector are selforthogonal to the right eigenvectors, i.e., X|Φ = 0, and the degenerate eigenvalues with energy λ = ih 0 thus correspond to exceptional points with order 2N − 2. This piling up of states can also be understood from studying the form of the OBC Hamiltonian at the EPs in more detail, where this argument is inspired by the argument in section SVIII in Ref. 33. In the following, we perform a unitary transformation such that σ x → σ z in H NH (k) in Eq. (6), and find the following Hamiltonian when taking OBC in x where we use the same conventions for c ( †) x,j as before. Firstly, we note that at the EPs, i.e., h 2 = h 3 = 0, this reduces to where λ is the eigenvalue of the EP as before. This Hamiltonian describes two decoupled chains, where there is only hopping to the left in the "A chain" and only hopping to the right in the "B chain" thus explaining the appearance of the EPs and the existence of only two eigenvectors. This model thus realizes two decoupled Hatano-Nelson chains [45] in the extreme limit, i.e., hopping in one direction only. Secondly, away from the EPs, H N,x can be interpreted as describing two coupled Hatano-Nelson chains, which still reside in the extreme limit. Therefore, one would naively expect that when h 2 and h 3 are small, the eigenstates still pile up albeit now with different eigenvalues. Interestingly, however, both from numerical checks as well as analytical computations (see below), we find that this is not the case. This is in agreement with the result in Ref. 61, where it is shown that the presence of a magnetic field coupling the two (pseudo)spin sectors, i.e., coupling A and B to each other, suppresses the NH skin effect. Thirdly, we stress that performing the unitary transformation (σ x → σ z ) does not alter the physics of our model upon considering OBC. Indeed, it is possible to show that U N = 1 N ⊗ U , where U † (σ x , σ y , σ z )U = (σ z , σ y , −σ x ) with U = (σ 0 −iσ y )/ √ 2, is a unitary matrix, which transforms H N,x into H N,x according to U † N H N,x U N = H N,x . This means that the spectra of H N,x and H N,x are equivalent, and U † N Φ and XU N are right and left eigenvectors of H N,x , respectively. Due to its specific form, U N only acts locally in the eigenvectors and thus does not change their overall behavior. Lastly, as one would expect, we find that H N,x equals the Hamiltonian for OBC in z, H N,z , with k x → k z in the latter. Indeed, all results in this paper are independent of whether one takes OBC in x, y or z as the Hamiltonian is equivalent in all three directions.
Eigenvectors away from EPs. Here, we follow the method presented in Ref. 62 to find an explicit form for the eigenvectors away from the EPs. We start by rewriting the equalities for the right eigenvectors in Eqs. (B7)-(B12) for general k, which yields the following and for the two boundaries, and for the bulk. Next, we make the following ansatz for the eigenfunction in unit cell n [62] Φ n = φ A,n φ B,n = j r n j φ (j) A φ (j) B .