# Surface acoustic waves in laterally periodic superlattices

A L Shuvalov

► To cite this version:

A L Shuvalov. Surface acoustic waves in laterally periodic superlattices. Wave Motion, 2024, 129, 10.1016/j.wavemoti.2024.103331. hal-04611203

HAL Id: hal-04611203

https://hal.science/hal-04611203

Submitted on 13 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Surface acoustic waves in laterally periodic superlattices

A L Shuvalov

► To cite this version:

A L Shuvalov. Surface acoustic waves in laterally periodic superlattices. Wave Motion, 2024, 129, 10.1016/j.wavemoti.2024.103331. hal-04611203

HAL Id: hal-04611203

https://hal.science/hal-04611203v1

Submitted on 13 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Surface acoustic waves in laterally periodic superlattices

A.L. Shuvalov

Univ. Bordeaux, I2M, UMR CNRS 5295, 351 cours de la Libération, Talence 33405, France

## ARTICLE INFO

**Keywords:**

Surface acoustic waves

Superlattice

Lateral periodicity

Plane wave expansion

Weak modulation

## ABSTRACT

The existence of surface acoustic waves (SAWs) is studied in laterally periodic superlattices, modelled as an anisotropic elastic half-space with an arbitrary periodic variation of its material properties along the stratification direction (call it X) parallel to the surface. Unlike a homogeneous half-space, such a structure allows for more than one (dispersive) SAW. Specifically, it is shown that any superlattice with a generic shape of periodicity profile admits at most three SAW dispersion branches $\omega(K_x)$, i.e., at most three different SAW frequencies at any fixed Bloch wavenumber $K_x$. Moreover, the total number of SAWs at fixed $K_x$ in a pair of superlattices with periodicity profiles obtained from one another by the inversion of the axis X cannot exceed three either. At least one SAW branch must exist in one of these two superlattices unless the bulk-wave threshold is the so-called exceptional (i.e., admits surface skimming wave). The SAW branch is unique in the particular case of a superlattice invariant to the inversion $X \to -X$. The above general results are illustrated by the perturbation theory derivations for the weakly modulated superlattices. Explicit leading-order formulas are obtained for the quasi-Rayleigh wave branch evolving from the Rayleigh wave in each of the mutually "inverse" superlattices and for the quasibulk wave branch evolving from the exceptional bulk-wave threshold in one of the superlattices.

## 1. Introduction

There is a growing interest in the laterally periodic superlattices, i.e. 1D-periodic multilayered substrates stratified normally to the surface. Such structures are highly promising for further increasing the generation frequency of hypersonic surface acoustic waves that lie at the heart of many modern applications ranging from spectroscopy for quality control to the products of information and communication technologies, see e.g. [1–4]. While interdigital or optoacoustic transducers based on depositing a surface grating allow the excitation of coherent surface waves with a frequency of tens GHz, reaching higher frequencies requires its period to be shortened to the nanometer scale, and this is about the limit for the surface patterning techniques. However, such a small period length is achievable by fabricating a 1D-periodic stack of “very long” layers, which can then be cleaved normally to the layer interfaces to produce a desirable laterally periodic substrate [5].

A fundamental issue in the general theory of surface waves is the possibility of a priori prediction of their existence and number regardless of the material and geometrical parameters of the medium. It is well-known that an anisotropic elastically homogeneous half-space (single-crystal substrate) can support only one subsonic surface (Rayleigh) wave, as established by the Barnett–Lothe theorems based on the Stroh formalism and the impedance matrix method [6–9]. In the case of a transversely periodic half-space (a superlattice with layering parallel to the surface), the above approach complemented by the transfer-matrix concept has led [10] to deduce that there may be at most three surface wave branches per a stopband of the dispersion spectrum. The present study addresses the subsonic surface waves in laterally periodic superlattices. This case is somewhat less amenable to the analysis and, therefore, less

<footnote>E-mail address: alexander.shuvalov@u-bordeaux.fr.</footnote>
<figure><img src="image_2.png" /></figure>

# Surface acoustic waves in laterally periodic superlattices

A.L. Shuvalov

Univ. Bordeaux, I2M, UMR CNRS 5295, 351 cours de la Libération, Talence 33405, France

## ARTICLE INFO

**Keywords:**

Surface acoustic waves

Superlattice

Lateral periodicity

Plane wave expansion

Weak modulation

## ABSTRACT

The existence of surface acoustic waves (SAWs) is studied in laterally periodic superlattices, modelled as an anisotropic elastic half-space with an arbitrary periodic variation of its material properties along the stratification direction (call it X) parallel to the surface. Unlike a homogeneous half-space, such a structure allows for more than one (dispersive) SAW. Specifically, it is shown that any superlattice with a generic shape of periodicity profile admits at most three SAW dispersion branches $\omega(K_x)$, i.e., at most three different SAW frequencies at any fixed Bloch wavenumber $K_x$. Moreover, the total number of SAWs at fixed $K_x$ in a pair of superlattices with periodicity profiles obtained from one another by the inversion of the axis X cannot exceed three either. At least one SAW branch must exist in one of these two superlattices unless the bulk-wave threshold is the so-called exceptional (i.e., admits surface skimming wave). The SAW branch is unique in the particular case of a superlattice invariant to the inversion $X \to -X$. The above general results are illustrated by the perturbation theory derivations for the weakly modulated superlattices. Explicit leading-order formulas are obtained for the quasi-Rayleigh wave branch evolving from the Rayleigh wave in each of the mutually "inverse" superlattices and for the quasibulk wave branch evolving from the exceptional bulk-wave threshold in one of the superlattices.

## 1. Introduction

There is a growing interest in the laterally periodic superlattices, i.e. 1D-periodic multilayered substrates stratified normally to the surface. Such structures are highly promising for further increasing the generation frequency of hypersonic surface acoustic waves that lie at the heart of many modern applications ranging from spectroscopy for quality control to the products of information and communication technologies, see e.g. [1–4]. While interdigital or optoacoustic transducers based on depositing a surface grating allow the excitation of coherent surface waves with a frequency of tens GHz, reaching higher frequencies requires its period to be shortened to the nanometer scale, and this is about the limit for the surface patterning techniques. However, such a small period length is achievable by fabricating a 1D-periodic stack of “very long” layers, which can then be cleaved normally to the layer interfaces to produce a desirable laterally periodic substrate [5].

A fundamental issue in the general theory of surface waves is the possibility of a priori prediction of their existence and number regardless of the material and geometrical parameters of the medium. It is well-known that an anisotropic elastically homogeneous half-space (single-crystal substrate) can support only one subsonic surface (Rayleigh) wave, as established by the Barnett–Lothe theorems based on the Stroh formalism and the impedance matrix method [6–9]. In the case of a transversely periodic half-space (a superlattice with layering parallel to the surface), the above approach complemented by the transfer-matrix concept has led [10] to deduce that there may be at most three surface wave branches per a stopband of the dispersion spectrum. The present study addresses the subsonic surface waves in laterally periodic superlattices. This case is somewhat less amenable to the analysis and, therefore, less

<footnote>E-mail address: alexander.shuvalov@u-bordeaux.fr.</footnote>
covered in the literature. Its standard theoretical treatment engages the plane wave expansion (PWE) operating with the matrices of formally infinite size [11]. Let us mention the pioneering paper [12], where the PWE was used for computing the quasi-Rayleigh dispersion branch in a superlattice of cubic bi-layers, and a relatively recent paper [13], where the PWE was integrated into the Stroh formalism of the boundary-value problem for an arbitrary laterally periodic half-space. However, the insight specifically into the surface wave existence in such superlattices amounts so far to the variational inequality sufficient criterion obtained in [14]. We intend to tackle the problem in the general way by applying the methodology of [9,10,13] to the PWE matrices and then transfer the derivation from the reciprocal space to the real space. By this means, considering a laterally periodic superlattice of unrestricted anisotropy and generic periodicity profile, we show that the maximum number of subsonic surface waves with different frequencies at a fixed Bloch wavenumber is three (but only one if the profile is invariant to spatial inversion) Moreover, it is established that the total number of surface waves that propagate in any superlattice and its "inverted" counterpart obtained by rotating 180° around the normal to the surface cannot exceed three either. The emergence of three surface waves in a pair of mutually "inverse" superlattices is exemplified in the paper under the assumption of weak laterally periodic modulation. Within this framework, the Rayleigh-wave gives rise to a surface wave branch in each of the two superlattices, and one more branch may split off the bulk-wave threshold in one of the superlattices but not in the other. Explicit analytical approximations for these branches are derived.

The paper is organised as follows. Part I deals with the general case of a semi-infinite laterally periodic superlattice. Section 2 presents the background by casting the PWE of the wave equation as a Stroh-like boundary-value problem. Section 3 introduces the impedance and admittance matrices in the real space (i.e. as functions of the coordinate) and in the PWE form and infers their algebraic and analytical properties. On this basis, Section 4 proves the statements on the existence and number of surface waves in laterally periodic superlattices. Part II considers the case of weak lateral modulation, for which the obtained analytical asymptotics both illustrate the general results of Part I and have an interest in themselves. Section 5 sets the context, and Section 6 formulates the zero-order PWE background. Section 7 contains the perturbation-theory derivation of the wave packet partial modes in the generic and degenerate cases. Section 8 and Section 9 present the leading-order formulas for the quasi-Rayleigh and quasibulk surface wave dispersion branches. Section 10 provides the conclusions and research outlooks.

Throughout the text, the superscripted symbols $^T$, * and + will indicate transpose, complex and Hermitian conjugation, respectively, and the same notation I will be used for the identity matrix regardless of its size. The large-size (formally infinite) PWE matrices and vectors are labelled by, respectively, big and small tilde or bar overlays, and their 3 × 3 blocks and 3-component subvectors are labelled by, respectively, big and small hat overlays (also, the hat over scalars marks the Fourier coefficients).

# Part I

## Arbitrary lateral periodicity

### 2. PWE in the Stroh-like format

Consider a laterally periodic superlattice understood as an elastic half-space 1D-periodically stratified along the surface, e.g. obtained by cutting an infinite periodic multilayer orthogonally to the layer interfaces. Introduce the spatial variables $x = \mathbf{r} \cdot \mathbf{e}_1$, $y = \mathbf{r} \cdot \mathbf{e}_2$, where $\mathbf{r}$ is the radius vector, and $\mathbf{e}_1$ and $\mathbf{e}_2$ are the unit vectors orthogonal to the layer interfaces and the boundary plane, respectively. The layers are made of homogeneous or functionally graded materials so that the density $\rho$ and the tensor of stiffness coefficients $\mathbf{c} = \{c_{ijkl}\}$, $i, j, k, l = 1, 2, 3$, characterising the superlattice are real T-periodic piecewise constant or piecewise continuous functions of $x$. Given the half-space $y \ge 0$, we shall seek two-dimensional steady surface waves with displacement $\mathbf{u}(x, y)e^{-i\omega t}$ maintaining zero traction at the boundary $y = 0$ and turning to zero at the infinite depth $y \to \infty$. The corresponding governing wave equation and the auxiliary conditions for the surface wave are

$$ \mathbf{t}_{1,x}(x, y) + \mathbf{t}_{2,y}(x, y) = -\omega^2 \rho(x) \mathbf{u}(x, y) \quad (1) $$

and

$$ \mathbf{t}_2(x, 0) = \mathbf{0}, \quad \lim_{y \to \infty} \mathbf{u}(x, y) = \mathbf{0}, \quad (2) $$

where $\mathbf{t}_j(x, y) = \mathbf{e}_j \cdot \sigma(x, y)$ are the traction vectors constituting the stress tensor $\sigma(x, y) = \{c_{ijkl}u_{k,l}\}$ and the comma in the subscript indicates derivative. Hereafter the factor $e^{-i\omega t}$ is omitted.

The PWE approach implies seeking the solution to Eq. (1) in the form

$$ \left( \begin{array}{c} \mathbf{u}(x, y) \\ i\mathbf{t}_2(x, y) \end{array} \right) = \left( \begin{array}{c} \mathbf{A}(x, y) \\ \mathbf{B}(x, y) \end{array} \right) e^{iK_x x} = \sum_{n=-\infty}^{\infty} \left( \begin{array}{c} \hat{\mathbf{A}}^{(n)}(y) \\ \hat{\mathbf{B}}^{(n)}(y) \end{array} \right) e^{ik_n x}, \quad (3) $$

where the vector functions $\mathbf{A}(x, y)$ and $\mathbf{B}(x, y)$ are T-periodic in $x$, $\hat{\mathbf{A}}^{(n)}(y)$ and $\hat{\mathbf{B}}^{(n)}(y)$ are their Fourier coefficients, and $K_x$ is the Bloch wavenumber taken within the 1D Brillouin zone so that

$$ k_n = K_x + gn, \quad |K_x| \le \frac{\pi}{T}, \quad g = \frac{2\pi}{T}, \quad n \in \mathbb{Z}. \quad (4) $$

The parametric dependence of the functions in (3) on $\omega$ and $K_x$ is understood and suppressed. On inserting expansion (3) truncated by the order $n = \pm N$ (i.e. by $2N + 1 \equiv M$ terms) along with the truncated Fourier series of $\rho(x)$ and $c_{ijkl}(x)$,

$$ \rho(x) = \sum_{n=-N}^{N} \hat{\rho}_n e^{i\lg nx}, \quad c_{ijkl}(x) = \sum_{n=-N}^{N} \hat{c}_{ijkl}^{(n)} e^{i\lg nx}, \quad (5) $$
covered in the literature. Its standard theoretical treatment engages the plane wave expansion (PWE) operating with the matrices of formally infinite size [11]. Let us mention the pioneering paper [12], where the PWE was used for computing the quasi-Rayleigh dispersion branch in a superlattice of cubic bi-layers, and a relatively recent paper [13], where the PWE was integrated into the Stroh formalism of the boundary-value problem for an arbitrary laterally periodic half-space. However, the insight specifically into the surface wave existence in such superlattices amounts so far to the variational inequality sufficient criterion obtained in [14]. We intend to tackle the problem in the general way by applying the methodology of [9,10,13] to the PWE matrices and then transfer the derivation from the reciprocal space to the real space. By this means, considering a laterally periodic superlattice of unrestricted anisotropy and generic periodicity profile, we show that the maximum number of subsonic surface waves with different frequencies at a fixed Bloch wavenumber is three (but only one if the profile is invariant to spatial inversion) Moreover, it is established that the total number of surface waves that propagate in any superlattice and its "inverted" counterpart obtained by rotating 180° around the normal to the surface cannot exceed three either. The emergence of three surface waves in a pair of mutually "inverse" superlattices is exemplified in the paper under the assumption of weak laterally periodic modulation. Within this framework, the Rayleigh-wave gives rise to a surface wave branch in each of the two superlattices, and one more branch may split off the bulk-wave threshold in one of the superlattices but not in the other. Explicit analytical approximations for these branches are derived.

The paper is organised as follows. Part I deals with the general case of a semi-infinite laterally periodic superlattice. Section 2 presents the background by casting the PWE of the wave equation as a Stroh-like boundary-value problem. Section 3 introduces the impedance and admittance matrices in the real space (i.e. as functions of the coordinate) and in the PWE form and infers their algebraic and analytical properties. On this basis, Section 4 proves the statements on the existence and number of surface waves in laterally periodic superlattices. Part II considers the case of weak lateral modulation, for which the obtained analytical asymptotics both illustrate the general results of Part I and have an interest in themselves. Section 5 sets the context, and Section 6 formulates the zero-order PWE background. Section 7 contains the perturbation-theory derivation of the wave packet partial modes in the generic and degenerate cases. Section 8 and Section 9 present the leading-order formulas for the quasi-Rayleigh and quasibulk surface wave dispersion branches. Section 10 provides the conclusions and research outlooks.

Throughout the text, the superscripted symbols $^T$, * and + will indicate transpose, complex and Hermitian conjugation, respectively, and the same notation I will be used for the identity matrix regardless of its size. The large-size (formally infinite) PWE matrices and vectors are labelled by, respectively, big and small tilde or bar overlays, and their 3 × 3 blocks and 3-component subvectors are labelled by, respectively, big and small hat overlays (also, the hat over scalars marks the Fourier coefficients).

# Part I

## Arbitrary lateral periodicity

### 2. PWE in the Stroh-like format

Consider a laterally periodic superlattice understood as an elastic half-space 1D-periodically stratified along the surface, e.g. obtained by cutting an infinite periodic multilayer orthogonally to the layer interfaces. Introduce the spatial variables $x = \mathbf{r} \cdot \mathbf{e}_1$, $y = \mathbf{r} \cdot \mathbf{e}_2$, where $\mathbf{r}$ is the radius vector, and $\mathbf{e}_1$ and $\mathbf{e}_2$ are the unit vectors orthogonal to the layer interfaces and the boundary plane, respectively. The layers are made of homogeneous or functionally graded materials so that the density $\rho$ and the tensor of stiffness coefficients $\mathbf{c} = \{c_{ijkl}\}$, $i, j, k, l = 1, 2, 3$, characterising the superlattice are real T-periodic piecewise constant or piecewise continuous functions of $x$. Given the half-space $y \ge 0$, we shall seek two-dimensional steady surface waves with displacement $\mathbf{u}(x, y)e^{-i\omega t}$ maintaining zero traction at the boundary $y = 0$ and turning to zero at the infinite depth $y \to \infty$. The corresponding governing wave equation and the auxiliary conditions for the surface wave are

$$ \mathbf{t}_{1,x}(x, y) + \mathbf{t}_{2,y}(x, y) = -\omega^2 \rho(x) \mathbf{u}(x, y) \quad (1) $$

and

$$ \mathbf{t}_2(x, 0) = \mathbf{0}, \quad \lim_{y \to \infty} \mathbf{u}(x, y) = \mathbf{0}, \quad (2) $$

where $\mathbf{t}_j(x, y) = \mathbf{e}_j \cdot \sigma(x, y)$ are the traction vectors constituting the stress tensor $\sigma(x, y) = \{c_{ijkl}u_{k,l}\}$ and the comma in the subscript indicates derivative. Hereafter the factor $e^{-i\omega t}$ is omitted.

The PWE approach implies seeking the solution to Eq. (1) in the form

$$ \left( \begin{array}{c} \mathbf{u}(x, y) \\ i\mathbf{t}_2(x, y) \end{array} \right) = \left( \begin{array}{c} \mathbf{A}(x, y) \\ \mathbf{B}(x, y) \end{array} \right) e^{iK_x x} = \sum_{n=-\infty}^{\infty} \left( \begin{array}{c} \hat{\mathbf{A}}^{(n)}(y) \\ \hat{\mathbf{B}}^{(n)}(y) \end{array} \right) e^{ik_n x}, \quad (3) $$

where the vector functions $\mathbf{A}(x, y)$ and $\mathbf{B}(x, y)$ are T-periodic in $x$, $\hat{\mathbf{A}}^{(n)}(y)$ and $\hat{\mathbf{B}}^{(n)}(y)$ are their Fourier coefficients, and $K_x$ is the Bloch wavenumber taken within the 1D Brillouin zone so that

$$ k_n = K_x + gn, \quad |K_x| \le \frac{\pi}{T}, \quad g = \frac{2\pi}{T}, \quad n \in \mathbb{Z}. \quad (4) $$

The parametric dependence of the functions in (3) on $\omega$ and $K_x$ is understood and suppressed. On inserting expansion (3) truncated by the order $n = \pm N$ (i.e. by $2N + 1 \equiv M$ terms) along with the truncated Fourier series of $\rho(x)$ and $c_{ijkl}(x)$,

$$ \rho(x) = \sum_{n=-N}^{N} \hat{\rho}_n e^{i\ln x}, \quad c_{ijkl}(x) = \sum_{n=-N}^{N} \hat{c}_{ijkl}^{(n)} e^{i\ln x}, \quad (5) $$
Eq. (1) can be transformed into the system of 6M equations

$$ \frac{d\tilde{\eta}(y)}{dy} = i\tilde{N}\tilde{\eta}(y), \qquad (6) $$

where $\tilde{\eta}(y)$ consists of 3M-component vectors $\tilde{\mathbf{A}}(y)$ and $\tilde{\mathbf{B}}(y)$, each composed of M subvectors $\tilde{\mathbf{A}}^{(n)}(y)$ and $\tilde{\mathbf{B}}^{(n)}(y)$ with n running from $-N$ to $N$, namely,

$$ \tilde{\eta}(y) = \begin{pmatrix} \tilde{\mathbf{A}}(y) \\ \tilde{\mathbf{B}}(y) \end{pmatrix}, \quad \tilde{\mathbf{A}}(y) = \{\tilde{\mathbf{A}}^{(n)}(y)\}, \\ \tilde{\mathbf{B}}(y) = \{\tilde{\mathbf{B}}^{(n)}(y)\}, \qquad (7) $$

and the $6M \times 6M$ Stroh-like matrix of coefficients is [13]

$$ \tilde{\mathbf{N}} \equiv \begin{pmatrix} \tilde{\mathbf{N}}_1 & \tilde{\mathbf{N}}_2 \\ \tilde{\mathbf{N}}_3 & \tilde{\mathbf{N}}_1^+ \end{pmatrix} = \begin{pmatrix} -\tilde{\mathbf{T}}^{-1}\tilde{\mathbf{R}}^+ & -\tilde{\mathbf{T}}^{-1} \\ \tilde{\mathbf{P}} - \tilde{\mathbf{R}}\tilde{\mathbf{T}}^{-1}\tilde{\mathbf{R}}^+ & -\tilde{\mathbf{R}}\tilde{\mathbf{T}}^{-1} \end{pmatrix}. \qquad (8) $$

The $3M \times 3M$ “tiled” submatrices of $\tilde{\mathbf{N}}$ are constructed of the $3 \times 3$ “hatted” blocks as follows

$$ \tilde{\mathbf{P}} = \{\tilde{\mathbf{P}}^{(mn)}\}, \quad \tilde{\mathbf{R}} = \{\tilde{\mathbf{R}}^{(mn)}\}, \quad \tilde{\mathbf{T}} = \{\tilde{\mathbf{T}}^{(mn)}\}, \quad m, n = -N, \dots, N, \\ \tilde{\mathbf{P}}^{(mn)} = k_m k_n \tilde{c}_{i_1 j_1}^{(m-n)} - \omega^2 \hat{\rho}_{m-n} \delta_{ij}, \quad \tilde{\mathbf{R}}^{(mn)}_{ij} = k_m \tilde{c}_{i_1 j_2}^{(m-n)}, \quad \tilde{\mathbf{T}}^{(mn)}_{ij} = \tilde{c}_{i_2 j_2}^{(m-n)}, \qquad (9) $$

where the block indices $m, n$ are superscripted, and the indices of their elements $i, j$ referred to the basis of vectors $e_1$ and $e_2$ are subscripted. Unless otherwise specified, both $\omega$ and $K_x$ are assumed to be real. Then the matrix $\tilde{\mathbf{N}}$ satisfies the symmetry relation

$$ \tilde{\mathbf{N}} = \tilde{\mathbf{T}}\tilde{\mathbf{N}}^+\tilde{\mathbf{T}}, \quad \tilde{\mathbf{T}} = \begin{pmatrix} \tilde{\mathbf{0}} & \mathbf{I} \\ \mathbf{I} & \tilde{\mathbf{0}} \end{pmatrix}, \qquad (10) $$

where $\tilde{\mathbf{0}}$ and $\mathbf{I}$ are the $3M \times 3M$ zero and identity matrices.

In the present case, where the material coefficients are independent of $y$ and hence so is the matrix $\tilde{\mathbf{N}} = \tilde{\mathbf{N}}(\omega, K_x)$, the general solution to Eq. (6) with semisimple (diagonalisable) matrix $\tilde{\mathbf{N}}$ is

$$ \tilde{\eta}(y) = \sum_{J=1}^{6M} c_J \tilde{\xi}_J e^{ik_{J,y}y}, \quad \tilde{\xi}_J = \begin{pmatrix} \tilde{\mathbf{A}}_J \\ \tilde{\mathbf{B}}_J \end{pmatrix}, \quad \tilde{\mathbf{A}}_J = \{\tilde{\mathbf{A}}_J^{(n)}\}, \\ \tilde{\mathbf{B}}_J = \{\tilde{\mathbf{B}}_J^{(n)}\}, \qquad (11) $$

where

$$ (\tilde{\mathbf{N}} - k_{yJ}\mathbf{I})\tilde{\xi}_J = \mathbf{0}, \quad J = 1, \dots, 6M, \qquad (12) $$

and $c_J$ are arbitrary constants. In Part I, we restrict our attention to the frequency interval $0 < \omega < \hat{\omega}(K_x)$, in which all the eigenvalues of $\tilde{\mathbf{N}}$ are complex and hence, due to (11), complex conjugated. Following the canonical terminology of the theory of surface waves in homogeneous materials [8,9], this range is called subsonic, and the threshold frequency $\hat{\omega}(K_x)$ transonic (its “hat” has indeed nothing to do with the Fourier series). The eigenvalues within the subsonic interval are conventionally grouped into two subsets as

$$ k_{yJ} = k_{y,3M+J}^* \quad \text{Im} k_{yJ} > 0, \quad J = 1, \dots, 3M. \qquad (13) $$

Denote the matrix of eigenvectors of $\tilde{\mathbf{N}}$ arranged according to (13) by

$$ \tilde{\Xi} = \|\tilde{\xi}_1 \cdots \tilde{\xi}_{6M}\| = \begin{pmatrix} \|\tilde{\mathbf{A}}_1 \cdots \tilde{\mathbf{A}}_{3M}\| & \|\tilde{\mathbf{A}}_{3M+1} \cdots \tilde{\mathbf{A}}_{6M}\| \\ \|\tilde{\mathbf{B}}_1 \cdots \tilde{\mathbf{B}}_{3M}\| & \|\tilde{\mathbf{B}}_{3M+1} \cdots \tilde{\mathbf{B}}_{6M}\| \end{pmatrix} = \begin{pmatrix} \tilde{\Xi}_1 & \tilde{\Xi}_2 \\ \tilde{\Xi}_3 & \tilde{\Xi}_4 \end{pmatrix}, \qquad (14) $$

where $\|\cdots\|$ implies a matrix consisting of the enclosed vector columns. By virtue of (10) and (13), the matrix $\tilde{\Xi}(\omega, K_x)$ taken in the subsonic domain satisfies the orthonormality relation

$$ \tilde{\Xi}^+ \tilde{\mathbf{T}} \tilde{\Xi} = \tilde{\mathbf{T}} \Leftrightarrow \tilde{\Xi}^T \tilde{\mathbf{T}} \tilde{\Xi}^+ = \tilde{\mathbf{T}}. \qquad (15) $$

In the particular case of a symmetric profile of periodicity, such that is described by the even functions $\rho(x) = \rho(-x)$ and $c_{ijkl}(x) = c_{ijkl}(-x)$ where $x = 0$ is fixed at the period edge or midpoint, the matrix $\tilde{\mathbf{N}}$ is real and hence satisfies $\tilde{\mathbf{N}} = \tilde{\mathbf{T}}\tilde{\mathbf{N}}^T\tilde{\mathbf{T}}$, so identity (15) splits in the following two independent ones

$$ \tilde{\Xi}^T \tilde{\mathbf{T}} = \tilde{\Xi}^*, \quad \tilde{\Xi}^T \tilde{\mathbf{T}} \tilde{\Xi} = \mathbf{I} \Leftrightarrow \tilde{\Xi}^T \tilde{\Xi}^T = \tilde{\mathbf{T}}, \qquad (16) $$

where the former is restricted to the subsonic domain, while the latter is not.

At the transonic threshold $\omega = \hat{\omega}(K_x)$, the (real) eigenvalue degeneracy $k_{y\hat{J}} = k_{y,3M+\hat{J}} \equiv k_{y,\text{deg}}$, $\hat{J} \in \{1, \dots, 3M\}$, renders matrix $\tilde{\mathbf{N}}$ non-semisimple (defective) so that the corresponding pair of eigenvectors $\tilde{\xi}_\hat{J}$ and $\tilde{\xi}_{3M+\hat{J}}$ merges into a single eigenvector $\tilde{\xi}_{\text{deg}} = (\tilde{\mathbf{A}}_{\text{deg}} \tilde{\mathbf{B}}_{\text{deg}})^T$ satisfying $\tilde{\xi}_{\text{deg}}^+ \tilde{\mathbf{T}} \tilde{\xi}_{\text{deg}} = 0$. If accidentally $\tilde{\mathbf{B}}_{\text{deg}} = \mathbf{0}$ (and hence $\mathbf{B}_{\text{deg}}(x, 0) = \mathbf{0}$, i.e. the traction-free boundary condition at $y = 0$ is fulfilled by a single partial mode $\tilde{\eta}(y) = \tilde{\xi}_{\text{deg}} e^{ik_{y,\text{deg}}y}$), then, by analogy with the above terminology [8,9], such transonic state is called exceptional while otherwise it is called normal.
Eq. (1) can be transformed into the system of 6M equations

$$ \frac{d\tilde{\eta}(y)}{dy} = i\tilde{N}\tilde{\eta}(y), \qquad (6) $$

where $\tilde{\eta}(y)$ consists of 3M-component vectors $\tilde{\mathbf{A}}(y)$ and $\tilde{\mathbf{B}}(y)$, each composed of M subvectors $\tilde{\mathbf{A}}^{(n)}(y)$ and $\tilde{\mathbf{B}}^{(n)}(y)$ with $n$ running from $-N$ to $N$, namely,

$$ \tilde{\eta}(y) = \begin{pmatrix} \tilde{\mathbf{A}}(y) \\ \tilde{\mathbf{B}}(y) \end{pmatrix}, \quad \tilde{\mathbf{A}}(y) = \{\tilde{\mathbf{A}}^{(n)}(y)\}, \\ \tilde{\mathbf{B}}(y) = \{\tilde{\mathbf{B}}^{(n)}(y)\}, \qquad (7) $$

and the $6M \times 6M$ Stroh-like matrix of coefficients is [13]

$$ \tilde{\mathbf{N}} \equiv \begin{pmatrix} \tilde{\mathbf{N}}_1 & \tilde{\mathbf{N}}_2 \\ \tilde{\mathbf{N}}_3 & \tilde{\mathbf{N}}_1^+ \end{pmatrix} = \begin{pmatrix} -\tilde{\mathbf{T}}^{-1}\tilde{\mathbf{R}}^+ & -\tilde{\mathbf{T}}^{-1} \\ \tilde{\mathbf{P}} - \tilde{\mathbf{R}}\tilde{\mathbf{T}}^{-1}\tilde{\mathbf{R}}^+ & -\tilde{\mathbf{R}}\tilde{\mathbf{T}}^{-1} \end{pmatrix}. \qquad (8) $$

The $3M \times 3M$ “tiled” submatrices of $\tilde{\mathbf{N}}$ are constructed of the $3 \times 3$ “hatted” blocks as follows

$$ \tilde{\mathbf{P}} = \{\tilde{\mathbf{P}}^{(mn)}\}, \quad \tilde{\mathbf{R}} = \{\tilde{\mathbf{R}}^{(mn)}\}, \quad \tilde{\mathbf{T}} = \{\tilde{\mathbf{T}}^{(mn)}\}, \quad m, n = -N, \dots, N, \\ \tilde{\mathbf{P}}^{(mn)} = k_m k_n \tilde{c}_{i_1 j_1}^{(m-n)} - \omega^2 \hat{\rho}_{m-n} \delta_{ij}, \quad \tilde{\mathbf{R}}^{(mn)}_{ij} = k_m \tilde{c}_{i_1 j_2}^{(m-n)}, \quad \tilde{\mathbf{T}}^{(mn)}_{ij} = \tilde{c}_{i_2 j_2}^{(m-n)}, \qquad (9) $$

where the block indices $m, n$ are superscripted, and the indices of their elements $i, j$ referred to the basis of vectors $\mathbf{e}_1$ and $\mathbf{e}_2$ are subscripted. Unless otherwise specified, both $\omega$ and $K_x$ are assumed to be real. Then the matrix $\tilde{\mathbf{N}}$ satisfies the symmetry relation

$$ \tilde{\mathbf{N}} = \tilde{\mathbf{T}}\tilde{\mathbf{N}}^+\tilde{\mathbf{T}}, \quad \tilde{\mathbf{T}} = \begin{pmatrix} \tilde{\mathbf{0}} & \mathbf{I} \\ \mathbf{I} & \tilde{\mathbf{0}} \end{pmatrix}, \qquad (10) $$

where $\tilde{\mathbf{0}}$ and $\mathbf{I}$ are the $3M \times 3M$ zero and identity matrices.

In the present case, where the material coefficients are independent of $y$ and hence so is the matrix $\tilde{\mathbf{N}} = \tilde{\mathbf{N}}(\omega, K_x)$, the general solution to Eq. (6) with semisimple (diagonalisable) matrix $\tilde{\mathbf{N}}$ is

$$ \tilde{\eta}(y) = \sum_{J=1}^{6M} c_J \tilde{\xi}_J e^{ik_{J,y}y}, \quad \tilde{\xi}_J = \begin{pmatrix} \tilde{\mathbf{A}}_J \\ \tilde{\mathbf{B}}_J \end{pmatrix}, \quad \tilde{\mathbf{A}}_J = \{\tilde{\mathbf{A}}_J^{(n)}\}, \\ \tilde{\mathbf{B}}_J = \{\tilde{\mathbf{B}}_J^{(n)}\}, \qquad (11) $$

where

$$ (\tilde{\mathbf{N}} - k_{yJ}\mathbf{I})\tilde{\xi}_J = \mathbf{0}, \quad J = 1, \dots, 6M, \qquad (12) $$

and $c_J$ are arbitrary constants. In Part I, we restrict our attention to the frequency interval $0 < \omega < \hat{\omega}(K_x)$, in which all the eigenvalues of $\tilde{\mathbf{N}}$ are complex and hence, due to (11), complex conjugated. Following the canonical terminology of the theory of surface waves in homogeneous materials [8,9], this range is called subsonic, and the threshold frequency $\hat{\omega}(K_x)$ transonic (its “hat” has indeed nothing to do with the Fourier series). The eigenvalues within the subsonic interval are conventionally grouped into two subsets as

$$ k_{yJ} = k_{y,3M+J}^*, \quad \mathrm{Im} \, k_{yJ} > 0, \quad J = 1, \dots, 3M. \qquad (13) $$

Denote the matrix of eigenvectors of $\tilde{\mathbf{N}}$ arranged according to (13) by

$$ \tilde{\Xi} = \|\tilde{\xi}_1 \cdots \tilde{\xi}_{6M}\| = \begin{pmatrix} \|\tilde{\mathbf{A}}_1 \cdots \tilde{\mathbf{A}}_{3M}\| & \|\tilde{\mathbf{A}}_{3M+1} \cdots \tilde{\mathbf{A}}_{6M}\| \\ \|\tilde{\mathbf{B}}_1 \cdots \tilde{\mathbf{B}}_{3M}\| & \|\tilde{\mathbf{B}}_{3M+1} \cdots \tilde{\mathbf{B}}_{6M}\| \end{pmatrix} = \begin{pmatrix} \tilde{\Xi}_1 & \tilde{\Xi}_2 \\ \tilde{\Xi}_3 & \tilde{\Xi}_4 \end{pmatrix}, \qquad (14) $$

where $\|\cdots\|$ implies a matrix consisting of the enclosed vector columns. By virtue of (10) and (13), the matrix $\tilde{\Xi}(\omega, K_x)$ taken in the subsonic domain satisfies the orthonormality relation

$$ \tilde{\Xi}^+ \tilde{\mathbf{T}} \tilde{\Xi} = \tilde{\mathbf{T}} \Leftrightarrow \tilde{\Xi}^T \tilde{\mathbf{T}} \tilde{\Xi}^+ = \tilde{\mathbf{T}}. \qquad (15) $$

In the particular case of a symmetric profile of periodicity, such that is described by the even functions $\rho(x) = \rho(-x)$ and $c_{ijkl}(x) = c_{ijkl}(-x)$ where $x = 0$ is fixed at the period edge or midpoint, the matrix $\tilde{\mathbf{N}}$ is real and hence satisfies $\tilde{\mathbf{N}} = \tilde{\mathbf{T}}\tilde{\mathbf{N}}^T\tilde{\mathbf{T}}$, so identity (15) splits in the following two independent ones

$$ \tilde{\Xi}^T \tilde{\mathbf{T}} = \tilde{\Xi}^*, \quad \tilde{\Xi}^T \tilde{\mathbf{T}} \tilde{\Xi} = \mathbf{I} \Leftrightarrow \tilde{\Xi}^T \tilde{\Xi}^T = \tilde{\mathbf{T}}, \qquad (16) $$

where the former is restricted to the subsonic domain, while the latter is not.

At the transonic threshold $\omega = \hat{\omega}(K_x)$, the (real) eigenvalue degeneracy $k_{y\hat{J}} = k_{y,3M+\hat{J}} \equiv k_{y,\mathrm{deg}}$, $\hat{J} \in \{1, \dots, 3M\}$, renders matrix $\tilde{\mathbf{N}}$ non-semisimple (defective) so that the corresponding pair of eigenvectors $\tilde{\xi}_\hat{J}$ and $\tilde{\xi}_{3M+\hat{J}}$ merges into a single eigenvector $\tilde{\xi}_{\mathrm{deg}} = (\tilde{\mathbf{A}}_{\mathrm{deg}}\tilde{\mathbf{B}}_{\mathrm{deg}})^T$ satisfying $\tilde{\xi}_{\mathrm{deg}}^+ \tilde{\mathbf{T}} \tilde{\xi}_{\mathrm{deg}} = 0$. If accidentally $\tilde{\mathbf{B}}_{\mathrm{deg}} = \mathbf{0}$ (and hence $\mathbf{B}_{\mathrm{deg}}(x, 0) = \mathbf{0}$, i.e. the traction-free boundary condition at $y = 0$ is fulfilled by a single partial mode $\tilde{\eta}(y) = \tilde{\xi}_{\mathrm{deg}} e^{ik_{y,\mathrm{deg}}y}$), then, by analogy with the above terminology [8,9], such transonic state is called exceptional while otherwise it is called normal.
We are interested in the surface-localised (evanescent) wave field, which falls off as the depth grows. In view of (7), (11) and (13), inserting

$$ \begin{pmatrix} \hat{\mathbf{A}}^{(n)}(y) \\ \hat{\mathbf{B}}^{(n)}(y) \end{pmatrix} = \sum_{J=1}^{3M} c_J \begin{pmatrix} \hat{\mathbf{A}}_{J}^{(n)} \\ \hat{\mathbf{B}}_{J}^{(n)} \end{pmatrix} e^{ik_{yJ}y} \qquad (17) $$

in Eq. (3) yields the PWE expression for such wave field in the half-space $y \ge 0$:

$$ \begin{pmatrix} \mathbf{u}(x, y) \\ it_2(x, y) \end{pmatrix} = e^{iK_x x} \sum_{J=1}^{3M} c_J \begin{pmatrix} \mathbf{A}_J(x) \\ \mathbf{B}_J(x) \end{pmatrix} e^{ik_{yJ}y} = \sum_{J=1}^{3M} \sum_{n=-N}^{N} c_J \begin{pmatrix} \hat{\mathbf{A}}_J^{(n)} \\ \hat{\mathbf{B}}_J^{(n)} \end{pmatrix} e^{i(k_n x + k_{yJ} y)}, \qquad (18) $$

where asymptotic vanishing at $y \to \infty$ is ensured by $\mathrm{Im} \, k_{yJ} > 0$ and the coefficients $c_J$ are to be determined from the boundary condition.

Note that the solution to (1) may be sought beyond PWE in the form $(\mathbf{u} \, it_2)^T = \sum_{J=1}^{6M} \eta_J(x) e^{ik_{yJ}y}$ with each $J$th $\eta_J(x)$ being the partial solution of the system of 6 equations

$$ \frac{d\eta(x)}{dx} = i\mathbf{N}(x; \omega, k_y) \eta(x). \qquad (19) $$

Here, $\mathbf{N}$ is the $6 \times 6$ Stroh matrix built from $3 \times 3$ blocks of the similar pattern as in (8), (9) but with $\rho(x)$ and $c_{ijkl}(x)$ replacing the Fourier coefficients $\hat{\rho}_n$ and $\hat{c}_{ijkl}^{(n)}$, the stiffness indices 1, 2 replacing 2, 1, respectively, and the wavenumber $k_y$ replacing $k_m$, $k_n$. The matricant $\mathbf{M}$ of Eq. (19) may be evaluated, numerically at least, in the following steps: first, its standard expression as a product of matrix exponentials or the Peano series is substituted into the characteristic equation det $[\mathbf{M}(T, 0; \omega, k_y) - e^{iK_x T}\mathbf{I}] = 0$, from which the values $k_{yJ}(K_x)$ with $\mathrm{Im} \, k_{yJ} > 0$ are calculated; then, for each $J = 1, \dots, 3M$, the eigenvector $\mathbf{w}$ of $\mathbf{M}(T, 0; \omega, k_{yJ})$ corresponding to its eigenvalue $e^{iK_x T}$ is found and used as an initial condition in the formula

$$ \eta_J(x) = \mathbf{M}(x, 0; \omega, k_{yJ}) \mathbf{w}. \qquad (20) $$

### 3. Impedance and admittance matrices

#### 3.1. Definitions

Introduce the $3 \times 3$ impedance $\mathbf{Z}(x)$ and admittance $\mathbf{Y}(x) = \mathbf{Z}^{-1}(x)$ as $T$-periodic $3 \times 3$ matrix functions connecting the displacement and traction vectors in (18), namely, $t_2(x, y) = -\mathbf{Z}(x)\mathbf{u}(x, y)$ and so

$$ \mathbf{B}_J(x) = -i\mathbf{Z}(x)\mathbf{A}_J(x), \quad \mathbf{A}_J(x) = i\mathbf{Y}(x)\mathbf{B}_J(x), \quad J = 1, \dots, 3M, \qquad (21) $$

where the dependence on the parameters $\omega$ and $K_x$ is understood. Substituting the (truncated) Fourier series of $\mathbf{B}_J(x)$, $\mathbf{A}_J(x)$ (see (3)) and of $\mathbf{Z}(x)$, $\mathbf{Y}(x)$,

$$ \mathbf{Z}(x) = \sum_{n=-N}^{N} \hat{\mathbf{Z}}^{(n)} e^{i\mathbf{g}nx}, \quad \mathbf{Y}(x) = \sum_{n=-N}^{N} \hat{\mathbf{Y}}^{(n)} e^{i\mathbf{g}nx}, \qquad (22) $$

into (21) yields the relations

$$ \hat{\mathbf{B}}_J^{(n)} = -i\hat{\mathbf{Z}}^{(n-m)}\hat{\mathbf{A}}_J^{(m)}, \quad \hat{\mathbf{A}}_J^{(n)} = i\hat{\mathbf{Y}}^{(n-m)}\hat{\mathbf{B}}_J^{(m)}. \qquad (23) $$

They may be collected into the form

$$ \tilde{\mathbf{B}}_J = -i\tilde{\mathbf{Z}}\tilde{\mathbf{A}}_J, \quad \tilde{\mathbf{A}}_J = i\tilde{\mathbf{Y}}\tilde{\mathbf{B}}_J \qquad (24) $$

with the $3M \times 3M$ block Toëplitz matrices

$$ \tilde{\mathbf{Z}} = \{\hat{\mathbf{Z}}^{(n-m)}\}, \quad \tilde{\mathbf{Y}} = \{\hat{\mathbf{Y}}^{(n-m)}\}, \quad n, m = -N, \dots, N. \qquad (25) $$

Rewriting $3M$ vector Eqs. (24) in the matrix form and invoking notations (14) yields $\tilde{\mathbf{\Xi}}_3 = -i\tilde{\mathbf{Z}}\tilde{\mathbf{\Xi}}_1$ and $\tilde{\mathbf{\Xi}}_1 = i\tilde{\mathbf{Y}}\tilde{\mathbf{\Xi}}_3$, i.e.

$$ \tilde{\mathbf{Z}}(\omega, K_x) = i\tilde{\mathbf{\Xi}}_3\tilde{\mathbf{\Xi}}_1^{-1}, \quad \tilde{\mathbf{Y}}(\omega, K_x) = -i\tilde{\mathbf{\Xi}}_1\tilde{\mathbf{\Xi}}_3^{-1}. \qquad (26) $$

Eq. (26) defines the Fourier coefficients of $\mathbf{Z}(x)$ and $\mathbf{Y}(x)$ via the eigenvectors of $\tilde{\mathbf{N}}$.

Next, consider another laterally periodic half-space whose periodicity profile is related to that of the given one by the inversion of the $X$-axis. Denote its density and stiffness coefficients by $\rho'(x)$ and $c'_{ijkl}(x)$ (the prime does not indeed mean a derivative); then, by definition,

$$ \rho'(x) = \rho(-x), \quad c'_{ijkl}(x) = c_{ijkl}(-x). \qquad (27) $$

The Fourier coefficients of mutually "inverted" profiles (27) are complex conjugated, hence so are their Stroh-like matrices (8), i.e. $\tilde{\mathbf{N}}' = \tilde{\mathbf{N}}^*$ (regarding the entry $\tilde{\mathbf{T}}^{-1}$, it is taken into account that the inverse of a complex conjugated matrix is complex conjugate of the inverse). Thus, the sets of eigenvalues and eigenvectors of $\tilde{\mathbf{N}}$ and $\tilde{\mathbf{N}}'$ are complex conjugates of each other. The former implies that
We are interested in the surface-localised (evanescent) wave field, which falls off as the depth grows. In view of (7), (11) and (13), inserting

$$ \begin{pmatrix} \hat{\mathbf{A}}^{(n)}(y) \\ \hat{\mathbf{B}}^{(n)}(y) \end{pmatrix} = \sum_{J=1}^{3M} c_J \begin{pmatrix} \hat{\mathbf{A}}_{J}^{(n)} \\ \hat{\mathbf{B}}_{J}^{(n)} \end{pmatrix} e^{ik_{yJ}y} \qquad (17) $$

in Eq. (3) yields the PWE expression for such wave field in the half-space $y \ge 0$:

$$ \begin{pmatrix} \mathbf{u}(x, y) \\ it_2(x, y) \end{pmatrix} = e^{iK_x x} \sum_{J=1}^{3M} c_J \begin{pmatrix} \mathbf{A}_J(x) \\ \mathbf{B}_J(x) \end{pmatrix} e^{ik_{yJ}y} = \sum_{J=1}^{3M} \sum_{n=-N}^{N} c_J \begin{pmatrix} \hat{\mathbf{A}}_{J}^{(n)} \\ \hat{\mathbf{B}}_{J}^{(n)} \end{pmatrix} e^{i(k_n x + k_{yJ} y)}, \qquad (18) $$

where asymptotic vanishing at $y \to \infty$ is ensured by $\mathrm{Im} \, k_{yJ} > 0$ and the coefficients $c_J$ are to be determined from the boundary condition.

Note that the solution to (1) may be sought beyond PWE in the form $(\mathbf{u} \, it_2)^T = \sum_{J=1}^{6M} \eta_J(x) e^{ik_{yJ}y}$ with each $J$th $\eta_J(x)$ being the partial solution of the system of 6 equations

$$ \frac{d\eta(x)}{dx} = i\mathbf{N}(x; \omega, k_y) \eta(x). \qquad (19) $$

Here, $\mathbf{N}$ is the $6 \times 6$ Stroh matrix built from $3 \times 3$ blocks of the similar pattern as in (8), (9) but with $\rho(x)$ and $c_{ijkl}(x)$ replacing the Fourier coefficients $\hat{\rho}_n$ and $\hat{c}_{ijkl}^{(n)}$, the stiffness indices 1, 2 replacing 2, 1, respectively, and the wavenumber $k_y$ replacing $k_m$, $k_n$. The matricant $\mathbf{M}$ of Eq. (19) may be evaluated, numerically at least, in the following steps: first, its standard expression as a product of matrix exponentials or the Peano series is substituted into the characteristic equation det $[\mathbf{M}(T, 0; \omega, k_y) - e^{iK_x T}\mathbf{I}] = 0$, from which the values $k_{yJ}(K_x)$ with $\mathrm{Im} \, k_{yJ} > 0$ are calculated; then, for each $J = 1, \dots, 3M$, the eigenvector $\mathbf{w}$ of $\mathbf{M}(T, 0; \omega, k_{yJ})$ corresponding to its eigenvalue $e^{iK_x T}$ is found and used as an initial condition in the formula

$$ \eta_J(x) = \mathbf{M}(x, 0; \omega, k_{yJ}) \mathbf{w}. \qquad (20) $$

# 3. Impedance and admittance matrices

## 3.1. Definitions

Introduce the $3 \times 3$ impedance $\mathbf{Z}(x)$ and admittance $\mathbf{Y}(x) = \mathbf{Z}^{-1}(x)$ as $T$-periodic $3 \times 3$ matrix functions connecting the displacement and traction vectors in (18), namely, $t_2(x, y) = -\mathbf{Z}(x)\mathbf{u}(x, y)$ and so

$$ \mathbf{B}_J(x) = -i\mathbf{Z}(x)\mathbf{A}_J(x), \quad \mathbf{A}_J(x) = i\mathbf{Y}(x)\mathbf{B}_J(x), \quad J = 1, \dots, 3M, \qquad (21) $$

where the dependence on the parameters $\omega$ and $K_x$ is understood. Substituting the (truncated) Fourier series of $\mathbf{B}_J(x)$, $\mathbf{A}_J(x)$ (see (3)) and of $\mathbf{Z}(x)$, $\mathbf{Y}(x)$,

$$ \mathbf{Z}(x) = \sum_{n=-N}^{N} \hat{\mathbf{Z}}^{(n)} e^{i\mathbf{g}nx}, \quad \mathbf{Y}(x) = \sum_{n=-N}^{N} \hat{\mathbf{Y}}^{(n)} e^{i\mathbf{g}nx}, \qquad (22) $$

into (21) yields the relations

$$ \hat{\mathbf{B}}_{J}^{(n)} = -i \hat{\mathbf{Z}}^{(n-m)} \hat{\mathbf{A}}_{J}^{(m)}, \quad \hat{\mathbf{A}}_{J}^{(n)} = i \hat{\mathbf{Y}}^{(n-m)} \hat{\mathbf{B}}_{J}^{(m)}. \qquad (23) $$

They may be collected into the form

$$ \tilde{\mathbf{B}}_J = -i \tilde{\mathbf{Z}} \tilde{\mathbf{A}}_J, \quad \tilde{\mathbf{A}}_J = i \tilde{\mathbf{Y}} \tilde{\mathbf{B}}_J \qquad (24) $$

with the $3M \times 3M$ block Toëplitz matrices

$$ \tilde{\mathbf{Z}} = \{\hat{\mathbf{Z}}^{(n-m)}\}, \quad \tilde{\mathbf{Y}} = \{\hat{\mathbf{Y}}^{(n-m)}\}, \quad n, m = -N, \dots, N. \qquad (25) $$

Rewriting $3M$ vector Eqs. (24) in the matrix form and invoking notations (14) yields $\tilde{\mathbf{\Xi}}_3 = -i\tilde{\mathbf{Z}}\tilde{\mathbf{\Xi}}_1$ and $\tilde{\mathbf{\Xi}}_1 = i\tilde{\mathbf{Y}}\tilde{\mathbf{\Xi}}_3$, i.e.

$$ \tilde{\mathbf{Z}}(\omega, K_x) = i\tilde{\mathbf{\Xi}}_3\tilde{\mathbf{\Xi}}_1^{-1}, \quad \tilde{\mathbf{Y}}(\omega, K_x) = -i\tilde{\mathbf{\Xi}}_1\tilde{\mathbf{\Xi}}_3^{-1}. \qquad (26) $$

Eq. (26) defines the Fourier coefficients of $\mathbf{Z}(x)$ and $\mathbf{Y}(x)$ via the eigenvectors of $\tilde{\mathbf{N}}$.

Next, consider another laterally periodic half-space whose periodicity profile is related to that of the given one by the inversion of the $X$-axis. Denote its density and stiffness coefficients by $\rho'(x)$ and $c'_{ijkl}(x)$ (the prime does not indeed mean a derivative); then, by definition,

$$ \rho'(x) = \rho(-x), \quad c'_{ijkl}(x) = c_{ijkl}(-x). \qquad (27) $$

The Fourier coefficients of mutually "inverted" profiles (27) are complex conjugated, hence so are their Stroh-like matrices (8), i.e. $\tilde{\mathbf{N}}' = \tilde{\mathbf{N}}^*$ (regarding the entry $\tilde{\mathbf{T}}^{-1}$, it is taken into account that the inverse of a complex conjugated matrix is complex conjugate of the inverse). Thus, the sets of eigenvalues and eigenvectors of $\tilde{\mathbf{N}}$ and $\tilde{\mathbf{N}}'$ are complex conjugates of each other. The former implies that
the transonic frequency $\omega$ is identical for both profile configurations; at the same time, once the eigenvalues $k_{y,J}$ and $k_{y,J+3M} = k_{y,J}^*$ of $\tilde{N}$ are arranged according to (13), the subset of eigenvalues $k'_{y,J}$ of $\tilde{N}'$ such that $\text{Im} k'_{y,J} > 0$ satisfies $k'_{y,J} = k_{y,3M+J}^*$ and so the corresponding subset of eigenvectors of $\tilde{N}'$ is $\tilde{\xi}'_J = \tilde{\xi}^*_{3M+J} = (\tilde{A}^*_{3M+J} \tilde{B}^*_{3M+J})^T$, i.e. $\tilde{\Xi}' = \tilde{\Xi}^* \tilde{T}$. In other words, the surface-localised wave field in the half-space $y \ge 0$ with an “inverted” profile may be written as

$$ \begin{pmatrix} \mathbf{u}(x, y) \\ i\mathbf{t}_2(x, y) \end{pmatrix} = e^{iK_x x} \sum_{J=1}^{3M} c'_J \begin{pmatrix} \mathbf{A}'_J(x) \\ \mathbf{B}'_J(x) \end{pmatrix} e^{ik'_{y,J} y} = \sum_{J=1}^{3M} \sum_{n=-N}^{N} c'_J \begin{pmatrix} \hat{\mathbf{A}}_{3M+J}^{(n)*} \\ \hat{\mathbf{B}}_{3M+J}^{(n)*} \end{pmatrix} e^{i(k_x x + k_{y,J} y)} \quad (28) $$

i.e. it is the complex conjugate of the “non-physical” (growing into depth) wave field in the half-space $y \ge 0$ with the “direct” profile. Introduce the impedance and admittance for the half-space $y \ge 0$ with the “inverted” profile by the definition

$$ \mathbf{B}'_J(x) = -i \tilde{Z}^{*}(x) \mathbf{A}'_J(x), \quad \mathbf{A}'_J(x) = i \tilde{Y}^{*}(x) \mathbf{B}'_J(x), \quad (29) $$

where

$$ \tilde{Z}^{*}(x) = \sum_{n=-N}^{N} \hat{Z}^{(n)*} e^{ignx}, \quad \tilde{Y}^{*}(x) = \sum_{n=-N}^{N} \hat{Y}^{(n)*} e^{ignx}, \quad (30) $$

and $\hat{Z}^{(n)}$ and $\hat{Y}^{(n)}$ are the Fourier coefficients of $Z'(x)$ and $Y'(x)$ relating $B'_J(-x)$ and $iA'_J(-x)$. Inserting (28) and (30) in (29) and taking the complex conjugate yields

$$ \hat{\mathbf{B}}_{3M+J}^{(n)} = i \hat{Z}^{(n-m)} \hat{\mathbf{A}}_{3M+J}^{(m)}, \quad \hat{\mathbf{A}}_{3M+J}^{(n)} = -i \hat{Y}^{(n-m)} \hat{\mathbf{B}}_{3M+J}^{(m)}. \quad (31) $$

Setting Eqs. (31) into the $3M \times 3M$ form leads to

$$ \tilde{Z}'(\omega, K_x) = -i\tilde{\Xi}_4\tilde{\Xi}_2^{-1}, \quad \tilde{Y}'(\omega, K_x) = i\tilde{\Xi}_2\tilde{\Xi}_4^{-1}, \quad (32) $$

where the $3M \times 3M$ block Toöplitz matrices $\tilde{Z}'$ and $\tilde{Y}'$ are arranged from $3 \times 3$ blocks $\hat{Z}^{(n-m)}$ and $\hat{Y}^{(n-m)}$. As a test, assume the case of a symmetric periodicity profile such that the "direct" and "inverted" half-spaces $y \ge 0$ are identical. Then, by (16)$_1$, $\tilde{\Xi}_1 = \tilde{\Xi}_2^*$, $\tilde{\Xi}_3 = \tilde{\Xi}_4^*$, hence Eqs. (26) and (32) yield $\tilde{Z} = \tilde{Z}^{*}$ and therefore $Z(x) = Z^{*}(x)$, i.e., the results for both above configurations coincide, as they must.

Note that the PWE version of the Barnett–Lothe integral formalism [13] allows direct evaluation of the impedances bypassing the eigenvalue problem of $\tilde{N}$. This may be done via the sign function of the matrix $i\tilde{N}$, which was shown to satisfy the equalities

$$ \operatorname{sign}i\tilde{N} = i\langle\tilde{N}_{\theta}\rangle = i\tilde{Y}, \quad (33) $$

where $\langle\tilde{N}_{\theta}\rangle = \frac{1}{\pi} \int_{0}^{\pi} \tilde{N}_{\theta} d\theta (= \tilde{T} \langle\tilde{N}_{\theta}\rangle + \tilde{T})$ with the matrix $\tilde{N}_{\theta}$ constructed similarly to (8), (9) but in the rotating frame of vectors $(e_1, e_2)$, and

$$ \tilde{Y} \equiv \begin{pmatrix} \tilde{Y}_1 & \tilde{Y}_2 \\ \tilde{Y}_3 & -\tilde{Y}_1 \end{pmatrix} = i\tilde{\Xi}\operatorname{diag}(\tilde{I}, -\tilde{I})\tilde{\Xi}^{-1} = i \begin{pmatrix} \tilde{I} - 2\tilde{\Xi}_2\tilde{\Xi}_3^+ & 2\tilde{\Xi}_1\tilde{\Xi}_2^+ \\ 2\tilde{\Xi}_3\tilde{\Xi}_4^+ & \tilde{I} - 2\tilde{\Xi}_4\tilde{\Xi}_1^+ \end{pmatrix}. \quad (34) $$

Hence

$$ \begin{aligned} \tilde{Z} &= -\tilde{Y}_2^{-1}(\mathbf{I} + i\tilde{Y}_1), \quad \tilde{Y} = \tilde{Y}_3^{-1}(\mathbf{I} + i\tilde{Y}_1^T); \\ \tilde{Z}' &= -\tilde{Y}_2^{-1}(\mathbf{I} - i\tilde{Y}_1), \quad \tilde{Y}' = \tilde{Y}_3^{-1}(\mathbf{I} - i\tilde{Y}_1^T). \end{aligned} \quad (35) $$

where the blocks of $\tilde{Y}$ may be calculated from the sign function either by integrating the matrix $\tilde{N}_{\theta}$ or else by invoking the corresponding projector functions, see [13].

## 3.2. Properties of the PWE impedance and admittance

According to the on-diagonal blocks of the orthonormality relation (15)$_1$, the impedance matrices (26) and (32) are Hermitian:

$$ \tilde{Z} = \tilde{Z}^+, \quad \tilde{Z}' = \tilde{Z}^{+}, \quad (36) $$

and hence so are their inverse, the admittances matrices. In turn, from the off-diagonal blocks of (15)$_1$ it follows that

$$ \tilde{Z} + \tilde{Z}' = i(\tilde{\Xi}_1\tilde{\Xi}_2^+)^{-1}(= -2\tilde{Y}_2^{-1}), \quad \tilde{Y} + \tilde{Y}' = -i(\tilde{\Xi}_3\tilde{\Xi}_4^+)^{-1}(= 2\tilde{Y}_3^{-1}), \quad (37) $$

where $\tilde{Y}_2 = \tilde{Y}_3^+$ and $\tilde{Y}_3 = \tilde{Y}_4^+$ are blocks of the matrix (34). If the periodicity profile is symmetric, then $\tilde{Z} = \tilde{Z}^{*}$ (see below (32)) and hence (37)$_1$ reduces to $\text{Re}\tilde{Z} = -\tilde{Y}_2^{-1}$.

The sign-definiteness of impedances and admittances in a laterally periodic half-space rests on their general link to energy quantities and may be proved along similar lines as for other half-space models. First, consider the cycle-averaged kinetic energy density $E(x, y) = \frac{1}{4}\rho(x)\omega^2 |\mathbf{u}(x, y)|^2$ for the evanescent wave field (18). Performing the Fourier-series expansion of $E(x, y)$ in $x$ and taking its zero-order coefficient $\hat{E}_0(y)$ yields

$$ \hat{E}_0(y) = \frac{1}{4}\omega^2 \tilde{A}^+(y) \tilde{p} \tilde{A}(y), \quad (38) $$
the transonic frequency $\omega$ is identical for both profile configurations; at the same time, once the eigenvalues $k_{y,J}$ and $k_{y,J+3M} = k_{y,J}^*$ of $\tilde{N}$ are arranged according to (13), the subset of eigenvalues $k'_{y,J}$ of $\tilde{N}'$ such that $\mathrm{Im} k'_{y,J} > 0$ satisfies $k'_{y,J} = k_{y,3M+J}^*$ and so the corresponding subset of eigenvectors of $\tilde{N}'$ is $\tilde{\xi}'_J = \tilde{\xi}^*_{3M+J} = (\tilde{\Lambda}^*_{3M+J} \tilde{\mathbf{B}}^*_{3M+J})^T$, i.e. $\tilde{\Xi}' = \tilde{\Xi}^* \tilde{\Pi}$. In other words, the surface-localised wave field in the half-space $y \ge 0$ with an “inverted” profile may be written as

$$ \begin{pmatrix} \mathbf{u}(x, y) \\ i\mathbf{t}_2(x, y) \end{pmatrix} = e^{iK_x x} \sum_{J=1}^{3M} c'_J \begin{pmatrix} \mathbf{A}'_J(x) \\ \mathbf{B}'_J(x) \end{pmatrix} e^{ik'_{y,J} y} = \sum_{J=1}^{3M} \sum_{n=-N}^{N} c'_J \begin{pmatrix} \hat{\mathbf{A}}_{3M+J}^{(n)*} \\ \hat{\mathbf{B}}_{3M+J}^{(n)*} \end{pmatrix} e_j^{i(k_j x + k_{yJ} y)}, \quad (28) $$

i.e. it is the complex conjugate of the “non-physical” (growing into depth) wave field in the half-space $y \ge 0$ with the “direct” profile. Introduce the impedance and admittance for the half-space $y \ge 0$ with the “inverted” profile by the definition

$$ \mathbf{B}'_J(x) = -i \mathbf{Z}^{*\dagger}(-x) \mathbf{A}'_J(x), \quad \mathbf{A}'_J(x) = i \mathbf{Y}^{*\dagger}(-x) \mathbf{B}'_J(x), \quad (29) $$

where

$$ \mathbf{Z}^{*\dagger}(-x) = \sum_{n=-N}^{N} \hat{\mathbf{Z}}^{(n)*} e^{ignx}, \quad \mathbf{Y}^{*\dagger}(-x) = \sum_{n=-N}^{N} \hat{\mathbf{Y}}^{(n)*} e^{ignx}, \quad (30) $$

and $\hat{\mathbf{Z}}^{(n)}$ and $\hat{\mathbf{Y}}^{(n)}$ are the Fourier coefficients of $\mathbf{Z}'(x)$ and $\mathbf{Y}'(x)$ relating $\mathbf{B}'_J(-x)$ and $i\mathbf{A}'_J(-x)$. Inserting (28) and (30) in (29) and taking the complex conjugate yields

$$ \hat{\mathbf{B}}_{3M+J}^{(n)} = i \hat{\mathbf{Z}}^{(n-m)} \hat{\mathbf{A}}_{3M+J}^{(m)}, \quad \hat{\mathbf{A}}_{3M+J}^{(n)} = -i \hat{\mathbf{Y}}^{(n-m)} \hat{\mathbf{B}}_{3M+J}^{(m)}. \quad (31) $$

Setting Eqs. (31) into the $3M \times 3M$ form leads to

$$ \tilde{\mathbf{Z}}'(\omega, K_x) = -i\tilde{\Xi}_4\tilde{\Xi}_2^{-1}, \quad \tilde{\mathbf{Y}}'(\omega, K_x) = i\tilde{\Xi}_2\tilde{\Xi}_4^{-1}, \quad (32) $$

where the $3M \times 3M$ block Toöplitz matrices $\tilde{\mathbf{Z}}'$ and $\tilde{\mathbf{Y}}'$ are arranged from $3 \times 3$ blocks $\hat{\mathbf{Z}}^{(n-m)}$ and $\hat{\mathbf{Y}}^{(n-m)}$. As a test, assume the case of a symmetric periodicity profile such that the "direct" and "inverted" half-spaces $y \ge 0$ are identical. Then, by (16)$_1$, $\tilde{\Xi}_1 = \tilde{\Xi}_2^*$, $\tilde{\Xi}_3 = \tilde{\Xi}_4^*$, hence Eqs. (26) and (32) yield $\tilde{\mathbf{Z}} = \tilde{\mathbf{Z}}^{*\dagger}$ and therefore $\mathbf{Z}(x) = \mathbf{Z}^{*\dagger}(x)$, i.e., the results for both above configurations coincide, as they must.

Note that the PWE version of the Barnett–Lothe integral formalism [13] allows direct evaluation of the impedances bypassing the eigenvalue problem of $\tilde{N}$. This may be done via the sign function of the matrix $i\tilde{N}$, which was shown to satisfy the equalities

$$ \mathrm{sign}i\tilde{N} = i\langle\tilde{N}_{\theta}\rangle = i\tilde{\mathbf{Y}}, \quad (33) $$

where $\langle\tilde{N}_{\theta}\rangle = \frac{1}{\pi} \int_{0}^{\pi} \tilde{N}_{\theta} d\theta (= \tilde{T}(\tilde{N}_{\theta})^{+} \tilde{T})$ with the matrix $\tilde{N}_{\theta}$ constructed similarly to (8), (9) but in the rotating frame of vectors $(\mathbf{e}_1, \mathbf{e}_2)$, and

$$ \tilde{\mathbf{Y}} \equiv \begin{pmatrix} \tilde{\mathbf{Y}}_1 & \tilde{\mathbf{Y}}_2 \\ \tilde{\mathbf{Y}}_3 & -\tilde{\mathbf{Y}}_1 \end{pmatrix} = i\tilde{\Xi}\mathrm{diag}(\tilde{\mathbf{I}}, -\tilde{\mathbf{I}})\tilde{\Xi}^{-1} = i \begin{pmatrix} \tilde{\mathbf{I}} - 2\tilde{\Xi}_2\tilde{\Xi}_3^+ & 2\tilde{\Xi}_1\tilde{\Xi}_2^+ \\ 2\tilde{\Xi}_3\tilde{\Xi}_4^+ & \tilde{\mathbf{I}} - 2\tilde{\Xi}_4\tilde{\Xi}_1^+ \end{pmatrix}. \quad (34) $$

Hence

$$ \begin{aligned} \tilde{\mathbf{Z}} &= -\tilde{\mathbf{Y}}_2^{-1}(\mathbf{I} + i\tilde{\mathbf{Y}}_1), \quad \tilde{\mathbf{Y}} = \tilde{\mathbf{Y}}_3^{-1}(\mathbf{I} + i\tilde{\mathbf{Y}}_1^T); \\ \tilde{\mathbf{Z}}' &= -\tilde{\mathbf{Y}}_2^{-1}(\mathbf{I} - i\tilde{\mathbf{Y}}_1), \quad \tilde{\mathbf{Y}}' = \tilde{\mathbf{Y}}_3^{-1}(\mathbf{I} - i\tilde{\mathbf{Y}}_1^T). \end{aligned} \quad (35) $$

where the blocks of $\tilde{\mathbf{Y}}$ may be calculated from the sign function either by integrating the matrix $\tilde{N}_{\theta}$ or else by invoking the corresponding projector functions, see [13].

## 3.2. Properties of the PWE impedance and admittance

According to the on-diagonal blocks of the orthonormality relation (15)$_1$, the impedance matrices (26) and (32) are Hermitian:

$$ \tilde{\mathbf{Z}} = \tilde{\mathbf{Z}}^+, \quad \tilde{\mathbf{Z}}' = \tilde{\mathbf{Z}}^{+\dagger}, \quad (36) $$

and hence so are their inverse, the admittances matrices. In turn, from the off-diagonal blocks of (15)$_1$ it follows that

$$ \tilde{\mathbf{Z}} + \tilde{\mathbf{Z}}' = i(\tilde{\Xi}_1\tilde{\Xi}_2^+)^{-1}(= -2\tilde{\mathbf{Y}}_2^{-1}), \quad \tilde{\mathbf{Y}} + \tilde{\mathbf{Y}}' = -i(\tilde{\Xi}_3\tilde{\Xi}_4^+)^{-1}(= 2\tilde{\mathbf{Y}}_3^{-1}), \quad (37) $$

where $\tilde{\mathbf{Y}}_2 = \tilde{\mathbf{Y}}_3^+$ and $\tilde{\mathbf{Y}}_3 = \tilde{\mathbf{Y}}_4^+$ are blocks of the matrix (34). If the periodicity profile is symmetric, then $\tilde{\mathbf{Z}} = \tilde{\mathbf{Z}}^{*\dagger}$ (see below (32)) and hence (37)$_1$ reduces to $\mathrm{Re}\,\tilde{\mathbf{Z}} = -\tilde{\mathbf{Y}}_2^{-1}$.

The sign-definiteness of impedances and admittances in a laterally periodic half-space rests on their general link to energy quantities and may be proved along similar lines as for other half-space models. First, consider the cycle-averaged kinetic energy density $E(x, y) = \frac{1}{4}\rho(x)\omega^2 |\mathbf{u}(x, y)|^2$ for the evanescent wave field (18). Performing the Fourier-series expansion of $E(x, y)$ in $x$ and taking its zero-order coefficient $\hat{E}_0(y)$ yields

$$ \hat{E}_0(y) = \frac{1}{4}\omega^2 \tilde{\mathbf{A}}^+(y) \tilde{\rho} \tilde{\mathbf{A}}(y), \quad (38) $$
where $\tilde{\rho} = \{\hat{\rho}_{n-m} \mathbf{I}\}$. Using (6)–(9), the value $\hat{E}_0(y)$ may be manipulated to the form

$$ \hat{E}_0(y) = -\frac{1}{8}\omega\tilde{\eta}^+(y) \frac{\partial\tilde{\Pi}\tilde{\mathbf{N}}}{\partial\omega}\tilde{\eta}(y) = \frac{1}{8}i\omega\frac{d}{dy}\tilde{\eta}^+(y) \tilde{\Pi}\frac{\partial\tilde{\eta}(y)}{\partial\omega}. \quad (39) $$

Consider (39) with $\tilde{\eta}(y)$ composed of the evanescent modes. It is described by the sum (11) taken up to $J = 3M$, which may also be written in terms of notations (14) as

$$ \tilde{\eta}(y) = \begin{pmatrix} \tilde{\Xi}_1 \\ \tilde{\Xi}_3 \end{pmatrix} \operatorname{diag}(e^{ik_{y1}y}, \dots, e^{ik_{y,3M}y}) \tilde{\mathbf{c}}, \quad (40) $$

where $\tilde{\mathbf{c}} = (c_1 \cdots c_{3M})^T$. Denote $\langle \hat{E}_0 \rangle \equiv \int_0^\infty \hat{E}_0(y) dy$. Bearing in mind that, by construction, $\hat{E}_0(y) \to 0$ at $y \to \infty$ and invoking (26)$_1$ and (36)$_2$, one finds

$$ \langle \hat{E}_0 \rangle = -\frac{1}{8}i\omega\tilde{\mathbf{c}}^+(\tilde{\Xi}_1^+ \frac{\partial\tilde{\Xi}_3}{\partial\omega} + \tilde{\Xi}_3^+ \frac{\partial\tilde{\Xi}_1}{\partial\omega})\tilde{\mathbf{c}} = -\frac{1}{8}i\omega\tilde{\mathbf{c}}^+\tilde{\Xi}_1^+ \frac{\partial\tilde{Z}}{\partial\omega}\tilde{\Xi}_1\tilde{\mathbf{c}}. \quad (41) $$

Now consider the static limit $\omega \to 0$ of the cycle-averaged internal energy density $W(x, y) = \frac{1}{8}(\sigma_{ij}^* u_{i,j} + \sigma_{ij} u_{i,j}^*)$. Taking into account Eq. (1) and integrating over the period $x \in [0, T]$ yields

$$ \hat{W}_0(y)|_{\omega=0} = \frac{1}{8} \int_0^T \operatorname{div} (\sigma_{ij}^* u_i + \sigma_{ij} u_i^*) dx = \frac{1}{8} \frac{d}{dy} \int_0^T (t_2^* \mathbf{u} + t_2 \mathbf{u}^*) dx. \quad (42) $$

Denote $\langle \hat{W}_0 \rangle \equiv \int_0^\infty \hat{W}_0(y) dy$. Inserting the wave field (18), using its vanishing at $y \to \infty$ and invoking (40) gives

$$ \langle \hat{W}_0 \rangle_{\omega=0} = -\frac{1}{8}i\tilde{\mathbf{c}}^+(\tilde{\Xi}_3^+ \tilde{\Xi}_1 - \tilde{\Xi}_1^+ \tilde{\Xi}_3)\tilde{\mathbf{c}} = \frac{1}{8}i\tilde{\mathbf{c}}^+ \tilde{\Xi}_1^+ \tilde{Z} \tilde{\Xi}_1 \tilde{\mathbf{c}}. \quad (43) $$

Using Eqs. (28), (32)$_1$ and (36)$_2$ leads to the same relations as (41) and (43) but with $\tilde{Z}'$ and $\tilde{\Xi}_2$, $\tilde{\Xi}_4$ replacing $\tilde{Z}$ and $\tilde{\Xi}_1$, $\tilde{\Xi}_3$, respectively.

Since $\langle \hat{W}_0 \rangle$ and $\langle \hat{E}_0 \rangle$ are positive and $\tilde{\mathbf{c}}$, $\tilde{\mathbf{c}}$ are arbitrary 3M-component vectors, the above expressions prove that the impedance matrices $\tilde{Z}(\omega, K_x)$ and $\tilde{Z}'(\omega, K_x)$ taken within the subsonic interval $\omega < \hat{\omega}$ at any fixed $K_x$ are positive definite in the limit $\omega \to 0$ and have negative-definite derivatives in $\omega$. Consequently, the admittances $\tilde{\mathbf{Y}}$ and $\tilde{\mathbf{Y}}'$ as the inverse matrices are positive definite at $\omega \to 0$ and have positive-definite derivatives in $\omega$. The latter, in turn, reveals, via (37), that $\tilde{\mathbf{Y}}_2(\omega, K_x)$ is negative definite and $\tilde{\mathbf{Y}}_3(\omega, K_x)$ is positive definite at $\omega \to 0$, and that they both have negative-definite derivatives in $\omega$.

The above sign-definiteness properties are seen to be similar to those of the counterpart 3 × 3 matrices in a homogeneous half-space [7,9] and transversely periodic half-space [10]. Another common feature is the guaranteed finiteness of the matrices $\tilde{\mathbf{Y}}_2$ and $\tilde{\mathbf{Y}}_3$ at $\omega < \hat{\omega}$ and their divergence at $\omega = \hat{\omega}$ unless, respectively, $\tilde{\mathbf{A}}_{\deg}$ or $\tilde{\mathbf{B}}_{\deg}$ in the eigenvector $\tilde{\xi}_{\deg}$ accidentally vanish, i.e. $\tilde{\Xi}_1$ or $\tilde{\Xi}_3$ is singular at $\hat{\omega}$. Both statements may be proved by way of expanding eigenvectors of the “almost non-semisimple” matrix in the same way as it was done for the transversely periodic case in [10]. Once they are established, the negative definiteness of (finite) $\tilde{\mathbf{Y}}_2$ at $\omega = 0$ and of its frequency derivative at $\omega < \hat{\omega}$ implies that $\tilde{\mathbf{Y}}_2$ and hence $\tilde{\Xi}_1$, $\tilde{\Xi}_2$ and $\tilde{\mathbf{Y}}$, $\tilde{\mathbf{Y}}'$ are non-singular, while the finiteness of $\tilde{\mathbf{Y}}_3$ means that so are $\tilde{Z}$ and $\tilde{Z}'$.

### 3.3. Properties of the x-dependent impedance and admittance

Let us pass to the ensuing properties of the 3 × 3 impedance matrices $\mathbf{Z}(x)$ and $\mathbf{Z}'(x)$ as T-periodic functions in the x-space. It is known that the matrix function whose Fourier coefficients form a Hermitian and sign-definite (or singular) block Toéplitz matrix is itself Hermitian and similarly sign-definite (or singular) for any $x \in [0, T]$, see e.g. [15]. Hence, from (36) and the above-listed properties of $\tilde{Z}$ and $\tilde{Z}'$ it follows that

$$ \mathbf{Z}(x) = \mathbf{Z}^+(x), \quad \mathbf{Z}'(x) = \mathbf{Z}^{\prime+}(x), \quad (44) $$

and that, for any fixed $K_x$ and $\omega < \hat{\omega}$,

$\mathbf{Z}(x)$ and $\mathbf{Z}'(x)$ are positive definite at $\omega = 0$,

$$ \frac{\partial \mathbf{Z}(x)}{\partial \omega} \text{ and } \frac{\partial \mathbf{Z}'(x)}{\partial \omega} \text{ are negative definite;} \quad (45) $$

$\mathbf{Y}(x)$ and $\mathbf{Y}'(x)$ are positive definite at $\omega = 0$,

$$ \frac{\partial \mathbf{Y}(x)}{\partial \omega} \text{ and } \frac{\partial \mathbf{Y}'(x)}{\partial \omega} \text{ are positive definite.} \quad (46) $$

Next, introduce the 3 × 3 matrices

$$ \mathbf{Y}_2(x) = 2i \{ \mathbf{A}_J(x) \} \{ \mathbf{A}'_J(-x) \}^T, \quad \mathbf{Y}_3(x) = 2i \{ \mathbf{B}_J(x) \} \{ \mathbf{B}'_J(-x) \}^T, \quad (46) $$

where the symbol $\{\cdot\}$ implies a 3 × 3M matrix composed of enclosed column vectors: e.g., $\{\mathbf{A}_J(x)\} = \|\mathbf{A}_1(x) \cdots \mathbf{A}_{3M}(x)\|$. Inserting $\{\mathbf{B}'_J(-x)\}^T = -i\mathbf{Z}^{\prime+}(x) \{\mathbf{A}'_J(-x)\}^T$ (see (29)$_1$ and (44)$_2$) yields

$$ \mathbf{Y}_3(x) = -2i\mathbf{Z}(x)\mathbf{Y}_2(x)\mathbf{Z}'(x). \quad (47) $$
where $\tilde{\rho} = \{\hat{\rho}_{n-m} I\}$. Using (6)–(9), the value $\hat{E}_0(y)$ may be manipulated to the form

$$ \hat{E}_0(y) = -\frac{1}{8}\omega\tilde{\eta}^+(y) \frac{\partial\tilde{\Pi}\tilde{N}}{\partial\omega}\tilde{\eta}(y) = \frac{1}{8}i\omega\frac{d}{dy}\tilde{\eta}^+(y) \tilde{\Pi}\frac{\partial\tilde{\eta}(y)}{\partial\omega}. \quad (39) $$

Consider (39) with $\tilde{\eta}(y)$ composed of the evanescent modes. It is described by the sum (11) taken up to $J = 3M$, which may also be written in terms of notations (14) as

$$ \tilde{\eta}(y) = \begin{pmatrix} \tilde{\Xi}_1 \\ \tilde{\Xi}_3 \end{pmatrix} \operatorname{diag}(e^{ik_{y1}y}, \dots, e^{ik_{y,3M}y}) \tilde{c}, \quad (40) $$

where $\tilde{c} = (c_1 \cdots c_{3M})^T$. Denote $\langle \hat{E}_0 \rangle \equiv \int_0^\infty \hat{E}_0(y) dy$. Bearing in mind that, by construction, $\hat{E}_0(y) \to 0$ at $y \to \infty$ and invoking (26)$_1$ and (36)$_2$, one finds

$$ \langle \hat{E}_0 \rangle = -\frac{1}{8}i\omega\tilde{c}^+(\tilde{\Xi}_1^+ \frac{\partial\tilde{\Xi}_3}{\partial\omega} + \tilde{\Xi}_3^+ \frac{\partial\tilde{\Xi}_1}{\partial\omega})\tilde{c} = -\frac{1}{8}i\omega\tilde{c}^+\tilde{\Xi}_1^+ \frac{\partial\tilde{Z}}{\partial\omega}\tilde{\Xi}_1\tilde{c}. \quad (41) $$

Now consider the static limit $\omega \to 0$ of the cycle-averaged internal energy density $W(x, y) = \frac{1}{8}(\sigma_{ij}^* u_{i,j} + \sigma_{ij} u_{i,j}^*)$. Taking into account Eq. (1) and integrating over the period $x \in [0, T]$ yields

$$ \hat{W}_0(y) \Big|_{\omega=0} = \frac{1}{8} \int_0^T \operatorname{div} (\sigma_{ij}^* u_i + \sigma_{ij} u_i^*) dx = \frac{1}{8} \frac{d}{dy} \int_0^T (t_2^* u + t_2 u^*) dx. \quad (42) $$

Denote $\langle \hat{W}_0 \rangle \equiv \int_0^\infty \hat{W}_0(y) dy$. Inserting the wave field (18), using its vanishing at $y \to \infty$ and invoking (40) gives

$$ \langle \hat{W}_0 \rangle_{\omega=0} = -\frac{1}{8}i\tilde{c}^+(\tilde{\Xi}_3^+ \tilde{\Xi}_1 - \tilde{\Xi}_1^+ \tilde{\Xi}_3)\tilde{c} = \frac{1}{8}i\tilde{c}^+ \tilde{\Xi}_1^+ \tilde{Z} \tilde{\Xi}_1 \tilde{c}. \quad (43) $$

Using Eqs. (28), (32)$_1$ and (36)$_2$ leads to the same relations as (41) and (43) but with $\tilde{Z}'$ and $\tilde{\Xi}_2$, $\tilde{\Xi}_4$ replacing $\tilde{Z}$ and $\tilde{\Xi}_1$, $\tilde{\Xi}_3$, respectively.

Since $\langle \hat{W}_0 \rangle$ and $\langle \hat{E}_0 \rangle$ are positive and $\tilde{c}$, $\tilde{c}$ are arbitrary 3M-component vectors, the above expressions prove that the impedance matrices $\tilde{Z}(\omega, K_x)$ and $\tilde{Z}'(\omega, K_x)$ taken within the subsonic interval $\omega < \hat{\omega}$ at any fixed $K_x$ are positive definite in the limit $\omega \to 0$ and have negative-definite derivatives in $\omega$. Consequently, the admittances $\tilde{Y}$ and $\tilde{Y}'$ as the inverse matrices are positive definite at $\omega \to 0$ and have positive-definite derivatives in $\omega$. The latter, in turn, reveals, via (37), that $\tilde{\mathbf{Y}}_2(\omega, K_x)$ is negative definite and $\tilde{\mathbf{Y}}_3(\omega, K_x)$ is positive definite at $\omega \to 0$, and that they both have negative-definite derivatives in $\omega$.

The above sign-definiteness properties are seen to be similar to those of the counterpart 3 × 3 matrices in a homogeneous half-space [7,9] and transversely periodic half-space [10]. Another common feature is the guaranteed finiteness of the matrices $\tilde{\mathbf{Y}}_2$ and $\tilde{\mathbf{Y}}_3$ at $\omega < \hat{\omega}$ and their divergence at $\omega = \hat{\omega}$ unless, respectively, $\tilde{\mathbf{A}}_{\deg}$ or $\tilde{\mathbf{B}}_{\deg}$ in the eigenvector $\tilde{\xi}_{\deg}$ accidentally vanish, i.e. $\tilde{\Xi}_1$ or $\tilde{\Xi}_3$ is singular at $\hat{\omega}$. Both statements may be proved by way of expanding eigenvectors of the “almost non-semisimple” matrix in the same way as it was done for the transversely periodic case in [10]. Once they are established, the negative definiteness of (finite) $\tilde{\mathbf{Y}}_2$ at $\omega = 0$ and of its frequency derivative at $\omega < \hat{\omega}$ implies that $\tilde{\mathbf{Y}}_2$ and hence $\tilde{\Xi}_1$, $\tilde{\Xi}_2$ and $\tilde{\mathbf{Y}}$, $\tilde{\mathbf{Y}}'$ are non-singular, while the finiteness of $\tilde{\mathbf{Y}}_3$ means that so are $\tilde{\mathbf{Z}}$ and $\tilde{\mathbf{Z}}'$.

### 3.3. Properties of the x-dependent impedance and admittance

Let us pass to the ensuing properties of the 3 × 3 impedance matrices $\mathbf{Z}(x)$ and $\mathbf{Z}'(x)$ as T-periodic functions in the x-space. It is known that the matrix function whose Fourier coefficients form a Hermitian and sign-definite (or singular) block Toéplitz matrix is itself Hermitian and similarly sign-definite (or singular) for any $x \in [0, T]$, see e.g. [15]. Hence, from (36) and the above-listed properties of $\tilde{\mathbf{Z}}$ and $\tilde{\mathbf{Z}}'$ it follows that

$$ \mathbf{Z}(x) = \mathbf{Z}^+(x), \quad \mathbf{Z}'(x) = \mathbf{Z}^{\prime+}(x), \quad (44) $$

and that, for any fixed $K_x$ and $\omega < \hat{\omega}$,

$\mathbf{Z}(x)$ and $\mathbf{Z}'(x)$ are positive definite at $\omega = 0$,

$$ \frac{\partial \mathbf{Z}(x)}{\partial \omega} \text{ and } \frac{\partial \mathbf{Z}'(x)}{\partial \omega} \text{ are negative definite;} \quad (45) $$

$\mathbf{Y}(x)$ and $\mathbf{Y}'(x)$ are positive definite at $\omega = 0$,

$$ \frac{\partial \mathbf{Y}(x)}{\partial \omega} \text{ and } \frac{\partial \mathbf{Y}'(x)}{\partial \omega} \text{ are positive definite.} \quad (46) $$

Next, introduce the 3 × 3 matrices

$$ \mathbf{Y}_2(x) = 2i \{ \mathbf{A}_J(x) \} \{ \mathbf{A}'_J(-x) \}^T, \quad \mathbf{Y}_3(x) = 2i \{ \mathbf{B}_J(x) \} \{ \mathbf{B}'_J(-x) \}^T, \quad (46) $$

where the symbol $\{\cdot\}$ implies a 3 × 3M matrix composed of enclosed column vectors: e.g., $\{\mathbf{A}_J(x)\} = \|\mathbf{A}_1(x) \cdots \mathbf{A}_{3M}(x)\|$. Inserting $\{\mathbf{B}'_J(-x)\}^T = -i\mathbf{Z}^{\prime+}(x) \{\mathbf{A}'_J(-x)\}^T$ (see (29)$_1$ and (44)$_2$) yields

$$ \mathbf{Y}_3(x) = -2i\mathbf{Z}(x)\mathbf{Y}_2(x)\mathbf{Z}'(x). \quad (47) $$
Bearing in mind notations (14), one can verify that matrices (46) have the same sign-definiteness properties as the 3M × 3M matrices $\tilde{\Upsilon}_2 = 2i\tilde{\Xi}_1\tilde{\Xi}_2^+$ and $\tilde{\Upsilon}_3 = 2i\tilde{\Xi}_3\tilde{\Xi}_4^+$ which were formulated at the end of §3.2. Thus,

$$ \begin{array}{l} \Upsilon_2(x) \text{ is negative definite at } \omega = 0, \\ \frac{\partial \Upsilon_2(x)}{\partial \omega} \text{ is negative definite;} \end{array} \qquad (48) $$

$$ \begin{array}{l} \Upsilon_3(x) \text{ is positive definite at } \omega = 0, \\ \frac{\partial \Upsilon_3(x)}{\partial \omega} \text{ is negative definite.} \end{array} $$

Note aside that non-singularity of the matrix $\Upsilon_2(x)$, which follows from (48)$_{1,2}$, rules out subsonic surface waves in a superlattice with clamped boundary condition $\mathbf{u}(x, 0) = \mathbf{0}$.

Other important properties taken over in the x-space are the finiteness of the matrices $\mathbf{Z}(x)$, $\mathbf{Z}'(x)$ and $\Upsilon_2(x)$, $\Upsilon_3(x)$ at $\omega < \hat{\omega}$, and divergence of $\Upsilon_2(x)$, $\Upsilon_3(x)$ at $\omega = \hat{\omega}$ unless, respectively, $\{\mathbf{A}_J(x)\}$ or $\{\mathbf{B}_J(x)\}$ is accidentally singular at the transonic state $\hat{\omega}$.

## 4. Surface waves on a free boundary

Let the surface of the half-space $y \ge 0$ with a laterally periodic profile of material properties $\rho(x)$ and $c_{ijkl}(x)$ or with its “inverted” version $\rho'(x) = \rho(-x)$ and $c'_{ijkl}(x) = c_{ijkl}(-x)$ maintain zero traction $t_2(x, 0) = \mathbf{0}$ at any $x$. By (18) and (28), this means the vanishing of, respectively, the following linear combinations of traction vectors:

$$ \begin{array}{l} \tilde{\Xi}_3 \tilde{c} = \mathbf{0} \Leftrightarrow \| \mathbf{B}_1(x) \dots \mathbf{B}_{3M}(x) \| \tilde{c} = \mathbf{0} \quad \forall x, \\ \tilde{\Xi}_4 \tilde{c}' = \mathbf{0} \Leftrightarrow \| \mathbf{B}'_1(x) \dots \mathbf{B}'_{3M}(x) \| \tilde{c}' = \mathbf{0} \quad \forall x, \end{array} \qquad (49) $$

where $\tilde{c} = (c_1 \cdots c_{3M})^T$ and $\tilde{c}' = (c'_1 \cdots c'_{3M})^T$ are constants. Hence, according to (26) and (32), the surface-wave dispersion branches $\omega(K_x)$ in the two above cases may be defined, respectively, from the equations

$$ \begin{array}{l} \det \tilde{\mathbf{Z}}(\omega, K_x) = 0 \Leftrightarrow \det \mathbf{Z}(x; \omega, K_x) = 0 \quad \forall x, \\ \det \tilde{\mathbf{Z}}'(\omega, K_x) = 0 \Leftrightarrow \det \mathbf{Z}'(x; \omega, K_x) = 0 \quad \forall x. \end{array} \qquad (50) $$

We consider the subsonic frequency interval $\omega \le \hat{\omega}(K_x)$ where the impedance matrices are Hermitian, hence Eq. (50) involves real functions (note that vanishing of determinants of $\tilde{\Xi}_3$ or $\tilde{\Xi}_4$ is also equivalent to a real equation, since, say, $\det \tilde{\Xi}_3 = 0$ implies $\det \tilde{\Xi}_3^* = 0$ and vice versa, as may be seen, e.g., from (15)$_2$). Specifically, the left-hand algebraic equations in (50) involve explicitly accessible 3M × 3M PWE matrices and thereby are appropriate for direct finding of surface-wave branches $\omega(K_x)$, see [13] for more details. On the other hand, the right-hand equations in (50) tell us about the possible number of these branches.

First, we note that, by (45)$_{1,2}$, the three eigenvalues of $\mathbf{Z}(x)$ and $\mathbf{Z}'(x)$ taken at any fixed $x$ and $K_x$ are continuously decreasing functions of $\omega \le \hat{\omega}$ positive at $\omega = 0$, so they may have at most three zeros at $\omega = \omega(K_x)$. In consequence, each of the right-hand equations in (50) considered at fixed $x$ and $K_x$ may have at most three solutions $\omega = \omega(K_x)$ and hence at most three ones $\omega = \omega(K_x)$ such that are the same for any $x \in [0, T]$. Thus, a half-space with a generic shape of lateral periodicity profile admits at most three surface waves at fixed $K_x$.

Now consider an arbitrary superlattice and its "inverted" version obtained from the given one by rotating 180° around the axis $Y$, which is equivalent to the inversion of the stratification axis $X$. Let us ask a question as to how many surface waves may propagate in these two superlattices in total. A somewhat counter-intuitive answer follows from Eqs. (47), (48)$_{3,4}$ and

$$ \det \Upsilon_3(x; \omega, K_x) = 0. \qquad (51) $$

Taking into account (48)$_{3,4}$ and finiteness of $\Upsilon_2(x)$ and $\Upsilon_3(x)$, we observe that Eq. (51) admits at most three solutions $\omega = \omega(K_x)$. Thus, for any fixed $K_x$, there may be at most three surface waves in total propagating in two mutually "inverse" superlattices. By (47), they are split, one way or another, between the solutions of equations (50)$_1$ and (50)$_2$. In other words, if, at a given $K_x$, one of these superlattices admits three or two or one wave, the other cannot support surface waves or admit at most one or two, respectively.

Furthermore, due to (48)$_{3,4}$ and the divergence of $\Upsilon_3(x)$ at the normal (not exceptional) transonic state $\omega = \hat{\omega}$ such that keeps $\{\mathbf{B}_J(x)\}$ non-singular, at least one of the eigenvalues of $\Upsilon_3(x)$ varies monotonically from positive value at $\omega = 0$ to $-\infty$ at $\omega = \hat{\omega}$ and hence turns to zero in between. Thus, with reference to (47), at least one surface wave is guaranteed to exist at any fixed $K_x$ for one of the mutually "inverse" superlattices, provided that the transonic state is normal.

The special case is when a superlattice has a symmetric periodicity profile such that it is invariant to the inversion $X \to -X$ about an appropriately chosen origin $x = 0$. Then $\mathbf{Z}(x) = \mathbf{Z}'^*(x)$, which implies that $\mathrm{Re}\,\mathbf{Z}(x)$ is positive definite at $\hat{\omega}$, and, independently, that a zero eigenvalue of $\Upsilon_3(x; \omega, K_x)$, if exists at $\omega < \hat{\omega}$, is a double one, see (47). Either of these arguments leads to the conclusion that Eq. (51) admits only one solution at $\omega < \hat{\omega}$. Thus, a superlattice with a symmetric profile admits only one surface wave at fixed $K_x$. In particular, this is the case for any periodically bilayered superlattice since an infinite sequence of pairs A-B can equally be seen as a sequence of symmetric patterns $\frac{1}{2}A-B-\frac{1}{2}A$.

While it follows from the above that at least one surface wave branch usually does exist in a laterally periodic superlattice, the likelihood of occurrence of two or three surface wave branches is yet to be explored. Certain insight may be gained via explicit
Bearing in mind notations (14), one can verify that matrices (46) have the same sign-definiteness properties as the 3M × 3M matrices $\tilde{\mathbf{\Upsilon}}_2 = 2i\tilde{\mathbf{\Xi}}_1\tilde{\mathbf{\Xi}}_2^+$ and $\tilde{\mathbf{\Upsilon}}_3 = 2i\tilde{\mathbf{\Xi}}_3\tilde{\mathbf{\Xi}}_4^+$ which were formulated at the end of §3.2. Thus,

$$ \mathbf{\Upsilon}_2(x) \text{ is negative definite at } \omega = 0, $$
$$ \frac{\partial \mathbf{\Upsilon}_2(x)}{\partial \omega} \text{ is negative definite;} \qquad (48) $$

$$ \mathbf{\Upsilon}_3(x) \text{ is positive definite at } \omega = 0, $$
$$ \frac{\partial \mathbf{\Upsilon}_3(x)}{\partial \omega} \text{ is negative definite.} $$

Note aside that non-singularity of the matrix $\mathbf{\Upsilon}_2(x)$, which follows from (48)$_{1,2}$, rules out subsonic surface waves in a superlattice with clamped boundary condition $\mathbf{u}(x, 0) = \mathbf{0}$.

Other important properties taken over in the x-space are the finiteness of the matrices $\mathbf{Z}(x)$, $\mathbf{Z}'(x)$ and $\mathbf{\Upsilon}_2(x)$, $\mathbf{\Upsilon}_3(x)$ at $\omega < \hat{\omega}$, and divergence of $\mathbf{\Upsilon}_2(x)$, $\mathbf{\Upsilon}_3(x)$ at $\omega = \hat{\omega}$ unless, respectively, $\{\mathbf{A}_J(x)\}$ or $\{\mathbf{B}_J(x)\}$ is accidentally singular at the transonic state $\hat{\omega}$.

## 4. Surface waves on a free boundary

Let the surface of the half-space $y \ge 0$ with a laterally periodic profile of material properties $\rho(x)$ and $c_{ijkl}(x)$ or with its “inverted” version $\rho'(x) = \rho(-x)$ and $c'_{ijkl}(x) = c_{ijkl}(-x)$ maintain zero traction $t_2(x, 0) = \mathbf{0}$ at any $x$. By (18) and (28), this means the vanishing of, respectively, the following linear combinations of traction vectors:

$$ \tilde{\mathbf{\Xi}}_3 \tilde{\mathbf{c}} = \mathbf{0} \Leftrightarrow \| \mathbf{B}_1(x) \dots \mathbf{B}_{3M}(x) \| \tilde{\mathbf{c}} = \mathbf{0} \quad \forall x, $$
$$ \tilde{\mathbf{\Xi}}_4 \tilde{\mathbf{c}}' = \mathbf{0} \Leftrightarrow \| \mathbf{B}'_1(x) \dots \mathbf{B}'_{3M}(x) \| \tilde{\mathbf{c}}' = \mathbf{0} \quad \forall x, \qquad (49) $$

where $\tilde{\mathbf{c}} = (c_1 \cdots c_{3M})^T$ and $\tilde{\mathbf{c}}' = (c'_1 \cdots c'_{3M})^T$ are constants. Hence, according to (26) and (32), the surface-wave dispersion branches $\omega(K_x)$ in the two above cases may be defined, respectively, from the equations

$$ \det \tilde{\mathbf{Z}}(\omega, K_x) = 0 \Leftrightarrow \det \mathbf{Z}(x; \omega, K_x) = 0 \quad \forall x, $$
$$ \det \tilde{\mathbf{Z}}'(\omega, K_x) = 0 \Leftrightarrow \det \mathbf{Z}'(x; \omega, K_x) = 0 \quad \forall x. \qquad (50) $$

We consider the subsonic frequency interval $\omega \le \hat{\omega}(K_x)$ where the impedance matrices are Hermitian, hence Eq. (50) involves real functions (note that vanishing of determinants of $\tilde{\mathbf{\Xi}}_3$ or $\tilde{\mathbf{\Xi}}_4$ is also equivalent to a real equation, since, say, $\det \tilde{\mathbf{\Xi}}_3 = 0$ implies $\det \tilde{\mathbf{\Xi}}_3^* = 0$ and vice versa, as may be seen, e.g., from (15)$_2$). Specifically, the left-hand algebraic equations in (50) involve explicitly accessible 3M × 3M PWE matrices and thereby are appropriate for direct finding of surface-wave branches $\omega(K_x)$, see [13] for more details. On the other hand, the right-hand equations in (50) tell us about the possible number of these branches.

First, we note that, by (45)$_{1,2}$, the three eigenvalues of $\mathbf{Z}(x)$ and $\mathbf{Z}'(x)$ taken at any fixed $x$ and $K_x$ are continuously decreasing functions of $\omega \le \hat{\omega}$ positive at $\omega = 0$, so they may have at most three zeros at $\omega = \omega(K_x)$. In consequence, each of the right-hand equations in (50) considered at fixed $x$ and $K_x$ may have at most three solutions $\omega = \omega(K_x)$ and hence at most three ones $\omega = \omega(K_x)$ such that are the same for any $x \in [0, T]$. Thus, a half-space with a generic shape of lateral periodicity profile admits at most three surface waves at fixed $K_x$.

Now consider an arbitrary superlattice and its "inverted" version obtained from the given one by rotating 180° around the axis $Y$, which is equivalent to the inversion of the stratification axis $X$. Let us ask a question as to how many surface waves may propagate in these two superlattices in total. A somewhat counter-intuitive answer follows from Eqs. (47), (48)$_{3,4}$ and

$$ \det \mathbf{\Upsilon}_3(x; \omega, K_x) = 0. \qquad (51) $$

Taking into account (48)$_{3,4}$ and finiteness of $\mathbf{\Upsilon}_2(x)$ and $\mathbf{\Upsilon}_3(x)$, we observe that Eq. (51) admits at most three solutions $\omega = \omega(K_x)$. Thus, for any fixed $K_x$, there may be at most three surface waves in total propagating in two mutually "inverse" superlattices. By (47), they are split, one way or another, between the solutions of equations (50)$_1$ and (50)$_2$. In other words, if, at a given $K_x$, one of these superlattices admits three or two or one wave, the other cannot support surface waves or admit at most one or two, respectively.

Furthermore, due to (48)$_{3,4}$ and the divergence of $\mathbf{\Upsilon}_3(x)$ at the normal (not exceptional) transonic state $\omega = \hat{\omega}$ such that keeps $\{\mathbf{B}_J(x)\}$ non-singular, at least one of the eigenvalues of $\mathbf{\Upsilon}_3(x)$ varies monotonically from positive value at $\omega = 0$ to $-\infty$ at $\omega = \hat{\omega}$ and hence turns to zero in between. Thus, with reference to (47), at least one surface wave is guaranteed to exist at any fixed $K_x$ for one of the mutually "inverse" superlattices, provided that the transonic state is normal.

The special case is when a superlattice has a symmetric periodicity profile such that it is invariant to the inversion $X \to -X$ about an appropriately chosen origin $x = 0$. Then $\mathbf{Z}(x) = \mathbf{Z}'^*(x)$, which implies that $\mathrm{Re}\,\mathbf{Z}(x)$ is positive definite at $\hat{\omega}$, and, independently, that a zero eigenvalue of $\mathbf{\Upsilon}_3(x; \omega, K_x)$, if exists at $\omega < \hat{\omega}$, is a double one, see (47). Either of these arguments leads to the conclusion that Eq. (51) admits only one solution at $\omega < \hat{\omega}$. Thus, a superlattice with a symmetric profile admits only one surface wave at fixed $K_x$. In particular, this is the case for any periodically bilayered superlattice since an infinite sequence of pairs A-B can equally be seen as a sequence of symmetric patterns $\frac{1}{2}A-B-\frac{1}{2}A$.

While it follows from the above that at least one surface wave branch usually does exist in a laterally periodic superlattice, the likelihood of occurrence of two or three surface wave branches is yet to be explored. Certain insight may be gained via explicit
derivation of the effect of weak periodic perturbation of properties of a homogeneous half-space. It may be expected that the "reference" (unperturbed) Rayleigh wave will give rise to two slightly different dispersion branches of surface "quasi-Rayleigh" waves associated, respectively, with the "direct" and "inverted" periodicity profiles. The third surface-wave branch may be conjectured to evolve for one of these profiles (but not for the other) from the transonic frequency in the unperturbed half-space, provided this transonic state is of the exceptional type. The perturbation theory intended to verify the above scenario is developed in Part II.

## Part II

### Weak periodic modulation

#### 5. Preamble

Let the periodic modulation of the density and stiffness coefficients be weak in the sense that the total variation $V_0^T(f)$ of the functions $f(x) \equiv \rho(x)$, $c_{ijkl}(x)$ and hence the norm of their nth Fourier coefficients $||f_n||$ ($n \in \mathbb{Z}/\{0\}$), both relative to the average $\langle|f|\rangle = \frac{1}{T} \int_0^T |f(x)| dx$, are small quantities:

$$ \frac{V_0^T(f)}{\langle|f|\rangle} \sim \epsilon, \quad \frac{\|f_n\|}{\langle|f|\rangle} \sim \frac{\epsilon}{|n|}, \quad \epsilon \ll 1, \qquad (52) $$

see [16]. Under this assumption, the wave solutions may be sought within the perturbation approach, where the unperturbed (zero-order) reference is provided by the solutions in a fictitious homogeneous medium characterised by the statically averaged values $\hat{\rho}_0 = \frac{1}{T} \int_0^T \rho(x) dx$ and $\hat{\varepsilon}_{ijkl}^{(0)} = \frac{1}{T} \int_0^T \hat{c}_{ijkl}(x) dx$.

As mentioned in Section 3.1, the numerical calculation of the surface wave spectrum $\omega(K_x)$ based on the sign function of the PWE system matrix $\tilde{N}(\omega, K_x)$ circumvents the partial eigenmodes [13]. On the other hand, they are involved in the subsequent perturbation theory derivation. In this context, it is expedient to proceed from the PWE eigenrelation (12) in the reshuffled form:

$$ (\bar{\mathbf{N}} - k_{yJ} \mathbf{I}) \bar{\xi}_J = \mathbf{0}, \quad J = 1, \dots, 6M, \qquad (53) $$

where the $6M \times 6M$ matrix $\bar{\mathbf{N}} = \{\hat{\mathbf{N}}^{(mn)}\}$, $m, n = -N, \dots, N$, consists of the $6 \times 6$ blocks

$$ \bar{\mathbf{N}} = \{\hat{\mathbf{N}}^{(mn)}\}, \quad \hat{\mathbf{N}}^{(mn)} = \begin{pmatrix} -\hat{(\mathbf{T}^{-1})^{(ml)}(\hat{\mathbf{R}}^{+})^{(ln)}} & -\hat{(\mathbf{T}^{-1})^{(mn)}} \\ \hat{\mathbf{P}}^{(mn)} - \hat{\mathbf{R}}^{(ml)}\hat{(\mathbf{T}^{-1})^{(lk)}}(\hat{\mathbf{R}}^{+})^{(kn)} & -\hat{\mathbf{R}}^{(ml)}\hat{(\mathbf{T}^{-1})^{(ln)}} \end{pmatrix}, \qquad (54) $$

expressed in matrices (9), and the 6-component vectors $\hat{\xi}_J^{(n)} = (\hat{\mathbf{A}}_J^{(n)} \hat{\mathbf{B}}_J^{(n)})^T$ consist of the 3-vectors $\hat{\mathbf{A}}_J^{(n)}$, $\hat{\mathbf{B}}_J^{(n)}$ appearing in (17), (18). Identity (10) yields $\bar{\mathbf{N}} = \bar{\mathbf{T}}\bar{\mathbf{N}}^+\bar{\mathbf{T}}$ where $\bar{\mathbf{T}} = \text{diag}\{\hat{\mathbf{T}}\}$ and $\hat{\mathbf{T}}$ has zero diagonal and identity off-diagonal $3 \times 3$ blocks.

#### 6. Zero-order PWE formulation

##### 6.1. Partial modes

Setting to zero the variation of material properties relative to their mean values specifies Eq. (53) to the form

$$ (\bar{\mathbf{N}}_0 - k_{y(n,a)}^{(0)}) \mathbf{I} \bar{\xi}_{0(n,a)} = \mathbf{0}, \qquad (55) $$

where $\bar{\mathbf{N}}_0(\omega, K_x) = \text{diag}\{\hat{\mathbf{N}}^{(n,n)}\}$ is a $6M \times 6M$ block-diagonal matrix, the double index $(n, \alpha)$, in which $n = -N, \dots, N$ and $\alpha = 1, \dots, 6$, replaces $J = 1, \dots, 6M$ used so far, and the additional index "0" is a reference to the zero-order state. Each eigenvector

$$ \bar{\xi}_{0(n,a)} = (0 \cdots \hat{\xi}_{0\alpha}^{(n)} \cdots 0)^T \text{ with } \hat{\xi}_{0\alpha}^{(n)} = (\hat{\mathbf{A}}_{0\alpha}^{(n)} \hat{\mathbf{B}}_{0\alpha}^{(n)})^T \qquad (56) $$

contains a single non-zero 6-tuple $\hat{\xi}_{0\alpha}^{(n)}$ at the nth position, which is the eigenvector of the nth on-diagonal block $\hat{\mathbf{N}}^{(nn)}$ so that

$$ (\hat{\mathbf{N}}^{(nn)} - k_{y(n,a)}^{(0)}) \mathbf{I} \hat{\xi}_{0\alpha}^{(n)} = \mathbf{0}, \quad \det(\hat{\mathbf{N}}^{(nn)} - k_{y(n,a)}^{(0)}) \mathbf{I} = \mathbf{0}, \quad \alpha = 1, \dots, 6. \qquad (57) $$

The blocks $\hat{\mathbf{N}}^{(nn)}$ represent each the Stroh matrix of the zero-order homogeneous medium while they differ by the value of the horizontal wavenumber $k_n$. The characteristic polynomial of $\hat{\mathbf{N}}^{(nn)}$ is homogeneous of degree 2 in $\omega$ and $k_n$, hence $k_{y(n,a)}^{(0)}(\omega, K_x) = p_a(\nu)k_n$ where $\nu = \omega/k_n$. The Stroh matrix is real, hence so is $\bar{\mathbf{N}}_0$; therefore, identity (54) reduces to the form $\bar{\mathbf{N}}_0 = \bar{\mathbf{T}}\bar{\mathbf{N}}_0^T\bar{\mathbf{T}}$. Given the common case of semisimple (diagonalisable) $\bar{\mathbf{N}}_0$, the matrix $\bar{\mathbf{\Xi}}_0$ of its eigenvectors satisfies the orthonormality relation

$$ \bar{\mathbf{\Xi}}_0^T \bar{\mathbf{T}} \bar{\mathbf{\Xi}}_0 = \mathbf{I} \Leftrightarrow \bar{\xi}_{0(n,a)}^T \bar{\mathbf{T}} \bar{\xi}_{0(m,\beta)} = \hat{\xi}_{0\alpha}^{(n)T} \hat{\mathbf{T}} \hat{\xi}_{0\beta}^{(m)} = \delta_{nm} \delta_{a\beta}, \qquad (58) $$

where $n, m = -N, \dots, N$ and $\alpha, \beta = 1, \dots, 6$.

For a homogeneous medium, the subsonic interval (the one where no bulk modes with real $k_y$ are admitted) is defined in terms of the trace velocity $v = \omega/k_x$ and is bounded by the first transonic state $v = \hat{v}$, at which the normal to the surface is tangent to the outmost bulk-wave slowness sheet, see [8]. Let $\hat{v}$ be such a value for the reference homogeneous material with the averaged density and stiffness $\hat{\rho}_0$ and $\hat{c}_{ijkl}^{(0)}$. Passing to the PWE formalism requires operating with frequency rather than velocity. Denote

$$ \hat{\omega}_n^{(0)} = \hat{v} |k_n|, \quad n = -N, \dots, N, \qquad (59) $$
derivation of the effect of weak periodic perturbation of properties of a homogeneous half-space. It may be expected that the "reference" (unperturbed) Rayleigh wave will give rise to two slightly different dispersion branches of surface "quasi-Rayleigh" waves associated, respectively, with the "direct" and "inverted" periodicity profiles. The third surface-wave branch may be conjectured to evolve for one of these profiles (but not for the other) from the transonic frequency in the unperturbed half-space, provided this transonic state is of the exceptional type. The perturbation theory intended to verify the above scenario is developed in Part II.

# Part II

## Weak periodic modulation

### 5. Preamble

Let the periodic modulation of the density and stiffness coefficients be weak in the sense that the total variation $V_0^T(f)$ of the functions $f(x) \equiv \rho(x)$, $c_{ijkl}(x)$ and hence the norm of their nth Fourier coefficients $||f_n|| (n \in \mathbb{Z}/\{0\})$, both relative to the average $\langle|f|\rangle = \frac{1}{T} \int_0^T |f(x)| dx$, are small quantities:

$$ \frac{V_0^T(f)}{\langle|f|\rangle} \sim \varepsilon, \quad \frac{\|f_n\|}{\langle|f|\rangle} \sim \frac{\varepsilon}{|n|}, \quad \varepsilon \ll 1, \qquad (52) $$

see [16]. Under this assumption, the wave solutions may be sought within the perturbation approach, where the unperturbed (zero-order) reference is provided by the solutions in a fictitious homogeneous medium characterised by the statically averaged values $\hat{\rho}_0 = \frac{1}{T} \int_0^T \rho(x) dx$ and $\hat{\varepsilon}_{ijkl}^{(0)} = \frac{1}{T} \int_0^T \hat{c}_{ijkl}(x) dx$.

As mentioned in Section 3.1, the numerical calculation of the surface wave spectrum $\omega(K_x)$ based on the sign function of the PWE system matrix $\tilde{N}(\omega, K_x)$ circumvents the partial eigenmodes [13]. On the other hand, they are involved in the subsequent perturbation theory derivation. In this context, it is expedient to proceed from the PWE eigenrelation (12) in the reshuffled form:

$$ (\bar{\mathbf{N}} - k_{yJ} \mathbf{I}) \bar{\xi}_J = \mathbf{0}, \quad J = 1, \dots, 6M, \qquad (53) $$

where the $6M \times 6M$ matrix $\bar{\mathbf{N}} = \{\hat{\mathbf{N}}^{(mn)}\}$, $m, n = -N, \dots, N$, consists of the $6 \times 6$ blocks

$$ \bar{\mathbf{N}} = \{\hat{\mathbf{N}}^{(mn)}\}, \quad \hat{\mathbf{N}}^{(mn)} = \begin{pmatrix} -(\tilde{\mathbf{T}}^{-1})^{(ml)}(\tilde{\mathbf{R}}^{+})^{(ln)} & -(\tilde{\mathbf{T}}^{-1})^{(mn)} \\ \hat{\mathbf{P}}^{(mn)} - \hat{\mathbf{R}}^{(ml)}(\tilde{\mathbf{T}}^{-1})^{(lk)}(\hat{\mathbf{R}}^{+})^{(kn)} & -\hat{\mathbf{R}}^{(ml)}(\tilde{\mathbf{T}}^{-1})^{(ln)} \end{pmatrix}, \qquad (54) $$

expressed in matrices (9), and the 6-component vectors $\hat{\xi}_j^{(n)} = (\hat{\mathbf{A}}_j^{(n)} \hat{\mathbf{B}}_j^{(n)})^T$ consist of the 3-vectors $\hat{\mathbf{A}}_j^{(n)}$, $\hat{\mathbf{B}}_j^{(n)}$ appearing in (17), (18). Identity (10) yields $\bar{\mathbf{N}} = \overline{\mathbf{T}\mathbf{N}^T\bar{\mathbf{T}}}$ where $\overline{\mathbf{T}} = \text{diag}\{\hat{\mathbf{T}}\}$ and $\hat{\mathbf{T}}$ has zero diagonal and identity off-diagonal $3 \times 3$ blocks.

### 6. Zero-order PWE formulation

#### 6.1. Partial modes

Setting to zero the variation of material properties relative to their mean values specifies Eq. (53) to the form

$$ (\bar{\mathbf{N}}_0 - k_{y(n,a)}^{(0)}) \mathbf{I} \bar{\xi}_{0(n,a)} = \mathbf{0}, \qquad (55) $$

where $\bar{\mathbf{N}}_0(\omega, K_x) = \text{diag}\{\hat{\mathbf{N}}^{(n,n)}\}$ is a $6M \times 6M$ block-diagonal matrix, the double index $(n, \alpha)$, in which $n = -N, \dots, N$ and $\alpha = 1, \dots, 6$, replaces $J = 1, \dots, 6M$ used so far, and the additional index "0" is a reference to the zero-order state. Each eigenvector

$$ \bar{\xi}_{0(n,a)} = (0 \cdots \hat{\xi}_{0\alpha}^{(n)} \cdots 0)^T \text{ with } \hat{\xi}_{0\alpha}^{(n)} = (\hat{\mathbf{A}}_{0\alpha}^{(n)} \hat{\mathbf{B}}_{0\alpha}^{(n)})^T \qquad (56) $$

contains a single non-zero 6-tuple $\hat{\xi}_{0\alpha}^{(n)}$ at the nth position, which is the eigenvector of the nth on-diagonal block $\hat{\mathbf{N}}^{(nn)}$ so that

$$ (\hat{\mathbf{N}}^{(nn)} - k_{y(n,a)}^{(0)}) \mathbf{I} \hat{\xi}_{0\alpha}^{(n)} = \mathbf{0}, \quad \det(\hat{\mathbf{N}}^{(nn)} - k_{y(n,a)}^{(0)}) \mathbf{I} = 0, \quad \alpha = 1, \dots, 6. \qquad (57) $$

The blocks $\hat{\mathbf{N}}^{(nn)}$ represent each the Stroh matrix of the zero-order homogeneous medium while they differ by the value of the horizontal wavenumber $k_n$. The characteristic polynomial of $\hat{\mathbf{N}}^{(nn)}$ is homogeneous of degree 2 in $\omega$ and $k_n$, hence $k_{y(n,a)}^{(0)}(\omega, K_x) = p_a(\nu)k_n$ where $\nu = \omega/k_n$. The Stroh matrix is real, hence so is $\bar{\mathbf{N}}_0$; therefore, identity (54) reduces to the form $\bar{\mathbf{N}}_0 = \overline{\mathbf{T}\mathbf{N}_0^T\bar{\mathbf{T}}}$. Given the common case of semisimple (diagonalisable) $\bar{\mathbf{N}}_0$, the matrix $\bar{\xi}_0$ of its eigenvectors satisfies the orthonormality relation

$$ \bar{\xi}_0^T \overline{\mathbf{T}} \bar{\xi}_0 = \mathbf{I} \Leftrightarrow \bar{\xi}_{0(n,a)}^T \overline{\mathbf{T}} \bar{\xi}_{0(m,\beta)} = \hat{\xi}_{0\alpha}^{(n)T} \hat{\mathbf{T}} \hat{\xi}_{0\beta}^{(m)} = \delta_{nm} \delta_{a\beta}, \qquad (58) $$

where $n, m = -N, \dots, N$ and $\alpha, \beta = 1, \dots, 6$.

For a homogeneous medium, the subsonic interval (the one where no bulk modes with real $k_y$ are admitted) is defined in terms of the trace velocity $v = \omega/k_x$ and is bounded by the first transonic state $v = \hat{v}$, at which the normal to the surface is tangent to the outmost bulk-wave slowness sheet, see [8]. Let $\hat{v}$ be such a value for the reference homogeneous material with the averaged density and stiffness $\hat{\rho}_0$ and $\hat{c}_{ijkl}^{(0)}$. Passing to the PWE formalism requires operating with frequency rather than velocity. Denote

$$ \hat{\omega}_n^{(0)} = \hat{v} |k_n|, \quad n = -N, \dots, N, \qquad (59) $$
so that the frequency range

$$ \omega < \hat{\omega}_0^{(0)} = \hat{v} |K_x| \qquad (60) $$

is the zero-order subsonic interval, where all eigenvalues of $\overline{\mathbf{N}}_0$ are complex and thus may be partitioned into two complex conjugated subsets

$$ k_{y(n,a)}^{(0)} = k_{y(n,a+3)}^{(0)*} (\mathrm{Im} \, k_{y(n,a)}^{(0)}) > 0, \quad \hat{\xi}_{0a}^{(n)} = \hat{\xi}_{0(a+3)}^{(n)*}, \quad \alpha = 1, 2, 3, \forall n \qquad (61) $$

(cf. (13)). The range $\omega > \hat{\omega}_0^{(0)}$ is supersonic in the sense that it allows for the bulk (propagative) partial modes with real-valued vertical wavenumbers $k_{y(n,a)}^{(0)}$. Given $K_x$ fixed and $\omega \le \hat{\omega}_m^{(0)}$, the values of $k_{y(n,a)}^{(0)}$ with $|n| > |m|$ are guaranteed to be complex, while one pair $k_{y(m,\beta)}^{(0)}$ and $k_{y(m,\beta+3)}^{(0)}$ with $\beta \in \{1, 2, 3\}$ becomes real-valued as the frequency exceeds $\hat{\omega}_m^{(0)}$. Thus, the number of bulk modes within the zero-order PWE packet grows with increasing frequency. In the following, we will refer to $\hat{\omega}_n^{(0)}$ as the transonic frequency and to $\hat{\omega}_0^{(0)}$ as the lowest one (both referred to zeroth order PWE as indicated by the subscript).

The union of eigenvalues $k_{y(n,a)}^{(0)}$ of $\overline{\mathbf{N}}_0$ admits degeneracies between the eigenvalues of the same block $\hat{\mathbf{N}}^{(nn)}$ or of the different blocks. The former is always the case at the transonic frequency $\hat{\omega}_n^{(0)}$, where the block $\hat{\mathbf{N}}^{(nn)}$ becomes non-semisimple (defective) due to the merging of a pair of its eigenvalues, which are real at the point $\hat{\omega}_n^{(0)}$ and on one side of it, and complex on the other. Commonly, there is a single degenerate pair, i.e.,

$$ k_{y(n,a)}^{(0)} = k_{y(n,a+3)}^{(0)} \equiv k_{y(n,\mathrm{deg})}^{(0)} \quad \text{at} \quad \omega = \hat{\omega}_n^{(0)} \quad (\mathrm{Im} \, k_{y(n,\mathrm{deg})}^{(0)}) = 0 \qquad (62) $$

for one of $\alpha \in \{1, 2, 3\}$ in (61). Let us set it to be $\alpha = 1$, so that

$$ (\hat{\mathbf{N}}^{(nn)} - k_{y(n,\mathrm{deg})}^{(0)}) \mathbf{I} \hat{\xi}_{01}^{(n)} = \mathbf{0}, \quad (\hat{\mathbf{N}}^{(nn)} - k_{y(n,\mathrm{deg})}^{(0)}) \mathbf{I} \hat{\xi}_{04}^{(n)} = \hat{\xi}_{01}^{(n)} \qquad (63) $$

and

$$ \hat{\xi}_{01}^{(n)T} \hat{\mathbb{P}} \hat{\xi}_{01}^{(n)} = 0, \quad \hat{\xi}_{04}^{(n)T} \hat{\mathbb{P}} \hat{\xi}_{04}^{(n)} = 0, \quad \hat{\xi}_{01}^{(n)T} \hat{\mathbb{P}} \hat{\xi}_{04}^{(n)} = 1, \qquad (64) $$

where the first equality follows from (54) and the two others result from the optional normalisation. By (63), both $\hat{\xi}_{01}^{(n)}$ and $\hat{\xi}_{04}^{(n)}$ are either real or pure imaginary once $k_{y(n,\mathrm{deg})}^{(0)}$ is real. A block $\hat{\mathbf{N}}^{(nn)}$ may also turn non-semisimple due to the degeneracy of complex eigenvalues $k_{y(n,a)}^{(0)}$ at $\omega/k_n < \hat{v}$, though this is a rather special case (for instance, it is ruled out in isotropic media; see [17] for more details). In turn, possible degeneracy between real eigenvalues $k_{y(n,a)}^{(0)}$ at some value of $\omega/k_n > \hat{v}$, which signals the so-called acoustic axis, keeps $\hat{\mathbf{N}}^{(nn)}$ and hence $\overline{\mathbf{N}}_0$ semisimple (diagonalisable). Finally, the matrix $\overline{\mathbf{N}}_0$ also remains semisimple in the relatively common case of degeneracy between the eigenvalues of two different blocks of $\overline{\mathbf{N}}_0$. Implications of the different types of degeneracy will be discussed in the Sections Sections 7.2 and 8.2.

## 6.2. Surface waves

According to (56), the $3M \times 3M$ blocks $\tilde{\Xi}_{0l}$ of the matrix $\tilde{\Xi}_0$ composed of the zeroth order eigenvectors $\hat{\xi}_{0a}^{(n)}$ by analogy with (14) have quasi-diagonal structure

$$ \tilde{\Xi}_{0l} = \mathrm{diag}\{\hat{\Xi}_i^{(nn)}\}, \quad i = 1, \dots, 4. \qquad (65) $$

Hence the equation for the zeroth order PWE surface waves at the traction-free boundary reads

$$ \det \tilde{\Xi}_{03} = \prod_{n=-N}^{N} \det \hat{\Xi}_3^{(nn)} = 0, \qquad (66) $$

where each $3 \times 3$ matrix $\hat{\Xi}_3^{(nn)}$ is composed of the columns $\hat{\mathbf{B}}_{0a}^{(n)}$, $\alpha = 1, 2, 3$, see (61). The union of roots of (66) considered within the Brillouin zone $|K_x| \le \pi/T$ may be visualised as a “ladder-like” sequence of segments

$$ \omega_{Rn}(K_x) = v_R |K_x + gn|, \quad n = -N, \dots, N, \qquad (67) $$

where $v_R < \hat{v}$ is the velocity of the subsonic Rayleigh wave (RW) propagating along the lateral axis X in the reference unperturbed medium. The displacement-traction field of RW is a superposition of three eigenvectors of $\hat{\mathbf{N}}^{(nn)}$ ($\omega, K_x$)

$$ \hat{\xi}_{0a}^{(n)} (\omega_{Rn}, K_x) = (\hat{\mathbf{A}}_{0(Rn,a)}^{(n)} \hat{\mathbf{B}}_{0(Rn,a)}^{(n)})^T, \quad \alpha = 1, 2, 3, \qquad (68) $$

which correspond to the triplet of eigenvalues with $\mathrm{Im} \, k_{y(n,a)}^{(0)} > 0$ (see (57)). The RW fulfils the traction-free boundary condition, hence

$$ \sum_{a=1}^{3} c_{(Rn,a)} \hat{\mathbf{B}}_{0(Rn,a)}^{(n)} = \hat{\Xi}_{3(Rn)}^{(nn)} \bar{c}_{Rn} = \mathbf{0}, \quad \det \hat{\Xi}_{3(Rn)}^{(nn)} = 0, \qquad (69) $$
so that the frequency range

$$ \omega < \hat{\omega}_0^{(0)} = \hat{v} |K_x| \qquad (60) $$

is the zero-order subsonic interval, where all eigenvalues of $\bar{\mathbf{N}}_0$ are complex and thus may be partitioned into two complex conjugated subsets

$$ k_{y(n,a)}^{(0)} = k_{y(n,a+3)}^{(0)*} (\mathrm{Im} \, k_{y(n,a)}^{(0)}) > 0, \quad \hat{\xi}_{0\alpha}^{(n)} = \hat{\xi}_{0(\alpha+3)}^{(n)*}, \quad \alpha = 1, 2, 3, \forall n \qquad (61) $$

(cf. (13)). The range $\omega > \hat{\omega}_0^{(0)}$ is supersonic in the sense that it allows for the bulk (propagative) partial modes with real-valued vertical wavenumbers $k_{y(n,\alpha)}^{(0)}$. Given $K_x$ fixed and $\omega \le \hat{\omega}_m^{(0)}$, the values of $k_{y(n,\alpha)}^{(0)}$ with $|n| > |m|$ are guaranteed to be complex, while one pair $k_{y(m,\beta)}^{(0)}$ and $k_{y(m,\beta+3)}^{(0)}$ with $\beta \in \{1, 2, 3\}$ becomes real-valued as the frequency exceeds $\hat{\omega}_m^{(0)}$. Thus, the number of bulk modes within the zero-order PWE packet grows with increasing frequency. In the following, we will refer to $\hat{\omega}_n^{(0)}$ as the transonic frequency and to $\hat{\omega}_0^{(0)}$ as the lowest one (both referred to zeroth order PWE as indicated by the subscript).

The union of eigenvalues $k_{y(n,a)}^{(0)}$ of $\bar{\mathbf{N}}_0$ admits degeneracies between the eigenvalues of the same block $\hat{\mathbf{N}}^{(nn)}$ or of the different blocks. The former is always the case at the transonic frequency $\hat{\omega}_n^{(0)}$, where the block $\hat{\mathbf{N}}^{(nn)}$ becomes non-semisimple (defective) due to the merging of a pair of its eigenvalues, which are real at the point $\hat{\omega}_n^{(0)}$ and on one side of it, and complex on the other. Commonly, there is a single degenerate pair, i.e.,

$$ k_{y(n,a)}^{(0)} = k_{y(n,a+3)}^{(0)} \equiv k_{y(n,\mathrm{deg})}^{(0)} \quad \text{at} \quad \omega = \hat{\omega}_n^{(0)} \quad (\mathrm{Im} \, k_{y(n,\mathrm{deg})}^{(0)}) = 0 \qquad (62) $$

for one of $\alpha \in \{1, 2, 3\}$ in (61). Let us set it to be $\alpha = 1$, so that

$$ (\hat{\mathbf{N}}^{(nn)} - k_{y(n,\mathrm{deg})}^{(0)}) \mathbf{I} \hat{\xi}_{01}^{(n)} = \mathbf{0}, \quad (\hat{\mathbf{N}}^{(nn)} - k_{y(n,\mathrm{deg})}^{(0)}) \mathbf{I} \hat{\xi}_{04}^{(n)} = \hat{\xi}_{01}^{(n)} \qquad (63) $$

and

$$ \hat{\xi}_{01}^{(n)T} \hat{\mathbb{P}} \hat{\xi}_{01}^{(n)} = 0, \quad \hat{\xi}_{04}^{(n)T} \hat{\mathbb{P}} \hat{\xi}_{04}^{(n)} = 0, \quad \hat{\xi}_{01}^{(n)T} \hat{\mathbb{P}} \hat{\xi}_{04}^{(n)} = 1, \qquad (64) $$

where the first equality follows from (54) and the two others result from the optional normalisation. By (63), both $\hat{\xi}_{01}^{(n)}$ and $\hat{\xi}_{04}^{(n)}$ are either real or pure imaginary once $k_{y(n,\mathrm{deg})}^{(0)}$ is real. A block $\hat{\mathbf{N}}^{(nn)}$ may also turn non-semisimple due to the degeneracy of complex eigenvalues $k_{y(n,a)}^{(0)}$ at $\omega/k_n < \hat{v}$, though this is a rather special case (for instance, it is ruled out in isotropic media; see [17] for more details). In turn, possible degeneracy between real eigenvalues $k_{y(n,a)}^{(0)}$ at some value of $\omega/k_n > \hat{v}$, which signals the so-called acoustic axis, keeps $\hat{\mathbf{N}}^{(nn)}$ and hence $\bar{\mathbf{N}}_0$ semisimple (diagonalisable). Finally, the matrix $\bar{\mathbf{N}}_0$ also remains semisimple in the relatively common case of degeneracy between the eigenvalues of two different blocks of $\bar{\mathbf{N}}_0$. Implications of the different types of degeneracy will be discussed in the Sections Sections 7.2 and 8.2.

## 6.2. Surface waves

According to (56), the $3M \times 3M$ blocks $\tilde{\Xi}_{0l}$ of the matrix $\tilde{\Xi}_0$ composed of the zeroth order eigenvectors $\hat{\xi}_{0\alpha}^{(n)}$ by analogy with (14) have quasi-diagonal structure

$$ \tilde{\Xi}_{0l} = \mathrm{diag}\{\hat{\Xi}_{l}^{(nn)}\}, \quad l = 1, \dots, 4. \qquad (65) $$

Hence the equation for the zeroth order PWE surface waves at the traction-free boundary reads

$$ \det \tilde{\Xi}_{03} = \prod_{n=-N}^{N} \det \hat{\Xi}_{3}^{(nn)} = 0, \qquad (66) $$

where each $3 \times 3$ matrix $\hat{\Xi}_3^{(nn)}$ is composed of the columns $\hat{\mathbf{B}}_{0\alpha}^{(n)}$, $\alpha = 1, 2, 3$, see (61). The union of roots of (66) considered within the Brillouin zone $|K_x| \le \pi/T$ may be visualised as a “ladder-like” sequence of segments

$$ \omega_{Rn}(K_x) = v_R |K_x + gn|, \quad n = -N, \dots, N, \qquad (67) $$

where $v_R < \hat{v}$ is the velocity of the subsonic Rayleigh wave (RW) propagating along the lateral axis X in the reference unperturbed medium. The displacement-traction field of RW is a superposition of three eigenvectors of $\hat{\mathbf{N}}^{(nn)}$ ($\omega, K_x$)

$$ \hat{\xi}_{0\alpha}^{(n)}(\omega_{Rn}, K_x) = (\hat{\mathbf{A}}_{0(Rn,\alpha)}^{(n)} \hat{\mathbf{B}}_{0(Rn,\alpha)}^{(n)})^T, \quad \alpha = 1, 2, 3, \qquad (68) $$

which correspond to the triplet of eigenvalues with $\mathrm{Im} \, k_{y(n,a)}^{(0)} > 0$ (see (57)). The RW fulfils the traction-free boundary condition, hence

$$ \sum_{a=1}^{3} c_{(Rn,a)} \hat{\mathbf{B}}_{0(Rn,a)}^{(n)} = \hat{\Xi}_{3(Rn)}^{(nn)} \bar{c}_{Rn} = \mathbf{0}, \quad \det \hat{\Xi}_{3(Rn)}^{(nn)} = 0, \qquad (69) $$
where

$$ \hat{\Xi}_{3(R_n)}^{(nn)} = \| \hat{\mathbf{B}}_{0(R_n,1)}^{(n)} \hat{\mathbf{B}}_{0(R_n,2)}^{(n)} \hat{\mathbf{B}}_{0(R_n,3)}^{(n)} \| \qquad (70) $$

denotes the nth on-diagonal block of the traction matrix $\tilde{\Xi}_{03}$ taken at $\omega_{Rn}$ and the vector $\tilde{c}_{(Rn)}$ consists of amplitudes $c_{(Rn,\alpha)}$ of the RW partial modes. From the orthonormality condition $\hat{\Xi}^{(nn)T} \overline{\mathbb{T}} \hat{\Xi}^{(nn)} = \mathbf{I}$ and the link $\hat{\Xi}_{1(Rn)}^{(nn)} = \hat{\Xi}_{2(Rn)}^{(nn)*}$, $\hat{\Xi}_{3(Rn)}^{(nn)} = \hat{\Xi}_{4(Rn)}^{(nn)*}$ (see (58) and (61)), it is seen that the overall displacement vector of the (subsonic) RW $\overline{\mathbf{A}}_{0(Rn)} = \hat{\Xi}_{1(Rn)}^{(nn)} \tilde{c}_{(Rn)}$ satisfies the equalities

$$ \overline{\mathbf{A}}_{0(Rn)}^{T} \hat{\Xi}_{3(Rn)}^{(nn)} = \tilde{c}_{Rn}^{T} (\Rightarrow \tilde{c}_{Rn}^{T} \tilde{c}_{Rn} = 0), \quad \overline{\mathbf{A}}_{0(Rn)}^{+} \hat{\Xi}_{3(Rn)}^{(nn)} = \mathbf{0}. \qquad (71) $$

Note that in the plane strain case, where the RW consists of two in-plane polarized modes (say, $\alpha = 1$ and 2), the vector $\tilde{c}_{Rn}$ takes the form $\tilde{c}_{Rn} = (1, \pm i, 0)^T$ [18].

The zeroth order boundary condition (65) may also be realised at the transonic frequency $\omega = \delta\omega_n^{(0)}$ provided that the transonic state $\hat{v} = \hat{\omega}_n^{(0)}/k_n$ is exceptional [8,9], i.e., that the eigenvector $\hat{\xi}_{01}^{(n)}$ fulfils (64)1 due to zero traction subvector $\hat{\mathbf{B}}_{01}^{(n)}$ and hence, by (56), the matrix $\tilde{\Xi}_{03}$ has a vanishing column $\hat{\mathbf{B}}_{0(n,1)} = \mathbf{0}$. According to the Barnett–Lothe theory [6–9], the subsonic Rayleigh wave in a homogeneous half-space may not exist iff the transonic state is exceptional.

# 7. Perturbation of partial modes

## 7.1. Non-degenerate case

Assume weak periodic variation (52) of material properties. Hereafter we denote the running PWE indices by $l$ and $m$, whereas $n$ is reserved for the Rayleigh-branch index (see (67)). All roots $k_{y(m,a)}^{(0)}$ of Eq. (57)2 are supposed to be distinct and their difference is much greater than $O(\epsilon)$. Consider Eq. (53) with the perturbed system matrix $\overline{\mathbb{N}} \approx \overline{\mathbb{N}}_0 + \delta\overline{\mathbb{N}}$, where $\delta\overline{\mathbb{N}} \sim O(\epsilon)$ is the leading-order truncation of $\Delta\overline{\mathbb{N}} \equiv \overline{\mathbb{N}} - \overline{\mathbb{N}}_0$ (the terms of the next orders in $\epsilon$ are due to $\tilde{T}^{-1}$). Since $\overline{\mathbb{N}}_0$ is block-diagonal and the perturbation $\delta\overline{\mathbb{N}} = \{\delta\tilde{\mathbb{N}}^{(lm)}\}$ is block-hollow (i.e. $\delta\tilde{\mathbb{N}}^{(mm)} = 0 \forall m$), it follows that the eigenvectors $\tilde{\xi}_{m,a} = \{\hat{\xi}_{m,a}^{(l)}\} = \tilde{\xi}_{0(m,a)} + \delta\tilde{\xi}_{m,a}$ have the following structure: the unperturbed part $\tilde{\xi}_{0(m,a)}$ contains a single non-zero sixtuple $\hat{\xi}_{m,a}^{(m)} = \hat{\xi}_{0a}^{(m)}$ (see (56)) and the perturbed part $\delta\tilde{\xi}_{m,a}$ consists of (sub)vectors $\hat{\xi}_{m,a}^{(l)} \sim O(\epsilon)$ with $l \neq m$. Each perturbation vector $\hat{\xi}_{m,a}^{(l)}$ may be decomposed into the linear superposition of the unperturbed vectors so that

$$ \hat{\xi}_{m,a}^{(l)} = \sum_{\beta=1}^{6} q_{\beta a}^{(lm)} \hat{\xi}_{0\beta}^{(l)}. \qquad (72) $$

The coefficients $q_{\beta a}^{(lm)}$ may be found from the first-order approximation of Eq. (53),

$$ (\overline{\mathbb{N}}_0 - k_{y(m,a)}^{(0)}) \mathbf{I} \delta \tilde{\xi}_{m,a} + (\delta \overline{\mathbb{N}} - \delta k_{y(m,a)} \mathbf{I}) \tilde{\xi}_{0(m,a)} = \mathbf{0}, \qquad (73) $$

via left-multiplying it by $\hat{\xi}_{0(l,\beta)}^T \overline{\mathbb{T}}$ successively for all $l \neq m$ and taking into account (58), whence

$$ q_{\beta a}^{(lm)}|_{l \neq m} = \frac{\delta N_{\beta a}^{(lm)}}{k_{y(m,a)}^{(0)} - k_{y(l,\beta)}^{(0)}} \sim O(\epsilon), \quad \alpha, \beta = 1, \dots, 6, \qquad (74) $$

$$ \delta N_{\beta a}^{(lm)} = \hat{\xi}_{0(l,\beta)}^T \overline{\mathbb{T}} (\delta \overline{\mathbb{N}}) \hat{\xi}_{0(m,a)}^{(m)} = \hat{\xi}_{0\beta}^{(l)T} \overline{\mathbb{T}} \delta \hat{\mathbb{N}}^{(lm)} \hat{\xi}_{0a}^{(m)} = \hat{\mathbf{A}}_{0\beta}^{(l)T} \delta \hat{\mathbf{A}}_{\beta a}^{(lm)} \hat{\mathbf{A}}_{0a}^{(m)} \\ \text{with } \delta \hat{\mathbf{A}}_{\beta a}^{(lm)} = \hat{\mathbf{A}}^{(lm)} + k_{y(m,a)}^{(0)} \hat{\mathbf{R}}^{(lm)} + k_{y(l,\beta)}^{(0)} \hat{\mathbf{R}}^{(ml)} + k_{y(l,\beta)}^{(0)} k_{y(m,a)}^{(0)} \hat{\mathbf{T}}^{(lm)}. $$

Note that, by (54), (61), $q_{\beta a}^{(lm)} = -q_{\gamma \sigma}^{(ml)*}$ for $|\beta - \sigma| = |\alpha - \gamma| = 3$ in the subsonic frequency interval $\omega < \hat{\omega}_0$. The coefficients $q_{\alpha \beta}^{(mm)}$ are indeterminate from Eq. (73); however, they will be shown below to not contribute to the sought surface wave evaluation (note that taking $q_{\alpha \beta}^{(mm)} = 0$ extends the normalisation of the type (58) to the first perturbation order).

Consider now the perturbed traction matrix $\tilde{\Xi}_3 = \tilde{\Xi}_{03} + \delta\tilde{\Xi}_3$, whose columns $\hat{\mathbf{B}}_{m,a}$, $\alpha = 1, 2, 3$, when taken in the subsonic interval, correspond to the evanescent partial modes. As discussed above, $\hat{\mathbf{B}}_{(m,a)} = \{\hat{\mathbf{B}}_{m,a}^{(l)}\}$ consists of subvectors

$$ \hat{\mathbf{B}}_{m,a}^{(m)} = \hat{\mathbf{B}}_{0a}^{(m)} \sim O(1), \quad \hat{\mathbf{B}}_{m,a}^{(l)}|_{l \neq m} = \sum_{\beta=1}^{6} q_{\beta a}^{(lm)} \hat{\mathbf{B}}_{0\beta}^{(l)} \sim O(\epsilon), \qquad (75) $$

i.e. $\tilde{\Xi}_{03}$ and $\delta\tilde{\Xi}_3$ are block-diagonal and block-hollow, respectively, and their blocks are as follows:

$$ \hat{\Xi}_{3}^{(mm)} = \| \hat{\mathbf{B}}_{01}^{(m)} \hat{\mathbf{B}}_{02}^{(m)} \hat{\mathbf{B}}_{03}^{(m)} \| \sim O(1), \quad \hat{\Xi}_{3}^{(lm)}|_{l \neq m} = \| \hat{\mathbf{B}}_{m,1}^{(l)} \hat{\mathbf{B}}_{m,2}^{(l)} \hat{\mathbf{B}}_{m,3}^{(l)} \| \sim O(\epsilon). \qquad (76) $$

Inserting (75)2 and splitting the sum in two yields

$$ \hat{\Xi}_{3}^{(lm)}|_{l \neq m} = \hat{\Xi}_{3}^{(ll)} \hat{\mathbf{q}}_{1}^{(lm)} + \hat{\Xi}_{4}^{(ll)} \hat{\mathbf{q}}_{2}^{(lm)}, \qquad (77) $$

where $\hat{\mathbf{q}}_{1}^{(lm)} = \{q_{\beta a}^{(lm)}\}$ and $\hat{\mathbf{q}}_{2}^{(lm)} = \{q_{\beta+3,\alpha}^{(lm)}\}$ are the 3 × 3 matrices with the elements $q_{\beta a}^{(lm)}$ and $q_{\beta+3,\alpha}^{(lm)}$, $\alpha, \beta = 1, 2, 3$, defined in (74). Recall that $\hat{\Xi}_{3}^{(ll)*} = \hat{\Xi}_{4}^{(ll)*} \forall l$ in the subsonic interval $\omega < \hat{\omega}_0^{(0)}$.
where

$$ \hat{\Xi}_{3(R_n)}^{(nn)} = \| \hat{\mathbf{B}}_{0(R_n,1)}^{(n)} \hat{\mathbf{B}}_{0(R_n,2)}^{(n)} \hat{\mathbf{B}}_{0(R_n,3)}^{(n)} \| \qquad (70) $$

denotes the nth on-diagonal block of the traction matrix $\tilde{\Xi}_{03}$ taken at $\omega_{Rn}$ and the vector $\tilde{c}_{(Rn)}$ consists of amplitudes $c_{(Rn,\alpha)}$ of the RW partial modes. From the orthonormality condition $\hat{\Xi}^{(nn)T} \overline{\mathbb{T}} \hat{\Xi}^{(nn)} = \mathbf{I}$ and the link $\hat{\Xi}_{1(Rn)}^{(nn)} = \hat{\Xi}_{2(Rn)}^{(nn)*}$, $\hat{\Xi}_{3(Rn)}^{(nn)} = \hat{\Xi}_{4(Rn)}^{(nn)*}$ (see (58) and (61)), it is seen that the overall displacement vector of the (subsonic) RW $\overline{\mathbf{A}}_{0(Rn)} = \hat{\Xi}_{1(Rn)}^{(nn)} \tilde{c}_{(Rn)}$ satisfies the equalities

$$ \overline{\mathbf{A}}_{0(Rn)}^{T} \hat{\Xi}_{3(Rn)}^{(nn)} = \tilde{c}_{Rn}^{T} (\Rightarrow \tilde{c}_{Rn}^{T} \tilde{c}_{Rn} = 0), \quad \overline{\mathbf{A}}_{0(Rn)}^{+} \hat{\Xi}_{3(Rn)}^{(nn)} = \mathbf{0}. \qquad (71) $$

Note that in the plane strain case, where the RW consists of two in-plane polarized modes (say, $\alpha = 1$ and 2), the vector $\tilde{c}_{Rn}$ takes the form $\tilde{c}_{Rn} = (1, \pm i, 0)^T$ [18].

The zeroth order boundary condition (65) may also be realised at the transonic frequency $\omega = \delta\omega_n^{(0)}$ provided that the transonic state $\hat{v} = \hat{\omega}_n^{(0)}/k_n$ is exceptional [8,9], i.e., that the eigenvector $\hat{\xi}_{01}^{(n)}$ fulfils (64)¹ due to zero traction subvector $\hat{\mathbf{B}}_{01}^{(n)}$ and hence, by (56), the matrix $\tilde{\Xi}_{03}$ has a vanishing column $\hat{\mathbf{B}}_{0(n,1)} = \mathbf{0}$. According to the Barnett–Lothe theory [6–9], the subsonic Rayleigh wave in a homogeneous half-space may not exist iff the transonic state is exceptional.

# 7. Perturbation of partial modes

## 7.1. Non-degenerate case

Assume weak periodic variation (52) of material properties. Hereafter we denote the running PWE indices by $l$ and $m$, whereas $n$ is reserved for the Rayleigh-branch index (see (67)). All roots $k_{y(m,a)}^{(0)}$ of Eq. (57)² are supposed to be distinct and their difference is much greater than $O(\epsilon)$. Consider Eq. (53) with the perturbed system matrix $\overline{\mathbb{N}} \approx \overline{\mathbb{N}}_0 + \delta\overline{\mathbb{N}}$, where $\delta\overline{\mathbb{N}} \sim O(\epsilon)$ is the leading-order truncation of $\Delta\overline{\mathbb{N}} \equiv \overline{\mathbb{N}} - \overline{\mathbb{N}}_0$ (the terms of the next orders in $\epsilon$ are due to $\overline{\mathbb{T}}^{-1}$). Since $\overline{\mathbb{N}}_0$ is block-diagonal and the perturbation $\delta\overline{\mathbb{N}} = \{\delta\overline{\mathbb{N}}^{(lm)}\}$ is block-hollow (i.e. $\delta\overline{\mathbb{N}}^{(mm)} = 0 \forall m$), it follows that the eigenvectors $\tilde{\xi}_{m,a} = \{\hat{\xi}_{m,a}^{(l)}\} = \tilde{\xi}_{0(m,a)} + \delta\tilde{\xi}_{m,a}$ have the following structure: the unperturbed part $\tilde{\xi}_{0(m,a)}$ contains a single non-zero sixtuple $\hat{\xi}_{m,a}^{(m)} = \hat{\xi}_{0a}^{(m)}$ (see (56)) and the perturbed part $\delta\tilde{\xi}_{m,a}$ consists of (sub)vectors $\hat{\xi}_{m,a}^{(l)} \sim O(\epsilon)$ with $l \neq m$. Each perturbation vector $\hat{\xi}_{m,a}^{(l)}$ may be decomposed into the linear superposition of the unperturbed vectors so that

$$ \hat{\xi}_{m,a}^{(l)} = \sum_{\beta=1}^{6} q_{\beta a}^{(lm)} \hat{\xi}_{0\beta}^{(l)}. \qquad (72) $$

The coefficients $q_{\beta a}^{(lm)}$ may be found from the first-order approximation of Eq. (53),

$$ (\overline{\mathbb{N}}_0 - k_{y(m,a)}^{(0)}) \mathbf{I} \delta \tilde{\xi}_{m,a} + (\delta \overline{\mathbb{N}} - \delta k_{y(m,a)} \mathbf{I}) \tilde{\xi}_{0(m,a)} = \mathbf{0}, \qquad (73) $$

via left-multiplying it by $\hat{\xi}_{0(l,\beta)}^T \overline{\mathbb{T}}$ successively for all $l \neq m$ and taking into account (58), whence

$$ q_{\beta a}^{(lm)}|_{l \neq m} = \frac{\delta N_{\beta a}^{(lm)}}{k_{y(m,a)}^{(0)} - k_{y(l,\beta)}^{(0)}} \sim O(\epsilon), \quad \alpha, \beta = 1, \dots, 6, \qquad (74) $$

$$ \delta N_{\beta a}^{(lm)} = \hat{\xi}_{0(l,\beta)}^T \overline{\mathbb{T}} (\delta \overline{\mathbb{N}}) \hat{\xi}_{0(m,a)} = \hat{\xi}_{0\beta}^{(l)T} \overline{\mathbb{T}} \delta \overline{\mathbb{N}}^{(lm)} \hat{\xi}_{0a}^{(m)} = \hat{\mathbf{A}}_{0\beta}^{(l)T} \delta \hat{\mathbf{A}}_{\beta a}^{(lm)} \hat{\mathbf{A}}_{0a}^{(m)} $$

$$ \text{with } \delta \hat{\mathbf{A}}_{\beta a}^{(lm)} = \hat{\mathbf{B}}^{(lm)} + k_{y(m,a)}^{(0)} \hat{\mathbf{R}}^{(lm)} + k_{y(l,\beta)}^{(0)} \hat{\mathbf{R}}^{(ml)} + k_{y(l,\beta)}^{(0)} k_{y(m,a)}^{(0)} \hat{\mathbf{T}}^{(lm)}. $$

Note that, by (54), (61), $q_{\beta a}^{(lm)} = -q_{\gamma\sigma}^{(ml)*}$ for $|\beta - \sigma| = |\alpha - \gamma| = 3$ in the subsonic frequency interval $\omega < \hat{\omega}_0$. The coefficients $q_{\alpha\beta}^{(mm)}$ are indeterminate from Eq. (73); however, they will be shown below to not contribute to the sought surface wave evaluation (note that taking $q_{\alpha\beta}^{(mm)} = 0$ extends the normalisation of the type (58) to the first perturbation order).

Consider now the perturbed traction matrix $\tilde{\Xi}_3 = \tilde{\Xi}_{03} + \delta\tilde{\Xi}_3$, whose columns $\hat{\mathbf{B}}_{m,a}$, $\alpha = 1, 2, 3$, when taken in the subsonic interval, correspond to the evanescent partial modes. As discussed above, $\hat{\mathbf{B}}_{(m,a)} = \{\hat{\mathbf{B}}_{m,a}^{(l)}\}$ consists of subvectors

$$ \hat{\mathbf{B}}_{m,a}^{(m)} = \hat{\mathbf{B}}_{0a}^{(m)} \sim O(1), \quad \hat{\mathbf{B}}_{m,a}^{(l)}|_{l \neq m} = \sum_{\beta=1}^{6} q_{\beta a}^{(lm)} \hat{\mathbf{B}}_{0\beta}^{(l)} \sim O(\epsilon), \qquad (75) $$

i.e. $\tilde{\Xi}_{03}$ and $\delta\tilde{\Xi}_3$ are block-diagonal and block-hollow, respectively, and their blocks are as follows:

$$ \hat{\Xi}_3^{(mm)} = \| \hat{\mathbf{B}}_{01}^{(m)} \hat{\mathbf{B}}_{02}^{(m)} \hat{\mathbf{B}}_{03}^{(m)} \| \sim O(1), \quad \hat{\Xi}_3^{(lm)}|_{l \neq m} = \| \hat{\mathbf{B}}_{m,1}^{(l)} \hat{\mathbf{B}}_{m,2}^{(l)} \hat{\mathbf{B}}_{m,3}^{(l)} \| \sim O(\epsilon). \qquad (76) $$

Inserting (75)² and splitting the sum in two yields

$$ \hat{\Xi}_3^{(lm)}|_{l \neq m} = \hat{\Xi}_3^{(ll)} \hat{\mathbf{q}}_1^{(lm)} + \hat{\Xi}_4^{(ll)} \hat{\mathbf{q}}_2^{(lm)}, \qquad (77) $$

where $\hat{\mathbf{q}}_1^{(lm)} = \{q_{\beta a}^{(lm)}\}$ and $\hat{\mathbf{q}}_2^{(lm)} = \{q_{\beta+3,\alpha}^{(lm)}\}$ are the 3 × 3 matrices with the elements $q_{\beta a}^{(lm)}$ and $q_{\beta+3,\alpha}^{(lm)}$, $\alpha, \beta = 1, 2, 3$, defined in (74). Recall that $\hat{\Xi}_3^{(ll)*} = \hat{\Xi}_4^{(ll)*} \forall l$ in the subsonic interval $\omega < \hat{\omega}_0^{(0)}$.
The wavenumber perturbation $\delta k_{y(m,a)}$ is of the order $O(\epsilon^2)$, which cannot be evaluated from (73). Since it does not contribute to the leading-order perturbation of the surface-wave dispersion but may be of interest otherwise, we omit the derivation and skip to the result:

$$ 
\begin{aligned}
\delta k_{y(n,a)} &= \sum_{m \neq n; \beta=1,...,6} \frac{\delta N_{\alpha\beta}^{(nm)} \delta N_{\beta\alpha}^{(mn)}}{k_{y(n,a)}^{(0)} - k_{y(m,\beta)}^{(0)}} = -\hat{\xi}_{0(n,a)}^{T} \left[ \sum_{m \neq n} \hat{\Pi} \delta \hat{N}^{(nm)} (\hat{N}^{(mm)} - k_{y(n,a)}^{(0)} \mathbf{I})^{-1} \delta \hat{N}^{(mn)} \right] \hat{\xi}_{0(n,a)} = \\
&= -\hat{\mathbf{A}}_{0(n,a)}^{T} \left\{ \sum_{m \neq n} \delta \hat{\mathbf{A}}^{(nm)} \left[ \hat{\mathbf{A}}^{(mm)} (k_{y(n,a)}^{(0)}) \right]^{-1} \delta \hat{\mathbf{A}}^{(mn)} \right\} \hat{\mathbf{A}}_{0(n,a)} \sim O(\epsilon^2),
\end{aligned}
\quad (78) $$

where $\hat{\mathbf{A}}^{(nm)}(k_y) = \hat{\mathbf{A}}^{(nm)} + k_y(\hat{\mathbf{R}}^{(mm)} + \hat{\mathbf{R}}^{(mm)T}) + k_y^2 \hat{\xi}_{0(n,a)}^{(nm)}$, which can be recognised as the zero-order acoustic tensor. Here and below the summation index $m \neq n$ is short for $m \in \{-N, ..., N\} \setminus \{n\}$. Note from (78) that if $k_{y(m,a)}^{(0)}$ is real, then its perturbation $\delta k_{y(m,a)}$ is also real; in other words, any $(m, \alpha)$th supersonic partial mode which is bulk (i.e. has real $k_y$) in the unperturbed state remains such under the perturbation. Thus, Eq. (78) confirms that the laterally periodic superlattices usually do not admit the "true" (non-leaky) surface waves with frequency in the supersonic range unless the extraordinary cases where the boundary condition is fulfilled by the incomplete set of partial modes. This is unlike the existence of surface waves in the upper stopbands of 2D-periodic structures. It is also noteworthy that Eq. (78) is invariant to replacing $\bar{N}$ with $\bar{N}^s (= \bar{\Pi}\bar{N}^T\bar{\Pi})$, which shows that the perturbation $\delta k_{y(m,a)}$ does not change on passing to the "inverted" periodicity profile (27).

## 7.2. Degeneracy with non-semisimple $\bar{N}_0$

Perturbation theory in the case of eigenvalue degeneracy is elaborated in many textbooks, and it comes down to the system of two linear homogeneous equations describing the splitting of the 2D eigensubspace corresponding to the (doubly) degenerate eigenvalue, see e.g. [19]. In particular, this is the case when two eigenvalues of different blocks of the zero-order quasi-diagonal matrix $\bar{N}_0$ fall into degeneracy; as a result, their perturbation $\delta k_y$ is of the order $O(\epsilon)$, which is by order of magnitude greater than (78).

The main interest for our purpose is the degeneracy at the transonic frequency $\hat{\omega}_n^{(0)}$ which renders the block $\hat{\mathbf{N}}^{(nm)}$ of $\bar{N}_0$ non-semisimple, see (63), (64). Consider the splitting of the degenerate eigensubspace of $\bar{N}_0$ under the combined perturbation by the block-hollow matrix $\delta \bar{N}$ of the order $\epsilon$ and the small frequency variation $\Delta\omega = \omega - \hat{\omega}_n^{(0)}$. Assume $|\Delta\omega| \sim \epsilon^2$ (while $K_x$ remains a free parameter). Confining to the leading order in this parameter, the result written in terms of the subsequently used traction vectors is

$$ \bar{\mathbf{B}}_{n,a} |_{\alpha=1,4} = C_{\alpha}(\bar{\mathbf{B}}_{0(n,1)} + \delta k_{y(n,a)} \bar{\mathbf{B}}_{0(n,4)} + \sum_{m \neq n, \beta=1,...,6} q_{\beta 1}^{(mn)} \bar{\mathbf{B}}_{0(m,\beta)}), \quad (79) $$

where the common factor $C_\alpha$ is $\pm 1$ or 1 as specified below, $\bar{\mathbf{B}}_{0(n,a)} = (0 \cdots \hat{\mathbf{B}}_{0a}^{(n)} \cdots 0)^T$ (see (56)), the sum in $m(\neq n)$ runs from $-N$ to $N$, and the coefficients $q_{\beta 1}^{(mn)}$ are defined in (74). The perturbation $\delta k_{y(n,a)}$ of the degenerate eigenvalue $k_{y(n,1)}^{(0)} = k_{y(n,4)}^{(0)}$ is approximated to the leading order as

$$ \delta k_{y(n,a)} = \pm \sqrt{\kappa k_n(\omega - \hat{\omega}_n^{(0)}) + \sum_{11}^{(n)}} \quad \alpha = 1, 4, \quad (80) $$

where

$$ \kappa = -2\hat{\rho}_0 \hat{\omega} \hat{\mathbf{A}}_{01}^{(n)T} \hat{\mathbf{A}}_{01}^{(n)} \quad (81) $$

is the curvature radius of the $(n, 1)$th slowness curve at the transonic state, and the quantity $\Sigma_{11}^{(n)} \sim O(\epsilon^2)$ is

$$ \Sigma_{11}^{(n)} = \sum_{m \neq n, \beta=1,...,6} \frac{\delta N_{1\beta}^{(nm)} \delta N_{\beta 1}^{(mn)}}{k_{y(n,1)}^{(0)} - k_{y(m,\beta)}^{(0)}} |_{\omega=\hat{\omega}_n^{(0)}} = -\hat{\xi}_{01}^{(n)T} \hat{\Pi} \sum_{m \neq n} [\delta \hat{\mathbf{N}}^{(nm)} (\hat{\mathbf{N}}^{(mm)} - k_{y(n,1)}^{(0)} \mathbf{I})^{-1} \delta \hat{\mathbf{N}}^{(mn)}]_{\xi_{01}^{(n)}}^{2}. \quad (82) $$

By (82), $\Sigma_{11}^{(n)}$ is real. Accordingly, setting $\delta k_{y(n,a)} = 0$ yields the perturbed value of the transonic frequency $\hat{\omega}_n$:

$$ \hat{\omega}_n = \hat{\omega}_n^{(0)} - \frac{\Sigma_{11}^{(n)}}{\kappa k_n} + O(\epsilon^4), \quad (83) $$

where $O(\epsilon^4)$ implies the next-order term omitted under the radical in (80). It is seen from (83) and (54) that $\Sigma_{11}^{(n)}$ is the same for the "direct" and "inverted" profiles and hence so is the perturbed transonic frequency, in agreement with the general statement of Part I (see Section 3).

Aiming at the non-leaky surface waves, we are mostly concerned with the subsonic frequency range $\omega < \hat{\omega}_0^{(0)}$ bounded by the first (lowest) transonic frequency, see (60). It corresponds to the outmost slowness curve, which is therefore convex, i.e. $\kappa > 0$, and hence, by (81), $\hat{\xi}_{01}^{(n)}$ is purely imaginary. In consequence, the common factor in (79) is $C_{1,4} = \pm 1$, which maintains the identity $\hat{\xi}_1^{(n)} = \hat{\xi}_4^{(n)*}$ at the subsonic frequency $\omega < \hat{\omega}_0^{(0)}$. Note that the right-hand side of Eq. (82) taken at $\omega < \hat{\omega}_0^{(0)}$ reduces to the form $\Sigma_{11}^{(0)} = 2 \operatorname{Re} \sum_{m \neq 0; \beta=1,2,3} \cdot$

The above equations equally apply to the case of subsequent transonic states, which admit a locally concave slowness curve with $\kappa < 0$ wherefrom $\hat{\xi}_{01}^{(n)}$ and $\hat{\xi}_{04}^{(n)}$ are real and $C_{1,4} = 1$. Moreover, these equations remain valid for any event of eigenvalue degeneracy
The wavenumber perturbation $\delta k_{y(m,a)}$ is of the order $O(\epsilon^2)$, which cannot be evaluated from (73). Since it does not contribute to the leading-order perturbation of the surface-wave dispersion but may be of interest otherwise, we omit the derivation and skip to the result:

$$ \delta k_{y(n,a)} = \sum_{m \neq n; \beta=1,...,6} \frac{\delta N_{\alpha\beta}^{(nm)} \delta N_{\beta\alpha}^{(mn)}}{k_{y(n,a)}^{(0)} - k_{y(m,\beta)}^{(0)}} = -\hat{\xi}_{0(n,a)}^T \left[ \sum_{m \neq n} \hat{\Pi} \delta \hat{N}^{(nm)}(\hat{N}^{(mm)} - k_{y(n,a)}^{(0)} \mathbf{I})^{-1} \delta \hat{N}^{(mn)} \right] \hat{\xi}_{0(n,a)} = -\hat{\mathbf{A}}_{0(n,a)}^T \left\{ \sum_{m \neq n} \delta \hat{\mathbf{A}}^{(nm)} \left[ \hat{\mathbf{A}}^{(mm)}(k_{y(n,a)}^{(0)}) \right]^{-1} \delta \hat{\mathbf{A}}^{(mn)} \right\} \hat{\mathbf{A}}_{0(n,a)} \sim O(\epsilon^2), \quad (78) $$

where $\hat{\mathbf{A}}^{(mm)}(k_y) = \hat{\mathbf{A}}^{(mm)} + k_y(\hat{\mathbf{R}}^{(mm)} + \hat{\mathbf{R}}^{(mm)T}) + k_y^2 \hat{\xi}_{0(m,n)}^{(mm)}$, which can be recognised as the zero-order acoustic tensor. Here and below the summation index $m \neq n$ is short for $m \in \{-N, ..., N\} \setminus \{n\}$. Note from (78) that if $k_{y(m,a)}^{(0)}$ is real, then its perturbation $\delta k_{y(m,a)}$ is also real; in other words, any $(m, \alpha)$th supersonic partial mode which is bulk (i.e. has real $k_y$) in the unperturbed state remains such under the perturbation. Thus, Eq. (78) confirms that the laterally periodic superlattices usually do not admit the "true" (non-leaky) surface waves with frequency in the supersonic range unless the extraordinary cases where the boundary condition is fulfilled by the incomplete set of partial modes. This is unlike the existence of surface waves in the upper stopbands of 2D-periodic structures. It is also noteworthy that Eq. (78) is invariant to replacing $\bar{N}$ with $\bar{N}^s (= \bar{\Pi}\bar{N}^T\bar{\Pi})$, which shows that the perturbation $\delta k_{y(m,a)}$ does not change on passing to the "inverted" periodicity profile (27).

## 7.2. Degeneracy with non-semisimple $\bar{N}_0$

Perturbation theory in the case of eigenvalue degeneracy is elaborated in many textbooks, and it comes down to the system of two linear homogeneous equations describing the splitting of the 2D eigensubspace corresponding to the (doubly) degenerate eigenvalue, see e.g. [19]. In particular, this is the case when two eigenvalues of different blocks of the zero-order quasi-diagonal matrix $\bar{N}_0$ fall into degeneracy; as a result, their perturbation $\delta k_y$ is of the order $O(\epsilon)$, which is by order of magnitude greater than (78).

The main interest for our purpose is the degeneracy at the transonic frequency $\hat{\omega}_n^{(0)}$ which renders the block $\hat{\mathbf{N}}^{(nm)}$ of $\bar{N}_0$ non-semisimple, see (63), (64). Consider the splitting of the degenerate eigensubspace of $\bar{N}_0$ under the combined perturbation by the block-hollow matrix $\delta \bar{N}$ of the order $\epsilon$ and the small frequency variation $\Delta\omega = \omega - \hat{\omega}_n^{(0)}$. Assume $|\Delta\omega| \sim \epsilon^2$ (while $K_x$ remains a free parameter). Confining to the leading order in this parameter, the result written in terms of the subsequently used traction vectors is

$$ \bar{\mathbf{B}}_{n,a} |_{\alpha=1,4} = C_{\alpha}(\bar{\mathbf{B}}_{0(n,1)} + \delta k_{y(n,a)} \bar{\mathbf{B}}_{0(n,4)} + \sum_{m \neq n, \beta=1,...,6} q_{\beta 1}^{(mn)} \bar{\mathbf{B}}_{0(m,\beta)}), \quad (79) $$

where the common factor $C_\alpha$ is $\pm 1$ or 1 as specified below, $\bar{\mathbf{B}}_{0(n,a)} = (\mathbf{0} \cdots \hat{\mathbf{B}}_{0a}^{(n)} \cdots \mathbf{0})^T$ (see (56)), the sum in $m(\neq n)$ runs from $-N$ to $N$, and the coefficients $q_{\beta 1}^{(mn)}$ are defined in (74). The perturbation $\delta k_{y(n,a)}$ of the degenerate eigenvalue $k_{y(n,1)}^{(0)} = k_{y(n,4)}^{(0)}$ is approximated to the leading order as

$$ \delta k_{y(n,a)} = \pm \sqrt{\kappa k_n(\omega - \hat{\omega}_n^{(0)}) + \sum_{11}^{(n)}} \quad \alpha = 1, 4, \quad (80) $$

where

$$ \kappa = -2\hat{\rho}_0 \hat{\omega} \hat{\mathbf{A}}_{01}^{(n)T} \hat{\mathbf{A}}_{01}^{(n)} \quad (81) $$

is the curvature radius of the $(n, 1)$th slowness curve at the transonic state, and the quantity $\Sigma_{11}^{(n)} \sim O(\epsilon^2)$ is

$$ \Sigma_{11}^{(n)} = \sum_{m \neq n, \beta=1,...,6} \frac{\delta N_{1\beta}^{(nm)} \delta N_{\beta 1}^{(mn)}}{k_{y(n,1)}^{(0)} - k_{y(m,\beta)}^{(0)}} |_{\omega=\hat{\omega}_n^{(0)}} = -\hat{\xi}_{01}^{(n)T} \hat{\Pi} \sum_{m \neq n} [\delta \hat{\mathbf{N}}^{(nm)}(\hat{\mathbf{N}}^{(mm)} - k_{y(n,1)}^{(0)} \mathbf{I})^{-1} \delta \hat{\mathbf{N}}^{(mn)}]_{\xi_{01}^{(n)}}^{\xi_{01}^{(n)}}. \quad (82) $$

By (82), $\Sigma_{11}^{(n)}$ is real. Accordingly, setting $\delta k_{y(n,a)} = 0$ yields the perturbed value of the transonic frequency $\hat{\omega}_n$:

$$ \hat{\omega}_n = \hat{\omega}_n^{(0)} - \frac{\Sigma_{11}^{(n)}}{\kappa k_n} + O(\epsilon^4), \quad (83) $$

where $O(\epsilon^4)$ implies the next-order term omitted under the radical in (80). It is seen from (83) and (54) that $\Sigma_{11}^{(n)}$ is the same for the "direct" and "inverted" profiles and hence so is the perturbed transonic frequency, in agreement with the general statement of Part I (see Section 3).

Aiming at the non-leaky surface waves, we are mostly concerned with the subsonic frequency range $\omega < \hat{\omega}_0^{(0)}$ bounded by the first (lowest) transonic frequency, see (60). It corresponds to the outmost slowness curve, which is therefore convex, i.e. $\kappa > 0$, and hence, by (81), $\hat{\xi}_{01}^{(n)}$ is purely imaginary. In consequence, the common factor in (79) is $C_{1,4} = \pm 1$, which maintains the identity $\hat{\xi}_1^{(n)} = \hat{\xi}_4^{(n)*}$ at the subsonic frequency $\omega < \hat{\omega}_0^{(0)}$. Note that the right-hand side of Eq. (82) taken at $\omega < \hat{\omega}_0^{(0)}$ reduces to the form $\Sigma_{11}^{(0)} = 2 \operatorname{Re} \sum_{m \neq 0; \beta=1,2,3} \cdot$

The above equations equally apply to the case of subsequent transonic states, which admit a locally concave slowness curve with $\kappa < 0$ wherefrom $\hat{\xi}_{01}^{(n)}$ and $\hat{\xi}_{04}^{(n)}$ are real and $C_{1,4} = 1$. Moreover, these equations remain valid for any event of eigenvalue degeneracy
such that renders the block $\hat{\mathbf{N}}^{(nm)}$ non-semisimple - including possibly at $\omega_0 < \hat{\omega}_n^{(0)}$, which is when the degenerate $k_{y(n,1)}^{(0)}$ and the eigenvector $\hat{\xi}_{01}^{(n)}$ appearing in (82) are complex, see Section 6.1.

Finally, a note on the shear horizontal (SH) wave is worthwhile. It can be shown that the perturbation of its traction is proportional to that of its vertical wavenumber and so they vanish simultaneously; hence the SH transonic state remains exceptional, i.e. laterally periodic modulation of properties cannot give rise to the SH surface (quasibulk) wave.

# 8. Quasi-Rayleigh surface waves

## 8.1. Non-degenerate case

Periodic variation of material properties couples the PWE partial modes and thus leads to the quasi-Rayleigh wave (qRW) dispersion branches. Given the case of weak modulation (52), the qRW frequency

$$ \omega_{qR}^{(n)}(K_x) = \omega_{Rn}(K_x) + \Delta\omega_{qR}^{(n)}(K_x) \quad (84) $$

is sought in the vicinity of the n-th reference branch of the RW “ladder” (67). Of major interest for us is the fundamental qRW dispersion branch with $n=0$, which arises from the lowest RW branch $\omega_{R0} = v_R |K_x|$. Provided that the difference between the RW velocity $v_R$ and the limiting velocity $\hat{v}$ is large enough relative to the perturbation (52), the whole branch $\omega_{qR}^{(0)}(K_x)$ is subsonic and hence real-valued, thus corresponding to a “true” surface wave. The onsets of the next-order branches $\omega_{qR}^{(\pm 1)}$ may also be real near the Brillouin zone edges; however, the upper branches $\omega_{qR}^{(n)}$ are generally complex-valued and describe leaky qRWs containing bulk (propagative) partial modes whose number increases with growing frequency, see the discussion below (61).

We proceed from a generic case where $\omega_{Rn}$ is a single root of the zero-order dispersion equation $\det \tilde{\Xi}_{03} = 0$, see (66). Imposing weak modulation (52) and taking $\omega = \omega_{Rn} + \Delta\omega$ close to $\omega_{Rn}$ casts the dispersion equation in the form $\det \tilde{\Xi}_3(\omega, K_x) = 0$, where the traction matrix $\tilde{\Xi}_3$ has the block structure described in (75)–(77) and its n-th diagonal block is nearly singular due to proximity of $\omega$ to $\omega_{Rn}$. Having the main focus on the fundamental qRW branch, it is advantageous to write the same dispersion equation via the $3M \times 3M$ impedance matrix $\tilde{Z} = i\tilde{\Xi}_3 \tilde{\Xi}_1^{-1}$ (26), i.e.

$$ \det \tilde{Z}(\omega, K_x) = 0, \quad (85) $$

since its left-hand side is assuredly real in the subsonic range. The perturbation theory derivation, employing Jacobi's formula and the fact that the adjugate of $\tilde{Z}$ is a scalar multiple of $\tilde{A}_{Rn}\tilde{A}_{Rn}^+$, yields

$$ \det \tilde{Z} = i D_0 \tilde{A}_{Rn}^+ \left[ \frac{\partial \tilde{\Xi}^{(nn)}}{\partial \omega} \bigg|_{\omega_{Rn}} \Delta\omega - \sum_{m \neq n} \tilde{\Xi}_{3Rn}^{(nn)} \tilde{\Xi}_{3Rn}^{(mm)-1} \tilde{\Xi}_{3Rn}^{(mn)} \right] \tilde{\Xi}_{1Rn}^{(nn)-1} \tilde{A}_{Rn} \equiv i D_0 (f_1 \Delta\omega - f_2), \quad (86) $$

where $D_0$ is a real positive constant (its explicit form is unimportant in the present context) and

$$ \tilde{\Xi}_{3(Rn)}^{(l,m)} = \| \hat{\mathbf{B}}_{m(Rn,1)}^{(l)}, \hat{\mathbf{B}}_{m(Rn,2)}^{(l)}, \hat{\mathbf{B}}_{m(Rn,3)}^{(l)} \|, \quad (87) $$

which extends notation (70) to the off-diagonal blocks of the matrix $\tilde{\Xi}_3(\omega_{Rn})$. Note that $f_1$ and $f_2$ must be purely imaginary in the subsonic range, where $\det \tilde{Z}$ is real. Following the same routine as in Section 7.1 and making use of (71), the first term on the right-hand side of (86) can be developed to the form

$$ f_1 = \bar{A}_{Rn}^+ \frac{\partial}{\partial \omega} \tilde{\Xi}_{03}^{(nn)} \bar{c}_{Rn} = \bar{A}_{Rn}^+ (\tilde{\Xi}_{3(Rn)}^{(nn)} \hat{s}_1^{(n)} + \tilde{\Xi}_{3(Rn)}^{(nn)*} \hat{s}_2^{(n)}) \bar{c}_{Rn} = \bar{c}_{Rn}^+ \hat{s}_2^{(n)} \bar{c}_{Rn}, \quad (88) $$

where $\hat{s}_2^{(n)}$ is the 3 × 3 matrix with the elements $s_{3+\beta,\alpha}^{(n)}$, $\alpha, \beta = 1, 2, 3$, defined as follows:

$$ \mathbf{s}_{3+\beta,\alpha}^{(n)} = \frac{S_{3+\beta,\alpha}^{(n)}}{k_{y(n,\alpha)}^{(0)} - k_{y(n,3+\beta)}^{(0)}}, \quad \mathbf{S}_{3+\beta,\alpha}^{(n)} = \hat{\xi}_{0,3+\beta}^{(n)T} \hat{\Pi} \frac{d\hat{\mathbf{N}}^{(nn)}}{d\omega} \hat{\xi}_{0\alpha}^{(n)} = -2\hat{\rho}_0 \omega_{Rn} \hat{\mathbf{A}}_{0(Rn,3+\beta)}^{(n)T} \hat{\mathbf{A}}_{0(Rn,\alpha)}^{(n)}, \quad (89) $$

By (58) and (61), $\hat{s}_2^{(n)} = -\hat{s}_2^{(n)+}$. Next, inserting (77) and again using (71) casts the second term as<fnref n="1" />

$$ \hat{\Omega}_n = \hat{\Xi}_{1Rn}^{(nn)+} \sum_{m \neq n} (\hat{\Xi}_{3Rn}^{(nm)} \hat{\Xi}_{3Rn}^{(mm)-1} \hat{\Xi}_{3Rn}^{(mn)}) = \sum_{m \neq n} \hat{\mathbf{q}}_2^{(nm)} (\hat{\mathbf{q}}_1^{(mn)} + \hat{\mathbf{\Gamma}}^{(m)} \hat{\mathbf{q}}_2^{(m,n)}), \quad (90) $$

where $\hat{\Gamma}^{(m)} = \hat{\Xi}_{3(Rn)}^{(mm)-1} \hat{\Xi}_{4(Rn)}^{(mm)} (= -\hat{\Gamma}^{(m)+} \forall m$ in the subsonic range). Thus, from (86), (88) and (90),

$$ \Delta\omega_{qR}^{(n)} = \frac{\bar{c}_{Rn}^+ \hat{\Omega}_n \bar{c}_{Rn}}{\bar{c}_{Rn}^+ \hat{s}_2^{(n)} \bar{c}_{Rn}} \sim O(\epsilon^2). \quad (91) $$

<footnote n="1">¹ Note that though an indeterminate projection of $\delta\hat{\xi}_{n,a}$ on $\hat{\xi}_{0(n,a)}$ mentioned in Section 7.1 may contribute to the on-diagonal blocks of the perturbed traction matrix $\tilde{\Xi}_3(\omega_{Rn})$, it keeps the $(nn)$th block singular and thus cannot affect the leading-order approximation formula (90).</footnote>
such that renders the block $\hat{\mathbf{N}}^{(nm)}$ non-semisimple - including possibly at $\omega_0 < \dot{\omega}_n^{(0)}$, which is when the degenerate $k_{y(n,1)}^{(0)}$ and the eigenvector $\hat{\xi}_{01}^{(n)}$ appearing in (82) are complex, see Section 6.1.

Finally, a note on the shear horizontal (SH) wave is worthwhile. It can be shown that the perturbation of its traction is proportional to that of its vertical wavenumber and so they vanish simultaneously; hence the SH transonic state remains exceptional, i.e. laterally periodic modulation of properties cannot give rise to the SH surface (quasibulk) wave.

# 8. Quasi-Rayleigh surface waves

## 8.1. Non-degenerate case

Periodic variation of material properties couples the PWE partial modes and thus leads to the quasi-Rayleigh wave (qRW) dispersion branches. Given the case of weak modulation (52), the qRW frequency

$$ \omega_{qR}^{(n)}(K_x) = \omega_{Rn}(K_x) + \Delta\omega_{qR}^{(n)}(K_x) \quad (84) $$

is sought in the vicinity of the nth reference branch of the RW “ladder” (67). Of major interest for us is the fundamental qRW dispersion branch with $n = 0$, which arises from the lowest RW branch $\omega_{R0} = v_R |K_x|$. Provided that the difference between the RW velocity $v_R$ and the limiting velocity $\hat{v}$ is large enough relative to the perturbation (52), the whole branch $\omega_{qR}^{(0)}(K_x)$ is subsonic and hence real-valued, thus corresponding to a “true” surface wave. The onsets of the next-order branches $\omega_{qR}^{(\pm 1)}$ may also be real near the Brillouin zone edges; however, the upper branches $\omega_{qR}^{(n)}$ are generally complex-valued and describe leaky qRWs containing bulk (propagative) partial modes whose number increases with growing frequency, see the discussion below (61).

We proceed from a generic case where $\omega_{Rn}$ is a single root of the zero-order dispersion equation $\det \tilde{\Xi}_{03} = 0$, see (66). Imposing weak modulation (52) and taking $\omega = \omega_{Rn} + \Delta\omega$ close to $\omega_{Rn}$ casts the dispersion equation in the form $\det \tilde{\Xi}_3(\omega, K_x) = 0$, where the traction matrix $\tilde{\Xi}_3$ has the block structure described in (75)–(77) and its nth diagonal block is nearly singular due to proximity of $\omega$ to $\omega_{Rn}$. Having the main focus on the fundamental qRW branch, it is advantageous to write the same dispersion equation via the $3M \times 3M$ impedance matrix $\tilde{Z} = i\tilde{\Xi}_3\tilde{\Xi}_1^{-1}$ (26), i.e.

$$ \det \tilde{Z}(\omega, K_x) = 0, \quad (85) $$

since its left-hand side is assuredly real in the subsonic range. The perturbation theory derivation, employing Jacobi's formula and the fact that the adjugate of $\tilde{Z}$ is a scalar multiple of $\tilde{A}_{Rn}\tilde{A}_{Rn}^+$, yields

$$ \det \tilde{Z} = i D_0 \tilde{A}_{Rn}^+ \left[ \frac{\partial \tilde{\Xi}^{(nn)}}{\partial \omega} \bigg|_{\omega_{Rn}} \Delta\omega - \sum_{m \neq n} \tilde{\Xi}_{3Rn}^{(nn)} \tilde{\Xi}_{3Rn}^{(mm)-1} \tilde{\Xi}_{3Rn}^{(mn)} \right] \tilde{\Xi}_{1Rn}^{(nn)-1} \tilde{A}_{Rn} \equiv i D_0 (f_1 \Delta\omega - f_2), \quad (86) $$

where $D_0$ is a real positive constant (its explicit form is unimportant in the present context) and

$$ \tilde{\Xi}_{3(Rn)}^{(l,m)} = \| \hat{\mathbf{B}}_{m(Rn,1)}^{(l)}, \hat{\mathbf{B}}_{m(Rn,2)}^{(l)}, \hat{\mathbf{B}}_{m(Rn,3)}^{(l)} \|, \quad (87) $$

which extends notation (70) to the off-diagonal blocks of the matrix $\tilde{\Xi}_3(\omega_{Rn})$. Note that $f_1$ and $f_2$ must be purely imaginary in the subsonic range, where $\det \tilde{Z}$ is real. Following the same routine as in Section 7.1 and making use of (71), the first term on the right-hand side of (86) can be developed to the form

$$ f_1 = \tilde{A}_{Rn}^+ \frac{\partial}{\partial \omega} \tilde{\Xi}_{03}^{(nn)} \tilde{c}_{Rn} = \tilde{A}_{Rn}^+ (\tilde{\Xi}_{3(Rn)}^{(nn)} \hat{s}_1^{(n)} + \tilde{\Xi}_{3(Rn)}^{(nn)*} \hat{s}_2^{(n)}) \tilde{c}_{Rn} = \tilde{c}_{Rn}^+ \hat{s}_2^{(n)} \tilde{c}_{Rn}, \quad (88) $$

where $\hat{s}_2^{(n)}$ is the 3 × 3 matrix with the elements $s_{3+\beta,\alpha}^{(n)}$, $\alpha, \beta = 1, 2, 3$, defined as follows:

$$ s_{3+\beta,\alpha}^{(n)} = \frac{S_{3+\beta,\alpha}^{(n)}}{k_{y(n,\alpha)}^{(0)} - k_{y(n,3+\beta)}^{(0)}}, \quad S_{3+\beta,\alpha}^{(n)} = \hat{\xi}_{0,3+\beta}^{(n)T} \hat{\Pi} \frac{d\hat{\mathbf{N}}^{(nn)}}{d\omega} \hat{\xi}_{0\alpha}^{(n)} = -2\hat{\rho}_0 \omega_{Rn} \hat{\mathbf{A}}_{0(Rn,3+\beta)}^{(n)T} \hat{\mathbf{A}}_{0(Rn,\alpha)}^{(n)}, \quad (89) $$

By (58) and (61), $\hat{s}_2^{(n)} = -\hat{s}_2^{(n)+}$. Next, inserting (77) and again using (71) casts the second term as<fnref n="1" />

$$ f_2 = \tilde{c}_{Rn}^+ \hat{\Omega}_n \tilde{c}_{Rn} \quad \text{with} \quad (90) $$

$$ \hat{\Omega}_n = \hat{\Xi}_{1Rn}^{(nn)+} \sum_{m \neq n} (\hat{\Xi}_{3Rn}^{(nm)} \hat{\Xi}_{3Rn}^{(mm)-1} \hat{\Xi}_{3Rn}^{(mn)}) = \sum_{m \neq n} \hat{\mathbf{q}}_2^{(nm)} (\hat{\mathbf{q}}_1^{(mn)} + \hat{\mathbf{r}}^{(m)} \hat{\mathbf{q}}_2^{(m,n)}), $$

where $\hat{\mathbf{r}}^{(m)} = \hat{\Xi}_{3(Rn)}^{(mm)-1} \hat{\Xi}_{4(Rn)}^{(mm)} (= - \hat{\mathbf{r}}^{(m)+} \forall m$ in the subsonic range). Thus, from (86), (88) and (90),

$$ \Delta\omega_{qR}^{(n)} = \frac{\tilde{c}_{Rn}^+ \hat{\Omega}_n \tilde{c}_{Rn}}{\tilde{c}_{Rn}^+ \hat{s}_2^{(n)} \tilde{c}_{Rn}} \sim O(\epsilon^2). \quad (91) $$

<footnote n="1">¹ Note that though an indeterminate projection of $\delta \tilde{\xi}_{n,a}$ on $\tilde{\xi}_{0,n,a}$ mentioned in Section 7.1 may contribute to the on-diagonal blocks of the perturbed traction matrix $\tilde{\Xi}_3(\omega_{Rn})$, it keeps the $(nn)$th block singular and thus cannot affect the leading-order approximation formula (90).</footnote>
The fundamental qRW frequency branch ($n=0$) is described by Eq. (91) with purely imaginary numerator and denominator.

Consider now the qRW branches

$$ \omega'_{qR} (K_x) = \omega_{Rn} (K_x) + \Delta\omega'_{qR}^{(n)} (K_x) \qquad (92) $$

evolving from the RW ones in the superlattice with the “inverted” profile (27) and hence with the PWE system matrix $\overline{\mathbf{N}}^*$. According to Section 4, the surface-wave dispersion equation in this case is

$$ \det \tilde{Z}'(\omega, K_x) = 0 \qquad (93) $$

with the impedance matrix $\tilde{Z}' = -i\tilde{\Xi}_4\tilde{\Xi}_2^{-1}$ (32). The leading-order formula for the frequency shift $\Delta\omega'_{qR}^{(n)}$ has the same form as (91),

$$ \Delta\omega'_{qR}^{(n)} = \frac{\bar{c}_{Rn}^{+} \hat{\Omega}'_{n} \bar{c}_{Rn}}{\bar{c}_{Rn}^{+} \hat{s}_{2}^{(n)} \bar{c}_{Rn}} \sim O(\epsilon^2), \qquad (94) $$

up to replacing the matrix $\hat{\Omega}_n$ with an appropriate $\hat{\Omega}'_n$.

There may be two ways to evaluate the matrix $\hat{\Omega}'_n$. First is to repeat the same line of the perturbation-theory derivation as above but with respect to $\tilde{\Xi}_4$ in place of $\tilde{\Xi}_3$. The second way, which we opt for, is to simply reshuffle the already obtained expression of $\hat{\Omega}'_n$ (90)$_2$ taking note that it is quadratic in the perturbation coefficients $q_{\beta a}^{(lm)}$ (74) and these, in turn, are proportional to the blocks of matrix $\delta\overline{\mathbf{N}} = \{\delta\hat{\mathbf{N}}^{(lm)}\}$, whereas the sought $\hat{\Omega}'_n$ is quadratic in the counterpart coefficients $(q_{\beta a}^{(lm)})'$ proportional to the complex conjugate $\delta\overline{\mathbf{N}}^* = \{\delta\hat{\mathbf{N}}^{(lm)*}\}$. This may be written as

$$ q_{\beta a}^{(lm)} = q_{\beta a}^{(s:lm)} + q_{\beta a}^{(a:lm)}, \quad (q_{\beta a}^{(lm)})' = q_{\beta a}^{(s:lm)} - q_{\beta a}^{(a:lm)}, \qquad (95) $$

where

$$ q_{\beta a}^{(s:lm)} = \frac{\hat{\xi}_{0\beta}^{(l)T} \hat{\Pi}(\text{Re} \delta\hat{\mathbf{N}}^{(lm)}) \hat{\xi}_{0a}^{(m)}}{k_{y(m,a)}^{(0)} - k_{y(l,\beta)}^{(0)}}, \quad q_{\beta a}^{(a:lm)} = \frac{\hat{\xi}_{0\beta}^{(l)T} \hat{\Pi}(i \text{Im} \delta\hat{\mathbf{N}}^{(lm)}) \hat{\xi}_{0a}^{(m)}}{k_{y(m,a)}^{(0)} - k_{y(l,\beta)}^{(0)}}, \quad \alpha, \beta = 1, \dots, 6, \qquad (96) $$

and the indices "s" and "a" highlight the fact that the blocks $\delta\hat{\mathbf{N}}^{(lm)}$ are linear in the Fourier coefficients $\hat{\rho}_{l-m}$ and $\hat{c}_{ijkl}^{(l-m)}$ (see (74)), therefore $\text{Re} \delta\hat{\mathbf{N}}^{(lm)}$ and $\text{Im} \delta\hat{\mathbf{N}}^{(lm)}$ are gauged, respectively, by the symmetric and antisymmetric parts of the functions $\rho(x)$ and $c_{ijkl}(x)$. In these terms, the matrices $\hat{\Omega}_n$ and $\hat{\Omega}'_n$ may be subdivided as

$$ \begin{aligned} \hat{\Omega}_n &= \hat{\Omega}_n^{(ss)} + \hat{\Omega}_n^{(aa)} + \hat{\Omega}_n^{(as)} + \hat{\Omega}_n^{(sa)}, \\ \hat{\Omega}'_n &= \hat{\Omega}_n^{(ss)} + \hat{\Omega}_n^{(aa)} - (\hat{\Omega}_n^{(as)} + \hat{\Omega}_n^{(sa)}), \end{aligned} \qquad (97) $$

where

$$ \hat{\Omega}_n^{(ij)} = \sum_{m \neq n} \hat{\mathbf{q}}_2^{(i:nm)} (\hat{\mathbf{q}}_1^{(j:mn)} + \hat{\mathbf{r}}_1^{(m)} \hat{\mathbf{q}}_2^{(j:mn)}), \quad i, j \equiv s, a, \qquad (98) $$

with $\hat{\mathbf{q}}_1^{(k:lm)} = \{q_{\beta a}^{(k:lm)}\}$ and $\hat{\mathbf{q}}_2^{(k:lm)} = \{q_{\beta+3,a}^{(k:lm)}\}$, $k \equiv s, a; \alpha, \beta = 1, 2, 3$. Thus, by (91), (94) and (97), the difference between the leading-order frequencies $\omega_{qR}^{(n)}$ and $\omega'_{qR}^{(n)}$ of the quasi-RW in the superlattices with mutually "inverse" profiles is proportional to the quantity

$$ \hat{\Omega}_n - \hat{\Omega}'_n = 2(\hat{\Omega}_n^{(as)} + \hat{\Omega}_n^{(sa)}), \qquad (99) $$

where the right-hand side is explicitly defined in (98) and (95), (96).

By definition, $\omega_{qR}^{(n)} \equiv \omega'_{qR}^{(n)}$ if two profiles coincide, i.e. $\rho(x)$ and $c_{ijkl}(x)$ may be taken to have zero antisymmetric parts on a period. If, otherwise, they have zero symmetric parts on the period, i.e., $\hat{\Omega}_n = \hat{\Omega}'_n = \hat{\Omega}_n^{(aa)}$, then the qRW frequency values $\omega_{qR}^{(n)}$ and $\omega'_{qR}^{(n)}$ coincide in the leading order of the weak modulation parameter $\epsilon$, though differ in the next order.

## 8.2. Implications of zero-order degeneracy

Consider first the case where the statically averaged material coefficients $\hat{\rho}^{(0)}$, $\hat{c}_{ijkl}^{(0)}$ and the orientation of the axes X, Y (the geometry of the problem) are such that there exists a secluded value of the subsonic trace velocity $\omega/k_n < \hat{\nu}$, denoted by $= v_{\text{deg}}$, at which the block $\overline{\mathbf{N}}^{(nn)}$ of the system matrix $\overline{\mathbf{N}}$ turns non-semisimple and hence one of the zero-order attenuated ($\text{Im} k_{y(n,a)}^{(0)} > 0$) partial modes has a linear-times-exponential dependence on $y$. Perturbation of this solution by weak lateral periodicity contributes to the qRW solution solely near isolated frequencies that come about as the intersections of the RW “ladder” (67) with another “ladder” $\omega_{\text{deg}}^{(m)}(K_x) = v_{\text{deg}}|K_x + gm|$. There is a theoretical possibility of the equality $v_{\text{deg}} = v_R$, in which case the qRW involves the perturbation of this particular “linear-times-exponential mode” all along the dispersion branch $\omega_{qR}^{(n)}(K_x) = \omega_{Rn}(K_x) + \Delta\omega_{qR}^{(n)}(K_x)$, where $\Delta\omega_{qR}^{(n)}(K_x)$ is no longer described by Eq. (91) and should take into account the perturbation formulas of Section 7.1. It is of interest to study such a solution; however, as outlined above, this is an extraordinary situation, and therefore we will not pursue its study here.
The fundamental qRW frequency branch ($n=0$) is described by Eq. (91) with purely imaginary numerator and denominator.

Consider now the qRW branches

$$ \omega'_{qR} (K_x) = \omega_{Rn} (K_x) + \Delta\omega'_{qR} (K_x) \quad (92) $$

evolving from the RW ones in the superlattice with the “inverted” profile (27) and hence with the PWE system matrix $\overline{\mathbf{N}}^*$. According to Section 4, the surface-wave dispersion equation in this case is

$$ \det \tilde{Z}'(\omega, K_x) = 0 \quad (93) $$

with the impedance matrix $\tilde{Z}' = -i\tilde{\Xi}_4\tilde{\Xi}_2^{-1}$ (32). The leading-order formula for the frequency shift $\Delta\omega'_{qR} (n)$ has the same form as (91),

$$ \Delta\omega'_{qR} = \frac{\bar{c}_{Rn}^{+} \hat{\Omega}'_{n} \bar{c}_{Rn}}{\bar{c}_{Rn}^{+} \hat{s}_{2}^{(n)} \bar{c}_{Rn}} \sim O(\epsilon^2), \quad (94) $$

up to replacing the matrix $\hat{\Omega}_n$ with an appropriate $\hat{\Omega}'_n$.

There may be two ways to evaluate the matrix $\hat{\Omega}'_n$. First is to repeat the same line of the perturbation-theory derivation as above but with respect to $\tilde{\Xi}_4$ in place of $\tilde{\Xi}_3$. The second way, which we opt for, is to simply reshuffle the already obtained expression of $\hat{\Omega}'_n$ (90)$_2$ taking note that it is quadratic in the perturbation coefficients $q_{\beta a}^{(lm)}$ (74) and these, in turn, are proportional to the blocks of matrix $\delta\overline{\mathbf{N}} = \{\delta\hat{\mathbf{N}}^{(lm)}\}$, whereas the sought $\hat{\Omega}'_n$ is quadratic in the counterpart coefficients $(q_{\beta a}^{(lm)})'$ proportional to the complex conjugate $\delta\overline{\mathbf{N}}^* = \{\delta\hat{\mathbf{N}}^{(lm)*}\}$. This may be written as

$$ q_{\beta a}^{(lm)} = q_{\beta a}^{(s:lm)} + q_{\beta a}^{(a:lm)}, \quad (q_{\beta a}^{(lm)})' = q_{\beta a}^{(s:lm)} - q_{\beta a}^{(a:lm)}, \quad (95) $$

where

$$ q_{\beta a}^{(s:lm)} = \frac{\hat{\xi}_{0\beta}^{(l)T} \hat{\Pi}(\text{Re} \delta\hat{\mathbf{N}}^{(lm)}) \hat{\xi}_{0a}^{(m)}}{k_{y(m,a)}^{(0)} - k_{y(l,\beta)}^{(0)}}, \quad q_{\beta a}^{(a:lm)} = \frac{\hat{\xi}_{0\beta}^{(l)T} \hat{\Pi}(i \text{Im} \delta\hat{\mathbf{N}}^{(lm)}) \hat{\xi}_{0a}^{(m)}}{k_{y(m,a)}^{(0)} - k_{y(l,\beta)}^{(0)}}, \quad \alpha, \beta = 1, \dots, 6, \quad (96) $$

and the indices "s" and "a" highlight the fact that the blocks $\delta\hat{\mathbf{N}}^{(lm)}$ are linear in the Fourier coefficients $\hat{\rho}_{l-m}$ and $\hat{c}_{ijkl}^{(l-m)}$ (see (74)), therefore $\text{Re} \delta\hat{\mathbf{N}}^{(lm)}$ and $\text{Im} \delta\hat{\mathbf{N}}^{(lm)}$ are gauged, respectively, by the symmetric and antisymmetric parts of the functions $\rho(x)$ and $c_{ijkl}(x)$. In these terms, the matrices $\hat{\Omega}_n$ and $\hat{\Omega}'_n$ may be subdivided as

$$ \begin{aligned} \hat{\Omega}_n &= \hat{\Omega}_n^{(ss)} + \hat{\Omega}_n^{(aa)} + \hat{\Omega}_n^{(as)} + \hat{\Omega}_n^{(sa)}, \\ \hat{\Omega}'_n &= \hat{\Omega}_n^{(ss)} + \hat{\Omega}_n^{(aa)} - (\hat{\Omega}_n^{(as)} + \hat{\Omega}_n^{(sa)}), \end{aligned} \quad (97) $$

where

$$ \hat{\Omega}_n^{(ij)} = \sum_{m \neq n} \hat{\mathbf{q}}_2^{(i:nm)} (\hat{\mathbf{q}}_1^{(j:mn)} + \hat{\mathbf{r}}_1^{(m)} \hat{\mathbf{q}}_2^{(j:mn)}), \quad i, j \equiv s, a, \quad (98) $$

with $\hat{\mathbf{q}}_1^{(k:lm)} = \{q_{\beta a}^{(k:lm)}\}$ and $\hat{\mathbf{q}}_2^{(k:lm)} = \{q_{\beta+3,a}^{(k:lm)}\}$, $k \equiv s, a; \alpha, \beta = 1, 2, 3$. Thus, by (91), (94) and (97), the difference between the leading-order frequencies $\omega_{qR}^{(n)}$ and $\omega'_{qR}^{(n)}$ of the quasi-RW in the superlattices with mutually "inverse" profiles is proportional to the quantity

$$ \hat{\Omega}_n - \hat{\Omega}'_n = 2(\hat{\Omega}_n^{(as)} + \hat{\Omega}_n^{(sa)}), \quad (99) $$

where the right-hand side is explicitly defined in (98) and (95), (96).

By definition, $\omega_{qR}^{(n)} \equiv \omega'_{qR}^{(n)}$ if two profiles coincide, i.e. $\rho(x)$ and $c_{ijkl}(x)$ may be taken to have zero antisymmetric parts on a period. If, otherwise, they have zero symmetric parts on the period, i.e., $\hat{\Omega}_n = \hat{\Omega}'_n = \hat{\Omega}_n^{(aa)}$, then the qRW frequency values $\omega_{qR}^{(n)}$ and $\omega'_{qR}^{(n)}$ coincide in the leading order of the weak modulation parameter $\epsilon$, though differ in the next order.

## 8.2. Implications of zero-order degeneracy

Consider first the case where the statically averaged material coefficients $\hat{\rho}^{(0)}$, $\hat{c}_{ijkl}^{(0)}$ and the orientation of the axes X, Y (the geometry of the problem) are such that there exists a secluded value of the subsonic trace velocity $\omega/k_n < \hat{\nu}$, denoted by $= v_{\text{deg}}$, at which the block $\overline{\mathbf{N}}^{(nn)}$ of the system matrix $\overline{\mathbf{N}}$ turns non-semisimple and hence one of the zero-order attenuated ($\text{Im} k_{y(n,a)}^{(0)} > 0$) partial modes has a linear-times-exponential dependence on $y$. Perturbation of this solution by weak lateral periodicity contributes to the qRW solution solely near isolated frequencies that come about as the intersections of the RW “ladder” (67) with another “ladder” $\omega_{\text{deg}}^{(m)}(K_x) = v_{\text{deg}}|K_x + gm|$. There is a theoretical possibility of the equality $v_{\text{deg}} = v_R$, in which case the qRW involves the perturbation of this particular “linear-times-exponential mode” all along the dispersion branch $\omega_{qR}^{(n)}(K_x) = \omega_{Rn}(K_x) + \Delta\omega_{qR}^{(n)}(K_x)$, where $\Delta\omega_{qR}^{(n)}(K_x)$ is no longer described by Eq. (91) and should take into account the perturbation formulas of Section 7.1. It is of interest to study such a solution; however, as outlined above, this is an extraordinary situation, and therefore we will not pursue its study here.
Another type of degeneracy mentioned in Section 7.2 is that between the eigenvalues of two different blocks of $\bar{N}_0$, i.e., $k_{y(m,\alpha)}^{(0)} = k_{y(l,\beta)}^{(0)}$ with $m \neq l$. It comes about along the curves $\omega_{\text{deg}}^{(m,l)}(K_x)$, each depending on $m, l$ as parameters. At the point $(\omega, K_x)$, where this curve intersects the RW “ladder” $\omega_{Rn}(K_x)$ (67) (normally $m, l \neq n$), the coupled $(m, \alpha)$th and $(l, \beta)$th modes remain of the same order under perturbation and hence the pair of blocks $\hat{\Xi}_3^{(ml)}$ and $\hat{\Xi}_3^{(lm)}$ of the perturbed PWE traction matrix $\tilde{\Xi}_3(\omega_{Rn})$ are of the order $O(1)$ (unlike the non-degenerate case, where all off-diagonal blocks of $\tilde{\Xi}_3$ are of the order $O(\epsilon)$). In consequence, formula (91) for the frequency shift $\Delta\omega_{\text{qR}}^{(n)}(K_x)$ near the secluded points $\omega_{Rn} = \omega_{\text{deg}}^{(m,l)}$ needs a slight modification by way of changing a few addends in the sum (90); however, this does not affect the order of magnitude of $\Delta\omega_{\text{qR}}^{(n)} \sim O(\epsilon^2)$. An important reservation, however, is that the above conclusions on the impact of the degeneracy of eigenvalues of different blocks of $\bar{N}_0$ on the qRW dispersion do not apply to the case where the structure possesses a symmetry plane orthogonal to the X axis and $K_x$ is taken at the centre and edges of the Brillouin zone, i.e. $K_x = 0, \pm \frac{1}{2}g$. Analysis of this case will be reported elsewhere.

## 9. Quasibulk waves

Let the transonic state $\hat{v} = \hat{\omega}_n^{(0)}/k_n$ in the reference homogeneous medium be of the exceptional type, which means that the traction vector of the degenerate partial mode $\alpha = 1$ is zero: $\hat{B}_{01} = 0$ (see Section 6.2). Assume weak periodic modulation (52). The perturbation theory derivation yields

$$ \det \tilde{Z} = i D_0 \hat{A}_{01}^{(n)T} [\delta \hat{\Xi}_3^{(nm)} - \sum_{m \neq n} \hat{\Xi}_3^{(nm)} \hat{\Xi}_3^{(mn)} - \hat{\Xi}_3^{(lm)}] \hat{\Xi}_1^{(nm)} \hat{A}_{01}^{(n)} \equiv i D_0 (g_1(\omega) - g_2), \quad (100) $$

where the real constant $D_0$ is specified below (see (108)) and

$$ g_1(\omega) = \delta k_{y(n,1)} = \sqrt{\kappa k_n(\omega - \hat{\omega}_n) + O(\epsilon^4)}, \quad g_2 = \sum_{m \neq n} (\overline{\mathbf{p}_1^{(nm)}} \cdot \overline{\mathbf{p}_2^{(mn)}} - \overline{\mathbf{p}_1^{(nm)T}} \hat{\Gamma}_m \overline{\mathbf{p}_1^{(nm)*}}) \sim O(\epsilon^2) \quad (101) \\ \text{with } \overline{\mathbf{p}_1^{(lk)}} \equiv \{q_{3\beta}^{(lk)}, \overline{\mathbf{p}_2^{(lk)}} \equiv \{q_{\beta 3}^{(lk)}\}, \beta = 1, 2, 3. $$

Here $\hat{\omega}_n$ is the frequency of the perturbed transonic state given by (83), the coefficients $q_{\beta\alpha}^{(lm)}|_{l \neq m}$ are defined in (74), and $\hat{\Gamma}_m \equiv \hat{\Xi}_3^{(nm)} - \hat{\Xi}_4^{(nm)}$ is taken at $\omega = \hat{\omega}_n^{(0)}$. Recall that we are primarily concerned with the subsonic range below the first transonic state $\omega < \hat{\omega}_0$ ($n=0$), which implies that $\det \tilde{Z}$ is real, $\kappa > 0$ and $\hat{A}_{01}^{(0)T} \hat{A}_{01}^{(0)} < 0$ (see (81)), $\text{Im} \delta k_{y(0,1)} > 0$, $\hat{\Xi}_4^{(mm)} = \hat{\Xi}_3^{(mm)*}$ (hence $\hat{\Gamma}_m = -\hat{\Gamma}_m^+$), and $q_{3\beta}^{(lk)} = q_{\beta+3,3}^{(kl)*}$, $q_{\beta 3}^{(lk)} = q_{3,\beta+3}^{(kl)*}$. Accordingly, $g_1$ and $g_2$ have pure imaginary values. We stick to this setting below.

By (100) and (101), provided that $\text{Im } g_2 > 0$, the dispersion equation $\det \tilde{Z} = 0$ defines the frequency $\omega_{qb} = \hat{\omega}_0 - O(\epsilon^4)$ of a weakly localized, or quasibulk, surface wave whose attenuation is measured by $\delta k_{y(n,1)} \sim \epsilon$. Note that the theory of quasibulk surface waves in homogeneous anisotropic media [20] gives the same asymptotic dependence on the perturbation parameter, which in that case is the small angle of deviation of the sagittal plane from its particular orientation in a crystal such that renders the transonic state exceptional.

Now consider the perturbed dispersion equation $\det \tilde{Z}' = 0$ for the "inverted" profile (27). By analogy with §8.1, let us use the linearity of the leading-order perturbation to the Fourier coefficients which allows presenting the coefficients (Eq. (101)) for the "direct" and "inverse" profiles of $\rho(x)$ and $c_{ijkl}(x)$ as superpositions of contributions of their symmetric and, with an opposite sign, antisymmetric parts, namely,

$$ \overline{\mathbf{p}}_{1,2}^{(lk)} = \overline{\mathbf{p}}_{1,2}^{(s:lk)} + \overline{\mathbf{p}}_{1,2}^{(a:lk)}, \quad (\overline{\mathbf{p}}_{1,2}^{(lk)})' = \overline{\mathbf{p}}_{1,2}^{(s:lk)} - \overline{\mathbf{p}}_{1,2}^{(a:lk)}. \quad (102) $$

From definition (74), it follows that $\overline{\mathbf{p}}_1^{(s:kl)} = -\overline{\mathbf{p}}_2^{(s:lk)}$ and $\overline{\mathbf{p}}_1^{(a:kl)} = \overline{\mathbf{p}}_2^{(a:lk)}$; hence

$$ \overline{\mathbf{p}}_1^{(kl)} = -(\overline{\mathbf{p}}_2^{(lk)})', \quad \overline{\mathbf{p}}_2^{(kl)} = -(\overline{\mathbf{p}}_1^{(lk)})'. \quad (103) $$

Substituting to (Eq. (101)) and bearing in mind $\hat{\Gamma}_m = -\hat{\Gamma}_m^+$ yields

$$ \det \tilde{Z}' = iD(g_1(\omega) - g_2'), \quad (104) $$

where $g_2'$ is related to $g_2$ of (100) as follows

$$ g_2' = g_2 + 4i \operatorname{Im} \sum_{m \neq n} \overline{\mathbf{p}_1^{(s:mn)T}} \hat{\Gamma}_m \overline{\mathbf{p}_1^{(a:mn)*}} \quad (105) $$

with $n=0$ for the subsonic range.

Thus, setting to zero (100) and (104) shows that the frequencies of the subsonic quasibulk wave $\omega_{qb} < \hat{\omega}_0$ and $\omega'_{qb} < \hat{\omega}_0$, which may come about for either the "direct" or "inverted" periodicity profile, satisfy the estimate

$$ \omega_{qb} - \hat{\omega}_0 = O(\epsilon^4) + \frac{g_2^2}{\kappa K_x}, \quad \omega'_{qb} - \hat{\omega}_0 = O(\epsilon^4) + \frac{g_2'^2}{\kappa K_x}, \quad (106) $$

where the term $O(\epsilon^4)$ is the same for both $\omega_{qb}$ and $\omega'_{qb}$, while $\frac{1}{\kappa K_x}(g_2^2 - g_2'^2) \sim \epsilon^4$ is their difference (more detailed evaluation requires next-order perturbation calculations).
Another type of degeneracy mentioned in Section 7.2 is that between the eigenvalues of two different blocks of $\bar{\mathbf{N}}_0$, i.e., $k_{y(m,\alpha)}^{(0)} = k_{y(l,\beta)}^{(0)}$ with $m \neq l$. It comes about along the curves $\omega_{\text{deg}}^{(m,l)}(K_x)$, each depending on $m, l$ as parameters. At the point $(\omega, K_x)$, where this curve intersects the RW “ladder” $\omega_{Rn}(K_x)$ (67) (normally $m, l \neq n$), the coupled $(m, \alpha)$th and $(l, \beta)$th modes remain of the same order under perturbation and hence the pair of blocks $\hat{\mathbf{\Xi}}_3^{(ml)}$ and $\hat{\mathbf{\Xi}}_3^{(lm)}$ of the perturbed PWE traction matrix $\hat{\mathbf{\Xi}}_3(\omega_{Rn})$ are of the order $O(1)$ (unlike the non-degenerate case, where all off-diagonal blocks of $\hat{\mathbf{\Xi}}_3$ are of the order $O(\epsilon)$). In consequence, formula (91) for the frequency shift $\Delta\omega_{\text{qR}}^{(n)}(K_x)$ near the secluded points $\omega_{Rn} = \omega_{\text{deg}}^{(m,l)}$ needs a slight modification by way of changing a few addends in the sum (90); however, this does not affect the order of magnitude of $\Delta\omega_{\text{qR}}^{(n)} \sim O(\epsilon^2)$. An important reservation, however, is that the above conclusions on the impact of the degeneracy of eigenvalues of different blocks of $\bar{\mathbf{N}}_0$ on the qRW dispersion do not apply to the case where the structure possesses a symmetry plane orthogonal to the X axis and $K_x$ is taken at the centre and edges of the Brillouin zone, i.e. $K_x = 0, \pm \frac{1}{2}g$. Analysis of this case will be reported elsewhere.

## 9. Quasibulk waves

Let the transonic state $\hat{v} = \hat{\omega}_n^{(0)}/k_n$ in the reference homogeneous medium be of the exceptional type, which means that the traction vector of the degenerate partial mode $\alpha = 1$ is zero: $\hat{\mathbf{B}}_{01} = \mathbf{0}$ (see Section 6.2). Assume weak periodic modulation (52). The perturbation theory derivation yields

$$ \det \tilde{\mathbf{Z}} = i D_0 \hat{\mathbf{A}}_{01}^{(n)T} [\delta \hat{\mathbf{\Xi}}_3^{(nm)} - \sum_{m \neq n} \hat{\mathbf{\Xi}}_3^{(nm)} \hat{\mathbf{\Xi}}_3^{(nm)-1} \hat{\mathbf{\Xi}}_3^{(nm)}] \hat{\mathbf{\Xi}}_1^{(nm)-1} \hat{\mathbf{A}}_{01}^{(n)} \equiv i D_0 (g_1(\omega) - g_2), \quad (100) $$

where the real constant $D_0$ is specified below (see (108)) and

$$ g_1(\omega) = \delta k_{y(n,1)} = \sqrt{\kappa k_n(\omega - \hat{\omega}_n) + O(\epsilon^4)}, \quad g_2 = \sum_{m \neq n} (\overline{\mathbf{p}}_1^{(nm)} \cdot \overline{\mathbf{p}}_2^{(nm)} - \overline{\mathbf{p}}_1^{(nm)T} \hat{\mathbf{\Gamma}}_m \overline{\mathbf{p}}_1^{(nm)*}) \quad (\sim O(\epsilon^2)) \\ \text{with } \overline{\mathbf{p}}_1^{(lk)} = \{q_{3\beta}^{(lk)}\}, \quad \overline{\mathbf{p}}_2^{(lk)} = \{q_{\beta 3}^{(lk)}\}, \quad \beta = 1, 2, 3. \quad (101) $$

Here $\hat{\omega}_n$ is the frequency of the perturbed transonic state given by (83), the coefficients $q_{\beta\alpha}^{(lm)}|_{l \neq m}$ are defined in (74), and $\hat{\mathbf{\Gamma}}_m \equiv \hat{\mathbf{\Xi}}_3^{(nm)-1} \hat{\mathbf{\Xi}}_4^{(nm)}$ is taken at $\omega = \hat{\omega}_n^{(0)}$. Recall that we are primarily concerned with the subsonic range below the first transonic state $\omega < \hat{\omega}_0$ ($n=0$), which implies that $\det \tilde{\mathbf{Z}}$ is real, $\kappa > 0$ and $\hat{\mathbf{A}}_{01}^{(0)T} \hat{\mathbf{A}}_{01}^{(0)} < 0$ (see (81)), $\text{Im} \delta k_{y(0,1)} > 0$, $\hat{\mathbf{\Xi}}_4^{(nm)} = \hat{\mathbf{\Xi}}_3^{(nm)*}$ (hence $\hat{\mathbf{\Gamma}}_m = -\hat{\mathbf{\Gamma}}_m^+$), and $q_{3\beta}^{(lk)} = q_{\beta+3,3}^{(kl)*}, q_{\beta 3}^{(lk)} = q_{3,\beta+3}^{(kl)*}$. Accordingly, $g_1$ and $g_2$ have pure imaginary values. We stick to this setting below.

By (100) and (101), provided that $\text{Im } g_2 > 0$, the dispersion equation $\det \tilde{\mathbf{Z}}' = 0$ defines the frequency $\omega_{qb} = \hat{\omega}_0 - O(\epsilon^4)$ of a weakly localized, or quasibulk, surface wave whose attenuation is measured by $\delta k_{y(n,1)} \sim \epsilon$. Note that the theory of quasibulk surface waves in homogeneous anisotropic media [20] gives the same asymptotic dependence on the perturbation parameter, which in that case is the small angle of deviation of the sagittal plane from its particular orientation in a crystal such that renders the transonic state exceptional.

Now consider the perturbed dispersion equation $\det \tilde{\mathbf{Z}}' = 0$ for the "inverted" profile (27). By analogy with §8.1, let us use the linearity of the leading-order perturbation to the Fourier coefficients which allows presenting the coefficients (Eq. (101)) for the "direct" and "inverse" profiles of $\rho(x)$ and $c_{ijkl}(x)$ as superpositions of contributions of their symmetric and, with an opposite sign, antisymmetric parts, namely,

$$ \overline{\mathbf{p}}_{1,2}^{(lk)} = \overline{\mathbf{p}}_{1,2}^{(s:lk)} + \overline{\mathbf{p}}_{1,2}^{(a:lk)}, \quad (\overline{\mathbf{p}}_{1,2}^{(lk)})' = \overline{\mathbf{p}}_{1,2}^{(s:lk)} - \overline{\mathbf{p}}_{1,2}^{(a:lk)}. \quad (102) $$

From definition (74), it follows that $\overline{\mathbf{p}}_1^{(s:kl)} = -\overline{\mathbf{p}}_2^{(s:lk)}$ and $\overline{\mathbf{p}}_1^{(a:kl)} = \overline{\mathbf{p}}_2^{(a:lk)}$; hence

$$ \overline{\mathbf{p}}_1^{(kl)} = -(\overline{\mathbf{p}}_2^{(lk)})', \quad \overline{\mathbf{p}}_2^{(kl)} = -(\overline{\mathbf{p}}_1^{(lk)})'. \quad (103) $$

Substituting to (Eq. (101)) and bearing in mind $\hat{\mathbf{\Gamma}}_m = -\hat{\mathbf{\Gamma}}_m^+$ yields

$$ \det \tilde{\mathbf{Z}}' = iD(g_1(\omega) - g_2'), \quad (104) $$

where $g_2'$ is related to $g_2$ of (100) as follows

$$ g_2' = g_2 + 4i \operatorname{Im} \sum_{m \neq n} \overline{\mathbf{p}}_1^{(s:nm)T} \hat{\mathbf{\Gamma}}_m \overline{\mathbf{p}}_1^{(a:nm)*} \quad (105) $$

with $n=0$ for the subsonic range.

Thus, setting to zero (100) and (104) shows that the frequencies of the subsonic quasibulk wave $\omega_{qb} < \hat{\omega}_0$ and $\omega'_{qb} < \hat{\omega}_0$, which may come about for either the “direct” or “inverted” periodicity profile, satisfy the estimate

$$ \omega_{qb} - \hat{\omega}_0 = O(\epsilon^4) + \frac{g_2^2}{\kappa K_x}, \quad \omega'_{qb} - \hat{\omega}_0 = O(\epsilon^4) + \frac{g_2'^2}{\kappa K_x}, \quad (106) $$

where the term $O(\epsilon^4)$ is the same for both $\omega_{qb}$ and $\omega'_{qb}$, while $\frac{1}{\kappa K_x}(g_2^2 - g_2'^2) \sim \epsilon^4$ is their difference (more detailed evaluation requires next-order perturbation calculations).
Let us show that the quasibulk wave may split off from the transonic state for one of the mutually inverse profiles but not for both simultaneously. The speculation is based on the sign analysis of $\det \tilde{Z}$ and $\det \tilde{Z}'$ taken at the (perturbed) transonic frequency $\omega = \dot{\omega}_0$ and on the general considerations of Part I. Appealing to (100) and (104), consider the leading-order evaluation

$$ \det \tilde{Z}(\dot{\omega}_0) = i D_0 g_2, \quad \det \tilde{Z}'(\dot{\omega}_0) = i D_0 g_2', \qquad (107) $$

where

$$ D_0 = \frac{1}{\hat{A}_{01}^{(0)T} \hat{A}_{01}^{(0)}} \prod_{(m,a) \neq (0,1)} \lambda_{0(m,a)} \qquad (108) $$

with $\hat{A}_{01}^{(0)T} \hat{A}_{01}^{(0)} < 0$ (see above) has the sign of the product of non-zero eigenvalues $\lambda_{0(m,a)}$ ($m = -N, \dots, N$; $\alpha = 2, 3$) of the zero-order quasi-diagonal impedance $\tilde{Z}_0(\dot{\omega}_0^{(0)}) = \text{diag}\{\hat{Z}^{(nn)}\}$. Recall that these eigenvalues of the block-diagonal matrices $\tilde{Z}_0 = (\hat{Z}^{(nn)})$ are monotonically decreasing functions of $\omega < \dot{\omega}_0^{(0)}$ positive at $\omega = 0$. One of the eigenvalues related to the block $\hat{Z}^{(00)}$ changes sign at the frequency $\omega_{R0} < \dot{\omega}_0^{(0)}$ of the subsonic Rayleigh wave (RW) frequency. Hence,

$$ \text{sign} D_0 \begin{cases} < 0, & \text{if the zero-order subsonic RW } \nexists. \\ > 0, & \text{if the zero-order subsonic RW } \exists. \end{cases} \qquad (109) $$

In turn, $g_2$ and $g_2'$ are quadratic in the perturbation coefficients (Eq. (102)) and hence may be partitioned accordingly as $\{g_2, g_2'\} = g_2^{(ss)} + g_2^{(aa)} \pm (g_2^{(sa)} + g_2^{(as)})$. Recall that they are pure imaginary. Consider available options for the signs of $\text{Im} g_2$ and $\text{Im} g_2'$.

First, assume that the periodicity profile is symmetric: $g_2 = g_2' = g_2^{(ss)}$. Then, as one of the results of Part I tells us, the analogue of the Barnett-Lothe uniqueness theorem holds true so that if qRW solution exists, then the quasibulk one is ruled out. Therefore, $\tilde{Z}$ and $\tilde{Z}'$ taken at $\dot{\omega}_0$ have one negative eigenvalue each, hence their determinants are negative, and so, according to (107) and (109), $\text{sign} \text{Im} g_2^{(ss)} > 0$ at $\dot{\omega}_0$. Next, we assume an asymmetric profile with pure antisymmetric modulation and invoke another result of Part I, which proves the impossibility of more than three surface waves for whichever profile. Regarding this case, we infer that if the unperturbed state admits the RW solution and hence the qRW one exists for both the "direct" and "inverted" profiles, then none of them admits a quasibulk subsonic wave (for otherwise, due to $g_2 = g_2' = g_2^{(aa)}$, both of them would admit it and then the overall number of subsonic waves would be four). Hence, once again $\det \tilde{Z} = \det \tilde{Z}' < 0$ and $\text{sign} \text{Im} g_2^{(aa)} > 0$ at $\dot{\omega}_0$. Finally, consider a generic pair of mutually "inverse" profiles with non-zero symmetric and antisymmetric parts. Let there exist a RW which gives rise to a perturbed qRW for each of the two profiles. In this case, as we have shown, $\text{sign}(\text{Im} g_2^{(ss)} + \text{Im} g_2^{(aa)})$ taken at $\dot{\omega}_0$ should be positive. Further, we note that had the quasibulk wave also existed for each of the two profiles, both $\det \tilde{Z}$ and $\det \tilde{Z}'$ would be positive at $\dot{\omega}_0$ (due to two negative eigenvalues). In view of (107) and (109), the latter means that $\text{Im}[(g_2^{(ss)} + g_2^{(aa)}) \pm (g_2^{(sa)} + g_2^{(as)})]$ should be negative at $\dot{\omega}_0$. However, as argued above, the term in the first parentheses should be positive. Therefore, the superposition of two terms within the brackets cannot become negative for both alternative signs. Thus, the foregoing assumption of two quasibulk surface waves occurring on top of the two quasi-Rayleigh waves is proven false. We have thereby exemplified the operation of the key result of Part I, which restricts the maximum total number of subsonic surface waves in a pair of mutually "inverse" profiles of lateral periodicity to three.

## 10. Conclusions

The problem of the existence and number of surface acoustic waves in laterally periodic superlattices of arbitrary anisotropy and period arrangement was treated. In Part I, it was shown that the upper bound on the possible number of surface wave branches (or, more precisely, of surface waves with different frequencies at the same Bloch wavenumber) admitted in an arbitrary superlattice is three and the same bound restricts their total number in two superlattices with mutually inverse patterns of periodic arrangements unless they coincide (the pattern is symmetric), in which case the surface wave solution is unique. The existence of at least one surface wave branch in either one or the other superlattice is guaranteed, provided there is no exceptional solution at the transonic frequency.

It is noteworthy that textually the same formulations apply to the number of surface waves per stopband in transversely periodic half-spaces with mutually inverse profiles such that $\rho'(y) = \rho(-y)$ and $c_{ijkl}(y) = c'_{ijkl}(y)$ [10], despite that the boundary value problems related to the lateral and transverse periodicity require dissimilar analytical techniques.

In Part II, the perturbation-theory derivation was performed for the case of weak periodic variation of properties. Explicit formulas were obtained for two quasi-Rayleigh branches, which evolve due to mutually "inverse" modulation profiles from the zero-order Rayleigh wave, and for one quasibulk surface wave, which evolves for one of the profiles from the transonic state provided it is of the so-called exceptional type. Given the modulation gauged by small parameter $|\epsilon| \ll 1$, the frequencies of the quasi-Rayleigh waves differ from the unperturbed Rayleigh wave frequency by the value of the order $O(\epsilon^2)$. The quasibulk wave attenuation is small of the order $O(\epsilon)$, and its frequency is by $O(\epsilon^4)$ less than the perturbed bulk-wave (transonic) threshold, which in turn differs by $O(\epsilon^2)$ from the unperturbed transonic frequency.

Further research may be envisaged in two principal directions. First, there is much room for numerical modelling: in particular, an example of three surface waves propagating in a given superlattice is yet to be provided; besides that, it is of interest to simulate different scenarios of how the arrival of the extra surface wave for one of the mutually "inverse" profiles eliminates the surface wave for the other profile so as to prevent the total number of waves exceeding three. Second, by analogy with the case of transversely periodic superlattices, see e.g. [21–23], the above methodology of the existence considerations for surface waves in laterally periodic elastic superlattices may straightforwardly be adapted to the case of interface waves and extended to the piezoelectric and optical superlattices.
Let us show that the quasibulk wave may split off from the transonic state for one of the mutually inverse profiles but not for both simultaneously. The speculation is based on the sign analysis of $\det \tilde{Z}$ and $\det \tilde{Z}'$ taken at the (perturbed) transonic frequency $\omega = \dot{\omega}_0$ and on the general considerations of Part I. Appealing to (100) and (104), consider the leading-order evaluation

$$ \det \tilde{Z}(\dot{\omega}_0) = i D_0 g_2, \quad \det \tilde{Z}'(\dot{\omega}_0) = i D_0 g_2', \qquad (107) $$

where

$$ D_0 = \frac{1}{\hat{A}_{01}^{(0)T} \hat{A}_{01}^{(0)}} \prod_{(m,a) \neq (0,1)} \lambda_{0(m,a)} \qquad (108) $$

with $\hat{A}_{01}^{(0)T} \hat{A}_{01}^{(0)} < 0$ (see above) has the sign of the product of non-zero eigenvalues $\lambda_{0(m,a)}$ ($m = -N, \dots, N$; $\alpha = 2, 3$) of the zero-order quasi-diagonal impedance $\tilde{Z}_0(\dot{\omega}_0^{(0)}) = \text{diag}\{\hat{Z}^{(nn)}\}$. Recall that these eigenvalues of the block-diagonal matrices $\tilde{Z}_0 = (\hat{Z}^{(nn)})$ are monotonically decreasing functions of $\omega < \dot{\omega}_0^{(0)}$ positive at $\omega = 0$. One of the eigenvalues related to the block $\hat{Z}^{(00)}$ changes sign at the frequency $\omega_{R0} < \dot{\omega}_0^{(0)}$ of the subsonic Rayleigh wave (RW) frequency. Hence,

$$ \text{sign} D_0 \begin{cases} < 0, & \text{if the zero-order subsonic RW } \nexists. \\ > 0, & \text{if the zero-order subsonic RW } \exists. \end{cases} \qquad (109) $$

In turn, $g_2$ and $g_2'$ are quadratic in the perturbation coefficients (Eq. (102)) and hence may be partitioned accordingly as $\{g_2, g_2'\} = g_2^{(ss)} + g_2^{(aa)} \pm (g_2^{(sa)} + g_2^{(as)})$. Recall that they are pure imaginary. Consider available options for the signs of $\text{Im} g_2$ and $\text{Im} g_2'$.

First, assume that the periodicity profile is symmetric: $g_2 = g_2' = g_2^{(ss)}$. Then, as one of the results of Part I tells us, the analogue of the Barnett-Lothe uniqueness theorem holds true so that if qRW solution exists, then the quasibulk one is ruled out. Therefore, $\tilde{Z}$ and $\tilde{Z}'$ taken at $\dot{\omega}_0$ have one negative eigenvalue each, hence their determinants are negative, and so, according to (107) and (109), sign $\text{Im}g_2^{(ss)} > 0$ at $\dot{\omega}_0$. Next, we assume an asymmetric profile with pure antisymmetric modulation and invoke another result of Part I, which proves the impossibility of more than three surface waves for whichever profile. Regarding this case, we infer that if the unperturbed state admits the RW solution and hence the qRW one exists for both the "direct" and "inverted" profiles, then none of them admits a quasibulk subsonic wave (for otherwise, due to $g_2 = g_2' = g_2^{(aa)}$, both of them would admit it and then the overall number of subsonic waves would be four). Hence, once again $\det \tilde{Z} = \det \tilde{Z}' < 0$ and $\text{sign} \text{Im}g_2^{(aa)} > 0$ at $\dot{\omega}_0$. Finally, consider a generic pair of mutually "inverse" profiles with non-zero symmetric and antisymmetric parts. Let there exist a RW which gives rise to a perturbed qRW for each of the two profiles. In this case, as we have shown, $\text{sign}(\text{Im}g_2^{(ss)} + \text{Im}g_2^{(aa)})$ taken at $\dot{\omega}_0$ should be positive. Further, we note that had the quasibulk wave also existed for each of the two profiles, both $\det \tilde{Z}$ and $\det \tilde{Z}'$ would be positive at $\dot{\omega}_0$ (due to two negative eigenvalues). In view of (107) and (109), the latter means that $\text{Im}[(g_2^{(ss)} + g_2^{(aa)}) \pm (g_2^{(sa)} + g_2^{(as)})]$ should be negative at $\dot{\omega}_0$. However, as argued above, the term in the first parentheses should be positive. Therefore, the superposition of two terms within the brackets cannot become negative for both alternative signs. Thus, the foregoing assumption of two quasibulk surface waves occurring on top of the two quasi-Rayleigh waves is proven false. We have thereby exemplified the operation of the key result of Part I, which restricts the maximum total number of subsonic surface waves in a pair of mutually "inverse" profiles of lateral periodicity to three.

## 10. Conclusions

The problem of the existence and number of surface acoustic waves in laterally periodic superlattices of arbitrary anisotropy and period arrangement was treated. In Part I, it was shown that the upper bound on the possible number of surface wave branches (or, more precisely, of surface waves with different frequencies at the same Bloch wavenumber) admitted in an arbitrary superlattice is three and the same bound restricts their total number in two superlattices with mutually inverse patterns of periodic arrangements unless they coincide (the pattern is symmetric), in which case the surface wave solution is unique. The existence of at least one surface wave branch in either one or the other superlattice is guaranteed, provided there is no exceptional solution at the transonic frequency.

It is noteworthy that textually the same formulations apply to the number of surface waves per stopband in transversely periodic half-spaces with mutually inverse profiles such that $\rho'(y) = \rho(-y)$ and $c_{ijkl}(y) = c'_{ijkl}(y)$ [10], despite that the boundary value problems related to the lateral and transverse periodicity require dissimilar analytical techniques.

In Part II, the perturbation-theory derivation was performed for the case of weak periodic variation of properties. Explicit formulas were obtained for two quasi-Rayleigh branches, which evolve due to mutually "inverse" modulation profiles from the zero-order Rayleigh wave, and for one quasibulk surface wave, which evolves for one of the profiles from the transonic state provided it is of the so-called exceptional type. Given the modulation gauged by small parameter $|\epsilon| \ll 1$, the frequencies of the quasi-Rayleigh waves differ from the unperturbed Rayleigh wave frequency by the value of the order $O(\epsilon^2)$. The quasibulk wave attenuation is small of the order $O(\epsilon)$, and its frequency is by $O(\epsilon^4)$ less than the perturbed bulk-wave (transonic) threshold, which in turn differs by $O(\epsilon^2)$ from the unperturbed transonic frequency.

Further research may be envisaged in two principal directions. First, there is much room for numerical modelling: in particular, an example of three surface waves propagating in a given superlattice is yet to be provided; besides that, it is of interest to simulate different scenarios of how the arrival of the extra surface wave for one of the mutually "inverse" profiles eliminates the surface wave for the other profile so as to prevent the total number of waves exceeding three. Second, by analogy with the case of transversely periodic superlattices, see e.g. [21–23], the above methodology of the existence considerations for surface waves in laterally periodic elastic superlattices may straightforwardly be adapted to the case of interface waves and extended to the piezoelectric and optical superlattices.
CRediT authorship contribution statement

A.L. Shuvalov: Conceptualization, Formal analysis, Writing – original draft.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

No data was used for the research described in the article.

Acknowledgements

The author is grateful to Prof. V. Gusev, Prof. A. Sarychev and Dr O. Poncelet for stimulating discussions.

References

[1] Y. Yang, C. Dejous, H. Hallil, Trends and applications of surface and bulk acoustic wave devices: A review, Micromachines 14 (2022) 43.

[2] A. Hagelauer, et al., From microwave acoustic filters to millimeter-wave operation and new applications, IEEE J. Microw. 3 (2022) 484–508.

[3] M.Z. Aslam, H. Zhang, V.S. Sreejith, M. Naghdi, S. Ju, Advances in the surface acoustic wave sensors for industrial applications: Potentials, challenges, and future directions: A review, Measurement 222 (2023) 113657.

[4] B.S. Damasceno, et al., Recent improvements on surface acoustic wave sensors based on graphenic nanomaterials, Mater. Sci. Semicond. Process. 167 (2023) 107811.

[5] C. Li, N. Chigarev, T. Thréard, et al., Generation and detection of coherent picosecond surface acousticwaves in superlattices, e-forum acusticum (2020), in: Generation and Detection of Coherent GHz Surface Phonons in Superlattices, APS March Meeting 2021, IEEE International Ultrasonics Symposium, 2021.

[6] D.M. Barnett, J. Lothe, Consideration of the existence of surface wave (Rayleigh wave) solutions in anisotropic elastic crystals, J. Phys. F: Met. Phys. 4 (1974) 671–686.

[7] J. Lothe, B.M. Barnett, On the existence of surface-wave solutions for anisotropic elastic half-spaces with free surface, J. Appl. Phys. 47 (1976) 428–433.

[8] P. Chadwick, G.D. Smith, Foundations of the theory of surface waves in anisotropic elastic materials, Adv. Appl. Mech. 17 (1977) 303–376.

[9] D.M. Barnett, J. Lothe, Free surface (rayleigh) waves in anisotropic elastic half-spaces: The surface impedance method, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 402 (1985) 135–152.

[10] A.N. Darinskii, A.L. Shuvalov, Surface acoustic waves on one-dimensional phononic crystals of general anisotropy: existence considerations, Phys. Rev. B 98 (2018) 024309.

[11] J. Vasseur, The plane wave expansion method, in: In the Fundamentals and Applications of Acoustic Metamaterials: From Seismic to Radio Frequency, vol. 1, 2019, Ch. 2.

[12] B. Djafari-Rouhani, A.A. Maradudin, R.F. Wallis, Rayleigh waves on a superlattice stratified normal to the surface, Phys. Rev. B 29 (1984) 6454.

[13] A.N. Norris, A.L. Shuvalov, A.A. Kutsenko, The matrix sign function for solving surface wave problems in homogeneous and laterally periodic elastic half-spaces, Wave Motion 50 (2013) 1239–1250.

[14] L.X. Hu, L.P. Liu, K. Bhattacharya, Existence of surface waves and band gaps in periodic heterogeneous half-spaces, J. Elasticity 107 (2012) 65–79.

[15] S. Serra, Spectral and computational analysis of block Toeplitz matrices having nonnegative definite matrix-valued generating functions, BIT 39 (1999) 152–175.

[16] W. Rudin, Principles of Mathematical Analysis, McGraw Hill, 1953.

[17] A.L. Shuvalov, On the theory of plane inhomogeneous waves in anisotropic elastic media, Wave Motion 34 (2001) 401–429; ibid. 36 (2002) 305, Erratum.

[18] V.I. Alshits, J. Lothe, Comments on the relation between surface wave theory and the theory of reflection, Wave Motion 3 (1981) 297–310.

[19] L. Brillouin, Wave Propagation in Periodic Structures, Dover Publications, 1953.

[20] A.N. Darinskii, Quasi-bulk Rayleigh waves in semi-infinite media of arbitrary anisotropy, Wave Motion 27 (1998) 79–93.

[21] A.N. Darinskii, A.L. Shuvalov, Existence of surface acoustic waves in one-dimensional piezoelectric phononic crystals of general anisotropy, Phys. Rev. B 99 (2019) 174305.

[22] A.N. Darinskii, A.L. Shuvalov, Surface electromagnetic waves in anisotropic superlattices, Phys. Rev. A 102 (2020) 033515.

[23] A.N. Darinskii, A.L. Shuvalov, Stoneley-type waves in anisotropic periodic superlattices, Ultrasonics 109 (2021) 106237.
CRediT authorship contribution statement

A.L. Shuvalov: Conceptualization, Formal analysis, Writing – original draft.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

No data was used for the research described in the article.

Acknowledgements

The author is grateful to Prof. V. Gusev, Prof. A. Sarychev and Dr O. Poncelet for stimulating discussions.

References

[1] Y. Yang, C. Dejous, H. Hallil, Trends and applications of surface and bulk acoustic wave devices: A review, Micromachines 14 (2022) 43.

[2] A. Hagelauer, et al., From microwave acoustic filters to millimeter-wave operation and new applications, IEEE J. Microw. 3 (2022) 484–508.

[3] M.Z. Aslam, H. Zhang, V.S. Sreejith, M. Naghdi, S. Ju, Advances in the surface acoustic wave sensors for industrial applications: Potentials, challenges, and future directions: A review, Measurement 222 (2023) 113657.

[4] B.S. Damasceno, et al., Recent improvements on surface acoustic wave sensors based on graphenic nanomaterials, Mater. Sci. Semicond. Process. 167 (2023) 107811.

[5] C. Li, N. Chigarev, T. Thréard, et al., Generation and detection of coherent picosecond surface acousticwaves in superlattices, e-forum acusticum (2020), in: Generation and Detection of Coherent GHz Surface Phonons in Superlattices, APS March Meeting 2021, IEEE International Ultrasonics Symposium, 2021.

[6] D.M. Barnett, J. Lothe, Consideration of the existence of surface wave (Rayleigh wave) solutions in anisotropic elastic crystals, J. Phys. F: Met. Phys. 4 (1974) 671–686.

[7] J. Lothe, B.M. Barnett, On the existence of surface-wave solutions for anisotropic elastic half-spaces with free surface, J. Appl. Phys. 47 (1976) 428–433.

[8] P. Chadwick, G.D. Smith, Foundations of the theory of surface waves in anisotropic elastic materials, Adv. Appl. Mech. 17 (1977) 303–376.

[9] D.M. Barnett, J. Lothe, Free surface (rayleigh) waves in anisotropic elastic half-spaces: The surface impedance method, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 402 (1985) 135–152.

[10] A.N. Darinskii, A.L. Shuvalov, Surface acoustic waves on one-dimensional phononic crystals of general anisotropy: existence considerations, Phys. Rev. B 98 (2018) 024309.

[11] J. Vasseur, The plane wave expansion method, in: In the Fundamentals and Applications of Acoustic Metamaterials: From Seismic to Radio Frequency, vol. 1, 2019, Ch. 2.

[12] B. Djafari-Rouhani, A.A. Maradudin, R.F. Wallis, Rayleigh waves on a superlattice stratified normal to the surface, Phys. Rev. B 29 (1984) 6454.

[13] A.N. Norris, A.L. Shuvalov, A.A. Kutsenko, The matrix sign function for solving surface wave problems in homogeneous and laterally periodic elastic half-spaces, Wave Motion 50 (2013) 1239–1250.

[14] L.X. Hu, L.P. Liu, K. Bhattacharya, Existence of surface waves and band gaps in periodic heterogeneous half-spaces, J. Elasticity 107 (2012) 65–79.

[15] S. Serra, Spectral and computational analysis of block Toeplitz matrices having nonnegative definite matrix-valued generating functions, BIT 39 (1999) 152–175.

[16] W. Rudin, Principles of Mathematical Analysis, McGraw Hill, 1953.

[17] A.L. Shuvalov, On the theory of plane inhomogeneous waves in anisotropic elastic media, Wave Motion 34 (2001) 401–429; ibid. 36 (2002) 305, Erratum.

[18] V.I. Alshits, J. Lothe, Comments on the relation between surface wave theory and the theory of reflection, Wave Motion 3 (1981) 297–310.

[19] L. Brillouin, Wave Propagation in Periodic Structures, Dover Publications, 1953.

[20] A.N. Darinskii, Quasi-bulk Rayleigh waves in semi-infinite media of arbitrary anisotropy, Wave Motion 27 (1998) 79–93.

[21] A.N. Darinskii, A.L. Shuvalov, Existence of surface acoustic waves in one-dimensional piezoelectric phononic crystals of general anisotropy, Phys. Rev. B 99 (2019) 174305.

[22] A.N. Darinskii, A.L. Shuvalov, Surface electromagnetic waves in anisotropic superlattices, Phys. Rev. A 102 (2020) 033515.

[23] A.N. Darinskii, A.L. Shuvalov, Stoneley-type waves in anisotropic periodic superlattices, Ultrasonics 109 (2021) 106237.