# ON THE MEAN FIRST ARRIVAL TIME OF BROWNIAN PARTICLES ON RIEMANNIAN MANIFOLDS

M. NURSULTANOV, L. TZOU, AND J.C. TZOU

**ABSTRACT.** We use geometric microlocal methods to compute an asymptotic expansion of mean first arrival time for Brownian particles on Riemannian manifolds. This approach provides a robust way to treat this problem, which has thus far been limited to very special geometries. This paper can be seen as the Riemannian 3-manifold version of the planar result of [1] and thus enable us to see the full effect of the local extrinsic boundary geometry on the mean arrival time of the Brownian particles. Our approach also connects this question to some of the recent progress on boundary rigidity and integral geometry [23, 20].

## 1. INTRODUCTION

Let $(M, g, \partial M)$ be a compact connected orientable Riemannian manifold with non-empty smooth boundary and without loss of generality we may assume that it is an open subset of an orientable Riemannian manifold $(\tilde{M}, g)$ without boundary oriented by the Riemannian volume form $dvol_g$. Let also $(X_t, \mathbb{P}_x)$ be the Brownian motion on $M$ with initial condition at $x$, that is, the stochastic process generated by the Laplace-Beltrami operator $\Delta_g$ (this article uses the convention $\Delta_g = -d^*d$ with negative spectrum, where $d$ is the exterior derivative). For any $\Gamma \subset \partial M$ open we denote by $\tau_\Gamma$ the first time the Brownian motion $X_t$ hits $\Gamma$, that is

$$ \tau_\Gamma := \inf\{t \ge 0 : X_t \in \Gamma\}. $$

In the case when $\Gamma = \Gamma_{\epsilon,a}$ is a small elliptic window of eccentricity $\sqrt{1-a^2}$ and size $\epsilon \to 0^+$ (to be made precise later), the narrow escape/mean first arrival time problem wishes to derive an asymptotic expansion as $\epsilon \to 0$ for the expected value $\mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0=x]$ of the first arrival time $\tau_{\Gamma_{\epsilon,a}}$ amongst all Brownian particles starting at $x$. Another quantity of interest is the *average* expected value over $M$:

$$ |M|^{-1} \int_M \mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0=x]dvol_g(x). $$

Here $|M|$ denotes the Riemannian volume of $M$ with respect to the metric $g$.

Many problems in cellular biology may be formulated as mean first arrival time problems; a collection of analysis methods, results, applications, and references may be found in [10]. For example, cells have been modelled as simply connected two-dimensional domains with small absorbing windows on the boundary representing ion channels or target binding sites; the quantity sought is then the mean time for a diffusing ion or receptor to exit through an ion channel or reach a binding site [28, 8, 24].

<footnote n="2010">2010 Mathematics Subject Classification. Primary: 58J65 Secondary: 60J65, 58G15, 92C37.</footnote>

<footnote>M. Nursultanov and L. Tzou are partially supported by ARC DP190103302 and ARC DP190103451 during this work. We thank Ben Goldys for the helpful discussion.</footnote>
There has been much progress for this problem in the setting of planar domains, and we refer the readers to [8, 24, 31, 1] and references therein for a complete bibliography. An important contribution was made in the planar case by [1] to introduce rigor into the computation of [24]. The use of layered potential in [1] also cast this problem in the mainstream language of elliptic PDE and facilitates some of the approach we use in this article.

Few results exists for three dimensional domains in $\mathbb{R}^n$ or Riemannian manifolds; see [3, 27, 30, 5] and references therein. The additional difficulties introduced by higher dimension are highlighted in the introduction of [1] and the challenges in geometry are outlined in [30]. In the case when $M$ is a domain in $\mathbb{R}^3$ with Euclidean metric and $\Gamma_{\epsilon,a}$ is a single small disk absorbing window, [27, 30] gave an expansion for the average of the expected first arrival time, averaged over $M$, up to an unspecified $O(1)$ term:

$$ (1.1) \quad |M|^{-1} \int_M \mathbb{E}[\tau_{\Gamma_{\epsilon,a}} | X_0 = x] \mathrm{dvol}_g(x) \sim \frac{|M|}{4\epsilon} \left[ 1 - \frac{\epsilon}{\pi} H \log \epsilon + \mathcal{O}(\epsilon) \right]. $$

Here, $H$ is the mean curvature of the boundary at the center of the absorbing window. The case when $\Gamma_{\epsilon,a}$ is a small elliptic window was also addressed in [27, 30].

When $M$ is a three dimensional ball with multiple circular absorbing windows on the boundary, an expansion capturing the explicit form of the $O(1)$ correction in equation (1.1) in terms of the Neumann Green's function and its regular part was done in [3]. The method of matched asymptotic used there required the explicit computation of the Neumann Green's function, which is only possible in special geometries with high degrees of symmetry/homogeneity. In these results one does not see the full effects of local geometry. This result was also rigorously proved in [2] but with a better estimate for the error term.

In this paper we outline an approach which allows one to derive all the main terms of $\mathbb{E}[\tau_{\Gamma_{\epsilon,a}} | X_0 = x]$ (up to a remainder vanishing as $\epsilon \to 0$) for Riemannian manifolds of dimension three with a multiple number of small absorbing windows which are boundary geodesic balls or ellipses. We will only demonstrate this approach for one absorbing window so as to not obscure the main idea. In the case when the window is a geodesic ball our approach also adapts naturally to Riemannian manifolds of any dimension as the proof of Proposition 1.1 as well as the analysis for inverting a key integral equation on the ball in Section 4 both carry through to higher dimensions.

We discuss briefly here on how to obtain a comprehensive singularity expansion at the boundary for the Neumann Green's function on a Riemannian manifold as the Euclidean case was of interest in [30] and [3]. We will define in Section 3 the Neumann Green's function $G(x, z)$ on $(M, g, \partial M)$ which satisfies

$$ \Delta_g G(x, z) = -\delta_x, \partial_{\nu_z} G(x, z) \Big|_{z \in \partial M} = \frac{-1}{|\partial M|}, \quad \int_{\partial M} G(x, z) \mathrm{dvol}_{\partial M}(z) = 0, $$

where $z \in \partial M \mapsto \nu_z$ is an outward pointing normal vector field and $|\partial M|$ is the area of the boundary.

Singer-Schuss-Holcman in [30] highlighted the difficulty in obtaining a comprehensive singularity expansion of $G(x, z) \mid_{x,z \in \partial M, x \neq z}$ in a neighbourhood of the diagonal $\{x = z\}$ when $M$ is a bounded domain in $\mathbb{R}^n$, but it turns out that even when $M$ is a general Riemannian manifold this question can be treated by the standard pseudodifferential operators approach. We only carry out this calculation in three dimensions
as it pertains to our application. Readers who are interested in the higher dimensional analogue can follow our treatment to carry out the (cumbersome) calculations for themselves:

**Proposition 1.1.** For $x, y \in \partial M$, set $H(x)$ to be the mean curvature of $\partial M$ at $x$, $d_h(x, y)$ the geodesic distance on the boundary given by metric $h := \iota_{\partial M}^*g$, $d_g(x, y)$ the geodesic distance given by the metric $g$, and

$$ \Pi_x(V) := \Pi_x(V, V), \quad V \in T_x \partial M $$

the scalar second fundamental quadratic form (see pages 235 and 381 of [16] for definitions).

i) The map

$$ f \in C^{\infty}(\partial M) \mapsto \left( \int_{\partial M} G(\cdot, y) f(y) d\mathrm{vol}_h(y) \right) \Big|_{\partial M} $$

is well defined and extends to a map from $H^k(\partial M) \to H^{k+1}(\partial M)$ for all $k \in \mathbb{R}$ whose Schwartz kernel we will denote by $G_{\partial M}(x, y) \in \mathcal{D}'(\partial M \times \partial M)$. Here the map $u \in H^1(M) \mapsto u|_{\partial M} \in H^{1/2}(\partial M)$ is the trace map.

ii) There exists an open neighbourhood of the diagonal

$$ \mathrm{Diag} := \{(x, y) \in \partial M \times \partial M \mid x = y\} $$

such that in this neighbourhood, the singularity structure of $G_{\partial M}(x, y)$ is given by:

$$ (1.2) \quad G_{\partial M}(x, y) = \frac{1}{2\pi} d_g(x, y)^{-1} - \frac{1}{4\pi} H(x) \log d_h(y, x) \\ + \frac{1}{16\pi} \left( \Pi_x \left( \frac{\exp_{x;h}^{-1}(y)}{|\exp_{x;h}^{-1}(y)|_h} \right) - \Pi_x \left( \frac{*\exp_{x;h}^{-1}(y)}{|\exp_{x;h}^{-1}(y)|_h} \right) \right) + R(x, y), $$

where $R(\cdot, \cdot) \in C^{0,\mu}(\partial M \times \partial M)$, for all $\mu < 1$, is called the regular part of the Green's function and $*$ is the Hodge-star operator (i.e. rotation by $\pi/2$ on the surface $\partial M$).

We recall the definition of the exponential map. Let $(X, g_0)$ be a geodesically complete manifold. For any $x \in X$ and $V \in T_x X$ there exists a unique geodesic $\gamma_{g_0}(t) = \gamma_{g_0}(t; V)$, defined on $[0, 1]$, such that $\gamma_{g_0}(0) = x$, $\gamma'_g(0) = V$. The exponential map based at $x$ is then a map taking $T_x X \to X$ defined by

$$ \exp_{x;g_0} : V \mapsto \gamma_{g_0}(1; V). $$

Observe that when $M$ is a Euclidean ball the singular term involving the second fundamental forms vanishes due to homogeneity and therefore does not show up. This is consistent with the explicit formula derived in [3].

An explicit formula for the regular part is only possible in special geometries such as the one considered in [3]. However, our approach in arriving at (1.2) also provides a way to numerically compute $R(x, y)$ via a Fredholm integral equation. See Remark 3.3.

We will use the formula in Proposition 1.1 to derive the mean first arrival time of a Brownian particle on a Riemannian manifold with a single absorbing window which is a small geodesic ellipse. As mentioned earlier, our method extends to multiple windows but we present the single window case to simplify notations. We first state the result when the window is a geodesic disk of the boundary $\partial M$ around a fixed point since the statement is cleaner:
**Theorem 1.2.** Let $(M, g, \partial M)$ be a smooth Riemannian manifold of dimension three with boundary and let $|M|$ be its volume.

i) Fix $x^* \in \partial M$ and let $\Gamma_\epsilon$ be a boundary geodesic ball centered at $x^*$ of geodesic radius $\epsilon > 0$. For each $x \notin \Gamma_\epsilon$,

$$ \mathbb{E}[\tau_{\Gamma_\epsilon}|X_0 = x] = F(x) + C_\epsilon - |M|G(x, x^*) + r_\epsilon(x), $$

with $\|r_\epsilon\|_{C^k(K)} \le C_{k,K}\epsilon$ for any integer $k$ and compact set $K \subset \overline{M}$ which does not contain $x^*$. The function $F$ is the unique solution to the boundary value problem

$$ \Delta_g F = -1, \quad \partial_\nu F = -|M|/|\partial M|, \quad \int_{\partial M} F = 0. $$

The constant $C_\epsilon$ is, modulo an error of $O(\epsilon \log \epsilon)$, given by

$$ C_{\epsilon} = \frac{|M|}{4\epsilon} - \frac{1}{4\pi}H(x^*)|M|\log \epsilon + R(x^*, x^*)|M| - F(x^*) - \frac{|M|H(x^*)}{4\pi} \left(2\log 2 - \frac{3}{2}\right), $$

where $R(x^*, x^*)$ is the evaluation at $(x, y) = (x^*, x^*)$ of the kernel $R(x, y)$ in (1.2).

ii) One has that the integral of $\mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0 = x]$ over $M$ satisfies

$$ \int_M \mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0 = x] \mathrm{dvol}_g(x) = \int_M F(x) \mathrm{dvol}_g(x) + C_\epsilon |M| - F(x^*) |M| + O(\epsilon). $$

Theorem 1.2 does not realize the full power of Proposition 1.1 as it does not see the non-homogeneity of the local geometry at $x^*$ (only the mean curvature $H(x^*)$ shows up). This is due to the fact that we are looking at windows which are geodesic balls. If we replace geodesic balls with geodesic *ellipses*, we see that the second fundamental form term in (1.2) contributes to a term in $\mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0 = x]$ which is the *difference* of principal curvatures.

To this end let $E_1(x^*), E_2(x^*) \in T_{x^*} \partial M$ be the unit eigenvectors of the shape operator at $x^*$ corresponding respectively to the principal curvatures $\lambda_1(x^*), \lambda_2(x^*)$. For $1 \ge a > 0$ fixed, let

$$ \Gamma_{\epsilon,a} := \{\exp_{x^*;h}(\epsilon t_1 E_1(x^*) + \epsilon t_2 E_2(x^*)) \mid t_1^2 + a^{-2}t_2^2 \le 1\} $$

be a small geodesic ellipse.

**Theorem 1.3.** Let $(M, g, \partial M)$ be a smooth Riemannian manifold of dimension three with boundary.

i) For each $x \in M \setminus \Gamma_{\epsilon,a}$,

$$ \mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0 = x] = F(x) + C_{\epsilon,a} - |M|G(x, x^*) + r_{\epsilon}(x) $$

with $\|r_\epsilon\|_{C^k(K)} \le C_{k,K}\epsilon$ for any integer $k$ and compact set $K \subset \overline{M}$ which does not contain $x^*$. The function $F$ is the unique solution to the boundary value problem

$$ \Delta_g F = -1, \quad \partial_\nu F = -|M|/|\partial M|, \quad \int_{\partial M} F = 0. $$
The constant $C_{\epsilon,a}$ is given by

$$ (1.3) \qquad \begin{aligned} C_{\epsilon,a} = & \frac{|M|K_a}{4a\epsilon\pi^2} - \frac{1}{4\pi}H(x^*)|M|\log\epsilon + aR(x^*, x^*)|M| - F(x^*) \\ & - \frac{|M|H(x^*)}{16\pi^3} \int_{\mathbb{D}} \frac{1}{(1-|s'|^2)^{1/2}} \int_{\mathbb{D}} \frac{\log((t_1-s_1)^2 + a^2(t_2-s_2)^2)^{1/2}}{(1-|t'|^2)^{1/2}} dt'ds' \\ & + \frac{|M|(\lambda_1 - \lambda_2)}{64\pi^3} \int_{\mathbb{D}} \frac{1}{(1-|s'|^2)^{1/2}} \int_{\mathbb{D}} \frac{(t_1-s_1)^2 - a^2(t_2-s_2)^2}{(t_1-s_1)^2 + a^2(t_2-s_2)^2} \frac{1}{(1-|t'|^2)^{1/2}} dt'ds' \\ & + O(\epsilon \log \epsilon), \end{aligned} $$

where $K_a = \frac{\pi}{2} \int_0^{2\pi} \left( \cos^2 \theta + \frac{\sin^2 \theta}{a^2} \right)^{-1/2} d\theta$ and $\mathbb{D}$ is the two dimensional unit disk centered at the origin.

ii) One has that the integral of $\mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0=x]$ over $M$ satisfies

$$ \int_M \mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0=x] \mathrm{dvol}_g(x) = \int_M F(x) \mathrm{dvol}_g + C_{\epsilon,a}|M| - F(x^*)|M| + O(\epsilon). $$

Note that while the dependence on the eccentricity of the ellipse is hidden in the integrals, the dependence on the difference of the principal curvatures $(\lambda_1 - \lambda_2)$ is easy to see in this formula. The integral which multiplies $(\lambda_1 - \lambda_2)$ turns out to vanish when $a=1$ which makes the above result consistent with Theorem 1.2.

The fact that our result is valid on general Riemannian three manifolds allows for the incorporation of spatial heterogeneity such as anisotropic diffusion. In contrast to [30], the fact that we are able to obtain explicitly an expression for the $O(1)$ term in (1.1) is due to the fact that in Proposition 1.1 we have the expansion of $G_{\partial M}(x, z)$ all the way to a remainder $R(x, y)$, which is Hölder continuous at the diagonal. We also appeal to some recent advances in integral geometry [29, 20, 23, 22, 12] to address the comment in [1] on the difficulty of treating this problem in higher dimensions.

The strategy and organization of this paper will be as follows. In Section 2 we will give a brief overview of pseudodifferential operators and their associated Schwartz kernels. The machinery of pseudodifferential operators serve as a bridge between the geometric and analytic objects appearing in Proposition 1.1 and we will compute their coordinate expression. In Section 3 we will use the tools we developed in Section 2 to prove Proposition 1.1. A singularity expansion for the Green's function such as Proposition 1.1 is the gateway for obtaining the asymptotic expansions of Theorems 1.2 and 1.3. However, there is an additional hurdle of inverting an integral transform as mentioned in [1]. Here we make use of some recent advancements in integral geometry and geometric rigidity [23, 12, 20] to overcome these difficulties. This approach is described in Section 4. Finally, in Section 5 we carry out the asymptotic calculation using the tools we have developed. The appendices characterize the expected first arrival time $\mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0=x]$ as the solution of an elliptic mixed boundary value problem. This is classical in the Euclidean case (see [26]) but we could not find a reference for the general case of a Riemannian manifold with boundary.

## 2. OVERVIEW OF PSEUDODIFFERENTIAL OPERATORS

### 2.1. Basic Definitions.
We give some basic definitions and properties of pseudodifferential operators. For a comprehensive treatment we refer the reader to Chapt 7 of
[34] or the book [32]. Readers who are already familiar with microlocal analysis can skip this section.

As usual, $C^\infty$ denotes the space of smooth functions. We use notation $C_c^\infty$ for compactly suported smooth functions and $D'$ for its dual. By $C^k$, we denote the space of $k$ time continuously differentiable functions. The spce of functions from $C^k$, whose $k$th derivatives are Hölder continuous with exponent $\mu \in (0, 1]$, is denoted by $C^{k, \mu}$

Let $a(x, \xi)$ be a smooth function on $T^*\mathbb{R}^n$ and for all $l \in \mathbb{R}$ we say that $a \in S_1^l(T^*\mathbb{R}^n)$ (or simply $S_1^l$) if for all multi-indices $\alpha, \beta$ there are constants $C_{\alpha, \beta}$ such that

$$ (2.1) \qquad |D_{\xi}^{\alpha} D_{x}^{\beta} a(x, \xi)| \leq C_{\alpha, \beta} \langle \xi \rangle^{l-|\alpha|} $$

where $D_{\xi}^{\alpha} = (-i)^{|\alpha|} \partial_{\xi}^{\alpha}$, $D_{x}^{\beta} = (-i)^{|\beta|} \partial_{x}^{\beta}$, and $\langle \xi \rangle := (1 + |\xi|^2)^{1/2}$. These are the Kohn-Nirenberg symbols. This class of symbols contain the classical symbols, denoted by $S_{cl}^l(T^*\mathbb{R}^n)$, which are defined by those $a(x, \xi) \in S_1^l(T^*\mathbb{R}^n)$ satisfying

$$ (2.2) \qquad a(x, \xi) \sim \sum_{m=0}^{\infty} a_{l-m}(x, \xi), $$

where each $a_{l-m}$ are homogeneous in the sense that $a_{l-m}(x, \tau\xi) = \tau^{l-m}a(x, \xi)$ for all $x \in \mathbb{R}^n$, $\tau > 1$ and $|\xi| > 1$. The expression (2.2) means that for all $N$,

$$ a(x, \xi) - \sum_{m=0}^{N} a_{l-m}(x, \xi) \in S_{1}^{l-N-1}(T^{*}\mathbb{R}^{n}). $$

If $a(x, \xi) \in S_1^l$ we can define an operator $a(x, D) : C_c^\infty(\mathbb{R}^n) \to D'(\mathbb{R}^n)$ by

$$ (2.3) \qquad a(x, D)u := \int_{\mathbb{R}^n} e^{i\xi \cdot x} a(x, \xi) \hat{u}(\xi) d\xi, $$

where $\hat{u}(\xi) := \mathcal{F}u := (2\pi)^n \int_{\mathbb{R}^n} e^{-ix \cdot \xi} u(x) dx$ is the Fourier transform. Recall that the absolutely convergent integral representation of the Fourier transform is well defined as an automorphism of the Schwartz class functions $S(\mathbb{R}^n)$ but extends to an automorphism of the tempered distributions $S'(\mathbb{R}^n)$. (See [4] for a comprehensive guide to distribution theory and definition of these spaces).

Operators taking $C_c^\infty(\mathbb{R}^n) \to D'(\mathbb{R}^n)$ which have the above representation are said to be in $\Psi_1^l(\mathbb{R}^n)$ and are called pseudodifferential operators. For the symbol class $S_1^l(T^*\mathbb{R}^n)$, Lemma 1.1 in Chapter 7.1 of [34] extends $a(x, D)$ to map $S'(\mathbb{R}^n) \to S'(\mathbb{R}^n)$.

The classical pseudodifferential operators $\Psi_{cl}^l(\mathbb{R}^n)$ are defined analogously by requiring that $a(x, \xi)$ belongs to $S_{cl}^l$. Note that knowing the operator $a(x, D) \in \Psi_1^l(\mathbb{R}^n)$ we can recover $a(x, \xi) \in S_1^l$ by the formula

$$ (2.4) \qquad a(x, \xi) = e_{-\xi}(x) a(x, D) e_{\xi}, $$

where $e_\xi(x) := e^{i\xi \cdot x}$. Note that if $A(x, y)$ is the Schwartz kernel of the operator $a(x, D)$ then

$$ (2.5) \qquad a(0, \xi) = \mathcal{F}_y^{-1}(A(0, y))(\xi) = \int_{\mathbb{R}^n} e^{i\xi \cdot y} A(0, y) dy. $$
Let $X$ be a compact manifold without boundary. An operator $\mathcal{A} : C^\infty(X) \to \mathcal{D}'(X)^1$ is said to be in $\Psi_1^l(X)$ if there exists coordinate covers $\{(O_j, \Phi_j) \mid \Phi_j : O_j \to U_j \subset \mathbb{R}^n\}$ and a partition of unity $\{\chi_j\}$ subordinate to $\{O_j\}$ such that the map

$$ (2.6) \qquad u \mapsto (\chi_k \mathcal{A} \chi_j \Phi_j^* u) \circ \Phi_k^{-1} $$

from $C^\infty(U_j) \to \mathcal{E}'(U_k)$ belongs to $\Psi_1^l(\mathbb{R}^n)$.

If $a \in C^\infty(T^*X)$ we say that it belongs to the symbol class $S_1^l(T^*X)$ if

$$ \chi_j \circ \Phi_j^{-1} a(\Phi_j^{-1}(\cdot), \Phi_j^*(\cdot)) \in S_1^l(T^*\mathbb{R}^n) $$

for all $j$. The classical pseudodifferential operators $\Psi_{cl}^l(X)$ and classical symbols $S_{cl}^l(T^*X)$ are defined analogously. These definitions depend a-priori on the choice of coordinate systems but turn out to be invariant (see Chapt 7 [34]).

There exists a linear isomorphism

$$ (2.7) \qquad \sigma_l : \Psi_{cl}^l(X)/\Psi_{cl}^{l-1}(X) \to S_{cl}^l(T^*X)/S_{cl}^{l-1}(T^*X) $$

called the *principal symbol map*. For each $\mathcal{A} \in \Psi_{cl}^l(X)$ it can be defined at each $x \in X$ by taking a coordinate neighborhood $O$ containing $x$ and a $\chi \in C_c^\infty(O)$ which is identically 1 near $x$ then considering the operator given in (2.6) for $\chi_j = \chi_k = \chi$ and $\Phi_j = \Phi_k = \Phi$. As the resulting operator in (2.6) is in $\Psi_{cl}^l(\mathbb{R}^n)$ with symbol $a \in S_{cl}^l(T^*\mathbb{R}^n)$, we may set

$$ \sigma_l(\mathcal{A})(x, \Phi^*\xi) := a_l(\Phi(x), \xi) $$

for all $x \in X$ and $\xi \in T_{\Phi_j(x)}^*\mathbb{R}^n$. This definition depends a-priori on the choice of coordinate systems but turns out to be invariant (see Chapt 5 of [32]). In practice these computations are often done in normal coordinates centered at the point of interest $x \in X$ then computing the inverse Fourier transform as in (2.5).

One important property of $\sigma_l$ which we will use is that it respects the product structure of $\Psi_{cl}^l$ and $S_{cl}^l$:

$$ (2.8) \qquad \sigma_l(\mathcal{A})\sigma_m(\mathcal{B}) = \sigma_{l+m}(\mathcal{A}\mathcal{B}) $$

for $\mathcal{A} \in \Psi_{cl}^l(X)$ and $\mathcal{B} \in \Psi_{cl}^m(X)$.

## 2.2. Coordinate Calculations.
We make some calculations for some geometric objects which will naturally appear in the singularity expansion for $G_{\partial M}$. These identities will be useful in proving Proposition 1.1.

Let $(M, g, \partial M)$ be a three dimensional Riemannian manifold with non-empty smooth boundary which inherits the metric $h := \iota_{\partial M}^*g$. Denote by $\Pi$ the scalar second fundamental form on the surface $\partial M$ and $H(x)$ be the mean curvature at $x \in \partial M$. Let $S\partial M$ denote the unit-sphere bundle over $\partial M$,

$$ S\partial M = \{v \in T\partial M \mid \|v\|_h = 1\}. $$

For any $x \in \partial M$, let $E_1(x), E_2(x) \in S_x \partial M$ be principal directions (i.e. unit eigenvectors) of the induced shape operator with eigenvalues $\lambda_1(x)$ and $\lambda_2(x)$. We will drop the dependence in $x$ from our notation when there is no ambiguity.

We choose $E_1$ and $E_2$ such that $E_1^\flat \wedge E_2^\flat \wedge \nu^\flat$ is a positive multiple of the volume form $dvol_g$ (see p.26 of [16] for the “musical isomorphism” notation of $^\flat$ and $^\flat$). Here

<footnote n="1">¹We shall avoid a discussion of distributional sections of density bundles by noting that in our setting of $X$ is always prescribed with a Riemannian volume form which provides a natural trivialization of density bundles. As such all distributional sections of density bundles are identified with sections of the trivial line bundle in this way.</footnote>
we use $\nu$ to denote the outward pointing normal vector field so that it is consistent with most PDE literature. However, in defining II and the shape operator we will follow geometry literature (e.g.[16]) and use the inward pointing normal so that the sphere embedded in $\mathbb{R}^3$ would have positive mean curvature in our convention.

For two points $x, y \in \partial M$ there are two distances to consider. The first is the shortest path amongst those that stay on the boundary which we denote by $d_h(x, y)$ and the other is the distance measured by paths allowing to enter $M$, which we denote by $d_g(x, y)$. Clearly, $d_h(x, y) \ge d_g(x, y)$.

For a fixed $x_0 \in \partial M$, we will denote by $B_h(\rho; x_0) \subset \partial M$ the geodesic disk of radius $\rho > 0$ (with respect to the metric $h$) centered at $x_0$ and $D_\rho$ to be the Euclidean disk in $\mathbb{R}^2$ of radius $\rho$ centered at the origin. In what follows $\rho$ will always be smaller than the injectivity radius of $(\partial M, h)$. Letting $t = (t_1, t_2, t_3) \in \mathbb{R}^3$, we will construct a coordinate system $x(t; x_0)$ by the following procedure:

Write $t \in \mathbb{R}^3$ near the origin as $t = (t', t_3)$ for $t' = (t_1, t_2) \in D_\rho$. Define first 2

$$x((t', 0); x_0) := \exp_{x_0;h}(t_1 E_1 + t_2 E_2),$$

where $\exp_{x_0;h}(V)$ denotes the time 1 map of $h$-geodesics with initial point $x_0$ and initial velocity $V \in T_{x_0}\partial M$. The coordinate $t' \in D_\rho \mapsto x((t', 0); x_0)$ is then an $h$-geodesic coordinate system for a neighborhood of $x_0$ on the boundary surface $\partial M$. We can extend this to become a coordinate system for points in $M$ near $x_0$ so that $t \mapsto x(t; x_0)$ is a boundary normal coordinate system with $t_3 > 0$ in $M$ as the boundary defining function. Readers wishing to know more about boundary normal coordinates can refer to [17] for a brief recollection of the basic properties we use here and Prop 5.26 of [16] for a detailed construction.

For convenience we will write $x(t'; x_0)$ in place of $x((t', 0); x_0)$. The boundary coordinate system $t \mapsto x(t; x_0)$ has the advantage that the metric tensor $g$ can be expressed as

$$(2.9) \qquad \sum_{j,k=1}^{3} g_{j,k}(t) dt_j dt_k = \sum_{\alpha,\beta=1}^{2} h_{\alpha,\beta}(t', t_3) dt_\alpha dt_\beta + dt_3^2,$$

where $h_{\alpha,\beta}(t', 0) = h_{\alpha,\beta}(t')$ is the expression of the boundary metric $h$ in the $h$-geodesic coordinate system $x(t'; x_0)$. Note that $(g_{j,k}(t))_{j,k=1}^3$ and $(h_{\alpha,\beta}(t', t_3))_{\alpha,\beta=1}^2$ are symmetric positive definite $3 \times 3$ and $2 \times 2$ matrices varying smoothly with respect to the variable $t = (t', t_3) = (t_1, t_2, t_3)$.

For $\epsilon > 0$ sufficiently small we define the (rescaled) $h$-geodesic coordinate by the following map

$$x^{\epsilon}(\cdot; x_0) : t' = (t_1, t_2) \in \mathbb{D} \mapsto x(\epsilon t'; x_0) \in B_h(\epsilon; x_0),$$

where $\mathbb{D}$ is the unit disk in $\mathbb{R}^2$. We derive some coordinate expressions for some of the geometric objects we will consider later.

**Lemma 2.1.** Let $h(s') = \sum_{\alpha,\beta=1}^2 h_{\alpha,\beta}(s') ds_\alpha ds_\beta$ be the pullback metric of $h$ by $s' \mapsto x(s'; x_0)$ on $D_\rho$. Denote by

$$r = |s' - t'|_{h(s')} := \left( \sum_{\alpha,\beta=1}^{2} h_{\alpha,\beta}(s')(s_{\alpha} - t_{\alpha})(s_{\beta} - t_{\beta}) \right)^{1/2}$$

<footnote n="2">²for example it is obvious below that $E_1$ and $E_2$ are elements of the tangent space over $x_0$ as they are inserted into the argument of $\exp_{x_0;h}(\cdot)$.</footnote>
so that $t' = s' + r\omega$ where $\omega \in S_{s'}\mathbb{D}_\rho$. We have that

$$d_h(x(s'; x_0), x(t'; x_0))^2 = r^2 \sum_{j,k=1}^{2} H_{j,k}(s', r, \omega) \omega_j \omega_k$$

for matrix $H_{j,k}(s', r, \omega)$ jointly smooth in $(s', r, \omega)$. It also satisfies $H_{j,k}(s', 0, \omega) = h_{j,k}(s')$ and

$$\sum_{j,k} \partial_r H_{j,k}(s', 0, \omega) \omega_k \omega_j = O(s').$$

*Proof.* Expressing $t'$ using $(s', r, \omega)$ we have that (see e.g. [35] Lemma 4.8)

$$d_h(x(s'; x_0), x(t'; x_0)) = r \left( \sum_{\alpha, \beta} H_{\alpha, \beta}(s', r, \omega) \omega_{\alpha} \omega_{\beta} \right)^{1/2}$$

with $H_{\alpha,\beta}$ symmetric, even under the map $(r, \omega) \mapsto (-r, -\omega)$, and $H_{\alpha,\beta}(s', 0, \omega) = h_{\alpha,\beta}(s')$. Setting $s' = 0$ and using the fact that we are using normal coordinates we obtain $r^2 = r^2 (\sum_{\alpha,\beta} H_{\alpha,\beta}(0, r, \omega) \omega_{\alpha} \omega_{\beta})$. Now Taylor expanding $H_{\alpha,\beta}(0, r, \omega)$ in $r$, we see that

$$r^2 = r^2 \left( 1 + \sum_{\alpha, \beta} \partial_r H_{\alpha, \beta}(0, 0, \omega) r \omega_{\alpha} \omega_{\beta} + O(r^2) \right).$$

The $r^2$ terms on left-hand and right-hand sides cancel, leaving

$$0 = \partial_r H_{\alpha,\beta}(0, 0, \omega) r^3 \omega_{\alpha} \omega_{\beta} + O(r^4).$$

Divide through by $r^3$ and take the limit as $r \to 0$ we see that $\sum_{\alpha,\beta} \partial_r H_{\alpha,\beta}(0, 0, \omega) \omega_{\alpha} \omega_{\beta} = 0$ for any $\omega \in S_0\mathbb{D}$. $\square$

**Lemma 2.2.** We use the same notation for $\omega$ and $r$ as in Lemma 2.1. One has that

$$d_h(x(s'; x_0), x(t'; x_0))^{-1} = r^{-1} + O(s') + O(r),$$

where $O(r)$ (respectively $O(s')$) denotes smooth functions of $(s', r, \omega) \in \mathbb{D}_\rho \times \mathbb{R} \times S_{s'}\mathbb{D}_\rho$ which vanish to first order as $r \to 0$ (respectively $s' \to 0$).

*Proof.* From Lemma 2.1 we have that

$$d_h(x(s'; x_0), x(t'; x_0))^{-1} = r^{-1} \left( \sum_{\alpha, \beta} h_{\alpha, \beta}(s') \omega_{\alpha} \omega_{\beta} + r \partial_r H_{\alpha, \beta}(s', 0, \omega) \omega_{\alpha} \omega_{\beta} + O(r^2) \right)^{-1/2}.$$

Using the fact that $\omega \in S_{s'}\mathbb{D}$ with respect to the metric given by $h_{\alpha,\beta}$ we have that

$$d_h(x(s'; x_0), x(t'; x_0))^{-1} = r^{-1} \left( \sum_{\alpha, \beta} 1 + r \partial_r H_{\alpha, \beta}(s', 0, \omega) \omega_{\alpha} \omega_{\beta} + O(r^2) \right)^{-1/2}.$$

For $r$ and $s'$ sufficiently small we may use Taylor's expansion to obtain the desired property. The fact that $\sum_{\alpha,\beta} \partial_r H_{\alpha,\beta}(s', 0, \omega) \omega_{\alpha} \omega_{\beta} = O(s')$ is stated in Lemma 2.1 $\square$

**Corollary 2.3.** For $\epsilon > 0$ sufficiently small we have that

$$d_h(x^\epsilon(s'; x_0), x^\epsilon(t'; x_0))^{-1} = \epsilon^{-1} r^{-1} + \epsilon r^{-1} A(\epsilon, s', r, \omega)$$

for some smooth function $A$ in the variables $(\epsilon, s', r, \omega) \in [0, \epsilon_0] \times \mathbb{D} \times \mathbb{R} \times S^1$. Here we use $r = |s' - t'|$ and $t' = s + r\omega$.
**Lemma 2.4.** Let $x(\cdot; x_0)$ be the coordinate system at the beginning of this section centered at $x_0$. For $s', t' \in \mathbb{R}^2$ sufficiently small, we have that

$$d_g(x(s'; x_0), x(t'; x_0))^2 = r^2 (1 + r \tilde{G}(s', \omega) + O(r^2)),$$

where $r = |s' - t'|_{h(s')}$ and $t' = s' + r\omega$. Here $\tilde{G}(s', \omega)$ is a smooth function of $(s', \omega)$ which vanishes at $s' = 0$. The $O(r^2)$ term is a smooth function in $(s', r, \omega) \in \mathbb{D}_\rho \times \mathbb{R} \times S_{s'}\mathbb{D}_\rho$ which vanishes to second order as $r \to 0$.

*Proof.* We begin with the identity that for any $s$ and $t$,

$$d_g(x(s; x_0), x(t; x_0))^2 = \sum_{j,k=1}^{3} \mathbf{G}_{j,k}(s, t)(s_j - t_j)(s_k - t_k),$$

with $\mathbf{G}_{j,k}(s, s) = g_{j,k}(s)$ given by (2.9). Now set $s = (s', 0)$ and $t = (t', 0)$ we have

$$d_g(x(s'; x_0), x(t'; x_0))^2 = r^2 \left( 1 + \sum_{j,k=1}^{2} \partial_r \mathbf{G}_{j,k}^0(s', 0, \omega) r \omega_j \omega_k + O(r^2) \right)$$

where $\mathbf{G}_{j,k}^0(s', r, \omega) := \mathbf{G}_{j,k}(s', s' + r\omega)$, is even in $(r, \omega) \mapsto -(r, \omega)$ and $O(r^2)$ is a smooth function of $(s', r, \omega)$ which is even and vanishes to second order as $r \to 0$. Observe that $\omega \mapsto \partial_r \mathbf{G}_{j,k}^0(s', 0, \omega)$ is odd in $\omega$.

We now need to argue that $\sum_{\alpha, \beta} \partial_r \mathbf{G}_{j,k}^0(0, 0, \omega) \omega_\alpha \omega_\beta = 0$. Setting $s' = 0$ in the above identity and using the fact that we are using boundary normal coordinates we have

$$r^2 = d_h(x(0; x_0), x(t'; x_0))^2 \geq d_g(x(0; x_0), x(t'; x_0))^2 = r^2 \left( 1 + \sum_{j,k=1}^{2} \partial_r \mathbf{G}_{j,k}^0(0, 0, \omega) r \omega_j \omega_k + O(r^2) \right).$$

Subtracting off the $r^2$ terms and dividing by $r^3$ we see that as $r \to 0$,

$$\sum_{j,k=1}^{2} \partial_r \mathbf{G}_{j,k}^0(0, 0, \omega) \omega_j \omega_k \leq 0$$

for all $\omega \in S^1$. We now use the fact that $\partial_r \mathbf{G}_{j,k}(0, 0, \omega)$ is odd to see that

$$\sum_{j,k=1}^{2} \partial_r \mathbf{G}_{j,k}^0(0, 0, \omega) \omega_j \omega_k = 0.$$

□

Just as how Lemma 2.2 and Corollary 2.3 followed from Lemma 2.1, we have the following:

**Corollary 2.5.** Using the expression $r = |s' - t'|_{h(s')}$ and $t' = s' + r\omega$, we have that

$$d_g(x(s'; x_0), x(t'; x_0))^{-1} = r^{-1} + O(s') + O(r),$$

where $O(r)$ (respectively $O(s')$) denotes smooth function of $(s', r, \omega) \in \mathbb{D}_\rho \times \mathbb{R} \times S_{s'}\mathbb{D}_\rho$ which vanishes to first order as $r \to 0$ (respectively $s' \to 0$).

**Corollary 2.6.** For $\epsilon > 0$ sufficiently small we have that

$$d_g(x^\epsilon(s'; x_0), x^\epsilon(t'; x_0))^{-1} = \epsilon^{-1} r^{-1} + \epsilon r^{-1} A(\epsilon, s', r, \omega)$$

for some smooth function $A$ in the variables $(\epsilon, s', r, \omega) \in [0, \epsilon_0] \times \mathbb{D} \times \mathbb{R} \times S^1$. Here we use $r = |s' - t'|$ and $t' = s + r\omega$.
**Lemma 2.7.** In the coordinate given by $y = x(s'; x_0)$,

$$ \partial_{\nu_y} d_g(x_0, y) = \frac{\lambda_1(x_0)s_1^2 + \lambda_2(x_0)s_2^2}{2|s'|} + O(|s'|^2) $$

where $O(|s'|^2)$ denotes a smooth function of the variables $(|s'|, \frac{s'}{|s'|})$ which vanishes to order 2 at the origin.

*Proof.* By Gauss Lemma, for $y$ near $x_0$, $\text{grad}_y d_g(x_0, y) = \dot{\gamma}(d_g(x_0, y))$ where $\gamma(\cdot)$ is the unit speed geodesic in $(M, g)$ from $x_0$ to $y$. Therefore we have that

$$ (2.10) \qquad \partial_{\nu_y} d_g(x_0, y) = - \left\langle \nu_y, \frac{\exp_{y;g}^{-1}(x_0)}{|\exp_{y;g}^{-1}(x_0)|_g} \right\rangle_g. $$

In the coordinates given by $s \mapsto x(s; x_0)$ the Christoffel symbols are

$$ (2.11) \qquad \Gamma^3_{3,3} = \Gamma^3_{\alpha,3} = \Gamma^3_{3,\alpha} = 0, \quad \Gamma^3_{\alpha,\beta} = -\frac{1}{2} \partial_3 h_{\alpha,\beta}. $$

Choose $(\hat{s}', 0) \in \mathbb{R}^3$ so that $x(\hat{s}', 0; x_0) = y$. By Lemma 2.4, $d_g(y, x_0)$ is a smooth function of $(|\hat{s}'|, \frac{\hat{s}'}{|\hat{s}'|})$ when we write $y = x(\hat{s}'; x_0)$ in these coordinates. Let $V(y) := \frac{\exp_{y;g}^{-1}(x_0)}{|\exp_{y;g}^{-1}(x_0)|_g} = \sum_{j=1}^3 V_j(\hat{s}') \partial_j$ be the unique unit vector over $y$ so that the $(M, g)$ unit velocity geodesic in these coordinates starting at $y$ with initial direction $V(y)$ reaches $x_0$ in time $d_g(y, x_0)$. In the coordinates given by $s' \mapsto x(s'; x_0)$, we want to argue that $V_j(\hat{s}')$ are smooth functions of $(|\hat{s}'|, \frac{\hat{s}'}{|\hat{s}'|})$. To do so, observe that in $g$-geodesic coordinates centered around $x_0$ this is of course the case. The result for any other coordinate system can then be obtained via a change of variable.

Note that the outward pointing normal is given by $-\partial_3$ since $s_3 > 0$ in $M$. Using (2.10) and the expression for the metric (2.9) we have that in the chosen boundary normal coordinate system

$$ (2.12) \qquad \partial_{\nu_y} d_g(x_0, y) = V_3. $$

After time $\tau$, the geodesic with initial position $(\hat{s}', 0)$ and initial unit velocity $V$ can be written in the $s = (s', s_3)$ coordinate as

$$ s(\tau) = (\hat{s}', 0) + \tau V + r(\tau; V) $$

for some remainder $r = (r_1, r_2, r_3)$ which has initial condition $r(0) = \dot{r}(0) = 0$. Taylor expanding $r(\cdot; V)$ we have that

$$ (2.13) \qquad s(\tau) = (\hat{s}', 0) + \tau V + \frac{\tau^2}{2} \ddot{r}(0; V) + \tau^3 r'(\tau; V) $$

for some $r'(\tau; V)$ depending smoothly on $\tau$ and $V$.

Due to (2.12), we are particularly interested in the evolution of the $r_3$ component which solves the ODE

$$ (2.14) \qquad \begin{aligned} \ddot{r}_3(\tau; V) &= -\sum_{j,k=1}^3 \Gamma^3_{j,k}(s(\tau))(V_j + \dot{r}_j)(V_k + \dot{r}_k) \\ &= \frac{1}{2} \sum_{\alpha,\beta=1}^2 \partial_3 h_{\alpha,\beta}(s(\tau))(V_\alpha + \dot{r}_\alpha)(V_\beta + \dot{r}_\beta). \end{aligned} $$

The last equality comes from (2.11).
Now set $\tau = d_g(y, x_0)$ so that $s(\tau) = 0$, we have from (2.13) that for $\alpha = 1, 2$,

$$V_{\alpha} = - \frac{\hat{s}_{\alpha}}{d_g(y, x_0)} + d_g(y, x_0) f_{\alpha}$$

for functions $f_{\alpha}$ which is smooth in $V$ and $d_g(y, x_0)$ and thus smooth functions of the $(\hat{s}', \hat{s}'/|\hat{s}'|)$ variable. Inserting this into (2.14) yields

$$\frac{d_g(x, y)^2}{2} \ddot{r}_3(0; V) = \frac{1}{4} \sum_{\alpha, \beta} \partial_3 h_{\alpha, \beta}(0) \hat{s}_{\alpha} \hat{s}_{\beta} + d_g(y, x)^3 f$$

for some function $f$ which is smooth in the variable $(|\hat{s}'|, \hat{s}'/|\hat{s}'|)$. Inserting this expression into (2.13) for $\tau = d_g(x_0, y)$ we have that

$$\begin{aligned} V_3(\hat{s}') &= -\frac{d_g(x, y)}{2} \ddot{r}_3(0; V) + d_g(x_0, y)^2 r_3'(d_g(x_0, y), V(s')) \\ &= -\frac{1}{4} \sum_{\alpha, \beta} \frac{\partial_3 h_{\alpha, \beta}(0)}{d_g(x_0, y)} \hat{s}_{\alpha} \hat{s}_{\beta} + d_g(y, x_0)^2 (f + r') \end{aligned}$$

Using Lemma 2.4 we have that $d_g(x_0, y)^{-1} = |\hat{s}'|^{-1} + F(\frac{\hat{s}'}{|\hat{s}'|}, s')$ for some smooth function $F(\cdot, \cdot)$. We may thus write

$$V_3(\hat{s}') = -\frac{1}{4} \sum_{\alpha, \beta} \frac{\partial_3 h_{\alpha, \beta}(0)}{|\hat{s}'|} \hat{s}_{\alpha} \hat{s}_{\beta} + O(|\hat{s}'|^2),$$

where $O(|\hat{s}'|^2)$ is a smooth function of $(\frac{\hat{s}'}{|\hat{s}'|}, \hat{s}')$ which vanishes to order 2 near the origin. Now use the fact that $-\frac{1}{2}\partial_3 h_{\alpha, \beta}(0)$ is the coordinate expression for the scalar second fundamental form with respect to the normal given by $\partial_3$ 3 at the point $x_0$ (see Proposition 8.17 of [16]) and our coordinate system $x(s'; x_0)$ is chosen so that the shape operator is diagonalized at $x_0$. Therefore,

$$V_3(\hat{s}') = \frac{1}{2} \sum_{\alpha, \beta} \frac{\lambda_1 \hat{s}_1^2 + \lambda_2 \hat{s}_2^2}{|\hat{s}'|} + O(|\hat{s}'|^2).$$

In view of (2.12) we have proven the required identity.

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

Let $x_0 \in \partial M$ and define $R_{II}(\cdot, \cdot), R_{II*}(\cdot, \cdot) \in L^{\infty}(B_h(\rho; x_0)B_h(\rho; x_0))$ by

$$\begin{aligned} (2.15) \qquad R_{II}(x, y) &:= \Pi_x \left( \frac{\exp_{x;h}^{-1} y}{|\exp_{x;h}^{-1} y|_h}, \frac{\exp_{x;h}^{-1} y}{|\exp_{x;h}^{-1} y|_h} \right), \\ R_{II*}(x, y) &:= \Pi_x \left( * \frac{\exp_{x;h}^{-1} y}{|\exp_{x;h}^{-1} y|_h}, * \frac{\exp_{x;h}^{-1} y}{|\exp_{x;h}^{-1} y|_h} \right). \end{aligned}$$

Here $*$ is the Hodge star operator associated to the metric $h$.

**Lemma 2.8.** Let $x(\cdot; x_0) : D_{\rho} \to \partial M$ be a normal coordinate system for $(\partial M, h)$ centred around $x_0 \in \partial M$ and let $h$ denote the pullback metric tensor on $D_{\rho}$ under this coordinate. Then for all $s', t' \in D_{\rho}$ sufficiently close to the origin,

$$\frac{\exp_{t';h}^{-1}(s')}{|\exp_{t';h}^{-1}(s')|_{h(t')}} = \sum_{j=1}^{2} \omega_j \partial_j + O(t') + O(r),$$

<footnote n="3">³recall that we are using the convention where $II$ and shape operator are defined with respect to the inward pointing normal</footnote>
where $\omega := \frac{s'-t'}{r} \in S^1$ with $r := |s' - t'|$ being the Euclidean distance between $s'$ and $t'$. The $O(t')$ (resp $O(r)$) term denotes a smooth map of $(t', \omega, r) \in \mathbb{D}_{\rho} \times S^2 \times [0, \rho]$ which vanishes to order 1 as $t' \to 0$ (resp $r \to 0$).

*Proof.* This comes from the fact that for some matrix $H_{j,k}(s', t')$ smooth in $(s', t')$

$$d_h(s', t')^2 = \sum_{j=1}^{2} H_{j,k}(s', t')(s_j - t_j)(s_k - t_k),$$

where $H_{j,k}(t', t') = h_{j,k}(t') = \delta_{j,k} + O(|t'|^2)$. Therefore

$$d_{h(s',t')}(s', t') = |s' - t'| + |s' - t'|F\left(t', \frac{s' - t'}{|s' - t'|}, |s' - t'|\right)$$

for some smooth function $F \in C^{\infty}(\mathbb{D}_{\rho} \times S^1 \times [0, r_0])$ which is $O(t') + O(r)$.

A coordinate calculation yields that $\frac{\exp_{t',h}^{-1}(s')}{|\exp_{t',h}^{-1}(s')|_{h(t')}} = \text{grad}_{t'} d_h(s', t')$ so

$$\frac{\exp_{t',h}^{-1}(s')}{|\exp_{t',h}^{-1}(s')|_{h(t')}} = \sum_{j,k=1}^{2} h_{j,k}(t') \partial_{t_j} d_{h^{\theta}(s',t')}(s', t') \partial_k = \sum_j \omega_j \partial_j + O(t') + O(r).$$

□

**Corollary 2.9.** Let $x_0 \in \partial M$ and $\lambda_1(x_0)$ and $\lambda_2(x_0)$ be the eigenvalues of the shape operator at $x_0$. Then,

i) in the $t' \mapsto x(t'; x_0)$ coordinate system prescribed at the beginning of this section,

$$R_{II}(x(t'; x_0), x(s'; x_0)) - \left( \lambda_1(x_0) \frac{(s_1 - t_1)^2}{|s' - t'|^2} + \lambda_2(x_0) \frac{(s_2 - t_2)^2}{|s' - t'|^2} \right) = O(r) + O(t')$$

$$R_{II*}(x(t'; x_0), x(s'; x_0)) - \left( \lambda_2(x_0) \frac{(s_1 - t_1)^2}{|s' - t'|^2} + \lambda_1(x_0) \frac{(s_2 - t_2)^2}{|s' - t'|^2} \right) = O(r) + O(t')$$

The $O(t')$ (resp $O(r)$) term denotes a smooth function of $(t', \omega, r) \in \mathbb{D}_{\rho} \times S^2 \times [0, \rho]$ which vanishes to order 1 as $t' \to 0$ (resp $r \to 0$).

ii) In the $t' \mapsto x^{\epsilon}(t'; x_0)$ coordinate systems prescribed at the beginning of this section,

$$R_{II}(x^{\epsilon}(t'; x_0), x^{\epsilon}(s'; x_0)) - \left( \lambda_1(x_0) \frac{(s_1 - t_1)^2}{|s' - t'|^2} + \lambda_2(x_0) \frac{(s_2 - t_2)^2}{|s' - t'|^2} \right) = \epsilon R_{\epsilon}(t', \omega, r),$$

$$R_{II*}(x^{\epsilon}(t'; x_0), x^{\epsilon}(s'; x_0)) - \left( \lambda_2(x_0) \frac{(s_1 - t_1)^2}{|s' - t'|^2} + \lambda_1(x_0) \frac{(s_2 - t_2)^2}{|s' - t'|^2} \right) = \epsilon R_{\epsilon}(t', \omega, r),$$

where $R_{\epsilon}(t', \omega, r)$ is smooth with derivatives of all orders uniformly bounded in $\epsilon$.

*Proof.* Since $x^{\epsilon}(t'; x_0) = x(\epsilon t'; x_0)$ we have that ii) is a consequence of i). For i), we will only prove this for $R_{II}$ since the statement for $R_{II*}$ can be obtained via a rotation.

Recall that the normal coordinate system $t' \mapsto x(t'; x_0)$ is chosen so that at $x_0$ the coordinate vectors $\{\partial_{t_1}, \partial_{t_2}\}$ pushes forward under $x(\cdot; x_0)$ to become of eigenvectors of the shape operator. Because of this we have that the pull-back of $\Pi_x(\cdot, \cdot)$ under this coordinate system is given by $\Pi_{x(t'; x_0)} = \sum_{j,k=1}^2 \Pi_{j,k} dt_j dt_k$ where $\Pi_{j,k} = \lambda_j(x_0) \delta_{j,k} + O(t')$. Here $O(t')$ denotes a smooth function of $t'$ which vanishes at the origin. Using the expression derived in Lemma 2.8 for the vector $\frac{\exp_{x;h}^{-1} y}{|\exp_{x;h}^{-1} y|_h}$ in the
coordinate given by $x = x(t'; x_0)$ and $y = x(s'; x_0)$, we have the desired expression for $R_{II}(x(t'; x_0), x(s'; x_0))$. $\square$

**2.3. Operator Estimates.** In this section we derive Sobolev estimates for some integral kernels we will encounter when obtaining the asymptotic expansions of Theorems 1.2 and 1.3. As these depend on a parameter $\epsilon > 0$ and do not immediately fit into the framework of semiclassical calculus, we need to keep track of the bounds by hand.

It is useful to take the Fourier transform with respect to only some variables. Let $u(s', t')$ be a family of tempered distributions in $t' \in \mathbb{R}^2$ depending smoothly on the parameter $s' \in \mathbb{R}^2$. That is, it is the distribution defined by $\phi \mapsto \int_{\mathbb{R}^2} u(s', t')\phi(t')dt'$ for all $\phi \in S(\mathbb{R}^2)$. We denote by $\mathcal{F}_{t'}(u(s', t'))(\xi')$ to be the Fourier transform with respect to the $t'$ variable only.

**Lemma 2.10.** Let $A(s', \omega)$ be a smooth function on $(s', \omega) \in \mathbb{R}^2 \times S^1$. For $j \ge 0$ and $\xi' \ne 0$ we have that for any multi-index $\alpha$,

$$D_{s'}^{\alpha} \mathcal{F}_{t'} \left( A \left( s', \frac{t'}{|t'|} \right) |t'|^{j-1} \right) (\xi') = \mathcal{F}_{t'} \left( D_{s'}^{\alpha} A \left( s', \frac{t'}{|t'|} \right) |t'|^{j-1} \right) (\xi').$$

Furthermore, $\mathcal{F}_{t'} \left( A \left( s', \frac{t'}{|t'|} \right) |t'|^{j-1} \right) (\xi')$ is jointly smooth in $(s', \xi') \in \mathbb{R}^2 \times \mathbb{R}^2 \setminus \{0\}$.

*Proof.* Let $\chi \in C_c^\infty(\mathbb{R}^2)$ be identically 1 near the origin. We can write for any positive integer $k$ and $\xi' \ne 0$,

$$\mathcal{F}_{t'} \left( A \left( s', \frac{t'}{|t'|} \right) |t'|^{j-1} \right) (\xi') = \mathcal{F}_{t'} \left( A \left( s', \frac{t'}{|t'|} \right) |t'|^{j-1} \chi(t') \right) (\xi') \\ + |\xi'|^{-2k} \mathcal{F}_{t'} \left( \Delta_{t'}^k \left( A \left( s', \frac{t'}{|t'|} \right) |t'|^{j-1} (1 - \chi(t')) \right) \right) (\xi').$$

The first integral is absolutely convergent by the $\chi(t')$ cut-off. The second Fourier transform is also absolutely convergent owing to the fact that the integrand is smooth and $\Delta_{t'}^k$ makes the integrand decay quickly as $t' \to \infty$ provided $k$ is chosen large enough. $\square$

**Lemma 2.11.** Let $A_\epsilon(s', r, \omega)$ be a family of $C_c^\infty(\mathbb{R}^2 \times \mathbb{R} \times S^1)$ whose support is uniformly bounded in $\epsilon \in [0, \epsilon_0]$ and whose derivatives are also uniformly bounded in $\epsilon$. Then for all $l \ge 0$,

$$A_\epsilon f := \int_{\mathbb{R}^2} A_\epsilon(s', r, \omega) r^l f(s' + r\omega) dr d\omega$$

is a map bounded uniformly in $\epsilon$ from $H^m(\mathbb{R}^2) \to H^{m+1+l}(\mathbb{R}^2)$.

*Proof.* We prove the estimates only for Schwartz functions $f \in S(\mathbb{R}^2)$. We first expand

$$(2.16) \qquad A_{\epsilon} = \sum_{j=1}^{N} \partial_{r}^{j} A_{\epsilon}(s', 0, \omega) r^{j} + r^{N+1} R_{\epsilon}(s', r, \omega).$$

All terms are smooth in its variables with derivatives uniformly bounded in $\epsilon$. Estimating the remainder term in (2.16) is easy:

$$A_{R}f := \int_{S^1} \int_{0}^{\infty} r^{N+1+l} R_{\epsilon}(s', r, \omega) f(s'+r\omega) dr d\omega = \int_{\mathbb{R}^2} |s'-t'|^{N+l} R_{\epsilon} \left( s', |s'-t'|, \frac{s'-t'}{|s'-t'|} \right) f(t') dt'$$
For any positive integer $m$ we write $f = \langle D \rangle^{2m} \langle D \rangle^{-2m} f$. Choose $N \gg m$ so that there is sufficient smoothness in the integral kernel to integrate by parts the formula

$$ \mathcal{A}_R f = \int_{\mathbb{R}^2} |s' - t'|^{N+l} R_\epsilon \left( s', |s' - t'|, \frac{s' - t'}{|s' - t'|} \right) \langle D_{t'} \rangle^{2m} \langle D_{t'} \rangle^{-2m} f(t') dt'. $$

We see from this that for a fixed positive integer $m$ we may choose $N$ large enough so that $\mathcal{A}_R : H^{-2m}(\mathbb{R}^2) \to H^{2m}(\mathbb{R}^2)$ is uniformly bounded in $\epsilon$.

For the integral involving the main term of (2.16), we write

$$ \mathcal{A}_j f := \int_{S^1} \int_0^\infty \partial_r^j A_\epsilon(s', 0, \omega) r^{j+l} f(s' + r\omega). $$

Let $\chi \in C_c^\infty(\mathbb{R}^2)$ be 1 near the origin and write

$$ (2.17) \qquad \mathcal{A}_j f = \mathcal{A}_j \chi(D)f + \mathcal{A}_j(1 - \chi(D))f. $$

To see the mapping property of the first term of (2.17) we write it out in Cartesian coordinates

$$ \mathcal{A}_j \chi(D) f = \int_{\mathbb{R}^2} \partial_r^j A_\epsilon \left( s', 0, \frac{t'}{|t'|} \right) |t'|^{j-1+l} (\tilde{\chi} * f)(s' - t') dt'. $$

Since $\partial_r^j A_\epsilon(s', 0, \frac{t'}{|t'|})$ is smooth in $s'$ with derivatives bounded uniformly in $\epsilon$ and $\chi(D)$ is smoothing, we have that $\mathcal{A}_j \chi(D) : H^{-m}(\mathbb{R}^2) \to H^m(\mathbb{R}^2)$ for any positive integer $m$ with bound uniform in $\epsilon$.

The second term of (2.17) is a pseudodifferential operator with full symbol

$$ a_j(s', \xi'; \epsilon) := (1 - \chi(\xi')) \mathcal{F}_{t'} \left( \partial_r^j A_\epsilon \left( s', 0, \frac{t'}{|t'|} \right) |t'|^{j-1+l} \right) $$

and away from $\xi = 0$ we can deduce from Lemma 2.10 that

$$ \mathcal{F}_{t'} \left( \partial_r^j A_\epsilon \left( s', 0, \frac{t'}{|t'|} \right) |t'|^{j-1+l} \right) = |\xi'|^{-j-l-1} \tilde{a}_j(s', \xi'|\xi'|; \epsilon) $$

for some $\tilde{a}_j(s', \omega; \epsilon) \in C_c^\infty(\mathbb{R}^2 \times S^1)$ with derivatives uniformly bounded in $\epsilon$. The operator $\langle D \rangle^{m+l+1} \mathcal{A}_j(1 - \chi(D)) \langle D \rangle^{-m}$ then has full symbol in $S^0$ with symbol semi-norms bounded uniformly in $\epsilon$. We can now apply Calderón-Vailancourt Theorem to deduce that $\mathcal{A}_j(1 - \chi(D))$ is bounded uniformly in $\epsilon$ from $H^m(\mathbb{R}^2) \to H^{m+l+1}(\mathbb{R}^2)$. $\square$

**2.4. Symbol Computations.** Compute symbol by taking the Fourier transform and multiply by $2\pi$. We compute the principal symbols of some of the main operators which we will encounter. The following list of inverse Fourier transforms will be useful for later computations and we will leave its proof to the reader:

**Lemma 2.12.** In $\mathbb{R}^2$ with $\xi = (\xi_1, \xi_2)$ and $x = (x_1, x_2)$ one has that for $|\xi| \ge 1$,

$$ i) \mathcal{F}^{-1}(\log|x|)(\xi) = -2\pi|\xi|^{-2}, \quad \mathcal{F}^{-1}(|x|^{-1})(\xi) = 2\pi|\xi|^{-1} $$

$$ ii) \mathcal{F}^{-1}(x_j|x|^{-1})(\xi) = 2\pi i \xi_j |\xi|^{-3} $$

$$ iii) \mathcal{F}^{-1}(x_j^2|x|^{-3})(\xi) = 2\pi \xi_k^2 |\xi|^{-3}, \quad k \neq j $$

$$ iv) \mathcal{F}^{-1}(x_j^2|x|^{-2})(\xi) = 2\pi(\xi_k^2 - \xi_j^2)|\xi|^{-4}, \quad k \neq j. $$

**Remark 2.13.** Note that we ignore the behaviour of $\mathcal{F}^{-1}(\cdot)$ near $\xi = 0$ as they are irrelevant to the principal symbol computations we are interested in.

**Lemma 2.14.** Let A be a pseudodifferential operator on $\partial M$ whose singularity along the diagonal is given by $d_g(x, y)^{-1}$. Then $\sigma(A)(x, \xi) \equiv 2\pi|\xi|_{h(x)}^{-1}$.
**Proof.** We compute the symbol at $x_0 \in \partial M$ using the normal coordinate $t' \mapsto x(t'; x_0)$. By Corollary 2.5 $d(x(t'; x_0), x_0)^{-1} = |t'|^{-1} + O(|t'|)$. By Lemma 2.12, the inverse Fourier transform of the leading order singularity is $2\pi|\xi|^{-1}$ for $\xi$ large and therefore $\sigma(A)(x, \xi) \equiv 2\pi|\xi|^{-1}_{h(x)}$. $\square$

**Lemma 2.15.** Let A be a pseudodifferential operator on $\partial M$ whose singularity near the diagonal is given by $\partial_{\nu_y} d_g(x, y)^{-1}|_{x,y \in \partial M \times \partial M}$. Then $\sigma(A)(x_0, \xi) \equiv -\pi \frac{\Pi_{x_0}(\xi^2)}{|\xi|_h^3}$ where $\xi \mapsto \xi^2$ denotes the musical isomorphism from $T^*\partial M$ to $T\partial M$ induced by the boundary metric h and $*$ is the Hodge star operator in this metric. Here we use $\Pi_x(V)$ to denote the quadratic form $\Pi_x(V, V)$ for $V \in T_x\partial M$.

**Proof.** Using Lemma 2.7 we see that in normal coordinates given by $x(s'; x_0)$ the leading order singularity of $\partial_{\nu_y} d_g(x_0, y = x(s'; x_0))$ is given by $-\frac{\lambda_1 s_1^2 + \lambda_2 s_2^2}{2|s'|^3}$. By Lemma 2.12 we have that

$$ \mathcal{F}^{-1} \left( -\frac{\lambda_1 s_1^2 + \lambda_2 s_2^2}{2|s'|^3} \right) = -\pi \frac{\lambda_1 \xi_2^2 + \lambda_2 \xi_1^2}{|\xi|^3}. $$

This is precisely the normal coordinate expression for $-\pi \frac{\Pi_{x_0}(\xi^2)}{|\xi|_h^3}$. $\square$

**Proposition 2.16.** Let $H(x)$ denote the mean curvature of $\partial M$ at $x$, $\Pi_x$ the second fundamental form of $\partial M$ at $x \in \partial M$, and $\Pi_x(V) := \Pi_x(V, V)$ for $V \in T_x\partial M$. Define

$$ (2.18) \quad K(x, y) := H(x)\pi \log d_h(y, x) - \frac{\pi}{4} \left( \Pi_x \left( \frac{\exp_x^{-1}(y)}{|\exp_x^{-1}(y)|_h} \right) - \Pi_x \left( \frac{\ast \exp_x^{-1}(y)}{|\exp_x^{-1}(y)|_h} \right) \right), $$

where $*$ is the Hodge star operator for the metric h. Let A : $C^\infty(\partial M) \to \mathcal{D}'(\partial M)$ be the operator defined by

$$ A : u \mapsto \int_{\partial M} K(x, y)u(y)dvol_h. $$

Then $A \in \Psi_{cl}^{-2}(\partial M)$ with principal symbol $a(x, \xi) \in C^\infty(T^*\partial M \setminus \{0\})$ given by

$$ a(x, \xi) = -\frac{2\pi^2}{|\xi|^4} \Pi_x(\ast\xi^2), $$

where $\xi \mapsto \xi^2 \in T\partial M$ denotes the raising of index with respect to the metric h on $\partial M$.

**Proof.** To see that A is a classical $\Psi$DO, we use Lemma 2.1 and Corollary 2.9 to see that the coordinate expression for the integral kernel $K(x, y)$ satisfies the polyhomo- geneous conditions of Prop 2.8 in Chapter 7 of [34]. Therefore $A \in \Psi_{cl}^{-2}(\partial M)$.

The principal symbol computation is done using normal coordinates. Fix $x_0 \in \partial M$ and denote by

$$ t' \mapsto x(t'; x_0) := \exp_{x_0; h}(t') $$

the normal coordinate system around $x_0$. By a rotation we can choose the coordi- nates so that $\partial_{t_j} x(0; x_0) \in T_{x_0}\partial M$ is an eigenvector of the shape operator at $x_0$ with eigenvalue $\lambda_j$.

According to Lemma 2.1 and Corollary 2.9, in these coordinates the terms of $K(x_0, y)$ can be expressed as

$$ (2.19) \quad K(x_0, x(t'; x_0)) = \frac{\lambda_1 + \lambda_2}{2} \pi \log |t'| - \frac{\pi}{4} \left( \frac{\lambda_1 t_1^2 + \lambda_2 t_2^2}{|t'|^2} - \frac{\lambda_1 t_2^2 + \lambda_2 t_1^2}{|t'|^2} \right) $$
for $t'$ close to the origin. Computing the principal symbol $a(x_0, \xi)$ amounts to taking the inverse Fourier transform of the above expression, and observe the behaviour as $|\xi| \to \infty$. Use the formula in Lemma 2.12, we obtain

$$a(x_0, \xi) = -\frac{2\pi^2}{|\xi|^4}(\lambda_1\xi_2^2 + \lambda_2\xi_1^2) = -\frac{2\pi^2}{|\xi|^4}\Pi_x(*\xi^\sharp).$$

The last equality holds due to the fact that we are using normal coordinates. $\square$

### 3. PROOF OF PROPOSITION 1.1

In this section we use layer potential methods to pick out the singularity structure of the Neumann Green's function at the boundary. Assume without loss of generality that $M$ is an open subset of a compact Riemannian manifold $(\tilde{M}, g)$ without boundary. Choose $M' \subset \tilde{M}$ a manifold with boundary which compactly contains $M$. For all $F \in C_0^\infty(M')$, standard elliptic theory shows that there exists a unique solution $U_F \in C^\infty(M')$ to

$$\Delta_g U_F = F, \quad U_F|_{\partial M'} = 0.$$

The map $F \mapsto U_F$ is a continuous linear operator from $C_0^\infty(M') \to \mathcal{D}'(M')$ and is therefore given by a Schwartz kernel $E(x, y) \in \mathcal{D}'(M' \times M')$ which we call the Green's function. Note that for any $u \in C_0^\infty(M')$, if we fix $x \in M'$ then by definition

$$u(x) = \int_{M'} \Delta_g u(y) E(x, y) d\mathrm{vol}_g(y) = \langle u(\cdot), \Delta_g E(x, \cdot) \rangle.$$

We formally write

$$(3.1) \qquad \Delta_g E(x, \cdot) = \delta_x(\cdot) \quad \text{on } M'$$

Note that if $M$ is a bounded domain in $\mathbb{R}^3$ then $\tilde{M}$ can be chosen to be the flat torus and $E(x, y)$ can be chosen to be $\frac{-1}{4\pi|x-y|}$ for the appropriate constant $c \in \mathbb{R}$.

Using standard elliptic parametrix construction in normal coordinates we express $E(x, y)$ in the following way.

**Lemma 3.1.** For all $x, y \in M'$

$$(3.2) \qquad E(x, y) = -\frac{1}{4\pi}d_g(x, y)^{-1} + \Psi_{cl}^{-4}(\tilde{M}).$$

Here $\Psi_{cl}^{-4}(\tilde{M})$ denotes the Schwartz kernel of an operator in $\Psi_{cl}^{-4}(\tilde{M})$.

*Proof.* Let $P \in \Psi_{cl}^{-2}(\tilde{M})$ be a parametrix for the $\Delta_g$ on the closed compact manifold $\tilde{M}$ without boundary meaning that

$$\Delta_g P = I + \Psi^{-\infty}(\tilde{M}).$$

By ellipticity, for any $\chi_0 \in C_c^\infty(M')$ we have

$$\chi_0(x)\chi_0(y)(E(x, y) - P(x, y)) \in \Psi^{-\infty}(\tilde{M}).$$

Therefore it suffices to show that

$$(3.3) \qquad P - P_{-2} \in \Psi_{cl}^{-4}(\tilde{M})$$

where $P_{-2} \in \Psi^{-2}(\tilde{M})$ is defined by

$$(P_{-2}u)(x) := \int_{\tilde{M}} \frac{-\chi(x, y)}{4\pi} d_g(x, y)^{-1} u(y) d\mathrm{vol}_g(y)$$
for some smooth function $\chi(\cdot, \cdot) \in C^\infty(\tilde{M} \times \tilde{M})$ satisfying $\chi(x, y) = \chi(y, x)$, $\chi(x, y) = 1$ if $d_g(x, y) < \text{Inj}_{\tilde{M}}/4$ and

$$ \text{supp}(\chi(\cdot, \cdot)) \subset \subset \{(x, y) \mid d_g(x, y) < \text{Inj}_{\tilde{M}}/2\}. $$

Here $\text{Inj}_{\tilde{M}}$ is the injectivity radius of the closed compact Riemannian manifold $(\tilde{M}, g)$.

By elliptic regularity (3.3) is equivalent to showing that

$$ (3.4) \qquad \Delta_g P_{-2} - I \in \Psi_{cl}^{-2}(\tilde{M}). $$

Taking the adjoint and use the self-adjointness of both $\Delta_g$ and $P_{-2}$, this is the same as

$$ (3.5) \qquad P_{-2}\Delta_g - I \in \Psi_{cl}^{-2}(\tilde{M}). $$

Using the principal symbol map defined in (2.7) it amounts to showing that

$$ (3.6) \qquad \sigma_{-1}(P_{-2}\Delta_g - I) = 0 $$

as an element of the quotient space $S_{cl}^{-1}/S_{cl}^{-2}$. In fact, since the symbol is classical, we now choose $\sigma_{-1}(P_{-2}\Delta_g - I)$ to be the representative in the equivalence class which is positively homogeneous of degree $-1$.

For each $y_0 \in \tilde{M}$ and covector $\eta_0 \in S_{y_0}^* \tilde{M}$ in the unit cosphere bundle we will show that

$$ (3.7) \qquad |\sigma_{-1}(P_{-2}\Delta_g - I)(y_0, \tau\eta_0)| \le C_{y_0, \eta_0}\tau^{-2} $$

as $\tau \to \infty$. Homogeneity would then ensure that $\sigma_{-1}(P_{-2}\Delta_g - I)(y, \eta) = 0$ for all $(y, \eta) \in T^*\tilde{M}$ which would then ensure (3.6).

To this end let $V_1, V_2, V_3 \in S_{y_0} \tilde{M}$ be three orthonormal vectors and choose normal coordinate given by $\Phi(t) := \exp_{y_0}(\sum_{j=1}^3 t_j V_j)$ for $|t|_{\mathbb{R}^3} < \text{Inj}_{\tilde{M}}/2$. Let $\chi_{\mathbb{R}^3} \in C_c^\infty(\mathbb{R}^3)$ take the value 1 in an open set containing the origin but $\text{supp}(\chi_{\mathbb{R}^3}) \subset \subset \{t \in \mathbb{R}^3 \mid |t|_{\mathbb{R}^3} < \text{Inj}_{\tilde{M}}/2\}$. Similarly let $\chi_{\tilde{M}} \in C^\infty(\mathbb{R}^3)$ take the value 1 in an open set containing $y_0$ but $\text{supp}(\chi_{\tilde{M}}) \subset \subset \{x \mid d_g(x, y_0) < \text{Inj}_{\tilde{M}}/2\}$.

Define the pullback operators $A, B: C^\infty(\mathbb{R}^3) \to C^\infty(\mathbb{R}^3)$ by

$$ A: u \mapsto \Phi_*(\chi_{\tilde{M}} P_{-2} \Phi^*(\chi_{\mathbb{R}^3} u)), \quad B: u \mapsto \Phi_*(\chi_{\tilde{M}} \Delta_g \Phi^*(\chi_{\mathbb{R}^3} u)) $$

where $\Phi^*$ and $\Phi_*$ are pullback by $\Phi$ and $\Phi^{-1}$ respectively.

Thanks to the invariance of the principal symbol map under symplectomorphism, we have

$$ (3.8) \qquad \sigma_{-1}(P_{-2}\Delta_g - I)(y_0, \Phi^*\xi) = \sigma_{-1}(AB - I)(0, \xi) $$

for all $\xi \in T_0^*\mathbb{R}^3$. We see then that (3.6) amounts to showing that $AB - I$ satisfies

$$ (3.9) \qquad \sigma_{-1}(AB - I)(0, \xi) = 0. $$

Let $a(t, \xi)$ and $b(t, \xi)$ be the full symbol of $A$ and $B$ respectively. The full symbol of $A$ can be computed by the formula

$$ a(t, \xi) = e_{-\xi}(t) (A e_{\xi})(t) $$

where $e_\xi(t) := e^{-it\cdot\xi}$. Since we are using the normal coordinate around $y_0$, $d_g(y_0, \Phi(t)) = |t|_{\mathbb{R}^3}$ and

$$ \Phi_* d\text{Vol}_g = \sqrt{|g|} dt = dt + H_2(t) dt $$

where $H_2(t)$ is a smooth function vanishing to order 2 at $t = 0$. So

$$ (A e_{\xi})(0) = \int_{\mathbb{R}^3} \frac{-4\pi}{|t|} e^{-it\cdot\xi} dt + \int_{\mathbb{R}^3} \frac{c}{|t|} e^{-it\cdot\xi} \tilde{H}_2(t) dt $$
for some smooth and compactly supported function $\tilde{H}_2(t)$ vanishing to order 2 at $t=0$. Computing the first term directly and treat the second term by expanding $H_2(t)$ using Taylor expansion we see that

$$ (3.10) \qquad a(0, \xi) = (Ae_{\xi})(0) = |\xi|^{-2} + S_{cl}^{-4}(T^*\mathbb{R}^3). $$

Since $B$ is the Laplace operator in the coordinate given by $\Phi$ we have that

$$ (3.11) \qquad b(t, \xi) = \sum_{j,k=1}^{3} g^{j,k}(t) \xi_j \xi_k + \frac{1}{\sqrt{|g|}} \sum_{j,k=1}^{3} \partial_{t_j} (\sqrt{|g|} g^{j,k})(t) \xi_k $$

Composition calculus gives that if $c(x, \xi)$ is the full symbol of the operator $AB$ then

$$ (3.12) \qquad c(0, \xi) = a(0, \xi)b(0, \xi) + i \sum_{j=1}^{3} \partial_{\xi_j} a(0, \xi) \partial_{t_j} b(t, \xi) \Big|_{t=0} + S_{cl}^{-2}. $$

Substituting into (3.12) the expression we have in (3.10), (3.11), and the fact that in normal coordinates $g^{j,k}(t) = \delta^{j,k} + O(|t|^2)$ for $t$ in a neighbourhood of the origin we have that the second term in (3.12) drops out. So the full symbol of $AB - I$ at the point $(0, \xi) \in T^*\mathbb{R}^3$ is

$$ c(0, \xi) - 1 \in S_{cl}^{-2}. $$

In light of (3.8), for each fixed $\xi \in T_{y_0}^* \tilde{M}$,

$$ |\sigma_{-1}(P_{-2}\Delta_g - I)(y_0, \tau\Phi^*\xi)| < C_\xi \tau^{-2} $$

as $\tau \to \infty$. Therefore (3.7) is verified.

□

For all $f \in C^\infty(\partial M)$ we define as in [34] the operators $S, N \in \Psi_{cl}^{-1}(\partial M)$ by the following

$$ (3.13) \quad Sf(x) := \int_{\partial M} E(x, y)f(y)\mathrm{vol}_h, \quad Nf(x) := 2 \int_{\partial M} \partial_{\nu_y} E(x, y)f(y)\mathrm{vol}_h $$

for $x \in \partial M$. Note that $Nf(x)$ is different from (see [34] Chapt 7 Sect 11)

$$ \lim_{x \to \partial M, x \in M} 2 \int_{\partial M} \partial_{\nu_y} E(x, y) f(y) dy = f(x) + Nf(x). $$

Modulo lower order pseudodifferential operator, $S$ and $N$ are given by the integral kernels $d_g(x, y)^{-1}$ and $\partial_{\nu_y} d_g(x, y)^{-1}$ respectively. Indeed, using (3.2) and equation (11.14) on page 38 of [34], we see that for $(x, y) \in \partial M \times \partial M$ in a neighbourhood of the diagonal,

$$ (3.14) \quad S = -\frac{1}{4\pi} d_g(x, y)^{-1} + \Psi_{cl}^{-3}(\partial M), \quad N = -\frac{1}{2\pi} \partial_{\nu_y} d_g(x, y)^{-1} + \Psi_{cl}^{-2}(\partial M). $$

Using (3.1) we can construct the so called *Neumann Green's function* on $M$ via the following procedure. For each fixed $x \in M^o$ we can solve the following Neumann boundary value problem to obtain the correction term $C(x, y)$ as a function of $y \in M$

$$ \Delta_g C(x, y) = 0, \quad \partial_{\nu_y} C(x, y) = \partial_{\nu_y} E(x, y) - \frac{1}{|\partial M|} \int_{\partial M} C(x, y) d\mathrm{vol}_h = \int_{\partial M} E(x, y) d\mathrm{vol}_h. $$
Setting $G(x, y) = -E(x, y) + C(x, y)$ we get, for each fixed $x \in M$ the unique solution (as a distribution in $z$) $G(x, z)$ to

$$ (3.15) \Delta_g G(x, z) = -\delta_x(z) , \partial_{\nu_z} G(x, z) \mid_{z \in \partial M} = \frac{-1}{|\partial M|} , \int_{\partial M} G(x, z) d\text{vol}_h = 0. $$

Fix for the time being $y \neq z$ in the interior of $M$ and observe that $x \mapsto G(z, x)$ is smooth in a neighbourhood of the singularity of the map $x \mapsto G(x, y)$ and vice versa. Therefore we can integrate by parts the expression $G(z, y) = -\int_M G(z, x) \Delta_x G(y, x) dx$ to obtain

$$ G(z, y) - G(y, z) = \int_{\partial M} G(y, x) \partial_{\nu_x} G(z, x) - G(z, x) \partial_{\nu_x} G(y, x) \text{vol}_h(x). $$

The boundary and orthogonality conditions in (3.15) ensures that the right side vanishes so we have

$$ (3.16) \qquad G(z, y) = G(y, z). $$

Let $\Lambda : H^k(\partial M) \to H^{k-1}(\partial M)$ denote the Dirichlet-to-Neumann map (see [17] for definition) whose range is precisely a codimension one subspace of $H^{k-1}$ which annihilates the constant function. By the orthogonality condition in (3.15), the behaviour of $G(x, y)$ is uniquely characterized by its action on the range of $\Lambda$. To this end, for $f \in C^\infty(\partial M)$, denote its harmonic extension by $u_f$. Integrating by parts the expression $0 = \int_M \Delta_g u_f(x) G(x, y) dx$ for $z$ in the interior of $M$ we have

$$ (3.17) \qquad \begin{aligned} u(y) &= \int_{\partial M} \partial_{\nu} u_f(x) G(x, y) \mathrm{dvol}_h(x) + \frac{1}{|\partial M|} \int_{\partial M} f \mathrm{dvol}_h \\ &= \int_{\partial M} (\Lambda f)(x) G(x, y) \mathrm{dvol}_h(x) + \frac{1}{|\partial M|} \int_{\partial M} f \mathrm{dvol}_h \end{aligned} $$

Observe that any $\tilde{f} \in C^\infty(\partial M)$ has a unique decomposition $\tilde{f} = c + \Lambda f$ for some constant function $c$ and $f \in C^\infty(\partial M)$ satisfying $\int_{\partial M} f = 0$. Therefore, using the orthogonality condition of (3.15) and taking the trace of (3.17) we see that the map

$$ \tilde{f} \mapsto \left( \int_{\partial M} \tilde{f}(x) G(x, y) \mathrm{dvol}_h(x) \right) \Big|_{y \in \partial M} $$

is well defined and takes $C^\infty(\partial M) \to C^\infty(\partial M)$. We denote this operator by $G_{\partial M}$ and its Schwartz kernel by $G_{\partial M}(x, y)$. Going back to (3.17) we see that

$$ f = G_{\partial M} \Lambda f + P f, $$

where

$$ (3.18) \qquad Pf := |\partial M|^{-1} \int_{\partial M} f $$

is smoothing. In operator form this is

$$ (3.19) \qquad I = G_{\partial M} \Lambda + P. $$

Since $\Lambda \in \Psi_{cl}^1(\partial M)$ is elliptic (see [17]) we can conclude that $G_{\partial M} \in \Psi_{cl}^{-1}(\partial M)$ which maps $H^k(\partial M) \to H^{k+1}(\partial M)$ for all $k \in \mathbb{R}$. This completes the proof of Proposition 1.1 part i).
**Remark 3.2.** A quick way to prove part ii) of Proposition 1.1 would be to observe that (3.19) implies $G_{\partial M}$ is a parametrix for $\Lambda$. The symbol expansion for $\Lambda$ has already been computed in [17] so constructing its parametrix follows from standard pseudodifferential calculus. However, we will choose instead to take the layer potential approach since it will be more conducive for future numerical implementations. See Remark 3.3 below.

Applying on the right the single-layered potential $S$ defined in(3.13) and using identity (11.58) of [34] we have

$$ (3.20) \qquad G_{\partial M} = -2S + G_{\partial M}N^* + 2PS. $$

Iterating this equation and using intertwining property (11.59) of [34] we get

$$ (3.21) \qquad G_{\partial M} = -2S - 2NS + \Psi_{cl}^{-3}(\partial M). $$

By (2.8) the principal symbol of the operator $NS$ is simply the product of the principal symbols of $S$ with the principal symbol of $N$. The leading singularities of the operators $S$ and $N$ are given in (3.14) and the principal symbols of these kernels are computed in Lemmas 2.14 and 2.15. Therefore, using Proposition 2.16, we see that modulo $\Psi_{cl}^{-3}(\partial M)$, the integral kernel of $NS$ is given by

$$ \frac{1}{8\pi} H(x) \log d_h(y, x) - \frac{1}{32\pi} \left( \Pi_x \left( \frac{\exp_x^{-1}(y)}{|\exp_x^{-1}(y)|_h} \right) - \Pi_x \left( \frac{*\exp_x^{-1}(y)}{|\exp_x^{-1}(y)|_h} \right) \right) $$

when $x, y \in \partial M$ are close to each other.

Inserting this into (3.21) we get that when $x, y \in \partial M$ are close to each other,

$$ (3.22) \qquad G_{\partial M}(x, y) = \frac{1}{2\pi} d_g(x, y)^{-1} - \frac{1}{4\pi} H(x) \log d_h(y, x) \\ \qquad \qquad + \frac{1}{16\pi} \left( \Pi_x \left( \frac{\exp_x^{-1}(y)}{|\exp_x^{-1}(y)|_h} \right) - \Pi_x \left( \frac{*\exp_x^{-1}(y)}{|\exp_x^{-1}(y)|_h} \right) \right) + R(x, y), $$

where $R(x, y)$ is the Schwartz kernel of an operator in $\Psi_{cl}^{-3}(\partial M)$ which we call the regular part of $G(x, y)$. Observe that since the principal symbol of $R$ is in $S_{cl}^{-3}(T^*\partial M)$ and $\partial M$ is dimension 2, Sobolev embedding yields that

$$ (3.23) \qquad R(\cdot, \cdot) \in C^{0,\alpha}(\partial M \times \partial M) $$

for all $\alpha < 1$. The proof of Proposition 1.1 is now complete.

**Remark 3.3.** Note that (1.2) peels off the "singular part" of the distribution $G_{\partial M}(x, y)$ and gives us the representation

$$ G_{\partial M}(x, y) = G_{\text{sing}}(x, y) + R(x, y) $$

with the singularity structure of $G_{\text{sing}}$ explicitly given by (1.2). Inserting this representation of $G_{\partial M}$ into (3.20) gives the following integral equation for the regular part $R(x, y)$:

$$ R(I - N^*) = -G_{\text{sing}} - 2S + G_{\text{sing}}N^* + 2PS $$

where the operators $P$, $S$, and $N$ are given by (3.18) and (3.13).

Since $N^* \in \Psi_{cl}^{-1}(\partial M)$, it is a compact operator which makes $I - N^*$ Fredholm with index zero. Therefore, numerically computing for $R(x, y)$ amounts to solving a Fredholm boundary integral equation subject to the orthogonality condition

$$ \int_{\partial M} G_{\partial M}(x, z) d\text{vol}_h(z) = 0. $$
## 4. INVERTING THE NORMAL OPERATOR

Let $\Omega \subset \mathbb{R}^n$ be a bounded convex domain with smooth boundary. We will analyze the mapping properties of the operator

$$ (4.1) \qquad L: f \mapsto \int_{\Omega} \frac{f(s)}{|s-t|^{n-1}} ds_1 \dots ds_n $$

and its inverse. Methods do exist [18],[7] for the explicit expression of the inverse of $L$ when $\Omega = \mathbb{D}$ (which is sufficient for our setting). When $\Omega$ is a two dimensional ellipse [27] computed explicitly the inverse of $L$ acting on the constant function.

The purpose of Section 4.1 is to provide a geometric perspective to the operator $L$ one of the advantages of which is that it provides an explicit formula for $L^{-1}(1)$ when $\Omega$ is a ball of any dimension. Our perspective is based on some of the recent progress on integral geometry (in particular [23], [20], [12]). Since the explicit formulas and estimates will be valid in all dimensions, this will provide the key ingredient in proving Theorem 1.3 in all dimensions. When $\Omega$ is not necessarily $\mathbb{B}$, this geometric point of view may also potentially provide ways to relate some of the quantities of interest to the geometry of $\Omega$.

Section 4.2 will provide some explicit formulas for the composition of $L^{-1}$ with other operators in the case when $\Omega = \mathbb{D}$, the two dimensional disk. Section 4.3 will do the same for when $\Omega$ is the two dimensional ellipse although the formulas will not be as explicit.

### 4.1. Mapping Properties of L. Denote by

$$ \partial_{+} S\Omega := \{(x, v) \in \partial\Omega \times S^{n-1} \mid v \cdot \nu(x) \le 0\} $$

to be the set of inward pointing unit vectors on $\partial\Omega$. Note that this is a closed submanifold of the sphere bundle $S\Omega$ and thus inherits its smooth structure. Define the X-Ray transform $I: C^\infty(\bar{\Omega}) \to C^\infty(\partial_+ S\Omega)$ by $I f(x, v) := \int_0^{\tau(x,v)} f(x + tv) dt$ where $\tau(x, v)$ is the time it takes for a ray of unit velocity $v$ starting at $x \in \bar{\Omega}$ to reach the boundary $\partial\Omega$. Note that because $\Omega$ is assumed to be convex, $\tau(x, v)$ is a smooth function on $\partial_+ SM$. Furthermore, $I$ is injective by [23].

By [29] Theorem 4.2.1 this operator extends to an operator $I: L^2(\Omega) \to L^2_\mu(\partial_+ S\Omega)$ where $\mu$ is the measure given by $\mu = |\nu(x) \cdot v| d\text{vol}_{\partial\Omega} d\text{vol}_{S^{n-1}}$. This $L^2$ space mapping property allows us to define the adjoint operator $I^*$ given by (see [23])

$$ (4.2) \qquad I^* \omega(x) = \int_{S^{n-1}} \omega(x + \tau(x, v)v) d\text{vol}_{S^{n-1}}(v) $$

when acting on smooth functions $\omega$. This allows us to define a self-adjoint normal operator $I^* I: L^2(\Omega) \to L^2(\Omega)$. It turns out that by [23] the Schwartz kernel of $I^* I$ is precisely $2|s - t|^{-n+1}$ and therefore $I^* I = 2L$. Let $d_\Omega(\cdot)$ be any smooth positive function on $\Omega$ which is equal to $\text{dist}(\cdot, \partial\Omega)$ near the boundary. By Theorem 2.2 and Theorem 4.4 of [20] respectively, we have that

$$ (4.3) \qquad I^* I: d_\Omega^{-1/2} C^\infty(\bar{\Omega}) \to C^\infty(\bar{\Omega}) $$

is a bijection and

$$ (4.4) \qquad 2L = I^* I: \{u \in H^{-1/2}(\mathbb{R}^n) \mid \text{supp}(u) \subset \Omega\} \simeq H^{1/2}(\Omega)^* \to H^{1/2}(\Omega) $$

is a self-adjoint homeomorphism. Thus there exists a unique function $u_0 \in d_\Omega^{-1/2} C^\infty(\bar{\Omega})$ such that $Lu_0 = 1$ which is equivalent to $I^* I u_0 = 2$. To compute $u_0$, observe that if
we find $u_0$ such that

$$ (4.5) \qquad I u_0(x, v) = \frac{2}{\mathrm{Vol}(S^{n-1})} $$

for all $(x, v) \in \partial_+ S\Omega$ then by (4.2) we would have $I^*Iu_0 = 2$. The solution of (4.5) is easy to compute explicitly when $\Omega = \mathbb{B}$. Indeed, direct computation shows that choosing

$$ (4.6) \qquad u_0(x) = \frac{2}{\pi \mathrm{Vol}(S^{n-1}) \sqrt{1 - |x|^2}} $$

one satisfies (4.5). In particular if $n = 2$ (which is the case we are interested in) we have that

$$ (4.7) \qquad L^{-1}(1) = u_0(x) = \frac{1}{\pi^2 \sqrt{1 - |x|^2}}, \quad \text{in dimension 2} $$

**Remark 4.1.** This process of computing solution to $I^*Iu_0 = \text{const}$ by solving for $Iu_0 = \text{const}$ unfortunately only works for $\Omega = \mathbb{B}$. In fact, thanks to the rigidity result of [12], we know that $Iu_0 = 1$ is solvable iff $\Omega = \mathbb{B}$. However, for more general domains the geometric view presented here could potentially allow one to apply the reconstruction formula for inverting $I$ [22] to solve $I^*Iu_0 = 1$ explicitly. To do so one must first invert $I^*$ (which has a large kernel but is surjective) into the range of $I$ and it is not clear how to do this when $\Omega \neq \mathbb{B}$.

4.2. **Integrals Involving L inverse of 1.** We also define $R_{\log}$ and $R_{\infty}$ to be operators with kernels $\log|s' - t'|$ and $\frac{(s_1-t_1)^2-(s_2-t_2)^2}{|s'-t'|^2}$ respectively.

**Lemma 4.2.** *The operators $R_{\infty}$ and $R_{\log}$ are bounded maps from $H^{1/2}(\mathbb{D})^*$ to $H^{3/2}(\mathbb{D})$.*

*Proof.* Observe that the integral kernels of both $R_{\log}$ and $R_{\infty}$ extends naturally to kernels representing operators in $\Psi^{\infty}(\mathbb{R}^2)$ which we denote by $\tilde{R}_{\log}$ and $\tilde{R}_{\infty}$ respectively. We denote by $E: H^{1/2}(\mathbb{D})^* \to \{u \in H^{-1/2}(\mathbb{R}^2) \mid \mathrm{supp}(u) \subset \Omega\}$ to be the isomorphism obtained by the trivial extension. Let $\chi \in C_0^{\infty}(\mathbb{R}^2)$ be identically 1 on $\mathbb{D}$. Then we have that

$$ (4.8) \qquad (R_{\log}u)|_{\mathbb{D}} = (\chi \tilde{R}_{\log} \chi Eu)|_{\mathbb{D}} $$

and the same holds for $R_{\infty}$.

By Proposition 7.2.8 of [34] we have that both $\chi \tilde{R}_{\log} \chi$ and $\chi \tilde{R}_{\infty} \chi$ are pseudodifferential operators of order $-2$. Therefore by (4.8) both $R_{\log}$ and $R_{\infty}$ are bounded operators from $H^{1/2}(\mathbb{D})^*$ to $H^{3/2}(\mathbb{D})$.

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

The following lemma was proved in [7, Theorem 4.2].

**Lemma 4.3.** For $u \in H^{\frac{1}{2}}(\mathbb{D})$, it follows

$$ (4.9) \qquad \langle L^{-1}u, 1 \rangle = \frac{1}{\pi^2} \int_{\mathbb{D}} \frac{u(t')}{(1 - |t'|^2)^{\frac{1}{2}}} dt_1 dt_2. $$

*Proof.* By (4.4) $L: H^{\frac{1}{2}}(\mathbb{D})^* \to H^{\frac{1}{2}}(\mathbb{D})$ is a self-adjoint homeomorphism. The result of this Lemma is a direct consequence.

<figure><img src="image_3.png" /></figure>
**Lemma 4.4.** Let

$$f(s') = R_{\log} L^{-1} 1 = \int_{\mathbb{D}} \log |t' - s'| [L^{-1} 1](t') dt_1 dt_2,$$

then

$$f(s') = \frac{2}{\pi} \log |s'| + \frac{2}{\pi} \left( \frac{1}{2} \log \left| (1 - |s'|^2)^{\frac{1}{2}} + 1 \right| - \frac{1}{2} \log \left| (1 - |s'|^2)^{\frac{1}{2}} - 1 \right| \right) - \frac{2}{\pi} (1 - |s'|^2)^{\frac{1}{2}}.$$

*Proof.* Note that $\frac{1}{2\pi} \log |t' - s'|$ is the fundamental solution for the Laplace operator in $\mathbb{R}^2$, therefore,

$$\Delta f = 2\pi[L^{-1}1] = \frac{2}{\pi} \frac{1}{(1 - |t'|^2)^{\frac{1}{2}}}$$

Since $[L^{-1}1](t')$ is radially symmetric, $f(t') = \tilde{f}(r)$ where $r = |t'|$. Writing the Laplace operator in polar coordinates, we get

$$(r \tilde{f}_r)_r = \frac{2}{\pi} \frac{r}{(1 - r^2)^{\frac{1}{2}}}.$$

Integration gives

$$\tilde{f}_r(r) = \frac{2}{\pi} \left( \frac{C_1}{r} - \frac{(1 - r^2)^{\frac{1}{2}}}{r} \right)$$

and

$$\tilde{f}(r) = \frac{2}{\pi} \left( C_1 \log r - \frac{1}{2} \log \left| (1 - r^2)^{\frac{1}{2}} - 1 \right| + \frac{1}{2} \log \left| (1 - r^2)^{\frac{1}{2}} + 1 \right| - (1 - r^2)^{\frac{1}{2}} + C_2 \right). \quad (4.10)$$

Let us find $C_1$ and $C_2$. Note that $\tilde{f}(r)$ does not have singularity at $r = 0$, namely,

$$\tilde{f}(0) = f(0) = \frac{1}{\pi^2} \int_{\mathbb{D}} \frac{\log |t'|}{(1 - |t'|^2)^{\frac{1}{2}}} = \frac{2}{\pi} \int_0^1 \frac{r \log r}{(1 - r^2)^{\frac{1}{2}}} dr = \frac{2}{\pi} (\log 2 - 1).$$

Therefore, the identities

$$C_1 \log r - \frac{1}{2} \log \left| (1 - r^2)^{\frac{1}{2}} - 1 \right| = \frac{1}{2} \log \left| \frac{r^{2C_1}}{(1 - r^2)^{\frac{1}{2}} - 1} \right| = \frac{1}{2} \log \left| \frac{r^{2C_1}}{-\frac{1}{2}r^2 + O(r^4)} \right|,$$

as $r \to 0$, implies that $C_1 = 1$. Hence, putting $r = 0$ into (4.10), gives

$$\frac{2}{\pi}(\log 2 - 1) = \frac{2}{\pi} \left( \frac{1}{2} \log 2 + \frac{1}{2} \log 2 - 1 + C_2 \right),$$

so that $C_2 = 0$.

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

Due to Lemmas 4.4 and 4.3, the following identity is a direct computation:

**Lemma 4.5.** The following identity holds

$$\langle L^{-1} R_{\log} L^{-1} 1, 1 \rangle = \frac{8}{\pi^2} \log 2 - \frac{6}{\pi^2}.$$

**Lemma 4.6.** The following identity holds

$$\langle L^{-1} R_{\infty} L^{-1} 1, 1 \rangle = 0.$$
*Proof.* By Lemmas (4.7) and 4.3, we know that

$$ \langle L^{-1} R_{\infty} L^{-1}, 1 \rangle = \frac{1}{\pi^2} \int_{\mathbb{D}} \int_{\mathbb{D}} \frac{(s_1 - t_1)^2 - (s_2 - t_2)^2}{|s' - t'|^2} \frac{1}{(1 - |s'|^2)^{\frac{1}{2}}} ds' \frac{1}{(1 - |t'|^2)^{\frac{1}{2}}} dt'. $$

Consider the following two changes of variables for the right-hand side

$$ (s_1, s_2, t_1, t_2) = (r \cos \phi, r \sin \phi, \rho \cos \theta, \rho \sin \theta), $$

$$ (s_1, s_2, t_1, t_2) = (r \sin \phi, r \cos \phi, \rho \sin \theta, \rho \cos \theta). $$

The results differ by multiplying by $-1$, which means that the right-hand side is 0. $\square$

**4.3. Explicit Formulas in 2 Dimensional Ellipse.** We now compute the inverse of the map $f \mapsto \int_{E_a} \frac{f(s')}{|s' - t'|} ds'$ where the domain of integration is the two dimensional ellipse $E_a := \{s_1^2 + \frac{s_2^2}{a^2} = 1\}$ instead of a ball. A change of variable leads us to consider the operator

$$ (4.11) \qquad L_a f = a \int_{\mathbb{D}} \frac{f(s')}{((t_1 - s_1)^2 + a^2(t_2 - s_2)^2)^{1/2}} ds' $$

acting on functions of the disk $\mathbb{D}$. By [27] we have that

$$ (4.12) \qquad L_a (K_a^{-1} (1 - |t'|^2)^{-1/2}) = 1, $$

on $\mathbb{D}$ where

$$ K_a = \frac{\pi}{2} \int_0^{2\pi} \frac{1}{\left(\cos^2 \theta + \frac{\sin^2 \theta}{a^2}\right)^{1/2}} d\theta. $$

By (4.4) this is the unique solution in $H^{1/2}(\mathbb{D})^*$ to $L_a u = 1$.

Next we denote

$$ R_{\log,a} f(t') := a \int_{\mathbb{D}} \log \left( (t_1 - s_1)^2 + a^2 (t_2 - s_2)^2 \right)^{1/2} f(s') ds', $$

$$ R_{\infty,a} f(t') := a \int_{\mathbb{D}} \frac{(t_1 - s_1)^2 - a^2(t_2 - s_2)^2}{(t_1 - s_1)^2 + a^2(t_2 - s_2)^2} f(s') ds'. $$

For general $a$, the quantities $\langle L_a^{-1}(1), \mathbb{R}_{\log,a} L_a^{-1}(a) \rangle$ and $\langle L_a^{-1}(1), \mathbb{R}_{\infty,a} L_a^{-1}(a) \rangle$ cannot be computed as explicitly as in the case when $a = 1$ in Section 4.2.

# 5. ASYMPTOTIC EXPANSION OF THE SINGULARLY PERTURBED PROBLEMS

## 5.1. Asymptotic Expansion of Mixed Boundary Value Problems.
We are now ready to compute the asymptotic expansion for the mixed boundary value problem $u_{\epsilon} \in H^1(M)$,

$$ (5.1) \qquad \Delta_g u_{\epsilon} = -1, \quad u_{\epsilon}|_{\Gamma_{\epsilon,a}} = 0, \quad \partial_{\nu} u_{\epsilon}|_{\partial M \setminus \Gamma_{\epsilon,a}} = 0 $$

which gives the compatibility condition

$$ (5.2) \qquad \int_{\partial M} \partial_{\nu} u_{\epsilon} d\text{vol}_h = -|M|. $$

All integrals on open subsets of $\partial M$ are with respect to the volume form given by the metric $h$.
Using Green's formula as in [9], also in [1], we can deduce that for points $x \in M^o$,
$u_\epsilon(x)$ satisfies the integral equation

$$
(5.3) \qquad u_{\epsilon}(x) = F(x) + C_{\epsilon,a} + \int_{\partial M} G(x, y) \partial_{\nu} u_{\epsilon}(y) \mathrm{dvol}_{h}(y),
$$

where $C_{\epsilon,a} = |\partial M|^{-1} \int_{\partial M} u_{\epsilon}$ and $F(x) = \int_{M} G(x, y)$ solves the boundary value problem

$$
(5.4) \qquad \Delta F = -1, \quad \partial_{\nu} F = -|M|/|\partial M|, \quad \int_{\partial M} F \, \mathrm{dvol}_h = 0.
$$

By Proposition 1.1 we can take the trace of (5.3) to the boundary and restrict to the
open subset $\Gamma^o_{\epsilon,a} \subset \partial M$. Using (5.1) we see that

$$
0 = F(x) + C_{\epsilon,a} + \int_{\Gamma_{\epsilon,a}} G_{\partial M}(x,y) \partial_{\nu} u_{\epsilon}(y) \mathrm{dvol}_h(y)
$$

for $x \in \Gamma^o_{\epsilon,a}$. We now replace $G_{\partial M}(x, y)$ for $x, y \in \Gamma^o_{\epsilon,a}$ with the expression in (1.2) to obtain

$$
-F(x) - C_{\epsilon,a} ={}& \frac{1}{2\pi} \int_{\Gamma_{\epsilon,a}} d_g(x,y)^{-1} \partial_\nu u_\epsilon(y) \mathrm{dvol}_h(y) - \frac{H(x)}{4\pi} \int_{\Gamma_{\epsilon,a}} \log d_h(x,y) \partial_\nu u_\epsilon(y) \mathrm{dvol}_h(y) \\
(5.5) \qquad & + \frac{1}{16\pi} \int_{\Gamma_{\epsilon,a}} \left( \Pi_x \left( \frac{\exp_x^{-1}(y)}{|\exp_x^{-1}(y)|_h} \right) - \Pi_x \left( \frac{*\exp_x^{-1}(y)}{|\exp_x^{-1}(y)|_h} \right) \right) \partial_\nu u_\epsilon(y) \mathrm{dvol}_h(y) \\
& + \int_{\Gamma_{\epsilon,a}} R(x,y) \partial_\nu u_\epsilon(y) \mathrm{dvol}_h(y).
$$

We will write this integral equation in the coordinate system given by

$$
(5.6) \qquad \mathbb{D} \ni (s_1, s_2) \mapsto x^{\epsilon}(s_1, as_2; x^{*}) \in \Gamma_{\epsilon,a},
$$

where $x^\epsilon(\cdot; x^*) : \mathcal{E}_a \to \Gamma_{\epsilon,a}$ is the coordinate defined in Section 2.2. To simplify notation we will drop the $x^*$ in the notation and denote $x^\epsilon(\cdot; x^*)$ by simply $x^\epsilon(\cdot)$.

Note that in these coordinates the volume form for $\partial M$ is given by

$$
(5.7) \qquad \mathrm{dvol}_h = a\epsilon^2 (1 + \epsilon^2 Q_\epsilon(s')) ds_1 \wedge ds_2, \quad s' \in \mathbb{D}
$$

for some smooth function $Q_{\epsilon}(s')$ whose derivatives of all orders are bounded uniformly
in $\epsilon$. We denote

$$
(5.8) \qquad \psi_{\epsilon}(s') := \partial_{\nu} u_{\epsilon}(x^{\epsilon}(s_{1}, a s_{2})).
$$

The compatibility condition (5.2) written using the expression for the volume form
(5.7) is then

$$
(5.9) \qquad \int_{\mathbb{D}} \psi_{\epsilon}(s')(1 + \epsilon^2 Q_{\epsilon}(s')) ds_1 ds_2 = - \frac{|M|}{a\epsilon^2}.
$$

Let us unwrap the right hand side of (5.5) term by term in the coordinate given
by $x^\epsilon(\cdot)$. Write out the integral of the first term using the expression of the volume
form (5.7) and the expression for $d_g(x, y)^{-1}$ in Corollary 2.6 and taking into account
that the coordinate system is scaled by a factor $a$ as in (5.6) gives
(5.10)

$$
\int_{\Gamma_{\epsilon,a}} d_g(x, y)^{-1} \partial_{\nu} u_{\epsilon}(y) \mathrm{dvol}_h(y) = a\epsilon \int_{\mathbb{D}} \frac{1}{((t_1 - s_1)^2 + a^2(t_2 - s_2)^2)^{1/2}} \psi_{\epsilon}(s') ds + \epsilon^3 \mathcal{A}_{\epsilon} \psi_{\epsilon},
$$
for some operator $\mathcal{A}_\epsilon$ whose Schwartz kernel is given by the second term of the expansion in Corollary 2.6. Here due to Lemma 2.11 we have that

$$ \mathcal{A}_{\epsilon} : H^{1/2}(\mathbb{D}; ds')^{*} \rightarrow H^{1/2}(\mathbb{D}; ds') $$

with operator norm bounded uniformly in $\epsilon$. From here on we will denote by $\mathcal{A}_\epsilon$ any operator which takes $H^{1/2}(\mathbb{D}; ds')^* \rightarrow H^{1/2}(\mathbb{D}; ds')$ whose operator norm is bounded uniformly in $\epsilon$.

Doing the same thing for the second term of (5.5) while using Lemma 2.2, Lemma 2.11, and (5.9) gives

$$ H(x) \int_{\Gamma_{\epsilon,a}} \log d_h(x, y) \partial_{\nu} u_{\epsilon}(y) \mathrm{dvol}_h(y) = -H(x^*) |M| \log \epsilon + \epsilon^2 H(x^*) R_{\log,a} \psi_{\epsilon} + \epsilon^3 \mathcal{A}_{\epsilon} \psi_{\epsilon} \\ (5.11) \qquad \qquad \qquad \qquad \qquad + O_{H^{1/2}(\mathbb{D})} (\epsilon \log \epsilon) $$

where $R_{\log,a}$ is defined at the very beginning of Section 4.2. Here $O_{H^{1/2}(\mathbb{D})}(\epsilon \log \epsilon)$ denotes a function with $H^{1/2}(\mathbb{D}; ds')$ norm vanishing to order $\epsilon \log \epsilon$. Note the volume for we use here is now the Euclidean one rather than $\mathrm{dvol}_h$ given by (5.7).

Finally, for the third term of (5.5) we get by using the coordinate expression derived in Corollary 2.9 and the estimate of Lemma 2.11:

$$ \int_{\Gamma_{\epsilon,a}} \left( \mathrm{II}_x \left( \frac{\exp_x^{-1}(y)}{|\exp_x^{-1}(y)|_h} \right) - \mathrm{II}_x \left( \frac{*\exp_x^{-1}(y)}{|\exp_x^{-1}(y)|_h} \right) \right) \partial_{\nu} u_{\epsilon}(y) = \epsilon^2 (\lambda_1 - \lambda_2) R_{\infty,a} \psi_{\epsilon} \\ (5.12) \qquad \qquad \qquad \qquad \qquad + \epsilon^3 \mathcal{A}_{\epsilon} \psi_{\epsilon}, $$

where $R_{\infty,a}$ is defined in Section 4. Inserting into (5.5) the identities (5.10), (5.11), and (5.12) we have

$$ \frac{-F(x^*) - C_{\epsilon,a} - H(x^*) |M| (4\pi)^{-1} \log \epsilon}{\epsilon} = \frac{1}{2\pi} L_a \psi_{\epsilon} - \epsilon \frac{H(x^*)}{4\pi} R_{\log,a} \psi_{\epsilon} + \epsilon \frac{\lambda_1 - \lambda_2}{16\pi} R_{\infty,a} \psi_{\epsilon} \\ (5.13) \qquad \qquad \qquad \qquad \qquad + a\epsilon \int_{\mathbb{D}} R(x^{\epsilon}(t'), x^{\epsilon}(s')) \psi_{\epsilon}(s') + \epsilon^2 \mathcal{A}_{\epsilon} \psi_{\epsilon} \\ \qquad \qquad \qquad \qquad \qquad + O_{H^{1/2}(\mathbb{D})} (\log \epsilon). $$

We would like to approximate the integral kernel $R(x^\epsilon(t'), x^\epsilon(s'))$ by the constant $R(x^*, x^*)$ and this is the content of

**Lemma 5.1.** Let $T_\epsilon : C_c^\infty(\mathbb{D}) \rightarrow \mathcal{D}'(\mathbb{D})$ be the operator defined by the integral kernel

$$ R(x^{\epsilon}(t'; x^{*}), x^{\epsilon}(s'; x^{*})) - R(x^{*}, x^{*}). $$

Then

$$ \|T_{\epsilon}\|_{(H^{1/2}(\mathbb{D}))^{*} \rightarrow H^{1/2}(\mathbb{D})} \leq C_{\epsilon} \log \epsilon. $$

*Proof.* Set $T(t', s') := R(x(t'; x^*), x(s'; x^*)) - R(x^*, x^*)$ for $t'$ and $s'$ small and extend it to be a smooth compactly supported kernel otherwise. The kernel for $T_\epsilon$ is then $T(\epsilon t', \epsilon s')$. Note that

$$ (5.14) \qquad T(0, 0) = 0. $$

Observe that if $\chi \in C_c^\infty(\mathbb{R}^2)$ which is identically 1 on $\mathbb{D}$ then the operator $T_\epsilon$ acting on distributions supported in $\mathbb{D}$ is given by the Schwartz kernel

$$ \chi(s')\chi(t')T(\epsilon t', \epsilon s') $$

for $t', s' \in \mathbb{D}$ and $\epsilon > 0$ small.
Observe that $T(t', s')$ is the integral kernel for an operator in $\Psi_{cl}^{-3}(\mathbb{R}^2)$. Applying Prop 2.8 in Chap 7 of [34] we can deduce that for all $k$ we may choose $N$ sufficiently large such that

$$T(t', s') = \sum_{l=0}^{N} (q_l(t', t' - s') + p_l(t', t' - s') \log |t' - s'|) + R_k(t', s'),$$

where for each integer $l$ and multi-index $\gamma$, $D_{t'}^{\gamma} q_l(t', \cdot)$ is a bounded continuous function of $t'$ with value in the space of smooth (away from the origin) homogeneous distributions of degree $l+1$, $p_l(t', \cdot)$ are homogenous polynomials of degree $l+1$ with coefficients which are smooth functions of $t'$, and for all multi-indices $\gamma$, $D_{t'}^{\gamma} R_k(t', \cdot)$ bounded continuous function of $t'$ with value in $C^k(\mathbb{R}^2)$.

Using (5.14) along with the homogenity degree of $q_l$ and $p_l$ we see that $R_k(0, 0) = 0$. Therefore, for $s', t' \in \mathbb{D}$ the integral kernel of $T_\epsilon$ is

$$T_\epsilon(t', s') = \sum_{l=0}^{N} \left( \epsilon^{l+1} q_l(\epsilon t', t' - s') + \epsilon^{l+1} p_l(\epsilon t', t' - s') \log \epsilon + \epsilon^{l+1} p_l(\epsilon t', t' - s') \log |t' - s'| \right) \\ (5.15) \qquad + R_k(\epsilon t', \epsilon s').$$

The kernel $R_k(s', t')$ is sufficiently smooth. Therefore, by doing a Taylor expansion and using $R_k(0, 0) = 0$ we see that the integral kernel $\chi(s')\chi(t')R_k(\epsilon t', \epsilon s')$ takes $H^{1/2}(\mathbb{D})^* \to H^{1/2}(\mathbb{D})$ with norm $\epsilon$. The worst term in the polyhomogeneous expansion part of (5.15) happens when $l = 0$ and this term is given by

$$\epsilon q_0(\epsilon t', z') + p_0(\epsilon t', z') \log |z'| + \epsilon \log \epsilon p_0(\epsilon t', z'),$$

where $z' = t' - s'$. Recall that both $q_0$ and $p_0$ are homogeneous of degree 1 in $z'$ so writing $z' = r\omega$ then applying Lemma 2.11 we have uniform estimates in $\epsilon$ for the kernels $\chi(s')\chi(t')q_0(\epsilon t', t' - s')$ and $\chi(s')\chi(t')p_0(\epsilon t', t' - s')$. For the term involving $\log|s' - t'|$, we use the fact that $p_0(t', z')$ is a linear function in $z'$ whose coefficients are smooth functions of $t'$. Therefore, if $u \in C_c^\infty(\mathbb{R}^2)$,

$$\int_{\mathbb{R}^2} \chi(t') p_0(\epsilon t', t' - s') \log |t' - s'| u(s') ds' = \chi(t') \int_{\mathbb{R}^2} \sum_j c_j(\epsilon t')(t_j - s_j) \log |s' - t'| u(s') ds' \\ = \chi(t') \sum_j c_j(\epsilon t') \int_{\mathbb{R}^2} a_j(\xi') \hat{u}(\xi') e^{it' \cdot \xi'} d\xi',$$

where $a_j(\cdot) \in S'(\mathbb{R}^2)$ for $j = 1, 2$ are derivatives of the distribution $\text{PF}|\xi|^{-2}$ with respect to $\partial_{\xi_j}$. We refer the reader to (8.31) in Chapt 3 of [33] for the definition of the distribution $\text{PF}|\xi|^{-2}$. From this we see that the integral kernel

$$\chi(t')\chi(s')p_0(\epsilon t', t' - s') \log |t' - s'|$$

also maps $H^{1/2}(\mathbb{D})^* \to H^{1/2}(\mathbb{D})$ with uniform bound in $\epsilon$. $\square$

Due to Lemma 5.1 we can write (5.13) as

$$\frac{2\pi}{\epsilon} \left( aR(x^*, x^*)|M| - F(x^*) - C_{\epsilon,a} - \frac{H(x^*)|M| \log \epsilon}{4\pi} \right) = \left( L_a - \frac{\epsilon H(x^*)}{2} R_{\log,a} + \frac{\epsilon(\lambda_1 - \lambda_2)}{8} R_{\infty,a} \right) \psi_\epsilon \\ + \epsilon T_\epsilon \psi_\epsilon + + \epsilon^2 A_\epsilon \psi_\epsilon + O_{H^{1/2}(\mathbb{D})}(\log \epsilon). \qquad (5.16)$$
We hit both sides with $L_a^{-1}$ and use (4.12) and (4.4) to get the identity

$$ (5.17) \quad \frac{2\pi}{\epsilon} \left( aR(x^*, x^*)|M| - F(x^*) - C_{\epsilon,a} - \frac{H(x^*)|M|\log\epsilon}{4\pi} \right) \frac{1}{K_a(1 - |t'|^2)^{1/2}} = \left( I - \frac{\epsilon H(x^*)}{2} L_a^{-1} R_{\log,a} + \frac{\epsilon(\lambda_1 - \lambda_2)}{8} L_a^{-1} R_{\infty,a} + T'_\epsilon \right) \psi_\epsilon + O_{H^{1/2}(\mathbb{D})^*}(\log \epsilon). $$

for some $T'_\epsilon: H^{1/2}(\mathbb{D})^* \to H^{1/2}(\mathbb{D})^*$ with operator norm $O(\epsilon^2 \log \epsilon)$. Use the mapping properties from Lemma 4.2 we see that the right side can be inverted by Neumann series to deduce

$$ (5.18) \quad \psi_{\epsilon} = - \frac{2\pi C_{\epsilon,a}}{\epsilon K_a (1 - |t'|^2)^{1/2}} + C_{\epsilon,a} O_{H^{1/2}(\mathbb{D})^*}(1) + O_{H^{1/2}(\mathbb{D})^*}(\epsilon^{-1} \log \epsilon). $$

Insert (5.18) into (5.9) we get that

$$ (5.19) \quad C_{\epsilon,a} = \frac{|M|K_a}{4a\epsilon\pi^2} + C'_{\epsilon,a} $$

with $C'_{\epsilon,a} = O(\log \epsilon)$. Insert (5.19) into (5.18)

$$ (5.20) \quad \psi_{\epsilon} = \frac{-|M|}{2a\pi\epsilon^2(1 - |t'|^2)^{1/2}} + \psi'_{\epsilon} = -\frac{|M|K_a}{2a\pi\epsilon^2} L_a^{-1}(1) + \psi'_{\epsilon} $$

with $\|\psi'_\epsilon\|_{H^{1/2}(\mathbb{D};ds')^*} \le C\epsilon^{-1}\log\epsilon$. Insert (5.19) and (5.20) into (5.17) we get

$$ (5.21) \quad \frac{2\pi}{\epsilon} \left( aR(x^*, x^*)|M| - F(x^*) - C'_{\epsilon,a} - \frac{H(x^*)|M|\log\epsilon}{4\pi} \right) \frac{1}{K_a(1 - |t'|^2)^{1/2}} = \psi'_\epsilon + \frac{|M|K_a}{2\pi\epsilon} \left( \frac{H(x^*)}{2} L_a^{-1} R_{\log,a} - \frac{(\lambda_1 - \lambda_2)}{8} L_a^{-1} R_{\infty,a} \right) L_a^{-1}(1) + O_{H^{1/2}(\mathbb{D})^*}(\log \epsilon). $$

Inserting the expression (5.20) into (5.9) we get that

$$ (5.22) \quad \int_{\mathbb{D}} \psi'_{\epsilon}(s')ds_1ds_2 = O(1). $$

Multiply (5.21) by $\epsilon$ then integrate over $\mathbb{D}$. Then (5.22) implies

$$ (5.23) \quad C'_{\epsilon,a} = -\frac{1}{4\pi}H(x^*)|M|\log\epsilon + aR(x^*, x^*)|M| - F(x^*) \\ - \frac{|M|H(x^*)K_a^2}{16\pi^3} \int_{\mathbb{D}} L_a^{-1}R_{\log,a}L_a^{-1}1(s')ds' \\ + \frac{|M|(\lambda_1 - \lambda_2)K_a^2}{64\pi^3} \int_{\mathbb{D}} L_a^{-1}R_{\infty,a}L_a^{-1}1(s')ds' + O(\epsilon \log \epsilon). $$

Since $L_a^{-1}$ is self-adjoint, we can express the last two integrals more explicitly:

$$ \int_{\mathbb{D}} L_a^{-1} R_{\log,a} L_a^{-1} 1(s') ds' = K_a^{-2} \langle (1 - |s'|^2)^{-1/2}, R_{\log,a} (1 - |s'|^2)^{-1/2} \rangle, $$

$$ \int_{\mathbb{D}} L_a^{-1} R_{\infty,a} L_a^{-1} 1(s') ds' = K_a^{-2} \langle (1 - |s'|^2)^{-1/2}, R_{\infty,a} (1 - |s'|^2)^{-1/2} \rangle. $$

We summarize this calculation into the following:
**Proposition 5.2.** We have that

$$ \psi_{\epsilon} = \frac{-|M|}{2a\pi\epsilon^2(1 - |t'|^2)^{1/2}} + \psi'_{\epsilon} $$

with $\psi'_{\epsilon} = O_{H^{1/2}(\mathbb{D})}(\epsilon^{-1} \log \epsilon)$. Furthermore

$$ (5.24) \qquad C_{\epsilon,a} = \frac{|M|K_a}{4a\epsilon\pi^2} - \frac{1}{4\pi}H(x^*)|M|\log\epsilon + aR(x^*, x^*)|M| - F(x^*) \\ - \frac{|M|H(x^*)}{16\pi^3} \langle (1 - |s'|^2)^{-1/2}, R_{\log,a}(1 - |s'|^2)^{-1/2} \rangle \\ + \frac{|M|(\lambda_1 - \lambda_2)}{64\pi^3} \langle (1 - |s'|^2)^{-1/2}, R_{\infty,a}(1 - |s'|^2)^{-1/2} \rangle \\ + O(\epsilon \log \epsilon), $$

where $F$ is the unique solution to (5.4) and $R(x^*, x^*)$ is the evaluation at $(x, y) = (x^*, x^*)$ of the kernel $R(x, y)$ in (1.2).

Observe that in the case of the disc (i.e. $a = 1$) we have that

$$ \langle (1 - |s'|^2)^{-1/2}, R_{\log,a}(1 - |s'|^2)^{-1/2} \rangle = \pi^2 (8 \log 2 - 6) $$

by Lemma 4.5 and

$$ \langle (1 - |s'|^2)^{-1/2}, R_{\infty,a}(1 - |s'|^2)^{-1/2} \rangle = 0 $$

by Lemma 4.6. Thus the formula (5.24) simplifies to

$$ (5.25) \qquad C_{\epsilon} := C_{\epsilon,1} = \frac{|M|K_a}{4a\epsilon\pi^2} - \frac{1}{4\pi}H(x^*)|M|\log\epsilon + aR(x^*, x^*)|M| - F(x^*) \\ - \frac{|M|H(x^*)}{16\pi}(8\log 2 - 6) + O(\epsilon \log \epsilon), $$

when $a = 1$.

**5.2. Proof of Theorems 1.2 and 1.3.** We now prove Theorem 1.3. Theorem 1.2 follows from the explicit expression for $C_{\epsilon}$ in (5.25).

By the result of Appendix A we have that $u_{\epsilon} = \mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0 = x]$ solves the mixed boundary value problem (5.1) so using Proposition 5.2, (5.3), and (5.7), the expansion for $\mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0 = x]$ is given by

$$ u_{\epsilon}(x) = \mathbb{E}[\tau_{\Gamma_{\epsilon,a}}|X_0 = x] = F(x) + C_{\epsilon,a} - |M|G(x, x^*) + r_{\epsilon}(x) $$

for each $x \in M \setminus \Gamma_{\epsilon,a}$. Here $F$ is the unique solution to (5.4) and the remainder $r_{\epsilon}$ is given by

$$ (5.26) \qquad r_{\epsilon}(x) = \int_{\partial M} (G(x, y) - G(x, x^*)) \partial_{\nu} u_{\epsilon}(y) d\mathrm{vol}_{h}(y). $$

Let $K \subset \bar{M}$ be a compact subset of $\bar{M}$ which has positive distance from $x^*$ and consider $x \in K$. Writing out this integral in the $x^\epsilon(\cdot; x^*)$ coordinate system and use (5.8), (5.7), and the expression of $\psi_{\epsilon}$ derived in Proposition 5.2 we get

$$ (5.27) \qquad r_{\epsilon}(x) = \epsilon \int_{\mathbb{D}} \frac{-|M|}{2\pi(1 - |s'|^2)^{1/2}} L(x, \epsilon s')(1 + \epsilon^2 Q_{\epsilon}(s')) ds' \\ + a\epsilon^3 \int_{\mathbb{D}} \psi'_{\epsilon}(s') L(x, \epsilon s')(1 + \epsilon^2 Q_{\epsilon}(s')) ds' $$
for some function $L(x, s')$ jointly smooth in $(x, s') \in K \times \mathbb{D}$. The second integral formally denotes the duality between $H^{1/2}(\mathbb{D})^*$ and $H^{1/2}(\mathbb{D})$. The estimate for $\psi'_\epsilon$ derived in Proposition 5.2 now gives for any integer $k$ and any compact set $K$ not containing $x^*$, $\|r_\epsilon\|_{C^k(K)} \le C_{k,K\epsilon}$.

Our pseudodifferential characterization of $G_{\partial M}$ also allows us to compute the asymptotic of the average $\int_M u_\epsilon$. Indeed, integrating (5.3) over $M$ we get (5.28)

$$ \int_M u_\epsilon \mathrm{dvol}_g = \int_M F(x) \mathrm{dvol}_g + C_{\epsilon,a} |M| + \int_M \int_{\Gamma_{\epsilon,a}} G(x, y) \partial_\nu u_\epsilon(y) \mathrm{dvol}_h(y) \mathrm{dvol}_g(x). $$

We compute the last integral by noting that

$$ v(x) := \int_{\Gamma_{\epsilon,a}} G(x, y) \partial_\nu u_\epsilon(y) \mathrm{dvol}_h(y) $$

is the unique solution to the Dirichlet boundary value problem:

$$ (5.29) \quad \Delta_g v = 0, \quad v(x) \Big|_{\partial M} = \int_{\Gamma_{\epsilon,a}} G_{\partial M}(x, y) u_\epsilon(y) \mathrm{dvol}_h(y) \in H^{1/2}(\partial M). $$

We concluded the boundary value is in $H^{1/2}$ because $G_{\partial M} \in \Psi_{cl}^{-1}(\partial M)$ by (1.2).

Let a sequence of smooth functions $f_j \to \partial_\nu u_\epsilon$ in $H^{1/2}(\partial M)$ and let $v_j$ solve

$$ \Delta_g v_j = 0, \quad v_j(x) \Big|_{\partial M} = \int_{\partial M} G_{\partial M}(x, y) f_j(y) \mathrm{dvol}_h(y). $$

Standard elliptic theory shows that $v_j \to v$ in $H^1(M)$. Therefore

$$ \begin{align*} \int_M \int_{\Gamma_{\epsilon,a}} G(x, y) \partial_\nu u_\epsilon(y) &= \lim_{j} \int_M \int_{\partial M} G(x, y) f_j(y) \\ &= \lim_{j} \int_{\partial M} f_j(y) \int_M G(x, y) \\ &= \lim_{j} \int_{\partial M} f_j(y) F(y) \\ &= \langle \partial_\nu u_\epsilon, F \rangle = -F(x^*)|M| + O(\epsilon), \end{align*} $$

where $F$ is the solution to the boundary value problem (5.4) and $\langle \cdot, \cdot \rangle$ denotes the pairing between $H^{-1/2}(\partial M)$ and $H^{1/2}(\partial M)$. The last equality comes from (5.2), smoothness of $F$, and $\mathrm{supp}(\partial_\nu u_\epsilon) \subset \Gamma_{\epsilon,a}$. Inserting this into (5.28) we have

$$ \int_M u_\epsilon \mathrm{dvol}_g = \int_M F(x) \mathrm{dvol}_g + C_{\epsilon,a} |M| - F(x^*)|M| + O(\epsilon). $$

The constant $C_{\epsilon,a}$ is given by Proposition 5.2.

## 6. APPENDIX A - ELLIPTIC EQUATION FOR THE FIRST PASSAGE TIME

In this appendix we show that $u(x) := \mathbb{E}[\tau_\Gamma | X_0 = x]$ satisfies the boundary value problem (5.1). This is standard material but we could not find a suitable reference which precisely addresses our setting. As such we are including this appendix for the convenience of the reader.

Let $(M, g, \partial M)$ be an orientable compact connected Riemannian manifold with non-empty smooth boundary oriented by $\mathrm{dvol}_g$. Let also $(X_t, \mathbb{P}_x)$ be the Brownian
motion on $M$ starting at $x$, that is, the stochastic process generated by the Laplace-Beltrami operator $\Delta_g$. Let $\Gamma$ be a geodesic ball on $\partial M$ with radius $\varepsilon > 0$. We denote by $\tau_\Gamma$ the first time the Brownian motion $X_t$ hits $\Gamma$, that is

$$ \tau_{\Gamma} := \inf\{t \ge 0 : X_t \in \Gamma\}. $$

We set

$$ \mathcal{P}_{\Gamma}(t, x) := \mathbb{P}[\tau_{\Gamma} \le t | X_0 = x]. $$

Let us note that $\mathcal{P}_\Gamma(t, x)$ is the probability that the Brownian motion hits $\Gamma$ before or at time $t$, and therefore, satisfies

$$ (6.1) \qquad \mathcal{P}_{\Gamma}(0, x) = 0, \quad x \in M \setminus \Gamma, $$

$$ (6.2) \qquad \mathcal{P}_{\Gamma}(t, x) = 1, \quad (t, x) \in [0, \infty) \times \Gamma. $$

Note that, for any compact subset $\Gamma \subset M$, it follows$^4$

$$ \mathrm{Cap}(\Gamma, M) := \inf_{u \in C^\infty(\bar{M}), u|_{\Gamma}=1} \int_M \|\nabla u\|^2 \mathrm{dvol}_g = 0. $$

Then, [13, Theorem 1.5] implies that $(M, g, \partial M)$ is parabolic, that is, the probability that the Brownian motion ever hits any compact set $F$ with non-empty interior is 1. Since $\Gamma \subset \partial M$ is connected with non-empty interior on $\partial M$, we can extend $M$ to a compact connected Riemannian manifold $\tilde{M}$ such that $\tilde{M} \setminus M$ is compact with non-empty interior and $\tilde{M} \setminus M \cap M = \Gamma$. Note that, the Brownian motion, starting at any point $M \setminus \Gamma$, hits $\tilde{M} \setminus M$ if and only if it hits $\Gamma$. Therefore, the parabolicity condition of $(M, g)$ gives

$$ (6.3) \qquad \lim_{t \to \infty} \mathcal{P}_{\Gamma}(t, x) = 1, \quad x \in M. $$

Further, let us define the mean first arrival time $u$, as

$$ (6.4) \qquad u(x) := \mathbb{E}[\tau_{\Gamma} | X_0 = x] := \int_0^{\infty} t d\mathcal{P}_{\Gamma}(t, x), $$

where the integral is a Riemann-Stieltjes integral. To investigate $u$, let us recall some properties of $\mathcal{P}_\Gamma$. By Remmark 2.1 in [6], it follows that

$$ 1 - \mathcal{P}_{\Gamma}(t, x) = (e^{t\Delta_{\mathrm{mix}}} 1)(x), $$

where $e^{t\Delta_{\mathrm{mix}}}$ is the semigroup with infinitesimal generator $\Delta_{\mathrm{mix}}$, and $\Delta_{\mathrm{mix}}$ is the Laplace operator $\Delta_g$ corresponding to the Dirichlet boundary condition on $\Gamma$ and Neumann boundary condition on $\partial M \setminus \Gamma$, which is defined as follows

$$ (6.5) \qquad \mathrm{D}(\Delta_{\mathrm{mix}}) := \{u \in H^1(M) : \Delta_g u \in L^2 \mid u|_{\Gamma} = 0, \partial_{\nu} u|_{\Gamma^c} = 0\} $$

$$ (6.6) \qquad \Delta_{\mathrm{mix}} u = \Delta_g u \quad u \in \mathrm{D}(\Delta_{\mathrm{mix}}). $$

In (6.5) we define $\partial_{\nu} u \in H^{-1/2}(\partial M)$ using the same method for defining the Dirichlet to Neumann map. That is, for $u \in H^1(M)$ such that $\Delta_g u \in L^2(M)$, the distribution $\partial_{\nu} u|_{\partial M} \in H^{-1/2}(\partial M)$ acts on $f \in H^{1/2}(\partial M)$ via

$$ \langle \partial_{\nu} u |_{\partial M}, f \rangle := \int_M \Delta u_g \overline{v_f} \, dvol_g + \int_M g(du, dv_f) \, dvol_g, $$

<footnote>$^4$Note that in [6] and [12], the authors consider the manifold together with its boundary, and $C_c^\infty(M)$, $C_0^\infty(M)$ denote the set of smooth (up to the boundary) functions with compact support. In case of compact manifold, these sets coincide with $C^\infty(\bar{M})$.</footnote>
where $v_f \in H^1(M)$ is the harmonic extension of $f$. We say that $\partial_\nu u|_{\bar{\omega}} = 0$, for non-empty open set $\omega \subset \partial M$, if $\langle \partial_\nu u|_{\partial M}, f|_{\partial M} \rangle = 0$ for all $f \in H^{1/2}(\partial M)$ such that $f|_{\partial M \setminus \bar{\omega}} = 0$.

Note that if $u$ sufficiently regular, for instance $u \in H^2(M)$, then $\langle \partial_\nu u|_{\partial M}, f \rangle$ is equal to the boundary integral of $\partial_\nu u|_{\partial M}$ and $f$.

In fact, $\Delta_{mix}$ can be equivalently defined by quadratic form; see Proposition 7.1 in Appendix B. Moreover, $\Delta_{mix}$ is the non-positive self-adjoint operator with the discrete spectrum, consisting of negative eigenvalues accumulating at $-\infty$; see Proposition 7.1 in Appendix. Hence, $\Delta_{mix}$ satisfies the quadratic estimate

$$ \int_0^\infty \|t\Delta_{mix}(1 + t^2\Delta_{mix}^2)^{-1}u\|_{L^2}^2 \frac{dt}{t} \le C\|u\|^2, $$

for some $C > 0$ and all $u \in L^2(M)$; see for instance [19, p. 221]. Therefore, $\Delta_{mix}$ admits the functional calculus defined in [21].

**Remark 6.1.** The functional calculus in [21] is defined for a concrete operator, which is denoted by $T$ in the notation used in that article. However, $\Delta_{mix}$ satisfy all necessary conditions to admit this functional calculus.

Therefore, the semigroup $e^{t\Delta_{mix}}$, which is contracting by Hille-Yosida theorem [14, Theorem 8.2.3], can be defined as follows

$$ e^{t\Delta_{mix}}u = \frac{1}{2\pi i} \int_{\gamma_{a,\alpha}} e^{t\zeta}(\zeta - \Delta_{mix})^{-1}ud\zeta, \quad u \in L^2(M), $$

where $a \in (\tau, 0)$, $\alpha \in (0, \frac{\pi}{2})$, and $\gamma_{a,\alpha}$ is the anti-clockwise oriented curve:

$$ \gamma_{a,\alpha} := \{\zeta \in \mathbb{C} : \mathrm{Re}\zeta \le a, |\mathrm{Im}\zeta| = |\mathrm{Re}\zeta - a| \tan \alpha\}. $$

Let $\varepsilon > 0$ such that $a + \varepsilon < 0$. Then $\Delta_{mix} + \varepsilon$ is also a negative self-adjoint operator, and hence generates contracting semigroup, $e^{t(\Delta_{mix}+\varepsilon)}$, as above.

By definition, we obtain, for $u \in L^2(M)$,

$$ \begin{align} (6.7) \qquad & e^{t\Delta_{mix}}u = \frac{1}{2\pi i} \int_{\gamma_{a,\alpha}} e^{t\zeta}(\zeta - \Delta_{mix})^{-1}ud\zeta \nonumber \\ &= \frac{e^{-t\varepsilon}}{2\pi i} \int_{\gamma_{a,\alpha}} e^{t(\zeta+\varepsilon)}(\zeta + \varepsilon - (\Delta_{mix} + \varepsilon))^{-1}ud\zeta \nonumber \\ &= \frac{e^{-t\varepsilon}}{2\pi i} \int_{\gamma_{a+\varepsilon,\alpha}} e^{t\zeta}(\zeta - (\Delta_{mix} + \varepsilon))^{-1}ud\zeta = e^{-t\varepsilon}e^{t(\Delta_{mix}+\varepsilon)}u, \end{align} $$

where $\gamma_{a+\varepsilon,\alpha} = \gamma_{a,\alpha} + \varepsilon \subset \{\mathrm{Re}\xi < 0\}$. Let $f_1$ the constant function on $M$ equals 1. By Theorem 8.2.2 in [14], we know, for $\lambda > 0$,

$$ (\lambda - (\Delta_{mix} + \varepsilon))^{-1}f_1 = \int_0^\infty e^{-\lambda t}e^{t(\Delta_{mix} + \varepsilon)}f_1 dt. $$

Let us choose $\lambda = \varepsilon$, then, by using (6.7), we obtain

$$ -\Delta_{mix}^{-1}f_1 = \int_0^\infty e^{-\varepsilon t}e^{t(\Delta_{mix}+\varepsilon)}f_1 dt = \int_0^\infty e^{t\Delta_{mix}}f_1 dt $$

and hence,

$$ (6.8) \qquad \int_0^\infty (1 - \mathcal{P}_\Gamma(t, x))dt = -(\Delta_{mix}^{-1}f_1)(x) < \infty. $$
Therefore, the dominated convergence theorem implies

$$ \lim_{b \to \infty} \int_{0}^{b} (\mathcal{P}_{\Gamma}(b, x) - \mathcal{P}_{\Gamma}(t, x)) dt = \int_{0}^{\infty} (1 - \mathcal{P}_{\Gamma}(t, x)) dt < \infty. $$

Hence, by using (6.4) and integration by parts, we obtain

$$ \begin{aligned} u(x) &= \lim_{b \to \infty} \left( \mathcal{P}_{\Gamma}(b, x)b - \int_{0}^{b} \mathcal{P}_{\Gamma}(t, x)dt \right) = \lim_{b \to \infty} \int_{0}^{b} (\mathcal{P}_{\Gamma}(b, x) - \mathcal{P}_{\Gamma}(t, x)) dt < \infty \\ &= \int_{0}^{\infty} (1 - \mathcal{P}_{\Gamma}(t, x))dt. \end{aligned} $$

Therefore, by (6.8), we obtain

$$ \Delta_{mix}u = -f_1 = -1. $$

In particular, $u \in D(\Delta_{mix})$, and hence,

$$ u|_{\Gamma} = 0, \quad \partial_{\nu}u|_{\partial M \setminus \Gamma} = 0. $$

We see that (5.1) is satisfied.

## 7. APPENDIX B - QUADRATIC FORM

Let $(M, g, \partial M)$ be a compact connected Riemannian manifold with non-empty smooth boundary. Let $\Gamma$ be a closed subset of $\partial M$ such that $\partial M \setminus \Gamma^c$ is a non-empty open set. Consider the quadratic form

$$ (7.1) \quad a[u, v] := \int_M g(du, d\bar{v}) \mathrm{dvol}_g, \quad u, v \in D(a) := \{u \in H^1(M) : u|_{\Gamma} = 0\}. $$

Note that $D(a)$ is closed subspace of $H^1(M)$ containing $H_0^1(M)$ and $a[\cdot, \cdot]$ is a non-negative, closed, densely defined form. Therefore, by Friedrichs Theorem 2.23 in [15], it generates a non-negative self-adjoint operator $-\Delta_a$ in $L^2(M)$ whose domain is contained in $D(a)$ such that $(-\Delta_a u, u)_{L^2(M)} = a[u, u]$ for $u \in D(-\Delta_a)$.

Let us show that the resolvents of $-\Delta_a$ are compact. Assume that $s$ belongs to the resolvent set of $-\Delta_a$. Since $-\Delta_a(-\Delta_a - s)^{-1} : L^2(M) \to L^2(M)$ is bounded, it is suffices to show that $D(-\Delta_a)$, endowed with the graph norm, compactly embedded into $L^2(M)$. Since, for $u \in D(-\Delta_a)$,

$$ \|du\|_{L^2(M)}^2 = (du, du)_{L^2(M)} = a[u, u] = (-\Delta_a u, u)_{L^2(M)} \le \frac{1}{4}(\| -\Delta_a u \|_{L^2(M)} + \|u\|_{L^2(M)})^2, $$

we see that any bounded sequence in $D(-\Delta_a)$, endowed with the graph norm, is bounded in $H^1(M)$, and hence, it contains a Cauchy subsequence in $L^2(M)$ by Rellich-Kondrachov theorem. This implies that the resolvents of $-\Delta_a$ are compact, and hence, the spectrum of $-\Delta_a$ is discrete, consisting of non-negative eigenvalues accumulating at $+\infty$. Assume that $\lambda = 0$ is an eigenvalue, and let $u_0 \in D(a)$ be a corresponding eigenfunction.

Then the Poincaré-Wirtinger inequality gives, for some $C > 0$,

$$ \left\| u - \frac{1}{|M|} \int_M u \mathrm{dvol}_g \right\|_{L^2} \le C \|du\|_{L^2} = C(-\Delta_a u, u)_{L^2} = 0, $$

so that $u_0$ is a constant in $L_2(M)$. Since $u_0 \in H^1(M)$, we conclude that $u_0 = \text{const}$ in $L^2(\partial M)$, and hence, $u_0 = 0$ by choice of $\Gamma$. Therefore $\lambda = 0$ is not an eigenvalue, and consequently, the spectrum of $-\Delta_a$ is positive. For sake of completeness, we prove the following well known result.
**Proposition 7.1.** Let $-\Delta_a$ be the operator defined above and $-\Delta_{mix}$ be the operator defined in Section 6, then $-\Delta_a = -\Delta_{mix}$. In particular, $-\Delta_{mix}$ is a self-adjoint operator with the positive discrete spectrum accumulating at infinity.

*Proof.* Assume that $u, v \in H^1(M)$ and $\Delta_g u, \Delta_g u \in L^2(M)$. Let $V$ be the harmonic extension of $v|_{\partial M}$, then $\omega := V - v \in H_0^1(M)$ and $\Delta_g \omega \in L^2(M)$. By $H^2$-regularity, $\omega \in H^2(M)$. The generalized Green's identity gives

$$ - \int_M u \Delta_g \omega \, dvol_g + \int_M \Delta_g u \omega \, dvol_g = \langle \partial_\nu u |_{\partial M}, \omega |_{\partial M} \rangle - \langle u |_{\partial M}, \partial_\nu \omega |_{\partial M} \rangle, $$

where the first term of the right hand side vanishes since $\omega \in H_0^1(M)$. Therefore we get

$$ \begin{aligned} 0 &= - \int_M u \Delta_g \omega \, dvol_g + \langle u |_{\partial M}, \partial_\nu \omega |_{\partial M} \rangle + \int_M \Delta_g u \omega \, dvol_g \\ &= \int_M g(du, d\bar{\omega}) \, dvol_g + \int_M \Delta_g u \omega \, dvol_g. \end{aligned} $$

Hence, we obtain

$$ \langle \partial_{\nu} u |_{\partial M}, v |_{\partial M} \rangle = \int_M \Delta_g u \bar{v} \, dvol_g + \int_M g(du, d\bar{v}) \, dvol_g $$

for $u, v \in H^1(M)$ and $\Delta_g u, \Delta_g u \in L^2(M)$.

Assume that $u \in D(\Delta_{mix}) \subset D(a)$ and $v \in D(\Delta_a)$, then, by above formula,

$$ a[u, v] = \int_M g(du, d\bar{v}) \, dvol_g = - \int_M \Delta_g u \bar{v} \, dvol_g + \langle \partial_\nu u |_{\partial M}, v |_{\partial M} \rangle $$

Note that the last term vanishes since $u \in D(\Delta_{mix})$ and $v \in D(a)$, so that

$$ a[u, v] = - \int_M \Delta_g u \bar{v} \, dvol_g = - \int_M \Delta_{mix} u \bar{v} \, dvol_g. $$

Since this holds for all $v \in D(\Delta_a)$, it follows from Theorem 2.1., in [15], that $u \in D(\Delta_a)$ and $\Delta_a u = \Delta_{mix} u$.

Conversely, assume that $u \in D(\Delta_a)$, then $-\Delta_g u = -\Delta_a u \in L^2(M)$. Then, it follows

$$ \langle \partial_{\nu} u |_{\partial M}, f \rangle = \int_M \Delta u_g \overline{V_f} \, dvol_g + \int_M g(du, dV_f) \, dvol_g = a[u, V_f] - a[u, V_f] = 0. $$

for any $f \in H^{1/2}(\partial M)$ such that $f|_{\Gamma} = 0$. This means that $\partial_\nu u|_{\Gamma_c} = 0$, so that $u \in D(\Delta_{mix})$ and $\Delta_a u = \Delta_{mix} u$.

□

## REFERENCES

[1] Habib Ammari, Kostis Kalimeris, Hyeonbae Kang, and Hyundae Lee. Layer potential techniques for the narrow escape problem. *J. Math. Pures Appl.* (9), 97(1):66–84, 2012.

[2] Xinfu Chen and Avner Friedman. Asymptotic analysis for the narrow escape problem. *SIAM J. Math. Anal.*, 43(6):2542–2563, 2011.

[3] Alexei F Cheviakov, Michael J Ward, and Ronny Straube. An asymptotic analysis of the mean first passage time for narrow escape problems: Part II: The sphere. *Multiscale Modeling & Simulation*, 8(3):836–870, 2010.

[4] F. G. Friedlander. *Introduction to the theory of distributions*. Cambridge University Press, Cambridge, 1982.
[5] Daniel Gomez and Alexei F. Cheviakov. Asymptotic analysis of narrow escape problems in nonspherical three-dimensional domains. *Phys. Rev. E*, 91:012137, Jan 2015.

[6] Alexander Grigor'yan and Laurent Saloff-Coste. Hitting probabilities for Brownian motion on Riemannian manifolds. *J. Math. Pures Appl. (9)*, 81(2):115-142, 2002.

[7] R. Hiptmair, C. Jerez-Hanckes, and C. Urzúa-Torres. Closed-form inverses of the weakly singular and hypersingular operators on disks. *Integral Equations Operator Theory*, 90(1):Art. 4, 14, 2018.

[8] D Holcman and Z Schuss. Escape through a small opening: receptor trafficking in a synaptic membrane. *Journal of Statistical Physics*, 117(5-6):975-1014, 2004.

[9] D. Holcman and Z. Schuss. Escape through a small opening: receptor trafficking in a synaptic membrane. *J. Statist. Phys.*, 117(5-6):975-1014, 2004.

[10] David Holcman and Zeev Schuss. Stochastic narrow escape in molecular and cellular biology. *Analysis and Applications*. Springer, New York, 2015.

[11] Lars Hörmander. Fourier Integral Operators. Acta Math., 127 (1971), 79-183

[12] Joonas Ilmavirta and Gabriel P. Paternain. Functions of constant geodesic x-ray transform. *Inverse Problems*, 35(6):065002, 11, 2019.

[13] Debora Impera, Stefano Pigola, and Alberto G. Setti. Potential theory for manifolds with boundary and applications to controlled mean curvature graphs. *J. Reine Angew. Math.*, 733:121-159, 2017.

[14] Jürgen Jost. *Partial differential equations*, volume 214 of *Graduate Texts in Mathematics*. Springer, New York, third edition, 2013.

[15] Tosio Kato. *Perturbation theory for linear operators*. Classics in Mathematics. Springer-Verlag, Berlin, 1995. Reprint of the 1980 edition.

[16] John M. Lee. *Introduction to Riemannian manifolds*, volume 176 of *Graduate Texts in Mathematics*. Springer, Cham, 2018. Second edition of [MR1468735].

[17] John M. Lee and Gunther Uhlmann. Determining anisotropic real-analytic conductivities by boundary measurements. *Comm. Pure Appl. Math.*, 42(8):1097-1112, 1989.

[18] P. A. Martin. Exact solution of some integral equations over a circular disc. *J. Integral Equations Appl.*, 18(1):39-58, 2006.

[19] Alan McIntosh. Operators which have an $H_\infty$ functional calculus. In *Miniconference on operator theory and partial differential equations (North Ryde, 1986)*, volume 14 of Proc. Centre Math. Anal. Austral. Nat. Univ., pages 210-231. Austral. Nat. Univ., Canberra, 1986.

[20] François Monard, Richard Nickl, and Gabriel P. Paternain. Efficient nonparametric Bayesian inference for X-ray transforms. *Ann. Statist.*, 47(2):1113-1147, 2019.

[21] Medet Nursultanov and Andreas Rosén. Evolution of time-harmonic electromagnetic and acoustic waves along waveguides. *Integral Equations Operator Theory*, 90(5):Paper No. 53, 32, 2018.

[22] Leonid Pestov and Gunther Uhlmann. On characterization of the range and inversion formulas for the geodesic X-ray transform. *Int. Math. Res. Not.*, (80):4331-4347, 2004.

[23] Leonid Pestov and Gunther Uhlmann. Two dimensional compact simple Riemannian manifolds are boundary distance rigid. *Ann. of Math. (2)*, 161(2):1093-1110, 2005.

[24] S Pillay, Michael J Ward, A Peirce, and Theodore Kolokolnikov. An asymptotic analysis of the mean first passage time for narrow escape problems: Part I: Two-dimensional domains. *Multiscale Modeling & Simulation*, 8(3):803-835, 2010.

[25] Xavier Saint Raymond *Elementary Introduction to the Theory of Pseudodifferential Operators*. Studies in Advanced Mathematics, 1991 CRC Press.

[26] Zeev Schuss. *Theory and applications of stochastic differential equations*. John Wiley & Sons, Inc., New York, 1980. Wiley Series in Probability and Statistics.

[27] Zeev Schuss, Amit Singer, and David Holcman. Narrow escape, part i. *Journal of Statistical Physics*, 122(41):437-463, 2006.

[28] Zeev Schuss, Amit Singer, and David Holcman. The narrow escape problem for diffusion in cellular microdomains. *Proceedings of the National Academy of Sciences*, 104(41):16098-16103, 2007.

[29] V. A. Sharafutdinov. *Integral geometry of tensor fields*. Inverse and Ill-posed Problems Series. VSP, Utrecht, 1994.

[30] Amit Singer, Z Schuss, and D Holcman. Narrow escape and leakage of Brownian particles. *Physical Review E*, 78(5):051111, 2008.
[31] Amit Singer, Zeev Schuss, and David Holcman. Narrow escape, part II: The circular disk. *Journal of statistical physics*, 122(3):465-489, 2006.

[32] Michael E. Taylor. *Pseudodifferential Operators*, volume 34 of PRINCETON MATHEMATICAL SERIES. Princeton University Press, first edition, 1981.

[33] Michael E. Taylor. *Partial differential equations I. Basic theory*, volume 115 of *Applied Mathematical Sciences*. Springer, New York, second edition, 2011.

[34] Michael E. Taylor. *Partial differential equations II. Qualitative studies of linear equations*, volume 116 of *Applied Mathematical Sciences*. Springer, New York, second edition, 2011.

[35] Lefeuvre Thibault. Sur la rigidité des variétés riemanniennes. 2019.

SCHOOL OF MATHEMATICS AND STATISTICS, UNIVERSITY OF SYDNEY
*Email address: medet.nursultanov@gmail.com*

SCHOOL OF MATHEMATICS AND STATISTICS, UNIVERSITY OF SYDNEY
*Email address: leo.tzou@gmail.com*

DEPARTMENT OF MATHEMATICS AND STATISTICS, MACQUARIE UNIVERSITY
*Email address: tzou.justin@gmail.com*