Wave Enhancement through Optimization of Boundary Conditions

TL;DR

This paper introduces an efficient method to optimize wave transmission in cavities by adjusting boundary conditions, utilizing recent theoretical insights and metasurface technology to enhance wave control at specific frequencies.

Contribution

It presents a novel approach combining boundary condition optimization with metasurface design, enabling precise control of wave transmission in cavities.

Findings

Effective boundary condition switching improves transmission signals.

Numerical experiments validate the accuracy of the proposed method.

The approach leverages eigenvalue monotonicity and Green's function sensitivity.

Abstract

It is well known that changing boundary conditions for the Laplacian from Dirichlet to Neumann can result in significant changes to the associated eigenmodes, while keeping the eigenvalues close. We present a new and efficient approach for optimizing the transmission signal between two points in a cavity at a given frequency, by changing boundary conditions. The proposed approach makes use of recent results on the monotonicity of the eigenvalues of the mixed boundary value problem and on the sensitivity of the Green s function to small changes in the boundary conditions. The switching of the boundary condition from Dirichlet to Neumann can be performed through the use of the recently modeled concept of metasurfaces which are comprised of coupled pairs of Helmholtz resonators. A variety of numerical experiments are presented to show the applicability and the accuracy of the proposed new methodology.

Figures

Tables

Table 1. Table 1 : We see Algorithm 4 performing on the unit circle with k ⋆ = 1 subscript 𝑘 ⋆ 1 k_{\star}=1 , x S = ( 0 , 0 ) T subscript 𝑥 S superscript 0 0 T x_{\textrm{S}}=(0,0)^{\mathrm{T}} , y ∈ { ( 0 , r ) T ∈ ℝ 2 ∣ r > 0 } 𝑦 conditional-set superscript 0 𝑟 T superscript ℝ 2 𝑟 0 y\in\{(0,r)^{\mathrm{T}}\in\mathbb{R}^{2}\mid r>0\} , C tol = 10 − 3 subscript 𝐶 tol superscript 10 3 C_{\mathrm{tol}}=10^{-3} and ε = 0.1 𝜀 0.1 \varepsilon=0.1 . Z Dirichlet k ​ ( x S , y ) superscript subscript Z Dirichlet 𝑘 subscript 𝑥 S 𝑦 \mathrm{Z}_{\mathrm{Dirichlet}}^{k}(x_{\textrm{S}},y) represent the Zaremba function on the partition ( ∂ Ω , ∅ ) Ω (\partial\Omega,\varnothing) of the boundary and Z End k ​ ( x S , y ) superscript subscript Z End 𝑘 subscript 𝑥 S 𝑦 \mathrm{Z}_{\mathrm{End}}^{k}(x_{\textrm{S}},y) represents the Zaremba function on the final partition, where the final partition is made out of two boundary intervals, one with Dirichlet boundary conditions and the other with Neumann boundary conditions. θ center ∈ [ 0 , 2 ⋅ p ​ i ) subscript 𝜃 center 0 ⋅ 2 𝑝 𝑖 \theta_{\mathrm{center}}\in[0,2\cdot pi) represents the angle of the center of the Neumann boundary intervals and l N subscript 𝑙 N l_{\mathrm{N}} its length. The shown values are the real, rounded values of the numerical results.

	$r=0.1$	$r=0.25$	$r=0.5$	$r=0.75$	$r=0.9$
$Z_{\text{Dirichlet}}^{k_{\star }}(x_{\text{S}},y)$	-0.412	-0.261	-0.138	-0.059	-0.022
$Z_{\text{End}}^{k_{\star }}(x_{\text{S}},y)$	-1288	-1438	-1634	-1754	-1788
$\left|\frac{Z_{\text{End}}^{k_{\star }}(x_{\text{S}},y)}{Z_{\text{Dirichlet}}^{k_{\star }}(x_{\text{S}},y)}\right|$	3123	5503	11824	29623	81687
${\theta }_{\mathrm{center}}$	$0.50\pi$	$0.50\pi$	$0.50\pi$	$0.50\pi$	$0.50\pi$
$l_N$	$1.32\pi$	$1.32\pi$	$1.32\pi$	$1.32\pi$	$1.32\pi$

Table 2. Table 2 : We see Algorithm 4 performing on the unit circle with k ⋆ = 15.4 subscript 𝑘 ⋆ 15.4 k_{\star}=15.4 , x S = ( 0 , 0 ) T subscript 𝑥 S superscript 0 0 T x_{\textrm{S}}=(0,0)^{\mathrm{T}} , y ∈ { ( 0 , r ) T ∈ ℝ 2 ∣ r > 0 } 𝑦 conditional-set superscript 0 𝑟 T superscript ℝ 2 𝑟 0 y\in\{(0,r)^{\mathrm{T}}\in\mathbb{R}^{2}\mid r>0\} , C tol = 10 − 3 subscript 𝐶 tol superscript 10 3 C_{\mathrm{tol}}=10^{-3} and ε = 0.05 𝜀 0.05 \varepsilon=0.05 . Z Dirichlet k ​ ( x S , y ) superscript subscript Z Dirichlet 𝑘 subscript 𝑥 S 𝑦 \mathrm{Z}_{\mathrm{Dirichlet}}^{k}(x_{\textrm{S}},y) , Z End k ​ ( x S , y ) superscript subscript Z End 𝑘 subscript 𝑥 S 𝑦 \mathrm{Z}_{\mathrm{End}}^{k}(x_{\textrm{S}},y) , θ center subscript 𝜃 center \theta_{\mathrm{center}} , and l N subscript 𝑙 N l_{\mathrm{N}} are defined as in Table 1 . The shown values are the real, rounded values of the numerical results.

	$r=0.1$	$r=0.25$	$r=0.5$	$r=0.75$	$r=0.9$
$Z_{\text{Dirichlet}}^{k_{\star }}(x_{\text{S}},y)$	0.341	-0.188	0.157	-0.085	0.118
$Z_{\text{End}}^{k_{\star }}(x_{\text{S}},y)$	36.341	-14.271	116.08	-15.811	232.28
$\left|\frac{Z_{\text{End}}^{k_{\star }}(x_{\text{S}},y)}{Z_{\text{Dirichlet}}^{k_{\star }}(x_{\text{S}},y)}\right|$	106.6	76.09	739.0	186.8	1962
${\theta }_{\mathrm{center}}$	$0.50\pi$	$1.90\pi$	$0.50\pi$	$0.46\pi$	$0.50\pi$
$l_N$	$0.064\pi$	$0.064\pi$	$0.064\pi$	$0.064\pi$	$0.064\pi$

Equations

$$\{\begin{array}{rlrlrl}\left(△+k^2\right){\mathrm{Z}}_{x_{\text{S}}}^k(x_{\text{S}},y)& ={\delta }_0(x_{\text{S}}-y)& & \text{for}& & y\in \Omega ,\\ {\mathrm{Z}}_{x_{\text{S}}}^k(x_{\text{S}},y)& =0& & \text{for}& & y\in {\Gamma }_{\text{D}},\\ {\partial }_{{\nu }_y}{\mathrm{Z}}_{x_{\text{S}}}^k(x_{\text{S}},y)& =0& & \text{for}& & y\in {\Gamma }_{\text{N}}.\end{array}$$

$$Z^k(x_s,\cdot )={\Gamma }^k(x_s,\cdot )+{\mathrm{R}}^k(x_s,\cdot )$$

$$\{\begin{array}{rlrl}\left(△+k^2\right){\mathrm{R}}^k(x_s,y)& =0,& & \text{mbox}in\Omega ,\\ {\mathrm{R}}^k(x_s,y)& =-{\Gamma }^k(x_s,y),& & \text{mbox}on{\Gamma }_D,\\ {\partial }_{{\nu }_y}{\mathrm{R}}^k(x_s,y)& =-{\partial }_{{\nu }_y}{\Gamma }^k(x_s,y)& & \text{mbox}on{\Gamma }_N.\end{array}$$

$$\{\begin{array}{rlrlrl}-△u& =k^{2}u& & \text{in}& & \Omega ,\\ u& =0& & \text{on}& & {\Gamma }_{\text{D}},\\ {\partial }_{{\nu }_y}u& =0& & \text{on}& & {\Gamma }_{\text{N}}.\end{array}$$

$$0$$

$$\{{\alpha }^2\in (0,\infty )\mid \exists n\in {\mathbb{N}}_0:\alpha \text{ is positive root of }{J}_n(x)\},$$

$$\{{\alpha }^2\in (0,\infty )\mid \exists n\in {\mathbb{N}}_0:\alpha \text{ is positive root of }{J}_n^{\prime }(x)\}.$$

$${\lambda }_j^{{\Gamma }_{\text{D}}}<{\lambda }_j^{{\Gamma }_{D^{\prime }}}\text{ for all }j\in \mathbb{N}.$$

$${\lambda }_j^{\varnothing }<{\lambda }_j^{{\Gamma }_{\text{D}}}<{\lambda }_j^{\partial \Omega },\text{ for all }j\in \mathbb{N}.$$

$$u(x)={\int }_{\partial \Omega }{\Gamma }^k(x,y)\psi (y)\mathrm{d}{\sigma }_y,$$

$${\mathcal{S}}_{{\Gamma }_{\text{D}}}^k[\psi ](x)$$

$$\partial {\mathcal{S}}_{{\Gamma }_{\text{D}}}^k[\psi ](x)$$

$$\mathcal{A}(k)\left[\begin{array}{c}\psi \mid {}_{{\Gamma }_{\text{D}}}\\ \psi \mid {}_{{\Gamma }_{\text{N}}}\end{array}\right]:=\left[\begin{array}{cc}{\mathcal{S}}_{{\Gamma }_{\text{D}}}^k& {\mathcal{S}}_{{\Gamma }_{\text{N}}}^k\\ \partial {\mathcal{S}}_{{\Gamma }_{\text{D}}}^k& -\frac{1}{2}{\mathrm{I}}_{L_{\omega }^2({\Gamma }_{\text{N}})}+({\mathcal{K}}_{{\Gamma }_{\text{N}}}^k{)}^{\ast }\end{array}\right]\left[\begin{array}{c}\psi \mid {}_{{\Gamma }_{\text{D}}}\\ \psi \mid {}_{{\Gamma }_{\text{N}}}\end{array}\right].$$

$$\displaystyle\mbox{ $k\mapsto\mathcal{A}(k)$ are the square roots of the Zaremba eigenvalues''}.$$

$$u(y)$$

$$\psi \mid {}_{{\Gamma }_{\text{D}}}(y)$$

$$\psi \mid {}_{{\Gamma }_{\text{N}}}(y)$$

$$k_j^{\epsilon }-k_j^0=\frac{1}{2\pi i}\mathrm{tr}{\int }_{\partial V}(\omega -k_j^0){\mathcal{A}}_{\epsilon }(\omega {)}^{-1}{\partial }_{\omega }{\mathcal{A}}_{\epsilon }(\omega )\mathrm{d}\omega .$$

$$k_j^{\epsilon }-k_j^0\approx \frac{-1}{2\pi i}\mathrm{tr}{\int }_{\partial V}{\mathcal{A}}_0(\omega {)}^{-1}({\mathcal{A}}_{\epsilon }(\omega )-{\mathcal{A}}_0(\omega ))\mathrm{d}\omega ,$$

$$k_j^{\epsilon }-k_j^0=\frac{-1}{2\pi i}\mathrm{tr}{\int }_{\partial V}(\mathrm{I}+(\omega -k_j^0){\mathcal{A}}_{\epsilon }(k_j^0{)}^{-1}{\partial }_{\omega }{\mathcal{A}}_{\epsilon }(k_j^0){)}^{-1}\mathrm{d}\omega$$

$$\displaystyle\times\bigg{[}1+\mathcal{O}(\mathrm{tr}\int_{\partial V}(\mathrm{I}+(\omega-k_{j}^{0})\mathcal{A}_{\varepsilon}(k_{j}^{0})^{-1}\partial_{\omega}\mathcal{A}_{\varepsilon}(k_{j}^{0}))^{-1}\,\mathrm{d}\omega)\bigg{]}.$$

$${\mathcal{A}}_{\epsilon }(\omega )={\mathcal{A}}_{\epsilon }(k_j^0)+(\omega -k_j^0){\partial }_{\omega }{\mathcal{A}}_{\epsilon }(k_j^0)+{\mathcal{B}}_{\epsilon }(\omega ),$$

$$\Vert {\left({\mathcal{A}}_{\epsilon }(k_j^0)+(\omega -k_j^0){\partial }_{\omega }{\mathcal{A}}_{\epsilon }(k_j^0)\right)}^{-1}{\mathcal{B}}_{\epsilon }(\omega )\Vert <1,$$

$$k_j^{♯}-k_j^0$$

$$\times {\partial }_{\omega }\left({\mathcal{A}}_{\epsilon }(k_j^0)+(\omega -k_j^0){\partial }_{\omega }{\mathcal{A}}_{\epsilon }(k_j^0)\right)\mathrm{d}\omega$$

$$=\frac{1}{2\pi i}\mathrm{tr}{\int }_{\partial V_{\epsilon }}(\omega -k_j^0){\left({\mathcal{A}}_{\epsilon }(k_j^0)+(\omega -k_j^0){\partial }_{\omega }{\mathcal{A}}_{\epsilon }(k_j^0)\right)}^{-1}{\partial }_{\omega }{\mathcal{A}}_{\epsilon }(k_j^0)\mathrm{d}\omega ,$$

$$k_j^{♯}-k_j^0$$

$$\times \left({\mathcal{A}}_{\epsilon }(k_j^0)+(\omega -k_j^0){\partial }_{\omega }{\mathcal{A}}_{\epsilon }(k_j^0)-{\mathcal{A}}_{\epsilon }(k_j^0)\right)\mathrm{d}\omega$$

$$=\frac{1}{2\pi i}\mathrm{tr}\left({\int }_{\partial V_{\epsilon }}\mathrm{I}\mathrm{d}\omega -{\int }_{\partial V_{\epsilon }}{\left({\mathcal{A}}_{\epsilon }(k_j^0)+(\omega -k_j^0){\partial }_{\omega }{\mathcal{A}}_{\epsilon }(k_j^0)\right)}^{-1}{\mathcal{A}}_{\epsilon }(k_j^0)\mathrm{d}\omega \right)$$

$$=-\frac{1}{2\pi i}\mathrm{tr}{\int }_{\partial V_{\epsilon }}{\left({\mathcal{A}}_{\epsilon }(k_j^0)+(\omega -k_j^0){\partial }_{\omega }{\mathcal{A}}_{\epsilon }(k_j^0)\right)}^{-1}{\mathcal{A}}_{\epsilon }(k_j^0)\mathrm{d}\omega .$$

$$k_j^{\epsilon }-k_j^{♯}=-({\mathcal{B}}_{\epsilon }(k_j^{♯}){\psi }_j^{♯},{\psi }_j^{♯}),$$

$$k_j^{\epsilon }-k_j^0=(k_j^{♯}-k_j^0)(1+\mathcal{O}(k_j^{♯}-k_j^0)),$$

$${\mathrm{Z}}_{\mathrm{N}}^k(x_{\text{S}},z)={\mathrm{Z}}_{\mathrm{D}}^k(x_{\text{S}},z)-{\epsilon }^2\frac{\pi }{2}{\partial }_{{\nu }_{y_{\star }}}{\mathrm{Z}}_{\mathrm{D}}^k(z,y_{\star }){\partial }_{{\nu }_{y_{\star }}}{\mathrm{Z}}_{\mathrm{D}}^k(x_{\text{S}},y_{\star })+\mathcal{O}\left(\frac{{\epsilon }^2}{|\log (\epsilon /2){|}^2}\right).$$

$${\mathrm{Z}}_{x_{\text{S}}}^k(y)=\sum _{j=1}^{\infty }\frac{u_j(x_S)u_j(y)}{k^2-{\lambda }_j^{{\Gamma }_{\text{D}}}}.$$

Full text

Mathematics Subject Classification (MSC2000).  35R30, 35C20.

Keywords.  Zaremba eigenvalue problem, boundary integral operators, mixed boundary conditions, metasurfaces.

1 Introduction

This paper develops a new and efficient approach for maximizing the transmission signal between two points at a chosen frequency through changes to specific eigenmodes of the cavity. These changes are achieved by changing the boundary conditions. The eigenmodes and the associated eigenfrequencies of a cavity are sensitively dependant on the geometric properties of the domains, as well as the location of Dirichlet and Neumann boundary conditions. Many recent works have been devoted to the understanding of the effect of changing the boundary condition on the eigenmodes and the eigenfrequencies [1, 2, 4, 5, 11, 13, 14, 16].

Through the use of a tunable reflecting metasurface, the boundary condition can be switched from Dirichlet to Neumann at some specific resonant frequencies [3]. In [3, Part I], the physical mechanism underlying the concept of tunable metasurfaces is modeled both mathematically and numerically. It is shown that an array of coupled pairs of Helmholtz resonators behaves as an equivalent surface with Neumann boundary condition at some specific subwavelength resonant frequencies, where the size of one pair of Helmholtz resonators is much smaller than the wavelengths at the resonant frequencies. The Green’s function of a cavity with mixed (Dirichlet and Neumann) boundary conditions (called also the Zaremba function) is also characterized. In [3, Part II], a one-shot optimization algorithm is proposed and used to obtain a good initial guess for the positions around which the boundary conditions should be switched from Dirichlet to Neumann.

In this paper, we present a new methodology for maximizing the Zaremba function between two points at a chosen frequency through specific eigemodes of the cavity. The paper is organized as follows. In Section 2, we first recall some useful results on the eigenvalues of the mixed boundary value problem (called Zaremba eigenvalue problem). Of particular interest is their monotonicity property with respect to the size of the Neumann part proven in [14]. Then we reformulate the eigenvalue problem using boundary integral operators. Based on this nonlinear formulation and the use of the generalized argument principle for the characterization of the characteristic values of finitely meromorphic operator-valued functions of Fredholm-type, we derive an accurate asymptotic formula of the changes of eigenfrequencies of a cavity with mixed boundary conditions in terms of the size of the part of the cavity boundary where the boundary condition is switched from Dirichlet to Neumann. Finally, we recall the asymptotic expansion of the Zaremba function in terms of the size of the Neumann part. The problem of changing a portion of a Dirichlet boundary to Neumann is more delicate than the converse. If a portion of the boundary is changed from having Neumann conditions to Dirichlet, the reverse consideration than in this paper, then an asymptotic expansion of the eigenvalues is easier to derive [5, 16]. The perturbation theory for the introduction of Neumann boundaries requires a careful consideration of the asymptotic behaviour of the Zaremba near the perturbation [3, Part II]. In Section 3, we derive a spectral decomposition of the Zaremba function. In Section 4, we consider the problem where we have a source in a bounded domain operating at a given frequency, and we want to determine, by exploiting the monotonicity property of the eigenvalues of the mixed boundary value problem, which part of the boundary to choose to be reflecting such that an eigenvalue of the mixed boundary value problem gets close enough to the operating frequency. In order to significantly enhance the signal at a given receiving point, both the emitter and the receiver should not belong to the nodal set corresponding to the eigenmode associated with the eigenvalue of the mixed boundary value problem.

There are two distinct issues: where to place the Neumann boundary condition, and how long it should be, to achieve the twin objectives of maximizing gain between a fixed source-receiver pair, and at a frequency close to a desired one.

Our main idea is to first nucleate the Neumann boundary conditions in order to maximize gain of the Zaremba function by making use of an asymptotic expansion of the Zaremba function in terms of the size of the Neumann part. Then the size of the Neumann part is changed in such away that an eigenvalue of the mixed boundary value problem gets close to the operating frequency by using the monotonicity property of the eigenvalues of the mixed eigenvalue problem. The optimization needs the high-accuracy evaluation of certain boundary integral operators, and this is done using techniques from [1, 2].

We present in Section 5 some numerical experiments to show the applicability and the accuracy of the proposed methodology.

2 Preliminaries

2.1 Laplace Eigenvalue with Mixed Boundary Conditions

Let $\Omega\subset\mathbb{R}^{2}$ be an open, bounded domain with a smooth boundary. We define $\overline{\Omega}$ as the topological closure of $\Omega$. We decompose the boundary $\partial\Omega\coloneqq\overline{\Omega}\setminus\Omega$ into two parts, \partial\Omega=\overline{{\Gamma_{\textrm{D}}}\mathbin{\mathchoice{\leavevmode\vtop{\halign{\hfil\m@th\displaystyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\textstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptscriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}}{\Gamma_{\textrm{N}}}}, where ${\Gamma_{\textrm{D}}}$ and ${\Gamma_{\textrm{N}}}$ are finite unions of open boundary sets. We define $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$ to be a partition of $\partial\Omega$. Let $x_{\textrm{S}}\in\Omega$ and $k\in(0,\infty)$. The Zaremba function $\mathrm{Z}^{k}_{x_{\textrm{S}}}(x_{s},\cdot):\Omega\setminus\{x_{\textrm{S}}\}\rightarrow\mathbb{R}$ is the Green’s function to the Zaremba problem, also known as the fundamental Helmholtz equation with mixed boundary conditions,

[TABLE]

Here $\nu_{y}$ denotes the outer normal at $y\in\partial\Omega$ and $\partial_{\nu_{y}}$ the normal derivative at $y\in\partial\Omega$. It is clear that we can write

[TABLE]

where $\Gamma^{k}(x,y)\coloneqq\frac{i}{4}H_{0}^{1}(k\,|x-y|)$ is the fundamental solution of the Helmholtz problem with wavenumber $k$, and $\mathrm{R}^{k}(x_{s},\cdot)$ is a smooth function satisfying the boundary value problem

[TABLE]

In Section 3, we will see that $\mathrm{Z}^{k}_{x_{\textrm{S}}}$ exists for all but countably many values of $k$, which are related to the unique solvability of the problem for $\mathrm{R}^{k}(x_{\textrm{S}},\cdot)$. These exceptional values of $k$ are the eigenvalues to the associated Laplace eigenvalue problem with mixed boundary conditions

[TABLE]

Equation (2.3) has a non-trivial solution $u\in\mathrm{H}^{1}(\Omega)$ for a countable set of real values of $k^{2}$ [15, Theorem 4.10], which we refer to as $\{\lambda^{{\Gamma_{\textrm{D}}}}_{j}\}_{j=1}^{\infty}$, so that $\lambda^{{\Gamma_{\textrm{D}}}}_{1}\leq\lambda^{{\Gamma_{\textrm{D}}}}_{2}\leq\lambda^{{\Gamma_{\textrm{D}}}}_{3}\leq\ldots\,$. We know that $\lambda^{{\Gamma_{\textrm{D}}}}_{1}\geq 0$ and that $\lim_{j\rightarrow\infty}\lambda^{{\Gamma_{\textrm{D}}}}_{j}=\infty$ for all partitions $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$ of $\partial\Omega$.

We denote by $\{\lambda^{\partial\Omega}_{j}\}_{j\in\mathbb{N}}$ the pure Dirichlet eigenvalues for $\Omega$, corresponding to the case ${\Gamma_{\textrm{D}}}=\partial\Omega$. We let $\{\lambda^{\varnothing}_{j}\}_{j\in\mathbb{N}}$ denote the Neumann eigenvalues associated to the case ${\Gamma_{\textrm{N}}}=\partial\Omega$. Then we have

[TABLE]

In [9], it is shown that $\lambda^{\varnothing}_{j+1}<\lambda^{\partial\Omega}_{j}$, for all $j\in\mathbb{N}$, for a very general class of domains $\Omega$.

Remark 2.1**.**

Let $\Omega$ be the unit circle, we have that $\{\lambda^{\partial\Omega}_{j}\}_{j=1}^{\infty}$ is (up to sorting) equal to

[TABLE]

where $J_{n}$ is a Bessel function of the first kind and order $n$. The eigenvalues corresponding to the roots of $J_{0}$ appear as simple Dirichlet eigenvalues; all others have multiplicity two. and $\{\lambda^{\varnothing}_{j}\}_{j=2}^{\infty}$ is (up to sorting) equal to

[TABLE]

*Again, the eigenvalues corresponding to the roots of $J^{\prime}_{0}$ appear as simple Neumann eigenvalues; all others have multiplicity two. We refer to [10]. *

Recently, Lotoreichik and Rohleder [14, Proposition 2.3] showed the following monotonicity statement.

Proposition 2.2**.**

Let $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}}),({\Gamma_{\textrm{D}}}^{\prime},{\Gamma_{\textrm{N}}}^{\prime})$ be two partitions of $\partial\Omega$, such that ${\Gamma_{\textrm{D}}}\subset{\Gamma_{\textrm{D}}}^{\prime}$. If ${\Gamma_{\textrm{D}}}^{\prime}\setminus{\Gamma_{\textrm{D}}}$ has a non-empty interior then

[TABLE]

With Proposition 2.2, we can readily infer that if $\varnothing\neq{\Gamma_{\textrm{D}}}$, and $\overline{{\Gamma_{\textrm{D}}}}\neq\partial\Omega$, then

[TABLE]

2.2 Boundary Integral Formulation of the Eigenvalue Problem

The solution $u$ of the eigenvalue(2.3) can be represented by a single layer potential

[TABLE]

with surface density $\psi\in L^{2}(\partial\Omega)$.

We define then the operators $\mathcal{S}_{\Gamma_{\textrm{D}}}^{k}:H^{-1/2}({\Gamma_{\textrm{D}}})\to H^{1/2}({\Gamma_{\textrm{D}}})$, $\mathcal{S}_{\Gamma_{\textrm{N}}}^{k}:H^{-1/2}({\Gamma_{\textrm{N}}})\to H^{-1/2}({\Gamma_{\textrm{D}}})$, $(\mathcal{K}^{k}_{\Gamma_{\textrm{N}}})^{\ast}:H^{-1/2}({\Gamma_{\textrm{N}}})\to H^{-1/2}({\Gamma_{\textrm{N}}})$ and $\partial\mathcal{S}^{k}_{{\Gamma_{\textrm{D}}}}:H^{-1/2}({\Gamma_{\textrm{D}}})\to H^{-1/2}({\Gamma_{\textrm{N}}})$ by

[TABLE]

where the ’p.v.’ stands for the principle value integral. This actually is the standard (Lebesgue-) integral for a smooth curved ${\Gamma_{\textrm{N}}}$, since $\partial_{\nu}\Gamma^{k}$ is a bounded and sufficiently smooth integral operator kernel. From [17, Chapter 11] we have that $\mathcal{S}_{\Gamma_{\textrm{D}}}^{k}$ is a Fredholm operator with index 0, we also readily infer that $(\mathcal{K}^{k}_{\Gamma_{\textrm{N}}})^{\ast}$, $\partial\mathcal{S}^{k}_{{\Gamma_{\textrm{D}}}}$, and $\mathcal{S}_{\Gamma_{\textrm{N}}}^{k}$ are compact operator.

We then define $\mathcal{A}(k):H^{-1/2}({\Gamma_{\textrm{D}}})\times H^{-1/2}({\Gamma_{\textrm{N}}})\to H^{1/2}({\Gamma_{\textrm{D}}})\times H^{-1/2}({\Gamma_{\textrm{N}}})$ in terms of these integral operators through

[TABLE]

We readily see that $\mathcal{A}(k)$ is an analytic Fredholm operator of index [math] in $\mathbb{C}\setminus i\mathbb{R}^{-}$.

To locate the Zaremba eigenvalues, we have the following statement:

[TABLE]

In [2, Section 3] and [1], it is shown that every square root of a Zaremba eigenvalue is a real positive characteristic value of $k\mapsto\mathcal{A}(k)$ and every real positive characteristic value of $k\mapsto\mathcal{A}(k)$ is the square root of a Zaremba eigenvalue.

We see that $\mathcal{A}(k)$ is invertible for $k\in(0,\infty)$ not a square root of a Zaremba eigenvalue.

We remark that the non-real characteristic values of $k\mapsto\mathcal{A}(k)$ cannot correspond to eigenvalues to the Laplace equation. This yields the undesirable, but avoidable, difficulty in choosing a neighbourhood $V$ to apply Proposition 2.5 in our algorithm, see also Section 4, comment on Line 13.

The Statement (2.2) allows for a discretization and thus a numerical approximation of the value $k$. We will use this further on. For these facts, we refer to [2, Sections 3 and 5].

Let us also consider the regularity of the solution $u$ and the density $\psi$ near a Dirichlet-Neumann junction. The following result can be found in [1, Theorems 4.2 and 4.3].

Proposition 2.3**.**

Let ${\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}}$ be non-empty. Let $k>0$ and $\psi$ satisfy the Statement (2.2). Let $y_{\star}\in\overline{{\Gamma_{\textrm{D}}}}\cap\overline{{\Gamma_{\textrm{N}}}}$. Then there exists a neighborhood $\mathcal{U}\subset\mathbb{R}^{2}$ around $y_{\star}$ such that for all $y\in\mathcal{U}$ and for all $n\in\mathbb{N}$

[TABLE]

*where $z\in\mathbb{C}$ is the complexification of $y-y_{\star}$, that is $z=\left(y_{1}-(y_{\star})_{1}\right)+i\left(y_{2}-(y_{\star})_{2}\right)$ with $i$ being the imaginary unit, and $\overline{z}$ being its conjugate value, and where $P^{n}_{y_{\star}},\mathrm{Q}_{\mathrm{D},y_{\star}}^{n},\mathrm{Q}_{\mathrm{N},y_{\star}}^{n}$ are polynomial functions of their respective arguments and of a degree such that none of their terms can be included in their respective error terms. *

2.3 Approximation of the Zaremba Eigenvalue using the Generalized Argument Principle

In this section we derive asymptotic expressions for the perturbation of the Zaremba eigenvalues, when a small portion of the boundary is changed from Dirichlet to Neumann.

Let ${\Gamma_{\Delta}}\subset\partial\Omega$ be a boundary interval of length $2\varepsilon$. Let ({\Gamma_{\textrm{D}}}\mathbin{\mathchoice{\leavevmode\vtop{\halign{\hfil\m@th\displaystyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\textstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptscriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}}{\Gamma_{\Delta}},{\Gamma_{\textrm{N}}}) be a partition of $\partial\Omega$. We associate the operator $\mathcal{A}_{0}(k)$, defined via (2.5), to that partition. This corresponds to ${\Gamma_{\Delta}}$ having a Dirichlet boundary condition. Then we define $\mathcal{A}_{\varepsilon}(k)$, also by obvious changes in the integrals in (2.5), to be the operator associated to the partition ({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}}\mathbin{\mathchoice{\leavevmode\vtop{\halign{\hfil\m@th\displaystyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\textstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptscriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}}{\Gamma_{\Delta}}). This in turn corresponds to ${\Gamma_{\Delta}}$ being a Neumann part. For ease of notation, we define k_{j}^{0}\coloneqq\sqrt{\lambda^{{\Gamma_{\textrm{D}}}\mathbin{\mathchoice{\leavevmode\vtop{\halign{\hfil\m@th\displaystyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\textstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptscriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}}{\Gamma_{\Delta}}}_{j}} and $k_{j}^{\varepsilon}\coloneqq\sqrt{\lambda^{{\Gamma_{\textrm{D}}}}_{j}}$ for all $j\in\mathbb{N}$, and call those characteristic values to their respective operators. From [6, Lemma 3.8] we then have the following Lemma:

Lemma 2.4**.**

Let $k_{j}^{0}$ be a simple characteristic value. Let $V\subset\mathbb{C}$ be a neighbourhood of $k_{j}^{0}$, such that $k_{j}^{\varepsilon}\in V$. Assume further that no other square root of Zaremba eigenvalue to the partition ({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}}\mathbin{\mathchoice{\leavevmode\vtop{\halign{\hfil\m@th\displaystyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\textstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptscriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}}{\Gamma_{\Delta}}) of $\partial\Omega$ is in $\overline{V}$. Then $k_{j}^{\varepsilon}$ is given by the contour integral

[TABLE]

Here $\partial_{\omega}$ denotes the variation of the operator in the wavenumber parameter $\omega$. This expression is exact. Unfortunately, its use in a practical algorithm is limited, since it would entail inverting the operator $\mathcal{A}_{\varepsilon}(\omega)$ for each $\varepsilon$ used in an optimization. It is useful, therefore, to locate an expression in which this inverse is approximated by $\mathcal{A}_{0}(\omega)$ instead.

From [6, Theorem 3.12] we get the approximation

[TABLE]

where we expect the error to be in $o\left(\frac{1}{|\log(\varepsilon)|}\right)$. We can, in fact, obtain a faster and even more accurate approximation, which we describe in the following proposition.

Proposition 2.5**.**

Let $k_{j}^{0}$ be the $jth$ (sorted) characteristic value of $\mathcal{A}_{0}(k)$ corresponding to the decomposition $\Gamma_{\textrm{D}},\Gamma_{\textrm{N}},$ and assume it is simple Then one can find a $\varepsilon>0$ and a neighbourhood $V\subset\mathbb{C}$ containing $k_{j}^{0}$ so that

• •

the $jth$ characteristic value $k_{j}^{\varepsilon}$ of the operator $\mathcal{A}_{\varepsilon}(k)$ (obtained by changing $\Gamma_{\Delta}$ to a Neumann boundary condition) is contained $\in V$;

• •

no other square root of the Laplace eigenvalues to the partition ({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}}\mathbin{\mathchoice{\leavevmode\vtop{\halign{\hfil\m@th\displaystyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\textstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptscriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}}{\Gamma_{\Delta}}) of $\partial\Omega$ are in $\overline{V}$.

• •

The characteristic value of the perturbed operator $k_{j}^{\varepsilon}$ is given by

[TABLE]

*Here $\mathrm{I}$ is the identity operator. *

Proof**.**

We first observe from Proposition 2.2 together with the fact that $\mathcal{A}_{\varepsilon}(k)$ is a Fredholm analytic operator of index [math] in $\mathbb{C}\setminus i\mathbb{R}^{-}$, we can see that $k_{j}^{\varepsilon}\nearrow k_{j}^{0}$ for $\varepsilon\searrow 0^{+}$. We now examine the perturbed operator $\mathcal{A}_{\varepsilon}$. Its characteristic value is $k_{j}^{\varepsilon}$. Provided $k_{j}^{0}$ is sufficiently close to $k_{j}^{\varepsilon}$, we have the following Taylor expansion:

[TABLE]

where $\mathcal{B}_{\varepsilon}(\omega)=\mathcal{O}((\omega-k_{j}^{0})^{2})$. This expansion holds only in a neighborhood $V^{0}_{\varepsilon}$ of $k_{j}^{0}$, and so $\varepsilon$ must be small enough such that $k_{j}^{0}\in V^{0}_{\varepsilon}$.

Then consider that we have in the operator-norm

[TABLE]

for $\omega\in V_{\varepsilon}\subset V^{0}_{\varepsilon}$ close enough to both $k_{j}^{\varepsilon}$ and $k_{j}^{0}$, because then the Taylor remainder $\mathcal{B}_{\varepsilon}(\omega)=\mathcal{O}((\omega-k_{j}^{0})^{2})$. If $\varepsilon$ is small enough, then $k_{j}^{\varepsilon}\in V_{\varepsilon}$. Then by the Generalization of Rouché’s Theorem [6, Theorem 1.15] we have that since $\mathcal{A}_{\varepsilon}(k_{j}^{0})+(\omega-k_{j}^{0})\partial_{\omega}\mathcal{A}_{\varepsilon}(k_{j}^{0})$ and $\mathcal{A}_{\varepsilon}(\omega)$ are close in operator norm, they both have the same number of characteristic values in $V_{\varepsilon}$. Thus $\mathcal{A}_{\varepsilon}(k_{j}^{0})+(\omega-k_{j}^{0})\partial_{\omega}\mathcal{A}_{\varepsilon}(k_{j}^{0})$ has a simple characteristic value $k^{\sharp}_{j}$ in $V_{\varepsilon}$. Now we can use Lemma 2.4, but replacing $\mathcal{A}_{0}(\omega)$ by $\left(\mathcal{A}_{\varepsilon}(k_{j}^{0})+(\omega-k_{j}^{0})\partial_{\omega}\mathcal{A}_{\varepsilon}(k_{j}^{0})\right)$:

to get

[TABLE]

and hence,

[TABLE]

Moreover, by a standard perturbation argument [6, Section 5.2.4], we have at the leading-order term

[TABLE]

where $\psi^{\sharp}_{j}$ is the root function associated with the characteristic value $k^{\sharp}_{j}$ evaluated at $k^{\sharp}_{j}$. Thus,

[TABLE]

and therefore, Proposition 2.5 holds.

$\blacksquare$

We remark on the significance of this result from the point of view of computation, and which makes it a key ingredient in our algorithm. If one seeks a high-accuracy approximation of the characteristic value $k_{j}^{\varepsilon}$ of $\mathcal{A}_{\varepsilon}(k)$, and one already has a good approximation of $k_{j}^{0}$, the approximation in Proposition 2.5 allows us to proceed by assembling only one matrix, that corresponding to $\mathcal{A}_{\varepsilon}(k_{j}^{0})$. The contour integrals can be effectively computed using the trapezoidal rule, making this an inexpensive but very accurate approximation of $k_{j}^{\varepsilon}$.

2.4 Approximation of the Zaremba Function

Let $\Omega\coloneqq\{z\in\mathbb{C}\mid|z|<1\}$, and let ${\Gamma_{\Delta}}\subset\partial\Omega$ be a boundary interval of length $2\varepsilon$ with center $y_{\star}\in{\Gamma_{\Delta}}$. Let $(\partial\Omega,\varnothing)$ be the partition of $\partial\Omega$, and with it we associate the Zaremba function $\mathrm{Z}_{\mathrm{D}}^{k}(x_{\textrm{S}},\cdot)$, for $x_{\textrm{S}}\in\Omega$, defined via (2.1). This corresponds to ${\Gamma_{\Delta}}$ having a Dirichlet boundary condition. Then we define $\mathrm{Z}_{\mathrm{N}}^{k}(x_{\textrm{S}},\cdot)\in\mathrm{L}^{2}(\Omega)$, $x_{\textrm{S}}\in\Omega$, also defined via (2.1), to be the Zaremba function associated to the partition $(\partial\Omega\setminus\overline{{\Gamma_{\Delta}}},{\Gamma_{\Delta}})$. This in turn corresponds to ${\Gamma_{\Delta}}$ having a Neumann boundary condition. We then have the following lemma.

Lemma 2.6**.**

Let $\Omega,y_{\star},{\Gamma_{\Delta}},\mathrm{Z}_{\mathrm{D}}^{k}(x_{\textrm{S}},\cdot)$ and $\mathrm{Z}_{\mathrm{N}}^{k}(x_{\textrm{S}},\cdot)$ be defined as described above. Let $\varepsilon>0$ be small enough. Let $k>0$, such that $k^{2}\neq\lambda_{j}^{\partial\Omega\setminus{\Gamma_{\Delta}}}$, and $k^{2}\neq\lambda_{j}^{\partial\Omega}$, for all $j\in\mathbb{N}$. Then for all $z\in\Omega$,

[TABLE]

Lemma 2.6 follows readily from combining the results in [3, Theorem 5.4] and [3, Equation (6.24)]

Numerical experiments confirm that $|\mathrm{Z}_{\mathrm{N},x_{\textrm{S}}}^{k}(y)-\mathrm{Z}_{\mathrm{D},x_{\textrm{S}}}^{k}(y)|$ is of order of $\varepsilon^{2}$, as long as $y$ is far enough away from the boundary.

3 Spectral Decomposition of the Zaremba Function

Let us again consider the more general setup at the beginning of Section 2, that is let $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$ be a partition of $\partial\Omega$, let $\{\lambda^{{\Gamma_{\textrm{D}}}}_{j}\}_{j=1}^{\infty}$ be the Zaremba eigenvalues and let $\{u_{j}\}_{j=1}^{\infty}$ be an $L^{2}$-orthonormal basis of associated eigenfunctions. Then we have the following statement about the Zaremba function $\mathrm{Z}^{k}_{x_{\textrm{S}}}$, $x_{\textrm{S}}\in\Omega$, defined by (2.1).

Theorem 3.1**.**

For all $y\in\Omega$, $y\neq x_{S}$ and for all $k>0$ which are not in the spectrum, ie, $k^{2}\neq\lambda^{{\Gamma_{\textrm{D}}}}_{j}$ of the Zaremba eigenvalue problem, the Zaremba function $\mathrm{Z}^{k}_{x_{\textrm{S}}}$, given by (2.1), exists and is in $\mathrm{L}^{2}_{\text{loc}}(\Omega)$. Furthermore, we can write it as

[TABLE]

Next, we will consider the proof of Theorem 3.1. To this end, we define

[TABLE]

Consider that the solution to the Laplace eigenvalue-equation $u$ is element of $\mathrm{H}_{0,{\Gamma_{\textrm{D}}}}^{1}(\Omega)$.

[TABLE]

The operator $-\mathop{}\!\mathbin{\bigtriangleup}$ is selfadjoint in $\mathrm{L}^{2}(\Omega)$, which we readily see using Green’s identity, and it has thus a discrete spectrum. Moreover, $-\mathop{}\!\mathbin{\bigtriangleup}$ corresponds to the sesquilinear form $\langle v_{1},v_{2}\rangle\mapsto(\nabla v_{1},\nabla v_{2})_{\mathrm{L}^{2}(\Omega)}$ with domain $\mathrm{H}_{0,{\Gamma_{\textrm{D}}}}^{1}$, since $(-\mathop{}\!\mathbin{\bigtriangleup}w_{1},w_{2})_{\mathrm{L}^{2}(\Omega)}=(\nabla w_{1},\nabla w_{2})_{\mathrm{L}^{2}(\Omega)}$ for all $w_{1},w_{2}\in\mathrm{dom}(-\mathop{}\!\mathbin{\bigtriangleup})$, see [7, 12, 18] for more details on semi-bounded self-adjoint operators and corresponding quadratic forms. And the form $\langle\cdot\,,\cdot\rangle$ is closed, non-negative and symmetric. This allows us to use the min-max principle. Thus we can write for all $j\in\mathbb{N}$,

[TABLE]

This leads us to the following lemma.

Lemma 3.2**.**

For all $f\in\mathrm{dom}(-\mathop{}\!\mathbin{\bigtriangleup})$, we have that

[TABLE]

*where $c_{j}\coloneqq(f\,,u_{j})_{\mathrm{L}^{2}(\Omega)}$, that is the linear subset spanned by eigenfunctions of the Laplace eigenvalue-equation with mixed boundary conditions (2.3) is dense in $\mathrm{dom}(-\mathop{}\!\mathbin{\bigtriangleup})$. *

Proof**.**

Let $\mathrm{r}_{N}\coloneqq f-\sum_{j=1}^{N}c_{j}\,u_{j}$. Then for all $i=1,\ldots,N$, we have that

[TABLE]

where we used Green’s identity and the fact that $f,u_{j}\in\mathrm{dom}(-\mathop{}\!\mathbin{\bigtriangleup})$. Next, we want to show that

[TABLE]

To this end, consider the min-max principle (3.1), it tells us that

[TABLE]

Thus, we can infer $\lambda^{{\Gamma_{\textrm{D}}}}_{N}\leq\frac{\left\lVert\nabla\mathrm{r}_{N}\right\rVert^{2}}{\left\lVert\mathrm{r}_{N}\right\rVert^{2}}$ from $\lambda^{{\Gamma_{\textrm{D}}}}_{n-1}\leq\frac{\left\lVert\nabla\mathrm{r}_{N}\right\rVert^{2}}{\left\lVert\mathrm{r}_{N}\right\rVert^{2}}$, which in turn is given by an induction argument, whose induction basis follows trivially from the min-max principle (3.1). Using the definition of $c_{j}$, we have that

[TABLE]

Thus, using (3.3), we have that

[TABLE]

Since $\left\lVert\nabla f\right\rVert^{2}_{\mathrm{L}^{2}}=(f,-\mathop{}\!\mathbin{\bigtriangleup}f)_{\mathrm{L}^{2}}\leq\left\lVert f\right\rVert_{\mathrm{L}^{2}}\left\lVert\mathop{}\!\mathbin{\bigtriangleup}f\right\rVert_{\mathrm{L}^{2}}<\infty$, $\left\lVert\nabla f\right\rVert^{2}_{\mathrm{L}^{2}}$ is bounded. Using the fact that $\lambda^{{\Gamma_{\textrm{D}}}}_{N}\xrightarrow{N\rightarrow\infty}\infty$, we have that $\left\lVert\mathrm{r}_{N}\right\rVert_{\mathrm{L}^{2}}^{2}\xrightarrow{N\rightarrow\infty}0$. This completes the proof of Lemma 3.2. $\blacksquare$

Proof** (Theorem 3.1).**

To show the existence of the Zaremba function $\mathrm{Z}^{k}_{x_{\textrm{S}}}$, we write $\mathrm{Z}^{k}_{x_{\textrm{S}}}(y)$, for all $y\in\Omega\,,y\neq x_{\textrm{S}}$ as

[TABLE]

where $\Gamma^{k}$ is the fundamental solution to the Helmholtz equation, and $\mathrm{R}^{k}$ satisfies

[TABLE]

The solution to (3.6) does exist, for those values of $k$ specified in the theorem, and it is in $\mathrm{H}^{1}(\Omega)$, see [15, Theorem 4.10]. Using that $\Gamma^{k}(x_{\textrm{S}},\cdot)\in L^{2}(\Omega)$, we have that $\mathrm{Z}^{k}_{x_{\textrm{S}}}(y)\in L^{2}(\Omega)$. Thus from Lemma 3.2 and the density of $\mathrm{dom}(-\mathop{}\!\mathbin{\bigtriangleup})$ in $L^{2}(\Omega)$, we have that for all $y\in\Omega$, $y\neq x_{\textrm{S}}$,

[TABLE]

for some $a_{j}\in\mathbb{R}$, depending on $x_{\textrm{S}}$. Let us give an expression for the $a_{j}$. Using Green’s identity, we have that

[TABLE]

where we used Fubini’s theorem to interchange summation and integration. With that we infer that for all $i\in\mathbb{N}$,

[TABLE]

and this concludes the proof. $\blacksquare$

4 The Algorithm

We next present our main algorithm for wave enhancement. We begin with a domain $\Omega$, the source point $x_{\textrm{S}}$ and the receiver point $y$, both in $\Omega$, and a predetermined target value ${k_{\star}}$ corresponding to a desired transmission frequency.

First, we determine the next higher Dirichlet eigenvalue to $k_{\star}^{2}$, which is done using a discretized version of the operator $\mathcal{A}(k)$ given in Section 2. The discretization follows the procedure developed in [2].

Second, we determine a location $y_{\star}$ on the boundary $\partial\Omega$, which yields a higher absolute value of $|\mathrm{Z}^{k}_{x_{\textrm{S}}}(x_{\textrm{S}},y)|$, when we insert a small enough Neumann boundary at that location. Finding the location is established using Lemma 2.6, that is we find the local maxima or minima of

[TABLE]

The computation of the Zaremba function is done by solving the problem 2.2 using the procedure described in [1], also uses the operator $\mathcal{A}(k)$.

Third, we successively increase the Neumann boundary until the characteristic value hits the target characteristic value. The computation of the new characteristic value after a small increase of the Neumann boundary is achieved using Proposition 2.5. It might be that we need to increase the boundary initially by a large amount, and the resulting characteristic value has to be computed with the time-expensive procedure described in [2].

A more detailed explanation is given in the comments after Algorithm 4. We note here that changing a boundary part from the Dirichlet boundary condition to the Neumann one, the associated Laplace eigenvalue $\lambda_{j}^{{\Gamma_{\textrm{D}}}}$ decreases, according to Proposition 2.2, and thus the characteristic value $\sqrt{\lambda_{j}^{{\Gamma_{\textrm{D}}}}}$ decreases as well. Moreover, $\lambda_{j}^{{\Gamma_{\textrm{D}}}}$ is between the Neumann and the Dirichlet eigenvalue, that is $\lambda_{j}^{\varnothing}\leq\lambda_{j}^{{\Gamma_{\textrm{D}}}}\leq\lambda_{j}^{\partial\Omega}$. Increasing boundary length enough, we eventually hit the target characteristic value $k_{\star}$, because \cup_{j=1}^{\infty}\Big{(}\lambda_{j}^{\varnothing},\lambda_{j}^{\partial\Omega}\Big{)}=(0,\infty), since $\lambda_{j+1}^{\varnothing}<\lambda_{j}^{\partial\Omega}$, proved in [9].

Algorithm 1 Finding an intensity maximizing partition of the boundary

Input: $\varepsilon>0$, $x_{\textrm{S}}\in\Omega$, $y\in\Omega$, $y\neq x_{\textrm{S}}$, $k_{\star}>0$, $C_{\mathrm{tol}}>0$.

Require: $\varepsilon$ is small enough, $C_{\mathrm{tol}}$ is big enough.

1:Let ${\Gamma_{\textrm{D}}}\coloneqq\partial\Omega$, ${\Gamma_{\textrm{N}}}\coloneqq\varnothing$.

2:Find the next higher square root of the Dirichlet eigenvalue $k$ to $k_{\star}$.

3:Compute the value $\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},y)$ and the normal derivative of the Zaremba functions $\partial_{\nu_{\cdot}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},\cdot),\partial_{\nu_{\cdot}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(y,\cdot)$ associated to the partition $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$ at the boundary.

4:if $\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},y)\geq 0$ then

5:     Let $\mathsf{S}$ be the location of a global minima of the function \partial\Omega\ni z\mapsto\big{(}\partial_{\nu_{z}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},z)\cdot\partial_{\nu_{z}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(y,z)\big{)}\in\mathbb{R}.

6:else if $\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},y)<0$ then

7:     Let $\mathsf{S}$ be the location of a global maxima of the function \partial\Omega\ni z\mapsto\big{(}\partial_{\nu_{z}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},z)\cdot\partial_{\nu_{z}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(y,z)\big{)}\in\mathbb{R}.

8:end if

9:$({\Gamma_{\textrm{D}}}^{0},{\Gamma_{\textrm{N}}}^{0})\coloneqq({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$.

10:while True do

11:     Define ${\Gamma_{\Delta}}$ to be a boundary interval of length $2\varepsilon$ with center $\mathsf{S}$.

12:     ({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})\coloneqq({\Gamma_{\textrm{D}}}\setminus\overline{{\Gamma_{\Delta}}},\;{\Gamma_{\textrm{N}}}\mathbin{\mathchoice{\leavevmode\vtop{\halign{\hfil\m@th\displaystyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\textstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}{\leavevmode\vtop{\halign{\hfil\m@th\scriptscriptstyle#\hfil\cr\cup\cr\cdot\crcr}}}}{\Gamma_{\Delta}}).

13:     Compute the perturbed characteristic value $k$ associated to the partition $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$ as described in Section 2.3 or with the procedure given in [2].

14:     if $|k-k_{\star}|\leq C_{\mathrm{tol}}$ then return $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$

15:     else if $k_{\star}+C_{\mathrm{tol}}<k$ then

16:         $\mathrm{BREAK}\;\mathrm{WHILE}$

17:     else

18:         $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})\coloneqq({\Gamma_{\textrm{D}}}^{0},{\Gamma_{\textrm{N}}}^{0})$

19:         $\varepsilon\coloneqq\frac{\varepsilon}{\sqrt{2}}$

20:     end if

21:end while

22:$({\Gamma_{\textrm{D}}}^{0},{\Gamma_{\textrm{N}}}^{0})\coloneqq({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$

23:while True do

24:     Define ${\Gamma_{\Delta}}$ to be the extension of the Neumann interval boundary with center $j$, extended on both sides by $\varepsilon/2$.

25:     $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})\coloneqq({\Gamma_{\textrm{D}}}\setminus\overline{{\Gamma_{\Delta}}},\;{\Gamma_{\textrm{N}}}\cup{\Gamma_{\Delta}})$

26:     Compute the perturbed characteristic value $k$ associated to the partition $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$ as described in Section 2.3.

27:     if $|k-k_{\star}|\leq C_{\mathrm{tol}}$ then return $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$

28:     else if $k_{\star}+C_{\mathrm{tol}}<k$ then

29:         $({\Gamma_{\textrm{D}}}^{0},{\Gamma_{\textrm{N}}}^{0})\coloneqq({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})$

30:     else

31:         $({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})\coloneqq({\Gamma_{\textrm{D}}}^{0},{\Gamma_{\textrm{N}}}^{0})$

32:         $\varepsilon\coloneqq{\varepsilon}\cdot 0.9$

33:     end if

34:end while

In the following we give an explanation for the choices.

• Line 2:

The reason we search for the next higher Dirichlet eigenvalue originates from the fact that, according to Proposition 2.2, when we insert Neumann boundaries, the corresponding eigenvalue decreases. The search for the next higher Dirichlet characteristic value an its multiplicity might be computationally expensive.

• Line 3:

Using the algorithm proposed in [1], we compute the Zaremba function using the decomposition $\mathrm{Z}^{k}(x_{\textrm{S}},y)=\Gamma^{k}(x_{\textrm{S}},y)+\mathrm{R}^{k}(x_{\textrm{S}},y)$, where $\Gamma^{k}$ is the fundamental solution to the Helmholtz equation, and $\mathrm{R}^{k}$ satisfies the partial differential equation (3.6). More exactly, we obtain a function $\varphi_{\mathrm{R}}$ on $\partial\Omega$, which is of the form in Proposition 2.3, with

[TABLE]

for $y\in\Omega$. Using the jump relations, see [6, Section 2.3.2], we get for $y\rightarrow{\Gamma_{\textrm{D}}}$ that

[TABLE]

where $\mathrm{I}_{\partial\Omega}$ denotes the identity operator.

Using a discretization to the operator $(\mathcal{K}^{k}_{\partial\Omega})^{\ast}$, which we also readily obtain from [2], we can calculate $\partial_{\nu_{y}}\mathrm{Z}^{k}(x_{\textrm{S}},\cdot)=\partial_{\nu_{y}}\Gamma^{k}(x_{\textrm{S}},\cdot)+\partial_{\nu_{y}}\mathrm{R}^{k}(x_{\textrm{S}},\cdot)$.

• Line 4-8:

In view of Lemma 2.6, we obtain that if $\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},y)\geq 0$ then we need a negative value of $\partial_{\nu_{z}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},z)\cdot\partial_{\nu_{z}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(y,z)$ to increase $\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},y)$ and vice-verca for $\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},y)\leq 0$. Taking the minima, respectively the maxima, we increase the absolute value of $\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},y)$.

We note that Lemma 2.6 only holds for the case where $\Omega$ is the unit circle, but we assume that it holds for all domains with smooth boundaries. We think, this can be established expanding the operator in [3, Theorem 5.4].

From Theorem 3.1 we know that the Zaremba function is real valued, but due to numerical cancellation errors, the Zaremba function might have a non-zero imaginary part.

In our numerical experiments, it always holds that a global minima is negative and a global maxima is positive, respectively. But we do not know if this holds true in general.

• Line 10:

In this while-loop we change a boundary interval with center $\mathsf{S}$ and length $2\varepsilon$ into a Neumann Boundary condition. Then we compute an approximation $k$ to the new characteristic value. If $|k-k_{\star}|<C_{\mathrm{tol}}$, we end the algorithm, if $k+C_{\mathrm{tol}}<k_{\star}$, we break the while-loop, and in the remaining case we decrease $\varepsilon$ and go through the loop again.

• Line 13:

To compute an approximation to the new characteristic value, which is smaller than $k$, we use the approximation stated in Proposition 2.5. To this end, we use as the complex domain $V$ encircling $k$ and $k_{\star}$ an ellipse with center $(k+k_{\star})/2$ and semi-major axis $(k-k_{\star})\cdot 0.55$ and semi-minor axis $(k-k_{\star})\cdot 0.1$, which is to avoid complex characteristic values. Those factors are chosen due to good numerical results. A discretization to the operator $\mathcal{A}(k)$ is computed using the algorithm described in [2]. For the complex derivative of $\mathcal{A}(k)$, we used the rough approximation $(\mathcal{A}(w+0.01)-\mathcal{A}(w))/0.01$. The integral is approximated with a inbuilt-process. The approximation may yield the same result as the former characteristic value, that is $k$. In that case, the new characteristic value is not within $V$, which happens when the new boundary interval with Neumann boundary conditions is too long, or cannot be detected by the approximation.

Here it might very well be that $k$ is not a simple eigenvalue, but instead for example a double eigenvalue, which occurs for $\Omega$ being the unit circle. Then we search for both new eigenvalues and pick the one closer to $k_{\star}$, but still larger than $k_{\star}$. This search costs more time than the approximation algorithm.

In numerical experiments it seems that the two eigenvalues of the double Dirichlet eigenvalue split such that one eigenvalue escapes subjectively faster from the double Dirichlet eigenvalue the longer the new boundary interval ${\Gamma_{\Delta}}$ is and the other eigenvalue subjectively slower. This is reminiscent of the behavior of the perturbation of a double eigenvalue in [8], where the eigenvalue splits in an eigenvalue with difference $\mathcal{O}(\varepsilon^{2})$ and an eigenvalue with difference $\mathcal{O}(1/|\log(\varepsilon)|)$, where $\varepsilon$ is a value associated to the perturbation.

• Line 23:

Next, we expand the boundary interval, which we established in Line 10-21. We expand it on both ends by a length $\varepsilon/2$, whose factor $1/2$ is again chosen due to good numerical approval for minimizing runtime. Then we compute an approximation $k$ to the new characteristic value. If $|k-k_{\star}|<C_{\mathrm{tol}}$, we end the algorithm, if $k+C_{\mathrm{tol}}<k_{\star}$, we extend the boundary interval once again, else decrease $\varepsilon$.

• Line 26:

To compute an approximation to the new characteristic value, we use the same setting as in Line 13: The complex domain $V$ encircling $k$ and $k_{\star}$ is an ellipse with center $(k+k_{\star})/2$ and semi-major axis $(k-k_{\star})\cdot 0.55$ and semi-minor axis $(k-k_{\star})\cdot 0.1$. A discretization to the operator $\mathcal{A}(k)$ is computed using the algorithm described in [2]. For the complex derivative of said operator we used the rough approximation $(\mathcal{A}(w+0.01)-\mathcal{A}(w))/0.01$. The integral is approximated with a inbuilt-process.

The approximation may again yield the same result as the former characteristic value, that is $k$, this happens when ${\Gamma_{\Delta}}$ is too long.

In this while-loop, it never happened that $k$ is not a simple eigenvalue.

Remark 4.1**.**

*If the function \partial\Omega\ni z\mapsto\big{(}\partial_{\nu_{z}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},z)\cdot\partial_{\nu_{z}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(y,z)\big{)}\in\mathbb{R} oscillates strongly on the boundary it might yield better results, when multiple, but smaller, boundary intervals are applied. The thought behind this is that using one long boundary interval might intersect the disadvantageous part of the function $\partial_{\nu_{z}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},z)\cdot\partial_{\nu_{z}}\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(y,z)$ and thus decrease the intensity of $\mathrm{Z}^{k}_{({\Gamma_{\textrm{D}}},{\Gamma_{\textrm{N}}})}(x_{\textrm{S}},y)$. *

5 Numerical Implementation and Tests

Our first numerical test shows the algorithm in the best case scenario. We have the domain $\Omega=\{x\in\mathbb{R}^{2}\mid\left\lVert x\right\rVert_{\mathbb{R}^{2}}<1\}$, the signal point $x_{\textrm{S}}=(0,0)^{\mathrm{T}}$, the target characteristic value $k_{\star}=1$ and $C_{\mathrm{tol}}=10^{-3}$ and $\varepsilon=0.1$. We remark here that the next higher Dirichlet characteristic value is a simple one at approximately $2.40482$. We let the receiving point $y\in\{(0,r)^{\mathrm{T}}\in\mathbb{R}^{2}\mid r>0\}$ vary. Here we want to mention that our implementation of the Zaremba function, as described in Section 4, comment on Line 3, yields a non-zero imaginary part for the Zaremba function, the same holds true for the approximation to the characteristic value $k$ as described in Section 4, comment on Line 13. We always choose the real part whenever in question. The number of discretization points for the operator $\mathcal{A}(k)$ was $3\cdot 64$. The results are displayed in Table 1. The Zaremba functions with Dirichlet boundary conditions and with final mixed boundary conditions, for the case $y=(0,0.5)^{\mathrm{T}}$, are displayed in Figure 5.

Our second numerical test shows the algorithm for a higher target characteristic value $k_{\star}$, namely $k_{\star}=15.4$. We have as the domain $\Omega$ the unit circle $\{x\in\mathbb{R}^{2}\mid\left\lVert x\right\rVert_{\mathbb{R}^{2}}<1\}$, as the signal point $x_{\textrm{S}}=(0,0)^{\mathrm{T}}$ and $C_{\mathrm{tol}}=10^{-3}$ and $\varepsilon=0.05$. We remark here that the next higher Dirichlet characteristic value has multiplicity two and is at approximately $15.5898$. We let the receiving point $y\in\{(0,r)^{\mathrm{T}}\in\mathbb{R}^{2}\mid r>0\}$ vary. The number of discretization points for the operator $\mathcal{A}(k)$ is $4\cdot 48$. The results are displayed in Table 2. The Zaremba functions with Dirichlet boundary conditions and with final mixed boundary conditions, for the case $y=(0,0.5)^{\mathrm{T}}$, are displayed in Figure 6.

Our third numerical test shows the algorithm for a different domain $\Omega$ namely a kite-shaped domain given by the following description for its boundary

[TABLE]

for $\tau\in[0,2\pi)$. The target characteristic value is $k_{\star}=1.5$. The signal point $x_{\textrm{S}}=(-1.25,1.25)^{\mathrm{T}}$ and receiving point $y=(-1.25,-1.25)^{\mathrm{T}}$. $C_{\mathrm{tol}}=10^{-2}$ and $\varepsilon=0.05$. We remark here that the next higher Dirichlet characteristic value has multiplicity one and is at approximately $2.2099$. The number of discretization points for the operator $\mathcal{A}(k)$ is $4\cdot 48$. The result is displayed in Figure 7. The center of the Neumann boundary condition ${\Gamma_{\textrm{N}}}$ is at $(-1.191,-1.493)^{\mathrm{T}}$ with length $\approx 3.119$. $\mathrm{Z}_{\text{Dirichlet}}^{k}(x_{\textrm{S}},y)\approx-4.05\,\cdot\,10^{-5}$ and $\mathrm{Z}_{\text{End}}^{k}(x_{\textrm{S}},y)\approx-39.38$

In Figure 8, we have the same set-up but for $k_{\star}=11.5$, with the next higher Dirichlet characteristic value around $11.6507$. Here, the center of the Neumann boundary condition ${\Gamma_{\textrm{N}}}$ is at $(-1.142,0.641)^{\mathrm{T}}$ with length $\approx 0.632$. $\mathrm{Z}_{\text{Dirichlet}}^{k}(x_{\textrm{S}},y)\approx 0.148$ and $\mathrm{Z}_{\text{End}}^{k}(x_{\textrm{S}},y)\approx 1.68$.