A new energy-gap cost functional approach for the exterior Bernoulli free boundary problem

The solution to a free boundary problem of Bernoulli type, also known as Alt-Caffarelli problem, is studied via shape optimization techniques. In particular, a novel energy-gap cost functional approach with a state constraint consisting of a Robin condition is proposed as a shape optimization reformulation of the problem. Accordingly, the shape derivative of the cost is explicitly determined, and using the gradient information, a Lagrangian-like method is used to formulate an efficient boundary variation algorithm to numerically solve the minimization problem. The second order shape derivative of the cost is also computed, and through its characterization at the solution of the Bernoulli problem, the ill-posedness of the shape optimization formulation is proved. The analysis of the proposed formulation is completed by addressing the existence of optimal solution of the shape optimization problem and is accomplished by proving the continuity of the solution of the state problems with respect to the domain. The feasibility of the newly proposed method and its comparison with the classical energy-gap type cost functional approach is then presented through various numerical results. The numerical exploration issued in the study also includes results from a second-order optimization procedure based on a Newton-type method for resolving such minimization problem. This computational scheme put forward in the paper utilizes the Hessian information at the optimal solution and thus offers a state-of-the-art numerical approach for solving such free boundary problem via shape optimization setting.


1.
Introduction. In this work, we are interested in the so-called exterior Bernoulli free boundary problem which is described as follows: given a bounded and connected domain A ⊂ R 2 with a fixed boundary Γ := ∂A and a constant λ < 0, find a bounded connected domain B ⊂ R 2 with a free boundary Σ := ∂B, B ⊃Ā, and an associated state function u := u(Ω), where Ω = B \Ā, such that the following conditions are satisfied − ∆u = 0 in Ω, u = 1 on Γ, u = 0 and ∂ n u = λ on Σ, where ∂ n u := ∇u · n denotes the normal derivative of u and n is the outward unit normal vector to Σ. The problem under consideration models various physical phenomena and is also known in the literature as the Alt-Caffarelli problem (see [1]). For a discussion of the physical background and an exhaustive bibliography regarding the problem we refer the reader to [23,24].
The presence of overdetermined conditions on Σ makes the problem ill-posed and difficult to solve. Nevertheless, the method known as shape optimization (see, e.g., [18,30,45]) is already an established tool to solve the free boundary problem. The point of departure of this approach is the observation that the original problem can be written as an optimization problem of the form min Ω J(Ω, u Ω ) subject to e(u Ω ) = 0, where J denotes a suitable objective functional that depends on a domain Ω as well as on a function u Ω , which is the solution of a partial differential equation e(u) = 0 posed on Ω.
There are different ways to reformulate (1) into a shape optimization setting and one possibility is to consider the minimization problem min ΩJ (Ω) = min where the state functions u D and u N are the solutions of the following PDE systems − ∆u D = 0 in Ω, u D = 1 on Γ, u D = 0 on Σ; − ∆u N = 0 in Ω, u N = 1 on Γ, respectively. Such shape optimization reformulation of (1) has already been studied and exploited in the literature (see, e.g., [9,11,19]). We mention here that the energy-gap cost functionalJ is often attributed to Kohn and Vogelius [33] since they were among the first to use it in the context of inverse problems. The minimization problem (2) can be carried out numerically using different strategies [41]. One particular choice is the application of a gradient type method for which gradient information is needed. Accordingly, sensitivity analysis is indispensable in any optimization problem such as (2). In a recent paper, Eppler and Harbrecht [19] computed the first-and second-order shape derivatives (shape gradient and shape Hessian, respectively) ofJ through formal differentiation (see [18,30,45]). These authors also proved algebraic ill-posedness of the Kohn-Vogelius formulation (2)-(4) through the investigation of sufficient second-order conditions. The technique used by Eppler and Harbrecht in [19], which was inspired by [21,22], however, restricts their findings to starlike domains. In [11], another method using only the Eulerian derivatives [18] of the states was used to characterize the shape gradient ofJ. By this approach, Ben Abda et al. were able to extend the results in [19] to more general C 1,α domains. It was pointed in [11] that an alternative technique in computing the shape gradient ofJ using the concept of shape derivatives of u D and u N combined with the calculus for domain integral cost functionals could be given. This suggested approach was delivered by Bacani in [9] wherein he also demonstrated how to apply yet another technique called rearrangement method, first used in [31], to establish the expression for the shape gradient ofJ.
Similar to [11,19], our present investigation is primarily devoted to the numerical solution of (1) by means of shape optimization methods. Our motivation roots from the fact that in most, if not all, shape optimization formulation of the Bernoulli problem, one is required to solve at least two systems of partial differential equations (PDEs) in order to evaluate the gradient of the associated cost function (see, for instance, [10,Thm. 33] or [11,Thm. 2] for the case ofJ, and also [19,Thm. 1] for the n-dimensional case of the problem). In this work, we show that the number of state constraints to be solved when considering an energy-gap cost functional approach to (1) (such as that of (2)) can be reduced. To accomplish such objective, we replace the pure Dirichlet problem (3) with a different state system consisting of a boundary condition of the third kind. More precisely, we let β > 0 be a fixed number, and consider the equivalent form of (1) with a Robin boundary condition: as one of the state system. With this new state constraint, the shape optimization problem now reads as min Ω J(Ω) = min where the state functions u N and u R are the only solutions to (4) and (5), respectively. Contrary to the classical Kohn-Vogelius formulation (2)-(4), we shall see in Section 3 that, with the right choice of β and an additional assumption on u R , the shape gradient of the cost function obtained from this new formulation depends only on the solution of the state system (4) (see Proposition 1). As a result, the proposed formulation is more attractive compared to that of the classical setting, especially in terms of numerical computations. Here, we numerically solve the minimization problem (6) using a Lagrangian-like method as opposed to [11] which applies an Eulerian type approach to solve the Bernoulli problem. We announce in advance that, in most of the numerical examples considered here, the proposed formulation provides less computing time per iteration than the classical Kohn-Vogelius formulation. The two formulations, however, are comparable in terms of mean over-all computing time. In addition, it seems that the proposed formulation provides a somewhat more accurate approximation of the optimal solution than the classical Kohn-Vogelius approach. We support these claims with various numerical examples that are reported in subsections 5.2.2 and 5.2.3. It is worth mentioning that, to the best of our knowledge, the shape optimization formulation (6), with state systems (4) and (5), of (1) has not yet been used in any previous investigations. Therefore, the proposed formulation is novel to our work. In addition, a numerical realization of the shape optimization problem (2)-(4) via a Lagrangian method has not yet been provided in any foregoing investigation. This fact further warrants the need to carry out a numerical exploration of the classical Kohn-Vogelius formulation via a different numerical technique to [11].
The rest of the paper is organized as follows. In Section 2, we describe the weak formulations of the state equations and briefly discuss the regularity of their solutions. In Section 3, we recall a few basic concepts from shape calculus which are used in subsequent parts of the paper. Then, we examine the sensitivity of the cost function with respect to the domain. We compute the shape gradient of the cost through chain rule approach and exhibit the explicit expression for the shape Hessian using Hadamard's domain differentiation formula. Also, we characterize the shape Hessian at the solution of the Bernoulli problem and prove the instability of the optimization problem. We show the latter result by proving that strict coercivity of the shape Hessian at the solution of the Bernoulli problem cannot be established in the energy space associated with the shape Hessian. In Section 4, we address the issue of existence of optimal solution of the proposed shape optimization formulation. We first introduce a suitable family of admissible domains for the shape optimization problem and then show the existence of optimal solution by proving in particular, the continuity of the state problems (4) and (5) with respect to the domain. We then complete the proof of existence of optimal solution by establishing the lower semi-continuity of the cost function J. In Section 5, we describe how the computed shape gradient can be utilized in a gradient-based scheme to obtain numerical solutions of the present problem and outline the algorithm to be implemented. In addition, we demonstrate how to use the Hessian information at the optimal solution to formulate a second-order optimization procedure for resolving the minimization problem. Afterwards, we present several numerical tests illustrating the feasibility of the newly proposed shape optimization formulation of (1). We also compare our results with the ones obtained via the classical Kohn-Vogelius formulation. Finally, we end the paper with a conclusion in Section 6.

2.
Preliminaries. Here let us review an essential quality of the state solutions which is vital in guaranteeing the existence of shape derivatives.
2.1. Weak formulation of the state equations. In this section, we state the variational formulations of the state equations (4) and (5) and briefly discuss the solvability of their solutions.
Consider the Hilbert space V 0 (Ω) , and the linear manifold defined by Note that · V0(Ω) is basically the H 1 (Ω)-seminorm. For the weak formulation of (4) and (5), we consider two fixed functions u N0 , u R0 ∈ H 1 (U ) such that u N0 = u R0 = 1 on Γ. Now, the weak formulation of (4) can be stated as follows: Similarly, the weak form of (5) can be expressed as follows: The variational equation (7) is known to admit a unique solution in H 1 (Ω), while it can easily be verified (for instance, by means of Lax-Milgram lemma) that (8) also have a unique solution in H 1 (Ω) (see [26]).

2.2.
Higher regularity of the states. The unique solution u N of the PDE system (4) actually possesses higher regularity property because of the regularity assumption on Ω. In fact, since Ω is of class C 2,1 , the weak solution u N ∈ H 1 (Ω) to (4) is also in H 3 (Ω), and in general, if Ω is of class C k+1,1 , where k is a non-negative integer, then u N is H k+2 regular. This assertion can easily be established since dist(Γ, Σ) > 0, as demonstrated in [10]. Similarly, the unique solution u R ∈ H 1 (Ω) of (5) also have higher regularity depending on the degree of smoothness of Ω. More precisely, if Ω is of class C k+1,1 (again k is a non-negative integer), then u R is also an element of H k+2 (Ω) (see, e.g., [32,Rem. 3.5]).
3. Shape sensitivity analysis of the states and cost function. Let us consider a fixed, bounded and connected C 2,1 domain U ⊃Ω and a family of deformation fields Θ := {V ∈ C 2,1 (U , R 2 ) : V = 0 on ∂U ∪ Γ}. (9) Given an element of Θ, we perturb Ω by means of the so-called perturbation of the identity operator (see, e.g., [18,Sec. 2.5.2,p. 147] or [10]): For sufficiently small t and for each V ∈ Θ, one can show that T t is a C 2,1 diffeomorphism from Ω onto its image (cf. [45]). Note that, in view of (9), Γ remains invariant after a deformation since V vanishes on Γ. Hence, Γ is a component of the boundary of Ω t for all t. Further, we note Ω 0 = Ω and Γ 0 = Γ. On the perturbed domain Ω t , the state solutions u Nt and u Rt satisfy respectively, where n t is the unit outward normal to Σ t . Here, we can actually drop t in Γ t since Γ t = Γ for all t.
In what follows, we recall some definitions from shape calculus. We say that the function u(Ω) has a material derivativeu and a shape derivative u at zero in the direction V if the limitṡ ). These expressions are related by u =u − (∇u · V) provided that ∇u · V exists in some appropriate function space [18,45]. In general, ifu and ∇u · V both exist in the Sobolev space W m,p (Ω), then u also exists in that space.
Given a functional J : Ω → R, we say that it has a directional Eulerian derivative at Ω in the direction V if the limit exists. In addition, if dJ(Ω)[V] exists for all V and the map V → dJ(Ω)[V] is linear and continuous, then J is shape differentiable at Ω, and this mapping will be referred to as the shape gradient of J at Ω.
3.1. Shape derivative of the states. In order to prove the existence of dJ(Ω)[V], one needs to show that the material and shape derivatives of the states u N and u R exist and, consequently, apply the chain rule. Deriving the expressions for these aforementioned quantities are quite laborious. Fortunately, the system of PDEs which the shape derivatives of u N and u R satisfy were already established in [9] and [46], respectively. To guarantee their existence, however, we require Ω to be at least C 2,1 regular.
Let Ω be a bounded C 2,1 domain and V ∈ Θ. Then, u N ∈ H 3 (Ω) is shape differentiable with respect to the domain, and its shape derivative where where κ denotes the mean curvature of Σ.

Lemma 3.2 ([46]).
Let Ω be a bounded C 2,1 domain and V ∈ Θ. Then, u R ∈ H 3 (Ω) is shape differentiable with respect to the domain, and its shape derivative where B(·)[V] is given by (13). If β = κ, then for the shape derivative u R of the solution of (5), it holds that u R ≡ 0 when Σ is the free boundary. Remark 1. The last statement in the previous lemma can easily be verified as follows. Indeed, from (14), we note that the shape derivative u R satisfies the equation Hence, if Σ is the free boundary, then u R = 0 on Σ, and so −∇ Σ u R · ∇ Σ ϕ = 0. This leaves us the equation Since ∂ n u R = λ on Σ, then by choosing β = κ, the right side of the above equation vanishes and we get u R ≡ 0. We emphasize that this result plays an important part in this investigation, particularly in simplifying the expression for the shape derivative of J which would make the minimization problem easier to solve.
3.2. First-order shape derivative of the cost function. Now, our objective here is to derive the shape derivative of the cost function J in the direction of the deformation field V ∈ Θ. For this purpose, we utilize Hadamard's differentiation formula (cf. [18,Thm. 4.3,p. 486], see also [30,45]) This formula holds when f ∈ C([0, ε], W 1,1 (U )) and d dt f (0) exists in L 1 (U ). In addition to (15), we shall also use the so-called tangential Green's formula (cf. [18,Eq. 5.27,p. 498]): let U be a bounded domain of class C 1,1 and Ω ⊂ U with boundary Γ. For V ∈ C 1,1 (U , R 2 ) and f ∈ W 2,1 (U ), we have where κ is the mean curvature of Γ and the tangential gradient ∇ Γ is given by We are now ready to prove the following result.

Proposition 1.
Let Ω be of class C 2,1 and V ∈ Θ. Then, the energy-gap cost functional J is shape differentiable with where Here, κ denotes the mean curvature of Σ. Moreover, τ represents the unit tangent vector on Σ and it is oriented in such a way that Σ is at the left of τ ; that is, if n = (n 1 , n 2 ) , then τ = (−n 2 , n 1 ) .
If, in addition, u R is the shape derivative of the solution of (5) where Σ is the free boundary, then, for β = κ, the shape gradient considerably simplifies to where Proof. Let Ω be of class C 2,1 and V ∈ Θ. Since u N and u R are sufficiently regular, we can apply (15) to obtain (21) where u N and u R satisfy (12) and (14), respectively. The second integral can be easily expanded as On the other hand, we can write the first integral as follows Equality 1 was obtained through Green's formula. Meanwhile, 2 follows from (12) and (14). Lastly, equality 3 was derived through the identity which holds since V · n∇ Σ ϕ · n = 0 (see (16)). Further, in 3 , we may write . Hence, combining the computed expressions for the first and second integral of (21), we get (17). Letting β = κ and u R be the shape derivative of the solution of (5) where Σ is the free boundary, we get u R = 0 and ∂ n u R = λ on Σ. By virtue of [46, Lem. 1], u R ≡ 0. In addition, we obtain the relation ∇u R = (∂ n u R )n on Σ from which we deduce that ∇u R · τ = 0 on Σ. Direct substitutions of these identities in (17) eventually lead to (19), completing the proof of the proposition.
As an immediate consequence of Proposition 1, we have the following. Corollary 1. Let the domain Ω * be such that u = u(Ω * ) satisfies the overdetermined boundary value problem (1); i.e., it holds that Then, the domain Ω * fulfils the necessary optimality condition Proof. At the shape solution Ω = Ω * of the Bernoulli problem (1), u N = 0 on Σ. Hence, ∇u N = (∂ n u N )n on Σ and it follows that ∇u N · τ = 0 on Σ. Thus, G is zero which implies the assertion.
Remark 3. We reiterate that, as opposed to the classical Kohn-Vogelius formulation (2)-(4) of the Bernoulli problem, one only needs to solve (4) in order to evaluate the shape gradient of J. We recall from [10,Thm. 33] (see also [11,Thm. 2] and [19,Thm. 1]) that under the classical setting, the first-order shape derivative of the Kohn-Vogelius cost functional in (2) is given by and the state functions u D and u N are the unique solutions to (3) and (4), respectively. Apparently, the shape gradient of J is obtained from (22) when ∂ n u D = λ. Of course, this condition holds in the case of (5) if one assumes that β = κ and u R is the shape derivative of the solution of (5) where Σ is the free boundary (cf. [46]).
3.3. Second-order shape derivative of the cost function. Let W be a velocity field in Θ. By virtue of Proposition 1 together with Remark 2, the derivative dJ(Ω s (W))[V] exists for all sufficiently small s. Our next objective now is to find an expression for the limit Here, Σ s := Σ s (W) denotes the free boundary of the deformed domain Ω s := Ω s (W) perturbed via the velocity field W ∈ Θ and u Ns ∈ H 3 (Ω s ) is the unique solution of the state system (4) onΩ =Ω s (cf. equation (10) with t replaced by s). In addition, κ s = div Σs n s ; and n s and τ s refer to the unit outward normal and unit tangent vectors on Σ s , respectively. If for all V and W, d 2 J(Ω)[V, W] exists and is bilinear and continuous with respect to V and W, then J is said to be twice shape differentiable at Ω. In this case, the map (V, W) → d 2 J(Ω)[V, W] is called the shape Hessian of J in the directions of V and W.
Before we give the characterization of the shape Hessian below, we first introduce some notations. For convenience, we use the notation V n := V · n for V ∈ Θ. Moreover, we let v = V| Σ and v = v Σ + v n n := (v · τ )τ + (v · n)n. Also, we recall from [18, Eq. (5.2), p. 495] the definition of the tangential differential operator D Σ called the tangential Jacobian matrix given as Furthermore, we denote by ϕ W the shape derivative of ϕ along a deformation field W ∈ Θ. Lastly, we recall the shape derivatives of the mean curvature κ and the tangential vector τ in the direction of W which are respectively given as (see [18,45]) We now state the shape Hessian of J at Ω in the following proposition.
Proposition 2. The second-order shape derivative d 2 J(Ω)[V, W] of the cost functional J at Ω in the directions of the deformation fields V, W ∈ Θ has the following structure: (26) where G is given by (20) and Proof. Using Stoke's formula, we write the gradient as From (15), we easily find that where G W given by (27) is simply obtained by differentiating G . Note that the shape Hessian (30) is consistent with the Hadamard-Zolésio structure theorem provided that u NW is linear and continuous function of V n on the boundary Σ (cf. [18]). Moreover, notice in equation (30) the non-symmetry of the second integral in V n and W n . Now we further extract from the second integral in (30) a few more symmetric expressions of the shape Hessian. The non-symmetric part will be obtained from (19) applied to the deformation field DVW (cf. [37]). Since Ω is a C 2,1 domain and we have enough regularity for the state u N , then in view of [18, Chap. 9, Sec. 6] (or, in particular, [18, Eq. 6.10, p. 505]), the following identity actually holds where K is given by (29). Here, the explicit form of ∂ n G is computed as follows Since ∂ n u N = λ on Σ, the last expression in (32) disappear and we immediately arrive at (28). Substituting this expression into (30) proves the expression for the shape Hessian given by (26).
We point out that we may in fact consider velocity fields V ∈ Θ that are normal to Σ, i.e., V| Σ = V n n. In this case, we observe that the expression in (29) vanishes. Moreover, the last integrand in (26) can be expressed as G (DV)W n = G ∂ n v n w n since (DV)W n n = D(v n n)w n n · n = w n n [n(∇(v n )) + v n Dn]n and Dn · n = 0 on Σ (cf. [45]). Thus, for deformation fields V, W ∈ Θ with null tangential part, the shape Hessian (26) simplifies to Note that, in the proof of Proposition 2, we managed not to use the second-order shape derivative of the states in characterizing the second-order Eulerian derivative of J contrary to [21]. Thanks to identity (31), the characterization of the shape Hessian was easily accomplished. (27) the term κ W which represents the shape derivative of the mean curvature κ along the deformation field W ∈ Θ. Note that the explicit form of κ W given in (24) essentially consists of a second order tangential derivative of the vector field W, and this derivative actually exists due to our assumption that Ω is of class C 2,1 (cf. [18,45]). Hence, we deduce that the shape Hessian defines a continuous bilinear form

Remark 4. Notice in
Let us now consider the shape solution Ω * of the Bernoulli problem (1), focusing on the case when β = κ. In order to study the stability of the shape optimization problem (6) at the solution of the Bernoulli problem, we first characterize the shape Hessian of J along the deformation fields V, W ∈ Θ at Ω * .
Corollary 2. At the solution of the Bernoulli problem (1), i.e., Ω = Ω * , we have Proof. Using Stoke's theorem, we write Then, in the direction of a velocity field W ∈ Θ, we have At the shape solution Ω = Ω * of the Bernoulli problem (1), we know that u N = 0 on Σ. Hence, the second and the third integral disappear, and we are left with For J 2 , we have Furthermore, at Ω = Ω * , we have that ∇u N · τ = 0 on Σ. Thus, dJ 2 (Ω * )[W] vanishes, and we conclude that Alternative Proof. Of course, equation (33) can be proven directly from Proposition 2 together with Corollary 1. Indeed, at Ω = Ω * , we get G = 0 and (27) reduces to G W = λκu NW . Moveover, we have ∂ n G = λ 2 κ. Inserting these identities to yield (33) as announced.

3.4.
Ill-posedness of the shape optimization problem. Here we prove the ill-posedness of the present shape optimization problem using the method already used in [20]. Particularly, we use a local regularity argument in order to prove that the shape optimization problem under consideration is not well-posed. To this end, we recall from [16,17] an important result about sufficient second order conditions: a local minimizer Ω * is stable if and only if the shape Hessian d 2 J(Ω * ) is strictly coercive in its corresponding energy space H 1 (Σ * ) (see also [19]). Throughout this section, for convenience, we use the standard notation " ", i.e., by B D we mean that there is some constant Proposition 3 (Coercivity of the Hessian at the solution of the Bernoulli problem). Let the mean curvature κ of Σ * be non-negative. Then, the shape Hessian at the solution of the Bernoulli problem (1) is strictly H 1/2 (Σ * )-positive, i.e., . To justify the above result, we first prove that multiplication by a Lipschitz function defines a bounded operator in H 1/2 (∂Ω).
Let φ be a Lipschitz function with Lipschitz constant L. Then, we get the inequality Hence, we can estimate |φv| 1/2,2,Γ as follows , then the assertion is proved.
With the aid of Lemma 3.3, the following multiplication operators are continuous: Here, of course, κ denotes the mean curvature of Σ * . We note that, in addition to being continuous, the map L : To complete our preparation for the proof of Proposition 3, let us define S as the Steklov-Poincaré operator on Σ * which is defined by (see [46]) where Ψ ∈ H 1 (Ω * ) satisfies The operator S, also called the Dirichlet-to-Neumann map (see, e.g., [19]), is . Its inverse R called the Neumann-to-Dirichlet map is defined by Proof of Proposition 3. Using R, we can write (33) as Using the continuity of the maps L and M, the bijectivity of R, and the fact that As we have just shown, the coercivity of the shape Hessian can only be attained with respect to a Sobolev space weaker than the energy space. Hence, the shape optimization problem (6), subject to (4) and (5), is algebraically ill-posed. This result also indicates that the gradient does not have a uniform sensitivity with respect to the deformation directions.

4.
Existence of optimal domains of the shape optimization problem. Before we proceed to the numerical realization of the proposed shape optimization formulation of (1), we first address the issue of existence of optimal solution to the shape optimization problem (6), subject to (4)- (5). Meanwhile, regarding the existence of solution to the exterior Bernoulli free boundary problem (1), we refer the readers to [1].
We emphasize that in (6), we are actually trying to minimize the cost functional (4) and (5), respectively. In view of their respective weak formulations (7) and (8), noting that each of these systems admits a unique solution in H 1 (Ω), we can then define the mapping Ω → (z N , z R ) = (z N (Ω), z R (Ω)) and denote its graph by So, in (6), we are actually looking for a solution (Ω, z N (Ω), z R (Ω)) that minimizes J(Ω) = J(Ω, z N (Ω), z R (Ω)) on F. This minimization problem is usually solved by endowing the set F with a topology for which F is compact and J is lower semicontinuous. In accomplishing the task, we shall follow the ideas developed in [27] and the ones used in [12,28]. To do this, we first need to define the set of admissible domains O ad and then give an appropriate topology on it. From Subsection 3.1, we recall that we used C 2,1 regularity of the free boundary Σ to guarantee the existence of the shape derivatives of the states. For the analysis of existence of optimal solution to (6), we may as well suppose that the free boundary Σ is C 2,1regular, however, as we have seen in [12], a C 1,1 regularity on the free boundary is enough to carry out the task. So, we assume that Σ is parametrized by where φ ∈ U ad and U ad is the set of vector functions φ ∈ C 1,1 (R, R 2 ) having the following properties: (P 1 ) φ is injective on (0, 1] and is 1-periodic; (P 2 ) there exist positive constants c 0 , c 1 , c 2 and c 3 such that The set U in assumption (P 3 ) and the one introduced in Section 3 are not necessarily the same set. However, we point out that in (P 3 ), we are assuming that all admissible domains Ω(φ) are contained in the hold-all domain U (in the same manner that the universal set U in equation (9) holds all the possible deformations of the reference domain Ω). Here, we are in fact requiring that dist(Σ(φ), ∂U ) > 0.
Remark 5. The function φ described above can actually be viewed as a mapping from the quotient space Remark 6. Note that, in assumption (P 2 ), we are already requiring the free boundary Σ to be represented by a closed C 1,1 -curve. In [28], a C 2 -regularity of the free boundary was used to construct a C 1 -diffeomorphism of a uniform tubular neighborhood of the free boundary. Hence, as a consequence of the definition of O ad , every admissible domain Ω(φ), φ ∈ U ad , is a uniformly open set in R 2 and therefore satisfy the uniform cone property (cf. [30]). We also mention that it is actually possible to construct a C 1 -diffeomorphism of a uniform tubular neighborhood of the boundary by using only C 1 -regularity of the boundary (see [12,Cor. 1]).
Given the definition of elements of U ad , the set of admissible domains O ad is then defined by In view of the Remark 6, we mention that the class of domains Ω(φ), φ ∈ U ad , being considered here possesses a very important extension property. In fact, for every k 1, p > 1 and domain Ω ∈ O ad , there exists an extension operator such that and C > 0 is independent of the domain Ω [15]. The uniform cone property of every admissible domain Ω(φ), together with the previously mentioned result, ensured that every function u ∈ H 1 (Ω) has a uniform extensionũ ∈ H 1 (U ) from Ω to U . With the above results now at our disposal, we now define the topology we shall work with. First, we define the convergence of a sequence {φ n } ⊂ U ad by i.e., if and only if φ n → φ in the C 1 -topology. Then, the convergence of a sequence of domains Meanwhile, the convergence of a sequence {z Nn } of solutions of (7) on Ω n to the solution of (7) on Ω can then be defined by Similarly, we define the convergence of a sequence {z Rn } of solutions of (8) on Ω n to the solution of (8) on Ω by Of course, in (44) and (45), the extensionsz i ,z in , i = N, R, are actually defined as Finally, the topology we introduce on F is the one induced by the convergence defined by We now state the main result of this section.

Proposition 4. The minimization problem
Find where z N and z R are weak solutions to (7) and (8), respectively, admits a solution in F.
As already mentioned, we established the above result by proving the compactness of F and the lower semi-continuity of J. Regarding the compactness of F with respect to the convergence (46), we note that the convergence φ n → φ easily follows from the compactness of U ad and the Arzelà-Ascoli theorem. Indeed, by the said theorem, there is a subsequence {φ k } that converges uniformly in the C 1 ([0, 1], R 2 )norm. Denoting its limit by φ, we know that |φ| c 0 . Since c 1 |φ k | c 2 a.e. in (0, 1) and we have the uniform convergence φ k → φ and φ k → φ in (0, 1), we deduce that c 1 |φ | c 2 a.e. in (0, 1). Hence, φ ∈ U ad . This implies that we only need to show the continuity of the state problems (4) and (5) with respect to the domain in order to complete the proof of compactness of F.

4.1.
Continuity of the state problems. In the next result, we state the continuity of the state problems.
Proof. The first convergence has already been verified to hold via the compactness of U ad and the Arzelà-Ascoli theorem. So, we proceed on showing the validity of the last two convergences. We only prove the convergencez Rk z R in H 1 (U ) since the proof for the other one uses the same argument. First, we prove the following claim.
Claim. There exists a uniform extensionz Rk of z Rk from Ω k to U and a constant C > 0 independent of k such that z Rk H 1 (U ) C.
By a result of D. Chenais [15], we know that a solution z Rk of (8) on Ω k admits an extensionz Rk in H 1 (U ) such that for some non-negative constantC (independent of k), the following inequality holds: Clearly, to verify our claim, we need to show that z Rk H 1 (Ω k ) is bounded. In view of (8), letting the test function be ϕ = z Rk ∈ V 0 (Ω k ), we have λz Rk dσ.
From this equality, we get the estimate where Σ k := Σ(φ k ). Next, we find an estimate for z Rk L 2 (Σ k ) in terms of z Rk V0(Ω k ) . To do this, we use a uniform Poincaré inequality proved in [14], apply the inequality (41), and utilize the uniform continuity of the trace operator with respect to the domain. The former result is precisely given by the following inequality: Note that, since every domain in O ad satisfies the uniform cone property, the above result actually follows from [14,Cor. 3(ii)]. On the other hand, concerning the trace operator, the following inequality holds for all real number q such that 1 2 < q 1, φ ∈ U ad and functions z ∈ H 1 (U ) (cf. [12,Thm. 4]): where, of course, · H q (U ) denotes the H q (U )−norm. (In fact, due to assumption (P 3 ) and the uniform cone property of the domain Ω(φ) ∈ O ad , the norm of the trace map tr : H 1 0 (U ) → L 2 (Σ(φ)) can be bounded uniformly with respect to Ω(φ) ∈ O ad ; see [35]. ) We first apply inequality (50) (with q = 1), making use of the extension of z Rk ∈ H 1 (Ω k ) to H 1 (U ), then apply (41), and finally employ (49) to obtain Hence, going back to (48), we have Thus, after some manipulations, we eventually get This proves our claim which means that the sequence { z Rk H 1 (U ) } is bounded.
Now, in what follows, we show that z R (φ) =z R | Ω(φ) is the solution of (8) on Ω(φ). In fact, we will show that the variational equation also holds for all test functions v ∈ V 0 := V 0 (U ) = {ṽ ∈ H 1 (U ) :ṽ = 0 on Γ}. Clearly, the restriction on Ω k of any element v of V 0 is in V 0 (Ω k ), for all k, which is exactly the test space of (8) on Ω k . Hence, we have Now, we prove that, by passing to the limit, we will obtain (51) from (52). To see this, we simply take the difference of equations (51) and (52) and then let k → ∞. As for the difference of the last two integrals, we have where in the last inequality, we applied assumption (P 2 ). Since v ∈ V 0 ⊂ H 1 (U ), then, according to [13, Cor. 1], the following limit actually holds for any sequence {φ n } ⊂ U ad and element φ ∈ U ad such that φ n → φ in the sense of (42). In addition, we have, in view of (50), the estimate v L 2 (Σ(φ)) K v H 1 (U ) . Finally, using the convergence φ k → φ in the C 1 ([0, 1], R 2 )-norm (cf. (42)), we get lim k→∞ I 4 = 0 as desired.
Similarly, we have Now, using the fact that z Rk =z Rk | Ω k ∈ H 1 (U ), together with the estimate (50), we deduce, via the application of [13,Cor. 1], that the first and the third summands in the last inequality above eventually vanished. Also, using (50), and the compactness of the injection of H 1 (U ) into H q (U ) for 1 2 < q < 1, the left side of the inequality also goes to zero as k → ∞. Likewise, the right side of the inequality also disappears because of (42). Hence, we also have lim k→∞ I 3 = 0. For the remaining two differences the desired limits lim k→∞ I 1 = lim k→∞ I 2 = 0 are obtained by applying the convergencez Rk z R in H 1 (U )-weak and the convergence of characteristic functions (see, e.g., [ together with the fact that the sequence { z Rk H 1 (U ) } is bounded. This proves that z R (φ) =z R | Ω(φ) is the solution of (8) on Ω(φ). Applying the same arguments used above, we can also show that there exists a sequence of uniform extensions {z Nk } of {z Nk } which is uniformly bounded in H 1 (U ), i.e., { z Nk H 1 (U ) } is bounded. Utilizing this result, we can also prove (again, following the lines of arguments used above) that z N (φ) =z N | Ω(φ) is in fact the solution of (7) on Ω(φ), finally completing the proof of Proposition 5.

4.2.
Lower semi-continuity of the cost function J. To complete the proof of Proposition 4, let us now establish the lower semi-continuity of J in the following result.
Proposition 6. The cost functional is lower semi-continuous on F in the topology induced by (46).
Proof. Let {(Ω n , u Nn , u Rn )} be a sequence in F, Ω n = Ω(φ n ), and assume that where Ω = Ω(φ) and (Ω, u N , u R ) ∈ F. For convenience, we let w n = u Rn − u Nn and w = u R − u N . Moreover, we denote their respective extensions byw n ,w ∈ H 1 (U ). In proving that J(Ω n , u Nn , u Rn ) → J(Ω, u N , u R ), we use the identity Using (53) and the fact that {w n } is bounded, we get lim n→∞ J 1 = 0. Meanwhile, we can write J 2 as Hence, we also have lim n→∞ J 2 = 0 becausew n w in H 1 (U )-weak. Therefore, lim n→∞ J = 0.
To end this section, let us formally provide the proof of Proposition 4 using the main results established in the last two subsections.
We apply Proposition 5 to obtain a subsequence (Ω k , z Nk , z Rk ) and an element Ω = Ω(φ) ∈ O ad such that Ω k → Ω (i.e., φ k → φ uniformly in the C 1 topology), z Nk z N ,z Rk z R in H 1 (U ), and the functionsz N | Ω andz R | Ω are the solutions to the variational equations (7) and (8) in Ω, respectively. Using these, together with the continuity of J proved in Proposition 6, we conclude that (by virtue of [27, Thm. 2.10])

Numerical algorithm and examples.
To numerically solve the minimization problem (6), we propose to apply a gradient based scheme by means of a Lagrangianlike method as oppose to [11] which employs an Eulerian-like type method known as level-set method (see [39]). Alternatively, one could also apply a variant of Newton's method which, in addition to the shape gradient, also requires the knowledge of the shape Hessian. This approach, however, is much more difficult to utilize and numerically implement (see, e.g., [38,44], and the references therein). Nevertheless, we shall present in Subsection 5.1.4 a second-order method based on a modified H 1 -Newton method (see [2]).

Numerical algorithm.
5.1.1. The Sobolev gradient method. Let us denote by Ω k the shape of the domain at the k th iteration. Then, at the (k + 1) th iteration, the shape Ω could be updated as Ω k+1 := Ω t k+1 = (I 2 + t k V)Ω, where t k 0 is some small step size parameter. In perturbing the domain Ω, we may take as the descent direction. This choice of the descent direction, however, may cause undesirable oscillations on the boundary of the approximate shape solution. To avoid such phenomena, we define the descent direction based from a variant of what we call Sobolev gradient method (cf. [4,36]); that is, we compute V as the unique solution in [V 0 (Ω)] 2 of the variational problem Here, of course, κ denotes the mean curvature of the boundary Σ. In this sense, the vector field V provides a smooth extension of G n over the entire domain Ω which not only smoothes the boundary [5] but also preconditions the descent direction. In (55), the boundary integral with the term κ can be viewed as a perimeter regularization, weighted with weights V n . Computing the descent direction V using (55) is actually inspired by the so-called H 1 gradient method [4] which, on the other hand, was based on the idea of the traction method [5,6,7,8]. Notice in equations (20) and (55) that, in addition to solving the state equations (4) and (5), we also need to evaluate the mean curvature κ of Σ to compute the descent direction V. We recall from [30,Prop. 5.4.8,p. 218] (see also [25,Lem. 16.1, p. 390]) that, for a domain Ω of C 2 class, the mean curvature can be defined as κ = div Σ n = div N, where N is any (unitary) extension of n that is of class C 1 . Having this idea in mind, we calculate κ by evaluating the expression div N, where N is the unique element in [H 1 (Ω)] 2 of the variational equation

5.1.2.
Step Size. The choice of the step size parameter t k is not a straightforward. Too large, the algorithm is unstable; too small, the rate of convergence is insignificant. In updating t k , one could follow a heuristic approach inspired by the Armijo-Goldstein line search strategy similar to [11] or use a back-tracking procedure discussed, for instance, in [40]. In this work, we follow the former procedure with a slight modification on the main formula for the step size. According to [11], the step size t k ∈ (0, ε], where ε > 0 is some small real number, can be updated as follows. Using the definition of Ω ε and V| Σ = −G n, we have . The requirement J(Ω ε ) = (1 − α)J(Ω 0 ) for some α ∈ (0, 1) then suggests the choice ε = αJ(Ω 0 )/ G 2 L 2 (Σ0) . Hence, at each iteration, we may choose for a fixed α, as the step size. However, in this investigation, we take into account the fact that the deformation field V is computed through equation (55). So, we replace the denominator G 2 L 2 (Σ k ) by V 2 L 2 (Σ k ) (see [41]) and take, for a fixed α ∈ (0, 1), the step size parameter t k as The above step size value is chosen whenever J(Ω k+1 ) J(Ω k ). Otherwise, we reduce the step size and go backward: the next iteration is initialized with the previous shape Ω k . The step t k is also decreased if reversed triangles are detected within the mesh update.

The Boundary Variation
Algorithm. Now we summarize the main steps required for the computation of the k th domain as follows: Step 1: Set the parameters α and η. Also, choose an initial shape Ω 0 .
Step 2: Compute the solutions u N of the state problem (4) on Ω k and evaluate the mean curvature κ = div N, where N is obtained from (56).
Step 3: Compute the descent direction V k using (55).
Step 4: Using V k and the current step size t k , perturb the current domain by Finally, to complete the above steps, we need to specify the stopping condition. A classical stopping criterion is to find that whether the shape gradients in some suitable norm are small enough, but here we terminate the algorithm as soon as the inequality condition |J(Ω k+1 ) − J(Ω k )| < η (58) is satisfied for some sufficiently small real number η > 0.

5.1.4.
Second order optimization method. We remark that, with the help of the shape Hessian information, a regularized Newton method could be used as a numerical procedure to solve the minimization problem (6) (see, e.g., [20]). However, we emphasize that our main purpose in this investigation is to reduce the number of associated PDE systems to be solved during the optimization procedure. Applying a second-order method will obviously lessen the number of iterations needed to reduce the cost at certain magnitude. The disadvantage, however, is the additional computational burden and time to carry out the task. Moreover, note that u NW depends on the velocity field W. Hence, applying a second-order method would require the solution of a system of PDEs for each velocity field W. This, in turn, will make the computation of the descent direction more involved. Nevertheless, the said issue can be resolved by introducing the adjoint method. Here, for the sake of comparison, we shall also formulate a second-order optimization algorithm to solve the minimization problem (6). To do this, we follow the idea about a second-order method for solving shape optimization problems proposed by Azegami in [3] (see also [2]). Particularly, we use a variant of the so-called H 1 Newton method which utilizes the Hessian information to compute the descent direction.
For our purpose, we use the Hessian information at the solution of (1). Using (12), we introduce the adjoint variable p N ∈ H 1 (Ω) which is the only solution to the PDE system Hence, we may write the shape Hessian d 2 J(Ω) [V, W] at Ω = Ω * as Having this expression at our disposal, we define the descent direction W ∈ [V 0 (Ω)] 2 as the unique solution of the variational equation Now, the main steps to compute the k th domain Ω k are basically the same as that given in Section 5.1.3. However, in order to take into account the procedure in computing W, we divide Step 3 of the original algorithm as follows: Step 3.1: Compute the descent direction V k using (55).
Step 3.2: Compute p N by solving the system (59) at Ω = Ω k .
Step 3.3: Compute the descent direction W k using (60).
Of course, in Step 4, we must replace V k with the new deformation field W k ; that is, we perturb the k th domain by Ω k+1 = (I 2 + t k W k )Ω k . Here, however, the step size t k is still chosen on the basis of the formula given by (57). We emphasize that the adjoint variable p N satisfying (59) is essentially the shape derivative u NV at Ω * .
In the next section we illustrate the feasibility of the above methods in solving the minimization problem (6). We first test the accuracy of the proposed method in Subsection 5.2.1 by examining a particular test case with known analytical solution. Also, we show in the same subsection that the proposed second-order method could improve the number of iterations required for solving the optimization process. Then, in Subsection 5.2.2, we consider an L-shaped fixed boundary with different values for λ and compare the results of the proposed method with that of the classical Kohn-Vogelius formulation. In Subsection 5.2.3, we feature a numerical example wherein the fixed boundary is a union of two disjoint curves. We mention here that by providing numerical results using the classical formulation using the boundary variation algorithm presented in Subsection 5.1.3, we are able to exhibit an alternative numerical scheme via a Lagrangian formulation to [11] which considers a level-set approach in numerically solving the optimization problem (2)-(4).

Numerical examples.
The numerical simulations exhibited here are performed in two-dimension using the programming software FreeFem++ (see [29]). All systems of partial differential equations are solved using P2 finite element discretization where the number of discretization points on the free and fixed boundaries are initially set to N ext and N int , respectively. In deforming the shape of the domain during the optimization procedure, we utilize the function movemesh of FreeFem++ and use the function adaptmesh to refine and avoid the degeneracy of the triangles in the meshes with a maximum edge size h max . In all test cases, the exterior and interior boundaries are discretized with N ext × N int = 120 × 100 discretization points. All computations are carried out on a 1.6 GHz Intel Core i5 Macintosh computer with 4GB RAM processors. 5.2.1. Example 1: Accuracy of the computed gradient. We begin by testing the accuracy of the computed gradient. To do this, we consider the exterior Bernoulli problem with where C(0, r) is the circle centered at the origin with radius r. In this case, the only solution is the circle C(0, R). We let r = 0.3 and R = 0.5. These give us λ = −3.9152. We consider three different initial guesses defined as follows (see Figure 1a (57) particular, the final values of the cost, the Hausdorff distances between the computed optimal shape and the exact shape, the average distancesR of points on the computed free boundaries to the origin, the relative errors betweenR and the exact radius of the free boundary R, the number of iterations until termination of the algorithm and over-all computing times. Clearly, as α increases in value, the number of iterations, as well as the computing time, decreases. Moreover, except for the case when the initial guess is Σ 3 0 and α = 0.1, the Hausdorff distance between the exact optimal shape Σ * = C(0, 0.5) and each of the computed optimal (final) shape Σ i f , i = 1, 2, 3, is (approx.) equal to 0.005. Meanwhile, we notice large number of iterations required to reach convergence when α is set to 0.1. These values can obviously be reduced by applying a second-order method. Indeed, by employing the modified H 1 -Newton method presented in subsection 5.1.4, we obtain a significant reduction in the number of iterations needed to reach convergence as evident in Table 2. The computed optimal free boundaries when α = 0.1, as well as the fixed    Figure 2. Looking at Figure 2a and Figure 2c, we observe that the first and second tests, where we respectively took Σ 1 0 and Σ 2 0 as initial guesses, have almost the same rate of convergence and are both faster compared to when taking Σ 3 0 as the initial profile for the free boundary Σ. However, in terms of convergence to the exact solution (measuring the Hausdorff distance between the k th approximation Σ k of the free boundary and its exact profile Σ * = C(0, 0.5)), the choice Σ 2 0 gives the fastest rate of convergence among the three choices for Σ 0 (with the third choice Σ 3 0 as the initial guess for the free boundary giving the slowest convergence rate to the exact solution) in case of using the first-order method is applied (see Figure 2b). On the other hand, it appears that the first and second test cases have almost the same convergence rate when using the second-order method (refer to Figure 2d). In testing the accuracy of the computed gradient, we also considered coarser mesh in solving the PDE systems involved in the formulations. It seems that the coarser the mesh is, the less computing time is needed to complete the optimization process (as expected). However, we obtained a more accurate final free boundary in all test cases when using finer mesh during the discretization process (also as expected).

Example 2:
An L-shaped fixed domain. Next, we consider the boundary Γ = ∂S of the L-shaped domain S = (−0.25, 0.25) 2 \ [0, 0.25] 2 and compute the optimal shape for all integers λ = −10, −9, . . . , −1. In all situations, we take α = 0.99 and choose C(0, 0.6) as the initial shape of the free boundary Σ. We compare our results with the ones obtained using the classical Kohn-Vogelius formulation (2)- (4). As in the previous example, we also set η = 10 −6 in the stopping condition (58) and set h max = 0.01. Table 3 summarizes the computational results for the present test  Table 3. Summary of computational results for an L-shaped fixed boundary Γ = ∂S with λ = −10, −9, . . . , −1 where α = 0.99 in (57) and η = 10 −6 in the stopping condition (58) cases obtained through the proposed shape optimization formulation (6), subject to (4) and (5), versus the classical Kohn-Vogelius formulation (2)-(4). The table shows, in particular, the initial step sizes, the final cost values, the number of iterations until termination of the algorithm and over-all computing times for each of the two formulations. Notice that when λ = −9, the computing time for the classical formulation is too large compared to other cases. Also, we observe that only in cases when λ = −4, −1 that the number of iterations of the proposed formulation is higher compared to the classical formulation. The resulting exterior boundaries are shown in Figure 3 (blue-colored lines) where the outermost boundary corresponds to λ = −1 and the innermost boundary to λ = −10. The fixed boundary Γ and the initial shape of the free boundary Σ are also depicted in the figure, and are respectively colored with black and magenta colors. We compared the computed optimal free boundaries from the two formulations and we noticed that, in all cases being considered, the results are indistinguishable from each other (see, e.g., Figure 4a for λ = −9). Also, we observed that the histories of cost values (as well as the L 2 -norms of the descent direction V, and hence the descent step sizes) from the proposed formulation exhibit an almost  Figures  4b-4d for the case λ = −9). We believe that this is due to large deformations of the domain caused by large values of descent step sizes (and therefore has to be reduced) during the optimization process. In addition to these observations, we also mention the following important remarks regarding the computational results summarized in Table 3. Firstly, for α = 0.99 and η = 10 −6 in (58), it seems that the proposed formulation requires less computing time to complete the optimization process than the classical formulation. Secondly, it appears that, for all values of λ ∈ {−10, −9, . . . , −1}, the initial step size for the classical formulation is larger (in fact, more than three times) than the magnitude of the initial step size for the proposed formulation. Lastly, we observed that the final cost values from the classical formulation are only of magnitude 10 −4 (or lower) while the proposed formulation produces final cost values that in the magnitude 10 −6 (or lower). Because of the last two remarks, it is actually difficult to say that the proposed formulation possesses faster convergence rate to the optimal solution than the classical Kohn-Vogelius formulation.
We further assess the quality of the two formulations in terms of numerically solving the exterior Bernoulli problem (1) by taking into account the above-mentioned key observations. To do this, instead of taking the same value of α for the two formulations, we choose α in such a way that the difference between the initial step sizes from the two formulations is small as possible. For simplicity, we take α = 0.99 for the proposed formulation and adjust the value of α in the classical scheme so that the initial step size for the two formulations are as close as possible. Also, for this purpose, we focus our attention to the case when λ = −9, −4, −1 since these are the cases where we see some sort of inconsistencies in the number of iterations and computing times shown in Table 3. Table 4 shows the results of the computations using the classical Kohn-Vogelius formulation with η = 10 −6 in (58). It seems that for λ = −9 (and possibly for smaller values of λ), the computing time when using   (57) the classical formulation varies greatly with respect to a small change in the value of α. In contrast, for λ = −1 (and possibly for values of λ < 0 near zero), the computing time is not sensitive to small change in α. Meanwhile, we notice that, still, the proposed formulation (refer to Table 3) requires less number of iterations and computing times to reach convergence to the optimal solution than the classical formulation.
We also examine the results of the two formulations when η is set to 10 −4 in the stopping condition (58) while taking the initial step size for the two formulations as close as possible. Among the three values of the parameter α listed in Table 4, we take the corresponding value of α 2 for each λ = −9, −4, −1 for the classical formulation. Table 5 shows the corresponding results for the given setup. Observe that, for the three cases considered, the proposed formulation requires less computing time than the classical one except for the case when λ = −4 where both formulations require 11 seconds to complete the optimization process. However, it seems that the classical formulation requires less number of iterations than the proposed scheme if we compare the number of iterations for the proposed formulation tabulated in Table 3 with that of the classical formulation shown in Table 5. Nevertheless, notice that, in most cases, the proposed formulation requires less computing time per iteration than the classical formulation. The histories of descent step sizes for the proposed and classical formulations are plotted in Figure 5a. Figure 5b, on the other hand, depicts the computed (optimal) exterior boundaries obtained through the proposed formulation versus the ones computed via the classical Kohn-Vogelius shape optimization formulation when η = 10 −4 . We observe that, in all cases examined, the computed optimal free boundaries from the two formulations are almost indistinguishable. We provide a few more numerical examples comparing the results between the proposed and classical formulation when λ = −10 for η = 10 −4 , 10 −5 , 10 −6 in (58). In this test case, we consider two different initial guesses given in Example 1. Particularly, we consider Σ 0 as the circle Σ 1 0 and as the ellipse Σ 3 0 . For the proposed formulation, we again take α = 0.99 which gives us t 0 = 0.2281500068 when Σ 0 = Σ 1 0 and t 0 = 0.1918396442 in case of taking Σ 0 = Σ 3 0 . On the other hand, by taking α = 0.2941418090 when Σ 0 = Σ 1 0 in the classical formulation, we get t 0 = 0.2281500069. Also, with α = 0.2904954592 when Σ 0 = Σ 3 0 in the classical formulation, we get t 0 = 0.1918396442. We examine the convergence of the approximate free boundaries to the optimal shape obtained through the two formulations. Because we do not know precisely the exact profile of the optimal domains corresponding to the solution of the test case in consideration, we use the computed optimal free boundary Σ in obtained using the improved Neumann-datatracking cost functional approach proposed in [41] as our benchmark. The results of the computations are summarized in Table 6 where we show the final cost values, the Hausdorff distance of the computed optimal free boundary Σ i f , i = 1, 3, with respect to Σ in , the number of iterations and the total computing times for each of the two formulations. Comparing the results of the proposed formulation when η = 10 −6 with that of the classical formulation when η = 10 −4 (see highlighted rows), we observe that the former formulation needs less computing time per iteration to complete the iteration process than the latter one. In addition, it seems that the computed optimal free boundary obtained through the proposed formulation is closer (in terms of the Hausdorff distance) to Σ in than the one obtained via the classical formulation when Σ 0 = Σ 3 0 . However, when Σ 0 = Σ 1 0 , the classical formulation produces smaller Hausdorff distance with respect to Σ in than the proposed formulation. The evolutions of the domains or histories of free boundaries computed for each of the cases considered through the two formulations are shown in Figure 6. Looking at the plots depicted in the said figures, we notice that the proposed formulation yields a more stable convergence behavior (in the sense that the shape evolution is almost monotone) to the optimal solution than the classical formulation. We also looked at the histories of the minimum, the maximum and the mean curvatures (respectively denoted by κ 2 , κ 1 and κ) of the free boundaries (see Figure 7) plotted in Figure 6, and we found out that the mean curvatures of the computed optimal free boundaries Σ f for all considered cases are positive. Hence, according to Proposition 3, the shape Hessian at the solution Ω * of the Bernoulli problem (1) when Γ = ∂S and λ = −10 (and therefore, for −10 < λ < 0,) is strictly H 1/2 -positive. Moreover, we noticed that, after a certain number of iterations, the history of minimum curvatures that corresponds to the histories of free boundaries obtained through the proposed formulation (refer to Figure 7a and Figure 7b) projects a decreasing trend. This is in contrast to the behavior of the graph of the minimum curvatures corresponding to the history of free boundaries obtained via the classical setting where we noticed an increase in the value of the minimum curvature from the 7 th to the 8 th iteration (see Figure 7c and Figure 7d).  Table 6. Comparison of computational results obtained through the proposed and classical formulations for an L-shaped fixed boundary Γ = ∂S when λ = −10 with almost the same initial step size t 0 for both formulations Before we go to our next and final set of examples, we reiterate the following key findings drawn in this subsection. For the same value of α, the classical formulation produces larger magnitude for t 0 than the proposed formulation. Moreover, it seems that the appropriate value of η in (58) is 10 −4 for the classical formulation and 10 −6 for the proposed formulation. One of the possible reason for this difference in the right choice for η is the fact that the classical formulation produces larger initial cost values than the proposed formulation. Furthermore, instead of simply comparing just the number of iterations or computing time to complete an iteration process, it seems reasonable to compare the mean computing time per iteration of the two formulations to evaluate their performance in numerically solving the optimization problem. This way of comparing the two formulations seems sensible because the two formulations almost have the same computing time. Moreover, we emphasize that we are actually applying the same algorithm for each of the two formulations, hence, considering the computing time per iteration as performance metrics in evaluating the two methods is justifiable. Lastly, we emphasize that the classical formulation requires less iteration number to complete an iteration process than the proposed formulation.
To end this section, let us complete Table 5 and compare the mean computing time per iteration and standard deviations for all λ = −10, −9, . . . , −1. Table 7 shows the final cost values, the number of iterations and computing times for all λ ∈ {−10, −8, −7, −6, −5, −3, −2} when using the classical approach. The values of α chosen for each cases and the corresponding sizes of the initial step t 0 are also shown in the table. Meanwhile, in Table 8, we show the means and standard deviations of the number of iterations, computing time and computing time per iteration from the two formulations for the present optimization problem. Notice  (58) that the two formulations are comparable in terms of the mean over-all computing time. However, it requires two additional iterations for the proposed formulation to complete the optimization process than the classical formulation. Nevertheless, the proposed formulation needs less computing time per iteration than the classical setting.

Example 3: A domain with fixed boundary having two disjoint components.
For the third and final example, we look at a similar test case studied in [34]. Particularly, we define the fixed boundary Γ as the union of two disjoint kite-shaped figures which are parametrically defined as follows: As in Example 2, we compute the optimal shape for all integers λ = −10, −9, . . . , −1 with C(0, 0.6) as the initial shape for the free boundary Σ. We set α = 0.99 for the proposed formulation and choose α for the classical formulation in such a way that it results to an initial step size having (almost) the same value as that of the  Table 9 summarizes the computational results obtained through the two formulations. It shows in particular the values of α used for each of the two formulations, the resulting initial step sizes, the final values of the costs, the number of iterations and over-all computing time to reach convergence to the optimal solution. Meanwhile, in Figure 9, we show a comparison between the results obtained through the proposed and classical formulation for the case λ = −10, −4, −1. In Figure 9a, we see that the computed optimal free boundaries obtained from the two formulations are indistinguishable from each other. The histories of the costs values, the L 2 -norms of V and the histories of the descent step sizes are plotted  Table 7. Computational results obtained through the classical formulation with η = 10 −4 in (58) for an L-shaped fixed boundary Γ = ∂S for λ = −10, −8, −7, −6, −5, −3, −2 with almost the same initial step size t 0 with respect to that of the proposed formulation shown in Table 3 formulation iteration cpu time  (58) in Figures 9b-9d, respectively. Looking at these figures, we observe that we get a more comparable convergence rate from the two formulations as λ < 0 closes to zero; that is, in case of λ = −1, the rate of convergence, for instance, of the cost for the proposed formulation is almost of the same value with that of the classical scheme. On the other hand, when λ = −10, the convergence rates of the two formulations differ greatly from each other. Meanwhile, in Figure 10, we show the evolutions or the histories of the free boundaries when λ = −10 for each of the two formulations (see, particularly, Figure 10a and Figure 10c for the respective results of the proposed and classical formulation). The figure also shows the evolutions of the free boundaries (refer to Figure 10b and Figure 10d) when the initial shape Σ 0 of the free boundary is taken to be the ellipse Υ := {(−0.1 + 0.4 cos θ, 0.05 + 0.5 sin θ) , 0 θ 2π}.
As in the previous example, we noticed that the proposed formulation yields a more stable convergence behavior (again, in the sense that the shape evolution is monotone) to the optimal solution than the classical formulation. The histories of the minimum, the maximum and the mean curvatures of the free boundaries depicted in Figure 10 are plotted in the graphs shown in Figure 11. Clearly, the optimal free boundary for the present test case has positive mean curvature. Finally, in Table  10, we show the means and standard deviations of the number of iterations, computing times and computing times per iterations for the present test case. Looking at the results shown in the table, we see that the numerical algorithm presented in Subsection 5.1.3 requires 12 iterations to complete the optimization process when using the proposed formulation. This is two iterations higher compared to when using the classical formulation. However, we get a comparable result for the two formulations in terms of computing time which means that the proposed formulation actually requires less computing time per iteration than the classical setting. Conclusion. In this work, we have presented a new shape optimization formulation of the exterior Bernoulli free boundary problem which was formulated by modifying the classical Kohn-Vogelius formulation of the original overdetermined problem. Using the first-order shape derivative of the energy-gap type cost functional studied in the paper, we have successfully devised an efficient iterative scheme based on a Lagrangian-like method to numerically solve the minimization problem. Based on the numerical results, we found that the proposed shape optimization formulation and the classical classical Kohn-Vogelius approach are comparable in terms of mean over-all computing time. The classical setting, however, requires less number of iterations to complete the optimization process than the proposed method. This fact, on the other hand, means that the proposed formulation actually demands less computing time per iteration to finish the computations. Furthermore, by inspecting the evolution of the free boundaries and the corresponding histories of their curvatures, we found out that the proposed method exhibits a more stable approximation of the optimal shape solution than the classical Kohn-Vogelius formulation, in the sense that the shape evolution of the free boundary is monotone. These observations lead us to conclude that the proposed energy-cost functional approach is somewhat more robust than the classical Kohn-Vogelius formulation. Meanwhile, using the shape Hessian information at the solution of the  Table 10. Means and standard deviations (std) of the number of iterations, computing time and computing time per iteration of the computational results shown in Table 9 free boundary problem, we have also formulated a state-of-the-art second-order optimization procedure based on a Newton-type method to numerically resolve the minimization problem. Although the second-order algorithm requires additional computing time, as expected, improvements in terms of number of iterations and accuracy of computed optimal shapes were observed from its numerical results. Consequently, the newly proposed shape optimization formulation of (1) provides an alternative computational strategy to [11] which, on the other hand, utilizes the classical Kohn-Vogelius formulation in a level-set approach to numerically solve the exterior Bernoulli free boundary problem. We expect that the computational technique offered in this investigation, especially the second-order optimization scheme, will also provide efficient numerical resolution to other related free boundary problems in the framework of shape optimization.