Locally optimal configurations for the two-phase torsion problem in the ball
Lorenzo Cavallina

Computationally Efficient Market Simulation Tool for Future Grid Scenario Analysis
    Shariq Riaz, Graduate Student Member, IEEE,
        Gregor Verbič, Senior Member, IEEE,
        and Archie C. Chapman, Member, IEEE
Shariq Riaz, Gregor Verbič and Archie C. Chapman are with the School of Electrical and Information Engineering, The University of Sydney, Sydney, New South Wales, Australia. 	e-mails: (shariq.riaz, gregor.verbic, archie.chapman@sydney.edu.au).
Shariq Riaz is also with the Department of Electrical Engineering, University of Engineering and Technology Lahore, Lahore, Pakistan.

    December 30, 2023
==================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




We consider the unit ball Ω⊂ (N≥2) filled with two materials with different conductivities. We perform shape derivatives up to the second order to find out precise information about locally optimal configurations with respect to the torsional rigidity functional. 
In particular we analyse the role played by the configuration obtained by putting a smaller concentric ball inside Ω.
In this case the stress function admits an explicit form which is radially symmetric: this allows us to compute the sign of the second order shape derivative of the torsional rigidity functional with the aid of spherical harmonics.  
Depending on the ratio of the conductivities a symmetry breaking phenomenon occurs. 







2010 Mathematics Subject classification. 49Q10



Keywords and phrases: torsion problem, optimization problem, elliptic PDE, shape derivative








§ INTRODUCTION

We will start by considering the following two-phase problem.
Let Ω⊂^N (N≥ 2) be the unit open ball centered at the origin.

Fix m∈ (0,	(Ω)), where here we denote the N-dimensional Lebesgue measure of a set by (·) . 
Let ω⊂⊂Ω be a sufficiently regular open set such that (ω)=m. 


Fix two positive constants σ_-, σ_+ and consider the following distribution of conductivities: 
    σ:=σ_ω:= σ_- 1_ω+σ_+ 1_Ω∖ω,

where by 1_A we denote the characteristic function of a set A (i.e. 1_A(x)=1 if x∈ A and vanishes otherwise).



Consider the following boundary value problem:

    -(σ_ω∇ u )= 1       in Ω, 
    
    u=0       on ∂Ω.



We recall the weak formulation of (<ref>): 

    ∫_Ωσ_ω∇ u ·∇φ = ∫_Ωφ     for all φ∈ H^1_0(Ω).


Moreover, since σ_ω is piecewise constant, we can rewrite (<ref>) as follows

    -σ_ωΔ u = 1       in ω∪(Ω∖ω), 
    σ_-∂_n u_-= σ_+∂_n u_+       on ∂ω,
    
    u=0       on ∂Ω,

where the following notation is used: the symbol  is reserved for the outward unit normal and :=∂/∂ denotes the usual normal derivative. Throughout the paper we will use the + and - subscripts to denote quantities in the two different phases (under this convention we have (_ω)_+=_+ and (_ω)_-=_- and our notations are “consistent" at least in this respect). The second equality of (<ref>) has to be intended in the sense of traces. In the sequel, the notation [f]:=f_+-f_- will be used to denote the jump of a function f through the interface ∂ω (for example, following this convention, the second equality in (<ref>) reads “u=0 on ∂ω”).  


We consider the following torsional rigidity functional: 

    E(ω):= ∫_Ωσ_ω |∇ u_ω|^2=∫_ω_- | u_ω|^2+∫_Ω∖ω̅_+ | u_ω|^2=
    ∫_Ω u_ω,

where u_ω is the unique (weak) solution of (<ref>).

Physically speaking, we imagine our ball Ω being filled up with two different materials and the constants _± represent how “hard" they are. The second equality of (<ref>), which can be obtained by integrating by parts after splitting the integrals in (<ref>), is usually referred to as transmission condition in the literature and can be interpreted as the continuity of the flux through the interface .

The functional E, then, represents the torsional rigidity of an infinitely long composite beam. Our aim is to study (locally) optimal shapes of ω with respect to the functional E under the fixed volume constraint.

The one-phase version of this problem was first studied by Pólya. 
Let D⊂ (N≥ 2) be a bounded Lipschitz domain. 
Define the following shape functional 

    ℰ(D):= ∫_D | u_D|^2,
 
where the function u_D (usually called stress function) is the unique solution to 

    - u = 1      in  D,
    
    u=0     on ∂ D.


The value ℰ(D) represents the torsional rigidity of an infinitely long beam whose cross section is given by D. 
The following theorem (see <cit.>) tells us that beams with a spherical section are the “most resistant”. 


The ball maximizes ℰ among all Lipschitz domains with fixed volume.


Inspired by the result of Theorem <ref> it is natural to expect radially symmetrical configurations to be optimizers of some kind for E (at least in the local sense).
From now on we will consider ω:= B_R (the open ball centered at the origin, whose radius, 0<R<1, is chosen to verify the volume constraint |B_R|=m) and use shape derivatives to analyze this configuration. 
This technique has already been used by Conca and Mahadevan in <cit.> and Dambrine and Kateb in <cit.> for the minimization of the first Dirichlet eigenvalue in a similar two-phase setting (Ω being a ball) and it can be applied with ease to our case as well. 

A direct calculation shows that the function u, solution to (<ref>) where ω=, has the following expression:

    u(x)= 1-R^2/2N_++R^2-|x|^2/2N_-     for |x|∈ [0,R],
    1-|x|^2/2N_+     for  |x|∈ [R,1].


In this paper we will use the following notation for Jacobian and Hessian matrix respectively.

    (D v)_ij:=∂ v_i/∂ x_j,    (D^2 f)_ij= ∂^2 f/∂ x_i∂ x_j,

for all smooth real valued function f and vector field v=(v_1,…, v_N) defined on Ω.
We will introduce some differential operators from tangential calculus that will be used in the sequel. For smooth f and v defined on ∂ω we set

    _τ f    := f-(f·)     ( tangential gradient),   
    _τ v   :=  v-·( D  v)     (tangential divergence),

where f and v are some smooth extensions on the whole Ω of f and v respectively. It is known that the differential operators defined in (<ref>) do not depend on the choice of the extensions.   
Moreover we denote by D_τ v the matrix whose i-th row is given by _τ v_i.
We define the (additive) mean curvature of ∂ω as H:=_τ (cf. <cit.>). According to this definition, the mean curvature H of ∂ is given by  (N-1)/R.  





A first key result of this paper is the following. 

For all suitable perturbations that fix the volume (at least at first order), the first shape derivative of E at B_R vanishes. 


An improvement of Theorem <ref> is given by the following precise result (obtained by studying second order shape derivatives).

Let _-,_+>0 and R∈(0,1). If _->_+ then B_R is a local maximizer for the functional E under the fixed volume constraint.

On the other hand, if _-<_+ then B_R is a saddle shape for the functional E under the fixed volume constraint. 



In section 2 we will give the precise definition of shape derivatives and introduce the famous Hadamard forumlas, a precious tool for computing them. In the end of section 2 a proof of Theorem <ref> will emerge as a natural consequence of our calculations. 
Section 3 will be devoted to the computation of the second order shape derivative of the functional E in the case ω=.
In Section 4 we will finally calculate the sign of the second order shape derivative of E by means of the spherical harmonics. The last section contains an analysis of the different behaviour that arises when volume preserving transformations are replaced by surface area preserving ones.  





§ COMPUTATION OF THE FIRST ORDER SHAPE DERIVATIVE:
 PROOF OF THEOREM <REF>


We consider the following class of perturbations with support compactly contained in Ω:

    𝒜:={Φ∈^∞([0,1)×,)   | Φ(0,·)=  Id, ∃ K⊂⊂Ω s.t. Φ(t,x)=x  ∀ t∈ [0,1) , ∀ x∈∖ K        }.


For Φ∈𝒜 we will frequently write Φ(t) to denote Φ(t,·) and, for all domain D in , we will denote by Φ(t)(D) the set of all Φ(t,x) for x∈ D. We will also use the notation D_t:=Φ(t)(D) when it does not create confusion.    
In the sequel the following notation for the first order approximation (in the “time” variable) of Φ will be used.

    Φ(t)= Id + t  +o(t)   as t→0,

where  is a smooth vector field. In particular we will write h_n:=· (the normal component of ) and _τ : = - on the interface. 
We are ready to introduce the definition of shape derivative of a shape functional J with respect to a deformation field Φ in 𝒜 as the following derivative along the path associated to Φ.

    d/dt J(D_t)t=0 = lim_t→ 0J(D_t)-J(D)/t.


This subject is very deep. Many different formulations of shape derivatives associated to various  kinds of deformation fields have been proposed over the years. We refer to <cit.> for a detailed analysis on the equivalence between the various methods. For the study of second (or even higher) order shape derivatives and their computation we refer to <cit.>.   

The structure theorem for first and second order shape derivatives (cf. <cit.>, Theorem 5.9.2, page 220 and the subsequent corollaries) yields the following expansion. 
For every shape functional J, domain D and pertubation field Φ in 𝒜, under suitable smoothness assumptions the following holds.


    J(D_t)=J(D)+t   l_1^J (D)(h_n)+t^2/2(  l_2^J(D)(,)+l_1^J(D)(Z)  	) + o(t^2)    as t→ 0,

for some linear l_1^J(D): ^∞ (∂ D)→ and bilinear form l_2^J(D): ^∞ (∂ D)×^∞ (∂ D)→ to be determined eventually. 
Moreover for the ease of notation we have set 

    Z:=( V'+D)·+ ((D_τ) )·-2_τ·,

where V(t,Φ(t)):=∂_t Φ(t) and V':=∂_t V(t,·).


According to the expansion (<ref>), the first order shape derivative of a shape functional depends only on its first order apporximation by means of its normal components.
On the other hand the second order derivative contains an “acceleration” term l_1^J(D)(Z).  
It is woth noticing that, (see Corollary 5.9.4, page 221 of <cit.>) Z vanishes in the special case when Φ=𝕀+t with _τ= on ∂ D (this will be a key observation to compute the bilinear form l_2^J).

We will state the following lemma, which will aid us in the computations of the linear and bilinear forms l_1^J(D) and l_2^J(D) for various shape functionals (cf.  <cit.>, formula (5.17), page 176 and formulas (5.110) and (5.111), page 227).

Take Φ∈ and let f=f(t,x)∈^2([0,T), L^1())∩^1([0,T), W^1,1()). For every smooth domain D in  define
J(D_t):=∫_D_t f(t)
(where we omit the space variable for the sake of readability). Then the following identities hold:

    l_1^J(D)()= ∫_D ∂_t ft=0 + ∫_∂ D f(0) h_n,


    l_2^J(D)(,)=∫_D ∂_tt^2 ft=0+∫_∂ D2∂_t ft=0h_n+(H f(0)+∂_n f(0) )h_n^2.



Since we are going to compute second order shape derivatives of a shape functional subject to a volume constraint, we will need to restric our attention to the class of perturbations in  that fix the volume of ω:

    ℬ(ω):= {Φ∈𝒜  | (Φ(t)(ω))=(ω)=m   for all  t∈ [0,1)  }.

We will simply write  in place of ().
Employing the use of Lemma <ref> for the volume functional  and of the expansion (<ref>), for all Φ∈ we get

    (ω_t)=(ω)+ t ∫_∂ω + t^2/2( ∫_∂ω H h_n^2 + ∫_∂ω Z     ) + o(t^2)  as t→0.

This yields the following two conditions: 

    ∫_∂ω =0,           (1^ st order volume preserving)
       ∫_∂ω H h_n^2+ ∫_∂ωZ=0.    (2^ nd order volume preserving)
 

For every admissible perturbation field Φ=𝕀+t in , with  satisfying (<ref>), we can find some perturbation field Φ∈ such that Φ=𝕀 + t  + o(t) as t→ 0.
For example, the following construction works just fine:

    Φ(t,x)=Φ(t,x)/η(x)( (Φ(t)(ω))/(ω))^1/N+(1-η(x)),

where η is a suitable smooth cutoff function compactly supported in Ω that attains the value 1 on a neighbourhood of ω. 




We will now introduce the concepts of “shape" and “material" derivative of a path of real valued functions defined on Ω.
Fix an admissible perturbation field Φ∈ and let u=u(t,x) be defined on [0,1)×Ω for some positive T.
Computing the partial derivative with respect to t at a fixed point x∈Ω is usually called shape derivative of u; we will write:

    u'(t_0,x):= ∂ u/∂ t (t_0,x),   for x∈Ω, t_0∈ [0,1).
 
On the other hand differentiating along the trajectories gives rise to the material derivative:

    u̇(t_0,x):= ∂ v/∂ t(t_0,x),  x∈Ω, t_0∈ [0,1);
 
where v(t,x):=u(t, Φ(t,x)).
From now on for the sake of brevity we will omit the dependency on the “time" variable unless strictly necessary and write u(x), u'(x) and u̇(x) foru(0,x), u'(0,x) and u̇(0,x). 
The following relationship between shape and material derivatives hold true:

    u'   =u̇- u ·.


We are interested in the case where u(t,·):=u_(B_R)_t i.e. it is the solution to problem (<ref>) when ω=Φ(t)(B_R).
In this case, since by symmetry we have u = ( u), the formula above admits the following simpler form on the interface ∂ B_R:

    u'=u̇-( u) .


It is natural to ask whether the shape derivatives of the functional E are well defined. Actually, by a standard argument using the implicit function theorem for Banach spaces (we refer to <cit.> for the details) it can be proven that the application mapping every smooth vector field h compactly supported in Ω to E((𝕀+)(ω)) is of class 𝒞^∞ in a neighbourhood of h = 0. This implies the shape differentiability of the functional E for any admissible deformation field Φ∈. As a byproduct we obtain the smoothness of the material derivative u̇.

As already remarked in <cit.> (Remark 2.1), in contrast to material derivatives, the shape derivative u' of the solution to our problem has a jump through the interface. This is due to the presence of the gradient term in formula (<ref>) (recall that the transmission condition provides only the continuity of the flux). On the other hand we will still be using shape derivatives because they are easier to handle in computations (and writing Hadamard formulas using them is simpler). 



For any given admissible Φ∈, the corresponding u' can be characterized as the (unique) solution to the following problem in the class of functions whose restriction to both  and Ω∖ is smooth:

    u' =0      in ∪ (Ω∖),
     u'=0      on ,
    u'=- u     on ,
    
    u'= 0      on ∂Ω.



Let us now prove that u' solves (<ref>).
First we take the shape derivative of both sides of the first equation in (<ref>) at points away from the interface:

    Δ u' = 0  in ω∪ (Ω∖).

In order to prove that u' vanishes on ∂ B_R we will proceed as follows.
We performing the change of variables y:=Φ(t,x) in (<ref>) and set φ(x)=:ψ( Φ(t,x) ). Taking the derivative with respect to t, bearing in mind the first order approximation of Φ given by (<ref>) yields the following.

    σ( - D∇ u + ∇u̇)·∇ψ -  u · Dψ +  u ·ψ = ψ.


Rearranging the terms yields:

    u̇·ψ +( -D  - D ^T + () I ) u ·ψ_⊛ = 
    -·ψ.


Let x and y be two sufficiently smooth vector fields in  such that D  x= ( D  x) ^T and D  y= ( D  y) ^T. It is easy to check that the following identity holds:

    ( -D  - D^T + () I )  x· y = (( x· y) ) - (· x)· y- (· y)· x.


We can apply this identity with x= u and y= ψ to rewrite ⊛ as follows:

    ( -D  - D ^T + () I ) u ·ψ=( u ·ψ)_1 -(· u)·ψ_2 -(·ψ)· u_3.

Thus 

    u̇·ψ -  u ·ψ- (· u)·ψ+_-  u_- ψ =0,

where we have split the integrals and integrated by parts to handle the terms coming from 1 and 3. 

Now, merging together the integrals on Ω in the left hand side by (<ref>) and exploiting the fact that u = ( u) on , the above simplifies to 

    u' ·ψ = 0.





Splitting the domain of integration and integrating by parts, we obtain

    0=-∫__- Δ u_-' ψ + ∫__-  u_-' ψ -∫_Ω∖_+ Δ u_+' ψ -∫__+ u_+' ψ-∫_∂Ω_+  u_+' ψ
    
    =  u'ψ,

where in the last equality we have used (<ref>) and the fact that ψ vanishes on ∂Ω.
By the arbitrariness of ψ∈ H_0^1(Ω) we can conclude that u'=0 on ∂.  
The remaining conditions of problem (<ref>) are a consequence of (<ref>).

To prove uniqueness for this problem in the class of functions whose restriction to both  and Ω∖ is smooth, just consider the difference between two solutions of such problem and call it w. Then w solves

    w =0      in ∪ (Ω∖),
     w=0      on ,
    w=0     on ,
    
    w=0      on ∂Ω;
 
in other words, w solves 

    - ( w)=0      in Ω,
    
    w=0      on ∂Ω.

Since the only solution to the problem above is the constant function 0, uniqueness for Problem (<ref>) is proven.


We emphasize that formulas (<ref>) and (<ref>) are valid only for f belonging at least to the class ^2([0,T), L^1())∩^1([0,T), W^1,1()). 
We would like to apply them to f(t)=u_t and f(t)=σ_t |∇ u_t|^2, where _t and u_t are the distribution of conductivities and the solution of problem (<ref>) respectively corresponding to the case ω=()_t. On the other hand, u_t is not regular enough in the entire domain Ω, despite being fairly smooth in both ω_t and Ω_t∖ω_t: therefore we need to split the integrals in order to apply (<ref>) and (<ref>) (this will give rise to interface integral terms by integration by parts).



For all Φ∈ we have

    l_1^E(B_R)()=- | u|^2.

In particular, for all Φ satisfying the first order volume preserving condition (<ref>) we get l_1^E(B_R)()=0.


We apply formula (<ref>) to 
E(ω_t)=∫_Ω u_t=∫_ω_t u_t+∫_ u_t
to get


    l_1^E(B_R)()= ∫_ u_-'+∫_u_- h_n + ∫_Ω∖B_R u_+'
    -∫_ u_+ h_n.

Using the jump notation we rewrite the previous expression as follows

    l_1^E(B_R)()= ∫_Ω u' - ∫_ [u h_n]=∫_Ω u';

notice that the surface integral in (<ref>) vanishes as both u and h_n are continuous through the interface.


Next we apply (<ref>) to E(ω_t)=∫_Ωσ_t |∇ u_t|^2.

    l_1^E(B_R)()=     2∫_σ_- ∇ u_-·∇ u_-' + ∫_σ_- |∇ u_-|^2 h_n+ 
        2∫_σ_+ ∇ u_+·∇ u_+' +∫_∂Ωσ_+ |∇ u_+|^2 h_n - ∫_σ_+ |∇ u_+|^2 h_n.





Thus we get the following:

    l_1^E(B_R)()= 2∫_Ωσ∇ u·∇ u' -σ |∇ u|^2h_n.
























Comparing (<ref>) (choose ψ=u) with 
(<ref>) gives

    l_1^E(B_R)()= - | u|^2.

By symmetry, the term σ |∇ u|^2 is constant on  and can be moved outside the integral sign. Therefore we have

    l_1^E(B_R)()= 0  for all Φ satisfying (<ref>).

This holds in particular for all Φ∈. 







§ COMPUTATION OF THE SECOND ORDER SHAPE DERIVATIVE


The result of the previous chapter tells us that the configuration corresponding to  is a critical shape for the functional E under the fixed volume constraint. In order to obtain more precise information, we will need an explicit formula for the second order shape derivative of E.
The main result of this chapter consists of the computation of the bilinear form l_2^E(B_R)(,).


For all Φ∈ we have 

    =-2_- u_-  u' -2_- u_- ∂_nn^2 uh_n^2-_-  u_- u H h_n^2.



Take Φ=𝕀+t in  with = on . As remarked after (<ref>), Z vanishes in this case. We get 

    =d^2/dt^2E(Φ(t)(B_R))t=0.

Hence, substituting the expression of the first order shape derivative obtained in Theorem <ref> yields

    = -d/dt( ∫_(∂)_t_t | u_t|^2∘Φ^-1(t)   ) t=0.

We unfold the jump in the surface integral above and apply the divergence theorem to obtain

    =d/dt( ∫_(B_R)_t(  _- | u_t|^2 ∘Φ^-1(t) )  -  ∫_Ω∖ (B_R)_t(  _+ | u_t|^2 ∘Φ^-1(t) )  )t=0

We will treat each integral individually.
By (<ref>) we have ∂_t (  Φ)t=0=-, therefore ∂_t(∘Φ^-1)t=0=-D. Now set f(t):=_- | u_t|^2.
By (<ref>) we have

    d/dt(∫_(B_R)_t( f(t) ∘Φ^-1(t) ))t=0 =   ∫_∂_t(  (  f(t) ∘Φ^-1(t)  )  )t=0_(A)
     +  (f(0))_(B).

We have

    (A)=∫_( ∂_t ft=0 + f(0)∂_t (∘Φ^-1(t))t=0) =    
      ∫_( ∂_t ft=0) - ∫_( f(0)D) =  
    ∫_∂_t ft=0 -  f(0) · D.

On the other hand

    (B)= (  f(0) ) = (   f(0)· + f(0)) .

Using the fact that = and -· D=: _τ ()=H (c.f. Equation (5.22), page 366 of <cit.>) we get 

    (A)+(B)=  f'  +  ( f+ Hf)h_n^2.

Substituting f(t)=_-| u_t|^2 yields

    (A)+(B)= 2 _-  u_- · u_-'  + 2_-  u_- ∂_nn^2 u_- h_n^2 + _- | u_-|^2 H h_n^2.

The calculation for the integral over Ω∖ ()_t in (<ref>) is analogous.
We conclude that 

    = -2 _-  u_-  u' -2_-  u_- ∂_nn^2 uh_n^2- _-  u_- u H h_n^2.


















The following is a first expression we obtain for E”()[,].

With the same notation as before, we have:

    l_2^E(B_R)(,)= 2  u h_n^2 -2 u  u'
    
    - 2  u  u· (D)   +  H u  u h_n^2.


In order to prove Proposition <ref> we will need the following identities concerning shape derivatives of u (whose proofs will be postponed for the sake of clarity).








    | u'|^2=  u  u'.




    u· u”=  u”+ 2  u  u'
    
    +  u ∂_nn u +  u  u· (D )  .





First of all we write the torsional rigidity functional as E(ω_t)= u_t and then apply formula (<ref>) to get

    E”=E”()[,]=  u” - 2 u' - Hu h_n^2- u h_n^2.

Since u is continuous through the interface, the term containing u vanishes. Moreover, rewriting u' using (<ref>) and recalling that the material derivative u̇ has no jump through the interface, we get

    E” =  u”+  u.

Next we apply (<ref>) to the other expression E(ω_t)=_t | u_t|^2:

    E”= 2 | u'|^2+ 2 u · u” - 4 u  u'
    
    -  H  u  u - 2 u ∂_nn u.

As usual our aim is to make all volume integral disappear by wise substitutions.
By Lemma <ref>, 

    E”= 2 u · u” - 2 u  u'
    
    -  H  u  u - 2 u ∂_nn u.

Then, substituting the term containing u· u ” by Lemma <ref> yields

    E”= 2  u” + 2  u  u· (D )   + 2  u  u'-  H  u  u.

We finally prove the claim by combining (<ref>) with (<ref>).


Let us now prove Lemmas <ref> to <ref>.

Let us multiply both sides of (<ref>) by u' and integrate on  and Ω∖.
We get

    0= ∫__-  u' u' +∫_Ω∖_+  u' u'.

Integrating by parts to get the boundary terms yields 

    0= - | u'|^2- u' u'.

Using the fact that u' has no jump (Proposition <ref>), we get


    0= - | u'|^2- u' u'.

By (<ref>) we have

    0= - | u'|^2+ u'  u.

Rearranging the terms above and using the transmission condition we conclude as follows

    | u'|^2 =  u'  u= u'  u=  u  u'.



Multiplying the first equation in (<ref>) by u” and integrating by parts on  and Ω∖ we get

    u” = -∫_ u u” - ∫_Ω∖ u u”=  u · u” + u u”.

We conclude by rewriting u” using (<ref>).


Having finally proven the two lemmas above concludes the proof of Proposition <ref>. What we got here is an expression for the shape Hessian of E at  made up of the sum of four surface integrals. There are a few reasons why this is still not enough for our purposes. First of all, we succeeded in removing the dependence on u” but still we are left with a term involving the normal derivative of u'. Unlike the function u, its shape derivative u' depends on the choice of  and fails to have an explicit formula to express it; we will have to deal with this problem later in the paper. 
Furthermore, the integral depending on the Jacobian of  can be simplified as follows. 

Employing the use of Eq.(5.19) of page 366 of <cit.> we get

    · (D ) =  - _ h  on ,

then, since we are working with a constrained problem, we can restrict ourselves to deformation fields that are tangent to the constraint, i.e. in our case we can suppose without loss of generality that =0 in a neighbourhood of  (cf. <cit.>).
Therefore we may write,

    · (D ) = - _ h= - _()=-H  on ,

where we have used Eq.(5.22) of page 366 of <cit.>. 

Thus, for any divergence free admissible =, the result of Proposition <ref> simplifies to 

    E” ()[,]= 2  u h_n^2 -2 u  u' + 3 H u  u h_n^2.






§ CLASSIFICATION OF THE CRITICAL SHAPE :
 PROOF OF THEOREM <REF>


In order to classify the critical shape  of the functional E under the volume constraint we will use the expansion shown in (<ref>).
For all Φ∈ and t>0 small, it reads

    E( Φ(t)())= E()+ t^2/2(   - | u|^2 Z  )+o(t^2).

Employing the use of the second order volume preserving condition (<ref>) and the fact that, by symmetry, the quantity | u|^2 is constant on the interface  we have

    - | u|^2 Z =  | u|^2 H h_n^2.

Combining this with the result of Theorem <ref> yields

    E( Φ(t)())= E()+ t^2{  -_- u_- u' -_- u_- ∂_nn^2 uh_n^2         }+o(t^2).

We will denote the expression between braces in the above by Q(). Since u' depends linearly on  (see (<ref>)), it follows immediately that Q() is a quadratic form in .
Since both u and u' verify the transmission condition (see (<ref>) and (<ref>)) we have 

    _- u_- u'= u u'= _-u'_- u on .

Using the explicit expression of u given in (<ref>), after some elementary calculation we write

    Q()=R/N(1/_--1/_+)( -_- u'_-  +1/N h_n^2  ).
 

In the following we will try to find an explicit expression for u'. To this end we will perform the spherical harmonic expansion of the function :→.
We set 

    (Rθ)=∑_k=1^∞∑_i=1^d_kα_k,i Y_k,i(θ)    for all θ∈∂ B_1 
    
    .

The functions Y_k,i are called spherical harmonics in the literature. They form a complete orthonormal system of L^2(∂ B_1) and are defined as the solutions of the following eigenvalue problem:

    -_τ Y_k,i=λ_k Y_k,i   on ∂ B_1,

where _τ:= _τ_τ is the Laplace-Beltrami operator on the unit sphere.
We impose the following normalization coniditon 

    ∫_∂ B_1 Y_k,i^2=R^1-N.

The following expressions for the eigenvalues λ_k and the corresponding multiplicities d_k are also known:

    λ_k= k(k+N-2),      d_k= N+k-1k-N+k-2k-1.

Notice that the value k=0 had to be excluded from the summation in (<ref>) because  verifies the first order volume preserving condition (<ref>).

Let us pick an arbitrary k∈{1,2,…} and i∈{1,…, d_k}. We will use the method of separation of variables to find the solution of problem (<ref>) in the particular case when (Rθ)=Y_k,i(θ), for all θ∈∂ B_1.


Set r:=|x| and, for x 0, θ:=x/|x|. 
We will be searching for solutions to (<ref>) of the form u'=u'(r,θ)=f(r)g(θ).
Using the well known decomposition formula for the Laplacian into its radial and angular components, the equation u'=0 in ∪ (Ω∖) can be rewritten as

    0= u'(x) = f_rr(r)g(θ)+N-1/rf_r(r)g(θ)+1/r^2f(r)_τg(θ)  for r∈(0,R)∪ (R,1), θ∈∂ B_1.

Take g=Y_k,i.
Under this assumption, we get the following equation for f:

    f_rr+N-1/rf_r-λ_k/r^2f=0   in  (0,R)∪ (R,1).

It can be easily checked that the solutions to the above consist of linear combinations of  r^η and r^ξ, where

    η   =η(k)=1/2( 2-N + √( (N-2)^2+4λ_k ))=k,
    ξ   =ξ(k)= 1/2( 2-N - √( (N-2)^2+4λ_k ))=2-N-k.

Since equation (<ref>) is defined for r∈ (0,R)∪ (R,1), we have that the following holds for some real constants A, B, C and D;

    f(r)= 
    Ar^2-N-k+Br^k      for  r∈(0,R),
    
    Cr^2-N-k+Dr^k      for  r∈(R,1).

Moreover, since 2-N-k is negative, A must vanish, otherwise a singularity would occur at r=0.
The other three constants can be obtained by the interface and boundary conditions of problem (<ref>) bearing in mind that u'(r,θ)=f(r)Y_k,i(θ)=f(r)h_n(R θ).
We get the following system:

    C R^2-N-k+ DR^k-BR^k= -R/N_- + R/N _+,
    _- kB  R^k-1= _+ (2-N-k) C R^2-N-k + _+ k D R^k-1,
    
    C+D=0.

Although this system of equations could be easily solved completely for all its indeterminates, we will just need to find the explicit value of B in order to go on with our computations.
We have 

    B=B_k=R^1-k/N_-·k(_–_+)R^k-(2-N-k)(_–_+)R^2-N-k/k(_–_+)R^k+((2-N-k)_+-k_-)R^2-N-k.

Therefore, in the particular case when (R ·)=Y_k,i, 

    u'_-=u'_-(r,θ)=B_k r^k Y_k,i(θ),    r∈[0,R),  θ∈∂ B_1,
 
where B_k is defined as in (<ref>).
By linearity, we can recover the expansion of u'_- in the general case (i.e. when (<ref>) holds):

    u'_-(r,θ)= ∑_k=1^∞∑_i=1^d_kα_k B_k r^k Y_k,i(θ),    r∈[0,R), θ∈∂ B_1,   and therefore
       u'_-(R,θ)=∑_k=1^∞∑_i=1^d_kα_k B_k k R^k-1 Y_k,i(θ),       θ∈∂ B_1.

We can now diagonalize the quadratic form Q, in other words we can consider only the case (R ·)=Y_k,i for all possible pairs (k,i). 

We can write Q as a function of k as follows:

    Q()=Q(k)=R/N(_+-_-/_+_-)(-_-B_k k R^k-1 +1/N)=
    R/N^2( _+-_-/_+_-)(1- k k(_–_+)R^k-(2-N-k)(_–_+)R^2-N-k/k(_–_+)R^k+((2-N-k)_+-k_-)R^2-N-k).


The following lemma will play a central role in determining the sign of Q(k) and hence proving Theorem <ref>.

For all R∈ (0,1) and _-,_+>0, the function k↦ Q(k) defined in (<ref>) is monotone decreasing for k≥ 1. 








Let us denote by ρ the ratio of the the conductivities, namely ρ:= _-/_+.
















We get 

    Q(k)=R/N^2( 1-ρ/_-)(1- k k(ρ-1)R^k-(2-N-k)(ρ-1)R^2-N-k/k(ρ-1)R^k+((2-N-k)-kρ)R^2-N-k).

In order to prove that the map k↦ Q(k) is monotone decreasing it will be sufficient to prove that the real function

    j(x):=xx-(2-N-x)R^2-N-2x/(1-ρ)x+ (-2+N+x+ρ x )R^2-N-2x

is monotone increasing in the interval (1,∞). Notice that this does not depend on the sign of ρ-1. 
From now on we will adopt the following notation:

    L:=R^-1>1,   M:=N-2≥ 0,   P=P(x):= L^2x+M.

Using the notation introduced above, j can be rewritten as follows 

    j(x)=x^2+(x^2+Mx)P/(1-ρ)x+(x+M+ρ x)P.

In order to prove the monotonicity of j, we will compute its first derivative and then study its sign.
We get 

    j'(x)=MP(MP+2Px+2x)+x^2(P+1)^2+ρ x^2P (P-1/P-4xlog(L)-2Mlog(L))/((1-ρ)x+(x+M+ρ x)P )^2.
 
The denominator in the above is positive and we claim that also the numerator is. To this end it suffices to show that the quantity multiplied by ρ x^2P in the numerator, namely P-1/P-4xlog(L)-2Mlog(L), is positive for x∈ (1,∞) (although, we will show a stronger fact, namely that it is positive for all x>0).

    d/dx(P-1/P-4xlog(L)-2Mlog(L))= 2log(L)(P+1/P-2)>0   for x>0,

where we used the fact that L>1 and that P↦ P+P^-1-2 is a non-negative function vanishing only at P=1 (which does not happen for positive x).
We now claim that 

    (P-1/P-4xlog(L)-2Mlog(L))x=0= L^M-1/L^M-2M log(L)≥ 0.

This can be proven by an analogous reasoning: treating M as a real variable and differentiating with respect to it, we obtain

    d/dM( L^M-1/L^M-2M log(L) ) = log(L)(L^M+1/L^M-2)≥ 0

(notice that the equality holds only when M=0), moreover,

    ( L^M-1/L^M-2M log(L) )M=0=0,

which proves the claim.








We are now ready to prove the main result of the paper.

Let _-,_+>0 and R∈(0,1). If _->_+
then 

    d^2/dt^2E(Φ(t)(B_R))t=0<0    for all Φ∈.





Hence, B_R is a local maximizer for the functional E under the fixed volume constraint. 
On the other hand, if _-<_+, then there exist some Φ_1 and Φ_2 in , such that 

    d^2/dt^2E(Φ_1(t)(B_R))t=0<0,  d^2/dt^2E(Φ_2(t)(B_R))t=0>0.

In other words, B_R is a saddle shape for the functional E under the fixed volume constraint. 



We have 

    Q(1)=R/N^2( 1-ρ/_-) Nρ/ρ(R^N+1)     +N-R^N-1.

Since N≥ 2, R∈ (0,1), we have N-R^N-1>0 and therefore it is immediate to see that Q(1) and 1-ρ have the same sign. 

If _->_+, then, by Lemma <ref>, we get in particular that Q(k) is negative for all values of k≥ 1. This implies that the second order shape derivative of E at B_R is negative for all Φ∈ and therefore B_R is a local maximizer for the functional E under the fixed volume constraint as claimed.

On the other hand, if _-<_+, by (<ref>) we have Q(1)>0. 
An elementary calculation shows that, for all _-,_+>0,    

    lim_k→∞ Q(k)=-∞.

Therefore, when _-<_+, B_R is a saddle shape for the functional E under the fixed volume constraint.




§ THE SURFACE AREA PRESERVING CASE

The method employed in this paper can be applied to other constraints without much effort. For instance, it might be interesting to see what happens when volume preserving perturbations are replaced by surface area preserving ones. Is B_R a critical shape for the functional E even in the class of domains of fixed surface area? If so, of what kind?
We set (D):=∫_∂ D 1 for all smooth bounded domain D⊂. 

The following expansion for the functional  can be obtained just as we did for (<ref>):

    (ω_t)=(ω)+ t ∫_∂ωH  + t^2/2( l_2^(ω) (,) + ∫_∂ωH Z     ) + o(t^2)  as t→0,

where, (cf. <cit.>, page 225)  

    l_2^(ω)(,)=∫_∂ω |_τ|^2 +∫_∂ω( H^2 -  tr((D_τ)^T D_τ) ) h_n^2.

We get the following first and second order surface area preserving conditions. 

    ∫_∂ωH  =0, 
       ∫_∂ω |_τ|^2 +∫_∂ω( H^2 -  tr((D_τ)^T D_τ) ) h_n^2+ ∫_∂ωH Z=0.


Notice that when ω=B_R, the first order surface area preserving condition is equivalent to the first order volume preserving condition (<ref>) and therefore, by Theorem <ref>, B_R is a critical shape for E under the fixed surface area constraint as well.   

The study of the second order shape derivative of E under this constraint is done as follows. Employing the use of (<ref>) together with the second order surface area preserving condition in (<ref>) we get

    d^2/dt^2E(Φ(t)(B_R))t=0= l_2^E(B_R)(,) +  | u|^2/H l_2^(B_R)(,).

In other words, we managed to write the shape Hessian of E as a quadratic form in . We can diagonalize it by considering (R·)=Y_k,i for all possible pairs (k,i), where we imposed again the normalization (<ref>). Under this assumption, by (<ref>) we get

    ∫_ |_τ|^2= λ_k/R^2=k(k+N-2)/R^2.
  
We finally combine the expression for l_2^E(B_R) of Theorem <ref> with that of l_2^ (<ref>) to obtain

    E(Φ(t)(B_R))= E(B_R) + t^2  Q(k) + o(t^2)  as t→ 0,

where

    Q(k)=R/N^2( 1-ρ/_-)(3/2     - k(k+N-2)/2(N-1) -k k(ρ-1)R^k-(2-N-k)(ρ-1)R^2-N-k/k(ρ-1)R^k+((2-N-k)-kρ)R^2-N-k).


It is immediate to check that Q(1)=Q(1) and therefore, Q(1) is negative for _->_+ and positive otherwise. On the other hand, lim_k→∞Q(k)=∞ for _->_+ and lim_k→∞Q(k)=-∞ for _-<_+. In other words, under the surface area preserving constraint B_R is always a saddle shape, independently of the relation between _- and _+. 


We can give the following geometric interpretation to this unexpected result. 
Since the case k=1 corresponds to deformations that coincide with translations at first order, it is natural to expect a similar behaviour under both volume and surface area preserving constraint. On the other hand, high frequency perturbations (i.e. those corresponding to a very large eigenvalue) lead to the formation of indentations in the surface of B_R. Hence, in order to prevent the surface area of B_R from expanding, its volume must inevitably shrink (this is due to the higher order terms in the expansion of Φ). This behaviour can be confirmed by looking at the second order expansion of the volume functional under the effect of a surface area preserving transformation Φ∈ on the ball:

    (Φ(t)())=()+t^2/2( 1/R-k(k+N-2)/(N-1)R) +o(t^2)  as t→ 0.

We see that the second order term vanishes when k=1, while getting arbitrarily large for k≫1.
Since this shrinking effect becomes stronger the larger k is, this suggests that the behaviour of E(Φ(t)(B_R)) for large k might be approximated by that of the extreme case ω=∅. For instance, when _->_+ we have that E(B_R)<E(∅) and this is coherent with what we found, namely Q(k)>0 for k≫ 1.  




§ ACKNOWLEDGMENTS
 
This paper is prepared as a partial fulfillment of the author's doctor's degree at Tohoku University.
The author would like to thank Professor Shigeru Sakaguchi (Tohoku University) for his precious help in finding interesting problems and for sharing his naturally optimistic belief that they can be solved.
Moreover we would like to thank the anonymous referee, who suggested to study the surface area preserving case and helped us find a mistake in our calculations. Their detailed analysis and comments on the previous version of this paper, contributed to make the new version shorter and more readable.  

10 

L. Ambrosio, G. Buttazzo, An optimal design problem with perimeter penalization, Calc. Var. Part. Diff. Eq. 1 (1993): 55-69.

conca C.Conca, R.Mahadevan, L.Sanz, Shape derivative for a two-phase eigenvalue problem and optimal configuration in a ball. In CANUM 2008, ESAIM Proceedings 27, EDP Sci., Les Ulis, France, (2009): 311-321.


bandle C. Bandle, A. Wagner, Second domain variation for problems with Robin boundary conditions. J. Optim. Theory Appl. 167 (2015), no. 2: 430-463.

sensitivity M. Dambrine, D. Kateb, On the shape sensitivity of the first Dirichlet eigenvalue for two-phase problems. Applied Mathematics and Optimization 63.1 (Feb 2011): 45-74. 

stability M. Dambrine, J. Lamboley, Stability in shape optimization with second variation. arXiv:1410.2586v1 [math.OC] (9 Oct 2014).

SG M.C. Delfour, Z.-P. Zolésio, Shapes and Geometries: Analysis, Differential Calculus, and Optimization. SIAM, Philadelphia (2001).

dephfig G. De Philippis, A. Figalli, A note on the dimension of the singular set in free interface problems, Differential Integral Equations Volume 28, Number 5/6 (2015): 523-536.

espfusc L. Esposito, N. Fusco, A remark on a free interface problem with volume constraint, J. Convex Anal. 18 (2011): 417-426.

GT D. Gilbarg, N.S. Trudinger, Elliptic Partial Differential Equation of Second Order, second edition, Springer.

henrot A. Henrot, M. Pierre, Variation et optimisation de formes. Mathématiques & Applications. Springer Verlag, Berlin (2005).

new  R. Hiptmair, J. Li, Shape derivatives in differential forms I: an intrinsic perspective, Ann. Mate. Pura Appl. 192(6) (2013): 1077-1098.

lar C.J. Larsen, Regularity of components in optimal design problems with perimeter penalization, Calc. Var. Part. Diff. Eq. 16 (2003): 17-29.

lin F.H. Lin, Variational problems with free interfaces, Calc. Var. Part. Diff. Eq. 1 (1993):149-168.

structure A. Novruzi, M. Pierre, Structure of shape derivatives. Journal of Evolution Equations 2 (2002): 365-382.

polya G. Pólya, Torsional rigidity, principal frequency, electrostatic capacity and symmetrization. Q. Appl. Math. 6 (1948): 267-277.

simon J. Simon, Second variations for domain optimization problems, International Series of Numerical Mathematics, vol. 91. Birkhauser, Basel (1989): 361-378. 










Research Center for Pure and Applied Mathematics, Graduate
School of
Information Sciences, Tohoku University, Sendai 980-8579
, Japan.
 

Electronic mail address:
cava@ims.is.tohoku.ac.jp




