This text is just a write-up of some fairly horrible calculations I have made recently. The aim is to find how an eigenvalue of the Laplace operator with Dirichlet boundary conditions, let us say the first to fix the ideas, varies when the boundary is deformed. The formula is very classical, and as far as I know was discovered by Jacques Hadamard. With some computations (quite a lot actually), it can be deduced from the Feynman-Hellman formula, which gives the first derivative of the eigenvalues of a self-adjoint operator.

Statement of the problem

Let $\Omega$ be a bounded open set in $\mathbb{R}^{d}$ ($d$ being the dimension) with a $C^{\infty}$ boundary $\partial \Omega$. One can think of the unit disk in $\mathbb{R}^2$ or the unit ball in $\mathbb{R}^3$ (in those cases the computations can actually be carried out). We consider the eigenvalue problem

$\displaystyle \left\{ \begin{array}{rcl} -\Delta u &=& \lambda u \mbox{ in }\Omega,\\ u &=&0 \mbox{ on }\partial \Omega.\\ \end{array} \right.$

We limit ourselves to the first eigenvalue (which is always simple according to a classical result). Now let us assume that $X$ is a $C^{\infty}$ vector field defined in $\mathbb{R}^d$, with compact support . For any real number $t$, we define the mapping $\Phi_t$ from $\mathbb{R}^d$ to $\mathbb{R}^d$ by

$\displaystyle \Phi_t(x)=x+tX(x).$

For $t$ small enough, this actually defines a smooth change of variable, as shown by the following lemma.

Lemma
There exists $\varepsilon>0$ such that for $-\varepsilon \le t \le \varepsilon$, the mapping $\Phi_t:\,\mathbb{R}\mapsto\mathbb{R}^d$ is a $C^{\infty}$-diffeomorphism.

Proof
Let us fix a matrix norm $\|\cdot\|$ on $M(d,\mathbb{R})$, for example the operator norm associated with the euclidean norm on $\mathbb{R}^d$.
For any $x \in \mathbb{R}^d$,

$\displaystyle D\Phi_t(x)=Id+tDX(x),$

where $Id$ is the $d\times d$ identity matrix. Therefore, if

$\displaystyle |t|<\left(\sup_{x\in \mathbb{R}^d}\|DX(x)\|\right)^{-1},$

$D\Phi_t(x)$ is invertible at any $x$. In that case, $\Phi_t$ is a local diffeomorphism, and in particular $\Phi_t(\mathbb{R}^d)$ is an open set.

It remains to show that the mapping $\Phi_t$ is a bijection. We have

$\displaystyle |\Phi_t(x)-\Phi_t(y)|=|x-y+t(X(x)-X(y))| \ge |x-y|-t|X(x)-X(y)|.$

Since

$\displaystyle |X(x)-X(y)|\le \sup_{1\le t \le 1}\|DX((1-t)x+ty)\||x-y|,$

we get

$\displaystyle |x-y|\le (1-t\sup \|DX\|)^{-1}|\Phi_t(x)-\Phi_t(y)|$

(recall that $|t|<\sup\|DX\|^{-1}$). This obviously implies that $\Phi_t$ is injective. It also implies that $\Phi_t(\mathbb{R}^d)$ is a closed set. Indeed, let $(z_n)$ be a sequence of points in $\Phi_t(\mathbb{R^d})$ converging to $z\in \mathbb{R}^d$. There exists a (unique) sequence $(x_n)$ of points in $\mathbb{R}^d$ such that $z_n=\Phi_t(x_n)$. As a convergent sequence, $(z_n)$ satisfies a Cauchy criterion, and therefore so does $(x_n)$ thanks to the preceding inequality. There exists a limit $x$ of $(x_n)$ and since $\Phi_t$ is continuous, $\Phi_t(x)=z$. The set $\Phi_t(\mathbb{R}^d)$ is both open and closed, and obviously nonempty. By connectedness, $\Phi_t(\mathbb{R}^d)=\mathbb{R}^d$.

Let us note $\Omega_t=\Phi_t(\Omega)$ and $\lambda(t)$ the first eigenvalue of the problem

$\displaystyle \left\{\begin{array}{rcl} -\Delta v_t &=& \lambda(t) v_t \mbox{ in }\Omega_t,\\ v_t &=&0 \mbox{ on }\partial \Omega_t.\\ \end{array} \right.$

The case $t=0$ correspond to the initial problem. We want to compute $\lambda'(0)$.

Pullback of the Laplacian

If $v$ is a fonction on $\Omega_t$, we define, using standard notations of differential geometry, the function $u$ on $\Omega$ wich is the \emph{pullback} of $v$:

$\displaystyle u(x)=\Phi_t^* v(x)= v\circ \Phi_t(x).$

This operation gives us an unitary mapping $U_t$ from $L^2(\Omega_t)$ to $L^2(\Omega)$:

$\displaystyle U_tv(x)=|J(\Phi_t)(x)|^{1/2}\Phi_t^*v(x)$

with $J(\Phi_t)(x)=\det(D\Phi_t(x))$ (this is easy to check with the change of variable formula). We then define the operator $T(t)$, depending on $t$, on the \emph{fixed} Hilbert space $L^2(\Omega)$ by

$\displaystyle T(t)=U_t\circ(-\Delta)\circ U_t^*.$

The operator $T(t)$ is selfadjoint, and by invariance of the Sobolev spaces under a change of variable, its domain is $H^1_0(\Omega)\cap H^2(\Omega)$.

Let us now give an explicit formula for $T(t)$. We will use a weak formulation (i.e. consider $T(t)$ as a distribution). Let us fix $v \in H^1_0(\Omega_t)\cap H^2(\omega_t)$ and $\psi \in C^{\infty}_{c}(\Omega_t)$ (i.e. $\psi$ smooth and compactly supported in $\Omega_t$). By the chain rule

$\displaystyle \nabla (v \circ \Phi_t)=(\nabla v \circ \Phi_t)D\Phi_t$

and therefore if $u=U_tv$,

$\displaystyle \nabla v \circ \Phi_t=\nabla (|J(\Phi_t)|^{-1/2}u)(D\Phi_t)^{-1}.$

We set $\varphi=U_t\psi$. We have then

$\displaystyle \langle T(t)u, \varphi \rangle= \langle -\Delta v, \psi \rangle = \langle \nabla v , \nabla \psi \rangle.$

By the change of variable $y=\Phi_t(x)$, we get

$\displaystyle \langle \nabla v, \nabla \psi \rangle=\int_{\Omega_t} \overline{\nabla v } \cdot \nabla \psi \,dy= \int_{\Omega} \overline{\nabla (J(\Phi_t)^{-1/2} u)} \left(D(\Phi_t)D(\Phi_t)^T\right)^{-1} (\nabla (J(\Phi_t)^{-1/2} \varphi))^TJ(\Phi_t)\,dx.$

Since $U_t(C^{\infty}_c(\Omega_t))=C^{\infty}_c(\Omega)$, we get the formula

$\displaystyle T(t)u=-J(\Phi_t)^{-1/2}\mbox{div} \left(J(\Phi_t)\nabla(J(\Phi_t)^{-1/2}u)\left(D(\Phi_t)D(\Phi_t)^T\right)^{-1}\right).$

We can then see $\lambda(t)$ as the first eigenvalue of the problem

$\displaystyle \left\{\begin{array}{rcl} T(t) u_t &=& \lambda(t) u_t \mbox{ in }\Omega,\\ u_t &=&0 \mbox{ on }\partial \Omega.\\ \end{array} \right.$

We have a variable selfadjoint operator acting on a fixed Hilbert space an we want to compute the derivative of its first eigenvalue.

The Feynman-Hellman formula

Formal approach

We first find the formula without worrying about proving the various convergences and regularities or determining the operator domains. In a general setting, let $\mathcal{H}$ be a Hilbert space and $T(t)$ a family of selfadjoint operators acting on $\mathcal{H}$. Assume that for each $t$, $(\lambda(t),u_t)$ is an eigenpair, i.e.

$\displaystyle T(t)u_t=\lambda(t)u_t.$

Differentiating the previous inequality at $t=0$, we get

$\displaystyle T'(0)u_0+T(t)u_0'=\lambda'(0)u_0+\lambda(0)u_0'.$

Let us now take the scalar product of both sides with $u_0$:

$\displaystyle \langle u_0, T(0)u_0' \rangle+\langle u_0, T'(0)u_0' \rangle= \lambda'(0)\|u_0\|^2+\lambda(0)\langle u_0, u_0'\rangle.$

Using the fact that $T(0)$ is selfadjoint and that $T(0)u_0=\lambda(0)u_0$, we find

$\displaystyle \lambda(0)\langle u_0, u_0' \rangle +\langle u_0, T'(0)u_0 \rangle= \lambda'(0)\|u_0\|^2+\lambda(0)\langle u_0, u_0'\rangle$

and finally

$\displaystyle \lambda'(0)=\frac{\langle u_0, T'(0)u_0 \rangle}{\|u_0\|^2}.$

A result

Let us now give a precise statement. It is surely not the most general possible, but it will suffice for our purpose.

Proposition
Let $\mathcal{H}$ be a Hilbert space, $\mathcal{D}$ a subspace of $\mathcal{H}$, and $T(t)$ a family of selfadjoint operators on $\mathcal{H}$. We assume that, for all $-\varepsilon \le t \le \varepsilon$, $T(t)$ as compact resolvent, the domain of $T(t)$ is $\mathcal{D}$, and $\lambda(t)$ is the first eigenvalue of $T(t)$, with $u_t$ an associated eigenvector. We assume additionally that $u_t \to u_0$ when $t\to 0$ and that there exists $v_0 \in \mathcal{H}$ such that

$\displaystyle \frac{(T(t)-T(0))u_0}{t}\to v_0$

when $t\to 0$. Then $t\mapsto \lambda(t)$ is differentiable at $0$ and

$\displaystyle \lambda'(0)=\frac{\langle u_0, v_0\rangle}{\|u_0\|^2}.$

Proof
The hypotheses are of course chosen so that we can mimic the formal computations. Let $t$ be smaller than $\varepsilon$. Then

$\displaystyle \frac{\lambda(t)u_t-\lambda(0)u_0}{t}=\frac{(T(t)-T(0))u_t}{t}+\frac{T(0)(u_t-u_0)}{t}$

and therefore

$\displaystyle \frac{\lambda(t)-\lambda(0)}{t}u_t+\frac{\lambda(0)}{t}(u_t-u_0)=\frac{(T(t)-T(0))u_t}{t}+\frac{T(0)(u_t-u_0)}{t}.$

We take the scalar product of both sides with $u_0$ and get

$\displaystyle \frac{\lambda(t)-\lambda(0)}{t}\langle u_0, u_t \rangle+\frac{\lambda(0)}{t}\langle u_0,u_t-u_0\rangle = \frac{1}{t}\langle u_0, (T(t)-T(0))u_t\rangle+\frac{1}{t}\langle u_0, T(0)(u_t-u_0)\rangle.$

We use the fact that $T(0)$ and $T(t)$ are selfadjoints and the equation $T(0)u_0=\lambda(0)u_0$ and get, after cancellation,

$\displaystyle \frac{\lambda(t)-\lambda(0)}{t}\langle u_0, u_t \rangle=\frac{1}{t}\langle (T(t)-T(0))u_0, u_t\rangle.$

Letting $t$ tend to $0$ gives the desired result.

Application to the original problem

$L^2$-convergence

Using the Courant minimax principle, one can prove (and we will admit here) that $\lambda(t)\to \lambda(0)$ when $t\to 0$.

Let us choose $u_0$ such that $T(0)u_0=\lambda(0)u_0$ and $\|u_0\|_{L^2(\Omega)}=1$. Then, for $-\varepsilon \le t \le \varepsilon$, we choose $u_t$ such that $T(t)u_t=\lambda(t)u_t$, $\|u_t\|=1$, and $\langle u_t, u_0 \rangle \ge 0$. Now let us pick $t_n \to 0$. Since $(u_{t_n})$ is a bounded sequence in $H^1_0(\Omega)$, $(u_{t_n})$ converges weakly in $H^1_0(\Omega)$ to some $u$ (up to a subsequence). Since $H^1(\Omega)$ is compactly embedded in $L^2(\Omega)$, $(u_n)$ converges strongly to $u$. Therefore $\|u\|=1$ and $\langle u, u_0 \rangle \ge 0$. It is easy to see that $u$ satisfies $-\Delta u =\lambda(0)u$ in the sense of distributions. Therefore $u$ is a normalized eigenfunction for $T(0)$ associated with $\lambda(0)$. Since we have assumed $\lambda(0)$ simple, and according to the conditions $\|u\|=1$ and $\langle u, u_0 \rangle \ge 0$, we find $u=u_0$. We have shown that $u_t\to u_0$ when $t\to 0$.

Strong operator convergence

The convergence result for $T(t)$ is obvious if we pay close attention to the form of $T(t)$. The operator $T(t)$ is shown by formula \eqref{eq.PullLap} to be a selfadjoint, elliptic, second order differential operator with smooth coefficients. It can be written

$\displaystyle T(t)=A+tB+t^2R(t)$

where $A$ and $B$ are selfadjoint second order differential operators with smooth coefficients depending only on $x$, and $R(t)$ is a selfadjoint second order with smooth coefficients, bounded in the operator norm $H^2\to L^2$ uniformly in $t$ ($R(t)$ is actually polynomial in $t$). Obviously $A=T(0)$ and therefore, for any $u \in H^2$,

$\displaystyle \frac{(T(t)-T(0))u}{t}\to Bu$

when $t\to 0$.

Asymptotic expansion

It remains to compute explicitly the operator $B$. It is again simpler to use the weak formulation. Let $u$ be in $H^2(\Omega)\cap H^1_0(\Omega)$ and $\varphi$ be in $C^{\infty}_c(\Omega)$. We fist note that

$\displaystyle D\Phi_t=Id+tDX$

and deduce

$\displaystyle J(\Phi_t)=1+t\mbox{Tr}(DX)+O(t^2)=1+t\mbox{div}(X)+O(t^2).$

We get

$\displaystyle \nabla (J(\Phi_t)^{-1/2}\varphi)=\nabla \varphi-\frac{t}{2} \nabla(\mbox{div}(X)\varphi)+O(t^2).$

and

$\displaystyle \left((D\Phi_t)(D\Phi_t)^T\right)^{-1}=\left((Id+t(DX+DX^T)+O(t^2)\right)^{-1}=Id-t(DX+DX^T)+O(t^2)$

Then

$\displaystyle \int_{\Omega}\overline{\nabla(J(\Phi_t)^{-1/2}u)}\left((D\Phi_t)(D\Phi_t)^T\right)^{-1}\nabla(J(\Phi_t)^{-1/2}\varphi)J(\Phi_t)\,dx= \int_{\Omega}\overline{\nabla u}\cdot \nabla \varphi \,dx+t\int_{\Omega}\mbox{div}(X)\overline{\nabla u}\cdot\nabla \varphi\,dx-t\int_{\Omega}\overline{\nabla u}\left(DX+DX^{T}\right)\nabla \varphi^{T}\,dx -\frac{t}{2}\int_{\Omega}\left(\overline{\nabla (\mbox{div}(X)u)}\cdot \varphi+\overline{\nabla u}\cdot \nabla (\mbox{div}(X)\varphi)\right)\,dx+O(t^2).$

We get the quadratic form associated with the operator $B$:

$\displaystyle \langle Bu, \varphi \rangle= -\int_{\Omega}\overline{\nabla u}\left(DX+DX^{T}\right)\nabla \varphi^{T}\,dx + \int_{\Omega}\mbox{div}(X)\overline{\nabla u}\cdot \nabla \varphi \,dx-\frac{1}{2}\int_{\Omega}\left(\overline{\nabla (\mbox{div}(X)u)}\cdot \varphi+\overline{\nabla u}\cdot \nabla (\mbox{div}(X)\varphi)\right)\,dx= -\int_{\Omega}\overline{\nabla u}\left(DX+DX^{T}\right)\nabla \varphi^{T}\,dx-\frac{1}{2}\int_{\Omega}\nabla(\mbox{div}(X))\cdot\left(\varphi\overline{\nabla u}+\overline{u}\nabla \varphi\right)\,dx.$

We get the second order differential operator

$\displaystyle Bu=\mbox{div}\left(\nabla u \left(DX+DX^T\right)\right)-\mbox{div}\left(\mbox{div}(X)\nabla u\right) +\nabla (\mbox{div}(X))\cdot \nabla u +\frac{1}{2}\mbox{div}(\nabla(\mbox{div}(X)))u.$

According to the proposition, the function $t\to \lambda(t)$ is differentiable and

$\displaystyle \lambda'(0)=-\int_{\Omega}\overline{\nabla u_0}\left(DX+DX^{T}\right)\nabla u_0^{T}\,dx -\frac{1}{2}\int_{\Omega}\nabla(\mbox{div}(X))\cdot\left(u_0\overline{\nabla u_0}+\overline{u_0}\nabla u_0\right)\,dx.$

We can recover the classical Hadamard formula from the previous one with an integration by part. We have

$\displaystyle \int_{\Omega} \overline{\nabla u_0}DX \nabla u_0^T\,dx=\sum_{j=1}^d\int_{\Omega} \overline{\partial_{x_j}u_0}\nabla X_j\cdot \nabla u_0 \,dx= \int_{\partial \Omega} (\overline{\nabla u_0}\cdot X)(\nabla u_0 n)\,ds-\sum_{j=1}^d\int_{\Omega} \nabla (\overline{\partial_{x_j}u_0})X_j\cdot \nabla u_0 \,dx-\sum_{j=1}^d\int_{\Omega}\overline{\partial_{x_j} u_0}X_j\Delta u_0\,dx$

and therefore

$\displaystyle -\int_{\Omega}\overline{\nabla u_0}\left(DX+DX^{T}\right)\nabla u_0^{T}\,dx = -2\int_{\partial \Omega}\left|\frac{\partial u_0}{\partial n}\right|^2(X\cdot n)\,ds+\int_{\Omega} \sum_{j=1}^d\left(X_j\partial_{x_j}(\overline{\nabla u_0})\cdot \nabla u_0+X_j\overline{\nabla u_0}\cdot \partial_{x_j}\nabla(u_0)\right)\,dx +\int_{\Omega}\left(\Delta u_0\overline{\nabla u_0}+\overline{\Delta u_0}\nabla u_0\right)\cdot X\,dx.$

We get

$\displaystyle -\int_{\Omega}\overline{\nabla u_0}\left(DX+DX^{T}\right)\nabla u_0^{T}\,dx = -2\int_{\partial \Omega}\left|\frac{\partial u_0}{\partial n}\right|^2(X\cdot n)\,ds+\int_{\Omega} \nabla\left(|\nabla u_0|^2\right)\cdot X\,dx +\int_{\Omega}\left(\Delta u_0\overline{\nabla u_0}+\overline{\Delta u_0}\nabla u_0\right)\cdot X\,dx= -\int_{\partial \Omega}\left|\frac{\partial u_0}{\partial n}\right|^2(X\cdot n)\,ds-\int_{\Omega}|\nabla u_0|^2\mbox{div}(X)\,ds+\int_{\Omega}\left(\Delta u_0\overline{\nabla u_0}+\overline{\Delta u_0}\nabla u_0\right)\cdot X\,dx$

On the other hand,

$\displaystyle -\frac{1}{2}\int_{\Omega}\nabla(\mbox{div}(X))\cdot\left(u_0\overline{\nabla u_0}+\overline{u_0}\nabla u_0\right)\,dx= \int_{\Omega}\mbox{div}(X)|\nabla u_0|^2\,dx+\frac{1}{2}\int_{\Omega}\mbox{div}(X)\left(u_0\overline{\Delta u_0}+\overline{u_0}\Delta u_0\right)\,dx= \int_{\Omega}\mbox{div}(X)|\nabla u_0|^2\,dx-\lambda(0)\int_{\Omega}\mbox{div}(X)|u_0|^2\,dx= \int_{\Omega}\mbox{div}(X)|\nabla u_0|^2\,dx+\lambda(0)\int_{\Omega}X\cdot\left(u_0\overline{\nabla u_0}+\overline{u_0}\nabla u_0\right)\,dx= \int_{\Omega}\mbox{div}(X)|\nabla u_0|^2\,dx-\int_{\Omega}X\cdot\left(\Delta u_0\overline{\nabla u_0}+\overline{\Delta u_0}\nabla u_0\right)\,dx.$

We obtain the classical Hadamard formula for the derivative:

$\displaystyle \lambda'(0)=-\int_{\partial\Omega}\left|\frac{\partial u_0}{\partial n}\right|^2(X\cdot n)\,ds.$