# Turbulence and Random Geometry

Yaron Oz  Raymond and Beverly Sackler School of Physics and Astronomy, Tel
Aviv University, Tel Aviv 69978, Israel

###### Abstract

We outline our proposal for a field theory description of steady state
incompressible fluid turbulence at the inertial range of scales in a general
number of space dimensions. The theory consists of a Kolmogorov linear scaling
mean field theory dressed by a Nambu-Goldstone dilaton mode that induces a
random measure on the inertial range. We derive a KPZ-type formula for the
anomalous scalings of the velocity structure functions, the velocity gradients
and the local energy dissipation, and relate the dimensionless intermittency
parameter to the boundary conformal anomaly.

Dedicated to the memory of Jacob Bekenstein

##  I Introduction

Most fluid motions in our universe are turbulent. However, despite centuries
of research we still lack an analytical description and understanding of fluid
flows in the non-linear regime. Turbulence is considered as the most important
unsolved problem of classical physics. The turbulent incompressible fluid
flows exhibit highly complex spatial and temporal intermittent structures.
When considering their statistical average properties a universal structure is
revealed in the inertial range of scales that exhibits homogeneity, isotropy
and anomalous scalings. Our aim is to gain insight to the nature of the
anomalous scalings and construct an analytical formula to calculate them.

In the following we will outline our proposal for a field theory description
of steady state incompressible fluid turbulence at the inertial range of
scales. We will consider a general number of space dimensions $d\geq 2$, with
a particular attention being given to the physical two and three space
dimensions. Physically, our theory consists of a Kolmogorov linear scaling
mean field theory (K41) Kolmogorov  dressed by a Nambu-Goldstone (NG) gapless
dilaton mode. The NG mode arises from a spontaneous breaking due to the energy
flux of the separate scale and time symmetries of the inviscid Navier-Stokes
(NS) equations to a K41 scaling with a dynamical exponent $z=\frac{2}{3}$
Oz:2017ihc . Mathematically, the dilaton mode induces a random measure on the
inertial range and the emerging random geometry structure is reminiscent of
two-dimensional quantum gravity Knizhnik:1988ak ; Distler:1988jt , with the
dilaton being the analog of the Liouville mode.

We will derive a KPZ-type Knizhnik:1988ak ; Distler:1988jt  analytical formula
for the non-linear anomalous scalings of the velocity structure functions, the
velocity gradients and the local energy dissipation:

$\Delta-\Delta_{K41}={\cal G}^{2}(d)\Delta(1-\Delta)\ ,$ (1)

where $\Delta_{K41}$ is the K41 linear scaling dimension of the fluid
observable and $\Delta$ is the exact anomalous dimension. ${\cal G}(d)$ is a
dimensionless parameter that depends on the number of space dimensions and
quantifies intermittency and the deviation from Kolomogorov linear scaling. It
can be measured numerically and experimentally and we will relate it to the
boundary conformal anomaly of the K41 mean field theory Oz:2017ihc . Formula
(1) can be viewed as a physically and mathematically consistent completion
Eling:2015mxa  of the random cascade Kolmogorov-Obukhov lognormal model K62 ;
Obukhov  to the strong intermittency regime. It is in agreement with the
experimental and numerical data and provides an infinite number of new
predictions Eling:2015mxa .

The inspiration for the proposal comes from classical black hole dynamics.
Black holes have a horizon: objects fall into the black hole through the
horizon but cannot emerge. This breaks time reversal symmetry and allows
Einstein equations to describe dissipative effects (see e.g. Kovtun:2004de ;
Bhattacharyya:2008jc ) that are crucial for turbulence. The evolution
equations of the horizon are in the non-relativistic limit the NS equations
that describe incompressible fluid flows Damour ; Eling:2009pb . The fluid
pressure and velocity denote the deviation from equilibrium of the horizon
hypersurface location and its normal, respectively. Since the horizon is
characterized by a geometry, the gravity variables suggest that there should
be a geometrical framework for studying the dynamics of fluids: geometrization
of turbulence (for a brief review see Eling:2010vr ).

We can map every solution of the incompressible NS equations to a geometrical
configuration of the black hole horizon and its embedding. Hence, studying the
statistical properties of fluid flows is mapped to studying the statistics of
geometries, i.e. random geometry. By a random geometry we will mean changing
the Euclidean measure $dx$ on $R^{d}$ to a random measure $d\mu(x)\sim
e^{\gamma\phi(x)}dx$, where the Gaussian random field $\phi(x)$ is log-
correlated (for a review see e.g. chaos ). We will identify the field $\phi$
with the dilaton and $\gamma$ with the intermittency parameter ${\cal G}$
Oz:2017ihc .

The paper is organized as follows. In section II we will construct the field
theory of turbulence, while section III will be devoted to a discussion and
outlook.

##  II Field Theory of Turbulence

In this section we construct a field theory of turbulence in the inertial
range of scales. The theory consists of a Kolmogorov linear scaling theory
dressed by a NG dilaton that induces a random measure on the inertial range of
scales.

###  II.1 Flux States

The incompressible NS equations provide a universal description of fluid flows
at low Mach number, that is $v\ll v_{s}$ where $v_{s}$ is the speed of sound.
In flat space they read :

$\partial_{t}v^{i}+v^{j}\partial_{j}v^{i}=-\partial_{i}p+\nu\partial_{jj}v^{i}+F^{i},~{}~{}~{}~{}~{}~{}\partial_{i}v^{i}=0,~{}~{}~{}i=1,...,d\
,$ (2)

where $d\geq 2$ is the number of space dimensions, $v^{i},i=1...,d$ is the
fluid velocity vector, $p$ is the fluid pressure, $\nu$ is the kinematic
viscosity and $F^{i}$ is an external random force.

In a steady state we have an equality between the incoming energy flux due to
the external force and the energy dissipation due to viscosity:

$\int d^{d}x\left(F^{i}v^{i}\right)=\int d^{d}x\epsilon(x)\ ,$ (3)

where

$\epsilon(x)=\frac{\nu}{2}\left(\partial_{i}v^{j}+\partial_{j}v^{i}\right)^{2}\
,$ (4)

is the local energy dissipation.

Fluid flows exhibit turbulence at large Reynolds numbers ${\cal
R}_{e}=\frac{l_{c}v}{\nu}\simeq 10^{3}$, where $l_{c}$ is a characteristic
length scale and $v$ is the velocity difference at that scale. In such a state
the viscous term $\nu\partial_{jj}v^{i}$ is negligible compared to the
nonlinear term $v^{j}\partial_{j}v^{i}$ in the NS equations (2), yet,
$\epsilon(x)$ is non-vanishing and balances the incoming flux (3). This is a
turbulent flux state. It is far from equilibrium and the Gibbs measure is
inappropriate for quantifying its statistics. We will essentially work in the
infinite Reynolds number limit.

One defines the inertial range to be the range of distance scales $l\ll r\ll
L$, where the scales $l$ and $L$ are determined by the viscosity and forcing,
respectively. Although the turbulent flow breaks all the symmetries of the NS
equations, turbulence at the inertial range of scales exhibits statistical
homogeneity, isotropy and anomalous scalings. We will be interested in the
energy cascade steady state statistics and the equal time correlation
functions.

###  II.2 Spontaneous Symmetry Breaking of Space and Time Dilatations

In the absence of a viscosity term, the (inviscid) NS equations (2) exhibit
two scale symmetries. These are independent scalings of space and time which
we denote by $R_{x}\times R_{t}$:

$R_{x}:x^{i}\rightarrow e^{\sigma}x^{i}\ ,~{}~{}~{}~{}~{}R_{t}:t\rightarrow
e^{z\sigma}t\ ,$ (5)

where $\sigma$ and $z$ are real parameters and an appropriate charge is being
assigned to the random force. For any value of $z$, one can define the total
dimension of an object that carries dimensions $(r_{x},r_{t})$ under
$R_{x}\times R_{t}$ as

$\Delta_{z}=r_{x}+zr_{t}\ .$ (6)

The dimension of $v^{i}$ is $\Delta_{z}=1-z$ and the dimension of the
different terms in the inviscid NS equations is $\Delta_{z}=1-2z$.

The local energy dissipation (alternatively the flux) breaks spontaneously the
symmetries of the inviscid NS equations to $R_{z=\frac{2}{3}}$, as follows
from the requirement that

$\Delta_{z}(\epsilon)=2-3z=0\ .$ (7)

This is the Kolmogorov linear scaling
$\Delta_{z=\frac{2}{3}}\equiv\Delta_{K41}$:

$\Delta_{K41}[v^{i}]=\frac{1}{3}\ .$ (8)

Thus, we interpret the turbulent flux state as a state that spontaneously
breaks the independent scalings of space and time to $z=\frac{2}{3}$.

When a scale symmetry is spontaneously broken one expects a Nambu-Goldston
gapless mode called the dilaton. We denote the expectation value (VEV) that
breaks $R_{x}\times R_{t}\rightarrow R_{z=\frac{2}{3}}$ by:

$\langle\epsilon\rangle=\bar{\epsilon}\ ,$ (9)

where $\langle...\rangle$ means integration over both space and the
probability distribution function of the force $F$. When considering
turbulence in flat space $\bar{\epsilon}$ is a constant. It will be space
dependent when studying turbulent flows in curved spaces.

The dilaton $\tau(x)$ is the fluctuation:

$\epsilon(x)=\bar{\epsilon}e^{\delta\tau(x)}\ ,$ (10)

where $\delta(d)$ is a dimension dependent constant that will be related to
intermittency and $\delta=0$ means zero intermittency.

The dilaton is invariant under the unbroken $R_{z=\frac{2}{3}}$ symmetry
Hason:2017lao :

$\tau(x)\rightarrow\tau(x)+\alpha\sigma+\frac{2}{3}\beta\sigma,~{}~{}~{}~{}\alpha+\frac{2}{3}\beta=0\
,$ (11)

and $\Delta_{z=\frac{2}{3}}[e^{\tau}]=0$. We will denote its charges under
$R_{x}\times R_{t}$ as $(-Q,\frac{3}{2}Q)$. We will relate $Q$ to the
intermittency.

###  II.3 The Dilaton Field Theory

The steady state dilaton effective action in the inertial range of scales
should be invariant under the symmetry (5) and under Galilean boosts which
forbid time derivative terms in the Lagrangian. The dilaton action reads
Oz:2017ihc :

$S_{D}(\tau,\hat{g})=\frac{d}{\Omega_{d}(d-1)!}\int_{M}d^{d}x\sqrt{\hat{g}}\left(\tau{\cal
P}_{\hat{g}}\tau+2Q{\cal Q}_{\hat{g}}\tau\right)\ ,$ (12)

where we use the notation of Levy:2018bdc . When considering turbulence in
flat space, $M=B_{d}(L)$ is the $d$-dimensional closed ball in Euclidean space
$R^{d}$ with radius $L$, whose boundary is the $(d-1)$-sphere $S^{d-1}$ of
radius $L$ and $\hat{g}$ is the Euclidean metric. This will be the case most
interesting to us.

$Q$ is a dimensionless parameter called a background charge and $-Q$ is the
dilaton charge under scaling of space as defined above.
$\Omega_{d}=\frac{2\pi^{\frac{d+1}{2}}}{\Gamma[\frac{d+1}{2}]}$ is the surface
volume of the $d$-dimensional sphere $S^{d}$.

${\cal P}_{\hat{g}}$ are the conformally covariant GJMS operators GJMS :

$\mathcal{P}_{\hat{g}}=(-\Box)^{\frac{d}{2}}+lower~{}order\ ,$ (13)

where $\Box={\hat{g}}^{ab}\nabla_{a}\nabla_{b}$ with $\nabla_{a}$ being the
covariant derivative on $M$ and is the partial derivative on $M=B_{d}(L)$
whose metric is flat. They are local in even number of dimensions, for
instance $\mathcal{P}_{d=2}=-\Box$, and non-local in odd number dimensions odd
. The dilaton is a log-correlated field reminiscent of the Liouville field in
two-dimensional quantum gravity. $\mathcal{P}_{\hat{g}}$ gets a contribution
from the boundary that depends on the boundary conditions.

$\mathcal{Q}_{\hat{g}}$ is the ${\cal Q}$-curvature scalar Q :

${\cal Q}_{\hat{g}}=\frac{1}{2(d-1)}(-\Box)^{\frac{d}{2}-1}R+...\ .$ (14)

where $R$ is the scalar curvature with respect to ${\hat{g}}$. It is local in
even number of dimensions, for instance ${\cal Q}_{d=2}=\frac{1}{2}R$, and
non-local in odd number of dimensions odd . $\mathcal{Q}_{\hat{g}}$ gets a
contribution from the boundary, for instance ${\cal Q}_{d=2}=\frac{1}{2}R+K$
where $K$ is the extrinsic curvature. The bulk terms in
$\mathcal{Q}_{\hat{g}}$ vanish on the ball $B_{d}(L)$ and it is determined by
the boundary terms.

On conformally flat manifolds, such as the ball, ${\cal Q}_{\hat{g}}$ is
related to the Euler density $E_{d}$:

$\int_{M}d^{d}x\sqrt{{\hat{g}}}{\cal
Q}_{\hat{g}}=\frac{1}{2}\Omega_{d}(d-1)!\chi(M)\ ,$ (15)

where $\chi(M)$ is the Euler characteristic of $M$, $\chi(B_{d})=1$, and we
normalized the Euler density as:

$\int_{S^{d}}\sqrt{g}d^{d}xE_{d}=\Omega_{d}d!\ .$ (16)

Under a Weyl transformation ${\hat{g}}_{ab}\rightarrow
e^{2\sigma}{\hat{g}}_{ab}$:

${\cal P}_{e^{2\sigma}{\hat{g}}}=e^{-d\sigma}{\cal
P}_{\hat{g}},~{}~{}~{}~{}~{}{\cal Q}_{e^{2\sigma}{\hat{g}}}=e^{-d\sigma}({\cal
Q}_{\hat{g}}+{\cal P}_{\hat{g}}\sigma)\ ,$ (17)

and this relation holds both for even and odd number of dimensions. We have:

$S_{D}(\tau-Q\sigma,e^{2\sigma}{\hat{g}})=S_{D}(\tau,{\hat{g}})-S_{D}(Q\sigma,{\hat{g}})=S_{D}(\tau,{\hat{g}})-\frac{dQ^{2}}{\Omega_{d}(d-1)!}S_{L}(\sigma,{\hat{g}})\
,$ (18)

where the Liouville action reads:

$S_{L}(\phi,{\hat{g}})=\int_{M}d^{d}x\sqrt{{\hat{g}}}\left(\phi{\cal
P}_{{\hat{g}}}\phi+2{\cal Q}_{{\hat{g}}}\phi\right)\ ,$ (19)

and

$S_{L}(\phi-\phi_{0},e^{2\phi_{0}}{\hat{g}})=S_{L}(\phi,{\hat{g}})-S_{L}(\phi_{0},{\hat{g}})\
.$ (20)

The linear shift in the dilaton under Weyl transformation (18) means that it
can be viewed as a NG boson of a spontaneous breaking of the Weyl symmetry by
the choice of the fixed metric ${\hat{g}}$.

The dilaton action is changed by a constant under a constant shift:

$S_{D}(\tau+\tau_{0},{\hat{g}})=S_{D}(\tau,{\hat{g}})+dQ\chi(M)\tau_{0}\ ,$
(21)

which imposes a selection rule in correlation functions.

The semiclassical limit of the dilaton field theory (12) is the large charge
limit $Q\rightarrow\infty$ and the convenient parametrization is in terms of
the rescaled dilaton field $\tau_{c}=\frac{\tau}{Q}$. This limit will
correspond to the zero intermittency limit of the turbulence field theory.

Consider $M=B_{d}(L)$. The two-point function of the dilaton reads:

$\langle\tau(x)\tau(y)\rangle=\frac{1}{d}\log\left(\frac{L}{x-y}\right)+...\
,$ (22)

where dots refer to a constant and terms that vanish in the limit
$L\rightarrow\infty$. The dimension of the operator
$V_{\alpha}(x)=e^{d\alpha\tau(x)}$ reads:

$\Delta_{\alpha}=\frac{d}{2}\alpha(Q-\alpha)\ .$ (23)

The dilaton theory (12) is a generalization of the two-dimensional Liouville
conformal field theory to higher dimensions Levy:2018bdc .

###  II.4 Dilaton Dressing and Random Geometry

Define the random measure on the inertial range of scales $d\mu(x)$:

$d\mu(x)=\sqrt{g}d^{d}x,~{}~{}~{}~{}\sqrt{g}=\sqrt{\hat{g}}e^{d\gamma\tau}\ ,$
(24)

where $\sqrt{\hat{g}}d^{d}x$ is the fixed measure, which when considering
turbulence in flat space is the Euclidean measure as noted above. Physically
we identify this random measure with the local energy dissipation (10), with
$\sqrt{\hat{g}}=\bar{\epsilon}$ and $\delta=d\gamma$ which vanishes when
$\gamma\rightarrow 0$.

Requiring that $\Delta_{\gamma}=d$ (23):

$\frac{d}{2}\gamma(Q-\gamma)=d\ ,$ (25)

we get:

$Q=\gamma+\frac{2}{\gamma}\ ,$ (26)

and

$\gamma_{\pm}=\frac{Q\pm\sqrt{Q^{2}-8}}{2}\ .$ (27)

There are two branches $\gamma_{\pm}$. We require $\gamma$ to be real, hence
$Q^{2}\geq 8$. Without loss of generality we take $Q\geq 0$. We have:

$0\leq\gamma_{-}^{2}\leq 2,~{}~{}~{}~{}~{}\gamma_{+}^{2}\geq
2,~{}~{}~{}~{}~{}\gamma_{-}\gamma_{+}=2\ .$ (28)

The critical value $\gamma=\sqrt{2}$ corresponds in two-dimensional Liouville
quantum gravity to the $c=1$ barrier Knizhnik:1988ak ; Distler:1988jt .

The operators in the theory are K41 operators $O_{K41}$ dressed by a dilaton
factor $e^{d\alpha\tau}$:

$O(x)=e^{d\alpha\tau}O_{K41}(x),~{}~{}~{}~{}~{}\alpha=\gamma(1-\Delta)\ ,$
(29)

where $d\Delta_{K41}$ is the undressed dimension of $O_{K41}$. Note, that both
$O_{K41}$ and its dilaton factor dressing are made of the fluid velocity and
pressure. Requiring the conformal weight of $O$ to be $d$ we get:

$\frac{d}{2}\alpha(Q-\alpha)+d\Delta_{K41}=d\ ,$ (30)

which using (26) and (29) reads:

$\Delta-\Delta_{K41}=\frac{\gamma^{2}}{2}\Delta(1-\Delta)\ .$ (31)

The partition of the theory in a fixed volume $V={\cal L}^{d}$ with respect to
the random metric scales as ${\cal L}^{-2d(1+\frac{1}{\gamma^{2}})}$. The
expectation value of the one-point function of the integrated dressed operator
(29) in a fixed volume scales as ${\cal L}^{d(1-\Delta)}$. $d(1-\Delta)$ is
the dimension of a set in the random geometry of the inertial range and
$d(1-\Delta_{K41})$ is the dimension with respect to the fixed Euclidean
metric on it. In general, the expectation value of the $n$-point function in a
fixed volume scales as ${\cal L}^{d\sum_{i=1}^{n}(1-\Delta_{i})}$.

Note, that in the field theory framework the natural units for the dimensions
(conformal weights) of operators are inverse length, while when we will be
applying (31) to turbulence we will associate dimensions of length to the
fluid observables. In the fluid setup there is an IR scale $L$, which can be
used in order to relate the two.

###  II.5 Anomalous Scaling

Formula (31) reads as (1) with ${\cal G}^{2}(d)=\frac{{\gamma}^{2}}{2}$ that
quantifies intermittency and the deviation from Kolomogorov linear scaling.
The limit of zero intermittency ${\cal G}^{2}\rightarrow 0$ is the
semiclassical limit of the dilaton theory $Q\rightarrow\infty$ with ${\cal
G}^{2}\sim\frac{2}{Q^{2}}$. $\Delta_{K41}$ is the K41 dimension of the
operator and $\Delta$ is the exact anomalous dimension. We have:

$0\leq{\cal G}^{2}_{-}\leq 1,~{}~{}~{}~{}~{}{\cal G}^{2}_{+}\geq
1,~{}~{}~{}~{}~{}{\cal G}_{+}{\cal G}_{-}=1\ .$ (32)

We will assume that the formula can be used when $\Delta_{K41}>0$.

Consider the longitudinal structure functions $S_{n}$:

$S_{n}(r)=\langle(\delta_{r}v)^{n}\rangle\sim r^{\xi_{n}}\ ,$ (33)

where $\delta_{r}v$ is the longitudinal velocity difference between points
separated by a fixed distance $r=|\vec{r}|$:

$\delta_{r}v=\left(\vec{v}(\vec{r},t)-\vec{v}(0,t)\right)\cdot\frac{\vec{r}}{r}\
.$ (34)

The K41 scaling dimension of $(\delta_{r}v)^{n}$ is
$\Delta_{K41}=\frac{n}{3}$, thus

$\xi_{n}-\frac{n}{3}={\cal G}^{2}(d)\xi_{n}(1-\xi_{n})\ .$ (35)

The numerical and experimental data in two (inverse cascade) and three space
dimensions suggest Eling:2015mxa :

${\cal G}^{2}_{d=2}\sim 0,~{}~{}~{}~{}{\cal G}^{2}_{d=3}\sim 0.16\ ,$ (36)

and implies that we should pick the $-$ branch (32). The corresponding charges
are $Q^{2}_{d=3}\sim 17$, and if we take ${\cal G}^{2}_{d=2}\sim 0.01$ we get
$Q^{2}_{d=2}\sim 205$. Both may be considered much larger than one, thus
intermittency in these cases should be explained to a good accuracy by the
semiclassical dynamics of the dilaton theory.

Formula (35) predicts a generic growth of the anomalous scalings at large $n$:

$\xi_{n}\sim\sqrt{n},~{}~{}~{}~{}n\gg 1\ .$ (37)

This is a rather different behaviour compared to other existing proposals for
the anomalous scalings. Unfortunately, at this point an accurate experimental
or numerical data of high enough moments is not available and we cannot check
this prediction.

The anomalous dimensions of the transverse structure functions of the velocity
vector field as well as the mixed cases is the same as that of the
longitudinal structure functions since the transverse and the longitudinal
velocity increments have the same K41 scaling dimension.

The normalized local space average of the energy dissipation (4) over a
$d$-dimensional ball of radius $r$, $B_{d}(r)$, centered at a point $x$ reads:

$\epsilon_{r}(x)=\frac{1}{Vol(B_{d}(r))}\int_{|x^{\prime}-x|\leq
r}d^{d}x^{\prime}\epsilon(x^{\prime})\ .$ (38)

Consider the quantity $(r\epsilon_{r})^{\frac{1}{3}}$. Its K41 scaling
dimension is $\frac{1}{3}$. Formula (1) suggests that its statistics is the
same as that of $\delta_{r}v$. Thus, if

$\langle\epsilon_{r}^{n}\rangle\sim r^{\zeta_{n}}\ ,$ (39)

we get:

$\zeta_{\frac{n}{3}}=\xi_{n}-\frac{n}{3}\ ,$ (40)

with $\xi_{n}$ given in (35).

Next, we define the quantity $D_{r}^{k}\delta_{r}v$ where
$D_{r}^{k}=(r\partial_{r})^{k}$. It has the same K41 scaling dimension as
$\delta_{r}v$, thus the anomalous dimensions of these operators are:

$\langle(D_{r}^{k}\delta_{r}v)^{n}\rangle\sim r^{\xi_{n}}\ .$ (41)

###  II.6 Comparison to Random Cascade Models

Random cascade models for the redistribution of fluid energy via splitting
into finer scale eddies allows a rich phenomenological framework for studying
fluid turbulence and intermittency (for a review see e.g. Frisch ; chaos ).
When the random variable is lognormally distributed one gets the Kolmogorov-
Obukhov lognormal model K62 ; Obukhov .

At leading order in the intermittency parameter ${\cal G}^{2}$, $\xi_{n}$ of
(35) coincides with the lognormal model:

$\displaystyle\xi_{n}-\frac{n}{3}={\cal
G}^{2}\frac{n}{3}\left(1-\frac{n}{3}\right)\ .$ (42)

They differ at higher orders in ${\cal G}^{2}$ and (35) may be viewed as a
completion of the Kolmogorov-Obukhov formula (42) to the strong intermittency
regime.

Unlike (42) which implies physically inconsistent supersonic velocities at
large $n$ and a violation of the convexity inequality, its completion (35) is
mathematically and physically consistent.

In the lognormal model, the moments of the local energy dissipation
$\epsilon_{r}$ scale as:

$\langle\epsilon_{r}^{n}\rangle\sim r^{{\cal G}^{2}n(1-n)}\ ,$ (43)

and

$\langle(\delta_{r}v)^{n}\rangle\sim\langle(r\epsilon_{r})^{\frac{n}{3}}\rangle\sim\langle\epsilon_{r}^{\frac{n}{3}}\rangle
r^{\frac{n}{3}}\sim r^{{\cal G}^{2}\frac{n}{3}(1-\frac{n}{3})+\frac{n}{3}}\ ,$
(44)

as in (42). In the random cascade models one separates the dynamics of the K41
theory from the that of the random measure and this is the reason for the
factorization in (44). However, this separation of the dynamics is only an
approximation that can be made at weak intermittency. In our language, the
dilaton field interacts with the K41 theory and we cannot seperate their
dynamics.

More generally, in the random cascade models with a random variable $W$ one
has:

$\langle\epsilon_{r}^{n}\rangle\sim
r^{\zeta_{n}},~{}~{}~{}~{}\zeta_{n}=-\log_{2}(W^{n})\ ,$ (45)

which coincides with (43) when $W$ is lognormally distributed, and

$\xi_{n}=\frac{n}{3}-\log_{2}(W^{\frac{n}{3}})\ .$ (46)

In the random geometry framework we get:

$\xi_{n}=\frac{n}{3}-\log_{2}(W^{\xi_{n}})\ ,$ (47)

which coincides with (35) when $W$ is lognormally distributed.

###  II.7 Conformal Anomaly and Intermittency

On the ball $B_{d}(L)$ the conformal anomaly coefficient $a$ is obtained by
varying the field theory partition function $Z_{B_{d}}$:

$a\sim\frac{\partial}{\partial\log L}Z_{B_{d}}(L)\ .$ (48)

Another way to see the anomaly is by the non-invariance of the partition
function under a Weyl transformation that results in a term of the form
$e^{-aS_{L}}$, where $S_{L}$ is the Liouville action (19) and we neglect terms
that are suppressed in the inertial range Oz:2017ihc . We will require that
the total conformal anomaly vanishes:

$a_{total}=a_{dilaton}+a_{K41}=0\ .$ (49)

We may justify this requirement by asking that the inertial range universal
structure and in particular the anomalous scalings should not depend on $L$.
In an even number of dimensions we have Dowker:2010bu ; Oz:2017ihc ;
Levy:2018bdc :

$a_{dilaton}=\frac{2}{\Omega_{d}(d!)^{2}}\int_{0}^{\frac{d}{2}}dt\prod_{i=0}^{\frac{d}{2}-1}\left(i^{2}-t^{2}\right)+\frac{(-1)^{\frac{d}{2}}}{\Omega_{d}(d-1)!}Q^{2}\
.$ (50)

Since $Q^{2}\geq 8$ we see that the first term in (50) can be neglected and

$a_{dilaton}\simeq\frac{(-1)^{\frac{d}{2}}}{\Omega_{d}(d-1)!}Q^{2}\ ,$ (51)

that is its classical value (18). Using (49) and the large $Q$ behaviour of
the intermittency parameter ${\cal G}^{2}\sim\frac{2}{Q^{2}}$ we get:

${\cal G}^{2}(d)\simeq\frac{2}{\Omega_{d}(d-1)!|a_{K41}(d)|}\ .$ (52)

This relation holds also for odd number of dimensions, where in both cases the
contribution to the anomaly comes from the boundary Graham:1999pm ; Oz:2017ihc
.

##  III Discussion and Outlook

We outlined our proposal for a field theory description of the universal
structure of incompressible fluid turbulence at the inertial range of scales,
derived an analytical formula for anomalous scalings and related the
intermittency parameter to the conformal anomaly coefficient. Note, that if
the K41 theory is not conformally invariant Oz:2018yaz  the same analysis as
above holds but $a_{K41}$ (52) gets contributions on the ball not only from
the Euler density term. Mathematically, the structure of the field theory is
that of a dynamical system (K41 theory) coupled to a random geometry.
Physically, the random geometry structure is induced by a NG dilaton that is a
consequence of a spontaneous symmetry breaking in the turbulent flux state and
the dilaton field is responsible for the intermittency. There are infinite
number of anomalous scaling dimensions predicted by the theory that may be
checked numerically or experimentally.

There are numerous directions to follow. A calculation of the anomaly
coefficient of the K41 mean field theory $a_{K41}$ will enable via the
relation (52) to calculate the intermittency parameter and compare with
numerical and experimental data. For this we need to gain a better
understanding of K41 theory, which is a non-local field theory. A calculation
of equal time correlations of the local energy dissipation such as
$\langle(\epsilon_{r}(x))^{p}(\epsilon_{R}(x))^{q}\rangle$,
$\langle(\epsilon_{r}(x_{1}))^{p}(\epsilon_{r}(x_{2}))^{q}\rangle$ and
velocity differences $\langle(\delta_{r}v)^{p}(\delta_{R}v)^{q}\rangle$ (see
e.g. biferale ), where $r\neq R$ are in the inertial range, can shed more
light on the theory and may be compared with numerical and experimental data.

The random cascade phenomenological framework for studying fluid turbulence
and intermittency is not directly related to the NS equations that govern the
motion of fluid flows. The spontaneous symmetry breaking mechanism that we
introduced provides a direct link between the two via the dilaton that
introduces the random structure. This topic deserves further study, in
particular since the spontaneous symmetry breaking is in the flux state which
is far from equilibrium.

It has been stressed by Polyakov Polyakov:1992er  that the IR boundary
conditions at the forcing scale $L$ are of much importance for turbulence and
that IR condensates break spontaneously conformal symmetry. Our dilaton theory
is a conformal field theory Levy:2018bdc  and it would be valuable to
understand the IR effects on it as well as on the K41 mean field theory. In
our framework the natural units for dimensions are inverse length and an
explicit introduction of the IR scale $L$ can explain the change to units of
length.

A study of the non-local dilaton theories in odd number of dimensions may shed
more light on odd-dimensional turbulence in general and on three-dimensional
turbulence in particular. Also, studying turbulence in curved space may offer
a new viewpoint on its dynamics. In our framework it amounts to introducing a
random dressing of a fixed curved metric on the inertial range.

Higher-dimensional turbulence is mathematically intriguing and may also
provide a $\frac{1}{d}$ expansion scheme. Denote by $u$ the velocity one-form,
then $du$ is the vorticity two-form and the incompressibility condition reads
$d^{*}u=0$. In addition to the energy one has a conserved helicity made from
the wedge products ${\cal H}=\int_{M}(du)^{m}$ when ${\rm dim}M=2m$ is even
and ${\cal H}=\int_{M}u\wedge(du)^{m}$ when $dimM=2m+1$ is odd arnold . The
well defined parity of the helicity in even dimensions suggests that as in two
space dimensions there may be two cascades in these cases: an inverse energy
cascade and a direct helicity cascade.

One may try to generalize our framework to other systems. An interesting case
is Burgers turbulence, which will correspond to a one-dimensional random
geometry. One may expect that the shock wave solutions of Burges equations
that determine the anomalous scaling correspond to singular one-dimensional
manifolds.

Finally, the structure of our proposed field theory of steady state inertial
range turbulence is reminiscent of the macroscopic structure of quantum
gravity as a sum over spacetime geometries. It will be interesting to obtain a
better understanding of this relationship.

## Acknowledgements

I would like to thank L. Biferale, V. Mukhanov and A. Polyakov for
discussions. This work is supported in part by the I-CORE program of Planning
and Budgeting Committee (grant number 1937/12), the US-Israel Binational
Science Foundation, GIF and the ISF Center of Excellence.

## References

  * (1) A. N. Kolmogorov, ”The local structure of turbulence in incompressible viscous fluid for very large Reynolds number”, Dokl. Akad. Nauk. SSSR 30, 299 (1941) and ”Dissipation of energy in in locally isotroptic turbulence” 32, 16 (1941). Reproduced in Proc. R. Soc. London, Ser. A 434 9 and 15 (1991).
  * (2) Y. Oz, “Spontaneous Symmetry Breaking, Conformal Anomaly and Incompressible Fluid Turbulence,” JHEP 1711, 040 (2017) [arXiv:1707.07855 [hep-th]].
  * (3) V. G. Knizhnik, A. M. Polyakov and A. B. Zamolodchikov, “Fractal Structure of 2D Quantum Gravity,” Mod. Phys. Lett. A 3, 819 (1988).
  * (4) J. Distler and H. Kawai, “Conformal Field Theory and 2D Quantum Gravity,” Nucl. Phys. B 321, 509 (1989).
  * (5) C. Eling and Y. Oz, “The Anomalous Scaling Exponents of Turbulence in General Dimension from Random Geometry,” JHEP 1509, 150 (2015) [arXiv:1502.03069 [nlin.CD]].
  * (6) A. N. Kolmogorov, “A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number,” J. Fluid Mech. 13 (1962) 82.
  * (7) A. M. Obukhov,“Some specific features of atmospheric turbulence,” J. Fluid Mech. 13 77 (1962).
  * (8) P. Kovtun, D. T. Son and A. O. Starinets, “Viscosity in strongly interacting quantum field theories from black hole physics,” Phys. Rev. Lett. 94, 111601 (2005) doi:10.1103/PhysRevLett.94.111601 [hep-th/0405231].
  * (9) S. Bhattacharyya, V. E. Hubeny, S. Minwalla and M. Rangamani, “Nonlinear Fluid Dynamics from Gravity,” JHEP 0802, 045 (2008) [arXiv:0712.2456 [hep-th]].
  * (10) T. Damour, ”Surface Effects in Black-Hole Physics”, Proceedings of the Second Marcel Grossmann Meeting on General Relativity, 587 (1982).
  * (11) C. Eling, I. Fouxon and Y. Oz, “The Incompressible Navier-Stokes Equations From Membrane Dynamics,” Phys. Lett. B 680, 496 (2009) [arXiv:0905.3638 [hep-th]].
  * (12) C. Eling, I. Fouxon and Y. Oz, “Gravity and a Geometrization of Turbulence: An Intriguing Correspondence,” Contemporary Physics 52, 43 (2011), [arXiv:1004.2632 [hep-th]].
  * (13) R. Rhodes and V. Vargas,”Gaussian multiplicative chaos and applications: A review”, Probab. Survey 11, 315 (2014), arXiv:1305.6221 [math.PR].
  * (14) I. Hason, “Dual Scale Invariance Spontaneous Symmetry Breaking and Turbulence,” arXiv:1708.08294 [hep-th].
  * (15) T. Levy and Y. Oz, “Liouville Conformal Field Theories in Higher Dimensions,” JHEP 1806, 119 (2018) [arXiv:1804.02283 [hep-th]].
  * (16) C. R. Graham, R. Jenne, L. J. Mason, G. A. J. Sparling, ”Conformally Invariant Powers of the Laplacian, I: Existence”, J. London Math. Soc. (2) 46, 557 (1992).
  * (17) C. Fefferman and C. R. Graham,”Q-Curvature and Poincare Metrics”, arXiv:math/0110271.
  * (18) T. Branson, ”Differential operators canonically associated to a conformal structure”, Math. Scand. 57 (1985) 293.
  * (19) U. Frisch, Turbulence: The Legacy of A. N. Kolmogorov, Cambridge University Press (1995).
  * (20) S. Deser and A. Schwimmer, “Geometric classification of conformal anomalies in arbitrary dimensions,” Phys. Lett. B 309, 279 (1993) [hep-th/9302047].
  * (21) J. S. Dowker, “Entanglement entropy for even spheres,” arXiv:1009.3854 [hep-th].
  * (22) C. R. Graham and E. Witten, “Conformal anomaly of submanifold observables in AdS / CFT correspondence,” Nucl. Phys. B 546, 52 (1999) [hep-th/9901021].
  * (23) Y. Oz, “On Scale Versus Conformal Symmetry in Turbulence,” Eur. Phys. J. C 78, no. 8, 655 (2018) doi:10.1140/epjc/s10052-018-6147-8 [arXiv:1801.04388 [hep-th]].
  * (24) R. Benzi, L. Biferale, and F. Toschi, ”Multiscale Velocity Correlations in Turbulence”, Phys. Rev. Lett. 80, 3244 (1998).
  * (25) A. M. Polyakov, “The Theory of turbulence in two-dimensions,” Nucl. Phys. B 396, 367 (1993) doi:10.1016/0550-3213(93)90656-A [hep-th/9212145].
  * (26) V. Arnold, Vladimir and B. Khesin, ”Topological Methods in Hydrodynamics”, Springer (1998).


term of the resulting equation separately. Analogously to (6.20) and (6.21),
we obtain

$\displaystyle\int_{-\infty}^{\infty}\boldsymbol{\Lambda}_{1}\cdot\partial_{z}\boldsymbol{\Phi}_{0}\mathrm{\,d}z$
$\displaystyle=0\quad\text{and}$ (6.32)
$\displaystyle\int_{-\infty}^{\infty}\boldsymbol{\vartheta}_{1}\cdot\partial_{z}\boldsymbol{\Phi}_{0}\mathrm{\,d}z$
$\displaystyle=\boldsymbol{\vartheta}_{1}\cdot\left(\boldsymbol{e}_{j}-\boldsymbol{e}_{i}\right).$
(6.33)

Considering the Lagrange multiplier $\boldsymbol{\mu}$, we recall (6.23)

$\displaystyle\mu_{i}^{\varepsilon}\,\nabla_{x}\varphi_{i}^{\varepsilon}=\boldsymbol{0}.$
(6.34)

Due to (4.17), its contribution of leading order $\mathcal{O}(1)$ in inner
coordinates is given by

$\displaystyle\mu_{1}^{i}\partial_{z}\Phi_{0}^{i}\boldsymbol{\nu}+\mu_{0}^{i}\nabla_{\Gamma}\Phi_{0}^{i}+\mu_{0}^{i}\partial_{z}\Phi_{1}^{i}\boldsymbol{\nu}=\boldsymbol{0}$

for all $i\in\\{1,\dots,N\\}$. Multiplying this identity by $\boldsymbol{\nu}$
and integrating the resulting equation with respect to $z$, we infer

$\displaystyle\int_{-\infty}^{\infty}\boldsymbol{\mu}_{1}\cdot\partial_{z}\boldsymbol{\Phi}_{0}\mathrm{\,d}z$
$\displaystyle=-\int_{-\infty}^{\infty}\boldsymbol{\mu}_{0}\cdot\partial_{z}\boldsymbol{\Phi}_{1}\mathrm{\,d}z.$
(6.35)

Furthermore, applying integration by parts twice and using that due to the
matching condition (5.2) all derivatives of $\boldsymbol{\Phi}_{0}$ with
respect to $z$ tend to $0$ as $z\to\pm\infty$, we obtain

$\displaystyle\int_{-\infty}^{\infty}\partial_{zz}\boldsymbol{\Phi}_{1}\cdot\partial_{z}\boldsymbol{\Phi}_{0}\mathrm{\,d}z=\int_{-\infty}^{\infty}\partial_{zz}\left(\partial_{z}\boldsymbol{\Phi}_{0}\right)\cdot\boldsymbol{\Phi}_{1}\mathrm{\,d}z.$
(6.36)

As $\partial_{z}\boldsymbol{\Phi}_{0}$ attains its values only in
$T\Sigma^{N}$, we deduce

$\displaystyle\int_{-\infty}^{\infty}P_{T\Sigma}\left[\psi_{0}^{\prime\prime}(\boldsymbol{\Phi}_{0})\boldsymbol{\Phi}_{1}\right]\cdot\partial_{z}\boldsymbol{\Phi}_{0}\mathrm{\,d}z$
$\displaystyle=\int_{-\infty}^{\infty}\psi_{0}^{\prime\prime}(\boldsymbol{\Phi}_{0})\,\partial_{z}\boldsymbol{\Phi}_{0}\cdot\boldsymbol{\Phi}_{1}\mathrm{\,d}z$
(6.37)

due to the symmetry of the Hessian matrix. Moreover, recalling that
$\mathbf{W}_{0}$ is independent of $z$ due to (6.3), a simple computation
yields

$\displaystyle\begin{aligned}
\int_{-\infty}^{\infty}\overline{\rho}^{\prime}(\boldsymbol{\Phi}_{0})\partial_{z}\boldsymbol{\Phi}_{0}\left|\mathbf{W}_{0}\right|^{2}\mathrm{\,d}z&=\int_{-\infty}^{\infty}\left[\frac{\text{d}}{\text{d}z}\overline{\rho}(\boldsymbol{\Phi}_{0})\right]\left|\mathbf{W}_{0}\right|^{2}\mathrm{\,d}z\\\
&=\int_{-\infty}^{\infty}\frac{\text{d}}{\text{d}z}\left[\overline{\rho}(\boldsymbol{\Phi}_{0})\left|\mathbf{W}_{0}\right|^{2}\right]\mathrm{\,d}z.\end{aligned}$
(6.38)

Furthermore, by the definition of the dyadic product, it holds

$\displaystyle\left(\nabla_{\Gamma}\mathbf{W}_{0}+\partial_{z}\mathbf{W}_{1}\otimes\boldsymbol{\nu}\right)^{\text{sym}}\boldsymbol{\nu}\cdot\partial_{zz}\mathbf{W}_{1}$
$\displaystyle\quad=\left(\nabla_{\Gamma}\mathbf{W}_{0}+\partial_{z}\mathbf{W}_{1}\otimes\boldsymbol{\nu}\right)^{\text{sym}}:\left(\partial_{zz}\mathbf{W}_{1}\otimes\boldsymbol{\nu}\right)^{\text{sym}}.$

Now we use (6.3) (which directly entails
$\partial_{z}\nabla_{\Gamma}\mathbf{W}_{0}=\boldsymbol{0}$), (6.4) and
$\partial_{z}\boldsymbol{\nu}=\boldsymbol{0}$ to deduce

$\displaystyle\begin{aligned}
&\int_{-\infty}^{\infty}\overline{\mathbb{C}}^{\prime}(\boldsymbol{\Phi}_{0})\partial_{z}\boldsymbol{\Phi}_{0}\left(\nabla_{\Gamma}\mathbf{W}_{0}+\partial_{z}\mathbf{W}_{1}\otimes\boldsymbol{\nu}\right)^{\text{sym}}:\left(\nabla_{\Gamma}\mathbf{W}_{0}+\partial_{z}\mathbf{W}_{1}\otimes\boldsymbol{\nu}\right)^{\text{sym}}\mathrm{\,d}z\\\
&\quad=\int_{-\infty}^{\infty}\left[\frac{\text{d}}{\text{d}z}\overline{\mathbb{C}}(\boldsymbol{\Phi}_{0})\right]\left(\nabla_{\Gamma}\mathbf{W}_{0}+\partial_{z}\mathbf{W}_{1}\otimes\boldsymbol{\nu}\right)^{\text{sym}}:\left(\nabla_{\Gamma}\mathbf{W}_{0}+\partial_{z}\mathbf{W}_{1}\otimes\boldsymbol{\nu}\right)^{\text{sym}}\mathrm{\,d}z\\\
&\quad=\int_{-\infty}^{\infty}\frac{\text{d}}{\text{d}z}\left[\overline{\mathbb{C}}(\boldsymbol{\Phi}_{0})\left(\nabla_{\Gamma}\mathbf{W}_{0}+\partial_{z}\mathbf{W}_{1}\otimes\boldsymbol{\nu}\right)^{\text{sym}}:\left(\nabla_{\Gamma}\mathbf{W}_{0}+\partial_{z}\mathbf{W}_{1}\otimes\boldsymbol{\nu}\right)^{\text{sym}}\right]\mathrm{\,d}z\\\
&\qquad-2\int_{-\infty}^{\infty}\frac{\text{d}}{\text{d}z}\left[\overline{\mathbb{C}}(\boldsymbol{\Phi}_{0})\left(\nabla_{\Gamma}\mathbf{W}_{0}+\partial_{z}\mathbf{W}_{1}\otimes\boldsymbol{\nu}\right)^{\text{sym}}\boldsymbol{\nu}\cdot\partial_{z}\mathbf{W}_{1}\right]\mathrm{\,d}z\end{aligned}$
(6.39)

by means of the product rule and integration by parts.

Collecting (6.32)–(6.39) and recalling (6.30), we eventually obtain

$\displaystyle\begin{aligned}
&\boldsymbol{\vartheta}_{1}\cdot\left(\boldsymbol{e}_{j}-\boldsymbol{e}_{i}\right)\\\
&\quad+\int_{-\infty}^{\infty}\left(\partial_{zz}\left(\partial_{z}\boldsymbol{\Phi}_{0}\right)-\psi_{0}^{\prime\prime}\left(\boldsymbol{\Phi}_{0}\right)\partial_{z}\boldsymbol{\Phi}_{0}\right)\cdot\boldsymbol{\Phi}_{1}\mathrm{\,d}z-\int_{-\infty}^{\infty}\boldsymbol{\mu}_{0}\cdot\partial_{z}\boldsymbol{\Phi}_{1}\mathrm{\,d}z\\\
&=-\frac{1}{\gamma}\sum_{r=1}^{l}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\lambda_{0,n_{r}}\left(\int_{-\infty}^{\infty}\frac{\text{d}}{\text{d}z}\left(\overline{\rho}(\boldsymbol{\Phi}_{0})\left|\mathbf{W}_{0,n_{r}}\right|^{2}\right)\mathrm{\,d}z\right)\\\
&\quad+\sigma_{ij}\hat{\kappa}+\frac{1}{\gamma}\sum_{r=1}^{l}\Bigg{\\{}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\\\
&\hskip
64.58313pt\cdot\Bigg{[}\int_{-\infty}^{\infty}\frac{\text{d}}{\text{d}z}\Big{(}\overline{\mathbb{C}}(\boldsymbol{\Phi}_{0})(...)^{\text{sym}}:(...)^{\text{sym}}\Big{)}\mathrm{\,d}z\\\
&\hskip
64.58313pt\quad-2\int_{-\infty}^{\infty}\frac{\text{d}}{\text{d}z}\left(\overline{\mathbb{C}}(\boldsymbol{\Phi}_{0})(...)^{\text{sym}}\boldsymbol{\nu}\cdot\partial_{z}\mathbf{W}_{1,n_{r}}\right)\mathrm{\,d}z\Bigg{]}\Bigg{\\}}\end{aligned}$
(6.40)

on $\Gamma_{ij}$, where $(...)^{\text{sym}}$ abbreviates
$\big{(}\nabla_{\Gamma}\mathbf{W}_{0,n_{r}}+\partial_{z}\mathbf{W}_{1,n_{r}}\otimes\boldsymbol{\nu}\big{)}^{\text{sym}}$.
Next, we want to show that

$\displaystyle\int_{-\infty}^{\infty}\big{(}\partial_{zz}(\partial_{z}\boldsymbol{\Phi}_{0})-\psi_{0}^{\prime\prime}(\boldsymbol{\Phi}_{0})\partial_{z}\boldsymbol{\Phi}_{0}\big{)}\cdot\boldsymbol{\Phi}_{1}\mathrm{\,d}z-\int_{-\infty}^{\infty}\boldsymbol{\mu}_{0}\cdot\partial_{z}\boldsymbol{\Phi}_{1}\mathrm{\,d}z=0.$
(6.41)

Differentiating (6.25) with respect to $z$, multiplying by
$\boldsymbol{\Phi}_{1}$ and integrating the resulting equation with respect to
$z$, we deduce

$\displaystyle\int_{-\infty}^{\infty}\big{(}\partial_{zz}(\partial_{z}\boldsymbol{\Phi}_{0})-\psi_{0}^{\prime\prime}(\boldsymbol{\Phi}_{0})\partial_{z}\boldsymbol{\Phi}_{0}\big{)}\cdot\boldsymbol{\Phi}_{1}\mathrm{\,d}z=-\int_{-\infty}^{\infty}\left[\partial_{z}(\boldsymbol{\Lambda}_{0}+\boldsymbol{\mu}_{0})\right]\cdot\boldsymbol{\Phi}_{1}\mathrm{\,d}z.$

Thus, in order to prove (6.41), it suffices to show

$\displaystyle\int_{-\infty}^{\infty}\left[\partial_{z}(\boldsymbol{\Lambda}_{0}+\boldsymbol{\mu}_{0})\right]\cdot\boldsymbol{\Phi}_{1}\mathrm{\,d}z+\int_{-\infty}^{\infty}\boldsymbol{\mu}_{0}\cdot\partial_{z}\boldsymbol{\Phi}_{1}\mathrm{\,d}z=0.$
(6.42)

By means of the product rule, the left-hand side can be reformulated as

$\displaystyle\int_{-\infty}^{\infty}\left[\partial_{z}(\boldsymbol{\Lambda}_{0}+\boldsymbol{\mu}_{0})\right]\cdot\boldsymbol{\Phi}_{1}\mathrm{\,d}z$
$\displaystyle=-\int_{-\infty}^{\infty}(\boldsymbol{\Lambda}_{0}+\boldsymbol{\mu}_{0})\cdot\partial_{z}\boldsymbol{\Phi}_{1}\mathrm{\,d}z$
$\displaystyle\qquad+\int_{-\infty}^{\infty}\frac{\text{d}}{\text{d}z}\left[(\boldsymbol{\Lambda}_{0}+\boldsymbol{\mu}_{0})\cdot\boldsymbol{\Phi}_{1}\right]\mathrm{\,d}z.$

Now as $\boldsymbol{\Phi}_{1},\partial_{z}\boldsymbol{\Phi}_{1}\in
T\Sigma^{N}$ pointwise, we know as in (6.20)

$\displaystyle\boldsymbol{\Lambda}_{0}\cdot\partial_{z}\boldsymbol{\Phi}_{1}=\boldsymbol{\Lambda}_{0}\cdot\boldsymbol{\Phi}_{1}=0.$

Thus, as we want to prove (6.42), it remains to show

$\displaystyle\int_{-\infty}^{+\infty}\frac{\text{d}}{\text{d}z}\left[\boldsymbol{\mu}_{0}\cdot\boldsymbol{\Phi}_{1}\right]\mathrm{\,d}z=0.$
(6.43)

Recalling once more formula (3.27), we infer

$\displaystyle[\mu^{\varepsilon}]^{i}[\varphi^{\varepsilon}]^{i}=0$

for all $i\in\\{1,\dots,N\\}$. Hence, for the $\mathcal{O}(1)$-contribution
and the $\mathcal{O}(\varepsilon)$-contribution of the inner expansions, we
obtain the relations

$\displaystyle\mu_{0}^{i}(z)\,\Phi^{i}_{0}(z)=0\quad\text{and}\quad\mu_{1}^{i}(z)\,\Phi^{i}_{0}(z)=-\mu_{0}^{i}(z)\,\Phi^{i}_{1}(z),$
(6.44)

respectively, for all $i\in\\{1,\dots,N\\}$ and $z\in\mathbb{R}$. Now, the
first equation in (6.44) implies that for any $z\in\mathbb{R}$ with
$\Phi_{0}^{i}(z)\neq 0$, we have $\mu_{0}^{i}(z)=0$ and thus also
$\mu_{0}^{i}(z)\Phi_{1}^{i}(z)=0$. On the other, for all $z\in\mathbb{R}$ with
$\Phi_{0}^{i}(z)=0$, we infer from the second equation in (6.44) that
$\mu_{0}^{i}(z)\Phi_{1}^{i}(z)=0$. Combining both statements, we conclude

$\displaystyle\mu_{0}^{i}(z)\,\Phi_{1}^{i}(z)=0\quad\text{for all
$i\in\\{1,\dots,N\\}$ and $z\in\mathbb{R}$.}$

This proves (6.43). By the above considerations, this verifies (6.42) which in
turn implies equation (6.41).

To conclude this section, we recall the definition of the jump, see (5.5).
Moreover, we recall from (4.23) that the mean curvature of $\Gamma_{ij}$ is
given by $\kappa_{ij}=-\nabla_{\Gamma_{ij}}\cdot\boldsymbol{n}_{\Gamma_{ij}}$.
Using the matching conditions (5.4), (6.5) and (6.6), we finally infer from
(6.40) that

$\displaystyle\begin{aligned}
\big{(}\vartheta_{1}^{j}-\vartheta_{1}^{i}\big{)}&=\sigma_{ij}\kappa_{ij}-\frac{1}{\gamma}\sum_{j=1}^{l}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\lambda_{0,n_{r}}\left[\overline{\rho\,}\left|\boldsymbol{w}_{0,n_{r}}\right|^{2}\right]_{i}^{j}\\\
&\quad+\frac{1}{\gamma}\sum_{r=1}^{l}\Bigg{\\{}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\\\
&\qquad\cdot\Big{(}\big{[}\overline{\mathbb{C}}\mathcal{E}(\boldsymbol{w}_{0,n_{r}}):\mathcal{E}(\boldsymbol{w}_{0,n_{r}})\big{]}_{i}^{j}-2\big{[}\overline{\mathbb{C}}\mathcal{E}(\boldsymbol{w}_{0,n_{r}})\boldsymbol{\nu}\cdot\nabla\boldsymbol{w}_{0,n_{r}}\boldsymbol{\nu}\big{]}_{i}^{j}\Big{)}\Bigg{\\}}\end{aligned}$
(6.45)

on $\Gamma_{ij}$ for all $i,j\in\\{1,\dots,N-1\\}$. In the case $j=N$ and
$i\neq N$, equation (6.40) simplifies to

$\displaystyle\begin{aligned}
\big{(}\vartheta_{1}^{j}-\vartheta_{1}^{i}\big{)}&=&&\sigma_{ij}\kappa_{ij}+\frac{1}{\gamma}\sum_{r=1}^{l}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\lambda_{0,n_{r}}\;\rho^{i}\left|(\boldsymbol{w}_{0,n_{r}})_{i}\right|^{2}\\\
&\phantom{=}&&-\frac{1}{\gamma}\sum_{r=1}^{l}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\mathbb{C}^{i}\mathcal{E}_{i}(\boldsymbol{w}_{0,n_{r}}):\mathcal{E}_{i}(\boldsymbol{w}_{0,n_{r}})\end{aligned}$
(6.46)

on $\Gamma_{iN}$ by the matching in (6.9) and (6.10).

To conclude this section we note that according to [20, Section 5.3] or [28,
Section 2.4] equation (6.25) induces a further solvability condition, namely
an angle condition for triple junctions. To see this, let us assume that the
regions $\Omega_{i}$, $\Omega_{j}$, $\Omega_{k}$ (with $i,j,k$ pairwise
different) meet at a triple point $m_{ijk}$ in the 2-dimensional case or on a
triple curve $m_{ijk}$ in the 3-dimensional case. Then the angle condition is
expressed via the normals of the three meeting interfaces as follows

$\displaystyle\sigma_{ij}\boldsymbol{n}_{\Gamma_{ij}}+\sigma_{jk}\boldsymbol{n}_{\Gamma_{jk}}+\sigma_{ki}\boldsymbol{n}_{\Gamma_{ki}}=0\quad\text{in
}m_{ijk}.$ (6.47)

This relation immediately implies the angle condition

$\displaystyle\frac{\sin(\theta_{ij})}{\sigma_{ij}}=\frac{\sin(\theta_{jk})}{\sigma_{jk}}=\frac{\sin(\theta_{ki})}{\sigma_{ki}},$

where $\theta_{ij}$ denotes the angle between $\boldsymbol{n}_{\Gamma_{jk}}$
and $\boldsymbol{n}_{\Gamma_{ki}}$. For the choice
$\psi_{0}({\boldsymbol{\varphi}})=\frac{1}{2}(1-{\boldsymbol{\varphi}}\cdot{\boldsymbol{\varphi}})$
we deduce that the transition energies denoted by
$\sigma_{ij},\sigma_{jk},\sigma_{ki}$ are always equal. This follows by
exploiting the symmetry of this potential in (6.28). Thus, in this case, we
know that triple junctions always occur at a $120^{\circ}$ contact angle.

##  7 The sharp-interface problem

Now we are in a position to state the complete problem that is obtained from
($SE^{\varepsilon}$) and ($GI^{\varepsilon}$) in the sharp-interface
situation.

###  7.1 The sharp-interface limit of the state equation

Therefore, we recall that the domain $\Omega$ is partitioned into $N$ regions
$\Omega_{i}$ for $i=1,\dots,N$ representing the presence of the $i$-th
material ($i<N$) or void ($i=N$) in its pure form. Those regions are separated
by interfaces $\Gamma_{ij}$. Furthermore we have chosen
$\boldsymbol{\eta}_{\Gamma_{ij}}$ to be the unit normal vector field on
$\Gamma_{ij}$ pointing from $\Omega_{i}$ into the region $\Omega_{j}$. This
means that

$\displaystyle
x+\delta\boldsymbol{\eta}_{\Gamma_{ij}}(x)\in\Omega_{j}\quad\text{and}\quad
x-\delta\boldsymbol{\eta}_{\Gamma_{ij}}(x)\in\Omega_{i}\quad\text{$x\in\Gamma_{ij}$
and $\delta>0$}.$

To capture the behavior of a function $\boldsymbol{v}$ across the interface
$\Gamma_{ij}$, we defined its jump by

$\displaystyle[\boldsymbol{v}]_{i}^{j}(x)\coloneqq\underset{\delta\searrow
0}{\lim}\;\Big{(}\boldsymbol{v}\big{(}x+\delta\boldsymbol{\eta}_{\Gamma_{ij}}(x)\big{)}-\boldsymbol{v}\big{(}x-\delta\boldsymbol{\eta}_{\Gamma_{ij}}(x)\big{)}\Big{)},$

for all $x\in\Gamma_{ij}$, see (5.5).

Combining the equations (${SE}^{i}_{0}$) derived in Claim 4.1 and the jump
conditions obtained in (6.5) and (6.10), we obtain the system

$\displaystyle\begin{cases}\begin{array}[]{rll}-\nabla\cdot\big{(}\mathbb{C}^{i}\mathcal{E}(\boldsymbol{w}_{0,n_{r}})\big{)}&=\lambda_{0,n_{r}}\rho^{i}\boldsymbol{w}_{0,n_{r}}&\quad\text{in
}\Omega_{i},\\\\[4.30554pt]
\left[\mathbb{C}\mathcal{E}(\boldsymbol{w}_{0,n_{r}})\boldsymbol{n}_{\Gamma_{ij}}\right]_{i}^{j}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{ij},\\\\[4.30554pt]
[\boldsymbol{w}_{0,n_{r}}]_{i}^{j}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{ij},\\\\[4.30554pt]
\mathbb{C}^{i}\mathcal{E}_{i}(\boldsymbol{w}_{0,n_{r}})\boldsymbol{n}_{\Gamma_{iN}}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{iN},\\\\[4.30554pt]
\mathbb{C}^{i}\mathcal{E}(\boldsymbol{w}_{0,n_{r}})\boldsymbol{n}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{0}\cap\partial\Omega_{i},\\\\[4.30554pt]
\boldsymbol{w}_{0,n_{r}}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{D}\cap\partial\Omega_{i},\end{array}\end{cases}$ (${SE}^{ij}_{r}$)

for $i,j=1,\dots,N-1$ and $r=1,\dots,l$, as the sharp-interface limit of the
state equation ($SE^{\varepsilon}$). Here, ${\boldsymbol{w}}_{0,n_{r}}$ is
normalized in the material regions, i.e.,

$\displaystyle
1=\sum_{i=1}^{N-1}\int_{\Omega_{i}}\rho^{i}\left|{\boldsymbol{w}}_{0,n_{r}}\right|^{2}\mathrm{\,d}x.$
(7.1)

Furthermore, we infer from (6.10) that

$\displaystyle[{\boldsymbol{w}}_{0,n_{r}}]_{i}^{N}=\boldsymbol{0}\quad\text{on
}\Gamma_{iN}$ (7.2)

for all $i\in\\{1,\dots,N-1\\}$ and each $r\in\\{1,\dots,l\\}$. However, this
condition does not provide any additional information as we do not know how
${\boldsymbol{w}}_{0,n_{r}}$ behaves in the void region. In particular, we see
that by interpreting (${SE}^{ij}_{r}$) as one system of PDEs in the material
region $\bigcup_{i=1}^{N-1}\Omega_{i}$, the homogeneous Neumann boundary
condition in the fourth line of (${SE}^{ij}_{r}$) is enough to obtain a closed
system.

Combining the Neumann type jump condition on $\Gamma_{ij}$ stated in the
second line of (${SE}^{ij}_{r}$) with the normality condition (7.1), we are
able to obtain the relation

$\displaystyle\int_{\Omega^{M}}\mathbb{C}^{M}\,\mathcal{E}({\boldsymbol{w}}_{0,n_{r}}):\mathcal{E}({\boldsymbol{w}}_{0,n_{r}})\mathrm{\,d}x=\lambda_{0,n_{r}},$
(7.3)

with

$\displaystyle\Omega^{M}:=\bigcup_{i=1}^{N-1}\Omega_{i}\quad\text{and}\quad\mathbb{C}^{M}:=\left(\sum_{i=1}^{N-1}\mathbb{C}^{i}\;\mathds{1}_{\Omega_{i}}\right),$

where $\mathds{1}_{\Omega_{i}}$ denotes the characteristic function on
$\Omega_{i}$. This means that the eigenvalue $\lambda_{0,n_{r}}$ in the sharp-
interface setting is indeed solely determined by an eigenvalue equation on the
material region $\Omega^{M}$ but does not have any contribution from the void
region.

To verify (7.3), we test (${SE}^{ij}_{r}$) with ${\boldsymbol{w}}_{0,n_{r}}$
and integrate by parts. This yields

$\displaystyle\begin{aligned}
&\int_{\Omega_{i}}\mathbb{C}^{i}\mathcal{E}({\boldsymbol{w}}_{0,n_{r}}):\mathcal{E}({\boldsymbol{w}}_{0,n_{r}})\mathrm{\,d}x-\int_{\partial\Omega_{i}}\mathbb{C}^{i}\mathcal{E}({\boldsymbol{w}}_{0,n_{r}})\boldsymbol{n}_{\Gamma_{i}}\cdot{\boldsymbol{w}}_{0,n_{r}}\mathrm{\,d}\Gamma=\lambda_{0,n_{r}}\\\
&\quad=\int_{\Omega_{i}}\rho^{i}\left|{\boldsymbol{w}}_{0,n_{r}}\right|^{2}\mathrm{\,d}x,\end{aligned}$
(7.4)

for all $i\in\\{1,\dots,N-1\\}$, where $\boldsymbol{n}_{\Gamma_{i}}$ stands
for the outer unit normal vector field of $\partial\Omega_{i}$. Noticing that
the outer unit normal vector simply switches its sign on neighboring
boundaries, we now use the second and the fourth line of (${SE}^{ij}_{r}$) to
infer

$\displaystyle\sum_{i=1}^{N-1}\int_{\partial\Omega_{i}}\mathbb{C}^{i}\mathcal{E}({\boldsymbol{w}}_{0,n_{r}})\boldsymbol{n}_{\Gamma_{i}}\cdot{\boldsymbol{w}}_{0,n_{r}}\mathrm{\,d}\Gamma=0.$

Thus, summing the equations (7.4) from $i=1$ to $N-1$ and using property
(7.1), we conclude

$\displaystyle\sum_{i=1}^{N-1}\int_{\Omega_{i}}\mathbb{C}^{i}\mathcal{E}({\boldsymbol{w}}_{0,n_{r}}):\mathcal{E}({\boldsymbol{w}}_{0,n_{r}})\mathrm{\,d}x=\lambda_{0,n_{r}}.$

By the linearity of the integral, this directly proves (7.3).

######  Remark 7.1.

As a refinement of Remark 4.2 (a), we now see that as long as at least one of
the material regions $\Omega_{1},\dots,\Omega_{N-1}$ shares a sufficiently
nice part of its boundary with $\Gamma_{D}$, we can apply Korn’s inequality in
order to deduce that all $\lambda_{0,n_{r}}$ are strictly positive. From a
physical point of view, this is reasonable since if the material region
$\Omega^{M}$ of the structure is not attached to some fixed boundary the shape
can freely move within the design domain just by translation without
exhibiting any vibrations.

###  7.2 The sharp-interface limit of the first-order optimality condition

Now let us turn to the limit of the gradient inequality ($GI^{\varepsilon}$).
For the sake of completeness, let us restate our final results from the
previous section, i.e., (6.45) and (6.46). We have

$\displaystyle
0=\gamma\sigma_{ij}\kappa_{ij}-\sum_{j=1}^{l}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\lambda_{0,n_{r}}\left[\overline{\rho\,}\left|\boldsymbol{w}_{0,n_{r}}\right|^{2}\right]_{i}^{j}+\gamma\big{(}\vartheta_{1}^{i}-\vartheta_{1}^{j}\big{)}$
$\displaystyle\qquad\quad\begin{aligned}
+\sum_{r=1}^{l}&\Big{\\{}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\\\
&\quad\cdot\Big{(}\big{[}\overline{\mathbb{C}}\mathcal{E}(\boldsymbol{w}_{0,n_{r}}):\mathcal{E}(\boldsymbol{w}_{0,n_{r}})\big{]}_{i}^{j}-2\big{[}\overline{\mathbb{C}}\mathcal{E}(\boldsymbol{w}_{0,n_{r}})\boldsymbol{\nu}\cdot\nabla\boldsymbol{w}_{0,n_{r}}\boldsymbol{\nu}\big{]}_{i}^{j}\Big{)}\Big{\\}}\end{aligned}$
(7.5)

on $\Gamma_{ij}$ for all $i,j=1\dots,N-1$, and

$\displaystyle\begin{aligned}
0&=\gamma\sigma_{iN}\kappa_{iN}+\sum_{r=1}^{l}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\lambda_{0,n_{r}}\;\rho^{i}\left|(\boldsymbol{w}_{0,n_{r}})_{i}\right|^{2}+\gamma\big{(}\vartheta_{1}^{i}-\vartheta_{1}^{N}\big{)}\\\
&\quad-\sum_{r=1}^{l}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\mathbb{C}^{i}\mathcal{E}_{i}(\boldsymbol{w}_{0,n_{r}}):\mathcal{E}_{i}(\boldsymbol{w}_{0,n_{r}})\end{aligned}$
(7.6)

on $\Gamma_{iN}$ for all $i=1\dots,N-1$ if $j=N$.

Here $\sigma_{ij}$ is defined as in (6.28) and stands for the total energy of
a transition across the interface $\Gamma_{ij}$. The vector
$\boldsymbol{\vartheta}_{1}\in\mathbb{R}^{N}$ denotes the
$\mathcal{O}(\varepsilon)$-contribution of the Lagrange-multiplier resulting
from the integral constraint
$\fint_{\Omega}{\boldsymbol{\varphi}}^{\varepsilon}\mathrm{\,d}x=\boldsymbol{m}$
that is hidden in the condition
${\boldsymbol{\varphi}}^{\varepsilon}\in\boldsymbol{\mathcal{G}}^{\boldsymbol{m}}$
(cf. Theorem 3.12).

Recalling (6.47), we additionally have the triple junction condition at any
junction $m_{ijk}$ with pairwise disjoint $i,j,k\\{1,\dots,N\\}$

$\displaystyle\sigma_{ij}\boldsymbol{n}_{\Gamma_{ij}}+\sigma_{jk}\boldsymbol{n}_{\Gamma_{jk}}+\sigma_{ki}\boldsymbol{n}_{\Gamma_{ki}}=0\quad\text{in
}m_{ijk}.$

###  7.3 The sharp-interface optimality system in the case of only one
material

We now want to state above equations for the simplest case of only one single
material (i.e., $N=2$) as this is the scenario we further study in the
subsequent sections.

In this case, we have $\Omega=\Omega^{M}\cup\Omega^{V}$, where $\Omega^{M}$
and $\Omega^{V}$ denote the material and the void parts of the domain,
respectively. We now denote the interface separating the two phases by
$\Gamma_{MV}$, its outer unit normal vector field by
$\boldsymbol{n}_{\Gamma_{MV}}$ and its mean curvature by
$\kappa_{MV}=-\nabla_{\Gamma_{MV}}\cdot\boldsymbol{n}_{\Gamma_{MV}}$. Using
the notation $\Gamma_{D}^{M}\coloneqq\Gamma_{D}\cap\partial\Omega^{M}$ and
$\Gamma_{0}^{M}\coloneqq\Gamma_{0}\cap\partial\Omega^{M}$, we obtain from
(${SE}^{ij}_{r}$), (7.1) and (7.2) the state equation

$\displaystyle\begin{cases}\begin{array}[]{rll}-\nabla\cdot\left(\mathbb{C}^{M}\mathcal{E}(\boldsymbol{w}_{0,n_{r}})\right)&=\lambda_{0,n_{r}}\rho^{M}\boldsymbol{w}_{0,n_{r}}&\quad\text{in
}\Omega^{M},\\\\[4.30554pt]
\mathbb{C}^{M}\mathcal{E}_{M}(\boldsymbol{w}_{0,n_{r}})\;\boldsymbol{n}_{\Gamma_{MV}}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{MV},\\\\[4.30554pt]
\boldsymbol{w}_{0,n_{r}}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{D}^{M},\\\\[4.30554pt]
\mathbb{C}^{M}\mathcal{E}(\boldsymbol{w}_{0,n_{r}})\;\boldsymbol{n}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{0}^{M},\end{array}\end{cases}$ (${SE}^{MV}_{r}$)

for $r=1,\dots,l$, along with the first-order necessary optimality condition

$\displaystyle\begin{aligned}
0&=\gamma\,\sigma_{MV}\,\kappa_{MV}+\sum_{r=1}^{l}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,{n_{1}}},\dots,\lambda_{0,{n_{l}}}\right)\lambda_{0}^{n_{r}}\rho^{M}\left|(\boldsymbol{w}_{0,n_{r}})_{M}\right|^{2}\\\
&\quad-\sum_{r=1}^{l}[\partial_{\lambda_{n_{r}}}\Psi]\left(\lambda_{0,n_{1}},\dots,\lambda_{0,n_{l}}\right)\mathbb{C}^{M}\mathcal{E}_{M}(\boldsymbol{w}_{0,n_{r}}):\mathcal{E}_{M}(\boldsymbol{w}_{0,n_{r}})+\gamma\left(\vartheta_{1}^{1}-\vartheta_{1}^{2}\right)\end{aligned}$
(7.7)

on $\Gamma_{MV}$. This means that the functions $\boldsymbol{w}_{0,n_{r}}$ are
eigenfunctions to the eigenvalues $\lambda_{0,n_{r}}$ which essentially solve
the eigenvalue problem for the elasticity equation subject to a homogeneous
Neumann boundary condition on the shape $\Omega^{M}$.

######  Remark 7.2.

Note that, in general, one cannot predict the behavior of solutions to
(${SE}^{MV}_{r}$). If $\Omega^{M}$ is merely a set of finite perimeter that
does not have a Lipschitz boundary or if
$\Gamma_{MV}\cap\Gamma_{D}^{M}=\emptyset$, the classical spectral theory (as
applied in Section 2.5) does not provide us with an infinite sequence of
positive eigenvalues. Nevertheless, as we want to consider a well posed
minimization problem and want to calculate shape derivatives associated to
this problem, we assume that these issues do not occur. In particular, we
always assume $\Omega^{M}$ to be sufficiently smooth and $\partial\Omega^{M}$
to have a suitably nice intersection with $\Gamma^{M}_{D}$ such that an
infinite sequence of positive eigenvalues actually exists (see also Remark
7.1).

##  8 Relating the first-order optimality condition
to classical shape calculus

We now want to compare the above results, especially (7.7), to the results in
[7], which were obtained using shape calculus. Our goal is to justify that the
gradient equality (7.7) is indeed the first-order condition of a sharp-
interface eigenvalue optimization problem, which is formally the limit of the
diffuse-interface problem we started with. Therefore, we need to fit the
notation of [7] to our setting.

As above consider the situation $N=2$, i.e.,
$\Omega=\Omega^{M}\cup\Omega^{V}$. Denote with $P_{\Omega}(\Omega^{M})$ the
perimeter of the shape $\Omega^{M}$ within the design domain $\Omega$, which
is given by the Hausdorff measure
$\mathcal{H}^{d-1}(\partial\Omega^{M}\cap\Omega)$ provided that $\Omega^{M}$
is non-empty and sufficiently smooth. Furthermore, we consider a prescribed
mass $m=\big{|}\Omega^{M}\big{|}<|\Omega|$. In order to be consistent with the
notation used in the previous chapters, we choose
$\boldsymbol{m}=(m_{1},m_{2})^{T}\in\Sigma^{2}$ with $m_{1}=m|\Omega|^{-1}$
and $m_{2}=1-m_{1}$. Then the sharp-interface structural optimization problem
that we intend to approximate via our diffuse-interface problem (2.17) reads
as

$\displaystyle\begin{cases}\begin{aligned}
&\min&&J(\Omega^{M})\coloneqq\Psi(\lambda_{n_{1}},\dots,\lambda_{n_{l}})+\gamma\,\sigma_{MV}\,P_{\Omega}(\Omega^{M}),\\\\[4.30554pt]
&\text{over}&&\mathcal{U}^{\text{ad}}=\left\\{\Omega^{M}\subset\Omega:\big{|}\Omega^{M}\big{|}=m\right\\},\\\\[4.30554pt]
&\text{s.t.}&&\eqref{SErMV}\begin{cases}\begin{array}[]{rll}-\nabla\cdot\left(\mathbb{C}^{M}\mathcal{E}(\boldsymbol{w}_{n_{r}})\right)&=\lambda_{n_{r}}\rho^{M}\boldsymbol{w}_{n_{r}}&\quad\text{in
}\Omega^{M},\\\\[4.30554pt]
\mathbb{C}^{M}\mathcal{E}_{M}(\boldsymbol{w}_{n_{r}})\;\boldsymbol{n}_{\Gamma_{MV}}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{MV},\\\\[4.30554pt]
\mathbb{C}^{M}\mathcal{E}(\boldsymbol{w}_{n_{r}})\;\boldsymbol{n}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{0}^{M},\\\\[4.30554pt]
\boldsymbol{w}_{n_{r}}&=\boldsymbol{0}&\quad\text{on
}\Gamma_{D}^{M},\end{array}\end{cases}\\\\[4.30554pt] &&&\text{for
all}\;r\in\\{1,\dots,l\\}.\end{aligned}\end{cases}$ ($\mathcal{P}_{l}^{0}$)

This system is the sharp-interface limit problem associated to the diffuse-
interface problem (2.17), where the side condition is exactly the sharp-
interface state equation (${SE}^{MV}_{r}$) and the perimeter
$\sigma_{MV}P_{\Omega}(\Omega^{M})$ is the rigorous $\Gamma$-limit of the
Ginzburg–Landau energy, see [14]. We recall that the constant $\sigma_{MV}$ we
obtained in (6.28) is exactly the one obtained in [14] in terms of the
rigorous $\Gamma$-limit, which is denoted by
$d(\boldsymbol{e}_{i},\boldsymbol{e}_{j})$ there. In particular, $\sigma_{MV}$
is independent of the shape $\Omega^{M}$.

In case an ambiguity might arise, we indicate the shape dependency explicitly
in the eigenfunctions and eigenvalues, i.e., we write
$(\lambda_{n_{r}}(\Omega^{M}),{\boldsymbol{w}}_{n_{r}}(\Omega^{M}))$ for
$r=1,\dots,l$. Now, we want to apply the calculus of shape derivatives from
[7, Thm. 2.5] to our situation. We obtain the following statement.

######  Theorem 8.1.

Let $\Omega^{M}$ be a smooth bounded open set and let $\boldsymbol{\theta}\in
W^{1,\infty}(\mathbb{R}^{d},\mathbb{R}^{d})$ with
${\boldsymbol{\theta}\cdot\boldsymbol{n}_{\Gamma_{\partial\Omega^{M}}}=0}$ on
$\partial\Omega^{M}\backslash\Gamma_{MV}$. We further assume that for
$r=1,\dots,l$, the eigenfunctions $\boldsymbol{w}_{n_{r}}(\Omega^{M})$ in
$({SE}^{MV}_{r})$ are sufficiently smooth, say
$\boldsymbol{w}_{n_{r}}(\Omega^{M})\in H^{2}(\Omega^{M};\mathbb{R}^{d})$.

Then, if the involved eigenvalues $\lambda_{n_{r}}$ for $r=1,\dots,l$ are all
simple, the shape derivative of $J$ at the shape $\Omega^{M}$ in the direction
$\boldsymbol{\theta}$ fulfills the equation

$\displaystyle\begin{aligned}
J^{\prime}(\Omega^{M})(\boldsymbol{\theta})&=\sum_{r=1}^{l}\Bigg{\\{}[\partial_{\lambda_{n_{r}}}\Psi](\lambda_{n_{1}}(\Omega^{M}),\dots,\lambda_{n_{l}}(\Omega^{M}))\\\\[-4.30554pt]
&\qquad\qquad\cdot\Bigg{[}\int_{\Gamma_{MV}}\mathbb{C}^{M}\mathcal{E}(\boldsymbol{w}_{n_{r}}(\Omega^{M})):\mathcal{E}(\boldsymbol{w}_{n_{r}}(\Omega^{M}))\boldsymbol{\theta}\cdot\boldsymbol{n}_{\Gamma_{MV}}\textup{\,d}\mathcal{H}^{d-1}\;\\\\[-4.30554pt]
&\qquad\qquad\qquad-\lambda_{n_{r}}(\Omega^{M})\int_{\Gamma_{MV}}\rho^{M}\big{|}\boldsymbol{w}_{n_{r}}(\Omega^{M})\big{|}^{2}\boldsymbol{\theta}\cdot\boldsymbol{n}_{\Gamma_{MV}}\textup{\,d}\mathcal{H}^{d-1}\Bigg{]}\Bigg{\\}}\\\
&\qquad-\int_{\Gamma_{MV}}\gamma\sigma_{MV}\,\kappa_{MV}\,\boldsymbol{\theta}\cdot\boldsymbol{n}_{\Gamma_{MV}}\textup{\,d}\mathcal{H}^{d-1}.\end{aligned}$
(8.1)

Here, the shape derivative of $J$ at a shape $\Omega^{M}$ is defined as the
Fréchet derivative of the functional

$\displaystyle
W^{1,\infty}(\mathbb{R}^{d};\mathbb{R}^{d})\to\mathbb{R},\quad\boldsymbol{\zeta}\mapsto
J\big{(}(\mathrm{Id}+\boldsymbol{\zeta})\Omega^{M}\big{)}$

evaluated at $\boldsymbol{\zeta}=\boldsymbol{0}$.

######  Remark 8.2.

  1. (a)

Note that the simplicity of eigenvalues is crucial here. Only then it is
guaranteed that the eigenvalues and eigenfunctions depend on the domain
$\Omega^{M}$ in a differentiable way. For a comprehensive overview over the
differentiablitiy of spectral quantities with repsect to the domain we refer
to [50, Section 5.7].

  2. (b)

For $\boldsymbol{\zeta}\in W^{1,\infty}(\mathbb{R}^{d};\mathbb{R}^{d})$ the
application

$\displaystyle T_{\boldsymbol{\zeta}}:\mathbb{R}^{d}\to\mathbb{R}^{d},\quad
x\mapsto(\mathrm{Id}+\boldsymbol{\zeta})(x),$

is invertible if $\left\|\boldsymbol{\zeta}\right\|_{W^{1,\infty}}<1$, and it
holds $(\mathrm{Id}+\boldsymbol{\zeta})^{-1}-\mathrm{Id}\in
W^{1,\infty}(\mathbb{R}^{d};\mathbb{R}^{d})$ with

$\displaystyle\left\|(\mathrm{Id}+\boldsymbol{\zeta})^{-1}-\mathrm{Id}\right\|_{W^{1,\infty}}\leq\left\|\boldsymbol{\zeta}\right\|_{W^{1,\infty}}(1-\left\|\boldsymbol{\zeta}\right\|_{W^{1,\infty}})^{-1}.$

This means the family $(T_{\boldsymbol{\zeta}})_{\boldsymbol{\zeta}\in
W^{1,\infty}}$ describes diffeomorphic perturbations of $\Omega^{M}$ “close”
to $\Omega^{M}$ if $\left\|\boldsymbol{\zeta}\right\|_{W^{1,\infty}}$ is
small, motivating the definition of the shape derivative above. For a detailed
discussion of this concept, we refer to [50, Section 5.2].

###### Proof.

We proceed analogously to [7, Theorem 2.5]. In the following,
$\Omega_{\boldsymbol{\zeta}}=(\mathrm{Id}+\boldsymbol{\zeta})(\Omega^{M})$
denotes the perturbation of $\Omega^{M}$ associated with a sufficiently small
$\boldsymbol{\zeta}\in W^{1,\infty}(\mathbb{R}^{d};\mathbb{R}^{d})$. First of
all, for $\boldsymbol{v}_{n_{r}}\in H^{1}(\mathbb{R}^{d};\mathbb{R}^{d})$ with
$r=1,\dots,l$, we introduce the Lagrangian

$\displaystyle\mathcal{L}(\Omega_{\boldsymbol{\zeta}},\boldsymbol{v}_{n_{1}},\dots,\boldsymbol{v}_{n_{l}})$
$\displaystyle\quad=\Psi\left(\frac{\int_{\Omega_{\boldsymbol{\zeta}}}\mathbb{C}^{M}\mathcal{E}({\boldsymbol{v}_{n_{1}}}):\mathcal{E}({\boldsymbol{v}_{n_{1}}})\mathrm{\,d}x}{\int_{\Omega_{\boldsymbol{\zeta}}}\rho^{M}\left|\boldsymbol{v}_{n_{1}}\right|^{2}\mathrm{\,d}x},\dots,\frac{\int_{\Omega_{\boldsymbol{\zeta}}}\mathbb{C}^{M}\mathcal{E}({\boldsymbol{v}_{n_{l}}}):\mathcal{E}({\boldsymbol{v}_{n_{l}}})\mathrm{\,d}x}{\int_{\Omega_{\boldsymbol{\zeta}}}\rho^{M}\left|\boldsymbol{v}_{n_{l}}\right|^{2}\mathrm{\,d}x}\right)$
$\displaystyle\qquad\quad+\gamma\sigma_{MV}P(\Omega_{\boldsymbol{\zeta}})\mathrm{\,d}s.$

For the partial Fréchet derivatives of the Lagrangian with respect to
$\boldsymbol{v}_{n_{r}}$ for $r=1,\dots,l$ at the point
$(\Omega_{\boldsymbol{\zeta}},\boldsymbol{w}_{n_{1}}(\Omega_{\boldsymbol{\zeta}}),\dots,\boldsymbol{w}_{n_{l}}(\Omega_{\boldsymbol{\zeta}}))$,
we obtain

$\displaystyle\partial_{\boldsymbol{v}_{n_{r}}}\mathcal{L}\big{(}\Omega_{\boldsymbol{\zeta}},\boldsymbol{w}_{n_{1}}(\Omega_{\boldsymbol{\zeta}}),\dots,\boldsymbol{w}_{n_{l}}(\Omega_{\boldsymbol{\zeta}})\big{)}=0.$
(8.2)

This is simply due to the fact, that the derivative of the Rayleigh quotient

$\displaystyle\mathcal{R}_{\boldsymbol{\zeta}}:H^{1}(\mathbb{R}^{d};\mathbb{R}^{d})\to\mathbb{R},\quad\boldsymbol{v}\mapsto\frac{\int_{\Omega_{\boldsymbol{\zeta}}}\mathbb{C}^{M}\mathcal{E}({\boldsymbol{v}}):\mathcal{E}({\boldsymbol{v}})\mathrm{\,d}x}{\int_{\Omega_{\boldsymbol{\zeta}}}\rho^{M}\left|\boldsymbol{v}\right|^{2}\mathrm{\,d}x},$

evaluated at an eigenfunction
$\boldsymbol{w}_{n}=\boldsymbol{w}_{n}(\Omega_{\boldsymbol{\zeta}})$ reads as

$\displaystyle\mathcal{R}^{\prime}_{\boldsymbol{\zeta}}({\boldsymbol{w}}_{n})\boldsymbol{v}$
$\displaystyle=\frac{2\int_{\Omega_{\boldsymbol{\zeta}}}\mathbb{C}^{M}\mathcal{E}(\boldsymbol{w}_{n}):\mathcal{E}(\boldsymbol{v})\mathrm{\,d}x\int_{\Omega_{\boldsymbol{\zeta}}}\rho^{M}\left|\boldsymbol{w}_{n}\right|^{2}\mathrm{\,d}x}{(\int_{\Omega_{\boldsymbol{\zeta}}}\rho^{M}\left|\boldsymbol{w}_{n}\right|^{2}\mathrm{\,d}x)^{2}}$
$\displaystyle\phantom{=}-\frac{2\int_{\Omega_{\boldsymbol{\zeta}}}\mathbb{C}^{M}\mathcal{E}(\boldsymbol{w}_{n}):\mathcal{E}(\boldsymbol{w}_{n})\mathrm{\,d}x\int_{\Omega_{\boldsymbol{\zeta}}}\rho^{M}\boldsymbol{w}_{n}\cdot\boldsymbol{v}\mathrm{\,d}x}{(\int_{\Omega_{\boldsymbol{\zeta}}}\rho^{M}\left|\boldsymbol{w}_{n}\right|^{2}\mathrm{\,d}x)^{2}}$

and this vanishes due to (${SE}^{MV}_{r}$).

On the other hand, recalling the definition of $J$ in ($\mathcal{P}_{l}^{0}$),
we obviously have

$\displaystyle
J(\Omega_{\boldsymbol{\zeta}})=\mathcal{L}(\Omega_{\boldsymbol{\zeta}},\boldsymbol{w}_{n_{1}}(\Omega_{\boldsymbol{\zeta}}),\dots,\boldsymbol{w}_{n_{l}}(\Omega_{\boldsymbol{\zeta}}))$

as the eigenvalues can be expressed by the corresponding Rayleigh quotients.
Note that due to the differentiability of eigenfunctions as discussed in
Remark 8.2, we can now apply the chain rule. Thus, using (8.2) we infer that
the shape derivative is given by

$\displaystyle J^{\prime}(\Omega^{M})$
$\displaystyle=\frac{\text{d}}{\text{d}\boldsymbol{\zeta}}[J((\mathrm{Id}+\boldsymbol{\zeta})(\Omega^{M}))]_{\boldsymbol{\zeta}=\boldsymbol{0}}$
$\displaystyle=\frac{\text{d}}{\text{d}\boldsymbol{\zeta}}[\mathcal{L}((\mathrm{Id}+\boldsymbol{\zeta})(\Omega^{M}),{\boldsymbol{w}}_{n_{1}}(\Omega^{M}),\dots,{\boldsymbol{w}}_{n_{l}}(\Omega^{M}))]_{\boldsymbol{\zeta}=\boldsymbol{0}}$

Applying the formulas for shape derivatives in [7, Lemma 2.3], we deduce

$\displaystyle\begin{aligned}
J^{\prime}(\Omega^{M})(\boldsymbol{\theta})&=\sum_{r=1}^{l}\Bigg{\\{}[\partial_{\lambda_{n_{r}}}\Psi](\lambda_{n_{1}}(\Omega^{M}),\dots,\lambda_{n_{l}}(\Omega^{M}))\\\\[-4.30554pt]
&\qquad\qquad\cdot\Bigg{[}\int_{\partial\Omega^{M}}\mathbb{C}^{M}\mathcal{E}(\boldsymbol{w}_{n_{r}}(\Omega^{M})):\mathcal{E}(\boldsymbol{w}_{n_{r}}(\Omega^{M}))\boldsymbol{\theta}\cdot\boldsymbol{n}_{\partial\Omega^{M}}\textup{\,d}\mathcal{H}^{d-1}\;\\\\[-4.30554pt]
&\qquad\qquad\qquad-\lambda_{n_{r}}(\Omega^{M})\int_{\partial\Omega^{M}}\rho^{M}\big{|}\boldsymbol{w}_{n_{r}}(\Omega^{M})\big{|}^{2}\boldsymbol{\theta}\cdot\boldsymbol{n}_{\partial\Omega^{M}}\textup{\,d}\mathcal{H}^{d-1}\Bigg{]}\Bigg{\\}}\\\
&\qquad-\int_{\partial\Omega^{M}}\gamma\sigma_{MV}\,\kappa_{M}\,\boldsymbol{\theta}\cdot\boldsymbol{n}_{\partial\Omega^{M}}\textup{\,d}\mathcal{H}^{d-1},\end{aligned}$

where $\kappa_{M}$ denotes the mean curvature of $\partial\Omega^{M}$. By the
assumption ${\boldsymbol{\theta}\cdot\boldsymbol{n}_{\partial\Omega^{M}}=0}$
on $\partial\Omega^{M}\backslash\Gamma_{MV}$, the boundary integrals vanish on
$\partial\Omega^{M}\backslash\Gamma_{MV}$ and we thus arrive at (8.1). Note
that in [7], the mean curvature is defined as
$\kappa=\nabla_{\partial\Omega^{M}}\cdot\boldsymbol{n}_{\partial\Omega^{M}}$,
whereas (in accordance with (4.23)) our mean curvature is given by
$\kappa=-\nabla_{\partial\Omega^{M}}\cdot\boldsymbol{n}_{\partial\Omega^{M}}$.
This explains the negative sign of our term involving $\kappa_{M}$. ∎

######  Remark 8.3.

The preceding theorem shows that using the approach of classical shape
calculus and additionally taking the volume constraint
$\big{|}\Omega^{M}\big{|}=m$ into account, we recover the gradient equality
(7.7) since the volume constraint produces a Lagrange multiplier as in our
previous analysis. This justifies our formal approach from the viewpoint of
classical shape calculus since (7.7) can be interpreted as the first-order
necessary optimality condition of the shape optimization problem
($\mathcal{P}_{l}^{0}$).

##  9 Numerical Examples

In the following, we present numerical results that illustrate the
applicability of our approach to find optimal topologies. After a brief
introduction of the numerical method, we investigate the dependence of
solutions on the parameter $\varepsilon$ in Section 9.1. Therefore, we study a
particular setting of an elastic beam that is known from literature (cf. [7]).
In Section 9.2, we consider a joint optimization of $\lambda_{1}$ and
$\lambda_{2}$ for this beam setup, and in Section 9.3, we investigate an
extended optimization problem to not only optimize the shape and topology of
this beam with respect to its first eigenvalue but also its compliance.

As in Subsection 7.3 and Section 8, we restrict ourselves to the case of only
two phases, i.e., material and void. In this situation, the vector-valued
phase-field ${\boldsymbol{\varphi}}=(\varphi^{1},\varphi^{2})$ can be
represented by a scalar order parameter

$\displaystyle\varphi:=\varphi^{1}-\varphi^{2}\in H^{1}(\Omega)\cap
L^{\infty}(\Omega),$

and in the Ginzburg–Landau energy
$E^{\varepsilon}(\varphi)=\int_{\Omega}\frac{\varepsilon}{2}|\nabla\varphi|^{2}+\psi(\varphi)\mathrm{\,d}x$,
we choose

$\displaystyle\psi(\varphi):=\psi_{0}(s)+I_{[-1,1]}(\varphi)\quad\text{with}\quad\psi_{0}(\varphi)=\tfrac{1}{2\varepsilon}(1-\varphi^{2}),$

where $I_{[-1,1]}$ is the indicator functional of the interval $[-1,1]$. This
means that $\varphi$ attains its values in $[-1,1]$, where “$1$” represents
the material and “$-1$” represents the void. The elastic tensor
$\mathbb{C}(\varphi)$ now is defined as

$\displaystyle\mathbb{C}(\varphi)\mathcal{E}(w):=\alpha(\varphi)\big{(}2\mu\,\mathcal{E}(w)+\ell\operatorname{tr}\big{(}\mathcal{E}(w)\big{)}\,\mathcal{I}\big{)}$
(9.1)

for Lamé parameters $\mu,\ell>0$ and the quadratic interpolation function
$\alpha(\varphi)$ satisfying $\alpha(1)=1$,
$\alpha(-1)=\underline{\alpha}\varepsilon$, and $\alpha^{\prime}(-1)=0$ for
some constant $\underline{\alpha}$. The eigenvalue equation is given by

$\displaystyle-\nabla\cdot\left[\mathbb{C}(\varphi)\mathcal{E}(w)\right]=\lambda\,\beta(\varphi)\rho\,w,$
(9.2)

with the quadratic interpolation function $\beta(\varphi)$ satisfying
$\beta(1)=1$, $\beta(-1)=\underline{\beta}\varepsilon^{2}$ and
${\beta^{\prime}(-1)=0}$ as well as an additional density function $\rho$ that
might depend on the spatial variable. If not stated differently, we use
$\underline{\alpha}=2\cdot 10^{-4}$ and $\underline{\beta}=10^{-4}$. Note that
this choice of interpolation functions is exactly reflected by the choice
(2.10) and the scaling choice of $k=1,l=2$ as discussed in Section 4.2. More
precisely we have

$\displaystyle\begin{aligned}
\beta(\varphi)&=\beta_{M}\left(1-\frac{1-\varphi}{2}\right)+\underline{\beta}\varepsilon^{2}\beta_{V}\left(\frac{1-\varphi}{2}\right)\\\
&=\left(1-\frac{1-\varphi}{2}\right)^{2}+\underline{\beta}\varepsilon^{2}\left[-\left(\frac{1-\varphi}{2}-1\right)^{2}+1\right],\end{aligned}$
(9.3)

and analogously for $\alpha$.

#### Numerical Solution Method.

The numerical implementation is based on linear finite elements for all
functions provided by the finite element package FEniCs [9, 57] together with
the PETSc linear algebra backend [12, 13]. For the eigenvalue problem, we use
the package SLEPc [51]. The optimization problem is solved by the VMPT method
that is proposed in [24]. In our case, it can be understood as an extension of
the projected gradient method into the space $H^{1}(\Omega)\cap
L^{\infty}(\Omega)$. We refer to [24, 43, 44] for more details.

###  9.1 Numerical investigation of the sharp-interface limit $\varepsilon\to
0$

In this section, to illustrate the sharp-interface limit, we present numerical
results for a sequence of decreasing values of $\varepsilon$.

We use the setup from [7, Sec. 7.1] to find a cantilever beam with maximal
first eigenvalue, i.e., we choose $\Psi(\lambda_{1})=-\lambda_{1}$. Our
computational domain is given by $\Omega=(0,2)\times(0,1)$. The Young’s
modulus is $E=1$ and Poisson’s ratio is $\nu=0.3$ leading to $\mu\approx 0.38$
and $\ell\approx 0.58$. We define the subset
$\Omega_{\rho}=(1.9,2.0)\times(0.45,0.55)$ and set $\rho(x)=1$ if
$x\not\in\Omega_{\rho}$ and $\rho(x)=100$ if $x\in\Omega_{\rho}$. We also fix
$\varphi(x)=1$ for all $x\in\Omega_{\rho}$. The beam is supposed to be
attached to the wall at the left boundary of $\Omega$, i.e., at
$\Gamma_{D}=\\{(0,\eta)\mid\eta\in(0,1)\\}\subset\partial\Omega$. This leads
to the boundary condition $w=0$ on $\Gamma_{D}$. We further set
$\Gamma_{0}=\partial\Omega\setminus\Gamma_{D}$ and we fix $\gamma=10^{-4}$ and
$\int_{\Omega}\varphi=0$.

Similar as in [7, Sec. 7.1], we start our optimization process with a
checkerboard type initial function given by
$\varphi_{0}(x)=\operatorname{sign}\big{(}v(x)\big{)}\,\big{|}v\big{(}x\big{)}\big{|}^{0.3}$
with $v(x)=\cos(3\pi x_{1})\cos(4\pi x_{2})$ for all $x\in\Omega$. We want to
emphasize that this problem is expected to have many local minima and thus,
the choice of initial function can significantly influence the shape and
topology of the local minimizer found by our numerical method.

We now solve the optimization problem for a decreasing sequence of values of
$\varepsilon$. In Table 1, we present the values of $\varepsilon$ together
with the corresponding value of the Ginzburg–Landau energy
$E^{\varepsilon}(\varphi)=\int_{\Omega}\frac{\varepsilon}{2}|\nabla\varphi|^{2}+\frac{1}{2\varepsilon}(1-\varphi^{2})\mathrm{\,d}x$
and the eigenvalue $\lambda_{1}$. Recall here that the values of the
Ginzburg–Landau energy converge to a weighted perimeter of the shape in the
sharp interface limit $\varepsilon\to 0$. In Figure 1, we present the zero
level lines of the (locally optimal) shapes we obtain for different values of
$\varepsilon$. Here we started with $\varepsilon=0.08$ and used the local
optimum as initial value for subsequent simulations.

$\varepsilon$ | $80\cdot 10^{-3}$ | $40\cdot 10^{-3}$ | $20\cdot 10^{-3}$ | $10\cdot 10^{-3}$ | $5\cdot 10^{-3}$ | $2.5\cdot 10^{-3}$ | $1.25\cdot 10^{-3}$
---|---|---|---|---|---|---|---
$\gamma E^{\varepsilon}(\varphi)$ | 0.00117 | 0.00120 | 0.00117 | 0.00115 | 0.00114 | 0.00114 | 0.00114
$\lambda_{1}$ | 0.01574 | 0.01626 | 0.01658 | 0.01678 | 0.01692 | 0.01699 | 0.01703
Table 1: Scaled Ginzburg–Landau energy $\gamma E^{\varepsilon}(\varphi)$ and
principal eigenvalue $\lambda_{1}$ of the optimal beam shape for decreasing
values of $\varepsilon$. This indicates that the values
$E^{\varepsilon}(\varphi)$ and $\lambda_{1}$ converge as $\varepsilon$
decreases. Figure 1: The zero level lines of the beam for the $\varepsilon\to
0$ test for all tested $\varepsilon$. The darker the line is, the smaller is
$\varepsilon$. We observe that the interface seems to stabilize with
decreasing values of $\varepsilon$ and that it only mildly depends on
$\varepsilon$.

###  9.2 Optimization of a beam

As a first test, we illustrate the influence of the regularization strength
$\gamma$ on the found structure. The parameter $\gamma$ acts as a weight for
the penalization of the length of the interface between void and material.
Thus a smaller value of $\gamma$ is expected to lead to thinner structures
which contain more braces. Using the same setup as before, we solve again the
optimization problem for the cantilever beam, but this time we fix
$\varepsilon=0.02$. We perform two simulations with
$\gamma\in\\{10^{-4},10^{-5}\\}$. The smaller $\gamma$ is chosen, the finer
structures we expect. We also expect that we reach a larger value for
$\lambda_{1}$, because less regularization is used.

In Figure 2, we present the found structures for these parameters. On the left
we present the result for $\gamma=10^{-4}$ and on the right for
$\gamma=10^{-5}$. As expected, it is clearly visible that the structure
obtained for the smaller value of $\gamma$ is finer and contains more braces.
Additionally, decreasing $\gamma$ also leads to sharper corners.

Figure 2: The optimal beam for $\Psi(\lambda_{1})=-\lambda_{1}$, i.e.,
maximization of the principal eigenvalue for $\gamma=10^{-4}$ (left) and
$\gamma=10^{-5}$ (right) with $\varepsilon=0.02$. We clearly observe finer
structures for smaller $\gamma$. We obtain $\lambda_{1}=1.68\cdot 10^{-2}$ for
$\gamma=10^{-4}$ and $\lambda_{1}=1.72\cdot 10^{-2}$ for $\gamma=10^{-5}$.
Thus, as expected, with less regularization we reach a larger value for
$\lambda_{1}$.

In a second test for the beam setup, we compare the numerical results for
different choices of $\Psi(\lambda_{1},\lambda_{2})$ as a linear combination
of $\lambda_{1}$ and $\lambda_{2}$. We set $\gamma=10^{-4}$ and use the
solution shown in Figure 2 as the initialization of the optimization method.
In Figure 3, we present numerical results for this setting with the choice
$\Psi(\lambda_{1},\lambda_{2})=-\lambda_{1}-\alpha\lambda_{2}$ for
$\alpha\in\\{10^{-2},2\cdot 10^{-2},6\cdot 10^{-2},10^{-1}\\}$. Moreover, in
Table 2 we list the corresponding values of $\lambda_{1}$ and $\lambda_{2}$.
Here, $\alpha=0$ corresponds to the result shown in Figure 2 on the left.

Figure 3: Optimization of the cantilever beam for $\gamma=10^{-4}$ and $\Psi(\lambda_{1},\lambda_{2})=-\lambda_{1}-\alpha\lambda_{2}$, where $\alpha\in\\{10^{-2},{2\cdot 10^{-2}},{6\cdot 10^{-2}},10^{-1}\\}$ (left to right). We observe that increasing the weight of $\lambda_{2}$ above a certain value reduces the amount of fine structures.  $\alpha$ | 0 | $1\cdot 10^{-2}$ | $2\cdot 10^{-2}$ | $6\cdot 10^{-2}$ | $1\cdot 10^{-1}$
---|---|---|---|---|---
$\lambda_{1}$ | $1.677\cdot 10^{-2}$ | $1.662\cdot 10^{-2}$ | $1.606\cdot 10^{-2}$ | $1.521\cdot 10^{-2}$ | $1.508\cdot 10^{-2}$
$\lambda_{2}$ | $9.181\cdot 10^{-2}$ | $11.874\cdot 10^{-2}$ | $15.178\cdot 10^{-2}$ | $17.663\cdot 10^{-2}$ | $18.047\cdot 10^{-2}$
Table 2: The first and second eigenvalue ($\lambda_{1}$ and $\lambda_{2}$) for
the optimal topologies for the beam example and
$\Psi(\lambda_{1},\lambda_{2})=-\lambda_{1}-\alpha\lambda_{2}$. As expected,
for larger weights $\alpha$ we reach a lower value for $\lambda_{1}$ and a
larger value for $\lambda_{2}$.

###  9.3 Joint optimization of compliance and principal eigenvalue

In this subsection, we extend the problem by using a linear combination of
compliance and the first eigenvalue as objective. For any given $\varphi\in
H^{1}(\Omega)\cap L^{\infty}(\Omega)$, the compliance problem is to find a
displacement field $\boldsymbol{u}_{c}^{\varphi}\in
H^{1}(\Omega;\mathbb{R}^{d})$ satisfying

$\displaystyle-\nabla\cdot(\mathbb{C}(\varphi)\mathcal{E}(\boldsymbol{u}_{c}^{\varphi}))$
$\displaystyle=\boldsymbol{0}\quad\mbox{ in }\Omega,$ (9.4)
$\displaystyle\boldsymbol{u}_{c}^{\varphi}$
$\displaystyle=\boldsymbol{0}\quad\mbox{ on }\Gamma_{D}\subset\partial\Omega,$
$\displaystyle\left[\mathbb{C}(\varphi)\mathcal{E}(\boldsymbol{u}_{c}^{\varphi})\right]\cdot\boldsymbol{n}$
$\displaystyle=\boldsymbol{g}\quad\mbox{ on }\Gamma_{g}\subset\partial\Omega,$
$\displaystyle\left[\mathbb{C}(\varphi)\mathcal{E}(\boldsymbol{u}_{c}^{\varphi})\right]\cdot\boldsymbol{n}$
$\displaystyle=\boldsymbol{0}\quad\mbox{ on }\Gamma_{0}\subset\partial\Omega,$

which minimizes the objective
$\int_{\Gamma_{g}}\boldsymbol{g}\cdot\boldsymbol{u}_{c}^{\varphi}$.

Combining this with our eigenvalue optimization problem for
$\Psi(\lambda_{1})=-\alpha\lambda_{1}$ for some $\alpha>0$, we arrive at

$\left\\{\begin{aligned}
&\min&&J(\varphi)=-\alpha\lambda_{1}^{\varphi}+\int_{\Gamma_{g}}\boldsymbol{g}\cdot\boldsymbol{u}_{\boldsymbol{c}}^{\varphi}+\gamma
E^{\varepsilon}(\varphi)\\\ &\mbox{over }&&\varphi\in H^{1}(\Omega)\cap
L^{\infty}(\Omega),\\\ &\mbox{s.t.
}&&\boldsymbol{u}_{\boldsymbol{c}}^{\varphi}\mbox{ solves the compliance
equation \eqref{eq:num:compliance}},\\\ &&&\lambda_{1}^{\varphi}\mbox{ is the
first eigenvalue of \eqref{state}}.\end{aligned}\right.$ (9.5)

This means that we are looking for a structure that simultaneously minimizes
the compliance with respect to a given force $\boldsymbol{g}$ and maximizes
the first eigenvalue $\lambda_{1}$. The optimization problem (9.5) is actually
a special case of the compliance and eigenvalue optimization problems studied
in [45] (with the quantities therein being chosen as $N=2$, $\alpha=1$,
$\beta=0$, $\boldsymbol{f}\equiv\mathbf{0}$, $J_{0}\equiv 0$ and
$U_{c}=H^{1}(\Omega)$, i.e., $S_{0}=S_{1}=\emptyset$). More details about the
formulation of the problem (9.5) can be found in [45, Section 2.7]. For the
optimality system of (9.5), which we are going to solve numerically, we refer
to [45, Theorem 7.2].

For the sharp-interface limit of compliance optimization problems without
eigenvalue optimization (such as (9.5), where $\alpha$ is set to zero), we
refer to [20]. As our sharp-interface analysis for eigenvalue optimization
problems without compliance optimization relies on the same expansions as in
[20], both approaches can be combined to formally derive the sharp-interface
limit of problem (9.5).

For our numerical computations, we use the same setup as in Section 9.2 for
the beam example and fix $\gamma=1\cdot 10^{-3}$. Moreover, the exterior force
is $\boldsymbol{g}=(0,-1)^{T}$ and acts on $\Gamma_{g}=\\{(2.0,y)\mid
y\in(0.45,0.55)\\}$. Note that $\Gamma_{g}$ belongs to the boundary of the
domain $\Omega_{\rho}$ on which we assume a higher value of the density
$\rho$.

In Figure 4, we show numerical result for this setting for different values of
$\alpha$. We observe that the structures become finer when we increase the
influence of the principal eigenvalue. In Table 3, we present the
corresponding values for compliance and $\lambda_{1}$ for these shapes. As
expected, we achieve a larger compliance when we increase the weight $\alpha$
of the principal eigenvalue. Simultaneously, we also obtain larger values for
the principal eigenvalue. It is worth mentioning that these results compare
very well with the ones obtained in [7], where a level-set method was used to
directly tackle the sharp-interface problem (see especially Fig. 2 and Fig. 5
in [7]).

Figure 4: Numerical results for joint optimization of compliance and principal eigenvalue with weight $\alpha\in\\{10,100,500\\}$ (left to right). We observe that increasing the weight $\alpha$ of the first eigenvalue leads locally to a finer structure. $\alpha$ | 10 | 100 | 200 | 500
---|---|---|---|---
Compliance | 0.5507 | 0.5629 | 0.5676 | 0.5769
$\lambda_{1}$ | 0.0164 | 0.0170 | 0.0172 | 0.0173
Table 3: Values of compliance and principal eigenvalue $\lambda_{1}$ for joint
optimization of compliance and principal eigenvalue with weight
$\alpha\in\\{10,100,200,500\\}$. We observe that increasing $\alpha$ leads, as
expected, to larger values of the principal eigenvalue and larger values for
the compliance.

## Acknowledgment

Harald Garcke, Paul Hüttl and Patrik Knopf were partially supported by the RTG
2339 “Interfaces, Complex Structures, and Singular Limits” of the German
Science Foundation (DFG). The support is gratefully acknowledged.

## References

  * [1] Abels, H., Garcke, H., and Grün, G. Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities.  Math. Models Methods Appl. Sci. 22, 3 (2012), 1150013, 40.
  * [2] Abels, H., and Liu, Y. Sharp interface limit for a Stokes/Allen-Cahn system.  Arch. Ration. Mech. Anal. 229, 1 (2018), 417–502.
  * [3] Allaire, G. Shape Optimization by the Homogenization Method.  Applied Mathematical Sciences. Springer-Verlag, New York, 2002.
  * [4] Allaire, G., Aubry, S., and Jouve, F. Eigenfrequency optimization in optimal design.  Comput. Methods Appl. Mech. Engrg. 190, 28 (2001), 3565–3579.
  * [5] Allaire, G., Dapogny, C., Delgado, G., and Michailidis, G. Multi-phase structural optimization via a level set method.  ESAIM Control Optim. Calc. Var. 20, 2 (2014), 576–611.
  * [6] Allaire, G., and Jakabčin, L. Taking into account thermal residual stresses in topology optimization of structures built by additive manufacturing.  Math. Models Methods Appl. Sci. 28, 12 (2018), 2313–2366.
  * [7] Allaire, G., and Jouve, F. A level-set method for vibration and multiple loads structural optimization.  Comput. Methods Appl. Mech. Engrg. 194, 30-33 (2005), 3269–3290.
  * [8] Almi, S., and Stefanelli, U. Topology optimization for incremental elastoplasticity: a phase-field approach.  SIAM J. Control Optim. 59, 1 (2021), 339–364.
  * [9] Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M., and Wells, G. The FEniCS project version 1.5.  Archive of Numerical Software 3, 100 (2015).
  * [10] Antunes, P. R. S., and Oudet, E. Numerical results for extremal problem for eigenvalues of the Laplacian.  In Shape Optimization and Spectral Theory. De Gruyter Open, Warsaw, 2017, pp. 398–411.
  * [11] Auricchio, F., Bonetti, E., Carraturo, M., Hömberg, D., Reali, A., and Rocca, E. A phase-field-based graded-material topology optimization with stress constraint.  Math. Models Methods Appl. Sci. 30, 8 (2020), 1461–1483.
  * [12] Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W., Kaushik, D., Knepley, M., May, D., McInnes, L. C., Mills, R. T., Munson, T., Rupp, K., Sanan, P., Smith, B., Zampini, S., Zhang, H., and Zhang, H. PETSc users manual.  Tech. Rep. ANL-95/11 - Revision 3.9, Argonne National Laboratory, 2018\.
  * [13] Balay, S., Gropp, W. D., McInnes, L. C., and Smith, B. F. Efficient management of parallelism in object-oriented numerical software libraries.  In Modern Software Tools for Scientific Computing (1997), E. Arge, A. M. Bruaset, and H. P. Langtangen, Eds., Birkhäuser Press, pp. 163–202.
  * [14] Baldo, S. Minimal interface criterion for phase transitions in mixtures of Cahn-Hilliard fluids.  Ann. Inst. H. Poincaré Anal. Non Linéaire 7, 2 (1990), 67–90.
  * [15] Barrett, J. W., and Blowey, J. F. An error bound for the finite element approximation of a model for phase separation of a multi-component alloy.  IMA J. Numer. Anal. 16, 2 (1996), 257–287.
  * [16] Barrett, J. W., Garcke, H., and Nürnberg, R. On sharp interface limits of Allen-Cahn/Cahn-Hilliard variational inequalities.  Discrete Contin. Dyn. Syst. Ser. S 1, 1 (2008), 1–14.
  * [17] Beck, L. Elliptic regularity theory, vol. 19 of Lecture Notes of the Unione Matematica Italiana.  Springer, Cham; Unione Matematica Italiana, Bologna, 2016.  A first course.
  * [18] Bendsøe, M., and Sigmund, O. Topology optimization.  Springer-Verlag, Berlin, 2003.  Theory, methods and applications.
  * [19] Blank, L., Farshbaf-Shaker, M., Garcke, H., Rupprecht, C., and Styles, V. Multi-material phase field approach to structural topology optimization.  In Trends in PDE constrained optimization, vol. 165 of Internat. Ser. Numer. Math. Birkhäuser/Springer, Cham, 2014, pp. 231–246.
  * [20] Blank, L., Garcke, H., Farshbaf-Shaker, M. H., and Styles, V. Relating phase field and sharp interface approaches to structural topology optimization.  ESAIM Control Optim. Calc. Var. 20, 4 (2014), 1025–1058.
  * [21] Blank, L., Garcke, H., Hecht, C., and Rupprecht, C. Sharp interface limit for a phase field model in structural optimization.  SIAM J. Control Optim. 54, 3 (2016), 1558–1584.
  * [22] Blank, L., Garcke, H., Sarbu, L., Srisupattarawanit, T., Styles, V., and Voigt, A. Phase-field approaches to structural topology optimization.  In Constrained optimization and optimal control for partial differential equations, vol. 160 of Internat. Ser. Numer. Math. Birkhäuser/Springer Basel AG, Basel, 2012, pp. 245–256.
  * [23] Blank, L., Garcke, H., Sarbu, L., and Styles, V. Nonlocal Allen-Cahn systems: analysis and a primal-dual active set method.  IMA J. Numer. Anal. 33, 4 (2013), 1126–1155.
  * [24] Blank, L., and Rupprecht, C. An extension of the projected gradient method to a Banach space setting with application in structural topology optimization.  SIAM J. Control Optim. 55, 3 (2017), 1481–1499.
  * [25] Blowey, J. F., and Elliott, C. M. The Cahn-Hilliard gradient theory for phase separation with nonsmooth free energy. I. Mathematical analysis.  European J. Appl. Math. 2, 3 (1991), 233–280.
  * [26] Bourdin, B., and Chambolle, A. Design-dependent loads in topology optimization.  ESAIM Control Optim. Calc. Var. 9 (2003), 19–48.
  * [27] Bourdin, B., and Chambolle, A. The Phase-Field Method in Optimal Design.  In IUTAM Symposium on Topological Design Optimization of Structures, Machines and Materials (Dordrecht, 2006), M. Bendsøe, N. Olhoff, and O. Sigmund, Eds., Springer Netherlands, pp. 207–215.
  * [28] Bronsard, L., Garcke, H., and Stoth, B. A multi-phase Mullins-Sekerka system: matched asymptotic expansions and an implicit time discretisation for the geometric evolution problem.  Proc. Roy. Soc. Edinburgh Sect. A 128, 3 (1998), 481–506.
  * [29] Bucur, D., Martinet, E., and Oudet, E. Maximization of Neumann Eigenvalues.  Arch. Ration. Mech. Anal. 247, 2 (2023), 19.
  * [30] Burger, M. A framework for the construction of level set methods for shape optimization and reconstruction.  Interfaces Free Bound. 5, 3 (2003), 301–329.
  * [31] Burger, M., and Stainko, R. Phase-field relaxation of topology optimization with local stress constraints.  SIAM J. Control Optim. 45, 4 (2006), 1447–1466.
  * [32] Carraturo, M., Rocca, E., Bonetti, E., Hömberg, D., Reali, A., and Auricchio, F. Graded-material design based on phase-field and topology optimization.  Comput. Mech. 64, 6 (2019), 1589–1600.
  * [33] Ciarlet, P. G. Linear and nonlinear functional analysis with applications.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2013\.
  * [34] Dedè, L., Borden, M., and Hughes, T. Isogeometric analysis for topology optimization with a phase field model.  Arch. Comput. Methods Eng. 19, 3 (2012), 427–465.
  * [35] Delfour, M. C., and Zolésio, J.-P. Shapes and geometries, second ed., vol. 22 of Advances in Design and Control.  Society for Industrial and Applied Mathematics (SIAM),r Philadelphia, PA, 2011.  Metrics, analysis, differential calculus, and optimization.
  * [36] Dondl, P., Poh, P. S. P., Rumpf, M., and Simon, S. Simultaneous elastic shape optimization for a domain splitting in bone tissue engineering.  Proc. R. Soc. A. 475, 2227 (2019), 20180718, 17.
  * [37] Ebeling-Rump, M., Hömberg, D., Lasarzik, R., and Petzold, T. Topology optimization subject to additive manufacturing constraints.  Journal of mathematics in industry 11 (2021), 19.
  * [38] Ebeling-Rump, M., Hömberg, D., Lasarzik, R., and Petzold, T. Topology optimization subject to additive manufacturing constraints.  J. Math. Ind. 11 (2021), 19.
  * [39] Eck, C., Garcke, H., and Knabner, P. Mathematical modeling.  Springer Undergraduate Mathematics Series. Springer, Cham, 2017.
  * [40] Elliott, C. M., and Luckhaus, S. A generalised diffusion equation for phase separation of a multi-component mixture with interfacial free energy.  In Preprint SFB 256 University Bonn (1991), vol. 195.
  * [41] Fife, P. C. Dynamics of internal layers and diffusive interfaces, vol. 53 of CBMS-NSF Regional Conference Series in Applied Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1988.
  * [42] Garcke, H., Hecht, C., Hinze, M., Kahle, C., and Lam, K. F. Shape optimization for surface functionals in Navier-Stokes flow using a phase field approach.  Interfaces Free Bound. 18, 2 (2016), 219–261.
  * [43] Garcke, H., Hinze, M., Kahle, C., and Lam, K. F. A phase field approach to shape optimization in Navier- Stokes flow with integral state constraints.  Adv. Comput. Math. 44, 5 (2018), 1345–1383.
  * [44] Garcke, H., Hüttl, P., Kahle, C., Knopf, P., and Laux, T. Phase-field methods for spectral shape and topology optimization.  ESAIM Control Optim. Calc. Var. 29 (2023), 57.
  * [45] Garcke, H., Hüttl, P., and Knopf, P. Shape and topology optimization involving the eigenvalues of an elastic structure: a multi-phase-field approach.  Adv. Nonlinear Anal. 11, 1 (2022), 159–197.
  * [46] Garcke, H., Lam, K. F., Nürnberg, R., and Sitka, E. A multiphase Cahn-Hilliard-Darcy model for tumour growth with necrosis.  Math. Models Methods Appl. Sci. 28, 3 (2018), 525–577.
  * [47] Garcke, H., Lam, K. F., Sitka, E., and Styles, V. A Cahn-Hilliard-Darcy model for tumour growth with chemotaxis and active transport.  Math. Models Methods Appl. Sci. 26, 6 (2016), 1095–1148.
  * [48] Gilbarg, D., and Trudinger, N. S. Elliptic partial differential equations of second order.  Classics in Mathematics. Springer-Verlag, Berlin, 2001.  Reprint of the 1998 edition.
  * [49] Grisvard, P. Elliptic problems in nonsmooth domains, vol. 69 of Classics in Applied Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2011.  Reprint of the 1985 original.
  * [50] Henrot, A., and Pierre, M. Shape variation and optimization, vol. 28 of EMS Tracts in Mathematics.  European Mathematical Society (EMS), Zürich, 2018.  A geometrical analysis, English version of the French publication with additions and updates.
  * [51] Hernandez, V., Roman, J. E., and Vidal, V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems.  ACM Trans. Math. Software 31, 3 (2005), 351–362.
  * [52] Herzog, R., Meyer, C., and Wachsmuth, G. Integrability of displacement and stresses in linear and nonlinear elasticity with mixed boundary conditions.  J. Math. Anal. Appl. 382, 2 (2011), 802–813.
  * [53] Hilhorst, D., Kampmann, J., Nguyen, T. N., and Van Der Zee, K. G. Formal asymptotic limit of a diffuse-interface tumor-growth model.  Math. Models Methods Appl. Sci. 25, 6 (2015), 1011–1043.
  * [54] Hüttl, P., Knopf, P., and Laux, T. A phase-field version of the Faber–Krahn theorem.  Preprint:  arXiv:2207.10946 [math.AP], 2022.
  * [55] Kao, C. Y., Osher, S., and Yablonovitch, E. Maximizing band gaps in two-dimensional photonic crystals by using level set methods.  Applied Physics B 81, 2 (2005), 235–244.
  * [56] Kevorkian, J., and Cole, J. D. Multiple scale and singular perturbation methods, vol. 114 of Applied Mathematical Sciences.  Springer-Verlag, New York, 1996.
  * [57] Logg, A., Mardal, K.-A., and Wells, G., Eds.  Automated Solution of Differential Equations by the Finite Element Method - The FEniCS Book, vol. 84 of Lecture Notes in Computational Science and Engineering.  Springer, 2012.
  * [58] Marino, M., Auricchio, F., Reali, A., Rocca, E., and Stefanelli, U. Mixed variational formulations for structural topology optimization based on the phase-field approach.  Struct. Multidiscip. Optim. 64, 4 (2021), 2627–2652.
  * [59] Murat, F., and Simon, S. Etudes de problèmes d’optimal design.  In Lecturenotes in Computer Science (1976), vol. 41, Springer Verlag, Berlin, pp. 54–62.
  * [60] Osher, S., and Santosa, F. Level Set Methods for Optimization Problems Involving Geometry and Constraints: I. Frequencies of a Two-Density Inhomogeneous Drum.  Journal of Computational Physics 171, 1 (2001), 272 – 288.
  * [61] Osher, S., and Sethian, J. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations.  Journal of Computational Physics 79, 1 (1988), 12 – 49.
  * [62] Oudet, E. Numerical minimization of eigenmodes of a membrane with respect to the domain.  ESAIM Control Optim. Calc. Var. 10, 3 (2004), 315–330.
  * [63] Pedersen, N. L. Maximization of eigenvalues using topology optimization.  Structural and Multidisciplinary Optimization 20, 1 (2000), 2–11.
  * [64] Penzler, P., Rumpf, M., and Wirth, B. A phase-field model for compliance shape optimization in nonlinear elasticity.  ESAIM Control Optim. Calc. Var. 18, 1 (2012), 229–258.
  * [65] Shi, P., and Wright, S. Higher integrability of the gradient in linear elasticity.  Math. Ann. 299, 3 (1994), 435–448.
  * [66] Simon, J. Differentiation with respect to the domain in boundary value problems.  Numer. Funct. Anal. Optim. 2, 7-8 (1980), 649–687.
  * [67] Sokolowski, J., and Zolesio, J.-P. Introduction to Shape Optimization: Shape Sensitivity Analysis, vol. 16 of Springer Series in Computational Mathematics.  Springer-Verlag Berlin Heidelberg, 1992.
  * [68] Sternberg, P. Vector-valued local minimizers of nonconvex variational problems.  Rocky Mountain J. Math. 21, 2 (1991), 799–807.
  * [69] Zeidler, E. Nonlinear functional analysis and its applications. II/B.  Springer-Verlag, New York, 1990.  Nonlinear monotone operators, Translated from the German by the author and Leo F. Boron.