Level Repulsion and Band Sorting in Phononic Crystals

TL;DR

This paper introduces a fast, accurate method for distinguishing level repulsion from normal crossings in phononic crystals, crucial for correct band sorting and property analysis, by analyzing eigenvector sign flips around exceptional points.

Contribution

The paper proposes a novel, efficient technique to identify exceptional points in phononic band structures, improving speed and accuracy over existing methods.

Findings

The method correctly identifies over 97% of cases.

Eigenvector sign flips are more complex in phononic crystals than in simple matrices.

The approach is applicable to general eigenvalue problems.

Abstract

In this paper we consider the problem of avoided crossings (level repulsion) in phononic crystals and suggest a computationally efficient strategy to distinguish them from normal cross points. This process is essential for the correct sorting of the phononic bands and, subsequently, for the accurate determination of mode continuation, group velocities, and emergent properties which depend on them such as thermal conductivity. Through explicit phononic calculations using generalized Rayleigh quotient, we identify exact locations of exceptional points in the complex wavenumber domain which results in level repulsion in the real domain. We show that in the vicinity of the exceptional point the relevant phononic eigenvalue surfaces resemble the surfaces of a 2 by 2 parameter-dependent matrix. Along a closed loop encircling the exceptional point we show that the phononic eigenvalues are exchanged, just as they are for the 2 by 2 matrix case. However, the behavior of the associated eigenvectors is shown to be more complex in the phononic case. Along a closed loop around an exceptional point, we show that the eigenvectors can flip signs multiple times unlike a 2 by 2 matrix where the flip of sign occurs only once. Finally, we exploit these eigenvector sign flips around exceptional points to propose a simple and efficient method of distinguishing them from normal crosses and of correctly sorting the band-structure. Our proposed method is roughly an order-of magnitude faster than the zoom-in method and correctly identifies > 97% of the cases considered. Both its speed and accuracy can be further improved and we suggest some ways of achieving this. Our method is general and, as such, would be directly applicable to other eigenvalue problems where the eigenspectrum needs to be correctly sorted.

Tables

Table 1. Table 1: The absolute values of inner products of p 1 − 6 subscript 𝑝 1 6 p_{1-6} and q 1 − 6 subscript 𝑞 1 6 q_{1-6} .

	$q_1$	$q_2$	$q_3$	$q_4$	$q_5$	$q_6$
$p_1$	0.9999	4.8268e-12	0.0033	1.6292e-13	0.0012	3.4603e-5
$p_2$	4.4115e-12	1.0000	5.3649e-13	0.0033	5.0022e-14	7.8003e-14
$p_3$	0.0033	4.4004e-13	0.9985	9.9670e-13	0.0453	0.0282
$p_4$	1.7215e-13	0.0034	5.0500e-13	0.9999	3.7333e-13	9.8231e-13
$p_5$	9.6713e-4	1.0209e-13	0.0432	2.5041e-13	0.9965	0.0676
$p_6$	2.3906e-6	5.1212e-14	0.0311	9.7758e-13	0.0662	0.9972

Equations

$$\mathbf{\Lambda }(\mathbf{u})+\mathbf{f}=\lambda \mathbf{u};\mathbf{\Lambda }(\mathbf{u})\equiv \frac{1}{\rho }{\left[C_{ijkl}{u}_{k,l}\right]}_{,j}$$

$$⟨\mathbf{v},\mathbf{w}⟩={\int }_{\Omega }\rho {\mathbf{v}}^{\ast }\cdot \mathbf{w}\mathrm{d}\Omega$$

$$C_{jkmn}(\mathbf{x}+n_i{\mathbf{h}}^i)=C_{jkmn}(\mathbf{x});\rho (\mathbf{x}+n_i{\mathbf{h}}^i)=\rho (\mathbf{x})$$

$${\mathbf{\Lambda }}^{(\mathbf{k})}({\mathbf{u}}^p)=\lambda {\mathbf{u}}^p$$

$${\mathbf{\Lambda }}^{(\mathbf{k})}({\mathbf{u}}^p)\equiv \frac{1}{\rho }{\left[C_{ijkl}{u}_{k,l}^p\right]}_{,j}+\frac{\mathrm{i}{q}_j{C}_{ijkl}}{\rho }{u}_{k,l}^p+\frac{\mathrm{i}{q}_l}{\rho }{\left[C_{ijkl}{u}_k^p\right]}_{,j}-\frac{q_l{q}_j{C}_{ijkl}}{\rho }{u}_k^p$$

$$⟨{\mathbf{v}}^p,\mathbf{\Lambda }({\mathbf{w}}^p)⟩=⟨\mathbf{\Lambda }({\mathbf{v}}^p),{\mathbf{w}}^p⟩$$

$$⟨u_{n;i}^p(\mathbf{x}),(u_{n;j}^p{)}^{\ast }({\mathbf{x}}^{\prime })⟩={\delta }_{ij}\delta (\mathbf{x}-{\mathbf{x}}^{\prime })$$

$$⟨{\mathbf{u}}_m^p,{\mathbf{u}}_n^p⟩={\delta }_{mn}$$

$$\left[{\mathbf{\Lambda }}^{(\mathbf{k})}\right]\{{\mathbf{U}}^{\mathbf{n}}\}=\lambda \{{\mathbf{U}}^{\mathbf{n}}\}$$

$$⟨{\mathbf{u}}_a^p,{\mathbf{u}}_b^p⟩=⟨U_i^{\mathbf{m};a}{e}^{\mathrm{i}{\mathbf{G}}^{\mathbf{m}}\cdot \mathbf{x}},U_i^{\mathbf{n};b}{e}^{\mathrm{i}{\mathbf{G}}^{\mathbf{n}}\cdot \mathbf{x}}⟩={\delta }_{ab}$$

$${\int }_{\Omega }{e}^{\mathrm{i}[{\mathbf{G}}^{\mathbf{m}}-{\mathbf{G}}^{\mathbf{n}}]\cdot \mathbf{x}}d\Omega =V{\int }_0^1{\int }_0^1{\int }_0^1{e}^{\mathrm{i}2\pi [\alpha H_1+\beta H_2+\gamma H_3]}d{H}_{1}d{H}_{2}d{H}_3=V{\delta }_{\mathbf{mn}}$$

$$⟨U_i^{\mathbf{m};a},U_i^{\mathbf{m};b}⟩=\frac{{\delta }_{ab}}{V}$$

$$⟨U_i^{\mathbf{m};a},U_i^{\mathbf{m};b}⟩={\delta }_{ab}$$

$$\sum _{i=1}^3{\left[U_i^{\mathbf{m};a}\right]}^{\ast }{U}_i^{\mathbf{m};b}\approx {\delta }_{ab}$$

$$\mathbf{A}=\left[\begin{array}{cc}{ϵ}_1& 0\\ 0& {ϵ}_2\end{array}\right]+\mu \mathbf{U}\left[\begin{array}{cc}{\omega }_1& 0\\ 0& {\omega }_2\end{array}\right]{\mathbf{U}}^{\mathrm{T}};\mathbf{U}=\left[\begin{array}{cc}\cos \varphi & -\sin \varphi \\ \sin \varphi & \cos \varphi \end{array}\right]$$

$${\mu }_c=-\frac{{ϵ}_1-{ϵ}_2}{{\omega }_1-{\omega }_2}\exp (\pm 2\mathrm{i}\varphi )$$

Full text

I Introduction

Phononic crystals are artificially designed periodic composite materials Pennec et al. (2010); Lee et al. (2012); Hussein et al. (2014); Srivastava (2015) which can significantly affect the propagation characteristics of acoustic and stress waves through them in 1- Nemat-Nasser (2015); Willis (2016); Srivastava (2016); Srivastava and Willis (2017); Shmuel and Band (2016); Chen and Elbanna (2016), 2- Shim et al. (2015); Mousanezhad et al. (2015), and 3-dimensions Babaee et al. (2015); Galich et al. (2016); Lu et al. (2017). Although the behavior of phononic crystals can be analyzed through full Finite Element simulations, a vitally important tool in phononics research is band-structure calculation. In its very basic form, the phononic band-structure is the dispersion relation of the phononic crystal derived by calculating the real frequency eigenvalues ($\omega$) of the equation of motion (acoustic or elastodynamic) for given values of real wavevector ($\mathbf{k}$) which span the boundaries of the crystal’s Irreducible Brillioun Zone (IBZ). It should be noted, however, that more complex bandstructures with complex values of frequency and/or wavevector are possible as wellTrainiti and Ruzzene (2016); Frazier and Hussein (2016); Mazzotti et al. (2017); Krattiger et al. (2016). The (real $\omega$-real $\mathbf{k}$) band-structure (phononic or photonic) represents a parametric and Hermitian eigenvalue problem where the frequency eigenvalue depends parametrically upon the wavevector. Any parametric eigenvalue problem exhibits two phenomena that are of interest to us Kato (1966); Heiss and Sannino (1990, 1991):

• •

Diabolic points: These are values of the parameter in the real domain which exhibit 2 degenerate eigenvalues whose associated eigenvectors are, nevertheless, linearly independent.

• •

Exceptional points: These are values of the parameter in the complex domain of the parameter which exhibit degenerate eigenvalues as well as degenerate associated eigenvectors.

It is well known that eigenvalues are repulsed (level repulsionWu and Huang (2004); Yeh et al. (2016)) from each other in the vicinity of exceptional points Heiss and Sannino (1990). Therefore, the existence of exceptional points in the vicinity of the real domain of the parameter will cause the associated eigenvalues to avoid each other in the real domain. Furthermore, in the real domain the eigenvectors of the associated eigenvalues will be rapidly exchanged in the level repulsion region. This rapid exchange of eigenvectors in the level repulsion region has been used to practical effect in the design of highly sensitive sensors in discrete and lower dimensional systemsManav et al. (2014); Foreman and Vollmer (2013). In the context of phononics and photonics, diabolic points correspond to real crosses of eigenvalue branches in the band-structure and exceptional points are associated with avoided crossings of the branchesDing et al. (2016); Lin et al. (2016). Branch continuations in the vicinity of crosses (real or apparent) become difficult to determine with coarse wavevector discretizations and it vitally affects some application areas of phononics. One area of such relevance is in employing phononic crystals for designing materials with extreme thermal propertiesDavis and Hussein (2014); Xiong et al. (2016). It has been shown that the correct sorting of the phononic band-structure is of prime importance in the accurate prediction of thermal properties of phononic crystalsZen et al. (2014); Anufriev and Nomura (2015). Since thermal properties directly depend upon the group velocities of phononic bands, incorrect band sorting can lead to more than 40% error in the predicted values of phonon emission power calculated over a coarse wave-vector gridZen et al. (2014).

The techniques which have been applied for phononics band-sorting in literature are essentially based upon eigenvector orthogonalityHuang et al. (2014); Zen et al. (2014) and polarizations of modesWu and Huang (2004); Achaoui et al. (2010). However, since eigenvectors and polarizations are exchanged in level repulsion regions, all such techniques are bound to fail in distinguishing real crosses from avoided crosses as we clearly show in this paper. The only fail-safe method proposed in literature to make this distinction is to calculate the eigenvalues of the system over very fine parametric discretizations in the crossing zones, thereby, explicitly distinguishing real crosses from apparent onesWu and Huang (2004); Achaoui et al. (2010). This process is highly computationally inefficient and simply infeasible in cases where thousands of such cross points might need to be considered, as is the case for thermal properties.

In this paper, we propose a new and highly computationally efficient method of distinguishing real crosses from avoided ones for the phononic bandstructure. The technique can be extended to any other system where eigenvalue bands need to be correctly sorted. Our method is based on the principal that eigenvectors suffer flips of signs along closed contours around exceptional pointsHeiss (1999, 2000). A large enough closed curve near a cross will either contain an exceptional point (if it exists) or not (if it doesn’t exist). Therefore, by calculating the eigenvectors at some points on this closed contour we can distinguish real crosses from apparent ones with near certainty. The method does not require us to actually find the location of the exceptional point as it only depends upon making the closed contour large enough to encircle it.

In the subsequent sections we first recast the phononic problem in a self-adjoint form ensuring that its associated finite dimensional representation will be Hermitian for real values of the wavevector. We then use the generalized Rayleigh quotient to explicitly identify an exceptional point in the complex wavevector domain. We show that in the vicinity of the exceptional point, the behavior of the eigenvalues resembles that of a $2\times 2$ matrix. As expected, the eigenvalues are exchanged along a closed curve around it. However, unlike for a $2\times 2$ matrix, for some exceptional points we show that the eigenvectors flip signs multiple times. We show that no such flips occur when the closed curve either does not encircle an exceptional point or pass a branch cut. Finally, we propose our method for distinguishing real crosses from apparent ones and use it to correctly sort an extended phononic band-structure.

II Problem Description

The frequency domain dynamics of a linear elastic medium with a spatially varying constitutive tensor $\mathbf{C}$ and density $\rho$ is given by:

[TABLE]

where $\lambda=-\omega^{2}$, $\mathbf{u}$ is the displacement field, $\mathbf{f}$ is the body force, and $\boldsymbol{\Lambda}$ is a linear differential operator. The operator operates on the Hilbert space of square integrable vector functions over which the inner product is defined with respect to the weight function $\rho$:

[TABLE]

It can be shown that the operator is self adjoint ($\langle\mathbf{v},\boldsymbol{\Lambda}(\mathbf{w})\rangle=\langle\boldsymbol{\Lambda}(\mathbf{v}),\mathbf{w}\rangle$). For phononic problems the problem domain is periodic and is defined by a repeating unit cell. For the general 3-dimensional case, the unit cell ($\Omega$) is characterized by 3 base vectors $\mathbf{h}^{i}$, $i=1,2,3$. Any point within the unit cell can be uniquely specified by the vector $\mathbf{x}=H_{i}\mathbf{h}^{i}=x_{i}\mathbf{e}^{i}$ where $\mathbf{e}^{i}$ are the orthogonal unit vectors and $0\leq H_{i}\leq 1,\forall i$. The unit cell is associated with a set of reciprocal base vectors $\mathbf{q}^{i}$ such that $\mathbf{q}^{i}\cdot\mathbf{h}^{j}=2\pi\delta_{ij}$. Reciprocal lattice vectors are represented as a linear combination of the reciprocal base vectors, $\mathbf{G}^{\mathbf{n}}=n_{i}\mathbf{q}^{i}$, where $n_{i}$ are integers. Fig. (1) shows the schematic of a 2-D unit cell, clearly indicating the unit cell basis vectors, the reciprocal basis vectors and the orthogonal basis vectors.

The material properties have the following periodicity:

[TABLE]

where $n_{i}(i=1,2,3)$ are integers. Due to the periodic nature of the problem, it accepts solutions of the form $\mathbf{u}(\mathbf{x})=\mathbf{u}^{p}(\mathbf{x})e^{\mathrm{i}\mathbf{k}.\mathbf{x}}$ where $\mathbf{k}$ is the Bloch wavevector with components $\mathbf{k}=Q_{i}\mathbf{q}^{i}$ where $0\leq Q_{i}\leq 1,\forall i$ and $\mathbf{u}^{p}$ is a periodic function. Under the substitution, the harmonic elastodynamic problem can be formally written as (neglecting the body force):

[TABLE]

where the superscript $(\mathbf{k})$ is now included to emphasize that the operator depends upon the Bloch wavevector. Explicitly we have:

[TABLE]

It can be shown that for real values of the wavevector the operator is self adjoint:

[TABLE]

It, therefore, has real eigenvalues $\lambda_{n}$ for real values of the wavevector. Furthermore, its eigenfunctions , $\tilde{\mathbf{u}}_{n}^{p}$, form a complete and orthonormal set:

[TABLE]

For a suitable span of the wavevector, the sets of eigenvalues $\lambda_{n}$ (and the corresponding frequencies $\omega_{n}$) constitute the phononic dispersion relation of the composite. There are several numerical techniques for calculating the eigenvalues but a common method is to expand the field variable $\mathbf{u}^{p}$ in an appropriate basis and then use the basis to convert the differential equation into a set of linear equations. Both the Plane Wave ExpansionKushwaha et al. (1993) method and Rayleigh quotientLu and Srivastava (2016) follow this strategy. As an example, trigonometric basis can be used to write $\mathbf{u}^{p}=\mathbf{U}^{\mathbf{n}}e^{\mathrm{i}\mathbf{G^{n}}\cdot\mathbf{x}}$ where Einstein summation is implied. $\mathbf{U}^{\mathbf{n}}$ is a constant vector with components $U^{\mathbf{n}}_{i};\;i=1,2,3$. Each component depends upon the appropriate Fourier component which, in 3-dimension, refers to the set $(n_{1}n_{2}n_{3})$. Either through a variational law or by expanding materials properties in the same basis, we arrive at a system of linear equations formally written as:

[TABLE]

which is a simple matrix eigenvalue problem. The $a^{\mathrm{th}}$ eigenvalue $\lambda_{a}$ has the associated eigenvector $\mathbf{U}^{\mathbf{n};a}$ which provides the $a^{\mathrm{th}}$ eigenfunction of the original problem through $\mathbf{u}^{p}_{a}=\mathbf{U}^{\mathbf{n};a}e^{\mathrm{i}\mathbf{G^{n}}\cdot\mathbf{x}}$. The orthogonality condition now reads:

[TABLE]

The above can be simplified further by noting that the basis functions themselves are orthogonal in the following sense:

[TABLE]

where $\alpha,\beta,\gamma$ are integers and $V$ is the volume of the unit cell. Employing this we have the following expression for the orthogonality of the eigenvectors:

[TABLE]

For the rest of the paper we have normalized the eigenvectors such that

[TABLE]

II.1 Example

As an example Wu and Huang (2004), a 2-D phononic composite is considered whose unit cell is shown in Fig. (2a). The unit cell is composed of Nickle ($C_{11}=324\text{ GPa},C44=80\text{ GPa},\rho=8905\text{kg/m}^{3}$) matrix with embedded circular PMMA ($C_{11}=7.11\text{ GPa},C44=2.03\text{ GPa},\rho=1200\text{kg/m}^{3}$) inclusions in a square lattice. The lattice constant is $1mm$ with an inclusion fill fraction of $50\%$. Fig. (2b) depicts the corresponding phononic dispersion relation for in-plane modes along the boundary of the Irreducible Brillouin Zone (IBZ) (Fig. 2c). This band structure, which is rendered from the self-adjoint form of elastodynamic eigenvalue problem using Rayleigh quotient, shows very good agreement with previously published plane wave expansion (PWE) results Wu and Huang (2004). The figure shows two vertical lines which denote adjacent wavenumbers over which the eigenvalue computations have been performed. Within the frequency range under consideration, the eigenvalue solutions at these wavenumbers include 6 solution points ($p_{i},q_{i};\;i=1,...6$). The corresponding eigenvectors obviously satisfy Eq. (12) when they all belong to the same wavenumber. However, if the wavenumbers P and Q are close enough to each other then it is expected that Eq. (12) will still be approximately satisfied. In other words, if $\mathbf{U}^{\mathbf{m};a}$ belongs to P and $\mathbf{U}^{\mathbf{m};b}$ belongs to Q then we should still have:

[TABLE]

if $P,Q$ are close to each other. This approximate orthogonality of the eigenvectors of two neighboring wavenumbers can be used to sort the modes of the phononic bandstructure in regions which are away from the regions of level repulsion. The logic behind this is that the modes which form a single phononic band over the span of the wavevector will have associated eigenvectors which produce nearly unitary inner products with their neighboring modes on the same band. They will, on the other hand, produce nearly zero values when their inner products are taken with neighboring modes on other bands. As an example we present the calculated absolute values of the inner products between modes $p_{1-6}$ and $q_{1-6}$ in Table 1. The table confirms that approximate orthogonality holds for the considered modes.

After sorting using this approximate orthogonality, the dispersion curves appear to cross each other at ten points. However, not all of these cross points are true cross points.

To clarify this further we focus on the two points marked by T and R in Fig. (2) (and Figs. 3a,b). The modes around these cross points are marked as $t_{i}$ and $r_{i}$ respectively, where $i=1,\cdots,4$. $t_{1},t_{2}$ and $t_{3},t_{4}$ belong to two different wavevectors which are separated from each other by a distance of $6.6667\times 10^{-3}$ in the reduced wavenumber domain. The same situation applies to $r_{i}$. The inner products between $t_{1},t_{4}$ and $t_{2},t_{3}$ are both 0.9999. Between $t_{1},t_{3}$ and $t_{2},t_{4}$ they are close to zero. Similarly, the inner products between $r_{1},r_{4}$ and $r_{2},r_{3}$ are both 1.0000 and between $t_{1},t_{3}$ and $t_{2},t_{4}$ they are close to zero. The inner products indicate that at this level of wavevector discretization, the two points should be determined as real cross points. Figs. (3c,d) show zoomed-in views of the T and R points with the wavenumber regions shown in Figs. (3a,b) further discretized in about 1000 intermediate points. The zoomed-in views clearly show that while at T the two curves indeed cross each other, at R the curves move apart. R is an apparent cross point when the IBZ is roughly discretized but is, in reality, an avoided crossing. There are several other apparent cross points in the bandstructure in Fig. (2b). These are marked off with green circles whereas the real cross points are marked off with red circles. The apparent cross points can always be distinguished from the real cross points by calculating the bandstructure over a very fine grid of wavevectors. This has been the dominant approach in literatureWu and Huang (2004); Achaoui et al. (2010). However, this approach is computationally wasteful and we would like to be able to distinguish apparent cross points from real cross points for accurate band sorting over a rough grid of wavevectors.

III Level Repulsion in Phononic Crystals

Consider an $N\times N$ matrix $\mathbf{A}(\mu)$ which depends upon a parameter $\mu$. For real values of $\mu$ consider that this matrix is self-adjoint. It, therefore, has $N$ real eigenvalues $\lambda_{n},\;n=1,...N$ which are determined by the solution of the equation $\mathrm{det}(\mathbf{A}(\mu)-\lambda\mathbf{I})=0$ for real $\mu$. Now extend the domain of $\mu$ to the complex plane. In this case, $\mathbf{A}$ generally ceases to be self-adjoint and the eigenvalues become complex. In fact, it can be shown that the solution of $\mathrm{det}(\mathbf{A}(\mu)-\lambda\mathbf{I})=0$ is a multi-valued function $\lambda(\mu)$ which is complex valued and is defined over the complex planeHeiss and Sannino (1990, 1991). It has a Riemannian sheet structure with $N$ Riemann sheets which are connected with each other through square root branch points Heiss and Sannino (1990, 1991). The branch points are termed exceptional points. At these points we have the coalescing of both a set of eigenvalues and their associated eigenvectors. This is in contrast with such points as the cross point T in Fig. (3) where there is a degeneracy of two eigenvalues but the associated eigenvectors are still linearly independent. The latter is called a diabolic point. The presence of exceptional points in the complex $\mu$ plane is associated with avoided crossings (level repulsion) of $\lambda_{n}$ on the real $\mu$ axis.

As an example, if $\mathbf{A}$ is a $2\times 2$ matrix of the form

[TABLE]

then it has an exceptional point at

[TABLE]

$\mathbf{A}$ has two eigenvalues $\lambda_{1},\lambda_{2}$ whose variations in the vicinity of the exceptional point are shown in Fig. (4). The exceptional point itself is marked with a red circle in the real and imaginary domains of the eigenvalues. The imaginary figure (Fig. 4b) is flipped with respect to the imaginary $\mu$ axis for easy visualization of the exceptional point. Fig. (4a) shows both the coalescing of the two surfaces at the exceptional point and the associated phenomenon of level repulsion which follows on the real $\mu$ axis. The latter is emphasized by the red curve. On the real $\mu$ axis the complex parts of the eigenvalues are zero as expected.

Although Fig. (4) pertains to the behavior of the eigenvalues near the exceptional point of a $2\times 2$ matrix, it can be shown that a $2\times 2$ matrix captures the essential characteristics of an $N\times N$ matrix near an exceptional point if the exceptional point is isolated Heiss and Sannino (1991); Heiss (1999). The discretized phononics eigenvalue matrix is such an $N\times N$ problem. To make a clear correspondence, note that phononic eigenvalue problem is self-adjoint for real values of the wavevector. The wavevector has the form $\mathbf{k}=0.5\mathbf{q}^{1}+Q_{2}\mathbf{q}^{2}$ where $0\leq Q_{2}\leq 0.5$ at all points along the $X-M$ direction in the bandstrucure in Fig. (2b). The eigenvalue matrix, along the $X-M$ direction, which includes the avoided cross point R, is, therefore, only dependent upon the parameter $Q_{2}$. In other words $\mathbf{\Lambda}^{(\mathbf{k})}$ is of the form $\mathbf{\Lambda}(Q_{2})$ and is Hermitian (self-adjoint) as long as $Q_{2}$ is real (or follows the boundary of the IBZ). It has $N$ real eigenvalues $\lambda_{n}(Q_{2})$ for a given value of real $Q_{2}$. However, $\mathbf{\Lambda}(Q_{2})$ is non Hermitian for complex $Q_{2}$. We are interested in identifying the exceptional point which is giving rise to the avoided crossing at point R.

Fig. (5) shows the behavior of the frequencies corresponding to the two avoided crossing branches when complex values of $Q_{2}$ are considered around point R. The required complex eigenvalues are calculated by using the generalized Rayleigh quotientLu and Srivastava (2016). As in Fig. (4), the imaginary $Q_{2}$ axis for the imaginary part of the frequency surfaces (Fig. 5b) is flipped for the easy visualization of the coalescing point. The exceptional point is identified as $Q_{2}=0.1854865\pm\mathrm{i}2.6130959\mathrm{\times}10^{-6}$. Even though Fig. (5) corresponds to the eigenvalues of a much larger matrix, its similarity with Fig. (4) shows that the behavior of the eigenvalues near the exceptional point is similar to those of a $2\times 2$ matrix. This similarity in appearance would clearly not have held if the larger matrix had another exceptional point in the vicinityHeiss and Sannino (1991).

The exceptional point forms a square root branch pointHeiss and Sannino (1990, 1991); Heiss (1999) and connects the two eigen-surfaces. A closed curve in the wavenumber domain around the exceptional point would require the two eigenfrequencies to switch if analytic continuation is to be maintained. Figs. (6b,c) show the movement of the real and imaginary parts of the two phononic eigenvalues as the wavenumber traces a closed curve around the exceptional point (Fig. 6a). The exceptional point is indicated by the bold black dot whereas red and black dots represent the individual eigenvalues. There is a branch cut in Fig. (6a) which starts at the exceptional point and increases vertically along the positive imaginary $Q_{2}$ direction. The real parts of the eigenvalues are non-differentiable across this branch cut (Fig. 6b) whereas the imaginary parts are discontinuous (Fig. 6c). However, they can be analytically continued by switching the eigenvalues across the branch cut. This is evident from the smooth switching of the red curves into the black curves and vice-versa. If the closed curve in the wavenumber domain excludes the exceptional point (Fig. 6d) then the corresponding eigenvalue curves remain separated and analytic (Figs. 6e,f).

In addition to the eigenvalue exchange phenomenon described above, the corresponding eigenvectors also behave in characteristic ways along closed curves which encircle an exceptional point Heiss (2000). It is known, for instance, that for a simple $2\times 2$ matrix, one of the two (and only one) eigenvectors will flip sign around such a closed curve. For the more complex and higher-dimensional phononic case the situation is slightly different. To understand this difference, first consider Figs. (7a-c). Fig. (7a) shows the behavior of the two relevant phononic eigenvalues as two different closed curves are traversed in the vicinity of point R in Fig. (2b). The smaller curve encircles one exceptional point and the larger curve encircles the same exceptional point as well as its complex conjugate exceptional point. Both curves are circles with the smaller circle centered at the exceptional point and the larger one centered on the real axis. We now calculate the eigenvector inner products over adjacent wavevector points as we traverse these circles. Fig. (7b) shows these inner products for the two eigenvectors over the smaller circle (as a function of the polar angle.) We note that while one of the two inner products (black curve) is very close to 1 over the entire circle, the other curve (red) achieves a value of -1 as we cross the branch cut. Over the closed curve, therefore, one of the eigenvectors has flipped its sign. Over the larger circle (Fig. 7c) there are two sign flips for one of the eigenvectors essentially nullifying the flip. The location of these sign flips are also denoted by green dots in Fig. (7a). The situation is more complex for another exceptional point which exists at a higher frequency (Fig. 7d corresponding to point L1 in Fig. 2b). In this case, even for the smaller circle (Fig. 7e), there are multiple flips of sign for the two involved eigenvectors. However, even in this case the essential effect of the flips is to ensure that only one of the eigenvectors flips sign over the close curve. For the larger circle (Fig. 7e), the essential effect of the sign flips is to ensure that over the closed curve there are no cumulative flips. Multiple flips, in these case, may be due to complex topology of the branch cut or other branch points (exceptional points) in the vicinity of the exceptional point considered. That there are more than one flips in the eigenvectors around an exceptional point does not necessarily defeat their use in our proposed band-sorting algorithm. What is important is that the flips do not happen at all when the closed curve does not encircle an exceptional point. This is shown in Figs. (7g-l). Figs. (7g,j) show four closed curves around two real cross points in the phononic bandstructure (T and L2 respectively in Fig. 2b). Figs. (7h,i,k,l) show the behavior of the inner products as these curves are traversed. There are no flips of sign as there are no exceptional points in the vicinity of real crosses. Numerical calculations show that the two eigenvalues in the case of point T coalesce on the real axis at $Q_{1}=0.37000376$, $Q_{2}=0$. For point L2 they again coalesce on the real axis at $Q_{1}=Q_{2}=0.41390153$.

IV Band Sorting for Phononic Crystals

Due to the presence of level repulsion, any band-sorting algorithm which depends solely upon eigenvector or modeshape continuation calculated on coarse wavenumber discretization is doomed from the start. This is because of the fact that in the vicinity of level repulsion not only do the eigenvalues avoid crossing each other but they also exchange their modeshapes. This ensures that any band-sorting algorithm which depends upon eigenvector continuation will mistake a level repulsion point as a real cross point if the modes of interest belong to two neighboring points on a rough grid of wavevectors. This phenomenon is explicitly shown in Fig. (8). Figs. (8a-d) are the displacement fields in the $x_{1}$ and $x_{2}$ directions respectively at points $t_{1}-t_{4}$ from Fig. (3a). The similarity of the modeshapes suggests that $t_{4}$ should be the continuation of $t_{1}$ and $t_{3}$ should be the continuation of $t_{2}$. This conclusion is correct since T is a real cross point. Figs. (8e-h) are the displacement fields in the $x_{1}$ and $x_{2}$ directions at points $r_{1}-r_{4}$ from Fig. (3b). In this case, the modeshape similarity suggests that $r_{4}$ should be the continuation of $r_{1}$ and $r_{3}$ should be the continuation of $r_{2}$. However, in this case the conclusion would be incorrect since point R is a level repulsion point as revealed by the fine discretization in Fig. (3d). Due to the exchange of eigenvectors, therefore, it is impossible to distinguish a level repulsion point from a real cross point purely from eigenvector orthogonality considerations.

Here we propose a more efficient method of sorting the band-structure which forgoes the need for the zoom-in method. The idea rests on the observation from Fig. (7) that flips of signs in eigenvectors should be present as one considers the eigenvalue problem at discrete locations on a closed curve which encircles an exceptional point. No such flips would occur if the closed curve does not encircle an exceptional point. The process of implementation of this idea is further explained in Fig. (9) with respect to the two points discussed in Fig. (3). Figs. (9a,b) correspond to Figs. (3a,b) respectively. We consider the location of the hypothetical cross point (real or apparent) between two curves of interest and as determined by the coarse discretization. This is determined by finding the intersection of the two lines which approximate the curves in the region of interestHeiss and Sannino (1990, 1991). We draw a circle in the complex domain of the parameter. Its lower extremity lies at the hypothetical location of the cross-point and its radius is large enough to encompass the exceptional point if there is one in the vicinity. In general, the distance of the exceptional point from the real axis of the parameter is small and dependent upon the distance between the two curves in the region of repulsionHeiss and Sannino (1990). In the case that level repulsion exists, such a circle will most likely encircle the exceptional point and exclude its complex conjugate exceptional point. The situation is shown in Fig. (9b) where the circle encircles the associated exceptional point. Fig. (9a) shows that the circle encircles no exceptional point since none exists in this case. The eigenvectors are now evaluated at some discrete points on the circle (marked in Fig. 9a,b). Since the radius of the circle is still very small compared to the domain of the parameter (it is smaller than the coarse discretization), the set of two relevant eigenvectors at the discrete points will be nearly the same unless there is a flip of sign. In other words, the inner products of the eigenvectors corresponding to different points on the closed curve would be either very close to +1, 0, or -1. In the case of a real cross, all mutual inner products will be either close to +1 or 0. In the case of level repulsion there will be some inner products which have a value close to -1.

To implement the above ideas into a band-sorting algorithm, the following processes are implemented: As a first step, we sort the eigenvalues in ascending order and rearrange the eigenvectors accordingly. Tentative mode continuation is determined by checking approximate orthogonality (13), at adjacent points on the coarse discretization. This process results, for example, in the band-structure shown in Fig. (10a) where individual colors indicate individual phononic bands. At this point the cross locations have been determined but level repulsion zones have not yet been distinguished from real crosses. We now form closed curves in the complex parametric domain near the vicinity of the crosses and analogous to the process shown in Fig. (9). By checking for flips of signs at discrete points over these closed curves and in accordance with the process described above, we distinguish real crosses from level repulsion zones. The results are shown in Fig. (10b) where the level repulsion zones have been taken into account. Comparing it to Fig. (10a), we can see that several mode continuations have been changed in the process. Band continuation of the first 6 curves agrees exactly with results in literatureWu and Huang (2004). For the above, the bandstrucuture is calculated by evaluating the eigenvalue problem at 256 wavevector points along the boundary of IBZ. To distinguish level repulsion zones from real crosses, circles of diameters $1\times 10^{-3}$ are constructed in the complex domain of the parameter. Flips of signs are then checked at 9 discrete points along each circle (Fig. 9).

In order to evaluate the overall performance of this method, we now present band-sorting results for larger number of phononic bands. The fully sorted bandstructutre for 683 phononic bands is shown in Fig. (11a). Sorting has been carried out using the zoom-in method here so there is no ambiguity in the results. There is a total of 1764 crosses (real or apparent) in the band-structure. Since the zoom-in method requires the evaluation of the eigenvalue problem at fine discretizations near the cross points, we had to evaluate the problem a total of 137592 times to achieve the sorting. Real cross points are marked by red markers and the level repulsion points are marked by green ones. It is interesting to note that all points between $X-M$, which would be determined as cross points using solely eigenvector continuity, turn out to be level repulsion points. The majority of points between $\Gamma-X$ and $M-\Gamma$ are real cross points. Fig. (11b) shows the sorted bandstructure based on examining the existence of sign flips. The latter process correctly identifies the crosses as real or apparent for 97.39% of the total 1764 crosses. In effect, only 46 out of the 1764 total points are incorrectly identified by the flip of sign method. This essentially means, for the current problem, that 92.69% of the 683 phononic bands have been assigned correct mode continuation through the flip of sign method. In contrast, if the bands were sorted in ascending order of magnitude, the accuracy of mode continuation would have been just 38.55%. If the bands are sorted based only on the inner product continuity on the coarse grid then the accuracy would have been 79.08%. The current implementation of the flip of sign method is roughly an order of magnitude faster than the zoom in method. Both the accuracy and the speed of the method can be improved as we discuss in the next section.

V Discussion and Conclusions

In this paper we have considered the problem of level repulsion in the context of phononics and suggested a computationally efficient strategy to distinguish them from normal cross points. This process is essential for the correct sorting of the phononic bands and, subsequently, for the accurate determination of mode continuation, group velocities, and emergent properties which depend on them. We have shown that in the vicinity of the exceptional point the relevant phononic eigenvalue surfaces resemble the surfaces of a 2 by 2 matrix. Along a closed loop encircling the exceptional point the phononic eigenvalues are exchanged, just as they are for the 2 by 2 matrix case. However, the behavior of the associated eigenvectors is shown to be more complex in the phononic case. Along a closed loop around an exceptional point, we show that the eigenvectors can flip signs multiple times unlike a 2 by 2 matrix where the flip of sign occurs only once. Based upon this sign-flip property near an exceptional point, we have proposed a band-sorting technique for the accurate sorting of phononic eigenvalues. The method is roughly an order of magnitude faster than the fail-safe zoom-in method. As implemented, the method correctly identified $>$97% of the total cross points in the example considered which translated into the correct mode continuation assignment for $>$90% of the phononic branches.

There are several ways of further increasing both the accuracy and the speed of the sorting method. First, it should be noted that the only reason that the flip of sign method, as currently implemented, incorrectly identifies a level repulsion point as a real cross point is that the considered closed curve (circle in our case) fails to encircle an exceptional point. The circle misses encircling an exceptional point if the exceptional point is too close to the real axis and if the hypothetical location of the intersection of the two curves (which is also the lower extremity of the considered circle) is such that the exceptional points sits just below the circle boundary. This issue can potentially be remedied by considering other closed curves such as a semicircle instead of a circle. It may even be possible to forgo the closed curve completely and just evaluate the inner products at randomly chosen locations in the complex domain of the parameter and near the hypothetical cross point of the two curves under consideration.

The current method never identifies a real cross point as a level repulsion point unless there is a real cross point in close vicinity of a level repulsion point. For example, Fig. (11c) shows that a real cross point interferes with level repulsion in the region marked by a small green square in Fig. (11b). The sorting algorithm draws two circles of the same diameter $1\times 10^{-3}$ in the complex domain with one located at the real cross point and the other at the level repulsion point. These circles are large enough such that the one at the real cross point ends up encircling the exceptional point which belongs to the level repulsion zone. Figs. (11d,e)show that both eigenvalue trajectories along the two curves pass the branch cut where the eigenvectors flip signs. Therefore, this real cross point is mistakenly identified as a level repulsion by the algorithm. Such a situation can be remedied by employing a smaller closed curve. However, it must be noted that a smaller curve will, in general, increase the chance that it may fail to encircle an exceptional point when it exists. Therefore, further strategies should take into account this situation. Finally, we have chosen to evaluate 9 discrete points along the closed curves to check for sign-flips. It is entirely possible that one can implement an algorithm which maintains the accuracy while evaluating the eigenvalue problem at far fewer locations. This would further speed-up the process of band-sorting.

Acknowledgment

A.S. acknowledges support from the NSF CAREER grant # 1554033 to the Illinois Institute of Technology.