Dispersion of waves in two and three-dimensional periodic media

TL;DR

This paper develops a new asymptotic method to analyze wave dispersion in 2D and 3D periodic media with small inclusions, providing explicit formulas for effective wave speed and dispersion relations.

Contribution

A novel approach to derive complete asymptotic expansions of dispersion relations for waves in periodic media with small inclusions, applicable to various boundary conditions.

Findings

Explicit asymptotic formulas for dispersion relations.

Effective wave speed as a function of frequency and material properties.

Illustrations of parameter dependence through graphs.

Abstract

We consider the propagation of acoustic time-harmonic waves in a homogeneous media containing periodic lattices of spherical or cylindrical inclusions. It is assumed that the wavelength has the order of the periods of the lattice while the radius $a$ of inclusions is small. A new approach is suggested to derive the complete asymptotic expansions of the dispersion relations in two and three-dimensional cases as $a \to 0$ and evaluate explicitly several first terms. Our method is based on the reduction of the original singularly perturbed (by inclusions) problem to the regular one. The Neumann, Dirichlet and transmission boundary conditions are considered. The effective wave speed is obtained as a function of the wave frequency, the filling fraction of the inclusions, and the physical properties of the constituents of the mixture. Dependence of asymptotic formulas obtained in the paper on geometric and material parameters is illustrated by graphs.

Equations

$$\nabla \cdot \frac{1}{\varrho (\mathbf{r})}\nabla u+\gamma (\mathbf{r}){\omega }^{2}u=0,\mathbf{r}\notin \partial \mathcal{B},$$

$$\Delta u+k_{\pm }^{2}u=0,\mathbf{r}\notin \partial \mathcal{B},$$

$$k_{-}=$$

$$k_{+}=$$

$$\displaystyle\left.\left\llbracket u(\bm{r})\right\rrbracket\right.$$

$$\displaystyle\left.\left\llbracket\frac{1}{\varrho(\bm{r})}\frac{\partial u(\bm{r})}{\partial n}\right\rrbracket\right.$$

$$u(\mathbf{r})=\Phi (\mathbf{r}){\mathrm{e}}^{-\mathrm{i}\mathbf{k}\cdot \mathbf{r}},$$

$$]]{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}}u(\mathbf{r})[[=0,$$

$$\begin{array}[]{l}\Delta u+k^{2}_{-}u=0,\quad r<a,\\[5.69054pt] \Delta u+k^{2}_{+}u=0,\quad\bm{r}\in\Pi\cap\{r>a\},\end{array}$$

$$\displaystyle\left\llbracket u(\bm{r})\right\rrbracket=0,\quad\left\llbracket\frac{1}{\varrho(\bm{r})}\frac{\partial u(\bm{r})}{\partial n}\right\rrbracket=0,\quad\rrbracket\mathrm{e}^{\mathrm{i}{\bm{k}}\cdot\bm{r}}u(\bm{r})\llbracket=0.$$

$$|\mathbf{k}{|}^2=G(\hat{\mathbf{k}},\omega ,a),\omega \ne 0,$$

$$|\mathbf{k}{|}^2$$

$$|\mathbf{k}{|}^2=\frac{{\omega }^2}{c^2}\left(1+\frac{1}{2}f+\frac{3}{20}{\left(\frac{\omega a}{c}\right)}^{2}f\right)+O\left(a^6\right),f=\frac{4\pi a^3}{3|\Pi |}.$$

$$c_1{a}^3$$

$$c_2{a}^5$$

$$+\frac{9}{5}{\left(\frac{\omega a}{c}\right)}^2\frac{f}{({\varrho }_{+}+2{\varrho }_{-}{)}^2}\left[{\varrho }_{+}^2-{\varrho }_{-}^2-{\varrho }_{-}{\varrho }_{+}\left(1-\frac{{\gamma }_{-}}{{\gamma }_{+}}\frac{{\varrho }_{-}}{{\varrho }_{+}}\right)\right].$$

$$|\mathbf{k}{|}^2=\sum _{i=0}^{\infty }\sum _{j=0}^{\infty }{\sigma }_{ij}(\hat{\mathbf{k}},\omega )a^{2i}{\left(a^2\ln a\right)}^j,a\to 0.$$

$$|\mathbf{k}{|}^2=\frac{{\omega }^2}{c^2}\left(1+c_1{a}^2-c_2{a}^4\ln \frac{\omega a}{c}\right)+g(\hat{\mathbf{k}},\omega ,a),|g|<C(\omega )a^4,$$

$$|\mathbf{k}{|}^2=\frac{{\omega }^2}{c^2}\left(1+f-\frac{3}{2}{\left(\frac{\omega a}{c}\right)}^{2}f\ln \frac{\omega a}{c}\right)+O\left(f^2\right),f=\frac{\pi a^2}{|\Pi |}.$$

$$c_1{a}^2=\left(2\alpha -1+\frac{{\gamma }_{-}}{{\gamma }_{+}}\right)f,c_2{a}^4=\left[\frac{1}{2}{\left(1-\frac{{\gamma }_{-}}{{\gamma }_{+}}\right)}^2+{\alpha }^2\right]{\left(\frac{\omega a}{c}\right)}^{2}f,$$

$$|\mathbf{k}{|}^2=\frac{{\omega }^2}{c^2}-\frac{4\pi a}{|\Pi |}+O\left(a^2\right).$$

$$|\mathbf{k}{|}^2=\frac{{\omega }^2}{c^2}-\frac{2\pi }{|\Pi |\ln \frac{\omega a}{c}}+O\left({\ln }^{-2}\left(\frac{\omega a}{c}\right)\right).$$

$$\left\{\begin{array}[]{rl}(\Delta+k^{2}_{+})u(\bm{r})&=\lambda u(\bm{r}),\\[5.69054pt] \left\rrbracket\mathrm{e}^{\mathrm{i}{\bm{k}}\cdot\bm{r}}u(\bm{r})\right\llbracket&=0\end{array}\right.$$

$$\left\{\begin{array}[]{rl}(\Delta+k^{2}_{\pm})u(\bm{r})&=\lambda u(\bm{r}),\quad\bm{r}\in\Pi,\quad r\gtrless a,\\[5.69054pt] \left\llbracket u(\bm{r})\right\rrbracket&=\displaystyle\left\llbracket\frac{1}{\varrho(\bm{r})}\frac{\partial u(\bm{r})}{\partial n}\right\rrbracket=0,\quad\left\rrbracket\mathrm{e}^{\mathrm{i}{\bm{k}}\cdot\bm{r}}u(\bm{r})\right\llbracket=0\end{array}\right.$$

$$\left\{\begin{array}[]{rl}(\Delta+k^{2}_{\pm})u(\bm{r})&=0,\quad\bm{r}\in\Pi,\quad r\gtrless a,\\[5.69054pt] \left\llbracket u(\bm{r})\right\rrbracket&=\displaystyle\left\llbracket\frac{1}{\varrho(\bm{r})}\frac{\partial u(\bm{r})}{\partial n}\right\rrbracket=0,\quad\left\rrbracket\mathrm{e}^{\mathrm{i}{\bm{k}}\cdot\bm{r}}u(\bm{r})\right\llbracket=0\end{array}\right.$$

$$(N_{\mathbf{k}}^{+}-N_a^{-})\psi =0$$

$$\left(\Delta +k_{+}^2\right)u^{+}(\mathbf{r})$$

$$\left(\Delta +k_{\pm }^2\right)v(\mathbf{r})$$

$${\mathcal{N}}_{\mathbf{k}}^{+},{\mathcal{N}}_a^{-}:H^1(\partial B_R)\to L^2(\partial B_R),{\mathcal{N}}_{\mathbf{k}}^{+}\psi ={\frac{\partial u^{+}}{\partial r}|}_{r=R},{\mathcal{N}}_a^{-}\psi ={\frac{\partial v}{\partial r}|}_{r=R},$$

$${\mathcal{N}}_{\mathbf{k}}^{+}-{\mathcal{N}}_a^{-}:H^1(\partial B_R)\to L^2(\partial B_R)$$

$${\int }_{r=R}\left({\mathcal{N}}_{\mathbf{k}}^{+}\psi \right)\overline{\varphi }\mathrm{d}S={\int }_{r=R}\psi \overline{\left({\mathcal{N}}_{\mathbf{k}}^{+}\varphi \right)}\mathrm{d}S,\psi ,\varphi \in H^1(\partial B_R),$$

$$n^{\pm }\psi ={\frac{\partial u^{\pm }}{\partial r}|}_{r=R}$$

$$\Delta u^{+}=0,R<r<1,{u^{+}|}_{r=R}=\psi ,{u^{+}|}_{r=1}=0,$$

Full text

1 Introduction

Periodic media offer a great deal of possibilities for manipulating wave propagation. These include electromagnetic waves in photonic crystals, where one can create bands and gaps in the wave spectrum, positive or negative group velocity [1], slowing down considerably the speed of light [2, 3, 4], nonreciprocal media [5], the self-collimation effect [6] and more that lead to the development of new devices [7].

Similar phenomena can be observed in phononic crystals [8] for elastic or acoustic waves [9] whose band gap structure is employed in sound filters, transducer design and acoustic mirrors. Research on breaking time-reversal symmetry in wave phenomena is a growing area of interest in the field of phononic crystals and metamaterials aimed at realizing one-way propagation devices which have many potential technological applications [10].

Deriving an explicit dispersion relation for the Floquet-Bloch waves in two and three-dimensional periodic media is an arduous problem and is usually performed numerically [11]. However, assuming that the wavelength is long compared to the period of the lattice or a characteristic size of the scatters one can obtain an asymptotic approximation. Typically it employs the method of matched asymptotic expansions as in [12, 13] for small scatterers with the Dirichlet or Neumann boundary conditions. The latter results were further developed in [14]. A semi-analytical approach using the multipole expansion method is described in [15]. A rigorous analysis of a sub-wavelength plasmonic crystal was presented in [16], where a solution of a nonlinear eigenvalue problem is given in terms of convergent high-contrast power series for the electromagnetic fields and the first branch of the dispersion relation. Explicit formulas for the effective dielectric tensor and the dispersion relation are obtained in [17] assuming that the cell size is small compared to the wavelength, but large compared to the size of the inclusions. Some other approached are presented in [18, 19, 20, 21, 22, 23].

We consider the propagation of waves, governed by the Helmholtz equation in a medium containing periodic lattices of spherical or cylindrical inclusions of radius $a$. We assume transmission boundary conditions on the inclusions’ interface and that $a$ is small relative to the period of the lattice while the wave length is comparable to (or larger than) the lattice size. The Dirichlet and Neumann boundary conditions are also discussed. We suggest a new method for determining the dispersion relations of the Floquet-Bloch waves. The method reduces the original singularly perturbed problem to the regular one and provides explicit formulas for the dispersion relations in two and three dimensional settings with rigorous estimates of the remainders. Our approach can be extended if small inclusions have arbitrary shape. The coefficients of asymptotic expansion in this case will be expressed through solutions of a regular boundary value problem. This will be published elsewhere.

2 Formulation of the problem

We consider the propagation of acoustic waves through an infinite medium containing a periodic array $\mathcal{B}$ of spherical obstacles. The periods of the lattice $\bm{\tau}_{1}$, $\bm{\tau}_{2}$ and $\bm{\tau}_{3}$ are normalized in such a way that $\ell=\min\{|\bm{\tau}_{1}|,|\bm{\tau}_{2}|,|\bm{\tau}_{3}|\}=1$, while the radius of the balls $r=a\ll 1$ (see Figure 1).

Assuming the excess pressure $p(\bm{r},t)$ to be time-harmonic $\displaystyle p(\bm{r},t)=u(\bm{r})\mathrm{e}^{-\mathrm{i}\omega t},$ the amplitude $u(\bm{r})$ should satisfy the equation

[TABLE]

where $\varrho(\bm{r})$ is the mass density and $\gamma(\bm{r})$ is the adiabatic bulk compressibility modulus of the media. Both $\varrho(\bm{r})$ and $\gamma(\bm{r})$ are periodic piecewise constant functions with the periods of the lattice. In what follows, we denote functions in the inclusions and the external medium using the subscripts $-$ and $+$, respectively. Thus,

[TABLE]

where

[TABLE]

We suppose that inclusions are penetrable and therefore impose the transmission conditions on their boundaries

[TABLE]

where the brackets $\llbracket\cdot\rrbracket$ denote the jump of the enclosed quantity across the interface $\partial\mathcal{B}$ of the inclusions. Solution of (2.2) is sought in the form of Floquet-Bloch waves

[TABLE]

where ${\bm{k}}=(k_{1},k_{2},k_{3})$ is the wave vector and $\Phi(\bm{r})$ is a periodic function with the periods of the lattice. The latter condition implies that function $\mathrm{e}^{-\mathrm{i}{\bm{k}}\cdot\bm{r}}u(\bm{r})$ is periodic. We write this as

[TABLE]

where the inverted brackets $\rrbracket\cdot\llbracket$ denote the jump of the enclosed expression and their first derivatives across the opposite sides of the cells of periodicity.

We reduce the above problem to the fundamental cell $\Pi$ centered at the origin:

[TABLE]

[TABLE]

In the inclusionless case there is a simple dispersion relation between the time frequency $\omega$ and the spacial frequency ${\bm{k}}$. Namely, there are waves propagating in any direction and $\omega=c|{\bm{k}}|$, where $c=1/\sqrt{\gamma_{+}\varrho_{+}}$ is the speed of waves in the host medium. The dispersion relation $\omega=H({\bm{k}},a)$ in the presence of inclusions is more complicated and our goal is to find it when $a$ is small. We will write the dispersion relation in the form $|{\bm{k}}|^{2}=G(\hat{\bm{k}},\omega,a),$ where $\hat{\bm{k}}={\bm{k}}/|{\bm{k}}|,$

Interaction of the incident wave with inclusions creates multiple Floquet-Bloch waves with the same direction. We consider the waves for which $|{\bm{k}}|$ is close to $\omega/c$ when $a$ is small. For other waves, $|{\bm{k}}+{\bm{L}}|$ is close to $\omega/c$, where ${\bm{L}}$ is an arbitrary period of the dual lattice. These waves can be studied similarly.

We consider both three and two dimensional cases ($d=3,2$). If $d=2$, the inclusions have the shape of infinite cylinders shown in Figure 2.

Theorem 1**.**

Let $d=3$. Then function $G$ in the dispersion equation

[TABLE]

*is infinitely smooth. Moreover, *

[TABLE]

where (2.12) has the following form in the case of the Neumann condition on inclusion boundaries

[TABLE]

Here $c=1/\sqrt{\gamma_{+}\varrho_{+}}$ is the speed of waves in the host medium, $|\Pi|$ is the volume of the cell $\Pi$ and $f$ is the filling fraction of the inclusions.

In the case of transmission conditions, we have

[TABLE]

Theorem 2**.**

Let $d=2$. Then the dispersion relation has the form

[TABLE]

Moreover,

[TABLE]

where (2.16) has the following form in the case of the Neumann condition on inclusion boundaries

[TABLE]

In the case of transmission conditions, we have

[TABLE]

where $\displaystyle\alpha=\frac{\varrho_{-}-\varrho_{+}}{\varrho_{+}+\varrho_{-}}$.

Note that the Neumann condition on inclusion boundaries can be obtained from the transmission conditions by passing to the limit $\lim\varrho_{-}\to\infty,\lim\gamma_{-}\varrho_{-}\to 0$. While formulas (2.13), (2.17) coincide with (2.14), (2.18), respectively, when $\varrho_{-}^{-1}=\gamma_{-}=0$, we do not justify these limiting transitions in the formulas. Instead, one can obtain the results for Neumann boundary conditions independently using the same approach.

Our approach can be also applied to the problem with the Dirichlet boundary condition. The answer in this case differs considerably from the transmission problem. Namely,

• •

$d=3$, Dirichlet boundary condition

[TABLE]

• •

$d=2$, Dirichlet boundary condition

[TABLE]

3 Outline of the proof

We begin with the observation that the unperturbed eigenvalue problem

[TABLE]

in the cell of periodicity $\Pi$ and a fixed ${\bm{k}}$ with $|{\bm{k}}|={k_{+}}$ has a simple eigenvalue $\lambda=0$ with the eigenfunction $u=\mathrm{e}^{-\mathrm{i}{\bm{k}}\cdot\bm{r}}$. This statement is independent of whether we solve (3.1) in the Sobolev space $H^{2}(\Pi)$ or in the space of infinitely smooth functions, since the ellipticity of the problem implies that solutions of (3.1) are infinitely smooth. The simplicity of a zero eigenvalue can be easily established using the substitution $\mathrm{e}^{-\mathrm{i}{\bm{k}}\cdot\bm{r}}u=v$ and expanding the periodic function $v$ into Fourier series.

Problem

[TABLE]

can be considered as a perturbation of problem (3.1). For arbitrary vector ${\bm{k}}$, we denote by ${{\bm{k}}_{+}}$ the vector with the same directions as ${\bm{k}}$ and with the magnitude ${k_{+}}$, i.e. $\displaystyle{{\bm{k}}_{+}}={k_{+}}\hat{\bm{k}}$. If $a$ and $\varepsilon={k_{+}}-|{\bm{k}}|$ are small, then the perturbation is small from the point of view of physics. Since $\lambda=0$ is a simple eigenvalue of the unperturbed problem, it follows that the eigenvalue $\lambda=\lambda(a,\varepsilon,\hat{{\bm{k}}})$ of problem (3.2) is smooth when $a$ and $|\varepsilon|$ are small. Hence one could find the dispersion relation by finding $\lambda(a,\varepsilon,\hat{{\bm{k}}})$ and solving the equation $\lambda(a,\varepsilon,\hat{{\bm{k}}})=0.$

It is not very easy to fulfill rigorously the approach discussed above, since problem (3.2) with $a\neq 0$ is a singular perturbation of problem (3.1), and the standard perturbation technique can not be applied. Thus we will consider problem (3.2) only with $\lambda=0$:

[TABLE]

and will reduce it to an equivalent operator equation

[TABLE]

for which the standard perturbation theory is valid. Here $\psi$ is a function on a sphere $r=R$ of a fixed radius $R<1$ (a circle if $d=2$), and operator $N^{+}_{{\bm{k}}}-N^{-}_{a}$ (which will be defined below) is symmetric and Fredholm. Moreover, $N^{+}_{{\bm{k}}}-N^{-}_{a}$ is infinitely smooth in ${\bm{k}},a$, and has a simple zero eigenvalue when ${\bm{k}}={{\bm{k}}_{+}},~{}a=0$. Thus, the dispersion relation can be found from (3.4) by the standard perturbation theory.

Let us define operators in (3.4). We split $\Pi$ into two domains $B_{R}=\{{\bm{r}}\mid r<R\}$, such that $B_{R}\subset\Pi$, and its compliment $\Omega=\Pi\smallsetminus B_{R}$ (see Figure 3). Constant $R$ will be fixed later, and $a\ll R$. Consider the following two separate problems in $\Omega$ and $B_{R}$, instead of problem (3.3).

[TABLE]

Both problems with $\psi\in C^{\infty}$ have unique $C^{\infty}$-solutions for all values of $R$, except, possibly, a discrete set $\{R_{i}\}$. We fix an $R\notin\{R_{i}\}$ and define operators ${\mathcal{N}}^{+}_{{\bm{k}}}$ and ${\mathcal{N}}^{-}_{a}$ as the Dirichlet-to Neumann operators with the derivatives in the direction of $r$:

[TABLE]

where $H^{s}(\partial B_{R})$ is the Sobolev space of functions on $\partial B_{R}$. Constant $a$ will be always assumed to be less than $R$.

4 Analysis of operators $\mathcal{N}^{+}_{{\bm{k}}},\mathcal{N}^{-}_{a}$

Lemma 1**.**

1. (1)

Operators ${\mathcal{N}}^{+}_{{\bm{k}}},~{}{\mathcal{N}}^{-}_{a}$ and

[TABLE]

are Fredholm. 2. (2)

Operators ${\mathcal{N}}^{+}_{{\bm{k}}}$ and ${\mathcal{N}}^{-}_{a}$ are symmetric in $L^{2}(\partial B_{R})$, i.e.

[TABLE]

and the same relation holds for ${\mathcal{N}}^{-}_{a}$.

Proof.

Dirichlet-to-Neumann operators (for elliptic equations of the second order) are elliptic pseudo-differential operators of the first order (see [24]). Thus the first statement for operators ${\mathcal{N}}^{+}_{{\bm{k}}},{\mathcal{N}}^{-}_{a}$ is obvious. If one proves that the difference ${\mathcal{N}}^{+}_{{\bm{k}}}-{\mathcal{N}}^{-}_{a}$ is still elliptic (i.e., the principal symbol of the difference is not vanishing), then the first statement for (4.1) becomes an immediate consequence of the ellipticity. This approach could be completed, since a calculation of the symbols for Dirichlet-to-Neumann operators can be found in [24]. We prefer to prove the ellipticity of the difference by considering a simplified version $n^{+}-n^{-}$ of operator ${\mathcal{N}}^{+}_{{\bm{k}}}-{\mathcal{N}}^{-}_{a}$, where

[TABLE]

are Dirichlet-to-Neumann operators related to the following problems (which are simplified versions of (3.5), (3.6)):

[TABLE]

[TABLE]

From local a priori estimates for elliptic operators it follows that a perturbation of the boundary value problem outside of a neighborhood of the boundary on which the Dirichlet-to-Neumann operator is defined does not change the symbol of a Dirichlet-to-Neumann operator. The principal symbols of Dirichlet-to-Neumann operators do not depend also on the lower order terms of the equations. Thus, the principal symbols of operators $n^{\pm}$ are the same as the ones for operators ${\mathcal{N}}^{+}_{{\bm{k}}},{\mathcal{N}}^{-}_{a}$, respectively. Hence, the first statement of the lemma will be proved if we prove that $n^{+}-n^{-}$ is an elliptic operator of the first order. For this purpose, it is enough to show that operator

[TABLE]

is an isomorphism. Consider the three-dimensional case (the two-dimensional case is similar). Denote by $\widetilde{\Delta}$ the Laplace-Beltrami operator on the sphere $\partial B_{R}$. Its spectrum consists of the eigenvalues $\lambda_{n}=n(n+1),~{}n\geqslant 0,$ of multiplicity $2n+1$. Let $\psi\in L^{2}(\partial B_{R})$ and let $\psi_{n}$ be the projection of $\psi$ into the eigenspace of $\widetilde{\Delta}$ with the eigenvalue $\lambda_{n}=n(n+1)$. Thus $\psi=\sum_{n=0}^{\infty}\psi_{n}$, and the norms of $\psi$ can be defined as follows:

[TABLE]

[TABLE]

Let us now estimate $(n^{+}-n^{-})\psi$. Obviously,

[TABLE]

and therefore

[TABLE]

Since $R<1$, the isomorphism of map (4.3) follows immediately from (4.4)-(4.6). This completes the proof of the first statement of the lemma.

Let us prove (4.2). From the symmetry of the problem (3.5) and Green’s second identity it follows that

[TABLE]

for solutions $u,v$ of (3.5) with data $\psi,\phi$ at $r=R$, respectively. The latter relation coincides with (4.2). ∎

Lemma 2**.**

• (1)

Relation $\psi=u|_{r=R}$ is a one-to-one correspondence between solutions $u\in C^{\infty}$ of (3.3) and solutions $\psi\in H^{1}(\partial B_{R})$ of (3.4).

• (2)

Zero is a simple eigenvalue of the operator $({\mathcal{N}}^{+}_{{\bm{k}}_{+}}-{\mathcal{N}}^{-}_{0})$ with the eigenfunction

[TABLE]

Proof.

Let $u\in C^{\infty}$ be a solution of (3.3). Then its restrictions to $\Omega$ and $B_{R}$ satisfy (3.5), (3.6), respectively, with $\psi=u|_{r=R}$. Moreover, relation (3.4) holds, since ${\mathcal{N}}^{+}_{{\bm{k}}}\psi={\mathcal{N}}^{-}_{a}\psi=u_{r}|_{r=R}$. Conversely, let $\psi\in H^{1}(\partial B_{R})$ satisfy (3.4). The ellipticity of the equation (3.4) implies that $\psi\in C^{\infty}$. Let $u$ coinside with solutions of (3.5), (3.6) in $\Omega$ and $B_{R}$, respectively. Then from (3.4) it follows that equation $(\Delta+k^{2}_{+})u=0$ holds also on $\partial R_{R}$, and therefore $u$ is a solution of (3.3). The first statement of the lemma is proved. The second statement follows from the first one applied to the unperturbed problem (${\bm{k}}={{\bm{k}}_{+}},a=0$) and the fact that zero is a simple eigenvalue of the problem (3.1) with the eigenfunction $\mathrm{e}^{-\mathrm{i}{{\bm{k}}_{+}}\cdot\bm{r}}|_{r=R}$. ∎

We will write each element $f$ in the domain $H^{1}(\partial B_{R})$ and range $L^{2}(\partial B_{R})$ of operator (4.1) in the vector form $f=(f_{1},f_{L}),$ where $f_{1}=c\psi_{0}$ is the projection of $f$ in $L^{2}(\partial B_{R})$ on element (4.7) and $f_{L}$ is orthogonal to $\hat{\psi}$ in $L^{2}(\partial B_{R})$. Denote by $L_{1},~{}L_{0}$ the subspaces of $H^{1}(\partial B_{R}),~{}L^{2}(\partial B_{R})$, respectively, that consist of functions orthogonal in $L^{2}(\partial B_{R})$ to $\psi_{0}$. Then, due to Lemmas 1,2, operator (4.1) has the following matrix form

[TABLE]

where $A:L_{1}\to L_{0}$ is an isomorphism.

Obviously, operator ${\mathcal{N}}^{+}_{{\bm{k}}}$ is infinitely smooth function of ${\bm{k}}$. Thus the operator ${\mathcal{N}}^{+}_{{\bm{k}}}-{\mathcal{N}}^{-}_{0}$ in the same basis chosen for ${\bm{k}}={{\bm{k}}_{+}}$ has a matrix representation

[TABLE]

where $C$ is a constant, $|\varepsilon|\ll 1$ and $D_{ij}=D_{ij}(\varepsilon,{k_{+}},\hat{{\bm{k}}}),~{}{k_{+}}>0,$ are infinitely smooth functions of all the arguments. We will evaluate constant $C$ in the next Lemma.

Lemma 3**.**

Constant $C$ in the matrix expansion (4.9) of the operator ${\mathcal{N}}^{+}_{\bm{k}}-{\mathcal{N}}^{-}_{0}$ is equal

[TABLE]

in both dimensions $d=2$ and $d=3$.

Proof.

Recall that $\varepsilon={k_{+}}-|{\bm{k}}|$ and ${{\bm{k}}_{+}}={k_{+}}\hat{{\bm{k}}}$. Hence ${\bm{k}}={{\bm{k}}_{+}}-\varepsilon\hat{{\bm{k}}}$, and the boundary condition on $\partial\Pi$ in (3.5) can be written in the form \left\rrbracket\mathrm{e}^{\mathrm{i}({{\bm{k}}_{+}}-\varepsilon\hat{{\bm{k}}})\cdot\bm{r}}u(\bm{r})\right\llbracket=0. Then problem (3.5) in $\Omega=\Pi\smallsetminus B_{R}$ becomes a regular perturbation of the same problem with ${\bm{k}}={{\bm{k}}_{+}}$. Thus constant $C$ in (4.9) is the coefficient in the linear term of the Taylor expansion of functions $q^{+}-q^{-}$ as $\varepsilon\to 0$ (i.e., ${\bm{k}}\to{{\bm{k}}_{+}}$), where

[TABLE]

Here $\hat{\psi}$ and $u^{+},v$ are solutions of (3.5), (3.6), respectively, with $\psi=\hat{\psi}$ and without the inclusion in problem (3.6).

Let us evaluate $q^{+}$. By the choice of $R$, problem (3.5) with ${\bm{k}}={{\bm{k}}_{+}}$ is uniquely solvable. If $\psi=\hat{\psi}$, then its solution is $\mathrm{e}^{-\mathrm{i}{{\bm{k}}_{+}}\cdot\bm{r}}$. Thus, solution $u^{+}$ of (3.5) can be expanded into the Taylor series in $\varepsilon$ with the zero order term being $\mathrm{e}^{-\mathrm{i}{{\bm{k}}_{+}}\cdot\bm{r}}$, i.e.,

[TABLE]

where the remainder term is small uniformly together with all its derivatives in $\bm{r}$, and $u_{1}$ is the solution of the problem

[TABLE]

We note that

[TABLE]

Thus from (4.11) and (4.13) it follows that

[TABLE]

In order to evaluate the integral above, we make the substitution $u_{1}(\bm{r})=w(\bm{r})+\mathrm{i}(\hat{{\bm{k}}}\cdot\bm{r})\mathrm{e}^{-\mathrm{i}{{\bm{k}}_{+}}\cdot\bm{r}}$, which reduces (4.13) to the following problem for $w$ with the homogeneous boundary condition on $\partial\Pi$:

[TABLE]

Denote the integral in (4.15) by $q^{+}_{1}$. Then

[TABLE]

We have

[TABLE]

The first integrand above is odd, and the corresponding integral is zero. Depending on the dimension $d=2,3$, the second integral above is $1/d$ times the Lebesgue measure of $\partial B_{R}$. Thus $\displaystyle q^{+}_{11}=\frac{{k_{+}}R}{d}|\partial B_{R}|$, where $|\partial B_{R}|=2\pi R$ in the dimension $d=2$ and $|\partial B_{R}|=4\pi R^{2}$ if $d=3$.

Using the Green formula and the symmetry of boundary problem (4.16), we obtain

[TABLE]

Hence

[TABLE]

Now we evaluate

[TABLE]

Function $v$ solves the problem

[TABLE]

Then

[TABLE]

The latter integrand is odd, and therefore $q^{-}=0$. This fact together with (4.17) imply (4.10).

∎

The next two lemmas are crucial in what follows, since they allow us to forget about the singularity of problem (3.3) as $a\to 0$.

Lemma 4**.**

Let $d=3$. Then the operator function ${\mathcal{N}}^{-}_{a}-{\mathcal{N}}^{-}_{0}:H^{1}(B_{R})\to L^{2}(B_{R})$ is infinitely smooth in $a\geqslant 0$, and there are bounded operators $P$ and $Q$ such that

[TABLE]

Moreover,

[TABLE]

where $c_{1}$ and $c_{2}$ are given by (2.14), (1) in Theorem 1.

Proof.

We prove this lemma using explicit construction of operator ${\mathcal{N}}^{-}_{a}-{\mathcal{N}}^{-}_{0}$. We denote the Laplace-Beltrami operator on the sphere $\partial B_{R}$ by $\widetilde{\Delta}$. Its spectrum consists of the eigenvalues $\lambda_{n}=n(n+1),~{}n\geqslant 0,$ of multiplicity $2n+1$ (and its eigenfunctions are the products of the associated Legendre polynomials of $\cos\theta$ and the exponents $e^{im\phi}$). Let $\psi_{n}$ be the projection in $L^{2}(\partial B_{R})$ of $\psi$ on the eigenspace of $\widetilde{\Delta}$ with the eigenvalue $\lambda_{n}=n(n+1)$. Thus $\psi=\sum_{n=0}^{\infty}\psi_{n}$. Then solution of problem (3.6) has the form

[TABLE]

where $j_{n},y_{n}$ are the spherical Bessel functions and constants $a_{n},b_{n},c_{n}$ are determined from the the boundary conditions in (3.6) . Solving the system above, we obtain that

[TABLE]

where

[TABLE]

Solution of the problem (3.6) without the inclusion has the form

[TABLE]

Hence, $\displaystyle{\mathcal{N}}^{-}_{0}\psi={k_{+}}\sum_{n=0}^{\infty}\frac{j_{n}^{\prime}({k_{+}}R)}{j_{n}({k_{+}}R)}\,\psi_{n}$, and

[TABLE]

where

[TABLE]

From the asymptotic behavior of Bessel functions $J_{n}(x),Y_{n}(x)$ [25] we have for $0<x\ll\sqrt{n}$ and $n\to\infty$

[TABLE]

Then

[TABLE]

where constants $C_{i}$ do not depend on $n$, and therefore $\displaystyle\left|\frac{f_{n}}{d_{n}\,j_{n}({k_{+}}R)}\right|\sim a^{2n+1}$ uniformly in $n\geqslant 1$ as $a\to 0$ and $\displaystyle\left|\frac{f_{0}}{d_{0}\,j_{0}({k_{+}}R)}\right|\sim a^{3}$. The latter relations and (4.27) justify the first statement of the lemma.

In order to obtain the exact values of coefficients $c_{1}$, $c_{2}$ in (4.19) and complete the proof of Lemma 4 we need asymptotic expansions of $f_{n},d_{n}$ for $n=0,1$ when $a\to 0$. They look as follows

[TABLE]

Then

[TABLE]

Let us substitute $\hat{\psi}$ for $\psi$ in (4.35) and evaluate quadratic form (4.19). We have

[TABLE]

Note that

[TABLE]

Differentiation of the first equality above with respect to $|{{\bm{k}}_{+}}|$ leads to

[TABLE]

Therefore

[TABLE]

and

[TABLE]

Hence

[TABLE]

Now formula (4.19) follows from the last two relations and (4.35). ∎

Lemma 5**.**

Let $d=2$. Then the operator function ${\mathcal{N}}^{-}_{a}-{\mathcal{N}}^{-}_{0}:H^{1}(B_{R})\to L^{2}(B_{R})$ has the following asymptotic expansion as $a\to 0$

[TABLE]

Moreover,

[TABLE]

where $c_{1}$ and $c_{2}$ are given in Theorem 2.

Proof.

The proof is similar to that in Lemma 4. Namely, solution $v$ of the problem (3.6) is represented in the form

[TABLE]

where $J_{n}(x),Y_{n}(x)$ are Bessel functions. Constants $A_{n},B_{n},C_{n}$ are determined from the boundary conditions in (3.6):

[TABLE]

where

[TABLE]

Then

[TABLE]

where

[TABLE]

Using the asymptotic behavior of the Bessel function as $x\to 0$

[TABLE]

one can derive the form of asymptotic expansion of $D_{n}$ and $F_{n}$ as $a\to 0$

[TABLE]

In particular, for $n=0,1$ we have

[TABLE]

As a result,

[TABLE]

Substitution of (4.49) into (4.48) provides first statement (4.39) of the Lemma.

We use the values of the integrals

[TABLE]

From here we calculate projections $\psi_{0}$, $\psi_{1}$ of the function $\mathrm{e}^{-\mathrm{i}{\bm{k}}_{+}\cdot\bm{r}}$:

[TABLE]

Hence

[TABLE]

and therefore

[TABLE]

Addition of the above formulas gives (4.40) in terms of notation (2.3),(2.4). ∎

5 Proof of Theorems 1, 2 and formulas (2.19), (2.20)

In order to prove Theorem 1, we use Lemma 2 and reduce problem (3.3) to the equivalent equation (3.4). Due to Lemma 1, the operator in the latter equation is Fredholm and symmetric. Let us rewrite (3.4) in the form

[TABLE]

Obviously, the first term on the left is infinitely smooth in ${\bm{k}}$ and has the matrix representation (4.9) in the decompositions $(\hat{\psi},L)$ of the domain and range of the operator. Recall that $\hat{\psi}$ is defined by (4.7) and $L$ is the subspace of the domain or range, respectively, which consists of functions orthogonal (in $L^{2}(\partial B_{R})$) to $\hat{\psi}$. Due to Lemma 4, the second operator in (5.1) is an infinitely smooth function of $a$. Hence equation (5.1) can be written in the form

[TABLE]

where $\phi_{1}=\sigma\hat{\psi}$ is the projection of $\psi$ on $\hat{\psi}$, component $\phi_{2}$ of $\psi$ is orthogonal to $\hat{\psi}$, operator functions $D_{ij}=D_{ij}(\varepsilon,{k_{+}},\hat{{\bm{k}}}),~{}{k_{+}}>0,$ and $B_{1j}(a)$ are infinitely smooth, and asymptotic expansions of operators $B_{ij}$ as $a\to 0$ are given in Lemma 4. In particular, $\|B_{ij}\|=O(a^{3}),~{}a\to 0$.

Since operator $A$ is invertible, the second equation of the system (5.2) can be solved for $\phi_{2}$, and the solution is proportional to the norm of $\phi_{1}$, i.e.,

[TABLE]

where $\hat{\phi_{2}}$ is a $\sigma$-independent element of $L$ which is infinitely smooth in $\varepsilon,{k_{+}},\hat{{\bm{k}}}$ and such that $\|\hat{\phi_{2}}\|=O(|\varepsilon|+a^{3}),~{}|\varepsilon|+a\to 0$. Then the first equation of (5.2) implies

[TABLE]

Hence, a non-trivial solution of (5.2) for small $\varepsilon,a$ exists if and only if

[TABLE]

where the left-hand side is an infinitely smooth function of $\varepsilon,{k_{+}},\hat{{\bm{k}}}$ (when ${k_{+}}>0)$ and operator $B_{11}(a)$ coincides with the left-hand side in (4.19). The dispersion relation is defined by solving the equation above for $\varepsilon$. We obtain that $|{\bm{k}}|-{k_{+}}$ is an infinitely smooth function of $a,{k_{+}},\hat{{\bm{k}}}$ and

[TABLE]

This justifies the infinite smoothness of the dispersion relation (2.11). In order to justify expansion (2.12), we use (4.19) for $B_{11}(a)$ and formula for $C$ from Lemma 3 and arrive at

[TABLE]

From here it follows that $|{\bm{k}}|-{k_{+}}=O(a^{3})$ and therefore $2{k_{+}}=|{\bm{k}}|+{k_{+}}+O(a^{3})$. The latter relation together with (5.4) implies (2.12) and completes the proof of Theorem 1. We also note that the first correction term in (2.12) agrees with that provided in [26].

Proof of Theorem 2 is absolutely similar. One need only to use Lemma 5 instead of Lemma 4. Hence all the functions of $a$ will have asymptotic expansions as in (4.39) instead of being infinitely smooth. Operators $B_{ij}$ will have order $O(a^{2}),~{}a\to 0$, instead of $O(a^{3})$, $B_{11}$ will coincide with the left-hand side in (4.40), and (5.3) will have the form

[TABLE]

In the case of the Dirichlet boundary condition, norms of operators $B_{ij}$ have order $O(a)$ and $O(1/\ln a)$ when $d=3,2$, respectively, and this leads to different dispersion relations.

6 Effective wave velocity

Asymptotic formulas (2.12), (2.16) allow us to calculate both the effective phase velocity $\displaystyle c_{\textrm{ph}}=\frac{\omega}{k}$ and the group velocity $c_{\ast}=\displaystyle\frac{\partial\omega}{\partial k}$. With the accuracy of the asymptotic formulas we have

• •

$d=2$, transmission boundary conditions:

[TABLE]

• •

$d=2$, Neumann boundary conditions:

[TABLE]

• •

$d=3$, transmission boundary conditions:

[TABLE]

• •

$d=3$, Neumann boundary conditions:

[TABLE]

Here $c_{1}$, $c_{2}$ are defined by (2.14),(1) in three dimensions and are given by (2.18) for $d=2$. The formulas are consistent with the results reported in [27], [28], [29]. If $\gamma_{-}/\gamma_{+}\gg 1$ as in the case of air bubbles in water where $\gamma_{-}/\gamma_{+}\sim 1.5\cdot 10^{4}$, then even a tiny concentration of the air bubbles causes a dramatic reduction of the wave speed. This effect was confirmed numerically in [30] and [28], however it is not attributed to the Minnaert resonance or wave localization. We rewrite first correction terms (2.14) and (2.18) in formulas (2.12) and (2.16), respectively, in terms of the wave speed $c_{\pm}$ in the constituent media

[TABLE]

The first terms in the numerators of the formulas indicate the mismatch in the characteristic impedances $\varrho c$ of the media. These formulas show that if the scatterers are slow ($c_{-}<c_{+}$) then the effective wave speed will always be slower than the wave speed in the host medium. On the other hand, the presence of fast scatterers $(c_{-}>c_{+})$ does not guarantee that the effective wave velocity will be larger than that in the matrix. It will be so only when the characteristic impedances of the two media are close enough so that their mismatch does not cause strong scattering.

Figure 4 depicts the comparison of the asymptotic formulas of the present paper with numerical calculations of the sound velocity for a hexagonal cluster of $151$ fluid cylinders embedded in a fluid or gas (left panel) reported in [28] when $\omega\to 0$. The right panel shows dependence of the long wavelength sound group velocity for a lattice of spherical inclusions of the same constituents calculated by formulas (6.9), (6.11). Compared with the two-dimensional case, the sound velocity decays about twice as slowly. We note that the only geometric parameter the sound velocity depends on in long wavelength approximation is the filling fraction, $f$.

Increasing the wave frequency $\omega$ enhances the scattering and makes the decay of the sound velocity even faster as shown in Figure 5 for $\omega/c=1$ and $|\Pi|=1$. This effect is also illustrated in Figure 6 when the group velocity of sound in an orthorhombic lattice of rigid spherical inclusions decreases as $\omega/c$ increases from [math] to $3$.

7 Conclusions

We have derived dispersion relations for the acoustic waves propagating in homogeneous medium containing a periodic lattice of spherical or cylindrical inclusions. The wavelength is assumed to be of the order of the periods of the lattice while the radius $a$ of the inclusion is small compared to the periods. We suggest a new approach to derive and justify asymptotic expansions of the dispersion relations in two and three-dimensional cases as $a\to 0$ and evaluate explicitly several first terms. The approach is based on the reduction of the original singularly perturbed (by inclusions) problem to an equivalent regular one. In order to get the latter problem we split the cell of periodicity into two parts by introducing an auxiliary spherical boundary of a fixed radius. We replace the original problem with the equality of Dirichlet-to-Neumann maps on the auxiliary boundary of the new two subdomains. The Neumann, Dirichlet and transmission boundary conditions on inclusion boundaries are considered. The effective wave speed is obtained as a function of the wave frequency, the filling fraction of the inclusions, and the physical properties of the constituents of the mixture. The method can also be extended to small inclusions of arbitrary shape.

Acknowledgments

The work of B.V. was partially supported by the NSF grant DMS-1714402 and the Simons Foundation grant 527180