On the Properties of Phononic Eigenvalue Problems

TL;DR

This paper investigates the fundamental properties of phononic eigenvalue problems, providing a unified theoretical framework, extending the Plane Wave Expansion method, and supporting findings with numerical calculations.

Contribution

It offers a unified understanding of phononic eigenvalues and eigenvectors, and extends the PWE method to broader classes of eigenvalue problems.

Findings

Eigenvalues can be real or complex, with eigenvectors forming complete bases under certain conditions.

Phononic problems can be reformulated into linear eigenvalue problems for direct computation of physical quantities.

The extended PWE method applies to new eigenvalue forms, broadening its applicability.

Abstract

In this paper, we consider the operator properties of various phononic eigenvalue problems. We aim to answer some fundamental questions about the eigenvalues and eigenvectors of phononic operators. These include questions about the potential real and complex nature of the eigenvalues, whether the eigenvectors form a complete basis, what are the right orthogonality relationships, and how to create a complete basis when none may exist at the outset. In doing so we present a unified understanding of the properties of the phononic eigenvalues and eigenvectors which would emerge from any numerical method employed to compute such quantities. We show that the phononic problem can be cast into linear eigenvalue forms from which such quantities as frequencies, wavenumbers, and desired components of wavevectors can be directly ascertained without resorting to searches or quadratic eigenvalue problems and that the relevant properties of such quantities can be determined apriori through the analysis of the associated operators. We further show how the Plane Wave Expansion (PWE) method may be extended to solve each of these eigenvalue forms, thus extending the applicability of the PWE method to cases beyond those which have been considered till now. The theoretical discussions are supplemented with supporting numerical calculations. The techniques and results presented here directly apply to wave propagation in other periodic systems such as photonics.

Tables

Table 1. Table 1: The inner products of 6 eigenmodes at β 1 = π / 6 subscript 𝛽 1 𝜋 6 \beta_{1}=\pi/6 . The results below are normalized via | u ( l , r ) i ⁣ † ​ B ​ u r j / | u ( l , r ) i ⁣ † ​ B ​ u r i | | superscript subscript 𝑢 𝑙 𝑟 𝑖 † 𝐵 superscript subscript 𝑢 𝑟 𝑗 superscript subscript 𝑢 𝑙 𝑟 𝑖 † 𝐵 superscript subscript 𝑢 𝑟 𝑖 \Big{|}u_{(l,r)}^{i\dagger}Bu_{r}^{j}/\big{|}u_{(l,r)}^{i\dagger}Bu_{r}^{i}\big{|}\Big{|} . Note that only the absolute values of the inner products are presented.

	$u_r^1$	$u_r^2$	$u_r^3$	$u_r^4$	$u_r^5$	$u_r^6$	$u_l^1$	$u_l^2$	$u_l^3$	$u_l^4$	$u_l^5$	$u_l^6$
$u_r^1$	1.00	0.77	0.05	0.01	0.01	0.01	1.00	0.00	0.00	0.00	0.00	0.00
$u_r^2$	0.19	1.00	0.17	0.01	0.05	0.03	0.00	1.00	0.00	0.00	0.00	0.00
$u_r^3$	0.03	0.41	1.00	0.15	0.33	0.05	0.00	0.00	1.00	0.00	0.00	0.00
$u_r^4$	0.12	0.31	0.13	1.00	0.20	0.31	0.00	0.00	0.00	1.00	0.00	0.00
$u_r^5$	0.27	0.66	0.35	0.29	1.00	0.04	0.00	0.00	0.00	0.00	1.00	0.00
$u_r^6$	0.05	0.28	0.19	0.37	0.05	1.00	0.00	0.00	0.00	0.00	0.00	1.00

Equations

$$Av=\lambda Bv$$

$$⟨u,v⟩=\int u{v}^{\ast }dx$$

$$⟨Au,v⟩=⟨u,A^{\ast }v⟩$$

$$A{A}^{\ast }=A^{\ast }A;|Av|=|A^{\ast }v|$$

$$|(A-\lambda B)v|=|(A-\lambda B{)}^{\ast }v|=0,$$

$${\lambda }_1⟨Bv,w⟩=⟨Av,w⟩=⟨v,A^{\ast }w⟩=⟨v,{\lambda }_2^{\ast }{B}^{\ast }w⟩={\lambda }_2⟨Bv,w⟩$$

$$[A]\{v\}=\lambda [B]\{v\}$$

$$[C]\{v\}=\lambda \{v\};[C]=[B{]}^{-1}[A]$$

$$[H][A]\{v\}=\lambda \{v\}$$

$$\{v{\}}^{†}[A][H][A]\{v\}=\lambda \{v{\}}^{†}[A]\{v\}$$

$$\lambda =\frac{\{v{\}}^{†}[A][H][A]\{v\}}{\{v{\}}^{†}[A]\{v\}}$$

$$u_l^{†}A=\lambda u_l^{†}B$$

$${\lambda }_1{u}_l^{†}B{v}_r=u_l^{†}A{v}_r={\lambda }_2{u}_l^{†}B{v}_r,$$

$$(A-\lambda B{)}^m{u}_r^{(m)}=0;u_l^{(m)†}(A-\lambda B{)}^m=0,$$

$$(A-\lambda B)u_r^{(m+1)}=B{u}_r^{(m)}$$

$$u_l^{(m+1)†}(A-\lambda B)=u_l^{(m)†}B.$$

$$u_{lj}^{(n)†}B{u}_{rk}^{(m)}={\delta }_{jk}{\delta }_{m+n,m_i+1}.$$

$$\mathbf{u}=\bar{\mathbf{u}}\exp \left[i(\omega t-\beta \mathbf{n}\cdot \mathbf{x})\right]$$

$$\mathbf{\sigma }=\bar{\mathbf{\sigma }}\exp \left[i(\omega t-\beta \mathbf{n}\cdot \mathbf{x})\right]$$

$$(C_{ijkl}{\bar{u}}_{k,l}-i\beta n_l{C}_{ijkl}{\bar{u}}_k{)}_{,j}-i\beta n_j{C}_{ijkl}{\bar{u}}_{k,l}-{\beta }^2{n}_j{C}_{ijkl}{n}_l{\bar{u}}_k=\lambda \rho {\bar{u}}_i$$

$$A\bar{v}=\lambda B\bar{v}$$

$$⟨A\bar{v},\bar{w}⟩=\int \left[(C_{ijkl}{\bar{v}}_{k,l}-i\beta n_l{C}_{ijkl}{\bar{v}}_k{)}_{,j}-i\beta n_j{C}_{ijkl}{\bar{v}}_{k,l}-{\beta }^2{n}_j{C}_{ijkl}{n}_l{\bar{v}}_k\right]{\bar{w}}_i^{\ast }\mathrm{d}\Omega$$

$$⟨A\bar{v},\bar{w}⟩=\int {\bar{v}}_k\left[(C_{ijkl}{\bar{w}}_{i,j}^{\ast }+i\beta n_j{C}_{ijkl}{\bar{w}}_i^{\ast }{)}_{,l}+i\beta n_l{C}_{ijkl}{\bar{w}}_{i,j}^{\ast }-{\beta }^2{n}_j{C}_{ijkl}{n}_l{\bar{w}}_i^{\ast }\right]\mathrm{d}\Omega \equiv ⟨\bar{v},A^{\ast }\bar{w}⟩$$

$$A^{\ast }\bar{w}=(C_{ijkl}^{\ast }{\bar{w}}_{i,j}-i{\beta }^{\ast }{n}_j{C}_{ijkl}^{\ast }{\bar{w}}_i{)}_{,l}-i{\beta }^{\ast }{n}_l{C}_{ijkl}^{\ast }{\bar{w}}_{i,j}-{\beta }^{\ast 2}{n}_j{C}_{ijkl}^{\ast }{n}_l{\bar{w}}_i$$

$$|A\bar{v}{|}^2=⟨A\bar{v},A\bar{v}⟩$$

$$=\int \left[\nabla \cdot (\mathbf{C}\nabla \bar{\mathbf{v}}-i\beta \mathbf{Cn}\bar{\mathbf{v}})-i\beta \mathbf{nC}\nabla \bar{\mathbf{v}}-{\beta }^2\mathbf{nCn}\bar{\mathbf{v}}\right]\cdot \left[\nabla \cdot ({\mathbf{C}}^{\ast }\nabla {\bar{\mathbf{v}}}^{\ast }+i{\beta }^{\ast }{\mathbf{C}}^{\ast }\mathbf{n}{\bar{\mathbf{v}}}^{\ast })+i{\beta }^{\ast }\mathbf{n}{\mathbf{C}}^{\ast }\nabla {\bar{\mathbf{v}}}^{\ast }-{\beta }^{\ast 2}\mathbf{n}{\mathbf{C}}^{\ast }\mathbf{n}{\bar{\mathbf{v}}}^{\ast }\right]\mathrm{d}\Omega$$

$$=\int \left[\nabla \cdot \mathbf{C}\nabla \bar{\mathbf{v}}-{\beta }^2\mathbf{nCn}\bar{\mathbf{v}}-i\beta (\nabla \cdot \mathbf{Cn}\bar{\mathbf{v}}+\mathbf{nC}\nabla \bar{\mathbf{v}})\right]\cdot \left[\nabla \cdot {\mathbf{C}}^{\ast }\nabla {\bar{\mathbf{v}}}^{\ast }-{\beta }^{\ast 2}\mathbf{n}{\mathbf{C}}^{\ast }\mathbf{n}{\bar{\mathbf{v}}}^{\ast }+i{\beta }^{\ast }(\nabla \cdot {\mathbf{C}}^{\ast }\mathbf{n}{\bar{\mathbf{v}}}^{\ast }+\mathbf{n}{\mathbf{C}}^{\ast }\nabla {\bar{\mathbf{v}}}^{\ast })\right]\mathrm{d}\Omega$$

$$=\int \left[\mathbf{a}\cdot {\mathbf{a}}^{\ast }+|\beta {|}^2\mathbf{b}\cdot {\mathbf{b}}^{\ast }-i\beta {\mathbf{a}}^{\ast }\cdot \mathbf{b}+i{\beta }^{\ast }\mathbf{a}\cdot {\mathbf{b}}^{\ast }\right]\mathrm{d}\Omega$$

$$|A^{\ast }\bar{v}{|}^2=⟨A^{\ast }\bar{v},A^{\ast }\bar{v}⟩$$

$$=\int \left[\nabla \cdot (\nabla \bar{\mathbf{v}}{\mathbf{C}}^{\ast }-i{\beta }^{\ast }\mathbf{n}\bar{\mathbf{v}}{\mathbf{C}}^{\ast })-i{\beta }^{\ast }\nabla \bar{\mathbf{v}}{\mathbf{C}}^{\ast }\mathbf{n}-{\beta }^{\ast 2}\mathbf{n}\bar{\mathbf{v}}{\mathbf{C}}^{\ast }\mathbf{n}\right]\cdot \left[\nabla \cdot (\nabla {\bar{\mathbf{v}}}^{\ast }\mathbf{C}+i\beta \mathbf{n}{\bar{\mathbf{v}}}^{\ast }\mathbf{C})+i\beta \nabla {\bar{\mathbf{v}}}^{\ast }\mathbf{Cn}-{\beta }^2\mathbf{n}{\bar{\mathbf{v}}}^{\ast }\mathbf{Cn}\right]\mathrm{d}\Omega$$

$$=\int \left[\mathbf{c}\cdot {\mathbf{c}}^{\ast }+|\beta {|}^2\mathbf{d}\cdot {\mathbf{d}}^{\ast }-i{\beta }^{\ast }{\mathbf{c}}^{\ast }\cdot \mathbf{d}+i\beta \mathbf{c}\cdot {\mathbf{d}}^{\ast }\right]\mathrm{d}\Omega$$

$$A^{\ast }\bar{w}=(C_{klij}{\bar{w}}_{i,j}-i\beta n_j{C}_{klij}{\bar{w}}_i{)}_{,l}-i\beta n_l{C}_{klij}{\bar{w}}_{i,j}-{\beta }^2{n}_j{C}_{klij}{n}_l{\bar{w}}_i$$

$$⟨A\bar{v},\bar{v}⟩=\int \left[(C_{ijkl}{\bar{v}}_{k,l}-i\beta n_l{C}_{ijkl}{\bar{v}}_k{)}_{,j}-i\beta n_j{C}_{ijkl}{\bar{v}}_{k,l}-{\beta }^2{n}_j{C}_{ijkl}{n}_l{\bar{v}}_k\right]{\bar{v}}_i^{\ast }\mathrm{d}\Omega$$

$$=\int \left[-{\bar{v}}_{i,j}^{\ast }{C}_{ijkl}{\bar{v}}_{k,l}-{\beta }^2{n}_j{\bar{v}}_i^{\ast }{C}_{ijkl}{n}_l{\bar{v}}_k+i\beta {\bar{v}}_{i,j}^{\ast }{C}_{ijkl}{n}_l{\bar{v}}_k-i\beta n_j{\bar{v}}_i^{\ast }{C}_{ijkl}{\bar{v}}_{k,l}\right]\mathrm{d}\Omega$$

$$=-\int \left[{\bar{v}}_{i,j}^{\ast }{C}_{ijkl}{\bar{v}}_{k,l}+{\beta }^2{n}_j{\bar{v}}_i^{\ast }{C}_{ijkl}{n}_l{\bar{v}}_k-i\beta {\bar{v}}_{i,j}^{\ast }{C}_{ijkl}{n}_l{\bar{v}}_k+i\beta n_j{\bar{v}}_i^{\ast }{C}_{ijkl}{\bar{v}}_{k,l}\right]\mathrm{d}\Omega$$

$$=-\int \left[{\bar{v}}_{i,j}^{\ast }-i\beta {\bar{v}}_i^{\ast }{n}_j\right]C_{ijkl}\left[{\bar{v}}_{k,l}+i\beta {\bar{v}}_k{n}_l\right]\mathrm{d}\Omega$$

Full text

Spectral theorem, Phononics, Metamaterial, Wave Propagation, Modeshape Orthogonality, Scattering

I Introduction

There has been considerable recent research interest in the field of wave propagation in periodic structures under the fields of photonics, phononics, and even metamaterials. Much of the progress in these fields depends upon the determination of wave propagation characteristics in such periodic systemsNemat-Nasser (2015); Shmuel and Band (2016); Chen and Elbanna (2016); Mokhtari et al. (2019). Historically, numerical efforts in this direction have been driven towards the calculation of the so-called bandstructureDeymier (2013) of the periodic system which is a graphical representation of the frequency-wavevector pairs which satisfy a certain kind of dispersion relationship for the system. Traditionally, such dispersion relations have been calculated through what we would call the conventional form of the eigenvalue problem – determining acceptable frequencies given a wavevector – termed $\omega(\boldsymbol{\beta})$ systems. A host of numerical techniques have been devised to solve this particular form of the eigenvalue problem. This includes the Plane Wave Expansion (PWE) method Ho et al. (1990); Deymier (2013); Kushwaha et al. (1994), the multiple scattering methodMei et al. (2005); Kafesaki and Economou (1999), variational techniquesSrivastava and Nemat-Nasser (2014); Lu and Srivastava (2016, 2017), FEM White et al. (1989); Veres and Berer (2012); Hladky-Hennion and Decarpigny (1991); Hussein (2009a), and Finite Difference Tanaka et al. (2000); Chan et al. (1995) techniques among others.

Of late, there has been growing interest in the solution of the eigenvalue problem of periodic systems when the problem is not in a traditional form, often involving complex frequencies and wavenumbers. Most commonly, complex wavenumbers and frequencies are a direct outcome of including dissipation in the system, with a complex frequency representing a temporally dissipating wave whereas a complex wavenumber representing a spatially dissipating wave. However, imaginary and/or complex wavenumbers appear even in the absence of dissipation and serve the very important purpose of satisfying interface continuity conditions in scattering problems. The simplest of the non-standard eigenvalue cases is the solution of the $\omega(\boldsymbol{\beta})$ eigenvalue form in the presence of dissipationHussein (2009b); Hussein and Frazier (2010); Frazier and Hussein (2015). In this case, it turns out that the resulting frequencies for real assumed wavevectors are complex. A further complication which has been considered in literature is the determination of the wavenumber when frequency is given (termed $\beta(\omega,\mathbf{n})$ form where $\mathbf{n}$ is the direction of wave propagation). The corresponding eigenvalue problem is most naturally quadratic and, therefore, more difficult to solve than linear eigenvalue problems. In Ref. Frazier and Hussein (2016), an algorithm for 1D systems was developed that provided dispersion curves for damped free wave motion based on frequencies and wavenumbers that are permitted to be simultaneously complex. The algorithm was applied to a viscously damped mass-in-mass metamaterial exhibiting local resonance. In their study, two eigenvalue problems were solved: Frequency solutions from linear eigenvalue problem, and wavenumber solutions from quadratic eigenvalue problem. As the latter problem is quadratic, a search algorithm was presented to find the wavenumber solutions for a given frequency. The problem can alternatively be converted into a linear eigenvalue form by using a state space representationLaude et al. (2009); Moiseyenko and Laude (2011); Andreassen and Jensen (2013); Krushynska et al. (2016). The resulting mixed-form of the elastodynamics problem has been considered in detail in the Finite Element literature (Least-Squares FEMGunzburger and Bochev (2009)) but its appearance in the phononics/photonics area is rare. Computational techniques used to solve such a problem in the area of phononics/photonics are, therefore, limited as well to what is called the Extended PWE methodHsue et al. (2005); Laude et al. (2009).

A further complication, rarely studied till now, which could be considered is the determination of one component of the wavevector when the other components and the frequency are given - termed $\beta_{3}(\omega,\beta_{\alpha})$ problems. Such problems naturally emerge in cases where scattering at an interface are being studied Srivastava and Willis (2017); Hsue et al. (2005) This is due to the fact that Snell’s law ensures that the components of the wavevector tangential to the interface are preserved and, therefore, there is a natural requirement to determine the remaining component when the preserved component is specified. Currently, there appears to exist no study in phononics which could directly solve this problem.

In this paper, we ask some basic questions pertaining to the phononic eigenvalue problems and propose some solutions. Foremost, is our attempt to analyze the phononic eigenvalue problems in the three different forms mentioned above through the lens of linear operators. This exercise reveals to us the basic properties of the eigenvalues and eigenvectors which can be expected from numerical calculations without actually doing any calculations. In all the eigenvalue forms, we are interested in determining the appropriate orthogonality conditions and whether the eigenvector basis is complete. In addition to the theoretical considerations, we also present extensions to the PWE method and give representative solutions for the eigenvalue forms considered. This paper is organized as follows: in section II a brief introduction of the spectral theorem and properties of self-adjoint and non-self-adjoint operators are presented. We then investigate the properties of $\omega(\beta)$ problem for different cases of material properties and also different types of wave number (real and complex) in section III. In section IV, the general problem of $\beta(\omega)$ and the corresponding properties of its eigenvalues and eigenvectors are studied. One special case of this problem is the $\beta_{3}(\omega,\beta_{\alpha})$ problem which is formulated in section IV.

II Some results on the eigenvalue properties of linear operators

Here we are concerned with eigenvalue problems of the following form:

[TABLE]

where $A,B$ are linear operators which act on a dense domain in the Hilbert space $H$ of complex functions over which an inner product has been defined Schmüdgen (2012); Davies (1995):

[TABLE]

where $u,v\in H$. The adjoint operator to $A$, denoted by $A^{*}$, is defined by:

[TABLE]

Consider the case where $A,B$ are normal operators:

[TABLE]

with similar relations for $B$. If $\lambda$ is an eigenvalue of $A$ with associated eigenvector $v$ ($Av=\lambda Bv$) then

[TABLE]

showing that $\lambda^{*}$ is an eigenvalue of $A^{*}$ with the same eigenvector ($A^{*}v=\lambda^{*}B^{*}v$). If $\lambda_{1},\lambda_{2}$ are two distinct eigenvalues with associated eigenvectors $v,w$ then

[TABLE]

showing that the eigenvectors are orthogonal ($\langle Bv,w\rangle=0$).

Most of the standard eigenvalue problems in phononics and photonics and other areas of physics are governed by linear operators that are self-adjoint or Hermitian in characterKostenbauder et al. (1997). The normal modes of such a system, referring to the eigenfunctions of the Hermitian operator, then form a complete and orthonormal set for expanding the field variables in the system like displacement, stress, etc. Morse P.M. (1953). However, there are several important physical problems in phononics (like scattering of the wave at the metamaterial interfaces) and other fields, for which the governing operator is not Hermitian or self-adjoint and cannot easily be made so. There is no universal way for unbounded differential operators to figure out if their eigenfunctions form a complete basis or not. But at a very general level, if a differential operator is acting on a finite domain with periodic boundary condition, then it has a point spectrum and it is possible to expand the state of the system using the normal eigenmodes Davies (1995). If the operator is not normal, it still has propagation eigenmode. However, these modes are no longer orthogonal. In fact, they are orthogonal to the eigenmodes of the adjoint operator, a case which is called bi-orthogonalityKostenbauder et al. (1997). Below we present a brief explanation on the spectral theory for Hermitian operators, non-Hermitian operators, the problem of bi-orthogonality and also the generalized eigenvectors and formation of a complete basis.

Hermitian operators:

Self-adjoint operators ($A=A^{*}$) form a subset of normal operators and if $A,B$ are self-adjoint and $B$ is positive definite then the associated eigenvalues are real. In both cases the eigenvectors form a complete basis. The matrix representations of $A,B$ are $[A],[B]$ respectively. If $A$ is positive definite ($\langle Av,v\rangle>0$) then $[A]$ is also positive definite ($\{x\}^{\dagger}[A]\{x\}>0$) where $\dagger$ represents conjugate transpose. Similar consideration holds for negative-definiteness as well. Using a finite set of basis functions, the generalized eigenvalue problem can be represented in its finite dimensional form:

[TABLE]

Now if $[B]$ is positive definite then all its eigenvalues are positive. Since the eigenvalues of $[B]^{-1}$ are the inverse of the eigenvalues of $[B]$, its eigenvalues are also positive. Therefore, positive-definiteness of $[B]$ implies the positive-definiteness of $[B]^{-1}$. Now converting the generalized eigenvalue problem above into a standard eigenvalue problem:

[TABLE]

Denoting $[B]^{-1}=[H]$ we have $(\{v\}[A])^{\dagger}[H][A]\{v\}=\{v\}^{\dagger}[A][H][A]\{v\}>0$ since $[H]$ is positive definite and $[A]$ is Hermitian. Now we have the generalized eigenvalue problem

[TABLE]

showing that all eigenvalues $\lambda$ of the generalized eigenvalue problem are negative since the numerator is positive and the denominator is negative. Therefore, in a eigenvalue problem $Av=\lambda Bv$ where $A$ is negative definite and self-adjoint and $B$ is positive definite, all $\lambda$ are real and negative.

Non-Hermitian Operators:

Non-Hermitian systems have drawn significant recent attention due to their intriguing physics revolving $\mathcal{PT}-$symmetry Bender and Boettcher (1998) and exceptional points Heiss (2000); Heiss and Sannino (1990); Lu and Srivastava (2018). A recent review on non-Hermitian physics is given by Ref El-Ganainy et al. (2018). Questions about orthogonality condition and completeness of basis of non-Hermitian operators are of particular interest in wave scattering settings because they require an appropriate set of wave functions to expand wave fields in Srivastava (2016). In connection to the earlier results, any self-adjoint (Hermitian) operator can be shown to be normal. Vice-versa, if an operator is not normal then it can be shown to be non-self-adjoint (non-Hermitian). It is well known that the right eigenvectors of a non-Hermitian operator are not orthogonal and an orthogonalization process such as the Gram-Schmidt method does not guarantee to work unless the eigenvectors are at least linearly independent Arfken (1985). To discuss this further, we need to consider two cases, diagonalizable and non-diagonalizable operators. In the discussions below, the square brackets denoting matrices are dropped and the matrix nature of the equations is assumed implicit.

First, for a generalized eigenvalue problem such as Eq. (1) we consider its left eigenvector $u_{l}$ which satisfies

[TABLE]

Note that $u_{l}$ can be solved as the right eigenvector of $A^{*}$ and $B^{*}$ with associated eigenvalue $\lambda^{*}$. If $\lambda_{1}$ is an eigenvalue with left eigenvector $u_{l}^{\dagger}$ and $\lambda_{2}\neq\lambda_{1}$ is another eigenvalue with right eigenvector $u_{r}$ then we have:

[TABLE]

The above shows that for two distinct eigenvalues $\lambda_{1},\lambda_{2}$ of a non-Hermitian eigenvalue problem, the bi-orthogonality relation $u_{l}^{\dagger}Bv_{r}=0$ is satisfied. In operator notation, this bi-orthogonality condition can be written as $\langle B(\mathbf{v_{r}}),\mathbf{u_{l}}\rangle=0$. A numerical example which verifies this bi-orthogonality for a relevant phononic case will be given later in Section IV.3.

This bi-orthogonality condition for non-Hermitian operators has been widely accepted in the physics community El-Ganainy et al. (2018), however, it should be noted that it only holds under the assumption that the eigenvalues are distinct (operator is diagonalizable). If the non-Hermitian operator has repeated eigenvalues (such as in the case of exceptional points), then it can be shown to be non-diagonalizableRiechers and Crutchfield (2018) in which case, the bi-orthogonality condition will not apply and enough linearly independent eigenvectors cannot be found through conventional eigenvalue computation. However, in such a case, the operator can be transformed into its Jordan canonical form, and generalized left and right eigenvectors corresponding to each Jordan block can be defined which then satisfy a bi-orthogonality condition. The generalized eigenvectors thus calculated will form a linearly independent set.

For example, if a repeated eigenvalue $\lambda$ corresponds to a Jordan block, $\mathbf{J}_{i}$, of size $m_{i}$, then it can be associated with one pair of conventional left and right eigenvectors $u_{l}^{(1)}$ and $u_{r}^{(1)}$, and $m_{i}-1$ pairs of generalized left and right eigenvectors defined through the modified eigenvalue equation

[TABLE]

where $0\leq m\leq m_{i}-1$ with $u_{r}^{(0)}=0$ and $u_{l}^{(0)}=0$. The generalized eigenvectors are obtained by first solving the conventional eigenvectors and then following the process shown below

[TABLE]

Ref. Riechers and Crutchfield (2018) has presented an excellent discussion on the bi-orthogonality properties of generalized eigenvectors. Here, we only note the primary result of immediate relevance to us. If $u_{lj}^{(m)},u_{rj}^{(n)}$ ($u_{lk}^{(m)},u_{rk}^{(n)}$) are the left and right generalized eigenvectors corresponding to the repeated eigenvalue $\lambda_{j}$ ($\lambda_{k}$) for which there exists the Jordan block $\mathbf{J}_{j}$ ($\mathbf{J}_{k}$) of size $m_{j}$ ($m_{k}$) then the following general bi-orthogonality relationship holds:

[TABLE]

This completes the discussion of the orthogonality conditions for non-Hermitian operators. It must be noted that for phononic eigenvalue problems, the need to create generalized left and right eigenvectors will only be in the case where an exceptional points is involved in the eigen-spectrum. In the sections below, we consider the phononic eigenvalue problem under various operator forms, elucidating their relevant properties. We also present simultaneous numerical examples which explicitly show how the operator equations may be solved. These numerical solutions are based upon modifications of the well studied Plane Wave Expansion method, however, other numerical methods are also possible. The theoretical discussions are presented in the form of stiffness tensors at some places and compliance tensors at other places. The choice is mainly directed by the simplicity of resulting expressions and has no bearing on the conclusions reached. The PWE scheme is always implemented based upon the stiffness tensor in this paper.

III $\omega(\beta)$ solutions for a wave with a given wavevector

Consider Bloch waves propagating in a phononic crystal in direction $\mathbf{n}$. The displacement and stress fields due to the wave will have the general form:

[TABLE]

where $\bar{\mathbf{u}},\bar{\boldsymbol{\sigma}}$ are $\Omega-$periodic. $\omega$ and $\beta$ can both potentially assume real, imaginary, or complex values, however, only certain combinations are physically meaningful. First consider the usual form of the eigenvalue problem:

[TABLE]

where $\lambda=-\omega^{2}$. which is in the form

[TABLE]

when one identifies $\bar{v}\equiv\{\bar{\mathbf{u}}\}$. Now consider two $\Omega$ periodic fields $\bar{v},\bar{w}$. The relevant inner product for the operator $A$ is:

[TABLE]

Since $C_{ijkl},\bar{v},\bar{w}$ are all $\Omega$ periodic terms, the above can be transformed using Gauss theorem into:

[TABLE]

showing that the adjoint operator $A^{*}$ is:

[TABLE]

In what follows we will use $\mathbf{C}$ for the tensor $C_{ijkl}$, $\mathbf{n}$ for the vector $n_{k}$, and $\bar{\mathbf{v}}$ for the vector $\bar{v}_{k}$. Tensor contraction to the right of $\mathbf{C}$ will represent contraction with respect to the last two indices and to the left will represent contraction with the the first two indices. Appropriate contractions are assumed without making them explicit. With this we have:

[TABLE]

where $\mathbf{a}=(\nabla\cdot\mathbf{C}\nabla\bar{\mathbf{v}}-\beta^{2}\mathbf{n}\mathbf{C}\mathbf{n}\bar{\mathbf{v}})$ and $\mathbf{b}=(\nabla\cdot\mathbf{C}\mathbf{n}\bar{\mathbf{v}}+\mathbf{n}\mathbf{C}\nabla\bar{\mathbf{v}})$. On the other hand we have:

[TABLE]

where $\mathbf{c}=(\nabla\cdot\nabla\bar{\mathbf{v}}\mathbf{C}^{*}-\beta^{*2}\mathbf{n}\bar{\mathbf{v}}\mathbf{C}^{*}\mathbf{n})$ and $\mathbf{d}=(\nabla\cdot\mathbf{n}\bar{\mathbf{v}}\mathbf{C}^{*}+\nabla\bar{\mathbf{v}}\mathbf{C}^{*}\mathbf{n})$. Therefore, in general, $|A\bar{v}|\neq|A^{*}\bar{v}|$, or $A$ is not a normal operator. However, for certain special cases $A$ becomes a normal operator and the spectral theorem applies to it. The stiffness tensor $\mathbf{C}$ can always be separated into its Hermitian and skew-Hermitian parts $\mathbf{C}=\mathbf{H}+\mathbf{N}$ where $H^{*}_{klij}=H_{ijkl}$ and $N^{*}_{klij}=-N_{ijkl}$.

III.1 Real $\beta$

$\mathbf{C}=\mathbf{H}$:

Consider the case when the skew-Hermitian part is zero ($\mathbf{N}=0$). In this case, $C_{klij}^{*}=C_{ijkl}$. and we have $\mathbf{c}=(\nabla\cdot\mathbf{C}\nabla\bar{\mathbf{v}}-\beta^{*2}\mathbf{n}\mathbf{C}\mathbf{n}\bar{\mathbf{v}})$ and $\mathbf{d}=(\nabla\cdot\mathbf{C}\mathbf{n}\bar{\mathbf{v}}+\mathbf{n}\mathbf{C}\nabla\bar{\mathbf{v}})$. If we further insist that $\beta$ is real ($\beta^{*}=\beta$) then $\mathbf{c}=\mathbf{a},\mathbf{d}=\mathbf{b}$ and, in fact, $|A\bar{v}|=|A^{*}\bar{v}|$ and $A$ is a normal operator. Moreover, in this case the adjoint operator is

[TABLE]

showing that $A^{*}=A$, or that $A$ is a self-adjoint operator. Since the operator $B$ is clearly self-adjoint (and, therefore, normal), the eigenvalue problem which emerges from assuming a real wavenumber in an elastodynamic system characterized by a stiffness tensor which respects the symmetry $C_{ijkl}^{*}=C_{klij}$ is self adjoint. The eigenvalues, $\lambda$, are, therefore, real and the corresponding frequencies $\omega=\pm\sqrt{-\lambda}$ can be either real or imaginary (but not complex with simultaneously nonzero real and imaginary parts). We further have:

[TABLE]

which is of the form $-\int s_{ij}C_{ijkl}s_{kl}^{*}\mathrm{d}\Omega$ when one identifies $s_{ij}\equiv\bar{v}_{i,j}^{*}-i\beta\bar{v}_{i}^{*}n_{j}$. If we assume that the stiffness tensor is positive-definite (as in conservative systems) then the integral is always less than zero showing that the operator $A$ is negative definite as well ($\langle A\bar{v},\bar{v}\rangle\leq 0$). Since the operator $B$ is clearly positive definite, this means that all the eigenvalues $\lambda$ of the generalized eigenproblem under real $\beta$ and a self-adjoint and positive definite $\mathbf{C}$ tensor will be negative. All the frequencies $\omega=\sqrt{-\lambda}$ will, therefore, be purely real. As an academic point which is a corollary of this analysis, if $\mathbf{C}$ is negative definite then all frequencies will be imaginary since all $\lambda$ will be positive and real.

$\mathbf{C}=\mathbf{N}$:

If the Hermitian part of the stiffness tensor is zero ($\mathbf{H}=0$) then $C^{*}_{klij}=-C_{ijkl}$. In this case, $\mathbf{c}=(-\nabla\cdot\mathbf{C}\nabla\bar{\mathbf{v}}+\beta^{*2}\mathbf{n}\mathbf{C}\mathbf{n}\bar{\mathbf{v}})$ and $\mathbf{d}=(-\nabla\cdot\mathbf{C}\mathbf{n}\bar{\mathbf{v}}-\mathbf{n}\mathbf{C}\nabla\bar{\mathbf{v}})$. If $\beta$ is such that $\beta^{*2}=\beta^{2}$ then we will have $\mathbf{c}=-\mathbf{a}$ and $\mathbf{d}=-\mathbf{b}$. Of course, in this situation we will have $\mathbf{c}\cdot\mathbf{c}^{*}=\mathbf{a}\cdot\mathbf{a}^{*},\mathbf{d}\cdot\mathbf{d}^{*}=\mathbf{b}\cdot\mathbf{b}^{*}$. Even in this case if $\beta$ is real then it is clear that $|A\bar{v}|=|A^{*}\bar{v}|$ and that $A$ is a normal operator. In this case, we have:

[TABLE]

showing that $A^{*}\neq A$. Therefore, unlike in the case where $\mathbf{N}=0$, $\mathbf{H}=0$ (for real $\beta$) leads to an eigenvalue problem which is not self-adjoint. It is, however, normal which means that the eigenvectors will be orthogonal and will form a complete basis. The eigenvalues, in this case, have no requirement of being real. Instead we have $A^{*}=-A$ and, therefore,

[TABLE]

showing that $\lambda^{*}=-\lambda$. This, in turn, means that all eigenvalues of the problem will be strictly imaginary. Since $\omega=\sqrt{-\lambda}$, the corresponding frequencies will be complex with equal nonzero real and imaginary parts.

$\mathbf{C}=\mathbf{N}+\mathbf{H}$:

In this case, the tensors $\mathbf{a},\mathbf{b}$ can be divided into two parts: one resulting from $\mathbf{H}$ and the other resulting from $\mathbf{N}$. We can write $\mathbf{a}=\mathbf{a}^{h}+\mathbf{a}^{n}$ and similarly for $\mathbf{b}$. Carrying out the same decomposition for $\mathbf{c},\mathbf{d}$, we can show that $\mathbf{c}=\mathbf{a}^{h}-\mathbf{a}^{n}$ and $\mathbf{d}=\mathbf{b}^{h}-\mathbf{b}^{n}$. As $\mathbf{c}\cdot\mathbf{c}^{*}\neq\mathbf{a}\cdot\mathbf{a}^{*}$ and $\mathbf{d}\cdot\mathbf{d}^{*}\neq\mathbf{b}\cdot\mathbf{b}^{*}$ then $|A\bar{v}|^{2}\neq|A^{*}\bar{v}|^{2}$ and the operator $A$ is not normal. In this case, the eigenvalues $\lambda$ are complex which result in complex frequencies.

III.2 Imaginary and complex $\beta$

Assuming $\beta=i\beta_{\Im}$, we have $\mathbf{c}=(\nabla\cdot\mathbf{C}\nabla\bar{\mathbf{v}}+\beta_{\Im}^{2}\mathbf{n}\mathbf{C}\mathbf{n}\bar{\mathbf{v}})$, $\mathbf{d}=(\nabla\cdot\mathbf{C}\mathbf{n}\bar{\mathbf{v}}+\mathbf{n}\mathbf{C}\nabla\bar{\mathbf{v}})$, $\mathbf{a}=(\nabla\cdot\mathbf{C}\nabla\bar{\mathbf{v}}+\beta_{\Im}^{2}\mathbf{n}\mathbf{C}\mathbf{n}\bar{\mathbf{v}})$, and $\mathbf{b}=(\nabla\cdot\mathbf{C}\mathbf{n}\bar{\mathbf{v}}+\mathbf{n}\mathbf{C}\nabla\bar{\mathbf{v}})$. then $\mathbf{c}=\mathbf{a},\mathbf{d}=\mathbf{b}$, however as $i\beta^{*}\neq i\beta$, then from Eqs.(20,21), we have $|A\bar{v}|^{2}\neq|A^{*}\bar{v}|^{2}$. When $\beta=\beta_{\Re}+i\beta_{\Im}$, then $\mathbf{c}\neq\mathbf{a},\mathbf{d}\neq\mathbf{b}$ and consequently the operator is not normal when $\beta$ is not real. In these cases, the eigenvalues are complex.

III.3 $\omega(\beta)$ solutions using PWE

The $\omega(\beta)$ solutions can be easily found by using the PWE method and the approach is standard. Consider a periodic structure with the reciprocal-lattice vectors $\mathbf{G}=(G_{1},G_{2},G_{3})$. The material properties and field variables can be expanded using Fourier series as followsHo et al. (1990); Deymier (2013):

[TABLE]

where $\alpha$ can be any of $\{\rho,\mathbf{C},\mu,\lambda\}$, $\mathbf{f}$ can be $\{\mathbf{u},\boldsymbol{\sigma}\}$, $\mathbf{r}=(x_{1},x_{2},x_{3})$ is the position vector, and $\mathbf{K}=\beta\mathbf{n}$ is the wave vector. Substituting the material properties and displacement field in the elastodynamics equation of motion we have:

[TABLE]

Assuming specific values of $\mathbf{K}$, the above can be solved as a generalized matrix eigenvalue problem. The calculated bandstructures for in-plane wave propagation in a square unit cell of length $1\ (m)$ with a circular hole of radius $0.25\ (m)$ is shown in Fig. (1). In all the sub-figures, we have assumed isotopic material properties. Fig. (1a) shows the bandstructure when $C=H$. In this case, $H$ has been defined as a self-adjoint elasticity tensor emerging from two independent material constants ($E=20e9\ (Pa),\nu=0.25$). The frequency eigenvalues are all real, as expected from the discussion above. Fig. (1b) shows the real and imaginary parts of the calculated frequency when $\mathbf{C}=\mathbf{N}$. In this case $\mathbf{N}$ has been defined ($E=20e9i\ (Pa),\nu=0.25$). In this case, the real and imaginary parts of the frequencies come out to be the same, as expected by the theoretical arguments. In Figs. (1c,d) the real and imaginary parts of the frequency are plotted when $10\%$ loss is added to the Young modulus thus giving rise to a $\mathbf{C}$ tensor which is of the form $\mathbf{H+N}$. In this case, we end up with both real and imaginary parts of the frequency eigenvalues which are not necessarily equal.

IV $\beta(\omega)$ solutions for a wave at a given frequency

As a slight modification of the problem, we can seek $\beta$ solutions given frequency $\omega$ and a direction $\mathbf{n}$. The relevant eigenvalue problem can be written in several ways but we will write it in a form in which the $B$ operator, at least, has certain desirable properties. Consider the following:

[TABLE]

It will be assumed that there is minor symmetry on the compliance tensor which automatically enforces the usual requirement that $u_{i,j}=u_{j,i}$. Now considering $\mathbf{u}=\bar{\mathbf{u}}\exp(-i\beta n_{i}x_{i})$ and $\boldsymbol{\sigma}=\bar{\mathbf{s}}\exp(-i\beta n_{i}x_{i})$:

[TABLE]

which after some rearrangement can be written as:

[TABLE]

where $\bar{{\phi}}\equiv\{\bar{\mathbf{u}}\;\;\bar{\mathbf{s}}\}^{T}$, $\lambda=i\beta$ and the linear operators are given by:

[TABLE]

or in matrix form:

[TABLE]

In this form, $B$ is clearly self-adjoint. The adjoint operator for $A$ is:

[TABLE]

showing that $A$ is not self-adjoint even if $D^{*}_{klij}=D_{ijkl}$. To check whether the operator $A$ is normal or not, we have:

[TABLE]

where $|a_{i}|^{2}=a_{i}a_{i}^{*}$ and $|b_{ij}|^{2}=b_{ij}b_{ij}^{*}$. Similarly we have

[TABLE]

showing that in general $|A\bar{\phi}|^{2}\neq|A^{*}\bar{\phi}|^{2}$ and thus the operator $A$ is not normal. However, it is possible to surmise one situation in which the operator becomes normal. This would require the symmetry $D_{klij}^{*}=-D_{ijkl}$ and would additionally require that density be a second order tensor as well with a similar major symmetry. Specifically, $\rho_{ji}^{*}=-\rho_{ij}$, in addition to the aforementioned symmetry on $\mathbf{D}$ will render $A$ to be a normal operator (but not self-adjoint). In such a case we will have $A^{*}=-A$.

IV.1 $\beta(\omega)$ solutions using PWE

This problem is generally a quadratic eigenvalue problem, but as shown above we can re-write it as a mixed linear eigenvalue problem. For numerical solutions we will use the following form of the eigenvalue problem:

[TABLE]

By substituting the Fourier expansion of stress and displacement into the above equation we have:

[TABLE]

Multiplying both sides of the equations by $e^{-i({\hat{\mathbf{G}}}+\mathbf{K}).\mathbf{r}}$, integrating over the unit cell, and separating the terms which contain $\beta$, we have the following:

[TABLE]

which is the matrix eigenvalue problem whose solutions give us the wavenumber values at assumed values of frequency and wave propagation direction.

IV.1.1 Out-of-plane waves in a 2-D phononic crystal

To use this formulation in an example, we calculate the complex bandsructure for out-of-plane waves traveling in a square unit cell in the $\Gamma-X$ direction. The unit cell has side of $1\ (m)$ with a central circular region of radius $0.3\ (m)$. The material properties of the matrix are $\rho_{0}=1100\ (kg/m^{3})$, $\mu_{0}=1.33e9\ (Pa)$ and for the inclusion $\rho_{1}=7630\ (kg/m^{3})$ and $\mu_{1}=82e9\ (Pa)$ which are taken from Ref.Laude et al. (2009). The complex bandstructure shown in Fig. (2a) is for the case of linear elastic materials. In this case, only two evanescent wave mode are plotted (among infinite number of evanescent waves which exist). The first branch has zero imaginary part in the pass band and a non zero part in the band gap (where the real part is equal to $\pi$ or 0). This branch is the traditionally studied branch in the phononics literature. The second branch is purely imaginary at all frequencies and has been noted in other papers Laude et al. (2009). It must further be noted that since this is an out-of-plane example, we do not expect complex wavenumber solutions due to the monocoupled nature of the problem and because the materials are linear elastic Mead (1975). However, the results are still consistent with the non-normality of $A$ since the eigenvalues are $i\beta$. If this was an in-plane problem, we would expect fully complex wavenumber solutions.

The bandstructure shown in Figs. (2b-c) are for the unit cell of length $1\ m$ and the inclusion of radius $0.41\ m$ with the following complex shear modulus for the matrix phase (the inclusion properties are unchanged):

[TABLE]

where $\eta=80\ (Pa.s)$ and $\mu_{0}=1.33e9\ (Pa)$. These properties are taken from Ref. Moiseyenko and Laude (2011) . As the stiffness tensor is frequency dependent in viscoelastic materials, the $\omega(\beta)$ formulation cannot solve this problem directly. However, it is easy to solve this problem using the linear eigenvalue formulation of $\beta(\omega)$. In some studiesZhu et al. (2016); Zhao and Wei (2009) only the real part of the shear modulus is used to compare the bandgap behavior, however in Fig. (2-b-c) both storage and loss moduli are used. As one can see, for the viscoelastic case, pure real solutions are not predicted and the eigenvalues are complex. These results are in good agreement with Refs. Moiseyenko and Laude (2011); Laude et al. (2009).

IV.2 $\beta_{3}(\omega,\beta_{\alpha})$ solutions

Now we consider another kind of the phononic eigenvalue problem which one of its applications is to find the scattered wave field when a wave is incident at the interface of metamaterials Srivastava and Willis (2017). In this case we are given the frequency $\omega$ and two of the three components of the wavevector $\boldsymbol{\beta}=\{\beta_{1},\beta_{2},\beta_{3}\}$. Without any loss of generality we assume that $\beta_{1},\beta_{2}$ are known. The fields are of the form $u=\bar{u}\exp(-i\beta_{i}x_{i})$ and $\sigma=\bar{s}\exp(-i\beta_{i}x_{i})$. Denoting by Greek letters the indices $1,2$ and by roman the indices $1,2,3$ we can write the equation of motion as:

[TABLE]

and the constitutive relation as:

[TABLE]

These equations can be cast in the generalized eigenvalue form by identifying two new vectors. Specifically we consider $\boldsymbol{\gamma}=\{\beta_{1},\beta_{2},0\}$ and $\beta_{3}\mathbf{n}=\beta_{3}\{0,0,1\}$. This renders:

[TABLE]

or

[TABLE]

where $\bar{\phi}\equiv\{\bar{\mathbf{u}},\bar{\mathbf{s}}\}$, $\lambda=i\beta_{3}$, and $A,B$ are defined as:

[TABLE]

in matrix form:

[TABLE]

$B$ is clearly self-adjoint. For $A$ the adjoint operator is:

[TABLE]

showing that $A\neq A^{*}$. In general the operator is not normal either as can be seen from the following:

[TABLE]

IV.3 PWE solution

For numerical solutions we will begin with the following form of the problem:

[TABLE]

By substituting Eqs. (26,27) into Eq. (IV.3) we have the following:

[TABLE]

Assuming that $\beta_{1}$ and $\beta_{2}$ are given, we can re-write the above equation in the following way to form a generalized eigenvalue with $\beta_{3}$ as the eigenvalues and $\{\mathbf{u},\boldsymbol{\sigma}\}$ as the eigenvectors:

[TABLE]

multiplying both sides of the equations by $e^{-i({\hat{\mathbf{G}}}+\mathbf{K}).\mathbf{r}}$ and integrating over the unit cell, Eqs. (55,IV.3) are written as follows:

[TABLE]

which is a generalized matrix eigenvalue problem whose solutions are $\beta_{3}$.

IV.3.1 In-plane waves in a 1-D phononic crystal

To show a numerical example, we consider an in-plane wave propagation problem in periodic layered composite as an example (Fig. 3). Anti-plane shear waves in similar structures have been shown to result in exotic wave phenomenon Willis (2016); Srivastava (2016) and in-plane waves have also been considered recently Shmuel and Srivastava (2018). Each repeated unit cell has the following material property and layer thickness configuration:

• •

$E_{1}=7\ (GPa),\;\nu_{1}=0.25,\;\rho_{1}=2700\ (kg/m^{3}),\;h_{1}=2/3\ (m),$

• •

$E_{2}=4.34\ (GPa),\;\nu_{2}=0.36,\;\rho_{2}=1180\ (kg/m^{3}),\;h_{2}=1/3\ (m)$

Studying some scattering problems in such a structure requires the knowledge of admissible $\beta_{2}$ values for given $\beta_{1},\omega$ values Srivastava (2016). The numerical method of choice for solving these cases is the Transfer Matrix Method (TMM) Haque and Shim (2016) which is an exact method. However, TMM solves for $\beta_{1}$ given $\beta_{2},\omega$ thus requiring further numerical searches. If the eigenvalue problem is cast in the $\beta_{2}(\omega,\beta_{1})$ form as we have done above, then it could be directly solved. Here, TMM is first employed to produce exact solutions at $\omega=3000\ (rad/s)$. At this frequency, TMM shows that there are two $\beta_{1}$ solutions for given $\beta_{2}$ values. Depending upon whether $\beta_{2}$ is real, imaginary, or complex, the corresponding $\beta_{1}$ solutions can also be complex. These two solutions are shown in Fig. (3b) in the form of two sheets of continuous surfaces (only the real part of $\beta_{1}$ are plotted). Out of these multitude of solutions, there appear certain $\beta_{2}$ values for which $\beta_{1}$ is real. These real $\beta_{1}$ solutions are extracted from the TMM solutions using a grid search and marked by red dots in Fig. (3b). These correspond to propagating modes in the $x_{1}$ direction. However, in the $x_{2}$ direction, these modes can be propagating (real $\beta_{2}$), evanescent (imaginary $\beta_{2}$), or non-propagating (complex $\beta_{2}$). This is in contrast with the anti-plane shear case Srivastava (2016), where only real or imaginary $\beta_{2}$ are associated with real $\beta_{1}$. Fig. (3c) shows that there are two propagating modes in the $x_{2}$ direction for the considered range of $\beta_{1}$ and Fig. (3d) shows the evanescent modes in the $x_{2}$ direction. In addition, there are non-propagating modes in the $x_{2}$ direction which are shown in Fig. (3b) but not plotted separately. The red dots in Fig. (3b), therefore, representing all modes which are propagating in the $x_{1}$ direction and can be extracted from TMM after considerable effort. These solutions are critical for solving interface problems Shmuel and Srivastava (2018). Fig. (3b) also shows the PWE solutions in black dots where the $\beta_{2}(\omega,\beta_{1})$ problem is directly solved by varying $\beta_{1}$ between 0 and $\pi$. It shows that all the $\beta_{2}$ solutions (real, imaginary, and complex) are computed directly, as is evidenced by the comparison between the locations of the black and red dots. In the PWE solution, 73 plane waves were used in the stress and displacement expansion.

As a further point of interest here and relating to the overall context of this work, we explicitly show the bi-orthogonality condition mentioned earlier. As shown by Eqs. (49), (50) and (51), the operator for the $\beta_{2}(\beta_{1},\omega)$ problem is not self-adjoint, therefore, there is no orthogonality for right or left eigenvectors themselves. However, since there is no exceptional point in the ranges of $\beta_{1},\beta_{2}$ considered, there is no need to create generalized eigenvectors either and the bi-orthogonality of the left and right eigenvectors should still hold. This is verified in Tab. 1. For $\beta_{1}=\pi/6$, $u_{(l,r)}^{1}$ and $u_{(l,r)}^{2}$ are the left and right eigenvectors, respectively, corresponding to the two $\Re\beta_{2}$ modes. Similarly, $u_{(l,r)}^{3}$ and $u_{(l,r)}^{4}$ correspond to the first two $\Im\beta_{2}$ modes, and $u_{(l,r)}^{5}$ and $u_{(l,r)}^{6}$ correspond to the first two complex $\beta_{2}$ modes. The table shows that while \Big{|}u_{(l)}^{i\dagger}Bu_{r}^{j}/\big{|}u_{(l,r)}^{i\dagger}Bu_{r}^{i}\big{|}\Big{|}=\delta_{ij}, \Big{|}u_{(r)}^{i\dagger}Bu_{r}^{j}/\big{|}u_{(l,r)}^{i\dagger}Bu_{r}^{i}\big{|}\Big{|}\neq\delta_{ij} thus showing that bi-orthogonality holds with respect to the $B$ operator.

V Conclusion

In this paper we consider the operator properties of various phononic eigenvalue problems and aim to answer some fundamental questions about the eigenvalues and eigenvectors of phononic operators. These include questions about the potential real and complex nature of the eigenvalues, whether the eigenvectors form a complete basis, what are the right orthogonality relationships, and how to create a complete basis when none may exist at the outset. In doing so we present a unified understanding of the properties of the eigenvalues and eigenvectors which would emerge from any numerical method employed to compute such quantities. We show that the phononic problem can be cast into linear eigenvalue forms from which such quantities as frequencies ($\omega(\boldsymbol{\beta})$), wavenumbers ($\beta(\omega,\mathbf{n})$), and desired components of wavevectors ($\beta_{3}(\omega,\beta_{\alpha})$) can be directly ascertained without resorting to searches or quadratic eigenvalue problems and that the relevant properties of such quantities can be determined apriori through the analysis of the associated operators. For the $\omega(\boldsymbol{\beta})$ problem we show that the associated operators may be normal under certain assumed symmetries of the stiffness tensor. With additional definiteness assumptions, the frequency eigenvalues can be shown to be purely real or imaginary. For the $\beta(\omega,\mathbf{n})$ and $\beta_{3}(\omega,\beta_{\alpha})$ problems we show that the associated operators are not normal thus giving rise to complex eigenvalues in general. In each case we show how the Plane Wave Expansion (PWE) method may be extended to solve the associated eigenvalue forms and we present associated numerical examples for each eigenvalue form.

Acknowledgements.

A.S. acknowledges support from the NSF CAREER grant #1554033 to the Illinois Institute of Technology and NSF grant #1825354 to the Illinois Institute of Technology. We would also like to thank Prof. Igor Cialenco at IIT for discussions and guidance on the spectral theorem.