From Feynman-Hellman to Hadamard

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.

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.

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)|.


\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.

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}.

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


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).


\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)


\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.

Hadamard formula

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.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s