{'fraction_non_alphanumeric': 0.05829988346497475, 'fraction_numerical': 0.01777159135051146, 'mean_word_length': 3.890454329586829, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatial transient behavior in waveguides with lossy impedance boundary conditions

Attenuation in acoustic waveguides with lossy impedance boundary conditions are associated with non-Hermitian and non-normal operators. This subject has been extensively studied in fundamental and engineering research, and it has been traditionally assumed that the attenuation behavior of total sound power can be totally captured by considering the decay of each transverse mode individually. One of the classical tools in this context is the Cremer optimum concept that aims to maximize the attenuation of the least attenuated mode. However, a typical sound field may be a superposition of a large number of transverse modes which are nonorthogonal, and the individual mode attenuation may have little to do with the total sound power attenuation. By using singular value decomposition, we link the least attenuated total sound power to the maximum singular value of the non-normal propagator. The behaviors of the least attenuated total sound power depend only on the lossy boundary conditions and frequency, but are independent of sources. The sound may be almost non-decaying along the waveguide transition region for any lossy impedance boundary conditions although all modes attenuate exponentially. This spatial transient appears particularly strongly if the impedance is close to an exceptional point of the propagator, at which a pair of adjacent modes achieve maximum attenuation predicted by Cremer optimum concept. These results are confirmed using non-modal numerical calculations and a two-by-two toy model.


Introduction
Attenuation in acoustic waveguides with lossy impedance at the walls has been extensively investigated since the pioneering work of Morse [1] because it has important applications such as noise mitigation in aircraft engines or duct facilities [2]. As in any waveguides, the sound pressure can be decomposed into a set of transverse modes, and, to maximize the global attenuation, the classical approach is to focus only on attenuation of each transverse mode individually. Following this line of thought, one of the fundamental results is Cremer optimum concept [3] that aims to maximize the attenuation of the mode able to propagate on the longest distance (i.e. the least attenuated transverse mode). According to this Cremer criterion, it appears that the maximum attenuation occurs when the wavenumbers of the two lowest transverse modes coalesce. Cremer optimum criterion has been one of the most important tool for duct liner design [2,4,5,6,7,8,9,10], and it has been revisited recently in references [11] [12]. It is true that if a source is dominated by a single transverse mode, considering the attenuation of each mode individually may describe the Total Sound Power (TSP) attenuation. However, typical sources are not so simple; usually, there is a mix of transverse modes with energy distributed among them in a manner that is generally unknown, and several studies have resulted in the observation that Cremer criterion can be inefficient if an arbitrary source is considered [13,14,15]. The reason of this difficulty is that the transverse modes are non-orthogonal due to the lossy boundary condition at the walls, and a typical sound field is a superposition of a large number of transverse modes. Therefore, the sum of power attenuations of individual mode may have little to do with the TSP attenuation, meaning that the superposition of modes may dominate.
A matrix or operator whose eigenfunctions are non-orthogonal is said to be non-normal, and its eigenvalues may not totally capture its behavior. This can be illustrated by the failure of eigenvalue analysis in the hydrodynamic transition from a laminar to turbulent flow of Poiseuille and Couette flow [16]. The critical Reynolds numbers (where the transition should occur) predicted by eigenvalue analysis are very different from those measured in the laboratories. Trefethen and Embree [16], Trefethen et al. [17], Reddyet al. [18], and Butler and Farrel [19] found that the transitions are induced by the time transient amplification of certain small perturbations although all modes decay monotonically along flow direction. This transient behavior appeared to be due to the non-normality of the operator governing the time evolution of the perturbations. Recently, similar transient behavior was found in a non-normal optical systems with an unbalanced distribution of loss and gain [20].
In general, this non-normal transient behavior is closely connected to Exceptional Points (EPs) that have received much attention in recent years for non-Hermitian systems. The experimental demonstration of EPs has been observed by Dembowski et al. [21] in a microwave cavity with dissipation. The important properties of EPs have been uncovered by Heiss [22,23,24,25], Rotter [26], and Berry [27] et al.. EPs have been found in different systems, such as, laser-induced ionization states of atoms [28], electronic circuits [29], a chaotic optical microcavity [30], PTsymmetric waveguides [31], and waveguides with lossy impedance boundary conditions [35][36] [37] etc. A range of extraordinary phenomena related with EPs have been illustrated, such as, loss-induced transparency [32], enhancing the sensitivity in applications of microcavity sensors [33], parametric instability [34].
In this paper, we study the spatial transient behavior of the least attenuated TSP G max (x) as defined in Sec. 3, where x is the axis coordinate, of a uniform waveguide governed by Helmholtz equation and lossy impedance boundary conditions which is described by the matrix operator L as defined in Eq. (3). We use singular value decomposition to link G max (x) to the maximum singular value of the propagator e ixL . We find that G max (x) shows almost no decay in transient region although all modes decay exponentially. By comparing the least attenuated TSP G max (x) with the attenuation of the least attenuated mode which dominates the modal attenuations, we find that only when x → ∞, the least attenuated mode has the same decay rate with G max (x). In spatial transient region, they are different. This difference tend to be infinite when the admittance is close to the EPs at which a pair of modes achieve maximum attenuation predicted by Cremer optimum criterion. For each x, we find a least attenuated source that is generally close to a Gaussian form. When this source is imposed, the sound pressure field can take a pattern that avoids the lossy impedance wall and the TSP transports almost without decay in transient region, although all modes decay exponentially. This almost non-decayed TSP can be realized for any complex admittance, even in the vicinity of the EPs. A two-by-two toy model is also proposed to give in-depth understandings of these non-normality effects.

Model
We consider a semi-infinite, two-dimensional (2D) waveguide with one wall uniformly lined with locally reacting materials and the other wall acoustically rigid. Linear and lossless wave propagation is assumed. The configuration is depicted in Fig. 1. The sound pressure satisfies the Helmholtz equation with the boundary conditions where Y = ik/Z and Z refers to dimensionless impedance, it is a complex number. All lengths are divided by the waveguide width h. k = ωh/c, where ω and c refer to angular frequency and the sound velocity, respectively. Time dependence e −iωt is omitted in the followings. An acoustic source is located at x = 0. The sound wave propagates along x positive direction without reflection. We use a spectral collocation method based on Chebyshev polynomials to discretize sound pressure p(x, y) along y, p T = [p 1 , p 2 , · · · , p n , · · · , p N ] T = [p(y 1 ), p(y 2 ), · · · , p(y n ), · · · , p(y N )] T , where p is a column vector, " T " refers to the transpose, N refers to truncation number, the interpolation points y n = cos[(n − 1)π/(N − 1)], n = 1, · · · N. Note that p still depends on x which has been omitted. The differential operator ∂ 2 /∂y 2 is approximated by D 2 calculated by the MATLAB program chebdif. By taking only the right-going wave, Eq. (1) is written as where L = D 2 + k 2 I is a matrix operator, I refers to identity matrix. A detailed derivation of Eq. (3) is shown in Appendix A. It needs to stress that the boundary conditions (2) have been included in the matrix operator L. L is non-normal, i.e., LL † L † L when Y is complex, where " † " refers to the complex conjugate transpose. The modes of this operator are not orthogonal. L is also non-Hermitian, i.e., L L † . For a non-Hermitian operator, the eigenvalues and the corresponding eigenfunctions of a pair of adjacent modes may coalesce at some special points in the complex parameter plane, e.g., admittance plane in our case. These points are called Exceptional points (EPs) as we have pointed out in the Introduction. In our case, there are infinite number of EPs [35] in the complex value admittance plane. In this work, we consider only the first one. The solution to Eq. (3) is p(x) = e ixL p(0).
The TSP can then be expressed in terms of sound pressure distribution vector p where we have used the Clenshaw-Curtis quadrature to numerically calculate the integral 1 0 |p(x, y)| 2 dy, and for simplicity, we have omitted the independent variable x for vector function p(x). W is the weight matrix, it is diagonal and positive definite. It can be decomposed into W = F † F using a Cholesky decomposition, where F is still a diagonal matrix. In the last equation, we have assumed q = Fp and q refers to L 2 -norm of the complex vector q.

The least attenuated TSP
We define the attenuation of TSP G(x) as where the TSP at x is defined as 1 0 |p(x, y)| 2 dy = 1 0 p * (x, y)p(x, y) dy. G(x) is different from the usual definition of energy flux attenuation in waveguides. We will show in Sec. 7 that this difference does not change our conclusions.  (6), the attenuation of TSP G(x) can be calculated as where q 0 = Fp 0 represents sources. In the whole paper, q 0 is normalized so that q 0 2 = q † 0 q 0 = 1 0 |p(0, y)| 2 dy = 1, i.e., we assume the sound power of the source is normalised to 1. The variations of G(x) along x depend strongly on the sources [15].
The least attenuated TSP, at each x, can be achieved by "optimizing" G(x) over all permissible sources [38], where we have used q(x) = Fe ixL F −1 q 0 and invoked the definition of a vector-induced matrix norm in above derivation. σ n are the singular values of matrix Fe ixL F −1 .
For each given x, G max (x) in Eq. (8) gives the least attenuated TSP by taking care of the optimization over all permissible sources for a given waveguide with lossy impedance boundary conditions and frequency which is described by matrix operator L. G max (x) is totally decided by exp(ixL) (Note that F is an integration weight matrix.), therefore it is independent of sources. However, it is important to realize that, for each given x = x f , the least attenuated TSP G max (x f ) is achieved by using a different source which is given by the principal right singular vector of matrix operator Fe ixL F −1 [16] [38]. We call this source the least attenuated source, expressed as q least (x f ). When this source inputs, the usually calculated TSP attenuation G(x) as defined in Eq. (7) will be maximum at x = x f by taking care of the optimization over all permissible sources, and G max = G at x = x f as shown in Fig. 2 , and x f = 5. Therefore, curve G max (x) is not a TSP attenuation curve G(x) as usually illustrated in the literature of lined ducts, but an envelop of G(x). From the point of view of noise reduction in lined ducts, G max (x) gives the worst case at each x.
G max (x) can also be interpreted in term of modes. In Eq. (8), using eigenvalue decomposition of matrix, L = VΛV −1 , where V is a matrix with eigenvectors as columns, Λ is a diagonal matrix with eigenvalues on the main diagonal. Eq. (8) is then rewritten as exp(ixΛ) 2 represents the contributions of eigenvalues. If the boundary conditions are lossy, it represents mode attenuations and is dominated by the attenuation of the least attenuated mode, where Im(k 1 ) denotes the imaginary part of the eigenvalue of the least attenuated mode. Only when the boundary conditions are non-lossy (Y is a pure real number), is L normal, V a unitary matrix, V 2 V −1 2 = 1, and G max (x) totally decided by exp(ixΛ) 2 . Otherwise L is non-normal, V is not a unitary, V 2 V −1 2 may be very large when L is far from normality [38]. G max (x) may be highly different from M max (x) when L is far from normality [16]. M max (x) cannot capture the behavior of G max (x). Therefore, we define to describe the departure of the least attenuated TSP G max (x) from the attenuation of the least attenuated mode M max (x). The farther away from normality of the matrix operator L is, the larger is the departure of G max (x) from M max (x). Typical G max (x) and R(x) as functions of x are plotted in Fig. 3 when k = 2 for Y = 1 + 2i ((a) and (e)) and k = 10 for Y = 1 + i ((b) and (f)), Y = 1 + 4i ((c) and (g)), and Y = 1.65 + 2.06i which is close to EP ((d) and (h)), respectively. For comparison, we also plot the attenuations of the least attenuated modes M max (x) in Fig. 3 (a)-(d).
R(x) tend to be constant when x goes beyond certain values x tran , e.g., x tran ≈ 3 in (e), x tran ≈ 200 in (f), x tran ≈ 10 in (g), and x tran ≈ 1000 in (h). This means that G max (x) decay exponentially at the same rate as M max (x) as shown, e.g., in (a) and (c). This is the asymptotic region. Only in this region, M max (x) dominates the TSP and capture the decayed rates of G max (x). However, G max (x) and M max (x) still differ by a constant. This constant is a function of Y, and when Y is close to the EP as shown in Fig. 3 (h), this constant becomes large.
When x < x tran , G max (x) and R(x) are in spatial transient region. In Fig. 3 (a), k = 2 and Y = 1 + 2i, for this very low frequency, it is often assumed that higher modes attenuate rapidly, the least attenuated mode dominate the TSP. However, it is easy to see that the difference between G max (x) and M max (x) is not negligible. When Y is close to the EP (not shown in this figure), the difference becomes large.
When k increases, transient region can be roughly divided into two sub-regions: (1) almost non-decayed transient regions in which G max (x) is almost non-decayed although all modes attenuate exponentially. This region is approximately in 0 < x < 2 for k = 10 as shown in Fig. 3 (c) and (d). This almost non-decayed transient region can be realized for any Y, even when Y is close to the EPs, as shown in Fig. 3 (d) for Y = 1.65 + 2.06i, at which Cremer optimum criterion predicts that a pair of modes achieve maximum attenuation. (2) damped-oscillation transient region in which the behavior of G max (x) may be analogous to that of the displacement of a damped harmonic oscillator. In this region, G max (x) may be decayed oscillation as shown in Fig. 3 (b) when it is "underdamping" or decay non-exponentially as shown in Fig. 3 (c) and (d) when it is "overdamping". The analogy of "underdamping" and "overdamping" can be seen more clearly in Fig. 3 (f)-(h) for R(x). As will be shown by a two-by-two toy model in Sec. 6, the periods of the decayed oscillations δ are decided by the difference of real parts of the axial wavenumbers. When k = 10, the range of damped-oscillation transient regions start roughly at x = 2 as shown in Fig. 3 (b)-(d) and stop roughly at x = x tran as shown in Fig. 3 (f)-(h).

Non-normality indicator
In Fig. 3(e), (f), (g), and (h) we observe that R(x), which describes the departure of G max from M max , always start at 1 and then grows up to its maximum. We use this maximum departure to define an indicator of non-normality for the matrix operator L H is a function of Y.
We plot the contours of non-normality indicator ln(H) in the complex Y plane as shown in the left panel of Fig. 4 for k = 10. It is shown that H = 1 is minimum, it is located on axis of Im(Y) = 0, corresponding to L is normal. 1 ≤ H < 3 for most Y. H increases when Y → Y EP = 1.650611293539765 + 2.059981457179885i. The increasing behavior of H as a function of Im(Y) when Re(Y) = Re(Y EP 1 ) is shown in Fig. 4. We find numerically that H ∝ 1/|Y − Y EP |. We will return to this point in Sec. 6.

The least attenuated sources and the corresponding sound fields
As we have indicated in Sec. 3, for each given x = x f , the least attenuated TSP G max (x f ) is achieved by using a least attenuated source q least (x f ). When this source inputs, the usually calculated TSP attenuation G(x) as defined in Eq. (7) will be maximum at x = x f by taking care of the optimization over all permissible sources. q least (x f ) is the principal right singular vector of matrix operator Fe ixL F −1 [38]. We find that when x in the almost non-decayed transient region, the amplitude |q least | is approximately a half of Gaussian function along y (If the boundary condition at y = 0 is same as that of y = 1, |q least | will be approximately a Gaussian function.).
In Fig. 5 (a), we plot the TSP attenuation ln(G(x)) when the least attenuated sources q least (x f = 1.8) and q least (x f = 4) are incident, respectively, for Y = 1 + 4i and k = 10. x = 1.8 is still in the almost non-decayed transient region (x < 2) as shown in Fig. 3 (c) and (g). The TSP transports almost non-decayed in the range of x < 1.8. The corresponding sound pressure |p(x, y)| is shown in Fig. 5 (b). It is highly interesting that the sound field can avoid the lossy impedance boundary, like a half of a Gaussian beam, in the almost non-decayed transient region. However, x = 4 is out of almost non-decayed transient region as shown in Fig. 3 (c) and (g). Although G(x) = G max (x) at x = 4 is still maximum by considering all possible sources, q least (x f = 4) deviates from Gaussian function; G decays in the transient region; |p(x, y)| gradually increases in the vicinity of the lossy impedance boundary as shown in Fig. 5 (c).
We stress that the almost non-attenuated TSP in transient region can be realized for any Y, even in the vicinity of the EPs at which the the rates of exponential attenuations of a pair of modes which dominate the TSP are maximum.

Toy model
In this section, we construct a two-by-two toy model and analytically analyze the behavior of G max , R(x), and H. To solve Eq. (1) with the boundary conditions (2), we expand the sound pressure as proposed in Ref. [39] p(x, y) = N n=0 p n (x)ψ n (y), where ψ n (y) are the normalized eigenfunctions of rigid modes which satisfy the Helmholtz equation (1) and rigid boundary conditions ∂p/∂y = 0 at y = 0 and y = 1 respectively, N refers to the truncation number. Using Eq. (13), we project Eqs. (1) and (2) onto the corresponding rigid modes ψ n (y), we obtain dp dx = iLp, where In the derivations of above equations, we have chosen the truncation number N = 1. α 0 = 0 and α 1 = π refer to transverse eigenvalues of the first two rigid modes. p = [p 0 (x), p 1 (x)] T , where " T " refers to the transpose. a 1 = k 2 and a 2 = k 2 − α 2 1 = k 2 − π 2 are real, therefore a 1 − a 2 = π 2 . The detailed derivations of Eqs. (14)(15)(16) are shown in Appendix A. Equation (14) is similar to Eq. (3) except the dimensions of p and L.
The two eigenvalues k n = √ γ n of matrix L coalesce, where γ n refer to the eigenvalues of matrix L 1 , n = 1, 2, when the parameter Y = Y EP (or Y * EP ), where From Eq. (A-6), the corresponding eigenvectors also coalesce. Therefore Y EP and Y * EP are two EPs. By keeping in mind that Im(Y) must be non-negative for a passive system, we take only Y EP for our lossy waveguide problem.
Note that to assure the truncation of the expansion (13) to N = 1 represents a true passive sound pressure field in a waveguide with lossy impedance boundary conditions, we have chosen the imaginary parts of eigenvalues k 1,2 to be non-negative (convention exp(−iωt)), and proven dW(x)/dx < 0 as shown in Appendix B.
To analytically calculate the least attenuated TSP G max , we use p = Tq to transform matrix L to a triangle matrix L , where T is a unitary matrix by which the transformation conserves the TSP. Eq. (14) can be written as The detailed derivations of Eq. (18), the elements of matrix T, and C 1 are shown in Appendix C. The solution of Eq. (18) is where q 0 refers to a source at x = 0, the exponential of matrix L is Similar to Eq. (8), the least attenuated TSP is where σ n refer to the singular values of exponential matrix e ixL , In the expression of P t , we have used P t refers to Petermann factor of mode 1. When Y → Y EP , P t → ∞. Petermann factor, which describes the self-overlap of a mode [35], plays important roles in excess-noise for unstable lasers [40]. Recently, its important roles have also been found in open quantum systems [41] and waveguides with lossy boundary conditions [35]. Equations (21)(22)(23) show clearly that G max is not dependent on the sources whose effects are described by p n in Eq. (13). Note that, without loss of the generality, we have assumed mode 1 to be the least attenuated mode in Eqs. (21) and (22), i.e., Im(k 1 ) < Im(k 2 ). If, on the other hand, Im(k 2 ) < Im(k 1 ), M max (x) = e −2Im(k 2 )x , the least attenuated mode is mode 2. Note that in this section, we will analytically analyze the behavior of R and G max as functions of x and complex Y. We will use alternatively x and Y as independent variable.
We can find from Eqs. (22) and (23) that three factors contribute to the departure of G max from M max : P t ∝ 1/|Y − Y EP | where |Y − Y EP | describes the departure from the EP in the complex Y plane; cos(k R x) represents the interferences between the modes; and exp(−k I x) represents the differences of attenuation rates between the modes. cos(k R x) exp(−k I x) in Eq. (23) dominates the variation of R(x) in most cases. It may be analogous to the displacement of a damped harmonic oscillator. The Y plane could be roughly divided into two parts: in the region Im(Y) < Im(Y EP ), k R = O(k I ) or |k R | > |k I |, the transient behavior of R(x) is in damped oscillation as a underdamped oscillator as shown in Fig. 3 (f); in the region Im(Y) > Im(Y EP ), k I becomes larger as Im(Y) increases. R(x) decays as a overdamped oscillator as shown in Fig. 3 (g) and (h).
When x is small such that k R x 1 and k I x 1 as derived in Appendix D, where C 4 = 4C 2 C 3 Im(Y) 2 is a positive real number. The coefficients of x and x 2 are positive as shown in Appendix D. Therefore, R(x) always increases when x is small. For any given Y and Y Y EP , when x → ∞, e −k I x → 0, R(x) turns to be a constant, where C 3 is a positive real number. As Im(Y) → 0, this constant is close to 1. This constant will be very large when Y is close to Y EP because of P t ∝ 1/|Y − Y EP |. G max (x) then tends to [1 + C 3 Im(Y) 2 P t ]M max (x), decreases exponentially with the same rate as the least attenuated mode M max . For any given x, when Y → Y EP , we have k R → 0, k I → 0, k R x 1, k I x 1, and P t → ∞, the approximation of R(Y) can be obtained as same as that given in Appendix D for x small, where the approximation of (k 2 R + k 2 I )/(|γ 1 − γ 2 | 2 ) when Y → Y EP is shown in Appendix E. C 4 /|a 1 + a 2 + 3Y| ≈ 0.447 at EP when k = 10. Equation (27) shows clearly that when Y → Y EP , for large x, R may be very large.
7. The influence of using TSP to calculate G max All the above results are based on the attenuation of TSP defined in Eq. (6). However, we usually use energy flux to describe the attenuation of energy in waveguides. To verify whether this difference will have important influence on our conclusions, we re-plot the G max (x) in Fig. 6 (a) and (b) by using the parameters Y = 1 + 4i and k = 10 as shown in Fig. 3(c). Using the same parameters, we plot G(x) when q least (x f = 1) and q least (x f = 3) are incident respectively. In the same figure, we also plot the attenuation of the energy flux defined as by incidence of the same least attenuated sources. It is shown that E f (x) is very close to G(x). At x = 1 and x = 3, the differences between E f (x) and G(x) are of the order O(10 −13 ) and therefore could be negligible. We have implemented many comparisons for different q least (x f ) to justify the above results. We then confidently conclude that all our conclusions about G max obtained from G(x) defined in Eq. (6) are valid for energy flux defined in Eq. (28).

Conclusions
We have studied the spatial transient behavior of the least attenuated total sound power G max (x) described by the largest singular value of exp(ixL) in a waveguide with lossy admittance boundary conditions which can be described by a non-Hermitian and non-normal operator L. G max (x) is not curves of total sound power attenuation along wave propagation when sources input as usually illustrated in Acoustic literature. It is the envelop of them. G max (x) can be roughly divided into three regions: (1) almost non-decayed transient region in which G max is almost non-decayed although all modes attenuated exponentially. This almost non-decayed transient behavior can be realized for any complex value admittance, even when it is close to the EPs, the spectral singularities of non-Hermitian operators, at which a pair of adjacent modes achieve maximum attenuation rates predicted by Cremer optimum criterion. (2) damped-oscillation transient region in which G max (x) may be analogous to the displacement of a damped harmonic oscillator in underdamping or overdamping. (3) asymptotic region.
We have conducted detail studies in the differences between G max and M max which describes the attenuation of the least attenuated mode. We have shown that the departures of G max from M max depend on x and the lossy admittance boundary conditions. The maximum departure, defined as non-normality indicator, tends to be infinite when the admittance is close to the EPs. Only in the asymptotic region, G max (x) decays with the same rates as M max (x). However, G max (x) and M max (x) still differ by a constant. This constant tends to be infinite when the admittance is close to the EPs.
We have shown that for each given x, the corresponding G max (x) is achieved by a different source. When this source inputs at the start of the waveguide x = 0, the sound field can avoid the lossy admittance boundary, like a half of a Gaussian beam, in the almost non-decayed transient region.
Our analysis based on singular-value decomposition and a two-by-two toy model is generic and capture not only the asymptotic behavior governed only by eigenvalues, but the spatial transient behavior in which the eigenvalues may be misleading and non-normality plays an important role. We believe that our analysis approach and conclusions would be interesting for other lossy acoustic systems and might shed new light on other lossy physical systems.
Appendix A Derivation of Eqs. (3), (14) and (16) We consider the semi-infinite, two-dimensional (2D) waveguide as defined by the configuration depicted in Fig. 1 and the dimensionless Helmholtz equation (1) and boundary conditions (2). All length variables are divided by the waveguide width h. Time dependence e −iωt is omitted in the followings.
To derive Eq. (3), we have used a spectral collocation method based on Chebyshev polynomials to discretize sound pressure p(x, y) along y, p T = [p 1 , p 2 , · · · , p n , · · · , p N ] T = [p(y 1 ), p(y 2 ), · · · , p(y n ), · · · , p(y N )] T , where p is a column vector, " T " refers to the transpose, N refers to truncation number, the interpolation points y n = cos[(n − 1)π/(N − 1)], n = 1, · · · N. The differential operator ∂ 2 /∂y 2 is approximated by D 2 calculated by the MATLAB program chebdif. The Eq. (1) is then written as In an infinite and uniform (range independent) waveguide, p can be decomposed as right-and left-going waves, where L = D 2 + k 2 I is a matrix operator, I refers to identity matrix. We consider only right-going wave from the source at x = 0 and suppress the "+", obtain dp dx = iLp.
To derive Eqs. (14) and (16) where ψ = [ψ 0 , ψ 1 , · · · , ψ N ] T , " T " refers to the transpose, N the truncation number of vector ψ. By using the expansion (13), we obtain ∂ 2 p ∂x 2 + (k 2 I − α + Yψ(y = 1)ψ T (y = 1))p = 0, where the elements on the main diagonal of the diagonal matrix α are α 2 n , n = 0, 1, · · · , N. The above equation depends only on x, it is then rewritten as Similar to the derivation above, we obtain Eq. (14) where L = √ L 1 . If we truncate matrix L 1 to N = 1, Eq. (16) is obtained. To calculate the elements of L, we first calculate the eigenvalues and their corresponding normalized eigenvectors of matrix L 1 defined in Eq. (16), which are where the normalization constants Λ 1,2 are defined as where " † " refers to the complex conjugate transpose. We can show that i.e., eigenvectors are not mutually orthogonal. This non-orthogonality is a general property of a non-normal matrix. Besides, we can find the basis that is biorthogonal to X 1 and X 2 . In our case, since the matrix L 1 is complex symmetric (L 1 T = L 1 ), this biorthognal basis is simply X * 1 and X * 2 such that X m T X n = 0, m n. By using the bi-orthogonal eigenfunctions X 1,2 of matrix L 1 , L can be expressed as where k 1 = √ γ 1 and k 2 = √ γ 2 are the eigenvalues of matrix L or axial wavenumbers of wave p(x). L 11 , L 12 , L 21 , and L 22 are then obtained as follows Appendix B Validity of the two-by-two toy model as a true passive system To assure the truncation of the expansion (13) to N = 1 represents a true passive sound pressure field in a waveguide with lossy impedance boundary conditions, we must choose the imaginary parts of eigenvalues k 1 and k 2 to be non-negative, and the axial variation of TSP, dW(x)/dx < 0, where where L 2 = i(L − L † ). In the derivation of Eq. (B-1), we have used the expansion Eq. (13) and the orthogonal properties of ψ n . It is easy to verify that L 2 is Hermitian so that its eigenvalues γ ± = −Im(k 1 + k 2 ) ± [Im(k 1 − k 2 )] 2 + |C 1 | 2 are real. To assure dW(x)/dx < 0, L 2 has to be negative semidefinite, i.e., γ ± are nonpositive. This is equivalent to the condition that 4Im(k 1 )Im(k 2 ) > |C 1 | 2 . (B-2) We have numerically verified this condition over complex Y plane. Note that Wiersig [33] has derived a similar condition for a two-by-two non-Hermitian passive system.
The coefficients of x and x 2 in Eq. (25) are both non-negative. Since C 4 ≥ 0, we can easily verify the coefficient of x is always non-negative. For the coefficient of x 2 , since To verify the right hand side of above equation is greater than 0, we need only to show 2k 2 I + C 4 A > 2k I k 2 I + C 4 A, by just squaring both sides of this equation because each term is non-negative.