# Stability of optimal shapes and convergence of thresholding algorithms in linear and spectral optimal control problems

Antonin Chambolle, Idriss Mazari-Fouquer, Yannick Privat

► To cite this version:

Antonin Chambolle, Idriss Mazari-Fouquer, Yannick Privat. Stability of optimal shapes and convergence of thresholding algorithms in linear and spectral optimal control problems. 2023. hal-04140177

HAL Id: hal-04140177

https://hal.science/hal-04140177

Preprint submitted on 24 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Stability of optimal shapes and convergence of thresholding algorithms in linear and spectral optimal control problems

Antonin Chambolle* Idriss Mazari-Fouquer† Yannick Privat‡§

June 24, 2023

## Abstract

We prove the convergence of the fixed-point (also called thresholding) algorithm in three optimal control problems under large volume constraints. This algorithm was introduced by Céa, Gioan and Michel, and is of constant use in the simulation of $L^\infty - L^1$ optimal control problems. In this paper we consider the optimisation of the Dirichlet energy, of Dirichlet eigenvalues and of certain non-energetic problems. Our proofs rely on new diagonalisation procedure for shape hessians in optimal control problems, which leads to local stability estimates.

**Keywords:** Convergence analysis, Numerical algorithms in optimal control, PDE constrained optimisation, Shape optimisation, Quantitative inequalities, Thresholding scheme.

**AMS classification (2020):** 49M05, 49M41, 49N05, 49Q10.

# 1 Introduction

## 1.1 Scope of the paper: the fixed point algorithm of Céa, Gioan and Michel for optimal control problems

The fixed-point algorithm, also dubbed thresholding algorithm, is ubiquitous in the numerical simulations of optimal control problems [6, 7, 10, 24, 27, 28, 29, 30, 43, 46]. It was first introduced by Céa, Gioan & Michel in the seminal [8]. This algorithm is designed to solve optimal control problems subject to constraints expressed in $L^1$ and $L^\infty$ norms. These problems are part of a broad class of optimal control problems writing:

$$ \max_{f \in \mathcal{F}} \left( J(f) := \int_{\Omega} j(u_f, \nabla u_f) \right) \quad \text{where} \quad \begin{cases} -\nabla \cdot (A\nabla u_f) = F(x, u_f) + H(f, u_f), \\ \mathcal{B}u_f = 0 \end{cases} \qquad (\mathbf{P}) $$

where $\mathcal{B}$ is a boundary conditions operator, $A$ is an elliptic matrix and $H$ accounts for the coupling between the state $u_f$ and the control $f$. The admissible class of controls is of the form

$$ \mathcal{F} := \left\{ f \in L^\infty(\Omega) : \ 0 \le f \le 1 \text{ a.e., } \int_{\Omega} f = V_0 \right\} $$

<footnote>*CEREMADE, UMR 7534, CNRS & Université Paris Dauphine PSL, 75775 Paris Cedex 16, France (chambolle@ceremade.dauphine.fr)</footnote>

<footnote>†CEREMADE, UMR 7534, Université Paris Dauphine PSL, 75775 Paris Cedex 16, France (mazari@ceremade.dauphine.fr)</footnote>

<footnote>‡IRMA, Université de Strasbourg, CNRS UMR 7501, Inria, 7 rue René Descartes, 67084 Strasbourg, France (yannick.privat@unistra.fr).</footnote>

<footnote>§Institut Universitaire de France (IUF).</footnote>
for some volume constraint $V_0 \in (0; 1)$. In most cases, no explicit description of the maximisers is available and the algorithm proposed in [8] can be recast in terms of the so-called *switch function*: provided all functionals at hand are regular enough we may represent, for a given $f \in \mathcal{F}$, the derivative of the criterion through a function $p_f$, called the *switch function of the optimal control problem*: for any admissible perturbation $h$ at $f$ (for the sake of simplicity, such that for any $\varepsilon > 0$ small enough $f + \varepsilon h \in \mathcal{F}$; see Definition 3 for the proper notion), the derivative of $J$ at $f$ in the direction $h$ writes

$$ \dot{J}(f)[h] = \int_{\Omega} hp_f. $$

The thresholding algorithm, in short, picks the direction $h$ that maximises $\dot{J}(f)[\cdot]$:

<caption>**Algorithm 1** Thresholding algorithm</caption>

1: Initialisation at $f_0 \in \mathcal{F}$
2: $k \leftarrow 0$
3: Compute $p_{f_k}$
4: Compute $c_k$ such that $\mathrm{Vol}(\{p_{f_k} > c_k\}) = V_0\mathrm{Vol}(\Omega)$.
5: $f_k \leftarrow \mathbb{1}_{\{p_{f_k} > c_k\}}$
6: $k \leftarrow k + 1$.

Some difficulties can arise. Indeed, we can only hope for this algorithm to be well-posed if the solutions to the optimal control problem (P) only has so-called bang-bang solutions i.e. any optimal $f^*$ writes $f^* = \mathbb{1}_{E^*}$ for some measurable $E^*$. Indeed, in general, one can only define a Lagrange multiplier $c_k$ such that

$$ \mathrm{Vol}(\{p_{f_k} > c_k\}) \leq V_0\mathrm{Vol}(\Omega) < \mathrm{Vol}(\{p_{f_k} \geq c_k\}). $$

For the algorithm to be properly defined, one needs to ensure that the level-set $\{p_{f_k} = c_k\}$ has zero measure, which “usually” implies that (P) indeed satisfies this bang-bang property. Observe that some simple linear control problems do not satisfy this property [42]. But even assuming this is not a problem, and despite the fact that this method has proved very efficient in past contributions, the theoretical convergence of the sequence of iterates was, to the best of our knowledge, never proved apart from cases where explicit computations are available; see [29] where a detailed study of the convergence rate of the method in the one-dimensional case is undertaken. In [8, Condition (2.4)] a sufficient condition for the convergence of the algorithm is given for very specific types of functionals; some of the conditions given in this paper that ensure convergence are similar and are akin to a coercivity condition in shape differentiation [16]. Finally, we mention that in [18] a study of convergence rate for discretised version of shape optimisation algorithms is obtained.

The purpose of this paper is to study Algorithm 1 for three optimal control problems and to establish its global convergence for large volume constraints $V_0$. This is done by building on recent progress in the study of quantitative inequalities in optimal control theory [37, 39], and by introducing a new diagonalisation procedure for shape hessians in optimal control theory. Such a procedure is of independent interest and relates to stability estimates in optimal control theory.

**Plan of the article** In section 1.2 we lay out our notations and introduce our admissible classes. We study three different problems; the introduction is divided accordingly.

• First, we consider the optimisation of the Dirichlet energy, that is, the optimisation problem

$$ \min_{f \in \mathcal{F}} \frac{1}{2} \int_{\Omega} |\nabla u_f|^2 - \int_{\Omega} f u_f \quad \text{subject to} \quad \begin{cases} -\Delta u_f = f & \text{in } \Omega, \\ u_f = 0 & \text{on } \partial \Omega. \end{cases} $$
We refer to section 1.3 of this introduction.

* Second, we consider the optimisation of weighted Dirichlet eigenvalues, that is, the optimisation problem

$$ \min_{f \in \mathcal{F}} \lambda(f) \text{ where } \lambda(f) \text{ is the first Dirichlet eigenvalue of } -\Delta - f $$

We underline that this is a bilinear control problem. We refer to Section 1.4.

* Third, we consider a non-energetic optimal control problem: for a given function $j = j(x, u)$ that is convex in $u$ we seek to solve

$$ \max_{f \in \mathcal{F}} J(f) := \int_{\Omega} j(x, u_f) \text{ subject to } \begin{cases} -\Delta u_f = f & \text{in } \Omega, \\ u_f = 0 & \text{on } \partial \Omega. \end{cases} $$

We refer to section 1.5.

For the sake of readability, we give in section 1.6 the general plan of the proof of convergence. The core of the paper is devoted to the proofs of the result for the Dirichlet energy. As the proofs for the other two problems are very similar we postpone them to appendices. In the conclusion of this article, we discuss several limits of our analysis and offer possible future research directions.

**Related algorithms.** Before we move on to the main part of the introduction, let us mention that Algorithm 1 is linked to the thresholding algorithm, studied in detail since its introduction in [40] and which, broadly speaking, seeks to approximate mean-curvature motions using a thresholding of the solutions to certain PDEs. The theoretical and numerical aspects of this scheme have been the subject of an intense research activity in the past years [4, 19, 25, 26, 31, 32, 44, 45]. Let us mention that most of the methods used in the aforementioned works can not be used here as they use the behaviour of solutions of PDEs defined in $\mathbb{R}^d$. The presence of boundary conditions in the models under consideration here prohibits relying on the tools these authors developed.

Finally, it should be mentioned that a closely related fixed-point type algorithm was used in [3, 2] to numerically solve topological optimization problems.

## 1.2 Notational conventions, preliminary definitions and setting

Throughout the paper, $\Omega$ is a fixed $\mathcal{C}^2$ bounded domain of $\mathbb{R}^d$ ($d \ge 2$) which is simply connected (in particular $\partial \Omega$ is a connected smooth submanifold of $\mathbb{R}^d$). For a fixed volume constraint $V_0 \in (0; 1)$ we define the set of admissible controls

$$ \mathcal{F}(V_0) := \left\{ f \in L^{\infty}(\Omega) : 0 \le f \le 1 \text{ a.e., } \int_{\Omega} f = V_0 \right\}. \qquad (\text{Adm}) $$

By [22, Proposition 7.2.17] the set $\mathcal{F}(V_0)$ is convex and closed for the weak $L^\infty$ – * topology. Moreover, its extreme points are characteristic functions of subsets: using $\text{Extr}(X)$ to denote the extreme points of a convex set $X$ we have

$$ \text{Extr}(\mathcal{F}(V_0)) = \{1_E, E \subset \Omega, \text{Vol}(E) = V_0 \text{Vol}(\Omega)\}. $$

It is convenient to introduce a bit of terminology:

**Definition 1.** A bang-bang function of $\mathcal{F}(V_0)$ is an extreme point of $\mathcal{F}(V_0)$.

In other words, $m$ is bang-bang in $\mathcal{F}(V_0)$ if, and only if, there exists a measurable subset $E \subset \Omega$ such that $\text{Vol}(E) = V_0 \text{Vol}(\Omega)$ and such that $m = 1_E$.
**Remark 2.** A salient feature of all the optimisation problems we consider in this paper is that their solutions are bang-bang functions. We explain in the conclusion why we can not at this point bypass this underlying assumption.

Finally, as our proofs rely on optimality conditions, we need to introduce the tangent cone to an element $f \in \mathcal{F}(V_0)$.

**Definition 3.** Let $f \in \mathcal{F}(V_0)$. The tangent cone to $\mathcal{F}(V_0)$ at $f$ is the set of functions $h \in L^\infty(\Omega)$ such that, for any sequence $\{t_k\}_{k \in \mathbb{N}}$ converging to 0, there exists a sequence $\{h_k\}_{k \in \mathbb{N}}$ that converges strongly to $h$ in $L^2(\Omega)$ and such that, for any $k \in \mathbb{N}$, $f + t_k h_k \in \mathcal{F}(V_0)$. Functions belonging to this tangent cone are called admissible perturbations at $f$.

Further characterisations of the tangent cone are available in [5, 15].

## 1.3 Optimisation of the Dirichlet energy

**The energy functional and the optimal control problem.** The first problem we study is the optimisation of the Dirichlet energy. For a given $f \in L^2(\Omega)$, let $u_f \in W_0^{1,2}(\Omega)$ be the unique solution of

$$ \begin{cases} -\Delta u_f = f & \text{in } \Omega, \\ u_f = 0 & \text{on } \partial\Omega. \end{cases} \qquad (1.1) $$

This equation has a variational formulation: $u_f$ is the unique minimiser in $W_0^{1,2}(\Omega)$ of the energy $\mathcal{E}_f$ defined as

$$ \mathcal{E}_f : W_0^{1,2}(\Omega) \ni u \mapsto \frac{1}{2} \int_{\Omega} |\nabla u|^2 - \int_{\Omega} fu. $$

The Dirichlet energy associated with $f$ is defined by

$$ \mathcal{E}(f) := \mathcal{E}_f(u_f) = \frac{1}{2} \int_{\Omega} |\nabla u_f|^2 - \int_{\Omega} f u_f. \qquad (1.2) $$

The optimal control problem is

$$ \min_{f \in \mathcal{F}(V_0)} \mathcal{E}(f). \qquad (\mathbf{P}_{\text{Dir}}) $$

The following result is standard.

**Lemma 4.** The problem $(\mathbf{P}_{\text{Dir}})$ has a solution. Furthermore, the map $\mathcal{E}$ is strictly concave. In particular any solution $f^*$ of $(\mathbf{P}_{\text{Dir}})$ is a bang-bang function in the sense of Definition 1: there exists $E^* \subset \Omega$ such that $f^* = \mathbb{1}_{E^*}$, with $\text{Vol}(E^*) = V_0 \text{Vol}(\Omega)$.

**The switch function and the thresholding algorithm.** To describe the thresholding algorithm we need to specify the switch function of $\mathcal{E}$. The following lemma is well-known.

**Lemma 5.** The map $f \mapsto \mathcal{E}(f)$ is Fréchet differentiable at any $f \in \mathcal{F}(V_0)$ and, for any $f \in \mathcal{F}(V_0)$, for any admissible perturbation $h$, there holds

$$ \dot{\mathcal{E}}(f)[h] = - \int_{\Omega} h u_f. $$

Consequently, for the Dirichlet energy, the switch function of the optimal control problem is $-u_f$. This allows to describe the thresholding algorithm:
**Algorithm 2** Thresholding algorithm for the Dirichlet energy

1: Initialisation at $f_0 \in \mathcal{F}$

2: $k \leftarrow 0$

3: Compute $u_{f_k}$

4: Compute $c_k$ such that $\mathrm{Vol}(\{u_{f_k} > c_k\}) = V_0\mathrm{Vol}(\Omega)$.

5: $f_k \leftarrow \mathbb{1}_{\{u_{f_k} > c_k\}}$

6: $k \leftarrow k + 1$.

**Remark 6.** Of course we have the same problem as in Algorithm 1, namely, that we would need to ensure that the algorithm is well-defined in the sense that for any index $k$ there indeed exists a $c_k \in \mathbb{R}$ such that $\mathrm{Vol}(\{u_{f_k} > c_k\}) = V_0\mathrm{Vol}(\Omega)$. This is also covered by our theorem.

**Some terminology.** We introduce the following definition:

**Definition 7 (Critical points).** For any $V_0 \in (0; 1)$, for any $f \in \mathcal{F}(V_0)$, we say $f$ is a **critical point of $\mathcal{E}$** if:

1. There exists a unique $\mu_f = \mu_{u_f, V_0}$ such that

$$ \mathrm{Vol}(\{u_f > \mu_f\}) = V_0\mathrm{Vol}(\Omega) = \mathrm{Vol}(\{u_f \ge \mu_f\}), $$

2. $f = \mathbb{1}_{\{u_f > \mu_f\}}$.

A **critical set** is a subset $E$ of $\Omega$ such that $\mathbb{1}_E \in \mathcal{F}(V_0)$ and such that $\mathbb{1}_E$ is a critical point of $\mathcal{E}$.

It is natural to expect that the algorithm will converge to a local minimiser of the functional. In general, one must clarify the meaning behind “local minimiser”; more than mere local minimisers we use the notion of “stable local minimisers”. Our definition here is the following:

**Definition 8.** An admissible control $f^* \in \mathcal{F}(V_0)$ is called a stable local minimiser of $\mathcal{E}$ in $\mathcal{F}(V_0)$ if there exist $C = C(f^*)$, $\delta = \delta(f^*) > 0$ such that

$$ \forall f \in \mathcal{F}(V_0), \|f - f^*\|_{L^1(\Omega)} \le \delta(f^*) \Rightarrow \mathcal{E}(f) \ge \mathcal{E}(f^*) + C(f^*) \|f - f^*\|_{L^1(\Omega)}^2. $$

**Main result** Our main result is the following:

**Theorem 1.** There exists $\varepsilon > 0$ such that, if $1 - \varepsilon \le V_0 < 1$, for any initialisation $f_0 \in \mathcal{F}(V_0)$, the sequence $\{f_k\}_{k \in \mathbb{N}}$ generated by Algorithm 2 converges strongly in $L^1(\Omega)$ to a stable local minimiser of $\mathcal{E}$ in $\mathcal{F}(V_0)$.

## 1.4 Optimisation of weighted eigenvalues

**The energy functional and the optimal control problem.** The second functional under consideration is the principal eigenvalue of a Dirichlet operator. For any $f \in \mathcal{F}(V_0)$, let $\lambda(f)$ be the first eigenvalue of the operator $-\Delta - f$ endowed with Dirichlet boundary conditions. The variational formulation of $\lambda(f)$ is

$$ \lambda(f) = \min_{u \in W_0^{1,2}(\Omega) \setminus \{0\}} \frac{\int_{\Omega} |\nabla u|^2 - \int_{\Omega} f u^2}{\int_{\Omega} u^2}. \quad (1.3) $$
It is standard to see that this eigenvalue is simple and that any associated eigenfunction has constant sign. Thus, up to normalisation, the eigenpair $(\eta_f, \lambda(f))$ satisfies

$$ \begin{cases} -\Delta\eta_f = \lambda(f)\eta_f + f\eta_f & \text{in } \Omega, \\ \eta_f \in W_0^{1,2}(\Omega), \\ \eta_f > 0 & \text{in } \Omega, \\ \int_{\Omega} \eta_f^2 = 1. \end{cases} \qquad (1.4) $$

The optimal control problem under consideration is

$$ \min_{f \in \mathcal{F}(V_0)} \lambda(f). \qquad (\mathbf{P}_{\text{Eigen}}) $$

This problem has been studied at length, see [21, Chapter 8] or the more recent [30]. Part of the recent interest in it is due to its applications in spatial ecology, as the sign of $\lambda(f)$ predicts extinction or survival of a species, see [36, Introduction] and the references therein. The thresholding algorithm was applied to this problem in [24, 27, 30].

Similar to Lemma 4, we have the following result.

**Lemma 9.** The problem $(\mathbf{P}_{\text{Eigen}})$ has a solution. Furthermore, the map $\lambda$ is strictly concave. In particular any solution $f^*$ of $(\mathbf{P}_{\text{Eigen}})$ is a bang-bang function in the sense of Definition 1: there exists $E^* \subset \Omega$ such that $f^* = \mathbb{1}_{E^*}$, with $\text{Vol}(E^*) = V_0\text{Vol}(\Omega)$.

**The switch function and the thresholding algorithm.** Let us now give the switch function of $(\mathbf{P}_{\text{Eigen}})$.

**Lemma 10.** The map $f \mapsto \lambda(f)$ is Fréchet differentiable at any $f \in \mathcal{F}(V_0)$ and, for any $f \in \mathcal{F}(V_0)$, for any admissible perturbation $h$ at $f$, there holds

$$ \dot{\lambda}(f)[h] = -\int_{\Omega} h\eta_f^2. $$

Consequently, for the Dirichlet eigenvalue, the switch function is $-\eta_f^2$. Note that as $\eta_f > 0$, computing a level-set of $\eta_f$ is the same as computing a level-set of $\eta_f^2$. This allows to describe the thresholding algorithm for $(\mathbf{P}_{\text{Eigen}})$.

<caption>**Algorithm 3** Thresholding algorithm for the Dirichlet eigenvalue</caption>

1: Initialization at $f_0 \in \mathcal{F}$

2: $k \leftarrow 0$

3: Compute $\eta_{f_k}$

4: Compute $c_k$ such that $\text{Vol}(\{\eta_{f_k} > c_k\}) = V_0\text{Vol}(\Omega)$.

5: $f_k \leftarrow \mathbb{1}_{\{\eta_{f_k} > c_k\}}$

6: $k \leftarrow k + 1$.

Remark 6 about the well-posedness of the algorithm applies here as well.

To alleviate the presentation, we give the precise definitions of “critical point” and of “local minimisers” at the beginning of the proof; they are identical to Definitions 7-8, up to replacing $\mathcal{E}$ with $\lambda$.
**Main result** The main theorem is the following:

**Theorem II.** There exists $\varepsilon > 0$ such that, if $1 - \varepsilon \le V_0 < 1$, for any initialisation $f_0 \in \mathcal{F}(V_0)$, the sequence $\{f_k\}_{k \in \mathbb{N}}$ generated by Algorithm 3 converges strongly in $L^1(\Omega)$ to a local minimiser of $\lambda$ in $\mathcal{F}(V_0)$.

The proof is very similar to that of Theorem I; it is given in [9, Section A].

We conclude this section with an illustration of the thresholding algorithm 3 for optimizing the functional $\lambda$. In Figures 2 and 3, we show the optimal domain, as well as the evolution of the quantity $\|f_{k+1} - f_k\|_{L^1(\Omega)}$, using the same notation as in algorithm 3, which is used as a stopping criterion for the algorithm, as a function of the iterations. We observe that the number of iterations required for the stopping criterion to be below the fixed tolerance (in this case, $10^{-6}$ for these figures) is very small. The numerical results also suggest that the statement of Theorem II is likely to be valid for large values of $\varepsilon$ close to 1.

<figure><img src="image_2.png" /><figcaption>Figure 1: Domain $E_0$ chosen to initialize algorithm 3, i.e., $f_0 = \mathbb{1}_{E_0}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: $V_0 = 0.3$. Left: representation of the optimal domain (black contour) and the eigenfunction (warmer colors indicate where the function takes its maximum values). Right: evolution of the quantity $\|f_{k+1} - f_k\|_{L^1(\Omega)}$ as a function of iterations, on a logarithmic scale.</figcaption></figure>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: $V_0 = 0.8$. Left: representation of the optimal domain (black contour) and the eigenfunction (warmer colors indicate where the function takes its maximum values). Right: evolution of the quantity $||f_{k+1} - f_k||_{L^1(\Omega)}$ as a function of iterations, on a logarithmic scale.</figcaption></figure>

## 1.5 Optimisation of non-energetic criteria

**The functional and the optimal control problem.** The third and final problem we consider is the optimisation of a non-energetic criterion. We consider, for any $f \in \mathcal{F}(V_0)$, the solution $u_f$ of the Dirichlet problem (1.1) (the same $u_f$ we used in the study of the Dirichlet energy). We fix a non-linearity $j = j(x, u)$ satisfying the following regularity assumptions:

$$ \left\{ \begin{array}{l} j \text{ is } \mathcal{C}^2 \text{ in } u, \mathcal{C}^1 \text{ in } x, \\ \text{For } k = 0, 1, 2, i = 0, 1, \text{ for any compact } K \subset \subset \mathbb{R} \text{ there holds } \sup_{\bar{\Omega} \times K} \|\partial_{x^i}^i \partial_{u^k}^k j\|_{L^\infty(\Omega)} < \infty, \\ j \text{ is strictly convex in } u \text{ and for any } u > 0 \inf_{\bar{\Omega}} \partial_{u^2}^2 j(x, u) > 0, \\ j \text{ is increasing in } u \text{ for any } x \in \bar{\Omega}, \text{ for any } u > 0, \partial_{u}j(x, u) > 0. \end{array} \right. \qquad (\mathbf{H}_j) $$

The problem under consideration is

$$ \max_{f \in \mathcal{F}(V_0)} \left( J(f) := \int_{\Omega} j(x, u_f(x)) dx \right). \qquad (\mathbf{P}_j) $$

We start with a standard result.

**Lemma 11.** The problem $(\mathbf{P}_j)$ has a solution. Furthermore, the map $J$ is strictly convex. In particular any solution $f^*$ of $(\mathbf{P}_{\text{Eigen}})$ is a bang-bang function in the sense of Definition 1: there exists $E^* \subset \Omega$ such that $f^* = \mathbb{1}_{E^*}$, with $\text{Vol}(E^*) = V_0 \text{Vol}(\Omega)$.

**The switch function and the thresholding algorithm.** We now describe the switch function in the following standard lemma:

**Lemma 12.** The map $J$ is Fréchet-differentiable and, for any $f \in \mathcal{F}(V_0)$, for any admissible perturbation $h$ at $f$ there holds

$$ \dot{j}(f)[h] = \int_{\Omega} hp_f $$
where $p_f$ is the unique solution of

$$ \begin{cases} -\Delta p_f = \partial_u j(x, u_f) & \text{in } \Omega, \\ p_f \in W_0^{1,2}(\Omega). \end{cases} \qquad (1.5) $$

Lemma 12 states that the function $p_f$ defined in (1.5) is the switch function of the optimal control problem ($P_j$). This allows to describe the thresholding algorithm for ($P_j$).

**Algorithm 4** Thresholding algorithm for non-energetic problems eigenvalue

1: Initialisation at $f_0 \in \mathcal{F}$

2: $k \leftarrow 0$

3: Compute $p_{f_k}$

4: Compute $c_k$ such that $\mathrm{Vol}(\{p_{f_k} > c_k\}) = V_0 \mathrm{Vol}(\Omega)$.

5: $f_k \leftarrow \mathbb{1}_{\{p_{f_k} > c_k\}}$

6: $k \leftarrow k + 1$.

Remark 6 applies here as well. The notion of “critical point” and of “local minimisers”, which were defined for the Dirichlet energy in Definitions 7-8, are identical for non-energetic criteria: it suffices to replace $\mathcal{E}$ with $J$.

**Main result.** The main result regarding non-energetic problems is the following convergence theorem:

**Theorem III.** Assume $j$ satisfies (H$_j$). There exists $\varepsilon > 0$ such that, if $1 - \varepsilon \le V_0 < 1$, for any initialisation $f_0 \in \mathcal{F}(V_0)$, the sequence $\{f_k\}_{k \in \mathbb{N}}$ generated by Algorithm 4 converges strongly in $L^1(\Omega)$ to a local minimiser of $J$ in $\mathcal{F}(V_0)$.

As it is very similar to the proof of Theorem I, the proof of Theorem III is given in [9, Section B].

## 1.6 The plan for the proof of convergence

As the schemes of proof for Theorems I, II and III are identical and as they are each quite long we now present the general plan. This allows use to single out the relevant elements for each step of the proof.

**Consequence of the convexity of the functional.** The first part of the proof of convergence is to use the convexity or concavity properties of the functional (depending on whether we are maximising or minimising). To be more precise, consider a generic functional $\mathcal{G}: \mathcal{F}(V_0) \to \mathbb{R}$ that we seek to minimise in $\mathcal{F}(V_0)$. Assume $\mathcal{G}$ is Fréchet-differentiable and denote by $-q_f$ its switch function at a given $f \in \mathcal{F}(V_0)$ so that, for any $f$, for any admissible perturbation $h$ at $f$, there holds

$$ \dot{\mathcal{G}}(f)[h] = - \int_{\Omega} q_f h. $$

Starting from $f_k$, we choose $f_{k+1} := \mathbb{1}_{\{q_{f_k} > c_k\}}$ where $\mathrm{Vol}(\{q_{f_k} > c_k\}) = V_0 \mathrm{Vol}(\Omega)$. If the functional $\mathcal{G}$ is concave, then we have the estimate

$$ \mathcal{G}(f_{k+1}) - \mathcal{G}(f_k) \le \dot{\mathcal{G}}(f_k)(f_{k+1} - f_k) = - \int_{\Omega} q_{f_k} (f_{k+1} - f_k). \qquad (1.6) $$

Thus we can be satisfied with gaining enough control of the first order derivative, which is a linear (in the perturbation $h$) functional.
**Consequence of large volume constraints I: quantitative bathtub principle and well-posedness of thresholding algorithms.** This is where the large volume constraint plays a first role. Keeping the same notations as in the previous paragraph, observe that from Lemmata 5-10-12 the switch function $q_f$ satisfies a PDE with Dirichlet boundary conditions; as the PDEs involved enjoy a maximum principle, it is expected that, if $V_0$ is close enough to 1, then the set $\{q_f > c\}$ such that $\text{Vol}(\{q_f > c\}) = V_0$ is close to the boundary, and even $\mathcal{C}^1$ close to $\partial\Omega$. Moreover, it is expected that by elliptic regularity $q_f$ is $\mathcal{C}^{1,\alpha}$ close to $q_{\bar{f}}\equiv 1$, the function associated to $f = 1$ (observe that when $V_0 = 1$ the set of admissible controls is reduced to a single point: $\mathcal{F}(V_0) = \{\bar{f} \equiv 1\}$). Applying the Hopf Lemma to $q_{\bar{f}}$ we may get that, for $V_0$ close to 1, the gradient of $q_f$ is always non-zero in a neighbourhood of $\partial\Omega$. With these two informations combined, we deduce that the level set $\{q_f = c_f\}$ has zero measure and, hence, that the thresholding algorithm is indeed well defined.

But there is another crucial information we can deduce, and it is that $|\nabla q_f| \neq 0$ on $\{q_f = c_f\}$. Combined with the aforementioned regularity of the boundary we are then in a position to apply a quantitative bathtub principle [39, Proposition 26] (see also the related [14]). One information we have not yet given is that the choice of the next iterate $f_{k+1}$ in the thresholding algorithm is that $f_{k+1}$ is actually the unique<fnref n="1" /> maximiser of the functional

$$T: \mathcal{F}(V_0) \ni f \mapsto \int_{\Omega} q_{f_k} f.$$

With this notation (1.6) rewrites $\mathcal{G}(f_{k+1}) - \mathcal{G}(f_k) \le T(f_k) - T(f_{k+1})$. The quantitative bathtub principle writes: under certain regularity assumptions, for a certain constant $C > 0$,

$$\forall f \in \mathcal{F}(V_0), T(f) - T(f_{k+1}) \le -C \|f_{k+1} - f\|_{L^1(\Omega)}^2.$$

Assuming that the constant $C$ is uniform in $k$ this would yield

$$G(f_{k+1}) - G(f_k) \le -C \|f_{k+1} - f_k\|_{L^1(\Omega)}^2$$

and, in turn, prove that the sequence $\{f_k\}_{k \in \mathbb{N}}$ satisfies $\sum_{k \in \mathbb{N}} \|f_{k+1} - f_k\|_{L^1(\Omega)}^2 < \infty$. Finally, we can prove that this implies that, either the sequence $\{f_k\}_{k \in \mathbb{N}}$ converges, or it has an infinite number of closure points. To rule the latter possibility out, we use a shape derivation argument.

**Consequence of large volume constraints II: shape derivative formalism and coercivity of shape hessians.** We use the shape derivative formalism of Hadamard; see [22, Chapter 5] for a full introduction. We can prove that any closure point of the sequence $\{f_k\}_{k \in \mathbb{N}}$ is a bang-bang function $f^* = \mathbb{1}_E$ that is also a critical point (in the sense that $\{q_{f^*} > c_{f^*}\} = E^*$). With an abuse of notation, for a bang-bang function $f = \mathbb{1}_E$, we may define the shape functional

$$\mathcal{G}(E) := \mathcal{G}(\mathbb{1}_E).$$

To show that possible closure points are isolated, which would allow to conclude, we prove that these closure points are stable in the sense of shapes [16]. Namely, considering a critical shape $E$ and smooth enough vector fields $\Phi$, we define shape Lagrangians (with a Lagrange multiplier $\mu = \mu(E^*)$ independent of $\Phi$ and encoded by first order optimality conditions)

$$g_{\Phi}(t) := \mathcal{G}((\text{Id} + t\Phi)E) + \mu \text{Vol}((\text{Id} + t\Phi)E)$$

so that $g'(0) = 0$. We then show a coercivity of second order shape derivatives or, in other words, that there exists a constant $C > 0$ such that for any $\Phi$

$$g_{\Phi}''(0) \ge C \|\langle \Phi, \nu \rangle\|_{X(\partial E)}^2$$

<footnote n="1">¹using the information about the measure of level sets that were already derived.</footnote>
for some norm $X$. In this case we say the second-order derivative is $X$-coercive. Our main contribution in this paper is to find a procedure to diagonalise these second order shape derivatives by using a new interior Steklov system of eigenvalues and eigenfunctions and to prove that, under large volume constraints, these second-order shape derivatives are indeed coercive with respect to the $L^2$-norm. We believe this is of independent interest. Once this is done we can apply the general procedure of [16], adapted in [35, 37, 39] to the setting of optimal control problems, to prove that critical points are isolated.

# 2 Proof of Theorem I

We begin with the study of $(P_{\text{Dir}})$.

## 2.1 First consequences of a large volume constraint

We gather here several crucial properties of the problem under large enough volume constraints. First, we show that when $V_0$ is close enough to 1 the level-sets generated by the thresholding algorithm enjoy strong regularity properties.

To state our results in a synthetic way we let, for any $f \in \mathcal{F}(V_0)$, $\mu_{f,V_0}$ be the unique real number such that

$$ \frac{\text{Vol}(\{u_f > \mu_{f,V_0}\})}{\text{Vol}(\Omega)} \le V_0 \le \frac{\text{Vol}(\{u_f \ge \mu_{f,V_0}\})}{\text{Vol}(\Omega)}. $$

We also define

$$ \omega_{f,V_0} := \{u_f > \mu_{f,V_0}\}. $$

**Lemma 13.** There exists $\varepsilon_0 > 0$ and $\delta_0 > 0$ such that, for any $V_0 \in (1 - \varepsilon_0; 1)$, for any $f \in \mathcal{F}(V_0)$, we have on the one hand

$$ \text{Vol}(\{u_f = \mu_{f,V_0}\}) = 0 $$

and, on the other hand, the following regularity properties:

1. $\omega_{f,V_0}$ has a $\mathcal{C}^{1,\alpha}$ boundary (for any $\alpha \in (0; 1)$),

2.

$$ \min_{\partial\omega_{f,V_0}} \left| \frac{\partial u_f}{\partial \nu} \right| \ge \delta_0, $$

3. $\partial\omega_{f,V_0}$ is locally a graph over $\partial\Omega$.

Finally, for any $\alpha \in (0; 1)$, $\partial\omega_{f,V_0}$ converges to $\partial\Omega$ in $\mathcal{C}^{1,\alpha}$ as $V_0 \to 1$, uniformly in $f \in \mathcal{F}(V_0)$.

**Proof of Lemma 13.** We will proceed in several steps.

**Convergence to the torsion function.** This lemma rests upon the study of the torsion function $w_\Omega$ of $\Omega$, defined as the solution of

$$ \begin{cases} -\Delta w_\Omega = 1 & \text{in } \Omega, \\ w_\Omega = 0 & \text{on } \partial\Omega. \end{cases} \qquad (2.1) $$

By the maximum principle of Hopf there exists $\underline{\delta} > 0$ such that

$$ \min_{\partial\Omega} \left| \frac{\partial w_\Omega}{\partial \nu} \right| \ge \underline{\delta}. \qquad (2.2) $$
Furthermore, observe that standard $L^p$-regularity estimates (see for instance [20, Corollary 9.10]) entail that, for any $p \in [1; +\infty)$, there exists a constant $C_p > 0$ such that for any $V_0 \in (0; 1)$ and any $f \in \mathcal{F}(V_0)$

$$ \|u_f\|_{W^{2,p}(\Omega)} \le C_p. $$

From Sobolev embeddings, we deduce that for any $\alpha \in (0; 1)$ there exists $H_\alpha$ such that for any $V_0 \in (0; 1)$ and any $f \in \mathcal{F}(V_0)$

$$ \|u_f\|_{\mathcal{C}^{1,\alpha}(\Omega)} \le H_{\alpha}. \qquad (2.3) $$

From these regularity estimates, we obtain at once the uniform convergence result

$$ \forall \alpha \in (0; 1), \lim_{V_0 \to 1} \sup_{f \in \mathcal{F}(V_0)} \|u_f - w_\Omega\|_{\mathcal{C}^{1,\alpha}(\Omega)} = 0. \qquad (2.4) $$

**Behaviour of $\mu_{f,V_0}$**. Let us now prove that

$$ \lim_{V_0 \to 1} \sup_{f \in \mathcal{F}(V_0)} |\mu_{f,V_0}| = 0. \qquad (2.5) $$

Argue by contradiction and assume (2.5) does not hold. Then there exist $\gamma > 0$, a sequence $\{V_k\}_{k \in \mathbb{N}}$ that satisfies $V_k \to 1$ and, for any $k \in \mathbb{N}$, a function $f_k \in \mathcal{F}(V_k)$ such that

$$ \forall k \in \mathbb{N}, |\mu_{f_k,V_0}| \ge \gamma > 0. \qquad (2.6) $$

To obtain a contradiction, we first prove that there exists $\beta_0, \varepsilon'_0 \in (0; 1)$ such that

$$ \forall V_0 \in (1 - \varepsilon'_0; 1), \forall f \in \mathcal{F}(V_0), \beta_0 w_\Omega \le u_f \le w_\Omega. $$

That $u_f \le w_\Omega$ simply follows from the Hopf maximum principle as, for any $V_0 \in (0; 1)$ and any $f \in \mathcal{F}(V_0)$ there holds $-\Delta u_f = f \le 1 = -\Delta w_\Omega$. Now fix $\beta_0 \in (0; 1)$. From (2.4) and (2.2) there exists $\varepsilon'_0 > 0$ such that for any $V_0 \in (1 - \varepsilon'_0; 1)$ and any $f \in \mathcal{F}(V_0)$ we have

$$ \beta_0 \max_{\partial \Omega} \left| \frac{\partial w_{\Omega}}{\partial \nu} \right| < \min_{\partial \Omega} \left| \frac{\partial u_f}{\partial \nu} \right|. \qquad (2.7) $$

As for any $\alpha \in (0; 1)$ we have

$$ w_{\Omega} \in \mathcal{C}^{1,\alpha}(\Omega), \quad \sup_{V_0 \in (0;1), f \in V_0} \|u_f\|_{\mathcal{C}^{1,\alpha}(\Omega)} < \infty $$

we deduce that there exists $r'_0 > 0$ small enough such that if we define

$$ \Omega(r'_0) := \{x \in \Omega, \mathrm{dist}(x, \partial\Omega) \le r'_0\} \qquad (2.8) $$

then (2.7) along with the Dirichlet boundary conditions implies

$$ \beta_0 w_{\Omega} \le u_f \text{ in } \Omega(r'_0). $$

From the maximum principle and (2.4) we also conclude that if $\varepsilon'_0 > 0$ is small enough ($r'_0 > 0$ being fixed) we have

$$ \beta_0 w_{\Omega} < u_f \text{ in } \Omega \setminus \Omega(r'_0). $$

Thus we obtain, if $\varepsilon'_0 > 0$ is small enough that, for any $V_0 \in (1 - \varepsilon'_0; 1)$, for any $f \in \mathcal{F}(V_0)$ we have

$$ \beta_0 w_{\Omega} \le u_f. \qquad (2.9) $$
Let us now prove that (2.6) can not hold. From (2.9) we deduce that for any $k \in \mathbb{N}$ we have

$$ \{u_f \ge \mu_{f_k, V_k}\} \subset \left\{w_{\Omega} \ge \frac{\mu_{f_k, V_k}}{\beta_0}\right\} \subset \left\{w_{\Omega} \ge \frac{\gamma}{\beta_0}\right\}. $$

Consequently we deduce that

$$ \forall k \in \mathbb{N}, \mathrm{Vol}(\{u_f \ge \mu_{f_k, V_k}\}) \le \mathrm{Vol}\left(\left\{w_{\Omega} \ge \frac{\gamma}{\beta_0}\right\}\right) \le 1 - s_0 $$

for a fixed $s_0 \in (0; 1)$, since $w_{\Omega} > 0$ in $\Omega$ and $\gamma > 0$. This is a contradiction with the definition of $\mu_{f_k, V_k}$, whence the conclusion.

**Study of $\omega_{f,V_0}$**. From (2.2) and the fact that $w_{\Omega} \in \mathcal{C}^{1,\alpha}$ there exists $r_0 > 0$ such that

$$ |\nabla w_{\Omega}| \ge \frac{\delta}{2} \quad \text{in } \Omega(r_0') $$

where $\Omega(r_0)$ is defined in (2.8). Fix such an $r_0 > 0$. Since $w_{\Omega} > 0$ in $\Omega$, (2.4)-(2.5) imply that there exists $\varepsilon_0'' > 0$ such that if $V_0 \in (1 - \varepsilon_0''; 1)$ then for any $f \in \mathcal{F}(V_0)$

$$ \{u_f = \mu_{f,V_0}\} \subset \Omega(r_0). $$

From (2.4) we also have, up to reducing $\varepsilon_0'' > 0$,

$$ \inf_{V_0 \in (1-\varepsilon_0''; 1), f \in \mathcal{F}(V_0)} \left( \inf_{\Omega(r_0)} |\nabla u_f| \right) \ge \frac{\delta}{4}. $$

From the implicit function theorem we deduce first that $\mathrm{Vol}(\{u_f = \mu_{f,V_0}\}) = 0$, second that for any $\alpha \in (0; 1)$, $\omega_{f,V_0}$ has a $\mathcal{C}^{1,\alpha}$ boundary and, finally, that there exists $\delta_0 > 0$ such that

$$ \min_{\partial\omega_{f,V_0}} \left| \frac{\partial u_f}{\partial \nu} \right| \ge \delta_0. $$

**Convergence of $\partial\omega_{f,V_0}$**. We first prove that $\partial\omega_{f,V_0}$ is (locally) a graph over $\partial\Omega$. Argue by contradiction and assume that there exists a sequence $\{V_k\}_{k \in \mathbb{N}}$ that converges to 1 and, for any $k \in \mathbb{N}$, a function $f_k \in \mathcal{F}(V_k)$, a sequence $\{x_k\}_{k \in \mathbb{N}} \in \partial\Omega^{\mathbb{N}}$ and two sequences $\{y_{k,0}\}_{k \in \mathbb{N}}$, $\{y_{k,1}\}_{k \in \mathbb{N}}$ of points of $\partial\omega_{f_k, V_k}$ such that there exist $(t_{k,i})_{k \in \mathbb{N}, i=0,1} \in (\mathbb{R}_+ \times \{0,1\})^{\mathbb{N}}$ with

$$ y_{k,i} = x_k - t_{k,i}\nu(x_k), \quad t_{k,i} \underset{k \to \infty}{\to} 0, \quad t_{k,0} < t_{k,1}. $$

We consider a closure point $x_{\infty}$ of $\{x_k\}_{k \in \mathbb{N}}$. By the intermediate value theorem, for any $k \in \mathbb{N}$ there exists $s_k \in (t_{k,0}; t_{k,1})$ such that, setting $z_k := x_k - s_k\nu(x_k)$ we have

$$ \langle \nabla u_f(z_k), \nu(x_k) \rangle = 0. $$

Clearly we have $z_k \underset{k \to \infty}{\to} x_{\infty}$ and $\nu(x_k) \underset{k \to \infty}{\to} \nu(x_{\infty})$. We obtain from (2.4)

$$ \langle \nabla w_{\Omega}(x_{\infty}), \nu(x_{\infty}) \rangle = 0, $$

in contradiction with (2.2). Consequently, as $V_0 \to 1$, for any $f \in \mathcal{F}(V_0)$, $\partial\omega_{f,V_0}$ is a local graph over $\partial\Omega$; the fact that it converges to $\partial\Omega$ uniformly in $f \in \mathcal{F}(V_0)$ is a simple consequence of (2.4).

□
We note the following consequence of Lemma 13: if we denote by Lip($\Sigma$) the Lipschitz-constant of a hypersurface $\Sigma$ then there exists $\varepsilon_0 > 0$ and a constant $M$ such that

$$ \forall V_0 \in (1 - \varepsilon_0; 1), \forall f \in \mathcal{F}(V_0), \mathrm{Lip}(\partial\omega_{f,V_0}) + \mathrm{Per}(\omega_{f,V_0}) \le M. \quad (2.10) $$

We will use below [11, Lemma 3] to prove that $\omega_{f,V_0}$ actually has an analytic boundary, see Lemma 19; this is not necessary for the time being as we do not need analytic regularity.

We conclude with the following information about $\omega_{f,V_0}$.

**Lemma 14.** There exists $\varepsilon'_0 > 0$ such that, for any $V_0 \in (1 - \varepsilon'_0; 1)$, for any $f \in \mathcal{F}(V_0)$ the boundary $\partial\omega_{f,V_0}$ is connected.

**Lemma 14.** We argue by contradiction. Then there exists a sequence $\{V_k\}_{k \in \mathbb{N}}$ converging to 1 and, for any $k \in \mathbb{N}$, a function $f_k \in \mathcal{F}(V_k)$ such that $\partial\omega_{f_k,V_0}$ has at least two connected components $\Sigma_{k,1}, \Sigma_{k,2}$. We pick, for any $k \in \mathbb{N}$, two points $x_{k,i} \in \Sigma_{k,i}$ ($i = 1, 2$). First of all, for any closure point $y_i$ of the sequence $\{x_{k,i}\}_{k \in \mathbb{N}}$ we have $y_i \in \partial\Omega$ ($i = 1, 2$). This is simply due to the fact that $u_{f_k}(x_{k,i}) = \mu_{f_k,V_0} \underset{k \to \infty}{\to} 0$. Hence, from (2.4) we have $w_\Omega(y_i) = 0$, and we conclude by the maximum principle applied to the torsion function. Consequently, for $k$ large enough, $\partial\omega_{k,i} \subset \Omega(r'_0)$ where $\Omega(r'_0)$ is defined in (2.8) and $r'_0$ is small enough to ensure that $\partial\omega_{f,V_0}$ is a graph over $\partial\Omega$. From the connectedness of $\partial\Omega$, we reach the desired contradiction. $\square$

From the proof of this lemma we see that the only thing required is for $\partial\omega_{f,V_0}$ to be a graph over $\partial\Omega$ whence we may choose $\varepsilon_0 = \varepsilon'_0$ where $\varepsilon_0$ is given in Lemma 13.

## 2.2 Applications of the quantitative bathtub principle

The first step of the proof is to apply the quantitative bathtub principle [39, Proposition 26] (see also [14]). Before we state this result in the form that will be used let us recall the bathtub principle [34, Theorem 1.14]: let $u \in \mathcal{C}^0(\Omega)$ and let $V_0 \in (0; 1)$ be such that there exists a unique $\mu_{u,V_0}$ such that

$$ \mathrm{Vol}(\{u > \mu_{u,V_0}\}) = V_0 \mathrm{Vol}(\Omega) = \mathrm{Vol}(\{u \ge \mu_{u,V_0}\}). \quad (2.11) $$

If there exists a unique $\mu_{u,V_0}$ such that (2.11) is satisfied we define

$$ \omega_{u,V_0} := \{u > \mu_{u,V_0}\}. \quad (2.12) $$

Then:

$$ f_u := \mathbb{1}_{\omega_{u,V_0}} \text{ is the unique solution of } \max_{f \in \mathcal{F}(V_0)} \int_{\Omega} fu. \quad (2.13) $$

The goal of the quantitative bathtub principle is to quantify the optimality of $f_u$.

**Proposition 15.** [39, Proposition 26] Let $\alpha \in (0; 1)$. For any triplet of constants $M, P, \delta_1 > 0$ there exists a constant $c_1 > 0$ such that the following holds: for any function $u$ satisfying

1. there exists a unique $\mu_{u,V_0}$ that satisfies (2.11),

2. $\omega_{u,V_0}$ has a $\mathcal{C}^1$ boundary and

$$ \inf_{\partial\omega_{u,V_0}} \left| \frac{\partial u}{\partial \nu} \right| \ge \delta_1 > 0, $$

3. $\|u\|_{\mathcal{C}^{1,\alpha}(\Omega)} \le M$ and $\mathrm{Per}(\omega_{u,V_0}) \le P$,

then, defining $f_u := \mathbb{1}_{\omega_{u,V_0}}$

$$ \forall f \in \mathcal{F}(V_0), \int_{\Omega} fu \le \int_{\Omega} f_u u - c_1 \|f - f_u\|_{L^1(\Omega)}^2. \quad (2.14) $$
**Application to the thresholding algorithm.** The goal of this paragraph is to apply Proposition 15 to the sequence generated by the thresholding algorithm.

**Lemma 16.** Let $\varepsilon_0, \delta_0$ be in the conditions of Lemma 13. For any $V_0 \in (1 - \varepsilon_0; 1)$, for any $f_0 \in \mathcal{F}(V_0)$, the sequence $\{f_k\}_{k \in \mathbb{N}}$ is uniquely defined (in the sense that for any $k \in \mathbb{N}$ there holds $\mathrm{Vol}(\{u_{f_k} = \mu_{f_k, V_0}\}) = 0$) and there holds

$$ \sum_{k=0}^{\infty} \|f_{k+1} - f_k\|_{L^1(\Omega)}^2 < \infty. \qquad (2.15) $$

*Proof of Lemma 16.* Let $\varepsilon_0, \delta_0$ in the conditions of Lemma 13. We fix $V_0 \in (1 - \varepsilon_0; 1)$. The fact that for any $k \in \mathbb{N}$ we have

$$ \mathrm{Vol}(\{u_{f_k} = \mu_{f_k, V_0}\}) = 0 $$

is a conclusion of Lemma 13.

Since $\Omega$ is $\mathcal{C}^2$, by elliptic regularity, we know that for any $p \in [1; +\infty)$

$$ \sup_{k \in \mathbb{N}} \|u_{f_k}\|_{W^{2,p}(\Omega)} = M_p < \infty $$

whence from Sobolev embeddings [33, Theorem 12.5], fixing $s \in (0; 1)$, we have

$$ \sup_{k \in \mathbb{N}} \|u_{f_k}\|_{\mathcal{C}^{1,s}(\Omega)} =: M < \infty. $$

Furthermore,

$$ \inf_{k \in \mathbb{N}} \min_{\partial \omega_{u_{f_k}, V_0}} \left| \frac{\partial u_{f_k}}{\partial \nu} \right| \ge \delta_0 > 0. $$

Finally, estimate (2.10) allows us to apply Proposition 15 to deduce that there exists a constant $c_1 > 0$ such that

$$ \forall k \in \mathbb{N}, \int_{\Omega} u_k f_k \le \int_{\Omega} u_k f_{k+1} - c_1 \|f_{k+1} - f_k\|_{L^1(\Omega)}^2 $$

or, equivalently,

$$ \forall k \in \mathbb{N}, c_1 \|f_{k+1} - f_k\|_{L^1(\Omega)}^2 \le \int_{\Omega} u_k (f_{k+1} - f_k). \qquad (2.16) $$

From the concavity of $\mathcal{E}$ (Lemma 4) we also know that

$$ \forall k \in \mathbb{N}, \mathcal{E}(f_{k+1}) - \mathcal{E}(f_k) \le \dot{\mathcal{E}}(f_k)[f_{k+1} - f_k] = \int_{\Omega} u_k (f_k - f_{k+1}). $$

Thus we obtain

$$ \forall k \in \mathbb{N}, c_1 \|f_{k+1} - f_k\|_{L^1(\Omega)}^2 \le \mathcal{E}(f_k) - \mathcal{E}(f_{k+1}). \qquad (2.17) $$

However, there exists $\underline{E}$ such that

$$ \forall V_0 \in (0; 1), \forall f \in \mathcal{F}(V_0), \underline{E} \le \mathcal{E}(f). \qquad (2.18) $$

This is a direct consequence of the variational formulation of the Dirichlet energy (1.2). Summing the estimates (2.17) for $k \in \mathbb{N}$ we get

$$ c_1 \sum_{k=0}^{\infty} \|f_{k+1} - f_k\|_{L^1(\Omega)}^2 \le \mathcal{E}(f_0) - \underline{E} < \infty. $$

As $c_1 > 0$, the proof is concluded. $\square$
The purpose of this lemma is to give us more insight into the possible asymptotic behaviours of {$f_k$}$_{k \in \mathbb{N}}$. We detail this in the next paragraph.

Let us start with a basic lemma:

**Lemma 17.** Let $\varepsilon_0$ be in the conditions of Lemma 13. For any $V_0 \in (1-\varepsilon_0; 1)$, for any initialisation $f_0 \in \mathcal{F}(V_0)$, any $L^\infty$- closure point $f_\infty$ of the sequence {$f_k$}$_{k \in \mathbb{N}}$ generated by the thresholding algorithm is a bang-bang function:

$$ \exists E_\infty \subset \Omega, f_\infty = 1_{E_\infty}. \qquad (2.19) $$

Furthermore $E_\infty = \{u_{f_\infty} > \mu_{u_{f_\infty}, V_0}\}$, which is uniquely defined. In other words, $E_\infty$ is a critical set in the sense of Definition 7.

*Proof of Lemma 17.* We do not relabel the $L^\infty$-converging subsequence and write it {$f_k$}$_{k \in \mathbb{N}}$; its closure point is called $f_\infty$. We recall that $\varepsilon_0$ is given by Lemma 13. By elliptic regularity, for any $p \in [1; +\infty)$,

$$ \sup_{k \in \mathbb{N}} \|u_{f_k}\|_{W^{2,p}(\Omega)} < \infty. $$

By Sobolev embeddings, for any $\alpha \in (0; 1)$,

$$ u_{f_k} \underset{k \to \infty}{\stackrel{\alpha}{\rightleftharpoons}} u_{f_\infty}. $$

Define, for any $k \in \mathbb{N}$, $\omega_k := \{u_{f_k} \ge \mu_{u_{f_k}, V_0}\}$. From (2.10) we deduce that {$\omega_k$}$_{k \in \mathbb{N}}$ satisfies a uniform $\varepsilon$-cone property ([22, Remark 2.4.8]). From [22, Theorem 2.4.10] the sequence {$\omega_k$}$_{k \in \mathbb{N}}$ converges in the $L^1$-topology to a measurable subset $\omega_\infty$ of $\Omega$, with $\mathrm{Vol}(\omega_\infty) = V_0\mathrm{Vol}(\Omega)$; as $f_{k+1} = 1_{\omega_k}$ we deduce that $f_\infty = 1_{\omega_\infty}$. Moreover, still from [22, Theorem 2.4.10] the sequence {$\partial\omega_k$}$_{k \in \mathbb{N}}$ converges to $\partial\omega_\infty$ in the Hausdorff distance. Consequently, passing to the limit in $\partial\omega_k = \{u_{f_k} = \mu_{f_k, V_0}\}$ we deduce that $u_{f_\infty}$ is constant on $\partial\omega_\infty$. $\omega_\infty$ is thus necessarily the superlevel set of $u_{f_\infty}$ with volume $V_0$, which concludes the proof. $\square$

Using this result we can provide a simplification of the behaviour of the sequence {$f_k$}$_{k \in \mathbb{N}}$.

**Lemma 18.** Let $\varepsilon_0$ be chosen as in Lemma 13. For any $V_0 \in (1 - \varepsilon_0; 1)$, for any initialisation $f_0 \in \mathcal{F}(V_0)$, let {$f_k$}$_{k \in \mathbb{N}}$ be the sequence generated by the thresholding algorithm. Then:

1. Either {$f_k$}$_{k \in \mathbb{N}}$ has a unique weak $L^\infty$- closure point $f_\infty$, in which case $f_k \underset{k \to \infty}{\to} f_\infty$ strongly in $L^1(\Omega)$,

2. Or {$f_k$}$_{k \in \mathbb{N}}$ has an infinite number of closure points.

*Proof of Lemma 18.* Observe that if {$f_k$}$_{k \in \mathbb{N}}$ has a weak $L^\infty$- closure point $f_\infty$, Lemma 17 ensures that $f_\infty$ is an extreme point of $\mathcal{F}(V_0)$, whereby $f_\infty$ is actually a strong $L^1$ closure point of the sequence.

We prove that if the sequence {$f_k$}$_{k \in \mathbb{N}}$ has two distinct closure points $f_{\infty,1} \ne f_{\infty,2}$, then it has infinitely many closure points. First of all, from Lemma 16,

$$ \|f_{k+1} - f_k\|_{L^1(\Omega)} \underset{k \to \infty}{\to} 0. \qquad (2.20) $$

By Lemma 17 there exist $F_\infty^1, F_\infty^2 \subset \Omega$ such that

$$ f_{\infty,i} = 1_{F_{\infty,i}} \quad (i = 1, 2). $$
Let $\delta_1 \in (0; \text{dist}_{L^1}(F_{\infty,1}, F_{\infty,2})/4)$ and define

$$ \mathbb{N}_1 := \{k \in \mathbb{N} : \min_{i=1,2} (\|f_k - f_{\infty,i}\|_{L^1(\Omega)}) \ge \delta_1\}. $$

From (2.20), $\mathbb{N}_1$ is infinite. We denote it as

$$ \mathbb{N}_1 = \{i_0, \dots, i_k, \dots\}. $$

The sequence $\{f_{i_k}\}_{k \in \mathbb{N}}$ contains an $L^\infty$- converging subsequence, still denoted by $\{f_{i_k}\}_{k \in \mathbb{N}}$. Let $f_{\infty,3}$ be its closure point and let us show that

$$ \min_{i=1,2} (\|f_{\infty,3} - f_{\infty,i}\|_{L^1(\Omega)}) \ge \delta_1. $$

We only show

$$ \|f_{\infty,3} - f_{\infty,1}\|_{L^1(\Omega)} \ge \delta_1 \quad (2.21) $$

as the same proof would yield that $\|f_{\infty,3} - f_{\infty,2}\|_{L^1(\Omega)} \ge \delta_1$. To prove (2.21) we use in a crucial manner that $f_{\infty,1}$ is a bang-bang function (see Definition 1). Consider $h_k := f_{i_k} - f_{\infty,1}$. Then we have

$$ h_k = h_{k,+} \mathbb{1}_{F_{\infty,1}^c} - h_{k,-} \mathbb{1}_{F_{\infty,1}} $$

with $h_{k,\pm} \ge 0$. This is a simple consequence of the fact that $f_{\infty,1}$ is bang-bang. Furthermore, we have

$$ \int_{F_{\infty,1}^c} h_{k,+} + \int_{F_{\infty,1}} h_{k,-} \ge \delta. $$

We consider a weak closure point $h_{\infty,+}$ of $\{h_{k,+} \mathbb{1}_{F_{\infty,1}^c}\}_{k \in \mathbb{N}}$ and a weak closure point $h_{\infty,-}$ of $\{h_{k,-} \mathbb{1}_{F_{\infty,1}^c}\}_{k \in \mathbb{N}}$. We obtain

$$ \int_{F_{\infty,1}^c} h_{\infty,+} + \int_{F_{\infty,1}} h_{\infty,-} \ge \delta. $$

Furthermore, by linearity of the weak convergence we have

$$ f_{\infty,3} - f_{\infty,1} = h_{\infty,+} \mathbb{1}_{F_{\infty,1}^c} - h_{\infty,-} \mathbb{1}_{F_{\infty,1}}. $$

Consequently

$$ \|f_{\infty,3} - f_{\infty,1}\|_{L^1(\Omega)} \ge \delta, $$

whence the conclusion. Thus we deduce that (2.21) holds and, adapting the reasoning for $f_{\infty,2}$, we obtain that

$$ \min_{i=1,2} (\|f_{\infty,3} - f_{\infty,i}\|_{L^1(\Omega)}) \ge \delta. $$

By Lemma 17, $f_{\infty,3}$ is a bang-bang function, and $\{f_k\}_{k \in \mathbb{N}}$ has three distinct closure points. It then suffices to iterate the procedure to construct an infinite sequence of closure points. This concludes the proof. □

**The goal of the next sections.** To obtain the uniqueness of the closure point and the strong convergence to a local minimiser, we need to show that every critical point (in the sense of Definition 7) is actually a local minimiser and that these critical points are isolated. This is done using shape derivatives and diagonalisation of shape hessians.
## 2.3 Qualitative study of critical points

In this section we give an in-depth analysis of critical points (in the sense of Definition 7). We begin with the analyticity of the boundaries of critical sets.

**Regularity of critical sets.** We will be using a shape derivative formalism and compute second order shape derivatives at critical shapes in order to conclude as to their minimality and to obtain a full convergence result for the thresholding algorithm. Doing so requires some regularity (at least $\mathcal{C}^3$ for second-order shape derivatives) of the boundary of critical sets. It should be noted that usually this type of regularity is proved for minimal sets (i.e. for $E \subset \Omega$ such that $1_E$ is a solution of $(P_{\text{Dir}})$). A paradigmatic result is the regularity of minimal sets in two dimensions [13]. However, the hard part in proving this regularity is usually obtaining an estimate of the gradient of the state function on the boundary of the optimal set. Working with a large volume constraint allows to bypass this regularity problem, as gradient estimates are readily provided by Lemma 13, and thus enable us to apply [11, Lemma 3, Theorem 8].

**Lemma 19.** Let $\varepsilon_0 > 0$ be as in Lemma 13. Any critical set $E$ (in the sense of Definition 7) has a compact analytic boundary and, furthermore, is uniformly bounded in the $\mathcal{C}^2$ topology.

*Proof of Lemma 19.* This proof is an adaptation of [11, Proof of Theorem 8]. Let us first recall the following simpler version of [11, Lemma 3]:

**Lemma 20** (Lemma 3, [11]). Let $h: \mathbb{R} \to \mathbb{R}$ be a locally bounded function and $w \in \mathcal{C}^1(\Omega)$ be a solution of

$$ \Delta w = h(w) \text{ in } \Omega. $$

Let $x_0 \in \Omega$ be such that $\nabla w(x_0) \neq 0$. There exists a ball $\mathbb{B}(x_0; r)$ ($r > 0$) such that the set $\{x \in \mathbb{B}: w(x) = w(x_0)\}$ is an analytic hypersurface of $\mathbb{R}^d$.

The original version of [11, Lemma 3] also features a dependency of $h$ on $\nabla w$; this is not necessary here. Now consider a critical set $E \subset \Omega$. Let $u_{1_E}$ be the solution of (1.1) with $f = 1_E$ and let $\mu_E$ be such that

$$ E = \{u_{1_E} > \mu_E\}. \qquad (2.22) $$

The fact that (2.22) holds follows from Lemma 13. Defining the (locally bounded) function $h$ : $\eta \mapsto 1_{(\mu_E;+\infty)}(\eta)$, we thus have

$$ -\Delta u_{1_E} = h(u_{1_E}). $$

From Lemma 13, for any $x_0 \in \partial E = \{x: u_{1_E}(x) = \mu_E\}$, we have

$$ \nabla u_{1_E}(x_0) \neq 0. $$

Applying Lemma 20 we deduce that $\partial E$ is locally analytic. In other words, we may write

$$ \partial E = \bigcup_{x \in \partial E} (\mathbb{B}(x; r_x) \cap \partial E) $$

with $r_x > 0$ for any $x \in \partial E$ and $\mathbb{B}(x; r_x) \cap \partial E$ is an analytic hypersurface of $\mathbb{R}^n$. As $\partial E$ is compact by continuity of $u_{1_E}$, we may extract a finite covering of $\partial E$ and hence conclude that $\partial E$ is a compact analytic hypersurface. To obtain the uniform $\mathcal{C}^2$ bounds it suffices to conclude as in [11, Lemma 3, Theorem 8], by using (2.10)$^2$. $\square$

<footnote>$^2$In the proof of [11, Lemma 3] the analytic bounds obtained only depend on the regularity of the local change of coordinates, which here is uniform in $E$ from Lemma 13 and the implicit function theorem.</footnote>
**Preliminary considerations about shape derivatives at critical shapes.** We first identify the functional $\mathcal{E}$ with a shape functional $\mathcal{F}$, by defining

$$ \mathcal{F}: \Omega \supset E \mapsto \mathcal{E}(1_E). $$

For any compact $\mathcal{C}^3$ subset $E$ of $\Omega$ and for any compactly supported vector field $\Phi \in W^{2,\infty}(\Omega; \mathbb{R}^d)$, we can define, for any $t \in (-1; 1)$ small enough, the function

$$ e_{E,\Phi}(t) := \mathcal{F}((\mathrm{Id} + t\Phi)E). $$

Provided they exists, the first (resp. second) order derivative of $\mathcal{F}$ at $E$ in the direction $\Phi$ is defined as

$$ \mathcal{F}'(E)[\Phi] = e'_{E,\Phi}(0), \text{ resp. } \mathcal{F}''(E)[\Phi, \Phi] = e''_{E,\Phi}(0). $$

We first check that $\mathcal{E}$ is indeed twice differentiable in the sense of shapes.

**Lemma 21.** For any $\mathcal{C}^3$ subset $E$ of $\Omega$, $\mathcal{F}$ is twice differentiable at $E$ in the direction of any compactly supported vector-field $\Phi \in W^{2,\infty}(\Omega; \mathbb{R}^d)$.

As this lemma is proved by a direct adaptation of [41] we omit it here and focus in a subsequent paragraph on the computation of first and second order shape derivatives.

**First order shape derivative and definition of the Lagrange multiplier.** In terms of optimality conditions, since we are working with volume constraints, we have to work with the shape derivative of the volume functional. We recall [22, Section 5.9.3, first example] that the map

$$ \Omega \supset E \mapsto \mathrm{Vol}(E) $$

is twice shape differentiable (i.e. differentiable in the direction of any compactly supported vector field $\Phi \in W^{2,\infty}(\Omega)$) at any $\mathcal{C}^3$ domain $E$ and that, for any compactly supported vector field $\Phi$, there holds

$$ \mathrm{Vol}'(E)[\Phi] = \int_{\partial E} \langle \Phi, \nu_E \rangle, \mathrm{Vol}''(E)[\Phi, \Phi] = \int_{\partial E} \mathcal{H}_E \langle \Phi, \nu \rangle^2 \quad (2.23) $$

where $\nu_E$ is the normal vector on $\partial E$ and $\mathcal{H}_E$ is the mean curvature of $\partial E$.

For semantical convenience, we introduce the following definition:

**Definition 22** (Critical shape). A $\mathcal{C}^2$ shape $E$ is a critical shape for $\mathcal{F}$ if, for any compactly supported $\Phi \in W^{2,\infty}(\Omega)$,

$$ \int_{\partial E} \langle \Phi, \nu_E \rangle = 0 \Rightarrow \mathcal{F}'(E)[\Phi] = 0. $$

In order to exploit second order optimality conditions it is more convenient to use a Lagrange multiplier. If a shape $E$ is critical, then the condition given in Definition 22 rewrites as: there exists a constant $\mu_E$ such that, for any compactly supported vector field $\Phi \in W^{2,\infty}(\Omega; \mathbb{R}^d)$,

$$ (\mathcal{F} + \mu_E \mathrm{Vol})'(E)[\Phi] = 0. \quad (2.24) $$

Computing this Lagrange multiplier is an important step in our subsequent analysis; to do it we need an expression of the first order derivative.

**Lemma 23.** For any $\mathcal{C}^3$ shape, for any compactly supported vector field $\Phi \in W^{2,\infty}(\Omega; \mathbb{R}^d)$,

$$ \mathcal{F}'(E)[\Phi] = - \int_{\partial E} u_{1_E} \langle \Phi, \nu_E \rangle. $$
Furthermore, the shape derivative of the map $E \mapsto u_{1_E}$ at $E$ in the direction $\Phi$ is the unique solution of

$$ \begin{cases} -\Delta u'_{\Phi} = 0 & \text{in } \Omega, \\ [\partial_{\nu} u'_{\Phi}] = -\langle \Phi, \nu_E \rangle & \text{on } \partial E \end{cases} \qquad (2.25) $$

where $[\cdot]$ denotes the jump of a function across a hypersurface.

**Proof of Lemma 23.** To prove this lemma we first need to compute the shape derivative of the map $u \mapsto u_{1_E}$. These computations were already carried out (for the shape derivative of the first eigenvalue of the operator $-\Delta - \mathbb{1}_E$) in details in [39], but we sketch them in the present case for the sake of completeness. Fix $E$ and a compactly supported vector field $\Phi$. Define, for any $t$ such that $|t|$ small enough, $E_t := \mathbb{1}_{(\text{Id}+t\Phi)_E}$ and $u_t$ as the solution of (1.1) with $f = \mathbb{1}_{E_t}$. The weak formulation of the equation on $u_t$ reads:

$$ \forall v \in W_0^{1,2}(\Omega), \int_{\Omega} \langle \nabla u_t, \nabla v \rangle = \int_{E_t} v. \qquad (2.26) $$

Taking the derivative of this formulation at $t=0$, it appears that the shape derivative $u'_{\Phi}$ satisfies the equation

$$ \forall v \in W_0^{1,2}(\Omega), \int_{\Omega} \langle \nabla u'_{\Phi}, \nabla v \rangle = \int_{\partial E} v \langle \Phi, \nu_E \rangle. \qquad (2.27) $$

Alternatively, this rewrites as: $u'_{\Phi}$ is the unique solution of the equation

$$ \begin{cases} -\Delta u'_{\Phi} = 0 & \text{in } \Omega, \\ [\partial_{\nu} u'_{\Phi}] = -\langle \Phi, \nu_E \rangle & \text{on } \partial E \end{cases} $$

exactly as claimed in the statement of the theorem. Now, going back to the definition of $\mathcal{F}$ we obtain

$$ \mathcal{F}'(E)[\Phi] = \int_{\Omega} \langle \nabla u_{1_E}, \nabla u'_{\Phi} \rangle - \int_{E} u'_{\Phi} - \int_{\partial E} u_{1_E} \langle \Phi, \nu_E \rangle. $$

Using the weak formulation of (1.1) this gives

$$ \mathcal{F}'(E)[\Phi] = - \int_{\partial E} u_{1_E} \langle \Phi, \nu_E \rangle. $$

□

This expression allows us to prove the following result

**Lemma 24.** A shape $E \subset \Omega$ is a critical set in the sense of Definition 7 if, and only if it is a critical shape in the sense of Definition 22.

**Proof of Lemma 24.** Let us first notice that from Lemma 19, critical shapes in the sense of Definition 7 are analytic and thus in particular $\mathcal{C}^2$.

Consider now a critical point $E$ in the sense of Definition 7. By definition there exists a constant $\mu_E$ such that $u_{1_E} = \mu_E$ on $\partial E$. As $E$ is $\mathcal{C}^2$ from Lemma 19 we can legitimately compute the first order shape derivative of $\mathcal{F}$ at $E$ in the direction of a fixed compactly supported vector field $\Phi$. By Lemma 23

$$ \mathcal{F}'(E)[\Phi] = - \int_{\partial E} u_{1_E} \langle \Phi, \nu_E \rangle = - \mu_E \int_{\partial E} \langle \Phi, \nu_E \rangle. $$

Hence, if $\int_{\partial E} \langle \Phi, \nu_E \rangle = 0$ we have $\mathcal{F}'(E)[\Phi] = 0$ and thus $E$ is a critical shape in the sense of Definition 22.
Conversely, consider a shape E of class $\mathcal{C}^3$ that is critical in the sense of Definition 22. As for any compactly supported vector field $\Phi$ such that $\int_{\partial E} \langle \Phi, \nu_E \rangle = 0$ we must have

$$ \int_{\partial E} u_{1_E} \langle \Phi, \nu_E \rangle = 0 $$

the fundamental lemma of the calculus of variations implies the existence of a constant $\mu_E > 0$ such that

$$ u_{1_E} = \mu_E \text{ on } \partial E. $$

From the maximum principle applied to (1.1) we already know that $u_{1_E} > 0$ in $\Omega$, hence $\mu_E > 0$. Consequently $\partial E \cap \partial \Omega = \emptyset$ and $u_{1_E}$ is a solution of

$$ \left\{ \begin{array}{ll} -\Delta u_{1_E} = 1 & \text{in } E, \\ u_E = \mu_E & \text{on } \partial E \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{ll} -\Delta u_{1_E} = 0 & \text{in } \Omega \setminus \bar{E}, \\ u_{1_E} = \mu_E & \text{on } \partial E, \\ u_{1_E} = 0 & \text{on } \partial \Omega. \end{array} \right. \qquad (2.28) $$

As $\partial E \cap \partial \Omega = \emptyset$, both E and $\Omega \setminus \bar{E}$ are smooth open sets. From the strong maximum principle we deduce that $u_{1_E} > \mu_E$ in E and that $u_{1_E} < \mu_E$ in $\Omega \setminus \bar{E}$. Thus E is a critical set in the sense of Definition 7. $\square$

From this analysis we deduce the following:

**Lemma 25.** Assume E is a critical shape in the sense of Definition 22. The associated Lagrange multiplier is $\mu_E$ where $E = \{u_{1_E} > \mu_E\}$ in the sense that

$$ \forall \Phi \in W^{2,\infty}(\Omega; \mathbb{R}^d), \Phi \text{ compactly supported}, (\mathcal{F} + \mu_E \text{Vol})'(E)[\Phi] = 0. $$

To alleviate notations we introduce a notation for the shape Lagrangian:

**Definition 26 (Shape Lagrangian).** Let E be a critical shape for $\mathcal{F}$. The associated shape Lagrangian is

$$ \mathcal{L}_E := \mathcal{F} + \mu_E \text{Vol} $$

$$ \text{where } E = \{u_{1_E} > \mu_E\}. $$

With the above ingredients at hand, we now move to second-order shape derivatives.

**First computations and elements for shape hessians.** In the subsequent paragraphs we use indifferently the expressions "shape hessians" and "second order shape derivatives". We begin with an expression of the shape hessian at a critical shape E in the direction $\Phi$.

**Lemma 27.** Let E be a critical shape in the sense of Definition 22 and let $\Phi \in W^{2,\infty}(\Omega; \mathbb{R}^d)$ be a compactly supported vector field. The second-order shape derivative of the shape Lagrangian $\mathcal{L}_E$ (defined in Definition 26) is given by the expression

$$ \mathcal{L}_E''(E)[\Phi, \Phi] = - \int_{\partial E} u'_{\Phi} \langle \Phi, \nu_E \rangle - \int_{\partial E} \frac{\partial u_E}{\partial \nu} \langle \Phi, \nu_E \rangle^2 $$

where $u'_{\Phi}$ solves (2.25).
Proof of Lemma 27. First of all, since E is a critical shape, it suffices to work with vector fields that are normal to ∂E (see [22, Remark on page 246]). We recall that the formula of Hadamard gives

$$ \left. \frac{d}{dt} \right|_{t=0} \int_{\Gamma(t)} f(t) = \int_{\Gamma(0)} f'(0) + \int_{\Gamma(0)} \left( \mathcal{H}_{\Gamma(0)} f + \frac{\partial f}{\partial \nu} \right) \langle \Phi, \nu \rangle $$

where $\mathcal{H}$ denotes the mean curvature of a hypersurface.

Let us start from the fact that, at a given shape E we have (Lemma 23)

$$ \mathcal{F}'(E)[\Phi] = - \int_{\partial E} u_{1_E} \langle \Phi, \nu_E \rangle. $$

Applying the formula of Hadamard we obtain

$$ \mathcal{F}''(E)[\Phi, \Phi] = - \int_{\partial E} u'_{\Phi} \langle \Phi, \nu \rangle - \int_{\partial E} \left( \mathcal{H}_E u_{1_E} + \frac{\partial u_{1_E}}{\partial \nu} \right) \langle \Phi, \nu_E \rangle^2. $$

Taking into account the definition of the Lagrange multiplier (Lemma 25) and (2.23) we finally derive the following expression for the second-order shape derivative of the Lagrangian (see Definition 26)

$$ \mathcal{L}_E''(E)[\Phi, \Phi] = - \int_{\partial E} u'_{\Phi} \langle \Phi, \nu \rangle - \int_{\partial E} \frac{\partial u_{1_E}}{\partial \nu} \langle \Phi, \nu_E \rangle^2. \qquad (2.29) $$

□

**The question of coercivity of shape Lagrangians.** The goal of upcoming sections is to provide coercivity results on shape Lagrangians at critical shapes. By coercivity we mean not only that the shape hessian is positive, i.e. that for any compactly supported vector field $\Phi \in W^{2,\infty}(\Omega; \mathbb{R}^d)$ we have

$$ \mathcal{L}_E''(E)[\Phi, \Phi] > 0 \text{ whenever } \langle \Phi, \nu_E \rangle \neq 0 \text{ on } \partial E $$

but that we can even quantify this positivity, in the sense that there exists a constant $c_E > 0$ (independent of $\Phi$) such that

$$ \mathcal{L}_E''(E)[\Phi, \Phi] \geq c_E \|\langle \Phi, \nu_E \rangle\|_{L^2(\partial E)}^2. $$

This leads to introducing the following definition:

**Definition 28 (L²-stability).** A critical shape E is said to be L²-stable if there exists a constant $c_E > 0$ such that for any compactly supported vector field $\Phi \in W^{2,\infty}(\Omega; \mathbb{R}^d)$,

$$ \mathcal{L}_E''(E)[\Phi, \Phi] \geq c_E \|\langle \Phi, \nu_E \rangle\|_{L^2(\partial E)}^2. $$

It should be noted that we can not expect a better coercivity norm than $L^2(\partial E)$; when diagonalising the shape hessian we actually will prove the optimality of this norm. While weaker than the usual $H^{\frac{1}{2}}(\partial E)$ norm [16, 17, 23] this is still enough to obtain local quantitative inequalities [39].

Our goal is henceforth to prove the following coercivity result:

**Proposition 29.** There exists $\varepsilon_1 > 0$, $\varepsilon_1 \le \varepsilon_0$, $c > 0$ such that, for any $V_0 \in (1 - \varepsilon_1; 1)$, for any critical set E (in the sense of Definition 22), for any compactly supported vector field $\Phi \in W^{2,\infty}(\Omega; \mathbb{R}^d)$,

$$ \mathcal{L}_E''(E)[\Phi, \Phi] \geq c \|\langle \Phi, \nu \rangle\|_{L^2(\partial E)}^2. $$
Before we prove this proposition we need to introduce the diagonalisation basis. We consider a fixed critical shape $E \subset \Omega$ and work with $V_0 \in (1 - \varepsilon_0; 1)$, where $\varepsilon_0$ is given by Lemma 13. Heuristically, given that $u'_\Phi$ satisfies (2.25), it is natural to solve the following eigenvalue problem: find $(\lambda, \phi) \in \mathbb{R} \times W^{1,2}(\Omega)$ such that

$$ egin{cases} -\Delta\phi = 0 & \text{in } E \cup (E)^c \\ [[\partial_\nu\phi]]_E = -\lambda\rho\phi & \text{on } \partial E, \\ \int_{\partial E} \phi^2 > 0, \\ \phi \in W_0^{1,2}(\Omega) \end{cases} \qquad (2.30) $$

where $\rho$ is a weight to be determined. Let us simply consider the eigenvalue problem (2.30) with a weight $\rho \in L^\infty(\partial E; [\delta, 1/\delta])$ for some small $\delta > 0$. To study this spectral problem we introduce the weighted space

$$ L_{\rho}^{2}(\partial E) := \{f \in L^{2}(\partial E), \int_{\partial E} \rho f^{2} < \infty\} = L^{2}(\partial E). $$

The last equality comes from the fact that the weight $\rho$ is uniformly bounded. The only difference here is of course the scalar product; we work with

$$ \langle \cdot, \cdot \rangle_{\rho} : (f,g) \mapsto \int_{\partial E} \rho fg. $$

Under these assumptions it is fairly standard to obtain the existence of a spectral basis associated with (2.30). Define the operator $T_\rho : L_\rho^2(\partial E) \to L^2(\partial E)$ as follows: for any $f \in L_\rho^2(\partial E)$ let $v_f \in W^{1,2}(\Omega)$ be the unique solution to

$$ egin{cases} -\Delta v_f = 0 & \text{in } \Omega, \\ [[\partial_\nu v_f]] = -\rho f & \text{on } \partial E, \\ v_f \in W_0^{1,2}(\Omega). \end{cases} \qquad (2.31) $$

The fact that for any $f \in L_\rho^2(\partial E)$ the function $v_f$ exists and is unique simply follows by minimising the functional

$$ W_{\rho,f} : W_0^{1,2}(\Omega) \ni v \mapsto \frac{1}{2} \int_{\Omega} |\nabla v|^2 - \int_{\partial E} \rho v f. $$

Let $\text{Tr}_E : W^{1,2}(\Omega) \to L^2(\partial E)$ be the trace operator. As $E$ is analytic (Lemma 19) this operator is well-defined. Set

$$ T_{\rho}(f) := \text{Tr}_{E}(v_{f}). \qquad (2.32) $$

In particular,

$$ \forall f \in L_{\rho}^{2}(\partial E), T_{\rho}(f) \in L_{\rho}^{2}(\partial E). $$

Clearly, $T_\rho$ is symmetric, positive and compact. In particular, by the spectral decomposition theorem, there exists a sequence of eigen-elements $\{(\sigma_{k,\rho}, \psi_{k,\rho})\}_{k \in \mathbb{N}} \in (\mathbb{R}_+^* \times L_\rho^2(\partial E))^\mathbb{N}$ such that

$$ \sigma_{k,\rho} \text{ is non-increasing in } k \text{ and } \sigma_{k,\rho} \underset{k \to \infty}{\longrightarrow} 0, $$

$$ \{\psi_{k,\rho}\}_{k \in \mathbb{N}} \text{ is an orthonormal basis of } L_{\rho}^{2}(\partial E) $$

and

$$ \forall k \in \mathbb{N}, \quad T_{\rho}(\psi_{k,\rho}) = \sigma_{k,\rho}\psi_{k,\rho}. $$
We set $\phi_{k,\rho} := v_{\psi_{k,\rho}}$ to extend it to a function on $\Omega$. Hence,

$$ egin{cases} -\Delta\phi_{k,\rho} = 0 & \text{in } \Omega, \\ \left[\partial_\nu\phi_{k,\rho}\right] = -\rho\psi_{k,\rho} = -\frac{1}{\sigma_{k,\rho}}\rho\phi_{k,\rho} & \text{on } \partial E, \\ \phi_{k,\rho} = 0 & \text{on } \partial\Omega. \end{cases} \tag{2.33} $$

We thus obtain, defining

$$ \lambda_{k,\rho} := \frac{1}{\sigma_{k,\rho}} $$

the system

$$ \forall k \in \mathbb{N}^*, \quad \left\{ \begin{array}{ll} -\Delta\phi_{k,\rho} = 0 & \text{in } \Omega, \\ \left[\partial_\nu\phi_{k,\rho}\right] = -\lambda_{k,\rho}\rho\phi_{k,\rho} & \text{on } \partial E, \\ \phi_{k,\rho} = 0 & \text{on } \partial\Omega \end{array} \right. \\ \text{and } \forall k, k' \in \mathbb{N}^*, \int_{\partial E} \rho\phi_{k,\rho}\phi_{k',\rho} = \delta_{k,k'}. \tag{2.34} $$

The family $\{\phi_{0,\rho}\} \cup \{\phi_{k,\rho}\}_{k \in \mathbb{N}^*} = \{\phi_{k,\rho}\}_{k \in \mathbb{N}}$ is an orthonormal basis of $L^2_\rho(\partial E)$ for the scalar product $\langle \cdot, \cdot \rangle_\rho$.

From the Courant-Fisher principle we have, furthermore, the following characterisation of the first eigenvalue:

$$ \lambda_{0,\rho} = \min_{u \in W_0^{1,2}(\Omega), \int_{\partial E} u^2 > 0} \frac{\int_{\Omega} |\nabla u|^2}{\int_{\partial E} \rho u^2}. \tag{2.35} $$

**Diagonalisation of $\mathcal{L}''_E(E)$**. We can now turn back to the expression of the second-order shape derivative of the Lagragian (2.29). We want to determine the weight $\rho$. We use the following notational convention: for any $\Phi \in W^{2,\infty}(\Omega; \mathbb{R}^d)$,

$$ \varphi := \langle \Phi, \nu \rangle, \tilde{\varphi} := \frac{1}{\rho}\varphi. $$

The first-order shape derivative $u'_\Phi$ satisfies

$$ \left\{ \begin{array}{ll} -\Delta u'_{\Phi} = 0 & \text{in } \Omega, \\ \left[ \frac{\partial u'_{\Phi}}{\partial \nu} \right] = -\rho \tilde{\varphi} & \text{on } \partial E, \\ u'_{\Phi} = 0 & \text{on } \partial \Omega. \end{array} \right. \tag{2.36} $$

We decompose $\tilde{\varphi}$ in the basis $\{\phi_{k,\rho}\}_{k \in \mathbb{N}}$ as

$$ \tilde{\varphi} = \sum_{k=0}^{\infty} \alpha_{k,\rho} \phi_{k,\rho}. \tag{2.37} $$

Thus,

$$ u'_{\Phi} = \sum_{k=0}^{\infty} \frac{\alpha_{k,\rho}}{\lambda_{k,\rho}} \phi_{k,\rho}. $$

As a conclusion

$$ \mathcal{L}''_E(E)[\Phi, \Phi] = - \int_{\partial E} u'_{\Phi} \langle \Phi, \nu_E \rangle - \int_{\partial E} \partial_{\nu} u_{1_E} \langle \Phi, \nu_E \rangle^2 $$
$$= - \int_{\partial E} \rho \tilde{\varphi} u'_{\Phi} - \int_{\partial E} \rho^2 \partial_{\nu} u_{1_E} (\tilde{\varphi})^2 .$$

Let us choose

$$\rho_E = -\frac{1}{\partial_{\nu} u_{1_E}}.$$

From Lemma 13, we have $\rho_E \in L^\infty([\delta; 1/\delta])$ for $\delta$ small enough. Finally, this yields

$$\begin{aligned} \mathcal{L}_E''(E)[\Phi, \Phi] &= - \int_{\partial E} u'_{\Phi} \langle \Phi, \nu_E \rangle - \int_{\partial E} \partial_{\nu} u_{1_E} \langle \Phi, \nu_E \rangle^2 \\ &= - \int_{\partial E} \rho_E \tilde{\varphi} u'_{\Phi} - \int_{\partial E} \rho_E^2 \partial_{\nu} u_{1_E} (\tilde{\varphi})^2 \\ &= \sum_{k=0}^{\infty} \alpha_{k, \rho_E}^2 \left( 1 - \frac{1}{\lambda_{k, \rho_E}} \right). \end{aligned}$$

In particular, given that the sequence $\{\lambda_{k, \rho_E}\}_{k \in \mathbb{N}}$ is non-decreasing, we have the estimate from below:

$$\mathcal{L}_E''(E)[\Phi, \Phi] \geq \left(1 - \frac{1}{\lambda_{0, \rho_E}}\right) \sum_{k=0}^{\infty} \alpha_{k, \rho_E}^2 = \left(1 - \frac{1}{\lambda_{0, \rho_E}}\right) \|\langle \Phi, \nu_E \rangle\|_{L^2(\partial E)}^2. \quad (2.38)$$

Hence we have the following sufficient condition for the stability of a critical shape:

**Lemma 30.** If $\lambda_{0, \rho_E} > 1$, then E is stable in the sense of Definition 28.

To prove Proposition 29 we study the asymptotic behaviour of $\lambda_{0, \rho_E}$ as $V_0 \to 1$.

**Asymptotic behaviour of $\lambda_{0, \rho_E}$ as $V_0 \to 1$**. We prove here the following lemma:

**Lemma 31.** There holds

$$\min_{E \text{ critical shape for the volume constraint } V_0} \lambda_{0,E} \xrightarrow{V_0 \to 1} \infty.$$

**Remark 32** (Heuristic approach to Lemma 31). Assume $\Omega = (-1; 1)$ and $V_0 = 1 - \varepsilon$, where $\varepsilon > 0$ is a small parameter. In this case, we may use the Schwarz rearrangement to obtain that the unique solution of $(P_{\text{Dir}})$ is $f^* = 1_{I_\varepsilon}$ where

$$I_{\varepsilon} = (-1 + \varepsilon/2; 1 - \varepsilon/2).$$

Let $u_\varepsilon$ be the associated solution of (1.1). Clearly,

$$u_{\varepsilon} \xrightarrow{\varepsilon \to 0} u_0$$

in $\mathcal{C}^1(\Omega)$, where $u_0$ is the solution of (1.1) with $f \equiv 1$. In particular,

$$|u'_{\varepsilon}(\pm 1 \mp \varepsilon/2)| = a_{\varepsilon} \xrightarrow{\varepsilon \to 0} 1.$$

Now let us study the lowest eigenvalue $\lambda_{0, \rho_E}$ which, for the sake of notational convenience, we rewrite as $\lambda_0(\varepsilon)$. But notice the following thing: if we define

$$\varphi_{0,\varepsilon}: x \mapsto \begin{cases} 1 & \text{if } |x| \le 1 - \varepsilon/2, \\ \frac{2}{\varepsilon a_{\varepsilon}}(1 - |x|) & \text{if } 1 - \varepsilon/2 \le |x| \le 1 \end{cases}.$$
it appears that $\varphi_{0,\varepsilon}$ is an eigenfunction associated with the eigenvalue $\frac{2}{\varepsilon}$. As it has a constant sign, and as the eigenfunction associated with the lowest eigenvalue is the only one that has a constant sign, we deduce that

$$ \lambda_{0,\rho_{I_{\varepsilon}}}(\varepsilon) = \frac{2}{\varepsilon}, $$

which indeed goes to $\infty$ as $\varepsilon \to 0^+$.

The situation is of course more complicated in higher dimensions. Here is a roadmap to prove the same type of results. We already observed in the proof of Lemma 13 that $u_f \xrightarrow{V_0 \to 1} w$ uniformly in $f$, as $\varepsilon \to 0$, where $w$ is the torsion function of the set $\Omega$. As $\Omega$ has a regular boundary, the Hopf maximum principle entails that

$$ \inf_{\partial \Omega} \left( -\frac{\partial w}{\partial \nu} \right) = \delta(\Omega) > 0. $$

Now consider the set $X_{\varepsilon}^*$ of critical shapes. From the definition of criticality, it follows that, for any $E^* \in X^*$, there exists $\mu(E^*)$ such that $E^* = \{u_{1E^*} > \mu(E^*)\}$. We want to prove that

$$ \liminf_{\varepsilon \to 0^+} \inf_{E^* \in X_{\varepsilon}^*} \lambda_{0,\rho_{E^*}} = +\infty. $$

We argue by contradiction and assume this is not the case. This gives a sequence $\varepsilon \to 0^+$ and a sequence $\{E_{\varepsilon}\}_{\varepsilon>0} \in \prod_{\varepsilon>0} X_{\varepsilon}^*$ such that

$$ \sup_{\varepsilon} \lambda_{0,\rho_{E^*}} \le M_0 < \infty. $$

By Lemma 13, $E_{\varepsilon}^*$ converges in $\mathcal{C}^{1,\alpha}$ to $\Omega$ (for any $\alpha < 1$). Now observe that, in fact, by a standard comparison principle, we may be satisfied to prove that $\lambda(\varepsilon) \to \infty$ where $\lambda(\varepsilon)$ is defined as

$$ \lambda(\varepsilon) := \inf_{\varphi \in W_0^{1,2}(\Omega)} \frac{\int_{\Omega} |\nabla \varphi|^2}{\int_{\partial E_{\varepsilon}} \varphi^2}. $$

This follows from the lower bound on $|\partial_{\nu} u_{\varepsilon}|$ given by Lemma 13. The idea is then that, if the eigenvalues are uniformly bounded, then we should have on the one-hand (up to renormalisation)

$$ \int_{E_{\varepsilon}} \varphi_{\varepsilon}^2 = 1 $$

and, on the other hand, we should have enough regularity on $\varphi_{\varepsilon}$ to write

$$ \int_{E_{\varepsilon}} \varphi_{\varepsilon}^2 \sim \int_{\partial \Omega} \varphi_{\varepsilon}^2 = 0. $$

This would be an obvious contradiction.

**Proof of Lemma 31.** We first observe that from Lemma 13 and the Rayleigh quotient formulation (2.35) of $\lambda_{0,\rho_E}$ we have, for $V_0 \in (1 - \varepsilon_0; \varepsilon)$,

$$ \lambda_{0,\rho_E} \ge \delta_0 \min_{u \in W_0^{1,2}(\Omega), \int_{\partial E} u^2 \ne 0} \frac{\int_{\Omega} |\nabla u|^2}{\int_{\partial E} u^2}. $$

Now, argue by contradiction and assume that there exists a sequence $\varepsilon \to 0$, a constant $M_0 > 0$, and, for any $\varepsilon > 0$, a critical set $E_{\varepsilon}$ such that $1_{E_{\varepsilon}} \in \mathcal{F}(V_0)$ with $V_0 \in (1 - \varepsilon; 1)$ that satisfies

$$ \lambda_{0,\rho_{E_{\varepsilon}}} \le M_0. $$
Define, for any $\varepsilon > 0$,

$$ \lambda(\varepsilon) := \min_{u \in W_0^{1,2}(\Omega), \int_{\partial E_\varepsilon} u^2 \neq 0} \frac{\int_{\Omega} |\nabla u|^2}{\int_{\partial E_\varepsilon} u^2}. $$

As already observed we have $\delta_0\lambda(\varepsilon) \le \lambda_{0,\rho_{E_\varepsilon}}$ so that the sequence $\{\lambda(\varepsilon)\}_{\varepsilon \to 0}$ is uniformly bounded. In other words we have

$$ 0 < \lambda(\varepsilon) \le M_1 < \infty \qquad (2.39) $$

for some suitable constant $M_1$. We define, for any $\varepsilon > 0$, $\varphi_\varepsilon$ as the first eigenfunction associated with $\lambda(\varepsilon)$; in other words, $\varphi_\varepsilon$ solves the equation

$$ \left\{ \begin{array}{ll} -\Delta\varphi_\varepsilon = 0 & \text{in } \Omega, \\ \left[ \frac{\partial\varphi_\varepsilon}{\partial\nu} \right] = -\lambda(\varepsilon)\varphi_\varepsilon & \text{on } \partial E_\varepsilon, \\ \varphi_\varepsilon \in W_0^{1,2}(\Omega), \\ \int_{\partial E_\varepsilon} \varphi_\varepsilon^2 = 1. \end{array} \right. \qquad (2.40) $$

Let us prove that if (2.39) holds then:

$$ \forall \beta \in (0; 1), \exists c_\beta, \sup_{\varepsilon>0} \|\varphi_\varepsilon\|_{\mathcal{L}^{0,\beta}} \le c_\beta. \qquad (2.41) $$

First of all, from the uniform regularity estimates of Lemma 13 and from [33, Theorem 18.34], for any $p \in (1; +\infty)$, there exists a constant $C_p$ such that

$$ \forall \varepsilon > 0, \|\mathrm{Tr}_{E_\varepsilon}\|_{\mathcal{L}(W_0^{1,p}(\Omega); W^{1-\frac{1}{p},p}(\partial E_\varepsilon))} \le C_p. $$

The idea is to use a bootstrap argument as well as the duality method of Stampacchia. Let show how to initialise the bootstrap. As

$$ \int_{\Omega} |\nabla \varphi_{\varepsilon}|^2 = \lambda(\varepsilon) $$

it follows that

$$ \sup_{\varepsilon} \|\varphi_{\varepsilon}\|_{W^{\frac{1}{2},2}(\partial E_{\varepsilon})} \le C_0 M_1 $$

where $C_0$ is a trace embedding constant which only depends on $\sup_\varepsilon \mathrm{Lip}(E_\varepsilon)$. Let $2' = \frac{2d}{2d-3}$ be the conjugate Sobolev exponent of 2 on the boundary. Then we have

$$ \sup_{\varepsilon} \|\varphi_{\varepsilon}\|_{L^{2'}(E_{\varepsilon})} < C_1. $$

By a bootstrapping argument, we see that it suffices to prove the following regularity Lemma:

**Lemma 33.** There exists a constant $C$ that only depends on $\mathrm{Lip}(E)$ and the covering number of $E$ such that the following holds: assume $\psi$ solves

$$ \left\{ \begin{array}{ll} -\Delta\psi = 0, & \partial E_\psi \\ \left[ \partial_\nu \psi \right] = -p \in L^q(\partial E) & \partial E_\psi \end{array} \right. $$

then

$$ \|\psi\|_{W^{1,q}(\Omega)} \le C \|p\|_{L^q(\partial E)}. $$
Proof. Recall that, by standard elliptic estimates (see for instance [1, Theorem 1.1]), for any smooth vector field $F \in \mathcal{C}^\infty(\Omega)$, the solution $v_F$ of

$$ \begin{cases} -\Delta v_F = \nabla \cdot F & \text{in } \Omega, \\ v_F = 0 & \text{on } \partial\Omega \end{cases} $$

satisfies

$$ \|v_F\|_{W^{1,r}(\Omega)} \le D_r \|F\|_{L^q(\partial E)} $$

for some constant $D_r = D(r, \Omega)$.

Now, consider a smooth vector field $F$, then there holds

$$ \int_{\Omega} \langle \nabla \psi, F \rangle = \int_{\Omega} \langle \nabla \psi, \nabla v_F \rangle = \int_{\partial E} v_F p \le C_E \|v_F\|_{W^{1,r}(\Omega)} \|p\|_{L^q(\partial E)} $$

Here, $C$ is the trace embedding constant; it only depends on Lip($E$) and on the covering number of $\partial E$. Thus we get

$$ \int_{\Omega} \langle \nabla \psi, F \rangle \le C_E D_r \|F\|_{L^r(\Omega)} \|p\|_{L^q(\partial E)} $$

Thus $\psi \in W^{1,q}(\Omega)$ and

$$ \|\psi\|_{W^{1,q}(\Omega)} \le C_E D_r \|p\|_{L^q(\partial E)} $$

□

Consequently we deduce by a bootstrap argument that, if (2.39) holds then, for any $p \in [1; +\infty)$,

$$ \sup_{\epsilon} \|\varphi_{\epsilon}\|_{W^{1,p}(\Omega)} < \infty. $$

We then conclude by the Sobolev embedding $W^{1,p}(\Omega) \hookrightarrow \mathcal{C}^{0,\beta}(\Omega)$. Thus (2.41) holds.

Finally, to derive a contradiction, we use once again the regularity results of Lemma 13: from (2.41) and the fact that $\partial E_\epsilon$ converges in $\mathcal{C}^1$ to $\partial \Omega$ we should have

$$ 1 = \int_{\partial E_{\epsilon}} \varphi_{\epsilon}^2 \to 0, \qquad (2.42) $$

an obvious contradiction. This concludes the proof of Lemma 31.

□

We can now prove the uniform stability result of Proposition 29.

Proof of Proposition 29. From the lower estimate (2.38) we have, for any critical shape $E$ and any compactly supported vector field

$$ \mathcal{L}_E''(E)[\Phi, \Phi] \ge \left(1 - \frac{1}{\lambda_{0,\rho_E}}\right) \|\langle \Phi, \nu_E \rangle\|_{L^2(\partial E)}^2. $$

From Lemma 31, there exists $\epsilon_1 > 0$ such that, for any $V_0 \in (1 - \epsilon_1; 1)$, for any critical shape $E$,

$$ \lambda_{0,\rho_E} > 2 $$

whence, fixing this $\epsilon_1$, for any $V_0 \in (1 - \epsilon_1; 1)$, for any critical shape $E$,

$$ \mathcal{L}_E''(E)[\Phi, \Phi] \ge 2 \|\langle \Phi, \nu_E \rangle\|_{L^2(\partial E)}^2. $$

This finishes the proof.

□
**Local quantitative inequalities around critical shapes.** We now come to the two consequences of the previous analysis. The first one is that critical shapes are $L^1$ isolated. The second one is that, whenever $V_0$ is close enough to 1, any critical shape is in fact a local $L^1$ minimiser. All the results of this section are straightforward adaptations of similar results in [39].

We begin with the first of these results.

**Lemma 34.** Let $\varepsilon_1 > 0$ be given by Proposition 29. For any $V_0 \in (1 - \varepsilon_1; 1)$ the critical shapes are isolated in the sense that

$$ \inf_{E, E' \text{ critical, } E' \neq E} |E' \Delta E| = \delta(V_0) > 0. $$

This is proved exactly as [39, Proposition 23].

**Lemma 35.** Let $\varepsilon_1$ be given by Proposition 29. For any $V_0 \in (1 - \varepsilon_1; 1)$, any critical shape $E$ is a local $L^1$ minimiser: there exist $c_E > 0$ and $r_E > 0$ such that, for any $f \in \mathcal{F}(V_0)$ with $\|f - 1_E\|_{L^1(\Omega)} \le r_E$, there holds

$$ \mathcal{E}(f) - \mathcal{E}(1_E) \ge c_E \|f - 1_E\|_{L^1(\Omega)}^2. $$

This is a straightforward adaptation of the proof of [39, Theorem 1].

## 2.4 Conclusion of the proof of Theorem I

*Proof of Theorem I.* Consider $\varepsilon \le \varepsilon_1$ where $\varepsilon_1$ is given by Proposition 29. For any $V_0 \in (1 - \varepsilon; 1)$, for any initialisation $f_0 \in \mathcal{F}(V_0)$, let $\{f_k\}_{k \in \mathbb{N}}$ be the sequence generated by the thresholding algorithm. From Lemma 18 the sequence $\{f_k\}_{k \in \mathbb{N}}$ has either one or infinitely many closure points. From Lemma 17 any closure point $f_\infty$ of $\{f_k\}_{k \in \mathbb{N}}$ is a critical point. From Lemma 34, critical points are isolated. Hence, there exists a unique closure point for the sequence $\{f_k\}_{k \in \mathbb{N}}$. Call it $f_\infty$. From Lemma 35, $f_\infty$ is a local minimiser. As any critical point is an extreme point of $\mathcal{F}(V_0)$ the sequence $\{f_k\}_{k \in \mathbb{N}}$ converges strongly in $L^1(\Omega)$ to $f_\infty$. The proof is concluded. $\square$

# 3 Conclusion, generalisations, obstructions

We have thus established the convergence of the thresholding algorithm in three different cases, but under several restrictive assumptions we do not yet know how to bypass. Let us list below certain research questions we plan on tackling in the future.

## 3.1 Lowering the volume constraint threshold

A first major step that would need to be taken next is the extension of the convergence results established in this paper to all possible volume constraints $V_0$. Two remarks are in order in this case. The first one is that one would need an *a priori* regularity theory. Although this is manageable, in the case of energetic functionals, in the two dimensional case thanks to [13], in higher dimensions, the solutions to the problem might no longer be regular enough to apply the shape derivative formalism put in place here. Second, the algorithm converges at best to a critical point. For low volume constraints, we can not guarantee that such critical points are local minimisers, let alone that they are regular enough (even in dimension 2) to use shape hessians, see [12, Remark 3.20].
## 3.2 Degenerate optimal control problems

The first major underlying assumption is that all of our optimal control problems have extremal points of the admissible set of controls as solutions. Put otherwise, this amounts to requiring that the switch function has no flat zone. Now, in linear control problems for semilinear equations it is often the case that the optimal controls are not extreme points: the switch function has flat zones, and this requires a fine tuning of the thresholding method used to obtain satisfactory numerical simulations see [42]. At this point it is unclear how one could tackle this question.

## 3.3 Non-energetic bilinear optimal control problems

There have been many contributions to the qualitative analysis of bilinear optimal control problems in recent years see [38] and the references therein. The problem with this type of queries is that, despite the fact that the thresholding algorithm provides satisfactory simulations, several basic questions would need to be answered: the first one is the regularity of optimal controls. A priori we can not hope for something better than a $C^{1,\alpha}$ boundary, which is not high enough to apply second order shape derivatives arguments. It may be possible to push the methods of [11] to cover this setting but this is at the moment unclear.

## References

[1] C. Amrouche, C. Conca, A. Ghosh, and T. Ghosh. Uniform $W^{1,p}$- estimates for an elliptic operator with Robin boundary condition in a $C^1$ domain. *Calculus of Variations and Partial Differential Equations*, 59(2), Mar. 2020.

[2] S. Amstutz. Analysis of a level set method for topology optimization. *Optim. Methods Softw.*, 26(4-5):555-573, 2011.

[3] S. Amstutz and H. Andrä. A new algorithm for topology optimization using a level-set method. *J. Comput. Phys.*, 216(2):573-588, 2006.

[4] G. Barles and C. Georgelin. A simple proof of convergence for an approximation scheme for computing motions by mean curvature. *SIAM Journal on Numerical Analysis*, 32(2):484-500, apr 1995.

[5] E. Bednarczuk, M. Pierre, E. Rouy, and J. Sokolowski. Calculating tangent sets to certain sets in functional spaces. Research Report RR-3190, INRIA, 1997.

[6] J. Bintz and S. Lenhart. Optimal resource allocation for a diffusive population model. *Journal of Biological Systems*, 28(04):945-976, dec 2020.

[7] B. Braida, J. Dalphin, C. Dapogny, P. Frey, and Y. Privat. Shape and topology optimization for maximum probability domains in quantum chemistry. *Numerische Math.*, To appear.

[8] J. Céa, A. Gioan, and J. Michel. Quelques résultats sur l'identification de domaines. *Calcolo*, 10(3-4):207-232, sep 1973.

[9] A. Chambolle, I. Mazari-Fouquer, and Y. Privat. Stability of optimal shapes and convergence of thresholding algorithms in linear and spectral optimal control problems. Supplementary material.

[10] S. Chanillo, D. Grieser, M. Imai, K. Kurata, and I. Ohnishi. Symmetry breaking and other phenomena in the optimization of eigenvalues for composite membranes. *Communications in Mathematical Physics*, 214(2):315-337, Nov. 2000.
[11] S. Chanillo, D. Grieser, and K. Kurata. The free boundary problem in the optimization of composite membranes, 2000.

[12] S. Chanillo and C. Kenig. Weak uniqueness and partial regularity for the composite membrane problem. *Journal of the European Mathematical Society*, pages 705–737, 2008.

[13] S. Chanillo, C. E. Kenig, and T. To. Regularity of the minimizers in the composite membrane problem in $\mathbb{R}^2$. *Journal of Functional Analysis*, 255(9):2299–2320, Nov. 2008.

[14] A. Cianchi and A. Ferone. A strengthened version of the Hardy-Littlewood inequality. *Journal of the London Mathematical Society*, 77(3):581–592, Feb. 2008.

[15] R. Cominetti and J.-P. Penot. Tangent sets of order one and two to the positive cones of some functional spaces. *Applied Mathematics and Optimization*, 36(3):291–312, nov 1997.

[16] M. Dambrine and J. Lamboley. Stability in shape optimization with second variation. *Journal of Differential Equations*, 267(5):3009–3045, Aug. 2019.

[17] M. Dambrine and M. Pierre. About stability of equilibrium shapes. *ESAIM: Mathematical Modelling and Numerical Analysis*, 34(4):811–834, July 2000.

[18] K. Eppler, H. Harbrecht, and R. Schneider. On convergence in elliptic shape optimization. *SIAM Journal on Control and Optimization*, 46(1):61–83, jan 2007.

[19] L. C. Evans. Convergence of an algorithm for mean curvature motion. *Indiana University mathematics journal*, pages 533–557, 1993.

[20] D. Gilbarg and N. Trudinger. *Elliptic Partial Differential Equations of Second Order*. Springer Berlin Heidelberg, 1983.

[21] A. Henrot. *Extremum problems for eigenvalues of elliptic operators*. Frontiers in Mathematics. Birkhäuser Verlag, Basel, 2006.

[22] A. Henrot and M. Pierre. *Shape Variation and Optimization*. European Mathematical Society Publishing House, Feb. 2018.

[23] A. Henrot, M. Pierre, and M. Rihani. Positivity of the shape Hessian and instability of some equilibrium shapes. *Mediterr. J. Math.*, 1(2):195–214, 2004.

[24] M. Hintermüller, C.-Y. Kao, and A. Laurain. Principal eigenvalue minimization for an elliptic problem with indefinite weight and Robin boundary conditions. *Applied Mathematics & Optimization*, 65(1):111–146, dec 2011.

[25] H. Ishii. A generalization of the Bence, Merriman and Osher algorithm for motion by mean curvature. *Curvature flows and related topics*, pages 111–127, 1995.

[26] H. Ishii, G. E. Pires, and P. E. Souganidis. Threshold dynamics type approximation schemes for propagating fronts. *Journal of the Mathematical Society of Japan*, 51(2), apr 1999.

[27] C.-Y. Kao, Y. Lou, and E. Yanagida. Principal eigenvalue for an elliptic problem with indefinite weight on cylindrical domains. *Math. Biosci. Eng.*, 5(2):315–335, 2008.

[28] C.-Y. Kao and S. A. Mohammadi. Extremal rearrangement problems involving Poisson's equation with Robin boundary conditions. *Journal of Scientific Computing*, 86(3), feb 2021.
[29] C.-Y. Kao, S. A. Mohammadi, and B. Osting. Linear convergence of a rearrangement method for the one-dimensional poisson equation. *Journal of Scientific Computing*, 86(1), jan 2021.

[30] J. Lamboley, A. Laurain, G. Nadin, and Y. Privat. Properties of optimizers of the principal eigenvalue with indefinite weight and Robin conditions. *Calculus of Variations and Partial Differential Equations*, 55(6), Dec. 2016.

[31] T. Laux and F. Otto. Convergence of the thresholding scheme for multi-phase mean-curvature flow. *Calculus of Variations and Partial Differential Equations*, 55(5):1–74, 2016.

[32] T. Laux and D. Swartz. Convergence of thresholding schemes incorporating bulk effects. *Interfaces and Free Boundaries*, 19(2):273–304, 2017.

[33] G. Leoni. *A First Course in Sobolev Spaces*. Graduate studies in mathematics. American Mathematical Society, 2nd. edition, 2017.

[34] E. Lieb and M. Loss. *Analysis*. American Mathematical Society, Providence, Rhode Island, 2001.

[35] I. Mazari. Quantitative inequality for the eigenvalue of a Schrödinger operator in the ball. *Journal of Differential Equations*, 269(11):10181–10238, Nov. 2020.

[36] I. Mazari. *Shape optimization and spatial heterogeneity in reaction-diffusion equations*. Theses, Sorbonne Université, July 2020.

[37] I. Mazari. Quantitative estimates for parabolic optimal control problems under $L^\infty$ and $L^1$ constraints in the ball: Quantifying parabolic isoperimetric inequalities. *Nonlinear Analysis*, 215:112649, Feb. 2022.

[38] I. Mazari, G. Nadin, and Y. Privat. Optimisation of the total population size for logistic diffusive equations: bang-bang property and fragmentation rate. *Communications in Partial Differential Equations*, 47(4):797–828, dec 2021.

[39] I. Mazari and D. Ruiz-Balet. Quantitative stability for eigenvalues of Schrödinger operator, quantitative bathtub principle, and application to the turnpike property for a bilinear optimal control problem. *SIAM Journal on Mathematical Analysis*, 54(3):3848–3883, jun 2022.

[40] B. Merriman, J. Bence, and S. Osher. Diffusion generated motion by mean curvature. In L. A. Department of Mathematics, University of California, editor, *CAM Report 92-33*, 1992.

[41] F. Mignot, J. Puel, and F. Murat. Variation d'un point de retournement par rapport au domaine. *Communications in Partial Differential Equations*, 4(11):1263–1297, 1979.

[42] G. Nadin and A. I. T. Marrero. On the maximization problem for solutions of reaction-diffusion equations with respect to their initial data. *Mathematical Modelling of Natural Phenomena*, 15:71, 2020.

[43] O. Pironneau. *Optimal Shape Design for Elliptic Systems*. Springer Berlin Heidelberg, 1984.

[44] S. J. Ruuth, B. Merriman, J. Xin, and S. Osher. Diffusion-generated motion by mean curvature for filaments. *Journal of Nonlinear Science*, 11(6):473–493, 2001.

[45] D. Swartz and N. K. Yip. Convergence of diffusion generated motion to motion by mean curvature. *Communications in Partial Differential Equations*, 42(10):1598–1643, sep 2017.

[46] M. Yousefnezhad, C.-Y. Kao, and S. A. Mohammadi. Optimal chemotherapy for brain tumor growth in a reaction-diffusion model. *SIAM Journal on Applied Mathematics*, 81(3):1077–1097, jan 2021.