---
abstract: |
    The results in [@O2] (see [@O1] for the quasistatics regime) consider the Helmholtz equation with fixed frequency $k$ and, in particular imply that, for $k$ outside a discrete set of resonant frequencies and given a source region $D_a\subset \RR^d$ ($d=\overline{2,3}$) and $u_0$, a solution of the homogeneous scalar Helmholtz equation in a set containing the control region $D_c\subset \RR^d$, there exists an infinite class of boundary data on $\partial D_a$ so that the radiating solution to the corresponding exterior scalar Helmholtz problem in $\RR^d\setminus D_a$ will closely approximate $u_0$ in $D_c$. Moreover, it will have vanishingly small values beyond a certain large enough “far-field" radius $R$ (see Figure \[fig:mainsetup\] for a geometric description).

    In this paper we study the minimal energy solution of the above problem (e.g. the solution obtained by using Tikhonov regularization with the Morozov discrepancy principle) and perform a detailed sensitivity analysis. In this regard we discuss the stability of the the minimal energy solution with respect to measurement errors as well as the feasibility of the active scheme (power budget and accuracy) depending on: the mutual distances between the antenna, control region and far field radius $R$, value of regularization parameter, frequency, location of the source.
author:
- Mark Hubenthal
- Daniel Onofrei
bibliography:
- 'controlpaper.bib'
title: Sensitivity analysis for active control of the Helmholtz equation
---

Introduction
============

During recent years, there has been a growing interest in the development of feasible strategies for the control of acoustic and electromagnetic fields with one possible application being the construction of robust schemes for sonar or radar cloaking.

One main approach controls fields in the regions of interest by changing the material properties of the medium in certain surrounding regions ([@Chan3; @Chan2; @Cummer; @Green2; @Green1; @Gunther-pr; @Pendry] and references therein). Several alternative techniques are proposed in the literature (other than transformation optics strategies) such as: plasmonic designs (see [@Alu] and references therein), strategies based on anomalous resonance phenomena (see [@Mil1; @Mil3; @Mil2]), conformal mapping techniques (see [@Ulf2; @Ulf1]), and complementary media strategies (see [@Chan]).

In the applied community, active designs for the manipulation of fields appear to have occurred initially in the context of low-frequency acoustics (or active noise cancellation). Especially notable are the pioneering works of Lueg [@Lueg] (feed-forward control of sound) and Olson & May [@Olson-May] (feedback control of sound). The reviews [@Elliot; @Fuller; @Tsynkov; @T1; @Peake; @Peterson], provide detailed accounts of past and recent developments in acoustic active control.

In the context of cloaking, the [**[*interior*]{}**]{} strategy proposed in [@Miller] employs a continuous active layer on the boundary of the control region while the [**[*exterior*]{}**]{} scheme discussed in [@OMV4; @OMV1; @OMV2; @OMV3] (see also [@CTchan-num]), uses a discrete number of active sources located in the exterior of the control region to manipulate the fields. The active exterior strategy for 2D quasistatics cloaking was introduced in [@OMV1], and, based on *a priori* information about the incoming field, the authors constructively described how one can create an almost zero field control region with very small effect in the far field. However, the proposed strategy did not work for control regions close to the active source. It “cloaked" large objects only when they are far enough from the source region (see [@OMV4]) and was not adaptable to three space dimensions. The finite frequency case was studied in the last section of [@OMV1] and in [@OMV3] (see also [@OMV4] for a recent review) where three (or four in 3D) active sources were needed to create a zero field region in the interior of their convex hull, while creating a very small scattering effect in the far field. The broadband character of the proposed scheme was numerically observed in [@OMV2]. All the above results were obtained assuming large amplitude and highly oscillatory currents on the active source regions. In this regard, in [@Norris] (see also [@Devaney0; @Miller]) the authors presented theoretical and numerical evidence that increasing the number of sources will decrease the power needed on each source and thus increase the feasibility of the scheme. Experimental designs and testing of active cloaking schemes in various regimes are reported in [@Du; @Ma; @Eleftheriades1; @Eleftheriades2].

In a recent development in [@O1], a general analytical approach based on the theory of boundary layer potentials is proposed for the active control problem in the quasi-static regime. By using the same integral equation approach, in [@O2] we extended the results presented in [@O1] to the active control problem for the exterior scalar Helmholtz equation. In particular, we characterized an infinite class of boundary functions on the source boundary $\partial D_a$ so that we achieve the desired manipulation effects in several mutually disjoint exterior regions. The method is novel in the sense that instead of using microstructures, exterior active sources modeled with the help of the above boundary controls are employed for the desired control effects. Such exterior active sources can represent velocity potential, pressure or currents.

In the current paper we study the active control problem in the context of cloaking, where one antenna $D_a$ protects a given control region $D_c$ from far field interrogation on $\partial B_{R}(\mathbf{0})$, with $R \gg 1$ (see Figure \[fig:mainsetup\]). We make use of the results in [@O2] and present a detailed sensitivity and feasibility study for the minimal norm solution of the problem.

The paper is organized as follows: In Section \[sec:background\] we recall the general result obtained in [@O2] in the context of exterior active cloaking. In Section \[sec:stability\] we present an $L^2$ conditional stability result for the minimal norm solution with respect to measurement errors of the incoming field. In Section \[sec:numerics\] we present the numerical details of the Tikhonov regularization algorithm with the Morozov discrepancy principle for the computation of the minimal norm solution of the exterior active cloaking problem in two dimensions. We will numerically observe the fact that the scheme requires large antenna powers in the far field and we will provide numerical support for our theoretical stability results. An important part of this section will be focused on the sensitivity analysis, where we will study: the dependence of the control results as a function of mutual distances between the antenna, control region and far field region; and the broadband character of our scheme in the near field region. Finally, in Section \[sec:conclusions\] we highlight the main results of the paper and discuss current and future challenges and extensions of our research.

Background {#sec:background}
==========

In this section we will recall the main result regarding the active exterior control problem for the Helmholtz equation obtained in [@O2]. We will focus only on the case where one active external source (antenna) $D_a$ protects a control region $D_c$ from an interrogating far field and maintains an overall small signature beyond a disk of large enough radius $R$.

The general setup for this question will be as follows. Let $B_{R} \subset \mathbb{R}^{d}$ be the ball of radius $R > 0$. We assume $\B0\in D_{a} \subset  B_{R}$ is the region inside a single antenna with sufficiently smooth boundary $\partial D_{a}$. We also let $D_{c} \Subset B_{R}$ be the control region, which is assumed to satisfy $\overline{D_{c}} \cap \overline{D_{a}} = \emptyset$ (see Figure \[fig:mainsetup\]). The numerical simulations in the current work are performed for the two dimensional case but the methods are adaptable to the three dimensional setting as well. Consider the function space $$\Xi = L^{2}(\partial D_{c}) \times L^{2}(\partial B_{R}),$$ endowed with the scalar product $$(\phi,\psi)_{\Xi} = \int_{\partial D_{c}}\phi_{1}(\mathbf{y})\overline{\psi}_{1}(\mathbf{y})\,dS_{\mathbf{y}} + \int_{\partial B_{R}} \phi_{2}(\mathbf{y})\overline{\psi}_{2}(\mathbf{y})\,dS_{\mathbf{y}},$$ which is a Hilbert space. For the remainder of the paper we will assume that every $L^2$ space of complex valued functions will be endowed with the usual inner product. As in [@O2] consider $K: L^{2}(\partial D_{a}) \to \Xi$, the double layer potential operator restricted to $\partial D_{c}$ and $\partial B_{R}$, respectively, defined by $$\label{EQ:K}
K\phi(\Bx,\Bz) =  (K_{1}\phi(\Bx), K_{2}\phi(\Bz)), \quad \phi \in L^{2}(\partial D_{a}),$$ where $$\begin{aligned}
K_1:L^{2}(\partial D_{a})\rightarrow L^{2}(\partial D_{c}), \, K_1\phi(\Bx) & = \int_{\partial D_{a}}\phi(\By)\frac{\partial \Phi(\Bx,\By)}{\partial {{\mbox{\boldmath ${\nu}$}}}_{\By}}ds_{\By}, \mbox{ for }\Bx\in
\partial D_c, \notag\\
K_2:L^{2}(\partial D_{a})\rightarrow L^{2}(\partial B_R), \, K_2\phi(\Bz) & = \int_{\partial D_{a}}\phi(\By)\frac{\partial \Phi(\Bz,\By)}{\partial {{\mbox{\boldmath ${\nu}$}}}_{\By}}ds_{\By}, \mbox{ for }\Bz\in
\partial B_R(\B0). \label{8}\end{aligned}$$ Here $\Phi(\Bx, \By)$ represents the fundamental solution of the relevant Helmholtz operator, i.e., $$\label{9}
\Phi(\Bx,\By)=\vspace{0.15cm}\left\{\begin{array}{ll}
\vspace{0.15cm}\displaystyle\frac{e^{ik|\Bx-\By|}}{4\pi|\Bx-\By|}, \mbox{ for } d=3\vspace{0.15cm}\\
\frac{i}{4}H_0^{(1)}(k|\Bx-\By|) , \mbox{ for } d=2\end{array}\right.$$ with $H_{0}^{(1)} = J_{0} + iY_{0}$ representing the Hankel function of first type. Note that in the integrals are to be understood as singular integrals defined through an operator extension from $C(\partial D_a)$. We will also consider $k$ such that $$\begin{aligned}
\label{ass-lambda-0}
                 && \mathbf{1)} \, -k^2 \mbox{ is not a Neumann eigenvalue for the Laplace operator in $D_a$ or $B_R(\B0)$},\nonumber\\
                 && \mathbf{2)} \, -k^2 \mbox{ is not a Dirichlet eigenvalue for the Laplace operator in $D_c$}.\end{aligned}$$ As in [@O2] we introduce the adjoint operator $K^{*}: \Xi \to L^{2}(\partial D_{a})$, which can be shown to satisfy $$K^{*}\psi(\mathbf{x}) = \int_{\partial D_{c}} \psi_{1}(\mathbf{y}) \overline{\frac{\partial \Phi(\mathbf{y},\mathbf{x})}{\partial \nu_{\mathbf{x}}}}\,dS_{\mathbf{y}} + \int_{\partial B_{R}} \psi_{2}(\mathbf{y}) \overline{\frac{\partial \Phi(\mathbf{y},\mathbf{x})}{\partial \nu_{\mathbf{x}}}}\, dS_{\mathbf{y}}, \quad \mathbf{x} \in \partial D_{a}. \label{eq:dlpotentialadjoint}$$

This paper proposes a sensitivity study for the following problem: Let $V\Subset D_c$ and $R'>R$. For a fixed wave number $k > 0$ and fixed $0 < \mu \ll 1$, find a function $h \in C(\partial D_{a})$ such that there exists $u \in C^{2}(\RR^{n} \setminus \overline{D_{a}}) \cap C^{1}(\RR^{n} \setminus D_{a})$ solving $$\left\{ \begin{array}{rl}
(\Delta + k^{2}) u(\mathbf{x}) & = 0 \quad \mathbf{x} \in \RR^{n} \setminus \overline{D_{a}}\\
u & = h \quad \textrm{ on $\partial D_{a}$}\\
\|u - f_{1}\|_{C(\overline{V})} & = \mathcal{O}(\mu)\mbox { and }\|u\|_{C(\RR^{n} \setminus B_{R'}(\mathbf{0}))} = \mathcal{O}(\mu),
\end{array}
\right. \label{eq:inverseproblem}$$ where $f_{1}$ is a solution of the Helmholtz equation in a neighborhood of the control region $D_{c}$. In fact, by using the operator $K$ and regularity arguments it is shown in [@O2] that a class of solutions for problem can be obtained by considering the following problem: for a fixed wave number $k>0$ satisfying conditions , a given function $f=(f_1,0) \in \Xi$ and $\mu > 0$, find a density function $\phi \in C(\partial D_{a})$ such that $$\|K\phi - f\|_{\Xi} \leq \mu. \label{eq:controlinequality}$$

Problem is in fact a Fredholm integral equation of the first kind, and it was studied in a very general setting in [@O2]. There the authors proved that the bounded and compact operator $K$ is also one-to-one and has a dense (but not closed) range, thus proving the existence of a class of solutions for (and thus for ). However, the fact that $K$ is compact and that its range is not closed also implies that problem is ill-posed. By using regularization, one can approximate a solution to problem with an arbitrary level of accuracy $\mu \ll 1$. There are several methods known in the literature, but we will use in this paper the Tikhonov regularization method [@reg1; @reg2]. This method, when applied to the operator $K:L^2(\partial D_a)\rightarrow \Xi$, proposes a solution $\phi_{\alpha}\in C(\partial D_{a})$ of the form $$\label{Tikhonov}
\phi_{\alpha}=(\alpha I+K^*K)^{-1}K^*f,  \textrm{ for } 0<\alpha \ll 1,$$ where $\alpha$ is a suitably chosen regularization parameter.

It is known that $\Vert K\phi_\alpha-f\Vert_\Xi\rightarrow 0$ as $\alpha\rightarrow 0$, (see [@Kirsch-Book96], Theorem 2.16), but the optimal choice of $\alpha$ is an essential step in designing a feasible method (e.g., finding a minimal norm solution), and there are various modalities to do this. In this paper we will use the Morozov discrepancy principle associated to the following weighted residual: $$\label{disc-fun}
E(\phi,h)= \displaystyle\frac{1}{\|h_{1}\|_{L^{2}(\partial D_{c})}^{2}}\| K_1\phi  - h_{1}\|_{L^{2}(\partial D_{c})}^{2} + \frac{1}{2\pi R} \|K_2\phi \|_{L^{2}(\partial B_{R})}^{2},$$ for every given $h=(h_1,0)\in \Xi$. The reasoning behind using the weighted residual discrepancy functional defined at \[disc-fun\] is as follows. Due to the asymptotic behavior of $\frac{\partial \Phi(\mathbf{x},\, \mathbf{y})}{\partial \nu_{\mathbf{y}}} = \mathcal{O}(|\mathbf{x}-\mathbf{y}|^{-1/2})$ as $|\mathbf{x}-\mathbf{y}| \to \infty$, we have that given a fixed density $\phi$, $\|K\phi\|_{L^{2}(\partial B_{R})} = \mathcal{O}(1)$ as $R \to \infty$. In other words, using the space $L^{2}(\partial B_{R})$ with the standard surface measure is not really suited to the decay properties of double layer potential solutions, because the decay of the normal derivative $\partial_{\nu}\Phi$ is too weak. Similarly, we use the relative norm $$\label{eq:F1}
\frac{\|K_1\phi - h_{1}\|_{L^{2}(\partial D_{c})}}{\|h_{1}\|_{L^{2}(\partial D_{c})}}$$ on $\partial D_{c}$ because this is a useful quantity for determining how good the control is, regardless of the norm of $h_{1}$. Thus our procedure for finding an approximate solution for problem is to first make use of the Tikhonov regularization for the operator $K:L^2(\partial D_a)\rightarrow \Xi$ as described in to obtain $\phi_\alpha$ and then apply the Morozov’s discrepancy principle for the unique choice of $\alpha$ ([@Kress-Book99]), i.e. such that $$\label{Morozov0}
E(\phi_{\alpha}, f)=\delta^2,$$ with $\displaystyle\delta^2\leq \mu^2\min\left\{\frac{1}{2\|f_{1}\|^2_{L^{2}(\partial D_{c})}}, \frac{1}{4\pi R}\right\}$.

In what follows, we will account for noise and measurement errors and will consider with $f=(f_1,0)\in \Xi$ replaced by $f_\epsilon=(f_{\epsilon,1}, f_{\epsilon,2})\in \Xi$, given by $$\label{random-f}
f_{\epsilon} = (f_{1} + \epsilon s, \, 0)\in \Xi,$$ where $s\in L^{2}(\partial D_{c})$ is a random perturbation with $\|s\|_{L^{2}(\partial D_{c})} \leq 2 \|f_{1}\|_{L^{2}(\partial D_{c})}$ and $f_1$ is a solution of the Helmholtz equation in a neighborhood of the control region $D_{c}$. We mention that in the numerical experiments of Section \[sec:numerics\], $f_{1}$ denotes the $k$ frequency component of the far field of a far field observer. Note that this assumption about the interrogating signal ensures that $f_1$ is a solution of the Helmholtz equation in $B_R$. In the noisy case (i.e. when $f$ is replaced by $f_\epsilon$) equation becomes $$\label{Morozov}
E(\phi_{\alpha}, f_\epsilon)=\delta^2,$$ where $\phi_{\alpha}=(\alpha I+K^*K)^{-1}K^*f_\epsilon$ is the Tikhonov regularization solution. From the definition of $E$ and classical results, [@Kirsch-Book96; @Kress-Book99], it follows that admits at least a solution $\alpha$. Moreover, as we will discuss in Section \[sec:stability\], motivated by numerical evidence, we formulate the hypothesis that there exists $\epsilon_0>0$ such that for each $\epsilon\in (0,\epsilon_0)$, problem has a unique solution $\alpha(\epsilon)$ which uniquely defines a differentiable function $\epsilon\rightsquigarrow\alpha(\epsilon)$. We will study the minimal norm solution uniquely determined by , discuss its stability for given noisy data in $\Xi$ and, in the case of data corresponding to a point source, analyze its sensitivity with respect to parameters such as: mutual distances between $D_a$, $D_c$ and $B_R(\B0)$; wave number $k$; and the location of the point source.

Stability estimate for the Tikhonov regularization {#sec:stability}
==================================================

In this section we present analytical and numerical arguments which indicate the stability of the minimum norm solution $\phi_{\alpha}$ with respect to noise level $\epsilon$ for a given fixed discrepancy level $\delta$. Next, we present below Lemma \[lemma:phif1constant\] which will provide bounds for $\|f_{1}\|_{L^{2}(\partial D_{C})}$ and $\alpha$ in terms of the operatorial norm of $K_{1}^{*}$.

Let $0<\delta<\frac{1}{\sqrt{2}}$ and $z=(z_1,0) \in \Xi'$ with $z_1\neq 0$. Consider the Tikhonov regularization solution $\phi_{\alpha}=(\alpha I + K^{*}K)^{-1}K^*z \in C(\partial D_a)$, with $\alpha$ such that $\|K\phi_{\alpha} - z\|_{\Xi'} \leq \delta$. Then we have $$\begin{aligned}
\|z_{1}\|_{L^2(\partial D_c)} & \leq 4\|K_1^*\|_{\mathcal O}\|\phi_{\alpha}\|_{L^{2}(\partial D_{a})},\label{lb-phi}\\
\alpha & \leq 4\delta \Vert K_1^*\Vert^2_{\mathcal O},\label{ub-alpha}\end{aligned}$$ where $K_1^*$ is the adjoint operator for $K_1$ defined by and $\|\cdot\|_{\mathcal O}$ denotes the operatorial norm. \[lemma:phif1constant\]

We will start with the proof of . Note that since $E(\phi_{\alpha}, z)=\delta^{2}$, we have $$\label{5.17}
\|K_{1}\phi_{\alpha} - z_{1}\|_{L^{2}(\partial D_{c})}^{2} \leq \delta^{2}\|z_{1}\|_{L^2(\partial D_c)}^{2}.$$ From we obtain $$\|K_{1}\phi_{\alpha}\|_{L^2(\partial D_c)}^{2} - 2 \, {\textrm{Re}}(K_{1}\phi_{\alpha}, \, z_{1})_{L^2(\partial D_c)} + \|z_{1}\|_{L^2(\partial D_c)}^{2} \leq \delta^{2}\|z_{1}\|_{L^2(\partial D_c)}^{2},$$ and this implies $$\|z_{1}\|_{L^2(\partial D_c)}^{2}(1 -\delta^{2}) \leq 2 \, {\textrm{Re}}(\phi_{\alpha}, \, K_{1}^{*}z_{1})_{L^2(\partial D_a)} \leq 2 \|\phi_{\alpha}\|_{L^2(\partial D_a)}\|K_{1}^{*} z_{1}\|_{L^2(\partial D_a)}. \label{pas1}$$ Then gives $$\begin{aligned}
\|\phi_{\alpha}\|_{L^2(\partial D_a)} & \geq \frac{\|z_{1}\|_{L^2(\partial D_c)}^{2}(1 - \delta^{2})}{2 \|K_{1}^{*}z_{1}\|_{L^2(\partial D_a)}}
= {\left ( \frac{1 - \delta^{2}}{2}\right )}{\left (  \frac{\|z_{1}\|_{L^2(\partial D_c)}}{\|K_{1}^{*}z_{1}\|_{L^2(\partial D_a)}} \right )}\|z_{1}\|_{L^2(\partial D_c)}\notag\\
& \geq \frac{1 - \delta^{2}}{2} {\left ( \frac{\|z_{1}\|_{L^2(\partial D_c)}}{\|K_{1}^{*}\|_{\mathcal O}}\right )}\notag\\
&\geq \frac{\|z_1\|_{L^2(\partial D_c)}}{4\|K_{1}^{*}\|_{\mathcal O}}.\label{margine}\end{aligned}$$

Next we proceed towards proving . From the definition of $\phi_{\alpha}$ we have $$\begin{aligned}
\alpha \phi_{\alpha} + K^{*}K\phi_{\alpha}&=& K^*z=K_1^*z_1,\nonumber\\
\alpha \phi_{\alpha}+K_2^{*}K_2\phi_{\alpha} &=& K_1^*(z_1-K_1\phi_{\alpha}).\label{5.18}\end{aligned}$$ Here we have used from and that $$\begin{aligned}
\label{adjunct-rel}
K^*\psi=K_1^*\psi_1+K^*_2\psi_2, \mbox{ for all }\psi\in\Xi,\nonumber\\
K^*Kv=K_1^*K_1v+K_2^*K_2v, \mbox{ for all }v\in L^2(\partial D_a).\end{aligned}$$ Multiplying with $\phi_{\alpha}$ in the sense of the usual scalar product in $L^2(\partial D_a)$, we obtain $$\label{rel11}
\alpha \|\phi_{\alpha}\|^2_{L^2(\partial D_a)} +\|K_2\phi_{\alpha}\|^2_{L^2(\partial D_a)}  = (K_1^*(z_1-K_1\phi_{\alpha}),\phi_{\alpha})_{L^2(\partial D_a)}.$$ Using , and in we then have $$\alpha \|\phi_{\alpha}\|_{L^2(\partial D_a)} \leq \delta\|K_1^*\|_{\mathcal O}\|z_1\|_{L^2(\partial D_c)}\Longrightarrow
\alpha\leq 4\delta\|K_1^*\|^2_{\mathcal O} \label{alfa}.$$

Next, before presenting the main stability result of this section, i.e., Proposition \[theorem:stabilityestimate\] below, we must understand the conditions on $\epsilon > 0$ under which admits a unique solution $\alpha(\epsilon)$ with the property that the resulting function $\epsilon\rightsquigarrow \alpha(\epsilon)$ is differentiable. For this, we consider the function $g:(0,\infty)\times(0,\infty)\rightarrow (0,\infty)$ defined by $$g(\alpha, \epsilon) = E(\phi_\alpha, f_\epsilon)\label{eq:galphaeps},$$ where $f_{\epsilon} \in \Xi$ was introduced in , and $\phi_\alpha$ is the Tikhonov regularization solution introduced in . With this notation, can be rewritten as $$\label{Morozov-g}
g(\alpha, \epsilon) = \delta^2,$$ where $\delta$ is the desired fixed discrepancy level. By using classical results (e.g., [@Kirsch-Book96; @Kress-Book99]) it can be observed that for every $\epsilon$, admits at least one solution in $(0,\infty)$ and that $g$ defined by is differentiable with respect to positive $\alpha$ and $\epsilon$, respectively. In fact, it follows from classical arguments that a maximum value of $\alpha$ for a given $\epsilon$ exists. This solution of corresponds to the $L^2$ minimal energy solution and we will further refer to it as the Morozov solution.

For the remainder of the paper, unless otherwise specified, $C$ will denote a generic constant which depends only on the operator $K$, $d_c=diam(D_c)$ and $d=dist(\partial D_c,\partial D_a)$. The next Proposition states a central stability result concerning the Morozov solution of . We have,

Let $0<\delta$ be as above, and $f_\epsilon$ and $f_1$ as defined in . For every $\epsilon \geq 0$ consider $\phi_{\alpha_{\epsilon}}=(\alpha_\epsilon I + K^{*}K)^{-1}K^{*}f_{\epsilon}\in C(\partial D_a)$ with $\alpha_\epsilon=\alpha(\epsilon)$ the Morozov solution of . Then we have, $$\displaystyle\frac{\| \phi_{\alpha_\epsilon} - \phi_{\alpha_0}\|_{L^{2}(\partial D_{a})}}{\|\phi_{\alpha_{\epsilon}}\|_{L^{2}(\partial D_{a})}} \leq \displaystyle\frac{ \left|\displaystyle\frac{\alpha_{\epsilon}}{\alpha_0} - 1\right| + \sqrt{ \left|\displaystyle\frac{\alpha_{\epsilon}}{\alpha_0} -1 \right|^{2} + 16\displaystyle \frac{\epsilon{\left (  2\delta + C\delta\epsilon +C\epsilon\right )}}{\alpha_0}\|K^*_{1}\|_{\mathcal O}}}{2}.$$\[theorem:stabilityestimate\]

Fix $\epsilon > 0$ and let $f=f_\epsilon\mbox{ for } \epsilon=0$. Let us recall that $\alpha(\epsilon)$ is uniquely implicitly defined by the equation $E((\alpha_{\epsilon}I+K^*K)^{-1}K^*f_{\epsilon}, f_{\epsilon})= \delta^2$ and by Lemma \[lemma-diff-alpha\] it will be differentiable in some interval $(0,\epsilon_0)$ for all wavenumbers $k$. Next consider $$\begin{aligned}
\alpha_{\epsilon}\phi_{\alpha_{\epsilon}} + K^{*}K\phi_{\alpha_{\epsilon}} & = K^{*}f_{\epsilon},\\
\alpha_{0} \phi_{\alpha_{0}} + K^{*}K\phi_{\alpha_0} & = K^{*}f.\end{aligned}$$ Subtracting, we obtain $$\begin{aligned}
 \alpha_{0}\phi_{\alpha_0}-\alpha_{\epsilon}\phi_{\alpha_{\epsilon}} + K^{*}K(\phi_{\alpha_0}-\phi_{\alpha_{\epsilon}}) & = K^{*}(f - f_{\epsilon}),\notag\\
\alpha_0(\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}}) + (\alpha_0 - \alpha_{\epsilon})\phi_{\alpha_{{\epsilon}}} + K^{*}K(\phi_{\alpha_0}-\phi_{\alpha_{\epsilon}}) & = K^{*}(f - f_{\epsilon}).\label{3.1.1}\end{aligned}$$ Integrating both sides of against $\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}}$ yields $$\begin{aligned}
& \displaystyle\alpha_0 \|\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}}\|_{L^2(\partial D_a)}^{2}  + (\alpha_0 - \alpha_{{\epsilon}})(\phi_{\alpha_{{\epsilon}}}, \, \phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}})_{L^2(\partial D_a)} + \|K(\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}})\|_\Xi^{2}\notag\\
& \quad\quad= \displaystyle(K(\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}}), \, f - f_{\epsilon})_{\Xi}= \displaystyle(K_{1}(\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}}), \, f_{1} - f_{\epsilon, 1})_{L^2(\partial D_c)}.\label{3.1.2}\end{aligned}$$ where, we have used and the fact that $f_2=f_{\epsilon,2}=0$ in the last equality above. Thus, $$\begin{aligned}
\label{ecuatie-1}
\!\!\!\!\!\!\!\!\!\alpha_0 \|\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}}\|_{L^2(\partial D_a)}^{2}  \leq & |\alpha_0 - \alpha_{{\epsilon}}|\|\phi_{\alpha_{{\epsilon}}}\|_{L^2(\partial D_a)}\|\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}}\|_{L^2(\partial D_a)} \nonumber\\
+ & \|f_{1} - f_{\epsilon, 1}\|_{L^2(\partial D_c)}\|K_{1}(\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}})\|_{L^2(\partial D_c)}\end{aligned}$$ Observe that $$\begin{aligned}
\|K_{1}(\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}})\|_{L^2(\partial D_c)} & \leq \|K_{1} \phi_{\alpha_0} - f_{ 1}\|_{L^2(\partial D_c)} + \|K_{1} \phi_{\alpha_{{\epsilon}}} - f_{{\epsilon},1}\|_{L^2(\partial D_c)} + \|f_{1} - f_{\epsilon, 1}\|_{L^2(\partial D_c)}\nonumber\\
&\leq {\delta}\|f_{1}\|_{L^2(\partial D_c)}+ {\delta}\|f_{{\epsilon},1}\|_{L^2(\partial D_c)} + C\epsilon \|f_{1}\|_{L^2(\partial D_c)}\nonumber\\
& = {\left (  2\delta + C\delta\epsilon +C\epsilon\right )}\|f_{1}\|_{L^2(\partial D_c)}. \label{ecuatie-2}\end{aligned}$$ where $f_{\epsilon,1}=f_{1} + \epsilon s$ with $\|s\|_{L^{2}(\partial D_{c})} \leq C \|f_{1}\|_{L^{2}(\partial D_{c})}$, and we have used the definition of $\phi_{\alpha_{\epsilon}}$ and in the inequalities above. Using in we obtain $$\begin{aligned}
\alpha_0 \|\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}}\|_{L^2(\partial D_a)}^{2} & \leq |\alpha_0 - \alpha_{{\epsilon}}| \|\phi_{\alpha_{{\epsilon}}}\|_{L^2(\partial D_a)}\|\phi_{\alpha_0} - \phi_{\alpha_{{\epsilon}}}\|_{L^2(\partial D_a)} \notag\\
& + \epsilon{\left (  2\delta + C\delta\epsilon +C\epsilon\right )}\|f_{1}\|^2_{L^2(\partial D_c)}. \label{ecuatie-3}\end{aligned}$$ If we define $A := \frac{\|\phi_{\alpha_{\epsilon}} - \phi_{\alpha_{0}} \|_{L^{2}(\partial D_{a})}}{\|\phi_{\alpha_{{\epsilon}}}\|_{L^{2}(\partial D_{a})}}$, inequality implies that $$\begin{aligned}
\alpha_0 A^{2} & \leq |\alpha_0 - \alpha_{{\epsilon}}| A + \frac{\epsilon {\left (  2\delta + C\delta\epsilon +C\epsilon\right )} \|f_{1}\|_{L^{2}(\partial D_{c})}^{2}}{\|\phi_{\alpha_{{\epsilon}}}\|_{L^{2}(\partial D_{a})}^{2}} \notag\\
& \leq |\alpha_0 - \alpha_{{\epsilon}}| A + 16\epsilon {\left (  2\delta + C\delta\epsilon +C\epsilon\right )}\|K^*_1\|^2_{\mathcal O}, \label{quadratic}\end{aligned}$$ where we have used of Lemma \[lemma:phif1constant\] in the last inequality above. Next, consider $$h(A) := A^{2} - \left|\displaystyle\frac{\alpha_{{\epsilon}}}{\alpha_0} - 1 \right| A - 16\|K^*_1\|^2_{\mathcal O}\displaystyle\frac{\epsilon {\left (  2\delta + C\delta\epsilon +C\epsilon\right )}}{\alpha_0}.$$ Then, from we have $h(A) \leq 0$ and this implies that $$A \leq \displaystyle\frac{ \left|\displaystyle\frac{\alpha_{{\epsilon}}}{\alpha_0} - 1 \right| + \sqrt{\left|\displaystyle\frac{\alpha_{{\epsilon}}}{\alpha_0} - 1 \right|^{2} + 64\|K^*_1\|^2_{\mathcal O}\displaystyle\frac{\epsilon {\left (  2\delta + C\delta\epsilon +C\epsilon\right )}}{\alpha_0}}}{2}$$ which completes the proof.

Regarding the monotonic character of $g$, we note that, as suggested by the numerics, $g$ is not in general *globally* monotonic with respect to $\alpha$ as can be seen in Figure \[fig:nonincreasing\_F\], which considers an antenna of radius $a = 0.01$, region of control characterized in polar coordinates by $r_{1} = 0.011$, $r_{2} = 0.015$, $\theta \in [3\pi/4, 5\pi/4]$, wave number $k=10$, $f_\epsilon$ given by with $f_1=\frac{1}{4}H^{(1)}_0(k|\Bx-\Bx_0|)$ with $\mathbf{x}_{0} = [10000,0]^{T}$, and noise level $\epsilon=0.005$. But on the other hand, for the same geometry and functional settings as in Figure \[fig:nonincreasing\_F\], we observe in Figure \[fig:plot\_g\_dg\_alpha\] that for each $\epsilon<0.015$, $g(\alpha, \epsilon) = E(\phi_{\alpha_{\epsilon}}, \, f_{\epsilon})$ is strictly increasing with respect to $\alpha$ in the interval $(10^{-4}, 1)$. Moreover, for every $\epsilon<0.015$ the Morozov solution $\alpha(\epsilon)$ is the unique solution of in $(10^{-4}, 1)$.

In fact, for the same geometry and functional settings as in Figure \[fig:plot\_g\_dg\_alpha\] and for $k \in [1,100]$ and $\epsilon = 0.005$, Table \[tab:pktable\] summarizes the values of $p_k>0$ for which $g(\alpha, \epsilon)$ is locally strictly monotonic with respect to $\alpha$ in an interval $(10^{-p_k}, 1)$, as well as the value of the Morozov solution for each $k$ (also see Figure \[fig:pk\_plot\]). Together with Figure \[fig:plot\_g\_dg\_alpha\] where $k=10$ and $\epsilon$ is varied in the interval $[0,0.015]$, this suggests that the Morozov solution $\alpha(\epsilon)$ of satisfies $\alpha(\epsilon)\in(10^{-p_k}, 1)$ at least for $\epsilon < 0.015$. This in turn together with the strict monotonicity $g$ implies the existence of a unique solution $\alpha(\epsilon)$ for . Then, uniqueness together with the fact that $\frac{\partial}{\partial \alpha}g(\alpha, \epsilon)\neq 0$ in $ (10^{-p_k}, 1)\times (0,0.015)$ (for $k=10$ shown in Figure \[fig:plot\_g\_dg\_alpha\]) implies the differentiability of the Morozov solution $\alpha(\epsilon)$ by using the implicit function theorem.

  --------------- ---------------- ------------------ ------ ---------------- ------------------
  $k$             $-p_{k}$         Morozov $\alpha$   $k$    $-p_{k}$         Morozov $\alpha$
  \[0.5ex\] 1.0   -5.74057337341   0.0021397          51.0   -5.63387601498   0.023987
  6.0             -5.15371857022   0.0022445          56.0   -5.5538513243    0.028226
  11.0            -4.61487654411   0.0027985          61.0   -5.58052339557   0.032734
  16.0            -4.43348001737   0.0037643          66.0   -7.93866063316   0.038053
  21.0            -7.22374205154   0.0052228          71.0   -7.78393971408   0.043891
  26.0            -6.73292218326   0.0072707          76.0   -7.60786606025   0.048884
  31.0            -6.41280555768   0.0098052          81.0   -7.41580451387   0.05425
  36.0            -6.18339417827   0.012607           86.0   -7.2077344143    0.060871
  41.0            -5.97530913764   0.015782           91.0   -7.00500135956   0.066982
  46.0            -5.78858575492   0.01971            96.0   -6.90364624001   0.07278
  --------------- ---------------- ------------------ ------ ---------------- ------------------

  : Table of values $-p_{k}$ such that $g(\alpha, \epsilon)$ is increasing with respect to $\alpha$ for $\alpha \geq 10^{-p_{k}}$. $\epsilon$ is fixed at $0.005$ in this case.[]{data-label="tab:pktable"}

![Threshold value $p_{k}$ for which $F(\alpha)$ is increasing when $\alpha > 10^{-p_{k}}$. Also shows the value of $\alpha$ for which $F(\alpha) = \delta^{2}$ with the same setting as in Figure \[fig:plot\_g\_dg\_alpha\], where $\delta = 0.02$ and the noise level $\epsilon = 0.005$.[]{data-label="fig:pk_plot"}](pk_plot){width="\textwidth"}

For simplicity of notation, in what follows we will write sometimes $\alpha$ instead of $\alpha_\epsilon=\alpha(\epsilon)$ and we will use $\alpha'$ and $f_{\epsilon,i}'$ to denote $\frac{d\alpha}{d\epsilon}$ and $\frac{df_{\epsilon,i}}{d\epsilon}$ respectively. Motivated by the above numerics, we formulate the following more general hypothesis:

\[hypothesis-I\] Assume the same geometrical setup as in Section \[sec:background\] and let $f_\epsilon, f_1$ be as in . Then there exists $p_0>0$ and $\epsilon_0>0$ such that $\frac{\partial}{\partial \alpha}g(\alpha, \epsilon)\neq 0$ in $ (10^{-p_0}, 1)\times (0,\epsilon_0)$ for all wave numbers $k$, and the Morozov solution $\alpha(\epsilon)$ is the unique solution of in $(10^{-p_0}, 1)$.

For example, as shown in Table \[tab:pktable\] and Figure \[fig:pk\_plot\], for $f_\epsilon, f_1$ as in with $s=\widehat{\nu}\Vert f_1\Vert_{ L^2(\partial D_c)}$ where $\widehat{\nu}\in L^2(\partial D_c)$ is a random perturbation with $\Vert \widehat{\nu}\Vert_{ L^2(\partial D_c)}=1$, and for the same geometry and data as in Figure \[fig:plot\_g\_dg\_alpha\], we have that Hypothesis \[hypothesis-I\] is satisfied for $p_0=10^{-4}$ and $\epsilon_0=0.015$ for all $k=\overline{1,100}$. Thus, whenever Hypothesis \[hypothesis-I\] is satisfied, the definition of $\alpha(\epsilon)$ and the implicit function theorem imply:

\[lemma-diff-alpha\] There exists $\epsilon_0>0$ such that for every $\epsilon\in(0,\epsilon_0)$, the function $\alpha:(0,\epsilon_0)\rightarrow(0,\infty)$, where for each $\epsilon\in(0,\epsilon_0)$, $\alpha(\epsilon)$ represents the Morozov solution of , will be differentiable for all wave numbers $k$.

The next Lemma is a technical result needed in the stability estimate obtained in Corollary \[corolar-1\].

\[teorema-int\] Let $f_{1}$ be a solution of the Helmholtz equation in a neighborhood of $D_c$ satisfying the following source type condition: $$\label{source-cond}
\displaystyle\left\Vert K_1\psi_0-\frac{f_1}{\Vert f_1\Vert_{L^{2}(\partial D_{c})}}\right\Vert_{L^{2}(\partial D_{c})}\!\!\!\!\leq C\delta \mbox{ for some }\psi_0\in L^2(\partial D_a) \mbox{ with } \Vert\psi_0\Vert_{L^2(\partial D_a)}\leq C \delta.$$ Assume $R$ (radius of $B_R(\B0)$) is such that, $$\label{R-bound}
\Vert f_1\Vert_{L^{2}(\partial D_{c})} \leq \sqrt{\pi R}.$$ Consider $s_\epsilon=f_{1} + \frac{\epsilon}{2} \widehat{\nu}\|f_{1}\|_{L^{2}(\partial D_{c})}$ where $\widehat{\nu}\in L^2(\partial D_c)$ is a random perturbation with $\|\widehat{\nu}\| \leq 1$. Assume the same functional framework as in Proposition \[theorem:stabilityestimate\] and that Hypothesis \[hypothesis-I\] holds true in the case when $f_\epsilon$ is given by $$\label{hyp-1}
f_{\epsilon} = (f_{\epsilon,1}, f_{\epsilon,2})= (f_{1} + \epsilon s_\epsilon, \, 0).$$ Then, there exists $\epsilon_0>0$ such that the Morozov solution of equation $\alpha=\alpha(\epsilon)$ satisfies $$\alpha|\alpha'| \leq C\frac{\delta^2}{\sqrt{\alpha}}, \mbox{ for all } \epsilon<\epsilon_0.$$

Define the weights $$\begin{aligned}
w_{1} & := \frac{1}{\|f_{\epsilon, 1}\|_{L^{2}(\partial D_{c})}^{2}}\\
w_{2} & := \frac{1}{2 \pi R}.\end{aligned}$$ and denote $T_{\alpha} := (K^{*}K + \alpha I)^{-1}$, $R_{\alpha} := T_{\alpha}K^{*}$. Then using the Einstein summation convention, we may write $$E(\phi_{\alpha}, f_{\epsilon})= w_{i}\|K_{i}\phi_{\alpha} - f_{\epsilon,i} \|_{L^{2}(W_{i})}^{2} = w^{i} {\left ( K_{i}\phi_{\alpha} - f_{\epsilon,i}, \, K_{i}\phi_{\alpha} - f_{\epsilon,i}\right )}_{L^{2}(W_{i})},$$ where $W_{1} = \partial D_{c}$, $W_{2} = \partial B_{R}$. Next, as in Lemma \[lemma-diff-alpha\], we observe that Hypothesis \[hypothesis-I\] together with the implicit function theorem imply the uniqueness and differentiability of $\alpha(\epsilon)$, on the interval $\epsilon \in (0, \epsilon_0)$ for some $\epsilon_0>0$, where $\alpha(\epsilon)$ is uniquely and implicitly defined by the equation $E((\alpha_{\epsilon}I+K^*K)^{-1}K^*f_{\epsilon}, \, f_{\epsilon})= \delta^2$. Differentiating the equation $E((\alpha_{\epsilon}I+K^*K)^{-1}K^*f_{\epsilon}, \, f_{\epsilon})= \delta^2$ with respect to $\epsilon$ and noting that $\delta$ is fixed, we obtain $$\begin{aligned}
0  = \partial_{\epsilon}E(\phi_{\alpha}, f_{\epsilon})& = 2w_{i}\,{\textrm{Re}}{ {\left ( K_{i}\phi_{\alpha}' - f_{\epsilon,i}', \, K_{i}\phi_{\alpha} - f_{\epsilon, i}\right )}}_{L^{2}(W_{i})} \notag\\
&  \quad - 2(w_{1})^{2}\,{\textrm{Re}}{{\left ( f_{\epsilon,1}', f_{\epsilon, 1}\right )}}_{L^{2}(\partial D_{c})}\|K_{1}\phi_{\alpha} - f_{\epsilon,1}\|_{L^{2}(\partial D_{c})}^{2}\label{eq:diffnorm1}.\end{aligned}$$

Next, from $(K^{*}K + \alpha I)\phi_{\alpha} = K^{*}f_{\epsilon}$ we observe that $$\phi_{\alpha}' = R_{\alpha}f_{\epsilon}' - \alpha' T_{\alpha}\phi_{\alpha}. \label{eq:phialphadiff}$$ Thus, we may write $$K_{i}\phi_{\alpha}' - f_{\epsilon, i}' = -\alpha' K_{i} T_{\alpha}\phi_{\alpha} + K_{i}R_{\alpha}f_{\epsilon}' - f_{\epsilon, i}'.\label{eq:sub1}$$ By using (\[eq:sub1\]) and we obtain that $$\begin{aligned}
\label{ec-11}
\!\!\!\!\!\!\!\!\!
2{\left ( K_{i}\phi_{\alpha}' - f_{\epsilon,i}', \, K_{i}\phi_{\alpha}\right )}_{L^{2}(W_{i})} &= -2\alpha' {\left ( T_{\alpha}\phi_{\alpha}, \, K_{i}^{*}K_{i}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})}\notag\\
& + 2{\left ( K_{i}R_{\alpha}f_{\epsilon}' - f_{\epsilon, i}', \, K_{i}\phi_{\alpha}\right )}_{L^{2}(W_{i})},\end{aligned}$$ and $$\begin{aligned}
\label{ec-12}
\!\!\!\!\!\!\!\!\!-2{\left (  f_{\epsilon,i}, \, K_{i}\phi_{\alpha}' - f_{\epsilon, i}'\right )}_{L^{2}(W_{i})} & = 2\alpha'{\left (  f_{\epsilon,i}, \, K_{i}T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(W_{i})} - 2{\left ( f_{\epsilon,i}, \, K_{i}R_{\alpha}f_{\epsilon}' - f_{\epsilon, i}'\right )}_{L^{2}(W_{i})}\notag\\
& = 2\alpha'{\left (  K_{i}^{*}f_{\epsilon,i}, \, T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})} - 2{\left ( f_{\epsilon,i}, \, K_{i}R_{\alpha}f_{\epsilon}' - f_{\epsilon, i}'\right )}_{L^{2}(W_{i})}.\end{aligned}$$

Let $P,Q$ be defined by $$\begin{aligned}
P & = 2w_{i}\left[ {\left ( K_{i}\phi_{\alpha} - f_{\epsilon, i}, \, K_{i}R_{\alpha}f_{\epsilon}' - f_{\epsilon,i}'\right )}_{L^{2}(W_{i})} + \alpha' {\left ( K_{i}^{*}f_{\epsilon,i} - K_{i}^{*}K_{i}\phi_{\alpha}, \, T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})} \right]\label{P}\\
Q & = 2(w_{1})^{2}\,{\left ( f_{\epsilon,1}', f_{\epsilon, 1}\right )}_{L^{2}(\partial D_{c})}\|K_{1}\phi_{\alpha} - f_{\epsilon,1}\|^{2}.\label{Q}\end{aligned}$$ Then from , used in we obtain $$\label{ec-15'}
0=\partial_{\epsilon}E(\phi_{\alpha}, f_{\epsilon})= {\textrm{Re}}(P) - {\textrm{Re}}(Q).$$ We focus first on $P$ introduced in . In this regard, let us define $$\label{P1}
L_{i} := {\left (  K_{i}\phi_{\alpha} - f_{\epsilon, i}, \, K_{i}R_{\alpha}f_{\epsilon}' - f_{\epsilon, i}'\right )}_{L^{2}(W_{i})}.$$ Observe that implies $$f_\epsilon'=(f'_{\epsilon,1}, f'_{\epsilon,2})=(f_1+\epsilon\widehat{\nu}\Vert f_1\Vert_{L^2(\partial D_c)},0)\label{P2}.$$ First note that by using classical arguments based on the singular value decomposition for $K:L^2(\partial D_a)\rightarrow\Xi$, one can adapt the results in [@Kirsch-Book96] (Theorem 2.7) and obtain, $$\label{regl-1}
\|KR_{\alpha}Kz - Kz \|_{\Xi} \leq C\sqrt{\alpha} \Vert z\Vert_{L^2(\partial D_a)}, \mbox{ for every } z\in L^2(\partial D_a).$$

Let $\displaystyle f=\left(\frac{f_1}{\Vert f_1\Vert_{L^2(\partial D_c)}},0\right)$ and $\displaystyle v=\frac{f_\epsilon'}{\Vert f_1\Vert_{L^2(\partial D_c)}}-f$. By using the definition of $E$ and $\Xi$, , , , , , Cauchy’s inequality and the triangle inequality in , we obtain $$\begin{aligned}
\label{P4}
2|w_iL_i|  &\leq& C\delta (\sqrt{w_1}\|K_1R_{\alpha}f_{\epsilon}' - f_{\epsilon,1}' \|_{L^2(\partial D_c)}+\sqrt{w_2}\|K_2R_{\alpha}f_{\epsilon}' - f_{\epsilon,2}' \|_{L^2(\partial B_R(\B0))})\nonumber\\
&\leq&C\delta (\|K_1R_{\alpha}(v+f) - (v_1+f_1) \|_{L^2(\partial D_c)}+\|K_2R_{\alpha}(v+f) - (v_2+f_2) \|_{L^2(\partial B_R(\B0))})\nonumber\\
&\leq& C\delta \|KR_{\alpha}(v+f) - (v+f) \|_{\Xi}\nonumber\\
&\leq& C\delta ( \|KR_{\alpha}f - f\|_{\Xi}\!+\!\|KR_{\alpha}v\! - \!v\|_{\Xi})\nonumber\\
&\leq& C\delta(\|KR_{\alpha}K\psi_0 - K\psi_0 \|_{\Xi}+\|KR_{\alpha}(K\psi_0 - f) \|_{\Xi}\!+\!\|K\psi_0 - f\|_{\Xi})+C\frac{\delta\epsilon}{\sqrt{\alpha}}\nonumber\\
&\leq& C\delta^2\sqrt{\alpha} +C\frac{\delta^2}{\sqrt{\alpha}}+C\delta^2+C\frac{\delta\epsilon}{\sqrt{\alpha}}\nonumber\\
&\leq& C\frac{\delta^2}{\sqrt{\alpha}},\end{aligned}$$ where Einstein summation convention was used and where, in the second inequality above we make use of to obtain $\sqrt{\frac{w_2}{w_1}}\leq 1$ and $\sqrt{w_1}\Vert f_1\Vert_{L^2(\partial D_c)}<C$ for small enough $\epsilon$, and in the fourth and fifth inequalities above we have used that $||R_\alpha||_{{\cal O}}\leq\frac{1}{2\sqrt{\alpha}}$ (e.g. see [@Kirsch-Book96]), and respectively, that $\epsilon < \delta$ and $\psi_0$ satisfies the source condition .

Expanding $P$ defined in and using the fact that $f_{\epsilon,2} = 0$ and $\phi_{\alpha} = T_{\alpha}K^{*}f_{\epsilon} = T_{\alpha}K_{1}^{*}f_{\epsilon,1}$, we obtain $$\begin{aligned}
\label{ec-13}
P & = \frac{2\alpha'}{\|f_{\epsilon,1}\|_{L^{2}(\partial D_{c})}^{2}} {\left (  K_{1}^{*}f_{\epsilon,1},\, T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})} - \frac{2\alpha'}{\|f_{\epsilon,1}\|_{L^{2}(\partial D_{c})}^{2}}{\left ( K_{1}^{*}K_{1}\phi_{\alpha}, \, T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})}\notag\\
& \quad - \frac{2 \alpha'}{2\pi R} {\left ( K_{2}^{*}K_{2}\phi_{\alpha},\, T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})} + 2w_{i}L_{i}\notag\\
& = \frac{2\alpha'}{\|f_{\epsilon,1}\|^{2}} {\left ( T_{\alpha}^{-1}\phi_{\alpha},\,T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})} - \frac{2\alpha'}{\|f_{\epsilon,1}\|^{2}}{\left ( K_{1}^{*}K_{1}\phi_{\alpha}, T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})}\notag\\
& \quad - \frac{2 \alpha'}{2\pi R} {\left ( K_{2}^{*}K_{2}\phi_{\alpha},\, T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})} + 2w_{i}L_{i}\notag\\
& = \frac{2\alpha'}{\|f_{\epsilon,1}\|^{2}} {\left ( (\alpha I + K_{2}^{*}K_{2})\phi_{\alpha},\,T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})} - \frac{2 \alpha'}{2\pi R} {\left ( K_{2}^{*}K_{2}\phi_{\alpha},\, T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})} + 2w_{i}L_{i}\notag\\
& = \frac{2\alpha \alpha'}{\|f_{\epsilon,1}\|^{2}} {\left ( \phi_{\alpha},\,T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})} + \frac{\alpha'}{\|f_{\epsilon,1}\|^{2}}B{\left ( K_{2}^{*}K_{2}\phi_{\alpha}, \, T_{\alpha}\phi_{\alpha}\right )}_{L^{2}(\partial D_{a})} + 2w_{i}L_{i},\end{aligned}$$ where $B = 2 - \displaystyle\frac{\|f_{\epsilon,1}\|_{L^{2}(\partial D_{c})}^{2}}{\pi R}$ and we have used that $T_\alpha^{-1}=\alpha I+K_1^*K_1+K_2^*K_2$ in the third equality above. Observe that implies $B\geq 0$. Introduce the notation $\widetilde{K}_{2} := \sqrt{B}K_{2}$, and denote $\displaystyle v_{\alpha} := \frac{\phi_{\alpha}}{\|f_{\epsilon,1}\|_{L^{2}(\partial D_{c})}}$. Then becomes $$\begin{aligned}
\label{ec-14}
P & = 2 \alpha \alpha' {\left ( v_{\alpha}, \, T_{\alpha}v_{\alpha}\right )}_{L^{2}(\partial D_{a})} + \alpha' {\left (  \widetilde{K}_{2}^{*}\widetilde{K}_{2}v_{\alpha}, \, T_{\alpha}v_{\alpha}\right )}_{L^{2}(\partial D_{a})} + 2w_{i}L_{i}\notag\\
& = \alpha \alpha' {\left ( v_{\alpha}, \, T_{\alpha}v_{\alpha}\right )}_{L^{2}(\partial D_{a})} + \alpha' {\left (  (\alpha I + \widetilde{K}_{2}^{*}\widetilde{K}_{2})v_{\alpha}, \, T_{\alpha}v_{\alpha}\right )}_{L^{2}(\partial D_{a})} + 2w_{i}L_{i}.\end{aligned}$$ From [@kato] (see Section V.3.10), for any self-adjoint linear operator $A:H\rightarrow H$, where $H$ is a given Hilbert space (real or complex), we have that: $$\label{Kato}
0<\gamma = \inf_{\lambda \in \mathrm{Sp}(A)}\lambda \Longrightarrow (Ax, x)_H \geq \gamma(x,x)_H,$$ where $(\cdot,\cdot)$ in denotes the usual Hilbert product and where $\mathrm{Sp}(A)$ denotes the real spectrum of the operator $A$. Then, by using for the operator $T_\alpha$ we obtain $$\begin{aligned}
\label{Kato-1}
{\left (  v_{\alpha}, \, T_{\alpha}v_{\alpha} \right )}_{L^{2}(\partial D_{a})} &\geq \frac{1}{\alpha + \mu_{1}^{2}}\|v_{\alpha}\|_{L^{2}(\partial D_{a})}^{2}\geq \frac{1}{1 + \mu_{1}^{2}}\|v_{\alpha}\|_{L^{2}(\partial D_{a})}^{2}\notag\\
&\geq C \|v_{\alpha}\|_{L^{2}(\partial D_{a})}^{2},\end{aligned}$$ where we have used that $\displaystyle \frac{1}{\alpha+\mu_1^2} = \inf_{\lambda \in \mathrm{Sp}(T_\alpha)}\lambda$ with $\mu_1$ denoting the largest singular value of $K$.

Next consider $\displaystyle D_{\alpha} := \alpha I + \widetilde{K}_{2}^{*}\widetilde{K}_{2}$. Then, because $D_\alpha$ and $T_\alpha$ are linear, bounded, self-adjoint, invertible and positive definite operators, we have that $D_\alpha T_{\alpha}$ will also be linear, bounded, self-adjoint, invertible and have strictly positive eigenvalues. Indeed, for any eigenvalue-eigenvector pair $(x,\lambda)$ of $D_\alpha T\alpha$ we have $$D_\alpha T_\alpha x=\lambda x\Rightarrow T_\alpha x=\lambda D_\alpha^{-1}x\Rightarrow \lambda=\frac{( T_\alpha x,x)}{( D^{-1}_\alpha x,x)}\geq 0.$$ Observing that $0\notin \mathrm{Sp}(D_\alpha T_\alpha)$ we have the claim, and the positive definiteness of $D_\alpha T_\alpha$ follows. Thus we have $$\label{ec-15}
{\left (  (\alpha I + \widetilde{K}_{2}^{*}\widetilde{K}_{2})v_{\alpha}, \, T_{\alpha}v_{\alpha}\right )}_{L^{2}(\partial D_{a})} = {\left (  D_{\alpha}v_{\alpha},\, T_{\alpha}v_{\alpha}\right )}_{L^{2}(\partial D_{a})}\geq 0$$ From , and used in , we obtain $$\begin{aligned}
\label{ec-17}
|2w_{i}L_{i}| + |Q| & \geq |\alpha'| \left| \alpha {\left ( v_{\alpha}, \, T_{\alpha} v_{\alpha}\right )}_{L^{2}(\partial D_{a})} + {\left (  D_{\alpha}v_{\alpha}, \, T_{\alpha}v_{\alpha}\right )}_{L^{2}(\partial D_{a})} \right|\notag\\
& \geq |\alpha'| C\alpha\|v_{\alpha}\|_{L^{2}(\partial D_{a})}^{2} \notag\\
& \geq C \alpha |\alpha'| \|v_{\alpha}\|_{L^{2}(\partial D_{a})}^{2}.\end{aligned}$$ From and elementary algebraic manipulations we obtain, $$\label{ec-18}
|Q| \leq \frac{2 \|f_{1}\|_{L^{2}(\partial D_{c})}}{\|f_{\epsilon, 1}\|_{L^{2}(\partial D_{c})}^{4}} \|f_{\epsilon, 1}\|_{L^{2}(\partial D_{c})}\cdot \delta^{2} \|f_{\epsilon,1}\|_{L^{2}(\partial D_{c})}^{2} = \frac{2 \|f_{1}\|_{L^{2}(\partial D_{c})} \delta^{2}}{\|f_{\epsilon, 1}\|_{L^{2}(\partial D_{c})}} \leq C\delta^{2}$$ Recalling that $\displaystyle v_{\alpha} := \frac{\phi_{\alpha}}{\|f_{\epsilon,1}\|_{L^{2}(\partial D_{c})}}$, Lemma \[lemma:phif1constant\] implies $$\label{ec-19}
\|v_{\alpha}\|_{L^{2}(\partial D_{a})} \geq C.$$ Then from , , and used in we finally obtain the statement of the Lemma: $$\label{ecuatie-4}
\alpha |\alpha'|
 \leq C \delta^{2}+ C\frac{\delta^2}{\sqrt{\alpha}}\leq C\frac{\delta^2}{\sqrt{\alpha}}, \mbox{ for } \epsilon<\epsilon_0.$$

Note that $\|s_\epsilon\|_{L^{2}(\partial D_{c})}\leq (1+\frac{\epsilon}{2})\|f_{1}\|_{L^{2}(\partial D_{c})}$ and thus $f_\epsilon$ as defined above satisfies .

The next result is a simple consequence of Proposition \[theorem:stabilityestimate\] and Lemma \[teorema-int\]. We have,

\[corolar-1\] Assume the same notation and the same framework as in Proposition \[theorem:stabilityestimate\]. Assume also that for $d_c=diam(D_c)$ small enough there exists $d=dist(\partial D_c,\partial D_a)$ small enough so that $\alpha(0)=\alpha_0\approx\delta$. Then we have $$\displaystyle\frac{\| \phi_{\alpha_\epsilon} - \phi_{\alpha_0}\|_{L^{2}(\partial D_{a})}}{\|\phi_{\alpha_{\epsilon}}\|_{L^{2}(\partial D_{a})}} \leq C\sqrt{\epsilon}.$$

From Lemma \[teorema-int\] we obtain that $$\label{estimare}
|\alpha_\epsilon'| \alpha_\epsilon^{\frac{3}{2}}\leq C\delta^2.$$ Estimate and Cauchy’s theorem implies that $$\label{ecuatie-5}
\left|\frac{\alpha_\epsilon^{\frac{5}{2}}}{\alpha_0^{\frac{5}{2}}}-1\right|= \frac{5}{2}\epsilon\left|\alpha'_{\epsilon_*}\alpha_{\epsilon_*}^{\frac{3}{2}}\alpha_0^{-\frac{5}{2}}\right|\leq C\epsilon\delta^2\alpha_0^{-\frac{5}{2}}\leq   C\sqrt{\epsilon},$$ where $\epsilon_*\in(0,\epsilon)$ and we used that $\epsilon<\delta$. Next, simple algebraic manipulation and imply $$\left|\frac{\alpha_\epsilon}{\alpha_0}-1\right|\leq \left|\frac{\alpha_\epsilon^5}{\alpha_0^5}-1\right|\leq C\sqrt{\epsilon}\left(\frac{\alpha_\epsilon^{\frac{5}{2}}}{\alpha_0^{\frac{5}{2}}}+1\right)\leq C\sqrt{\epsilon}(2+ C\sqrt{\epsilon})\leq C\sqrt{\epsilon}.$$

This together with Proposition \[theorem:stabilityestimate\] imply the statement of the Corollary.

We note that all the above results can be adapted to three dimensions. The proof follows exactly the same steps by considering the natural extension of the definition of the discrepancy function $E$ in three dimensions.

The assumption made in Corollary \[corolar-1\] that $\alpha_0\approx\delta$ for $d_c=diam(D_c)$ and $dist(\partial D_c,\partial D_a)=d$ small enough is based on of Lemma \[lemma:phif1constant\] and on the numerical results presented in Section \[sec:numerics\]. In particular, for the same settings as in Figure \[fig:nonincreasing\_F\], Figure \[fig:besselFS\_alpha0\_vary\_k\_mu\] shows that given small $d_c$ for small enough $d$, we have roughly that $10^{-2.5} \leq \alpha_{0} \leq 10^{-1}$. Since $\delta = 2\cdot10^{-2}$ we have in this case that $\frac{1}{2\sqrt{5}} \delta \lesssim \alpha_{0} \lesssim 5 \delta$.

\[constant\] As suggested by our numerical results in Section \[sec:numerics\] we beleive that the constants (denoted by $C$) in Proposition \[theorem:stabilityestimate\], Lemma \[teorema-int\] and Corollary \[corolar-1\] will only have small values for $d_c=diam(D_c)$ and $dist(\partial D_c,\partial D_a)=d$ small enough.

Numerics {#sec:numerics}
========

In this section we proceed with the numerical study of the minimal norm solution for obtained through Tikhonov regularization with the Morozov discrepancy principle for the choice of the regularization parameter in two dimensions. First we focus on the general setup of our numerical approach, and then in Section \[subsec:parametersetup\] we discuss more specifically the parameters used in our numerical examples. In Sections \[subsec:nearfieldstability\] and \[subsec:sphericalsource\] we present numerical data which demonstrates how stably $\phi$ depends on $f$ and various control statistics for a spherical point source. All figures generally display their respective parameters in an offset legend.

For all of the numerical computations done in this section, we discretize the integral operator $K$ via the method of moment collocation. We refer to ([@Kress-Book99], §17.4) for more details on the method. First we choose an approximate basis of functions for $L^{2}(\partial D_{a})$. To do this we suppose the domain $D_{a}$ can be parametrized in polar coordinates by points $$(s(\tau)\cos{\tau}, \, s(\tau)\sin{\tau})), \quad \tau \in [0,2\pi],$$ where $s:\mathbb{R} \to \mathbb{R}_{+}$ is a $2\pi$-periodic smooth function. Using these coordinates, any function $\phi$ defined on $\partial D_{a}$ can be realized via the pullback as a function of $\tau$: $$\phi(s(\tau)\cos{\tau}, \, s(\tau)\sin{\tau}).$$ For convenience, let us use the notation $\widehat{\tau} = (\cos{\tau},\sin{\tau})$ and $\widehat{\tau}^{\perp} = (-\sin{\tau}, \cos{\tau})$.

Now let $n_{a} \in \mathbb{N}$ and let $\tau_{j} = \frac{2\pi j}{n_{a}}$, $0 \leq j \leq n_{a}-1$ be $n_{a}$ equally spaced points on the interval $[0,2\pi)$. We then use the exponential basis functions $\{e^{il\tau}\}_{l=0}^{n_{a}-1}$ for $L^{2}([0,2\pi])$ and approximate a given $\phi \in L^{2}(\partial D_{a})$ via interpolation at the points $\{\tau_{j}\}_{j=0}^{n_{a}-1} \subset [0,2\pi]$. Note that $$\begin{aligned}
\int_{\partial D_{a}}\phi(\mathbf{y}) \frac{\partial \Phi}{\partial \nu_{\mathbf{y}}}(\mathbf{x}, \, \mathbf{y})\,dS_{\mathbf{y}} & = \int_{0}^{2\pi}\phi(s(\tau)\cos{\tau}, \, s(\tau)\sin{\tau})\frac{\partial \Phi}{\partial \nu_{y}}(\mathbf{x}, (s(\tau)\cos{\tau}, \, s(\tau)\sin{\tau})) \notag\\
& \quad \cdot \sqrt{s(\tau)^{2} + s'(\tau)^{2}}\,d\theta.\end{aligned}$$ Furthermore, since $\left( s'(\tau)\cos{\tau} - s(\tau)\sin{\tau}, \, s(\tau)\cos{\tau} + s'(\tau)\sin{\tau}\right)$ is a tangent vector to $\partial D_{a}$, we have that $$\begin{aligned}
\nu(\mathbf{y}) = \nu(\tau) & = \frac{(s(\tau)\cos{\tau} + s'(\tau)\sin{\tau}, \, s(\tau)\sin{\tau} - s'(\tau)\cos{\tau})}{\sqrt{ s(\tau)^{2} + s'(\tau)^{2}}}\\
& = \frac{s(\tau)\widehat{\tau} - s'(\tau)\widehat{\tau}^{\perp}}{\sqrt{ s(\tau)^{2} + s'(\tau)^{2}}}\end{aligned}$$ is the unit outward normal vector to $\partial D_{a}$. It is then straightforward to compute in the case of the Helmholtz equation in 2D that $$\begin{aligned}
& \quad \frac{\partial \Phi}{\partial \nu_{\mathbf{y}}}(\mathbf{x}, \, (s(\tau)\cos{\tau}, \, s(\tau)\sin{\tau}))\\
& = \nabla_{y}\Phi(\mathbf{x}, (s(\tau)\cos{\tau}, \, s(\tau)\sin{\tau})) \cdot \nu(\tau)\\
& = \frac{ik}{4}H_{0}^{(1)'}(k|\mathbf{x} - s(\tau)\widehat{\tau}|)\frac{s(\tau)\widehat{\tau} - \mathbf{x}}{\sqrt{s(\tau)^{2} + |\mathbf{x}|^{2} - 2s(\tau)\mathbf{x}\cdot \widehat{\tau}}} \cdot \frac{s(\tau)\widehat{\tau} - s'(\tau)\widehat{\tau}^{\perp}}{\sqrt{s(\tau)^{2} + s'(\tau)^{2}}}.\end{aligned}$$

Let $n_{a} \in \mathbb{N}$ be the number of sample points on the antenna, $\partial D_{a}$, and let $n_{c} \in \mathbb{N}$ be the total number of sample points on $\partial D_{c}$. Also let $n_{R}$ be the total number of sample points on $\partial B_{R}$. We write the $2 \times (n_{c} + n_{R})$ matrix of points $$\mathbf{X} := [x_{1}, \ldots, x_{n_{c}+n_{R}}],$$ where each $x_{j}$ is a $2$-vector, $\{x_{j}\}_{j=1}^{n_{c}} \subset \partial D_{c}$ and $\{x_{j}\}_{j=n_{c}+1}^{n_{c} + n_{R}} \subset \partial B_{R}$. Approximations of all the integrals involved are then computed using a standard left endpoint sum with the appropriate quadrature weights. All the numerical examples presented herein take $D_{c}$ to be an annular sector parametrized by $r \in [r_{1}, r_{2}]$ and $\theta \in [\theta_{1}, \theta_{2}]$. See Figure \[fig:numericssetup\] for details.

For each $1 \leq j \leq n_{c} + n_{R}$ and each $0 \leq l \leq n_{a}-1$, we compute $K[e^{il\tau}](x_{j})$ via the approximation $$\frac{2\pi}{n_{a}}\sum_{m=0}^{n_{a}-1} \frac{\partial \Phi( x_{j}, [s(\tau_{m})\cos{\left ( \tau_{m}\right )}, \, s(\tau_{m})\sin{\tau_{m} }]^{T})}{\partial \nu_{\mathbf{y}}} e^{il \tau_{m}}\sqrt{s(\tau_{m})^{2} + s'(\tau_{m})^{2}}.$$ If we fix $j$ and vary $l$, we see that the above sum is equivalent to computing the discrete Fourier transform of the $n_{a}$-vector $$\mathbf{v}_{j} := \left[\frac{\partial \Phi( x_{j}, [s(\tau_{m})\cos{\tau_{m}}, \, s(\tau_{m})\sin{\tau_{m} }]^{T})}{\partial \nu_{\mathbf{y}}}\sqrt{s(\tau_{m})^{2} + s'(\tau_{m})^{2}}\right]_{m=0}^{n_{a}-1}, \label{eq:vj}$$ which can be computed efficiently using the Fast Fourier Transform algorithm. In particular, for the vector $\mathbf{v}_{j}$ in $$[K\{e^{il\tau}\}(x_{j})]_{l=0}^{n_{a}-1} \approx 2\pi \mathsf{FFT}(\mathbf{v}_{j}),$$ where $\mathsf{FFT}$ is defined on $n_{a}$-vectors $\mathbf{w} = [w_{1},\ldots, w_{n_{a}}]^{T} \in \mathbb{C}^{n_{a}}$ by $$\mathsf{FFT}(\mathbf{w}) = \left[ \frac{1}{n_{a}}\sum_{j=1}^{n_{a}}w_{j}e^{\frac{2\pi i (j-1)(l-1)}{n_{a}}} \right]_{l=1}^{n_{a}}.
\label{eq:fft}$$ So the matrix representation of $K$ is then the $n_{a} \times (n_{c} + n_{R})$ matrix $$A := 2\pi [\mathsf{FFT}(\mathbf{v}_{1}), \cdots, \mathsf{FFT}(\mathbf{v}_{n_{c}+n_{R}})].$$

Now, in order to approximately solve the ill-posed problem $K\phi = (f_{1}, f_{2})$, we attempt to solve the linear system $$\begin{aligned}
K_{1}\phi(x_{j}) & = f_{1}(x_{j}), \quad 1 \leq j \leq n_{c}\\
K_{2}\phi(x_{j}) & = f_{2}(x_{j}), \quad n_c + 1 \leq j \leq n_c + n_{R}.\end{aligned}$$ Since $A$ is computed with respect to the functions $e^{il\theta}$, we first consider the approximate coefficients of $\phi$ with respect to the finite basis $\{e^{il\tau}\}_{l=0}^{n_{a}-1}$, given by $$c_{l} := \frac{1}{n_{a}}\sum_{m=0}^{n_{a}-1}e^{-il\tau_{m}}\phi(s(\tau_{m})\cos(\tau_{m}), \, s(\tau_{m})\sin(\tau_{m}))\,d\tau \approx \frac{1}{2\pi}\int_{0}^{2\pi} e^{-il\tau}\phi(s(\tau)\cos(\tau), \, s(\tau)\sin(\tau))\,d\tau.$$ Let $$h = [c_{0}, c_{1}, \ldots, c_{n_{a} - 1}]^{T} \in \mathbb{C}^{n_{a}}.$$

We then numerically compute the Tikhonov regularized solution $$h_{\alpha} := (A^{*}A + \alpha I)^{-1}A^{*}f,$$ with $\alpha > 0$. The solution vector $h_{\alpha}$ yields the approximate coefficients $c_{l}$ of the desired density $\phi$ with respect to the functions $\{e^{il\tau}\}_{l=0}^{n_{a}-1}$. We obtain the density $\phi_{\alpha}$ corresponding to $h_{\alpha}$ sampled at the angles $\tau_{m}$ on $\partial D_{a}$ by the formula $$\phi_{\alpha}(\tau_{m}) := \sum_{l=0}^{n_{a}-1}[h_{\alpha}]_{l}e^{il\tau_{m}}.$$

After computing the residual $K\phi - f$ (e.g. for $\phi = \phi_{\alpha}$), we will then need to compute $$E(\phi_\alpha, f) = \frac{1}{\|f_{1}\|_{L^{2}(\partial D_{c})}^{2}}\| K_{1}\phi  - f_{1}\|_{L^{2}(\partial D_{c})}^{2} + \frac{1}{2\pi R} \|K_{2}\phi - f_{2}\|_{L^{2}(\partial B_{R})}^{2}.$$ Recall that the discrepancy function $F(\alpha)$ was defined by $$F(\alpha) = E(\phi_\alpha, f) - \delta^{2}, \label{eq:F}$$ where $\delta > 0$ is a fixed error parameter. As discussed in Section \[sec:stability\], the mapping $$\alpha \mapsto E(\phi_\alpha, f)$$ is not globally increasing, as can be numerically demonstrated. However, for certain feasible regions of $(\alpha, \epsilon)$, $F$ is increasing. And in this case, there is a unique $\alpha_{\delta}$ such that $F(\alpha_{\delta}) = 0$. To find $\alpha_{\delta}$, we use Newton’s method combined with an initial coarse line search to identify a good starting point.

First note that if we split the matrix $A$ into two blocks $A_{near}$ ($n_{c}$ by $n_{a}$) and $A_{far}$ ($n_{R}$ by $n_{a}$) so that $$A = \left[ \begin{array}{c}
A_{near}\\
A_{far}
\end{array}\right],$$ then $[A\phi]_{1} = A_{near}\phi$, $[A\phi]_{2} = A_{far}\phi$, and $A^{*}A = A_{near}^{*}A_{near} + A_{far}^{*}A_{far}$. In the discretized setting, instead of (\[eq:F1\]) we take $$F(\alpha) =\frac{1}{\|f_{1}\|^{2}} \|A_{near}h_{\alpha} - f_{1}\|_{L^{2}(\partial D_{c})}^{2} + \frac{1}{2\pi R}\|A_{far}h_{\alpha} - f_{2}\|_{L^{2}(\partial B_{R})}^{2} - \delta^{2}$$ with $$h_{\alpha} = (A^{*}A + \alpha I)^{-1}A^{*}f = (A_{near}^{*}A_{near} + A_{far}^{*}A_{far} + \alpha I)^{-1}{\left ( A_{near}^{*}f_{1} + A_{far}^{*}f_{2}\right )}.$$ Then in the same spirit as that presented in [@CoKr-Book98] for Tikhonov regularization with respect to the standard $L^2$ norm, we compute $$\begin{aligned}
F'(\alpha) & = \frac{-2\alpha}{\|f_{1}\|_{L^{2}(\partial D_{c})}^{2}}{\textrm{Re}}{\left ( \frac{\partial h_{\alpha}}{\partial \alpha}, \, h_{\alpha}\right )} \notag\\
& \quad + {\left (  \frac{1}{\pi R} - \frac{2}{\|f_{1}\|_{L^{2}(\partial D_{c})}^{2}}\right )}{\textrm{Re}}{\left ( \frac{\partial h_{\alpha}}{\partial \alpha}, \, A_{far}^{*}A_{far}h_{\alpha} - A_{far}^{*}f_{2}\right )} \label{eq:dF1}\\
\frac{\partial h_{\alpha}}{\partial \alpha} & = -(A^{*}A + \alpha I)^{-1}h_{\alpha}, \label{eq:dF2}\end{aligned}$$ where $( \cdot, \cdot )$ denotes the $L^{2}$ inner product on $\partial D_{a}$.

The function $f_{1}$ defined on $\partial D_{c}$ could be, for example, the trace of a plane wave, or of the fundamental solution to the Helmholtz equation based at some fixed point $\mathbf{x}_{0}$, i.e. a point source. For this paper, we focus on the case where $f_{1}$ is a point source and where $f_{2} \equiv 0$ on $\partial B_{R}$. A spherical point source is represented as $$\frac{i}{4}H_{0}^{(1)}(k|\mathbf{x} - \mathbf{x}_{0}|), \label{eq:sphericalsource}$$ where $\mathbf{x}_{0}$ is the source point (typically outside of $B_{R}$).

For such an $f_{1}$, there are some quantities in which we will be interested so as to determine the effectiveness of a given density $\phi$ in solving the problem $K\phi =f$. These are: the relative error of $K\phi$ on $\partial D_{c}$; the $L^{2}$ average of $K\phi$ on $\partial B_{R}$; the relative and absolute stability of $\phi$ when applying a small perturbation to $f_{1}$; the norm of $\phi$ on $\partial D_{a}$. In other words, we will measure $$\frac{\|K_{1} \phi - f_{1}\|_{L^{2}(\partial D_{c})}}{\|f_{1}\|_{L^{2}(\partial D_{c})}}, \quad \frac{1}{\sqrt{2\pi R}} \|K_{2}\phi\|_{L^{2}(\partial B_{R})}, \label{eq:computevar1}$$ $$\frac{\|\phi_{\alpha_{\epsilon}}- \phi_{\alpha_0}\|_{L^{2}(\partial D_{a})}}{\|\phi_{\alpha_0}\|_{L^{2}(\partial D_{a})}}, \quad \|\phi_{\alpha_\epsilon} - \phi_{\alpha_0}\|_{L^{2}(\partial D_{a})}, \label{eq:computevar2}$$ and $$\|\phi\|_{L^{2}(\partial D_{a})}, \label{eq:phinorm}$$ where $\phi_{\alpha_\epsilon}$ is the Tikhonov regularization solution to $K\phi = (f_{1,\epsilon},0)$ with $\|f_{1} - f_{1,\epsilon}\|_{L^{2}(\partial D_{c})} = \epsilon \|f_{1}\|_{L^{2}(\partial D_{c})}$, and $\phi_{\alpha_0}$ is the solution with unperturbed $f_{1}$. The Morozov solution $\alpha_0$ and $\alpha_\epsilon$ are computed via Newton’s Method using and such that $$\begin{aligned}
E(\phi_{\alpha_0}, f)  & = \delta^2 \notag\\
E(\phi_{\alpha_\epsilon}, f_{\epsilon}) & = \delta^2.\end{aligned}$$ See also the discussion following . Recall from (\[random-f\]), that when adding noise to the data $(f_{1}, 0)$, we choose a random perturbation $\eta \in L^{2}(\partial D_{c})$ and set $$f_{1,\epsilon} = f_{1} + \epsilon \widehat{\eta} \|f_{1}\|_{L^{2}(\partial D_{c})}, \label{eq:noiseterm}$$ where $\epsilon > 0$ represents the relative percentage of noise added. In the discrete case, the noise is chosen to be a complex $n_c$-vector $\nu$ whose real and imaginary components are pseudorandom numbers (we used uniformly distributed noise, but any distribution would yield similar results) on the interval $(-1,1)$. Furthermore, for reproducibility, whenever generating $\nu$ using a pseudorandom number generator, we always reset the seed to the same value.

Parameters Used for Numerical Experiments\[subsec:parametersetup\]
------------------------------------------------------------------

Here we describe some of the parameters used for the various numerical experiments presented. In Section \[subsec:sphericalsource\] we always assume that $\partial D_{a}$ is a circle with radius given by $a = 0.01$, and that $\partial D_{c}$ is a sector of an annulus with $\theta_{1} = 3\pi/4$ and $\theta_{2} = 5\pi/4$. We also take $R = 10$ in all computational examples. We remark that we always restrict the distance from $D_{c}$ to $D_{a}$ to be no smaller than $10^{-3}$ due to the numerical limitations of our approach. This is due to the fact that $K\phi$ is a singular integral when evaluating at points very near to $\partial D_{a}$. Therefore, at points on $\partial D_{c}$ that are near $\partial D_{a}$ the limitations of machine precision become more and more apparent. Numerically, we observed that our direct approach starts to break down near $d = \mathrm{dist}(\partial D_{c}, \partial D_{a}) = 10^{-4}$. However, we stress that one could most likely obtain high accuracy in computing $K\phi$ for $d \leq 10^{-4}$ by using the Nyström method as discussed in [@Kress-Book99].

For the collocation method, we use $n_{a} = 256$ sample points on $\partial D_{a}$, and $n_{\mathrm{arc}_{1}} = 256$ points on the inner arc of $\partial D_{c}$, with the remaining points chosen so as to keep the quadrature weights approximately constant. Thus for a very thin region, $n_{c} \approx 512$. We also take $n_{R} = 256$ (number of sample points on $\partial B_{R}$). Note that increasing $n_{c}$ or $n_{R}$ will put more emphasis on matching $f$ on $\partial D_{c}$ or $\partial B_{R}$, respectively. The discrepancy parameter $\delta$ used for Tikhonov regularization will typically be fixed at $0.02$. The key variables under consideration are $d = r_{1} - a$, $k$, and $\epsilon$ (perturbation parameter for adding noise to $f_{1}$). All of the plots presented in the following sections involve varying two of the aforementioned parameters and plotting different quantities of interest, as stated in (\[eq:computevar1\])-(\[eq:phinorm\]).

When evaluating the relative change in $\phi$ given a perturbation of $f_{1}$, denoted by $f_{\epsilon, 1}$, we remark that for the parameter choices we used, a $0.5\%$ change ($\epsilon = 0.005$) in $f_{1}$ yielded a roughly $5\%$ change in $\phi$. However, one must keep in mind that this depends quite a lot on the parameters used. In particular, setting the discrepancy $\delta = 0.02$ in all the simulations had an important effect on the numerical results. If we had used $\delta = 0.05$ instead (which leads to approximately a $5 \%$ mismatch on the region $\partial D_{c}$), then the relative change in $\phi$ given $\epsilon = 0.005$ would be noticeably smaller. So ultimately there is a tradeoff between $\mu$, $R$, $k$, $\delta$, and $\epsilon$ which is not entirely trivial, but this still can be examined experimentally as we have done.

Near field stability\[subsec:nearfieldstability\]
-------------------------------------------------

We present below Figure \[fig:singularvalues\], which shows how the first $50$ singular values of the operator $K = (K_{1}, K_{2})$ vary with $d$. It is clear that for $d$ small (i.e. for control in the nearfield of the antenna), the rate of decay of the singular values of $K$ is considerably slower than for larger $d$. This in turn provides some experimental evidence for the fact that nearfield control seems to be more feasible in terms of stable dependence of the solution $\phi$ on $f$. We also show Figure \[fig:singularvaluesdiff\], which shows the behavior of the first and sixth singular value of $K$ with respect to $d$ and $k$.

![Plot of first $50$ singular values of $K:L^{2}(\partial D_{a}) \to \Xi$ for $\partial D_{a}$ a circular antenna of radius $a = 0.01$ and $\partial D_{c}$ an annular region of varying distance from $\partial D_{a}$.[]{data-label="fig:singularvalues"}](50_singular_values){width="\textwidth"}

![Plot of first singular value of $K:L^{2}(\partial D_{a}) \to \Xi$ as well as the relative difference of the first and sixth singular values with respect to $d$ and $k$. Again, $\partial D_{a}$ is a circular antenna of radius $a = 0.01$ and $\partial D_{c}$ an annular region.[]{data-label="fig:singularvaluesdiff"}](singular_values_diff){width="\textwidth"}

Control for a Spherical Point Source\[subsec:sphericalsource\]
--------------------------------------------------------------

We now consider the case that $$f_{1}(\mathbf{x}) = \frac{i}{4}H_{0}^{(1)}(k|\mathbf{x} - \mathbf{x}_{0}|),$$ where $\mathbf{x}_{0}$ is the source point. In all examples presented in this section, we have $R = 10$ unless otherwise noted, and $\mathbf{x}_{0} = [20, 0]^{T}$ or $\mathbf{x}_{0} = [10000, 0]^{T}$ (to approximate a source at infinity).

First we observe how the frequency $k$ and distance $d$ from $\partial D_{c}$ to $\partial D_{a}$ affects the various control criteria. In figures \[fig:besselFS\_epsilon\_0005\_delta\_002\_vary\_k\_01\_to\_100\_mu\_0001\_to\_003\_data\] and \[fig:besselNS\_epsilon\_0005\_delta\_002\_vary\_k\_01\_to\_100\_mu\_0001\_to\_003\_data\] we vary $k$ from $0.1$ to $100$ and $d$ from $0.001$ to $0.003$ with $a = 0.01$. With the error discrepancy set at $\delta = 0.02$, we have in both figures that relative error on $\partial D_{c}$ is very close to $2\%$ for all data points. Moreover, with $0.5\%$ noise added to $f$, roughly a $5\%$ change in $\phi$ is observed at all frequencies when $d$ is near its lower limit. A bit more sensitivity is observed for frequencies $k < 20$ when $d$ increases beyond $0.01$. Interestingly, for $k > 20$ the optimal parameter $\alpha$ is larger and corresponding power budget smaller in order to achieve discrepancy $\delta$.

In figures \[fig:besselFS\_epsilon\_0005\_delta\_002\_vary\_k\_01\_to\_100\_mu\_0001\_to\_003\_data\] and \[fig:besselNS\_epsilon\_0005\_delta\_002\_vary\_k\_01\_to\_100\_mu\_0001\_to\_003\_data\] it is clear that for smaller $d$ values the sensitivity of $\phi$ to $0.5\%$ noise added to $f$ is close to $5\%$. As $d$ increases, sensitivity of $\phi$ to noise increases as expected. Having $\mathbf{x}_{0}$ nearer or farther from $\partial B_{R}$ does not have a very significant effect on the overall shape of each subplot.

Figures \[fig:besselFS\_vary\_mu\_epsilon\] and \[fig:besselNS\_vary\_mu\_epsilon\] show how the quantities of interest change with $d$ and the noise factor $\epsilon$, both with $k = 10$. The reason for choosing $k=10$ instead of, e.g., $k=1$ is that from figure \[fig:besselFS\_epsilon\_0005\_delta\_002\_vary\_k\_01\_to\_100\_mu\_0001\_to\_003\_data\] we see a slightly higher sensitivity of $\phi$ to noise for approximately $1 < k < 20$ when $d$ starts to increase. So the goal was to capture the worst case scenario for the control stability. For smaller values of $d$ we see as before that a roughly $0.5\%$ change in $f_{1}$ yields about a $5\%$ change in $\phi$. Moreover, the dependence on $\epsilon$ for fixed $d$ is superlinear, consistent with the illposedness of the problem. Interestingly, sensitivity of $\phi$ at $d \approx 0.015$ is better than at nearby values, but of course such a value depends on the other parameters of the problem setup.

Finally, we consider Figure \[fig:besselFS\_vary\_k\_R\], which shows the dependence on $R$ and $k$ for a source at $\mathbf{x}_{0} = [10000,\, 0]^{T}$. Overall, one can see that $R$ can be decreased to around $R = 3$ at any frequency between $0.1$ and $100$ and still achieve the same approximate level of sensitivity for $\phi$ as in the previous plots with $R = 10$.

Conclusions and Future Work {#sec:conclusions}
===========================

In this paper we studied the feasibility of the active control scheme for the scalar Helmholtz equation. In the $L^2$ setting, we presented analytic conditional stability results as well as detailed numerical sensitivity studies for the minimal energy solution. We provided several analytic and numerical arguments for the scheme’s feasibility and broadband character in the near field when the interrogating field is a far field of a far field observer.

We focused our discussion in this paper only on the case of an interrogating far field point source (i.e. similar to a plane wave with corresponding decay) because we believe that this situation is relevant in usual radar or sonar detection problems. In contrast, the case of an interrogating plane wave corresponds to a different problem, where the observer is close to the source and control region and thus the interrogating signal does not have sufficient decay.

In fact, we have numerically studied the case when the interrogating field is a plane wave without decay or a given uniform field. We observed the scheme does not behave well for the uniform field and that although the stability and accuracy of the near field scheme are essentially independent of the plane wave direction, the overall performance of the scheme is not as good when compared to the case of an interrogating signal coming from a far field observer presented above. In fact, for the same settings as in Figure \[fig:nonincreasing\_F\] when comparing the case of an interrogating far field point source with an interrogating plane wave, we obtained $5\%$ versus $8\%$ stability error and power budget levels of $\approx 10^{-1}$ versus $\approx 10$. We conclude that the scheme performance depends not only on the location of the control region with respect to the source region but also on the amplitude and oscillatory pattern of the incoming field.

Currently we are considering a more localized basis for $L^2(\partial D_a)$ (e.g. delta function basis, or splines) in order to better observe the field characteristics in the control region $D_c$ and around the antenna $D_a$. We also plan to study the active control scheme for linear arrays and for large elongated antennas. Then, as a next step in our research efforts, we will work on the extension of the current numerical sensitivity study to three dimensions and full Maxwell system and on the study of near field control with planar and conformal arrays.
