# Stability of eigenfunctions

In [CZ1994] some bounds for the Neumann heat kernels $P_t(x,y)$ on a general domain are given; in particular, one can bound this kernel by a multiple of the Euclidean heat kernel for t small enough. Using this bound, [BP2008] showed uniform stability of the second eigenfunction (in the uniform topology) and eigenvalue with respect to uniform perturbations of the domain.

## Formal theory

Suppose one has a one-parameter family $t \mapsto L(t)$ of self-adjoint operators (on some Hilbert space, e.g. $L^2(\Omega)$, though for this formal computation the domain will not be important), and one-parameter families $t \mapsto u(t)$, $t \mapsto \lambda(t)$ to the eigenfunction equation

$L(t) u(t) = \lambda(t) u(t)$. (1)

We normalise $u(t)$ to have norm 1:

$\langle u(t), u(t) \rangle = 1$. (2)

Formally, if we differentiate the norm equation (2) at time zero, we get

$\langle \dot u(0), u(0) \rangle = 0$ (3)

while if we differentiate (1) at time zero we obtain

$\dot L(0) u(0) + L(0) \dot u(0) = \dot \lambda(0) u(0) + \lambda(0) \dot u(0)$. (4)

Taking the inner product of (4) with u(0) and using (2), (3) we conclude the Hadamard first variation formula for eigenvalues:

$\dot \lambda(0) = \langle \dot L(0) u(0), u(0) \rangle$. (5)

If we instead take the orthogonal projection $\pi_{u(0)}^\perp$ onto the orthogonal complement of $u(0)$, we obtain the Hadamard first variation formula for eigenfunctions:

$(L(0)-\lambda) \dot u(0) = - \pi_{u(0)}^\perp( \dot L(0) u(0) )$. (6)

Formally, if $\lambda$ is a simple eigenvalue, then $L(0)-\lambda$ is invertible on the orthogonal complement of $u(0)$, and so (6) and (3) allow one to solve for $\dot u(0)$.

## L^2 and H^1 theory

Suppose that one has a domain $\Omega$ with second Neumann eigenvalue $\lambda_2$, and third Neumann eigenvalue (not counting multiplicity) $\lambda_3 \gt \lambda_2$. Thus, one has

$\int_\Omega |\nabla u|^2 \geq \lambda_2 \int_\Omega |u|^2$ (7)

for all mean zero u, with equality when u lies in the second eigenspace $V_2$, and one can improve this to

$\int_\Omega |\nabla u|^2 \geq \lambda_3 \int_\Omega |u|^2$ (8)

when u is orthogonal to $V_2$.

Now consider a perturbation $B\Omega$ of $\Omega$, where B is an invertible linear transformation. Then a second eigenfunction of $B\Omega$, after change of variables, becomes a function u on $\Omega$ that minimizes the modified Rayleigh quotient

$\int_\Omega M \nabla u \cdot \nabla u / \int_\Omega |u|^2$

where $M := (B^{-1}) (B^{-1})^T$. We may normalize this eigenfunction as $u = u_2 + v$, where u_2 is a unit eigenfunction in V_2 and v is orthogonal to V_2, so that $\|u\|_2^2 = 1 + \|v\|_2^2$. Then the modified Rayleigh quotient of u is less than or equal to that of u_2, and hence

$\int_\Omega M \nabla u \cdot \nabla u \leq (\int_\Omega M \nabla u_2 \cdot \nabla u_2) (1 + \|v\|_2^2 ).$

Expanding out u as u_2+v and rearranging, we end up at

$\int_\Omega M \nabla v \cdot \nabla v \leq \|v\|_2^2 \int_\Omega M \nabla u_2 \cdot \nabla u_2 - 2 \int_\Omega M \nabla u_2 \cdot \nabla v.$

Note that $\nabla v$ is orthogonal to $\nabla V_2$ by integration by parts, and so we may replace $M \nabla u_2$ on the RHS by the orthogonal projection $\pi_{\nabla V_2}^\perp(M \nabla u_2)$. Letting $\sigma_1(M) = \sigma_2(B)^{-2}$ be the least singular value of M, we now apply Cauchy-Schwarz and conclude that

$\sigma_1(M) \|\nabla v \|_2^2 \leq \|v\|_2^2 \int_\Omega M \nabla u_2 \cdot \nabla u_2 + 2 \| \pi_{\nabla V_2}^\perp(M \nabla u_2) \|_2 \|\nabla v\|_2$.

By (8) we may bound $\|v\|_2^2 \leq \frac{1}{\lambda_3} \|\nabla v\|_2^2$, and so we conclude that

$(\sigma_1(M) - \frac{\int_\Omega M \nabla u_2 \cdot u_2}{\lambda_3}) \|\nabla v \|_2 \leq 2 \| \pi_{\nabla V_2}^\perp(M \nabla u_2) \|_2$. (9)

This gives an H^1 bound on the error v between the perturbed eigenfunction u and the original eigenfunction u_2. To understand this bound, note that we may upper bound

$\int_\Omega M \nabla u_2 \cdot \nabla u_2 \leq \sigma_2(M) \int_\Omega |\nabla u_2|^2 = \lambda_2 \sigma_2(M)$

and also $\pi_{\nabla V_2}^\perp(M \nabla u_2) = \pi_{\nabla V_2^\perp}((M-\frac{\sigma_1(M)+\sigma_2(M)}{2}) \nabla u_2)$ so that

$\|\pi_{\nabla V_2^\perp}((M-\frac{\sigma_1(M)+\sigma_2(M)}{2}) \nabla u_2)\|_2 \leq \frac{\sigma_2(M)-\sigma_1(M)}{2} \|\nabla u_2\|_2 = \frac{\sigma_2(M)-\sigma_1(M)}{2} \lambda_2^{1/2}$

and so one has

$\|\nabla v \|_2 \leq \frac{(\sigma_2(M)-\sigma_1(M)) \lambda_2^{1/2}}{\sigma_1(M) - \sigma_2(M) \frac{\lambda_2}{\lambda_3}}$

provided that the denominator is positive. In terms of the condition number $\kappa := \sigma_2(B)/\sigma_1(B)$ of B, this becomes

$\|\nabla v \|_2 \leq \frac{(\kappa^2-1) \lambda_2^{1/2}}{1 - \kappa^2 \frac{\lambda_2}{\lambda_3}}$ (15)

which is a non-trivial bound when $\kappa \lt (\lambda_3/\lambda_2)^{1/2}$, and is of the order of $O(\lambda_2^{1/2} (\kappa-1))$ when $\kappa$ is close to 1. Using (8), we conclude in particular that

$\|v \|_2 \leq \frac{(\kappa^2-1) (\lambda_2/\lambda_3)^{1/2}}{1 - \kappa^2 \frac{\lambda_2}{\lambda_3}}$ (16)

For these calculations performed on a third reference triangle, http://www.math.sfu.ca/~nigam/polymath-figures/Perturbation.pdf

## An alternate approach to the L^2 and H^1 theory

Let us keep the reference triangle $\Omega$ fixed, and view the matrix M=M(t) as being smoothly time dependent, so that the eigenvalues $\lambda_k = \lambda_k(t)$ and L^2-normalised eigenfunctions $u_k = u_k(t)$ are also time dependent. Assume for the sake of argument that eigenvalues stay simple and all functions depend smoothly on t. We have the eigenfunction equation

$-\nabla \cdot M \nabla u_k = \lambda_k u_k$

and the Neumann boundary condition

$-n \cdot M \nabla u_k = 0$

and the L^2 normalisation

$\langle u_k, u_k \rangle = 1$.

Differentiating, we conclude that

$(-\nabla \cdot M \nabla - \lambda_k) \dot u_k = \nabla \cdot \dot M \nabla u_k + \dot \lambda_k u_k$ (10)

and

$-n \cdot M \nabla \dot u_k = n \cdot \dot M \nabla u_k$ (11)

and

$\langle \dot u_k, u_k \rangle = 0$. (12)

Taking the inner product of (10) with u_k and using (12) yields

$\langle -\nabla \cdot M \nabla \dot u_k, u_k \rangle = \langle \nabla \cdot \dot M \nabla u_k, u_k \rangle + \dot \lambda_k$.

Integrating by parts using the Neumann condition and (11) yields

$\langle -\nabla \cdot M \nabla \dot u_k, u_k \rangle = \langle \dot u_k, -\nabla \cdot M \nabla u_k \rangle + \int_{\partial \Omega} (n \cdot \dot M \nabla u_k) u_k$.

By the eigenfunction equation and (12), the first inner product on the RHS vanishes. By Stokes theorem one has

$\int_{\partial \Omega} (n \cdot \dot M \nabla u_k) u_k = \langle \nabla \cdot \dot M \nabla u_k, u_k \rangle + \langle \dot M \nabla u_k, \nabla u_k \rangle$

and thus we have the variation formula

$\dot \lambda_k = \langle \dot M \nabla u_k, \nabla u_k \rangle$. (14)

Next, if we take the inner product of (10) with u_l for some l distinct from k, one has

$\langle (-\nabla \cdot M \nabla - \lambda_k) \dot u_k, u_l \rangle = \langle \nabla \cdot \dot M \nabla u_k, u_l \rangle$.

Integrating by parts as before, we have

$\langle -\nabla \cdot M \nabla \dot u_k, u_l \rangle = \langle \dot u_k, -\nabla \cdot M \nabla u_l \rangle + \int_{\partial \Omega} (n \cdot \dot M \nabla u_k) u_l$.

By the eigenfunction equation, the first inner product on the RHS is $\lambda_l \langle \dot u_k, u_l \rangle$. By Stokes theorem we have

$\int_{\partial \Omega} (n \cdot \dot M \nabla u_k) u_l = \langle \nabla \cdot \dot M \nabla u_k, u_l \rangle + \langle \dot M \nabla u_k, \nabla u_l \rangle$

and thus

$(\lambda_l - \lambda_k) \langle \dot u_k, u_l \rangle = \langle \dot M \nabla u_k, \nabla u_l \rangle$

and thus by eigenfunction expansion (and (12))

$\dot u_k = \sum_{l \neq k} \frac{1}{\lambda_l-\lambda_k} \langle \dot M \nabla u_k, \nabla u_l \rangle u_l$

where the convergence is in an unconditional L^2 sense. (Note that $\nabla u_l / \lambda_l^{1/2}$ is an orthonormal system and so from Bessel's inequality we know that $\sum_l \langle \dot M \nabla u_k, \nabla u_l \rangle^2 / \lambda_l \lt \infty$, which is enough decay to justify the L^2 convergence of the RHS.) In fact we may differentiate and conclude that

$\dot \nabla u_k = \sum_{l \neq k} \frac{1}{\lambda_l-\lambda_k} \langle \dot M \nabla u_k, \nabla u_l \rangle \nabla u_l$

where the convergence is again in the unconditional L^2 sense (i.e. the previous convergence was in the unconditional H^1 sense). From the Bessel inequality we see in particular that

$\| \dot \nabla u_k \|_{L^2}^2 = \sum_{l \neq k} \frac{\lambda_l^2}{(\lambda_l-\lambda_k)^2} |\langle \dot M \nabla u_k, \nabla u_l \rangle|/\lambda_l$
$\leq \| \dot M \nabla u_k \|_{L^2}^2 \sup_{l \neq k} \frac{\lambda_l^2}{(\lambda_l - \lambda_k)^2}$;

in particular, we have

$\| \dot \nabla u_2 \|_{L^2}^2 \leq \frac{\lambda_3}{\lambda_3 - \lambda_2} \| \dot M \nabla u_2 \|_{L^2}$

and thus

$\| \dot u_2 \|_{L^2}^2 \leq \frac{1}{\lambda_3 - \lambda_2} \| \dot M \nabla u_2 \|_{L^2}$

which can be viewed as infinitesimal variants of the estimates in the previous section.

## A Sobolev inequality

Lemma Suppose that $u: B(0,R) \to {\mathbf R}$ obeys the inhomogeneous eigenfunction equation $-\Delta u = \lambda u + F$. Suppose also that $\eta: {\mathbf R} \to {\mathbf R}$ is a smooth function that equals 1 at 0 and vanishes on $[R,+\infty)$. Then
$|u(0)| \leq \frac{\sqrt{\pi}}{2\sqrt{2}} (\|F\|_{L^2(B(0,R)} (\int_0^\infty \eta(r)^2 Y_0(\sqrt{\lambda} r)^2\ r dr)^{1/2} + \lambda^{1/2} \|u\|_{L^2(B(0,R))} (\int_0^R \eta'(r)^2 Y'_0(\sqrt{\lambda}r)^2\ r dr)^{1/2}$
$+ \|\nabla u\|_{L^2(B(0,R))} (\int_0^R \eta'(r)^2 Y_0(\sqrt{\lambda}r)^2\ r dr)^{1/2} ).$

Proof We introduce the averaged functions

$\bar u(r) := \frac{1}{2\pi} \int_0^{2\pi} u(r,\theta)\ d\theta$
$\bar F(r) := \frac{1}{2\pi} \int_0^{2\pi} F(r,\theta)\ d\theta$

and observe that $\bar u$ obeys the inhomogeneous Bessel equation

$-\bar u''(r) - \frac{1}{r} \bar u'(r) = \lambda \bar u(r) + \bar F(r)$

with initial condition $\bar u(0) = u(0)$.

We can solve this equation by separation of variables, writing

$\bar{u}(r) = a(r) J_0(\sqrt{\lambda} r) + b(r) Y_0(\sqrt{\lambda} r)$
$\bar{u}'(r) = \sqrt{\lambda} a(r) J_0'(\sqrt{\lambda} r) + \sqrt{\lambda} b(r) Y_0'(\sqrt{\lambda} r)$.

Using the standard Wronskian identity

$J_0(\sqrt{\lambda} r) Y'_0(\sqrt{\lambda} r) - Y_0(\sqrt{\lambda} r) J'_0(\sqrt{\lambda} r) = \frac{2}{\sqrt{\lambda} \pi r}$

one soon arrives at the equations

$a(r) = \frac{\pi r}{2} ( \sqrt{\lambda} \bar u(r) Y'_0(\sqrt{\lambda} r) - \bar u'(r) Y_0(\sqrt{\lambda} r) )$

and

$a'(r) = -\frac{\pi r}{2} \bar F(r) Y_0(\sqrt{\lambda} r)$.

Meanwhile, by integration by parts

$a(0) = -\int_0^R a'(r) \eta(r)\ dr - \int_0^R a(r) \eta'(r)\ dr$

and thus

$|a(0)| \leq \frac{\pi}{2} ( \int_0^R |\bar F(r)| |Y_0(\sqrt{\lambda} r)|\ r dr + \int_0^R \sqrt{\lambda} |\bar u(r)| |Y'_0(\sqrt{\lambda} r)|\ r dr + \int_0^R |\bar u'(r)| |Y_0(\sqrt{\lambda} r)|\ r dr).$

From Cauchy-Schwarz one has

$2\pi \int_0^R |\bar u(r)|^2\ r dr \leq \| u \|_{L^2(B(0,R))}^2$

and

$2\pi \int_0^R |\bar u'(r)|^2\ r dr \leq \| \nabla u \|_{L^2(B(0,R))}^2$

and similarly

$2\pi \int_0^R |\bar F(r)|^2\ r dr \leq \| F \|_{L^2(B(0,R))}^2.$

The claim then follows from further applications of the Cauchy-Schwarz inequality. []

In principle, the above inequality, when combined with the previous L^2 and H^1 estimates, gives L^infty control on $\dot u$ in the interior of the triangle. One should also be able to get control near the interior of an edge by reflection. One needs to modify the argument a bit though to handle what is going on at vertices.

Here is an alternate approach that avoids the use of Bessel functions of the first and second kinds. Suppose one has a solution to the inhomogeneous eigenfunction equation $-\Delta u = \lambda u + F$ on a ball B(0,R). Let $\eta$ be a smooth (or at least C^2) cutoff on this ball which equals 1 at the origin. Then $u\eta$ is compactly supported with Laplacian

$-\Delta(\eta u) = \lambda \eta u + \eta F - 2 \nabla \eta \cdot \nabla u -\Delta \eta u$

and hence by the fundamental solution to the Laplacian

$u(0) = \eta u(0) = \frac{1}{2\pi} \int_{B(0,R)} \log |x| (\lambda \eta u + \eta F - 2 \nabla \eta \cdot \nabla u -\Delta \eta u)(x)\ dx$

and hence by Cauchy-Schwarz

$|u(0)| \leq \frac{1}{2\pi} (\|u\|_{L^2(B(0,R))} \| \log|x| (\lambda \eta - \Delta \eta) \|_{L^2(B(0,R)} + 2 \|\nabla u \|_{L^2(B(0,R)} \| \nabla \eta \|_{L^2(B(0,R))} + \|F\|_{L^2(B(0,R)} \| \log|x| \eta \|_{L^2(B(0,R))}.$

Because $\Delta(\eta u)$ has mean zero, one can also replace $\log |x|$ here by $\log |x|-C$ for any constant C, for instance one can use $\log(|x|/R)$.

## Sobolev inequalities in logarithmic coordinates

The computations appear to become cleaner if we work in logarithmic coordinates

$(s,\theta) := (\log r, \theta)$.

Thus, for instance, a sector $\{ (r,\theta): 0 \lt r \lt R; 0 \leq \theta \leq \alpha \}$ in polar coordinates, becomes a half-infinite strip $\{ (s,\theta): s \lt \log R; 0 \leq \theta \leq \alpha \}$ in logarithmic coordinates. A triangle $\Omega$ which contains this sector as one of its corners, then becomes a slightly enlarged version $\log \Omega$ of this half-infinite strip, with a concave boundary connecting the two half-infinite edges.

Consider an eigenfunction $-\Delta u = \lambda u$ in the original domain $\Omega$. In polar coordinates, the eigenfunction equation is

$-\partial_{rr} u - \frac{1}{r} \partial_r u - \frac{1}{r^2} \partial_{\theta\theta} u = \lambda u$

the Rayleigh quotient is

$\int_\Omega (\partial_r u)^2 + \frac{1}{r^2} (\partial_\theta u)^2\ r dr d\theta / \int_\Omega u^2\ r dr d\theta$

and the Neumann boundary condition on the two radial edges are

$\partial_\theta u = 0$.

If we convert to logarithmic coordinates $v: \log \Omega \to \R$ by the change of variables

$v(s,\theta) := u(e^s, \theta)$

then one can compute that the eigenfunction equation now becomes

$- \partial_{ss} v - \partial_{\theta\theta} v = \lambda e^{2s} v$ (A.1)

the Rayleigh quotient becomes

$\int_{\log \Omega} (\partial_s v)^2 + (\partial_\theta v)^2\ ds d\theta / \int_{\log \Omega} |v|^2 e^{2s}\ ds d\theta$

and the Neumann boundary conditions on the two half-infinite edges are

$\partial_\theta v = 0$. (A.2)

We can reflect $\log \Omega$ infinitely often across the half-infinite edges, extending v to a half-infinite domain that contains the half-space $\{ (s,\theta): s \lt \log R, \theta \in \R \}$. Thanks to the Neumann condition (A.2), v will remain smooth and obey the eigenfunction equation (A.1) in the entirety of this half-space.

Let's suppose that we have normalised our original eigenfunction u to have an L^2 norm of 1, then after transforming into logarithmic coordinates we have

$\int_{\log \Omega} |v|^2 e^{2s}\ ds = 1$ (A.3)

and hence by the Rayleigh quotient we have H^1 control:

$\int_{\log \Omega} (\partial_s v)^2 + (\partial_\theta v)^2\ ds d\theta = \lambda$. (A.4)

Now we can start moving towards L^infty control of v. We first look at the mean

$\bar{v}(s) := \frac{1}{\alpha} \int_0^\alpha v(s,\theta)\ d\theta$.

From the eigenfunction equation and the Neumann boundary condition, the mean obeys the equation

$\partial_{ss} \bar{v} = \lambda e^{2s} \bar{v}$. (A.5)

On the other hand, the L^2 normalisation of v and Cauchy-Schwarz gives

$\int_{-\infty}^{\log R} |\bar v|^2 e^{2s}\ ds \leq 1/\alpha$ (A.6)

and the H^1 bounds similarly gives

$\int_{-\infty}^{\log R} |\partial_s \bar v|^2\ ds \leq \lambda/\alpha$. (A.7)

Among other things, (A.7) tells us that $\partial_s \bar v$ must go to zero along some subsequence as $s \to -\infty$. (Indeed, a Bessel expansion of the original eigenfunction shows that it goes to zero uniformly, but we won't need this fact here.) Integrating (A.5) and using (A.6) and Cauchy-Schwarz then gives us the bound

$|\partial_s \bar{v}(s)| \leq 2^{-1/2} \alpha^{-1/2} \lambda e^s$ (A.8)

and hence on a second integration

$|\bar{v}(s)-\bar{v}(-\infty)| \leq 2^{-1/2} \alpha^{-1/2} \lambda e^s$. (A.9)

From the triangle inequality, we conclude that

$(\int_{-\infty}^{\log R} |\bar v|^2 e^{2s}\ ds)^{1/2} \geq |\bar v(-\infty)| 2^{-1/2} R - 2^{-3/2} \alpha^{-1/2} \lambda R^2$ (A.10)

and hence by (A.6)

$|\bar v(-\infty)| \leq 2^{1/2} \alpha R^{-1} + 2^{-1} \alpha^{-1/2} \lambda R$ (A.11)

(which already bounds the original eigenfunction at the vertex) and hence by (A.9)

$|\bar v(s)| \leq 2^{1/2} \alpha R^{-1} + 2^{-1} \alpha^{-1/2} \lambda R + 2^{-1/2} \alpha^{-1/2} \lambda e^s$. (A.12)

This gives L^infty bounds on the mean $\bar v$. Now we move on to bounding v itself. We have the Sobolev inequality

$\sup_\theta |v(s,\theta)-\bar v(s)| \leq \int_0^\alpha |\partial_\theta v(s,\theta)|\ d\theta \leq \alpha^{1/2} (\int_0^\alpha |\partial_\theta v(s,\theta)|^2\ d\theta)^{1/2}$ (A.13)

so it will suffice to get pointwise bounds on the angular energy

$E(s) := (\int_0^\alpha |\partial_\theta v(s,\theta)|^2\ d\theta)^{1/2}$ (A.14)

since we then have

$|v(s,\theta)| \leq |\bar v(s)| + \alpha^{1/2} E(s)$ (A.15)

from the triangle inequality.

Note that (A.4) gives an integrated bound

$\int_{-\infty}^{\log R} E(s)^2\ ds \leq \lambda$. (A.16)

On the other hand, differentiating (A.1) in $\theta$ gives

$- \Delta \partial_\theta v = \lambda e^{2s} \partial_\theta v$

and hence by integration by parts

$\int_{\log \Omega} |\nabla \partial_\theta v|^2 \eta^2 = \int_{\log \Omega} \lambda e^{2s} |\partial_\theta v|^2 \eta^2 - 2 \int_{\log \Omega} \partial_\theta v \nabla \eta \cdot \nabla \partial_\theta v \eta$

for any test function $\eta$. Bounding

$-2 \int_{\log \Omega} \partial_\theta v \nabla \eta \cdot \nabla \partial_\theta v \eta \leq \frac{1}{2} \int_{\log \Omega} |\nabla \partial_\theta v|^2 \eta^2 + 2 \int_{\log \Omega} |\partial_\theta v|^2 |\nabla \eta|^2$

we conclude that

$\int_{\log \Omega} |\nabla \partial_\theta v|^2 \eta^2 \leq \int_{\log \Omega} |\partial_\theta v|^2 (2 \eta^2 \lambda e^{2s} + 4 |\nabla \eta|^2)$.

In particular, if $\eta$ is a function just of s and is supported in $(-\infty,\log R]$, we conclude that

$\int_{-\infty}^{\log R} |\partial_s E(s)|^2 \eta^2(s)\ ds \leq \int_{-\infty}^{\log R} E(s)^2 (2 \eta^2(s) \lambda e^{2s} + 4 \eta'(s)^2)\ ds.$. (A.17)

Combining this with (A.16) for various choices of cutoff $\eta$ gives integral control on $\partial_s E(s)$, and this combined with (A.16) gives pointwise control on E(s) and hence on v(s). This gives L^infty control on the eigenfunction u near the vertex. (The computations are somewhat ugly, though, particularly if one wants to optimise in eta.)

## H^2 theory

Formally, if on a triangle $\Omega$, a function v obeys the inhomogeneous Laplace equation

$-\Delta v = F$ (B.1)

with homogeneous Neumann boundary data $n \cdot \nabla v = 0$, then after two integrations by parts we have the H^2 bound

$\int_\Omega |\nabla^2 v|^2 = \int_\Omega |F|^2$

(Bochner-Weitzenbock identity.) The integrations by parts can be justified if v is smooth on Omega, with first derivative uniformly bounded and second derivative blowing up by at most $O(r^{-1+\varepsilon})$ at each vertex for some $\varepsilon\gt0$, by the usual trick of inserting a cutoff function to smoothly cut out the integral at the vertices, performing the integration by parts, and then sending the cutoff parameter to zero. I think this sort of regularity for acute-angled triangles can be achieved in practice without difficulty.

Now suppose that one is solving the inhomogeneous Laplace equation (B.1) but with inhomogeneous Neumann data

$n \cdot \nabla v = f$. (B.2)

Then the integration by parts becomes much less favorable, as there are now a number of nasty boundary terms to take care of. But suppose that one can find a reasonably smooth function w (of the same regularity needed for v, i.e. first derivative uniformly bounded and second derivative blowing up slower than 1/r) which solves (B.2) but not (B.1), i.e.

$n \cdot \nabla w = f$

on the boundary. Then the difference v-w solves (B.1) with forcing term $F + \Delta w$ and thus

$\| \nabla^2(v-w) \|_{L^2(\Omega)} = \| F + \Delta w \|_{L^2(\Omega)}.$

From the pointwise estimate $|\Delta w| \leq \sqrt{2} |\nabla^2 w|$ and the triangle inequality we conclude that

$\| \nabla^2 v \|_{L^2(\Omega)} \leq \| F \|_{L^2(\Omega)} + (1+\sqrt{2}) \| w \|_{L^2(\Omega)}.$ (B.3)

We can apply this to the first variation $\dot u = \dot u_2$ of the second eigenfunction at M=I. From the L^2 and H^1 stability theory one has

$-\Delta \dot u = \lambda_2 \dot u + \nabla \cdot \dot M \nabla u + \dot \lambda_2 u$

with boundary condition

$n \cdot \nabla \dot u = - n \cdot \dot M \nabla u$ (B.4)

and thus we have the H^2 bound

$\| \nabla^2 \dot u \|_{L^2(\Omega)} \leq \lambda_2 \| \dot u \|_{L^2(\Omega)} + \| \nabla \cdot \dot M \nabla u \|_{L^2(\Omega)} + |\dot \lambda_2| \|u \|_{L^2(\Omega)} + (1+\sqrt{2}) \|w\|_{L^2(\Omega)}$ (B.5)

where w is any function with boundary value (B.4). From the L^2 stability theory one has

$\| \dot u \|_{L^2(\Omega)}\leq \frac{1}{\lambda_3 - \lambda_2} \| \dot M \nabla u \|_{L^2}$

and so one has bounds on everything on the RHS of (B.5) except for the last term.

To control this term, we divide the triangle $\Omega = ABC$ into three overlapping regions $\Omega_A, \Omega_B, \Omega_C$, each one being a region around one of the three vertices A,B,C that avoids the other two vertices, as well as the edge joining those vertices. In each of these regions, we find a function w that solves the boundary condition (B.5) in that region, which obeys good H^2 bounds; combining these functions together by a partition of unity will then give a global w in $H^2(\Omega)$ that obeys (B.5) on the entire boundary.

Let's look at the region $\Omega_A$ near the vertex A, which we set to be the origin. In polar coordinates $(r,\theta)$, we know that $\partial_\theta u = 0$ on the edges AB, AC, and so the boundary condition (B.4) becomes

$n \cdot \nabla \dot u = c_{AB} \partial_r u$

on AB and

$n \cdot \nabla \dot u = c_{AC} \partial_r u$

on AC for some explicit scalar constants $c_{AB}, c_{AC}$ depending only on geometry of the triangle ABC and on $\dot M$. On $\Omega_A$, we can build a function w_A that exhibits this data by the formula

$w_A(x) := u(\Phi(x))$

where $\Phi: \Omega_A \to \Omega$ is a smooth function that is the identity on the edges AB, BC inside $\Omega_A$, and whose Hessian maps the outward normal n to a vector with radial component $c_{AB}$ on AB and $c_{AC}$ on AC. It is not difficult to see that one can choose such a function on $\Omega_A$ so that $\nabla \Phi, (\nabla \Phi)^{-1} = O(1)$ and $\nabla^2 \Phi = O(1/r)$ (one can choose a function of the form $\Phi: (r,\theta) \mapsto (f(\theta) r, \theta)$ for a suitable function f). The function w_A defined as above will solve the Neumann conditions, and from the chain rule one has

$\nabla^2 w_A(x) = O( |\nabla \Phi(x)|^2 |\nabla^2 u(\Phi(x))| ) + O( |\nabla^2 \Phi(x)| |\nabla u(\Phi(x))| )$

and so

$\| \nabla^2 w_A \|_{L^2(\Omega_A)} \ll \| \nabla^2 u \|_{L^2(\Omega)} + \| \frac{\nabla u}{|x|} \|_{L^2(\Omega)}.$

A similar estimate (with some additional lower order terms) obtains when one glues together the three partial solutions w_A, w_B, w_C to create a solution w to the boundary problem with

$\| \nabla^2 w \|_{L^2(\Omega)} \ll \| u \|_{H^2(\Omega)} + \sum_{v=A,B,C} \| \frac{\nabla u}{|x-v|} \|_{L^2(\Omega)}.$

To finish the task of estimating $\dot u$ in H^2, one therefore needs bounds on $\| \frac{\nabla u}{|x|} \|_{L^2(\Omega)}$ (and similarly for the other two vertices of the triangle). I think such bounds can be obtained as follows. Let's work in $\Omega_A$ with some suitable smooth cutoff and ignore terms coming from the derivatives of the cutoff (which don't involve the singularity 1/|x| and so are easy to estimate). If we subtract off $u(0) J_0(\sqrt{\lambda} r)$ then u now has mean zero in angular directions. From integration by parts we have

$\int_{\Omega_A} \frac{|\nabla u|^2}{|x|^2} =\int_{\Omega_A} \frac{\lambda |u|^2}{|x|^2} + \int_{\Omega_A} \frac{u \partial_r u}{|x|^3}$
$= \int_{\Omega_A} \frac{\lambda |u|^2}{|x|^2} + \int_{\Omega_A} \frac{|u|^2}{|x|^4}.$

But if A has angle $\alpha$, then the Poincare inequality in the angular direction gives

$\int_{\Omega_A} \frac{|u|^2}{|x|^4} \leq \frac{\alpha^2}{\pi^2} \int_{\Omega_A} \frac{|\nabla u|^2}{|x|^2}$

and so as $\alpha \lt \pi$, we can bound $\int_{\Omega_A} \frac{|\nabla u|^2}{|x|^2}$ as required.

## An alternate approach to the H^2 theory

Yet another approach to the H^2 theory, which seems to give better results, comes from using Schwarz-Christoffel transformations, with the upper half-plane H being the reference domain.

More precisely, given three angles $\alpha+\beta+\gamma=\pi$ for a triangle, we can find a Schwarz-Christoffel transformation $\Phi_{\alpha,\beta}: H \to \Omega(\alpha,\beta)$ from the half-plane $H = \{z: Im(z) \gt 0\}$ to a triangle $\Omega_{\alpha,\beta}$ with angles $\alpha,\beta,\gamma$ by solving the differential equation

$\frac{d}{dz} \Phi_{\alpha,\beta}(z) = \frac{1}{z^{1-\alpha/\pi} (z-1)^{1-\beta/\pi}}$

with initial condition $\Phi_{\alpha,\beta}(0)=0$, keeping $\Phi_{\alpha,\beta}$ conformal. In particular, $\Phi_{\alpha,\beta}(z)$ behaves like $A + \frac{\pi}{\alpha} z^{\alpha/\pi}$ near the origin, $B + \frac{\pi}{\beta} (z-1)^{\beta/\pi}$ near 1, and $C -\frac{\pi}{\gamma} z^{-\gamma/\pi}$ at infinity, for some complex numbers A,B,C with A=0. The image $\Omega_{\alpha,\beta}$ is a triangle ABC with angles $\alpha,\beta,\gamma$.

Write $w = \Phi_{\alpha,\beta}(z)$, then

$dw = \frac{dz}{z^{1-\alpha/\pi} (z-1)^{1-\beta/\pi}}.$

The Lebesgue measure dm_w on the triangle is related to the Lebesgue measure dm_z on the half-plane by the change of variables formula

$\int_{\Omega_{\alpha,\beta}} F(w) dm_w = \int_H \frac{F(\Phi_{\alpha,\beta}(z))}{|z|^{2(1-\alpha/\pi)} |z-1|^{2(1-\beta/\pi)}} dm_z.$

We write the expression $\frac{1}{|z|^{2(1-\alpha/\pi)} |z-1|^{2(1-\beta/\pi)}}$ as $e^{2\omega}$, where

$\omega = (\frac{\alpha}{\pi}-1) \log |z| + (\frac{\beta}{\pi}-1) \log |z-1|.$

The Rayleigh quotient, pulled back to the reference half-plane H, is then

$\int_H |\nabla u|^2\ dm_z / \int_H |u|^2 e^{2\omega}\ dm_z$

and so the second eigenfunction u_2, in the reference half-plane will obey the equation

$-\Delta u_2 = \lambda_2 e^{2\omega} u_2$ (C.1)

with Neumann boundary condition

$n \cdot \nabla u_2 = 0$ (C.2)

on the boundary of the half-plane and the mean zero condition

$\int_H u_2 e^{2\omega}\ dm_z = 0$. (C.3)

We can normalise u_2 in L^2 norm,

$\int_H |u_2|^2 e^{2\omega}\ dm_z = 1$ (C.4)

which by the Rayleigh quotient immediately gives an H^1 bound

$\int_H |\nabla u_2|^2\ dm_z = \lambda$ (C.5)

Now suppose that we vary $\alpha,\beta$ smoothly with respect to a time parameter; this means that $\omega$ will also vary by the formula

$\dot \omega = \frac{\dot \alpha}{\pi} \log |z| + \frac{\dot \beta}{\pi} \log |z-1|.$ (C.6)

The eigenvalue $\dot \lambda_2$ can also vary at a rate which can be controlled by the previous perturbation theory (e.g. equation (14)). As for the variation of the eigenfunction $\dot u_2$, we can compute this by differentiating (C.1) to obtain

$- \Delta \dot u_2 = \dot \lambda_2 e^{2\omega} u_2 + 2 \lambda_2 \dot \omega e^{2 \omega} u_2 + \lambda_2 e^{2\omega} \dot u_2$. (C.7)

Also, on differentating (C.2) we see that $\dot u_2$ conveniently satisfies a homogeneous Neumann condition:

$n \cdot \nabla \dot u_2 = 0$ (C.8)

Using the orthonormal basis $u_k e^\omega$ of $L^2(H, dm_z)$ as in the previous notes on L^2 and H^1 stability, we can obtain an explicit formula for $\dot u_2$. Indeed, integrating (C.7) against another eigenfunction u_l and integrating by parts using (C.1), (C.8) gives

$\lambda_l \int_H \dot u_2 u_l e^{2\omega}\ dm_z= 2 \lambda_2 \int_H \dot \omega e^{2 \omega} u_2 u_l\ dm_z+ \lambda_2 \int_H e^{2\omega} \dot u_2 e^{2\omega} u_l\ dm_z$

and thus

$\int_H \dot u_2 u_l e^{2\omega}\ dm_z= 2 \frac{\lambda_2}{\lambda_l-\lambda_2} \int_H \dot \omega e^{2 \omega} u_2 u_l\ dm_z$

and hence by the eigenfunction expansion

$\dot u_2 = \sum_{l \neq 2} 2 (\frac{\lambda_2}{\lambda_l-\lambda_2} \int_H \dot \omega e^{2 \omega} u_2 u_l\ dm_z) u_l$

which among other things gives the L^2 identity

$\int_H (\dot u_2)^2 e^{2\omega}\ dm_z = \sum_{l \neq 2} 4 |\frac{\lambda_2}{\lambda_l-\lambda_2} \int_H \dot \omega e^{2 \omega} u_2 u_l\ dm_z|^2$

and hence by the Bessel inequality

$\int_H (\dot u_2)^2 e^{2\omega}\ dm_z \leq \frac{4 \lambda_2^2}{|\lambda_3-\lambda_2|^2} X$ (C.9)

where X is the quantity

$X : =\int_H |\dot\omega|^2 e^{2\omega} u_2^2\ dm_z$. (C.10)

The quantity cannot quite be bounded directly by (C.4) because of the logarithmic divergence in the weight $|\dot \omega|^2$, but it should be controllable by also employing (C.5) and (C.3), perhaps after transferring back to the triangle ABC.

The quantity X also controls $\dot \lambda_2$. Indeed, on integrating (C.7) against u_2 and using (C.1), (C.2), (C.4) we obtain that

$0 = \dot \lambda_2 + 2 \lambda_2 \int_H \dot \omega e^{2\omega} u_2\ dm_z$

and hence by Cauchy-Schwarz and (C.4)

$|\dot \lambda_2| \leq 2 \lambda_2 X^{1/2}.$ (C.11)

If we now multiply (C.7) by $e^{-\omega}$ and then take L^2 norms using the triangle inequality and the previous estimates (C.9), (C.10), (C.11) we conclude that

$(\int_H |\Delta \dot u_2|^2 e^{-2\omega}\ dm_z)^{1/2} \leq 2\lambda_2 X^{1/2} + 2 \lambda_2 X^{1/2} + \lambda_2 (\frac{4\lambda_2^2}{|\lambda_3-\lambda_2|^2} X)^{1/2}$

which simplifies to

$(\int_H |\Delta \dot u_2|^2 e^{-2\omega}\ dm_z)^{1/2} \leq \frac{4 \lambda_2 \lambda_3}{\lambda_3-\lambda_2} X^{1/2}.$

We can transfer this back to the triangle $\Omega_{\alpha,\beta}=ABC$ by the change of variables

$\dot v(w) := \dot u( \Phi_{\alpha,\beta}^{-1}(w))$

giving

$(\int_{\Omega_{\alpha,\beta}} |\Delta v|^2\ dm_w)^{1/2} \leq \frac{4\lambda_2 \lambda_3}{\lambda_3-\lambda_2} X^{1/2}$

which after two integration by parts (Bochner-Weitzenbock inequality) gives an H^2 bound on v:

$(\int_{\Omega_{\alpha,\beta}} |\nabla^2 v|^2\ dm_w)^{1/2} \leq \frac{4\lambda_2 \lambda_3}{\lambda_3-\lambda_2} X^{1/2}.$

By Sobolev embedding, (noting also that we have L^2 and H^1 control on v) we can thus obtain L^infty control on v (or equivalently, u) in terms of the quantity X. So the main task is to understand how to bound X, which in the ABC coordinates is basically the L^2 norm of v weighted by a reasonably explicit logarithmic weight.

This is done more explicitly at this set of notes.

The Bessel function $J_\nu(x)$ for $\nu \geq 0$ has the Weierstrass factorisation

$J_\nu(x) = \frac{x^\nu}{2^\nu \Gamma(\nu+1)} \prod_{i=1}^\infty (1 - \frac{x^2}{j_{\nu,i}^2})$

where $j_{\nu,i}$ are the positive zeroes of $J_\nu$; see e.g. this link.

In particular, we have

$J_\nu(x)/J_\nu(y) = (x/y)^\nu \prod_{i=1}^\infty (1 - \frac{x^2}{j_{\nu,i}^2})/(1 - \frac{y^2}{j_{\nu,i}^2}).$

It is known that $j_{\nu,i}$ are non-decreasing in $\nu$ (page 508 of WATSON, G.N., A Treatise on the Theory of Bessel Functions, 2nd ed., Cambridge University Press, 1944; one can also use here the Sturm comparison theorem). Thus if $0 \leq x \leq y \leq j_{\nu_0,1}$ and $\nu \geq \nu_0$, then $(1 - \frac{x^2}{j_{\nu,i}^2})/(1 - \frac{y^2}{j_{\nu,i}^2})$ is non-increasing in $\nu$, and thus

$J_\nu(x)/J_\nu(y) \leq (x/y)^{\nu-\nu_0} J_{\nu_0}(x)/J_{\nu_0}(y)$. (1)

Now, consider a solution to the eigenfunction equation $-\Delta u = \lambda u$ in the sector $\{ (r,\theta): 0 \leq r \leq R; 0 \leq \theta \leq \alpha\}$ for some acute $0 \lt \alpha \lt \pi/2$ and some radius $R\gt0$, obeying the Neumann conditions $\partial_\theta u = 0$ for $\theta=0,\alpha$. For technical reasons we will also need the upper bound

$\sqrt{\lambda} R \leq j_{\pi/\alpha,1}$ (2)

which seems to hold in practice.

By separation of variables, we can write u in polar coordinates as

$u(r,\theta) = \sum_{k=0}^\infty c_k J_{\pi k/\alpha}(\sqrt{\lambda} r) \cos(\pi k \theta/\alpha)$

for some coefficients $c_k$. From (2) and the triangle inequality we thus have the bound

$|u(r_0,\theta) - u(0) J_0(\sqrt{\lambda} r_0)| \leq \sum_{k=1}^\infty |c_k| J_{\pi k/\alpha}(\sqrt{\lambda} r_0)$

for any $0 \leq r_0 \leq R$. To bound the coefficients, we assume an L^2 normalisation

$\int_0^R \int_0^{\alpha} |u(r,\theta)|^2\ rdr \leq 1.$

The left-hand side may be expanded as

$\alpha |u(0)|^2 \int_0^R J_0(\sqrt{\lambda} r)^2 r dr + \frac{\alpha}{2} \sum_{k=1}^\infty |c_k|^2 \int_0^R J_{\pi k/\alpha}(\sqrt{\lambda} r)^2 r dr$

and thus

$\sum_{k=1}^\infty |c_k|^2 \int_0^R J_{\pi k/\alpha}(\sqrt{\lambda} r)^2 r dr \leq A$

where A is the quantity

$A := \frac{2}{\alpha} - 2 |u(0)|^2 \int_0^R J_0(\sqrt{\lambda} r)^2 r dr.$ (3)

Thus by the Cauchy-Schwarz inequality, we have

$|u(r_0,\theta) - u(0) J_0(\sqrt{\lambda} r_0)| \leq r_0^{\pi/\alpha} A^{1/2} B(r_0)^{1/2}$ (4)

where B(r_0) is the quantity

$B(r_0) := \sum_{k=1}^\infty \frac{1}{ \int_0^R (J_{\pi k/\alpha}(\sqrt{\lambda} r)/J_{\pi k/\alpha}(\sqrt{\lambda} r_0))^2 r dr }.$

We may upper bound B(r_0) as follows. We may reduce the range of integration from [0,R] to $[r_1,R]$ for some intermediate radius $r_0 \leq r_1 \leq R$. By (1) and (2), we may replace $\tilde J_{\pi k/\alpha}(\sqrt{\lambda} r)/J_{\pi k/\alpha}(\sqrt{\lambda} r_0)$ by $(r_1/r_0)^{\pi(k-1)/\alpha} \tilde J_{\pi/\alpha}(\sqrt{\lambda} r)/J_{\pi/\alpha}(\sqrt{\lambda} r_0)$. Thus

$B(r_0) \leq \sum_{k=1}^\infty \frac{(r_0/r_1)^{2(k-1)\pi/\alpha} J_{\pi/\alpha}(\sqrt{\lambda} r_0) }{ \int_{r_1}^R J_{\pi/\alpha}(\sqrt{\lambda} r)^2 r dr }.$

One may then sum the geometric series and conclude that

$|u(r_0,\theta) - u(0) J_0(\sqrt{\lambda} r_0)| \leq J_{\pi/\alpha}(\sqrt{\lambda} r_0) A^{1/2} ((r_1/r_0)^{2\pi/\alpha}-1)^{-1/2} C(r_1)^{-1/2}$ (5)

for any $r_0 \leq r_1 \leq R$, where

$C(r_1) := \int_{r_1}^R J_{\pi/\alpha}(\sqrt{\lambda} r)^2 r dr.$

One can use (5) to exclude extrema for small values of r_0. For instance, (5) shows that $|u(r_0,\theta)|$ cannot exceed $|u(0)|$ unless

$|u(0)| (1 - |J_0(\sqrt{\lambda} r_0)|) \geq J_{\pi/\alpha}(r_0) A^{1/2} ((r_1/r_0)^{2\pi/\alpha}-1)^{-1/2} C(r_1)^{-1/2}.$

The left hand side behaves like r_0^2, and the right-hand side like $r_0^{\pi/\alpha}$, so for acute $\alpha$ this inequality becomes impossible for small enough $r_0$ (and "small enough" can be made numerically explicit if one is willing to compute with Bessel functions). Things become bad when $\alpha$ is close to a right angle, but it seems in those cases that the extremum is obtained elsewhere.