E Bretin

S Masnou

E Oudet

Phase-field approximations of the Willmore functional and flow

We discuss in this paper phase-field approximations of the Willmore functional and the associated L 2 -flow. After recollecting known results on the approximation of the Willmore energy and its L 1 relaxation, we derive the expression of the flows associated with various approximations, and we show their behavior by formal arguments based on matched asymptotic expansions. We introduce an accurate numerical scheme, whose local convergence can be proved, to describe with more details the behavior of two flows, the classical and the flow associated with an approximation model due to Mugnai. We propose a series of numerical simulations in 2D and 3D to illustrate their behavior in both smooth and singular situations.

Introduction

Phase-field approximations of the Willmore functional have raised quite a lot of interest in recent years, both from the theoretical and the numerical viewpoints. In particular, attention has been given to understanding the continuous and numerical approximations of both smooth and singular sets with finite relaxed Willmore energy. Various approximation models have been proposed so far, whose properties are known only partially. Our main motivation in this paper is a better understanding of these models, and more precisely:

1. Exhibiting algebraic differences/similarities between the various approximations; 2. Deriving the L 2 -flows associated with these models; 3. Studying the asymptotic behavior of the flows, at least in smooth situations; 4. Simulating numerically these flows, and observing whether and how singularities may appear.

We focus on four models due, respectively, to De Giorgi, Bellettini, and Paolini [START_REF] Giorgi | Some remarks on Γ-convergence and least square methods[END_REF][START_REF] Bellettini | Approssimazione variazionale di funzionali con curvatura[END_REF], Bellettini [START_REF] Bellettini | Variational approximation of functionals with curvatures and related properties[END_REF], Mugnai [START_REF] Mugnai | Gamma-convergence results for phase-field approximations of the 2D-Euler elastica functional[END_REF], and Esedoglu, Rätz, and Röger [START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF]. The paper is organized as follows: Section 2 is an introductory section where we collect known results on the diffuse approximation of the perimeter, the diffuse approximation of the Willmore energy, and the critical issue of approximating singular sets with finite relaxed Willmore energy. We also recall the definitions of the above mentioned approximations. In Section 2.5, we make new observations on the differences between these different diffuse energies. Section 3 is devoted to the derivation of the L 2 -flows associated with, respectively, Bellettini's, Mugnai's, and Esedoglu-Rätz-Röger's models (actually a variant of the latter), and, for every flow, we use the formal method of matched asymptotic expansions to derive the asymptotic velocity of the limit 1 Université de Lyon, CNRS UMR 5208, INSA de Lyon, Institut Camille Jordan, 20, avenue Albert Einstein, F-69621 Villeurbanne Cedex, France. Email: bretin@cmap.polytechnique.fr 2 Université de Lyon, CNRS UMR 5208, Université Lyon 1, Institut Camille Jordan, 43 boulevard du 11 novembre 1918, F-69622 Villeurbanne-Cedex, France. Email: masnou@math.univ-lyon1.fr 3 Lab. Jean Kuntzmann, Université Joseph Fourier, Tour IRMA, BP 53, 51, rue des Mathématiques, F-38041 Grenoble Cedex 9, France. Email: edouard.oudet@imag.fr interface as the diffuse approximation becomes asymptotically sharp. We show in particular that, in dimensions 2 and 3 for all flows, and in any dimension for some of them, they correspond asymptotically to the continuous Willmore flow as long as the interface is smooth. In Section 4, we focus on the numerical simulation of De Giorgi-Bellettini-Paolini's flow (which we shall refer to as the classical flow) and Mugnai's flow, and we propose a fixed-point algorithm whose local convergence can be proved. We illustrate with various numerical examples the behavior of both flows in space dimensions 2 and 3, both in smooth and singular situations. We show in particular that our scheme can capture with good accuracy well-known singular configurations yielded by the classical flow, and that these configurations evolve as if the parametric Willmore flow were used. We also illustrate with several simulations that, in contrast, Mugnai's flow prevents the creation of singularities.

2 What is known?

Genesis : the van der Waals-Cahn-Hilliard interface model and the diffuse approximation of perimeter

In his 1893 paper on the thermodynamic theory of capillary (see an English translation, with interesting comments, in [START_REF] Translation | The thermodynamik of capillarity under the hypothesis of a continuous variation of density[END_REF]), van der Waals studied the free energy of a liquid-gas interface. Arguing that the density of molecules at the interface can be modelled as a continuous function of space u, he used thermodynamic and variational arguments to derive an expression of the free energy, in a small volume V enclosing the interface, as V ( f 0 (u) + λ|∇u| 2 )dx, where f 0 (u) denotes the energy of a homogeneous phase at density u and λ is the capillarity coefficient. The same expression was derived by Cahn and Hilliard in 1958 in their paper [START_REF] Cahn | Free energy of a nonuniform system. I. Interfacial free energy[END_REF] on the interface energy, to a first approximation, of a binary alloy with u denoting the mole fraction of one component. Cahn and Hilliard argued that both terms in the energy have opposite contributions: if the transition layer's size increases, then the gradient term diminishes, but this is possible only by introducing more material of nonequilibrium composition, and thus at the expense of increasing V f 0 (u)dx. Rescaling the energy, and changing the notations in the obvious way, yields the general form

F ε (u) = V ( ε 2 |∇u| 2 + W(u) ε )dx (1)
In the original papers of van der Waals, Cahn and Hilliard, f 0 was a smooth double-well function, yet with a slope between the local minima. For simplicity, since it does not modify the mathematical analysis, W will denote in the sequel a smooth double-well function with no slope (we will take in general W(s) = 1 2 s 2 (1s) 2 ). Two equations are usually associated with the van der Waals-Cahn-Hilliard energy, and will be used in this paper: the Allen-Cahn and the Cahn-Hilliard equations. The evolution Allen-Cahn equation is the L 2 -gradient descent associated with [START_REF] Ambrosio | Geometric evolution problems, distance function and viscosity solutions[END_REF] and is written

u t = W (u) ε -ε∆u.
We shall also refer to the stationary Allen-Cahn equation

W (u) -ε 2 ∆u = 0.
The Cahn-Hillard evolution equation is derived in a different manner: from a mathematical viewpoint, it is the H -1 -gradient flow associated with the van der Waals-Cahn-Hilliard energy [START_REF] Cahn | Overview no. 113 -surface motion by surface diffusion[END_REF][START_REF] Fife | Models for phase separation and their mathematics[END_REF]. The physical derivation of the equation is also instructive [START_REF] Cahn | On spinodal decomposition[END_REF][START_REF] Cahn | Free energy of a nonuniform system. I. Interfacial free energy[END_REF]: since ∇ u F ε (u) = -ε∆u + W (u) ε quantifies how the energy changes when molecules change position, it coincides with the chemical potential µ. Fick's first law states that the flux of particles is proportional to the gradient of µ, i.e. J = -α∇µ. Finally, the conservation law u t + div J = 0 yields the Cahn-Hilliard evolution equation

u t = α∆(-ε∆u + W (u) ε ).
To summarize, the Allen-Cahn equation describes the motion of phase boundaries driven by surface tension, whereas the Cahn-Hilliard equation is a conservation law that characterizes the motion induced by the chemical potential, which is the gradient of the surface tension.

Let us now recall the asymptotic behavior of F ε (u), as ε → 0 + , that has been exhibited by Modica and Mortola [START_REF] Modica | Il limite nella Γ-convergenza di una famiglia di funzionali ellittici[END_REF] following a conjecture of De Giorgi. We first fix some notations. Let Ω ⊂ R n be open, bounded and with Lipschitz boundary. W ∈ C 3 (R, R + ) is a double-well potential with two equal minima (in the sequel we will work, unless specified, with W(s) = 1 2 s 2 (1s) 2 ). Modica and Mortola have shown that the Γ-limit in L 1 (Ω) of the family of functionals

P ε (u) =    Ω ε 2 |∇u| 2 + W(u) ε dx if u ∈ W 1,2 (Ω) +∞ otherwise in L 1 (Ω)
is c 0 P(u) where

P(u) = |Du|(Ω) if u ∈ BV(Ω, {0, 1}) +∞ otherwise in L 1 (Ω)
and c 0 = 1 0 2W(s)ds. In particular, if E ⊂ R n has finite perimeter in Ω and u := 1 E ∈ BV(Ω, {0, 1}), then one can build a sequence of functions (u ε ) ∈ W 1,2 (Ω) such that u ε → u ∈ L 1 (Ω) and

P ε (u ε ) → c 0 |Du|(Ω) = c 0 P(E, Ω) with P(E, Ω) the perimeter of E in Ω.
To prove it, it is enough by density to restrict to smooth sets. Being E smooth, a good approximating sequence is given by u ε = q( d(x) ε ) (actually a variant of this expression, but we shall skip the details for the moment) where d is the signed distance function at ∂E, i.e. d(x) = -d(x, ∂E) if x ∈ E and d(x, ∂E) else, and

q(t) = 1-tanh(t) 2 is the unique decreasing minimizer of R ( |ϕ (t)| 2 2 + W(ϕ(t)))dt (2)
under the assumptions lim t→-∞ ϕ(t) = 1, lim t→∞ ϕ(t) = 0 and ϕ(0) = 1 2 . In other words, the approximation of u = 1 E is done by a suitable rescaling of the level lines of the distance function to ∂E. Such rescaling is optimal, in the sense that it minimizes the transversal energy (2) and forces the concentration as ε → 0 + .

Observe now that

Ω ε 2 |∇u| 2 + W(u) ε dx ≥ Ω √ 2|∇u| W(u)
and the equality holds if ε 2 |∇u| 2 = W(u) ε . Therefore, by lower semicontinuity arguments, the quality of an approximation depends on the so-called discrepancy measure

ξ ε = ( ε 2 |∇u| 2 - W(u) ε )L 2 ,
that will play an even more essential role for some diffuse approximations of the Willmore functional.

De Giorgi-Bellettini-Paolini's approximation of the Willmore energy

Based on a conjecture of De Giorgi [START_REF] Giorgi | Some remarks on Γ-convergence and least square methods[END_REF], several authors [START_REF] Bellettini | Approssimazione variazionale di funzionali con curvatura[END_REF][START_REF] Tonegawa | Phase field model with a variable chemical potential[END_REF][START_REF] Bellettini | On the approximation of the elastica functional in radial symmetry[END_REF][START_REF] Moser | A higher order asymptotic problem related to phase transitions[END_REF][START_REF] Röger | On a modified conjecture of De Giorgi[END_REF][START_REF] Nagase | A singular perturbation problem with integral curvature bound[END_REF] have investigated the diffuse approximation of the Willmore functional, which is for a set E ⊂ R N with smooth boundary in Ω:

W(E, Ω) = 1 2 ∂E∩Ω |H ∂E (x)| 2 dH N-1
where H ∂E (x) is the classical mean curvature vector at x ∈ ∂E. The approximation functionals are defined as

W ε (u) =    1 2ε Ω ε∆u - W (u) ε 2 dx if u ∈ L 1 (Ω) ∩ W 2,2 (Ω) +∞ otherwise in L 1 (Ω)
Introduced by Bellettini and Paolini in [START_REF] Bellettini | Approssimazione variazionale di funzionali con curvatura[END_REF], they differ from the original De Giorgi's conjecture in the sense that the perimeter is not explicitly encoded in the expression. They have however the advantage to be directly related to the Cahn-Hilliard equation, whose good properties [START_REF] Chen | Global asymptotic limit of solutions of the Cahn-Hilliard equation[END_REF] play a key role in the approximation. In the sequel, we shall refer to these functionals as the classical approximation model.

The reason why ε∆u -W (u) ε is related to the mean curvature can be simply understood at a formal level: it suffices to observe that the mean curvature of a smooth surface is associated with the first variation of its area, and that -ε∆u

+ W (u) ε is the L 2 gradient of ε 2 |∇u| 2 + W(u)
ε that appears in the approximation of the surface area.

The results on the asymptotic behavior of W ε as ε → 0 + have started with the proof by Bellettini and Paolini [START_REF] Bellettini | Approssimazione variazionale di funzionali con curvatura[END_REF] of a Γlim sup property, i.e. the Willmore energy of a smooth hypersurface E is the limit of W ε (u ε ), up to a multiplicative constant, where u ε is defined exactly as for the approximation of the perimeter.

The Γlim inf property is much harder to prove. The contributions on this point [START_REF] Tonegawa | Phase field model with a variable chemical potential[END_REF][START_REF] Bellettini | On the approximation of the elastica functional in radial symmetry[END_REF][START_REF] Moser | A higher order asymptotic problem related to phase transitions[END_REF][START_REF] Röger | On a modified conjecture of De Giorgi[END_REF], 66] culminated with the proof by Röger and Schätzle [START_REF] Röger | On a modified conjecture of De Giorgi[END_REF] in space dimensions N = 2, 3 and, independently, by Nagase and Tonegawa [START_REF] Nagase | A singular perturbation problem with integral curvature bound[END_REF] in dimension N = 2, that the result holds true for smooth sets. More precisely, given u = 1 E with E ∈ C 2 (Ω), and u ε converging to u in L 1 with a uniform control of

P ε (u ε ), then c 0 W(E, Ω) ≤ lim inf ε→0 + W ε (u ε ).
The proof is based on a careful control of the discrepancy measure

ξ ε = ( ε 2 |∇u| 2 -W(u) ε )L 2 that guar- antees good concentration properties, i.e. the varifolds v ε = |∇u ε |L 2 ⊗ δ ∇u ⊥ ε
(whose mass is related naturally to the variations of the approximating functions u ε ) concentrate to a limit integer varifold that has generalized mean curvature in L 2 and that is supported on a supset of ∂E. Then, in Röger and Schätzle's proof, a lower semicontinuity argument and the locality of integer varifolds' mean curvature yields the result. It holds in dimensions N = 2, 3 at most due to dimensional requirements for Sobolev embeddings and for the control of singular terms used in the proof. The result in higher dimension is still open.

What about unsmooth sets? Can the approximation results be extended to the relaxed Willmore functional? The answer is negative in general, as discussed below.

The approximation does not hold in general for unsmooth limit sets

Define for any set E of finite perimeter in Ω its relaxed Willmore functional

W(E, Ω) = inf{lim inf W(E h , Ω), ∂E h ∈ Ω ∈ C 2 , E h → E in L 1 (Ω)}.
The properties of this relaxation are fully known in dimension 2 [START_REF] Bellettini | Semicontinuity and relaxation properties of curvature depending functional in 2D[END_REF][START_REF] Bellettini | Characterization and representation of the lower semicontinuous envelope of the elastica functional[END_REF][START_REF] Bellettini | A varifold representation of the relaxed elastica functional[END_REF][START_REF] Masnou | A coarea-type formula for the relaxation of a generalized Willmore functional[END_REF]] and partially known in higher dimension [START_REF] Ambrosio | A direct variational approach to a problem arising in image reconstruction[END_REF][START_REF] Leonardi | Locality of the mean of rectifiable varifolds[END_REF][START_REF] Masnou | Gradient Young measures, varifolds, and a generalized Willmore functional[END_REF]. It is natural to ask whether the Γ-convergence of W ε to W can be extended to W. Unfortunately, this is not the case as it follows from the following observations (for simplicity we denote Γlim W ε (E) = Γlim W ε (1 E )) that are illustrated in Figure 1: 1. there exists a bounded set E 1 ⊂ R 2 of finite perimeter such that

Γ -lim W ε (E 1 ) < ∞ and W(E 1 ) = +∞ 2. there exists a bounded set E 2 ⊂ R 2 of finite perimeter such that Γ -lim W ε (E 2 ) < W(E 2 ) < +∞
The reason why W(E 1 ) = +∞ is a result by Bellettini, Dal Maso and Paolini [START_REF] Bellettini | Semicontinuity and relaxation properties of curvature depending functional in 2D[END_REF] according to which a non oriented tangent must exist everywhere on the boundary. Besides, W(E 2 ) < +∞ because, still by a result of Bellettini, Dal Maso and Paolini, the boundary is smooth out of evenly many cusps. Let us now explain why, in both cases, Γlim W ε (E 1,2 ) < +∞. The reason for this is the existence of smooth solutions with singular nodal sets for the Allen-Cahn equation

∆u -W (u) = 0 Figure 1: Top: a set E 1 such that Γ -lim W ε (E 1 ) < ∞ and W(E 1 ) = +∞.
Bottom, from left to right, a set E 2 , the limit configuration whose energy coincides with Γlim W ε (E 2 ), a configuration whose energy coincides with W(E 2 ).

According to Dang, Fife and Peletier [START_REF] Dang | Saddle solutions of the bistable diffusion equation[END_REF], there exists for such equation in R 2 a unique saddle solution u with values in (-1, 1). By saddle solution, it is meant that u(x, y) > 0 in quadrants I and III, and u(x, y) < 0 in quadrants II and IV, in particular u(x, y) = 0 on the nodal set xy = 0. Considering u ε (x) = u(εx), we immediately get that

ε 2 ∆u ε -W (u ε ) = 0
thus the second term in W ε (u ε ) vanishes, and the first term being obviously bounded, it follows from the lower semicontinuity of the Γ-limit that Γlim W ε (E 1,2 ) < +∞. Furthermore, the approximation of E 2 can be made so as to create a cross in the limit, as in bottom-middle figure. The limit energy is therefore lower than the energy obtained by pairwise connection without crossing of the cusps (bottom-right figure). Thus, Γlim W ε (E 2 ) < W(E 2 ) < +∞.

For the reader not familiar with varifolds, it must be emphasized that this is not in contradiction with the results described in the previous section, and more precisely with the fact that the discrepancy measure guarantees the concentration of the diffuse varifolds at a limit integer varifold with generalized curvature in L 2 . Indeed, the boundary curves of E 1 and of the bottom-middle set can be canonically associated with a varifold having L 2 generalized curvature because, by compensation between the tangents associated with each branch meeting at the cross, there is no singularity.

We end this section with the question that follows naturally from the discussion above: is it possible to find a diffuse approximation that Γ-converges to W (up to a multiplicative constant) whenever W(E) < +∞ ?

Diffuse approximations of the relaxed

W Be ε (u) =    1 2 Ω\{|∇u|=0} (| div ∇u |∇u| | 2 )( ε 2 |∇u| 2 + W(u) ε )dx if u ∈ C ∞ (Ω) +∞ otherwise in L 1 (Ω)
P ε + W Be ε )(E) = c 0 (P(E) + W(E))
for every E of finite perimeter such that W(E) < +∞

The constructive part of the proof is based, as usual, on using approximating functions of the form u ε = q( d E ε ). As for the lower semicontinuity part, it is facilitated by the explicit appearance of the mean curvature in the expression. Recall indeed that, u being smooth, for almost every t, H u (x) := (div ∇u |∇u| ) ∇u |∇u| (x) is the mean curvature at a point x of the isolevel {y, u(y) = t}. Let (u ε ) be a sequence of smooth functions that approximate u = 1 E in L 1 (R N ) and has uniformly bounded total variation. Then, by the coarea formula,

W Be ε (u ε ) ≥ 1 2 {|∇u h | =0} |∇u ε | 2W(u ε )|H u ε | 2 dx = 1 2 1 0 2W(t) {u ε =t}∩{|∇u ε | =0} |H u ε | 2 dH N-1 dt.
The last inequality is important: it guarantees a control of the Willmore energy of the isolevel surfaces of u ε . This is a major difference with the classical approximation, for which such control does not hold. Then, it suffices to observe that, by the Cavalieri formula and for a suitable subsequence, |{u ε ≥ t}∆{u ≥ t}| → 0 for almost every t. In addition, {u ≥ t} = E for almost every t, and by the lower semicontinuity of the relaxation F(E) ≤ lim inf ε→0 F(u ε ). Fatou's Lemma finally gives

lim inf ε→0 W Be ε (u ε ) ≥ F(E) 1 0 2W(t)dt = c 0 F(E).
Bellettini's approximation has however a drawback that will be explained with more details later: when one computes the flow associated with the functional, the 4th order term is nonlinear, which raises difficulties at the numerical level since it cannot be treated implicitly.

Mugnai's approximation in dimension N = 2

In the regular case and in dimensions 2,3, it follows from the results of Bellettini and Mugnai [START_REF] Bellettini | Approximation of Helfrich's via diffuse interfaces[END_REF] that, up to a uniform control of the perimeter, the Γ-limit of the functionals defined by

W Mu ε (u) =    1 2ε Ω ε∇ 2 u - W (u) ε ν u ⊗ ν u 2 dx if u ∈ C 2 (Ω) +∞ otherwise in L 1 (Ω),
where ν u = ∇u |∇u| when |∇u| = 0, and ν u =constant unit vector on {|∇u| = 0}, coincides with

c 0 Ω∩∂E |A ∂E (x)| 2 dx
for every smooth E, with A ∂E (x) the second fundamental form of ∂E at x. Again, this approximation allows a control of the mean curvature of the isolevel surfaces of an approximating sequence u ε , thus prevents from the creation of saddle solutions to the Allen-Cahn equation since, by [START_REF] Bellettini | Approximation of Helfrich's via diffuse interfaces[END_REF]Lemma 5.3] and [65, Lemma 5.2],

|∇u|| div ∇u |∇u| | ≤ 1 ε |ε∇ 2 u - W (u) ε ν u ⊗ ν u |
In dimension 2, the second fundamental form along a curve coincides with the curvature. Therefore, by identifying the limit varifold obtained when u ε converges to u = 1 E , and using the representation results of [START_REF] Bellettini | A varifold representation of the relaxed elastica functional[END_REF], Mugnai was able to prove in [START_REF] Mugnai | Gamma-convergence results for phase-field approximations of the 2D-Euler elastica functional[END_REF] that, in dimension 2, the Γ-limit of W Mu ε (with uniform control of the perimeter) coincides with F(E) for any E with finite perimeter.

Esedoglu-Rätz-Röger's approximation in dimension N ≥ 2

The model of Esedoglu, Rätz, and Röger in [START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF] is a modification of the classical energy that aims to preserve the "parallelity" of the level lines of the approximating functions, and avoids the formation of saddle points, by constraining the level lines' mean curvature using a term à la Bellettini. More precisely, one can calculate that ε∆u -W (u)

ε = ε|∇u| div ∇u |∇u| -∇ξ ε • ∇u |∇u| 2 . with ξ ε = ε 2 |∇u| 2 -W(u)
ε the discrepancy function (with a small abuse of notation, we use the same notation for the discrepancy measure and its density).

Therefore, ε∆u -W (u) ε approximates correctly the mean curvature (up to a multiplicative constant) if the projection of ∇ξ ε on the orthogonal direction ∇u is small. Equivalently, it can be required that

ε∆u - W (u) ε -ε|∇u| div ∇u |∇u|
be small, therefore a natural profile-forcing approximation model is (with α ≥ 0 a parameter):

W EsRäRö ε (u) =            1 2ε Ω ε∆u - W (u) ε 2 dx+ 1 2ε 1+α Ω (ε∆u - W (u) ε -ε|∇u| div ∇u |∇u| ) 2 dx if u ∈ C ∞ (Ω) +∞ otherwise in L 1 (Ω)
To simplify the theoretical analysis, the model proposed by Esedoglu, Rätz, and Röger is slightly different. It uses the fact that, if a phase field u ε resembles q( d ε ), one has ε|∇u| ∼ 2W(u), which leads Esedolu, Rätz, and Röger to penalize

ε∆u - W (u) ε -(ε|∇u|(2W(u)) 1 2
)

1 2 div ∇u |∇u| .
Finally, they propose the following approximating functional

W EsRäRö ε (u) =            1 2ε Ω ε∆u - W (u) ε 2 dx+ 1 2ε 1+α Ω (ε∆u - W (u) ε -(ε|∇u| 2W(u)) 1 2 div ∇u |∇u| ) 2 dx if u ∈ C ∞ (Ω) +∞ otherwise in L 1 (Ω)
This energy controls the mean curvature of the level lines of an approximating function since (see [START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF])

W EsRäRö ε (u) ≥ ε -α 2 + 2ε -α 1 0 2W(t) {u=t}∩{∇u =0} (div ∇u |∇u| ) 2 dH N-1 dt,
which, once again, excludes Allen-Cahn solutions. With the control above, the authors prove with the same argument as Bellettini [START_REF] Bellettini | Variational approximation of functionals with curvatures and related properties[END_REF] that, for any α > 0,

Γ -lim ε→0 P ε + W EsRäRö ε = c 0 P + W in L 1 (Ω).
With α = 0 the Γ-convergence result does not hold anymore, but instead, with a uniform control of the perimeter,

Γ -lim ε→0 W EsRäRö ε ≥ c 0 2 W.
which still guarantees a control of W.

For the sake of numerical simplicity, another version is tackled numerically in [START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF], based again on the approximation ε|∇u| ∼ 2W(u):

W EsRäRö ε (u) =            1 2ε Ω ε∆u - W (u) ε 2 dx+ 1 2ε 1+α Ω (ε∆u - W (u) ε -2W(u) div ∇u |∇u| ) 2 dx if u ∈ C ∞ (Ω) +∞ otherwise in L 1 (Ω)
We will focus in the sequel on W EsRäRö ε , whose flow will be derived, as well as its asymptotic behavior as ε goes to 0.

Few remarks on the connections between the different approximations 2.5.1 From Mugnai's model to Esedoglu-Rätz-Röger's

We saw previously that the phase-field approximations W Be ε and W EsRäRö ε Γ-converge, up to a uniform control of perimeter, to c 0 W in any dimension, and the same holds true in dimension 2 for W Mu ε . We will now emphasize the connections between these approximations. More precisely, we will see that Mugnai's approximation W Mu ε can be viewed as the sum of a geometric-type approximation of the Willmore energy plus a profile penalization term of the same kind as in Esedoglu, Rätz, Röger's model (or, more precisely, the initial model W EsRäRö ε ). Indeed we have, denoting ν = ∇u |∇u| when |∇u| = 0, and ν =constant unit vector on {|∇u| = 0},

W Mu ε (u) = 1 2ε Ω ε∇ 2 u - W (u) ε ν ⊗ ν 2 dx = 1 2ε Ω\{|∇u|=0} ε∇ 2 u : ν ⊗ ν - W (u) ε 2 dx + Ω ε 2 |∇ 2 u| 2 -(∇ 2 u : ν ⊗ ν) 2 dx.
where, being A, B two matrices, we denote as A : B = ∑ i,j A ij B ij the usual matrix scalar product. Using A : e 1 ⊗ e 2 =< Ae 2 , e 1 >, we observe that ∇ 2 u : ν ⊗ ν = ∆u -|∇u| div ∇u |∇u| , therefore the first term of W Mu ε coincides with the second term of W EsRäRö ε for α = 0. The second term of W Mu ε can be splitted as

Ω ε 2 |∇ 2 u| 2 -(∇ 2 u : ν ⊗ ν) 2 dx = 1 2 Ω ∇ ∇u |∇u| 2 ε|∇u| 2 dx, + Ω ε 2 |∇ 2 u ν| 2 -|∇ 2 u : ν ⊗ ν| 2 dx. Note that Ω ε |∇ 2 u ν| 2 -∇ 2 u : ν ⊗ ν 2 dx ≥ 0,
is positive and vanishes for all functions u of the general form u = η (d(x)) with η smooth. It is therefore a soft profile-penalization term that forces the approximating function to be a profile, yet not necessarily the optimal profile q. As for the term

Ω ∇ ∇u |∇u| 2 ε|∇u| 2 dx, it is purely geometric
and constrains the approximating function's level lines mean curvature. It would therefore be worth addressing the Γ-convergence of the new functional

W New ε (u) = 1 2ε Ω (ε∆u - 1 ε W (u)) 2 dx + 1 2ε α ε |∇ 2 u ν| 2 -∇u 2 : ν ⊗ ν 2 dx.
The reason why such approximation would be interesting is that, if it indeeds Γ-converges, the associated flow would not be influenced by the asymptotic behavior of the penalization term, since it vanishes for approximating functions that are profiles. More precisely, the Willmore flow could be captured at low order of ε, and not at the numerically challenging order ε 3 as for the Esedoglu-Rätz-Röger model with α = 0 or 1.

Towards a modification of Mugnai's energy that forces the Γ-convergence in dimension ≥ 3

Obviously, we cannot expect that Mugnai's energy Γ-converges to the Willmore energy in dimension greater than 2, since for E smooth

|A ∂E | 2 = |H ∂E | 2 -∑ i =j κ i κ j
where κ 1 , κ 2 . . . κ N-1 are the principal curvatures. This identity suggests however that a suitable correction could force the Γ-convergence, i.e. by subtracting to W Mu ε an approximation of

J(E, Ω) = ∂E∩Ω ∑ i =j κ i κ j dH N-1 .
Recalling our assumptions that d < 0 in E, and as an easy consequence of Lemma 14.17 in [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF] (see also [START_REF] Ambrosio | Geometric evolution problems, distance function and viscosity solutions[END_REF]), we obtain in a small tubular neighborhood of ∂E:

div (∆d(x)∇d(x)) = (∆d(x)) 2 + ∇∆d(x) • ∇d(x) = ∑ i κ i (π(x)) 1 + d(x)κ i (π(x)) 2 -∑ i κ i (π(x)) 2 (1 + d(x)κ i (π(x))) 2 ∑ i =j κ i κ j on ∂E,
where π(x) is the projection of x on Γ. Thus, a possible approximation of c 0 J(E, Ω) is

J 1 ε (u) = - 2 ε Ω ε∆u - 1 ε W (u) W (u) ε dx.
Indeed, with u = q(d/ε) and with a suitable truncation of q so that q (d/ε) vanishes on ∂Ω (which is always possible if E ⊂⊂ Ω), integrating by parts yields:

J 1 ε (u) = - 2 ε 2 Ω ∆dq d ε q d ε dx = - 1 ε Ω ∆d∇ q d ε 2 .∇d dx = 1 ε Ω div (∆d∇d) q d ε 2 dx c 0 ∂E∩Ω ∑ i =j κ i κ j dH N-1 .
Remark also that since a profile function

u = q(d/ε) satisfies 1 ε W (u) = ε∇ 2 u : N(u) where N(u) = ν ⊗ ν = ∇u |∇u| ⊗ ∇u |∇u| , the following energies        J 2 ε (u) = -2 ε 2 Ω ε∆u -ε∇ 2 u : N(u) W (u)dx J 3 ε (u) = -2 Ω ε∆u -1 ε W (u) ∇ 2 u : N(u)dx J 4 ε (u) = -2 Ω ε∆u -ε∇ 2 u : N(u) ∇ 2 u : N(u)dx
approximate also c 0 J(E, Ω). In particular, as a modified version of Mugnai's energy, we can consider

W Mu ε = W Mu ε + 1 2 (J 1 ε (u) + J 3 ε (u) -J 4 ε (u))
We have indeed

J 1 ε (u) + J 3 ε (u) -J 2 ε (u) -J 4 ε (u) = 2 ε Ω ε∇ 2 u : ν ⊗ ν - W (u) ε 2 dx. and W Mu ε = W ε (u) - 1 2 J 2 ε (u) since |∇ 2 u| 2 = (∆u) 2 , N(u) : N(u) = 1, |ε∇ 2 u -ε -1 W (u)N(u)| 2 = ε 2 |∇ 2 u| 2 -2W (u)∇ 2 u : N(u) + ε -2 W (u) 2 N(u) : N(u) and (ε∆u -ε -1 W (u)) 2 = ε 2 (∆u) 2 -2W (u)∆u + ε -2 (W (u)) 2 .
Therefore

W Mu ε = W ε (u) + 1 ε Ω ε∇ 2 u : ν ⊗ ν - W (u) ε 2 dx,
which resembles Esedoglu-Rätz-Röger's approximation with α = 0 since, in both approximations, the second term forces u to be a "profile" function, and vanishes at the limit. In view of the approximation result of Esedoglu, Rätz, and Röger, it is reasonable to expect that W Mu ε Γ-converges to the relaxed Willmore energy in any dimension.

The Willmore flow and its approximation by the evolution of a diffuse interface

This section is devoted to the approximation of the Willmore flow by L 2 -gradient flows associated with the approximating energies introduced above. In particular, we shall derive explicitly each approximating gradient flow and, using the matched asymptotic expansion method [START_REF] Caginalp | Dynamics of layered interfaces arising from phase boundaries[END_REF][START_REF] Pego | Front migration in the nonlinear Cahn-Hilliard equation[END_REF][START_REF] Bellettini | Quasi-optimal error estimates for the mean curvature flow with a forcing term[END_REF][START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF], we will show that, at least formally and for smooth interfaces, there is convergence to the Willmore flow, at least in dimensions 2 and 3 for all flows, and in any dimension for some of them. The general question "if a sequence of functionals Γ-converges to a limit functional, is there also convergence of the associated flows?" is rather natural, since Γ-convergence implies convergence of minimizers, up to the extraction of a subsequence. However, the question is difficult and remains open for the Willmore functional. Our results below give formal indications that the convergence holds. Serfaty discussed in [START_REF] Serfaty | Gamma-convergence of gradient flows on hilbert and metric spaces and applications[END_REF] a general theorem on the Γ-convergence of gradient flows, provided that the generalized gradient of the associated functional can be controled (see in particular the discussion on the Cahn-Hilliard flow). Such control is so far out or reach for the Willmore functional.

On the Willmore flow

Let E(t), 0 ≤ t ≤ T, denote the evolution by the Willmore flow of smooth domains, i.e. the outer normal velocity V(t) is given at x ∈ ∂E(t) by

V = ∆ S H - 1 2 H 3 + H A 2 ,
where ∆ S is the Laplace-Beltrami operator on ∂E(t), H the scalar mean curvature, A the second fundamental form, and A 2 is the sum of the squared coefficients of A.

In the plane, the Willmore flow coincides with the flow of curves associated with the Bernoulli-Euler elastica energy, i.e., denoting by κ the scalar curvature

V = ∆ S κ + 1 2 κ 3 .
The long time existence of a single curve evolving by this flow is established in [START_REF] Dziuk | Evolution of elastic curves in R n : Existence and computation[END_REF], and any curve with fixed length converges to an elastica. In higher dimension, Kuwert and Schätzle give in [START_REF] Kuwert | The Willmore flow with small initial energy[END_REF][START_REF] Kuwert | Removability of point singularities of Willmore surfaces[END_REF] a long time existence proof of the Willmore flow and the convergence to a round sphere for sufficiently small initial energy. Singularities may appear for larger initial energies, as indicated by numerical simulations [START_REF] Mayer | A numerical scheme for axisymmetric solutions of curvature driven free boundary problems, with applications to the Willmore flow[END_REF].

Approximating the Willmore flow with the classical De Giorgi-Bellettini-Paolini approach

The L 2 -gradient flow of the approximating energy

W ε (u) = 1 2ε Ω ε∆u - W (u) ε 2 dx,
is equivalent to the evolution equation

∂ t u = -∆ ∆u - 1 ε 2 W (u) + 1 ε 2 W (u) ∆u - 1 ε 2 W (u) ,
that can be rewritten as the phase field system

ε 2 ∂ t u = ∆µ -1 ε 2 W (u)µ µ = W (u) -ε 2 ∆u. (3)
Existence and well-posedness The well-posedness of the phase field model (3) at fixed parameter ε has been studied in [START_REF] Colli | A phase-field approximation of the Willmore flow with volume constraint[END_REF] with a volume constraint fixing the average of u, and in [START_REF] Colli | A phase-field approximation of the Willmore flow with volume and area constraints[END_REF] with both volume and area constraints.

Convergence to the Willmore flow Loreti and March showed in [START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF], by using the formal method of matched asymptotic expansions, that if ∂E is smooth and evolves by Willmore flow, it can be approximated by level lines of the solution u ε of the phase field system (3) as ε goes to 0. In addition, u ε and µ ε are expected to take the form

   u ε (x, t) = q d(x,E(t)) ε + ε 2 A 2 -1 2 H 2 η 1 d(x,E(t)) ε + O(ε 3 ) µ ε (x, t) = -εHq d(x,E(t)) ε + ε 2 H 2 η 2 d(x,E(t)) ε + O(ε 3 ) ,
where η 1 and η 2 are two functions depending only of the double well potential W, and defined as the solutions of η 1 (s) -W (q(s))η 1 (s) = sq (s), with lim s→±∞ η 1 (s) = 0, η 2 (s) -W (q(s))η 2 (s) = q (s), with lim s→±∞ η 2 (s) = 0.

An important point is that the second-order term in the asymptotic expansion of u ε has an influence on the limit law as ε goes to zero [START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF]. This is a major difference with the Allen-Cahn equation, for which the velocity law follows from the expansion at zero and first orders only [START_REF] Bellettini | Quasi-optimal error estimates for the mean curvature flow with a forcing term[END_REF]. As a consequence, addressing numerically the Willmore flow is more delicate and requires using a high accuracy approximation in space to guarantee a sufficiently good approximation of the expansion of u ε .

Approximating the Willmore flow with Bellettini's model

We focus now on the approximation model

W Be ε (u) = 1 2 Ω div ∇u |∇u| 2 ε 2 |∇u| 2 + 1 ε W(u) dx
We will prove in the next section that its L 2 -gradient flow is equivalent to the evolution equation

∂ t u = K(u) 2 2 ∆u - 1 ε 2 W (u) + 1 2 ∇[K(u) 2 ].∇u - 1 ε div P u ∇ [K(u)h ε (u)] |∇u| , (4)
where

P u = I d -∇u |∇u| ⊗ ∇u |∇u| = I d -N(u), h ε (u) = ε 2 |∇u| 2 + 1 ε W(u) and K(u) = div ∇u |∇u| .
Existence and well-posedness of this equation are open questions. Numerical simulations performed with this flow are shown in [START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF]. Note that the fourth-order nonlinear term makes numerics harder.

Using the formal method of matched asymptotic expansions, we show below that the phase field model (4) converges in any dimension, at least formally, to the Willmore flow as ε goes to 0. More precisely, we observe an asymptotic expansion of u ε of the form

u ε (x, t) = q d(x, E(t)) ε + O(ε 2 ),
where the second-order term does not have any influence on the limit velocity law as ε goes to zero, in contrast with the classical approximation of the previous section.

∂ t u = K(u) 2 2 ∆u - 1 ε 2 W (u) + 1 2 ∇[K(u) 2 ].∇u - 1 ε div P u ∇ [K(u)h ε (u)] |∇u| , where P u = I d -∇u |∇u| ⊗ ∇u |∇u| , h ε (u) = ε 2 |∇u| 2 + 1 ε W(u) and K(u) = div ∇u |∇u| .
Proof: The differential of K at u satisfies

K (u)(w) = lim t→0 K(u + tw) -K(u) t = div ∇w |∇u| - ∇u.∇w∇u |∇u| 3 , therefore (W Be ε (u)) (w) = Ω [K(u)h ε (u)] div ∇w |∇u| - ∇u.∇w∇u |∇u| 3 dx + 1 2 Ω K(u) 2 ε∇u∇w + 1 ε W (u)w dx
It follows that the L 2 -gradient of W Be ε reads as

∇W Be ε (u) = div ∇ [K(u)h ε (u)] |∇u| -div ∇ [K(u)h ε (u)] |∇u| • ∇u |∇u| ∇u |∇u| - 1 2 ε div K(u) 2 ∇u - 1 ε K(u) 2 W (u) .
whence the L 2 -gradient flow of W Be ε (u) follows.

Asymptotic analysis

In this section, we compute the formal expansions of the solution u ε (x, t) to the phase field model (4).

Preliminaries

We assume without loss of generality that the isolevel set

Γ(t) = {u ε = 1 2 } is a smooth n -1 dimensional boundary Γ(t) = ∂E(t) = ∂{x ∈ R d ; u ε (x, t) ≥ 1/2}.
We follow the method of matched asymptotic expansions proposed in [START_REF] Caginalp | Dynamics of layered interfaces arising from phase boundaries[END_REF][START_REF] Pego | Front migration in the nonlinear Cahn-Hilliard equation[END_REF][START_REF] Bellettini | Quasi-optimal error estimates for the mean curvature flow with a forcing term[END_REF][START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF]. We assume that the so-called outer expansion of u ε , i.e. the expansion far from the front Γ, is of the form

u ε (x, t) = u 0 (x, t) + εu 1 (x, t) + ε 2 u 2 (x, t) + O(ε 3 )
In a small neighborhood of Γ, we define the stretched normal distance to the front,

z = d(x, t) ε ,
were d(x, t) denotes the signed distance to E(t) such that d(x, t) < 0 in E(t). We then focus on inner expansions of u ε (x, t), i.e. expansions close to the front, of the form

u ε (x, t) = U(z, x, t) = U 0 (z, x, t) + εU 1 (z, x, t) + ε 2 U 2 (z, x, t) + O(ε 3 )
Let us define a unit normal m to Γ and the normal velocity V to the front as

V = -∂ t d(x, t), m = ∇d(x, t), x ∈ Γ,
where ∇ refers to spatial derivation only (the same holds for further derivation operators used in the sequel). Following [START_REF] Pego | Front migration in the nonlinear Cahn-Hilliard equation[END_REF][START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF] we assume that U(z, x, t) does not change when x varies normal to Γ with z held fixed, or equivalently ∇ x U.m = 0. This amounts to requiring that the blow-up with respect to the parameter ε is coherent with the flow.

Claim 3.2

In a suitable regime provided by the method of matched asymptotic expansions, the normal velocity of the 1 2 -front Γ(t) = ∂E(t) associated with a solution u ε (x, t) to Bellettini's phase field model (4) is the Willmore velocity

V = ∆ Γ H + A 2 H - H 3 2 ,
and

u ε (x, t) = q( d(x, E(t)) ε ) + O(ε 2 )
Following [START_REF] Pego | Front migration in the nonlinear Cahn-Hilliard equation[END_REF][START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF], it is easily seen that

     ∇u = ∇ x U + ε -1 m∂ z U ∆u = ∆ x U + ε -1 ∆d∂ z U + ε -2 ∂ 2 zz U ∂ t u = ∂ t U -ε -1 V∂ z U.
Recall also that in a sufficiently small neighborhood of Γ, according to Lemma 14.17 in [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF] (see also [START_REF] Ambrosio | Geometric evolution problems, distance function and viscosity solutions[END_REF]), we have

∆d(x, t) = n-1 ∑ i=1 κ i (π(x)) 1 + κ i (π(x))d(x, t) = n-1 ∑ i=1 κ i (π(x)) 1 + κ i (π(x))εz
where π(x) is the projection of x on Γ, and κ i are the principal curvatures on Γ. In particular this implies that ∆d

(x, t) = H -εz A 2 + O(ε 2 ),
where H and A 2 denote, respectively, the mean curvature and the squared 2-norm of the second fundamental form on Γ at π(x).

Outer solution:

We now compute the solution u ε in the outer region. By equation ( 4), u 0 satisfies W (u 0 ) = 0 and

u 0 (x, t) = 1 if x ∈ E(t) 0 otherwise
We also see that u 1 = 0 is a possible solution at the first order.

Matching condition :

The inner and outer expansions are related by the following matching condition

u 0 (x, t) + εu 1 (x, t) + • • • = U 0 (z, x, t) + εU 1 (z, x, t) + • • •
with x near the front Γ and εz between O(ε) and •(1). With the notation

u ± i (x, t) = lim s→0± u i (x + sm, t), one has that u ± 0 (x, t) = lim z→±∞ U 0 (z, x, t) lim z→±∞ u ± 1 (x, t) + zm • ∇u ± 0 (x, t) = lim z→±∞ U 1 (z, x, t)
In particular, for the phase field model ( 4), it follows that 2 , therefore, using the orthogonality condition ∇ x U.m = 0:

lim z→+∞ U 0 (z, x, t) = 0, lim z→-∞ U 0 (z, x, t) = 1 and lim z→±∞ U 1 (z, x, t) = 0, Inner solution: Note that ∇u |∇u| = m -ε∇ x U/∂ z U 1 + ε 2 |∇ x U| 2 /(∂ z U)
               K(u) = ∆d + O(ε) h ε (u) = 1 ε 1 2 (∂ z U) 2 + W(U) + O(1) 1 2 ∇ K(u) 2 .∇u = 1 ε (∆d∇(∆d).∇d) ∂ z U + O(1) 1
ε div P u ∇[K(u)h ε (u)] |∇u| = 1 ε div (∇(∆d) -∇(∆d) • ∇d∇d) 1 2 (∂ z U) 2 +W(U) |∂ z U| + O(1)
Recall also that in a sufficiently small neighborhood of Γ,

∆d(x, t) = n-1 ∑ i=1 κ i (π(x)) 1 + κ i (π(x))εz , thus (∆d∇(∆d).∇d) = -A 2 H + O(ε) div (∇(∆d) -∇(∆d) • ∇d∇d) = ∆ Γ H + O(ε).
Then, the first order in ε -2 of Equation ( 4) reads

H 2 2 ∂ 2 zz U 0 -W (U 0 ) = 0.
Adding the boundary condition obtained from the matching condition, and using U 0 (0, x, t) = 1/2 leads to U 0 = q(z).

Moreover, the second order in ε -1 of (4) shows that

-V∂ z U 0 = H 2 2 ∂ 2 zz U 1 -W (U 0 )U 1 + H 3 2 ∂ z U 0 -A 2 H∂ z U 0 -∆ Γ H 1 2 (∂ z U 0 ) 2 + W(U 0 ) |∂ z U 0 | .
As U 0 (z, x, t) = q(z) and q = -2W(q), we obtain

-Vq = H 2 2 ∂ 2 zz U 1 -W (q)U 1 + H 3 2 -A 2 H -∆ Γ H q .
Then, multiplying by q and integrating over R, it follows that

V = ∆ Γ H + A 2 H - H 3 2 ,
thus the sharp interface limit ∂E(t) as ε goes to zero evolves, at least formally, as the Willmore flow. In addition, we have U 1 = 0, therefore

u ε (x, t) = q( d(x, E(t)) ε ) + O(ε 2 )
and the second-order term does not appear in the expression of V. This explains the numerical stability, despite the use of an explicit Euler scheme, observed by Esedoglu, Rätz and Röger in [START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF].

Approximating the Willmore flow with Mugnai's model

The aim of this section is the derivation and the study of the L 2 -gradient flow associated with Mugnai's energy

W Mu ε (u) = 1 2ε Ω εD 2 u - 1 ε W (u) ∇u |∇u| ⊗ ∇u |∇u| 2 dx.
We will show that the flow is equivalent to the phase field system

ε 2 ∂ t u = ∆µ -1 ε 2 W (u)µ + W (u)B(u) µ = 1 ε 2 W (u) -∆u, (5)
where

B(u) = div div ∇u |∇u| ∇u |∇u| -div ∇ ∇u |∇u| ∇u |∇u| .
Note that this system coincides with the classical one, up to the addition of a penalty term L(u) = W (u)B(u).

The well-posedness of the phase field model (5) at fixed parameter ε is open, and requires presumably a regularization of the term B(u) as done numerically in the next section.

Claim 3.3

In a suitable regime provided by the method of matched asymptotic expansions, the normal velocity of the 1 2 -front Γ(t) = ∂E(t) associated with a solution (u ε , µ ε ) to Mugnai's phase field model (5) is

V = ∆ Γ H + ∑ i κ 3 i - 1 2 A 2 H. and    u ε (x, t) = q d(x,E(t)) ε + ε 2 A 2 2 η 1 d(x,E(t) ε + O(ε 3 ) µ ε (x, t) = -εHq d(x,E(t) ε + A 2 ε 2 η 2 d(x,E(t) ε O(ε 3 ),
where η 1 and η 2 are profile functions.

Remark 3.4

The front velocity associated with Mugnai's phase field model coincides, up to a multiplicative constant, with the velocity of the L 2 -flow of the squared second fundamental form energy Γ A 2 dH N-1 . Indeed, according to [2, Section 5.3], the latter is

Ṽ = 2∆ Γ H + 2H A 2 -H 3 + 6 ∑ i<j< κ i κ j κ Observing that H A 2 = ∑ i κ 3 i + ∑ i =j κ i κ 2 j and
H 3 = ∑ κ 3 i + 3 ∑ i =j κ i κ 2 j + 6 ∑ i<j< κ i κ j κ , one has H A 2 - 1 2
H 3 + 3 ∑ i<j< κ i κ j κ = 1 2 ( ∑ i κ 3 i -∑ i =j κ i κ 2 j ). Since ∑ i κ 3 i - 1 2 H A 2 = 1 2 ( ∑ κ 3 i -∑ i =j κ i κ 2 j )
we finally get that Ṽ = 2V.

Remark 3.5 It is easily seen that, in dimensions 2 and 3, Mugnai's flow coincides with the Willmore flow. It is obvious in dimension 2, whereas in dimension 3 one has

∑ κ 3 i - 1 2 H A 2 = κ 3 1 + κ 3 2 - 1 2 (k 1 + k 2 )(k 2 1 + k 2 2 ) = A 2 H - H 3 2 .
Another explanation involves Gauss-Bonnet Theorem. In Mugnai's model, the energy associated with the squared 2-norm of the second fundamental form prevents from topological changes. By Gauss-Bonnet Theorem, this energy coincides with the Willmore energy up to a topological additive constant, and thus both associated flows coincide.

Derivation of the gradient flow

Proposition 3. [START_REF] Barrett | Parametric approximation of Willmore flow and related geometric evolution equations[END_REF] The L 2 -gradient flow of Mugnai's model is equivalent to

ε 2 ∂ t u = ∆µ -1 ε 2 W (u)µ + W (u)B(u) µ = 1 ε 2 W (u) -∆u, where B(u) = div div ∇u |∇u| ∇u |∇u| -div ∇ ∇u |∇u| ∇u |∇u| . Proof: Let V(u) = εD 2 u - 1 ε W (u) ∇u |∇u| ⊗ ∇u |∇u| .
The differential of V in the direction w is

V (u)(w) = lim t→0 (V(u + tw) -V(u))/t = εD 2 w - 1 ε W (u)w ∇u |∇u| ⊗ ∇u |∇u| - 1 ε W (u) ∇u ⊗ ∇w + ∇w ⊗ ∇u |∇u| 2 + 2 ε W (u) ∇u ⊗ ∇u |∇u| 4 < ∇u, ∇w > Denoting N(u) = ∇u⊗∇u |∇u| 2 , we have ε∇W Mu ε (u) = εD 2 : V(u) - 1 ε W (u)N(u) : V(u) + 2 ε div W (u) V(u)∇u |∇u| 2 - 2 ε div W (u) (V(u) : N(u)) ∇u |∇u| 2 ,
where

D 2 : V(u) = ∇ ⊗ ∇ : V(u) = ∑ ij ∂ 2 ij V ij (u) = div(div V(u))
with the same abuse of notation as when one writes div w = ∇ • w. The gradient of W Mu ε can be also expressed as

ε∇W Mu ε (u) = εD 2 : V(u) - 1 ε W (u)N(u) : V(u)+ 2 ε div W (u) V(u)∇u |∇u| 2 -< V(u)∇u/|∇u| 2 , ∇u/|∇u| > ∇u |∇u|
We now give an explicit expression of each previous term.

Evaluation of εD 2 : V(u)

For any operator Λ and any real-valued function u → ρ(u), we have

D 2 : (ρ(u)Λ(u)) = Λ(u) : D 2 ρ(u) + 2 < ∇ρ(u), div(Λ(u)) > +ρ(u)D 2 : Λ(u).
In particular, applying to Λ(u) = V(u)

εD 2 : V(u) = ε 2 D 2 : D 2 u -D 2 : W (u)N(u) = ε 2 ∆ 2 u -W (u)D 2 : N(u) + 2W (u) < ∇u, div(N(u)) > +(D 2 (W (u))) : N(u) .
The last term reads as follows

D 2 (W (u)) : N(u) = W (3) (u)∇u ⊗ ∇u + W (u)∇ 2 u : N(u) = W (3) (u)|∇u| 2 + W (u) < D 2 u∇u, ∇u > |∇u| 2 = ∆ W (u) -W (u) ∆u - < D 2 u∇u, ∇u > |∇u| 2 ,
where we used

∆W (u) = W (u)∆u + W (3) (u)|∇u| 2 .
Recalling that for all vector fields w 1 , w 2 ,

div(w 1 ⊗ w 2 ) = div(w 2 )w 1 + (∇w 1 )w 2 ,
and applying to the estimation of div(N(u)), one gets that

div(N(u)) = div ∇u |∇u| ∇u |∇u| + ∇ ∇u |∇u| ∇u |∇u| . Note that ∇ ∇u |∇u| ∇u |∇u| • ∇u = D 2 u∇u |∇u| 2 -< D 2 u∇u |∇u| 2 , ∇u |∇u| > ∇u |∇u| • ∇u = 0. Therefore 2W (u) < ∇u, div(N(u)) >= 2W (u)|∇u| div ∇u |∇u| = 2W (u) ∆u - < D 2 u∇u, ∇u > |∇u| 2
Lastly,

W (u)D 2 : N(u) = W (u) div (div (N(u))) = W (u) div div ∇u |∇u| ∇u |∇u| + ∇ ∇u |∇u| ∇u |∇u| = W (u) div div ∇u |∇u| ∇u |∇u| + div ∇ ∇u |∇u| ∇u |∇u| therefore εD 2 : V(u) = ε 2 ∆ 2 u -∆W (u) -W (u) ∆u - < D 2 u∇u, ∇u > |∇u| 2 -W (u) div div ∇u |∇u| ∇u |∇u| + div ∇ ∇u |∇u| ∇u |∇u| .

Sum of the first two terms of ε∇W

Mu ε (u) Let I 1 = εD 2 : V(u) - 1 ε W (u)N(u) : V(u). Remark that 1 ε W (u)N(u) : V(u) = 1 ε W (u)N(u) : εD 2 u - 1 ε W (u)N(u) = W (u) < D 2 u∇u, ∇u > |∇u| 2 - 1 ε 2 W (u) = W (u) ∆u - 1 ε 2 W (u) -W (u) ∆u - < ∇ 2 u∇u, ∇u > |∇u| 2
Combining with the previous estimation of εD 2 : V(u), we obtain

I 1 = εD 2 : V(u) - 1 ε W (u)N(u) : V(u) = ε∆ ε∆u - 1 ε W (u) - 1 ε W (u) ε∆u - 1 ε W (u) -W (u) div div ∇u |∇u| ∇u |∇u| + div ∇ ∇u |∇u| ∇u |∇u| .

Estimation of the divergence term

Let I 2 = 2 ε div W (u) V(u)∇u |∇u| 2 -W (u)(V(u) : N) ∇u |∇u| 2 .
On the one hand, with

V(u) ∇u |∇u| 2 = εD 2 u ∇u |∇u| 2 - 1 ε W (u) ∇u |∇u| 2 , we see that div W (u)V(u) ∇u |∇u| 2 = W (u) ε div D 2 u∇u |∇u| 2 - 1 ε div W (u) ∇u |∇u| 2 + W (u) ε < D 2 u∇u, ∇u > |∇u| 2 - 1 ε W (u)
On the other hand,

W (u)(V(u) : N) ∇u |∇u| 2 = W (u) ε < D 2 u∇u, ∇u > |∇u| 2 ∇u |∇u| 2 - 1 ε W (u) ∇u |∇u| 2 and div W (u)(V(u) : N) ∇u |∇u| 2 = W (u) ε < D 2 u∇u, ∇u > |∇u| 2 - 1 ε W (u) + W (u) div ε < D 2 u∇u, ∇u > |∇u| 2 ∇u |∇u| 2 - 1 ε W (u) ∇u |∇u| 2 .
Finally,

I 2 = 2W (u) div D 2 u∇u |∇u| 2 -< D 2 u∇u |∇u| , ∇u |∇u| > ∇u |∇u| 2 = 2W (u) div ∇ ∇u |∇u| ∇u |∇u|
Evaluation of the energy gradient

ε∇W Mu ε (u) = I 1 + I 2 = ε∆µ - 1 ε W (u)µ -W (u) div div ∇u |∇u| ∇u |∇u| -div ∇ ∇u |∇u| ∇u |∇u| , where µ = ε∆u - 1 ε W (u).
whence the L 2 -gradient flow associated with Mugnai's model follows.

Formal asymptotic expansions

We apply in this section the formal method of matched asymptotic expansions to the solution (u ε , µ ε ) of ( 5)

ε 2 ∂ t u = ∆µ -1 ε 2 W (u)µ + W (u)B(u) µ = W (u) -ε 2 ∆u.
. Again, we assume without loss of generality that the isolevel set

Γ(t) = {u ε = 1 2 } is a smooth n -1 dimensional boundary Γ(t) = ∂E(t) = ∂{x ∈ R d ; u ε (x, t) ≥ 1/2}.
In addition, we assume that there exist outer expansions of u ε and µ ε far from the front Γ of the form

u ε (x, t) = u 0 (x, t) + εu 1 (x, t) + ε 2 u 2 (x, t) + O(ε 3 ) µ ε (x, t) = µ 0 (x, t) + εµ 1 (x, t) + ε 2 µ 2 (x, t) + O(ε 3 )
Considering the stretched variable z = d(x,t) ε on a small neighborhood of Γ, we also look for inner expansions of u ε (x, t) and µ ε (x, t) of the form

u ε (x, t) = U(z, x, t) = U 0 (z, x, t) + εU 1 (z, x, t) + ε 2 U 2 (z, x, t) + O(ε 3 ) µ ε (x, t) = W(z, x, t) = W 0 (z, x, t) + εW 1 (z, x, t) + ε 2 W 2 (z, x, t) + O(ε 3 )
As before, we define a unit normal m to Γ and the normal velocity V to the front as

V = -∂ t d(x, t), m = ∇d(x, t), x ∈ Γ.
Let us now expand u ε and µ ε .

Outer expansion: Analogously to [START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF], we obtain

u 0 (x, t) = 1 if x ∈ E(t) 0 otherwise , and
u 1 = u 2 = u 3 = µ 0 = µ 1 = µ 2 = 0.

Matching conditions:

The matching conditions (see [START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF] for more details) imply in particular that 2 , and using ∇ x U.m = 0, it follows that

lim z→+∞ U 0 (z, x, t) = 0, lim z→-∞ U 0 (z, x, t) = 1, lim z→±∞ U 1 (z, x, t) = 0 and lim z→±∞ U 2 (z, x, t) = 0 and lim z→±∞ W 0 (z, x, t) = 0, lim z→±∞ W 1 (z, x, t) = 0 and lim z→±∞ W 2 (z, x, t) = 0 Penalization term B(u) : With ∇u |∇u| = m -ε∇ x U/∂ z U 1 + ε 2 |∇ x U| 2 /(∂ z U)
B(u) = div div ∇u |∇u| ∇u |∇u| -div ∇ ∇u |∇u| ∇u |∇u| = div (∆d∇d) -div ∇ 2 d∇d + O(ε) = ∑ i κ i (π(x)) 1 + zεκ i (π(x) 2 -∑ i κ i (π(x)) 2 (1 + zεκ i (π(x))) 2 + O(ε) = H 2 -A 2 + O(ε).

Inner expansion:

We can derive the asymptotics of the second equation of the system ( 5)

µ = W (u) -ε 2 ∆u as follows      W 0 = W (U 0 ) -∂ 2 z U 0 W 1 = W (U 0 )U 1 -∂ 2 z U 1 -κ∂ z U 0 W 2 = W (U 0 )U 2 -∂ 2 z U 2 + 1 2 W (3) (U 0 )U 2 1 -H∂ z U 1 + z A 2 ∂ z U 0 -∆ x U 0 ,
As for the first equation

ε 2 ∂ t u = ∆µ - 1 ε 2 W (u)µ + W (u)B(u) therefore 0 = ∂ 2 z W 0 -W (U 0 )W 0 , 0 = ∂ 2 z W 1 + H∂ z W 0 -W (U 0 )W 1 + W (3) (U 0 )U 1 W 0 , . and 0 = ∂ 2 z W 2 + H∂ z μ1 -A 2 z∂ z W 0 + ∆ x W 0 + W (U 0 ) H 2 -A 2 - W (U 0 )W 2 + W (3) (U 0 )U 1 W 1 + W (3) (U 0 )U 2 W 0 + 1 2 W (4) (U 0 )U 2 1 W 0

First order:

The two following equations

∂ 2 z W 0 -W (U 0 )W 0 = 0, and W 0 = W (U 0 ) -∂ 2 z U 0 , associated with the boundary conditions lim z→-∞ U 0 (z, x, t) = 1, lim z→+∞ U 0 (z, x, t) = 0,

and lim

z→±∞ W 0 (z, x, t) = 0, admit as solution pair U 0 (z, x) = q(z), and W 0 = 0.

Second order:

The second order gives

∂ 2 z W 1 -W (q)W 1 = 0, W 1 = W (q)U 1 -∂ 2 z U 1 -Hq (z),
which has the solution U 1 = 0, and W 1 = -Hq (z).

Third order:

Using U 0 = q, W 0 = U 1 = 0 and W 1 = -Hq , the first equation can be rewritten as

∂ 2 z W 2 -W (q)W 2 = -H∂ z W 1 -H 2 -A 2 W (q) = H 2 q (z) -H 2 -A 2 q (z) = A 2 q (z),
and implies that

W 2 = A 2 η 2 (z) + c(x, t)q (z),
where η 2 is defined as the solution of

η 2 (z) -W (q(z))η 2 (z) = q (z), with lim z→±∞ η 2 (z) = 0.
Remark that η 2 can be expressed as

η 2 (z) = 1 2 zq (z).
Note that the second equation also reads as

∂ 2 z U 2 -W (q)U 2 = z A 2 q (z) -W 2 = 1 2 z A 2 q (z) -c(x, t)q (z).
In particular, multiplying by q and integrating over R in z shows that c(x, t) = 0. We then deduce that

U 2 = 1 2 A 2 η 1 (z),
where η 1 is defined as the solution of

η 1 (s) -W (q(s))η 1 (s) = sq (s), with lim s→±∞ η 1 (s) = 0.
In conclusion, we have

W 2 = 1 2 A 2 zq (z) and U 2 = 1 2 A 2 η 1 (z).
Fourth order and estimation of the velocity V:

We can now explicit the term of order 1 in ε of B(u). Indeed we have U(z, x, t) = q(z) + O(ε 2 ), and as ∇ x q(z) = 0, we have

B(u) = div div ∇u |∇u| ∇u |∇u| -div ∇ ∇u |∇u| ∇u |∇u| = div (∆d∇d) -div ∇ 2 d∇d + O(ε 2 ) = ∑ i κ i (π(x)) 1 -zεκ i (π(x)) 2 -∑ i κ i (π(x)) 2 (1 -zεκ i (π(x))) 2 + O(ε 2 ) = H 2 -A 2 -ε2z H A 2 -Θ 3 + O(ε 2 ), where Θ 3 = ∑ i k i (π(x)) 3 .
The fourth order of the first equation now reads

-Vq = ∂ 2 z W 3 -W (q)W 3 -W (3) (q)U 2 W 1 + H∂ z W 2 -A 2 z∂ z W 1 + ∆ x W 1 -zW (q)2 H A 2 -Θ 3 = ∂ 2 z W 3 -W (q)W 3 + 1 2 W (3) (q)H A 2 η 1 q + 1 2 A 2 H 3zq + q -(∆ Γ H) q -2 H A 2 -Θ 3 zq = ∂ 2 z W 3 -W (q)W 3 + 1 2 W (3) (q)H A 2 η 1 q + - 1 2 A 2 H + 2B 3 zq + 1 2 A 2 Hq -(∆ Γ H) q
Multiplying by q and integrating over R leads to

V = - 1 S 1 2 A 2 HS + - 1 2 A 2 H + 2Θ 3 R zq q dz + 1 2 A 2 H R W (3) (q)η 1 (q ) 2 dz -∆ Γ HS ,
where S = R q (z) 2 dz. Remark also that

R zq q dz = 1 2 R z((q ) 2 ) dz = - 1 2 R (q ) 2 dz = - 1 2 S.
Moreover, recall that η 1 satisfies

η 1 -W (q)η 1 = zq η 1 -W (q)η 1 -W (3) (q)q η 1 = (zq ) , then we have R W (3) (q)η 1 (q ) 2 dz = R η 1 -W (q)η 1 q dz - R (zq ) q dz = - 1 2 S,
and we conclude that

V = ∆ Γ H + Θ 3 - 1 2 A 2 H.

Approximating the Willmore flow with Esedoglu-Rätz-Röger's energy

We now consider the following variant of the Esedoglu-Rätz-Röger's energy, which we introduced in Section 2.4.3

W EsRäRö ε (u) = 1 2ε Ω ε∆u - W (u) ε 2 dx + βJ ε (u)
where the penalization term J ε (u) reads The aim of this section is to derive and study the PDE obtained as the L 2 -gradient flow of W EsRäRö ε (u). We will show that the flow is equivalent to the phase field system [START_REF] Barrett | Parametric approximation of Willmore flow and related geometric evolution equations[END_REF] Note that this system coincides with the classical one, up to the addition of a penalty term -β L(u).

J ε (u) = 1 ε 1+α Ω ε∇ 2 u : N(u) - W (u)
                 ε 2 ∂ t u = ∆µ -1 ε 2 W (u)µ -β L(u) µ = W (u) -ε 2 ∆u. ξ ε = ε∇ 2 u : N(u) -W (u) ε L(u) = 2ε 1-α N(u) : ∇ 2 ξ ε -1 ε 2 W (u)ξ ε + 2 div ∇u |∇u| ∇u |∇u| • ∇ξ ε + B(u)ξ ε B(u) = div div ∇u |∇u| ∇u |∇u| -div ∇ ∇u |∇u| ∇u |∇u|
The well-posedness of the phase field model ( 6) at fixed parameter ε is open, and requires presumably a regularization as done numerically in [START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF].

By formal arguments involving matched asymptotic expansions again, we will show that this approximating flow is expected to converge, as ε goes to zero, to the Willmore flow in dimension N ≥ 2, at least whenever α = 0 or α = 1. More precisely we will show the Claim 3.7 In a suitable regime provided by the method of matched asymptotic expansions, the normal velocity of the 1 2 -front Γ(t) = ∂E(t) associated with a solution (u α ε , µ α ε , ξ α ε ) to Esedoglu-Rätz-Röger's phase field model (6) in both cases α = 0 and α = 1 is the Willmore velocity

V = ∆ Γ H + A 2 H - H 3 2 .
In addition, for α = 0:

         u 0 ε (x, t) = q d(x,E(t)) ε + ε 2 A 2 -H 2 1+2β η 1 d(x,E(t) ε + O(ε 3 ) µ 0 ε (x, t) = -εHq d(x,E(t) ε + ε 2 H 2 -2β 2 A 2 -H 2 1+2β η 2 d(x,E(t) ε + O(ε 3 ), ξ 0 ε (x, t) = ε 2 A 2 -H 2 1+2β η 2 d(x,E(t) ε + O(ε 2 )
where η 2 (z) = zq (z) is a profile function. For α = 1:

         u 1 ε (x, t) = q d(x,E(t)) ε + O(ε 3 ) µ 1 ε (x, t) = -εHq d(x,E(t) ε + 2ε 2 A 2 η 2 d(x,E(t) ε + O(ε 3 ), ξ 1 ε (x, t) = ε 2 (2 A 2 -H 2 ) 4β η 2 d(x,E(t) ε + O(ε 3 ) Remark 3.8
The previous claim gives indications on the design of a numerical scheme for simulating the Esedoglu-Rätz-Röger's flow in the cases α = 0, 1. Clearly, the flow acts at the second order for u in the case α = 0, and not less than at the third order (at least) whenever α = 1. This implies that capturing with accuracy the motion of the interface should be much more delicate when α = 1.

Derivation of the

                 ε 2 ∂ t u = ∆µ -1 ε 2 W (u)µ -β L(u) µ = W (u) -ε 2 ∆u. ξ ε = ε∇ 2 u : N(u) -W (u) ε L(u) = 2ε 1-α N(u) : ∇ 2 ξ ε -1 ε 2 W (u)ξ ε + 2 div ∇u |∇u| ∇u |∇u| • ∇ξ ε + B(u)ξ ε B(u) = div div ∇u |∇u| ∇u |∇u| -div ∇ ∇u |∇u| ∇u |∇u| Proof: Since ξ ε (u) = ε∇ 2 u : N(u) - W (u) ε = ε ∇ 2 u ∇u |∇u| • ∇u |∇u| - 1 ε W (u), one has that ξ ε (u)(w) = ε ∇ 2 w : N(u) + 2 ∇ 2 u : ∇w ⊗ ∇u |∇u| 2 -2∇ 2 u : N(u) ∇u • ∇w |∇u| 2 - 1 ε W (u)w.
The gradient of J ε (u) follows, recalling that

P u = I d -N(u): ∇J ε (u) = 2 ε α ∇ 2 : [N(u)ξ ε ] - 1 ε 2 W (u)ξ ε -2 div ξ ε (u)∇ 2 u∇u |∇u| 2 + 2 div ∇ 2 u : N(u)ξ ε ∇u |∇u| 2 = 2 ε α ∇ 2 : [N(u)ξ ε ] - 1 ε 2 W (u)ξ ε -2 div ξ ε (u)P u ∇ 2 u∇u |∇u| 2
More precisely, using

∇ 2 : [N(u)ξ ε ] = (∇ 2 : N(u))ξ ε + 2 div(N(u)) • ∇ξ ε + N(u) : ∇ 2 ξ ε , div(N(u)) = div ∇u
∇J ε (u) = 2 ε α N(u) : ∇ 2 ξ ε - 1 ε 2 W (u)ξ ε + 2 div ∇u |∇u| ∇u |∇u| • ∇ξ ε + B(u)ξ ε .
from which the L 2 -gradient flow of Esedoglu-Rätz-Röger's 's model follows.

Asymptotic analysis

In this section, we perform the formal method of matched asymptotic expansions to the system solution

(u ε , µ ε , ξ ε ) of            ε 2 ∂ t u = ∆µ -1 ε 2 W (u)µ -β L(u) µ = W (u) -ε 2 ∆u. ξ ε = ε∇ 2 u : N(u) -W (u) ε L(u) = 2ε 1-α N(u) : ∇ 2 ξ ε -1 ε 2 W (u)ξ ε + 2 div ∇u |∇u| ∇u |∇u| • ∇ξ ε + B(u)ξ ε ,
in both cases α = 0 and α = 1.

As previously, we assume that the 1/2-isolevel set of u ε is a smooth (n -1)-dimensional interfaces Γ(t) defined as the boundary of a set

E(t) = x ∈ R d ; u ε (x, t) ≥ 1/2 .
We assume that there exist outer expansions of u ε , µ ε and ξ ε far from the front Γ of the form

     u ε (x, t) = u 0 (x, t) + εu 1 (x, t) + ε 2 u 2 (x, t) + O(ε 3 ) µ ε (x, t) = µ 0 (x, t) + εµ 1 (x, t) + ε 2 µ 2 (x, t) + O(ε 3 ) ξ ε (x, t) = 1 ε ξ -1 (x, t) + ξ 0 (x, t) + εξ 1 (x, t) + ε 2 ξ 2 (x, t) + O(ε 3 )
Considering the stretched variable z = d(x,t) ε on a small neighborhood of Γ, we also look for inner expansions of u ε (x, t), µ ε (x, t) and ξ ε (x, t) of the form

     u ε (x, t) = U(z, x, t) = U 0 (z, x, t) + εU 1 (z, x, t) + ε 2 U 2 (z, x, t) + O(ε 3 ) µ ε (x, t) = W(z, x, t) = W 0 (z, x, t) + εW 1 (z, x, t) + ε 2 W 2 (z, x, t) + O(ε 3 ) ξ ε (x, t) = Φ(z, x, t) = ε -1 Φ -1 (z, x, t) + Φ 0 (z, x, t) + εΦ 1 (z, x, t) + ε 2 Φ 2 (z, x, t) + O(ε 3 )
In particular, remark that the third equation of ( 6) yields

Φ(z, x, t) = 1 ε   ∂ 2 zz U 1 + ε 2 |∇ x U| 2 (∂ z U) 2 -W (U)   + ε ∂ z |∇ x U| 2 ∂ z U 1 + ε 2 |∇ x U| 2 (∂ z U) 2 -1
.

As before, it can be observed for the outer expansions that

u 0 (x, t) = 1 if x ∈ E(t) 0 otherwise , and
u 1 = u 2 = u 3 = µ 0 = µ 1 = µ 2 = ξ -1 = ξ 0 = ξ 2 = 0.
The matching conditions imply the following boundary conditions on the inner expansions:

lim z→+∞ U 0 (z, x, t) = 0 lim z→-∞ U 0 (z, x, t) = 1 , lim z→±∞ U i (z, x, t) = 0 for i ∈ {1, 2},
and lim z→±∞ W i (z, x, t) = 0, for i ∈ {0, 1, 2}, and lim

z→±∞ Φ i (z, x, t) = 0, for i ∈ {-1, 0, 1, 2}.
Inner expansion with α = 0 : This paragraph is devoted to the derivation of the expression of the inner expansion in the special case α = 0.

First order:

We have the following system

     0 = ∂ 2 z W 0 -W (U 0 )W 0 -2β [∂ zz Φ -1 -W (U 0 )Φ -1 ] W 0 = W (U 0 ) -∂ 2 z U 0 Φ -1 = ∂ 2 z U 0 -W (U 0 )
, which admits the solution triplet U 0 = q(z), W 0 = 0, and Φ -1 = 0.

Second order:

At second order, we obtain

     0 = ∂ 2 z W 1 -W (q)W 1 -2β [∂ zz Φ 0 -W (q)Φ 0 ] W 1 = W (q)U 1 -∂ 2 z U 1 -Hq Φ 0 = ∂ 2 z U 1 -W (q)U 1 ,
whose solution is given by U 1 = 0, W 1 = -Hq , and Φ 0 = 0.

Third order:

At third order,

     0 = ∂ 2 z W 2 -W (U 0 )W 2 + H∂ z W 1 -2β [∂ zz Φ 1 -W (q)Φ 1 ] W 2 = W (q)U 2 -∂ 2 z U 2 + A 2 zq Φ 1 = ∂ 2 z U 2 -W (q)U 2 ,
and we are now looking for a system of solutions of the form

W 2 = c W (x, t)η 2 (z), U 2 = c U (x, t)η 1 (z), and Φ 1 = c Φ (x, t)η 2 (z),
where the two profiles η 1 and η 2 are solutions, respectively, of

η 1 -W (q)η 1 = zq and η 2 -W (q)η 2 = q .
Furthermore, the first equation gives

0 = c W η 2 -W (q)η 2 -H 2 q -2c Φ η 2 -W (q)η 2 = c W -H 2 -2βc Φ q ,
the second equation implies that

1 2 c W zq = -c U η 1 -W (q)η 1 + A 2 zq = (-c U + A 2 )zq .
and the third equation shows that

1 2 c Φ zq = c U η 1 -W (q)η 1 = (c U )zq .
This provides a linear system

c W -β2c Φ = H 2 , c W + 2c U = 2 A 2 and c Φ = 2c U ,
which admits as solutions

c W = H 2 -2β 2 A 2 -H 2 1 + 2β , c U = A 2 -H 2 /2 1 + 2β , and c Φ = 2 A 2 -H 2 1 + 2β . Therefore,          W 2 = H 2 -2β 2 A 2 -H 2 1+2β η 2 , U 2 = A 2 -H 2 /2 1+2β η 1 , Φ 1 = 2 A 2 -H 2 1+2β η 2 .
Fourth order and estimation of the velocity V:

The fourth order reads as follows

-Vq = ∂ 2 z W 3 -W (q)W 3 -W (3) (q)U 2 W 1 + H∂ z W 2 -A 2 z∂ z W 1 + ∆ x W 1 -2β ∂ zz Φ 2 -W (q)Φ 2 + 2H∂ z Φ 1 = ∂ 2 z (W 3 -2βΦ 2 ) -W (q)(W 3 -2βΦ 2 ) + H A 2 -H 3 /2 1 + 2β W (3) (q)η 1 q -∆ Γ Hq + H 3 /2 + A 2 H -Hβ 2 A 2 -H 2 1 + 2β zq + H 3 /2 -Hβ 2 A 2 -H 2 1 + 2β q Recalling that      R (q (z)) 2 dz = S R zq q dz = -1 2 S R W (3) (q)η 1 (q ) 2 dz = -1 2 S
, multiplying the last equation by q and integrating over R leads to

V = ∆ Γ H + A 2 H 2 - H 3 4 + H A 2 -H 3 /2 2(1 + 2β) (1 + 2β) = ∆ Γ H + A 2 H - 1 2 H 3
Inner expansion with α = 1 :

We are now looking for the inner expansion of the PDE system (6) in the special case α = 1 :

First order:

We have the following system

     0 = -2β [∂ zz Φ -1 -W (U 0 )Φ -1 ] W 0 = W (U 0 ) -∂ 2 z U 0 Φ -1 = ∂ 2 z U 0 -W (U 0 )
, whose solution is given by U 0 = q(z), W 0 = 0, and Φ -1 = 0.

Second order:

At second order

     0 = ∂ 2 zz W 0 -W (q)W 0 -2β ∂ 2 zz Φ 0 -W (q)Φ 0 W 1 = W (q)U 1 -∂ 2 zz U 1 -Hq Φ 0 = ∂ 2 zz U 1 -W (q)U 1 ,
which admits as solution U 1 = 0, W 1 = -Hq , and Φ 0 = 0.

Third order:

At third order

     0 = ∂ 2 zz W 1 -W (q)W 1 -2β [∂ zz Φ 1 -W (q)Φ 1 ] W 2 = W (q)U 2 -∂ 2 zz U 2 + A 2 zq Φ 1 = ∂ 2 z U 2 -W (q)U 2 ,
whose solution triplet is

U 2 = 0, W 2 = A 2 zq , and Φ 1 = 0.
Fourth order:

From 0 = ∂ 2 z W 2 -W (q)W 2 + H∂ z W 1 -2β ∂ zz Φ 2 -W (q)Φ 2 , we deduce that ∂ zz Φ 2 -W (q)Φ 2 = (2 A 2 -H 2 ) 2β q ,
and then

Φ 2 = (2 A 2 -H 2 ) 2β η 2 = (2 A 2 -H 2 ) 4β zq .
Last order and estimation of the velocity V:

We have

-Vq = ∂ 2 z W 3 -W (q)W 3 -W (3) (q)U 2 W 1 + H∂ z W 2 -A 2 z∂ z W 1 + ∆ x W 1 -2β ∂ zz Φ 3 -W (q)Φ 3 + 2H∂ z Φ 2 = ∂ 2 z (W 3 -2βΦ 3 ) -W (q)(W 3 -2βΦ 3 ) -∆ Γ Hq + 2 A 2 H -Hβ 2 A 2 -H 2 β zq + A 2 H -Hβ 2 A 2 -H 2 β q = ∂ 2 z (W 3 -2βΦ 3 ) -W (q)(W 3 -2βΦ 3 ) -∆ Γ Hq + H 3 zq + -A 2 H + H 3 q .
As previously, this shows that the velocity V of the interfaces equals the Willmore velocity

V = ∆ Γ H + A 2 H - H 3 2 .

Remark 3.10

The analysis of the asymptotic behavior for α non integer is more delicate because it requires studying non integer orders of ε and it is far from being clear how integer and non integer scales may combine. As for integer values of α > 1, a careful study at higher orders of ε should be possible but is out of the scope of the present paper.

2D and 3D numerical simulations for the classical and Mugnai's diffuse flows

There is an important literature on numerical methods for the approximation of interfaces evolving by a geometric law. They can be roughly classified into three categories: parametric methods [START_REF] Deckelnick | Discrete anisotropic curvature flow of graphs[END_REF][START_REF] Rusu | An algorithm for the elastic flow of surfaces[END_REF][START_REF] Deckelnick | Computation of geometric partial differential equations and mean curvature flow[END_REF][START_REF] Barrett | On the parametric finite element approximation of evolving hypersurfaces in R 3[END_REF][START_REF] Barrett | A variational formulation of anisotropic geometric evolution equations in higher dimensions[END_REF][START_REF] Polthier | Polyhedral surfaces of constant mean curvature[END_REF][START_REF] Hsu | Minimizing the squared mean curvature integral for surfaces in space forms[END_REF], level-set formulations [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF][START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF][START_REF] Osher | Geometric Level Set Methods in Imaging, Vision and Graphics[END_REF][START_REF] Evans | Motion of level sets by mean curvature[END_REF][START_REF] Chen | Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations[END_REF], phase-field approaches [START_REF] Modica | Un esempio di Γ--convergenza[END_REF][START_REF] Chen | Generation and propagation of interfaces for reaction-diffusion equations[END_REF][START_REF] Bellettini | Anisotropic motion by mean curvature in the context of Finsler geometry[END_REF][START_REF] Paolini | An efficient algorithm for computing anisotropic evolution by mean curvature[END_REF]. See for instance [START_REF] Deckelnick | Computation of geometric partial differential equations and mean curvature flow[END_REF] for a complete review (with a particular emphasis on the mean curvature flow, but fourthorder flows are also addressed) and a comparison between the different strategies. In the context of fourth order geometric evolution equations, in particular the Willmore flow, parametric approaches have been proposed in [START_REF] Dziuk | Computational parametric Willmore flow[END_REF][START_REF] Barrett | A parametric finite element method for fourth order geometric evolution equations[END_REF][START_REF] Barrett | Parametric approximation of Willmore flow and related geometric evolution equations[END_REF] for curves and surfaces using a semi-implicit finite element method.

In [START_REF] Olischläger | A nested variational time discretization for parametric Willmore flow[END_REF], a fully implicit approach via a variational formulation is also analyzed for the approximation of anisotropic Willmore flow. The level set methods have been applied for the first time in [START_REF] Droske | A level set formulation for Willmore flow[END_REF]. Concerning the phase field approach, semi-implicit schemes including standard finite element differences, finite elements, and Fourier spectral methods are developed in [START_REF] Du | A phase field approach in the numerical study of the elastic bending energy for vesicle membranes[END_REF][START_REF] Du | Simulating the deformation of vesicle membranes under elastic bending energy in three dimensions[END_REF][START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF] and analyzed in [START_REF] Du | Convergence of numerical approximations to a phase field bending elasticity model of membrane deformations[END_REF]. A fully implicit scheme coupled with a finite element method has been more recently introduced in [START_REF] Franken | A phase field based PDE constraint optimization approach to time discrete Willmore flow[END_REF] via a variational formulation. An adaptation to fourth order geometric evolution equations of the Bence-Merriman-Osher algorithm [START_REF] Bence | Diffusion generated motion by mean curvature[END_REF] is also proposed in [START_REF] Esedoglu | Threshold dynamics for high order geometric motions[END_REF]. Let us finally mention the discrete methods involving surface triangulations and discrete curvature operators [START_REF] Hsu | Minimizing the squared mean curvature integral for surfaces in space forms[END_REF][START_REF] Bobenko | Discrete Willmore Flow[END_REF][START_REF] Wardetzky | Discrete quadratic curvature energies[END_REF].

In this paper, we will consider a quite different and new scheme to solve both the classical and Mugnai's phase field systems (3) and ( 5). The simulations can be compared with those obtained by Esedoglu, Rätz and Röger in [START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF] for their phase field system (6) (actually a variant of it, see Section 2.4.3), and for Bellettini's phase field system (4).

Here, we use an implicit scheme to ensure the decreasing of the diffuse Willmore energy, and a Fourier spectral method in order to get high accuracy approximation in space. At each step time, it is necessary to solve a nonlinear equation. A Newton algorithm like in [START_REF] Franken | A phase field based PDE constraint optimization approach to time discrete Willmore flow[END_REF] appears to be very efficient in practice, but not in accordance with a Fourier spectral discretization, so we opted for a fixed point approach.

New numerical schemes for the approximation of classical and Mugnai's flows

Classical diffuse approximation flow

We introduce a new scheme to approximate numerically some solutions of the phase field system

   ∂ t u = 1 αε 2 ∆µ - 1 αε 4 W (u)µ µ = αW (u) -αε 2 ∆u,
where α is a positive constant. The particular phase field system (3) corresponding to the classical diffuse Willmore flow is obtained for α = 1.

We compute the solution for any time t ∈ [0, T] in a box Ω = [0, 1] N with periodic boundary conditions. We use a Euler semi-implicit discretization in time:

u n+1 = δ t 1 αε 2 ∆µ n+1 -1 αε 4 W (u n+1 )µ n+1 + u n µ n+1 = αW (u n+1 ) -αε 2 ∆u n+1 ,
where δ t is the time step, u n and µ n are the approximations of the solutions u and µ, respectively, evaluated at time t n = n δ t . The system can be written as

u n+1 -δ t αε 2 ∆µ n+1 = E µ n+1 + αε 2 ∆u n+1 = F, with E = u n -δ t αε 4 W (u n+1 )µ n+1 , F = αW (u n+1 ). Therefore, u n+1 + δ t ∆ 2 u n+1 = E + δ t αε 2 ∆F µ n+1 + δ t ∆ 2 µ n+1 = F -αε 2 ∆E.
Thus, (u n+1 , µ n+1 ) is the solution of the nonlinear equation

u n+1 µ n+1 = φ u n+1 µ n+1 , (7)
where

φ u n+1 µ n+1 = I d + δ t ∆ 2 -1 I d δ t αε 2 ∆ -αε 2 ∆ I d u n -δ t αε 4 W (u n+1 )µ n+1 αW (u n+1 )
A natural way to approximate the solution (u n+1 , µ n+1 ) to ( 7) is a fixed point iterative method. The space discretization is built with Fourier series. It has the advantage of preserving a high order approximation in space while allowing a fast and simple processing of the homogeneous operator

G = I d + δ t ∆ 2 -1 I d δ t αε 2 ∆ -αε 2 ∆ I d = I d -δ t αε 2 ∆ αε 2 ∆ I d -1
.

In practice, the solutions u(x, t n ) and µ(x, t n ) at time t n = nδ t are approximated by the truncated Fourier series :

u n P (x) = ∑ p ∞ ≤P
u n p e 2iπx•p , and µ n

P (x) = ∑ p ∞ ≤P µ n p e 2iπx•p
where p ∞ = max 1≤i≤N |p i |, P is the maximal number of Fourier modes in each direction, and the coefficients u n p , µ n p are derived from a prior Fourier decomposition of

u n -δ t αε 4 W (u n+1 )µ n+1 αW (u n+1 ) com-
bined with an application in the Fourier domain of the operator G. More precisely, the fixed point algorithm that we propose reads as follows:

Algorithm 1 Initialization : u n+1 0 = u n , µ n+1 0 = µ n While u n+1 k+1 -u n+1 k + µ n+1 k+1 -µ n+1 k > 10 -8
, perform the loop on k:

1) Compute

h n P = u n - δ t αε 4 W (u n+1 k )µ n+1 k
, and hn P = αW (u n+1 k )

2) Using the Fast Fourier Transform, compute the truncated Fourier series of h n P and hn P :

h n P (x) = ∑ p ∞ ≤P
h n p e 2iπx•p , and hn

P (x) = ∑ p ∞ ≤P hn p e 2iπx•p 3) Compute u n+1 k+1 (x) = ∑ p ∞ ≤P (u k+1 ) n p e

End

Note that the semi-implicit scheme

u n+1 = δ t 1 ε 2 α ∆µ n+1 -1 ε 4 α W (u n )µ n + u n µ n+1 = αW (u n ) -αε 2 ∆u n+1 , implies the following scheme on u n u n+1 = I d + δ t ∆ 2 -1 u n + δ t ε 2 ∆W (u n ) + δ t ε 2 W (u n ) ∆u n - 1 ε 2 W (u n-1 ) ,
which is expected to be stable under a condition of the form

δ t ≤ C min ε 2 δx 2 , ε 4
where δx = 1/(2P ) and C is a constant depending only the double-well potential W.

In practice, using fixed point iterations instead of a semi-implicit scheme appears more accurate numerically. This can be justified with the following proposition:

Proposition 4.1 Algorithm 1 converges locally under the assumptions

max [αM 2 ] 2 + 2[ δ t ε 4 M 3 (M 1 + N 3/2 π 2 ε 2 δ 5/2 x )] 2 , 2[ δ t αε 4 M 2 ] 2 < 1 ( 8
)
where M i = sup s∈[0,1] |W (i) (s)|.

Proof: We look for the conditions such that

Dφ(u n+1 , µ n+1 )(δ u , δ µ ) 2 < (δ u , h µ ) 2 ,
where the differential of φ is such that

Dφ(u n+1 , µ n+1 )(δ u , δ µ ) = I d -δ t αε 2 ∆ αε 2 ∆ I d -1 -δ t αε 4 W (3) (u n+1 )µ n+1 δ u + W (2) (u n+1 )δ µ αW (2) (u n+1 )δ u .
Note that the eigenvalues of the operator

I d -δ t αε 2 ∆ αε 2 ∆ I d are λ p = 1 ± 4π 2 i δ t |p| 2 , for p ∞ ∈ [0, P ].
In particular, this implies that

I d -δ t αε 2 ∆ αε 2 ∆ I d -1 ≤ 1.
Moreover, remark also that

|µ n+1 | = |αW (u n+1 ) -αε 2 ∆u n+1 | ≤ α M 1 + N 3/2 π 2 ε 2 δ 5/2 x . It follows that    |αW (2) (u n+1 )δ u | ≤ αM 2 |δ u | W (3) (u n+1 )µ n+1 δ u + W (2) (u n+1 )δ µ ≤ αM 3 (M 1 + N 3/2 π 2 ε 2 δ 5/2 x )|δ u | + M 2 |δ µ | and then -δ t αε 4 W (3) (u n+1 )µ n+1 δ u + W (2) (u n+1 )δ µ αW (2) (u n+1 )δ u 2 ≤ max [αM 2 ] 2 + 2[ δ t ε 4 M 3 (M 1 + N 3/2 π 2 ε 2 δ 5/2 x )] 2 , 2[ δ t αε 4 M 2 ] 2 (δ u , δ µ ) 2 ,
which concludes the convergence proof of the fixed-point iteration procedure if conditions (8) above are fulfilled.

Mugnai's flow

We now use a similar scheme for the following generalization of Mugnai's phase field system:

   ∂ t u = 1 ε 2 α ∆µ - 1 ε 4 α W (u)µ + B(u) µ = αW (u) -αε 2 ∆u, with B(u) = W (u)B(u) αε 4
. The exact Mugnai's flow (5) corresponds to the choice α = 1 (up to a time rescaling). We use again a semi-implicit scheme:

u n+1 = δ t 1 ε 2 α ∆µ n+1 -1 αε 4 W (u n+1 )µ n+1 + B(u n ) + u n µ n+1 = αW (u n+1 ) -ε 2 α∆u n+1 ,
where the penalization term B(•) is treated explicitly. We use a fixed point iteration to approximate the solution pair (u n+1 , µ n+1 ) to the system:

(u n+1 , µ n+1 ) = φ(u n+1 , µ n+1 ) = I d + δ t ∆ 2 -1 I d δ t αε 2 ∆ -αε 2 ∆ I d u n -δ t αε 4 W (u n+1 )µ n+1 + δ t B(u n ) αW (u n+1
) .

For it is highly singular, the penalization term B(•) needs to be regularized to avoid numerical errors.

Observing that

B(u) = W (u) ∇ ∇u |∇u| 2 -div ∇u |∇u| 2 -curl curl ∇u |∇u| • ∇u |∇u|
we consider the regularized penalization term

B σ (u) = W (u) |∇ν u,σ | 2 -|div ν u,σ | 2 -curl (curl (ν u,σ )) • ν u,σ
where ν u,σ = ∇u √ |∇u| 2 +σ 2 with σ a small regularization parameter. In particular, the positivity

|∇ν u,σ | 2 -|div ν u,σ | 2 ≥ 0,
is ensured, which is in accordance with the continuous case. In practice, finite differences are used for the numerical evaluation of B σ (u). Finally, we propose the following algorithm:

Algorithm 2 Initialization : u n+1 0 = u n , µ n+1 0 = µ n While u n+1 k+1 -u n+1 k + µ n+1 k+1 -µ n+1 k > 10 -8 ,
[αM 2 ] 2 + 2[ δ t ε 4 M 3 (M 1 + N 3/2 π 2 ε 2 δ 5/2 x )] 2 , 2[ δ t αε 4 M 2 ] 2 < 1 ( 9
)
where M i = sup s∈[0,1] |W (i) (s)|.

Proof: Being the penalty term B σ (u) treated explicitly, it does not appear in the expression of the differential of φ, thus the same proof as for Proposition 4.1 can be used.

Numerical simulations of the classical flow

The following simulations have been realized using Matlab. The isolevel sets Γ(t) = {x : u(x, t) = 1 2 } are computed and drawn using the Matlab functions contour in 2D and isosurface in 3D. We use the double-well potential W(s) = 1 2 s 2 (1s) 2 and consider the PDE system

∂ t u = ∆µ -1 ε 2 W (u)µ µ = 1 ε 2 W (u) -∆u.
with initial conditions u(x, 0) and µ(x, 0) of the form

   u(x, 0) = γ d(x,E) ε µ(x, 0) = -1 ε ∆d(x, E)γ d(x,E) ε

Evolution of a circle

The first test plotted in Figure 2 illustrates the good behavior of our scheme with respect to the exact solution. The initial set E is a circle of radius R 0 = 0.15. The continuous Willmore flow preserves the circle yet increases the radius according the law

R(t) = R 4 0 + 2t 1/4 .
The left picture in Figure 2 represents the interfaces Γ(t) at different times t obtained with the following numerical parameters: P = 2 7 , ε = 2/P and δ t = ε 2 2P 2 . The right picture in Figure 2 depicts the error between the numerical radius R ε (t) and the theoretical radius R(t), at different times and for two different values of ε (the other parameters are kept unchanged). It is reasonable to believe that this experiment illustrates the numerical convergence of R ε (t) to R(t) as ε goes to zero. Evolution of two disjoint circles and formation of singularities One of our motivations in this study is to understand and observe the behavior of the diffuse Willmore solution in the situations where singularities appear. As it was discussed in Section 2.3, this may happen for instance with the classical approximation flow. We consider as initial set Ω 0 the union of two disjoint circles of radius R = 0.15. Each circle should evolve as a circle with increasing radius, up to the contact occurs. To the best of our knowledge, the theoretical Willmore flow is not clearly defined after this critical collision time. Therefore, the asymptotic limit of the solution t → u ε (•, t) as ε goes to 0 could be a good candidate for the definition of a weak Willmore flow. However, different behaviors of t → u ε (•, t) have been observed in the literature. For instance, the two circle merge in [START_REF] Franken | A phase field based PDE constraint optimization approach to time discrete Willmore flow[END_REF] whereas a crossing of interfaces appears at collision time in [START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF].

We plot on Figure 3 the graph of t → u ε (•, t) computed for different values of ε. We choose for the other parameters: P = 2 7 , δ t = 1/P -4 . In the first experiment obtained with ε = 5/P, the two circles merge. In contrast, a crossing of interfaces appear for the cases ε = 3/P and ε = 1.5. More precisely, we can distinguish three different periods in the last two experiments: in the first period, both circles evolve independently one from the other. The second period begins when the distance between the two circles is about the size of the diffuse interfaces and the formation of a crossing is observed. This corresponds to a solution of the Allen-Cahn equation with unsmooth nodal set. After contact, the interfaces continues to evolve while the crossing seems to be numerically stable and does not influence the interface evolution. More precisely, the interface Γ(t) seems to converge to a growing eight, which is one of the closed planar elasticae described in Langer and Singer's work [START_REF] Langer | Curve straightening and a minimax argument for closed elastic curves[END_REF].

Numerical examples of saddle-shaped solutions of Allen-Cahn equation

We already mentioned in Section 2.3 the existence result due to Dang, Fife, and Peletier [START_REF] Dang | Saddle solutions of the bistable diffusion equation[END_REF] of an entire solution in the plane to the Allen-Cahn equation whose nodal set coincides with {(x, y), xy = 0} (it can be generalized to every even dimension [START_REF] Cabré | Saddle-shaped solutions of bistable diffusion equations in all of R 2m[END_REF]). By restricting to a sector and using consecutive reflections, it is possible to build solutions whose nodal set has an arbitrary number of branches with the property of dihedral symmetry, i.e. of equal angle between two consecutive branches [START_REF] Gui | Hamiltonian identities for elliptic partial differential equations[END_REF]. Actually, by a result of Hartman and Wintner [START_REF] Hartman | On the local behavior of solutions of non-parabolic partial differential equations[END_REF], saddle-shaped solutions must satisfy the equal angle property. We illustrate numerically in Figure 4, left, two solutions with different branching degrees. These examples are classical, and have been previously obtained by various authors [40,[START_REF] Lowengrub | Phase-field modeling of the dynamics of multicomponent vesicles: Spinodal decomposition, coarsening, budding, and fission[END_REF][START_REF] Esedoglu | Colliding Interfaces in Old and New Diffuse-interface Approximations of Willmore-flow[END_REF].

To be complete, let us mention that 2k-ended solutions, i.e. solutions whose nodal set coincides outside any compact set with the union of 2k straight lines which cross at the origin, do exist without the dihedral symmetry requirement [START_REF] Del Pino | Multiple-end to the Allen-Cahn equation in R 2[END_REF]. In the particular case of 4-ended solutions, the result can even be proved for arbitrary angles between the lines [START_REF] Kowalczyk | Four ended solutions to the Allen-Cahn equation on the plane[END_REF]. Of course, by Hartman and Wintner's result, the nodal set itself cannot self-intersect at the origin if the dihedral symmetry does not hold, but remains smooth instead. This is illustrated in Figure 4, right, where the 1/2-isolevel line has been represented

Comparison between phase field and parametric approaches

We observed previously that the evolution of two disjoint circles after contact and creation of a crossing is similar to the evolution of a eight-like single curve. To highlight this point, we tested on the same configurations both the evolution provided by the classical diffuse flow and the evolution with respect to a discrete parametric Willmore flow. In particular, we consider two different initial conditions corresponding, respectively, to the union of two or three contiguous circles. The parametric choice of Γ(0) is illustrated on Figure 5 and corresponds to using a single smooth C 1,1 curve that covers two or three circles, respectively. The discrete parametric Willmore flow is computed with a finite element method as proposed and analyzed by Dziuk in [START_REF] Dziuk | Computational parametric Willmore flow[END_REF]. The phase field simulations are done with the set of parameters: P = 2 7 , ε = 1/P and δ t = εP -2 /10. Numerical results are shown on Figure 6. As expected, these two different approaches give very similar results. This suggests that the interface obtained by a phase field approximation converges, after apparition of a singularity, to an interface which evolves as a regular parametric Willmore flow. This is actually very much in favor of a varifold interpretation, at least on this example, of both flows. What really cares is the support, and its geometry, and not the fact that it is seen either as an isolevel set or as a parametrized set.

The experiments on Figure 7 are in the same spirit. On the first line, we illustrate the evolution of two phases forming a circle cut by a straight line. Note that the circle seems to evolve independently of the line. The second situation is quite similar with two disjoint circles cut by a line, and the same conclusion holds.

Evolution by the classical diffuse Willmore flow of two contiguous circles

We now compare in Figure 8 the evolution by a discrete parametric Willmore flow of the different interfaces obtained from three different initial parameterizations of two contiguous circles. In the first parametrization (in blue), both circles are parametrized independently. The second parameterization (in magenta) corresponds to the covering of the two circles by a unique smooth parametric curve that self-crosses at the origin. The third parameterization (in red) corresponds to the singular curve that does not cross the horizontal axis at the origin (thus forming a double cusp point). The first two pictures in Figure 8 show the three interfaces obtained at different times. The third one depicts the evolution of the Willmore energy associated to each evolving interface. It is interesting to compare the second and the third parameterization. The energy of the second parameterization (i.e. passing from two circles to the eight-type curve after contact) decreases smoothly, therefore the parameterization seems to relate naturally to a continuous evolution. In contrast, the energy of the third parameterization explodes at contact, and then decreases strongly to become the lowest (after time t > 10 -5 ) with respect to the other parameterizations. This experiment illustrates clearly the bifurcation at contact, and justifies why different configurations have been observed in the literature. For instance, there is no crossing observed in [START_REF] Franken | A phase field based PDE constraint optimization approach to time discrete Willmore flow[END_REF] because the authors used a large time step δ t and therefore ignored the contact. However, after a while, the energy of the non-crossing configuration is indeed the best. It may be argued that a continuous flow should go up to contact, and therefore a numerical flow that is accurate enough to capture the singularity should be the best. On the other hand, once the crossing configuration has been chosen, there is no way to have an energy as low as the non crossing configuration's energy.

Experiments in space dimension 3 For the 3D simulations presented hereafter, we used the parameters : P = 2 7 , ε = 1.5/P and δ t = 1/10 P -2 ε 2 .

The first simulation illustrates the evolution of a torus. According to the Willmore conjecture, which seems to have been proved in [START_REF] Marques | Min-Max theory and the Willmore conjecture[END_REF], the torus that minimizes the Willmore energy is Clifford's, whose ratio between both radii equals √ 2. In the first line of Figure 9, we plot for different values of t, a Clifford torus (in blue) and the interface Γ(t) (in red) obtained numerically. As expected, the interface Γ(t) converges to the Clifford torus. The second line of Figure 9 shows the evolution of a parallelepiped with two holes. The interface converges to a Lawson-Kusner surface of genus 2, that is conjectured to minimize the Willmore energy among surfaces with genus 2 [START_REF] Kusner | Comparison surfaces for the Willmore problem[END_REF][START_REF] Hsu | Minimizing the squared mean curvature integral for surfaces in space forms[END_REF]. The same experiment is done for a genus 4 surface on the last line, and there is again convergence to a Lawson-Kusner surface. We believe that these simulations illustrate the good quality of our numerical scheme and its ability to recover We present additional experiments in Figure 10 which illustrate the formation of singularities in dimension 3. On the first line, two spheres evolve by the classical diffuse flow. As the distance between the two spheres is about ε, they merge. We take in the second experiment the initial set Γ(0) as the union of two parallel cylinders. The two cylinders grow up until collision time, at which a crossing arises. The last example shows the evolution of a cube cut by a plane (more precisely, both the plane and the cube's boundary separate the two phases, as in the 2D situation of Figure 7). The cube seems to evolve to a sphere without being disturbed by the presence of the plane. All these experiments show that the classical diffuse flow may yield singularities, although the comprehension of singular solutions to the Allen-Cahn equation in dimension 3 remains uncomplete. • It is possible to simulate the crossings of more than two interfaces;

• The evolution, by the classical diffuse flow, of two interfaces after crossing seems to be similar to the evolution by a smooth parametric approach. This is in favor of a varifold interpretation of the Willmore flow.

Numerical simulations of Mugnai's flow

We now consider the PDE system associated with Mugnai's flow in the form that we introduced in Section 4. We set the approximation parameter σ = 10 -3 and we solve numerically the system using Algorithm 2.

Convergence of Mugnai's approximation

The first example illustrated in Figure 11 shows the evolution of a circle taken as initial set Γ 0 , and the comparison with the exact solution. The numerical parameters are P = 2 7 , ε = 2/P or 3//P, and δ t = 1/2ε 2 1/P 2 . The smaller is ε, the closer the numerical flow is with respect to the continuous flow. This may indicate that the penalization term B σ (u) does not influence the evolution of smooth interfaces. This is also illustrated in Figure 12 where an initial torus in 3D evolves to the Clifford torus. We present two experiments in Figure 13 obtained with the set of parameters P = 2 7 , ε = 2/P and δ t = 1/8ε 2 P -2 . The simulations indicate that the additional penalization term B σ (u) prevents the interfaces from colliding. This is coherent with what we argued in Section 2.5, i.e. that Mugnai's energy equals the classical energy plus a functional that penalizes non profile functions.

The same observation is illustrated in 3D on Figure 14. Both cylinders grow up, but deform themselves rather than colliding.

Conclusion To conclude this experimental section on Mugnai's flow, let us observe that

As long as the interfaces are smooth, Mugnai's and the classical flow behave in the same way, which was of course expected from the theoretical properties of the associated functionals. In particular, the penalization term B σ (u) has no critical influence on the evolution of a smooth interface, as long as the evolution remains smooth as well with the classical flow.

• Since Mugnai's energy W Mu ε Γ-converges in dimension 2 to the relaxation of the Willmore energy, the associated flow prevents from crossing, which is confirmed by the simulations. In 3D as well, our simulations indicate that no crossing should occur. This indicates that the Γ-convergence property should also be true in 3D for Mugnai's energy, which is so far an open question that requires a better understanding of the diffuse approximation of the genus (having in mind the Gauss-Bonnet Theorem).

3. 3 . 1 Proposition 3 . 1

 3131 Derivation of the L 2 -gradient flow of W Be ε (u) The L 2 -gradient flow of Bellettini's model is equivalent to

L 2 -Proposition 3 . 9

 239 gradient flow of W EsRäRö ε The L 2 -gradient flow of Esedoglu-Rätz-Röger's 's model is equivalent to

Figure 2 :

 2 Figure 2: Left : sampling of Γ(t) at different times t ; Right : The graphs of t → R ε (t) for ε = 2 P and ε = 3 P , compared with the exact solution.

Figure 3 :

 3 Figure 3: Evolution by the classical approximation of the Willmore flow of two disjoint circles, for various values of ε. First line: ε = 5/P ; second line: ε = 3/P ; third line: ε = 1.5/P; The curve Γ(t) is observed at times: t = 0 (left), t = 0.0004 (middle), t = 0.0008 (right).

  in black.

Figure 4 :

 4 Figure 4: Left : two examples of saddle-shaped Allen-Cahn solutions with 4 and 8 ends. Right: Allen-Cahn solution without dihedral symmetry (thus without saddle point, the 1/2-isolevel line is shown in black).

Figure 5 :

 5 Figure 5: Two different choices of a smooth parametric initial curve Γ(0) forming either two or three circles.

Figure 6 :

 6 Figure 6: Comparison between the classical phase field flow and the parametric Willmore flow (black line) starting from the initial curves of Figure 5; Left : t = 0 ; Middle : t = 5.10 -5 ; Right : t = 5.10 -4 . Both flows yield the same numerical solution.

Figure 7 :

 7 Figure 7: Two examples where the evolution of either one or two circles is not altered by an additional separating line.

Figure 8 :

 8 Figure 8: Parametric evolution of two contiguous circles associated to three different initial parametrizations; Left : Interfaces at t = 10 -5 ; Middle : Interfaces at t = 5.10 -5 ; Right : Evolution of the Willmore energy of each curve.

Figure 9 :

 9 Figure 9: Smooth evolution of Γ(t) by the classical diffuse Willmore flow in 3 D. First line: convergence to a Clifford torus (in blue). Second and third lines: convergence to Lawson-Kusner's surfaces of genus 2 and 4, respectively.

Figure 10 :

 10 Figure 10: 3D-examples of evolutions by the classical diffuse flow yielding singularities.

1 . 2 :

 12 ∂ t u = ∆µ -1 ε 2 W (u)µ + B σ (u) µ = 1 ε 2 W (u) -∆u. where B σ (u) = W (u) |∇ν u,σ | 2 -|div ν u,σ | 2curl (curl (ν u,σ )) • ν u,σwith ν u,σ = ∇u √ |∇u| 2 +σ . The initial conditions u(x, 0) and µ(x, 0) have the form   u(x, 0) = γ d(Γ 0 ) ε µ(x, 0) = -1 ε ∆d(Γ 0 )γ d(Γ 0 ) ε

Figure 11 :

 11 Figure 11: Left: evolution of a circle by Mugnai's flow; Right : Graphs of t → R ε (t) for different values of ε

Figure 12 :

 12 Figure 12: Smooth evolution by Mugnai's flow of a torus in 3D. The blue torus is the target Clifford torus.

Figure 13 :

 13 Figure 13: Illustrations in 2D that Mugnai's flow prevents colliding.

Figure 14 :

 14 Figure 14: Illustrations in 3D that Mugnai's flow prevents from colliding. The interfaces preferably deform themselves rather than merging.

Willmore functional 2.4.1 Bellettini's approximation in dimension N ≥ 2

  (x, ∇d E , ∇ 2 d E ))dH N-1 where E is smooth and d E is the signed distance function from ∂E. Such functionals include the Willmore energy since, on ∂E, H = (∆d E )∇d E = tr(∇ 2 d E )∇d E thus |H| 2 = | tr(∇ 2 d E )| 2 . Particularizing Bellettini's approximation model to this case yields the smooth functionals

	In [8], G. Bellettini proposed a diffuse model for approximating the relaxations of geometric functionals
	of the form	∂E (1 + f

Acknowledgements

The authors thank Luca Mugnai, Selim Esedoglu, Petru Mironescu, and Giovanni Bellettini for fruitful discussions.