# SPHERICAL CAPS DO NOT ALWAYS MAXIMIZE NEUMANN EIGENVALUES ON THE SPHERE

DORIN BUCUR, RICHARD S. LAUGESEN, ELOI MARTINET, AND MICKAËL NAHON

**ABSTRACT.** We prove the existence of an open set $\Omega \subset \mathbb{S}^2$ for which the first positive eigenvalue of the Laplacian with Neumann boundary condition exceeds that of the geodesic disk having the same area. This example holds for large areas and contrasts with results by Bandle and later authors proving maximality of the disk under additional topological or geometric conditions, thereby revealing such conditions to be necessary.

## 1. INTRODUCTION

On a smooth subdomain $\Omega$ of the standard unit sphere $\mathbb{S}^2$ in $\mathbb{R}^3$, the Laplace-Beltrami operator with Neumann boundary conditions has discrete spectrum consisting of eigenvalues, denoted

$$0 = \mu_0(\Omega) < \mu_1(\Omega) \le \dots \le \mu_k(\Omega) \to +\infty$$

and counted with respect to multiplicity. The eigenfunction corresponding to $\mu_0 = 0$ is constant. The first and second positive eigenvalues are

$$\mu_1(\Omega) = \min_{u \in H^1(\Omega) \setminus \{0\}, \int_{\Omega} u = 0} \frac{\int_{\Omega} |\nabla u|^2}{\int_{\Omega} u^2}, \quad \mu_2(\Omega) = \min_{S \in \mathcal{S}_2} \max_{u \in S \setminus \{0\}} \frac{\int_{\Omega} |\nabla u|^2}{\int_{\Omega} u^2},$$

where integrals are with respect to surface area measure and $\mathcal{S}_2$ is the family of 2-dimensional subspaces of $\{u \in H^1(\Omega) : \int_{\Omega} u = 0\}$. The eigenfunctions $(u_k)_{k \ge 0}$ satisfy

$$\begin{cases} -\Delta_{\mathbb{S}^2} u_k = \mu_k(\Omega) u_k & \text{in } \Omega, \\ \partial_{\nu} u_k = 0 & \text{on } \partial \Omega, \end{cases}$$

where $\partial_{\nu}$ is the inward normal derivative.

The question of maximizing $\mu_1(\Omega)$ under an area constraint $|\Omega| = m$ for given “mass” $m \in (0, 4\pi)$ has been studied intensively in the last fifty years, and yet intriguing open problems remain unsolved.

To briefly summarize this history, Ashbaugh and Benguria [2] proved that if $m \le 2\pi$ and $\Omega$ is contained in a hemisphere then

$$\mu_1(\Omega) \le \mu_1(B_{\theta_m}) \quad (1)$$

where $B_{\theta_m}$ is a spherical cap of area $m$ and $\theta_m$ is the corresponding aperture angle or geodesic radius, measured from the center of the cap. Whether this hemisphere inclusion condition is merely technical or instead indicates a genuine mathematical obstruction has remained an open question, stated for example in [10, Conjecture 1.3]. A step toward understanding this point was taken by Bucur, Martinet and Nahon [4], who relaxed the hemisphere condition to inclusion in the complement of a spherical cap of mass $m$. The smaller $m$ is, the larger the region in which $\Omega$ may be placed. However, in the same paper, the authors numerically exhibited a density or weight whose first eigenvalue is larger than that of the spherical cap having the same mass, thus leading to the conclusion that even if inequality (1) holds true without the inclusion restriction, the way to prove it cannot follow existing mass transplantation arguments. Further, for domains having area somewhat greater than $2\pi$, Martinet [12] found an example whose first eigenvalue, computed numerically, exceeds that of the spherical cap having the same mass. This counterexample domain is not simply connected and indeed has rather complicated topology.

<footnote n="2020">2020 Mathematics Subject Classification. Primary 35P15. Secondary 58J50.</footnote>

<footnote>Key words and phrases. Vibrating free membrane, isoperimetric inequality, Laplacian eigenfunction, Laplace-Beltrami.</footnote>
<figure><img src="image_1.png" /><figcaption>FIGURE 1. Different views of one punctured domain used in the proof of Theorem 1.</figcaption></figure>

Results by Bandle [3] preceded those of Ashbaugh and Benguria, although with different hypotheses. She proved that inequality (1) holds provided one restricts to simply connected domains $\Omega$ having at most half the area of the sphere, $|\Omega| \le 0.50|\mathbb{S}^2|$. This approach imposes no restriction on the parts of the sphere the set is permitted to occupy. Langford and Laugesen [10] recently improved from 50% to 94% of the sphere, proving (1) for simply connected domains with area $|\Omega| \le 0.94|\mathbb{S}^2|$.

The main purpose of this paper is to prove analytically that inequality (1) can fail for multiply connected domains, by constructing a rigorous counterexample whose topology is relatively simple. Precisely, for a large mass $m$, e.g. above $0.81|\mathbb{S}^2|$, we prove that there exists a non-simply connected set of mass $m$ whose first eigenvalue is larger than that of the corresponding spherical cap.

Write $\Theta$ for the unique aperture defined by the relation $\mu_1(B_\Theta) \sin^2 \Theta = 1$; it is shown in [10, Propositions 3.1 and 4.2] that $\mu_1(B_\Theta) \sin^2 \theta > 1$ on $(0, \Theta)$ and $\mu_1(B_\Theta) \sin^2 \theta < 1$ on $(\Theta, \pi)$ and that the value of $\Theta$ is approximately $0.7\pi$:

$$0.70\pi < \Theta < 0.71\pi.$$

The counterexample to (1) is provided by:

**Theorem 1.** There exists $\delta > 0$ such that for each $\theta \in (\Theta - \delta, \pi)$, there exists an open set $\Omega_\theta \subseteq \mathbb{S}^2$ with $|\Omega_\theta| = |B_\theta|$ and

$$\mu_1(\Omega_\theta) > \mu_1(B_\theta).$$

By no means do we claim optimality of the lower bound $\Theta - \delta$. On the contrary, we provide numerical evidence supporting the idea that such a counterexample can exist for lower areas, although still above $2\pi$. Apertures in the range $(\Theta, \pi)$ correspond to the area or mass $m$ lying between about $0.81|\mathbb{S}^2|$ and $|\mathbb{S}^2|$.

The key point underlying the theorem is that when $\theta > \Theta$, the geodesic ball $B_\theta$ possesses hot spots, meaning interior local extrema of the first eigenfunctions. The strategy for proving the theorem is to start with such a spherical cap of aperture $\theta$ and remove four small ellipses of large eccentricity that are symmetrically placed at the latitude of the hot spots (see Figure 1). The theorem then follows from an asymptotic analysis of the first (simple) eigenvalue of this modified set, building on work of Ozawa [15], however with the extra difficulty of dealing with a geometry where the first eigenvalue is multiple.

A different kind of counterexample is pursued in Section 5, where we provide a non-rigorous argument for “helmet” or “basket” counterexamples, as depicted in Figure 2. We support our construction also with numerical evidence. The approach suggests that the lower bound $\Theta - \delta$ in Theorem 1 might be significantly lowered.

**Open problems.** Some significant problems remain to be resolved.

* When $\theta \le \frac{1}{2}\pi$, is it true that the first positive eigenvalue $\mu_1$ is maximal on the spherical cap $B_\theta$ among domains of measure $|B_\theta|$, with no further restriction on the domain?

* Among simply connected domains, is there a counterexample when $|B_\theta| > 0.94|\mathbb{S}^2|$?
**Overview.** The paper is structured as follows. In the next section we recall known facts about the first and second eigenfunctions of spherical caps. In Section 3 we give a precise asymptotic for Neumann eigenvalues in the presence of several holes, following in the footsteps of Ozawa [15]. Section 4 estimates the second term of the asymptotic development of $\mu_1$ on the perturbed spherical cap, and hence proves the main theorem. In Section 5 we present the non-rigorous argument for helmet counterexamples. Section 6 proves that our rigorous method cannot provide counterexamples for simply connected sets having mass close to $4\pi$. This is done by comparing the first term asymptotic of the whole sphere minus a single hole with the capacity of that hole. Finally, Section 7 establishes various technical estimates used in proving the main theorem.

## 2. STRUCTURE OF THE FIRST AND SECOND EIGENFUNCTIONS ON THE SPHERE

Recall that $\mu_0 = 0$. The separated form of the eigenspace for $\mu_1(\theta)(:= \mu_1(B_\theta))$ of the spherical cap $B_\theta$ is given by the next proposition, which is drawn from [10, Proposition 4.1]. Parts (a) and (b) are distinguished by whether the aperture is smaller or larger than the value $\Theta \simeq 0.7\pi$ that was specified earlier.

Throughout the paper we adopt spherical coordinates, parametrizing the sphere $S^2$ by

$$ (t, \phi) \in [0, \pi] \times \mathbb{R}/2\pi\mathbb{Z} \mapsto (\sin t \cos \phi, \sin t \sin \phi, \cos t). $$

**Proposition 2.** If $0 < \theta < \pi$ then the eigenspace of $-\Delta_{S^2} u = \mu_1(\theta)u$ on the spherical cap $B_\theta$ with Neumann boundary condition is spanned by two functions of the form

$$ g(t) \cos \phi \quad \text{and} \quad g(t) \sin \phi, $$

where the radial part $g$ is smooth and satisfies the Neumann boundary condition $g'(\theta) = 0$, the differential equation

$$ -\frac{1}{\sin t} ((\sin t)g'(t))' + \frac{1}{\sin^2 t} g(t) = \mu_1(B_\theta)g(t) \quad (2) $$

and has $0 = g(0) < g(t)$ whenever $t \in (0, \theta]$.

(a) If $0 < \theta \le \Theta$ then $g(t)$ is strictly increasing for $t \in (0, \theta)$, with $g' > 0$ there.

(b) If $\Theta < \theta < \pi$ then $g(t)$ first strictly increases and then strictly decreases; more precisely, a maximum point $t_{max}(\theta) \in (\frac{\pi}{2}, \Theta)$ exists such that $g' > 0$ on $(0, t_{max}(\theta))$ and $g' < 0$ on $(t_{max}(\theta), \theta)$.

The proposition is illustrated in [10, Figure 2], which we recommend consulting.

## 3. MULTIPLE NEUMANN EIGENVALUES ON PERTURBED DOMAINS

In this section we develop some convergence results for eigenfunctions and multiple eigenvalues of the Neumann Laplacian on singularly perturbed domains, specifically, on domains enjoying simultaneous perforation with several small holes and boundary variation. The goal of this section is to establish Proposition 3, which gives a precise asymptotic for the eigenvalues in such a situation.

The literature is rich in these type of estimates. For instance, we refer to [13, 14, 6, 11, 9] and references therein. However, for the convenience of the reader and since we did not find a reference where all our specific needs are met, we prefer to produce the computation in a self-contained manner, based on the approach of Nazarov and Sokolowski [14].

**3.1. Statement of the result.** Throughout this section, we fix $M$ to be a compact closed 2-dimensional Riemannian manifold and take an open subset $\Omega \subseteq M$ with smooth boundary (which allows also for the case $\partial\Omega = \emptyset$). Consistent with the notation used in the introduction for subdomains of the sphere, we write $\mu_i(\Omega)$ for the Neumann eigenvalues of the Laplace-Beltrami operator on $\Omega$, starting with $\mu_0(\Omega) = 0$. Let $(u_i)_{i \in \mathbb{N}}$ be a corresponding choice of $L^2$-normalized eigenfunctions. Take:

* $N(\ge 1)$ distinct points $O_1, \dots, O_N \in \Omega$;
• for each $n = 1, \dots, N$, a conformal local chart $\psi_n$ from a neighborhood of 0 (say $B_r$ for some small radius $r > 0$) in $\mathbb{R}^2$ into $\Omega$ with $\psi_n(0) = O_n$. We suppose $r$ is small enough that the sets $(\psi_n(B_r))_{n=1}^N$ are disjoint to each other and do not intersect the boundary. Denote by $\rho_n : B_r \to \mathbb{R}_+$ the conformal factor associated to the pullback of the metric of $M$ through $\psi_n$ (meaning, the pullback of the metric of $M$ by $\psi_n$ is $\rho_n g_0$, where $g_0$ is the euclidean metric on $B_r$).

• for each $n = 1, \dots, N$, a smoothly bounded closed set $\omega_n$ in $\mathbb{R}^2$ whose interior contains the origin, $0 \in \stackrel{\circ}{\omega}_n$.

• a bounded integer interval $I \subset \mathbb{N}^*$ such that

$$ \mu_{\min(I)-1}(\Omega) < \mu_{\min(I)}(\Omega) = \mu_{\max(I)}(\Omega) < \mu_{\max(I)+1}(\Omega). $$

We denote by $\mu$ the common value of the $\mu_i(\Omega)$ for $i \in I$. In particular, note that $\mu > 0$.

• a smooth vector field $\phi$ on $M$ that vanishes outside a small neighborhood of $\partial\Omega$ which does not intersect the sets $\psi_n(B_r)$. Write $(\phi^t)_{t \in \mathbb{R}}$ for its associated flow.

We define

$$ \Omega^{\epsilon} = \phi^{\epsilon^2} \left( \Omega \setminus \bigcup_{n=1}^{N} \psi_n(\epsilon\omega_n) \right). $$

We will denote by $\nu$ a normal vector on $\partial\Omega^\epsilon$ with the following convention: it points outward of $\psi_n(\epsilon\omega_n)$ (and so points into $\Omega^\epsilon$), and on the exterior boundary $\partial\Omega^\epsilon \setminus \cup_n \psi_n(\partial(\epsilon\omega_n)) = \partial\phi^{\epsilon^2}(\Omega)$ it points outward from $\Omega^\epsilon$. Then Stokes' theorem gives for arbitrary vector fields $f$ that $\int_{\Omega} \nabla \cdot f = \int_{\partial\phi^{\epsilon^2}(\Omega)} \langle f, \nu \rangle - \sum_{n=1}^N \int_{\partial\psi_n(\epsilon\omega_n)} \langle f, \nu \rangle$.

We consider two exterior problems, for each $\omega_n$.

1) There exists a unique smooth vector field $W_n : \mathbb{R}^2 \setminus \omega_n \to \mathbb{R}^2$ such that

$$ \left\{ \begin{array}{ll} \Delta W_n = 0 & (\text{in } \mathbb{R}^2 \setminus \omega_n) \\ W_n(x) = \mathcal{O}_{|x| \to \infty} \left( \frac{1}{|x|} \right) & \\ \partial_\nu W_n(x) = -\nu(x) & (\text{on } \partial\omega_n) \end{array} \right. $$

where $\nu$ is the exterior normal on $\partial\omega_n$. This $W_n$ is easily shown to exist; it tends to a finite limit at infinity (no logarithmic term) since $\int_{\partial\omega_n} \partial_\nu W_n = -\int_{\partial\omega_n} \nu = 0$. After subtracting a suitable constant vector, one obtains that $W_n$ vanishes at infinity: $W_n(x) = \mathcal{O}_{|x| \to \infty}(1/|x|)$.

For $p, q \in \{1, 2\}$, let

$$ (M_n)_{p,q} = |\omega_n| \delta_{p,q} + \int_{\mathbb{R}^2 \setminus \omega_n} \nabla(W_n)_p \cdot \nabla(W_n)_q $$

where $(W_n)_p$ is the $p$-th component of $W_n$. Note $M_n$ is a $2 \times 2$ symmetric matrix that, as a sum of two positive matrices, is itself positive. Further, $M_n$ determines the behavior of $W_n$ near infinity:

$$ W_n(x) = \frac{M_n x}{2\pi|x|^2} + \mathcal{O}_{|x| \to \infty} \left( \frac{1}{|x|^2} \right) \quad (3) $$

by [14, formula (15)]; we give a proof of this fact in Section 7. (Incidentally, $M_n$ is linked to the so-called *virtual mass*, which appears in [16, Note G].)

2) Similarly, there exists a unique $2 \times 2$ matrix-valued function $L_n : \mathbb{R}^2 \setminus \omega_n \to M_2(\mathbb{R})$ such that

$$ \left\{ \begin{array}{ll} \Delta L_n = 0 & (\text{in } \mathbb{R}^2 \setminus \omega_n) \\ \partial_\nu L_n(x) = x \otimes \nu(x) & (\text{on } \partial\omega_n) \\ L_n(x) = \left(\frac{|\omega_n|}{2\pi} \log|x|\right) I_2 + \mathcal{O}_{|x| \to \infty} \left(\frac{1}{|x|}\right). & \end{array} \right. $$

This time the existence proof uses the divergence theorem to show that

$$ \int_{\partial\omega_n} \partial_\nu L_n(x) = |\omega| I_2 = \int_{\partial B_R} \partial_\nu \left( \frac{|\omega_n|}{2\pi} \log|x| \right) I_2 $$
for each large $R > 0$.

When these exterior problems are defined from a domain $\omega$, we may also denote them by $W_\omega$ (associated to the virtual mass matrix $M_\omega$) and $L_\omega$.

The next proposition quantifies how the eigenvalues of $\Omega^\epsilon$ perturb away from the eigenvalues of $\Omega$. The formula involves the vector field $\phi$ that was chosen above on a neighborhood of $\partial\Omega$.

**Proposition 3.** If the $(\kappa_i)_{i \in I}$ are the (increasingly ordered) eigenvalues of the real symmetric matrix

$$ \left( \int_{\partial\Omega} (\langle \nabla u_i, \nabla u_j \rangle - \mu u_i u_j) \langle \phi, \nu \rangle + \sum_{n=1}^{N} (\rho_n(0)\mu|\omega_n|u_i(O_n)u_j(O_n) - \langle M_n \nabla(u_i \circ \psi_n)(0), \nabla(u_j \circ \psi_n)(0) \rangle) \right)_{i,j \in I \times I} $$

then for some $\epsilon_0 > 0$ we have

$$ |\mu_i(\Omega^\epsilon) - \mu - \epsilon^2 \kappa_i| \le (\text{const.}) \epsilon^{2+1/4} $$

for all $\epsilon \in (0, \epsilon_0)$.

In the case where $\omega$ is simply connected and smooth, the matrix $M_\omega$ (associated to the domain $\omega$) may be computed as follows: write $D$ for the unit disk in $\mathbb{R}^2 \simeq \mathbb{C}$ and denote by $T_\omega : \mathbb{R}^2 \setminus D \to \mathbb{R}^2 \setminus \omega$ the unique biholomorphic conformal mapping that maps $\infty$ to $\infty$ with the normalization

$$ T_{\omega}(z) = \mathrm{cap}[\omega] \left( z + \sum_{k \ge 0} c_k[\omega] z^{-k} \right) $$

where $\mathrm{cap}[\omega]$ is the logarithmic capacity of $\omega$. Let $e_\theta = \begin{pmatrix} \cos \theta \\ \sin \theta \end{pmatrix}$. Then

$$ \langle M_{\omega} e_{\theta}, e_{\theta} \rangle = 2\pi \mathrm{cap}[\omega]^2 (1 - \mathrm{Re}(c_1[\omega]e^{-2i\theta})) \quad (4) $$

by using $\varphi = W \cdot e_\theta$ and $\vec{h} = e_\theta$ in [16, eq. (11) p. 244]. In particular, $\mathrm{Tr}(M_\omega) = 4\pi \mathrm{cap}[\omega]^2$.

**3.2. Definition of the ansatz.** First we choose a $L^2(\Omega)$-orthonormal basis of the eigenspace associated to $\mu$ (denoted $(u_i)_{i \in I}$) such that the symmetric matrix of Proposition 3 is diagonal, with the $(i, i)$ coefficient being $\kappa_i$.

We let $\zeta \in C_c^\infty(B_r, \mathbb{R})$ be a function that is equal to 1 in a neighborhood of 0. We then let $F_{n,i} : \mathbb{R}^2 \setminus \{0\} \to \mathbb{R}$ be the function defined by

$$ F_{n,i}(x) = (\rho_n^{-1}\Delta + \mu) \left[ \zeta(x) \left( \frac{\langle x, M_n \nabla(u_i \circ \psi_n)(0) \rangle}{2\pi|x|^2} + \frac{\rho_n(0)\mu|\omega_n|u_i(O_n)}{2\pi} \log|x| \right) \right]. \quad (5) $$

Since $x/|x|^2$ and $\log|x|$ are harmonic on $\mathbb{R}^2 \setminus \{0\}$, their Laplacians vanish. Thus the singularity of $F_{n,i}(x)$ near the origin arises only from the factor of $\mu$ times $x/|x|^2$ and $\log|x|$, and so $F_{n,i}$ is dominated there by $|x/|x|^2| = 1/|x|$. Hence $F_{n,i} \in L^p(\mathbb{R}^2)$ for all $p \in [1, 2)$.

Let $\Delta_t$ be the laplacian of the metric $g_t$ obtained by pulling back the metric of $M$ through the flow $\phi^t$, and $\nu_t$ be the exterior normal vector of $\partial\Omega$ for this same metric (meaning $\nu_t = (\phi^t)^*\nu_{\partial\phi^t(\Omega)}$ where $\nu_{\partial\phi^t(\Omega)}$ is the exterior normal vector of $\partial\phi^t(\Omega)$), and denote

$$ Kf := \lim_{t \to 0} \frac{\Delta_t f - \Delta f}{t}, \quad \eta = \lim_{t \to 0} \frac{\nu_t - \nu}{t} \quad \text{on } \partial\Omega. $$

Notice $Kf$ vanishes outside the neighborhood of $\partial\Omega$ on which the flow acts nontrivially.

We assert that for every $i$, the following equation defines a unique function $v_i \in H^1(\Omega)$:

$$ \begin{gathered} -\Delta v_i - \mu v_i = \kappa_i u_i + Ku_i + \sum_{n=1}^{N} F_{n,i} \circ \psi_n^{-1} \quad (\text{in } \Omega), \\ \partial_{\nu} v_i = -\langle \eta, \nabla u_i \rangle \quad (\text{on } \partial\Omega). \end{gathered} \quad (6) $$

$$ \int_{\Omega} v_i u_j = 0, \quad \forall j \in I $$
Indeed, to justify the solvability it suffices to observe that the righthand side of the equation in $\Omega$ belongs to $L^p$ for some $p > 1$ (in fact, for all $p \in (1, 2)$) and then to show for each $j \in I$ that

$$ \int_{\Omega} \left( \kappa_i u_i + Ku_i + \sum_{n=1}^{N} F_{n,i} \circ \psi_n^{-1} \right) u_j = - \int_{\partial \Omega} (\partial_{\nu} v_i) u_j. $$

This last assertion is equivalent to

$$ \kappa_i \delta_{i,j} = - \int_{\Omega} \left( Ku_i + \sum_{n=1}^{N} F_{n,i} \circ \psi_n^{-1} \right) u_j + \int_{\partial \Omega} \langle \eta, \nabla u_i \rangle u_j. $$

It may be checked that the right-hand side above is exactly the matrix of Proposition 3, and therefore this equation can be verified by our definitions of the $(\kappa_i)_{i \in I}$: we do not give the full computation, but instead point out that

$$ - \int_{\Omega} (K u_i) u_j + \int_{\partial \Omega} \langle \eta, \nabla u_i \rangle u_j = \int_{\partial \Omega} (\langle \nabla u_i, \nabla u_j \rangle - \mu u_i u_j) \langle \phi, \nu \rangle \quad (7) $$

by shape derivative computations with Neumann boundary conditions (see Section 7 for details that are essentially based on the computations in Fall and Weth [5]) and

$$ - \int_{\Omega} (F_{n,i} \circ \psi_n^{-1}) u_j = \rho_n(0) \mu |\omega_n| u_i(O_n) u_j(O_n) - \langle M_n \nabla (u_i \circ \psi_n)(0), \nabla (u_j \circ \psi_n)(0) \rangle \quad (8) $$

by an application of Green's theorem (for details, see Section 7).

We may then define the ansatz: for $i \in I$ and $\epsilon > 0$, let

$$ \begin{aligned} \hat{u}_{i,\epsilon}(x) = u_i(x) + \epsilon^2 v_i(x) + \epsilon \sum_{n=1}^{N} \zeta(\psi_n^{-1}(x)) \langle W_n(\psi_n^{-1}(x)/\epsilon), \nabla (u_i \circ \psi_n)(0) \rangle \\ + \epsilon^2 \sum_{n=1}^{N} \zeta(\psi_n^{-1}(x)) \left( -L_n(\psi_n^{-1}(x)/\epsilon) : D^2(u_i \circ \psi_n)(0) + \frac{\rho_n(0)\mu|\omega_n|u_i(O_n)}{2\pi} \log \epsilon \right) \end{aligned} \quad (9) $$

for $x \in \Omega \setminus \cup_{n=1}^N \psi_n(\epsilon\omega_n)$, where $D^2$ is the Hessian operator and $A: B \stackrel{\text{def}}{=} \sum_{k,l=1}^2 A_{k,l}B_{k,l}$,

$u_{i,\epsilon} = \hat{u}_{i,\epsilon} \circ \phi^{-\epsilon^2}$ (defined on $\Omega^\epsilon$),

$\mu_{i,\epsilon} = \mu + \epsilon^2 \kappa_i$.

(In the definition of $\hat{u}_{i,\epsilon}(x)$, one takes $\zeta(\psi_n^{-1}(x)) = 0$ when $x \notin \psi_n(B_r)$.)

**Lemma 4.** There exists some constant $C > 0$ such that for every small enough $\epsilon > 0$ and every $z \in H^1(\Omega^\epsilon)$ with unit norm, we have

$$ \int_{\Omega^\epsilon} (\mu_{i,\epsilon} u_{i,\epsilon} z - \langle \nabla u_{i,\epsilon}, \nabla z \rangle) \le C \epsilon^{2+\frac{1}{2}}. $$

We defer the proof of this lemma to Section 7.

3.3. **Proof of Proposition 3.** For the proof it is more convenient to work with a compact operator on $\Omega^\epsilon$ rather than an unbounded one. Denote by

$$ \langle u, v \rangle_{H^1(\Omega^\epsilon)} = \int_{\Omega^\epsilon} (\langle \nabla u, \nabla v \rangle + uv) $$

the inner product that defines the norm on $H^1(\Omega^\epsilon)$.

We first define, for a domain $V \subset M$, the operator $S_V: L^2(M) \to L^2(M)$ such that for any $f \in L^2(M)$ we have $S_V f = 0$ in $M \setminus V$ and $S_V f|_V \in H^1(V)$, and

$$ \begin{cases} S_V f - \Delta(S_V f) = f & (\text{in } V) \\ \partial_{\nu_V} S_V f = 0 & (\text{on } \partial V) \end{cases} $$
and such that the associated eigenvalues of $S_V$ are exactly the numbers $\left(\frac{1}{1+\mu_i(V)}\right)_{i \in \mathbb{N}}$ (and 0).

These eigenvalues are all bounded by 1 and so the operator norm is likewise bounded: $\|S_V\| \le 1$.

We claim $S_{\Omega^\epsilon}$ converges to $S_\Omega$ in the $L^2(M)$-operator norm as $\epsilon \to 0$. Indeed $\Omega^\epsilon$ (and its complement) converges in the Hausdorff sense as $\epsilon \to 0$ to $\Omega \setminus \{O_1, \dots, O_N\}$ (and its complement), which is equivalent to $\Omega$ with regard to the Laplace problem. Moreover, the set

$$ \bigcup_{0<\epsilon<1} \{u \in L^2(M) : \mathrm{Spt}(u) \subset \Omega^\epsilon, \|u\|_{H^1(\Omega^\epsilon)} \le 1\} \qquad (10) $$

is compact in $L^2(M)$, since there is a uniform $H^1$-extension operator to a neighborhood of $\bigcup_{0<\epsilon<1} \Omega^\epsilon$ (achieved by extending harmonically into the holes using Lemma 7 in Section 7, and taking a standard $H^1$-extension across the outer boundary of $\Omega^\epsilon$ to obtain an extended function with compact support in the neighborhood). To prove the operator norm convergence, consider a bounded sequence $(f_\epsilon)_\epsilon$ in $L^2(M)$. It is sufficient to prove that $\|S_\Omega f_\epsilon - S_{\Omega^\epsilon} f_\epsilon\|_{L^2(M)} \to 0$. In particular, it suffices to do so along every subsequence. From any subsequence we may extract another subsequence (denoted still $f_\epsilon$, for simplicity) such that $f_\epsilon$ converges $L^2(M)$-weakly to some $f \in L^2(M)$. Since $S_\Omega$ is compact (by Rellich's theorem), $S_\Omega f_\epsilon$ converges in $L^2(M)$ to its limit $v = S_\Omega f$ that is characterized by the weak form of the equation:

$$ \langle \nabla v, \nabla \varphi \rangle_{L^2(M)} + \langle v, \varphi \rangle_{L^2(M)} = \langle f \chi_{\Omega^\epsilon}, \varphi \rangle_{L^2(M)}, \quad \varphi \in C^\infty(M), \qquad (11) $$

where $\nabla v = 0$ outside $\Omega$ by convention. Now let $v_\epsilon = S_{\Omega^\epsilon} f_\epsilon$, so that $\|v_\epsilon\|_{H^1(\Omega^\epsilon)}$ is bounded independently of $\epsilon$. The compactness in (10) above yields a subsequence $v_{\epsilon_k}$ that converges in $L^2(M)$ to some $w \in L^2(M)$. Similarly, $\nabla v_{\epsilon_k}$ (which we recall is extended by 0 outside $\Omega^{\epsilon_k}$) is bounded in $L^2(M)$, and so some subsequence converges $L^2(M)$-weakly to a limiting vector field $G \in L^2(M)$. By the Hausdorff convergence of the domains, $w$ and $G$ are supported in $\Omega$. Using test functions supported in $\Omega$, one finds $w$ is weakly differentiable there with $\nabla w = G$. Next, the partial differential equation satisfied by $v_{\epsilon_k}$ says in its weak form that

$$ \langle \nabla v_{\epsilon_k}, \nabla \varphi \rangle_{L^2(M)} + \langle v_{\epsilon_k}, \varphi \rangle_{L^2(M)} = \langle f \chi_{\Omega^{\epsilon_k}}, \varphi \rangle_{L^2(M)}, \quad \varphi \in C^\infty(M). $$

Passing to the limit as $\epsilon_k \to 0$ implies that

$$ \langle \nabla w, \nabla \varphi \rangle_{L^2(M)} + \langle w, \varphi \rangle_{L^2(M)} = \langle f \chi_{\Omega^\epsilon}, \varphi \rangle_{L^2(M)}, \quad \varphi \in C^\infty(M). \qquad (12) $$

Thus the functions $v$ in (11) and $w$ in (12) satisfy the same variational characterization, so that $w = v$ by uniqueness. Hence $S_{\Omega^{\epsilon_k}} f_{\epsilon_k} = v_{\epsilon_k} \to w = v = S_\Omega f$. That concludes the proof that $S_{\Omega^\epsilon} \to S_\Omega$ in the operator norm.

As a consequence, for every $i$ we have

$$ \mu_i(\Omega^\epsilon) \xrightarrow{\epsilon \to 0} \mu_i(\Omega). \qquad (13) $$

Now define a compact resolvent operator $R_\epsilon : H^1(\Omega^\epsilon) \to H^1(\Omega^\epsilon)$ by

$$ \forall u, v \in H^1(\Omega^\epsilon) : \quad \langle R_\epsilon u, v \rangle_{H^1(\Omega^\epsilon)} = \langle u, v \rangle_{L^2(\Omega^\epsilon)}. \qquad (14) $$

In other words, $R_\epsilon u$ verifies $R_\epsilon u - \Delta R_\epsilon u = u$ in $\Omega^\epsilon$ and $\partial_\nu R_\epsilon u = 0$ on $\partial\Omega^\epsilon$.

Let $(\mu_i(\Omega^\epsilon), u_i(\Omega^\epsilon))_{i \in \mathbb{N}}$ be the eigenvalues and eigenfunctions of $-\Delta$ on $\Omega^\epsilon$ with Neumann conditions, normalized in $L^2(\Omega^\epsilon)$. Then the normalized eigenpairs associated to $R_\epsilon$ are exactly

$$ \left( \frac{1}{1 + \mu_i(\Omega^\epsilon)}, \frac{u_i(\Omega^\epsilon)}{\|u_i(\Omega^\epsilon)\|_{H^1(\Omega^\epsilon)}} \right). $$

where the eigenvalues are this time enumerated in a decreasing way. For an interval $A \subset \mathbb{R}$ and $u \in H^1(\Omega^\epsilon)$, we define the projection operator

$$ P_{\epsilon}^{A} u = \sum_{i:(1+\mu_i(\Omega^\epsilon))^{-1} \in A} \left\langle u, \frac{u_i(\Omega^\epsilon)}{\|u_i(\Omega^\epsilon)\|_{H^1(\Omega^\epsilon)}} \right\rangle_{H^1(\Omega^\epsilon)} \frac{u_i(\Omega^\epsilon)}{\|u_i(\Omega^\epsilon)\|_{H^1(\Omega^\epsilon)}}. $$
Let now $\lambda_{i,\epsilon} = \frac{1}{1+\mu_{i,\epsilon}}$, $U_{i,\epsilon} = \frac{u_{i,\epsilon}}{\|u_{i,\epsilon}\|_{H^1(\Omega^\epsilon)}}$. Note that $\|u_{i,\epsilon}\|_{H^1(\Omega^\epsilon)} = \sqrt{\mu_i(\Omega^\epsilon) + 1} \xrightarrow{\epsilon \to 0} \sqrt{\mu + 1}$. Then

$$
\begin{align*}
\|R_{\epsilon}U_{i,\epsilon} - \lambda_{i,\epsilon}U_{i,\epsilon}\|_{H^1(\Omega^{\epsilon})} &= \frac{1}{(1+\mu_{i,\epsilon})\|u_{i,\epsilon}\|_{H^1(\Omega^{\epsilon})}} \|(1+\mu_{i,\epsilon})R_{\epsilon}u_{i,\epsilon} - u_{i,\epsilon}\|_{H^1(\Omega^{\epsilon})} \\
&= \mathcal{O}(1) \sup_{\|z\|_{H^1(\Omega^{\epsilon})} \le 1} \langle (1+\mu_{i,\epsilon})R_{\epsilon}u_{i,\epsilon}, z \rangle_{H^1(\Omega^{\epsilon})} - \langle u_{i,\epsilon}, z \rangle_{H^1(\Omega^{\epsilon})} \\
&= \mathcal{O}(1) \sup_{\|z\|_{H^1(\Omega^{\epsilon})} \le 1} \langle (1+\mu_{i,\epsilon})u_{i,\epsilon}, z \rangle_{L^2(\Omega^{\epsilon})} - \langle u_{i,\epsilon}, z \rangle_{H^1(\Omega^{\epsilon})} \quad \text{by (14)} \\
&= \mathcal{O}(1) \sup_{\|z\|_{H^1(\Omega^{\epsilon})} \le 1} \int_{\Omega^{\epsilon}} (\mu_{i,\epsilon}u_{i,\epsilon}z - \langle \nabla u_{i,\epsilon}, \nabla z \rangle) \\
&\le C\epsilon^{2+\frac{1}{2}} \tag{15}
\end{align*}
$$

by Lemma 4. Consequently, defining the interval $A_{i,\epsilon} = [\lambda_{i,\epsilon} - \epsilon^{2+1/4}, \lambda_{i,\epsilon} + \epsilon^{2+1/4}]$, we have

$$
\begin{align*}
& \| U_{i,\epsilon} - P_{\epsilon}^{A_{i,\epsilon}} U_{i,\epsilon} \|_{H^1(\Omega^{\epsilon})}^2 \\
&= \sum_{j:(1+\mu_j(\Omega^{\epsilon}))^{-1} \notin A_{i,\epsilon}} \left\langle U_{i,\epsilon}, \frac{u_j(\Omega^{\epsilon})}{\|u_j(\Omega^{\epsilon})\|_{H^1(\Omega^{\epsilon})}} \right\rangle_{H^1(\Omega^{\epsilon})}^2 \\
&\le \epsilon^{-4-1/2} \sum_{j:(1+\mu_j(\Omega^{\epsilon}))^{-1} \notin A_{i,\epsilon}} \left( (1+\mu_j(\Omega^{\epsilon}))^{-1} - \lambda_{i,\epsilon} \right)^2 \left\langle U_{i,\epsilon}, \frac{u_j(\Omega^{\epsilon})}{\|u_j(\Omega^{\epsilon})\|_{H^1(\Omega^{\epsilon})}} \right\rangle_{H^1(\Omega^{\epsilon})}^2 \\
&\le \epsilon^{-4-1/2} \|(R_{\epsilon} - \lambda_{i,\epsilon})U_{i,\epsilon}\|_{H^1(\Omega^{\epsilon})}^2 \\
&\le C\epsilon^{1/2} \tag{16}
\end{align*}
$$

by squaring (15).

When $\epsilon < C^{-2}$, this inequality already gives that $P_{\epsilon}^{A_{i,\epsilon}} U_{i,\epsilon} \neq 0$, and so $A_{i,\epsilon}$ contains some eigenvalue of $\Omega^{\epsilon}$. To finish proving the proposition, it is necessary to be a bit more attentive about the multiplicity of the eigenvalues $(\kappa_i)_{i \in I}$. By the convergence of the eigenvalues in (13), there exists $\delta > 0$ such that for all sufficiently small $\epsilon > 0$,

$$ \mu_i(\Omega^\epsilon) \in [\mu - \delta, \mu + \delta] \text{ if and only if } i \in I. $$

This insures that the eigenvalues we obtain in the intervals $A_{i,\epsilon}$ have the correct index (belonging to $I$). We partition $I$ into $\cup_{r=1}^s I_r$ such that $i, j$ are in the same subinterval of the partition if and only if $\kappa_i = \kappa_j$. Let $B_{r,\epsilon} = A_{i,\epsilon}$ for any $i \in I_r$. Then for small enough $\epsilon$, the intervals $(B_{r,\epsilon})_{r=1,...,s}$ are disjoint. We also argue that for $\epsilon$ small enough, the functions

$$ (P_{\epsilon}^{B_{r,\epsilon}} U_{i,\epsilon})_{i \in I_r}, \quad r = 1, \dots, s, $$

are linearly independent. Independence is clear between different $r$-values due to disjointness of the $B_{r,\epsilon}$. For functions with the same $r$-value, the matrix $\langle P_{\epsilon}^{B_{r,\epsilon}} U_{i,\epsilon}, P_{\epsilon}^{B_{r,\epsilon}} U_{j,\epsilon} \rangle_{H^1(\Omega^{\epsilon})}$ of inner products converges by (16) to the same limit (as $\epsilon \to 0$) as the matrix $\langle U_{i,\epsilon}, U_{j,\epsilon} \rangle_{H^1(\Omega^{\epsilon})}$, which converges to the identity.

Since the total number of independent eigenfunctions on $\Omega^{\epsilon}$ associated to eigenvalues in the interval $[\mu - \delta, \mu + \delta]$ is exactly $|I|$, and since

$$ (1 + \mu_i(\Omega^\epsilon))^{-1} \in A_{i,\epsilon} \implies |\mu_i(\Omega^\epsilon) - (\mu + \epsilon^2 \kappa_i)| \le C\epsilon^{2+1/4}, $$

we see the proof of Proposition 3 is complete.
4. PROOF OF THE MAIN RESULT, THEOREM 1

Let us first treat the case $\theta \in (\Theta, \pi)$. Take $t = t_{max}(\theta) \in (\pi/2, \Theta)$ as in Proposition 2(b), noting that $g'(t) = 0$ and hence $\mu_1(B_t) = \mu_1(B_\theta)$. We consider the four equally spaced points $x_0, x_1, x_2, x_3$ at latitude $t$ on the unit sphere with spherical coordinates

$$ \left( t, \frac{j\pi}{2} \right), \quad j = 0, 1, 2, 3. $$

Let $\lambda = \lambda(\theta) > 1$ be a parameter that will later be fixed close to 1, define

$$ T(z) = \frac{\lambda z + \lambda^{-1} z^{-1}}{\sqrt{\lambda^2 - \lambda^{-2}}}, $$

and consider the domain $\omega = \mathbb{R}^2 \setminus T(\mathbb{R}^2 \setminus D)$, where $D$ is the unit disk in $\mathbb{R}^2 \simeq \mathbb{C}$. In other words $\omega$ is an ellipse of area $\pi$ whose horizontal and vertical semi-axes are respectively

$$ \sqrt{\frac{\lambda + \lambda^{-1}}{\lambda - \lambda^{-1}}} \quad \text{and} \quad \sqrt{\frac{\lambda - \lambda^{-1}}{\lambda + \lambda^{-1}}}. $$

Now take a small number $\epsilon > 0$, let $\psi_j : \mathbb{R}^2 \to \mathbb{S}^2$ be the stereographic projection such that $\psi_j(0) = x_j$ and $d\psi(0)e_2$ is a unit vector pointing upward (in the direction of the line of longitude) in the tangent plane to the sphere, and let $\Omega = B_\theta$ and

$$ \Omega^\epsilon := B_{\theta(\epsilon)} \setminus \bigcup_{j=0}^{3} \psi_j(\epsilon\omega), $$

where $\theta(\epsilon)$ is chosen such that $|\Omega^\epsilon| = |B_{\theta(\epsilon)}|$: this means that $\theta(\epsilon) = \theta + \frac{4\pi\epsilon^2}{|\partial B_\theta|} + o(\epsilon^2)$. Thus $\Omega^\epsilon$ is a spherical cap with four (approximately) elliptical holes removed near aperture $t$, with the longer axis of each hole being (approximately) horizontal.

Let us now compute, for a fixed $\lambda$, the asymptotic of $\mu_1(\Omega^\epsilon)$ as $\epsilon \to 0$. We denote by

$$ u_1(\cdot, \phi) = g(\cdot) \cos \phi \quad \text{and} \quad u_2(\cdot, \phi) = g(\cdot) \sin \phi \qquad (17) $$

the eigenfunctions corresponding to the eigenvalue $\mu = \mu_1(B_\theta) = \mu_2(B_\theta)$, where $g$ is normalized such that $u_1, u_2$ are orthonormal in $L^2(B_\theta)$. Since $\Omega$ is invariant by rotation of the sphere by 90 degrees around the vertical axis and the eigenfunctions associated to $\mu_1(\Omega^\epsilon)$ and $\mu_2(\Omega^\epsilon)$ are close (in the $L^2$ sense) to linear combinations of $u_1, u_2$ that are not invariant under the rotation, it follows that $\mu_1(\Omega^\epsilon) = \mu_2(\Omega^\epsilon)$.

Proposition 3 with $N = 4$ and $I = \{1, 2\}$ gives the convergence of $\mu_1(\Omega^\epsilon)$ to $\mu_1(B_\theta)$, with

$$ \begin{aligned} \mu_1(\Omega^\epsilon) - \mu_1(B_\theta) &= \frac{1}{2} \sum_{i=1}^{2} [\mu_i(\Omega^\epsilon) - \mu] = \frac{\epsilon^2}{2}(\kappa_1 + \kappa_2) + o(\epsilon^2) \\ &= \frac{\epsilon^2}{2}(\text{trace of the matrix in Proposition 3}) + o(\epsilon^2). \end{aligned} $$

Evaluating the diagonal entries of the matrix in Proposition 3 gives that

$$ \begin{aligned} \mu_1(\Omega^\epsilon) - \mu_1(B_\theta) &= \frac{\epsilon^2}{2} \sum_{i=1}^2 \sum_{j=0}^3 [\mu_i(B_\theta) |\omega| |u_i(x_j)|^2 - \langle M_\omega \psi_j^* \nabla u_i(x_j), \psi_j^* \nabla u_i(x_j) \rangle] \\ &\quad + \frac{1}{2} \frac{4\pi\epsilon^2}{|\partial B_\theta|} \sum_{i=1}^2 \int_{\partial B_\theta} (|\nabla u_i|^2 - \mu_i(B_\theta) |u_i|^2) + o(\epsilon^2) \\ &= 2\pi\epsilon^2 \left[ \mu_1(B_\theta) g(t)^2 - \frac{(M_\omega)_{1,1}}{\pi} \frac{g(t)^2}{\sin^2 t} - \frac{(M_\omega)_{2,2}}{\pi} g'(t)^2 + \left( \frac{1}{\sin^2 \theta} - \mu_1(B_\theta) \right) g(\theta)^2 \right] + o(\epsilon^2), \end{aligned} \qquad (18) $$
as we now explain in detail. Based on the form of $u_1$ and $u_2$ in equation (17), one has $u_1(\theta, \phi)^2+u_2(\theta, \phi)^2 = g(\theta)^2$ and

$$|\nabla u_1(\theta, \phi)|^2 + |\nabla u_2(\theta, \phi)|^2 = g'(\theta)^2 + \frac{g(\theta)^2}{\sin(\theta)^2} = \frac{g(\theta)^2}{\sin(\theta)^2}$$

since $g'(\theta) = 0$. For fixed $j$, the values taken by $u_1 \circ \psi_j$ and $u_2 \circ \psi_j$ at the origin are $\{\pm g(t), 0\}$ and the gradients are $\{\pm g'(t)e_2, \pm \frac{g(t)}{\sin(t)}e_1\}$. Here each $\pm$ sign depends on $j$, but fortunately, the choice of sign makes no difference to the following calculation. Hence

$$\begin{aligned} & \sum_{i=1}^{2} [\mu_i(B_\theta)|\omega| |u_i(x_j)|^2 - \langle M_\omega \psi_j^* \nabla u_i(x_j), \psi_j^* \nabla u_i(x_j) \rangle] \\ & = \mu_1(B_\theta)|\omega|g(t)^2 - (M_\omega)_{1,1} \frac{g(t)^2}{\sin^2 t} - (M_\omega)_{2,2} g'(t)^2. \end{aligned}$$

Summing over $j = 0, 1, 2, 3$ and recalling that $|\omega| = \pi$ yields the final line of (18).

The third and fourth terms in (18) vanish because our choice $t = t_{max}(\theta)$ gives $g'(t) = 0$. For the fifth term, since $\theta > \Theta$ we have $\frac{1}{\sin^2 \theta} - \mu_1(B_\theta) > 0$. For the second term, we now compute $(M_\omega)_{1,1}$ using formula (4). In this case

$$\begin{aligned} \operatorname{cap}[\omega] &= \frac{\lambda}{\sqrt{\lambda^2 - \lambda^{-2}}}, & c_1[\omega] &= \lambda^{-2}, \\ (M_\omega)_{1,1} &= \frac{2\pi\lambda^2}{\lambda^2 - \lambda^{-2}} \left(1 - \frac{1}{\lambda^2}\right) = \frac{2\pi}{1 + \lambda^{-2}}. \end{aligned}$$

Coming back to (18), we obtain

$$\begin{aligned} \mu_1(\Omega^\epsilon) - \mu_1(B_\theta) &\ge 2\pi\epsilon^2 \left[ \mu_1(B_\theta)g(t)^2 - \frac{1}{\pi}\frac{2\pi}{1+\lambda^{-2}}\frac{g(t)^2}{\sin^2 t} \right] + o(\epsilon^2) \\ &= 2\pi\epsilon^2 \frac{g(t)^2}{\sin^2 t} \left[ \mu_1(B_\theta)\sin^2 t - \frac{2}{1+\lambda^{-2}} \right] + o(\epsilon^2). \end{aligned}$$

Our choice of $t$ ensures that $\mu_1(B_\theta) = \mu_1(B_t)$, and since $t < \Theta$ we know that $\mu_1(B_t) \sin^2 t > 1$. Thus we may choose $\lambda$ close enough to 1 that $\mu_1(B_\theta) \sin^2 t > \frac{2}{1+\lambda^{-2}}$. Notice this choice of $\lambda$ depends only on $\theta$. Hence for all small enough $\epsilon$, we have $\mu_1(\Omega^\epsilon) > \mu_1(B_\theta)$. This concludes the proof of the theorem in the case $\theta > \Theta$.

**Remark 5.** If we sought instead counterexamples to the maximality of the geodesic ball for the arithmetic mean $(\mu_1 + \mu_2)/2$ or harmonic mean $2/(\frac{1}{\mu_1} + \frac{1}{\mu_2})$, then it would suffice to employ a single elliptical hole instead of four holes as used above. Piercing by multiple holes is done simply to ensure the double multiplicity of $\mu_1$. Thus while the harmonic mean is maximal at a geodesic ball among simply connected domains of fixed area smaller than $0.94|\mathbb{S}^2|$, by [10], for each area larger than $|B_\Theta|(\approx 0.8|\mathbb{S}^2|)$ a doubly connected domain exists whose harmonic mean of $\mu_1$ and $\mu_2$ is larger than for the geodesic ball of the same area. $\square$

The remaining case in the theorem, for $\theta \in (\Theta - \delta, \Theta]$, follows from the next proposition. Recall that $\omega = \omega_\lambda$ is an ellipse of area $\pi$. As above, $\Omega^\epsilon = B_{\theta(\epsilon)} \cup \bigcup_{j=0}^3 \psi_j(\epsilon\omega)$ is a cap with four holes such that $|\Omega^\epsilon| = |B_\theta|$.

**Proposition 6.** There exists $\delta > 0, \lambda > 1$ and $t < \Theta - \delta$ such that if $\theta \in (\Theta - \delta, \Theta]$ then $\mu_1(\Omega^\epsilon) > \mu_1(B_{\theta(\epsilon)})$ for all small $\epsilon > 0$.

*Proof.* Let $t \in (0, \Theta)$ and compute the asymptotic of $\mu_1(\Omega^\epsilon) - \mu_1(B_\theta)$ as follows. With the same computation as previously, applying Proposition 3 to $\Omega^\epsilon$ gives, following equation (18):

$$\mu_1(\Omega^\epsilon) - \mu_1(B_\theta) = 2\pi\epsilon^2 \left( \mu_1(B_\theta)g(t)^2 - \frac{(M_{\omega_\lambda})_{1,1}}{\pi} \frac{g(t)^2}{\sin^2 t} - \frac{(M_{\omega_\lambda})_{2,2}}{\pi} g'(t)^2 + (1 - \mu_1(B_\theta)\sin^2 \theta) \frac{g(\theta)^2}{\sin^2 \theta} \right) + o(\epsilon^2).$$
The factor in parentheses multiplying $2\pi\epsilon^2$ is continuous with respect to $\theta, t, \lambda$, using here that the radial part $g$ of the eigenfunction on $B_\theta$ varies continuously with the radius $\theta$ of that ball, and so to finish the proof of the theorem it is enough to prove that the factor is positive at $\theta = \Theta$ for some suitable choice of $t$ and $\lambda$; then $\delta$ can be chosen small enough to ensure that the factor in parentheses remains positive when $\theta \in (\Theta - \delta, \Theta]$. When $\theta = \Theta$ we have $\mu_1(B_\Theta) \sin^2(\Theta) = 1$ and so the task simplifies to finding $t$ and $\lambda$ such that

$$ \mu_1(B_\Theta)g(t)^2 - \frac{(M_{\omega_\lambda})_{1,1}}{\pi} \frac{g(t)^2}{\sin^2 t} - \frac{(M_{\omega_\lambda})_{2,2}}{\pi} g'(t)^2 > 0, \quad (19) $$

where $g$ is the radial part of the eigenfunction for $B_\Theta$, with $\theta = \Theta$.

Let $\lambda = 1 + \eta$ for some small $\eta > 0$ to be fixed later. Using formula (4) as previously, we have

$$ \frac{(M_{\omega_{\lambda}})_{1,1}}{\pi} = \frac{2}{1 + \lambda^{-2}} = 1 + \eta + \mathcal{O}(\eta^2) $$

and similarly

$$ \frac{(M_{\omega_{\lambda}})_{2,2}}{\pi} = \frac{2}{1 - \lambda^{-2}} = \frac{1}{\eta} + \mathcal{O}(1). $$

Note also that $g'(\Theta) = 0$ by the Neumann boundary condition (recall we have chosen $\theta = \Theta$), while the equation (2) satisfied by $g$ implies $g''(\Theta) = 0$ (a special property at the radius $\Theta$, where $\mu_1(B_\Theta) \sin^2(\Theta) = 1$) and $g'''(\Theta) = -2(\cos\Theta/\sin^3\Theta)g(\Theta) > 0$. Note $\cos\Theta < 0$ since $\Theta$ is larger than $\pi/2$. Consequently, writing $t = \Theta - \tau$ where $\tau > 0$ is to be determined, we find

$$ g(\Theta - \tau) = g(\Theta) + \frac{1}{3} \frac{\cos \Theta}{\sin^3 \Theta} g(\Theta) \tau^3 + \mathcal{O}(\tau^4), $$

$$ g'(\Theta - \tau) = -\frac{\cos \Theta}{\sin^3 \Theta} g(\Theta) \tau^2 + \mathcal{O}(\tau^3), $$

$$ \frac{g(\Theta - \tau)^2}{\sin^2(\Theta - \tau)} = \frac{g(\Theta)^2}{\sin^2 \Theta} (1 + 2\tau \cot \Theta + \mathcal{O}(\tau^2)). $$

Hence the left side of (19) equals

$$ \begin{aligned} & \mu_1(B_\Theta)g(\Theta - \tau)^2 - \frac{(M_{\omega_{1+\eta}})_{1,1}}{\pi} \frac{g(\Theta - \tau)^2}{\sin^2(\Theta - \tau)} - \frac{(M_{\omega_{1+\eta}})_{2,2}}{\pi} g'(\Theta - \tau)^2 \\ &= \mu_1(B_\Theta)(g(\Theta)^2 + \mathcal{O}(\tau^3)) - (1 + \eta + \mathcal{O}(\eta^2)) \frac{g(\Theta)^2}{\sin^2 \Theta} (1 + 2\tau \cot \Theta + \mathcal{O}(\tau^2)) \\ & \quad - (\eta^{-1} + \mathcal{O}(1)) \left( -\frac{\cos \Theta}{\sin^3 \Theta} g(\Theta) \tau^2 + \mathcal{O}(\tau^3) \right)^2 \\ &= (-2\tau \cot \Theta - \eta) \frac{g(\Theta)^2}{\sin^2 \Theta} + \mathcal{O} \left( \tau^2 + \eta^2 + \frac{\tau^4}{\eta} \right) \end{aligned} $$

by substituting $\mu_1(B_\Theta) = 1/\sin^2\Theta$ and simplifying. We now specify $\tau$ in terms of $\eta$ by requiring $\eta = -\tau \cot\Theta > 0$. Then (19) holds for any small enough choice of $\eta$, because its left side equals

$$ \eta \frac{g(\Theta)^2}{\sin^2 \Theta} + \mathcal{O}(\eta^2) > 0. $$

This inequality finishes the proof of Proposition 6.

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

## 5. NUMERICAL AND HEURISTIC ANALYSIS OF A NEW CLASS OF DOMAINS

In order to understand how low the value of $\Theta - \delta$ in Theorem 1 could be, in this section we are interested in a different kind of geometric perturbation of a spherical cap. This kind of counterexample is loosely inspired by the numerical optima found by Martinet [12]. The arguments we bring in this section are not rigorous enough to produce a fully mathematically justified lower value for $\Theta - \delta$. However, our analysis is supported by numerical computations, gives a good hint of a new class of competitive geometries, and could possibly be made rigorous with extra effort. The sets we deal with consist of the union of a spherical
cap (geodesic ball) and a graph-like domain. For a geodesic cap $B_{\theta}$ having its center at $e_3$, define the domain

$$ \Omega_{\varepsilon}(\theta) = B_{\theta\varepsilon} \cup \{|x_1| < \varepsilon\} \cup \{|x_2| < \varepsilon\} $$

to be the union of a cap and two strips, where the aperture $\theta^{\varepsilon} < \theta$ is chosen such that $|\Omega_{\varepsilon}(\theta)| = |B_{\theta}|$. Such domains, for several values of $\theta$ and $\varepsilon$, are depicted Figure 2. We call them “helmet” domains, although if they were turned upside down they could also be called “basket” domains.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Sets $\Omega_{\varepsilon}(\theta)$ for several values of the aperture $\theta$ and for strip thickness $\varepsilon = 0.05$ (top) and $\varepsilon = 0.1$ (bottom). Here $e_3$ is oriented vertically.</figcaption></figure>

Using the finite element software FreeFem++ [7], we can generate meshes on such domains and compute their eigenvalues. Figure 3 plots the value of $\mu_1(\Omega_{\varepsilon}(\theta))$ as a function of $\theta$, for different values of $\varepsilon > 0$, each of which provides counterexamples to maximality of the cap ($\varepsilon = 0$) over slightly different intervals of $\theta$.

This figure deserves several comments. Most importantly, the helmet counterexamples work for certain aperture values $\theta$ that are substantially smaller than $\Theta$, which Theorem 1 cannot do due to the very nature of the earlier construction. On the other hand, for a fixed $\varepsilon$ the helmet domains fail to provide counterexamples when $\theta$ is close to $\pi$, whereas Theorem 1 holds all the way up to $\pi$: in other words, as $\varepsilon$ gets small we obtain two new critical angles $\Theta'$, $\Theta''$ with

$$ \frac{\pi}{2} < \Theta' < \Theta < \Theta'' \le \pi $$

such that a helmet domain gives a counterexample for $\theta \in (\Theta', \Theta'')$.

Next, the two eigenvalue curves in the figure exhibit downward-pointing corners around the aperture $\Theta'(< \Theta)$, leading the current constructions to stop being counterexamples as the aperture gets smaller. It is shown in Figure 4 that this corner is linked to the change of nature of the eigenfunction associated to $\mu_1$ (with the eigenfunction transitioning to one that is concentrated on the strips) and is a crossing point for which we have $\mu_1 = \mu_2 = \mu_3$. Lastly, it seems that decreasing the width $\varepsilon$ of the strip allows one to get counterexamples for smaller values of $\theta$. A question of interest would be to provide the precise asymptotic behaviour of $\mu_1(\Omega_{\varepsilon}(\theta))$ when $\varepsilon$ goes to 0, with the purpose to determine the smallest $\theta$ that can yield a counterexample of this type.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Value of $\mu_1(\Omega_\varepsilon(\theta))$ as a function of $\theta/\pi$, for $\varepsilon = 0.05$ and $\varepsilon = 0.1$. The red curve $\varepsilon = 0$ corresponds to the spherical cap $B_\theta$. These graphs provide numerical counterexamples to maximality of $\mu_1$ at the cap, for cap apertures $\theta \in (0.623\pi, 0.925\pi)$ for $\varepsilon = 0.05$ and $\theta \in (0.646\pi, 0.912\pi)$ for $\varepsilon = 0.1$. The vertical dotted line is at $\Theta \approx 0.7\pi$; recall that Theorem 1 provides rigorous counterexamples to the right of that aperture.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Eigenfunctions associated to $\mu_1$ for $\varepsilon = 0.05$, for $\theta = 1.96$ (left) and $\theta = 1.97$ (right). The red color corresponds to the value 1 while deep blue corresponds to -1.</figcaption></figure>
For replicability purposes, the examples used to generate Figure 3 are available at https://github.com/EloiMartinet/PuncturedSpheres.

Let us now give a heuristic (and non-rigorous) justification of these counterexamples, based on the results of Arrieta [1] proved in the case of simple eigenvalues in “dumbbell” domains. To reduce the number of leaps of faith, let us proceed slightly differently to above by no longer attempting to preserve the original area: simply add the straps to the cap $B_\theta$ by letting

$$H_{\theta,\epsilon} = B_{\theta} \cup \{x : |x_1| \le \epsilon \text{ or } |x_2| \le \epsilon\}$$

where $H$ stands for “Helmet”, and compare with the cap

$$B_{\theta(\epsilon)} \quad \text{where } \theta(\epsilon) \text{ is chosen to satisfy } |B_{\theta(\epsilon)}| = |H_{\theta,\epsilon}|.$$

Let us first compute an asymptotic for $\mu_1(B_{\theta(\epsilon)})$, with a shape derivative. Since $|H_{\theta,\epsilon}| = |B_\theta| + 8(\pi - \theta)\epsilon + o(\epsilon)$, we have $\theta(\epsilon) = \theta + \frac{8(\pi-\theta)\epsilon}{|\partial B_\theta|} + o(\epsilon)$. As a consequence,

$$\begin{aligned} \mu_1(B_{\theta(\epsilon)}) &= \mu_1(B_\theta) + \frac{8(\pi - \theta)\epsilon}{|\partial B_\theta|} \int_{\partial B_\theta} (|\nabla u_1|^2 - \mu_1(B_\theta)u_1^2) \\ &= \mu_1(B_\theta) + 4(\pi - \theta)\epsilon \left( \frac{1}{\sin^2 \theta} - \mu_1(B_\theta) \right) g(\theta)^2 \end{aligned}$$

where $u_1$ is an eigenfunction associated to $\mu_1(B_\theta)$, of the form $g(t) \sin \phi$ or $g(t) \cos \phi$.

Let us now compute an asymptotic for $\mu_1(H_{\theta,\epsilon})$: following Arrieta [1, Theorem 2.5] we expect the Neumann spectrum of $H_{\theta,\epsilon}$ to be given by the union (counting multiplicities) of the spectra $(\nu_i^\epsilon)_{i \ge 1}$ and $(\mu_i^\epsilon)_{i \ge 0}$, where

* $\nu_i^\epsilon$ corresponds to an eigenfunction for which the $L^2$ norm is concentrated on the “helmet straps” $H_{\theta,\epsilon} \setminus B_\theta$. One expects that $\nu_i^\epsilon = \nu_i + o(1)$ as $\epsilon \to 0$, where the $(\nu_i)_{i \ge 1}$ are the solutions to the 1-dimensional Dirichlet eigenvalue problem on the graph

$$G := (\{x_1 = 0\} \cup \{x_2 = 0\}) \setminus B_\theta$$

on the sphere, with boundary $\partial G := (\{x_1 = 0\} \cup \{x_2 = 0\}) \cap \partial B_\theta$, and where the eigenvalue problem is

$$\begin{cases} -\Delta_G w_i = \nu_i w_i & \text{in } G, \\ w_i = 0 & \text{on } \partial G, \end{cases}$$

and the first equation is understood as $-\partial_\tau^2 w_i = \nu_i w_i$ along the edges of the graph, with a Kirchhoff condition at the (only) vertex. Explicitly, one can show the first eigenvalue is

$$\nu_1 = \frac{\pi^2}{4(\pi - \theta)^2}$$

associated to the eigenfunction

$$w_1 \left( \pi + \tau, \frac{k\pi}{2} \right) = \cos \left( \frac{\pi\tau}{2(\pi - \theta)} \right), \quad \tau \in [-(\pi - \theta), \pi - \theta], \quad k = 0, 1.$$

In particular $\nu_1 > \mu_1(B_\theta)$ if and only if $(\pi - \theta)\sqrt{\mu_1(B_\theta)} < \frac{\pi}{2}$. Letting $\Theta'$ be the value for which $(\pi - \theta)\sqrt{\mu_1(B_\theta)} = \frac{\pi}{2}$, it may be checked that $\Theta' \simeq 0.61\pi$ and that the condition $(\pi - \theta)\sqrt{\mu_1(B_\theta)} < \frac{\pi}{2}$ holds for all $\theta > \Theta'$. Accordingly, from now on we suppose $\theta \in (\Theta', \pi)$.

* $\mu_i^\epsilon$ corresponds to a perturbation of $\mu_i(B_\theta)$, and is associated to an eigenfunction $u_i^\epsilon$ of $H_{\theta,\epsilon}$ that is an $L^2$-perturbation of an eigenfunction associated to $\mu_i(B_\theta)$, denoted $u_i$. In this case an asymptotic of $\mu_i^\epsilon$ is

$$\mu_i^\epsilon = \mu_i(B_\theta) + \epsilon \int_G (|\nabla_G v_i|^2 - \mu_i(B_\theta)|v_i|^2) + o(\epsilon),$$
where $v_i$ is defined as the solution to

$$ \begin{cases} -\Delta_G v_i = \mu_i(B_\theta)v_i & \text{in } G, \\ v_i = u_i & \text{on } \partial G, \end{cases} $$

with again a Kirchhoff condition at the vertex. This $v_i$ exists and is unique only when $\mu_i(\theta) \neq \nu_j$ for all $j$, which is true for $i = 1$ when we suppose $\theta > \Theta'$.

Under the hypothesis that $\theta > \Theta'$, we know that for a small enough $\epsilon$ the first eigenvalue of $H_{\theta,\epsilon}$ is of the form $\mu_1^\epsilon$ as described above. For $i = 1$, the Neumann eigenvalue $\mu_1(B_\theta)$ is associated to the eigenfunctions

$$ u_1(t, \phi) = g(t) \cos \phi, \quad u_2(t, \phi) = g(t) \sin \phi, $$

where $t$ is the angle or geodesic distance from the north pole and $\phi$ is the angle around the line of latitude. We may compute directly that the first two eigenfunctions on the graph concentrate on the two “straps”, respectively:

$$ v_1(\pi + \tau, 0) = - \frac{\sin(\sqrt{\mu_1(B_\theta)} \tau)}{\sin(\sqrt{\mu_1(B_\theta)}(\pi - \theta))}, \quad v_1(\pi + \tau, \pm \frac{\pi}{2}) = 0, $$

for $\tau \in [-(\pi - \theta), \pi - \theta]$, and similarly

$$ v_2(\pi + \tau, 0) = 0, \quad v_2(\pi + \tau, \pm \frac{\pi}{2}) = - \frac{\sin(\sqrt{\mu_1(B_\theta)} \tau)}{\sin(\sqrt{\mu_1(B_\theta)}(\pi - \theta))} g(\theta). $$

So for $i = 1, 2$ we have

$$ \begin{aligned} & \int_G (|\nabla_G v_i|^2 - \mu_i(B_\theta) v_i^2) d\mathcal{H}^1 \\ &= \frac{g(\theta)^2 \mu_1(B_\theta)}{\sin^2(\sqrt{\mu_1(B_\theta)}(\pi - \theta))} \int_{-(\pi-\theta)}^{\pi-\theta} \left( \cos^2(\sqrt{\mu_1(B_\theta)}\tau) - \sin^2(\sqrt{\mu_1(B_\theta)}\tau) \right) d\tau \\ &= 2g(\theta)^2 \sqrt{\mu_1(B_\theta)} \cot(\sqrt{\mu_1(B_\theta)}(\pi - \theta)). \end{aligned} $$

Hence

$$ \mu_1(H_{\theta,\epsilon}) = \mu_1(B_\theta) + 4\epsilon g(\theta)^2 \sqrt{\mu_1(B_\theta)} \cot(\sqrt{\mu_1(B_\theta)}(\pi - \theta)) + o(\epsilon). $$

As a consequence, for any $\theta > \Theta'$ we obtain $\mu_1(H_{\theta,\epsilon}) > \mu_1(B_{\theta(\epsilon)})$ for a small enough $\epsilon$ as soon as the inequality

$$ \sqrt{\mu_1(B_\theta)} \cot \left( \sqrt{\mu_1(B_\theta)} (\pi - \theta) \right) > (\pi - \theta) \left( \frac{1}{\sin^2 \theta} - \mu_1(B_\theta) \right) $$

is verified. When $\theta \in (\Theta', \Theta]$ the left-hand side is positive and the right-hand side is nonpositive, so that the inequality holds on $(\Theta', \Theta'')$, for some $\Theta''$ slightly larger than $\Theta$. Numerics suggest that the inequality holds for all $\theta \in (\Theta', \pi)$, so that helmet counterexamples can be constructed for all such $\theta$, but we do not prove this claim.

## 6. NO ASYMPTOTIC COUNTEREXAMPLE WITH A SINGLE HOLE

Here we are interested in the possibility of finding a simply connected counterexample. Precisely, consider a set $\Omega^\epsilon$ of the form

$$ \Omega^\epsilon = \mathbb{S}^2 \setminus \psi(\epsilon\omega), $$

where $\omega \subset \mathbb{R}^2$ is a smooth simply connected set with finite measure, where $\epsilon > 0$ is small and $\psi: \mathbb{R}^2 \to \mathbb{S}^2$ is the stereographic projection sending the origin 0 to the south pole $-e_3$, with Jacobian 1 at 0.

Is it possible to obtain, for a fixed $\omega$ and arbitrarily small $\epsilon$, a counterexample to the maximality of the spherical cap among domains of the same area? We prove that such a counterexample works neither
for the functional $\mu_1$ nor for $\mu_1 + \mu_2$ or $(\frac{1}{\mu_1} + \frac{1}{\mu_2})^{-1}$. Our result does not imply that the spherical cap maximizes $\mu_1$ among “large” spherical domains, since we do not explore general sequences of domains expanding to fill the sphere, merely domains that are given by the complement of a dilation of a fixed simply connected set $\omega$.

Since $\mu_1 = \mu_2 = \mu_3 = 2$ on the full sphere $\mathbb{S}^2$, with eigenfunctions given by the $L^2$-normalized coordinate functions $u_i = \sqrt{3/4\pi} x_i$ for $i = 1, 2, 3$, Proposition 3 with $M = \Omega = \mathbb{S}^2$ implies that the eigenvalues of the sphere with hole $\psi(\epsilon\omega)$ are of the form

$$ \mu_i(\Omega^\epsilon) = 2 + \epsilon^2 \kappa_i + o(\epsilon^2) $$

where $\kappa_i$ is the $i$-th eigenvalue of the matrix

$$ (2|\omega|v_i(0)v_j(0) - \langle M_\omega \nabla v_i(0), \nabla v_j(0) \rangle)_{1\le i,j\le 3} $$

and where we wrote $v_i = u_i \circ \psi$. In this case we have $v_1(0) = v_2(0) = 0$, $v_3(0) = -\sqrt{3/4\pi}$, and $\nabla v_3(0) = 0$, $\nabla v_i(0) = \sqrt{3/4\pi} e_i$ for $i = 1, 2$. Thus this matrix has the form

$$ \frac{3}{4\pi} \begin{pmatrix} -(M_\omega)_{1,1} & -(M_\omega)_{1,2} & 0 \\ -(M_\omega)_{2,1} & -(M_\omega)_{2,2} & 0 \\ 0 & 0 & 2|\omega| \end{pmatrix}. $$

So, we obtain

$$ \mu_1(\Omega^\epsilon) + \mu_2(\Omega^\epsilon) = 4 - \frac{3\mathrm{Tr}(M_\omega)}{4\pi}\epsilon^2 + o(\epsilon^2), \\ \mu_3(\Omega^\epsilon) = 2 + \frac{3|\omega|}{2\pi}\epsilon^2 + o(\epsilon^2). $$

We recall that using the formula (4), we have $\mathrm{Tr}(M_\omega) = 4\pi \mathrm{cap}[\omega]^2$. As a consequence

$$ \mu_1(\Omega^\epsilon) + \mu_2(\Omega^\epsilon) = 4 - 3 \mathrm{cap}[\omega]^2 \epsilon^2 + o(\epsilon^2). $$

As observed above, the disk has the smallest capacity among sets of given measure. Hence in the limit as $\epsilon \to 0$, the sum $\mu_1(\Omega^\epsilon) + \mu_2(\Omega^\epsilon)$ of the first two eigenvalues is smaller than the corresponding sum when $\omega$ is a disk, that is, when $\Omega^\epsilon$ is a spherical cap. The same conclusion follows for the first eigenvalue $\mu_1(\Omega^\epsilon)$, since $\mu_1 = \mu_2$ for the spherical cap. Analogous computations may be done for $(\frac{1}{\mu_1} + \frac{1}{\mu_2})^{-1}$.

## 7. PROOF OF FORMULAS (3), (7), (8), AND LEMMA 4, AND HARMONIC EXTENSION TO THE HOLES

Here we prove some results stated in Section 3. Consider the setting and notations of Proposition 3.

**Proof of formula (3).** The harmonic vector field $W_n$ vanishes at infinity and, after transforming infinity to the origin by the Kelvin inversion, we know $W_n$ is locally linear. That is,

$$ W_n(x) = \frac{C_n x}{2\pi |x|^2} + \mathcal{O}_{|x|\to\infty} \left( \frac{1}{|x|^2} \right) \qquad (20) $$

for some $2 \times 2$ matrix $C_n$. We will show for formula (3) that $C_n = M_n$.

From now on, we drop the subscript $n$ on the vector field $W$, the matrices $C$ and $M$, and the set $\omega$. The goal is to relate the matrix entries by $M_{p,q} = C_{p,q}$. Let $R \gg 1$. Writing $W_p$ for the $p$-th component of
$W$, we find

$$
\begin{align*}
& \delta_{p,q} |\omega| + \int_{D_R \setminus \omega} \nabla W_p \cdot \nabla W_q \\
&= \int_{\omega} \nabla \cdot (x_q e_p) + \int_{D_R \setminus \omega} \nabla \cdot (W_p \nabla W_q) \\
&= \int_{\partial \omega} (x_q \nu_p - W_p \partial_{\nu} W_q) + \int_{\partial D_R} W_p \partial_{\nu} W_q \quad \text{by the divergence theorem} \\
&= \int_{\partial \omega} (-x_q \partial_{\nu} W_p + W_p \nu_q) + \mathcal{O}_{|x| \to \infty} \left( \frac{1}{R^2} \right)
\end{align*}
$$

by the boundary condition $\partial_\nu W_n = -\nu$ on $\partial\omega$. Applying Green's theorem and recalling the convention that the normal vector points outward from both $\omega$ and the disk $D_R$, we see

$$
\begin{align*}
& \delta_{p,q} |\omega| + \int_{D_R \setminus \omega} \nabla W_p \cdot \nabla W_q \\
&= \int_{\partial D_R} (-x_q \partial_\nu W_p + W_p \nu_q) + \int_{D_R \setminus \omega} (x_q \Delta W_p - W_p \Delta x_q) + \mathcal{O}_{|x| \to \infty} \left( \frac{1}{R^2} \right).
\end{align*}
$$

The second integral on the right vanishes since both $W_p$ and $x_q$ are harmonic. For the first integral on the right, observe that the normal vector $\nu = x/|x|$ is radially outward on the circle, and use (20) to compute that

$$
W_p = \frac{C_p \cdot x}{2\pi|x|^2} + \mathcal{O}_{|x| \to \infty} \left( \frac{1}{|x|^2} \right), \quad -x_q \partial_\nu W_p = x_q \frac{C_p \cdot x}{2\pi|x|^3} + \mathcal{O}_{|x| \to \infty} \left( \frac{1}{|x|^2} \right),
$$

where $C_p$ is the $p$-th row of $C$. It is now straightforward to integrate over the circle $\partial D_R$, deducing that

$$
\delta_{p,q} |\omega| + \int_{D_R \setminus \omega} \nabla W_p \cdot \nabla W_q = C_{p,q} + \mathcal{O}_{|x| \to \infty} \left( \frac{1}{R} \right).
$$

Letting $R \to \infty$ shows that $M_{p,q} = C_{p,q}$.

**Proof of formula (7).** The goal is to show

$$
- \int_{\Omega} (K u_i) u_j + \int_{\partial \Omega} \langle \eta, \nabla u_i \rangle u_j = \int_{\partial \Omega} (\langle \nabla u_i, \nabla u_j \rangle - \mu u_i u_j) \langle \phi, \nu \rangle
$$

where $K$ and $\eta$ are defined at the beginning of Section 3.2. Similar computations may be found in [8, Th. 5.7.2.]. In our case we rely on the computations in Fall and Weth [5, Sections 2 and 3].

Since the metric varies in this section, we make its dependence appear explicitly: $\langle \cdot, \cdot \rangle_t$ is the scalar product associated to $g_t$ and $\nabla_t$ is the gradient, so that $\nabla_t \cdot$ is the divergence. Also, $\Delta_t = \nabla_t \cdot \nabla_t$ is the Laplacian, $dm_t$ is the area measure associated to $g_t$, and $\nu_t$ is the outward normal unit vector on $\partial\Omega$, and $ds_t$ is the boundary measure on $\partial\Omega$. We also denote by $D$ the Levi-Civita connection associated to $g$. Green's theorem says

$$
- \int_{\Omega} u_j \Delta_t u_i \, dm_t + \int_{\partial \Omega} u_j \, du_i(\nu_t) \, ds_t = \int_{\Omega} \langle \nabla_t u_i, \nabla_t u_j \rangle_t \, dm_t.
$$

A commutator identity will be useful. For vector fields $\phi$ and $\psi$ and functions $h$, one has

$$
\langle [\psi, \phi], \nabla h \rangle = \psi \langle \phi, \nabla h \rangle - \phi \langle \psi, \nabla h \rangle.
$$

When $\psi$ is chosen to be a gradient field, say $\psi = \nabla f$, the identity implies that

$$
\langle [\nabla f, \phi], \nabla h \rangle = \langle \nabla f, \nabla \langle \phi, \nabla h \rangle \rangle - \phi \langle \nabla f, \nabla h \rangle. \quad (21)
$$
We have

$$
\begin{align*}
& \left. \frac{d}{dt} \right|_{t=0} \int_{\Omega} \langle \nabla_t u_i, \nabla_t u_j \rangle_t dm_t \\
&= \int_{\Omega} \left( \langle \nabla u_i, \nabla u_j \rangle \nabla \cdot \phi - \left\{ \langle D_{\nabla u_i} \phi, \nabla u_j \rangle + \langle D_{\nabla u_j} \phi, \nabla u_i \rangle \right\} \right) dm \\
&\qquad \text{by using Fall and Weth [5, formulas (7) and (14)]} \\
&= \int_{\partial \Omega} \langle \nabla u_i, \nabla u_j \rangle \langle \phi, \nu \rangle ds - \int_{\Omega} \langle \phi, \nabla \langle \nabla u_i, \nabla u_j \rangle \rangle dm \quad \text{by parts} \\
&\qquad - \int_{\Omega} \left\{ \langle D_{\phi} \nabla u_i, \nabla u_j \rangle + \langle D_{\phi} \nabla u_j, \nabla u_i \rangle + \langle [\nabla u_i, \phi], \nabla u_j \rangle + \langle [\nabla u_j, \phi], \nabla u_i \rangle \right\} dm \\
&\qquad \text{by the torsion-free property of the connection} \\
&= \int_{\partial \Omega} \langle \nabla u_i, \nabla u_j \rangle \langle \phi, \nu \rangle ds - \int_{\Omega} \phi \langle \nabla u_i, \nabla u_j \rangle dm \\
&\qquad - \int_{\Omega} \left\{ \phi \langle \nabla u_i, \nabla u_j \rangle + \langle \nabla u_i, \nabla \langle \phi, \nabla u_j \rangle \rangle + \langle \nabla u_j, \nabla \langle \phi, \nabla u_i \rangle \rangle - 2\phi \langle \nabla u_i, \nabla u_j \rangle \right\} dm
\end{align*}
$$

by the Leibniz property of the connection, and the commutator identity (21). Cancelling terms, we find

$$
\begin{align*}
& \left. \frac{d}{dt} \right|_{t=0} \int_{\Omega} \langle \nabla_t u_i, \nabla_t u_j \rangle_t dm_t \\
&= \int_{\partial \Omega} \langle \nabla u_i, \nabla u_j \rangle \langle \phi, \nu \rangle ds - \int_{\Omega} \{ \langle \nabla u_i, \nabla \langle \phi, \nabla u_j \rangle \rangle + \langle \nabla u_j, \nabla \langle \phi, \nabla u_i \rangle \rangle \} dm \\
&= \int_{\partial \Omega} \langle \nabla u_i, \nabla u_j \rangle \langle \phi, \nu \rangle ds + \int_{\Omega} \{ \langle \phi, \nabla u_j \rangle \Delta u_i + \langle \phi, \nabla u_i \rangle \Delta u_j \} dm \\
&\qquad \text{by Green's theorem since } \langle \nabla u_i, \nu \rangle = 0 \text{ and } \langle \nabla u_j, \nu \rangle = 0 \text{ on } \partial \Omega \\
&= \int_{\partial \Omega} \langle \nabla u_i, \nabla u_j \rangle \langle \phi, \nu \rangle ds - \mu \int_{\Omega} \langle \phi, \nabla (u_i u_j) \rangle dm
\end{align*}
$$

since $\Delta u_i = -\mu u_i$ and $\Delta u_j = -\mu u_j$.

Hence from the divergence theorem we conclude

$$
\left. \frac{d}{dt} \right|_{t=0} \int_{\Omega} \langle \nabla_t u_i, \nabla_t u_j \rangle_t dm_t = \int_{\partial \Omega} \{ \langle \nabla u_i, \nabla u_j \rangle - \mu u_i u_j \} \langle \phi, \nu \rangle ds + \mu \int_{\Omega} u_i u_j (\nabla \cdot \phi) dm. \quad (22)
$$

Since $Kf$ and $\eta$ are respectively the derivatives at $t=0$ of $\Delta_t f$ and $\nu_t$, we have

$$
\begin{align*}
& - \int_{\Omega} u_j K u_i \, dm + \int_{\partial \Omega} u_j du_i(\eta) \, ds \\
&= \left. \frac{d}{dt} \right|_{t=0} \left[ - \int_{\Omega} u_j \Delta_t u_i \, dm_t + \int_{\partial \Omega} u_j du_i(\nu_t) \, ds_t + \int_{\Omega} u_j \Delta u_i \, dm_t - \int_{\partial \Omega} u_j du_i(\nu) \, ds_t \right] \\
&= \left. \frac{d}{dt} \right|_{t=0} \left[ \int_{\Omega} \langle \nabla_t u_i, \nabla_t u_j \rangle_t \, dm_t + \int_{\Omega} u_j \Delta u_i \, dm_t \right] && \text{since } du_i(\nu) = 0 \\
&= \int_{\partial \Omega} \{ \langle \nabla u_i, \nabla u_j \rangle - \mu u_i u_j \} \langle \phi, \nu \rangle \, ds
\end{align*}
$$

by (22), where in the last line we used again that $\Delta u_i = -\mu u_i$, and that $\left.\frac{d}{dt}\right|_{t=0} dm_t = (\nabla \cdot \phi) dm$ (see [5, formula (7)]).
**Proof of formula (8).** Recall $F_{n,i}$ was defined in (5). Write $\tilde{u}_j = u_j \circ \psi_n$ in the proof below. Since $F_{n,i}$ is supported in $B_r$, the left side of (8) equals

$$ 
\begin{aligned}
-\int_{\psi_n(B_r)} (F_{n,i} \circ \psi_n^{-1}) u_j &= -\int_{B_r} \rho_n F_{n,i} \tilde{u}_j \\
&= -\int_{B_r \setminus \{0\}} \tilde{u}_j(\Delta + \mu \rho_n)g && \text{by (5), with } g \text{ as defined below,} \\
&= -\int_{B_r \setminus \{0\}} (\tilde{u}_j \Delta g - g \Delta \tilde{u}_j) && \text{since } \Delta \tilde{u}_j = -\mu \rho_n \tilde{u}_j, \\
&= \lim_{s \to 0} \int_{\partial B_s} (\tilde{u}_j \partial_r g - g \partial_r \tilde{u}_j) && \text{by Green's theorem,}
\end{aligned}
$$

where the compactly supported function $g$ on $B_r$ is

$$ g(x) = \zeta(x) \left( \frac{\langle x, M_n \nabla \tilde{u}_i(0) \rangle}{2\pi|x|^2} + \frac{\rho_n(0)\mu|\omega_n|\tilde{u}_i(0)}{2\pi} \log|x| \right). $$

We consider each boundary integral in turn. First, the definition of $g$ implies that

$$ 
\begin{aligned}
\int_{\partial B_s} g \partial_r \tilde{u}_j &= \int_{\partial B_s} \left( \frac{1}{2\pi s} \langle x/|x|, M_n \nabla \tilde{u}_i(0) \rangle \langle x/|x|, \nabla \tilde{u}_j(0) \rangle + o\left(\frac{1}{s}\right) \right) \\
&\to \frac{1}{2} \langle M_n \nabla \tilde{u}_i(0), \nabla \tilde{u}_j(0) \rangle
\end{aligned}
$$

as $s \to 0$, by parameterizing the circle with $x = (s \cos \phi, s \sin \phi)$. Next,

$$ 
\begin{aligned}
\tilde{u}_j \partial_r g &= (\tilde{u}_j(0) + s \langle x/|x|, \nabla \tilde{u}_j(0) \rangle + o(s)) \\
&\quad \times \left( -\frac{1}{2\pi s^2} \langle x/|x|, M_n \nabla \tilde{u}_i(0) \rangle + \frac{1}{2\pi s} \rho_n(0) \mu |\omega_n| \tilde{u}_i(0) \right) \\
&= -\frac{\tilde{u}_j(0)}{2\pi s^2} \langle x/|x|, M_n \nabla \tilde{u}_i(0) \rangle - \frac{1}{2\pi s} \langle x/|x|, \nabla \tilde{u}_j(0) \rangle \langle x/|x|, M_n \nabla \tilde{u}_i(0) \rangle \\
&\quad + \frac{1}{2\pi s} \rho_n(0) \mu |\omega_n| \tilde{u}_i(0) \tilde{u}_j(0) + o\left(\frac{1}{s}\right).
\end{aligned}
$$

Parameterizing the circle and integrating yields that

$$ \int_{\partial B_s} \tilde{u}_j \partial_r g \to 0 - \frac{1}{2} \langle \nabla \tilde{u}_j(0), M_n \nabla \tilde{u}_i(0) \rangle + \rho_n(0) \mu |\omega_n| \tilde{u}_i(0) \tilde{u}_j(0). $$

Combining the two boundary limits completes the proof of formula (8).

**Proof of Lemma 4.** To lighten the notation, we drop the index “i” in all that follows. Consider $z \in H^1(\Omega^\epsilon)$ with unit norm.

First we separate the estimate in Lemma 4 into two parts, by Green's theorem:

$$ \int_{\Omega^\epsilon} (\mu_\epsilon u_\epsilon z - \langle \nabla u_\epsilon, \nabla z \rangle) = \underbrace{\int_{\Omega^\epsilon} z(\Delta + \mu_\epsilon)u_\epsilon}_{\text{regular term}} - \int_{\partial\phi^{\epsilon^2}(\Omega)} z\partial_\nu u_\epsilon + \sum_{n=1}^N \underbrace{\int_{\partial(\epsilon\omega_n)} (z \circ \psi_n) \partial_\nu(u_\epsilon \circ \psi_n)}_{\text{boundary term}}, \quad (23) $$

noting in the boundary term that $\psi_n$ is conformal and so the change of variable factor for the normal derivative compensates exactly for the tangential change of variable factor due to the arclength element.
*Boundary term.* Fix $n \in \{1, \dots, N\}$ and write $\tilde{f} = f \circ \psi_n$ for $f = u, u_\epsilon, v, z$, so that the “tilde” version of each function is defined either on $B_r$ or $B_r \setminus \epsilon\omega_n$, as the case may be. The boundary term is $\int_{\partial(\epsilon\omega_n)} \tilde{z} \partial_\nu \tilde{u}_\epsilon$.

In terms of these functions, on $B_r \setminus \epsilon\omega_n$ the ansatz (9) becomes

$$ 
\begin{aligned}
\tilde{u}_\epsilon(x) = \tilde{u}(x) &+ \epsilon^2 \tilde{v}(x) + \epsilon\zeta(x) \langle W_n(x/\epsilon), \nabla \tilde{u}(0) \rangle \\
&+ \epsilon^2 \zeta(x) \left( -L_n(x/\epsilon) : D^2 \tilde{u}(0) + \frac{\rho_n(0)\mu|\omega_n|\tilde{u}(0)}{2\pi} \log \epsilon \right) \\
\mu_\epsilon = \mu &+ \epsilon^2 \kappa
\end{aligned}
$$

where $-\Delta \tilde{u} = \rho_n \mu \tilde{u}$ and $\tilde{v}$ verifies

$$ -\Delta \tilde{v} - \rho_n \mu \tilde{v} = \rho_n \kappa \tilde{u} + \rho_n F_n $$

in $B_r$, where we recall that $F_n$ is defined in equation (5). The normal derivative of $\tilde{u}_\epsilon$ can be computed from the ansatz, for $x \in \partial(\epsilon\omega_n)$:

$$ 
\begin{aligned}
\partial_\nu \tilde{u}_\epsilon(x) &= \langle \nu(x), \nabla \tilde{u}(x) \rangle + \langle (\partial_\nu W_n)(x/\epsilon), \nabla \tilde{u}(0) \rangle - \epsilon \partial_\nu L_n(x/\epsilon) : D^2 \tilde{u}(0) + \epsilon^2 \partial_\nu \tilde{v}(x) \\
&= \langle \nu(x), \nabla \tilde{u}(x) - \nabla \tilde{u}(0) - D^2 \tilde{u}(0)x \rangle + \epsilon^2 \partial_\nu \tilde{v}(x)
\end{aligned}
$$

by using the boundary conditions $\partial_\nu W_n(x) = -\nu(x)$ and $\partial_\nu L_n(x) = x \otimes \nu(x)$ for $x \in \partial\omega_n$, together with the fact that the unit normal vector remains invariant under scaling $(\nu_{\omega_n}(x/\epsilon) = \nu_{\epsilon\omega_n}(x)$ when $x \in \partial(\epsilon\omega_n)$).

Regard $\tilde{z} = z \circ \psi_n$ as being extended harmonically into $\epsilon\omega_n$. The extended $\tilde{z}$ is still bounded (independent of $\epsilon$) in $H^1(B_r)$, by Lemma 7 below and using the scale-invariance of the Dirichlet energy. As a consequence, $\tilde{z}$ is also bounded in $L^p(B_r)$ for all $p < \infty$, and in particular is bounded in $L^{12}(B_r)$.

From the definition of $v$ in (6), we see $\Delta \tilde{v} \in L^{3/2}(B_r)$ and so elliptic regularity implies $\tilde{v} \in W^{2,3/2}(B_{r/2})$ and so $\nabla \tilde{v} \in L^4(B_{r/2})$. (In fact, $\Delta \tilde{v} \in L^p$ for any $p \in [1, 2)$ and so $\tilde{v} \in W^{2,p}(B_{r/2})$ and hence $\nabla v \in L^q(B_{r/2})$ for all $q < \infty$.) We decompose

$$ 
\begin{aligned}
\int_{\partial(\epsilon\omega_n)} \tilde{z} \partial_\nu \tilde{u}_\epsilon &= \int_{\partial(\epsilon\omega_n)} \tilde{z} \langle \nabla \tilde{u} - \nabla \tilde{u}(0) - D^2 \tilde{u}(0)x + \epsilon^2 \nabla \tilde{v}, \nu \rangle \\
&= \int_{\epsilon\omega_n} \nabla \tilde{z} \cdot (\nabla \tilde{u} - \nabla \tilde{u}(0) - D^2 \tilde{u}(0)x + \epsilon^2 \nabla \tilde{v}) \\
&\quad + \int_{\epsilon\omega_n} \tilde{z} \nabla \cdot (\nabla \tilde{u} - \nabla \tilde{u}(0) - D^2 \tilde{u}(0)x + \epsilon^2 \nabla \tilde{v})
\end{aligned}
$$

by the divergence theorem. In the last line, note that $\nabla \cdot D^2 \tilde{u}(0)x = \Delta \tilde{u}(0)$. We bound terms separately as follows, using Hölder's inequality and that the area of $\epsilon\omega_n$ equals $\epsilon^2|\omega_n|$:

$$ 
\begin{aligned}
\left| \int_{\epsilon\omega_n} \nabla \tilde{z} \cdot (\nabla \tilde{u} - \nabla \tilde{u}(0) - D^2 \tilde{u}(0)x) \right| &\le C \|\tilde{z}\|_{H^1(B_r)} |\epsilon\omega_n|^{\frac{1}{2}} \sup_{\epsilon\omega_n} |\nabla \tilde{u} - \nabla \tilde{u}(0) - D^2 \tilde{u}(0)x| \le C' \epsilon^3 \\
\left| \int_{\epsilon\omega_n} \tilde{z} (\Delta \tilde{u} - \Delta \tilde{u}(0)) \right| &\le C \|\tilde{z}\|_{L^4(B_r)} |\epsilon\omega_n|^{1-\frac{1}{4}} \sup_{\epsilon\omega_n} |\Delta \tilde{u} - \Delta \tilde{u}(0)| \le C' \epsilon^{2+\frac{1}{2}} \\
\left| \int_{\epsilon\omega_n} \nabla \tilde{z} \cdot \epsilon^2 \nabla \tilde{v} \right| &\le C \epsilon^2 \|\tilde{z}\|_{H^1(B_r)} \|\nabla \tilde{v}\|_{L^4(B_{r/2})} |\epsilon\omega_n|^{\frac{1}{4}} \le C' \epsilon^{2+\frac{1}{2}} \\
\left| \int_{\epsilon\omega_n} \epsilon^2 \tilde{z} \Delta \tilde{v} \right| &\le C \epsilon^2 \|\tilde{z}\|_{L^{12}(B_r)} \|\Delta \tilde{v}\|_{L^{3/2}(B_r)} |\epsilon\omega_n|^{\frac{1}{4}} \le C' \epsilon^{2+\frac{1}{2}}
\end{aligned}
$$

Combining the estimates above, we see that the boundary term (near the n-th hole) is bounded by $C\epsilon^{2+\frac{1}{2}}$, as desired for Lemma 4.

**Regular term 1.** We start by considering the first regular term $\int_{\Omega^\epsilon} z(\Delta + \mu_\epsilon)u_\epsilon$, splitting it into a bulk term and terms near the holes:

$$ \int_{\Omega^\epsilon} z(\Delta + \mu_\epsilon)u_\epsilon = \int_{\Omega^\epsilon \setminus \cup_n \psi_n(B_r)} z(\Delta + \mu_\epsilon)u_\epsilon + \sum_{n=1}^N \int_{B_r \setminus \epsilon\omega_n} \tilde{z}(\Delta \tilde{u}_\epsilon + \rho\mu_\epsilon \tilde{u}_\epsilon). \quad (24) $$
We start with the term near the n-th hole. Let us give a few preliminary estimates on the exterior problems for $W_n$ and $L_n$ in Section 3: there is a constant $C > 0$ such that for all $x \in \mathbb{R}^2 \setminus \omega_n$,

$$ \left| W_n(x) - \frac{M_n x}{2\pi |x|^2} \right| \le \frac{C}{|x|^2}, \quad \left| \partial_{x_p} W_n(x) - \partial_{x_p} \frac{M_n x}{2\pi |x|^2} \right| \le \frac{C}{|x|^3}, $$

and similarly

$$ \left| L_n(x) - \frac{|\omega_n|}{2\pi} (\log|x|) I_2 \right| \le \frac{C}{|x|}, \quad \left| \partial_{x_p} L_n(x) - \frac{|\omega_n|}{2\pi} (\partial_{x_p} \log|x|) I_2 \right| \le \frac{C}{|x|^2}, $$

for $p = 1, 2$.

Formulas (5), (6) and (9) give equations satisfied by the functions $F_n$, $\tilde{v} = v \circ \psi_n$ and $\tilde{u}_\epsilon = u_\epsilon \circ \psi_n$ when $x \in B_r \setminus \epsilon\omega_n$, and by using those equations we may compute

$$ \begin{aligned} & (\Delta + \rho\mu_\epsilon)\tilde{u}_\epsilon(x) \\ &= (\Delta + \rho\mu)\tilde{u} \\ & \quad + \epsilon(\Delta + \rho\mu) \left[ \zeta(x) \left( W_n(x/\epsilon) \cdot \nabla\tilde{u}(0) - \frac{\langle x/\epsilon, M_n \nabla\tilde{u}(0) \rangle}{2\pi|x/\epsilon|^2} \right) \right] \\ & \quad + \epsilon^2(\Delta + \rho\mu) \left[ \zeta(x) \left( -L_n(x/\epsilon) : D^2\tilde{u}(0) - \frac{\rho_n(0)\mu|\omega_n|\tilde{u}(0)}{2\pi} \log|x/\epsilon| \right) \right] \\ & \quad + \epsilon^3\rho\kappa \left[ \zeta(x)W_n(x/\epsilon) \cdot \nabla\tilde{u}(0) + \epsilon\tilde{v}(x) + \epsilon\zeta(x) \left( -L_n(x/\epsilon) : D^2\tilde{u}(0) + \frac{\rho_n(0)\mu|\omega_n|\tilde{u}(0)}{2\pi} \log\epsilon \right) \right]. \end{aligned} $$

Since $(\Delta + \rho\mu)\tilde{u} = 0$ and $W_n(x), L_n(x), x/|x|^2$ and $\log|x|$ are harmonic, the last formula implies that

$$ \begin{aligned} (\Delta + \rho\mu_\epsilon)\tilde{u}_\epsilon(x) ={}& \epsilon(\Delta\zeta + \rho\mu\zeta) \left( W_n(x/\epsilon) \cdot \nabla\tilde{u}(0) - \frac{\langle x/\epsilon, M_n \nabla\tilde{u}(0) \rangle}{2\pi|x/\epsilon|^2} \right) \\ & + 2 \sum_{p=1}^2 \partial_{x_p} \zeta \left( (\partial_{x_p} W_n)(x/\epsilon) \cdot \nabla\tilde{u}(0) - \left. \partial_{x_p} \frac{\langle \cdot, M_n \nabla\tilde{u}(0) \rangle}{2\pi|\cdot|^2} \right|_{x/\epsilon} \right) \\ & - \epsilon^2(\Delta\zeta + \rho\mu\zeta) \left( L_n(x/\epsilon) : D^2\tilde{u}(0) + \frac{\rho_n(0)\mu|\omega_n|\tilde{u}(0)}{2\pi} \log|x/\epsilon| \right) \\ & - 2\epsilon \sum_{p=1}^2 \partial_{x_p} \zeta \left( (\partial_{x_p} L_n)(x/\epsilon) : D^2\tilde{u}(0) + \frac{\rho_n(0)\mu|\omega_n|\tilde{u}(0)}{2\pi} (\partial_{x_p} \log|\cdot|)(x/\epsilon) \right) + \mathcal{O}(\epsilon^3). \end{aligned} $$

Using the error bounds on the exterior problems and noting that $I_2 : D^2\tilde{u}(0) = \Delta\tilde{u}(0)$, we find

$$ |(\Delta + \rho\mu_\epsilon)\tilde{u}_\epsilon(x)| \le C \left( \epsilon^3 + \frac{\epsilon^3}{|x|^2} \right) \le C\epsilon^3 \frac{1}{|x|^2}, \quad x \in B_r \setminus \epsilon\omega_n. $$

As a consequence, and using the fact that $\tilde{z}$ is bounded in $L^4(B_r)$, we get that

$$ \begin{aligned} \int_{B_r \setminus \epsilon\omega_n} \tilde{z}(\Delta \tilde{u}_\epsilon + \rho \mu_\epsilon \tilde{u}_\epsilon) &\le C \epsilon^3 \|\tilde{z}\|_{L^4(B_r)} \left( \int_{B_r \setminus \epsilon\omega_n} \left( \frac{1}{|x|^2} \right)^{4/3} \right)^{3/4} \\ &\le C \epsilon^{2+\frac{1}{2}}. \end{aligned} $$

Next we treat the bulk term in (24), which is $\int_{\Omega^\epsilon \setminus \cup_n \psi_n(B_r)} z(\Delta + \mu_\epsilon)u_\epsilon$. First note that $u_\epsilon$ is bounded in $C^p(\Omega^\epsilon \setminus \cup_n \psi_n(B_r))$ for any $p \in \mathbb{N}$, by the definition (9) of $u_\epsilon$ and elliptic regularity for $v$ (which is valid up to the boundary since $\partial\Omega$ is smooth and the Neumann condition is satisfied).
Denoting by $g_t$ the metric on $\Omega$ obtained by pulling back the area measure of $\phi^t(\Omega)$, and by $\Delta_t$ and $m_t$ the associated Laplacian and area measure, we have

$$ \int_{\Omega^{\epsilon} \setminus \cup_n \psi_n(B_r)} z(\Delta + \mu_{\epsilon}) u_{\epsilon} = \int_{\Omega \setminus \cup_n \psi_n(B_r)} (z \circ \phi^{\epsilon^2})(\Delta_{\epsilon^2} + \mu_{\epsilon}) \hat{u}_{\epsilon} dm_{\epsilon^2} $$

where $\hat{u}_{\epsilon} = u_{\epsilon} \circ \phi^{\epsilon^2}$ was defined in the ansatz (9). Then, in $\Omega \setminus \cup_n \psi_n(B_r)$:

$$ \begin{aligned} (\Delta_{\epsilon^2} + \mu_{\epsilon}) \hat{u}_{\epsilon} &= (\Delta + \mu + \epsilon^2(K + \kappa) + (\Delta_{\epsilon^2} - \Delta - \epsilon^2 K))(u + \epsilon^2 v) && \text{by ansatz (9)} \\ &= (\Delta + \mu)u + \epsilon^2 (\Delta v + \mu v + \kappa u + Ku) \\ &\quad + (\Delta_{\epsilon^2} - \Delta - \epsilon^2 K)(u + \epsilon^2 v) + \epsilon^4(K + \kappa)v. \end{aligned} $$

The first of these two lines vanishes by the definitions of $u$ and $v$, and the second line is of order $O(\epsilon^4)$ (by definition of $K$, $\Delta_{\epsilon^2} - \Delta - \epsilon^2 K$ is a second order differential operator with coefficient bounded uniformly by $C\epsilon^4$).

**Regular term 2.** Finally, we treat similarly the second regular term of equation (23), which is $\int_{\partial\phi^{\epsilon^2}(\Omega)} z \partial_{\nu} u_{\epsilon}$.

Recall from earlier in the paper that $\nu_t = (\phi^t)^* \nu_{\partial\phi^t(\Omega)}$ is the outward normal vector with respect to the metric $g_t$ on $\Omega$. Taking $t = \epsilon^2$, the formula says $\nu_{\epsilon^2} = (\phi^{\epsilon^2})^* \nu$ where $\nu$ is the outward normal vector on $\partial\phi^{\epsilon^2}(\Omega)$. Hence the normal derivative of $u_{\epsilon}$ is $\partial_{\nu} u_{\epsilon} = du_{\epsilon}(\nu) = d\hat{u}_{\epsilon}(\nu_{\epsilon^2})$, using that $u_{\epsilon} \circ \phi^{\epsilon^2} = \hat{u}_{\epsilon}$. Thus

$$ \begin{aligned} \partial_{\nu} u_{\epsilon} &= d\hat{u}_{\epsilon}(\nu_{\epsilon^2}) \\ &= \langle \nabla \hat{u}_{\epsilon}, \nu_{\epsilon^2} \rangle \\ &= \langle \nabla(u + \epsilon^2 v), \nu + \epsilon^2 \eta + (\nu_{\epsilon^2} - \nu - \epsilon^2 \eta) \rangle \\ &= \langle \nabla u, \nu \rangle + \epsilon^2 (\langle \nabla v, \nu \rangle + \langle \nabla u, \eta \rangle) \\ &\quad + \langle \nabla(u + \epsilon^2 v), \nu_{\epsilon^2} - \nu - \epsilon^2 \eta \rangle + \epsilon^4 \langle \nabla v, \eta \rangle. \end{aligned} $$

The first line vanishes by the Neumann boundary conditions on $u$ and $v$ (that is, $\langle \nabla u, \nu \rangle = 0$ and $\langle \nabla v, \nu \rangle = -\langle \nabla u, \eta \rangle$), and the second line is bounded by $O(\epsilon^4)$.

This uniform estimate along with a trace inequality for $z$ on $\partial\Omega^{\epsilon}$ completes the proof that the second regular term of equation (23) is bounded by $C\epsilon^4$.

**Lemma 7** (Harmonic extension across a hole). Let $\omega \subset \mathbb{R}^2$ be a smoothly bounded open set. If $h \in H^1(\mathbb{R}^2)$ is compactly supported and harmonic on $\omega$ then $\int_{\omega} |\nabla h|^2 \le C_{\omega} \int_{\mathbb{R}^2 \setminus \omega} |\nabla h|^2$ for some constant $C_{\omega} > 0$.

**Proof.** Without loss of generality, suppose $\omega$ contains the origin. Write $i(z) = z/|z|^2$ for the inversion map on $\mathbb{R}^2 \cup \{\infty\}$, so that $\omega^* = i(\mathbb{R}^2 \cup \{\infty\} \setminus \bar{\omega})$ is also a smoothly bounded domain that contains the origin. Denote the fractional semi-norm by $[h]_{H^{1/2}(\partial\omega)} = \left( \int_{\partial\omega} \int_{\partial\omega} \frac{|h(x)-h(y)|^2}{|x-y|^2} ds(x)ds(y) \right)^{1/2}$ where $ds$ is the arclength element. Then for some constants $C_1, C_2, C_3 > 0$, we have

$$ \int_{\omega} |\nabla h|^2 \le C_1 [h]_{H^{1/2}(\partial\omega)}^2 \le C_2 [h \circ i]_{H^{1/2}(\partial\omega^*)}^2 \le C_3 \int_{\omega^*} |\nabla (h \circ i)|^2 = C_3 \int_{\mathbb{R}^2 \setminus \omega} |\nabla h|^2 $$

where the first “inverse trace” inequality comes from harmonicity of $h$, the second inequality uses the definition of the semi-norm, the third is the $H^1 \to H^{1/2}$ trace theorem applied to $u - \text{avg}_{\Omega}(u)$ and followed by Poincaré’s inequality, and the final equality relies on invariance of the Dirichlet integral under the conformal map $i$. $\square$

## ACKNOWLEDGEMENTS

Laugesen’s research was supported by grants from the Simons Foundation (964018) and the U.S. National Science Foundation (2246537).
REFERENCES

[1] J. M. Arrieta, Rates of eigenvalues on a dumbbell domain. Simple eigenvalue case. Trans. Amer. Math. Soc. 347 (1995), no. 9, 3503-3531.

[2] M. Ashbaugh and R. Benguria, *Sharp upper bound to the first nonzero Neumann eigenvalue for bounded domains in spaces of constant curvature*, J. London Math. Soc. (2), 52(2) (1995), 402–416.

[3] C. Bandle, *Isoperimetric inequality for some eigenvalues of an inhomogeneous, free membrane*, SIAM J. Appl. Math., 22 (1972), 142–147.

[4] D. Bucur, E. Martinet and M. Nahon, *Sharp inequalities for Neumann eigenvalues on the sphere*, J. Differential Geom., to appear. doi:arXiv:2208.11413

[5] M. M. Fall and T. Weth, *Critical domains for the first nonzero Neumann eigenvalue in Riemannian manifolds*, J. Geom. Anal. 29 (2019), 3221–3247.

[6] V. Felli, L. Liverani and R. Ognibene, *Quantitative spectral stability for the Neumann Laplacian in domains with small holes*, J. Funct. Anal. 288 (2025), Paper No. 110817.

[7] F. Hecht, *New development in freefem++*, J. Numer. Math. 20 (2012), 251–265.

[8] Henrot, A., Pierre, M., *Shape variation and optimization*, 2018.

[9] Y. Kabeya, *Eigenvalues of the Laplace-Beltrami operator under the homogeneous Neumann condition on a large zonal domain in the unit sphere*, Discrete Contin. Dyn. Syst. 40.6 (2020), pp. 3529–3559.

[10] J. J. Langford and R. S. Laugesen, *Maximizers beyond the hemisphere for the second Neumann eigenvalue*, Math. Ann. 386 (2023), 2255–2281.

[11] M. Lanza de Cristoforis, *Simple Neumann eigenvalues for the Laplace operator in a domain with a small hole. A functional analytic approach*. Rev. Mat. Complut. 25 (2012), 369–412.

[12] E. Martinet, *Numerical optimization of Neumann eigenvalues of domains in the sphere*, J. Comput. Phys. 508 (2024), Paper No. 113002, 26 pp.

[13] V. Maz'ya, S. Nazarov and B. Plamenewskij, *Asymptotic Theory of Elliptic Boundary Value Problems in Singularly Perturbed Domains*. Volumes 1 and 2, Oper. Theory Adv. Appl., 111 and 112, Birkhäuser Verlag, Basel, 2000.

[14] S. Nazarov and J. Sokolowski, *Shape sensitivity analysis of eigenvalues revisited*, Control Cybernet. 37 (2008), 999–1012.

[15] S. Ozawa, *Spectra of domains with small spherical Neumann boundary*, J. Fac. Sci. Univ. Tokyo Sect. IA Math. 30 (1983), 259–277.

[16] G. Pólya and G. Szegő, *Isoperimetric Inequalities in Mathematical Physics*, No. 27. Princeton University Press, 1951.

(Dorin Bucur) UNIV. SAVOIE MONT BLANC, CNRS, LAMA, 73000 CHAMBÉRY, FRANCE

Email address, D. Bucur: dorin.bucur@univ-savoie.fr

(Richard Laugesen) DEPARTMENT OF MATHEMATICS, UNIVERSITY OF ILLINOIS, URBANA, IL 61801, U.S.A.

Email address: Laugesen@illinois.edu

(Eloi Martinet) INSTITUT FÜR MATHEMATIK, JMU WÜRZBURG, 97074 WÜRZBURG, GERMANY

Email address: eloi.martinet@uni-wuerzburg.de

(Mickaël Nahon) UNIV. GRENOBLE ALPES, CNRS, GRENOBLE INP, LJK, 38000 GRENOBLE, FRANCE.

Email address: mickael.nahon@univ-grenoble-alpes.fr