Universal non-Hermitian skin effect in two and higher dimensions

Skin effect, experimentally discovered in one dimension, describes the physical phenomenon that on an open chain, an extensive number of eigenstates of a non-Hermitian Hamiltonian are localized at the end(s) of the chain. Here in two and higher dimensions, we establish a theorem that the skin effect exists, if and only if periodic-boundary spectrum of the Hamiltonian covers a finite area on the complex plane. This theorem establishes the universality of the effect, because the above condition is satisfied in almost every generic non-Hermitian Hamiltonian, and, unlike in one dimension, is compatible with all point-group symmetries. We propose two new types of skin effect in two and higher dimensions: the corner-skin effect where all eigenstates are localized at corners of the system, and the geometry-dependent-skin effect where skin modes disappear for systems of a particular shape, but appear on generic polygons. An immediate corollary of our theorem is that any non-Hermitian system having exceptional points (lines) in two (three) dimensions exhibits skin effect, making this phenomenon accessible to experiments in photonic crystals, Weyl semimetals, and Kondo insulators.


A. Some numerical examples of the theorem
In this subsection, we provide some numerical examples of the theorem. For simplicity, we consider the following single-band model where t ij represents the hopping parameter between the site (m, n) and the site (m, n) + (i, j). The two rows in Fig. 1 represent two different models. In the first one, the hopping parameters are shown in Fig. 1  (a), and the periodic boundary spectrum, i.e. E(k) = H(k) with k ∈ BZ, is shown in Fig. 1 (b). Here the k 200 × 200 is used. One can notice that the spectral area of the first model is nonzero. As a result, the open boundary eigenstates show the localization behaviors, namely, the emergence of skin effect, as shown in Fig. 1 (c-d). In order to illustrate the localization properties, is plotted in Fig. 1 (c-d), , where ψ n (x) is the n-th normalized eigenstate of the Hamiltonian with open-boundary condition, and N is the number of eigenstates. For the second example, since the spectral area is zero, as shown in Fig. 1 (f), there is no skin effect. Indeed, W (x) shown in Fig. 1 (g-h) are uniform distributed on the lattice. The results of the numerical calculation support our theorem, and then we will formally prove the theorem.

B. The proof of the theorem in two-dimensions
In this subsection, we will prove the theorem in two-dimensional systems. As illustrated in Fig. 2(a), the proof is divided into three steps. Therefore, we partition this subsection into three parts, and each part corresponds to one step of the proof.

1.
Step I : spectral area and spectral winding We begin with a two-dimensional single-band non-Hermitian model with periodic boundary in both x and y directions where u and v are real functions about k = (k x , k y ). For any k r ∈ BZ, one can define the following winding number where Γ kr represents the infinitesimal counterclockwise loop enclosing k r . Here E r = H(k r ) represents the reference energy, which is shown in Fig. 2 (b) with red point. We note that for different k r , the reference energy is different. This topological invariant describes the spectral winding on the complex plane. As shown in Fig. 2 (b), if ν(k r ) is nonzero, the image of Γ kr , i.e. H(Γ kr ), forms a closed loop that encloses E r . First we prove the following two statements, 1. if there are some generic k r points in the BZ with nonzero topological charge, the spectral area must be nonzero; 2. if all the k r points in the BZ have zero topological charge, the spectral area must be zero.
The above two statements establishes the equivalence relation between "spectral area" and "topological charge of k point". Based on the definition of winding number, the statement 1 is obvious. Therefore, we only need to prove the statement 2. In order to show this, we expand the Hamiltonian at the point k r ≡ (k r x , k r y ) as follows H(k) − H(k r ) ≈ ∂ x H(k r )q x + ∂ y H(k r )q y ,

Nonzero Spectral Winding of Straight Lines in the BZ
Step I Step II Step III where q x , q y are the displacements k − k r in x, y directions respectively. For a generic k r point, the first derivative of H does not vanish, i.e., the coefficients of q x and q y in Eq. (5) cannot be zero at the same time. The reason is that in order to make ∂ x H(k r ) = ∂ y H(k r ) = 0, four independent real conditions, Re[∂ x H(k r )] = Im[∂ x H(k r )] = Re[∂ y H(k r )] = Im[∂ y H(k r )] = 0, need to be satisfied. However, in two dimensions, there are only two free parameters (k x , k y ), which cannot satisfy the above four equations generally.

Conjecture
In order to calculate the topological charge of a generic k r point, according to Eq. (3), one can define where the notation ∂ x/y refers to ∂/∂ k x/y . When det[C(k r )] = 0, the topological charge of k r is the sign of the determinant of C(k r ), expressed as ν(k r ) = sign[det[C(k r )]].
Therefore, a sufficient and necessary condition for the zero charge of every k ∈ BZ (the statement 2) is A theorem (the corollary of theorem 13.2) in Ref.
[1] tells us that if C(k) = 0 and det[C(k)] = 0 for an open set S, then u(k) and v(k) have a functional dependent relation with k ∈ S ⊂ BZ, that is, u(k) = g(v(k)) with g being a complex function. Applying the theorem to the entire BZ (except for some isolated points where the first derivative vanishes), one can reexpress the single-band Hamiltonian that satisfies Eq. (8) as where h(k) is a real and periodic function of k, and P is a complex polynomial of h. Since h(k) is a real periodic function, its image must be an arc on the real axis, e.g. h(k) ∈ [h 1 , h 2 ], where h 1/2 are real numbers. Therefore, the image of P [h(k)] must be an arc on the complex plane, that is to say, the spectral area is zero. This completes the proof of the statement 2 for single-band models. Generalizing the above discussion to the multi-band case, for each k r ∈ BZ, the topological charge defined for the m-th band is where E m (k r ) is the energy of the m-th band with the momentum k r . For the second equal sign in Eq. (10), we have used det[H(k) − E m (k r )] = n [E n (k) − E m (k r )]. For Eq. (10), if k r is not the degeneracy point, only n = m term in the summation has contributions to the topological charge. Therefore, the Eq. (10) further becomes Using the similar approaches in the single-band case, one can conclude that, the real and imaginary parts of E m (k) are locally functional dependent on the neighborhood of k ∈ BZ. As a result, the spectrum of E m (k) must be an arc. The above conclusion applies for each band. So far, we have proven the equivalence relation between "spectral area" and "topological charge". In following contents, we will prove that if there are some generic k points in the BZ with nonzero topological charge, then there must be nonzero spectral winding number of the straight lines in the BZ.
We first notice that the BZ in two-dimensional systems can be covered by a set of straight lines of any slope, labeled as {L s }. Here, the subscript s indicates the slope of the set {L s }, and L s represents a generic straight line belonging to {L s }. For example, if we fix k y (k x ) and change k x (k y ) from 0 to 2π, we obtain a horizontal (vertical) straight line with the slope s being 0 (∞) in BZ, and the set of all horizontal or vertical straight lines ({L 0 } or {L ∞ }) covers the entire BZ. Particularly, an inclined straight line goes out from one side of BZ and again enters from another side as shown in Fig. 2(b). Since the straight lines on the BZ are periodic, one can define the spectral winding number for each straight lines with respect to the prescribed reference energy E r .
Obviously, if all the k points on the BZ have zero topological charge, ν(L s , E r ) must be zero for arbitrary L s and E r . Otherwise, one can always find some L s , such that ν(L s , E r ) is nonzero. Next we briefly prove the latter statement.
Assuming that each straight line in {L 0 } and {L ∞ } has zero spectral winding, and there is a k r point carrying nonzero topological charge on the BZ. For a generic inclined straight line, one can always find corresponding horizontal and vertical straight lines, such that together with the inclined straight line to form a closed path enclosing k r in BZ. Therefore, the closed path has nonzero winding number with respect to E r . Due to the zero spectral winding of the horizontal and vertical straight lines as we assumed, hence a generic inclined straight line must have nonzero spectral winding number. Let's take an example to show this. Consider a simple single-band model with the Hamiltonian H(k x , k y ) = 2 cos k x + 2i sin k y .
As shown in Fig. 2, the spectrum of the Hamiltonian along each horizontal or vertical straight line (gray lines) in the BZ has zero spectral winding number with respect to any reference energy. However, if we choose the straight line in the dark cyan color, these three straight lines (two gray lines and the darker cyan line) together form a closed path that encloses k r . Therefore, the closed path has nonzero spectral winding number with respect to E r . Due to the zero spectral winding of two gray lines, the dark cyan straight line must have nonzero winding number regarding E r as illustrated in Fig. 2(b).
In conclusion, we have proven the equivalence between "spectral area" and "spectral winding number of straight lines in the BZ". Up to this point, we have completed the first step of the proof shown in Fig. 2(a).

2.
Step II : spectral winding and skin effect on the stripe geometry Before, we have proved that if the spectral area is zero, the spectral winding number along each direction in the BZ is zero, and vice versa; if the spectral area is nonzero, there is at least one direction along which the spectral winding number is nonzero. Now we transform the momentum basis from k = (k x , k y ) T to q = (q , q y ) T by where S is a 2 * 2 integer matrix due to the lattice momentum in the BZ is discrete. For any direction in the BZ, one can always choose the appropriate q basis such that under this q basis, each straight line in this direction can be obtained by fixing q and running q y from 0 to 2π. Accordingly, the Hamiltonian H(k x , k y ) can be transformed intõ H(q , q y ). For any fixed q , the spectral winding can be expressed as where E b is the reference energy. According to what we have proved, if the spectral area is zero, under any transformation S, the spectral winding w E b (q ) = 0 for any q and E b . It means that the Hamiltonian on any stripe geometry does not exhibit skin effect if spectral area is zero. Here the stripe geometry refers to the geometry with open boundary in only one direction and periodic boundary in other direction (preserves the momentum q ). If the spectral area is nonzero, there is at least one stripe geometry on which the Hamiltonian establishes skin effect. Therefore, we have completed the second step of the proof.

3.
Step Note that the above is the equivalent statement of the conjecture in step III of Fig. 2(a), that is, skin effect on stripe geometry implies skin effect on fully open-boundary geometry (i.e., the universal skin effect).
Next we explain why the conjecture makes sense, and take some numerical results to support our conjecture.
• First, we believe that whether a given edge exhibits a skin effect or not depends only on the bulk spectral topology in the direction perpendicular to this edge, having nothing to do with other boundary conditions, which accords with the spirit of bulk-boundary correspondence.
• Second, for a given edge, no skin effect means that the wave packet obeys the conventional law of reflection on this edge (as shown in Fig.4(c) of the main text). Otherwise, the skin effect leads to the anomalous dynamical behavior on this edge (as shown in Fig.4(d)). This local physical consequence caused by the skin effect can not be affected by any change to other edges in the thermodynamic limit. Therefore, we conjecture that if a given edge under stripe geometry does not exhibit the skin effect, it still does not show the skin effect under any open boundary geometry, and vice versa.

No Skin Effect
Supplementary • Third, we show the numerical results in Fig. 3 to support the conjecture. The bulk Hamiltonian of this model reads The mirror-y symmetry of the bulk Hamiltonian makes the zero spectral winding for any a k x -subsystem. Based on the first two steps of the proof and the above conjecture, the theorem of the universal skin effect can be obtained.

C. The proof of the theorem in three-dimensions
In this section, we extend the above proof of two-dimensional systems into three dimensions, and obtain the conclusion that nonzero spectral area means the existence of the universal skin effect.
Consider a general three-dimensional single-band tight-binding Hamiltonian, which consists of real-and imaginarypart functions We choose a generic k r point and use its energy H(k r ) as the reference energy. For a given reference energy E r , we can obtain a one-dimensional curve in the there-dimensional BZ by solving the following two real equations, Each equation determines a surface, and the intersection of two surfaces is one-dimensional curve in there-dimensional BZ. The tangent direction of the curve at k r is perpendicular to the normal vector of the two surfaces at this k r point. The tangent vector at k r is expressed as where ∇u(k r ) represents the gradient of u. We choose the local coordinate system (R 3 space) with k r as the origin, and the gradient is reexpressed as where q i ≡ ( k−k0 |k−k0| ) i , k and k 0 represent two vectors in the global coordinate system. Next we expand the Hamiltonian into Taylor series around the origin of the local coordinate system, where the subscription i represents the partial differential to x, y, z. And q i represents the deviation of k i from k r,i and the last term is the infinitesimal of higher order of |q|. Obvious, the zero winding condition requires or equivalently, Next, we prove that if all k points in three-dimensional BZ satisfy T k = 0, then the entire 3D periodic-boundary spectrum must be an arc in the complex plane. We define a two-tuple function W (k) := [u(k) v(k)] t with three variables, the exterior derivative of the vector-valued function is expressed as Eq. (23) implies that the rank of dW less than 2 (the number of components of W ). To be precise, there are the following cases. (i.) Both the gradients of u and v are not zero vector, and they are linearly dependent on each other. (ii.) One of the gradients of u and v is zero vector. (iii.) Both the gradients of u and v are zero vector. In all these cases, we can obtain the final conclusion that u and v are linearly functional dependent on each other. Therefore, the spectrum must be arcs on the complex plane. A 3D BZ can be divided into a series of plane systems, each plane corresponds to a two-dimensional subsystem. If the spectral area of a three-dimensional system is nonzero, then for each reference energy on the spectral area, its preimage (1D ring) has nonzero topological charge. Equivalently, the two-dimensional subsystem, of which the BZ (2D plane) has intersections with the ring, also has nonzero topological charge for the intersecting k points. Hence, the 2D subsystem has nonzero spectral area, and has the universal skin effect. Correspondingly, we come to the same conclusion in 3D systems that nonzero spectral area implies the presence of the universal skin effect. Figure 4.
(a) shows the periodic-boundary spectrum of Eq. (28) with gray color, and the pre-images of E0 = 1 + i (red point in (a)) are the four red points in (b). The periodic-boundary spectrum of Eq. (29) is the gray line in (c), and k(E0 = 1 + i) is plotted by the red lines in (d).

SUPPLEMENTARY NOTE 2: A PHYSICAL EXPLANATION FOR THE THEOREM
Here, we use some examples to illustrate the intuition that motivates the theorem. Consider the following onedimensional model placed on a chain of length L. Under the periodic boundary condition, the two Bloch waves e ik0x and e −ik0x have the same energy E(k 0 ) = 2 cos k 0 . When the system has open boundary condition, the Bloch wave e ik0x will be reflected to e −ik0x with a π-phase shift. Their linear superposition e ik0x − e −ik0x is an eigenstate with energy 2 cos k 0 that satisfies the zero boundary condition at x = 0, L, thus being an open-boundary eigenstate. When the system is added a momentum-dependent dissipation, the spectrum E(k) becomes complex and forms an ellipse in the complex plane. In this case, the degeneracy is broken, e.g. E(k) = E(−k), which implies the open-boundary eigenstates are no longer the linear superposition of the extended Bloch waves. This implies the emergence of skin effect. Extend the above arguments to two dimensions, and we can provide a physical explanation for the theorem proved in the Supplementary Note 1.B.
Formally, we consider a single-band model When the real and imaginary parts of which are functionally independent, the Hamiltonian will have a non-zero spectral area. For a given eigenvalue E 0 of the Bloch Hamiltonian, by solving H 0 (k) = Re E 0 and Γ(k) = Im E 0 , one can obtain a finite set of preimage of E 0 , i.e., K(E 0 ) = {k 1 , ..., k m }, which includes all Bloch waves having energy E 0 . Now suppose that one of the Bloch waves k i ∈ K(E 0 ) is incident on the boundary, depending on the normal direction of the boundary, the corresponding momentum of the reflected wave can be arbitrary. However, the number of elements of K(E 0 ) is finite, and as such cannot support so many reflection channels. This failure of reflection mechanism at a generic boundary means the failure in forming an open boundary eigenstate from Bloch waves, which implies the emergence of skin effect under a generic open-boundary geometry. However, the spectrum collapses into an arc (zero spectral area) if the real and imaginary parts of the Hamiltonian are functionally dependent, and the number of the corresponding solutions of H(k) = E 0 is infinite. It means that there are infinite reflection channels to satisfy the open boundary of any shape, and an open boundary eigenstate can be formed from superimposing all Bloch-wave channels. Concretely, we choose two examples to demonstrate the above arguments. The first example is of which the spectral area is nonzero shown in Fig. 4(a). For a given eigenvalue E 0 = 1 + i, by solving 2 cos k x = 1 and 2 cos k y = 1, we can obtain a finite set of pre-images of E 0 , that is, K(E 0 ) = {k 1 , k 2 , k 3 , k 4 } [red points in Fig. 4 The finite solutions of H(k) = E 0 cannot support so many reflection channels, that is to say, cannot form an openboundary eigenstate on a generic geometry by superimposing these Bloch waves specified by k i=1,2,3,4 . Therefore, the Hamiltonian Eq. (28) has skin effect under open-boundary geometry of a generic shape. The second example reads the periodic-boundary spectrum of which is an arc [the gray line in Fig. 4(c)]. The set of pre-images of E 0 = 1 + i has infinite elements [the red lines in Fig. 4(d)], which means that there are infinite ways of superimposing these Bloch waves to satisfy the open boundary condition of any shape. Therefore, the Hamiltonian Eq. (29) has no skin effect under any open-boundary geometry.

SUPPLEMENTARY NOTE 3: THE CURRENT FUNCTIONAL AND SKIN EFFECT
We divided this section into three parts. First, we give a general mathematical definition of the current functional (including zero and nonzero current functional). Second, we can completely classify the skin effect in two and higher dimensions into two types, that is, non-reciprocal skin effect (NRSE) and generalized reciprocal skin effect (GRSE), according to current functional. Finally, we discuss the restriction of symmetry on the current functional and the compatibility of skin effect with symmetry.

A. The definition of the current functional
In d dimensions, generally, the current functional is defined as under the periodic-boundary condition, where µ represents the energy band index and ∇ α is the directional derivative along certain direction α = d i=1 α iêki in d-dimensional momentum space (Hereê ki represents the i-th basis in momentum space).
Here, n(E, E * ) represents a distribution function when the system is in a steady state and only depends explicitly on the complex energy of the system state. Therefore, the current functional J α [n] is defined as the function of the distribution function n(E, E * ), and different input n(E, E * ) gives different output J α [n] with fixed α. In Hermitian case, the current functional becomes is the group velocity along direction α in the µ-th energy band. In non-Hermitian systems, the energy E is generally a complex number, therefore, the directional derivative of E µ (k) along α in Eq.(30) becomes which represents the generalized (complex) velocity in d dimensions.
Specially, in one-dimensional non-Hermitian system, the current functional reduces to where µ represents the band index. The generalized velocity for µ-th band Eq.(31) becomes which corresponds to the tangent vector of the µ-th energy band on the complex plane. Based on the above definition of current functional, we further define the "nonzero current functional" as, which is simply labeled as J = 0. As the complementary set, the "zero current functional" (labeled by J = 0) is defined as which means the current functional is zero regardless of the choice of α and n(E, E * ). For example, a Hermitian system always has zero current functional [2].

B. Current functional and the classification of skin effect
The nonzero current functional and zero current functional together constitute a complete set mathematically, and they are mutually exclusive. Based on this, we can completely classify the skin effect according to the zero and nonzero current functional.

One-dimensional skin effect
In one dimension, the current functional reduces to Eq.(32). In Ref. [2], the authors claim that for a one-dimensional system without any symmetry, if it has nonzero current functional (J = 0) in Eq.(34), then the system has skin effect (which is called Z skin effect), and vice versa. Another type of 1D skin effect -Z 2 skin effect has been reported in Ref. [3]. In this case, the system need to respect the spinful anomalous time-reversal symmetry, namely, In fact, we can prove that for the Z 2 skin effect, the system always has zero current functional (J = 0) in Eq.(35). A simple proof is present as follows.
This anomalous time-reversal symmetry requires the energy bands always come in pair and satisfy E ↑ (k) = E ↓ (−k). Correspondingly, for this pair of energy bands, the current functional satisfies J ↑ [n] = −J ↓ [n]. Here n(E, E * ) is invariant under the anomalous time-reversal symmetry, because n only depends on the energy and the complex energy is invariant under this symmetry. Therefore, the current functional for all energy bands Eq.(32) must be sum up to zero, regardless of the distribution function n(E, E * ).
To sum up, the 1D skin effect can be completely classified into Z skin effect with nonzero current functional (J = 0) and Z 2 skin effect with zero current functional (J = 0). (Note that this classification of skin effect according to the current functional is different from the classification of intrinsic point-gap topology for symmetry class [3].)

Two-and higher-dimensional skin effect
The theorem tells us that in two and higher dimensions the system has the universal skin effect, if and only if the spectral area is nonzero. According to the current functional, the universal skin effect can be further classified into non-reciprocal skin effect (J = 0) and generalized reciprocal skin effect (J = 0), as shown in Fig. 3 in the main text. By definition in Eq.(34) and Eq.(35), these two types of skin effect are complete and mutually exclusive.
We call the skin effect with nonzero current functional as non-reciprocal skin effect (NRSE), because this type of skin effect incompatible with the inversion or anomalous time-reversal symmetry [4,5]. Equivalently, the inversion and anomalous time-reversal symmetry always requires zero current functional (the rigorous proof is provided in the next subsection). Therefore, the NRSE is similar to the Z skin effect in one dimension, because the Z skin effect in 1D is forbidden by the inversion or anomalous time-reversal symmetry.
We name the skin effect with zero current functional as generalized reciprocal skin effect (GRSE) based on two reasons. The first reason is that this type of skin effect is compatible with the inversion or anomalous time-reversal symmetry. Note that an analog of GRSE in one dimension is the Z 2 skin effect [3], which can appear when the bulk Hamiltonian respects the spinful anomalous time-reversal symmetry (T 2 = −1). The second reason is stated as follows.
A reciprocal system requires the Hamiltonian to satisfy H t = H (H t (k) = H(−k) in momentum space) [6], which results in zero current functional. But generally, the zero current functional means more. For example, a Hamiltonian satisfying H t (k) = H(k + k θ ) (k θ = 0 and k θ = −2k) also gives rise to zero current functional. A simple model Hamiltonian reads H(k) = 2 cos k x + i sin k y , which satisfies H t (k x , k y ) = H(−k x , π − k y ) and has zero current functional. It can be seen that the zero current functional includes but is not limited to the reciprocal skin effect [6]. Therefore, we term the type of skin effect with zero current functional as the generalized reciprocal skin effect.

C. Current functional and symmetry
In this subsection, we first investigate the restriction of all point-group symmetries and time-reversal symmetry on the current functional. Meanwhile, we discuss the compatibility of all types of skin effect with these symmetries. Finally, we discuss the relationship between the universal skin effect and symmetry.

Current functional under point-group symmetry
In this part, we discuss the restriction of all point groups in three dimensions on the current functional. Here, if a point-group symmetry requires zero current functional, then the non-reciprocal skin effect (J = 0) can not appear under this symmetry. That is to say, the non-reciprocal skin effect is incompatible with this point-group symmetry. Obviously, the generalized reciprocal skin effect (J = 0) is compatible with all point-group symmetries.
In what follows, we demonstrate that the non-reciprocal skin effect is only compatible with point groups C m and C 2,3,4,6,2v,3v,4v,6v .
As the above definition, the current functional for µ-th band can be expressed as where ∇ α = d j=1 α j ∂ kj and d is the dimension of the system. Note that the distribution function n(E, E * ) only depends on the energy of the system state. Therefore, n(E, E * ) is invariant under the point-group operation due to the complex energy being invariant.
Inversion: Consider a system that only has inversion symmetry I, and each band satisfies . With the inversion symmetry, the current functional for µ-th band becomes It means that if the Hamiltonian has only inversion symmetry, the current functional for each band must be zero regardless of the choice of n(E, E * ). Equivalently, the non-reciprocal skin effect must vanish in the system with inversion symmetry. Therefore, the non-reciprocal skin effect is incompatible with the point groups including inversion symmetry, that is, C i,3i,2h,4h,6h , D 3d,2h,4h,6h , T h and O h . Rotation: Consider a system that is invariant under a point group including rotation operator R, then for mth energy band E µ (k) = E µ (R k). Under this rotation symmetry, the current functional becomes (the following derivation temporarily ignores the band index µ for simplicity) where det [J k,q ] in the first row is the determinant of the Jacobian J k,q that measures the change of differential volume element under different representations and the sign of det [J k,q ] is positive because the rotational operator preserves orientation. One can always choose an appropriate basis transformation such that det [J k,q ] = 1. In addition, since the Brillouin zone has the same symmetry group as the Hamiltonian and rotational operator R does not change the orientation, the integral region BZ is invariant under the point group.
The last equation of Eq. (38) requires α i = d j=1 α j R t ji = d j=1 R ij α j , which can be represented as which means that the direction of α is parallel to the rotational axis of R. Meanwhile, the component of α perpendicular to the rotational axis must be zero. For example, assume that R is a rotation operator that rotates θ along k z axis (θ = 0). The matrix representation of R in momentum space is, then Eq. (39) can be represented as the following matrix form, which requires α x = α y = 0. Therefore, only J z [n] component is allowed under the rotational symmetry along the z axis.
We conclude that if a point group contains two or more rotations with non-parallel rotational axes, the current functional for each band must be zero. If the point group contains only one rotation, a nonzero current functional along the rotational axis is allowed, thus non-reciprocal skin effect is compatible with the point groups including only one rotation, namely C 2,3,4,6 .
Mirror: Similarly, a mirror symmetry requires which means the nonzero current functional in the mirror-invariant plane is allowed. Therefore, the non-reciprocal skin effect is compatible with the point groups having only one mirror symmetry, that is, C m . Consider the point groups with both of rotation symmetry and mirror symmetry. In this case, the α need to satisfy Therefore, the current functional J α [n] can be nonzero when the rotational axis lies on the mirror plane. The point groups satisfying the above condition include C 2v,3v,4v,6v . Therefore, the non-reciprocal skin effect is compatible with these point groups C 2v,3v,4v,6v . So far, we discussed all three-dimensional point groups that allow the current functional to be nonzero. The same procedure can be easily done for the point groups in two dimensions. Equivalently, we can conclude that the non-reciprocal skin effect can appear under (or is compatible with) the following point groups,

Current functional with time-reversal symmetry
Before, we examine the restriction of all unitary point groups on the current functional and all types of skin effect. In this part, we simply discuss the current functional and skin effect in the system with time-reversal symmetry.
Here, the time-reversal symmetry refers to the collection of the two different time-reversal symmetries in Ref. [4], namely, the conventional (complex-conjugate type) one T and the anomalous (transpose type) oneT . Note that in one dimension, the non-Hermitian skin effect is compatible with T but incompatible with the spinless anomalous time-reversal symmetry (T 2 = +1) [5]. Therefore, onlyT need to be considered for the higher-dimensional skin effect.
A Hamiltonian respecting the anomalous time-reversal symmetry (aTRS) satisfies The spinless aTRS requires each energy band to satisfy E µ (k) = E µ (−k). As a result, the current functional for each band follows However, under the spinful aTRS, the energy bands come in pair and have the restriction, E µ (k) = Eμ(−k), which results in In this case, the current functional summing over all bands must be zero. We conclude that the non-reciprocal skin effect is incompatible with the aTRS, while the generalized reciprocal skin effect can appear under this symmetry.

The universal skin effect and symmetries
In the previous parts, we analyze the restrictions of all point groups and time-reversal symmetry. If a symmetry requires zero current functional, then the non-reciprocal skin effect is incompatible with this symmetry, which means that the non-reciprocal skin effect must vanish under this symmetry. Instead, as another type of universal skin effect, the generalized reciprocal skin effect can appear under these symmetries. It implies that the universal skin effect can appear under the above discussed symmetries, including all point groups and time-reversal symmetry.
In fact, according to our theorem, the universal skin effect disappears if and only if the spectral area is zero. Therefore, in order to discuss the role of symmetry on the universal skin effect, we need only to consider whether under the symmetry the spectral area is zero or not. However, there are no such internal (or point-group) symmetries to restrict the spectral area to be zero. Assume that the system has pseudo Hermiticity ηH † (k)η −1 = H(k) [4]. The eigenvalues of the Hamiltonian must be real or come in complex-conjugate pairs. It is possible for the system to have entire real Bloch spectrum, under which the spectral area is zero and the universal skin effect disappears. A particular example is η = 1, which reduces to the Hermitian case. On the other hand, once not all the Bloch spectrum are real but come in complex-conjugate pairs to cover a finite area on the complex plane, the spectral area becomes nonzero, and consequently, the universal skin effect appears.

SUPPLEMENTARY NOTE 4: CORNER-SKIN EFFECT AND GEOMETRY-DEPENDENT-SKIN EFFECT
In this section, we first define the corner-skin effect (CSE) and geometry-dependent-skin effect (GDSE), and discuss the localization of open-boundary eigenstates in CSE and GDSE. In addition, we show some features of the spectrum in GDSE example in the main text, and numerically verify that GDSE obeys the volume law, which differentiates from normal boundary states or higher-order skin modes [7,8].

A. The corner-skin effect
In this subsection, we define the CSE, and calculate the generalized Brillouin zone of the CSE example (Eq.(2) in the main text) to demonstrate the localization of the open-boundary eigenstates.

Definition
We define the CSE as a type of the non-reciprocal skin effect (J = 0) that exhibits the particular phenomenon that almost all eigenstates are localized at corners of the open-boundary geometry.
Therefore, the CSE inherits the features of the non-reciprocal skin effect, including nonzero current functional in Eq.(34) and incompatibility with all point groups except for {C m , C 2 , C 3 , C 4 , C 6 , C 2v , C 3v , C 4v , C 6v }.

The localization of eigenstates in CSE
The precise position where the eigenstates are localized depends on the generalized Brillouin zone (GBZ). As will be explained in the following contents, in the model shown in Fig. 3 Fig. 3(c). We first write down the Hamiltonian of the model Notice that under the open boundary of square geometry , the system can be solved by the separation of variables method. By replacing e ikx and e iky with complex variables β x and β y , respectively, the Hamiltonian can be rewritten as a polynomial where For the x-and y-directions, the corresponding 1D GBZs are determined by the following characteristic equations respectively. Using the approach developed in previous Refs. [9][10][11], the results are shown in Fig. 5 (a) and (b), where T x and T y denote the GBZs in the x-and y-directions, respectively. Having obtained T x and T y , the corresponding asymptotic energy spectrum of the 2D system in the thermodynamic limit becomes which is shown as the light blue region in Fig. 5 (c). Now the precise localization position of the eigenstates can be analyzed the corresponding 2D GBZ, which is T x ×T y in this model. Since the radius of T y is less than 1, all eigenstates are localized in the negative y-direction. In addition, the radius of T x greater than 1 indicates that all eigenstates concentrate in the positive x-direction. Therefore, the open-boundary eigenstates are localized in the lower right corner of the square geometry, as shown in Fig. 5(d). we define the quantity

Volume Law in CSE
where ψ n (x) is the n-th normalized right eigenstate and N is the number of these eigenstates [N = 60 * 60 = 3600 in Fig. 6(a)]. We immediately know that the sum of W (x) over the entire lattice sites must be 1, or symbolically, where "Lat" symbolically refers to the set of all lattice sites. For each given Hamiltonian, we can plot the spatial distribution of W (x) and mark the value of W (x) in different colors. This map is reflected in the colorbar, as shown in Fig. 6(a). If the system has no skin effect, the value of W (x) at each site is basically uniform and approximately equal to 1/N . However, in the system with CSE [Eq.(2) in the main text], W (x) is inhomogeneous in the square lattice, and, as shown in the colorbar of Fig. 6(a), the maximum of W (x) is about 0.0912. Here, we sum W (x) over the right lower corner region R C and obtain Note that the sum of W (x) over the right lower corner region is close to 1, which reveals that almost all eigenstates in CSE concentrate in the corner region R C . In addition, we take a scaling analysis for the CSE model. First, we judge an eigenstate ψ i (x) as a corner mode when it satisfies the following condition, which means that we can find this mode with 90% probability in the corner region R C . Here, the first point is that R C always refers to the 9 * 9 sites at the right lower corner of the square lattice regardless of the lattice size. The second point is that 90% here is not a rigid criterion for judging the corner skin mode, while what truly matters is the volume law. Second, we count the number of corner skin modes as increasing the lattice size from L = 18 to L = 70 (V = L 2 ). The data has been plotted in Fig. 6(b), where the red dots are the raw data and the blue line is the fitting curve with which tells us the corner skin modes follow the volume law.

B. The geometry-dependent-skin effect
In this subsection, we define the GDSE and give a criterion to determine where the open-boundary eigenstates are localized in the system with GDSE. Then we show some features in the spectra of the GDSE example (Eq.(3) in the main text). Finally, we take a numerical verification of the volume law.

Definition
Similar to the definition of CSE, the GDSE is one type of generalized reciprocal skin effect (J = 0) showing the unique phenomenon that there is at least one fully open boundary geometry under which the skin effect does not appear. Equivalently, the system exhibits GDSE when the system satisfies the following two points at the same time, (i) the spectral area is nonzero; (ii) there is at least one fully open boundary geometry under which the skin effect does not appear.
Therefore, GDSE inherits the features of the generalized reciprocal skin effect, including the zero current functional in Eq.(35) and compatibility with all point-group symmetries.
According to the definition of GDSE, the skin effect will disappear in certain open-boundary geometries. In what follows, we take an example with mirror symmetry to introduce an efficient method to find the geometry not showing the skin effect.
If the bulk Hamiltonian H(k) has one mirror symmetry, e.g., the mirror-x symmetry M x H(k x , k y )M −1 x = H(−k x , k y ), the (vertical) boundaries parallel to the mirror line does not exhibit the skin effect due to the spectral winding number being zero, regardless of the reference energy E b . If the bulk Hamiltonian has another mirror symmetry, for example, the mirror-y symmetry M y H(k x , k y )M −1 y = H(k x , −k y ), the (horizontal) boundaries parallel to the mirror-y line does not show the skin effect for the same reason. Therefore, we conclude that if the bulk Hamiltonian has two or more mirror symmetries, the skin effect does not appear on the open-boundary geometry with each boundary parallel to one of these mirror lines (one example is Fig. 3(g) in the main text).

The localization of eigenstates in GDSE
In this part, we provide a method to determine which edge (or surface) the eigenstates are localized on for a given open-boundary geometry.  Now we define the winding number. In two dimensions, for each edge of the open-boundary geometry, we first transform the momentum basis from k = (k x , k y ) T to q = (q , q y ) T by q = Sk; det S = 1. (60) Here, q is chosen as the momentum parallel to the edge, and S is a 2 * 2 matrix with unitary determinant, which is the element of group SL 2 (Z). Note that the choice of q y is not unique, it only needs to satisfy Eq.(60). Accordingly, the Hamiltonian transforms from H(k x , k y ) toH(q , q y ). Based on the above definition, we state that if the following condition is satisfied, then, no eigenstates are localized at this edge parallel to q . On the contrary, if the above condition is destroyed, some eigenstates must be localized at the corresponding edge. Note that this statement is strictly true when the other directions are periodic, and it is our conjecture in Fig. 2 For the "edge 1" and "edge 2" of the triangle geometry in Fig. 7(a)(c), the transformation matrix S equal to {{1, 0}, {0, 1}} and {{1, 1}, {0, 1}}, respectively. Obviously, for any fixed k x , the spectral winding number w E b (k x ) = 0, as shown in Fig. 7(b). The zero spectral winding number is guaranteed by the mirror symmetry in y direction of the bulk Hamiltonian, H(k x , k y ) = H(k x , −k y ). For "edge 2" of the geometry, the Hamiltonian can be transformed intõ For the fixed q , the spectrum ofH(q) forms an closed loop on the complex plane, as shown in Fig. 7(d). The above process can be taken for every edge of the open-boundary geometry. Therefore, the eigenstates are localized on the "edge 2" instead of "edge 1" of the triangle geometry, as shown in Fig. 3(h) of the main text.

The feature of spectrum in the GDSE example
We calculate the eigenstates and spectra of Hamiltonian Eq.(62) on the triangle geometry with different system size. From Fig. 8 we can observe that as the system size increases, the area covered by the open-boundary spectrum on the complex plane becomes gradually larger. It can be excepted that the area of the open-boundary spectrum under the thermodynamic limit tend to be the same as the spectral area (the light blue region in Fig. 8).
Note that even if the area of the energy spectrum under some open-boundary geometry seems to be the same as the spectral area, the system still has the skin effect due to the different density of states on the complex-energy plane, which is a unique feature in two-and higher-dimensional skin effect distinct from the one-dimensional case. We show the density of states plot to support the statement about GDSE in the main text. In Fig. 9, the openboundary spectra under periodic boundary, square open-boundary geometry and triangle open-boundary geometry are shown in Fig. 9(a)(b)(c), respectively. The corresponding density of states on the complex energy plane are plotted in Fig. 9(d-f), where the z-axis D(E) indicates the probability density of eigenvalues lying in the unit energy interval at E on the complex plane. The spectra under periodic boundary condition and square geometry have the same density of states on the complex energy plane as shown in Fig. 9(d)(e), where the spectral distribution at the corner is denser than the center. It is not the case on the triangle geometry. Even though the same region is covered by their spectra, the eigenvalues under the triangle geometry are more densely distributed at the center of the spectra as shown in Fig. 9(f).

Volume law in GDSE
We numerically verify that the geometry-dependent-skin effect of the model in the main text follows the volume law, that is, the increase in the number of skin modes is proportional to the increase in volume of the system, Here, we judge an eigenstate ψ(x) as a skin mode when it satisfies where B represents the boundary region that we specify. The Hamiltonian of the tight-binding model for GDSE is H(k x , k y ) = 2 cos k x + i cos k y . The spatial distribution W (x) of the eigenstates of this Hamiltonian under different open boundaries are plotted in Fig. 10(a)(b). One can observe that the skin modes disappear under the square geometry, but reappear under the parallelogram geometry, which is the feature of GDSE. For the parallelogram geometry, we specify the thickness of the boundary to be the width of three unit cells, and use the blue dashed lines to distinguish the boundary from the bulk in Fig. 10(b). If an eigenstate of the Hamiltonian under the parallelogram geometry satisfies Eq.(65), we count it as a skin mode. We count the number of skin modes N skin for different system size (V = L x × L y ), and fit the data (the red dots) into the blue line in Fig. 10(c). The results show that the GDSE satisfies the volume law, specifically, δN skin ≈ 0.494 δV in this model.

SUPPLEMENTARY NOTE 5: THE COROLLARY OF THE THEOREM
In this section, we will prove the corollary of our theorem, that is, all stable exceptional semimetals imply the universal skin effect. We review the topological charge of non-Hermitian band degeneracy in the first two subsection. Finally, we obtain the corollary of the theorem and prove the statement in the main text that the stable exceptional points for two-band model indicate nonzero spectral area A ± .

A. Non-Hermitian band degeneracy
Consider a general m-band non-Hermitian Bloch Hamiltonian (with periodic boundary condition), where Γ s are the generators of Lei algebra su(m) and h r s (k) and h i s (k) are real functions of k. When m = 2, 3, 4, Γ s refer to the Pauli, Gell-Mann, and Dirac matrices, respectively. The eigenvalues of H(k) can be obtained by solving the following characteristic polynomial where E i (k) is the ith eigenvalue of the non-Hermitian Hamiltonian H(k). At the degeneracy point k D , two bands must have the same energy, i.e.
for some i = j. In Ref. [12,13], the authors have shown that the above condition is equivalent to the vanishing of the discriminant of f E (k), i.e. where is the discriminant of f E (k). Although the discriminant is defined by the roots of f E (k) = 0, it can be computed directly from the determinant of the Sylvester matrix of f E (k) and ∂ E f E (k), which can be expressed by the coefficients of f E (k). Now we take a generic two-band model as an example to demonstrate this.
Example: Consider a generic two-band model where h µ (k) = h r µ (k) + ih i µ (k) are complex functions of k. The characteristic polynomial of the two-band model can be written as Computing the discriminant of polynomial (72) with respect to the energy E, we obtain the following condition for the existence of degeneracy points This condition can also be obtained from the energy spectrum, that is, the two bands The solutions of the above equation correspond to the non-Hermitian degeneracy points in 2D and lines in 3D.

B. Topological charge of non-Hermitian band degeneracy
In this subsection, we will review the topological charge of the non-Hermitian band degeneracies. Based on the discriminant of the characteristic polynomial, one can define the topological charge of the degeneracy point k D , i.e.
where Γ (k D ) is a loop encircling the degeneracy point k D . Since Disc E [H](k) is single valued, this invariant is quantized, which is called the discriminant number in Ref. [13]. Putting into ν(k D ), one can obtain Therefore, for a two-band system, which describes the winding of the complex energy between two bands. Now we show a concrete example of the winding number.
The topological charge ν(k D ) can be used to classify the non-Hermitian degeneracy points. However, the classification is not complete. As a comparison with H 1 (δk), we consider the following low energy Hamiltonian, Obviously, δk = 0 is a degeneracy point. One can further prove that its topological charge is +1, which is equal to the charge of δk = 0 in H 1 (δk). However, these two degeneracy points have different properties, for example, One can notice that H 1 (δk = 0) is non-diagonal. This type of non-Hermitian degeneracy points are called exceptional points (EPs). In Ref. [13], the authors have shown that only the exceptional points with ν(k D ) = ±1 are robust in 2D, while any other non-Hermitian band degeneracy points are unstable against non-Hermitian perturbations.

C. Stable exceptional points necessitate nonzero spectral area
Our corollary of the theorem is that all lattice Hamiltonians having stable exceptional points have universal skin effect. Here the stable exceptional point refers to the exceptional point of which the topological charge (discriminant number in Eq.(75)) is ±1. For the stable exceptional point, the spectral area A i must be nonzero. This fact is guaranteed by their topological properties. Therefore, the corollary can be obtained according to our theorem.
For simplicity, we here only consider a two-band model, and the Bloch Hamiltonian can be written as where Here we have omitted the h 0 (k)σ 0 term, which is irrelevant to the discussion of exceptional points. The eigenvalues of the Hamiltonian is The emergence of the exceptional point k EP requires that and there exists an invertible matrix P such that with a = 0 a general complex number. The topological charge of the exceptional point is which describes the winding number of the Bloch spectrum around E 0 = E ± (k EP ) = 0 and Γ(k EP ) represents a counterclockwise path enclosing k EP . Thus one can imagine that when ν (k EP ) = 0, the corresponding spectral area must be nonzero. Here we note that the above derivation does not require that the degeneracy point is an exceptional point. Actually, any degeneracy point with nonzero topological charge means the nonzero spectral area.
Now we consider the two-dimensional system with stable exceptional points (that is ν(k EP ) = ±1). The dispersion around the exceptional point k EP can be expanded [14] as where c x , c y are nonzero complex numbers and c x/y = c r x/y +ic i x/y with the superscript r/i indicates the real/imaginary part. The above equation implies that with c 1 = c y /c x = c r 1 + i c i 1 , where∆(q) = Re∆(q) + i Im∆(q) = (q x + c r 1 q y ) + i c i 1 q y . Putting this equation into the topological charge formula Eq.(88), equivalently, one can obtain that ν (k EP ) = sgn(det ∂ qx Re∆(q) ∂ qy Re∆(q) ∂ qx Im∆(q) ∂ qy Im∆(q) ) = sgn(det 1 c r where sgn(c i 1 ) represents the sign of the c i 1 . Note that in our statement in the "Corollary" section of the main text, c 0 = √ c x and c 1 = c y /c x . Therefore, a stable EP (ν(k EP ) = ±1) requires a nonzero imaginary part of c 1 , which further ensures a nonzero spectral area A ± .

SUPPLEMENTARY NOTE 6: THE PHOTONIC CRYSTAL MODEL
In this section, we supplement the spectra and eigenstates of the photonic crystal model [Eq.(4) in the main text]. Then we stack the 2D photonic crystal model along z axis and obtain a three-dimensional example that shows the GDSE. Finally, we explain the anomalous wave-packet dynamics shown in Fig 4(d) of the main text, and comment on the time evolution under the non-Hermitian setting.
A. The spectra and eigenstates of the photonic crystal model The Hamiltonian of the photonic crystal model reads where σ = (σ 0 , σ x , σ y , σ z ) is a vector of the Pauli matrices and d(k) is a vector with four components, that is, Here, we take the non-Hermitian parameter γ as 1/4. We numerically calculate the spectra in Fig. 11(a)(c) and spatial distribution of the wave function in Fig. 11 Fig. 11(b), and reappears under diamond geometry in Fig. 11(d), which is the signature of GDSE. Note that the spectrum under square geometry (the red dots in Fig. 11(a)) coincides with the spectrum under periodic boundary condition (the gray dots in Fig. 11(a)(c)), while the spectrum under diamond geometry (the red dots in Fig. 11(a)) does not.

B. An example of GDSE in three dimensions
In this subsection, we construct a three-dimensional model with exceptional lines, and demonstrate that it exhibits the GDSE.
Here, we construct a three-dimensional model by coupling the two-dimensional photonic crystal model Eq.(92) along the z axis. In this model, we can observe four exceptional lines crossing the Brillouin zone along k z direction, meanwhile, the system exhibits GDSE. The Hamiltonian of this model reads Here, we choose γ = 1, the only non-Hermitian parameter. The other parameters are the same as that in the model Eq.(92). In this 3D model, the coupling term along z direction is cos k z σ z . The Hamiltonian possesses four exceptional lines in 3D Brillouin zone, which are plotted as the blue lines in Fig. 12(a). When k z = ±π/2, the Hamiltonian Eq. (93) reduces to two-dimensional photonic crystal model Eq.(92) with four exceptional points. Next we demonstrate that the system exhibits GDSE. By definition, the sufficient and necessary condition for the existence of GDSE is that both of the following two points are satisfied: (i) the spectral area is nonzero; (ii) there is at least one geometry under which the skin effect disappear. The first point is satisfied due to the presence of the stable exceptional lines. The second point can be satisfied when we put the Hamiltonian on the cube geometry in Fig. 12(b), which is explained as follows.
First, the presence of the mirror symmetry, H(k x , k y , k z ) = H(k x , k y , −k z ), forbids the skin effect along z axis of the cube geometry. Second, for each fixed k z , cos k z can be absorbed into µ z , and the two-dimensional subsystem has no skin effect with square geometry. Therefore, with the cube geometry, the three-dimensional Hamiltonian does not exhibit skin effect, which is also verified by the numerical calculation. As shown in Fig. 12(b), the spatial distribution of eigenstates W (x) = 1 N n |ψ n (x)| 2 is uniform on the cube geometry, with system size L x = L y = L z = 12. The two points have been satisfied, therefore, the model has GDSE.
As a consequence of GDSE, the skin effect will reappear under other geometries. For example, if we cut the cube into two right-angle triangular prisms along z axis, the eigenstates will concentrate on the cross-section. Numerically, we calculate several two-dimensional subsystems with some generic k z , and observe that skin modes will reappear under the triangular geometry, as shown in Fig. 12(c). It implies the appearance of the skin modes on the cross-sections of two triangular prisms.

C. Wave-packet dynamics
In this subsection, we explain the transverse motion of wave packet shown in Fig. 4(d) in the main text, and comment on the time evolution of system state in classical wave system with gain and/or loss and dissipative open quantum system.

The transverse motion of wave packet induced by skin effect
We state that the skin effect causes some components of the wave packet to drift parallel to the edge in red color in Fig. 13(a), and suppresses the reflection components perpendicular to this edge. Therefore, after several bounces between the edges, the wave packet finally appears as a transverse motion shown in Fig. 4(d) in the main text.
In our example, the wave packet has Gaussian form centered at k c = (−2, −2) in momentum space, and evolves according to |ψ t = N (t)e −iHOBCt |ψ 0 .
Therefore, the wave packet will hit the edge in red color in Fig. 13(a) and scattering off it. In this process, the momentum component of the wave packet parallel to the edge is conserved. Hence, we transform the momentum basis from (k x , k y ) to (q x , q y ), as shown in Fig. 13(a), where k = Sq with S = {{−1, −1}, {1, −1}}. In q basis, the Gaussian wave packet is centered at q c = (0, 2). Accordingly, the Hamiltonian Eq.(92) can be rewritten in q basis, For a given q y (here we select q y = q 0 y = 3/2), the spectral winding number w E b (q 0 y ) is nonzero, which implies the skin effect on the edge 1 in Fig. 13(a). Next, we will show that this skin effect results in the transverse drift of some components of the wave packet.
• The transverse drift of the wave packet from the skin effect on edge 1 We focus on the component that has plane-wave form in q y direction and Gaussian form in q x direction. Note that the actual two-dimensional Gaussian wave packet is the coherent superposition of these components with a certain weight. Here we select the component with q 0 y = 3/2. In q x direction, the component is a Gaussian wave packet centered at q x = 0 as shown in Fig. 13(b)(c). In Hermitian case (γ = 0), the component will disperse with zero group velocity in q x direction. In Fig. 4(c) of the main text, therefore, the Gaussian wave packet always slowly disperses with time in q x direction.
However, in non-Hermitian case (γ = 1/4), the component will have a drift along the positive q x direction. As shown in Fig. 13(c), around q x = 0 each energy band has symmetric real part but asymmetric imaginary part, which corresponds to the presence of skin effect. In the range from q x = 0 to 1, the energy band in blue color has a larger imaginary part and positive group velocity (∂ Re E/∂q x > 0), while in the range from q x = −1 to 0, the red-color energy band has a larger imaginary part and positive group velocity. The components with larger imaginary part decay more slowly and therefore dominate the evolution of the wave packet as time goes on. Consequently, some components of the Gaussian wave packet have a transverse drift along the positive q x direction as shown in Fig. 4(d) of the main text, which is different from the Hermitian case.
Next, we will simply state that the suppression of the reflection wave components off the edge 2 ascribes to the skin effect on this edge.
• The suppression of the reflection wave components from the skin effect on edge 2 In the process of wave packet hitting the boundary and scattering off, q x is preserved. Therefore, we consider the wave-packet component that has plane-wave form in q x direction and Gaussian form in q y direction, labeled by q x -component. The actual two-dimensional Gaussian wave packet is composed of these q x -components with different weight. The skin effect on edge 2 means that for a given q x , the spectral winding number is nonzero, which suppresses the reflection wave of the q x -component.
To sum up, due to the skin effect on the edge 1 and edge 2, the Gaussian wave packet, going through several bounces between the two edges (edge 1 and the other parallel edge), finally appears as a transverse motion.

A comment on the time evolution in non-Hermitian system
In classical wave system, the evolution of a system state is governed by wave equation formally analogous to the Schrödinger equation [15]. In the cases with gain and/or loss, like the photonic crystal model in Ref. [16], the time evolution of a localized excitation (wave packet) is nonunitary, which can be captured by an effective non-Hermitian matrix [18]. Therefore, we believe that the phenomena in the simulation of wave-packet dynamics, e.g., the anomalous dynamical behavior shown in Fig. 4(d) of the main text, can be observed in a realistic classical wave system.
In the (driven or dissipative) open quantum system, a system state is described by the density matrixρ, whose time evolution is governed by the Lindblad master equation [19], HereĤ = xy H xyĉ † xĉy is the system Hamiltonian, andL x = y D xyĉy is the Lindblad dissipators describing quantum jumps due to coupling to the environment. In this setting, the dynamics of the density matrixρ(t) can be capture by the single-particle correlation C xy (t) = tr[ρ(t)ĉ † xĉy ]. After some tedious calculations, one obtains that the time evolution of the correlation follows [20,21] C(t) = G(t)C(0)G † (t); G(t) = e iH eff t , where H eff = (H − i D † D) t is the effective non-Hermitian matrix, the superscript "t" representing the transpose of the matrix. In the closed (Hermitian) quantum system, the system has no coupling with the external environment, and the density matrix evolves according to dρ(t)/dt = −i[Ĥ,ρ]. As a result, the dynamics of correlation satisfies where H t is transpose of the Hamiltonian matrix H. Comparing Eq.(97) with Eq.(98), we conclude that the time evolution of the correlation C(t) in dissipative open quantum system has the same form as the Hermitian case.