Benoît Braida

Jérémy Dalphin

Charles Dapogny

Pascal Frey

Yannick Privat

SHAPE AND TOPOLOGY OPTIMIZATION FOR MAXIMUM PROBABILITY DOMAINS IN QUANTUM CHEMISTRY

Keywords: Quantum Chemistry, Maximum Probability Domains, Shape Optimization, Level Set Method, Fixed Point Algorithm AMS classification: 35P20, 93B07, 58J51, 49K20. Contents

Shape and topology optimization for maximum probability domains in quantum chemistry

1. Introduction 2. Presentation of the chemical setting and of the considered shape optimization problem 2.1. The wave function of an n-electron system and the probability functionals P ν (Ω) 2.2. The Hartree-Fock approximation for the ground state of an n-electron system 2.3. The shape optimization problem and its relaxation 3. Existence of Maximum Probability Domains The traditional vision of the electron cloud of a molecule used in chemistry rests on the intuition whereby electrons are "localized" in well-defined regions of space around the nuclei, and arranged, for instance, as pairs or shells. In particular, this assumption is the basic principle of the celebrated Lewis model [START_REF] Lewis | The atom and the molecule[END_REF] and of additional concepts, such as resonance and arrow pushing formalism [START_REF] Pauling | The nature of the chemical bond. application of results obtained from the quantum mechanics and from a theory of paramagnetic susceptibility to the structure of molecules[END_REF][START_REF] Shaik | The charge-shift bonding concept: electron-pair bonds with very large ionic-covalent resonance energies[END_REF]. Although this "localized" theory of the organization of electrons shows remarkable agreement with experimental observations, in particular in the field of organic chemistry, it has found severe limitations: in particular, its inability to account for the two-slit diffraction phenomenon of electrons rather vouches for a wave-like description.

Along this line of thinking, the modern theory of quantum chemistry [START_REF] Basdevant | Mécanique quantique[END_REF][START_REF] Goudsmit | La découverte du spin de l'électron[END_REF][START_REF] Rivail | Eléments de chimie quantique à l'usage des chimistes[END_REF][START_REF] Szabo | Modern Quantum Chemistry: introduction to advanced electronic structure theory[END_REF] is based on a "delocalized" point of view about the electrons of a molecule. These are described in terms of a wave function, often denoted by Ψ, encoding the density of probability that the total system be in any possible state (notably in terms of the spatial positions of its elements). One drawback of this vision is that it does not naturally bring into play the historic, key concept of chemical bonds between atoms.

To reconcile this dilemma, much effort has been devoted lately to extract "localized" chemical properties of electrons from the knowledge of the wave function of the system, thus bridging the gap between chemical intuition and accurate predictions from quantum mechanics. Among these so-called "interpretative methods", let us mention the Valence Bond theory [START_REF] Pauling | The nature of the chemical bond. application of results obtained from the quantum mechanics and from a theory of paramagnetic susceptibility to the structure of molecules[END_REF][START_REF]The nature of the chemical bond. ii. the one-electron bond and the three-electron bond[END_REF][START_REF]The Nature of the Chemical Bond[END_REF][START_REF] Shaik | Chemist's Guide to Valence Bond Theory, Chemist's Guide to Valence Bond Theory[END_REF], the Electron Localization Function (ELF) method [START_REF] Becke | A simple measure of electron localization in atomic and molecular-systems[END_REF][START_REF] Savin | Elf: The electron localization function[END_REF][START_REF] Silvi | Classification of chemical-bonds based on topological analysis of electron localization functions[END_REF], the Atoms in Molecules (AIM) theory [START_REF] Bader | Atoms in molecules[END_REF][START_REF]A bond path: A universal indicator of bonded interactions[END_REF], and the Non Covalent Interaction (NCI) framework [START_REF] Contreras-Garcia | Nciplot: A program for plotting noncovalent interaction regions[END_REF][START_REF] Lane | Are bond critical points really critical for hydrogen bonding?[END_REF]; see also [START_REF] Grabowsky | Complementary Bonding Analysis[END_REF]. These methods are still the subject of intense research, stimulated by the tremendous rise in computational power and the development of efficient numerical methods, which allow to consider larger and larger chemical systems. Unfortunately, they all suffer from conceptual obstacles, to the point that their predictions may be difficult to interprete. The featured objects (for instance, the ELF basins, or the bond critical points in AIM) do not always correspond directly to the Lewis electron pairs and chemical bonds which are familiar to experimentalists [START_REF] Ayers | Six questions on topology in theoretical chemistry[END_REF][START_REF]Bond paths are not chemical bonds[END_REF][START_REF] Causa | Atoms and bonds in molecules and chemical explanations[END_REF][START_REF] Poater | Hydrogen-hydrogen bonding in planar biphenyl, predicted by atoms-inmolecules theory, does not exist[END_REF][START_REF]A model of the chemical bond must be rooted in quantum mechanics, provide insight, and possess predictive power[END_REF][START_REF] Savin | The electron localization function (elf ) and its relatives: interpretations and difficulties[END_REF][START_REF]On the significance of elf basins[END_REF].

An alternative interpretative method was proposed about twenty years ago by A. Savin in the seminal article [START_REF] Savin | Probability distributions and valence shells in atoms[END_REF]. This paradigm of Maximal Probability Domains (MPDs) consists in reconstructing localized information about an n-electron system from the wave function Ψ by looking at the probability P ν (Ω) that ν among these n electrons belong to a domain Ω of the Euclidean space R 3 , and that the remaining n -ν electrons lie in the complement R 3 \Ω. Of particular interest are those domains Ω for which the probability P ν (Ω) is maximal, a feature which is naturally expressed in the language of shape optimization. Formally, Maximal Probability Domains (MPDs) are defined as the solutions (whenever they exist) to the shape optimization problem:

(1. [START_REF] Acke | Maximum probability domains for hubbard models[END_REF] sup

Ω⊆R 3 P ν (Ω) .
This formalism is reminiscent of the seminal works [START_REF] Daudel | Sur la localisabilité des corpuscules dans les noyaux et les cortèges électroniques des atomes et des molécules[END_REF][START_REF] Aslangul | Aspects of the localizability of electrons in atoms and molecules: Loge theory and related methods[END_REF] devoted to the loge theory, which brings into play similar probability functionals with the aim to optimize a partition of the molecular space into multiple subsets called "loges" (and not one single domain as in (1.1)). Multiple variants of the problem (1.1) have been considered in the framework of the MPD method, with equal interest in chemical applications. For instance, in [START_REF] Cances | How electrons guard the space: shape optimization with probability distribution criteria[END_REF][START_REF] Savin | Probability distributions and valence shells in atoms[END_REF], the difference ∆P ν (Ω) := P ν (Ω) -P ind ν (Ω) is introduced as an indicator of the valence shell of atoms, where P ind ν (Ω) stands for the probability of finding ν electrons in Ω, were the n electrons in the system completely independent.

From the chemical point of view, the MPD method relies on a simple mathematical formulation, and leads to a geometrical and spatial characterization of the electronic molecular structures and their interactions. It indeed provides unambiguous images of the core and valence domains of atoms [START_REF] Cances | How electrons guard the space: shape optimization with probability distribution criteria[END_REF][START_REF] Savin | Probability distributions and valence shells in atoms[END_REF], of the lone and bonding pairs [START_REF] Menendez | A view of covalent and ionic bonding from maximum probability domains[END_REF], and of the regions where the electrons may move inside a simple molecule [START_REF] Causà | The bond analysis techniques (ELF and maximum probability domains): application to a family of models relevant to bio-inorganic chemistry[END_REF][START_REF] Gallegos | Maximal probability domains in linear molecules[END_REF][START_REF] Lopes | Understanding Maximum Probability Domains with Simple Models[END_REF][START_REF] Menéndez | On the stability of some analytically solvable maximum probability domains[END_REF], a liquid [START_REF] Agostini | Maximum probability domains for the analysis of the microscopic structure of liquids[END_REF], or a crystal [START_REF] Causà | Maximum probability domains in crystals: the rock-salt structure[END_REF][START_REF]Maximum probability domains in the solid-state structures of the elements: the diamond structure[END_REF]. The instance of (1.1) where ν = 2 is especially appealing in chemistry, since the local or global maximizers of the probability of finding exactly 2 electrons can be directly interpreted as the locations of the Lewis pairs from traditional chemistry, or in any event, as stable electronic configurations. We refer to Section 2.1 below for a brief introduction to the method of Maximum Probability Domains, and to [START_REF] Lopes | Understanding Maximum Probability Domains with Simple Models[END_REF] for a more exhaustive presentation.

Maximum Probability Domains have been the focus of relatively few mathematical and algorithmic investigations in the literature, to the best of our knowledge. Shape optimization problems of the form (1.1) present several peculiarities, despite their simple appearance. In particular, the unknown sets Ω may be unbounded and the functional P ν (Ω), although purely geometric in the sense that it does not involve the solution to a partial differential equation depending on Ω, features multiple integrals on Ω and its complement R 3 \Ω, so that it is difficult to evaluate numerically.

In the case of single determinant wave functions, an efficient recursive formula for computing the probability P ν (Ω) was derived in [START_REF] Cances | How electrons guard the space: shape optimization with probability distribution criteria[END_REF], leading to a preliminary implementation of a shape optimization algorithm for the calculation of MPDs in 2d. In the three-dimensional context, a simple algorithm was proposed in [START_REF] Scemama | Investigating the volume maximizing the probability of finding v electrons from variational monte carlo data[END_REF][START_REF] Scemama | Maximum probability domains from quantum monte carlo calculations[END_REF], relying on a Monte Carlo sampling of the wave function and preliminary chemical explorations. Yet another numerical program was proposed in [START_REF] Lopes | Understanding Maximum Probability Domains with Simple Models[END_REF], featuring a discretization of the 3d space with uniform cubic grids and a standard use of shape gradients. The model investigations conducted with these numerical resources have revealed the great potential of the MPD method to unambiguously validate some key concepts in chemistry and suggest new ones [START_REF] Acke | Maximum probability domains for hubbard models[END_REF][START_REF] Agostini | Maximum probability domains for the analysis of the microscopic structure of liquids[END_REF][START_REF] Causa | Electron localization function and maximum probability domains analysis of semi-ionic oxides crystals, surfaces and surface defects[END_REF][START_REF] Hiberty | Pleading for a dual molecular-orbital/valence-bond culture[END_REF][START_REF] Menendez | A view of covalent and ionic bonding from maximum probability domains[END_REF][START_REF] Turek | Bonding in heavier group 14 zero-valent complexes-a combined maximum probability domain and valence bond theory approach[END_REF]. Unfortunately, in spite of their achievements, these algorithms suffer from severe technical limitations, which are prejudicial to the chemical interpretation: notably, the use of uniform Cartesian computational grids inevitably causes large approximations errors, and boundary variation algorithms employed for the numerical resolution of shape optimization problems are well-known to be very sensitive to the initial design. These serious obstructions call for the device of efficient numerical techniques for computing MPDs in the context of large electron systems.

The present article is devoted to the mathematical and numerical analyses of (1.1); our contributions are twofold. From the theoretical point of view, we investigate the existence of Maximum Probability Domains as global solutions to the shape optimization problem (1.1). This question is certainly of theoretical interest; nevertheless, the proof also reveals a new characterization of maximizing domains, suggesting an interesting improvement of the boundary variation algorithms applied hitherto in this context. In Theorem 2.5, roughly speaking, we show the existence of a solution to Problem (1.1), under slight and non-restrictive assumptions on the quantum chemistry model. Furthermore, we prove that any optimal domain Ω * necessarily satisfies the identity

(1.2) Ω * = {F Ω * ν > 0}
where the explicit expression of the function F Ω * ν : R 3 → R is given in Definition 2.4 below. Our second topic of attention in this paper concerns the design of an efficient numerical algorithm for the calculation of MPDs, via the resolution of the shape optimization problem (1.1). "Classical" techniques from shape and topology optimization must be adapted in order to handle several features which are quite specific to the calculation of MPDs. In particular, as we have mentioned, MPDs are often unbounded domains, which is not a common situation in the practice of shape optimization. Also, the strong dependency of the optimized shape on the initialization choice makes it difficult to catch regions far from nuclei, which typically feature very small values of the shape gradient, but which can nevertheless significantly contribute to the value of the optimized probability functional. Moreover, the regions of the 3d space related to chemical interactions between different nuclei often show sharp variations of the shape gradient, which may cause instabilities in the optimization process. To overcome these difficulties, we introduce a hybrid algorithm, combining a very accurate shape gradient algorithm, based on a level set based mesh evolution strategy, with a fixed point method inspired from the identity (1. 2). An open-source implementation of this method, as well as a short manual, can be downloaded at the following address: https://github.com/ISCDtoolbox/MPD Admittedly, multiple contributions in the literature have been devoted to the calculation of MPDs, see for instance [START_REF] Cances | How electrons guard the space: shape optimization with probability distribution criteria[END_REF][START_REF] Scemama | Investigating the volume maximizing the probability of finding v electrons from variational monte carlo data[END_REF][START_REF] Scemama | Maximum probability domains from quantum monte carlo calculations[END_REF][START_REF] Lopes | Understanding Maximum Probability Domains with Simple Models[END_REF][START_REF] Acke | Maximum probability domains for hubbard models[END_REF][START_REF] Agostini | Maximum probability domains for the analysis of the microscopic structure of liquids[END_REF][START_REF] Causa | Electron localization function and maximum probability domains analysis of semi-ionic oxides crystals, surfaces and surface defects[END_REF][START_REF] Hiberty | Pleading for a dual molecular-orbital/valence-bond culture[END_REF][START_REF] Menendez | A view of covalent and ionic bonding from maximum probability domains[END_REF][START_REF] Turek | Bonding in heavier group 14 zero-valent complexes-a combined maximum probability domain and valence bond theory approach[END_REF]). Nevertheless, one unique asset of our method is its ability to capture unbounded domains. To the best of our knowledge, this is the first method for computing MPDs with this property.

The remainder of this article is organized as follows. Section 2 collects background material about quantum mechanics and the Hartree-Fock model for the wave function of the considered electron system. The considered shape optimization problem (1.1) for the calculation of MPDs is then introduced, and the main theoretical result of the article is stated. The next Section 3 is devoted to the proof of this result. The strategy of the proof is constructive, and of interest for numerical purposes. It suggests an efficient computation algorithm, which is the core of Section 4. In the latter, we present the main ingredients of our numerical strategy for the resolution of (1.1): the evaluations of wave functions, orbitals, and shape derivatives are discussed, as well as our combined level set based mesh evolution approach with a new fixed-point method. We finally propose several numerical illustrations of our method in Section 5, ranging from validation examples in the context of simple and well-understood molecules to more challenging computations raising interesting chemical interpretations.

Presentation of the chemical setting and of the considered shape optimization problem

In this section, we outline the physical context of our investigations. A few classical notions about the quantum representation of electron clouds are recalled in Section 2.1; meanwhile, we set the main notation used throughout the article. The next Section 2.2 focuses on the Hartree-Fock approximation for the wave function of a many-electron system, allowing for an explicit expression of the latter. Finally, in Section 2.3, we discuss the defining shape optimization problem (1.1) for MPDs and we state the related theoretical findings of this article.

2.1. The wave function of an n-electron system and the probability functionals P ν (Ω)

The basic material of this section is mainly excerpted from the reference books [START_REF] Basdevant | Mécanique quantique[END_REF][START_REF] Rivail | Eléments de chimie quantique à l'usage des chimistes[END_REF][START_REF] Szabo | Modern Quantum Chemistry: introduction to advanced electronic structure theory[END_REF]. Throughout this article, we consider a quantum system composed of the n ≥ 2 electrons of a molecule; the coordinates of the i th electron are described in terms of its spatial location x i ∈ R 3 and its spin σ i ∈ - 1 2 , 1 2 . The chemical state of the n-electron system is described by a wave function (2.1)

Ψ : R 3 × -1 2 , 1 2 n -→ C x1 σ1 , . . . , xn σn -→ Ψ x1 σ1 , . . . , xn σn , which is normalized in the space L 2 (R 3 × -1 2 , 1 2 ) n , C , that is: (2.2) (σ1,...,σn)∈{-1 2 , 1 2 } n R 3n Ψ x 1 σ 1 , . . . , x n σ n 2 dx 1 ...dx n = 1.
This function Ψ encodes all the available information about the considered electronic system, in the following sense: the joint probability that for all i = 1, ..., n, the i th electron occupy the region Ω i ⊂ R 3 with given spin value σ i = ± 1 2 equals:

Ω1×...×Ωn Ψ x 1 σ 1 , . . . , x n σ n 2 dx 1 ...dx n .
In addition to (2.2), Ψ satisfies the antisymmetry principle (a strong form of the Pauli exclusion principle): for all x1 σ1 , . . . , xn σn ∈ R 3 × - The mathematical assumption encodes the facts that the electrons of the system cannot be distinguished from one another by an external observer -(2.3) indeed implies that the probability density function |Ψ| 2 is unaltered when the coordinates xi σi , xj σj of any two electrons are interchanged -, and that no two electrons can occupy the same spatial position with identical spins.

In this setting, for all ν = 0, . . . , n, the probability P ν (Ω) that exactly ν electrons of the considered system belong to a subset Ω ⊂ R 3 (and that n -ν electrons belong to the complement R 3 \ Ω), regardless of their spins and of their numbering, reads:

(2.4) P ν (Ω) = Ω ν ×(R 3 \Ω) n-ν S ν (x 1 , . . . , x n )dx 1 ...dx n ,
where the function

S ν (x 1 , ..., x n ) ∈ L 1 (R 3n
) is defined by:

(2.5) S ν (x 1 , . . . , x n ) = n ν σ∈{-1 2 , 1 2 } n Ψ x1 σ1 , . . . , xn σn 2 .
Note that in (2.4), we have used the convention A 0 × B = B × A 0 = B for arbitrary sets A and B. Besides, since the total number n of electrons in the system is fixed once and for all, we omit the dependence of P ν and S ν with respect to n.

2.2.

The Hartree-Fock approximation for the ground state of an n-electron system

In this section, we provide a few more details about the wave function Ψ of the considered n-electron system, whose structure is inherited from the so-called Hartree-Fock approximation; see for instance [START_REF] Szabo | Modern Quantum Chemistry: introduction to advanced electronic structure theory[END_REF]Chapter 2 & 3]. This construction relies on the datum of n approximate wave functions φ 1 x σ , . . . , φ n x σ , called spin orbitals, accounting for the individual, independent behavior of each electron in the system. Mathematically, these functions belong to the space L 2 (R 3 × {- 1 2 , 1 2 }, C), and form an orthonormal family of this space:

(2.6) ∀i, j = 1, . . . , n,

σ∈{-1 2 , 1 2 } R 3 φ i x σ φ j x σ dx = 1 if i = j, 0 otherwise.
The wave function Ψ of the total n-electron system is defined from φ 1 , . . . , φ n as a single Slater determinant, that is:

(2.7) Ψ x 1 σ 1 , . . . , x n σ n := 1 √ n! det    φ 1 x1 σ1 • • • φ 1 xn σn . . . . . . φ n x1 σ1 • • • φ n xn σn    ,
a structure under which the antisymmetry principle (2.3) is obviously fulfilled.

As far as the nature of the functions φ 1 , . . . , φ n is concerned, we only consider pure spin orbitals, that is, each function φ i x σ is of the form:

φ i x 1/2 = ψ i (x) φ i x -1/2 = 0 or φ i x 1/2 = 0 φ i x -1/2 = ψ i (x)
, where the spatial orbitals ψ i : R 3 → C are functions of the spatial position x ∈ R 3 only. When the total number n of electrons is even, one may impose that the n spin orbitals φ 1 , . . . , φ n be grouped by pairs, being associated to the same spatial orbital, with opposite spins: this is the so-called restricted (or closed shell ) Hartree-Fock framework. On the contrary, in the unrestricted (or open shell ) Hartree-Fock framework, the spatial orbitals ψ i of the different spin orbitals φ i are free to be independent functions.

Usually, in computational chemistry, each spatial orbital ψ i is expanded as a linear combination of Gaussian basis functions, namely;

(2.8) ∀x = (x, y, z) ∈ R 3 , ψ i (x) = R r=1 c r (x -x r ) αr (y -y r ) βr (z -z r ) γr e -ar|x-xr| 2 ,
for some real numbers c r ∈ R and a r > 0, r = 1, . . . , R; see again [START_REF] Szabo | Modern Quantum Chemistry: introduction to advanced electronic structure theory[END_REF]. From the physical point of view, x r := (x r , y r , z r ) ∈ R 3 is the center of the r th nucleus of the considered molecule, and α r , β r , γ r are integer exponents, accounting for the type of the spatial orbital ψ i . The values of the coefficients c r featured in (2.8) are characterized as those minimizing the energy of the wave function Ψ of the total system. In practice, the description (2.8) of the optimized spatial orbitals ψ i (and thereby, that of the Hartree-Fock wave function Ψ of the total system) is supplied by the computational chemistry package Gaussian [START_REF] Frisch | Gaussian~16 Revision B.01[END_REF].

Let us mention that the above Hartree-Fock framework is only one particular, single-determinant approximation of the wave function Ψ of the system, which arises in general as the sum of several Slater determinants of the form (2.7), constructed from different sets of spin orbitals. Other single-determinant approximations of Ψ (all originating from single-particle models) may be considered, such as those supplied by the Kohn-Sham Density Functional Theory; see for instance [START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF].

In any event, the theoretical developments of this article do not require such a precise knowledge as (2.8) about the nature of the wavefunction Ψ. Often, we shall solely assume that the following Assumption 2.1 is fulfilled.

Assumption 2.1.

• The function Ψ is antisymmetric with respect to the interchange between the spatial and spin coordinates of any two electrons i and j, i.e. (2.3) holds.

• The function Ψ is normalized: for σ = (σ 1 , ..., σ n ) ∈ {-1 2 , 1 2 } n , let us introduce the reduced functional (2.9) Ψ σ : R 3 × . . . × R 3 -→ C (x 1 , . . . , x n ) -→ Ψ x1 σ1 , . . . , xn σn ;
we then impose that:

(2.10)

σ∈{-1 2 , 1 2 } n (R 3 ) n |Ψ σ (x 1 , . . . , x n )| 2 dx 1 . . . dx n = 1.
• The function Ψ σ is real analytic on R 3n .

Remark 2.2.

• Our main theoretical result, Theorem 2.5 below, actually holds under the slightly less restrictive assumption that Ψ be real analytic on R 3n \Z, where Z is a finite set of points. This observation is of particular interest when different structures from (2.8) are assumed about the spatial orbitals ψ i ; for instance, another popular set of basis functions in computational chemistry is that of Slater orbitals, of the form Q(x -x r )e -|x-xr| where Q(x) is a trivariate polynomial; these functions are not differentiable at x r , thus reflecting the singularity of the wave function Ψ near the center of the nuclei of the considered molecule. • Our theoretical framework below relies only on Assumption 2.1, leaving the room for multi-determinant wave functions Ψ. However, the numerical algorithms detailed in Section 4, used to calculate P ν and its derivative, are not easily generalized to this multi-determinant case. For this reason, the applications of Section 5 stick to single-determinant wave functions as presented in this section. Note however that these could actually be supplied by different chemical approximations from the Hartree-Fock framework, as we have already hinted at.

The shape optimization problem and its relaxation

Let 0 ≤ ν ≤ n be given; the shape optimization problem characterizing the Maximum Probability Domains of interest reads:

(P s.o. ν ) sup Ω∈M P ν (Ω),
where the probability P ν (Ω) that exactly ν electrons belong to Ω is defined by (2.4), and M is the set of all Lebesgue measurable sets in R 3 .

Remark 2.3. The value of the functional P ν (Ω) is obviously unaltered if the considered shape Ω is perturbed by a set with null Lebesgue measure. Hence, (P s.o. ν ) could be formulated in terms of the corresponding equivalence classes of measurable subsets of R 3 . For simplicity, we prefer to stick to the formulation (P s.o. ν ) of the problem, omitting at times that its solutions ought to be understood up to negligible sets.

Our theoretical analysis of (P s.o.

ν ) relies on a fairly classical strategy involving a relaxed optimization problem; the latter features an enlarged set of admissible designs and a corresponding extended definition of the optimized functional P ν , with better mathematical properties; see e.g. [START_REF] Henrot | Shape variation and optimization: a geometrical analysis[END_REF]. We trade "black-andwhite" designs Ω ∈ M (or, equivalently their characteristic functions 1 Ω ∈ L ∞ (R 3 , {0, 1})) for more general "density" functions m ∈ L ∞ (R 3 , [0, 1]) which are allowed to take intermediate, "grayscale" values in (0, 1). This is motivated by the well-known fact that [START_REF] Henrot | Shape variation and optimization: a geometrical analysis[END_REF] Prop. 7.2.17). The extended version of the objective

L ∞ (R 3 , [0, 1]) is the weak- * closure of L ∞ (R 3 , {0, 1}) in L ∞ (R 3 ) (see
P ν (Ω) to density functions m ∈ L ∞ (R 3 ) is the functional J ν : L ∞ (R 3 ) → R defined by: (2.11) J ν (m) = R 3n P ν (m)(x 1 , . . . , x n )S ν (x 1 , . . . , x n ) dx 1 . . . dx n ,
where P ν : L ∞ R 3 -→ L ∞ R 3n is given by:

(2.12)

P ν (m) : (x 1 , . . . , x n ) → ν i=1 m (x i ) n j=ν+1
(1 -m (x j )) .

Obviously, it holds P ν (Ω) = J ν (1 Ω ) for any Ω ∈ M. The relaxed version of our problem (P s.o. ν ) now reads:

(P relax ν ) sup m∈L ∞ (R 3 ,[0,1]) J ν (m).
The main theoretical result of this article ensures that, under mild assumptions, global solutions to (P s.o. ν ) and (P relax ν ) coincide. Before giving a precise statement and a complete proof of this fact in the next Section 3, a definition is in order: Definition 2.4. For an arbitrary density m ∈ L ∞ R 3 and for ν = 0, ..., n, the function F m ν ∈ L 1 (R 3 ) is defined by:

(2.13) F m ν (x) = R 3(n-1) (ν -nm(x ν+1 )) ν i=2 m(x i ) n j=ν+2 (1 -m(x j ))S ν (x, x 2 , . . . , x n )dx 2 . . . dx n ,
with the conventions that products over reversed indices equal 1 and that m(x n+1 ) = 0. When m = 1 Ω is the characteristic function of a shape Ω ∈ M, we use the shorthand

F Ω ν ≡ F 1Ω ν .
Theorem 2.5. Let the wave function Ψ fulfill Assumption 2.1, and assume in addition that there exists

m 0 ∈ L ∞ (R 3 , [0, 1]
) such that:

(2.14) J ν (m 0 ) 1 + n -ν ν + ν n -ν -1
, with the convention that the above right-hand side equals 1 if ν = 0 or ν = n. Then, ) is an open subset of R 3 which satisfies the identity:

(2.15)

Ω * = {F Ω * ν > 0} := {x ∈ R 3 , F Ω * ν (x) > 0}, up to a subset of R 3 with null Lebesgue measure.
Remark 2.6. The characterization (2.15) of an optimal domain Ω * for (P s.o.

ν ) actually implies that it is much more regular than a mere open subset of R 3 . It indeed entails that Ω * is semi-analytic, i.e. for every point x 0 ∈ Ω * , there exists a neighborhood U of x 0 in R 3 as well as 2pq real analytic functions g ij , h ij , i = 1, . . . , p, j = 1, . . . , q, such that:

Ω * ∩ U = p i=1 {x ∈ U | g ij (x) = 0 and h ij (x) > 0, j = 1, . . . , q}.
Semi-analytic sets enjoy quite convenient properties, such as being stratifiable in the sense of Whitney; see [START_REF] Hardt | Stratification of real analytic mappings and images[END_REF][START_REF] Hironaka | Subanalytic sets, in Number theory, algebraic geometry and commutative algebra[END_REF]. Since they are not needed in the following, we do not insist further on these features.

In spite of its theoretical aspect, Theorem 2.5 has interesting practical consequences. As we shall see in Section 4.5 below, the characterization (2.15) of optimal domains paves the way to an efficient fixed-point iteration procedure for the calculation of Ω * .

Let us end this section with a few comments about the assumption (2.14) in Theorem 2.5. Since the latter is rather technical and not very intuitive, we point out three simple and already quite general sufficient conditions under which it is satisfied.

(1) A calculation reveals that (2.14) is implied by the following assumption:

There exists Ω 0 ∈ M such that P ν (Ω 0 ) min 1 3 , ν n , n -ν n .
This is a simple consequence of the inequality:

∀n ∈ IN \ {0} , ∀ν = 0, . . . , n, 1 + n -ν ν + ν n -ν -1 ≤ min 1 3 , ν n , n -ν n .
(2) Testing whether (2.14) is satisfied by a constant density m(x) = m for a.e. x ∈ R 3 , and recalling that R 3n S ν (x 1 , . . . , x n ) dx 1 . . . dx n = n ν , it follows that (2.14) holds true as soon as:

max m∈[0,1] n ν m ν (1 -m) n-ν 1 + n -ν ν + ν n -ν -1
, or, after calculating explicitly the maximum in the above left-hand side:

(2.16) n ν ν n ν 1 - ν n n-ν 1 + n -ν ν + ν n -ν -1 .
This condition is easy to check numerically; for instance, we identify in Table 1 the values of ν for which it is satisfied, for small electron numbers n.

(3) Regardless of the total number n of electrons, (2.16) (and so (2.14)) is always satisfied if ν = 0, 1 or if ν = 2, a situation of much interest in chemistry, as we have already hinted at. This point is indeed obvious for ν = 0 and 1; in the case ν = 2, (2.16) is actually equivalent to:

f (n) < 0 where f (s) := s s -s(s -1)(s -2) s-3 (s 2 -2s + 4);
The latter fact holds true, as follows from an elementary study of the function f : R → R.

Remark 2.7. In addition to the previous sufficient conditions, (2.16) is asymptotically satisfied in the following sense: for a fixed integer ν, let n tend to ∞. Using the classical equivalents n ν ∼ n ν /ν! and (n -ν) n-ν-1 ∼ e -ν n n-ν-1 as n → ∞, we see that, in this limit, (2.16) is equivalent to the (true) inequality:

ν ν-1 n 3 ≥ e ν ν!.
Remark 2.8. The assumption (2.14), or the weaker version (2.16), are technical sufficient conditions for the conclusions of Theorem 2.5 to hold, which are naturally brought into play in the proof of the latter, and do not lend themselves to a satisfactory physical interpretation.

Indices ν = 0, . . . , n for which (2.16) holds Indices ν = 0, . . . , n for which (2.16) ) for various values of the number n of electrons and of the index ν = 0, . . . , n.

fails n = 1 0, 1 n = 2 0, 1, 2 n = 3 0, 1, 2, 3 n = 4 0, 1, 2, 3, 4 n = 5 0, 1, 2, 3, 4, 5 n = 6 0, 1,

Existence of Maximum Probability Domains

This section is dedicated to the proof of Theorem 2.5 about the existence of MPDs; unless stated otherwise, we assume that n ≥ 2, and that the wave function Ψ satisfies Assumption 2.1. The proof follows a quite classical outline when it comes to proving the existence of a solution to a shape optimization problem. At first, we consider the relaxed version (P relax ν ) of the shape optimization problem (P s.o. ν ) of interest. Preliminary properties of this problem are investigated in Section 3.1; then, we prove in Section 3.2 that (P

relax ν ) has a global maximizer m * ∈ L ∞ (R 3 , [0, 1]
). The next Section 3.3 is devoted to the calculation of the firstorder optimality conditions associated to the relaxed problem (P relax ν ). Eventually, in Section 3.4, we prove that any global maximizer m * of (P relax ν ) is actually the characteristic function 1 Ω * of a measurable subset Ω * ∈ M which is optimal for (P s.o. ν ).

Preliminary remarks

We start with simple and useful properties about the function S ν (x 1 , . . . , x n ) in (2.5) and the relaxed objective function J ν (m) defined by (2.11).

Lemma 3.1. For n ≥ 2 and 0 ≤ ν ≤ n, the following properties hold:

(i) The function S ν is symmetric with respect to the permutation of any two variables:

∀i < j = 1, ..., n, ∀(x 1 , . . . x n ) ∈ R 3n , S ν (x 1 , . . . , x i , . . . , x j , . . . , x n ) = S ν (x 1 , . . . , x j , . . . , x i , . . . , x n ). (ii) It holds J n (1) = J 0 (0) = 1. More generally, for any 0 ≤ ν ≤ n and any density m ∈ L ∞ (R 3 ), J ν (m) = J n-ν (1 -m). Consequently, if m * is one maximizer of J ν over L ∞ (R 3 , [0, 1]), then 1 -m * is also one maximizer of J n-ν over L ∞ (R 3 , [0, 1]). (iii) For arbitrary m ∈ L ∞ (R 3 ), one has n ν=0 J ν (m) = 1.
Proof. (i): This is an immediate consequence of the antisymmetry principle (2.3).

(ii): This follows readily from the definitions (2.11) and (2.12) of J ν and P ν .

(iii): Using the definitions of P ν and S ν , we calculate:

(3.1) n ν=0 J ν (m) = R 3n n ν=0 n ν ν i=1 m (x i ) n j=ν+1 (1 -m (x j ))    σ∈{-1 2 , 1 2 } n Ψ x1 σ1 , . . . , xn σn 2    dx 1 . . . dx n = R 3n n ν=0 Iν ⊂{0,...,n} #Iν =ν i∈Iν m (x i ) j / ∈Iν (1 -m (x j ))    σ∈{-1 2 , 1 2 } n Ψ x1 σ1 , . . . , xn σn 2    dx 1 . . . dx n = R 3n I⊂{0,...,n} i∈I m (x i ) j / ∈I (1 -m (x j ))    σ∈{-1 2 , 1 2 } n Ψ x1 σ1 , . . . , xn σn 2    dx 1 . . . dx n ,
where we have used the antisymmetry principle (2.3) to pass from the first line to the second one. We now remark that:

1 = (m(x 1 ) + 1 -m(x 1 )) . . . (m(x n ) + 1 -m(x n )) = I⊂{0,...,n} i∈I m (x i ) j / ∈I (1 -m (x j )) ,
where we used a direct expansion of the product in the first line of the above formula. Combining this with (3.1) and the normalization condition (2.10), the desired result follows.

We now state two easy technical identities about the function F m ν introduced in Definition 2.4. Lemma 3.2. Let m be an arbitrary element in L ∞ (R 3 , [0, 1]); the following identities hold:

R 3 m(x)F m ν (x)dx = νP ν (m) -(n -ν) P ν+1 (m) and R 3 (1 -m(x)) F m ν (x)dx = νP ν-1 (m) -(n -ν) P ν (m) . (3.2)
with the convention that P n+1 (m) = P -1 (m) = 0.

Proof. For brevity, we solely verify these formulas in the case where 0 < ν < n. Inserting the identity

∀ x ν+1 ∈ R 3 , ν -nm (x ν+1 ) = ν(1 -m * (x ν+1 )) -(n -ν) m (x ν+1 ) ,
into the definition (2.13) of F m ν , we obtain, for arbitrary x 1 ∈ R 3 :

F m ν (x 1 ) = ν (R 3 ) n-1 ν i=2 m (x i ) n j=ν+1 (1 -m(x j ))S ν (x 1 , . . . , x n ) dx 2 . . . dx n -(n -ν) (R 3 ) n-1 ν+1 i=2 m (x i ) n j=ν+2 (1 -m(x j ))S ν (x 1 , . . . , x n ) dx 2 . . . dx n .
Both formulas in (3.2) are then obtained by multiplying the above identity by m(x 1 ) or (1 -m(x 1 )), integrating over R 3 , and using the symmetry of S ν (x 1 , . . . , x n ) with respect to the permutation of any two of the variables x i , x j .

3.2.

Existence of a solution to the relaxed problem (P relax ν )

Let us start with an elementary lemma from measure theory.

Lemma 3.3. Let (X, µ) and (Y, ν) be two complete, σ-finite measured spaces, and let the product space X × Y be equipped with the completed product measure. Let f k (x, y) and ϕ k (y) be two sequences of functions in L 1 (X × Y ) and L ∞ (Y ), respectively; we define the sequence g k ∈ L 1 (X) by:

g k (x) = Y ϕ k (y)f k (x, y) dν(y).
We assume that

f k k→∞ ----→ f ∞ strongly in L 1 (X × Y ), and ϕ k k→∞ ----→ ϕ ∞ weakly - * in L ∞ (Y ).
Then,

g k k→∞ ----→ g ∞ , strongly in L 1 (X), where g ∞ (x) := Y ϕ ∞ (y)f ∞ (x, y) dν(y).
Proof. Let us first observe that the Fubini theorem (see e.g. [START_REF] Rudin | Real and complex analysis[END_REF]) implies that the functions g k (x) and g ∞ (x) are indeed defined µ-almost everywhere on X, and that they belong to L 1 (X). We then decompose:

g k (x) -g ∞ (x) = h (1) k (x) + h (2) k (x), where h (1) k (x) := Y ϕ k (y)(f k (x, y) -f ∞ (x, y))dν(y), and h (2) k (x) = Y (ϕ k (y) -ϕ ∞ (y))f ∞ (x, y)dν(y).

The function h

(1) k (x) can be estimated as:

X |h (1) k (x)|dµ(x) ≤ ϕ k L ∞ (Y ) f k -f ∞ L 1 (X×Y ) , so that h (1) k → 0 strongly in L 1 (X) as k → ∞.
On the other hand, the Fubini theorem ensures the existence of a subset A ⊂ X with µ(A) = 0 such that for any x ∈ X \A, the partial mapping y → f ∞ (x, y) belongs to L 1 (Y ). From the L ∞ (Y )-weak- * convergence of ϕ k , we infer that:

∀x ∈ X \ A, h (2) k (x) k→∞ ----→ 0. Since X |h (2) k (x)|dµ(x) ≤ sup k ϕ k L ∞ (Y ) + ϕ ∞ L ∞ (Y ) f ∞ L 1 (X×Y ) ,
it follows from the Lebesgue dominated convergence theorem that h Proof. Without loss of generality, we limit ourselves with proving that J n satisfies the required property, the adaptation of the argument to the case where 0 ≤ ν < n being straightforward.

Let m k be a sequence of density functions in L ∞ R 3 , [0, 1] which converges to m ∞ for the weak- * convergence. By definition,

J n (m k ) = R 3 m k (x 1 ) R 3 m k (x 2 ) • • • R 3 m k (x n )S ν (x 1 , • • • , x n ) dx n • • • dx 2 dx 1 .
Using Lemma 3.3, the functions

R 3(n-1) (x 1 , • • • , x n-1 ) → R 3 m k (x n )S ν (x 1 , • • • , x n ) dx n belong to L 1 (R 3n-3
) and form a sequence which converges strongly in this space to the function:

(x 1 , • • • , x n-1 ) → R 3 m ∞ (x n )S ν (x 1 , • • • , x n ) dx n .
Repeating this argument with respect to the variables x n-1 , . . . , x 1 successively, the desired result follows.

The existence of a maximizer for (P relax Proof. We apply the direct method of calculus of variations for the existence of solutions to a maximization problem. Let (m k ) k∈IN be a maximizing sequence for (P relax ν ). Since this sequence is uniformly bounded (by 1) in L ∞ (R 3 ), the Banach-Alaoglu theorem ensures that it converges (up to a subsequence still denoted by

m k ) weakly- * in L ∞ (R 3 ) to an element m ∞ ∈ L ∞ (R 3 ). Besides, one has: R 3 m k (x)φ(x)dx k→∞ ----→ R 3 m ∞ (x)φ(x)dx, and R 3 (1 -m k (x))φ(x)dx k→∞ ----→ R 3 (1 -m ∞ )(x)φ(x)dx, for any positive function φ ∈ C ∞ c (R 3
), and so it easily follows that the function m ∞ actually belongs to

L ∞ (R 3 , [0, 1]).
Using now the continuity of J ν supplied by Proposition 3.4, we see that:

sup m∈L ∞ (R 3 ,[0,1]) J ν (m) = lim k→∞ J ν (m k ) = J ν (m ∞ ),
which completes the proof.

3.3.

Optimality conditions for the relaxed problem (P relax ν )

Let us first calculate the derivative of the functional J ν (m).

Proposition 3.6. The functional

J ν (m) is Fréchet differentiable at any density m ∈ L ∞ (R 3 , [0, 1]). Its
Fréchet derivative, which we denote by h → J ν (m)(h), reads:

(3.3) ∀h ∈ L ∞ (R 3 ), J ν (m)(h) = R 3 F m ν (x)h(x)dx,
where

F m ν (x) is defined in (2.13). Proof. At first, the mapping m → P ν (m) is Fréchet differentiable from L ∞ (R 3 ) into L ∞ (R 3n
) and its Fréchet derivative reads:

∀h ∈ L ∞ (R 3 ), P ν (m)(h) =   n j=ν+1 (1 -m(x j ))   ν i=1 m(x 1 ) . . . m(x i-1 )h(x i )m(x i+1 ) . . . m(x ν ) - ν i=1 m(x i ) n j=ν+1 (1 -m(x ν+1 )) . . . (1 -m(x j-1 ))h(x j )(1 -m(x j+1 )) . . . (1 -m(x n )). Since S ν (x 1 , . . . , x n ) belongs to L 1 (R 3n ), J ν (m) is Fréchet differentiable at any m ∈ L ∞ (R 3
), and its Fréchet derivative reads, for arbitrary h ∈ L ∞ (R 3 ):

J ν (m)(h) = R 3n   n j=ν+1 (1 -m(x j ))   ν i=1 m(x 1 ) . . . m(x i-1 )h(x i )m(x i+1 ) . . . m(x ν ) S ν (x 1 , . . . , x n )dx 1 . . . dx n - R 3n ν i=1 m(x i )   n j=ν+1 (1 -m(x ν+1 )) . . . (1 -m(x j-1 ))h(x j )(1 -m(x j+1 )) . . . (1 -m(x n ))   S ν (x 1 , . . . , x n )dx 1 . . . dx n .
Using the symmetry of S ν with respect to the exchange of any two of the x i and changing variables in the above integrals, it follows readily that:

J ν (m)(h) = R 3n νh(x 1 ) ν i=2 m(x i )   n j=ν+1 (1 -m(x j ))   S ν (x 1 , . . . , x n )dx 1 . . . dx n - R 3n (n -ν)h(x 1 ) ν i=2 m(x i ) m(x ν+1 )   n j=ν+2 (1 -m(x j ))   S ν (x 1 , . . . , x n )dx 1 . . . dx n .
Rearranging both integrals and using the Fubini theorem, we finally arrive at:

J ν (m)(h) = R 3 F m ν (x)h(x)dx,
which is the desired result.

We now infer from the previous result the first-order necessary conditions of optimality for the relaxed problem (P relax ν ).

Proposition 3.7. Let m * ∈ L ∞ (R 3 , [0, 1]) be one local maximizer of (P relax ν ). Then, the following inclusions hold, up to subsets of R 3 with null Lebesgue measure:

(i) {F m * ν < 0} ⊆ {m * = 0} ⊆ {F m * ν 0}, (ii) {F m * ν > 0} ⊆ {m * = 1} ⊆ {F m * ν 0}, (iii) {0 < m * < 1} ⊆ {F m * ν = 0}. Proof. Since F m * ν is in L 1 (R 3 ) we know from the Lebesgue differentiation theorem that almost every point x 0 ∈ R 3 is a Lebesgue point for F m *
ν , that is:

(3.4) 1 |B(x 0 , r)| B(x0,r) |F m * ν (x) -F m * ν (x 0 )|dx r→0 ---→ 0.
For the same reason, we know that almost every point x 0 in {m * = 1} has density one:

(3.5) |B(x 0 , r) ∩ {m * = 1} | |B(x 0 , r)| r→0 ---→ 1;
see [39, §1.7] about these facts.

Let us then consider any point x 0 ∈ {m * = 1} which is a Lebesgue point of F m * ν and a point with density one of {m * = 1}. The strategy of the proof is to use both properties to construct an admissible perturbation of m * and to extract information about the integrand F m * ν of the derivative of J ν (m * ) around x 0 . More precisely, introducing a sequence r k of positive numbers decreasing to 0, let us introduce the sets G k := B(x 0 , r k ) ∩ {m * = 1}; it follows from the definitions that the sequence of functions g k := 1 G k is such that m * -ηg k belongs to L ∞ (R 3 , [0, 1]) for any 0 < η < 1. Since m * is a global maximizer of (P relax ν ), it also holds:

J ν (m * -ηg k ) -J ν (m * ) 0; dividing the previous inequality by η and letting η → 0, we infer that:

-J ν (m * )(g k ) = - G k F m * ν (x)dx 0.
Now dividing both sides by |G k |, letting k → +∞ and using that x 0 satisfies (3.4) and (3.5), we obtain that F m * ν (x 0 ) 0. Hence, we have proved that: For a.e.

x ∈ R 3 s.t. m * (x) = 1, F m * ν (x) 0.
A similar argument based on perturbations of the form m * + ηh k of m * , where ν ), one has necessarily:

h k = 1 H k is the charac- teristic function of a subset H k = {m * = 0} ∩ B(x 0 , r k ) centered
• For a.e. x ∈ Ω, F Ω ν (x) ≥ 0; • For a.e. x ∈ R 3 \ Ω, F Ω ν (x) ≤ 0.
3.4. End of the proof of Theorem 2.5

Proof of (i): This is exactly the statement of Corollary 3.5 above.

Proof of (ii): For simplicity, we assume that 0 < ν < n (the cases where ν = 0 and ν = n being handled in a similar way). Let m * ∈ L ∞ (R 3 , [0, 1]) be any global maximizer of the relaxed problem (P relax ν ); we prove that the set {0 < m * < 1} has null Lebesgue measure.

To this end, we assume that {0 < m * < 1} has positive measure, and work towards a contradiction with the assumption (2.14). From Proposition 3.7, the set {F m * ν = 0} has positive Lebesgue measure. Since S ν is a real analytic function over R 3n , standard results about the analyticity under the integral sign (see e.g. [START_REF] Tauvel | Analyse complexe pour la licence 3: cours et exercices corrigés[END_REF], §7.5) imply that F m * ν is real analytic over R 3 . That this function vanish on a subset of R 3 with positive Lebesgue measure then implies that it must vanish identically on R 3 ; hence {F m * ν = 0} = R 3 . We now infer from the identities of Lemma 3.2 that:

(3.6) νJ ν (m * ) = (n -ν)J ν+1 (m * ) and νJ ν-1 (m * ) = (n -ν)J ν (m * ),
and so:

J ν-1 (m * ) + J ν+1 (m * ) = J ν (m * ) n -ν ν + ν n -ν .
Now using Lemma 3.1 yields:

1 J ν-1 (m * ) + J ν (m * ) + J ν+1 (m * ) 1 + n -ν ν + ν n -ν J ν (m * ).
Therefore, one has necessarily

J ν (m * ) 1 + n -ν ν + ν n -ν -1
.

The last inequality is actually strict. Indeed, assume that

J ν (m * ) = 1 + n -ν ν + ν n -ν -1 .
Then, according to the above computations, one has J k (m * ) = 0 for every k ∈ {1, . . . , ν -2, ν + 2, . . . , n} (recall that we have assumed ν = 0). Since P k (m * ) and S ν are nonnegative functions, this entails P k (m * )S ν = 0 in R 3n . Since {0 < m * < 1} has positive measure, we infer that S ν necessarily vanishes on a non-zero Lebesgue measure set of R 3n , which is impossible by analyticity.

All things considered, we have proved that:

J ν (m * ) < 1 + n -ν ν + ν n -ν -1
, which is the sought contradiction with (2.14). This proves that both sets {F m * ν = 0} and {0 < m * < 1} have null Lebesgue measure, and so m * is the characteristic function of a measurable set Ω * ∈ M, as desired.

Proof of (iii): As an immediate consequence of the previous point, any global maximizer m * of the relaxed problem (P relax ν ) (the existence of which being guaranteed by (i)) is the characteristic function 1 Ω * of a measurable set Ω * ∈ M, which is necessarily a global maximizer of (P s.o. ν ).

Proof of (iv): Let Ω * ∈ M be any global solution to (P s.o. ν ); it follows from the previous analyses that m * := 1 Ω * is one global solution to (P s.o.

ν ). The fact proved in (ii) that F m * ν = 0 has zero Lebesgue measure, together with the necessary optimality conditions of Proposition 3.7 imply in this case:

{F m * ν < 0} = {m * = 0}, {F m * ν > 0} = {m * = 1},
these identities being understood up to subsets of R 3 with null Lebesgue measure. This yields (2.15), and the set

Ω * is open since F m * ν is continuous over R 3 .

A shape and topology optimization algorithm for the calculation of Maximum Probability Domains

We now turn to the numerical resolution of our shape optimization problem (P s.o. ν ), emphasizing on the specific treatments imposed by the context of quantum chemistry.

In the first Section 4.1, we recall a few basic notions about shape derivatives. Then, in Sections 4.2 and 4.3, we discuss the practical calculation of the probability functionals P ν (Ω) and their shape derivatives. We then outline in Section 4.4 the level set based mesh evolution strategy used for the numerical representation of shapes and their evolution, and we sketch the shape gradient algorithm inspired by the aforementioned ingredients. In Section 4.5, we present a fixed point algorithm based on the characterization (2.15) of optimal shapes, whose combination with our shape gradient algorithm yields an optimization strategy which conveniently handles topological changes, as well as the capture of regions of space where the shape gradient takes small values. Finally, the construction of the fixed box D where our computations take place in practice is not trivial when it comes to the calculation of MPDs, since the latter may be unbounded domains; we discuss this task in Section 4.6.

Computation and use of shape derivatives

The numerical resolution of a shape optimization problem of the form (P s.o.

ν ) rests upon the knowledge of the derivative of the optimized functional with respect to the domain. This notion can be understood in a variety of manners, one of them relying on Hadamard's boundary variation method which we sketch briefly for the convenience of the reader; see e.g. [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Henrot | Shape variation and optimization: a geometrical analysis[END_REF][START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF] and the recent survey [START_REF] Allaire | Shape and topology optimization[END_REF] for more exhaustive presentations.

In this framework, variations of a Lipschitz domain Ω are considered under the form (4.1)

Ω θ = (Id + θ)(Ω), where θ ∈ W 1,∞ (R 3 , R 3
) is a "small" vector field and Id : R 3 → R 3 is the identity mapping. One function J(Ω) of the domain is said to be shape differentiable at a particular shape Ω if the underlying mapping

θ → J(Ω θ ) from W 1,∞ (R 3 , R 3 ) into R is Fréchet differentiable at θ = 0.
The corresponding derivative J (Ω)(θ) is the shape derivative of J at Ω and the following expansion holds:

(4.2) J(Ω θ ) = J(Ω) + J (Ω)(θ) + o(θ), where o(θ) θ W 1,∞ (R 3 ,R 3 ) θ→0 ---→ 0.
The following proposition supplies the shape derivative of the probability functionals P ν (Ω): Proposition 4.1. For n ≥ 2 and 0 ≤ ν ≤ n, the function of the domain P ν (Ω) is shape differentiable at any (possibly unbounded) Lipschitz domain Ω, and the corresponding shape derivative reads:

(4.3) ∀θ ∈ W 1,∞ (R 3 , R 3 ), P ν (Ω)(θ) = ∂Ω F Ω ν θ • n ds,
where n is the unit normal vector to ∂Ω, pointing outward Ω, and F Ω ν is defined in (2.13).

Proof. For θ ∈ W 1,∞ (R 3 , R 3 ), a change of variables yields (see [START_REF] Evans | Measure theory and fine properties of functions[END_REF] in the present Lipschitz context):

P ν (Ω θ ) = Ω ν θ (R 3 \Ω θ ) n-ν S ν (x 1 , . . . , x n )dx 1 . . . dx n , = Ω ν (R 3 \Ω) n-ν | det(I + ∇θ(x 1 ))| . . . | det(I + ∇θ(x n ))|S ν (x 1 + θ(x 1 ), . . . , x n + θ(x n ))dx 1 . . . dx n .
Classical shape calculus then entails that P ν (Ω) is shape differentiable, with derivative:

P ν (Ω)(θ) = ν Ω ν (R 3 \Ω) n-ν (divθ(x 1 )S ν (x 1 , . . . , x n ) + ∇ x1 S ν (x 1 , . . . , x n ) • θ(x 1 )) dx 1 . . . dx n + (n -ν) Ω ν (R 3 \Ω) n-ν divθ(x ν+1 )S ν (x 1 , . . . , x n ) + ∇ xν+1 S ν (x 1 , . . . , x n ) • θ(x ν+1 ) dx 1 . . . dx n ,
where we have used the symmetry of S ν (x 1 , . . . , x n ) with respect to the exchange of any two of the variables x i . Now integrating by parts, we obtain the surface expression:

P ν (Ω)(θ) = ν ∂Ω Ω ν-1 (R 3 \Ω) n-ν S ν (x 1 , . . . , x n )dx ν+1 . . . dx n dx 2 . . . dx ν θ • n(x 1 )ds(x 1 ) -(n -ν) ∂Ω Ω ν (R 3 \Ω) n-ν-1 S ν (x 1 , . . . , x n )dx ν+2 . . . dx n dx 1 . . . dx ν θ • n(x ν+1 )ds(x ν+1 ).
After inspection, and using once more the symmetry of S ν (x 1 , . . . , x n ), this is the desired formula.

Note that the expression (4.3) for P ν (Ω)(θ) complies with the so-called structure theorem for shape derivatives (see e.g. [51, §5.9]): P ν (Ω)(θ) depends only on the values of the normal component θ • n of the considered deformation θ on the boundary ∂Ω of the optimized shape.

Taking advantage of this expression, it is easy to infer an ascent direction for P ν (Ω) from a given shape Ω; for instance, selecting a deformation θ with the property (4.4) θ = F Ω ν n on ∂Ω, it stems from (4.2) that, for a small enough (pseudo-)time step τ > 0, the following expansion holds:

(4.5) P ν (Ω τ θ ) = P ν (Ω) + τ ∂Ω F Ω ν 2 ds + o(τ ), where o(τ ) τ τ →0 ---→ 0,
that is, Ω τ θ enjoys an improved performance with respect to P ν (Ω), provided F Ω ν does not vanish identically on ∂Ω (in which case Ω is already a critical point of P ν (Ω)).

Remark 4.2. In (4.4), an ascent direction for P ν (Ω) is identified as the gradient associated to the derivative P ν (Ω)(θ) via the L 2 (∂Ω) inner product; however natural, this choice suffers from at least two drawbacks. On the one hand, θ as in (4.4) only makes sense on the boundary ∂Ω of the considered shape Ω, while our numerical method for tracking the evolution of Ω requires that it is defined on the whole computational domain D (see Section 4.4 below). On the other hand, such a deformation θ may be quite irregular on ∂Ω (mathematically, it only belongs to L 2 (∂Ω)). This may induce numerical artifacts such as parasitic oscillations of the boundary of the optimized shape; see for instance [66, §6.2.4] for a description of this phenomenon.

A popular idea to overcome both issues consists in identifying a shape gradient for P ν (Ω) via another inner product a(•, •), acting on a (Hilbert) space V of vector fields defined on D as a whole and enjoying higher regularity than that of functions in L 2 (∂Ω). For instance, one may take

(4.6) V = H 1 (D) d , and a(u, v) = D (α 2 ∇u • ∇v + uv) dx,
where α stands for a regularization length scale. The identification of the shape gradient g Ω ∈ V of P ν (Ω) then amounts to the resolution of the variational problem:

Search for g Ω ∈ V s.t. ∀w ∈ V, a(g Ω , w) = P ν (Ω)(w),
which is achieved in practice by using a standard finite element method. An alternative choice to (4.6) is to use the inner product induced by the bilinear form of the linear elasticity system, as in e.g. [START_REF] Dapogny | Geometrical shape optimization in fluid mechanics using FreeFem++[END_REF]; we refer to [START_REF] Allaire | Shape and topology optimization[END_REF][START_REF] Burger | A framework for the construction of level set methods for shape optimization and reconstruction[END_REF][START_REF] Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF][START_REF] Dogan | Discrete gradient flows for shape optimization and applications[END_REF] about this procedure.

Numerical evaluation of the probabilities P ν (Ω)

The defining formulas (2.4) and (2.5) are unfortunately awkward when it comes to the numerical evaluation of the probabilities P ν (Ω). For instance, the straightforward use of (2.4) requires the calculation, at each quadrature point of the subset Ω ν × (R 3 \ Ω) n-ν of the 3n-dimensional space, of 2 n determinants of n × n matrices, which is totally impractical, even for moderately large values of n.

To circumvent this drawback, we rely on the algorithm proposed in [START_REF] Cances | How electrons guard the space: shape optimization with probability distribution criteria[END_REF] in the present case of a single determinant wave function. We briefly summarize this procedure in the context of the relaxed functional J ν (m), involving a general density function m ∈ L ∞ (R 3 , [0, 1]), which will prove useful in the next section.

Let S(m) ∈ C n×n be the overlap matrix, defined by: (4.7) ∀i, j = 1, . . . , n, S (m

) ij := σ∈{-1 2 , 1 2 } R 3 m(x)φ i x σ φ j x σ dx.
This matrix is obviously Hermitian, and it can be proved that its (real-valued) eigenvalues (λ i (m)) i=1,...,n completely determine J ν (m) via the relation:

(4.8) J ν (m) = Iν ⊆{1,...,n} # Iν =ν i∈Iν λ i (m)   j∈{1,...,n}\Iν (1 -λ j (m))   .
The calculation of J ν (m) can then be carried out along the following lines, for any integer ν = 0, . . . , n:

(1) calculate the entries of the overlap matrix S(m) in (4.7);

(2) diagonalize S(m) to obtain the eigenvalues λ i (m), i = 1, . . . , n;

(3) compute J ν (m) via the following iterative procedure: let the sequence of numbers p i j i=0,...,n j=0,...,i be defined by p 0 0 = 1, and:

(4.9) ∀i = 1, . . . , n,    p i 0 = (1 -λ i (m)) p i-1 0 , p i j = λ i (m) p i-1 j-1 + (1 -λ i (m)) p i-1 j , for j = 1, . . . , i -1, p i i = λ i (m) p i-1 i-1 .

It then turns out that: (4.10)

J ν (m) = p n ν . This computation procedure for J ν (m) is much more efficient than the aforementioned naive evaluation of the integral (2.4); the first and third steps can indeed be realized in O(n 2 ) operations, while the second one requires no more than O(n 3 ) operations if a standard diagonalization algorithm is used. Remark 4.3. Unfortunately, this algorithm is restricted to the case where Ψ is defined as a single Slater determinant, the generalization to the case of multi-determinant wave functions being far from obvious.

Numerical evaluation of F Ω ν

The numerical realizations of the shape gradient algorithm proposed in Section 4.4 and of the fixed point iteration procedure of Section 4.5 require the efficient evaluation of the quantity F m ν (x) in (2.13) at all points x in (a discretization of) the computational domain D. This is achieved thanks to a variation of the procedure described in Section 4.2, which was also originally proposed in [25, §2.2].

Let m ∈ L ∞ (R 3 , [0, 1]) be a given density function, and let S(m) be the associated overlap matrix, see (4.7). We denote by Λ(m) := diag(λ i (m)) i=1,...,n and X(m) ∈ C n×n the diagonal matrix made from its eigenvalues, and the unitary matrix gathering its eigenvectors, respectively:

(4.11) S(m)X(m) = X(m)Λ(m), and X T (m)X(m) = I n ,
where I n is the identity n × n matrix.

Assuming that all the eigenvalues of S(m) are simple, the implicit function theorem implies that the mappings m → λ i (m) and m → X(m) are Fréchet differentiable at m, and that:

(4.12) ∀h ∈ L ∞ (R 3 ), Λ (m)(h) = X(m) T S (m)(h)X(m),
where the derivative of S(m) is given by:

(4.13) S (m) ij (h) = σ∈{-1 2 , 1 2 } R 3 h(x)φ i x σ φ j x σ dx.
As we have seen in Section 4.2, the functional J ν (m) can be evaluated via the recursive procedure (4.9) and (4.10), in which all the quantities p i j are multivariate polynomial expressions in terms of λ 1 (m), . . . , λ n (m).

Hence, the chain rule yields:

(4.14) J ν (m)(h) = n l=1 q n ν,l λ l (m)(h),
where an elementary calculation reveals that the sensitivities q i ν,l := ∂p i ν ∂λ l can be computed recursively by setting q 0 0,l = 1 for all l = 1, . . . , n, and:

∀i = 1, . . . , n,    q i 0,l = α i,l (m) q i-1 0,l , q i j,l = β i,l (m) q i-1 j-1,l + α i,l (m) q i-1 j,l for j = 1, . . . , i -1, q i i,l = β i,l (m) q i-1 i-1,l , with the coefficients α l,l (m) = -1, β l,l (m) = 1, and for k = l, α k,l (m) = 1 -λ k (m), β k,l (m) = λ k (m).
Combining (4.12) to (4.14), we now obtain, for arbitrary h ∈ L ∞ (R 3 ),

J ν (m)(h) = R 3    n i,j=1 n l=1 q n ν,l X il (m)X jl (m)    σ∈{-1 2 , 1 2 } φ i x σ φ j x σ       h(x) dx;
comparing this expression with (3.3) and using the continuity of F m ν on R 3 , we eventually infer the practical expression of F m ν (x):

∀x ∈ R 3 , F m ν (x) = n i,j=1 n l=1 q n ν,l X il (m)X jl (m)    σ∈{-1 2 , 1 2 } φ i x σ φ j x σ    .
Thanks to this expression, once the overlap matrix (4.7) and its eigendecomposition (4.11) have been computed (both operations being already required for the evaluation of J ν (m)), the value of F m ν (x) can be calculated at any point x ∈ D in constant time.

Maximization of P ν (Ω) using the method of Hadamard and a level-set based mesh evolution method

The chief numerical technique involved in our practical resolution of the shape optimization problem (P s.o. ν ), is the level set based mesh evolution method from our previous work [START_REF] Allaire | Topology and geometry optimization of elastic structures by exact deformation of simplicial mesh[END_REF][START_REF]A mesh evolution algorithm based on the level-set method for geometry and topology optimization[END_REF][START_REF] Allaire | Shape optimization with a level set based mesh evolution method[END_REF] (see also [START_REF] Feppon | Shape optimization of a coupled thermal fluidstructure problem in a level set mesh evolution framework[END_REF][START_REF] Feppon | Topology optimization of thermal fluid-structure systems using body-fitted meshes and parallel computing[END_REF] for recent applications). The latter features an explicit (meshed) description of any shape arising in the course of the optimization process, and still allows for arbitrarily large deformations of the latter (including topological changes). In a nutshell, we rely on two complementary descriptions of Ω:

• A meshed description: Ω is meshed exactly. More precisely, the mesh T of the computational domain D is arranged so that (i) it is conforming in the sense of finite element computations;

(ii) it contains an explicit discretization of the actual shape Ω as a submesh;

(iii) the regions of D where sharp variations of the spin orbitals φ i occur contain smaller elements. This meshed representation of shapes allows for an accurate calculation of the coefficients of the overlap matrix (4.7) and of the probabilities P ν (Ω); see Fig. 1 (left) for an illustration. • A level set description: following the seminal work [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF] dedicated to the level set method (see also [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF][START_REF] Wang | A level set method for structural topology optimization[END_REF] about its introduction in the context of shape optimization), Ω is the negative subdomain of a "level set function" ϕ : D → R, that is:

(4.15)    ϕ(x) < 0 if x ∈ Ω, ϕ(x) = 0 if x ∈ ∂Ω, ϕ(x) > 0 if x ∈ D \ Ω;
see Fig. 1 (right). This description is particularly convenient when it comes to updating the shape; indeed, if Ω(t) is a domain evolving according to a velocity field with normal component v(t, x) over a time period (0, T ), it is well-known that an associated level set function ϕ(t, x) (i.e. (4.15) holds at any time t) is obtained as the solution to the following Hamilton-Jacobi equation: In practice, (4.16) is solved on the simplicial mesh of D by means of a semi-implicit scheme very much similar to that of [START_REF] Strain | Semi-lagrangian methods for level set equations[END_REF]; see [START_REF] Bui | An accurate anisotropic adaptation method for solving the level set advection equation[END_REF]. We refer to [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF][START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics[END_REF] for more exhaustive presentations of the stakes and capabilities of the level set method. Efficient algorithms are used in order to switch from one representation of a shape Ω ⊂ D to the other:

• When Ω is explicitly discretized in the mesh T of D, one level set function ϕ : D → R for Ω is computed as the signed distance function d Ω to Ω, that is:

∀x ∈ D, d Ω (x) =    -d(x, ∂Ω) if x ∈ Ω, 0 if x ∈ ∂Ω, d(x, ∂Ω) if x ∈ D \ Ω,
where d(x, ∂Ω) = min y∈∂Ω |x -y| is the usual Euclidean distance from x to ∂Ω. The calculation of d Ω is realized thanks to the open-source library mshdist [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF].

• From the datum of a level set function ϕ : D → R for Ω, in practice defined at the vertices of an arbitrary mesh T of D, a new mesh T of D is created in which Ω appears explicitly as a submesh, and which is adapted to a desired local size feature, by using the open-source remeshing library mmg from [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF]; see [START_REF] Dapogny | Mmg[END_REF]. The shape optimization algorithm combining the ingredients of this section and the previous ones is summarized in Algorithm 1.

Remark 4.4. The Algorithm 1 is local, insofar as it is only guaranteed to converge to a local maximizer of (P s.o. ν ), which (strongly) depends on the initial shape Ω 0 . As we have mentioned, and as we shall see throughout Section 5, these local solutions (and not only the global maximizer of (P s.o. ν )) are of much interest from the chemical point of view, since they account for stable, plausible electronic configurations.

Noting in addition that if ϕ is one level set function for a domain Ω, then so is cϕ for any constant c > 0, the optimality conditions (4.17) can be enforced by seeking a function ϕ : D → R such that:

ϕ L 2 (D) = 1, ϕ = g Ω ν , where Ω = {x ∈ D s.t. ϕ(x) < 0} and g Ω ν := - 1 F Ω ν L 2 (D) F Ω ν .
We now apply a fixed-point algorithm with relaxation to construct such a function ϕ. Starting from an initial shape Ω 0 , and a corresponding level set function ϕ 0 with ϕ 0 L 2 (D) = 1, a sequence of shapes Ω k and associated level set functions ϕ k is calculated in the following way:

(4.18) ϕ k+1 = 1 ϕ k+1 L 2 (D) ϕ k+1 , where ϕ k+1 := (1 -γ k )ϕ k + γ k g k , g k ≡ g Ω k ν
, and γ k ∈ (0, 1) is a suitable relaxation step. In practice, we rely on a slightly more efficient version of this algorithm, based on the so-called spherical linear interpolation procedure [START_REF] Shoemake | Animating rotation with quaternion curves[END_REF], which conveniently alleviates the need to impose explicitly the normalization of ϕ k+1 in (4.18). Introducing a k := arccos(ϕ k , g k ) L 2 (D) , we observe that:

ϕ k+1 2 L 2 (D) = (1 -γ k ) 2 + γ 2 k + 2γ k (1 -γ k )(ϕ k , g k ) L 2 (D) . Defining the "new" pseudo time-step τ k > 0 via the relation sin(τ k a k ) sin a k := γ k (1 -γ k ) 2 + γ 2 k + 2γ k (1 -γ k ) cos a k 1/2 ,
an elementary calculation yields:

sin((1 -τ k )a k ) sin a k = (1 -γ k ) (1 -γ k ) 2 + γ 2 k + 2γ k (1 -γ k ) cos a k 1/2 .
Hence, the updated level set function ϕ k+1 in (4.18) is obtained directly as:

ϕ k+1 = 1 sin a k sin((1 -τ k )a k )ϕ k + sin(τ k a k )g k .
This fixed-point iteration procedure for the resolution of (P s.o. ν ) is summarized in Algorithm 2.

Algorithm 2 Fixed-point iteration algorithm for the resolution of (P s.o. ν ). Initialization: Level set function ϕ 0 for the initial shape Ω 0 such that ϕ 0 L 2 (D) = 1.

for k = 0, 1, 2, . . . until convergence do

(1) Calculate the gradient g k := -

1 F Ω k ν L 2 (D) F Ω k ν . (2) Calculate a k := arccos(φ k , g k ) L 2 (D) ; (3) Choose a relaxation step 0 < τ k < 1; (4) The updated level set function ϕ k+1 is ϕ k+1 = 1 sin a k sin((1 -τ k )a k )ϕ k + sin(τ k a k )g k ,
corresponding to the updated shape Ω k+1 = {ϕ k+1 < 0}. end for return Optimized shape Ω k .

In principle, Algorithm 2 could be used as is to solve (P s.o. ν ). Unfortunately, like most fixed-point strategies, it may become unstable if the initial guess Ω 0 is "too far" from the expected local maximizer (see Section 5 and especially Section 5.1). For this reason, in practice, we combine Algorithm 2 with our main shape optimization Algorithm 1, by periodically interrupting the process of the latter (say, every 10 iterations) to perform one iteration of the procedure in Algorithm 2.

Remark 4.5.

• The fixed-point Algorithm 2 is very similar to that proposed in [START_REF] Amstutz | Analysis of a level set method for topology optimization[END_REF][START_REF] Amstutz | A new algorithm for topology optimization using a level-set method[END_REF] for the optimization of a shape functional using its topological derivative, once we have observed that F Ω ν is exactly the topological derivative of the shape functional P ν (Ω); indeed, a variation of the proof of (3.3) reveals that for any x ∈ Ω,

P ν (Ω \ B(x, r)) = P ν (Ω) -r 3 |B(0, 1)|F Ω ν (x) + o(r 3 ), where |B(0, 1)| = 4π
3 is the measure of the unit three-dimensional ball. We refer to [START_REF] Novotny | Topological derivatives in shape optimization[END_REF] for more details about topological derivatives and to [START_REF] Allaire | Structural optimization using topological and shape sensitivity via a level set method[END_REF][START_REF]Topological sensitivity analysis and applications in shape optimization[END_REF][START_REF] Burger | Incorporating topological derivatives into level set methods[END_REF] about their use in shape and topology optimization. It follows from this remark that the use of Algorithm 2, and particularly its combination with the boundary variation Algorithm 1, can be understood as a complementary means to appraise the sensitivity of the optimized function P ν (Ω) with respect to topological changes, which echoes to the "hybrid" strategies in shape and topology optimization based on shape and topological derivatives, see again [START_REF] Allaire | Structural optimization using topological and shape sensitivity via a level set method[END_REF].

• This similarity between the shape derivative, the topological derivative and the derivative with respect to density functions for a shape functional and its relaxed counterpart is not a surprise, and a connection between those three concepts was established in a more general context in [START_REF] Amstutz | A consistent relaxation of optimal design problems for coupling shape and topological derivatives[END_REF].

Construction of the computational domain D

The selection of the bounded domain D where all computations take place in numerical practice is not a completely straightforward task in the present setting, where the sought optimized shapes may be unbounded.

It is therefore of utmost importance to ensure that D enclose the n electrons of the considered molecule with a probability close to 1.

To achieve this, we start from a very large cube around the nuclei of the considered molecule, which is equipped with a coarse mesh. Since the problem of maximizing the probability P n (Ω) to find all n electrons inside Ω has the full space R 3 for global maximizer, a first use of Algorithm 1 for this problem yields an increasing sequence of shapes Ω k , k = 1, . . ., with P n (Ω k ) → 1 as k → ∞. We interrupt the procedure at the stage k 0 when the probability P n (Ω k0 ) exceeds 1 -ε hold-all (in practice, ε hold-all = 0.005). The resulting domain D := Ω k0 is then meshed thanks to the numerical tools presented in Section 4.4, and yields a suitable computational box for the resolution of the problem (P s.o. ν ) of interest.

Numerical results

In this section, we provide several calculation examples of Maximum Probability Domains which illustrate the capabilities of the proposed numerical methods, and also suggest interesting comments from the chemical point of view. Most of our numerical examples focus on the calculation of MPDs associated to the probability P 2 (Ω) of finding ν = 2 electrons. Indeed, as we have hinted at in the introduction, the historic Lewis model from classical chemistry predicts that the valence electrons of a molecule tend to arrange in pairs; three categories of these structures are usually distinguished, corresponding to different locations within the molecule and different chemical properties:

• the core electron pairs are attached to one specific atom; they are chemically inert.

• The valence lone pairs are also attached to a particular atom, but they may be chemically reactive and responsible for specific properties of the molecule. • The chemical bonds are electron pairs shared by two atoms (and sometimes more in the case of delocalised bonding). Hence, MPDs associated to the value ν = 2 are the quantum counterparts of the electron pairs featured by the historic Lewis model. Note however that the instances of the optimization problem (P s.o.

ν ) where ν = 1 or ν > 2 are also pertinent from the chemical viewpoint, since the corresponding MPDs may be interpreted as the locations of radical electrons, and multiple bonds with delocalized electrons, respectively.

In the following examples, we aim to identify these Lewis pairs as local solutions to the shape optimization problem (P s.o. ν ) for ν = 2. Even though Theorem 2.5 and the subsequent remarks guarantee the existence of global maximizers in the cases under scrutiny, all local maximizers are plausible, stable configurations for electron pairs, and are therefore of great physical relevance. Our numerical strategy relies on the coupling between the boundary variation Algorithm 1 and the fixed-point Algorithm 2 described in Sections 4.4 and 4.5: using initial shapes with different locations within the considered molecule, the numerical resolution of (P s.o.

ν ) is driven towards MPDs corresponding to as many different pairs of the considered electron cloud. Our first two examples in Sections 5.1 and 5.2 are devoted to the computation of several MPDs in the context of quite simple and well-understood molecules from the chemical point of view, namely HF and H 2 O, respectively; we then consider in Section 5.3 the ethylene molecule C 2 H 4 , which displays a double carbon (carbon bond). Finally, in Sections 5.4 and 5.5, we apply our methods in the contexts of the dicarbon C 2 and of the larger propellane molecule C 5 H 6 (in terms of the total number of electrons), respectively, with the aim to try and gain insight about the locations of particular electron pairs which are the subject of active debates within the chemical community.

Before proceeding, let us warn the reader that the chemical interpretations raised by our computations should somehow be tempered, since they rely on drastic simplifications of the realistic physical model, notably because of the use of single determinant wave functions; see Remark 2.2 and the perspectives outlined in Section 6.

The hydrogen-fluorine molecule HF

This first example aims to illustrate the main features of the shape and topology optimization Algorithms 1 and 2 in the context of the simple hydrogen-fluorine molecule H-F. This molecule possesses n = 10 electrons in total, which are organized in five pairs, according to the Lewis model: one of them corresponds to the core electrons of the fluorine atom F, three are valence lone pairs of the latter, and the fifth pair is the chemical bond between the hydrogen and fluorine atoms.

We expect that each of these pairs be associated to a Maximum Probability Domain, characterized as a local maximizer of (P s.o.

ν ) for ν = 2. Following chemical intuition, the MPD associated with the core electron pair should be a small, almost spherical domain centered around the fluorine nucleus. Since the wave function Ψ is that of an isolated molecule in vacuum, we expect that the other four valence MPDs be unbounded in R 3 .

Here, and throughout the numerical examples of this article, Ψ is built as a single Slater determinant of the form (2.7) in the restricted Hartee-Fock framework (see Section 2.2); the constituent molecular orbitals φ i in (2.8) are supplied by the standard def2-TZVPP basis sets [START_REF] Weigend | Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for h to rn: Design and assessment of accuracy[END_REF] taken from the Basis Set Exchange database [START_REF] Pritchard | New basis set exchange: An open, up-to-date resource for the molecular sciences community[END_REF]. In the present context of the HF molecule, Ψ involves n = 10 molecular orbitals φ i , each of them being described by means of R = 59 Gaussian primitives.

We focus our attention on the calculation of the MPD associated with the bond between the H and F atoms, which is the electron pair closest to H. In this context, we select the sphere whose radius equals 0.58 times the interatomic distance, centered on the hydrogen atom by way of the initial guess Ω 0 (see Fig. 2 top left). We then experiment several algorithmic strategies to solve (P s.o. ν ), based on the numerical methods presented in Section 4:

• Case 1: the boundary variation Algorithm 1 is used as is;

• Case 2: the fixed-point Algorithm 2 is used in a standalone fashion;

• Case 3: the combination of Algorithms 1 and 2 is applied. In this section, and throughout this article, our numerical simulations are conducted on a standard HP Desktop computer tower equipped with Ubuntu 16.04, 16G of RAM, and Intel Xeon E5 -2665 2.4/3.1 GHz with 8 physical threads (16 in hyper-threading).

The optimized shapes resulting from these experiments are displayed in Fig. 2; the details of the optimization processes, and the associated convergence histories are reported in Table 2 and Fig. 3, respectively. In all three cases, the average number of vertices (resp. tetrahedra) in the featured meshes of the computational domain D is 35,000 (resp. 200,000). Obviously, the optimization procedure has not converged in Case 1; as can be seen from Fig. 3 and the values of the probability P 2 (Ω) reported in Table 2, which are significantly lower than those found in our next experiments. The sole application of a shape gradient algorithm is not able to capture regions far from the nuclei of the molecule. This can be explained by the very large contrast between the "reasonably large" values taken by the shape gradient of P 2 (Ω) inside the interatomic regions (of the order of 10 -4 ), when compared to its very small magnitude (of the order of 10 -8 , i.e. below the orders of magnitude of the roundoff and numerical errors) far away from it. The procedure conducted in Case 2 blows up very rapidly (the resulting shape is not reported): as often in the practice of fixed-point iteration strategies, the choice of an initial guess Ω 0 which is "too far" from the expected result makes the procedure totally unstable. On the contrary, the strategy of Case 3, resulting in the shape displayed in Fig. 2 (right) reveals a satisfying behavior: the addition of fixed point steps to the boundary variation Algorithm 1 allows to circumvent the difficulties encountered in Cases 1 and 2. It enables the rapid capture of an unbounded domain (only 14 iterations are needed for the whole computation), resulting in a shape with better performance P 2 (Ω) (see again Table 2). As a complementary assessment of the physical relevance of the obtained MPDs Ω, we calculate their population, that is, the expectation of the number of electrons contained in them:

(5.1) n ν=0 νP ν (Ω).

In the case of the MPD obtained in Case 3, the value of this population is very close to 2, even though no particular effort is paid in enforcing this feature during our computations. This is further confirmation that this domain is a plausible location for the chemical bond between both atoms of the HF molecule.

Case

The water molecule H 2 O

This section is devoted to the water molecule H 2 O, whose chemical behavior is also well understood. According to the Lewis theory, the n = 10 electrons of this molecule are arranged in five electron pairs: one core pair and two valence lone pairs are attached to the oxygen nucleus, and each hydrogen atom is connected to the oxygen via a chemical bond. We search for the Maximum Probability Domains of the H 2 O molecule accounting for these five pairs. Again, this goes through the numerical resolution of the shape optimization problem (1.1) for ν = 2, thanks to our combination of Algorithms 1 and 2, using different initial shapes Ω 0 . The restricted Hartee-Fock wave function Ψ describing the present electron cloud is a single Slater determinant of the form (2.7), involving n = 10 spin orbitals φ i , each of them being of the form (2.8) with R = 67 Gaussian primitives.

We first look for the MPD associated to the core electron pair of the molecule, and we investigate two different scenarii as regards the initial shape Ω 0 :

• Case 1: Ω 0 is a ball with radius 0.55 times the distance O-H, centered at the oxygen atom; see Fig. 4 (top, left). • Case 2: Ω 0 is a ball with radius 0.55 times the distance O-H, whose center is shifted away from the oxygen atom, in order to investigate the stability of the centered configuration; see Fig. 4 (bottom, left). In both cases, the average number of vertices of the meshes of the computational domain D is 35, 000 (corresponding to, roughly, 200, 000 tetrahedra). The optimized shapes obtained in both situations are displayed on Fig. 4: as expected, in Case 1, where Ω 0 is very close to the expected MPD, only very few iterations (8) of our numerical procedure are required for convergence; the total computation takes about 10 mn. On the contrary, in Case 2, where Ω 0 lies a little farther, the calculation requires a few more iterations (19 in total), corresponding to a computational time of about 45 mn.

We next turn to the identification of the four valence electron pairs of the water molecule. The four initial shapes used in this context are depicted on Fig. 5 (upper row), and the corresponding results are displayed in Fig. 5 (bottom row), which may be grouped into two sets of symmetrically related solutions: the optimized shapes on the two leftmost columns in Fig. 5 correspond to the oxygen lone pairs, and the optimized shapes on the two rightmost columns are the O-H bonds of the molecule. The final values of the probability function P 2 (Ω) agree with those obtained by other numerical strategies in [START_REF] Lopes | Understanding Maximum Probability Domains with Simple Models[END_REF][START_REF] Menendez | A view of covalent and ionic bonding from maximum probability domains[END_REF][START_REF] Scemama | Maximum probability domains from quantum monte carlo calculations[END_REF][START_REF] Turek | Bonding in heavier group 14 zero-valent complexes-a combined maximum probability domain and valence bond theory approach[END_REF]: we obtain P 2 (Ω) = 0.416 for the MPDs associated to the oxygen lone pair and P 2 (Ω) = 0.453 for those corresponding to the bonds. As in the case of the H-F molecule tackled in Section 5.1, our numerical strategy manages to capture unbounded optimized shapes, with several connected components, even when a bounded and connected initial guess is used. Let us emphasize that this is a salient feature of the proposed coupling of the boundary variation Algorithm 1 (which does not allow such dramatic topological changes) with the fixed point iteration strategy of Algorithm 2. Note that the smaller connected components featured by these MPDs are not numerical artifacts from the calculation. Indeed, we have verified that their removal from the obtained optimized shape leads to a slight decrease in the value of the probability P 2 (Ω). In a similar spirit, solving (P s.o. ν ) anew, starting from these versions of the optimized shapes where the small connected components are removed, leads to their restoration after a few iterations.

Interestingly, these calculations suggest a possible answer about an old, but still lively debate between two rival visions for the water lone pairs (i.e. the two valence pairs attached to the oxygen atom) [START_REF] Clauss | Rabbit-ears hybrids, vsepr sterics, and other orbital anachronisms[END_REF][START_REF] Hiberty | Pleading for a dual molecular-orbital/valence-bond culture[END_REF][START_REF] Hiberty | Comment on "rabbit-ears hybrids, vsepr sterics, and other orbital anachronisms". a reply to a criticism[END_REF]. On the one hand, the prevailing representation within the organic chemistry community, provided by the so-called "hybridization theory", features two equivalent "rabbit-ear" lone pairs, with identical shapes; see Fig. 6 (left). On the other hand, the Molecular Orbital reading of the H 2 O molecule predicts two different lone pairs, one of so-called σ type, and the other of so-called π type; see Fig. 6 (right). In this perspective,

The ethylene molecule C 2 H 4

We now study the ethylene molecule C 2 H 4 , which features n = 16 electrons, and we more particularly focus on the double C=C bond lying between the two carbon atoms. Similarly to what happens for multiply bonded hydrocarbons in general, it was originally suggested in [START_REF] Pauling | The nature of the chemical bond. application of results obtained from the quantum mechanics and from a theory of paramagnetic susceptibility to the structure of molecules[END_REF] that this double bond corresponds to two bent bonds with similar shapes, also called banana bonds [START_REF]The Nature of the Chemical Bond[END_REF]. Yet, another description featuring two different bonds of σ and π types has become increasingly popular within the chemical community [START_REF] Hueckel | Quantum contributions to the benzene problem[END_REF][START_REF]Quantum contributions to the problem of aromatic and unsaturated compounds. 3[END_REF]; see Fig. 7 for a schematic representation of both orbital sets.

In order to put both hypotheses to the test, we use our combination of Algorithms 1 and 2 to search for the Maximum Probability Domains of the ethylene molecule in the C-C interatomic region, leaving aside the core pairs. Again, this means that we perform several resolutions of (P s.o.

ν ) for ν = 2, starting from multiple initial guesses situated in this region. In all the conducted experiments, the average number of vertices in the computational mesh of D is 70, 000 (and so the average number of tetrahedra is about 375, 000).

Our first experiment features a sphere located in the C-C interatomic axis or slightly above, with radius 0.79 times the distance between carbon atoms, by way of initial guess Ω 0 ; see Fig. 8 (left). The optimized domain Ω resulting from our combination of Algorithms 1 and 2 is that displayed on Fig. 8 (right). Judging from its symmetry with respect to the plane of the molecule, it is tempting to associate the MPD Ω to a C-C bond of σ type (again, see Fig. 7, (right)). However, despite many attempts, we did not find any corresponding complementary MPD that would correspond to a π-bond between these two carbon atoms.

In a second set of experiments, using initial domains Ω 0 situated from either side of the plane of symmetry of the molecule (see Fig. 9, left column) leads to two MPDs in the C-C interatomic region, which are represented on Fig. 9 (right column). These could correspond to two equivalent "banana bonds", which are symmetric with respect to the plane of the molecule. Besides, the final value 0.404 of the probability P 2 (Ω) for either of these two MPDS Ω is slightly larger than that 0.401 obtained in the first experiment, in the case of the potential σ-bond MPD of Fig. 8.

All things considered, both arguments suggest that the "banana bond" model may be the most plausible structure for the two bonds between the carbon atoms of the ethylene.

The quadruple bond of the dicarbon dimer C 2

In this section, we conduct investigations in the context of the dicarbon molecule C 2 . It contains n = 12 electrons, arranged in two core pairs (one attached to each carbon atom), and four bonds, whose configuration is the main focus of this section, since it has aroused lively debates within the chemical community lately [START_REF] Frenking | Critical comments on "one molecule, two atoms, three views, four bonds?[END_REF][START_REF]Comment on "the quadruple bonding in c-2 reproduces the properties of the molecule[END_REF][START_REF] Hermann | The chemical bond in c-2[END_REF][START_REF] Shaik | The quadruple bonding in c-2 reproduces the properties of the molecule[END_REF][START_REF]A response to a comment by g. frenking and m. hermann on: "the quadruple bonding in c-2 reproduces the properties of the molecule[END_REF][START_REF] Shaik | Quadruple bonding in c-2 and analogous eight-valence electron species[END_REF][START_REF] Shaik | One molecule, two atoms, three views, four bonds?[END_REF]:

(1) Contrary to the common belief that two atoms within the same group of the periodic table are incapable of forming more than three bonds together, chemical arguments exposed in [START_REF] Shaik | Quadruple bonding in c-2 and analogous eight-valence electron species[END_REF] vouch for a quadruple bonding situation in the case of the dicarbon dimer C 2 , a setting which we shall refer to as the C≡C model in the sequel. A tentative explanation for this peculiar bonding was proposed in [START_REF] Shaik | The quadruple bonding in c-2 reproduces the properties of the molecule[END_REF], relying on Valence Bond analysis: the dimer C 2 would feature three standard covalent bonds, and a fourth, weaker C-C bond, of particular "inverted" type, resulting from the interaction of two atomic orbitals located near the two different carbon atoms, oriented in opposite directions, and nevertheless interacting with each other. Contrary to traditional covalent bonds, this type of bonding can occur without significant overlap between the orbitals at stake: it is enabled by the so-called Taking advantage of our numerical algorithms, we aim to gain insight about the relevance of these three visions. As usual, we apply our combination of Algorithms 1 and 2 to identify the local maximizers of the shape optimization problem (P s.o.

ν ) for ν = 2, starting from different initial shapes. In all the forthcoming experiments, the computational mesh of D contains about 70, 000 vertices (for about 375, 000 tetrahedra).

We start by assessing the scenarii featuring three distinct bonds between the two carbon atoms. Since the C 2 molecule has rotational symmetry with respect to the C-C axis, we rely on three different initial guess, which are as many balls Ω 0 with common radius 0.85 times the interatomic distance, whose centers are slightly shifted from the C-C axis, with a rotational symmetry of ±120 • . This results in the three MPDs displayed in Figs. 10 and11, which have the form of "banana bonds", like those found in Section 5.3, in the case of the ethylene molecule. In order to verify that the collection of these three MPDs does correspond to three distinct bonds, we calculate the population (5.1) contained in their reunion. The resulting number is 5.6, which is fairly close to the expected value 6. This seems to support the fact that there are three classical covalent bonds between the two carbon atoms, thus ruling out out the |C=C| model.

To proceed, we then search for an MPD associated to ν = 2 electrons outside the interatomic CC region. To this end, we rely on the initial shape Ω 0 depicted on Fig. 11 (top, left), which is "far" from this region. At this stage, we have identified four disjoint MPDs for the C 2 molecule (leaving aside the core pairs). This collection is a priori compatible with both the quadruple bond C≡C and the diradical •C≡C• representations of the four bonds. In order to come to a decision about the relevance of both models, we perform one last experiment: we search for an MPD associated to ν = 1 electron, outside the interatomic region, using one of the two connected components of the domain Ω 2 ext as initial guess. The rationale behind this test is that, if according to the •C≡C• model, the remaining two valence electrons of C 2 were independent, each being ext , as depicted in Fig. 11: it also has two connected components, with a lower total volume. From the chemical point of view, this result suggests that the two electrons of C 2 lying outside the interatomic region cannot be really separated as two independent electrons; rather, they ought to be seen as the constituent of a genuine electron pair. All things considered, our MPD analysis would comfort the quadruple bond vision C ≡ C between the two carbon atoms of the C 2 dimer, with three central standard "banana" bonds, and a fourth "inverted" bond.

The propellane molecule C 5 H 6

We conclude this numerical section with a more challenging calculation, arising in the context of the [1,1,1] propellane molecule C 5 H 6 , which contains n = 36 electrons; see Fig. 12 for an illustration. This investigation is motivated by a conjecture from the chemical community regarding the chemical bond between the two bridgehead carbon atoms of this polycyclic molecule (i.e. the two axial carbons which do not belong to one of the CH 2 groups); see [START_REF] Feller | Abinitio studies of 1.1.1 propellane and 2.2.2 propellane[END_REF][START_REF] Jackson | The c1-c3 bond in 1.1.1 propellane[END_REF]. For this molecule, the standard Def2-SVP basis set has been used for the molecular orbitals φ i in the definition (2.7) of the wave function Ψ. This central bond is certainly not conventional in nature: practical experiments have indeed evidenced its strength, which far exceeds that of a classical bond between two atoms lying at such distance from one another. The mechanism of "inverted" bond has been proposed in [START_REF] Wu | The inverted bond in 1.1.1 propellane is a charge-shift bond[END_REF] to account for this phenomenon, like in the case of the C ≡ C model for the C 2 dimer discussed in Section 5.4. This interpretation has been challenged recently, leading to a revival of the controversy [START_REF] Braida | Comment on "the 'inverted bonds' revisited. analysis of 'in silico' models and of 1.1.1 propellane using orbital forces[END_REF][START_REF] Laplaza | The "inverted bonds" revisited: Analysis of "in silico" models and of 1.1.1 propellane by using orbital forces[END_REF].

Thanks to our numerical Algorithms 1 and 2, we search for an MPD for ν = 2 electrons, related to an "inverted" bond between the two central carbons: inspired from the case of the C 2 dimer in Section 5.4, we start from an initial guess composed of two disjoint spheres, one near each carbon atom; see Fig. 12 (left). On average, the meshes of the computational domain D involved in the process contain 170, 000 vertices (corresponding to about 1, 000, 000 tetrahedra), and the total calculation takes about 8h.

As displayed in Fig. 12 (right), the resulting optimized shape resembles much the MPD Ω 2 ext found in Section 5.4. Like in there, the outlook of this optimized shape could be explained either by the existence of a two-electron "inverted" bond between the bridgehead carbon atoms of the propellane molecule, or alternatively by the presence of two radical electrons, one being attached to each of these carbon atoms.

In order to explore further both hypotheses, we rely on the same strategy as in Section 5.4: we search for MPDs associated to ν = 1 electron in the vicinity of these two carbon atoms. Contrary to the case of the C 2 molecule, this search results in two types of MPDs, depending on the initial guess used in the practice of our algorithmic strategy. On the one hand, two MPDs for ν = 1 were found, localized near either carbon atom; on the other hand, another, "delocalized" MPD was found, which very much resembles that obtained in the case of the C 2 dimer; see Fig. 13. Both types of solutions are associated with fairly large values of the probability functional P 1 (Ω): this quantity equals 0.450 for each of the two aforementioned "localized" MPDs, and 0.429 for the "delocalized" MPD.

The existence of a "delocalized" MPD for ν = 1, associated with a large value of the quantity P 1 (Ω) tends to support the presence of a genuine inverted bond in the propellane molecule. Let us however emphasize once more that, like in the case of the C 2 dimer dealt with in Section 5.4, these conclusions would have to

3 . 1 . Preliminary remarks 3 . 2 .) 3 . 3 .) 3 . 4 .

 31323334 Existence of a solution to the relaxed problem (P relax ν Optimality conditions for the relaxed problem (P relax ν End of the proof of Theorem 2.5 4. A shape and topology optimization algorithm for the calculation of Maximum Probability Domains 4.1. Computation and use of shape derivatives 1 Institute of Engineering Univ. Grenoble Alpes

1 4. 2 . 4 . 3 . 4 . 4 . 4 . 5 . 6 .

 124344456 Numerical evaluation of the probabilities P ν (Ω) 16 Numerical evaluation of F Ω ν 17 Maximization of P ν (Ω) using the method of Hadamard and a level-set based mesh evolution method 18 Addition of a fixed point strategy 20 4.6. Construction of the computational domain D 22 5. Numerical results 22 5.1. The hydrogen-fluorine molecule HF 23 5.2. The water molecule H 2 O 25 5.3. The ethylene molecule C 2 H 4 27 5.4. The quadruple bond of the dicarbon dimer C 2 28 5.5. The propellane molecule C 5 H 6 31 Conclusion and open problems 33 References 34 1. Introduction

  (i) The relaxed problem (P relax ν ) has a (possibly non unique) global solution in L ∞ (R 3 , [0, 1]). (ii) Any global solution m * to the relaxed problem (P relax ν ) is of the form m * = 1 Ω * where Ω * ∈ M. (iii) The shape optimization problem (P s.o. ν ) has a (possibly non unique) global solution in M. (iv) Any global solution Ω * to (P s.o. ν

( 2 )

 2 k → 0 strongly in L 1 (X) as k → ∞, which terminates the proof. Proposition 3.4. Assume that the wave function Ψ given by (2.1) satisfies Assumption 2.1. Then, the functional J ν (m) defined by (2.11) is sequentially continuous for the weak- * convergence of L ∞ (R 3 , [0, 1]).

ν

  ) is now a simple consequence of Proposition 3.4. Corollary 3.5. The relaxed problem (P relax ν ) has one global solution.

  at a point x 0 which is at the same time a Lebesgue point of F m * ν and a point in {m * = 0} with density one, we obtain: For a.e. x ∈ R 3 s.t. m * (x) = 0, F m * ν (x) 0. Finally, repeating again the same procedure, it follows: For a.e. x ∈ R 3 s.t. 0 < m * (x) < 1, F m * ν (x) = 0. The complete statement of the proposition results from these facts since the sets {F m * ν < 0}, {F m * ν > 0} and {F m * ν = 0} are a partition of R 3 , as well as the sets {m * = 0}, {m * = 1} and {0 < m * < 1} (up to a set with null Lebesgue measure). Remark 3.8. A similar argument shows that if Ω ∈ M is one local maximizer for the shape optimization problem (P s.o.

  x) + v(t, x)|∇ϕ(t, x)| = 0 for (t, x) ∈ (0, T ) × D, ϕ(0, •) is one level set function for Ω(0).

Figure 1 .

 1 Figure 1. (Left) One mesh of D, featuring an explicit discretization of the shape Ω (in black); (right) Level set representation ϕ of the shape Ω at the vertices of the mesh T .

Figure 2 .

 2 Figure 2. Maximum Probability Domains for ν = 2 associated with the chemical bond pair of the H-F molecule obtained using different algorithmic strategies. The H-F molecule is represented by the usual "ball and stick" model, with green/white ball representing the location of the fluorine/hydrogen atoms respectively, while the bicolour stick stands for the chemical bond between them. (Left) initial shape; optimized domains in (middle) Case 1, and (right) Case 3.

Figure 3 .

 3 Figure 3. Convergence histories of the numerical experiments conducted in Section 5.1: the dashed and dotted lines are respectively associated to the standard boundary variation Algorithm 1. The continuous line is associated to the combination of Algorithms 1 and 2.

Figure 4 .

 4 Figure 4. Maximum Probability Domains attached to the core electron pair of the oxygen atom in the H 2 O molecule. The oxygen atom is represented in red, and the two hydrogen atoms are in white, at the ends of the bicolor stick representing the two O-H bonds. The same final domain is obtained starting from two different initializations Ω 0 : (upper row, from left to right) iterations 0, 2 and 8 when Ω 0 is centered on the oxygen atom; (lower row, from left to right) iterations 0, 11, and 19 when Ω 0 is shifted away from the oxygen atom.

Figure 5 .

 5 Figure 5. Maximum Probability Domains associated with the valence electron pairs in the H 2 O molecule; the initial domains are displayed on the upper row, and the optimized ones are on the bottom row. The first two examples on the left correspond to MPDs associated with the oxygen lone pairs, while the other two on the right are the MPDs associated with the O-H chemical bond.

Figure 6 .

 6 Figure 6. (Left) In the H 2 O molecule, the "rabbit-ear" description of the water lone pairs predicts that the oxygen atom has two lone pairs with equivalent shapes, which are symmetric with respect to the plane of the molecule; (right) the σ-π vision features two different electron pairs, one of σ-type (comprised in the plane of the molecule), the other of π-type (which is orthogonal to it).

Figure 7 .

 7 Figure 7. (Left) the "bent-bond" model for the double C=C bond of the ethylene molecule features two electron pairs which are symmetric with respect to the plane of the molecule; (right) the σ -π description of this bond involves one σ electron pair within the molecule plane, and an orthogonal π electron pair.

Figure 8 .

 8 Figure 8. Maximum Probability Domain associated to the potential σ bond pair between the two carbon atoms of the ethylene molecule, as discussed in Section 5.3 (no corresponding π MPD has been found). The carbon atoms are in grey, and the hydrogen ones are in white: (left) initial guess; (right) optimized shape.

Figure 9 .

 9 Figure 9. Maximum Probability Domains associated with the C-C bond pairs in the ethylene molecule C 2 H 4 , investigated in Section 5.3. The carbon atoms are represented in grey, and the hydrogens are in white: (left) initial guess; (right) MPDs interpreted as banana bonds.

Figure 10 .

 10 Figure 10. Representation of the three MPDs for ν = 2 associated to the "banana bonds" of the C 2 dimer considered in Section 5.4.

Figure 11 .

 11 Figure 11. (Upper row) Maximum Probability Domain for ν = 2, associated with the "inverted" C-C bond at stake in Section 5.4; (left) initial guess and (right) optimized shape. (Lower row) maximum Probability Domain for ν = 1; (left) initial guess and (right) optimized shape.

Figure 13 .

 13 Figure 13. Calculation of Maximum Probability Domains for ν = 1 in the case of the propellane molecule C 5 H 6 investigated in Section 5.5. (Upper row) case where the initial guess is made of one connected component of the MPD for ν = 2 of Fig. 12; (left) initial guess; (right) resulting, "localized" MPD for ν = 1. (Lower row) case where the initial guess is made of two spheres, located near the two central carbons; (left) initial guess; (right) resulting, "delocalized" MPD for ν = 1.

Table 1 .

 1 Fulfillment of condition(2.16

		2, 4, 5, 6	3
	n = 7	0, 1, 2, 5, 6, 7	3, 4
	n = 8	0, 1, 2, 6, 7, 8	3, 4, 5
	n = 9	0, 1, 2, 7, 8, 9	3, 4, 5, 6
	n = 10	0, 1, 2, 3, 7, 8, 9, 10	4, 5, 6

Table 2 .

 2 Details of the results obtained in the search for the valence electron pairs of the H-F molecule in Section 5.1.

	P 2 (Ω MPD ) Volume Population Residual	Cpu time Status Iter.
	Case 1 0.40445	7.302	1.90208	4.43 × 10 -4 3h12mn	CV	76
	Case 3 0.42908	65.51	1.94882	1.01 × 10 -4 31mn11s CV	14
	Case 2 0.00000	277.4	9.99520	4.40 × 10 -54 2mn49s	Failed 5

Acknowledgements. The work of C.D. and Y.P. is partly supported by the project ANR-18-CE40-0013 SHAPO, financed by the French Agence Nationale de la Recherche (ANR).

Algorithm 1 Level set mesh evolution algorithm for the resolution of (P s.o. ν ). Initialization:

• Computational domain D (see Section 4.6 below about its practical determination).

• Analytical expression of the wave function Ψ of the considered electron system; see Section 2.2.

• Initial shape Ω 0 , and mesh T 0 of D where Ω 0 explicitly appears as a submesh.

for k = 0, 1, 2, . . . until convergence do (1) Calculation of an ascent direction. Calculate F Ω k ν (x) and P ν (Ω k )(θ) as in Sections 4.1 and 4.3. A smooth ascent direction θ k : D → R 3 is inferred by using the extension-regularization procedure detailed in Remark 4.2.

(2) Update of the shape.

(a) Generate one level set function ϕ k for Ω k on the actual mesh T k of D;

(b) Solve the level set evolution equation (4.16) with velocity field θ k and initial data ϕ k over a suitably short time period (0, τ k ). A level set function for the new shape Ω k+1 is obtained at the vertices of T k . (c) Create a new mesh T k+1 of D where Ω k+1 explicitly appears as a submesh; the mesh T k+1 is adapted to the geometry of Ω k+1 and to the variations of the spin orbitals φ i . (3) Stop the procedure when max

< ε stop (in practice ε stop = 0.005). end for return Optimized shape Ω k .

Addition of a fixed point strategy

In spite of its efficiency, the shape gradient Algorithm 1 is not devoid of flaws when it comes to the resolution of (P s.o.

ν ); more precisely: • (Local) maximizers of (P s.o. ν ) are often unbounded domains, which are not easily captured, especially when the chosen initial shape is bounded. Indeed, the spin orbitals φ i in (2.8) show exponential decay away from the nuclei of the considered molecule, and so do the function S ν in (2.5) and the derivative F Ω ν in (2.13). As a result, significant updates of the shape in regions of R 3 which are "far" from the nuclei are unlikely.

• In principle, shape optimization algorithms based on the method of Hadamard do not leave the room for changes in the topology of shapes: all the variations Ω θ in (4.1) of a given shape Ω are indeed homeomorphic to Ω. In practice, up to some abuse of the theoretical framework, certain topological changes may occur in the course of the optimization process; separate regions of the boundary of the shape (e.g. distant holes) may merge, but under no circumstances could holes appear inside the bulk of the shape, or "islands" could emerge within void.

To alleviate these issues, insofar as possible, we propose to couple the shape gradient Algorithm 1 with a fixed-point iteration strategy, inspired by that in [START_REF] Amstutz | Analysis of a level set method for topology optimization[END_REF][START_REF] Amstutz | A new algorithm for topology optimization using a level-set method[END_REF]; the latter appraises the sensitivity of the optimized functionals P ν (Ω) with respect to a different type of perturbations, of a more topological nature.

Let us first present this technique as a standalone resolution method for (P s.o. ν ). From Proposition 3.7 and Remark 3.8, the optimality conditions for (P s.o.

ν ) read as follows: if Ω is a global maximizer of P ν (Ω), then: (4.17) for a.e. x ∈ Ω, F Ω ν (x) ≥ 0, and for a.e.

In other words, if Ω is a global maximizer of P ν (Ω), then ϕ := -F Ω ν is one level set function for Ω, in the sense that (4.15) holds. be confirmed by the use of a more realistic chemical description of the considered electronic system, notably involving multi-determinant wavefunctions; see Remark 2.2 about this point.

Conclusion and open problems

Elaborating on the landmark model of [START_REF] Savin | Probability distributions and valence shells in atoms[END_REF], we have analyzed the shape optimization problem (P s.o.

ν ) characterizing Maximum Probability Domains, which are the delocalized, quantum counterparts of the Lewis structures from classical chemistry.

From the theoretical standpoint, a new characterization of MPDs was obtained, involving a relaxed version (P relax ν ) of the problem (P s.o.

ν ), which allows to prove existence of a global maximizer to (P s.o. ν ) in a wide variety of concrete situations. From the practical point of view, we have introduced a shape gradient algorithm for (P s.o. ν ), enriched with a fixed point iteration strategy suggested by our theoretical findings, which offers significant improvements with respect to prior resolution techniques dedicated to (P s.o. ν ). The numerical experiments conducted owing to these techniques have interesting interpretations, paving the way to promising further applications in computational chemistry.

The present work opens the way to multiple directions for future work. The most crucial challenge towards unleashing the full potential of the MPD paradigm for chemical interpretations is certainly the efficient treatment of multi-determinant wave functions for the description of the quantum arrangement of electrons within the considered molecules. This task is highly non trivial, since the intricate structure of these general wave functions does not lend itself to the convenient algorithmic simplifications brought to light in [START_REF] Cances | How electrons guard the space: shape optimization with probability distribution criteria[END_REF], which are key in making numerical computations affordable.

From the algorithmic viewpoint, it would be desirable to improve the computational efficiency of our shape optimization algorithm, which remains computationally expensive when large chemical structures are concerned. For instance, the (costly) remeshing stages involved in Algorithm 1 could be parallelized thanks to the new capabilities of the library mmg.

Eventually, it would be interesting to investigate the second-order shape sensitivity of the probability functionals P ν (Ω) featured in (P s.o. ν ). On the one hand, such information would make it possible to appraise the stability of the calculated MPD; on the other hand, it could be exploited to devise more efficient secondorder algorithms.