A boundary integral equation approach to computing eigenvalues of the Stokes operator

TL;DR

This paper introduces a boundary integral equation method for efficiently computing eigenvalues and eigenfunctions of the Stokes operator, especially in complex geometries, with theoretical justification and numerical validation.

Contribution

It provides a new BIE-based framework for Stokes eigenvalue problems, including theoretical analysis and practical numerical implementation for complex domains.

Findings

The BIE approach reduces problem dimension and avoids high-frequency pollution.

Numerical results confirm the method's accuracy and efficiency.

Theoretical analysis ensures the uniqueness and validity of the approach.

Abstract

The eigenvalues and eigenfunctions of the Stokes operator have been the subject of intense analytical investigation and have applications in the study and simulation of the Navier-Stokes equations. As the Stokes operator is a fourth-order operator, computing these eigenvalues and the corresponding eigenfunctions is a challenging task, particularly in complex geometries and at high frequencies. The boundary integral equation (BIE) framework provides robust and scalable eigenvalue computations due to (a) the reduction in the dimension of the problem to be discretized and (b) the absence of high frequency "pollution" when using a Green's function to represent propagating waves. In this paper, we detail the theoretical justification for a BIE approach to the Stokes eigenvalue problem on simply and multiply-connected planar domains, which entails a treatment of the uniqueness theory for oscillatory Stokes equations on exterior domains. Then, using well-established techniques for discretizing BIEs, we present numerical results which confirm the analytical claims of the paper and demonstrate the efficiency of the overall approach.

Figures

Tables

Table 1. Table 1 : Roots of 67 .

k
5.135622301840683
8.417244140399865
11.61984117214906
14.79595178235126

Equations

$$-\Delta \mathbf{u}+\nabla p$$

$$\nabla \cdot \mathbf{u}$$

$$\mathbf{u}=0\text{on}\Gamma .$$

$$\mathbf{\sigma }(\mathbf{u},p)=-p\mathbf{I}+2\mathbf{e}(\mathbf{u}),$$

$$e_{ij}(\mathbf{u})=\frac{1}{2}\left({\partial }_{x_j}{u}_i+{\partial }_{x_i}{u}_j\right).$$

$$\mathbf{t}=\mathbf{\sigma }\cdot \mathbf{\nu }.$$

$$A(\Omega )=\{(\mathbf{u},p)\text{ s.t. }\mathbf{u}\in \left(C^2(\Omega )\times C^2(\Omega )\right)\cap \left(C(\bar{\Omega })\times C(\bar{\Omega })\right),p\in C^1(\Omega )\cap C(\bar{\Omega })\},$$

$${\mathcal{L}}_x{\int }_{{\mathbb{R}}^2}G(\mathbf{x},\mathbf{y})f(\mathbf{y})d\mathbf{y}=f(\mathbf{x}).$$

$$\Delta (\Delta +k^2)u=0,$$

$$G^{\text{BH}}(\mathbf{x},\mathbf{y})=\frac{1}{k^2}\left(\frac{1}{2\pi }\log |\mathbf{x}-\mathbf{y}|+\frac{i}{4}{H}_0^{(1)}(k|\mathbf{x}-\mathbf{y}|)\right),$$

$$G^{\text{L}}(\mathbf{x},\mathbf{y})=\frac{1}{2\pi }\log |\mathbf{x}-\mathbf{y}|,$$

$$G^{\text{H}}(\mathbf{x},\mathbf{y})=-\frac{i}{4}{H}_0^{(1)}(k|\mathbf{x}-\mathbf{y}|).$$

$$\nabla p-\Delta \mathbf{u}-k^2\mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$\Delta G^{\text{L}}(\mathbf{x},\mathbf{y})={\delta }_{\mathbf{y}}(\mathbf{x}).$$

$$p=\nabla G^{\text{L}}(\mathbf{x},\mathbf{y})\cdot \mathbf{f}.$$

$$\mathbf{u}$$

$$=\left(-\Delta +\nabla \otimes \nabla \right)G^{\text{BH}}\mathbf{f}.$$

$$\mathbf{G}=-\mathbf{I}\Delta G^{\text{BH}}+\nabla \otimes \nabla G^{\text{BH}}$$

$${\sigma }_{ij}=-p{\delta }_{ij}+\left({\partial }_{x_j}{u}_i+{\partial }_{x_i}{u}_j\right),$$

$$T_{ij\ell }$$

$$+{\partial }_{x_i}\left(-\Delta G^{\text{BH}}{\delta }_{\ell j}+{\partial }_{x_{\ell }}\left({\partial }_{x_j}{G}^{\text{BH}}\right)\right).$$

$$\mathbf{S}[\mathbf{\mu }](\mathbf{x})={\int }_{\Gamma }\mathbf{G}(\mathbf{x},\mathbf{y})\mathbf{\mu }(\mathbf{y})dS(\mathbf{y}).$$

$$\mathbf{D}[\mathbf{\mu }](\mathbf{x})={\int }_{\Gamma }{\left({\mathbf{T}}_{\cdot ,\cdot ,\ell }(\mathbf{x},\mathbf{y}){\nu }_{\ell }(\mathbf{y})\right)}^{⊺}\mathbf{\mu }(\mathbf{y})dS(\mathbf{y}),$$

$${\left({\mathbf{T}}_{\cdot ,\cdot ,\ell }(\mathbf{x},\mathbf{y}){\nu }_{\ell }(\mathbf{y})\right)}^{⊺}\mathbf{\mu }(\mathbf{y})$$

$$+{\nabla }^{\perp }\left({\partial }_{\tau \tau }-{\partial }_{\nu \nu }\right)G^{\text{BH}}(\mathbf{x},\mathbf{y}){\mu }_{\tau }(\mathbf{y}).$$

$$\mathcal{S}[\mathbf{\mu }](\mathbf{x})={\int }_{\Gamma }\mathbf{G}(\mathbf{x},\mathbf{y})\mathbf{\mu }(\mathbf{y})dS(\mathbf{y})$$

$$\mathcal{D}[\mathbf{\mu }](\mathbf{x})=\text{P.V.}{\int }_{\Gamma }{\left({\mathbf{T}}_{\cdot ,\cdot ,\ell }(\mathbf{x},\mathbf{y}){\nu }_{\ell }(\mathbf{y})\right)}^{⊺}\mathbf{\mu }(\mathbf{y})dS(\mathbf{y}),$$

$$⟨\mathbf{f},\mathbf{g}⟩={\int }_{\Gamma }\mathbf{f}\cdot \mathbf{g}dS.$$

$${\mathcal{D}}^{⊺}[\mathbf{\mu }](\mathbf{x})=\text{P.V.}{\int }_{\Gamma }\left({\mathbf{T}}_{\cdot ,\cdot ,\ell }(\mathbf{x},\mathbf{y}){\nu }_{\ell }(\mathbf{x})\right)\mathbf{\mu }(\mathbf{y})dS(\mathbf{y}).$$

$$\underset{h\downarrow 0^{+}}{\lim }{\mathbf{\sigma }}_{\mathbf{S}}[\mathbf{\mu }]({\mathbf{x}}_0\pm h\mathbf{\nu }({\mathbf{x}}_0))\cdot \mathbf{\nu }({\mathbf{x}}_0)$$

$$\underset{h\downarrow 0^{+}}{\lim }\mathbf{D}[\mathbf{\mu }]({\mathbf{x}}_0\pm h\mathbf{\nu }({\mathbf{x}}_0))$$

$$\mathbf{S}[\mathbf{t}](\mathbf{x})-\mathbf{D}[\mathbf{u}](\mathbf{x})=\{\begin{array}{ll}\mathbf{u}(\mathbf{x})& \mathbf{x}\in \Omega \\ 0& \mathbf{x}\in E\end{array},$$

Full text

1 Introduction

The planar incompressible Stokes equations describe creeping flows in two dimensions. Let $\Omega\subset\mathbb{R}^{2}$ be a bounded domain with $C^{2}$ boundary denoted by $\Gamma$. The Stokes eigenvalue problem is to find values $k^{2}$ such that

[TABLE]

subject to boundary conditions, has a non-trivial solution $({\boldsymbol{u}},p)$. In this work, we consider the eigenvalue problem subject to the Dirichlet boundary condition,

[TABLE]

It is well known that the values $k^{2}$ are necessarily real and positive and that there is a countable collection of such values $0<k_{1}^{2}\leq k_{2}^{2}\leq\ldots\uparrow\infty$, counting multiplicities.

Remark 1**.**

When $k=i\alpha$, the differential equation 1 is known as the modified Stokes equation. As there appears to be no preferred name for the equation with real-valued $k$, we will refer to 1 as the oscillatory Stokes equation.

The eigenvalues (and eigenfunctions) of the Stokes operator have applications in the stability analysis of stationary solutions of the Navier–Stokes equations [40], in the study of decaying two dimensional turbulence [48], and as a trial basis for numerical simulations of the Navier–Stokes equations [6]. The eigenvalues and eigenfunctions of the Stokes operator are also the subject of intense analytical investigation [52, 51, 42, 10, 2, 35, 29, 1], especially as they relate to the eigenvalues and eigenfunctions of the Laplacian.

Historically, the Stokes eigenvalue problem serves as a common model problem for numerical eigenvalue analysis with a fourth order operator (here, the bi-Laplacian). Further, numerical simulation has long played an important role in the analyses cited above — both for computing the eigenvalues and eigenfunctions on domains of practical interest and in forming new conjectures.

Borrowing the language of [54], which concerns the eigenvalues of the Laplacian (also known as the membrane or “drum” problem), the numerical treatment of the Stokes eigenvalue problem can be divided into two basic approaches. The first class of methods directly discretize the differential operator, typically with a finite element basis, and the eigenvalues are found as the eigenvalues of the discrete system. The second class of methods reformulate the oscillatory Stokes equations as a boundary integral equation (BIE) which is discretized. The eigenvalues are then found by a nonlinear search for the values of $k$ where the BIE is not invertible.

There is a large body of research on the first class of methods for the Stokes eigenvalue problem. We do not seek to review this literature here, but point to [28, 46, 38, 8, 26, 13, 37, 25, 12] for some representative examples.

As noted in [54], integral equation based methods provide several advantages. Because the BIE is defined on the boundary alone, there is a reduction in the dimension of the domain to be discretized. This approach reduces the number of unknowns over finite element discretizations, and does not suffer from high-frequency “pollution” due to the large number of unknowns required to resolve the computation when $k$ is large [4] (which typically affects finite element approaches to determining the eigenvalues).

Further, Zhao and Barnett [54] show how to alleviate some of the costliness of the nonlinear optimization introduced by formulating the problem as an integral equation. The standard approach searches for “V”-shaped minima of the singular values of the BIE; see, for instance, [53]. Instead, Zhao and Barnett utilize the Fredholm determinant (see section 3.6) which, for certain BIEs, is an analytic function of $k$ with roots precisely when $k^{2}$ is an eigenvalue. The Fredholm determinant can be estimated using a Nyström discretization of the BIE [9, 54]. Then, the eigenvalues can be estimated efficiently by using high order root finding methods applied to the discretized determinant.

With the efficiency of the approach of [54] for the drum problem in mind, we develop an integral equation based method for the Stokes eigenvalue problem. This requires that a layer potential representation of the solution of 1 be given and that the resulting BIE is not invertible precisely when $k^{2}$ is an eigenvalue. The first requirement is straightforward to satisfy because well-known layer potential representations for the modified Stokes equation [43, 7, 27, 34] are directly applicable. Proving the invertibility of the associated operators away from the eigenvalues is a more involved task and forms the bulk of the theoretical component of this paper.

1.1 Relation to other work

While integral equation based methods for the related “buckling” eigenvalue problem (which is equivalent on simply connected domains [29]) have been considered previously, these typically relied on first-kind integral equation formulations of the underlying PDE, i.e. formulations in which the BIE operator is compact [31, 1]. This is unsatisfying from a numerical perspective, because the spectrum of a compact operator either contains zero or has zero as a limit point (by design, the representations in [31, 1] are not injective precisely when $k^{2}$ is an eigenvalue). This obscures the relation between the non-invertibility of discrete approximations of the operator and the eigenvalues; in particular, common measures of the “non-invertibility” of a matrix, like the smallest singular value or the determinant, converge rapidly to zero for all values of $k^{2}$ as the boundary is refined. The measure of whether $k^{2}$ is an approximate eigenvalue is then relative to the current grid for first kind formulations.

The classical single and double layer representations for oscillatory Stokes considered in this paper result in second kind equations, i.e. integral equations of the form $\mathcal{I}-\mathcal{K}_{k}$ where $\mathcal{K}_{k}$ is compact. Such equations have a more satisfying theory [47, 15, 33], which translates well to numerical implementation [3, 9, 19, 54]. The use of a second kind representation is standard for the drum problem [5, 54] and was used recently to compute the vibrating modes of thin, clamped plates [36].

1.2 Paper outline and contributions

The rest of this paper proceeds as follows. In section 2, we set the notation, provide some mathematical preliminaries, and review properties of single and double layer potentials for the oscillatory Stokes equations. Then, in section 3, we develop the necessary theory for proving the main results of this work (theorems 8 and 10), which show that the BIEs resulting from these layer potential representations are not invertible precisely when $k^{2}$ is an eigenvalue. These theoretical developments include a detailed discussion of the uniqueness of oscillatory Stokes boundary value problems in exterior domains. To the best of our knowledge, the invertibility and uniqueness results are new to the literature. Section 3.6 then outlines how the Fredholm determinant can be used in the oscillatory Stokes context. In section 4, we describe the numerical methods we use to discretize the BIEs and to perform determinant calculations. While the underlying methods are well-established, the combination of a high-order singular quadrature rule and a fast-direct method for determinant evaluations in a BIE framework is novel. At the moderate frequencies considered in this paper, we find that standard fast-direct solvers provide a reasonably efficient determinant evaluation. We also present numerical experiments which demonstrate some of the paper’s analytical claims as well as the effectiveness of the overall framework. Finally, we provide some concluding thoughts, describe plans for future research, and outline some open questions in section 5.

2 Mathematical Preliminaries

In this paper, vector- and tensor-valued quantities are denoted by bold letters (e.g. ${\boldsymbol{h}}$ and $\mathbf{T}$). Subscript indices of non-bold characters (e.g. $h_{j}$ or $T_{ij\ell}$) are used to denote the entries within a vector or tensor. We use the standard Einstein summation convention; i.e., there is an implied sum taken over the repeated indices of any term (e.g. the symbol $a_{j}b_{j}$ is used to represent the sum $\sum_{j}a_{j}b_{j}$). If ${\boldsymbol{x}}=(x_{1},x_{2})^{\intercal}$, then ${\boldsymbol{x}}^{\bot}=(-x_{2},x_{1})^{\intercal}$. Similarly, $\nabla^{\bot}=(-\partial_{x_{2}},\partial_{x_{1}})^{\intercal}$. Upper-case script characters (e.g. $\mathcal{K}$) are reserved for operators on Banach spaces, with $\mathcal{I}$ denoting the identity. Given a set $X$, we denote the closure of $X$ by $\overline{X}$.

For a velocity field ${\boldsymbol{u}}$ and pressure $p$, let ${\boldsymbol{\sigma}}({\boldsymbol{u}},p)$ denote the Cauchy stress tensor is given by

[TABLE]

where ${\boldsymbol{e}}({\boldsymbol{u}})$ is the strain tensor given by

[TABLE]

When it is clear from context, we will drop the dependence of ${\boldsymbol{\sigma}}$ on ${\boldsymbol{u}}$ and $p$. If $\Gamma$ is the boundary of a region $\Omega$ and ${\boldsymbol{\nu}}$ is the outward normal to $\Gamma$, the surface traction ${\boldsymbol{t}}$ on $\Gamma$ is the Neumann data, i.e.

[TABLE]

We seek solutions of 1 in the space

[TABLE]

where $\Omega$ is an open domain.

2.1 Green’s functions

Let $\mathcal{L}_{x}$ denote a linear differential operator. A fundamental solution $G({\boldsymbol{x}},{\boldsymbol{y}})$ of $\mathcal{L}_{x}$ satisfies the equation $\mathcal{L}_{x}G({\boldsymbol{x}},{\boldsymbol{y}})=\delta_{y}({\boldsymbol{x}})$ in the distributional sense, i.e. for sufficiently smooth $f$

[TABLE]

We consider here free-space Green’s functions, i.e. fundamental solutions which satisfy appropriate radiation conditions as $|{\boldsymbol{x}}-{\boldsymbol{y}}|\to\infty$. The Green’s function of the oscillatory biharmonic equation,

[TABLE]

is given by

[TABLE]

where $k$ is the Helmholtz parameter in the oscillatory biharmonic equation, and $H_{0}^{1}(r)$ is the Hankel function of the first kind of order zero. Note that this is a scaled difference of the Green’s function for Laplace, i.e.

[TABLE]

and the Green’s function for the Helmholtz equation

[TABLE]

2.2 The Fredholm Alternative

We require some standard results from the theory of Fredholm integral equations. Interested readers may consult [47, 15, 33], among others, for the relevant background.

We first recall some definitions. Let $X$ and $Y$ be Banach spaces with a non-degenerate bilinear form $\langle\cdot,\cdot\rangle:X\times Y\to\mathbb{C}$.

• •

Two operators ${\mathcal{A}}:X\to X$ and ${\mathcal{B}}:Y\to Y$ are adjoint operators if $\langle A\phi,\psi\rangle=\langle\phi,B\psi\rangle$ for every $\phi\in X$ and $\psi\in Y$.

• •

For an operator ${\mathcal{M}}:X\to X$, we can define the range $R({\mathcal{M}})$ as the set $\{\phi\in X:\exists\phi_{0}\textrm{ with }{\mathcal{M}}\phi_{0}=\phi\}$ and the null space $N({\mathcal{M}})$ as the set $\{\phi\in X:{\mathcal{M}}\phi=0\}$.

• •

An operator ${\mathcal{A}}$ is said to be compact if $\overline{{\mathcal{A}}V}$ is a compact set for any bounded subset $V\subset X$.

• •

Given a subspace $V\subset X$, we can define the subspace $V^{\perp}\subset Y$ as the set $V^{\perp}=\{\psi\in Y:\langle\phi,\psi\rangle=0\textrm{ for each }\phi\in V\}$, with the analogous definition for subspaces of $Y$.

Operators of the form $\mathcal{I}-{\mathcal{A}}$ have existence and uniqueness properties analogous to matrices. This is known as the Fredholm Alternative; we present the version provided in [15].

Theorem 1** (Fredholm Alternative [15]).**

Let $X$ and $Y$ be Banach spaces and $\langle\cdot,\cdot\rangle:X\times Y\to\mathbb{C}$ be a bilinear form. Suppose that ${\mathcal{A}}:X\to X$ and ${\mathcal{B}}:Y\to Y$ are compact adjoint operators. Then $\dim N(\mathcal{I}-{\mathcal{A}})=\dim N(\mathcal{I}-{\mathcal{B}})\in\mathbb{N}$, $R(\mathcal{I}-{\mathcal{A}})=N(\mathcal{I}-{\mathcal{B}})^{\perp}$, and $R(\mathcal{I}-{\mathcal{B}})=N(\mathcal{I}-{\mathcal{A}})^{\perp}$.

2.3 Properties of the oscillatory Stokes layer

potentials

Recall that, in the case $k=i\alpha$ for some real-valued $\alpha$, the oscillatory Stokes equations 1 are known as the modified Stokes equations and are of particular interest for their application to the analysis and numerical simulation of unsteady flow [43, 7, 27, 34]. The equations are well-studied in that setting and integral representations which lead to second kind integral equations have been developed. We review some of the relevant results here, translating to the oscillatory setting.

2.3.1 Oscillatory Stokeslets and stresslets

Consider the solution of 1 where a $\delta$-mass centered at ${\boldsymbol{y}}$ with strength ${\boldsymbol{f}}$ has been added to the right-hand side of 1, i.e.

[TABLE]

Recall that

[TABLE]

If we substitute (9) into (8) and take the divergence, we obtain

[TABLE]

We then have, formally,

[TABLE]

The tensor

[TABLE]

is then the analog of a Stokeslet [43] for (1).

A related object is the stresslet, which is defined in terms of the stress tensor of the velocity, pressure pair induced by a Stokeslet. For these tensors, we find that it is more convenient to express them in index notation with the Einstein index summing convention. Recall that the stress tensor ${\boldsymbol{\sigma}}$ is defined as

[TABLE]

where $\delta_{ij}$ is the standard Kronecker delta notation. The stresslet ${\bf T}$ is defined to be

[TABLE]

Let $u_{i}=G_{ij}f_{j}$ and $p=\partial_{x_{i}}G^{\textrm{L}}f_{i}$ be a solution of the Stokes equations induced by a Stokeslet. Then the corresponding stress tensor is given by $\sigma_{i\ell}=T_{ij\ell}f_{j}$.

2.3.2 Layer potentials

We now use the Stokeslet and stresslet to define the single and double layer potentials for the oscillatory Stokes problem. For ${\boldsymbol{x}}\in\mathbb{R}^{2}$, the single layer potential with density ${\boldsymbol{\mu}}$ is defined to be

[TABLE]

We use the notation ${\boldsymbol{\sigma}}_{\bf S}[{\boldsymbol{\mu}}]$ to denote the stress tensor of the single layer at any given point ${\boldsymbol{x}}\in\mathbb{R}^{2}\setminus\Gamma$.

For ${\boldsymbol{x}}\in\mathbb{R}^{2}\setminus\Gamma$, the double layer potential with density ${\boldsymbol{\mu}}$ is defined to be

[TABLE]

where ${\boldsymbol{\nu}}$ denotes the outward unit normal to the boundary. If we write ${\boldsymbol{\mu}}={\boldsymbol{\nu}}\mu_{\nu}+{\boldsymbol{\tau}}\mu_{\tau}$, where ${\boldsymbol{\tau}}={\boldsymbol{\nu}}^{\bot}$ is the positively oriented unit tangent to the curve, then we have

[TABLE]

Let ${\mathcal{S}}[{\boldsymbol{\mu}}]:C(\Gamma)\to C(\Gamma)$, and ${\mathcal{D}}[{\boldsymbol{\mu}}]:C(\Gamma)\to C(\Gamma)$ denote the restrictions of the layer potentials ${\bf S}[{\boldsymbol{\mu}}]$ and ${\bf D}[{\boldsymbol{\mu}}]$ on the boundary $\Gamma$, i.e. for ${\boldsymbol{x}}\in\Gamma$,

[TABLE]

and

[TABLE]

where the P.V. indicates that the integral is to be evaluated in the principal value sense.

For two vector valued functions ${\boldsymbol{f}}$ and ${\boldsymbol{g}}$ defined on $\Gamma$, consider the bilinear form

[TABLE]

The definition of the adjoint used throughout the paper will be the one induced by this form.

The adjoint of ${\mathcal{D}}$ with respect to the above bilinear form is of particular interest; and is given by

[TABLE]

In the following lemma, we review the limiting values of the layer potentials ${\bf S}$ and ${\bf D}$ on the boundary $\Gamma$.

Lemma 1** (Jump conditions).**

Suppose that $\Omega$ is a bounded region with a $C^{2}$ boundary $\Gamma$. Let ${\boldsymbol{\nu}}({\boldsymbol{x}})$ denote the outward pointing normal at ${\boldsymbol{x}}\in\Gamma$. Suppose that ${\boldsymbol{\mu}}\in C(\Gamma)$. Then ${\bf S}[{\boldsymbol{\mu}}]$ is continuous across $\Gamma$, and the exterior and interior limits of the surface traction of ${\bf D}[{\boldsymbol{\mu}}]$ are equal. Furthermore, for ${\boldsymbol{x}}_{0}\in\Gamma$,

[TABLE]

The above expressions are derived by noting that the leading order singularity of these integral kernels is the same as for the original Stokes case, so that the standard jump conditions for Stokes [30, 43] apply.

Lemma 2**.**

Suppose that $\Omega$ is a bounded region with a $C^{2}$ boundary $\Gamma$. Then the operators ${\mathcal{S}}$ and ${\mathcal{D}}$ defined above are compact operators on $C(\Gamma)\times C(\Gamma)$ and $\mathbb{L}^{2}(\Gamma)\times\mathbb{L}^{2}(\Gamma)$.

Compactness is proved by considering the asymptotic expansion of each kernel about ${\boldsymbol{x}}={\boldsymbol{y}}$ and noting that each is at most weakly singular.

2.3.3 Representation Theorem

In the following theorem, we sketch the proof of the equivalent of the Green’s identity for oscillatory Stokes setting, which is well-known in the Stokes and modified Stokes settings [43, 7, 34].

Theorem 2**.**

Let $\Omega$ be a bounded domain with $C^{2}$ boundary and let the pair $({\boldsymbol{u}},p)$ satisfy the oscillatory Stokes equations 1 in $\Omega$. Let ${\boldsymbol{t}}$ denote the surface traction associated with $({\boldsymbol{u}},p)$. Then

[TABLE]

where $E=\mathbb{R}^{2}\setminus\bar{\Omega}$ is the exterior of the domain.

Proof.

Suppose that ${\boldsymbol{x}}\in\Omega$. By the definitions of ${\bf G}$ and $G^{\textrm{L}}$, we have

[TABLE]

Applying Green’s identity and the divergence theorem, we obtain

[TABLE]

Substituting the definition of the PDE and applying the divergence theorem again, we obtain

[TABLE]

From the divergence theorem and the divergence-free properties of ${\boldsymbol{u}}$ and ${\bf G}$, we then get

[TABLE]

Adding 22 and 23, we get the desired result. The argument for the case ${\boldsymbol{x}}\in E$ is similar. ∎

A consequence of theorem 2 and the analyticity of $G^{\textrm{BH}}$ is

Corollary 1**.**

Let $({\boldsymbol{u}},p)\in A(\Omega)$ be a solution of 1. Then each component of ${\boldsymbol{u}}({\boldsymbol{x}})$ is an analytic function of the coordinates ${\boldsymbol{x}}$ in $\Omega$.

The proof of corollary 1 follows the same reasoning as that for the Helmholtz case; see [15, Theorem 3.5].

2.3.4 Null-space correction

Without modification, the standard layer potentials can result in rank-deficient representations for the boundary value problems. The nature of this deficiency is treated in section 3 but for now we introduce a standard operator used to correct this. For any integrable density ${\boldsymbol{\mu}}$, let ${\mathcal{W}}[{\boldsymbol{\mu}}]$ be defined by

[TABLE]

for any ${\boldsymbol{x}}\in\Gamma$. We have

Lemma 3**.**

Let $\Omega$ be a domain with $C^{2}$ boundary and ${\boldsymbol{\mu}}$ be an integrable function defined on $\Gamma$. Then

• •

${\mathcal{W}}[{\mathcal{W}}[{\boldsymbol{\mu}}]]={\mathcal{W}}[{\boldsymbol{\mu}}]$,

• •

${\mathcal{W}}^{\intercal}={\mathcal{W}}$,

• •

${\mathcal{W}}[{\boldsymbol{\mu}}-2{\mathcal{D}}[{\boldsymbol{\mu}}]]=0$,

• •

${\mathcal{W}}[{\mathcal{S}}[{\boldsymbol{\mu}}]]=0$,

where the transpose is induced by the bilinear form 17.

Proof.

The first two results follow from the definitions of ${\mathcal{W}}$ and the normal vector. The other two follow from the fact that ${\bf S}[{\boldsymbol{\mu}}]$ and ${\bf D}[{\boldsymbol{\mu}}]$ are divergence-free and an application of lemma 1. ∎

Remark 2**.**

The operators ${\bf S},{\bf D},{\mathcal{D}}$ and ${\mathcal{D}}^{\intercal}$ depend on the Helmholtz parameter $k$ of the oscillatory Stokes equation. In places where it is essential to highlight this dependence, in a slight abuse of notation, we will use the symbols ${\bf S}_{k},{\bf D}_{k},{\mathcal{D}}_{k}$ and ${\mathcal{D}}^{\intercal}_{k}$ to denote this dependence. Similarly, we will use $A^{\Gamma}$ instead of the operator $A$ to highlight the dependence of the operator $A$ on the boundary of the region $\Gamma$.

3 Fredholm analysis of the integral representations

In this section, we establish how layer potential representations of oscillatory Stokes velocity fields can be used to compute Stokes eigenvalues. The main results show that for certain representations the resulting integral equation is not invertible precisely when $k^{2}$ is an eigenvalue. For the interior Dirichlet eigenvalue problem, this is proved separately for a double layer representation on simply connected domains in theorem 8 and for a combined-field representation on multiply connected domains in theorem 10.

Before proving the main theorems, we require a number of uniqueness results for oscillatory Stokes boundary value problems. To prove the uniqueness results, we follow the structure presented in Colton and Kress [15, Ch. 3] for the scalar Helmholtz equation. While uniqueness results for interior Dirichlet, Neumann, and impedance problems follow from energy considerations and compactness arguments, the proofs for the uniqueness of exterior problems are more involved. In particular, the exterior problems are only well-posed after imposing an appropriate radiation condition. We formulate well-posed boundary value problems for both interior and exterior domains with Dirichlet, Neumann and impedance boundary conditions and present uniqueness results for each.

Along with the Fredholm alternative, these uniqueness results are sufficient to prove theorems 8 and 10. Once the main theorems are established, the details of how to use the Fredholm determinant as a numerical tool for computing the Stokes eigenvalues follow in a straightforward manner from the results in [54]. We reproduce these results in the present context for completeness.

3.1 Boundary value problems — interior

Let $\Omega$ be a bounded domain with a $C^{2}$ boundary denoted by $\Gamma$. We summarize the interior Dirichlet, Neumann and impedance boundary value problems in definitions 1, 2 and 3 below.

Definition 1** (Interior Dirichlet problem).**

Let ${\boldsymbol{f}}\in C(\Omega)$ be given. Find $({\boldsymbol{u}},p)\in A(\Omega)$ such that

[TABLE]

Note that the divergence-free constraint for the oscillatory Stokes equations implies a compatibility condition on the Dirichlet data ${\boldsymbol{f}}$, namely that

[TABLE]

Definition 2** (Interior Neumann problem).**

Let ${\boldsymbol{g}}\in C(\Omega)$ be given. Find $({\boldsymbol{u}},p)\in A(\Omega)$ such that

[TABLE]

Definition 3** (Interior impedance problem).**

Let ${\boldsymbol{h}}\in C(\Omega)$ be given and suppose $\eta\in\mathbb{C}$ with $\real{(\eta)}>0$ and ${\textrm{{Im}}}{(\eta)}\geq 0$. Find $({\boldsymbol{u}},p)\in A(\Omega)$ such that

[TABLE]

3.2 A radiation condition for the oscillatory Stokes

equation

Let $\Omega$ be the union of a finite collection of simply connected domains, i.e. $\Omega=\bigcup_{i=1}^{m}\Omega_{i}$ for some $m\in\mathbb{N}$, and let $E=\mathbb{R}^{2}\setminus\bar{\Omega}$ denote its exterior; see fig. 1 for an example with $m=4$. Let $\Gamma=\partial E$ denote the boundary of $E$ and ${\boldsymbol{\nu}}({\boldsymbol{y}})$ denote the exterior normal to the point ${\boldsymbol{y}}$ on $\Gamma$, i.e. the normal vector pointing out of $E$ into $\Omega$. For a given function ${\boldsymbol{f}}$ defined on $\Gamma$, the exterior Dirichlet boundary value problem is to find a pair $({\boldsymbol{u}},p)$ which satisfies:

[TABLE]

In addition to the boundary condition on $\Gamma$, we must impose radiation conditions at $\infty$, analogous to the Helmholtz equation.

Let $B_{r}(0)$ denote the disc of radius $r$ centered at the origin and $\partial B_{r}(0)$ its boundary. We propose the following radiation condition.

Definition 4**.**

Let $({\boldsymbol{u}},p)$ satisfy the oscillatory Stokes equations in the exterior of a bounded domain. We say that the pair $({\boldsymbol{u}},p)$ is radiating if

[TABLE]

uniformly in direction where ${\boldsymbol{t}}={\boldsymbol{\sigma}}\cdot{\boldsymbol{\nu}}$ with ${\boldsymbol{\nu}}={\boldsymbol{x}}/|{\boldsymbol{x}}|$, i.e. ${\boldsymbol{t}}$ is the surface traction on $\partial B_{r}(0)$.

In the following lemma, we show that the oscillatory Stokeslet satisfies the radiation condition.

Lemma 4**.**

The oscillatory Stokeslet, as defined in (10), satisfies the radiation condition in Definition 4.

Proof.

Consider the Stokeslet induced by an arbitrary charge $k^{2}{\boldsymbol{\psi}}$ at the origin where ${\boldsymbol{\psi}}\in\mathbb{C}^{2}$ is a constant. Let $r=|{\boldsymbol{x}}|$, ${\boldsymbol{\nu}}({\boldsymbol{x}})={\boldsymbol{x}}/|{\boldsymbol{x}}|$, and ${\boldsymbol{\tau}}({\boldsymbol{x}})={\boldsymbol{\nu}}({\boldsymbol{x}})^{\perp}$. We have

[TABLE]

Note that derivatives of $\log r$ are $o(1/\sqrt{r})$ and that the pressure associated with the Stokeslet is $p=\nabla G^{\textrm{L}}({\boldsymbol{x}})\cdot{\boldsymbol{\psi}}$. We then have

[TABLE]

Because $H_{0}^{(1)}(kr)$ has the asymptotic expansion

[TABLE]

as $r\to\infty$, we have

[TABLE]

Finally, since $H_{0}^{(1)}(kr)$ is radially symmetric, we have

[TABLE]

so that the Stokeslet satisfies the radiation condition. ∎

A consequence of the above lemma is that the oscillatory Stokes single layer potential satisfies the radiation condition.

Corollary 2**.**

Suppose that $\Gamma$ is the boundary of a region $\Omega$ and is $C^{2}$. Suppose that ${\boldsymbol{\mu}}\in C(\Gamma)\times C(\Gamma)$, then the oscillatory Stokes single layer potential ${\bf S}[{\boldsymbol{\mu}}]$, as defined in 12, satisfies the radiation condition.

Unfortunately, the stresslet, as defined in 11, does not necessarily satisfy the radiation condition. The reason for failure is the logarithmic growth of the pressure at $\infty$. However, this turns out to be a rank-one issue and the oscillatory Stokes double layer potential does satisfy the radiation condition if the density satisfies an integral constraint. The following lemma proves this result.

Lemma 5**.**

Suppose that $\Gamma$ is the boundary of a region $\Omega$ and is $C^{2}$. Suppose that ${\boldsymbol{\mu}}\in C(\Gamma)\times C(\Gamma)$ and satisfies $\int_{\Gamma}{\boldsymbol{\mu}}\cdot{\boldsymbol{\nu}}dS=0$, where ${\boldsymbol{\nu}}$ denotes the outward normal to the curve $\Gamma$. Then, the oscillatory Stokes double layer potential ${\bf D}[{\boldsymbol{\mu}}]$, as defined in (13), also satisfies the radiation condition.

Proof.

We only establish the decay of the pressure, which we will denote by $p^{\bf D}$; the rest of the terms in (29) can be bounded using an argument like that for the Stokeslet above. Because $p^{\bf D}$ is harmonic in the exterior of any disc containing $\Gamma$, it is sufficient to show that $|\nabla p^{\bf D}|=\mathcal{O}(1/r^{2})$. Let ${\boldsymbol{\tau}}$ denote the positively oriented tangent to the curve $\Gamma$ and $\mu_{\nu}({\boldsymbol{y}})$ and $\mu_{\tau}({\boldsymbol{y}})$ denote ${\boldsymbol{\mu}}({\boldsymbol{y}})\cdot{\boldsymbol{\nu}}({\boldsymbol{y}})$ and ${\boldsymbol{\mu}}({\boldsymbol{y}})\cdot{\boldsymbol{\tau}}({\boldsymbol{y}})$, respectively. Substituting ${\bf D}{\boldsymbol{\mu}}$ into (1), we obtain

[TABLE]

The other terms are higher-order derivatives of $G^{\textrm{L}}$, so it is sufficient to show that the term

[TABLE]

is $\mathcal{O}(1/r^{2})$. In the following, let $z=x_{1}+ix_{2}$ be the point corresponding to ${\boldsymbol{x}}$ in the complex plane and let $R$ be the radius of some disc containing $\Gamma$. If $|{\boldsymbol{x}}|>2R$, we can use the standard multipole expansion of $\log(z-(y_{1}+iy_{2}))$ and the assumption that $\int_{\Gamma}{\boldsymbol{\mu}}\cdot{\boldsymbol{\nu}}=0$ to obtain

[TABLE]

∎

3.3 Boundary value problems — exterior

Let $E$ and $\Gamma$ be as in the previous subsection. The radiation condition allows for a well-posed formulation of the exterior boundary value problems, which we summarize in definitions 5, 6 and 7 below.

Definition 5** (Exterior Dirichlet problem).**

Let ${\boldsymbol{f}}\in C(\Gamma)$ be given. Find $({\boldsymbol{u}},p)\in A(E)$ such that

[TABLE]

and $({\boldsymbol{u}},p)$ satisfies the radiation condition in definition 4.

Definition 6** (Exterior Neumann problem).**

Let ${\boldsymbol{g}}\in C(\Gamma)$ be given. Find $({\boldsymbol{u}},p)\in A(E)$ such that

[TABLE]

and $({\boldsymbol{u}},p)$ satisfies the radiation condition in definition 4.

Definition 7** (Exterior impedance problem).**

Let ${\boldsymbol{h}}\in C(\Gamma)$ be given and suppose that $\eta\in\mathbb{C}$ with $\real{(\eta)}>0$ and ${\textrm{{Im}}}{(\eta)}\geq 0$. Find $({\boldsymbol{u}},p)\in A(E)$ such that

[TABLE]

and $({\boldsymbol{u}},p)$ satisfies the radiation condition in definition 4.

3.4 Uniqueness results

Before moving on to the exterior uniqueness theorems, we establish the well-known result that the $k$ corresponding to interior eigenvalues, $k^{2}$, are real-valued.

Theorem 3**.**

Let $\Omega$ be a bounded domain and suppose that ${\textrm{{Im}}}(k)\neq 0$. Then both the interior Dirichlet and Neumann boundary value problems have unique solutions.

Proof.

A couple applications of the divergence theorem establish that

[TABLE]

Suppose that either ${\boldsymbol{u}}=0$ or ${\boldsymbol{t}}=0$ on $\Gamma$. Then, the right hand side of 33 is zero. Taking the real and imaginary parts of 33, it is clear that ${\boldsymbol{u}}\equiv 0$, if $\text{Im}(k)\neq 0$. ∎

In the following lemma, we prove uniqueness for the interior impedance problem.

Theorem 4**.**

Let $\Omega$ be a bounded domain and suppose that $\real(k),{\textrm{{Im}}}(k)>0$. Then the interior impedance problem has a unique solution.

Proof.

Plugging ${\boldsymbol{t}}=i\eta{\boldsymbol{u}}$ in 33 and taking the imaginary part, we get

[TABLE]

from which it is is clear that ${\boldsymbol{u}}\equiv 0$, since $\real{(\eta)},\real{(k)},{\textrm{{Im}}}{(k)}>0$. ∎

We now turn our attention to the proofs of the exterior boundary value problems for the oscillatory Stokes equation. The following lemmas are useful for proving these results.

Lemma 6**.**

Let the unbounded region $E$ be given as the exterior of a finite collection of bounded domains. Suppose that $({\boldsymbol{u}},p)$ satisfies the oscillatory Stokes equation in $E$ as well as the radiation condition 29. Then

[TABLE]

Proof.

Since $({\boldsymbol{u}},p)$ satisfies the radiation condition, we have that

[TABLE]

Since ${\boldsymbol{u}}$ satisfies the oscillatory Stokes equation $E\cap B_{r}(0)$, using a couple of applications of the divergence theorem, we have that

[TABLE]

Combining 36 and 37, we get

[TABLE]

∎

In the next lemma, we prove the analogue of Rellich’s lemma for the oscillatory Stokes equation.

Lemma 7**.**

Let the unbounded region $E$ be given as the exterior of a finite collection of bounded domains. Suppose $k$ is real, ${\boldsymbol{u}}$ satisfies the oscillatory Stokes equation in $E$, and that

[TABLE]

Then each component of ${\boldsymbol{u}}$ is harmonic in $E$.

Proof.

We first note that each component of ${\boldsymbol{u}}=(u_{1},u_{2})$ satisfies the oscillatory biharmonic equation in $E$, i.e.

[TABLE]

For $r$ sufficiently large, we can express $u_{j}$ in the Fourier basis as

[TABLE]

Using Parseval’s identity then

[TABLE]

Since ${\boldsymbol{u}}$ satisfies 38, we conclude that

[TABLE]

Since $u_{j}$, $j=1,2$ satisfies the oscillatory biharmonic equation, the functions $a_{j,n}$ are linear combinations of

[TABLE]

and

[TABLE]

where $H_{n}^{(1),(2)}(\cdot)$ are the Hankel functions of the first and second kind of order $n$. Since $a_{j,n}(r)$ satisfy 39, and using the asymptotic expansion of $H_{n}^{(1),(2)}(kr)$ when $k$ and $r$ are real-valued, we note that the projection of $a_{j,n}$ on $r^{|n|}$, and $H_{n}^{1,2}(kr)$ must be zero. Thus, for sufficiently large $r$,

[TABLE]

i.e. $u_{j}$ is harmonic in $B_{r}(0)^{c}$. Finally, by corollary 1, ${\boldsymbol{u}}$ is analytic in $E$. Therefore, each $u_{j}$ is harmonic throughout $E$. ∎

Remark 3**.**

Note that if ${\boldsymbol{u}}$ satisfies the assumptions of lemma 7, then each component is harmonic and thus ${\boldsymbol{u}}$ satisfies

[TABLE]

Remark 4**.**

It should be noted that in lemma 7, ${\boldsymbol{u}}$ need not be a radiating solution. All that is assumed of ${\boldsymbol{u}}$ is that it satisfies the oscillatory Stokes equations in $E$.

We now have the results needed to establish the uniqueness of exterior boundary value problems.

Theorem 5** (Uniqueness of the Exterior Dirichlet Problem).**

Let the unbounded region $E$ be given as the exterior of a finite collection of bounded domains. Suppose that ${\textrm{{Im}}}(k)\geq 0$ and that $({\boldsymbol{u}},p)$ is a radiating solution to the oscillatory Stokes equation in $E$ with ${\boldsymbol{u}}=0$ on the boundary $\Gamma$, then ${\boldsymbol{u}}\equiv 0$ in $E$.

Proof.

Since ${\boldsymbol{u}}=0$ on $\Gamma$, it follows from lemma 6 that

[TABLE]

Suppose that ${\textrm{{Im}}}(k)>0$. Then, it is immediate that ${\boldsymbol{u}}\equiv 0$ in $E$.

Suppose that $k$ is real valued. It is clear that

[TABLE]

Thus, the conditions on ${\boldsymbol{u}}$ and $k$ in lemma 7 are satisfied, and each component of ${\boldsymbol{u}}$ is a harmonic function with ${\boldsymbol{u}}\to 0$ as $r\to\infty$. Furthermore, since ${\boldsymbol{u}}=0$ on $\Gamma$, by the uniqueness of solutions to the Dirichlet problem for Laplace’s equation on exterior domains, we conclude that ${\boldsymbol{u}}\equiv 0$ in $E$. ∎

Theorem 6** (Uniqueness of the Exterior Neumann Problem).**

Suppose that $\Omega$ is the union of a finite collection of simply connected domains, i.e. $\Omega=\bigcup_{i=1}^{m}\Omega_{i}$ for some $m\in\mathbb{N}$, with $C^{2}$ boundaries, and let $E=\mathbb{R}^{2}\setminus\bar{\Omega}$ denote its exterior; see fig. 1 for an example with $m=4$. Let $\Gamma_{i}$ denote the boundary of $\Omega_{i}$, and $\Gamma=\bigcup_{i=1}^{m}\Gamma_{i}$ denote the boundary of $\Omega$. Suppose that ${\textrm{{Im}}}(k)\geq 0$ and that $({\boldsymbol{u}},p)$ is a radiating solution to the oscillatory Stokes equation in $E$ with ${\boldsymbol{t}}=0$ on the boundary $\Gamma$, then ${\boldsymbol{u}}\equiv 0$ in $E$.

Proof.

Since ${\boldsymbol{t}}=0$ on $\Gamma$, it follows from 35 that

[TABLE]

Suppose that ${\textrm{{Im}}}(k)>0$. It is then immediate that ${\boldsymbol{u}}\equiv 0$ in $E$.

Suppose that $k$ is real. It is clear that

[TABLE]

Thus, the conditions on ${\boldsymbol{u}}$ and $k$ in lemma 7 are satisfied, and each component of ${\boldsymbol{u}}$ is a harmonic function with ${\boldsymbol{u}}\to 0$ as $r\to\infty$. Furthermore, as observed in remark 3, $k^{2}{\boldsymbol{u}}=\nabla p$. Then, the boundary condition becomes $0={\boldsymbol{t}}=-p{\boldsymbol{\nu}}+2\nabla\partial_{\nu}p/k^{2}$. Because $0={\boldsymbol{\tau}}\cdot{\boldsymbol{t}}=2\partial_{\tau\nu}p/k^{2}$, $\partial_{\nu}p=c_{i}$ on $\Gamma_{i}$ for each $\Gamma_{i}$, where $c_{i}$ is a constant. Observe that $|{\boldsymbol{u}}|$ and $|{\boldsymbol{e}}({\boldsymbol{u}})|$ must be $O(1/r)$ as $r\to\infty$. Thus, for a radiating pair $({\boldsymbol{u}},p)$ with $p$ harmonic, we have that $|p|=O(1/r)$ and $|\nabla p|=O(1/r^{2})$ as $r\to\infty$. Since the boundary is $C^{2}$ and the boundary data for $p$ is analytic, we conclude that $p$ is $C^{2}$ in $\overline{E}$. Furthermore, ${\boldsymbol{t}}=0$ implies $p{\boldsymbol{\nu}}=2\nabla\partial_{\nu}p/k^{2}$, and taking the dot product with ${\boldsymbol{\nu}}$, we get $p=2\partial_{\nu\nu}p/k^{2}$. It then follows that $p=-2\partial_{\tau\tau}p/k^{2}$ on $\Gamma$ since $p$ is harmonic in $E$ and $C^{2}$ in $\overline{E}$. Since $p$ satisfies the radiation condition at $\infty$, we have

[TABLE]

Thus, $p$ is a constant in $E$. Furthermore, since $p\to 0$ at $\infty$, we conclude that $p\equiv 0$ in $E$. Finally, since $k^{2}{\boldsymbol{u}}=\nabla p$, we conclude that ${\boldsymbol{u}}\equiv 0$ in $E$. ∎

Theorem 7** (Uniqueness of the Exterior Impedance Problem).**

Let the unbounded region $E$ be given as the exterior of a finite collection of bounded domains. Suppose that the complex numbers $\eta$ and $k$ satisfy that $\real(\eta),\real(k)>0$ and ${\textrm{{Im}}}(\eta),{\textrm{{Im}}}(k)\geq 0$. Suppose further that $({\boldsymbol{u}},p)$ is a radiating solution of the oscillatory Stokes equation in $E$ which satisfies the homogeneous impedance boundary condition

[TABLE]

Then ${\boldsymbol{u}}\equiv 0$ for ${\boldsymbol{x}}\in E$.

Proof.

Since ${\boldsymbol{u}}$ satisfies the radiation condition at $\infty$ and ${\boldsymbol{t}}=-i\eta{\boldsymbol{u}}$ on $\Gamma$, it follows from 35 that

[TABLE]

Because all of the quantities in the last expression above are nonnegative, we have that

[TABLE]

The result then follows from the uniqueness of solutions to the exterior Dirichlet problem. ∎

3.5 The integral equations and their null-spaces

In this section, we establish the correspondence between Stokes eigenvalues and the invertibility of certain integral equations arising from layer potential representations of solutions to the oscillatory Stokes equation.

Let $\Omega$ be a bounded domain given as the intersection of a simply connected domain $\Omega_{0}$ and the exteriors of a finite collection of bounded, simply connected domains $\{\Omega_{i}\}_{i=1}^{m}$ whose closures are contained in $\Omega_{0}$; see fig. 2 for an example with four inclusions. Note that the exterior of $\Omega$, which we denote by $E$, is the disjoint union of the exterior of $\Omega_{0}$, which we denote $E_{0}$, with the sets $\{\Omega_{i}\}_{i=1}^{m}$. Let $\Gamma$ denote the boundary of $\Omega$ with the normal ${\boldsymbol{\nu}}$ pointing out of $\Omega$. We will use superscript $+$ and $-$ signs to indicate the limit values of a function on $\Gamma$ as approached from the exterior and interior, respectively.

For the sake of brevity, we consider only the Dirichlet eigenvalue problem but a similar analysis could be applied to the Neumann eigenvalue problem. We analyze two different representations for the Dirichlet problem: a double layer potential, i.e. setting ${\boldsymbol{u}}={\bf D}_{k}{\boldsymbol{\mu}}$, and a combined-field layer potential, i.e. setting ${\boldsymbol{u}}=(i\eta{\bf S}_{k}+{\bf D}_{k}){\boldsymbol{\mu}}$. While both representations result in a second kind integral equation for the oscillatory Stokes Dirichlet problem, the double layer potential has spurious non-trivial nullspaces on domains with positive genus, as explained below.

Remark 5**.**

The application of the Fredholm alternative here again follows the structure used for the Laplace eigenvalue problem in [15, Ch. 3].

3.5.1 Dirichlet eigenvalues — double layer

representation

Suppose that the solution to the oscillatory Stokes Dirichlet problem, 25, is represented using a double layer potential defined on $\Gamma$, i.e. setting ${\boldsymbol{u}}={\bf D}_{k}{\boldsymbol{\mu}}$ where ${\boldsymbol{\mu}}$ is an unknown density. Substituting this expression into the boundary condition and applying lemma 1, we obtain

[TABLE]

The rank deficiency of $\mathcal{I}-2{\mathcal{D}}_{k}$ is well-known [7] and we summarize it in the lemma below.

Lemma 8**.**

In the notation above, ${\boldsymbol{\nu}}\in{\mathcal{N}}(\mathcal{I}-2{\mathcal{D}}_{k}^{\intercal})$.

Proof.

From lemma 3, we note that $W[(\mathcal{I}-2{\mathcal{D}}_{k}){\boldsymbol{\mu}}]=0$ implies that $\left<(\mathcal{I}-2{\mathcal{D}}_{k}){\boldsymbol{\mu}},{\boldsymbol{\nu}}\right>=0$ for all ${\boldsymbol{\mu}}$, i.e. ${\boldsymbol{\nu}}\in R(\mathcal{I}-2{\mathcal{D}}_{k})^{\perp}$, where $R(A)$ denotes the range of the operator $A$. By the Fredholm alternative, the result then follows. ∎

Thus, we instead analyze the equation

[TABLE]

Note that if ${\boldsymbol{f}}$ satisfies the compatibility condition $\int_{\Gamma}{\boldsymbol{f}}\cdot{\boldsymbol{\nu}}\,dS=0$, then 43 implies 42.

On simply connected domains, there is a one-to-one correspondence between the eigenvalues of the Dirichlet problem for the Stokes equation and the values of $k$ for which the operator $(\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}})$ is not invertible. To prove this, we also need the following lemma:

Lemma 9**.**

If ${\boldsymbol{t}}^{-}$ is the surface traction associated with an interior Dirichlet Stokes eigenfunction ${\boldsymbol{u}}$, then ${\boldsymbol{t}}^{-}$ and ${\boldsymbol{\nu}}$ are linearly independent.

Proof.

We first note that ${\bf S}_{k}[{\boldsymbol{\nu}}]({\boldsymbol{x}})=0$ for all ${\boldsymbol{x}}\in\Omega$. This follows from an application of the divergence theorem and the fact that oscillatory Stokeslet is divergence free in $\Omega$. If ${\boldsymbol{t}}^{-}$ is a surface traction associated with a Stokes eigenvalue, then using Green’s theorem, it follows that ${\bf S}_{k}[{\boldsymbol{t}}^{-}]({\boldsymbol{x}})={\boldsymbol{u}}({\boldsymbol{x}})\neq 0$ for ${\boldsymbol{x}}\in\Omega$ and thus ${\boldsymbol{t}}^{-}$ and ${\boldsymbol{\nu}}$ are linearly independent. ∎

Theorem 8**.**

Suppose that $\Omega$ is a bounded, simply connected domain. Then, the operator $(\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}})$ is not invertible if and only if $k^{2}$ is a Dirichlet eigenvalue for the Stokes equation on $\Omega$.

Proof.

Suppose that $k^{2}$ is not a Dirichlet eigenvalue for the Stokes equation on $\Omega$. Suppose further that ${\boldsymbol{\mu}}$ satisfies

[TABLE]

i.e. ${\boldsymbol{\mu}}$ is in the null-space of $(\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}})$. Applying the operator ${\mathcal{W}}$ to 44 and using lemma 3, we get

[TABLE]

Thus 44 reduces to

[TABLE]

Suppose now ${\boldsymbol{u}}=-2{\bf D}_{k}[{\boldsymbol{\mu}}]$ in $\Omega$. Then ${\boldsymbol{u}}$ is a solution to the oscillatory Stokes equation in $\Omega$, and applying lemma 1, we get that the interior limit of the velocity ${\boldsymbol{u}}^{-}=(\mathcal{I}-2{\mathcal{D}}_{k}){\boldsymbol{\mu}}=0$ on $\Gamma$. Since $k^{2}$ is not a Dirichlet eigenvalue for the Stokes equation on $\Omega$, we conclude that ${\boldsymbol{u}}\equiv 0$ in $\Omega$. In particular, this implies that the interior limit of the surface traction denoted by ${\boldsymbol{t}}^{-}=0$ on $\Gamma$. Using lemma 1 again, we conclude that the exterior limit of the surface traction, ${\boldsymbol{t}}^{+}$, is [math] on $\Gamma$. Note that ${\boldsymbol{u}}$ is a radiating solution of the oscillatory Stokes equation in the exterior $E$, as ${\mathcal{W}}[{\boldsymbol{\mu}}]=0$ implies $\int_{\Gamma}{\boldsymbol{\mu}}\cdot{\boldsymbol{\nu}}=0$. From the uniqueness of solutions to the exterior Neumann problem, we conclude that ${\boldsymbol{u}}\equiv 0$ in $E$ as well, which in particular implies that the exterior limit of the velocity, ${\boldsymbol{u}}^{+}$, is [math] on $\Gamma$. Using the jump conditions in lemma 1 again, we get that $2{\boldsymbol{\mu}}={\boldsymbol{u}}^{-}-{\boldsymbol{u}}^{+}=0$. Thus $(\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}})$ is invertible if $k^{2}$ is not a Dirichlet eigenvalue for the Stokes equation on $\Omega$.

To prove the converse, note that from theorem 2 we have

[TABLE]

Suppose that $k^{2}$ is a Dirichlet eigenvalue for the Stokes equation on $\Omega$ and let ${\boldsymbol{u}}$ denote the corresponding eigenfunction with ${\boldsymbol{t}}^{-}$ the corresponding surface traction on the boundary $\Gamma$. Since ${\boldsymbol{u}}$ is a Dirichlet eigenfunction, the velocity restricted to the boundary, ${\boldsymbol{u}}^{-}$, is [math]. Using the Green’s theorem representation for the pair $({\boldsymbol{u}},{\boldsymbol{t}}^{-})$ and evaluating the surface traction using lemma 1, we get

[TABLE]

From the above and lemmas 8 and 9, we know that ${\boldsymbol{t}}^{-},{\boldsymbol{\nu}}\in{\mathcal{N}}(\mathcal{I}-2{\mathcal{D}}^{\intercal}_{k})$ are two linearly independent vectors in the null space. Let $c=\left<{\boldsymbol{t}}^{-},{\boldsymbol{\nu}}\right>$. Then, it follows that $\left<{\boldsymbol{t}}^{-}-c{\boldsymbol{\nu}},{\boldsymbol{\nu}}\right>=0$ and thus ${\mathcal{W}}[{\boldsymbol{t}}^{-}-c{\boldsymbol{\nu}}]=0$. Since ${\boldsymbol{t}}^{-}$ and ${\boldsymbol{\nu}}$ are linearly independent, we note that ${\boldsymbol{t}}^{-}-c{\boldsymbol{\nu}}\neq 0$. Combining these results, we get that $(\mathcal{I}-2{\mathcal{D}}^{\intercal}_{k}-2{\mathcal{W}})({\boldsymbol{t}}^{-}-c{\boldsymbol{\nu}})=0$. Since ${\boldsymbol{t}}^{-}-c{\boldsymbol{\nu}}$ is non-trivial, and ${\mathcal{W}}$ is self-adjoint, it follows from the Fredholm alternative that the operator $\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ is also not invertible. ∎

The correspondence result above does not hold on multiply connected domains. In particular, while the operator $\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ is indeed not invertible when $k^{2}$ is a Dirichlet eigenvalue, it turns out that the operator is also not invertible when $k^{2}$ is a Neumann eigenvalue corresponding to the interior of one of the obstacle regions, i.e. one of the $\Omega_{i}$ with $i>0$. The following theorem proves this result for a region with one obstacle; the extension to the general case is straightforward.

Theorem 9**.**

Suppose that $\Omega$ is a multiply connected domain given by the intersection of a domain $\Omega_{0}$ and the exterior of a single domain $\Omega_{1}$ with $\bar{\Omega}_{1}\subset\Omega_{0}$. Then, the operator $\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ is not invertible if $k^{2}$ is a Neumann eigenvalue of $\Omega_{1}$.

Proof.

Suppose that $k^{2}$ is a Neumann eigenvalue of $\Omega_{1}$, and let $\tilde{{\boldsymbol{u}}}$ denote the corresponding eigenfunction. Note that $\tilde{{\boldsymbol{u}}}$ is not identically [math] on the boundary of $\Omega_{1}$, which we denote by $\Gamma_{1}$. Since $\tilde{{\boldsymbol{u}}}$ is an interior Neumann eigenfunction, we note that the surface traction corresponding to the solution, ${\boldsymbol{t}}^{-}$, is [math] on the boundary. Applying theorem 2 to the solution $\tilde{{\boldsymbol{u}}}$ in the interior and taking the interior limit we get

[TABLE]

Note that the sign of the ${\bf D}$ operator in the representation theorem is switched since the normal is pointing inwards for the boundary $\Omega_{1}$. Thus, $\tilde{{\boldsymbol{u}}}$ is a non-trivial null vector of the operator $\frac{1}{2}\mathcal{I}-{\mathcal{D}}^{\Gamma_{1}}_{k}$. Furthermore since $\tilde{{\boldsymbol{u}}}$ is the boundary data of the solution of the oscillatory Stokes equation in $\Omega_{1}$, we get that ${\mathcal{W}}^{\Gamma_{1}}[\tilde{{\boldsymbol{u}}}]=0$. Setting ${\boldsymbol{\mu}}=\tilde{{\boldsymbol{u}}}$ on $\Gamma_{1}$, and ${\boldsymbol{\mu}}=0$ on $\Gamma_{0}$, we obtain a non-trivial null vector for the operator $\mathcal{I}-2{\mathcal{D}}^{\Gamma}_{k}-2{\mathcal{W}}^{\Gamma}$ on the boundary $\Gamma=\Gamma_{0}\cup\Gamma_{1}$. ∎

The spurious eigenvalues of the operator $\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ are demonstrated in section 4.2.2 on an annulus, where both the true and spurious eigenvalues can be determined analytically. Analogous with the observation in [54], this lack of one-to-one correspondence between the invertibility of the integral operator $\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ and the Dirichlet eigenvalues of the Stokes operator on multiply connected domains also causes non-robustness and introduces near-resonances for simply-connected domains which are almost multiply-connected.

3.5.2 Dirichlet eigenvalues — combined-field representation

The non-robustness in using the double layer potential representation can be remedied by using a combined-field, or mixed layer potential, representation, i.e. setting ${\boldsymbol{u}}=({\bf D}_{k}+i\eta{\bf S}_{k}){\boldsymbol{\mu}}$, with $\eta$ real and positive. Imposing the Dirichlet boundary condition and using lemma 1, we obtain

[TABLE]

As with the double layer representation, this integral equation is rank deficient for any $k$. Instead, we consider

[TABLE]

We now prove that for any bounded region $\Omega$ (simply or multiply connected) with $C^{2}$ boundaries, there exists a one-to-one correspondence between the invertibility of the operator $\mathcal{I}-2{\mathcal{D}}_{k}-2i\eta{\mathcal{S}}_{k}-2{\mathcal{W}}$ and the Dirichlet eigenvalues.

Theorem 10**.**

Suppose $\Omega$ is a bounded region defined by the intersection of a simply connected domain $\Omega_{0}$ and the exteriors of a finite collection bounded simply connected domains $\{\Omega_{i}\}_{i=1}^{m}$. As above, let $\Gamma_{i}$ denote the boundary of $\Omega_{i}$ and let $\Gamma=\cup_{i=0}^{m}\Gamma_{i}$ denote the boundary of $\Omega$. Then, the operator $\mathcal{I}-2{\mathcal{D}}_{k}-2i\eta{\mathcal{S}}_{k}-2{\mathcal{W}}$ is invertible if and only if $k^{2}$ is not a Dirichlet eigenvalue for the Stokes operator on $\Omega$.

Proof.

Suppose that $k^{2}$ is not a Dirichlet eigenvalue for the Stokes equation on $\Omega$. Suppose further that ${\boldsymbol{\mu}}$ satisfies

[TABLE]

i.e. ${\boldsymbol{\mu}}$ is in the null-space of $(\mathcal{I}-2{\mathcal{D}}_{k}-2i\eta{\mathcal{S}}_{k}-2{\mathcal{W}})$. Applying the operator ${\mathcal{W}}$ to 52 and using lemma 3, we get

[TABLE]

Thus 52 reduces to

[TABLE]

Suppose now ${\boldsymbol{u}}=-2{\bf D}_{k}[{\boldsymbol{\mu}}]-2i\eta{\bf S}_{k}[{\boldsymbol{\mu}}]$ in $\Omega$. Then ${\boldsymbol{u}}$ is a solution to the oscillatory Stokes equation in $\Omega$, and applying lemma 1, we get that the interior limit of the velocity ${\boldsymbol{u}}^{-}=(\mathcal{I}-2{\mathcal{D}}_{k}-2i\eta{\mathcal{S}}_{k}){\boldsymbol{\mu}}=0$ on $\Gamma$. Since $k^{2}$ is not a Dirichlet eigenvalue for the Stokes equation on $\Omega$, we conclude that ${\boldsymbol{u}}\equiv 0$ in $\Omega$. This in particular implies that the interior limit of the surface traction, denoted by ${\boldsymbol{t}}^{-}$, is [math] on $\Gamma$.

Using lemma 1 we observe that ${\boldsymbol{t}}^{+}=2i\eta{\boldsymbol{\mu}}({\boldsymbol{x}})$ and ${\boldsymbol{u}}^{+}=-2{\boldsymbol{\mu}}({\boldsymbol{x}})$ on $\Gamma$, i.e. ${\boldsymbol{t}}^{+}+i\eta{\boldsymbol{u}}^{+}=0$ and ${\boldsymbol{u}}^{+}$ satisfies the homogeneous exterior impedance problem. We first show that ${\boldsymbol{\mu}}=0$ on $\Gamma_{0}$. To this end, note that ${\boldsymbol{u}}$ is a radiating solution of the oscillatory Stokes equation in the exterior $E$, since ${\mathcal{W}}[{\boldsymbol{\mu}}]=0$ implies that $\int_{\Gamma}{\boldsymbol{\mu}}\cdot{\boldsymbol{\nu}}=0$. From the uniqueness of the impedance problem in the exterior $E_{0}$ of $\Omega_{0}$, we conclude that ${\boldsymbol{u}}\equiv 0$ in $E_{0}$ as well, which in particular implies that ${\boldsymbol{u}}^{+}=0$ on $\Gamma_{0}$. Using the jump conditions in lemma 1 again, we get that $2{\boldsymbol{\mu}}={\boldsymbol{u}}^{-}-{\boldsymbol{u}}^{+}=0$ on $\Gamma_{0}$.

Remark 6**.**

Note that there is potential for confusion here in that the exterior limit with respect to $\Omega$ for the boundary $\Gamma_{j}$ is the traditional interior limit with respect to the obstacle region $\Omega_{j}$.

To show that ${\boldsymbol{\mu}}=0$ on $\Gamma_{j}$, we observe that ${\boldsymbol{u}}$ is also a solution to the oscillatory Stokes equation in each of the obstacles $\Omega_{j}$. Using the jump conditions in lemma 1, we get that ${\boldsymbol{t}}^{+}=2i\eta{\boldsymbol{\mu}}$ and ${\boldsymbol{u}}^{+}=-2{\boldsymbol{\mu}}$. However, the normal is inward pointing inside $\Omega_{j}$ on the boundary $\Gamma_{j}$. If we revert back to the normal being defined as an outward normal to $\Omega_{j}$, then the boundary conditions on $\Gamma_{j}$ is ${\boldsymbol{t}}-i\eta{\boldsymbol{u}}=0$. From the uniqueness of solutions to the interior impedance problem, we conclude that ${\boldsymbol{u}}\equiv 0$ in $\Omega_{j}$, which in particular implies that $2{\boldsymbol{\mu}}={\boldsymbol{u}}^{-}-{\boldsymbol{u}}^{+}=0$ for ${\boldsymbol{x}}\in\Gamma_{j}$, $j=1,2,\ldots m$. Thus, $\mathcal{I}-2{\mathcal{D}}_{k}-2i\eta{\mathcal{S}}_{k}-2{\mathcal{W}}$ is invertible when $k^{2}$ is not a Dirichlet eigenvalue for the Stokes equation on $\Omega$.

From theorem 2, we have

[TABLE]

Suppose that $k^{2}$ is Dirichlet eigenvalue for the Stokes equation on $\Omega$ and let ${\boldsymbol{u}}$ denote the corresponding eigenfunction and ${\boldsymbol{t}}$ denote its surface traction. Note that theorem 2 implies that ${\mathcal{S}}_{k}[{\boldsymbol{t}}^{-}]=0$, since ${\boldsymbol{u}}^{-}=0$ on $\Gamma$. Applying theorem 2 to the pair ${\boldsymbol{t}}^{-},{\boldsymbol{u}}^{-}$ and evaluating the traction on $\Gamma$ using lemma 1, we get

[TABLE]

Combining these two identities, we get that

[TABLE]

As in the proof of theorem 8, letting $c=\left<{\boldsymbol{t}}^{-},{\boldsymbol{\nu}}\right>$, it follows that

[TABLE]

where ${\boldsymbol{t}}^{-}-c{\boldsymbol{\nu}}\neq 0$. Since ${\boldsymbol{t}}^{-}-c{\boldsymbol{\nu}}$ is non-trivial and both $i\eta{\mathcal{S}}_{k}$ and ${\mathcal{W}}$ are self-adjoint with respect to the bilinear form 17, it follows from the Fredholm alternative that the operator $\mathcal{I}-2{\mathcal{D}}_{k}-2i\eta{\mathcal{S}}_{k}-2{\mathcal{W}}$ is also not invertible. ∎

3.6 Fredholm determinants

In this section, we show how the Fredholm determinant can be used as a computational tool for detecting the non-invertibility of $\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}}$. The arguments here follow the structure of the analogous arguments in [54] for Laplace eigenvalues.

Let ${\mathcal{J}}_{1}(X)$ denote the space of trace class operators on $X$, where $X$ is a Hilbert space, which is a subspace of the space of compact operators on $X$. A compact operator ${\mathcal{A}}$ with eigenvalues $\lambda_{i},i\in\mathbb{N}$ is in ${\mathcal{J}}_{1}(X)$ if $\sum_{i}|\lambda_{i}|<\infty$. If ${\mathcal{A}}$, is a trace class operator, then the Fredholm determinant of the operator $\mathcal{I}+{\mathcal{A}}$ is defined by

[TABLE]

So far, we have discussed the Fredholm theory in the space $C(\Gamma)\times C(\Gamma)$ equipped with the bilinear form 17. However, it is more convenient to discuss the theory of Fredholm determinants on Hilbert spaces. We note that both the operators ${\mathcal{D}}_{k}$ and ${\mathcal{S}}_{k}$ are also compact operators mapping $Y\to Y$ where $Y=\mathbb{L}^{2}(\Gamma)\times\mathbb{L}^{2}(\Gamma)$. Furthermore, it is well-known that the spectrum of compact operators with weakly singular kernels coincide on $C(\Gamma)\times C(\Gamma)$ and $Y$ (see [33], for example). So for the rest of the section, we present the discussion of the relevant operators on $Y$ instead of $C(\Gamma)\times C(\Gamma)$.

The operator $-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ is trace class:

Lemma 10**.**

Suppose that $\Gamma$ is a $C^{2}$ curve. Then $-2{\mathcal{D}}_{k}-2{\mathcal{W}}\in{\mathcal{J}}_{1}(Y)$ for all $k\in\mathbb{C}\setminus\{0\}$

Proof.

Using Bessel function asymptotics, we note that the kernel of ${\mathcal{D}}_{k}$ given by ${\bf T}_{\cdot,\cdot,\ell}\nu_{\ell}({\boldsymbol{x}},{\boldsymbol{y}})$ has a leading order singularity of $|{\boldsymbol{x}}-{\boldsymbol{y}}|^{2}\log{|{\boldsymbol{x}}-{\boldsymbol{y}}|^{2}}$ as ${\boldsymbol{x}}\to{\boldsymbol{y}}$ for all $k\in\mathbb{C}\setminus\{0\}$. It follows from the criteria listed in [9, Sec. 2] that ${\mathcal{D}}_{k}$ is a trace-class operator. Since ${\mathcal{W}}$ is a rank-one perturbation independent of $k$, and trace-class operators are a vector space, we conclude that $-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ is also a trace-class operator. ∎

Let $f(k)=\text{det}(\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}})$. First, note that $f(k)$ is an analytic function of $k$ for $k\in\mathbb{C}\setminus\{0\}$, since the kernel of ${\mathcal{D}}_{k}$ is an analytic function of $k$ on that domain, and the Fredholm determinant of an analytic operator is analytic on the domain of analyticity of the operator (see [54], for example).

The zeros of the Fredholm determinant indicate when the operator is not invertible. The following lemma summarizes this result.

Lemma 11**.**

With $f(k)$ defined as above, $f(k)=0$ if and only if $\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ is not invertible.

Proof.

The proof is standard; see, for example, [50, p. 34]. ∎

When $\Omega$ is simply connected, lemma 11 and theorem 8 together imply that $f(k)=0$ if and only if $k^{2}$ is a Dirichlet eigenvalue of the Stokes equation. This reduces the problem of finding eigenvalues to finding the roots of an analytic function.

We now show how this fact can be used to numerically estimate the Dirichlet eigenvalues. Suppose that $D_{k}^{N}$ is a Nyström discretization of the operator $-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ when the boundary $\Gamma$ is discretized with $N$ points. Let $f^{N}(k)=\text{det}(I+D_{k}^{N})$ where here det is the standard matrix determinant. Note that the discretized matrix also depends on the choice of quadrature rule used in the Nyström discretization of the operator.

In [54], the authors prove that for computing the Laplace eigenvalues on regions with analytic boundaries, when the integral operators are discretized using Kress quadrature — a spectrally accurate quadrature rule for such kernels, see [32] — the determinant of the Nyström discretized operators at the true eigenvalues converge to [math] exponentially in $N$. Thus, if the eigenvalues have multiplicity $1$, i.e. the derivative of the determinant is non-zero at the true-eigenvalues, then the analyticity of the discretized determinant implies that the zeros of the determinant of the Nyström discretization of the linear operator converge exponentially to the true Dirichlet eigenvalues for Laplace’s equation.

The proof presented in [54] applies to the BIE approach for computing the Dirichlet eigenvalues of the Stokes operator as well. The result is summarized below.

Theorem 11**.**

Suppose that $\Omega$ is a simply connected domain with an analytic boundary. Let $k_{j}^{2}$, $j=1,2,\ldots M$ denote all the Dirichlet eigenvalues of Stokes equation on $\Omega$ contained in the interval $[a,b]$. Suppose further that all the eigenvalues have multiplicity $1$. Let $f^{N}(k)=\text{det}(I+D^{N}_{k})$, where $D^{N}_{k}$ is the Nyström discretization of $-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ with Kress quadrature. Then there exists $N_{0}\in\mathbb{N}$ such that for all $N>N_{0}$, $f^{N}(k)$ has exactly $M$ zeros on the interval $[a,b]$. Let $\omega_{j}$, $j=1,2\ldots M$ denote the zeros of $f^{N}$. Furthermore, there exist constants $a>0$ and $C$, such that $\sup_{j=1}^{M}|\omega_{j}-k_{j}|<Ce^{-aN}$.

Proof.

The proof follows from small modifications of the proofs contained in [54]. ∎

Remark 7**.**

In practice, using Kress quadrature for large problems is problematic owing to the global nature of the quadrature rule. First, the use of a global rule does not allow for adaptive refinement at a complicated, local feature of the boundary. Second, the integration weight in each entry of the matrix depends on both the column and the row in a non-separable way. As a result, the fast multipole method is not directly applicable to the resulting matrix and many fast-direct methods for computing the determinant lose efficiency (for instance, the reasoning behind the use of a proxy surface [14] no longer holds). Over the last two decades, many high-order quadrature methods which are compatible with the fast multipole method and fast-direct methods have been developed. Our numerical experiments, see section 4.2.1, suggest that the zeros of the determinants of linear systems discretized using these quadrature methods are also high order approximations of Dirichlet eigenvalues for the Stokes operator — the error is observed to be proportional to the quadrature error for the eigenfunction ${\boldsymbol{t}}^{-}$ associated with the eigenvalue. We leave a proof of this to future work.

Remark 8**.**

The same analysis does not carry through for the operator $\mathcal{I}-2{\mathcal{D}}_{k}-2i\eta{\mathcal{S}}_{k}-2{\mathcal{W}}$, since ${\mathcal{S}}_{k}$ is not a trace class operator. For brevity, let $\mathcal{C}_{k}=-2{\mathcal{D}}_{k}-2i\eta{\mathcal{S}}_{k}-2{\mathcal{W}}$. The operator $\mathcal{C}_{k}$ is in ${\mathcal{J}}_{2}(Y)$ where ${\mathcal{J}}_{2}(Y)$ is the space of Hilbert-Schmidt operators on $Y$ (the singular values of the operator are square summable, as opposed to being summable). Thus the Fredholm determinant of $\mathcal{I}+\mathcal{C}_{k}$ is not necessarily finite. However, as noted in [54], the convergence result theorem 11 should be true up to a logarithmic factor in the rate of convergence, since the singular values of the operator $\mathcal{C}_{k}$ decay like $\frac{1}{n}$, and the Fredholm determinant diverges logarithmically. In section 4.2.1, we demonstrate this fact numerically on the annulus, where the eigenvalues are analytically known.

4 Numerical Results

In this section, we demonstrate the analytical claims above with numerical examples and highlight the performance of the BIE approach with demonstrations on domains of analytical and practical interest. The software used to generate the figures is available online 111https://doi.org/10.5281/zenodo.2641296.

4.1 Numerical methods

First, we describe the numerical tools needed to compute Stokes eigenvalues in a BIE framework.

4.1.1 Discretizing the BIE

In order to turn the BIEs analyzed above into discrete linear systems, we require some standard techniques from the BIE literature.

Let the boundary be divided into $N_{p}$ panels. We parameterize panel $j$ as ${\boldsymbol{x}}_{j}(t)$, with $t$ ranging over the interval $[-1,1]$. Each component of ${\boldsymbol{x}}_{j}$ is taken to be a polynomial interpolant over the standard 16th-order Legendre nodes on $[-1,1]$, denoted by $t_{n}$, so that the total number of discretization points is $N=16N_{p}$. See fig. 3 for an example discretization.

Another important quantity below is the arc-length density of a panel, which we denote by $s_{j}(t):=|{\boldsymbol{x}}^{\prime}_{j}(t)|$. Finally, we denote the set of panels which are adjacent to panel $j$ by $A(j)$. On a closed curve, $A(j)$ contains two integers.

The integral kernels of the single and double layer potentials have weak singularities of the form $|{\boldsymbol{x}}-{\boldsymbol{y}}|^{p}\log|{\boldsymbol{x}}-{\boldsymbol{y}}|$ for some $p\in\mathbb{N}_{0}$, which require special quadrature rules to achieve high-order accuracy. In the examples below, we use generalized Gaussian quadrature (GGQ) [11]. To demonstrate the idea, we consider evaluating the convolution of a kernel $K$ with a density $\sigma$ at the boundary node ${\boldsymbol{x}}_{j}(t_{l})$. GGQ is a Nyström-type discretization — the density is approximated by its values at the discretization nodes, which we denote by $\sigma_{qp}:=\sigma({\boldsymbol{x}}_{q}(t_{p}))$. The basis of a GGQ rule is a set of support nodes and weights for the contribution to the integral from the “self” panel (panel $j$) and the adjacent panels (with index in $A(j)$).

For the self panel, there is a special set of nodes and weights for each interpolation point. Denote the nodes and weights for interpolation point $l$ by $t^{(l)}_{n}$ and $w^{(l)}_{n}$, respectively, with $1\leq l\leq 16$ and $1\leq n\leq N_{s}$. The adjacent panels are handled by a single set of over-sampled support nodes and weights. We denote these nodes and weights by $\tilde{t}_{n}$ and $\tilde{w}_{n}$, respectively, for $1\leq n\leq N_{a}$. For the 20th-order rule we used, $N_{s}=16$ and $N_{a}=48$. The contribution of other panels is assumed to be given to high accuracy by the standard Gauss-Legendre weights, which we denote by $w_{n}$. Adding these contributions together, we obtain the quadrature

[TABLE]

where ${\bf B}^{(l)}$ and ${\bf C}$ are interpolation matrices from the standard Legendre nodes to the self and adjacent panel support nodes, respectively. Observe that the quadrature is linear in $\sigma_{qp}$. In practice, we pre-compute and store the self and adjacent matrix entries for each interpolation point, which is a parallelizable $O(N)$ calculation. The “far” interactions are computed on-the-fly.

Remark 9**.**

We ensure that “far” interactions are handled to high precision by requiring that no two adjacent panels differ in length by more than a factor of 2. On a domain which does not nearly self-intersect this guarantees that no “far” interactions occur which are much closer than 1/2 of a panel away (assuming panels are relatively flat). Because the location of the singularity is bounded away from the panel and the smooth rule is of high order, we obtain a quadrature rule with sufficient precision.

The overall order of accuracy of the GGQ we use is 20th-order, up to the precision of the “far” interactions.

4.1.2 Fast determinant method

Once the discretization is set, we can form a compressed representation of the system matrix using recursive skeletonization [24]. We use the implementation of this procedure included in the fast linear algebra in MATLAB (FLAM) package [23]. At low-to-medium frequencies, the scaling of the recursive skeletonization algorithm is $O(N\log N)$ in operation count and storage and, by using a generalization of the Sylvester determinant formula, allows for a fast determinant calculation in $O(N\log N)$ time as a follow-up step. At higher-frequencies, the recursive skeletonization procedure, which is based on the assumption that off-diagonal blocks of the matrix are of low rank, breaks down and does not offer a speed advantage. These algorithms take a precision parameter ${\epsilon_{\texttt{FLAM}}}$ which determines the accuracy to which any sub-blocks of the matrix should be compressed. In all experiments, we set ${\epsilon_{\texttt{FLAM}}}=10^{-14}$.

The compressed representation also allows for fast applications of the system matrix, its transpose, the inverse of the system matrix, and the inverse transpose to vectors. In particular, this allows us to estimate the smallest singular values by performing randomized subspace iteration, see [20, Algorithm 4.4], on the inverse operator. Below, we use the smallest singular value as a measure of the quality of the eigenvalues found by approximating the roots of the determinant. We also evaluate the second smallest singular value if the root finding procedure suggests a possible double root.

4.1.3 Interpolation and root-finding

To estimate the eigenvalues, we fit a Chebyshev interpolant to the discretized determinant as a function of $k$ on intervals. This is done adaptively so that the Chebyshev coefficients of the determinant have decayed to the point that the ratio of the last coefficient to the largest coefficient is below some threshold. In all experiments, we set this threshold as ${\epsilon_{\texttt{cheb}}}=10^{-13}$. We perform this fit using the chebfun utility in the package of the same name [16] so that we can make use of the roots utility to approximate the roots of the determinant.

The roots utility returns the roots of the polynomial in the complex plane, with some minimal internal processing to remove spurious roots. Because our numerical determinant evaluation is somewhat noisy and we fit the function up to precision ${\epsilon_{\texttt{cheb}}}$, we perform some further post-processing to eliminate remaining spurious roots. Let $k^{(l)}_{\texttt{cheb}}$ denote the roots of the interpolants. We ignore any of the returned roots with $|\operatorname{Im}(k^{(l)}_{\texttt{cheb}})|>\sqrt{{\epsilon_{\texttt{cheb}}}}$, as these are too far from real-valued to be non-spurious. For the remaining roots, we consider the properties of $\real(k^{(l)}_{\texttt{cheb}})$. We inspect any pairs of roots $(k^{(p)}_{\texttt{cheb}},k^{(q)}_{\texttt{cheb}})$ for which $|\real(k^{(p)}_{\texttt{cheb}}-k^{(q)}_{\texttt{cheb}})|<\sqrt{{\epsilon_{\texttt{cheb}}}}$, as these are possibly spurious double roots. For these pairs, we compare the right singular vector of the appropriate BIE operator corresponding to the smallest singular value for each of $\real(k^{(p)}_{\texttt{cheb}})$ and $\real(k^{(q)}_{\texttt{cheb}})$, which we denote by ${\boldsymbol{v}}_{p}$ and ${\boldsymbol{v}}_{q}$. If $\|{\boldsymbol{v}}_{p}-{\boldsymbol{v}}_{q}{\boldsymbol{v}}_{q}^{*}{\boldsymbol{v}}_{p}\|<10^{-5}$, then we consider the pair to be spurious. For these near double roots, we also check that there is no two dimensional null-space corresponding to the root by estimating the second smallest singular value of the BIE operator. If this is larger than $10^{-5}$, then we declare it to be a simple root.

We can obtain an a posteriori estimate of the error in a computed root as follows. Let $f$ denote an analytic function, $P$ be the polynomial interpolant of that function over some interval, $\delta f=f-P$ be the difference, and $k_{\texttt{cheb}}$ denote a computed root of $P$ which is simple (i.e. assume that $P^{\prime}(k_{\texttt{cheb}})\neq 0$). The algorithm used by chebfun to approximate the roots of $P$ is backward stable [39]. Therefore the error in the roots will be small relative to the error of the fit and we set $P(k_{\texttt{cheb}})=0$ below. Suppose that $f(k_{\texttt{cheb}}+\delta k)=0$ for some small $\delta k$. Then

[TABLE]

In practice, we can obtain an approximate upper bound for $|\delta f(k_{\texttt{cheb}}+\delta k)|$ as ${\epsilon_{\texttt{cheb}}}\|P\|_{\infty}$ so that ${\epsilon_{\texttt{cheb}}}\|P\|_{\infty}/|P^{\prime}(k_{\texttt{cheb}})|$ provides an approximate upper bound for the error in the root.

4.2 Eigenvalues of an annulus

We test our numerical machinery and validate our analytical and numerical claims by comparing the results to the true eigenvalues on the annulus which are known analytically (see appendix A). In all of the examples below, we work on the annulus $r_{1}<r<r_{2}$ with $r_{1}=1$ and $r_{2}=1.7$. If the inner boundary is discretized using $N_{1}$ panels, then the outer boundary is discretized using $N_{2}=\lceil r_{2}/r_{1}N_{1}\rceil+1$ panels to ensure that the panels are approximately the same length on both the boundaries. The total number of discretization points is then given by $N=16(N_{1}+N_{2})$. Let $D^{N}_{k}$ denote the linear system corresponding to the Nyström discretization of $-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ using generalized Gaussian quadrature, and let $C^{N}_{k}$ denote the linear system corresponding to the Nyström discretization of $-2{\mathcal{D}}_{k}-2i{\mathcal{S}}_{k}-2{\mathcal{W}}$. Let $f_{D}^{N}(k)=\text{det}(I+D^{N}_{k})$, and $f_{C}^{N}(k)=\text{det}(I+C^{N}_{k})$.

4.2.1 Convergence study

We demonstrate that for sufficiently large $N$, if $k_{0}$ is a Dirichlet eigenvalue of the annulus, then $f_{D}^{N}(k_{D})=0$ and $f_{C}^{N}(k_{C})=0$ where $|k_{D}-k_{0}|\lesssim N^{-20}$, and $|k_{C}-k_{0}|\lesssim N^{-20}$. Recall that the GGQ we use has an expected order of convergence of $N^{-20}$ for evaluating convolutions with such integral kernels, so that the error of the roots is observed to have the same order as the quadrature rule. In fig. 4, we show this result for $k_{0}=13.48025717955055$ and plot the errors $|k_{D}-k_{0}|$ and $|k_{C}-k_{0}|$ as a function of $N$.

4.2.2 Spurious eigenvalues

As noted in section 3.5.1, if $k_{0}$ is a Neumann eigenvalue corresponding to the interior inclusion, which in our case is the disk $r\leq r_{1}$, then $f_{D}^{N}(k_{D})=0$ with $|k_{D}-k_{0}|=O(\varepsilon)$, even though $k_{0}$ is not a Dirichlet eigenvalue of the annulus, i.e. the integral equation $-\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ has a spurious eigenvalue. In fig. 5, we demonstrate this result and also show that $f_{C}^{N}(k_{0})\neq 0$, i.e., the combined field representation is robust and invertible at all values of $k$ which are not the Dirichlet eigenvalues of the annulus.

4.2.3 Speed

In this section, we demonstrate the $O(N\log{N})$ scaling of evaluating $f^{N}_{C}(k)$ as long as $N$ is large enough to resolve the interactions at the Helmholtz parameter $k$. When $N$ is smaller than that, we observe a worse scaling since the assumption that far-interactions are low-rank is no longer valid at the tolerance of FLAM. We plot the timing results corresponding to three different values of $k$ in fig. 6. The times are as recorded for a laptop with 16Gb of RAM and an Intel Core i7-6600U CPU at 2.60GHz with 4 cores.

4.3 Eigenvalues of a barbell-shaped domain

We consider the barbell-shaped domain in fig. 7(a). This domain is the union of a square of side-length 6, a square of side-length 3, and a “bridge” connecting them of height 1 and width 5/2. For the sake of simplicity, we round the corners of the domain to obtain a smooth object. Applying the approach described in [17], the corners of the domain are rounded by convolving with the Gaussian kernel

[TABLE]

with $h\approx 0.06$. This leaves the domain unperturbed to high precision outside of a radius of $0.1$ around each corner. The eigenfunctions of such a domain display the well-known localization property [53]: many of the eigenfunctions are approximately supported within one of the squares. We compute these eigenfunctions corresponding to eigenvalues $k^{2}$ with $k$ in the range $0.5\leq k\leq 6.5$.

The panels are divided adaptively so that the smallest panels in the rounded corners are smaller than $10^{-2}$, which keeps the panels relatively flat. This results in $N_{p}=412$ after enforcing the level-restriction property described in remark 9 and enforcing that no panel is larger than one wavelength for the largest $k$ (here $\lambda=2\pi/6.5$).

As this is a simply-connected domain, the eigenvalues are estimated by finding the values $k$ for which $\mathcal{I}-2{\mathcal{D}}_{k}-2{\mathcal{W}}$ is non-invertible. Let $f^{N}(k)=\det(\mathcal{I}^{N}-2{\mathcal{D}}_{k}^{N}-2{\mathcal{W}}^{N})$. To find the roots of $f^{N}(k)$, we fit a chebfun representation of $f^{N}(k)$ on each of the intervals $[j/2,(j+1)/2]$ for $j=1,\ldots,12$. We plot the absolute value of the Chebyshev coefficients (normalized by the absolute value of the first coefficient) of $f^{N}(k)$ on the intervals $[0.5,1.0]$, $[3.0,3.5]$, and $[6.0,6.5]$ in fig. 9(b). As expected, the coefficients decay exponentially to zero, with more terms required at higher frequencies.

We compute the roots of these Chebyshev interpolants and apply the post-processing described above. There were 135 total roots: 3 were removed because the imaginary part was too large and 13 pairs were found with values within $\sqrt{{\epsilon_{\texttt{cheb}}}}$ of each other. For these 13 pairs, none represented two distinct eigenvalues or a double root. This leaves 119 roots in the range $0\leq k\leq 6.5$. We plot the smallest singular value of $\mathcal{I}^{N}-{\mathcal{D}}_{k}^{N}-2{\mathcal{W}}^{N}$ for each of these roots in fig. 9(a) and plot the vorticity of the eigenfunctions in fig. 8. The singular values suggest that the quality of the eigenvalues is good. From the plots, we see that localization occurs until about the 100th eigenvalue.

4.4 Eigenvalues of a domain with several inclusions

We now consider the multiply-connected domain in fig. 7(b). The domain is defined by a smooth rectangular region of width 3 and height 2, with an array of randomly rotated “starfish” shapes removed. Such shapes are of interest in materials design, see, for instance, [41]. We compute the eigenfunctions corresponding to eigenvalues $k^{2}$ with $k$ in the range $3\leq k\leq 9$ (this range includes the smallest eigenvalue).

For this smooth shape, ensuring that no panel is larger than one wavelength for the largest $k$ (here $\lambda=2\pi/9$) is sufficient to resolve the object to high precision. After enforcing the level-restriction property described in remark 9, we end up with $N_{p}=224$.

As this is a multiply-connected domain, the eigenvalues are estimated by finding the values $k$ for which $\mathcal{I}-2{\mathcal{D}}_{k}-2i{\mathcal{S}}_{k}-2{\mathcal{W}}$ is non-invertible. Let $f^{N}(k)=\det(\mathcal{I}^{N}-2{\mathcal{D}}_{k}^{N}-2i{\mathcal{S}}_{k}^{N}-2{\mathcal{W}}^{N})$. To find the roots of $f^{N}(k)$, we fit a chebfun representation of $f^{N}(k)$ on each of the intervals $[j/2,(j+1)/2]$ for $j=6,\ldots,13$ and the intervals $[j/4,(j+1)/4]$ for $j=28,\ldots,35$. It should be noted that, due to the relative sizes of the domains, this represents a lower frequency problem than that for the barbell when measured in the number of wavelengths across the object. Thus, the use of a finer grid in frequency results from the difficulty in resolving the Fredholm determinant for this problem, which has a larger dynamical range than that for the barbell. We plot the absolute value of the Chebyshev coefficients of $f^{N}(k)$ on the intervals $[3,3.5]$, $[6.5,7]$, and $[8.75,9]$ in fig. 9(b). As expected, the coefficients decay exponentially to zero, with more terms required at higher frequencies (note that the interval $[8.75,9]$ is smaller than the others).

We compute the roots of these Chebyshev interpolants and apply the post-processing described above. There were 103 total roots: 21 were removed because the imaginary part was too large and 6 pairs were found with values within $\sqrt{{\epsilon_{\texttt{cheb}}}}$ of each other. For these 6 pairs, none represented two distinct eigenvalues or a double root. This leaves 76 roots in the range $0\leq k\leq 9$. We plot the smallest singular value of $\mathcal{I}^{N}-{\mathcal{D}}_{k}^{N}-2i{\mathcal{S}}_{k}^{N}-2{\mathcal{W}}^{N}$ for each of these roots in fig. 11(a) and plot the vorticity of the eigenfunctions in fig. 10. In the vorticity plots, we observe a different type of localization property than that seen in the barbell, with many of the eigenfunctions approximately supported in a small, connected subset of the domain. This is consistent with other studies [18, 36].

The singular values suggest that the quality of the eigenvalues is good, with a few outliers. To explain these outliers, we consider two quantities which affect the singular value at a computed root. As described above, we can approximate the error in the computed root at $k_{\texttt{cheb}}$ by ${\epsilon_{\texttt{cheb}}}\|P\|_{\infty}/|P^{\prime}(k_{\texttt{cheb}})|$, where $P$ is the interpolating polynomial. The singular value estimate itself is affected by the error incurred in applying the inverse of the compressed BIE matrix, which can be hard to quantify [24]. We approximate this error by $O(\sqrt{N}){\epsilon_{\texttt{FLAM}}}$ and assume this is the order of the error in the singular value estimate. We plot the maximum of these two estimates along with the computed singular values in fig. 12(a). There is a reasonably good correlation between the maximum of the error estimates and the observed smallest singular value for the BIE, especially for larger errors.

The worst outliers are from the left half of the interval $[4.5,5]$. Because the determinant is much larger on the right half than the left half of $[4.5,5]$, we can improve the estimate for the error in the roots by subdividing the interval. We plot the smallest singular value of the BIE for the roots obtained by fitting a polynomial on $[4.5,4.75]$ in fig. 12(b); the roots on the refined interval are of significantly higher quality.

Remark 10**.**

The above experience suggests that the ratio $\|P\|_{\infty}/|P^{\prime}(k_{\texttt{cheb}})|$ is a useful diagnostic for performing automated eigenvalue estimation. Note that at a multiple root, this ratio will be more difficult to bound.

5 Conclusion

In the preceding, we have demonstrated a BIE framework for computing the eigenvalues of the Stokes operator in the plane which is robust and scalable. To justify the approach, we developed a uniqueness theory for the oscillatory Stokes equations in exterior domains in analogy with the discussion of the Helmholtz equation in [15]. This lead to the primary theoretical results of the paper which show that the BIEs resulting from double layer and combined-field representations of the velocity field are not invertible precisely when $k^{2}$ is an eigenvalue on simply connected and multiply connected domains, respectively. As in [54], the costliness of performing the nonlinear minimization associated with computing these eigenvalues can be alleviated by computing instead the approximate zeros of the discrete Fredholm determinant.

The results of this paper can be extended in a number of ways. The theory extends directly to three dimensions, where computational efficiency and numerical implementation will be the primary concern. In the numerical examples above, we consider only domains with differentiable boundaries for simplicity. When using the eigenfunctions as a trial basis for simulating the Navier–Stokes equations [6], it is necessary to handle domains with corners because the domains of interest arise from domain decomposition, e.g. dividing a larger domain into quadrilaterals. Fortunately, there has been recent progress toward efficient discretization of the layer potentials of elliptic operators on domains with corners [22, 49, 45, 21] which makes the solution of such problems tractable. Of course, the theoretical considerations are different for such domains. Computing the Stokes eigenvalues of regions with corners is the topic of ongoing research.

As noted in the introduction, the eigenvalues of the Stokes operator are equivalent to the so-called “buckling” eigenvalues of the biharmonic operator on simply connected domains [29]. This can be seen through the stream function formulation of the oscillatory Stokes equation, i.e. setting ${\boldsymbol{u}}=\nabla^{\perp}\psi$ where $\psi$ now satisfies

[TABLE]

Note that the buckling problem enforces the clamped, or first Dirichlet, boundary condition on $\psi$

[TABLE]

On a multiply connected domain, there are Stokes eigenfunctions which do not have a corresponding clamped stream function, so that the buckling eigenvalues are a subset of the Stokes eigenvalues. By adapting the approach of [44], a suitable layer potential representation of the buckling problem can be derived based on the oscillatory Stokes layer potentials. This is the subject of a follow-up paper which is in preparation.

There are some interesting questions to answer on the use of Fredholm determinants in numerical calculations. As observed in [54], the combined-field representation causes some difficulty in that the Fredholm determinant is not defined when the single layer, which is not trace-class, is included. Zhao and Barnett [54] suggest looking into representations of the form $\mathcal{I}-2{\mathcal{D}}-2i\eta{\mathcal{S}}^{2}-2{\mathcal{W}}$, which would have a well-defined Fredholm determinant. The relative performance of such an approach should be explored. Further, as discussed above, the convergence of the determinant of integral equations discretized with panel-corrected schemes (as described in section 4) is yet to be proved. When addressing high frequency problems, the fast-direct solver used to evaluate determinants in this paper will no longer have near-linear scaling [24]. The design of fast-direct solvers in this regime is the subject of ongoing research, and to the best of our knowledge, fast determinant computation at high frequency is an open question. Finally, it is worth exploring alternatives to the Fredholm determinant which perform well for layer potentials that are not trace class on the boundary or for problems at higher frequencies.

6 Acknowledgements

The authors would like to thank Alex Barnett, Leslie Greengard, and Shidong Jiang for many useful discussions.

T. Askham was supported by the Air Force Office of Scientific Research under Grant FA9550-17-1-0329.

Appendix A Dirichlet eigenvalues and eigenfunctions on the annulus

In this section, we compute some of the Dirichlet eigenvalues corresponding to a subset of the radially symmetric eigenfunctions on the annulus. In polar coordinates $(r,\theta)$, consider the annulus defined by $R_{1}<r<R_{2}$. Suppose that ${\boldsymbol{u}}$ is of the form

[TABLE]

and $p=0$.

Clearly, this pair satisfies the osciallatory Stokes equation with parameter $k$, since $J_{0}(kr)$ and $H_{0}(kr)$ satisfy the Helmholtz equation on the annulus.

Let $\hat{r},\hat{\theta}$ denote the unit vectors in polar coordinates. A simple calculation shows that

[TABLE]

This in particular implies that on $r=R_{1}$, $u_{\theta}$ takes on the constant value, $k(\alpha H_{0}^{\prime}(kR_{1})+\beta J_{0}^{\prime}(kR_{1}))$. Similarly, on $r=R_{2}$, $u_{\theta}$ takes on the constant value, $k(\alpha H_{0}^{\prime}(kR_{2})+\beta J_{0}^{\prime}(kR_{2}))$.

Thus, if $k$ satisfies,

[TABLE]

and for those values of $k$ if $\alpha,\beta$ are non-zero solutions to system of equations

[TABLE]

then $k$ is a Dirichlet eigenvalue and ${\boldsymbol{u}}$ defined by 61 is the corresponding eigenfunction.

Appendix B Neumann eigenvalues and eigenfunctions on the unit disk

In this section, we derive an analytical expression which can be used to compute some of the radially symmetric Neumann eigenvalues on the unit disk for the Stokes operator.

Suppose that ${\boldsymbol{u}}$ is of the form

[TABLE]

and the pressure is given by $p=0$, as ${\boldsymbol{u}}$ satisfies $(\Delta+k^{2}){\boldsymbol{u}}=0$.

Then, the surface traction ${\boldsymbol{t}}$ on the disk of radius $r$ is given by

[TABLE]

Thus, $k$ which satisfies

[TABLE]

are Neumann eigenvalues on the unit disk. The first $4$ roots of the 67 are in table 1.