Magnetic Flux Response of Non-Hermitian Topological Phases

We derive the response of non-Hermitian topological phases with intrinsic point gap topology to localized magnetic flux insertions. In two spatial dimensions, we identify the necessary and sufficient conditions for a flux skin effect that localizes an extensive number of in-gap modes at a flux core. In three dimensions, we furthermore establish the existence of: a flux spectral jump, where flux tube insertion fills up the entire point gap only at a single parallel crystal momentum; a higher-order flux skin effect, which occurs at the ends of flux tubes in presence of pseudo-inversion symmetry; and a flux Majorana mode that represents a spectrally isolated mid-gap state in the complex energy plane. We uniquely associate each non-Hermitian symmetry class with intrinsic point gap topology with one of these cases or a trivial flux response, and discuss possible experimental realizations.


I. INTRODUCTION
Gapped Hermitian topological phases can be differentiated from trivial phases, as well as between each other, by a quantized response to local magnetic fluxes [1][2][3][4][5][6][7][8][9]. Such fluxes take the form of zero-dimensional (0D) flux cores (vortices) in two dimensions (2D), or one-dimensional (1D) flux tubes in three dimensions (3D). In particular, in presence of a symmetry that quantizes the magnetic flux ϕ through a plaquette (stack of plaquettes) of a 2D (3D) lattice -most commonly to values ϕ = 0, π -topological insulators generically host flux-localized bound states [10]. In 2D, these states may be constrained to occur at zero energy by a spectral symmetry such as the particle-hole symmetry of (mean-field) superconductors. For example, magnetic vortices in a 2D p + ip superconductor are known to host unpaired Majorana zeromodes [11,12]. Without spectral symmetry, flux bound states in 2D systems can be moved out of the gap, but may still contribute to a filling anomaly of the ground state that cannot be trivialized without breaking a symmetry or closing a gap [13]. In 3D, flux bound states fall into two categories: in the first case, they form a gapless state localized along the 1D flux tube, as is the case for π-flux tubes in 3D time-reversal symmetric topological insulators [3][4][5][6][7][8][9]. In the recently discovered second case, the flux tube ends bind 0D states that give rise to a filling anomaly or are pinned to zero energy by a spectral symmetry [14]. This latter case is realized in higher-order topological phases protected by crystalline symmetries.
The goal of our present work is to generalize the theory of flux bound states to non-Hermitian (NH) topological phases that break energy conservation. Such phases have generated increased interest recently due to their unconventional bulk-boundary correspondence [15][16][17][18][19][20][21][22][23][24][25]. Fundamentally, there exist two kinds of NH systems [15][16][17][18]: Hamiltonians H with a line gap in their complex energy spectrum can be adiabatically (without closing the line gap) deformed to purely Hermitian (H † = H) or anti-Hermitian (H † = −H) Hamiltonians. On the other hand, Hamiltonians without any line gap may still host a point gap -a region of complex energy that is devoid of, but surrounded by, eigenstates. Point-gapped systems without a line gap cannot be adiabatically deformed to any (anti-)Hermitian limit, and therefore realize intrinsically NH topology. A prime example is the 1D Hatano-Nelson chain [26][27][28] in NH symmetry class A (no symmetries) whose bulk winding number invariant W ∈ Z in periodic boundary conditions (PBC) results in the NH skin effect under open boundary conditions (OBC): when W ̸ = 0, an extensive number of (almost all) OBC eigenstates accumulate at only one edge of the system [29][30][31][32][33][34][35][36][37]. While this observation is reminiscent of a bulk-boundary correspondence, it is unclear how the OBC spectrum may crisply differentiate between systems with different nonzero W . To alleviate this issue, we here employ the NH pseudospectrum that consists of the collection of spectra associated with all O(ϵ)deformed Hamiltonians (see Sec. II A for details) [38][39][40]. The pseudospectrum reduces to the spectrum in the Hermitian case, but represents the more physically adequate quantity in the NH case: if a state fails to be an eigenstate only by O(ϵ) terms, then it behaves as an eigenstate with respect to realistic measurements. Moreover, the pseudospectrum restores the full bulk-boundary correspondence of NH point-gapped systems: for instance, a Hamiltonian with nonzero bulk winding number W has a pseudospectrum that fills the point gap and is W -fold degenerate in presence of boundaries [39]. In the thermodynamic limit, the pseudospectrum is equivalent to the spectrum in semi-infinite boundary conditions (SIBC) [38][39][40]. Correspondingly, we characterize the flux response of NH systems by their SIBC spectrum in presence of flux defects. We exclusively focus on NH phases with intrinsic point gap topology, which were classified in all 38 NH symmetry classes in Refs. 16 and 39: in addition to stability under manipulations preserving point gap and symmetry class, these phases cannot be trivialized even when coupling with arbitrary NH linegapped phases is allowed.
We find that, depending on the NH symmetry class, the flux response of intrinsically NH point-gapped phases falls into one of 5 classes: arXiv:2208.11712v2 [cond-mat.mes-hall] 1 Sep 2023 (1) No flux response: Some NH phases have a trivial flux response even when their bulk is topologically nontrivial.
(2) Flux skin effect: Inserting a π-flux core in a 2D NH system can induce a skin effect response where a macroscopic number of OBC eigenstates localizes at the flux core [39]. Concomitantly, the SIBC spectrum in presence of flux fills up the entire point gap in the complex energy plane. This effect is similar to the dislocation skin effect of Refs. 41 and 42, with the important distinction that it does not rely on (discrete) translational symmetry and therefore probes strong instead of weak NH topology [14].
(3) Flux spectral jump: Threading a 1D π-flux tube through a 3D NH system with nontrivial point gap topology can result in a spectral jump as the momentum coordinate k ∥ along the flux tube is varied. In particular, the SIBC spectrum in presence of flux remains gapped at generic values of k ∥ , but completely fills up with an extensive number of states at one of the two special points k ∥ = 0, π. Such a gapless NH dispersion cannot be realized in any purely 1D lattice system, and therefore represents an intrinsically NH anomalous 1D state. We note that a purely 1D system cannot realize a filled point gap at a single momentum, as there is only a finite Hilbert space available at any given momentum. Instead, the discontinuity in the 1D flux tube dispersion described here capitalizes on a topologically nontrivial 2D bulk.
(4) Higher-order flux skin effect: Even when the SIBC spectrum of a 3D NH phase in presence of a 1D flux tube is fully gapped (does not exhibit a spectral jump), preserving crystalline pseudo-inversion symmetry [43][44][45][46][47][48][49][50][51] may result in a nontrivial higherorder response. In this case, an extensive number of eigenstates accumulates at the ends of flux tubes when these terminate at sample surfaces. Concomitantly, the SIBC spectral point gap remains empty with PBC along the flux tubes, but is completely filled up in OBC when a surface termination perpendicular to the flux tubes is introduced.
(5) Flux Majorana modes: In just two NH symmetry classes, intrinsic 3D point gap topology furthermore manifests itself in an unpaired flux Majorana mode pinned to the center of the SIBC spectral point gap.
For the case of NH class D, this mode is localized along the entire 1D length of the flux tube. Such localization is impossible in Hermitian systems where energetically isolated bound states must be pointlike.
In Tab. I, we identify the flux response of each NH symmetry class that allows for intrinsic point gap topology with one of the 5 possible NH flux responses. In Sec. II we outline the formalism that we use to derive the flux responses, before discussing the individual effects: in Sec. III (1)], induces a spectral collapse of the point gap under OBC. However, the resulting OBC spectrum does not uniquely specify the nonzero value of W (E). This ambiguity is resolved in the SIBC spectrum, which restores the NH bulk-boundary correspondence: surface states fill the point gap, and their degeneracy corresponds to W (E).
the flux spectral jump, in Sec. V we introduce the higherorder flux skin effect, and finally Sec. VI unveils the flux Majorana modes. We close by discussing potential experimental realizations in Sec. VII. The appendices contain exhaustive auxiliary derivations and toy model Hamiltonians.

II. FORMALISM
We begin by pedagogically reviewing the concepts needed to characterize flux responses in the NH context.
A. NH Bulk-boundary correspondence NH systems can differ dramatically in their spectra under PBC and OBC. One of the prime examples of this feature is the NH skin effect, whereby a point-gapped bulk collapses to a line under OBC (see Fig. 1). Connected with this is a pile-up of all states at a single boundary. In 1D, this effect is associated with a winding number W (E) [39], where the Bloch Hamiltonian H(k) is point-gapped about the reference energy E. In particular, if W (E) ̸ = 0, the NH skin effect must occur. However, it is important to note that W (E) ∈ Z is an integer-valued invariant, while the presence or absence of a skin effect only provides a Z 2 quantifier. This observation raises the question of how the OBC spectrum of two systems with nonzero but different W (E) can be crisply distinguished. The problem is further compounded by the dramatic sensitivity of the OBC spectrum [52] to small changes of the Hamiltonian.
To solve both issues, the ϵ-pseudospectrum of NH systems was introduced in Refs. [38][39][40]. It describes the change in the spectrum under a small perturbation ϵ, σ ϵ (H) = {E ∈ C| ||(H − E) |v⟩ || < ϵ for at least one |v⟩ with ⟨v|v⟩ = 1}. (2) This definition is practical: if a state is just O(ϵ) away from being an eigenstate, it will still behave as one with regards to realistic measurements, which always include a small error. In contrast to other approaches reestablishing the NH bulk-boundary correspondence, for instance the OBC treatment of Ref. 53, the pseudospectrum, therefore, does not suffer from a sensitivity to infinitesimal errors. The pseudospectrum is particularly useful because it can be related to the spectrum in SIBC σ SIBC (H), which is the spectrum in presence of just a single boundary in the thermodynamic limit (system size L → ∞). In particular, it holds that [38][39][40] This correspondence between pseudospectrum and SIBC spectrum allows for a precise definition of a NH bulk-boundary correspondence. In particular, in the SIBC spectrum of a 1D system with nonzero W (E), the point gap fills completely with boundary localized states whose degeneracy equals W (E) (see Fig. 1).

B. Extended Hermitian Hamiltonian
To find the flux response of NH systems, we rely on the topological equivalence between a NH Hamiltonian H and an extended Hermitian Hamiltonian (EHH)H [54,55], defined asH By construction,H enjoys a chiral symmetry, The presence of a point gap of H around E = E 0 translates into a gapped spectrum ofH about zero energy. Conversely, exact topological zero-energy eigenvalues of H correspond to protected E = E 0 states within the NH point gap, because Consequently, we can predict the presence of protected in-gap modes from the EHH spectrum.
Any topological flux response relies on a quantization of the admissible values of flux ϕ. Restricting to local NH symmetries corresponding to one of the 38 NH symmetry classes [16] -potentially taken together with a crystalline symmetry like pseudo-inversion -we find that, depending on the symmetry class, flux is either unquantized or must take values ϕ = 0, π in PBC. In the cases where ϕ = 0, π, the NH symmetry class must either contain a time-reversal (TRS) or particle-hole (PHS) symmetry, or include crystalline pseudo-inversion symmetry.
Here, we focus on PBC in order to cleanly separate the NH flux response from boundary states or skin effects. The PBC geometry requires (at least) two flux cores/tubes with strength ±ϕ [14]. Such a pair of fluxes can be introduced by the Peierls substitution: the hoppings encircling each flux must accumulate a phase e ±iϕ (see Fig. 2a). A convenient electromagnetic gauge choice is then to multiply all hoppings across the line (plane) connecting the two fluxes in a 2D (3D) system by a factor of e iϕ in one direction (s 1 → s 2 , see Fig. 2b) and by e −iϕ in the other direction (s 2 → s 1 ), where we denote the sites above the line (plane) as s 1 and the ones beneath as s 2 . We choose the orientation of this line (plane) along the x-(x-and z-)direction in 2D (3D), respectively. SIBC corresponds to the idealized limit where only one flux (= only one "boundary") is present, while the second flux is infinitely far away. In our conventions, a system with SIBC then has a single flux core/tube at x = 0, is infinitely extended in the x-direction, and is finite with PBC in all remaining directions (but potentially still thermodynamically large).
We derive the flux response for all NH symmetry classes with intrinsic point gap topology from their respective EHH. The EHH experiences the same flux defect as the NH Hamiltonian: in the gauge described above, The energy spectrum of the EHH for NH class AII † [Eq. (E1)] shows four topological zero energy modes within the bulk energy gap, caused by two π-flux cores. c The flux skin effect localizes an extensive number of eigenmodes at the fluxes, indicated by two peaks of the summed density ρ(r) = α,i |⟨ri|ψα⟩| 2 , where α ranges over all eigenstates ψα of the Hamiltonian with flux defect, r denotes the lattice site, and the summation i runs over sublattice degrees of freedom.

Peierls substitution implies
The corresponding matrix element for the EHH reads ⟨s α 1 |H ϕ |s β 2 ⟩, where α, β = 1, 2 label the two sublattices introduced by the Hermitian extension. The only nonzero contributions are and ⟨s α=2 where, importantly, the flux enters in both cases as e iϕ . The reason is that Hermitian conjugation in Eq. (9) not only flips the sign of the flux, but also transposes the matrix elements s 1 → s 2 to s 2 → s 1 . In Eq. (7) we have used that E 0 multiplies an identity matrix in Eq. (4) and therefore remains unaffected under flux insertion. In summary, our strategy for determining the flux response of a given NH symmetry class X is to: (1) Find the corresponding Hermitian symmetry classX of the EHH.
(2) Derive the flux response ofX using the Dirac formalism detailed in the appendices.
(3) Infer the topological properties of the NH SIBC spectrum that result from exact EHH zeromodes.

III. NH FLUX SKIN EFFECT
We begin our survey of NH flux responses with 2D systems. Here we first consider a specific example in NH class AII † (see App. E 1 a for details). Point-gapped systems in NH class AII † are classified by a Z 2 invariant ν ∈ {0, 1}, defined in Ref. 16. In order to derive the flux response of the nontrivial phase where ν = 1, we rely on the EHH for a given energy E 0 inside the point gap (see Sec. II B). Irrespective of the choice of E 0 , the EHH for NH class AII † enjoys the symmetries of Hermitian class DIII (see App. B 4 a and Tab. II therein). Moreover, since Hermitian class DIII is likewise Z 2 -classified in 2D, the EHH associated with a nontrivial NH phase in class AII † must itself realize a nontrivial topological insulator phase in Hermitian class DIII. For this phase, it was shown in Ref. 10 that flux cores bind two degenerate zero-energy states (a Kramers pair). In the NH SIBC spectrum, these modes then correspond to a single flux localized state at complex energy E 0 in the point gap [39]. Since we can perform this construction for all E 0 inside the point gap, we obtain an extensive number of modes localized at the flux core (see Fig. 3a,b). Such an extensive accumulation of states defines the NH flux skin effect [39].
In a PBC geometry with two flux cores, this response is topologically equivalent to the Z 2 skin effect [39] of a 1D model in NH class AII † situated on the line terminated by the two flux cores (see Fig. 3c). Indeed, the number of flux-localized modes scales with the length of the 1D edge connecting the two defects, denoted by L ⊥ (see App. C 1 for details on the finite-size scaling). The equivalence of the flux skin effect and the 1D NH bulk-boundary correspondence is consistent with the fact that point gap topology in 1D and NH class AII † is Z 2 -classified [16]. This observation leads us to the following generalizing hypothesis, which we prove in App. B by exhaustion: (I) A 2D system with nontrivial intrinsic point gap topology in NH symmetry class X, where X pins ϕ = 0, π, exhibits a flux skin effect for ϕ = π iff X also has a 1 6 -1 Flux-localized states fill the SIBC spectral point gap at a single momentum k ∥ = 0 along the flux tube. b At this momentum, the EHH spectrum exhibits a helical crossing at zero energy for any E0 located in the central point gap. c The flux spectral jump localizes an extensive number of eigenmodes at the flux tubes, indicated by two peaks of the summed density ρ(r) = α,i |⟨ri|ψα⟩| 2 , where α ranges over all eigenstates ψα of the Hamiltonian with flux defect, r denotes the lattice site and the summation i runs over sublattice degrees of freedom.
nontrivial point gap classification in 1D.
The flux response of the complete set of NH symmetry classes with nontrivial intrinsic point gap topology [39] is summarized in Tab. I. Only the classes satisfying the above condition, including NH class AII † , realize a flux skin effect. An exhaustive collection of model Hamiltonians can be found in App. E 1.
We note in closing that there is a fundamental difference between crystal dislocations and flux defects [14]: Whereas lattice Burgers vectors are integer-valued, flux is only defined mod 2π. Consequently, since ϕ = +π and ϕ = −π are physically equivalent, both fluxes ϕ = ±π of a PBC geometry have to collect skin modes, which is compatible only with a Z 2 -skin effect but not a Z-skin effect [39]. In contrast, dislocations can give rise to both the Z and the Z 2 -skin effect [41].

IV. NH FLUX SPECTRAL JUMP
We now study the π-flux response of intrinsically NH topological phases in 3D [39] in presence of TRS ( †) or PHS ( †) . We begin with the example of a point-gapped NH system in class AII † , which is classified by a nontrivial 3D winding number W 3D (E) ∈ Z [16] (see App. E 2 a for details). We derive the flux response of the simplest nontrivial case, where W 3D (E 0 ) = 1 for a given energy E 0 inside the point gap, from the corresponding EHH (see Sec. II B and Fig. 4a). Irrespective of the choice of E 0 , this EHH enjoys the symmetries of Hermitian class DIII (see App. B 5 b and Tab. III therein), which is also Z-classified in 3D. Hence, the EHH associated with W 3D (E) = 1 corresponds to a nontrivial topological insulator phase [59]. For this phase, it was shown in Ref. 10 that flux tubes with ϕ = π bind a 1D helical Majorana mode. This mode crosses zero energy at a single fixed TRS-invariant momentum, k ∥ ∈ {0, π}, contributing two exact zero-energy states (see Fig. 4b). Importantly, no symmetry allowed perturbation is able to move these flux localized states away from zero energy at k ∥ = 0, π. Correspondingly, there is a single flux-localized state at complex energy E 0 in the NH SIBC spectrum [39]. Repeating this construction for all E 0 inside the NH point gap yields an extensive number of states localized along the 1D flux tube for a single momentum k ∥ ∈ {0, π} (see Fig. 4c). Away from this momentum, the EHH remains gapped, corresponding to the absence of in-gap modes in the NH SIBC spectrum. The extensive pile-up of flux-localized states only for a single value of k ∥ constitutes a novel type of NH bulk-defect correspondence, which we term NH flux spectral jump.
In a PBC geometry with two flux tubes, this response is equivalent to that of a nontrivial 2D NH phase situated in the plane separating the flux tubes: NH class AII † is Z 2 classified in 2D [16] and protects an infernal point, i.e. the extensive accumulation of edge states only at a single momentum [60]. The number of flux localized modes scales with L ⊥ , the distance between the two defects (see App. C 2 for details on the finite-size scaling). For OBC in the z−direction, where k ∥ is not conserved, the NH flux spectral jump still implies an accumulation of states at the flux tube. There may potentially also be surface states contributed by a surface skin effect. Next to their different real-space localization, the two effects can be distinguished in terms of the scaling with system size: whereas the surface skin effect is expected to localize O(L ∥ = L z ) modes, the flux spectral jump localizes O(L ⊥ ) modes.
We note that all NH phases with nontrivial point gap topology are Z 2 -classified in 2D. On the other hand, they may be Z 2 -or Z-classified in 3D. In the latter case, the flux spectral jump only arises when the Z-valued invariant is odd. In our example of NH class AII † , this means that only systems where W 3D (E) mod 2 ̸ = 0 exhibit a flux response (see Tab. I). Additional responses would be needed to probe the Z nature of W 3D (E). The above observations lead us to the following hypothesis, which we prove in App. B by exhaustion: (II) A 3D system with nontrivial intrinsic point gap topology in NH symmetry class X and invariant w ∈ Z 2 or Z, where X pins ϕ = 0, π, shows a flux spectral jump for ϕ = π iff class X has a nontrivial point gap classification in 2D and w mod 2 ̸ = 0.
The flux response of the complete set of NH symmetry classes with nontrivial intrinsic point gap topology [39] is summarized in Tab. I. Only classes fulfilling the above condition, including NH class AII † , realize a flux spectral jump. An exhaustive collection of model Hamiltonians can be found in App. E 2.

V. NH HIGHER-ORDER FLUX SKIN EFFECT
We next study the π-flux response of intrinsically NH topological phases in 3D [39] in the absence of TRS ( †) or PHS ( †) . In order to still obtain a quantized flux response, we consider the presence of pseudo-inversion symmetry (see also App. A 1), which, in contrast to normal inversion symmetry, entails a Hermitian conjugation of the Hamiltonian H(k).
As explained in App. D 2, pseudo-inversion represents the simplest crystalline symmetry that is compatible with nontrivial point-gap topology. In a PBC geometry with two inversion-related flux tubes, this symmetry implies that the magnitude of the two fluxes is the same: ϕ 1 = ϕ 2 . At the same time, PBC leads to the constraint ϕ 1 + ϕ 2 mod 2π = 0, thereby recovering the quantization of magnetic flux to ϕ = 0, π. This argument only requires PBC in the direction separating the flux tubes (i.e. in our case in the x-direction); in fact, we will consider OBC in the direction along the flux tubes later on [61]. We note that the addition of pseudo-inversion symmetry may change the classification of NH point gap topology. However, in this work, we are only concerned with the flux response of point-gapped systems with local NH symmetries that are enriched by pseudo-inversion (see App. D 2). We begin by considering the NH system discussed in Sec. IV and App. E 2 a in absence of TRS † . Recall that, in that section, we studied NH symmetry class AII † which exhibits a flux spectral jump. Upon the relaxation of TRS † , we obtain NH symmetry class A, which is still classified by the same topological invariant W 3D (E) ∈ Z [16] (see App. E 2 b for details). Our present example therefore again yields W 3D (E 0 ) = 1 for any given energy E 0 inside the point gap. As explained above, we furthermore assume pseudo-inversion symmetry.
We derive the flux response for this phase from the corresponding EHH (see Sec. II B). The EHH for NH class A retains the symmetries of Hermitian class AIII, irrespective of the choice of E 0 (see App. B 5 c and Tab. III therein). Additionally, NH pseudo-inversion symmetry induces a Hermitian inversion symmetry of the EHH. Hermitian class AIII is also Z-classified in 3D, such that the EHH associated with a nontrivial point-gapped phase in NH class A must itself realize a nontrivial topological insulator phase in Hermitian class AIII [59]. Importantly, this class does not protect gapless modes along 1D flux tubes [10]. Hence, the PBC spectrum of the EHH in presence of two π-flux tubes is gapped (see Fig. 5a), where here we assume PBC both along and perpendicular to the flux tube directions. Correspondingly, the NH point gap shows no SIBC spectral in-gap modes. However, inversion symmetry can protect a flux response as soon as OBC are introduced in the direction of the flux tubes.
To derive this response, we first note that Ref. 44 showed that W 3D (E 0 ) = 1 is indicated by a double band inversion in the corresponding EHH when pseudoinversion symmetry is present: to obtain a trivial (atomic limit) state, two occupied EHH inversion eigenvalues must flip sign at high-symmetry momenta in the 3D Brillouin zone. Without chiral symmetry (Hermitian class A), a double band inversion induces a higher-order topological (axion insulator) phase protected solely by inversion symmetry [62,63]. Such axion insulators exhibit a higher-order flux response: Ref. [14] derived that localized states appear at alternating ends of the two flux tubes, resulting in one state per flux tube for the EHH. In Hermitian class A, these states are generally not zeromodes, but instead contribute to a filling anomaly of the full inversion-symmetric system. Returning to the present case of Hermitian class AIII, chiral symmetry enforces that such flux end states arise at zero energy.
Besides pinning flux end states to zero energy, Hermitian class AIII differs in another important aspect from class A: it protects surface Dirac cones [59] that render the EHH OBC spectrum gapless even in absence of flux tubes (see Fig. 5b, left inset). However, the presence of such 2D surface states does not, in general, result in a skin effect in NH class A: The EHH surface Dirac cones are located at arbitrary surface momenta, which are in general not sampled over in the discrete surface Brillouin zone associated with any finite system size. Hence, there are no exact zero-energy states in the EHH spectrum [64]. Instead, each EHH Dirac cone implies a single exceptional point or single sheet of eigenvalues in the NH surface dispersion [21]. On the other hand, each flux tube introduces a single exact zero-mode per flux tube in the EHH. This constitutes a fractional flux response, since there exists no 1D Hermitian model with a single end state under OBC (see also App. D 1). The zero-energy flux end state of the EHH then corresponds to a fluxlocalized SIBC spectral in-gap state at complex energy E 0 . Repeating this construction for all E 0 inside the NH point gap yields an extensive number of modes localized at the end of a flux tube (see Fig. 5c). In a finite system, the number of localized modes scales with the length of the flux tube, denoted by L ∥ (see App. C 2 for details on the finite-size scaling). This extensive pile-up of states on  [16,39]. Class labels refer to the Altland-Zirnbauer (AZ) classes, their AZ † counterparts, and additional sublattice symmetry (SLS) [16]. The subscript of SLS, S±, determines whether SLS commutes (+) or anticommutes (-) with time-reversal symmetry (TRS) and/or particle-hole symmetry (PHS). For the NH symmetry classes with both TRS and PHS (BDI, DIII, CII, and CI), the first subscript refers to TRS and the second one to PHS. All classes are identified by the square of TRS ( †) , PHS ( †) and chiral symmetry CS. The topological classification is reproduced from Ref. 39, with d ℓ,1D denoting the 1D line gap classification for the cases with trivial point-gapped phases in both 1D and 2D. I † denotes the additional presence and form of pseudo-inversion symmetry that, in NH symmetry classes without an anti-unitary operation, is needed to protect a topological flux response. The subscript of SLS/TRS next to the pseudo-inversion symmetry I † , S±/T±, determines whether SLS/TRS commutes (+) or anticommutes (-) with I † . L ⊥ denotes the distance between the flux tubes in a finite system, and L ∥ their length. We distinguish the flux skin effect (FSE) in 2D, the flux spectral jump (FSJ), higher-order flux skin effect (HO-FSE), the flux Majorana mode (FMM) and the higher-order flux Majorana mode (HO-FMM) in 3D. dim class symmetry classification a region of dimension d − 3 represents a NH higher-order flux skin effect. Interestingly, the flux skin modes appear at the same end for both flux tubes, highlighting a uniquely NH feature (see Consequently, when the plane spanned between the two π-flux tubes is interpreted as a 1D NH system aligned along the flux tube direction, and with an extensively large unit cell along the perpendicular direction [65], the higher-order flux skin effect is equivalent to the skin effect of a nontrivial 1D system [39] in NH class A. Due to pseudo-inversion symmetry, which maps between the two flux tubes, the resulting extensive number of end-localized modes is evenly distributed between the two flux tube ends, so that each flux tube on its own realizes half of a conventional 1D skin effect. On the other hand, the 2D point-gap classification of NH symmetry class A is trivial. This observation leads us to the following hypothesis, which we prove in App. B by exhaustion: (III) A 3D system with nontrivial point gap topology in NH symmetry class X and with pseudo-inversion symmetry exhibits a higher-order flux skin effect for ϕ = π iff the point gap classification of X is trivial in 2D but nontrivial in 1D.
The flux response of the complete set of NH symmetry classes with nontrivial intrinsic point gap topology [39] is summarized in Tab. I. Only classes fulfilling the above condition, including NH class A, show a NH higher-order flux skin effect. An exhaustive collection of model Hamiltonians can be found in App.  5. NH higher-order flux skin effect in 3D. a The energy spectrum of the EHH for a system with higher-order flux skin effect shows a gap for PBC along the flux tubes, here using a model in NH class A [Eq. (E8)]. b Under OBC in the flux tube direction, the presence of a nontrivial flux ϕ = π introduces two exact zero-energy states in addition to the gapless Dirac cone surface dispersion, one for each flux tube. c In the NH system, the resulting flux skin modes are localized at the same ends of the two flux tubes. As explained in Sec. V, this is a consequence of pseudo-inversion symmetry I † , which maps states between the two flux tubes, but also changes their chirality. Only states with more than 85% support at the flux tubes are shown. All panels are generated for a model of size 20 × 20 × 20 unit cells in NH class A with additional pseudo-inversion symmetry (see Sec. E 2 b).
We note in closing that the higher-order flux skin effect realizes a novel flavour of NH topology that is fundamentally different from the Hermitian case: In a Hermitian system with gapless surface states, it is impossible to resolve a higher-order response to flux defects, because the nontrivial surface state would obscure any fluxlocalized modes. On the contrary the NH higher-order flux skin effect is observable: In addition to O(L ⊥ × L ⊥ ) NH surface modes deriving from gapless surface states, flux tubes induce a skin effect in z-direction, localizing another O(L ∥ ) modes on the surface, only at the flux tube ends (see App. C 2 for details on the finite-size scaling and App. D 3 for the implications on the EHH spectrum). Alternatively, one could reduce the extent of the flux tubes in the z−direction, thereby separating the flux response from the surface states.

VI. ISOLATED FIRST-AND HIGHER-ORDER FLUX MAJORANA MODES
All NH flux responses discussed so far have led to a SIBC spectrum with a completely filled point gap in presence of π-flux defects. Surprisingly, it is also possible to obtain isolated flux modes at specific energies within the point gap. These are similar to Hermitian bound states, but associated with NH point gap topology.
To understand this phenomenon, we consider NH systems in class D, which are classified by the Z-valued 3D winding number W 3D (E) [16] (see App. E 2 d for details). We derive the flux response for a given energy E 0 inside the nontrivial point gap, in the simplest case with W 3D (E 0 ) = 1, from the corresponding EHH (see Sec. II B). The EHH for E 0 = 0 of NH class D retains the symmetries of Hermitian class DIII (see App. B 5 e and Tab. III therein). Hermitian class DIII is also Z-classified in 3D, such that the EHH associated with a nontrivial NH phase in class D must itself realize a nontrivial topologi-cal insulator phase in Hermitian class DIII [59]. For this phase, it was shown in Ref. 10 that flux tubes bind two degenerate zero-energy states (see Fig. 6a), corresponding to a single flux-localized state at E 0 = 0 in the NH SIBC spectrum (see Fig. 6b). Away from E 0 = 0, the EHH retains only a chiral symmetry, resulting in Hermitian class AIII (see App. B 5 e). The flux response of the EHH in Hermitian class AIII is trivial, such that the point gap remains empty for E 0 ̸ = 0 (see Fig. 6c).
We diagnose the nature of the unpaired NH flux state at E 0 = 0 by investigating (in a gedankenexperiment) its stability under coupling with a 1D line-gapped phase spanned between the two flux tubes. 1D line gap topology in NH class D is Z 2 -classified (see Tab. I). The nontrivial phase is nothing but the NH generalization of the 1D Kitaev chain [39,66], which hosts a single Majorana zero-mode at each end. Two zero-modes -the end state of the 1D line gap phase combined with the flux localized mode -are not protected. Consequently, the flux state must be a Majorana mode as well, and flux defects in NH class D can only probe W 3D (E) mod 2. An unpaired flux Majorana mode remains pinned at zero complex energy. It is localized along the entire length of the flux tube (see Fig. 6d), highlighting another NH peculiarity: in Hermitian systems, spectrally isolated bound states are pointlike.
We next study NH class D S+ , which is also classified by a Z invariant [16], the 3D winding number W 3D (E) restricted to even values (see App. E 2 e for details). Choosing an energy E 0 inside the nontrivial point gap with W 3D (E 0 ) = 2 allows to construct an EHH: for E 0 = 0 the EHH contains two interdependent unitary subspaces in Hermitian class AIII (see App. B 5 f and Tab. III therein). Due to the Z-classification of Hermitian class AIII in 3D, the EHH associated with a nontrivial phase in NH class D S+ must itself realize a nontrivial topological insulator phase per unitary subspace in Hermitian class AIII. However, without additional crystalline sym- The energy spectrum of the EHH for NH class D realizes a helical metal for E0 = 0. b This corresponds to a single Majorana mode at E0 = 0 in the NH point gap. c Away from E0 = 0, the EHH is gapped. This implies that the NH point gap does not fill up with flux-localized states, but only pins the Majorana mode. d The Majorana mode is localized along the entire flux defect (k ∥ = 0), a scenario that is not possible for a single mode in a Hermitian system. e The energy spectrum of the EHH for NH class D S + with pseudo-inversion symmetry is gapped in presence of two PBC-preserving π-fluxes. f This results in an empty point gap in the NH spectrum. g When terminating the system along a surface perpendicular to the flux tubes, the EHH features four exact zero modes at E0 = 0, which are absent when the flux tubes are removed. (Note that in this limit, approximate zero-modes remain due to the gapless surface states in this symmetry class, highlighted in the inset for ϕ = 0.) h The two resulting NH higher-order Majorana modes are then localized at the ends of the flux defects. Panels a-d are generated for a model of size 20 × 20 (×20 in OBC) unit cells in NH class D (see Sec. E 2 d), panels e-h for a model of size 20 × 20 (×20 in OBC) unit cells in NH class D S + (see Sec. E 2 e). metries, the zero energy modes arising in the flux Dirac theory of systems in Hermitian class AIII are not protected [10] (see Fig. 6e). Adding an inversion symmetry (corresponding to pseudo-inversion in the NH Hamiltonian, App. D 2) quantizes the flux and allows for a stable flux response [14]: under OBC, each flux tube now hosts one zero energy mode per unitary subspace of the EHH. These flux end states then correspond to flux-localized NH SIBC spectral in-gap states at E 0 = 0 (see Fig. 6f and g). Away from E 0 = 0, the EHH is situated in Hermitian class CI: due to the absence of Kramers theorem, pairs of zero-energy states are no longer protected (see App. B 5 f). Consequently, the flux modes in each unitary subspace can be gapped, so that the higher order flux response persists only at E 0 = 0. The resulting NH higher-order Majorana mode localizes at the ends of the flux tubes (see Fig. 6h). Detecting its presence in a real system is, however, intricate, as it will in general be obscured by the NH surface state that is already present for ϕ = 0.
We can understand the existence of (higher-order) flux Majorana modes by noting that, among all NH symmetry classes that do not exhibit a flux spectral jump or higher-order flux skin effect in 3D, NH class D and D S+ are the only ones with a nontrivial line gap classification (see Tab. I) in 1D. As line-gapped phases are adiabatically deformable to Hermitian systems, flux tubes in these classes can only cause the presence of isolated modes (similar to the end states of 1D Hermitian topological insulators) instead of NH skin effects.

VII. EXPERIMENTAL REALIZATION OF NH FLUX RESPONSE
In a solid-state setting, we usually assume isolated systems that are governed by a Hermitian Hamiltonian. Most meta-material platforms are however either accidentally or tunably lossy, such that the description with an effective Hamiltonian involves NH terms. The same holds for interacting electronic quantum systems in which quasiparticles acquire a finite lifetime. Such a scenario equips the single-electron Green's function with a complex self-energy, leading to an effective NH Hamilto-nian [67][68][69][70][71][72]. Consequently, NH systems appear quite naturally in experimental settings.
We now outline the NH platforms in which flux defects can be studied experimentally. Superconductors are natural starting points for studying flux defects, as magnetic fields penetrate samples as flux vortices. Consequently, superconducting NH systems in 2D or those in direct proximity to a superconductor may experience a flux skin effect in presence of vortices. The size of vortices should be matched with the lattice scale by the use of moire substrates/potentials, allowing for precisely localized defects and quantized responses. However, in general we expect flux localized states to not disappear immediately with increasing vortex size, but rather to spread out over the vortex region. In 3D, the prototypical NH point-gapped phase is the exceptional topological insulator (ETI) [21]. The ETI emerges naturally from a Hermitian 3D topological insulator or Weyl semimetal, if quasi-particles acquire a finite lifetime. This could for instance be caused by electron or electron-phonon interactions [21]. As the ETI does not require any symmetry to be stabilized, it is naturally situated in NH class A. With additional pseudo-inversion symmetry, it can thus give rise to a higher-order flux skin effect upon flux insertion. Additionally, the ETI phase was shown to arise in NH class AII † [21], thereby allowing for the flux realization of the NH flux spectral jump. Consequently, the ETI serves as the ideal platform to investigate 3D flux defects. Due to their versatility, meta-materials provide convenient classical analogs to quantum mechanical topological states. A flux defect can for instance be realized in electrical meta-materials [35,73] by introducing operational amplifiers. These equip arbitrary hoppings t with a tuneable phase t → te iϕ , thereby allowing to implement ϕ flux tubes [74]. Alternatively, one might realize effective magnetic fluxes in photonic [75], mechanical [76], and ultra-cold atom systems [77,78]. Additionally, one can envision a realization in NH superconductors [79][80][81]. Such systems, for instance in NH class C S+ and DIII S+− , can give rise to both higher-order flux skin effect and flux spectral jump upon flux insertion, respectively. Similarly, the NH flux Majorana mode appears in a NH superconductor in class D.

VIII. DISCUSSION
We have derived the flux response of all NH symmetry classes with intrinsic point gap topology in 2D and 3D. In 2D, we found the NH flux skin effect, in which a macroscopic number of states localizes at the flux core. We identified the necessary conditions for its emergence and predict its occurrence for four NH symmetry classes. In 3D, we discovered the flux spectral jump, where a π-flux tube causes a NH skin effect only at a single momentum along the flux tube direction. This response derives from a corresponding novel 2D phase exhibiting an infernal point, where an extensive number of states collapses to the boundary at a single momentum. Developing a physical understanding for this phenomenon, including a field-theoretical description, is an important future endeavor. Additionally, we found NH symmetry classes giving rise to a higher-order flux skin effect, which occurs at the ends of flux tubes. Being a higher-order response, it relies on the presence of a crystalline symmetry. We here considered the case of pseudo-inversion symmetry, but expect further exciting responses for other crystalline symmetries. Finally, we identified the presence of flux Majorana modes, forming spectrally isolated mid-gap states in the NH point gap. Realizing such modes in open quantum systems is an interesting question for future theoretical and experimental research. This section defines symmetries in NH systems and their relation to the corresponding symmetries of the EHH. Hermitian quantities are denoted with an overline.

NH symmetries
The Hermitian time-reversal symmetry given as is generalized to a NH time-reversal symmetry TRS as well as a pseudo time-reversal symmetry TRS † The Hermitian particle-hole symmetry given as is generalized to a NH particle-hole symmetry PHS as well as a pseudo particle-hole symmetry PHS † The Hermitian chiral symmetry given as is generalized to a NH chiral symmetry CS Additonally, one can have sublattice symmetry SLS, defined by as well as pseudo-Hermiticity Finally, we investigate the presence of inversion and pseudo-inversion symmetry

Appearance of NH symmetries in the EHH
The presence of symmetries for the NH Hamiltonian H(k) imposes constraints on the EHHH(k), Specifically, for E 0 = 0, it holds: for TRS andŪ for TRS † of the NH Hamiltonian. withŪ for PHS andŪ for PHS † of the NH Hamiltonian. Additionally we havē By construction,H(k) enjoys an additional chiral (sublattice) symmetry for arbitrary E 0 : Besides the chiral symmetry introduced in (A23), sublattice symmetry S and PHS U P allow a TRS and a PHSŪ For daggered NH symmetry classes, TRS † can be combined with (A23) to a PHS withŪ Appendix B: Classification of flux response in all NH point-gapped symmetry classes We classify the flux response of all NH systems with nontrivial intrinsic point gap topology. Since topological zero energy eigenvalues of the EHHH(k) correspond to states at E 0 within the NH point gap, we first investigate the flux response of Hermitian systems. This procedure relies on a Dirac treatment of the flux defect, introduced in App. B 1. We employ this procedure for all relevant Hermitian symmetry classes in 2D (App. B 2) and 3D (App. B 3). The corresponding NH flux response follows from protected EHH zero energy modes, derived in App. B 4 for 2D and App. B 5 for 3D.

Dirac theory of Hermitian flux response
For a Hermitian topological insulator in one of the 10 Altland-Zirnbauer (AZ) symmetry classes [82,83], we can derive the response for a ϕ = π flux using the Dirac theory of its edge states [84].
To derive the flux response of 2D (3D) topological insulators, we first cut the 2D (3D) bulk in half (see Fig. 7a) to create two sets of edge (surface) states, one for the top layer and one for the bottom layer. In general, these will be described by a 1D (2D) Dirac Hamiltonian where k is the momentum along the edge and k the momenta on the surface, with τ z a Pauli matrix acting on layer space. In the presence of N edge (surface) states, h(k) is a N ×N matrix capturing the edge (surface) state dispersion, and M is a 2N × 2N mass matrix that couples the edge (surface) states to yield a gapped bulk. The insertion of a flux core (flux tube) ϕ can be viewed as a  [14]. c In 3D, π-fluxes similarly bind 1D states.
modification of all hoppings crossing a line (plane) emanating from the flux core (flux tube), multiplying all such hopping amplitudes with a Peierls phase t → te iϕ . For a flux of ϕ = π, this corresponds to a sign flip [14]. The mass term M coupling the two layers therefore changes sign at the flux tube, forming a domain wall binding point-like (line-like) states (see Fig. 7b,c). In the following, we therefore consider a flux-Dirac theory for all relevant Altland-Zirnbauer classes.

2D systems, symmetry classes ofH(k)
In 2D, the EHHH(k) for nontrivial intrinsic pointgapped NH models resides in Hermitian classes A, AII, D, DIII or C (see Tab. II). In the following, we investigate whether a flux induced mass domain wall leads to zero energy modes ofH(k).

a. Class A
The prototypical model in Hermitian class A in 2D is a Chern insulator, whose edge Hamiltonian is described where k is the momentum along the edge. The 1D Dirac Hamiltonian can then be written as where mx implements the flux core induced domain wall in the mass term at x = 0, τ µ are Pauli matrices (µ = 0, x, y, z) and δ multiplies a symmetry allowed perturbation. In order to derive the presence of topological zero energy states, we consider OBC along the edge (x−direction) and solve for boundary-localized zeroenergy states determined by the Hamiltonian For δ = 0 we can find one normalizable zero energy solution |ψ⟩,H x |ψ⟩ = 0, given by To study the effect of the perturbation δ, we rely on firstorder perturbation theory Consequently, the EHH for flux cores in Hermitian class A (d = 2) has no zero energy mode, as any finite δ is able to remove the in gap state. This means the flux response of 2D NH systems that lead to Hermitian class A forH(k) is trivial. Note that additional crystalline symmetries do not change this statement, contrary to the case of 3D.

b. Class AII
Hermitian class AII has a TRS withŪ TŪ * T = −1, describing 2D topological insulators. As such, the surface Hamiltonian of the nontrivial phase is gapless, where k is the momentum along the edge and mass terms are forbidden by the presence of TRSŪ T = σ y . The 1D Dirac Hamiltonian can then be written as where mx implements the flux core induced domain wall in the mass term at x = 0, τ µ , σ µ are Pauli matrices (µ = 0, x, y, z) and δ multiplies a collection of symmetry allowed perturbation terms O. In order to derive the presence of topological zero energy states, we consider OBC along the edge (x−direction) and solve for boundary-localized zero-energy states determined by the HamiltonianH For δ = 0 we can find two normalizable zero energy solutions |ψ⟩,H x |ψ⟩ = 0, given by In order to reveal that symmetry allowed perturbations O are able to gap out these modes, we rely on degenerate first-order perturbation theory in δ. Fixing O = τ 0 σ 0 as one symmetry preserving choice under TRSŪ T = τ 0 σ y , we obtain Consequently, the EHH for flux cores in Hermitian class AII (d = 2) is gapped and hosts no zero energy modes. This means the flux response of 2D NH systems that lead to Hermitian class AII forH(k) is trivial.
c. Class D Hermitian class D describes the thermal quantum Hall effect, with a PHS withŪ PŪ * P = +1. The surface Hamiltonian for the nontrivial phase is gapless, where k is the momentum along the edge. The 1D Dirac Hamiltonian follows as where mx implements the flux core induced domain wall in the mass term at x = 0, τ µ are Pauli matrices (µ = 0, x, y, z) and further symmetry allowed perturbation terms O are not allowed by the presence ofŪ P = τ z . In order to derive the presence of topological zero energy states, we consider OBC along the edge (x−direction) and solve for boundary-localized zero-energy states determined by the Hamiltonian We can find one normalizable zero energy solution |ψ⟩, H x |ψ⟩ = 0, given by Consequently, the EHH for flux cores in Hermitian class D (d = 2) is gapless, and hosts one topologically protected zero mode localized at the flux core.

d. Class DIII
In Hermitian symmetry class DIII, the surface Hamiltonian along the edgeh(k) is gapless in the nontrivial phase and assumes the form where k is the momentum along the edge and mass terms are forbidden by the presence of TRS, PHS and chiral symmetry. Choosing chiral symmetry asŪ C = τ z σ z , the 1D Dirac Hamiltonian can be written as where mx implements the flux core induced domain wall in the mass term at x = 0, τ µ , σ µ are Pauli matrices (µ = 0, x, y, z) and δ multiplies a collection of symmetry allowed perturbation terms O. In order to derive the presence of topological zero energy states, we consider OBC along the edge (x−direction) and solve for boundary-localized zero-energy states determined by the HamiltonianH For δ = 0 we can find two normalizable zero energy solutions |ψ⟩,H x |ψ⟩ = 0, given by In order to reveal that symmetry allowed perturbations O are not able to gap out these modes, we rely on degenerate first-order perturbation theory in δ. Fixing O = τ y σ z as the only symmetry preserving choice under TRSŪ T = τ 0 σ y and PHSŪ P = τ z σ x , we obtain H flux mn = ⟨ψ m |H x |ψ n ⟩ = 0. (B21) Consequently, the EHH for flux cores in Hermitian class DIII (d = 2) is gapless, and hosts two topologically protected zero modes localized at the flux cores.
e. Class C In Hermitian symmetry class C withŪ PŪ * P = −1, the surface Hamiltonian along the edgeh(k) is gapless in the nontrivial phase and assumes the form where k is the momentum along the edge and the terms multiplied bym 1,2,3 only shift the zero energy crossing. UsingŪ P = τ z σ y , the 1D Dirac Hamiltonian can then be written as where mx implements the flux core induced domain wall in the mass term at x = 0, τ µ , σ µ are Pauli matrices (µ = 0, x, y, z) and δ multiplies a collection of symmetry allowed perturbation terms O. In order to derive the presence of topological zero energy states, we consider OBC along the edge (x−direction) and solve for boundary-localized zero-energy states determined by the HamiltonianH For δ = 0 we can find two normalizable zero energy solutions |ψ⟩,H x |ψ⟩ = 0, given by To study the influence of symmetry allowed perturbations O, we rely on degenerate first-order perturbation theory in δ. Fixing O = τ 0 σ z as one symmetry preserving choice under PHSŪ P = τ z σ y , we obtain Consequently, the EHH for flux cores in Hermitian class C (d = 2) is gapped, and hosts no zero modes localized at the flux cores. This is a consequence of PHS being able to protect only a single state, not two.

3D systems, symmetry classes ofH(k)
In 3D, the EHHH(k) for nontrivial point-gapped NH models resides in Hermitian classes AII, AIII, DIII, CI or CII (see Tab. III). In the following, we investigate whether a flux induced mass domain wall leads to zero energy modes ofH(k).

a. Class AII
In App. B 2 b we derived the absence of zero energy modes at flux defects in 2D Hermitian class AII. In 3D, the edge Hamiltonian in the nontrivial phase is given as where k 1 is the momentum along the edge, k 2 the momentum along the flux tube, we obtain the 1D Dirac Hamiltonian as with the Pauli matrices τ µ , σ µ (µ = 0, x, y, z). For k 2 = δ = 0 we still find the two normalizable zero energy solutions given by (B11). In order to derive the effective Hamiltonian along the flux tube, we rely on degenerate first-order perturbation theory in k 2 and δ. Using δO = i=x,y,z δ i τ y σ i + δ 4 τ 0 σ 0 + δ 5 τ z σ 0 as the symmetry preserving choices under TRSŪ T = τ 0 σ y , we obtain Since the identity matrix constitutes just a trivial spectral shift, the Hamiltonian along the flux tube is still gapless. Consequently, flux tubes in Hermitian class AII (d=3) form a helical metal, with zero modes ofH(k) appearing at k 2 = 0.

b. Class AIII
Hermitian class AIII contains a chiral (sublattice) symmetry, such that the edge Hamiltonian for the nontrivial phase appears in (B2) as where k 1 is the momentum along the edge, k 2 the momentum along the flux tube, M the mass matrix that couples the edge states to yield a gapped bulk, τ µ , σ µ Pauli matrices (µ = 0, x, y, z) and δ multiplies a collection of symmetry-allowed perturbation terms O. The insertion of a π−flux tube forms a mass domain wall, causing a sign change in M . Note however that Hermitian class AIII does not contain a symmetry fixing the flux to zero or π. In order to derive the presence of topological zero energy states, we consider OBC along the edge (x−direction). Without further symmetries, there exists an ambiguity in the suitable mass domain wall, which can equivalently be realized as τ y σ 0 or τ x σ 0 . Therefore there exists no bound state without additional (crystalline) symmetries. Adding for instance inversion symmetryĪ, realized here asĪ = τ x σ 0 , removes the ambiguity by rendering τ y σ 0 an incompatible choice. However, our current single mass domain wall inherently breaks inversion symmetry. We therefore have to introduce a second flux tube, representing the opposite mass domain wall. A corresponding 2D-Dirac theory is discussed in Ref. 14, showing that each flux tube localizes one zero energy mode under OBC.

c. Class DIII
Since the DIII 3D topological insulator can be understood as a pumping cycle of the corresponding 2D system along a third momentum direction, the flux response follows directly from the respective case in 2D. In App. B 2 d we derived the presence of two zero modes under π-flux insertion. Accounting for the additional dimension by usingh where k 1 is the momentum along the edge, k 2 the momentum along the flux tube, we obtain the 2D Dirac Hamiltonian as with the Pauli matrices τ µ , σ µ (µ = 0, x, y, z). For k 2 = δ = 0 we still find the two normalizable zero energy solutions given by (B20). In order to derive the effective Hamiltonian along the flux tube, we rely on degenerate first-order perturbation theory in k 2 and δ. Using again O = τ y σ z as the only symmetry preserving choice under TRSŪ T = τ 0 σ y and PHSŪ P = τ z σ x , we obtain Consequently, flux tubes in Hermitian class DIII (d=3) form a gapless helical metal, with zero modes ofH(k) appearing at k 2 = 0.

d. Class CI
Hamiltonians in Hermitian symmetry class CI possess TRS, chiral and PHS, where only the latter squares to minus one. The corresponding surface Hamiltonian along the edgeh(k 1 , k 2 ) is gapless in the nontrivial phase and assumes the form where k 1 is the momentum along the edge, k 2 the momentum along the flux tube, ρ µ , τ µ , σ µ Pauli matrices (µ = 0, x, y, z) and allowed mass terms only constitute trivial shifts of the location of the zero energy mode in momentum space. UsingŪ T = τ 0 ρ 0 σ 0 ,Ū P = τ z ρ y σ z , U C = τ z ρ y σ z the 1D Dirac Hamiltonian can then be written asH where mx implements the flux tube induced domain wall in the mass term at x = 0 and δ multiplies a collection of symmetry allowed perturbation terms O. In order to derive the presence of topological zero energy states, we consider OBC along the edge (x−direction) and solve for boundary-localized zero-energy states determined by the Hamiltonian For k 2 = δ = 0 we can find four normalizable zero energy solutions |ψ⟩,H x |ψ⟩ = 0, given by (B37) In order to derive the effective Hamiltonian along the flux tube, we rely on degenerate first-order perturbation theory in k 2 and δ. Using O = τ y ρ y σ 0 as one symmetry preserving choice, we obtain Consequently, flux tubes in Hermitian class CI (d=3) are gapped, and also host no zero energy end states under OBC.

e. Class CII
In Hermitian symmetry class CII, the surface Hamiltonian along the edgeh(k 1 , k 2 ) is gapless for the nontrivial phase and assumes the form where k 1 is the momentum along the edge, k 2 the momentum along the flux tube, ρ µ , σ µ are Pauli matrices (µ = 0, x, y, z) and mass terms are forbidden by the presence of TRS, PHS and chiral symmetry. Choosing chiral symmetry asŪ C = τ z ρ z σ z , the 1D Dirac Hamiltonian can be written as where M implements the π−flux tube induced domain wall in the mass term at x = 0, τ µ are Pauli matrices (µ = 0, x, y, z) and δ multiplies a collection of symmetry allowed perturbation terms O. Note however that Hermitian class AIII does not contain a symmetry fixing the flux to zero or π. Requiring TRS asŪ T = iτ 0 ρ x σ y and PHS asŪ P = iτ z ρ y σ x , one can find two possible mass terms M , τ x ρ 0 σ 0 or τ y ρ z σ 0 . Due to this ambiguity there exists no bound state without additional (crystalline) symmetries. Adding for instance inversion symmetryĪ, realized here asĪ = τ x ρ 0 σ 0 , removes the ambiguity by rendering τ y ρ z σ 0 an incompatible choice. We therefore choose the mass domain wall to be realized as M = mxτ x ρ 0 σ 0 . However, the single mass domain wall again breaks inversion symmetry. We therefore have to introduce a second flux tube, representing the opposite mass domain wall. Using again the results of the 2D-Dirac theory discussed in Ref. 14, one obtains two zero energy modes per flux tube, localized at the flux tube ends.

NH flux response in 2D
Zero modes of the EHH and NH point gap states are in one-to-one correspondence. Using the Hermitian flux Dirac theory, we here derive the flux response and the corresponding topological properties of the SIBC spectrum of all intrinsically nontrivial 2D point-gapped NH symmetry classes, summarized in Tab. II. a. Class AII † NH class AII † has a TRS † squaring to minus one, which quantizes the flux ϕ = 0, π. As outlined in A 2, the corresponding EHH introduces a chiral symmetrȳ Σ C , which combines with TRS to form a PHS with U PŪ * P = −Ū TŪ * T = +1. Consequently, the EHH is in Hermitian class DIII, which is Z 2 classified for both 1D and 2D [83].
The flux response of the EHH in Hermitian class DIII shows two zero energy modes for each flux core (see App. B 2 d). As the construction of the EHH can be repeated for every complex energy inside the point gap, the point gap of a corresponding NH system fills with flux core localized modes. Consequently, models in NH class AII † show a flux skin effect for ϕ = π. A periodic system contains two flux cores ϕ = ±π, each localizing an extensive number of modes. This number scales with the extension of the system along the direction containing the two defects, denoted by L ⊥ (see App. C 1).

b. Class DIII †
The NH class DIII † possesses TRS, PHS and chiral symmetry with (U T U * T , U P U * P ) = (−1, 1), where U C = U T U P . As outlined in App. A 2, the corresponding EHH at E 0 = 0 obtains TRSŪ T , PHSŪ P and two chiral symmetries {Ū C ,Σ C } = 0. Both chiral symmetries can be combined to a unitary symmetryŪ =Ū CΣC with U 2 = −1.
SinceŪ has an imaginary spectrum and [Ū T ,Ū ] = 0, the eigenspaces ofŪ are exchanged by anti-unitary TRS. Since {Ū P ,Ū } = 0, we retain anti-unitary PHS in each eigenspace, corresponding to Hermitian class D which has a Z classification in 2D [83]. After modding out line-gap phases, this is reduced to a Z 2 classification [39]. Since 1D Hermitian class D is also Z 2 classified, we expect that a π-flux induces a single bound state perŪ eigensector (see App. B 2 c). Hence, the full Hermitian spectrum hosts a zero-energy Kramers pair per flux.
Away from E 0 = 0, we retain onlyŪ TŪ * T = −1 wherē U T andΣ C anti-commute, resulting in Hermitian class DIII. Due to TRS and chiral symmetry, the flux-bound Kramers pair cannot move away from zero energy. Hence the entire point gap of a corresponding NH system fills with flux core localized modes. Consequently, models in NH class DIII † quantize the flux to ϕ = 0, π and show a flux skin effect for ϕ = π.
In App. B 2 a, we derived the flux response of EHHs in Hermitian class A, which did not yield protected zero energy modes pinned to the flux defects. Consequently, NH systems in class AIII S− do not show flux induced states within the NH point gap. Therefore, the flux response for this NH symmetry class is trivial.
Away from E 0 = 0, we can only form the TRS . This results in Hermitian class DIII, such that due to TRS and chiral symmetry, the flux-bound Kramers pair cannot move away from zero energy. Hence the entire point gap of a corresponding NH system fills with flux core localized modes. Consequently, models in NH class BDI S+− quantize the flux to ϕ = 0, π and show a flux skin effect for ϕ = π. e. Class D S − Models in NH class D S− possess PHS U P (U P U * P = +1) and sublattice symmetry S, with {U P , S} = 0. As outlined in App. A 2, the corresponding EHH introduces a chiral symmetryΣ C , which combines with sublattice symmetry to form a unitary symmetryŪ =SΣ C with U 2 = +1. The unitaryŪ satisifies [Ū P ,Ū ] = [S,Ū ] = [Σ C ,Ū ] = 0. We may furthermore define a TRSŪ T = U PΣC which satisfiesŪ TŪ * T = −1. Hence we obtain Hermitian class DIII in eachŪ subspace, giving a Z 2 ⊕ Z 2 classification in 2D [83]. By modding out line-gap phases, this is reduced to a Z 2 classification where the nontrivial element corresponds to having only a singleŪ subspace being nontrivial [39]. Hence the full Hermitian spectrum hosts a single zero-energy Kramers pair under π-flux insertion (see App. B 2 d).
Away from E 0 = 0, we can only form the TRS . This results in Hermitian class DIII, such that the flux-bound Kramers pair cannot move away from zero energy. Hence the entire point gap of a corresponding NH system fills with flux core localized modes. Consequently, models in NH class D S− quantize the flux to ϕ = 0, π and show a flux skin effect for ϕ = π.

f. Class DIII S +−
The NH class DIII S+− possesses TRS, PHS and chiral symmetry with (U T U * T , U P U * P ) = (−1, 1), where U C = U T U P , as well as sublattice symmetry S. As outlined in App. A 2, the corresponding EHH at E 0 = 0 obtains TRS U T , PHSŪ P and three chiral symmetriesŪ C ,Σ C ,S. The chiral symmetries can be combined to two commuting unitary symmetriesŪ 1 =Ū CΣC withŪ 2 1 = −1 andŪ 2 = SΣ C withŪ 2 2 = +1. SinceŪ 1 has an imaginary spectrum and the eigenspaces ofŪ 1 are not independent and individually enjoyŪ T andŪ 2 symmetry. Moreover, we have so that theŪ 2 eigenspaces are independent and individually preserveŪ T symmetry. They therefore lie in Hermitian class AII and yield a Z 2 ⊕ Z 2 classification in 2D [83] that is reduced to Z 2 by line-gap phases [39]. The nontrivial element corresponds to having only a singleŪ 2 subspace nontrivial. In App. B 2 b, we derived the flux response of EHHs in Hermitian class AII, which did not yield protected zero energy modes pinned to the flux defects. Consequently, systems in NH class DIII S+− do not show flux induced states within the NH point gap. Therefore, the flux response for this NH symmetry class is trivial.

g. Class CII S −+
The NH class CII S−+ possesses TRS, PHS and chiral symmetry with (U T U * T , U P U * P ) = (−1, −1), where U C = U T U P , as well as sublattice symmetry S. As outlined in App. A 2, the corresponding EHH at E 0 = 0 obtains TRS U T , PHSŪ P and three chiral symmetriesŪ C ,Σ C ,S. The chiral symmetries can be combined to two commuting unitary symmetriesŪ 1 =Ū CΣC withŪ 2 1 = −1 andŪ 2 = SΣ C withŪ 2 2 = +1. SinceŪ 1 has an imaginary spectrum and the eigenspaces ofŪ 1 are not independent and individually enjoyŪ P andŪ 2 symmetry. Moreover, we have {Ū T ,Ū 2 } = {Ū P ,Ū 2 } = [S,Ū 2 ] = [Ū C ,Ū 2 ] = [Σ C ,Ū 2 ] = 0, so that theŪ 2 eigenspaces are exchanged byŪ P symmetry. This leaves us with Hermitian class A, which has a Z classification in 2D [83] that is reduced to Z 2 by line-gap phases [39]. In App. B 2 a, we derived the flux response of EHHs in Hermitian class A, which did not yield protected zero energy modes pinned to the flux defects. Consequently, systems in NH class CII S−+ do not show flux induced states within the NH point gap. Therefore, the flux response for this NH symmetry class is trivial.

h. Class CII S +−
The NH class CII S+− possesses TRS, PHS and chiral symmetry with (U T U * T , U P U * P ) = (−1, −1), where U C = U T U P , as well as sublattice symmetry S. As outlined in App. A 2, the corresponding EHH at E 0 = 0 obtains TRS U T , PHSŪ P and three chiral symmetriesŪ C ,Σ C ,S. The chiral symmetries can be combined to two commuting unitary symmetriesŪ 1 =Ū CΣC withŪ 2 1 = −1 andŪ 2 = SΣ C withŪ 2 2 = +1. SinceŪ 1 has an imaginary spectrum and the eigenspaces ofŪ 1 are not independent and individually enjoyŪ P andŪ 2 symmetry. Moreover, we have [Ū T ,Ū 2 ] = [Ū P ,Ū 2 ] = [S,Ū 2 ] = [Ū C ,Ū 2 ] = [Σ C ,Ū 2 ] = 0, so that theŪ 2 eigenspaces are independent and individually preserveŪ P symmetry. We therefore obtain Hermitian class C, giving a Z ⊕ Z classification in 2D [83]. This is reduced to Z 2 under the addition of line-gap phases, and the nontrivial element corresponds to only having oneŪ 2 eigenspace nontrivial [39]. Since the 1D classification of Hermitian class C is trivial and App. B 2 e showed no protected zero energy modes, we do not expect a stable flux response. Therefore, the flux response for systems in NH class CII S+− is trivial.

i. Class CI S −+
The NH class CI S−+ possesses TRS, PHS and chiral symmetry with (U T U * T , U P U * P ) = (1, −1), where U C = U T U P , as well as sublattice symmetry S. As outlined in App. A 2, the corresponding EHH at E 0 = 0 obtains TRS U T , PHSŪ P and three chiral symmetriesŪ C ,Σ C ,S. The chiral symmetries can be combined to two commuting unitary symmetriesŪ 1 =Ū CΣC withŪ 2 1 = −1 andŪ 2 = SΣ C withŪ 2 2 = +1. SinceŪ 1 has an imaginary spectrum and the eigenspaces ofŪ 1 are not independent and individually enjoyŪ T andŪ 2 symmetry. Moreover, we have {Ū T ,Ū 2 } = {Ū P ,Ū 2 } = [S,Ū 2 ] = [Ū C ,Ū 2 ] = [Σ C ,Ū 2 ] = 0, so that theŪ 2 eigenspaces are exchanged byŪ T symmetry. This leaves us with Hermitian class A, which has a Z classification in 2D [83] that is reduced to Z 2 by line-gap phases [39]. In App. B 2 a, we derived the flux response of EHHs in Hermitian class A, which did not yield protected zero energy modes pinned to the flux defects. Consequently, systems in NH class CI S−+ do not show flux induced states within the NH point gap. Therefore, the flux response for this NH symmetry class is trivial.

NH flux response in 3D
Similar to the case of 2D, we rely on the Hermitian flux Dirac theory to derive the flux response of all nontrivial 3D point-gapped NH symmetry classes, summarized in Tab. III.
a. Class AI † NH class AI † has a TRS † squaring to plus one. As outlined in App. A 2, the corresponding EHH introduces a chiral symmetryΣ C , which combines with TRS to form a PHS withŪ PŪ EHH for all energies inside the point gap E 0 is in Hermitian class CI.
The flux response of the EHH in Hermitian class CI shows no protected zero energy modes (see App. B 3 d). Consequently, systems in NH class AI † do not show flux induced states within the NH point gap. Therefore, the flux response for this NH symmetry class is trivial. b. Class AII † As for the case of 2D, NH class AII † has a TRS † squaring to minus one, which quantizes the flux ϕ = 0, π. As outlined in App. A 2, the corresponding EHH introduces a chiral symmetryΣ C , which combines with TRS to form a PHS withŪ PŪ * P = −Ū TŪ * T = +1. Consequently, the EHH is in Hermitian class DIII, which is Z classified for 3D and hosts a Z 2 index for 2D [83].
The flux response of the EHH in Hermitian class DIII forms a gapless helical metal along the flux tube (see App. B 3 c). As the construction of the EHH can be repeated for every complex energy inside the point gap, the point gap of a corresponding NH system fills with an extensive number of states at a single momentum k ∥ . This is the defining signature of the NH flux spectral jump, discussed in Sec. IV. Therefore, systems in NH class AII † show a NH flux spectral jump when introducing a flux defect. Note that due to the Z 2 -classification in 2D, flux defects only probe the Z 2 nature of the 3D Z invariant in NH class AII † . c. Class A+I † As NH class A does not contain any symmetries, the corresponding EHH only introduces a chiral symmetrȳ Σ C (see App. A 2). Consequently, the EHH is in Hermitian class AIII, which is Z classified for 3D and 1D while being trivial in 2D [83].
Without additional crystalline symmetries, the zero energy modes arising in the flux Dirac theory of systems in Hermitian class AIII are not protected (see App. B 3 b). This indicates the absence of a NH flux spectral jump, but allows for a different, higher order phase. Adding a crystalline symmetry I † serves two purposes: first, it should fix the flux to be either zero or π, to allow for a topological flux response. Second, it should protect zero energy states from being gapped by perturbations. A crystalline symmetry fixing the flux requires the presence of at least two defects, resulting in a 2D Hermitian flux Dirac theory. The results derived in App. B 3 b show the presence of zero energy states under OBC along the flux tubes, which can be distinguished from the surface signature for ϕ = 0. Consequently, these zero modes appear extensively in the point gap of the NH model. This constitutes a higher order skin effect, with skin modes appearing at the ends of the flux tubes in a finite geometry. The precise localization of skin modes is then determined by the present pseudo inversion symmetry I † . This is the defining signature of the NH higher-order flux skin effect, discussed in Sec. V. d. Class A S + I †,S − NH class A S possesses a sublattice symmetry S, which combines with the chiral symmetryΣ C in the corresponding EHH (see App. A 2) to form a unitary symmetrȳ U =SΣ C withŪ 2 = +1. The unitaryŪ satisifies [S,Ū ] = [Σ C ,Ū ] = 0. Hence we obtain Hermitian class AIII in eachŪ subspace, giving a Z ⊕ Z classification in 3D and 1D, while 2D is trivial [83]. By modding out linegap phases, this is reduced to a Z classification where the nontrivial element corresponds to having only a singleŪ subspace being nontrivial [39].
Away from E 0 = 0, we are only left with the chiral symmetryΣ C , resulting in Hermitian class AIII. Without additional crystalline symmetries, the zero energy modes arising in the flux Dirac theory of systems in Hermitian class AIII are not protected (see App. B 3 b). Adding a crystalline symmetry I †,S− , however, quantizes the flux and allows for a stable flux response. The results derived in App. B 3 b show the presence of zero energy states under OBC along the flux tubes, which can be distinguished from the surface signature for ϕ = 0. Consequently, these zero modes appear extensively in the point gap of the NH model, constituting a higherorder flux skin effect (see Sec. V) for NH class A S +I †,S− . e. Class D Models in NH class D possess a PHS U P with U P U * P = +1, yielding a Z classification in 3D, while lower dimensions are trivial. As outlined in App. A 2, the corresponding EHH introduces a chiral symmetryΣ C , which combines with PHS to a TRSŪ T =Ū PΣC which satis-fiesŪ TŪ * T = −1. Hence we obtain Hermitian class DIII, giving a Z classification in 3D [83]. The flux response of the EHH in Hermitian class DIII forms a gapless helical metal along the flux tube (see App. B 3 c).
Away from E 0 = 0, we are only left with the chiral symmetryΣ C , resulting in Hermitian class AIII. Hermitian class AIII allows to gap the helical metal by introducing suitable mass terms. Consequently, only at E 0 = 0 we expect flux localized modes to appear in the NH point gap. The number of modes does not scale extensively with the system size, the NH system only pins the Hermitian surface state of the extended model. This surface state is a Majorana, models in NH class D therefore show a flux Majorana mode, as introduced in Sec. VI. Note that due to the Z 2 line gap classification in 1D, flux defects only probe the Z 2 nature of the 3D Z invariant in NH class D.
f. Class D S + + I †,S − Models in NH class D S+ possess PHS U P (U P U * P = +1) and sublattice symmetry S, with [U P , S] = 0. As outlined in App. A 2, the corresponding EHH introduces a chiral symmetryΣ C , which combines with sublattice symmetry to form a unitary symmetryŪ =SΣ C with U 2 = +1. The unitaryŪ satisifies {Ū P ,Ū } = [S,Ū ] = [Σ C ,Ū ] = 0. We may furthermore define a TRSŪ T = U PΣC which satisfiesŪ TŪ * T = −1 and anti-commutes withŪ . Hence the eigenspaces ofŪ are not independent and individually enjoy a chiral symmetry, leading to Hermitian class AIII with a Z classification in 3D [83]. Without additional crystalline symmetries, the zero energy modes arising in the flux Dirac theory of systems in Hermitian class AIII are not protected (see App. B 3 b). Adding a crystalline symmetry I †,S− , however, quantizes the flux and allows for a stable flux response. The results derived in App. B 3 b show the presence of one zero energy mode per flux under OBC, occurring in eachŪ subspace.
Away from E 0 = 0, we can only form the TRS . This results in Hermitian class CI which does not protect a Kramers degeneracy. Consequently, the pair of flux modes can be now gapped, yielding a higher order response only at E 0 = 0. Models in NH class D S+ + I †,S− therefore show a higher-order Hermitian flux response, the higher-order flux Majorana mode, as introduced in Sec. VI.
Hence we obtain Hermitian class DIII in eachŪ subspace, giving a Z ⊕ Z classification in 3D [83]. By modding out line-gap phases, this is reduced to a Z classification where the nontrivial element corresponds to having only a singlē U subspace being nontrivial [39]. Hence the full Hermitian spectrum hosts a helical metal under π-flux insertion (see App. B 3 c).
Away from E 0 = 0, we can only form the TRS . This results in Hermitian class DIII, such that the flux-bound helical metal cannot move away from zero energy. Hence the entire point gap of a corresponding NH system fills with flux tube localized modes. Consequently, models in NH class D S− quantize the flux to ϕ = 0, π and show a NH flux spectral jump, discussed in Sec. IV. Note that due to the Z 2 -classification of NH class D S− in 2D, flux defects only probe the Z 2 nature of the 3D Z invariant.
h. Class DIII S +− As for the case of 2D, models in NH class DIII S+− possess TRS, PHS and chiral symmetry with (U T U * T , U P U * P ) = (−1, 1), where U C = U T U P , as well as sublattice symmetry S. As outlined in App. A 2, the corresponding EHH at E 0 = 0 obtains TRSŪ T , PHS U P and three chiral symmetriesŪ C ,Σ C ,S. The chiral symmetries can be combined to two commuting unitary symmetriesŪ 1 =Ū CΣC withŪ 2 1 = −1 andŪ 2 =SΣ C withŪ 2 2 = +1. SinceŪ 1 has an imaginary spectrum and {Ū T , the eigenspaces ofŪ 1 are not independent and individually enjoyŪ T andŪ 2 symmetry. Moreover, we have so that theŪ 2 eigenspaces are independent and individually preserveŪ T symmetry. They therefore lie in Hermitian class AII and yield a Z 2 ⊕ Z 2 classification in 3D [83] that is reduced to Z 2 by line-gap phases [39]. The nontrivial element corresponds to having only a singleŪ 2 subspace nontrivial. In App. B 3 a, we derived the flux response of EHHs in Hermitian class AII, which shows a helical metal pinned to the flux defects.
Away from E 0 = 0, we can only form the TRS . This results in Hermitian class DIII, such that the flux-bound helical metal cannot move away from zero energy. Hence the entire point gap of a corresponding NH system fills with an extensive number of states at a single momentum k ∥ . Consequently, models in NH class DIII S+− quantize the flux to ϕ = 0, π and show a NH flux spectral jump, discussed in Sec. IV. We may furthermore define a PHSŪ P =Ū TΣC which satisfiesŪ PŪ * P = −1 and commutes withŪ . Hence the eigenspaces ofŪ are independent and individually enjoy TRS, PHS and chiral symmetry, leading to Hermitian class CII with a Z 2 classification in each sector in 3D [83]. Without additional crystalline symmetries, the zero energy modes arising in the flux Dirac theory of systems in Hermitian class CII are not protected (see App. B 3 e). Adding a crystalline symmetry I †,S−,T− , however, quantizes the flux and allows for a stable flux response. Introducing a π-flux results in the presence of zero energy states under OBC along the flux tubes, in addition to the surface signature for ϕ = 0.
Away from E 0 = 0, we are only left with the chiral symmetryΣ C , resulting in Hermitian class AIII. As Her-mitian class AIII with additional crystalline symmetry I †,S−,T− protects flux modes localized at the end of the flux tubes, these zero modes appear extensively in the point gap of the NH model, constituting a higher-order flux skin effect (see Sec. V) for NH class AII S+ +I †,S−,T− .
j. Class C Models in NH class C possess a PHS U P with U P U * P = −1, yielding a 2Z classification in 3D, while lower dimensions are trivial. As outlined in App. A 2, the corresponding EHH introduces a chiral symmetrȳ Σ C , which combines with PHS to a TRSŪ T =Ū PΣC which satisfiesŪ TŪ * T = +1. Hence we obtain Hermitian class CI, giving a 2Z classification in 3D [83]. The flux response of the EHH in Hermitian class CI is trivial, without protected flux localized zero modes (see App. B 3 d). Consequently, systems in NH class C do not show flux induced states within the NH point gap. Therefore, the flux response for this NH symmetry class is trivial.
k. Class C S + + I †,S − Models in NH class C S+ possess PHS U P (U P U * P = −1) and sublattice symmetry S, with [U P , S] = 0. As outlined in App. A 2, the corresponding EHH introduces a chiral symmetryΣ C , which combines with sublattice symmetry to form a unitary symmetryŪ =SΣ C with U 2 = +1. The unitaryŪ satisifies {Ū P ,Ū } = [S,Ū ] = [Σ C ,Ū ] = 0. We may furthermore define a TRSŪ T = U PΣC which satisfiesŪ TŪ * T = +1 and anti-commutes withŪ . Hence the eigenspaces ofŪ are not independent and individually enjoy a chiral symmetry, leading to Hermitian class AIII with a Z classification in 3D [83]. Without additional crystalline symmetries, the zero energy modes arising in the flux Dirac theory of systems in Hermitian class AIII are not protected (see App. B 3 b). Adding a crystalline symmetry I †,S− , however, quantizes the flux and allows for a stable higher order flux response, with a single mode at the end of each flux tube perŪ eigenspace.
Away from E 0 = 0, we can only form the TRS . This results in Hermitian class DIII, such that the two flux modes in eachŪ eigenspace are still protected. Consequently, these zero modes appear extensively in the point gap of the NH model. This constitutes a higher order skin effect, with skin modes appearing at the ends of the flux tubes in a finite geometry. The precise localization of skin modes is then determined by the present crystalline symmetry I †,S− . This is the defining signature of the NH higher-order flux skin effect, discussed in Sec. V. l. Class C S − NH class C S− possesses PHS U P (U P U * P = −1) and sublattice symmetry S, with {U P , S} = 0. As outlined in App. A 2, the corresponding EHH introduces a chiral symmetryΣ C , which combines with sublattice symmetry to form a unitary symmetryŪ =SΣ C withŪ 2 = +1. The unitaryŪ satisifies [Ū P ,Ū ] = [S,Ū ] = [Σ C ,Ū ] = 0. We may furthermore define a TRSŪ T =Ū PΣC which satisfiesŪ TŪ * T = +1, commuting withŪ . Hence we obtain Hermitian class CI in eachŪ subspace, giving a 2Z ⊕ 2Z classification in 3D [83]. By modding out linegap phases, this is reduced to a Z classification where the nontrivial element corresponds to having only a singleŪ subspace being nontrivial [39].
Away from E 0 = 0, we can only form the TRS

Scaling of flux core-localized modes in 2D
Sec. III predicts a skin effect localizing modes at the points where flux cores with ϕ = π pierce a 2D system. Their number scales with the system size, precisely with the extension of the system parallel to the line connecting the two flux cores L ⊥ , as shown in Fig. 8. All states in this 1D slice collapse towards the flux cores. For small separations d, states at the two flux cores can hybridize, yielding a smaller number of localized states. The topologically relevant regime has both flux cores well separated, for which the number of skin modes does not change with d. Consequently, all L ⊥ states experience the flux skin effect, irrespective of the extension of the system to the left and right of the flux defects.

Scaling of flux tube-localized modes in 3D
In 3D, NH flux defects probe 4 unique topological phases: the flux spectral jump, the higher-order flux skin effect, the NH flux Majorana mode and its higher-order cousin. For a π-flux, the flux spectral jump causes an extensive number of defect localized modes. Their number scales with the extension of the system parallel to the line connecting the two flux tubes L ⊥ , as shown in Fig. 9a. Conversely, the higher-order flux skin effect appears along π-flux tubes, thereby localizing O(L ∥ ) modes, where L ∥ is the extension of the system along the defect as shown in Fig. 9b. In contrast, flux Majorana modes appear only as isolated modes, without a connected skin effect. Their number then derives from the EHH Dirac theory, which predicts two zero-energy modes per flux tube, localized at the flux defects (see App. B 2 d). Correspondingly, the model in NH class D hosts one mode per flux tube (see Fig. 9c  This appendix provides details on the fractional nature of the higher-order flux skin effect and the role of (pseudo-)inversion symmetry in it. Additionally, we provide details for NH symmetry classes AII S+ and C S+ .

Fractional nature of the higher-order flux skin effect
The NH higher-order flux skin effect results from the presence of one flux-end localized zero energy mode per flux tube in the corresponding EHH. The single flux bound state is fractional, as there exists no 1D system with a single end state. In this appendix, we investigate the fractional nature of this state by considering the weight of a given chirality on each surface, which reveals an imbalance only present for finite flux ϕ = π.
The eigenstates of the EHH are given as where L ⊥ , L ∥ are the system dimensions, with L i describing internal (sublattice, spin) degrees of freedom. We can then consider an interval around zero energy −ϵ < E j < ϵ: requiring that we do not separate degenerate states, the set of states in the interval 2ϵ forms the subset E ⊂ Ξ, for instance as highlighted in the inset of Fig. 5b. Expressing the chiral symmetry of the EHH Σ C in this subspace, allows to diagonalize Σ Σv ± l = ±v ± l , l ∈ {1, . . . , n ± }, and obtain eigenstates of the chiral symmetryΣ C . We denote eigenvectors with positive (negative) chiral eigenvalue as v + l (v − l ) and their number as n + (n − ). We are interested in the chirality of the original eigenstates of the EHHH, which we expand in the obtained chiral basis: where |ψ ± ⟩ are now eigenstates of positive (+) or negative (-) chirality. For a system without flux defects, for instance in Hermitian class AIII as considered in Sec. V, each surface has zero net-chirality [83]. We calculate the net chirality by summing all states of (+)/(-) chirality for the top/bottom surface, ρ ± top,bottom , and investigating their differences. Specifically, we use where we restrict the summation for bottom (top) surface, denoted by r ∥ = 0 (r ∥ = L ∥ ), once the states ψ ± decay in the bulk (|ψ ± (r ∥ )| 2 < δ, δ ≪ 1). Forming the differences between top and bottom surface for a fixed chirality, e.g. (+), yields In the absence of a flux defect, ϕ = 0, each chirality has equal weight one each surface, ∆ ± = 0. Conversely, for a nontrivial flux ϕ = π, each chirality has predominant weight on one surface. For a PBC system with two flux tubes, ∆ ± is equal to one: there is a single state per chirality per flux tube (see Fig. 10). Viewing a flux tube as an effective 1D model highlights the fractional nature of this response: A 1D model always contains end states of both chiralities, for instance as in the Su-Schrieffer-Heeger model [85]. The fact that we obtain a single state per chirality shows the fractional nature in the EHH, causing the higher-order flux skin effect in the NH system.

Role of inversion symmetry in the higher-order flux skin effect
This appendix investigates the effect of (pseudo-)inversion symmetry on NH symmetry classes with a higher-order flux skin effect under π-flux insertion. The relevant NH symmetry classes are A, A S , AII S+ and C S+ . In order to obtain a nontrivial response to flux defects in these classes, we require the presence of additional crystalline symmetries, restricting the magnetic flux ϕ to 0, π in PBC, while still allowing for nontrivial topology. Specifically, we investigate inversion symmetry appearing in the EHH as and pseudo-inversion symmetry appearing in the EHH as We have to identify choices of (pseudo-)inversion symmetry that leave the intrinsic point gap classification invariant, i.e. do not trivialize the flux response.

a. Class A
In NH symmetry class A, the EHH can be formed as outlined in Eq. (4) of the main text, which adds a chiral symmetryΣ C , resulting in Hermitian class AIII. We can either use an inversion symmetry (D8) or a pseudoinversion symmetry (D10) to fix the flux ϕ to 0, π in PBC in the NH Hamiltonian. However, these two choices differ in their commutation with the chiral symmetryΣ C when considering the EHH: Whereas "normal" inversion symmetry commutes with chiral symmetry pseudo-inversion symmetry anticommutes, This means both chiral subspaces have the same inversion eigenvalues. Consequently, at every inversionsymmetric momenta in the Brillouin zone, the occupied and unoccupied subspace have the same number of positive and negative inversion eigenvalues. Therefore every Hermitian model with these symmetry properties cannot have a band inversion and so is topologically trivial [86].
As topological zero modes of Hermitian extended models stand in one-to-one correspondence with NH modes within a point-gapped bulk (see Sec II B), also the corresponding NH model is trivial. Accordingly, Ref. 87 outlines that with commuting inversion symmetry, the resulting Hermitian class AIII I+ is trivial, in agreement with the vanishing 3D winding number of the NH system. Conversely, anticommuting pseudo-inversion symmetry yields a EHH in Hermitian class AIII I− which is classified by a Z index. Contrary to normal inversion, pseudo-inversion hence allows for a nontrivial flux response, the higher-order flux skin effect.

b. Class A S
In NH symmetry class A S , the EHH can be formed for each sublattice symmetry block (see App. B 5 d) as outlined in Eq. (4). This construction adds a chiral symmetryΣ C , resulting in Hermitian class AIII in each block, of which only one has to be nontrivial. Again, we can either use an inversion symmetry (D8) or a pseudo-inversion symmetry (D10) to fix the flux ϕ to 0, π in PBC in the NH Hamiltonian. However, these two choices differ in their commutation with the chiral symmetryŪ C when considering the EHH: Whereas "normal" inversion symmetry commutes with chiral symmetry pseudo-inversion symmetry anticommutes (see Sec. D 2 a). Additionally, the sublattice symmetryS of the NH Hamiltonian acts as a chiral symmetry ofH, too. Depending on the commutation or anti-commutation of bothŪ C andS with (pseudo-)inversion symmetry the topological classification differs. Ref. 87 outlines that both chiral symmetries have to anti-commute with inversion symmetry to form a nontrivial phase classified by a Z index. Correspondingly, an inversion symmetry in the NH Hamiltonian always results in a trivial phase. Pseudo-inversion symmetry on the other hand only results in a nontrivial flux response for an anti-commuting sublattice symmetry S. Contrary to normal inversion, pseudo-inversion hence allows for a nontrivial flux response, the higher-order flux skin effect.
c. Class AII S + In NH symmetry class AII S− , the EHH follows as outlined in Eq. (4) (see Sec. B 5 i). The EHH possesses a unitary symmetryŪ , whose eigenspaces are independent and individually enjoy TRS, PHS and chiral symmetry, leading to Hermitian class CII in each sector. A stable flux response in this class requires the presence of a crystalline symmetry: we can either use an inversion symmetry (D8) or a pseudo-inversion symmetry (D10). We have to analyse these two choices with respect to their commutation with TRS and PHS, to derive their topological classification [87]. In the EHH we can form two TRS:Ū Similarly, we can form two PHS, A 4Z classified phase is obtained if both TRSs commute while the PHSs anticommute with inversion. Conversely, if both TRSs anticommute while the PHSs commute with inversion, we obtain a Z 2 classification [87]. Only the latter is compatible with the Z 2 -point gap classification without crystalline symmetries. All other cases, also including inversion symmetry, do not yield a nontrivial response. Therefore, we should choose a pseudo-inversion symmetry commuting with sublattice symmetry S, but anticommuting with TRS U T . Contrary to normal inversion, pseudo-inversion hence allows for a nontrivial flux response, the higher-order flux skin effect.
d. Class C S + In NH symmetry class C S+ , the EHH possesses chiral symmetriesS,Σ C in each independent unitary eigenspace (see App. B 5 k). Consequently, the EHH sits in Hermitian class AIII for each block, of which only one has to be nontrivial. We can either use an inversion symmetry (D8) or a pseudo-inversion symmetry (D10) 11. Higher-order flux skin effect in NH class AII S + and C S + . a The EHH for models in NH class AII S + [Eq. E17] shows a gapless surface even for ϕ = 0. b Under OBC in NH class AII S + , the presence of a nontrivial flux ϕ = π introduces four additional flux states. c Similarly, the EHH for models in NH class C S + [Eq. E19] shows a gapless surface for ϕ = 0. d Under OBC, however, the presence of a nontrivial flux ϕ = π introduces two additional flux states per unitary subspace, yielding in total four additional modes. All panels are generated for systems of size 20 × 20 × 20 unit cells.
to fix the flux ϕ to 0, π in the NH Hamiltonian. However, these two choices differ in their commutation with the chiral symmetryΣ C when considering the EHH: Whereas "normal" inversion symmetry commutes with chiral symmetry, pseudo-inversion symmetry anticommutes (see Sec. D 2 a). Additionally, the sublattice symmetryS of the NH Hamiltonian acts as a chiral symmetry ofH, too. Depending on the commutation or anti-commutation of bothΣ C andS with (pseudo-)inversion symmetry the topological classification differs. Ref. 87 outlines that both chiral symmetries have to anti-commute with inversion symmetry to form a nontrivial phase classified by a Z index. Correspondingly, an inversion symmetry in the NH Hamiltonian always results in a trivial phase. Pseudo-inversion symmetry on the other hand only results in a nontrivial flux response for an anti-commuting sublattice symmetry S. Contrary to normal inversion, pseudo-inversion hence allows for a nontrivial flux response, the higher-order flux skin effect. 3. Resolving the higher-order flux skin effect in the EHH spectrum In Sec. V, we derived the presence of a higher-order flux skin effect in NH symmetry classes A, A S , AII S+ , and C S+ . While being observable in the NH system due to the extensive localization of states, the occurrence in the corresponding EHH is more intricate: the inherent nontrivial surface state of the EHH obscures flux-localized modes. We can nevertheless cleanly identify the presence of flux induced zero-energy modes in the EHH by considering an interval around zero energy −ϵ < E < ϵ, within the bulk energy gap (|ϵ| < E bulk [88]). We denote the set of states within this interval by E (see Fig. 5b and Fig. 11 for examples). The number of states |E| can be used to resolve the higher-order flux skin effect.
NH class A and A S map to Hermitian class AIII, where the total number of states |E| has to be a multiple of 4 in the absence of a magnetic flux ϕ = 0, |E| mod 4 = 0. (D18) This can be understood from the fact that the EHH in Hermitian class AIII hosts an integer number of Dirac cones on each surface [83]. Since there is no net chirality per surface [89], energy eigenvalues appear in pairs (E, −E) on each surface. Inversion symmetry maps between the two surfaces, such that E contains a multiple of 4 states. On the other hand, for a finite flux ϕ = π, we obtain one additional zero-energy mode per flux tube (see App. B 3 b), and so we find for a PBC system hosting two flux tubes that |E| mod 4 = 2. (D19) Calculating |E| mod 4 for our example system in NH class A (Fig. 5b) yields |E| ϕ=0 mod 4 = 0 (upper left inset), while we find |E| ϕ=π mod 4 = 2 as predicted (lower right inset). This relation is modified for the NH classes AII S+ and C S+ . NH class AII S+ maps to Hermitian class CII⊗CII, of which only one subspace is nontrivial (see App. B 5 i). Hermitian class CII shows a 4-fold degenerate Dirac cone per surface, as chiral symmetry combined with TRS yields pairs of doubly degenerate states. As inversion symmetry maps between the two surfaces, states appear as multiples of eight in E. The introduction of two flux tubes amends this by 4 additional states (see App. B 3 e), hence we obtain for the number of states in E, denoted by |E|, that |E| mod 8 = 0 (4), for flux ϕ = 0 (π). Fig. 11a and b show the comparison of both cases and the validity of the above index. NH class C S+ maps to Hermitian class AIII in two interdependent unitary subspaces (see App. B 5 k). As outlined in Sec. V, each subspace comes with a multiple of 4 states. The full system hence has E containing multiples of 8 states. The introduction of two flux tubes amends this by 2 additional states per unitary subspace (see App. B 3 b), hence |E| mod 8 = 0 (4), for flux ϕ = 0 (π). Fig. 11c and d show the comparison of both cases and the validity of the above index.
Appendix E: Toy models for all NH symmetry classes with nontrivial flux response

2D models
In this appendix we present models for all intrinsically nontrivial 2D point-gapped NH classes showing a nontrivial flux response. As outlined in the main text, we consider flux defects oriented along the x-direction.

a. Class AII †
The model in NH class AII † is given by the Hamiltonian with the Pauli matrices σ µ (µ = 0, x, y, z), possessing pseudo TRS as U T = iσ y and ∆ ̸ = 0 is a real parameter breaking residual symmetries.

b. Class DIII †
The nontrivial point-gapped model in NH class DIII † is given by the Hamiltonian with the Pauli matrices σ µ (µ = 0, x, y, z). We obtain TRS † U T = σ y , PHS † U P = σ x and chiral symmetry U C = σ z , fulfilling the required commutation relations.

d. Class D S −
A nontrivial point-gapped model in NH class D S− is realized by the Hamiltonian with Q(k; µ) =i sin(k x )σ z − sin(k y )σ 0 with the Pauli matrices σ µ and τ µ (µ = 0, x, y, z). The Hamiltonian possesses PHS U P = τ x σ 0 and sublattice symmetry S = τ z σ 0 , fulfilling the required commutation relations. Note that ∆ multiplies a term to remove unwanted residual symmetries.

3D models
In this appendix we present models for all intrinsically nontrivial 3D point-gapped NH classes showing a nontrivial flux response. As outlined in the main text, we consider flux defects oriented along the x-and zdirection.
where the Pauli matrices σ µ and τ µ act on the spin and orbital degrees of freedom, respectively, with µ = 0, x, y, z and the 0-th Pauli matrix as the 2 × 2 identity matrix. ∆ ̸ = 0 is a real parameter breaking residual symmetries. The invariant is w 3D = 1 for 3 − δ/λ < M < 3+δ/λ. By changing M , we transition to the trivial phase at M = 3 ± δ/λ. The Hamiltonian has pseudo-inversion symmetry with I = τ z σ 0 and pseudo TRS represented by U T = τ 0 σ y .

b. Class A+I †
The prototypical phase in 3D NH class A is formed by the exceptional exceptional topological insulator introduced in Ref. 21: where the Pauli matrices σ µ and τ µ act on the spin and orbital degrees of freedom, respectively, with µ = 0, x, y, z and the 0-th Pauli matrix as the 2 × 2 identity matrix. ∆ ̸ = 0 is a real parameter breaking residual symmetries. The Hamiltonian is only left with pseudoinversion symmetry with I = τ z σ 0 . c. Class A S + I † The nontrivial model in NH class A S is given by with the Pauli matrices σ µ and τ µ (µ = 0, x, y, z), ∆ ̸ = 0 a real parameter breaking residual symmetries and Q(k; µ) =i sin(k x )σ x + i sin(k y )σ y + i sin(k z )σ z +   i=x,y,z H A S has a sublattice symmetry S = τ z σ 0 and pseudoinversion symmetry I = τ x σ 0 .
with the Pauli matrices σ µ (µ = 0, x, y, z), possessing a PHS U P = σ 0 . e. Class D S + + I † The nontrivial point-gapped phase in NH class D S− is based on the Hamiltonian in NH class D, formed by with the Pauli matrices σ µ and τ µ (µ = 0, x, y, z) and ∆ ̸ = 0 a real parameter breaking residual symmetries. We obtain PHS U P = τ 0 σ 0 and sublattice symmetry S = τ z σ 0 , fulfilling the required commutation relations. Additionally, H D S + has pseudo-inversion symmetry with I = τ y σ y .

f. Class D S −
The nontrivial point-gapped phase in NH class D S− is formed by with the Pauli matrices σ µ and τ µ (µ = 0, x, y, z) and Q(k; µ) =i sin(k x )σ x + i sin(k y )σ y + i sin(k z )σ z +   i=x,y,z cos(k i ) − µ   σ y .