Complex-Eigenfrequency Band Structure of Viscoelastic Phononic Crystals

Featured Application

The method in this paper applies to the prediction of temporal damping in phononic crystal structures composed of materials that can be described by a viscoelastic model such that the imaginary part of elastic constants is proportional to frequency.

Abstract

The consideration of material losses in phononic crystals leads naturally to the introduction of complex valued eigenwavevectors or eigenfrequencies representing the attenuation of elastic waves in space or in time, respectively. Here, we propose a new technique to obtain phononic band structures with complex eigenfrequencies but real wavevectors, in the case of viscoelastic materials, whenever elastic losses are proportional to frequency. Complex-eigenfrequency band structures are obtained for a sonic crystal in air, and steel/epoxy and silicon/void phononic crystals, with realistic viscous losses taken into account. It is further found that the imaginary part of eigenfrequencies are well predicted by perturbation theory and are mostly independent of periodicity, i.e., they do not account for propagation losses but for temporal damping of Bloch waves.

1. Introduction

Phononic crystals for elastic waves are classical analogs to crystal lattices for phonons [1,2]. They are described with continuum mechanics and can also be viewed as periodic composites. At the frequencies involved, that go from audible sound (a few Hz), to ultrasound (up to a few tens of GHz), it is usual to neglect material losses in a first approximation. In any solid material, however, material loss is present and ultimately increases with frequency. The Kelvin–Voigt viscoelastic model for isotropic solid materials, that turns Young’s modulus E into a complex and dispersive number $E+ı\omega \eta$ with $\omega$ being the angular frequency and $\eta$ the viscosity, applies well at ultrasonic frequencies. In crystals such as silicon or quartz, material loss can often be modeled by adding an imaginary part to the elastic tensor, in a first approximation proportional to the frequency [3,4]. As a result, the elastic tensor $c_{ijkl}$ becomes the complex valued tensor $c_{ijkl}+ı\omega {\eta }_{ijkl}$, with ${\eta }_{ijkl}$ the phonon viscosity tensor. Sonic crystals for acoustic waves also contain viscous fluids that can be described by a dynamical viscosity leading to a bulk modulus whose imaginary part increases proportionally to frequency. This applies to common fluids supporting sound propagation, such as air and water.

In the literature on phononic crystals, having a reliable measure of the effect of material loss on wave propagation has been a standing salient problem, with obvious practical implications [5,6,7,8,9,10,11,12,13,14]. Since the band structure gives the dispersion relation of Bloch waves, i.e., of the eigenmodes of phononic crystals, it is natural to try to generalize it to complex quantities describing attenuation in space or damping in time. The complex band structure gives the complex-valued Bloch wavevector as a function of real frequency [7,15,16]. It is in essence well suited to the description of the attenuation in space of monochromatic waves originating from a source of finite extent excited at a given frequency. There are other situations where the wavevector can be imposed, such as with laser generated ultrasound, but waves will be damped as time goes by after the excitation process. Enhanced energy damping in metamaterials has also been proposed [17,18]. It is the purpose of this paper to propose a complex-eigenfrequency band structure method that is suited for viscoelastic materials. In the following, an augmented complex eigenvalue problem is defined using an auxilliary field technique, giving direct damping information for every Bloch wave in the band structure. Examples are given for a sonic crystal of rigid scatterers in air, for steel-epoxy phononic crystals, and for a holey silicon phononic crystal. It is found that numerical results follow closely the perturbation theory in Ref. [11], according to which the damping of a given Bloch wave can be estimated at any frequency from a volume average of the viscoelastic constituents weighted by the Bloch wave distribution. In particular, the complex-eigenfrequency band structure $\omega \left(k\right)$ does not capture any spatial propagation information but provides one with information complementary to the complex band structure $k\left(\omega \right)$.

2. Theory

2.1. Complex Eigenvalue Problem

We consider a primitive unit-cell of an artificial crystal, such as those depicted in Figure 1, Figure 2 and Figure 3. Elastic Bloch waves have a displacement vector of the form $\mathbf{U}(\mathbf{r},t)=\mathbf{u}\left(\mathbf{r}\right)exp\left(ı\right(\omega t-\mathbf{k}·\mathbf{r}\left)\right)$ where $\mathbf{u}\left(\mathbf{r}\right)$ is the periodic part of the Bloch wave, $\omega$ is the angular frequency, t is the time variable, and $\mathbf{k}$ is the Bloch wavevector. In elastic solids, the stress and strain tensors are related by Hooke’s law via the order-four elastic tensor $c_{ijkl}$. Bloch waves of phononic crystals are obtained by solving an eigenfrequency problem under periodic boundary conditions. In viscoelastic solids, Hooke’s law is generalized to a complex-valued elastic tensor $c_{ijkl}+ı\omega {\eta }_{ijkl}$, with ${\eta }_{ijkl}$ the order-four phonon viscosity tensor. The imaginary part is thus explicitly proportional to the angular frequency. As a result, the eigenfrequency problem becomes non trivial, since the matrix coefficients depend on frequency. A complex eigenfrequency band structure can be obtained anyway using the method described in the following and with the assumption that viscoelastic losses are proportional to $\omega$. We start with the formal eigenproblem with a viscoelastic term

$$\begin{array}{c}\hfill (\mathbf{K}\left(\mathbf{k}\right)+ı\omega \mathbf{V}\left(\mathbf{k}\right)-{\omega }^2\mathbf{M})\mathbf{u}=0,\end{array}$$

(1)

with $\mathbf{u}$ a vector of the degrees of freedom (d.o.f), $\mathbf{K}$ a stiffness matrix, $\mathbf{M}$ a mass matrix, and $\mathbf{V}$ a viscosity matrix. Note that both the stiffness matrix and the viscosity matrix are a function of the Bloch wavevector $\mathbf{k}$; we omit the dependence on $\mathbf{k}$ in the rest of this sub-section for simplicity of the presentation. The mass matrix generally has constant coefficients. All matrices are square and the number of lines equals the number of d.o.f of the system. We set $\lambda =ı\omega$ and rewrite Equation (1)

$$\begin{array}{c}\hfill (\mathbf{K}+\lambda \mathbf{V}+{\lambda }^2\mathbf{M})\mathbf{u}=0.\end{array}$$

(2)

This second-degree polynomial eigen-equation is equivalent to the first-degree system of equations

$$\begin{array}{c}\hfill \mathbf{v}=\lambda \mathbf{u},\end{array}$$

(3)

$$\begin{array}{c}\hfill \mathbf{K}\mathbf{u}+\lambda \mathbf{V}\mathbf{u}+\lambda \mathbf{M}\mathbf{v}=0,\end{array}$$

(4)

or finally equivalent to the double-size eigenvalue problem

$$\begin{array}{c}\hfill \left(\begin{array}{cc}\mathbf{K}& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}\mathbf{u}\\ \mathbf{v}\end{array}\right)=\lambda \left(\begin{array}{cc}-\mathbf{V}& -\mathbf{M}\\ 1& 0\end{array}\right)\left(\begin{array}{c}\mathbf{u}\\ \mathbf{v}\end{array}\right).\end{array}$$

(5)

This asymmetric eigenvalue problem yields complex eigenvalues, and hence the complex-eigenfrequency band structure. The wavevector $\mathbf{k}$ enters via the dependence of matrices $\mathbf{K}$ and $\mathbf{V}$. Note that the eigensystem can be written in an Hermitian symmetric form, providing the mass matrix can be inverted and is symmetric, and both K and V are Hermitian symmetric, e.g.,

$$\begin{array}{c}\hfill \left(\begin{array}{cc}\mathbf{K}& 0\\ 0& -{\mathbf{M}}^{-1}\end{array}\right)\left(\begin{array}{c}\mathbf{u}\\ -\mathbf{M}\mathbf{v}\end{array}\right)=\lambda \left(\begin{array}{cc}-\mathbf{V}& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}\mathbf{u}\\ -\mathbf{M}\mathbf{v}\end{array}\right)\end{array}$$

(6)

As a result, eigenvalues come in complex conjugate pairs $(\lambda ,{\lambda }^{*})$. Each pair of eigenfrequencies have the same real part (and hence propagation constant) but the opposite imaginary part. The eigenfrequency with positive imaginary part is a mode of vibration that attenuates in time, whereas the eigenfrequency with the negative imaginary part is a mode of vibration that amplifies in time. Note that their respective excitation is dictated by initial boundary conditions. For instance, specifying that the energy in the mode cannot grow to infinity with increasing time disqualifies the amplifying mode of vibration.

2.2. Phononic Crystal Containing Viscoelastic Materials

We now specify how to obtain an equation of the form (1) for a phononic crystal composed of viscoelastic materials. For elastic waves in solids, the elastodynamic equation for the displacement vector components $U_i$ can be written as the partial differential equation [2]

$$\begin{array}{c}\hfill -{\left[(c_{ijkl}\left(\mathbf{r}\right)+ı\omega {\eta }_{ijkl}\left(\mathbf{r}\right))U_{k,l}\left(\mathbf{r}\right)\right]}_{,j}={\omega }^2\rho \left(\mathbf{r}\right)U_i\left(\mathbf{r}\right).\end{array}$$

(7)

A variational formulation of the problem of Bloch wave propagation allows one to write the eigenvalue problem in the form (1). Note that we consider the periodic part of the elastic Bloch waves as the variables for the variational formulation. A weak form suitable for finite element implementation can be obtained directly by considering a mixed finite element space with variables $(\mathbf{u},\mathbf{v})$ and test functions $({\mathbf{u}}^{\prime },{\mathbf{v}}^{\prime })$ and reads (where the dependence of functions on the space coordinates is implicit)

$$\begin{array}{c}\hfill {\int }_{\Omega }{\left(\begin{array}{c}{\mathbf{u}}^{\prime }\\ {\mathbf{v}}^{\prime }\end{array}\right)}^T\left(\begin{array}{cc}(\nabla +ı\mathbf{k})c(\nabla -ı\mathbf{k})& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}\mathbf{u}\\ \mathbf{v}\end{array}\right)=\lambda {\int }_{\Omega }{\left(\begin{array}{c}{\mathbf{u}}^{\prime }\\ {\mathbf{v}}^{\prime }\end{array}\right)}^T\left(\begin{array}{cc}-(\nabla +ı\mathbf{k})\eta (\nabla -ı\mathbf{k})& -\rho \\ 1& 0\end{array}\right)\left(\begin{array}{c}\mathbf{u}\\ \mathbf{v}\end{array}\right)\end{array}$$

(8)

under periodic boundary conditions. This is the formulation we use for implementation of the complex-eigenfrequency band structures in Section 3.2 and Section 3.3. The different matrices thus have the following formal expressions

$$\begin{array}{cc}\hfill K\left(\mathbf{k}\right)& ={\int }_{\Omega }(u_{i,j}^{\prime }+ık_j{u}_i^{\prime })c_{ijkl}(u_{k,l}-ık_l{u}_k),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill V\left(\mathbf{k}\right)& ={\int }_{\Omega }(u_{i,j}^{\prime }+ık_j{u}_i^{\prime }){\eta }_{ijkl}(u_{k,l}-ık_l{u}_k),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill M& ={\int }_{\Omega }{u}_i^{\prime }\rho v_i.\hfill \end{array}$$

(11)

2.3. Extension to Sonic Crystals

Acoustic Bloch waves have a pressure field of the form $P(\mathbf{r},t)=p\left(\mathbf{r}\right)exp\left(ı\right(\omega t-\mathbf{k}·\mathbf{r}\left)\right)$ with $p\left(\mathbf{r}\right)$ the periodic part of the pressure field. Wave propagation in a sonic crystal can be written as the partial differential equation

$$\begin{array}{c}\hfill -{[\rho {\left(\mathbf{r}\right)}^{-1}{P}_{,j}\left(\mathbf{r}\right)]}_{,j}={\omega }^2{(B\left(\mathbf{r}\right)+ı\omega \eta \left(\mathbf{r}\right))}^{-1}P\left(\mathbf{r}\right).\end{array}$$

(12)

Here the elastic modulus B and the viscosity $\eta$ are scalar functions of position. Generally speaking, we have a problem to obtain an equation of the form (1) since the kinetic part of the equation involves the inverse function ${(B\left(\mathbf{r}\right)+ı\omega \eta \left(\mathbf{r}\right))}^{-1}$. Restricting the problem to a single homogeneous viscous fluid, however, the material constants are independent of space coordinates and (12) can be multiplied by $1+ı\omega \eta /B$ to get

$$\begin{array}{c}\hfill -{\left[\frac{1}{\rho }{P}_{,j}\left(\mathbf{r}\right)\right]}_{,j}-{\left[\frac{ı\omega \eta }{\rho B}{P}_{,j}\left(\mathbf{r}\right)\right]}_{,j}={\omega }^2\frac{1}{B}P\left(\mathbf{r}\right).\end{array}$$

(13)

As a result, we can introduce $q=ı\omega p$ and the test functions $(p^{\prime },q^{\prime })$ to obtain the following variational formulation

$${\int }_{\Omega }{\left(\begin{array}{c}{p}^{\prime }\\ q^{\prime }\end{array}\right)}^T\left(\begin{array}{cc}(\nabla +ı\mathbf{k})\frac{1}{\rho }(\nabla -ı\mathbf{k})& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}p\\ q\end{array}\right)=\lambda {\int }_{\Omega }{\left(\begin{array}{c}{p}^{\prime }\\ q^{\prime }\end{array}\right)}^T\left(\begin{array}{cc}-(\nabla +ı\mathbf{k})\frac{\eta }{\rho B}(\nabla -ı\mathbf{k})& -\frac{1}{B}\\ 1& 0\end{array}\right)\left(\begin{array}{c}p\\ q\end{array}\right)$$

(14)

under periodic boundary conditions. The different matrices thus have the following formal expressions

$$\begin{array}{cc}\hfill K\left(\mathbf{k}\right)& ={\int }_{\Omega }(\nabla q+ı\mathbf{k}q)\frac{1}{\rho }(\nabla p-ı\mathbf{k}p),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill V\left(\mathbf{k}\right)& ={\int }_{\Omega }(\nabla q+ı\mathbf{k}q)\frac{\eta }{\rho B}(\nabla p-ı\mathbf{k}p),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill M& ={\int }_{\Omega }{p}^{\prime }\frac{1}{B}q.\hfill \end{array}$$

(17)

3. Results

The variational formulations of the previous section were implemented with the finite element method with language FreeFem$++$ [19]. Obviously, this choice is not unique and other finite element codes could be employed for the same purpose. In the following of this section, we consider representative two-dimensional artificial crystals. Their unit cells are represented on a finite element mesh enclosed by boundaries. Periodic boundary conditions are applied on the external boundaries. Lagrange finite elements of degree two are used for the approximation of all unknown and test functions.

3.1. Sonic Crystal of Rigid Rods in Air

As a first example, we consider a two-dimensional sonic crystal of rigid rods in air, as was for instance considered experimentally by Miyashita for cylindrical rods [20,21] and by Lu et al. for triangular rods [22]. The solid rods are supposed to be sufficiently heavy and stiff that only negligible elastic waves can be excited inside them for sound waves incident from air. As a result, the unit-cell of the crystal depicted in Figure 1a (square lattice) or Figure 1c (hexagonal lattice) is composed of a hollow air square or hexagon. With only a single fluid constituent, the theory of Section 2.3 can be applied. The viscosity of air is conventionally $0.0183$ mPa.s at a temperature of $20{}^{\circ }$C, as reported in

Table 1. Independent material constants used in this work.

Air	$\rho =1.2041$$\mathrm{kg}$/${\mathrm{m}}^3$
	$B=142$$\mathrm{k}$$\mathrm{Pa}$
	$\eta =0.0183$$\mathrm{m}$$\mathrm{Pa}$.$\mathrm{s}$
Epoxy	$\rho =7780$$\mathrm{k}\mathrm{g}$/${\mathrm{m}}^3$
(isotropic)	$c_{11}=7.54$$\mathrm{G}$$\mathrm{Pa}$	$c_{66}=1.48$$\mathrm{G}$$\mathrm{Pa}$
	${\eta }_{11}=1$$\mathrm{Pa}$.$\mathrm{s}$	${\eta }_{66}=0.25$$\mathrm{Pa}$.$\mathrm{s}$
Steel	$\rho =1142$$\mathrm{k}\mathrm{g}$/${\mathrm{m}}^3$
(isotropic)	$c_{11}=264$$\mathrm{G}$$\mathrm{Pa}$	$c_{66}=84$$\mathrm{G}$$\mathrm{Pa}$
Silicon	$\rho =2331$$\mathrm{k}\mathrm{g}$/${\mathrm{m}}^3$
(cubic)	$c_{11}=165.7$$\mathrm{G}$$\mathrm{Pa}$	$c_{66}=79.62$$\mathrm{G}$$\mathrm{Pa}$	$c_{12}=63.9$$\mathrm{G}$$\mathrm{Pa}$
	${\eta }_{11}=1.505$$\mathrm{m}$$\mathrm{Pa}$.$\mathrm{s}$	${\eta }_{66}=0.553$$\mathrm{m}$$\mathrm{Pa}$.$\mathrm{s}$	${\eta }_{12}=-0.532$$\mathrm{m}$$\mathrm{Pa}$.$\mathrm{s}$

. Combined with an elastic modulus $B=142$ kPa, one would expect theoretically for longitudinal plane waves a quality factor $Q=\frac{B}{\omega \eta }$ = 12,350 at a frequency of 100 kHz [3]. For the m-th band of the crystal, the quality factor is estimated from the complex eigenfrequencies as

$$\begin{array}{cc}\hfill Q_m\left(\mathbf{k}\right)& =\frac{\mathrm{Re}\left[{\omega }_m\left(\mathbf{k}\right)\right]}{2\mathrm{Im}\left[{\omega }_m\left(\mathbf{k}\right)\right]}.\hfill \end{array}$$

(18)

Figure 1 shows the complex-eigenfrequency band structures in the square and the hexagonal lattice cases, for lattice constant $a=1$ mm. Also plotted are the empty-lattice models obtained when there is no rod in the unit-cell. Though the band structures are very different if the rod is present or not, the quality factor appears to be the same in both case. Actually, Q is simply inversely proportional to frequency in this single-fluid case (see the discussion of Section 4).

3.2. Epoxy–Steel Phononic Crystal

The width of the Bragg band gaps of solid–solid phononic crystals is known to be favored by a high contrast between the elastic constants and the mass density of the contituent solids, in case of hard and heavy inclusions in a soft and light matrix. Phononic crystals of steel rods or spheres in an epoxy matrix are an archetypal example, in both 2D [23] and 3D [24,25]. With lattice constants of the order of millimeters, attenuation was not observed to be problematic for experiments, though epoxy materials are presumably viscoelastic. The viscosity of epoxy materials is quite variable, depending on the composition and on the processes used. We consider in Table 1 rather arbitrary values for the viscosity tensor of epoxy, for illustration purposes. Here, the viscosity of steel is neglected with respect to that of epoxy. For the 2D hexagonal lattice phononic crystal of cylindrical steel rods in an epoxy matrix depicted in Figure 2a, the complex-eigenfrequency band structure shown in Figure 2b as a large Bragg band gap [23] and the quality factor Q decreases again hyperbolically with increasing frequency.

If the roles of steel and epoxy are reversed, for the case of epoxy rods in a steel matrix depicted in Figure 2c, then it is known that Bragg band gaps are not favored anymore but that local resonances are induced [26,27,28]. In the complex-eigenfrequency band structure shown in Figure 2d, such a local resonance appears for instance around a frequency of $1.25$ $\mathrm{M}$$\mathrm{Hz}$. It can be observed that the flatter bands have a much smaller quality factor compared to steeper bands. Actually, and despite appearances, this phenomenon is not related with band steepness but with the spatial distribution of Bloch waves within the unit-cell, as discussed in Section 4.

3.3. Phononic Crystal of Holes in Silicon

Another system that has been often considered for microsonic applications is the case of holes in a crystal matrix. In the case of silicon phononic crystals [29,30], cylindrical holes can be obtained by etching at the microscale. The phonon viscosity tensor for silicon has been discussed e.g., in Ref. [31], from which the values in Table 1 are taken. Alternative values are given by Helme and King [32]. For a 2D square lattice phononic crystal of holes in silicon depicted in Figure 3a, the complex-eigenfrequency band structure shown in Figure 3b has a clear Bragg band gap, whereas in the case of the 2D hexagonal lattice depicted in Figure 3c the band gap in Figure 3d closes. More significantly, if it is again observed that the quality factor varies as the inverse of frequency, there are different damping rates for the different bands. As discussed in Section 4, this variability is anisotropic in nature and the quality factor depends on the vector polarization of Bloch waves.

4. Discussion

In all numerical results presented in Section 3, the imaginary parts of the viscoelastic constants were much smaller than their real parts. These results are thus mainly relevant to the case where one is interested in engineering wave propagation properties through a sonic or phononic crystal with small damping. Of course, the equations defining the complex-eigenfrequency band structure also apply beyond this small attenuation regime, as could be interesting for applications to the mitigation of sound or vibrations. In the small damping limit, the first-order perturbation theory of Ref. [11] applies. The main result of this theory in our context is that the relative variation of the frequency of a band at a fixed value of the Bloch wavevector is given by

$$\begin{array}{cc}\hfill \frac{{\left.\delta \omega \right|}_{\mathbf{k}}}{\omega }& =\frac{ı}{2}\omega \frac{\langle {\nabla }_{\mathbf{k}}\mathbf{u}\left|\eta \right|{\nabla }_{\mathbf{k}}\mathbf{u}\rangle }{\langle {\nabla }_{\mathbf{k}}\mathbf{u}\left|c\right|{\nabla }_{\mathbf{k}}\mathbf{u}\rangle },\hfill \end{array}$$

(19)

where

$$\begin{array}{cc}\hfill \langle {\nabla }_{\mathbf{k}}\mathbf{u}\left|c\right|{\nabla }_{\mathbf{k}}\mathbf{u}\rangle & ={\int }_{\Omega }(u_{i,j}^{*}+ık_j{u}_i^{*})c_{ijkl}(u_{k,l}-ık_l{u}_k),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \langle {\nabla }_{\mathbf{k}}\mathbf{u}\left|\eta \right|{\nabla }_{\mathbf{k}}\mathbf{u}\rangle & ={\int }_{\Omega }(u_{i,j}^{*}+ık_j{u}_i^{*}){\eta }_{ijkl}(u_{k,l}-ık_l{u}_k).\hfill \end{array}$$

(21)

The latter equations are the spatial averages of the real and imaginary part of the viscoelastic tensor taken with respect to the particular Bloch wave considered. They include all the influence of the details of the crystal on the appearance of an imaginary part for the eigenfrequency. The quality factor for band m can then be estimed as

$$\begin{array}{cc}\hfill Q_m\left(\mathbf{k}\right)& =\frac{\langle {\nabla }_{\mathbf{k}}\mathbf{u}\left|c\right|{\nabla }_{\mathbf{k}}\mathbf{u}\rangle }{{\omega }_m\left(\mathbf{k}\right)\langle {\nabla }_{\mathbf{k}}\mathbf{u}\left|\eta \right|{\nabla }_{\mathbf{k}}\mathbf{u}\rangle }.\hfill \end{array}$$

(22)

A practical algorithm, and an alternative to the theory of Section 3 in the limit of small viscoelasticity, is to solve for the lossless band structure first, thus obtaining the real eigenvalues ${\omega }_m\left(\mathbf{k}\right)$ and their corresponding Bloch waves. Equation (22) then gives a first-order perturbation theory estimate of the quality factor for any band.

When there is a single homogeneous material in a sonic crystal, e.g., as is the case for the air crystal with rigid inclusions of Figure 1, the relative frequency variation simplifies to

$$\begin{array}{cc}\hfill \frac{{\left.\delta \omega \right|}_{\mathbf{k}}}{\omega }& =\frac{ı}{2}\omega \frac{\eta }{B},\hfill \end{array}$$

(23)

i.e., the distribution of the Bloch wave has no incidence on the result and the quality factor has the same expression as for a plane wave in an homogeneous medium, $Q=\frac{B}{\omega \eta }$. This is the reason why the sonic crystal and the empty lattice model show the same variation of the quality factor with frequency in Figure 1. Note that the quality factor is also independent of the filling fraction of the crystal in this scalar homogeneous case. When there is a single homogeneous but anisotropic material in a phononic crystal, as is the case for the holey silicon crystal of Figure 3, similar conclusions can be drawn but the quality factor is not independent of the filling fraction and it further depends on the direction of propagation and on the polarization of the particular Bloch wave, i.e., the ratio of the material averages $\langle C\rangle /\langle \eta \rangle$ is different for quasi longitudinal and quasi shear Bloch waves. As a consequence, considering the empty lattice model in this case would not be meaningful.

In the case of a phononic crystal with several constituents, as in the case of steel-epoxy crystals in Figure 2, the material averages can be decomposed over matrix and inclusion regions according to

$$\begin{array}{cc}\hfill \frac{{\left.\delta \omega \right|}_{\mathbf{k}}}{\omega }& =\frac{ı}{2}\omega \frac{\langle {\nabla }_{\mathbf{k}}\mathbf{u}\left|\eta \right|{\nabla }_{\mathbf{k}}{\mathbf{u}\rangle }_{\mathrm{matrix}}+\langle {\nabla }_{\mathbf{k}}\mathbf{u}\left|\eta \right|{\nabla }_{\mathbf{k}}{\mathbf{u}\rangle }_{\mathrm{incl}.}}{\langle {\nabla }_{\mathbf{k}}\mathbf{u}\left|c\right|{\nabla }_{\mathbf{k}}{\mathbf{u}\rangle }_{\mathrm{matrix}}+\langle {\nabla }_{\mathbf{k}}\mathbf{u}\left|c\right|{\nabla }_{\mathbf{k}}{\mathbf{u}\rangle }_{\mathrm{incl}.}}.\hfill \end{array}$$

(24)

As a result, damping further depends on the spatial distribution of the Bloch wave. For instance, in the case of a local resonance that results from the hybridization of a resonance of an inclusion (i.e., a flat band) with propagating waves in the matrix (i.e., a band with a given slope), the Bloch waves of the flat band are generally localized inside the inclusion and suffer larger damping if the viscosity of the inclusion is larger than the viscosity of the matrix. That is the situation observed for the local resonances in Figure 2d.

We note that in comparison with the complex band structure that can be used for an arbitrary dependence of the material constants with frequency [7,10,12], the complex eigenfrequency band structure introduced in the present paper does not capture any effect linked with the spatial attenuation of the Bloch waves as they propagate away from a source. Hence, the results we have obtained are independent of the local group velocity and do not reveal anything regarding frustrated evanescent Bloch waves or evanescent Bloch waves inside band gaps [7]. They essentially account for temporal damping of Bloch waves in crystals and include the contribution of local resonances.

There are further situations where the complex-eigenfrequency band structure could be useful and we indicate some of them as perspectives. In the case of a defect cavity embedded in a very large crystal, the method would give the quality factor of the resonance, if combined with a super-cell technique. Beyond time-harmonic excitations, transient or short pulse excitations are equally important for experiments. Since damped eigenmodes provide a complete basis for analysis of the frequency response function (FRF) of structures, the general methodology of the linear combination of damped eigenmodes can be combined with the complex-eigenfrequency band structure.