Elastostatics with multi-layer metamaterial structures and an algebraic framework for polariton resonances

TL;DR

This paper develops a mathematical framework for elastostatics in multi-layer metamaterial structures, analyzing their spectral properties and establishing an algebraic approach to study polariton resonances relevant for advanced metamaterial applications.

Contribution

It introduces a general algebraic framework for analyzing elastostatics and polariton resonances in multi-layer elastic metamaterials with arbitrary layers and materials.

Findings

Derived exact perturbed field expressions using elastic momentum matrices.

Analyzed spectral properties of the elastic momentum matrix.

Established a comprehensive algebraic framework for polariton resonances.

Abstract

Multi-layer structures are ubiquitous in constructing metamaterial devices to realise various frontier applications including super-resolution imaging and invisibility cloaking. In this paper, we develop a general mathematical framework for studying elastostatics within multi-layer material structures in $\mathbb{R}^d$, $d=2,3$. The multi-layer structure is formed by concentric balls and each layer is filled by either a regular elastic material or an elastic metamaterial. The number of layers can be arbitrary and the material parameters in each layer may be different from one another. In practice, the multi-layer structure can serve as the building block for various material devices. Considering the impingement of an incident field on the multi-layer structure, we first derive the exact perturbed field in terms of an elastic momentum matrix, whose dimension is the same as the number of layers. By highly intricate and delicate analysis, we derive a comprehensive study of the spectral properties of the elastic momentum matrix. This enables us to establishe a handy algebraic framework for studying polariton resonances associated with multi-layer metamaterial structures, which forms the fundamental basis for many metamaterial applications.

Tables

Table 1. Table 1. Roots to the characteristic polynomial with layers which are chosen by ( 6.1 ).

				$N=19$ and $\xi =3$
$q$		2.9691		2.9496		2.9185
$\beta$	0	2.9923	-0.9923	2.9874	-0.9874	2.9795	-0.9795
${\epsilon }_i$	-3.0000	-0.0019	-517.0116	-0.0032	-315.5802	-0.0051	-194.3029
$q$		2.8668		2.7735		2.5941
$\beta$		2.9664	-0.9664	2.9425	-0.9425	2.8958	-0.8958
${\epsilon }_i$		-0.0085	-118.0785	-0.0146	-68.6177	-0.0267	-37.3961
$q$		2.2295		1.5099		0.4949
$\beta$		2.7971	-0.7971	2.5843	-0.5843	2.2227	-0.2227
${\epsilon }_i$		-0.0534	-18.7110	-0.1160	-8.6217	-0.2412	-4.1458
				$N=19$ and $\xi =0.8$
$q$		0.7918		0.7866		0.7783
$\beta$	0	0.7954	-0.9954	0.7925	-0.9925	0.7878	-0.9878
${\epsilon }_i$	-0.8000	-0.0026	-392.1546	-0.0042	-239.1326	-0.0068	-147.0737
$q$		0.7645		0.7396		0.6918
$\beta$		0.7800	-0.9800	0.7658	-0.9658	0.7377	-0.9377
${\epsilon }_i$		-0.0112	-89.1786	-0.0194	-51.6170	-0.0358	-27.9008
$q$		0.5945		0.4026		0.1320
$\beta$		0.6775	-0.8775	0.5424	-0.7424	0.2768	-0.4768
${\epsilon }_i$		-0.0730	-13.6950	-0.1670	-5.9869	-0.4098	-2.4404

Table 2. Table 2. Roots to the characteristic polynomial with layers which are chosen by ( 6.2 ).

				$N=19$ and $\xi =3$
$q$		2.2677		2.2227		2.1412
$\beta$	0	2.8077	-0.8077	2.7952	-0.7952	2.7724	-0.7724
${\epsilon }_i$	-3.0000	-0.0505	-19.7986	-0.0540	-18.5314	-0.0603	-16.5712
$q$		2.0119		1.8171		1.5334
$\beta$		2.7355	-0.7355	2.6784	-0.6784	2.5917	-0.5917
${\epsilon }_i$		-0.0708	-14.1215	-0.0874	-11.4391	-0.1137	-8.7959
$q$		1.1403		0.6549		0.1979
$\beta$		2.4630	-0.4630	2.2864	-0.2864	2.0945	-0.0945
${\epsilon }_i$		-0.1551	-6.4484	-0.2171	-4.6057	-0.2926	-3.4174
				$N=19$ and $\xi =0.8$
$q$		0.6047		0.5927		0.5710
$\beta$	0	0.6840	-0.8840	0.6764	-0.8764	0.6622	-0.8622
${\epsilon }_i$	-0.8000	-0.0689	-14.5228	-0.0738	-13.5583	-0.0829	-12.0654
$q$		0.5365		0.4846		0.4089
$\beta$		0.6393	-0.8393	0.6033	-0.8033	0.5472	-0.7472
${\epsilon }_i$		-0.0981	-10.1980	-0.1227	-8.1490	-0.1634	-6.1211
$q$		0.3041		0.1746		0.0527
$\beta$		0.4604	-0.6604	0.3297	-0.5297	0.1505	-0.3505
${\epsilon }_i$		-0.2325	-4.3008	-0.3537	-2.8274	-0.5645	-1.7715

Equations

$$D_0:=\{r>r_1\},D_j:=\{r_{j+1}<r⩽r_j\},j=1,2,\dots ,N-1,D_N:=\{r⩽r_N\},$$

$$S_j:=\left\{|x|=r_j\right\},j=1,2,\dots ,N,$$

$$(\tilde{\lambda },\tilde{\mu })=\epsilon (\mathbf{x})\cdot (\lambda ,\mu )\text{mbox}and\epsilon (\mathbf{x})={\epsilon }_c(\mathbf{x})\chi (D)+{\epsilon }_0\chi ({\mathbb{R}}^d\backslash \overline{D}),$$

$$\mu >0\text{mbox}andd\lambda +2\mu >0.$$

$${\epsilon }_c(\mathbf{x})={\epsilon }_j,\mathbf{x}\in D_j,j=1,2,\dots ,N.$$

$${\epsilon }_j\in \mathbb{C}\text{mbox}with\Im {\epsilon }_j\ge 0,j=1,2,\dots ,N.$$

$${\mathbf{C}}_{\tilde{\lambda },\tilde{\mu }}=({\mathrm{C}}_{ijkl}),{\mathrm{C}}_{ijkl}(\mathbf{x}):=\tilde{\lambda }(\mathbf{x}){\mathbf{\delta }}_{ij}{\mathbf{\delta }}_{kl}+\tilde{\mu }(\mathbf{x})({\mathbf{\delta }}_{ik}{\mathbf{\delta }}_{jl}+{\mathbf{\delta }}_{il}{\mathbf{\delta }}_{jk}),\mathbf{x}=({\mathbf{x}}_1,\dots ,{\mathbf{x}}_d)\in {\mathbb{R}}^d,$$

$${\mathcal{L}}_{\tilde{\lambda },\tilde{\mu }}\mathbf{u}:=\nabla \cdot {\mathbf{C}}_{\tilde{\lambda },\tilde{\mu }}{\nabla }^s\mathbf{u}=\tilde{\mu }\Delta \mathbf{u}+(\tilde{\lambda }+\tilde{\mu })\nabla \nabla \cdot \mathbf{u},$$

$${\nabla }^s\mathbf{u}=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T),$$

$$\left\{\begin{array}[]{ll}\mathcal{L}_{\tilde{\lambda},\tilde{\mu}}\mathbf{u}=0,&\mbox{in}\quad\mathbb{R}^{d},\\ \mathbf{u}-\mathbf{H}=\mathcal{O}(|\mathbf{x}|^{1-d}),&|\mathbf{x}|\rightarrow\infty,\end{array}\right.$$

$${\mathcal{T}}_n^m(\mathbf{x})=\nabla (r^n{Y}_n^m(\hat{\mathbf{x}}))\times \mathbf{x},n⩾1,-n⩽m⩽n,$$

$${\mathcal{M}}_n^m(\mathbf{x})=\nabla (r^n{Y}_n^m(\hat{\mathbf{x}})),n⩾1,-n⩽m⩽n,$$

$${\mathcal{N}}_n^m(\mathbf{x})={\gamma }_n^m{r}^{n-1}{Y}_{n-1}^m(\hat{\mathbf{x}})\mathbf{x}+(1-\frac{{\gamma }_n^m}{2n-1}-r^2)\nabla (r^{n-1}{Y}_{n-1}^m(\hat{\mathbf{x}})),$$

$${\gamma }_n^m=\frac{2(n-1)\lambda +2(3n-2)\mu }{(n+2)\lambda +(n+4)\mu },n⩾1,-(n-1)⩽m⩽n-1,$$

$${\tilde{\mathcal{T}}}_n^m(\mathbf{x})=\nabla (r^{-n-1}{Y}_n^m(\hat{\mathbf{x}}))\times \mathbf{x},n⩾1,-n⩽m⩽n,$$

$${\tilde{\mathcal{M}}}_n^m(\mathbf{x})={\tilde{\gamma }}_n^m{r}^{-n-2}{Y}_{n-1}^m(\hat{\mathbf{x}})\mathbf{x}+(1+\frac{{\tilde{\gamma }}_n^m}{2n+1}-r^2)\nabla (r^{-n-2}{Y}_{n-1}^m(\hat{\mathbf{x}})),n⩾1,-(n-1)⩽m⩽n-1,$$

$${\tilde{\gamma }}_n^m=\frac{2(n+1)\lambda +(6n+11)\mu }{(n-1)\lambda +(n-5)\mu },$$

$${\tilde{\mathcal{N}}}_n^m(\mathbf{x})=\nabla (r^{-n-1}{Y}_n^m(\hat{\mathbf{x}})),n⩾1,-n⩽m⩽n,$$

$$\mathbf{H}=\sum _{n=1}^{\infty }\sum _{m=-n}^n({\mathfrak{t}}_{n,m}^0{\mathcal{T}}_n^m+{\mathfrak{m}}_{n,m}^0{\mathcal{M}}_n^m+{\mathfrak{n}}_{n,m}^0{\mathcal{N}}_n^m).$$

$$\mathbf{u}=\{\begin{array}{rlrl}& \sum _{n=1}^{\infty }\sum _{m=-n}^n({\mathfrak{t}}_{n,m}^N{\mathcal{T}}_n^m+{\mathfrak{m}}_{n,m}^N{\mathcal{M}}_n^m+{\mathfrak{n}}_{n,m}^N{\mathcal{N}}_n^m),& \mathbf{x}& \in D_N,\\ & \sum _{n=1}^{\infty }\sum _{m=-n}^n({\mathfrak{t}}_{n,m}^j{\mathcal{T}}_n^m+{\mathfrak{m}}_{n,m}^j{\mathcal{M}}_n^m+{\mathfrak{n}}_{n,m}^j{\mathcal{N}}_n^m)& & \\ & +\sum _{n=1}^{\infty }\sum _{m=-n}^n({\tilde{\mathfrak{t}}}_{n,m}^j{\tilde{\mathcal{T}}}_n^m+{\tilde{\mathfrak{m}}}_{n,m}^j{\tilde{\mathcal{M}}}_n^m+{\tilde{\mathfrak{n}}}_{n,m}^j{\tilde{\mathcal{N}}}_n^m),& \mathbf{x}& \in D_j,j=N-1,N-2,\dots ,1,\\ & \mathbf{H}+\sum _{n=1}^{\infty }\sum _{m=-n}^n({\tilde{\mathfrak{t}}}_{n,m}^0{\tilde{\mathcal{T}}}_n^m+{\tilde{\mathfrak{m}}}_{n,m}^0{\tilde{\mathcal{M}}}_n^m+{\tilde{\mathfrak{n}}}_{n,m}^0{\tilde{\mathcal{N}}}_n^m),& \mathbf{x}& \in D_0.\end{array}$$

$$\mathbf{H}=a_0{\mathcal{N}}_1^0(\mathbf{x}).$$

$$\mathbf{u}=\{\begin{array}{rlrl}& a_N{\mathcal{N}}_1^0,& \mathbf{x}& \in D_N,\\ & a_j{\mathcal{N}}_1^0+r^{-3}{b}_j{\mathcal{N}}_1^0,& \mathbf{x}& \in D_j,j=N-1,N-2,\dots ,1.\\ & a_0{\mathcal{N}}_1^0+r^{-3}{b}_0{\mathcal{N}}_1^0,& \mathbf{x}& \in D_0,\end{array}$$

$${\mathbf{u}|}_{-}={\mathbf{u}|}_{+}\text{ and }{{\epsilon }_j\frac{\partial \mathbf{u}}{\partial \mathbf{\nu }}|}_{-}={{\epsilon }_{j-1}\frac{\partial \mathbf{u}}{\partial \mathbf{\nu }}|}_{+}\text{ on }{S}_j,$$

$$\frac{\partial \mathbf{u}}{\partial \mathbf{\nu }}=\lambda (\nabla \cdot \mathbf{u})\mathbf{\nu }+\mu \left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right)\mathbf{\nu }$$

$${v|}_{\pm }(\mathbf{x})=\underset{h\to 0^{+}}{\lim }v(\mathbf{x}\pm h\mathbf{\nu }),\mathbf{x}\in S_j,$$

$${\beta }_j=\frac{\frac{4\mu }{3\lambda +2\mu }{\epsilon }_{j-1}+{\epsilon }_j}{{\epsilon }_{j-1}-{\epsilon }_j},j=1,2,\dots ,N.$$

$${\mathbb{P}}_N:=\left[\begin{array}{ccccc}{\beta }_1& -1& -1& \cdots & -1\\ \frac{4\mu }{3\lambda +2\mu }(r_2/r_1{)}^3& {\beta }_2& -1& \cdots & -1\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ \frac{4\mu }{3\lambda +2\mu }(r_{N-1}/r_1{)}^3& \frac{4\mu }{3\lambda +2\mu }(r_{N-1}/r_2{)}^3& \frac{4\mu }{3\lambda +2\mu }(r_{N-1}/r_3{)}^3& \cdots & -1\\ \frac{4\mu }{3\lambda +2\mu }(r_N/r_1{)}^3& \frac{4\mu }{3\lambda +2\mu }(r_N/r_2{)}^3& \frac{4\mu }{3\lambda +2\mu }(r_N/r_3{)}^3& \cdots & {\beta }_N\end{array}\right]$$

$$\mathbf{u}-\mathbf{H}=r^{-3}{\mathbf{e}}^T{{\rm Y}}_N({\mathbb{P}}_N^T{)}^{-1}\mathbf{e}\mathbf{H},$$

$${{\rm Y}}_N:=\left[\begin{array}{ccccc}{r}_1^3& 0& 0& \cdots & 0\\ 0& r_2^3& 0& \cdots & 0\\ 0& 0& r_3^3& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & r_N^3\end{array}\right].$$

$$\mathbf{u}=\{\begin{array}{rlrl}& {\mathfrak{t}}_{n,m}^N{\mathcal{T}}_n^m,& \mathbf{x}& \in D_N,\\ & {\mathfrak{t}}_{n,m}^j{\mathcal{T}}_n^m+{\tilde{\mathfrak{t}}}_{n,m}^j{\tilde{\mathcal{T}}}_n^m,& \mathbf{x}& \in D_j,j=N-1,N-2,\dots ,1,\\ & {\mathfrak{t}}_{n,m}^0{\mathcal{T}}_n^m+{\tilde{\mathfrak{t}}}_{n,m}^0{\tilde{\mathcal{T}}}_n^m,& \mathbf{x}& \in D_0.\end{array}$$

Full text

Keywords:  multi-layer solid structure, negative elastic materials, polariton resonance, elastic momentum matrix, spectral analysis, characteristic polynomial

2020 Mathematics Subject Classification:   35B34, 35L05, 35P15, 74B05, 74J20

1. Introduction

We are concerned with the mathematical study of elastic deformation within multi-layer solid structures and the peculiar resonance phenomena it induces. Specifically, we allow the presence of negative elastic materials, i.e. elastic materials with negative bulk moduli which form an important class of the so-called elastic metamaterials. In fact, exotic elastic materials with negative stiffness have been artificially structured [19, 20], and have been proposed for many revolutionary applications including super-resolution elastic imaging and invisibility cloaking. In constructing metamaterial devices to realise those frontier applications, multi-layer structures are frequently used as the building blocks. Moreover, we would like to point out that the so-called polariton resonances induced by those negative structures form the fundamental basis of those applications.

In [11], the authors studied the surface polariton resonance associated with a single elastic metamaterial particle in linear elasticity and considered its application in elastic wave imaging. This can be regarded as a single-layer metamaterial structure. In [8, 9, 22, 26], the core-shell structures were proposed for realising invisibility cloaking of elastic deformation based on the so-called anomalous localised resonance. Here, the core is filled with a regular elastic material whereas the shell is filled with a negative elastic material, and both the regular and negative elastic materials are uniform. If the material parameters in the core and the shell match in a delicate way involving the geometric configurations as well, then anomalous localised resonance can occur which in turn induces the cloaking effect due to an impinging field. The core-shell device can be regarded as a two-layer metamaterial structure. The aforementioned studies have also been extended from the micro-scale to the macro-scale in [12, 13, 24]. Similarly, one-layer or two-layer metamaterial structures have also been extensively used and theoretically investigated in optics and acoustics; see e.g. [3, 4, 10, 14, 15, 17, 18, 21, 23, 25, 28, 29] and the references cited therein; and in particular in [16], multi-layer metamaterial structures have been studied in electrostatics. Motivated by those studies, we consider elastostatic deformation within much more general layered metamaterial structures in $\mathbb{R}^{d}$, $d=2,3$. In fact, we consider the case that the multi-layer structure is formed by concentric balls and each layer is filled by either a regular elastic material or an elastic metamaterial. The number of layers can be arbitrary and the material parameters in each layer may be different, though uniform. This can cover many of the existing studies mentioned above with the number of layers being 1 or 2. Then we consider the polariton resonance that can be induced by such general layered elastic metamaterial structures. To that end, we first derive the elastostatic field within the multi-layer structure due to an impinging field in terms of the so-called elastic momentum matrix, whose dimension is the same as the number of the layers. By highly intricate and delicate analysis, we derive a comprehensive study of the spectral properties of the elastic momentum matrix. It in turns establishes a convenient and handy algebraic framework for studying polariton resonances associated with multi-layer metamaterial structures. In practice, the multi-layer structure can serve as the building block for various material devices, and our algebraic framework can facilitate the proper selection of the metamaterial parameters in order to induce the desired resonance and in turn the realisation of customised applications. We shall investigate along this direction in our forthcoming work.

There is one more motivation and physical relevance of our study which can be described as follows. In fact, as mentioned earlier, the multi-layer structure in our study can be completely occupied by regular elastic materials, i.e. there is no presence of metamaterials. Such regular multi-layer structures have been proposed for achieving the so-called GPT-vanishing (Generalised Polarisation Tensors) or ESC-vanishing (Elastic Scattering Coefficients) structures and hence cloaking devices with enhanced invisibility effects via the transformation approach; see [1, 27] for the related study in the elastic case, [7] in the electromagnetic case, [5] in the electrostatic case and [6] in the acoustic case. The key technical ingredient in those studies is an algebraic system defined by certain integral operators whose coefficient matrix is exactly the momentum matrix derive in this paper in the elastic case. The explicit formula of the characteristic polynomial and the estimation of the corresponding roots of the momentum matrix are set to be open problems in [5, Page 258] and [6, Page 497]:

1. (1)

• However, as the number of layers gets larger, solving analytically the characteristic polynomial seems too complicated, and even proving existence of solutions to the characteristic polynomial seems to be quite challenging. These numerical evidences show us that the characteristic polynomial has solutions, even though we are not able to prove it.

2. (2)

*As in the conductivity case [5], it should be emphasized that one does not know if a solution exists for any number of layers. *

Due to such a reason, in the works mentioned above, the authors only considered the numerical verification of the existence of such multi-layer structures with a small number of layers. Clearly, in the current article, we can completely address these open problems, at least in the elastostatic case. This paves the way for constructing much more general GPT-vanishing structures for enhanced invisibility cloaking devices by following a similar spirit as the works mentioned above. However, we choose to explore along that direction somewhere else and to focus on the development of the algebraic framework for the polariton resonance in the current paper.

The remainder of this paper is organized as follows. In Section 2, we first present the elastostatic scattering problem with multi-layer structures. Second, we give the representation of the perturbed field in terms of the elastic moment matrix. In Section 3, we establish an algebraic framework for polariton resonances and give the main results on the explicit formula of the characteristic polynomial and the estimation of the roots. Section 4 is devoted to the proofs of the main results in Section 3. In Section 5, we give some similar results in the two dimensional setting. In Section 6, numerical examples are presented in finding all the polariton resonance modes in a fixed multi-layer structure and polariton resonance is simulated. Some conclusions are made in Section 7.

2. Elastostatics with multi-layer metamaterial structures

2.1. Multi-layer metamaterial structures

We give the general $N$-layer elastic structure as shown in FIGURE 1.

Precisely, we construct a sequence of layers, $D_{0},D_{1},\ldots,D_{N}$, by

[TABLE]

and denote the interfaces between the adjacent layers by

[TABLE]

where $N\in\mathbb{N}$ and $r_{j}\in\mathbb{R}_{+}$. The elastic structure configuration in $\mathbb{R}^{d}$, $d=2,3$, is characterized by the Lamé parameters $\tilde{\lambda}$ and $\tilde{\mu}$ which are given by

[TABLE]

where $\chi$ denotes the characteristic function; $\mathbb{R}^{d}\backslash\overline{D}$ is the matrix, namely the background space; and $D=\cup_{j=1}^{N}D_{j}$ is the $N$-layer elastic structure. In (2.2), $\lambda$ and $\mu$ signify the bulk moduli of a regular elastic material which are real constants satisfying the following strong convexity conditions,

[TABLE]

$(\lambda,\mu)$ serves as the reference elastic configuration in our study. $\varepsilon_{0}$ is set to be a positive constant and hence $\varepsilon_{0}\cdot(\lambda,\mu)$ defines a regular elastic material in the matrix $\mathbb{R}^{d}\backslash\overline{D}$. For the $N$-layer structure, we set $\varepsilon_{c}$ to be of the following form

[TABLE]

where

[TABLE]

That is, $\varepsilon_{c}$ is of a layered-piecewise constant structure. If $\Re\varepsilon_{j}>0$, then $\varepsilon_{j}\cdot(\lambda,\mu)$ defines a regular elastic material in $D_{j}$ with the bulk moduli being $\Re\varepsilon_{j}\cdot(\lambda,\mu)$ which clearly fulfils the strong convexity conditions in (2.3), and $\Im\varepsilon_{j}\cdot(\lambda,\mu)$ signifying the lossy/dampling parameters of the material. If $\Re\varepsilon_{j}\leq 0$, then $\varepsilon_{j}\cdot(\lambda,\mu)$ defines an elastic metamaterial in $D_{j}$ with the bulk moduli $\Re\varepsilon_{j}\cdot(\lambda,\mu)$ clearly broken the convexity conditions (2.3). In this paper, we simply refer to $\varepsilon_{j}\cdot(\lambda,\mu)$ with $\Re\varepsilon_{j}\leq 0$ as a negative elastic material in $D_{j}$.

In summary, we consider a rather general multi-layer structure where the number of layers can be arbitrarily given and the material in each layer can either be regular or negative. Moreover, the material parameters in each layer can be different from one another.

2.2. Elastostatic scattering and momentum matrix

Let $\mathbf{C}_{\lambda,\mu}(\mathbf{x})=(\mathrm{C}_{ijkl}(\mathbf{x}))_{i,j,k,l=1}^{d}$ be a four-rank tensor such that

[TABLE]

where $\bm{\delta}$ is the Kronecker delta, and $\tilde{\lambda}$ and $\tilde{\mu}$ are the elastic moduli defined in (2.2). The corresponding Lamé operator $\mathcal{L}_{\tilde{\lambda},\tilde{\mu}}$ is defined by

[TABLE]

with the strain tensor

[TABLE]

where $\mathbf{u}$ is a $\mathbb{C}^{d}$-valued function, which signifies the elastic deformation field, and $T$ denotes the the transpose of a matrix.

Associated with the multi-layer elastic structure configuration described by (2.2)–(2.4), the elastostatic scattering is governed by the following Lamé system:

[TABLE]

where the displacement field $\mathbf{u}\in H^{1}_{loc}(\mathbb{R}^{d})^{d}$, and the background field $\mathbf{H}$ satisfies $\mathcal{L}_{{\lambda},{\mu}}\mathbf{H}=0$ in $\mathbb{R}^{d}$ with $\mathcal{L}_{{\lambda},{\mu}}$ defined in (2.6) with the subscripts replaced by ${\lambda}$ and ${\mu}$.

In what follows, we shall confine our study of the electrostatic system (2.7) within the three dimensional setting, and shall present the corresponding two-dimensional extension in Section 5. Recall that, in [12, Lemma 3.2], three types of vectorial polynomials are introduced:

[TABLE]

[TABLE]

and

[TABLE]

where

[TABLE]

which are solutions to the Lamé $\mathcal{L}_{{\lambda},{\mu}}\mathbf{u}=0$ in $\mathbb{R}^{d}$. By straight forward computations, one can also verify that

[TABLE]

[TABLE]

where

[TABLE]

and

[TABLE]

are solutions to the elastic system $\mathcal{L}_{{\lambda},{\mu}}\mathbf{u}=0$ in $\mathbb{R}^{d}$ excluding the origin.

Since $\{\mathcal{T}_{n}^{m},\mathcal{M}_{n}^{m},\mathcal{N}_{n}^{m}\}$ forms a complete system, we can suppose that the background field $\mathbf{H}$ is represented by

[TABLE]

Then the displacement field $\mathbf{u}$ can be written as the following form:

[TABLE]

By using transmission conditions on each layer and making use of the orthogonality of the vectorial spherical harmonic functions, one can separate each term above and form different transmission equations. In order to ease the exposition, we next mainly consider the case that the background field contains only one type of modes with $n=1$ and $m=0$; see Remark 2.2 for more related discussion about this point.

Suppose the background field is uniformly distributed in $\mathbb{R}^{3}$, i.e. the field $\mathbf{H}$ can be represented as

[TABLE]

Based on the symmetric properties of the multi-layer structure, we assume that the displacement field $\mathbf{u}$ takes the form

[TABLE]

The solution $\mathbf{u}$ satisfies the transmission conditions

[TABLE]

for $j=1,2,\ldots,N$, where the conormtal derivative (or traction) on $S_{j}$ is defined by

[TABLE]

where $\bm{\nu}$ is the exterior unit normal vector to $S_{j}$ and

[TABLE]

for an arbitrary function $v$. In what follows, we define

[TABLE]

Now we give the representation of the perturbed field $\mathbf{u}-\mathbf{H}$ outside the multi-layer structure.

Theorem 2.1**.**

Suppose $\mathbf{u}$ is the solution to the Lamé system (2.1)–(2.7) in $\mathbb{R}^{3}$, with the parameter $\varepsilon_{c}$ given by (2.4). Let $\mathbf{H}$ be given by (2.10). Define the elastic momentum matrix $\mathbb{P}_{N}$ as follows.

[TABLE]

If $\mathbb{P}_{N}$ is invertible, then the Lamé system (2.7) is uniquely solvable with the solution given by the following formula:

[TABLE]

where $\bm{e}:=(1,1,\ldots,1)^{T}$, and

[TABLE]

Remark 2.1**.**

It is clear that the invertibility of the momentum matrix $\mathbb{P}_{N}$ is equivalent to the well-posedness of the Lamé system (2.7). If $\Re\varepsilon>0$, namely $\mathbf{C}_{\tilde{\lambda},\tilde{\mu}}$ is a regular elastic medium, then $\mathcal{L}_{\tilde{\lambda},\tilde{\mu}}$ is an elliptic PDO (Partial Differential Operator), which guarantees the well-posedness of (2.7) and hence the invertibility of $\mathbb{P}_{N}$. On the other hand, if for some layers the elastic materials are negative, the ellipticity of $\mathbb{P}_{N}$ might be broken. Nevertheless, if for those negative materials, say $\varepsilon_{j}$, if one has $\Im\varepsilon_{j}>0$, then the ellipticity of $\mathcal{L}_{\tilde{\lambda},\tilde{\mu}}$ is retained and hence the invertibility of $\mathbb{P}_{N}$. This can be seen in our subsequent analysis.

Remark 2.2**.**

We remark that we only consider the uniform background field case to derive the exact perturbation formula (2.15). However, even if the background field contains more modes, one should be able to derive the related perturbation formula by using the form (2.9). For an illustration, if $\mathbf{H}=\mathfrak{t}_{n,m}^{0}\mathcal{T}_{n}^{m}$ for some specific $m$ and $n$, then the displacement field $\mathbf{u}$ can be written by

[TABLE]

By using a similar strategy as the proof of Theorem 2.1 in what follows, one can derive the exact formula for the coefficient ${\mathfrak{\tilde{t}}}_{n,m}^{0}$. If the background field contain other modes, say e.g. $\mathcal{M}_{n}^{m}$, the representation formula can be more complicated, but it will not bring any essential difficulty. For the sake of simplicity and to better express our main idea, we only study the case for a uniformly distributed background field.

2.3. Proof of Theorem 2.1

By using the transmission conditions (2.11) across the interface $S_{j}$, $j=1,2,\ldots N$, we can deduce the following equations

[TABLE]

where we set $b_{N}=0$. By using (2.13) and some basic arrangements to the equations (2.17), we obtain that

[TABLE]

The matrix $\mathbb{P}_{N}$ is invertible in the case that the elastic material parameters $\varepsilon_{j}$, $j=1,2,\ldots N$, are all positive, and then there holds that

[TABLE]

where $\bm{a}:=(a_{1},a_{2},\ldots,a_{N})^{T}$ and $\Xi$ is defined by

[TABLE]

Furthermore, by the first equation in (2.17) and $b_{N}=0$, we have

[TABLE]

Combining (2.18) with (2.20), there holds

[TABLE]

where $\bm{b}:=(b_{0},b_{1},\ldots,b_{N-1})^{T}$. By extracting the first element in $\bm{b}$, we derive that

[TABLE]

The proof is complete.

3. An algebraic framework for polariton resonances

In this section, we shall establish an algebraic framework for polariton resonances. We mention that in order to analyze the polariton resonance phenomena, one should consider the situation that the elastic material is not regular, that is (2.3) does not hold any more.

3.1. Polariton resonance

It is known that polariton resonance is usually relevant to some eigenvalue problem arising from the PDE systems. We shall also consider the related eigenvalue problem associated with elastostatics system for multi-layer structure. According to our earlier discussion in Remark 2.1, if the parameters $\varepsilon_{j}$, $j=1,2,\ldots,N$, are all positive real-valued, the elastostatic system (2.1)–(2.7) has only trivial solution if $\mathbf{H}=0$. We seek non-trivial solutions to (2.1)–(2.7) when the parameters are allowed to be negative valued, i.e., $\varepsilon_{j}<0$ for some $j=1,2,\ldots,N$. In order to simplify the analysis, in our subsequent study, we always assume that

[TABLE]

and

[TABLE]

where $\varepsilon^{*}$ is a positive number to be chosen and $\delta>0$ is sufficiently small, signifying a lossy parameter. $\mathrm{i}=\sqrt{-1}$. Define

[TABLE]

Then one can readily see that

[TABLE]

for $j=1,2,\ldots,N.$ Hence, the perturbed field (2.15) can be rewritten by

[TABLE]

where $\tilde{\bm{e}}:=(1,-1,1,\dots,(-1)^{N-1})^{T}$ and the matrix $\mathbb{K}_{N}$ is given by

[TABLE]

Next, we give the definition of the polariton resonance for our subsequent study.

Definition 3.1**.**

Consider the system (2.7) associated with the $N$-layer structure $D$, where the elastic material configuration is described in (3.1). Polariton resonance occurs if there exists elastic material configuration such that, in the limiting case as the loss parameter $\delta$ goes to zero, the corresponding PDO $\mathcal{L}_{\tilde{\lambda},\tilde{\mu}}$ possesses non-trivial kernel which in turn induces resonance.

Theorem 3.1**.**

Suppose $\mathbf{u}$ is the solution to the Lamé system (2.1)–(2.7) in $\mathbb{R}^{3}$, with the parameter $\varepsilon_{c}$ be given in (3.1). If there holds

[TABLE]

where $\beta^{*}$ is an eigenvalue of $\mathbb{K}_{N}$, then the polariton resonance occurs.

Proof.

According to Definition 3.1, one only need to verify that the corresponding PDO $\mathcal{L}_{\tilde{\lambda},\tilde{\mu}}$ possesses a non-trivial kernel as $\delta\rightarrow 0$. To this end, let us consider the following Lamé system:

[TABLE]

We shall search for the solution of the following form:

[TABLE]

By following a similar treatment as the proof of Theorem 2.1, one can show that

[TABLE]

That is

[TABLE]

where $\Sigma$ and $\bm{a}$ are defined in proof of Theorem 2.1. By letting $\delta\rightarrow 0$, one has $|\beta I-\mathbb{K}_{N}^{T}|=0$ and thus the homogeneous linear equations (3.6) has non-trivial solutions. Note that $\Sigma^{-1}\bm{a}=(a_{1},a_{2}-a_{1},\ldots,a_{N}-a_{N-1})$, one immediately has that $a_{j}$, $j=1,\ldots,N$ are not identically zero and thus the system (3.5) has non-trivial solutions. $\Box$

In the following, we shall find out the suitable elastic material configurations to make the PDO $\mathcal{L}_{\tilde{\lambda},\tilde{\mu}}$ possess a non-trivial kernel by establishing an algebraic framework.

3.2. Resonant fields and energy blowup

For $\mathbf{u},\mathbf{v}\in H^{1}_{\text{loc}}(\mathbb{R}^{3})^{3}$ , we introduce a bilinear form associated with the Lamé parameters $\lambda,\mu$ by

[TABLE]

where and also in what follows, $\mathbf{A}:\mathbf{B}=\sum_{i,j=1}^{3}a_{ij}b_{ij}$ for two matrices $\mathbf{A}=(a_{ij})_{i,j=1}^{3}$ and $\mathbf{B}=(b_{ij})_{i,j=1}^{3}$. The corresponding quadratic form is defined by

[TABLE]

The resonant field demonstrates an energy blowup. Next, similar to [14, 22], we give the formal definition of the polariton resonance.

Definition 3.2**.**

Consider the system (2.7) associated with the $N$-layer structure $D$, where the elastic material configuration is described in (3.1). Then polariton resonance occurs if the following condition is fulfilled:

[TABLE]

where $\mathbf{E}_{\mathbb{R}^{3}\setminus\overline{D}}^{\lambda,\mu}(\mathbf{u}-\mathbf{H})$ is given by (3.7)–(3.8) with the integration domain replaced by $\mathbb{R}^{3}\setminus\overline{D}$.

In view of Definition 3.2, the polariton resonance mode is parallel to the elastic material configurations which make the parameter $\beta$ the eigenvalue of $\mathbb{K}_{N}$ as $\delta\rightarrow 0$. For this reason, we shall focus on the explicit computation of eigenvalues to the matrix $\mathbb{K}_{N}$, or the exact formula to the roots of $|\mathbb{P}_{N}|$.

3.3. An algebraic framework

In this subsection, we shall derive a precise connection between the polariton resonance and the choice of parameters $\varepsilon_{c}$ in the multi-layer structure $D$. By observing the explicit formula of $|\mathbb{P}_{N}(\beta)|$, $4\leqslant N\leqslant 6$ (see Subsection 4.1 below for details), we give the delicate result about the explicit formula of $|\mathbb{P}_{N}(\beta)|$ for all $N\in\mathbb{N}$ in the following theorem. For the convenience of description, we denote by $\mathbf{i}_{m}$ the multi-index $(i_{1},i_{2},\ldots,i_{m})$ and $\tau_{\mathbf{i}_{m}}$ its sign given by

[TABLE]

Moreover, we denote by $C^{i,m}_{n}$ the set of all combinations of $m$ out $n$, $m\leqslant n$, say e.g., for one combination

[TABLE]

we order them in ascending way, that is, $i_{1}<i_{2}<\cdots<i_{m}$.

Theorem 3.2**.**

Define

[TABLE]

Then for $N\in\mathbb{N}$, it holds that

[TABLE]

Remark 3.1**.**

Theorem 3.2 shows that if the number of layers, $N$, is odd then there exists a zero root $\beta=0$ to the equation $|\mathbb{P}_{N}(\beta)|=0$. The other roots are contained in a quadratic polynomial, whose constant terms are the roots of a polynomial equation of order $\lfloor N/2\rfloor$. Here we denote by $\lfloor t\rfloor$ the integer part of $t\in\mathbb{R}$. In other words, to solve the equation $|\mathbb{P}_{N}(\beta)|=0$, we first solve the equation $|\mathbb{P}_{N}(q)|=0$, which is a polynomial equation of order $\lfloor N/2\rfloor$, with respect to $q$. We then solve the following quadratic equation

[TABLE]

to find the roots of $|\mathbb{P}_{N}(\beta)|=0$.

We shall analyze all the roots to the polynomial equation $|\mathbb{P}_{N}(\beta)|=0$. Note that the explicit formula (3.9) of $|\mathbb{P}_{N}(\beta)|$ can also be written by the following compact form:

[TABLE]

As already mentioned before, in order to ensure that the polariton resonance occurs, it is essential to find the roots of the polynomial

[TABLE]

To this end, we have the following elementary result on the roots of $f_{N}(q)$:

Theorem 3.3**.**

Let $f_{N}(q)$ be defined in (3.11) and $q^{*}$ be the root to $f_{N}(q)$. Then there exists $\lfloor N/2\rfloor$ real values of roots to $f_{N}(q)$. Moreover, we have

[TABLE]

Remark 3.2**.**

1. (i)

It follows from Theorem 3.3 that the roots of $|\mathbb{P}_{N}(\beta)|$ are all real values. In fact, for any real solution $q^{*}$ to (3.11), one can derive that the polynomial equation $|\mathbb{P}_{N}(\beta)|=0$ has two roots,

[TABLE]

this, togrther with (3.12), implies that

[TABLE] 2. (ii)

We also remark that if the polariton resonance occurs, the parameters $\varepsilon_{2j-1}$, $j=1,2,\ldots,N$, defined in (3.1) should have a negative real part, i.e., $\varepsilon^{*}>0$. Indeed, if $\varepsilon^{*}<0$, by the definition of $\beta$ in (3.2) we have that

[TABLE]

Thus in such a case, we cannot have the polariton resonance since the spectrum of the matrix $\mathbb{K}_{N}$ lies in $[-1,\xi]$.

Next, we shall consider that the radius of the layers are extreme large.

Theorem 3.4**.**

Suppose $r_{i}=R+c_{i}$, where $R\gg 1$ and $c_{i}$ are regular constants, $i=1,2,\ldots,N$. Then the polynomial (3.11) can be rewritten as

[TABLE]

Remark 3.3**.**

This indicates that the possible polariton modes for such a set of multi-layer structure are only $\beta=-1,\xi$ if the number of layers $N$ is even and $\beta=0,-1,\xi$ if $N$ is odd.

4. Proofs of Theorems 3.2–3.4

In this section, we study the polariton resonances for the elastostatic system, which in turn casts some light on deriving the explicit formula in Theorem 3.2, and hence facilitates the proofs of Theorems 3.2–3.4.

4.1. Eigenvalue problem

As mentioned above, the polariton mode is parallel to the elastic material configurations which make the parameter $\beta$ the eigenvalue of $\mathbb{K}_{N}$ as $\delta\rightarrow 0$. So, we focus on the eigenvalues of the matrix $\mathbb{K}_{N}$, or the explicit formula of $|\mathbb{P}_{N}|$ defined in (2.14) in this subsection. For this, we first define

[TABLE]

and set $P^{3}_{2}=1$.

Next, we give the recursion formula for $|\mathbb{P}_{N}|$ in the following lemma.

Lemma 4.1**.**

Let $N\geqslant 4$. Then there holds the following recursion formula:

[TABLE]

Proof.

Denote by $\mathbf{E}_{i,j*k}$ the elementary matrix which is transform of identity matrix by multiplying each element of the $j$-th row of the identity matrix by $k$, and then adding it to the $i$-th row. Thus, by using some elementary translation, we can obtain

[TABLE]

where we only changed the first row and column, together with the last row and column. By virtue of the notation (4.1) and the Laplace expansion theorem for determinant one thus has that

[TABLE]

holds for all $N\geqslant 4$, $N\in\mathbb{N}$. $\Box$

Next, we shall study the dependence of the determinant of $P^{i}_{m}$ defined in (4.1), $i\leqslant m$, on $\beta$. By direct computations, we first have that for $m\geqslant 1$,

[TABLE]

Combining the recursion formula (4.2) with (4.5) we obtain that

[TABLE]

In a similar manner, by straightforward (but lengthy and tedious) we have that

[TABLE]

Similarly,

[TABLE]

Observing the above results (4.6)–(4.8) gives us a clue to derive the exact formula in Theorem 3.2. Next, we prove it by mathematical induction.

Proof of Theorem 3.2.

We shall use induction to prove Theorem 3.2. In view of (3.9) we can apply (4.5) to deduce that (3.9) holds for $1\leqslant N\leqslant 3$. Suppose that (3.9) holds for all $N\leqslant N_{0}$, $N_{0}\geqslant 4$, we show that it also holds for $N=N_{0}+1$. Note that $q(\xi-1-\beta)=q(\beta)$, in what follows, we shall denote $q(\beta)$ by $q$ for simplicity.

1. Case i

$N_{0}$ is even. Since (3.9) and (4.2) hold for all $N\leqslant N_{0}$ and the fact that $N=N_{0}+1$ is odd, we obtain

[TABLE]

where we use the notations

[TABLE]

and

[TABLE]

From (4.10), we set

[TABLE]

By direct computations, we deduce that

[TABLE]

where it can be seen that

[TABLE]

and

[TABLE]

where the last equality follows by using the relations

[TABLE]

We can also deduce that

[TABLE]

Furthermore, one can readily obtain that

[TABLE]

for all $k=2,3,\ldots,N_{0}/2-1$. By some proper arrangements to the equation (4.14), we have

[TABLE]

for $k=2,3,\ldots,N_{0}/2-1$. By direct computations, there holds the following relations:

[TABLE]

and

[TABLE]

By substituting (4.16)–(4.17) into (4.15) and straight forward computations, we can deduce that

[TABLE]

which competes the proof for the case that $N_{0}$ is even. 2. Case ii

$N_{0}$ is odd. The proof follows from a similar argument to the case that $N_{0}$ is even. We shall only briefly sketch it. First, it can be derived that

[TABLE]

Similarly to (4.11)–(4.13), we can obtain that $\mathrm{F}_{0}=1$ and

[TABLE]

We can also deduce that $\mathrm{F}_{k}$, $k=2,3,\ldots,N/2-1$, have the same form with (4.15), thus there holds

[TABLE]

The proof is complete. $\Box$

Next, we proceed with the proof of Theorem 3.3.

Proof of Theorem 3.3.

Arguing by contradiction, suppose that

[TABLE]

Combining (4.18) with the fact that $\beta^{2}-(\xi-1)\beta=q^{*}$, we have

[TABLE]

Since $\beta$ satisfies the polynomial equation $|\mathbb{P}_{N}(\beta)|=0$, we can obtain that there exists non-trivial solution to the following linear equations

[TABLE]

where $\mathbf{y}\in\mathbb{R}^{N}$ and $\mathbf{y}\neq\mathbf{0}$. Next, for $\mathbf{x}=(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3})\in\mathbb{R}^{3}$, we define

[TABLE]

where $\bm{a}=(a_{1},a_{2},\ldots,a_{N})^{T}$ and $\bm{b}=(b_{0},b_{1},\ldots,b_{N-1})^{T}$ are determined by

[TABLE]

where $\Upsilon_{N}$ is defined in (2.16) and $\Xi$ is defined by (2.19). Then it can be verified that $\mathbf{u}$ defined in (4.20) is the solution to

[TABLE]

where

[TABLE]

and the parameter $\varepsilon_{j}$ is given by

[TABLE]

Now we claim that $\mathbf{u}$ is not identically zero vector in $\mathbb{R}^{3}$. Indeed, if $\mathbf{u}\equiv\mathbf{0}$, it then follows from (4.20) that

[TABLE]

where

[TABLE]

Furthermore, by using the transmission conditions across the boundary of each layer, we have that

[TABLE]

Combining (4.24) with (4.26), we obtain that

[TABLE]

which implies that $\bm{a}=\mathbf{0}$. This is a contradiction since $\bm{a}=\Xi\mathbf{y}$ and $\mathbf{y}\neq\mathbf{0}$. Thus the claim is validated.

On the one hand, by using Green’s formula [2] for Lamé operator, it holds that

[TABLE]

where the last equality follows by using the transmission conditions across the interfaces $S_{j}$, $j=1,2,\ldots,N$. On the other hand, by using (4.23) it also holds that

[TABLE]

We next distinguish the proof into two cases.

• Case i

$\beta\not\in[-1,\xi]$. Then one can readily see that

[TABLE]

which means that the elastic material parameter $\varepsilon_{j},$ $j=1,2,\ldots,N$, are all positive valued. Then by the well-posedness of the elastostatic system (4.21), one must have that $\mathbf{u}\equiv\mathbf{0}$ in $\mathbb{R}^{3}$. Indeed, in view of (4.29), and using that $\nabla^{s}\mathbf{u}(\mathbf{x}):{\mathbf{C}}_{\lambda,\mu}\nabla^{s}\mathbf{u}(\mathbf{x})\geqslant 0$, we have that $\nabla^{s}\mathbf{u}(\mathbf{x}):{\mathbf{C}}_{\lambda,\mu}\nabla^{s}\mathbf{u}(\mathbf{x})=0$ in $\mathbb{R}^{3}$. Thus, from [2], one then has

[TABLE]

this, together with the decay behavior of $\mathbf{u}$ at infinity, implies that $\mathbf{u}\equiv\mathbf{0}$ in $\mathbb{R}^{3}$.

• Case ii

$\Im{\beta}\neq 0$. It then follows from (4.29) that

[TABLE]

Thus we have

[TABLE]

By using Green’s formula [2] for Lamé operator in $D_{2j-1}$ ,$j=1,2,\ldots,\lfloor(N+1)/2\rfloor$, and the transmission conditions across the interface $S_{{2j}}$, $j=1,2,\ldots,\lfloor N/2\rfloor$, one can deduce that

[TABLE]

Thus we also have

[TABLE]

Combining (4.31) with (4.32), we obtain that $\nabla^{s}\mathbf{u}(\mathbf{x}):{\mathbf{C}}_{\lambda,\mu}\nabla^{s}\mathbf{u}(\mathbf{x})=0$ in $\mathbb{R}^{3}$. Thus $\mathbf{u}\equiv\mathbf{0}$ in $\mathbb{R}^{3}$.

We have shown that $\mathbf{u}\equiv\mathbf{0}$ holds for either $\beta\not\in[-1,\xi]$ or $\Im{\beta}\neq 0$, which contradicts with our assumptions. The proof is complete. $\Box$

4.2. Extreme case

In this subsection, we consider that the radius of the layers are extreme large. It is worth mentioning that the structure of the Earth or other celestial bodies can be treated in this case. To describe this case mathematically, we set $r_{i}=R+c_{i}$, where $R\gg 1$ and $c_{i}$ are regular constants, $i=1,2,\ldots,N$. Then we have

[TABLE]

Let

[TABLE]

By direct computations, we have

[TABLE]

and

[TABLE]

Moreover, we can also observe that

[TABLE]

and

[TABLE]

Next, we give some recursion formulae for $h_{N,k}$ in the following lemma.

Lemma 4.2**.**

Let $h_{N,k}$ be given by (4.34). Then there holds

[TABLE]

and

[TABLE]

Proof.

First, by using (4.35) and (4.36) we have

[TABLE]

which verifies (4.37). Next, by using again (4.35) and (4.36), we obtain

[TABLE]

where the last equality follows the fact that

[TABLE]

Then the second equation in (4.38) holds. In order to show the first equation in (4.38), we shall make use of induction. It can be readily verified that the first equation in (4.38) holds for $N\leqslant 4$. Suppose it holds for $N\leqslant N_{0}-1$, $N_{0}\geqslant 5$, we show that it also holds for $N=N_{0}$. Combining (4.35)–(4.37) with the first equation in (4.38), one finally has

[TABLE]

which completes the proof. $\Box$

In what follows, we begin with the proof of Theorem 3.4.

proof of Theorem 3.4.

By using (3.11) and (4.33), the associated polynomial $f_{N}(q)$ is then given by

[TABLE]

where $h_{N,2k}$ can be calculated recursively by using Lemma 4.2. By using elementary combination theory, we can deduce that

[TABLE]

where $C_{\lfloor N/2\rfloor}^{k}$ denotes the number of combinations for $k$ out of $\lfloor N/2\rfloor$. It then follows from (4.39) and (4.40) that

[TABLE]

The proof is complete. $\Box$

5. Two dimensional case

In this section, for the sake of completeness, we consider the elastostatic scattering problem in the two dimensional space $\mathbb{R}^{2}$.

5.1. Representation of the perturbed field

Let us keep the notations for the multi-layer structure $D$. Suppose the background field $\mathbf{H}$ is represented as

[TABLE]

Then the displacement field $\mathbf{u}$ takes the following form

[TABLE]

By using the transmission conditions across the interface $\{|x|=r_{j}\}$, $j=1,2,\ldots N$, we deduce

[TABLE]

where we set $\tilde{b}_{N}=0$. By similar argument as in the proof of Theorem 2.1, we can obtain the following result:

Theorem 5.1**.**

Suppose $\mathbf{u}$ is the solution to the Lamé system

[TABLE]

where the Lamé parameters $\tilde{\lambda}$ and $\tilde{\mu}$ are given by

[TABLE]

Let $\mathbf{H}$ be given by (5.1) and

[TABLE]

Then we have

[TABLE]

where $\bm{e}:=(1,1,\ldots,1)^{T}$, the matrix $\tilde{\mathbb{P}}_{N}$ and $\tilde{\Upsilon}_{N}$ are defined by

[TABLE]

and

[TABLE]

5.2. Polariton resonance modes

In this subsection, we study the polariton resonances for the elastostatic system in $\mathbb{R}^{2}$. As in the 3D case, we shall introduce a similar structural design, i.e., we assume that $\varepsilon_{j}$, $j=1,2,\ldots,N$, is given by (3.1). Then one can readily see that the matrix $\tilde{P}_{N}$ in (3.6) is reduced to

[TABLE]

where

[TABLE]

For such a setup, all possible polariton resonance modes are contained in the solution to

[TABLE]

which is equivalent to

[TABLE]

By similar argument as in the proof of Theorem 3.3, we can obtain that

Theorem 5.2**.**

Let $\tilde{\mathbb{P}}_{N}(\tilde{\beta})$ be defined in (5.10) and $\tilde{\beta}^{*}$ be the root to (5.11). Then there exists $N$ real values of roots to (5.11). Moreover, we have

[TABLE]

6. Numerical examples

In this section, we present some numerical examples to corroborate our theoretical findings in the previous sections. From a theoretical point of view (see Theorem 3.2), in order to find the polariton modes for any multi-layer structure, it suffices to find the roots of the characteristic polynomial (3.9). Moreover, from Definition 3.2, to show the polariton resonance phenomena, one only needs to observe the elastic energy outside the structure, determined by the term

[TABLE]

We next investigate these two scenarios numerically.

6.1. Roots to characteristic polynomials

First, we consider the radius of layers are equidistance. Let $\varepsilon_{0}=1$. For $N$-layer structure, set

[TABLE]

TABLE 1 shows the roots to the characteristic polynomial $f_{19}(q)$ with $\xi=3$ and $\xi=0.8$, respectively. In FIGURE 2, we plot the values of the polynomial $f_{19}(q)$ with $\xi=3$ in the span $[2.59,3.0]$ (approximately the span between the fourth eigenvalue and last eigenvalue) and the values of the polynomial $f_{19}(q)$ with $\xi=0.8$ in the span $[0.13,0.8]$ (approximately the span between the first eigenvalue and last eigenvalue), it is interesting that the values of the characteristic polynomial in the span are all very close to zero (of the order $10^{-5}$). These spans may have great potential applications such as wave guide for seismic wave.

Next, we consider the radius of layers are decreasing with the same scale $s$, that is

[TABLE]

Let $r_{1}=1$ and $s=0.7$. TABLE 2 presents all the roots to the characteristic polynomial $f_{19}(q)$ with $\xi=3$ and $\xi=0.8$, respectively. Similarly, in FIGURE 3, one can also find out that the values of the polynomial $f_{19}(q)$ with $\xi=3$ in the span $[1.53,2.27]$ ( approximately the span between the fourth eigenvalue and last eigenvalue) and the values of the polynomial $f_{19}(q)$ with $\xi=0.8$ in the span $[0.05,0.61]$ (approximately the span between the first eigenvalue and last eigenvalue), are all very close to zero. Besides, it is worth mentioning that in both set up of structures, the roots $q$ are all positive values. Similar results can be found in FIGURE 5 and 5 for $N=13$ and $N=16$, respectively.

6.2. Illustration of polariton resonances

In (3.1), $\varepsilon_{i}$ are replaced by

[TABLE]

We define the elastic moment tensor here by

[TABLE]

The functionality of $r_{1}^{-3}$ appearing in (6.4) is to reduce the scale of the structure. In FIGURE 6, we show the norm of the elastic moment tensor defined in (6.4) with the multi-layer structure designed by (6.1), where $\xi=3$ and $N=14$. It can be seen that the peaks of the norm of the elastic moment tensor are in accordance with the polariton modes in the setup (6.3).

7. Concluding remarks

In this paper, we studied the elastostatic scattering from a rather general multi-layer metamaterial structures and derived the exact scattering field in terms of the elastic momentum matrix. By highly intricate and delicate analysis of the momentum matrix, we established a handy algebraic framework for analysing the polariton resonances associated with such material structures which yields explicit relationships between the polariton resonance and the geometric and material configurations of the structure. This facilitates the design of metamaterial structures to induce customised resonances. Our study opens up an intriguing direction of further development with many possible extensions as well as applications. In practice, the multi-layer structure can serve as the building block for various material devices. For example, it may serve for constructing elastic/seismic waveguide with enlarged SPR-like (Surface Polariton Resonance) band by employing more layers of metamaterials. As also discussed in the introduction, it can be used for the quantitative design of ESC-vanishing structures to produce elastic invisibility cloaking devices with enhanced cloaking effects with a rigorous basis. Moreover, our study can be extended to the other physical systems including the electromagnetism and acoustics. We shall study these and other extensions in our forthcoming works.

Acknowledgement

The work of Y. Deng was supported by NSFC-RGC Joint Research Grant No. 12161160314 and NSF grant of China No. 11971487. The work of H. Liu was supported by NSFC/RGC Joint Research Scheme, N_CityU101/21, ANR/RGC Joint Research Scheme, A-CityU203/19, and the Hong Kong RGC General Research Funds (projects 12302919, 12301420 and 11300821).

Data availability statement

This is a mathematical paper containing all the necessary theoretical proofs. There are no data to be reported concerning this work.