Shape and Topology Optimization involving the eigenvalues of an elastic structure: A multi-phase-field approach

A cost functional involving the eigenvalues of an elastic structure, that is described by a multi-phase-field equation, is optimized. This allows us to handle topology changes and multiple materials. We prove continuity and differentiability of the eigenvalues and we establish the existence of a global minimizer to our optimization problem. We further derive first-order necessary optimality conditions for local minimizers. Moreover, an optimization problem combining eigenvalue and compliance optimization is also discussed.


Introduction
The main goal in structural topology optimization is to find the optimal distribution of materials in a so called design domain. In contrast to shape optimization, topological changes including the design and distribution of holes in the structure are also allowed in topology optimization. The shape and the topology of the structure are initially unknown and are to be determined by minimizing a suitable objective functional. In many applications (e.g., design engineering), support conditions, volume restrictions, prescribed solid regions or voids, as well as applied loads have to be taken into account.
Several mathematical techniques to deal with shape or topology optimization problems can be found in the literature. The traditional approach is the method of boundary variations to compute shape derivatives. In this way the value of the objective functional can be decreased by deforming the boundary in a certain descent direction (see, e.g., [16,22,27,28] and the references cited therein). The drawbacks of this technique are its high computational costs and that topological changes are not allowed. In some situations, it is also possible to deal with changes of the topology by homogenization methods (see, e.g., [1]) or variants of this approach such as the SIMP method (see, e.g., [9]). Especially in recent times, the level-set method has been a popular tool to approach topology optimization problems. It was originally developed in [24] and was afterwards frequently used in the literature (see, e.g., [11,23]). Although this method can handle topological changes, difficulties can arise if voids are to be created.
In this paper, however, we pursue a different ansatz. We describe an elastic structure by a vector valued phase-field variable ϕ representing the distribution of materials. Here, with respect to the space variable the phase-field ϕ does not change its values abruptly but exhibits continuous phase transitions. The phase-field method for topology optimization was first introduced in [9] and was subsequently used, e.g., in [6][7][8]12,15,25,29,[31][32][33]. The main advantage of this approach is that topological changes can be handled directly without having to switch the framework. In particular, the creation of voids does not impose any problems. Moreover, these models are very well suited to be treated with methods from mathematical analysis.
In many technical applications (especially from engineering sciences), it is not only desired to optimize the material distribution of an elastic structure but also to minimize (or maximize) its dynamical response to a given driving frequency of a specific range. A classical example is a machine whose running engine generates vibrations that affect other parts of the device. Mostly, it is essential that the vibrations do not match certain eigenfrequencies of other components of the machine to avoid a resonance disaster.
A typical concrete example is discussed in [5] where the engine of an airplane is considered. While the engine is running it creates vibrations that affect other components of the aircraft, especially its wings. Therefore, they must be designed in such a way that these vibrations are not amplified. Otherwise this could lead to fatal flight instabilities and might even damage or break the material. From a mathematical point of view, the eigenmodes of the engine and the wings should be taken into account within the optimization process to keep them as different from the modes of the engine as possible. One possibility to guarantee this behavior would be to maximize the smallest eigenmode of the wings, as consequently also all larger eigenmodes will be separated from those of the engine (which are generally rather small).
Models for eigenvalue problems and their analysis have already been discussed in the literature, see [2,10,13,19,21,26,30]. In [2,30] the authors investigate models similar to the one we intend to study. In these papers, the density distribution ρ is assumed to depend only on the spatial variable x ∈ Ω meaning that the dependence on the structure (represented by the phase-field ϕ) is neglected. However, as the material distribution of the structure is actually to be optimized, the optimal density distribution is initially unknown.
In this paper, we study the following eigenvalue problem to describe an elastic structure: with the disjoint splitting ∂Ω = Γ D ∪ Γ 0 , where we additionally demand Γ D to have strictly positive Hausdorff measure. Here, C denotes the elasticity tensor, E(w) stands for the symmetrized gradient of w, λ ϕ is the eigenvalue (depending on ϕ) and w = w ϕ denotes a corresponding eigenfunction. In contrast, to the similar models studied in [2,30], the density distribution ρ(ϕ) is now allowed to depend on the phase-field ϕ. For more details about the notation and the motivation of this model, we refer the reader to Section 2.
We prove the existence of eigenvalues and eigenfunctions for the problem (1.1) and we establish essential properties needed for the theory of calculus of variations such as suitable continuity statements. This allows us to investigate an optimal control problem where an objective functional J ε l (ϕ) = Ψ(λ ϕ i 1 , . . . , λ ϕ i l ) + γE ε (ϕ), (1.2) is to be minimized under the constraint that ϕ and λ ϕ i j satisfy the state equation (1.1). Here, the function Ψ : (R >0 ) l → R is continuously differentiable and bounded from below, and penalizes the eigenvalues. The expression E ε (ϕ) stands for the Ginzburg-Landau energy where ε > 0 corresponds to the thickness of the diffuse interface and ψ stands for the bulk potential that usually has a double-well structure. A typical choice is ψ(s) = (s 2 − 1) 2 , s ∈ R. As the energy E ε is an approximation of the perimeter of the material boundaries, minimizing (1.2) can be related to a shape and topology optimization problem with a perimeter penalization (see, e.g., [4]). The phase-field ϕ represents the control and is supposed to satisfy suitable restrictions. For reasons of mathematical analysis we use a diffuse interface approach, i.e. the components of ϕ do not change their values abruptly but continuously at interfacial regions between the materials. As in [7], the sharp interface limit could be considered to describe a discrete material distribution and would allow it to become a formulation with a perimeter regularization. The optimization problem will be introduced in more detail in Section 2.
To derive first-order necessary conditions for locally optimal controls we need to show that the considered eigenvalues λ ϕ are differentiable with respect to ϕ. Based on the theory developed in [26], we can show that the smallest eigenvalue is semi-differentiable with respect to the phase-field. In addition, we prove Fréchet differentiability of simple eigenvalues by means of the implicit function theorem after introducing a proper sign convention for the eigenfunctions. A positive side benefit of the implicit function theorem is that we also obtain the Fréchet derivatives of the corresponding eigenfunctions.
With the approach described in this paper, it is also possible to approach classical shape and topology optimization problems also for other elliptic operators, i.e., the Laplacian. The idea is to send C and ρ to 0 in the void phase. In a forthcoming paper we plan to study this limit in more detail. To describe the idea of our approach in such a setting we consider a phase-field approximation for a spectral optimization problem for the Neumann-Laplace operator. Here, a domain D ⊂ R n is to be optimized such that the eigenvalues of the Neumann problem for the Laplace operator −∆u = λu in D, are maximal. Choosing a scalar phase-field variable ϕ and a function a such that a(−1) = δ, where δ denotes a small parameter, a(1) = 1 such that a is a smooth, positive interpolation in between. We then solve −∇ · (a(ϕ)∇u) = λa(ϕ)u in D, with a given function Ψ. We conjecture that for ε, δ → 0 classical spectral optimization problems for the Laplace operator are recovered. As´Ω ε 2 |∇ϕ| 2 + 1 ε ψ(ϕ) converges in the sense of a Γ-limit to the perimeter functional the limits contain a perimeter regularization. In a similar way (1.2) under the constraint (1.1) can be related to a sharp interface problem for eigenvalues of the elasticity operator.
Our paper is structured as follows. First, we precisely formulate the mathematical model for the problem with a special emphasis on the eigenvalue problem and its analytic difficulties. After the first continuity results for eigenvalues and eigenfunctions with respect to the phase-field ϕ, we are able to show existence of a minimizer of the objective functional. Here we do not need to assume anything about simplicity of eigenvalues. The most elaborate part is then dedicated to deriving differentiability results and to improve the aforementioned continuity statements, which will yield the desired variational inequality. In this context, it is crucial to assume simplicity of the considered eigenvalues. In the last part we want to combine the eigenvalue problem with compliance minimization problems.

Formulation of the problem
This section is devoted to the introduction of the mathematical model and the structural optimization problem.

The design domain and the phase-field variable
We fix a bounded Lipschitz design domain Ω ⊂ R d with d ∈ N whose boundary is split into two disjoint parts: A homogeneous Dirichlet boundary Γ D with strictly positive (d − 1)dimensional Hausdorff measure and a homogeneous Neumann boundary Γ 0 . The distribution of N ∈ N materials is represented by the vector valued phase-field ϕ : Ω → R N . This means, for any i ∈ {1, ..., N }, the component ϕ i can be interpreted as the concentration of the i-th material. In this regard, ϕ i = 0 describes the absence of the i-th material, whereas ϕ i = 1 means that only the i-th material is present. We use the convention that voids are also interpreted as a sort of material, whose distribution is given by the N -th component of the vector ϕ.
For reasons of mathematical analysis we use a diffuse interface approach, i.e. the components of ϕ do not change their values abruptly but continuously at interfacial regions between the materials. Furthermore, we want to prescribe the total spatial amount of each phase. To this end, we impose the mean value constraint ffl is a fixed given number for any i ∈ {1, . . . , N }. In addition, we want the vector m to be an element of the set This constraint is a plausible consequence of the physical fact that at each point in space the volume fractions of the materials should sum up to 1. Furthermore, being a volume fraction, each component clearly has to be non-negative. For the upcoming analysis, we additionally want to prescribe a suitable regularity for the phase-field, namely H 1 (Ω; R N ). All these constraints are expressed in the set Here, G is given by The set G is referred to as the Gibbs-Simplex.

The Ginzburg-Landau energy
For the objective functional and especially the well-posedness of the minimization problem the following so called Ginzburg-Landau energy is crucial. In our model the function ψ : R N → R ∪ {∞} should attain exactly N global minima of value zero attained at the unit vectors e i ∈ R N , i.e. for all i ∈ {1, ..., N }, Furthermore, ψ is assumed to exhibit the decomposition ψ(ϕ) = ψ 0 (ϕ) + I G (ϕ) with ψ 0 ∈ C 1,1 (R N , R) and I G being the indicator functional This type of obstacle functional is used to enforce that ϕ attains its values only in G as this set is not penalized by the indicator functional. We refer to Elliot and Luckhaus [18] who first introduced this obstacle formulation of the energy E ε .

The density function
The density distribution ρ depends directly on the phase-field ϕ and this way, ρ is not just a given function but represents the density of the actual structure we want to optimize.
To this end, we assume that the density function ρ belongs to C 1,1 loc (R N ; R) and is uniformly positive, i.e., there exists a constant ρ 0 > 0 such that ρ(ϕ) ≥ ρ 0 for all ϕ ∈ R N . This directly yields for all ϕ, u ∈ R N . For any fixed ϕ ∈ R N , there exist constants C ϕ , C ′ ϕ > 0 (that may locally depend on ϕ, i.e., C ϕ and C ′ ϕ can be chosen uniformly on bounded sets), such that Due to the above assumptions, we can use this notation to define a family of scalar products on L 2 (Ω; R d ) depending on ϕ ∈ L ∞ (Ω; R N ) by These scalar products canonically induce norms that are all equivalent to the standard norm on L 2 (Ω; R d ). To indicate the norm we consider L 2 (Ω; R d ) to be equipped with, we will use the notation L 2 ϕ (Ω; R d ). A reasonable choice of ρ would be Here, for any i ∈ {1, ..., N − 1}, the coefficient ̺ i > 0 stands for the density of the i-th material which is assumed to be constant. In our model we interpret the void as a material of very low density. Hence, we chose ̺ N = ̺ N (ε) = ε 2̺ N as corresponding density, wherẽ ̺ N > 0 is a fixed constant. The scaling with ε 2 would guarantee the desired behaviour of ρ in the sharp interface limit, see [7] who treat the sharp interface limit for a related problem.
However, for the sake of mathematical analysis, we have to extend the definition of ρ onto the hyperplane Σ N . To this end, we define the cut-off function for any δ > 0 which will be specified later. Here, a δ and b δ are monotonically increasing C 1,1 -functions such that σ δ ∈ C 1,1 (R; R). We now define the function ρ by Obviously, it holds that ρ ∈ C 1,1 (R N ; R) and since σ δ (ϕ i ) = ϕ i for all i ∈ {1, ..., N } as long as ϕ ∈ G, the relation (2.4) holds true for this definition.
It remains to show that ρ is uniformly positive, at least if δ is chosen sufficiently small. To this end, we fix an arbitrary vector ϕ ∈ Σ N and define the index sets Recalling the definition of Σ N , we infer that I ≥0 ϕ i ≥ 1, and thus also we conclude the estimate Since ϕ ∈ Σ N was arbitrary, this estimate holds for all ϕ ∈ Σ N . We point out that ρ 0 does not depend on ϕ and thus, this estimate is uniform. This means that the function ρ defined in (2.6) is admissible as it exhibits all demanded properties.

The elasticity tensor
Another important tool in linear elasticity are the tensors appearing in Hooke's Law (see, e.g., [17,20]), namely the strain and the elasticity tensor which describe the stress tensor. To introduce the strain tensor we consider the displacement vector u : Ω → R d that describes the deformation of the structure under applied forces or vibrations. Now, the strain tensor of u can be defined as where A sym := 1 2 A + A T for any matrix A ∈ R d×d . The elasticity tensor C is a fourth order tensor whose components are demanded to fulfill C ijkl ∈ C 1,1 loc (R N , R) as well as the symmetry properties for all i, j, k, l ∈ {1, . . . , d}. From the regularity property we conclude that for any ϕ ∈ R N , there exist constants Λ ϕ , Λ ′ ϕ > 0 locally depending on ϕ such that . denotes the derivative of C(ϕ) in the direction h. Furthermore, we demand that there exists a positive constant θ such that for all symmetric matrices A ∈ R d×d \ {0} and for all ϕ, h ∈ R N it holds Recall that the application of a fourth order tensor onto a quadratic matrix is given by A concrete choice of the elasticity tensor in analogy to the construction of ρ is where for i = 1, . . . , N − 1, C i andC N denote constant material specific elasticity tensors.
To guarantee (2.9) we need to assume the existence of positive constantsθ i , ϑ i such that for all A ∈ R d×d \ {0}, it holds that for all i = 1, . . . , N . Now, proceeding similarly as for the density ρ, we can construct an extension to R N taking (2.9) into account. For more details we refer to [7, Sect. 2.2].

The system of PDEs describing the elastic structure
We now introduce the system of equations describing the elastic structure: Endowed with the standard inner product and norm given by is a Hilbert space. For any matrices A, B ∈ R d×d and any fourth-order tensor C ∈ R d×d×d×d , we introduce the notation Then the mapping defines a scalar product on H 1 D (Ω; R d ). By Korn's inequality (see, e.g., [35]), the norm induced by this inner product is equivalent to the standard norm on H 1 In what follows, we will always choose for a given ϕ this inner product and induced norm on H 1 D (Ω; R d ). Using this notation and invoking the symmetry property (2.7), the weak formulation of (2.10) can be expressed as In Section 3, we will see that for any ϕ ∈ L ∞ (Ω, R N ), there exists a sequence of eigenvalues and corresponding eigenfunctions {w ϕ 1 , w ϕ 2 , ...} ⊂ H 1 D (Ω; R d ) which form an orthonormal basis of L 2 ϕ (Ω; R d ). We also want to mention that in [26] such eigenvalue problems are analyzed on a general, rather abstract level. The eigenvalue problem therin is given by the relation where λ ϕ and w(ϕ) stand for the corresponding eigenvalue and eigenfunction, respectively, H is a Hilbert space and η ∈ H can be interpreted as a test function. For the analysis, it is assumed that a and b are of the form where A and B are linear, continuous operators, B is compact, and (·, ·) H denotes the inner product on V . For this problem, continuity and (semi-)differentiability of λ ϕ and w(ϕ) with respect to ϕ is established using an approach involving inverse operators that differs from the one discussed in this paper. However, optimization problems are not addressed in [26]. The concept of "semi-differentiability" applied to our setting will be explained in Section 5.1 in more detail.
Next, we introduce the structural optimization problem in which the system (2.10) can be regarded as the state equation.

The structural optimization problem
We also want to introduce constraints on the structure which prescribe void or material in certain parts of Ω. Mathematically speaking, we fix two disjoint measurable sets S i ⊂ Ω with i ∈ {0, 1} and define the set to fix material on S 0 and complete void on S 1 .
For l ∈ N and i 1 , . . . , i l ∈ N, the eigenvalues λ i 1 , ..., λ i l are to be penalized via a function which is assumed to be C 1 and bounded from below, i.e., we find a constant c Ψ > 0 such that Ψ(x) ≥ −c Ψ for all x ∈ (R >0 ) l . As mentioned above, we have to include the Ginzburg-Landau energy into our minimization problem. Hence we define the objective functional as with γ > 0. Consequently, the overall optimization problem reads as . . , λ ϕ i l are eigenvalues of (2.12).
(P ε l ) To investigate this optimal control problem, we first have to establish the existence of eigenvalues along with suitable associated eigenfunctions. This topic is addressed in the next section.

A combination of compliance and eigenvalue optimization
In [7], the problem of minimizing the mean compliance with f ∈ L 2 (Ω; R d ) and g ∈ L 2 (Γ g , R d ), and the deviation with respect to a target displacement u Ω ∈ L 2 (Ω; R d ) given by is also considered. Here, c ∈ L ∞ (Ω) denotes a function with |supp c| > 0, where |supp c| stands for the Lebesgue measure of the support. The boundary ∂Ω is split into two relatively open, disjoint subsets Γ C , Γ g ⊂ ∂Ω such that ∂Ω = Γ C ∪ Γ g . Moreover, the state equation is determined by the mean compliance in order to obtain as the displacement vector under the given forces. It reads as Combining this problem with the one discussed in Subsection 2.6., we obtain a structure that is on the one hand as stiff as possible (i.e., it has small compliance) and on the other hand realizes the desired vibration properties (e.g., a large first eigenvalue). In Section 7, we will present an existence result as well as the variational inequality for this combined problem. Then the combination of (P ε ) in [7] and (P ε l ) reads as and λ ϕ i 1 , . . . , λ ϕ i l are eigenvalues of (2.12),

Analysis of the state equation
In this case, the function w ϕ is called an eigenfunction to the eigenvalue λ ϕ .
The assumptions of the previous section allow us to prove two classical functional analytic results in our setting.
(a) There exists a sequence possessing the following properties: • For all k ∈ N, w ϕ k is an eigenfunction to the eigenvalue λ ϕ k in the sense of Definition 3.1.
• The eigenvalues λ ϕ k (which are repeated according to their multiplicity) can be ordered in the following way: Moreover, it holds that λ ϕ k → ∞ as k → ∞, and there exist no further eigenvalues of the state equation (3.1).

Moreover, the maximum is attained at the subspace
Proof. Using the Lax-Milgram theorem and the fact that This allows us to define a solution operator we can easily show that T is a compact, self-adjoint, and bounded linear operator. Thus, the assertions in (a) directly follow from the spectral theorem for compact self-adjoint operators (see, e.g., [3,Sect. 12.12]).
To prove (b), we first infer from (3.1) that the sequence forms an orthonormal basis of H 1 D (Ω; R d ) when taking the inner product (2.11). For any where the series on the right-hand side converges in H 1 D (Ω; R d ). In the following, we will sometimes omit the exponent ϕ for a more convenient depiction.
To establish the Courant-Fischer characterization we now fix an arbitrary subspace V ∈ S k−1 . Let us denote the orthogonal projection from L 2 ϕ (Ω; R d ) to V with respect to the scalar product on L 2 ϕ (Ω; R d ) by must be linearly dependent. Hence, for every i ∈ {1, . . . , k}, we find coefficients α i ∈ R that are not all equal to zero such that Per construction of the orthogonal projection this is equivalent to As not all of the coefficients vanish, and since the eigenfunctions {w 1 , . . . , w k } are linearly independent, we know that v = 0. Using the orthogonality of eigenfunctions and the fact, that the sequence of eigenvalues increases, we conclude that As the infimum in (3.2) obviously exists, we can find a minimizing sequence Hence, due to the Banach-Alaoglu theorem and the compact embedding along a non-relabeled subsequence. In particular, since all members of the sequence (ũ l ) l∈N are normalized with respect to the L 2 Hence, it is also weakly (sequentially) closed and we thus know that Using the fact that norms are always weakly lower semi-continuous we infer thatũ is a minimizer of the expression in (3.2) via the direct method in the calculus of variations.
Since this holds for any arbitrary We now select a special (k − 1)-dimensional subspace defined by Then the definition of the orthogonal complement yields where the series converges in H 1 D (Ω; R d ). Consequently, we obtain Altogether, we conclude that This means that the maximum is attained at the subspace , which proves the claim.

Weak sequential continuity of the eigenvalues
First of all we only consider the first eigenvalue λ 1 to establish continuity results with respect to the phase-field ϕ. Afterwards, we proceed inductively to obtain these results also for all the other eigenvalues.
We consider the mapping associated with the first eigenvalue.
The first continuity result for the eigenvalue λ 1 is obtained by proving lower and upper semi-continuity. Lower semi-continuity is established by the following lemma.
Then it holds that along a non-relabeled subsequence.
Proof. We first notice that the assumptions of Lemma 4.1 imply that ϕ ∈ L ∞ (Ω; R N ). Let denote an orthonormal basis of eigenfunctions corresponding to the sequence of eigenvalues λ ϕ i i∈N from Theorem 3.2.
Now, for any k ∈ N, we choose an arbitrary L 2 ϕ k (Ω; R d )-normalized eigenfunction u k that fulfills (3.1) for λ ϕ k 1 . This choice is not necessarily unique up to multiplication with ±1, as we do not assume simplicity of λ ϕ k 1 or λ ϕ 1 yet. Using the Courant-Fischer representation from Theorem 3.2(b) for the first eigenvalue and the continuity of C and ρ, we see that the sequence (u k ) k∈N ⊂ H 1 D (Ω; R d ) is bounded. By the Banach-Alaoglu theorem and the compact embedding as k → ∞, up to a subsequence. With the help of Lebesgue's theorem and the assumptions on the sequence ϕ k , this yields as k → ∞ after another subsequence extraction. This implies u L 2 ϕ (Ω;R d ) = 1 since the members u k were chosen as L 2 ϕ k (Ω; R d )-normalized eigenfunctions. In particular, this implies that , and invoking the increasing order of the sequence λ ϕ i i∈N , we conclude that If we can now show that the proof would be complete. Using the convergence results we have just established, the Cauchy-Schwarz inequality and the weak formulation (3.1), we infer that We further know that Using Lebesgue's convergence theorem, the boundedness of (u k ) k∈N , the local Lipschitz continuity of C and the assumptions on (ϕ k ) k∈N , we conclude that the first summand converges to zero along a non-relabeled subsequence. The second summand converges to zero as a direct consequence of (4.1).
In summary, we obtain that which completes the proof. Now, we establish the corresponding result for weak upper semi-continuity.
Then it holds that along a non-relabeled subsequence.
Proof. As C and ρ satisfy suitable continuity properties we can proceed as in [21,Thm. 8.1.3] and use once more the Courant-Fischer representation for the first eigenvalue to prove the claim.
Combining both lemmata we can conclude that λ 1 is weakly sequentially continuous.
Then it holds that

4)
i.e., the whole sequence of eigenvalues converges and not just a subsequence.
Proof. The convergence λ ϕ k 1 → λ ϕ 1 as k → ∞ follows from Lemma 4.1 and Lemma 4.2 after extraction of a subsequence. Moreover, as the limit λ ϕ 1 does not depend on the choice of the subsequence, we conclude by a standard contradiction argument that the convergence remains true for the whole sequence.
The convergence properties of (u k ) k∈N ⊂ H 1 D (Ω; R d ) and the fact that the weak limit for all test functions η ∈ H 1 D (Ω; R d ). Using the convergence of eigenvalues (4.4) and proceeding as in the proof of Lemma 4.1, we infer that for any any η ∈ H 1 as k → ∞, after extraction of a subsequence. Hence, we can pass to the limit in equation (4.5) to obtain which proves that u ∈ H 1 D (Ω; R d ) is indeed an eigenfunction corresponding to λ ϕ 1 .
Corollary 4.3 now serves as initial case for the following inductive proof which yields convergence of all eigenvalues.
as k → ∞ along a non-relabeled subsequence.
Proof. As mentioned before we proceed by induction. The initial step has already been established in Corollary 4.3.
Now, we assume that the statement is already verified for the index (j − 1) ∈ N. Our task is to prove that the assertion is true for the j-th eigenvectors and the associated eigenfunctions.
In this regard, the Courant-Fischer representation of Theorem 3.2(b) will be a helpful tool.
For k ∈ N, we fix the (j − 1)-dimensional subspace of H 1 D (Ω; R d ) that realizes the maximum in the Courant-Fischer representation discussed in Theorem 3.2(b), namely Analogously, we define Then by the induction hypothesis we know that for every i = 1, . . . , j − 1, there exists a for m = l, using the convergence properties of the sequence (ϕ k ) k∈N along with Lebesgue's convergence theorem. In particular, the family {u 1 , . . . , u j−1 } ⊂ L 2 ϕ (Ω; R d ) is linearly independent, which yields that all eigenfunctions to eigenvalues strictly smaller than λ ϕ j are contained in its span W := u 1 , . . . , u j−1 span . Hence, we conclude that As the minimum is attained, we infer that we find a non-trivial Otherwise the inequality in (4.8) would be strict, which would be a contradiction to Theorem 3.2. This means we have shown the existence of a function fulfilling (4.10).
Let now the sequence (v k ) k∈N be defined by for all k ∈ N. By this construction, we immediately observe that We now intend to show that the convergences as k → ∞, hold along a non-relabeled subsequence.
To verify (4.13), we consider the decomposition .

(4.15)
For the first product on the right-hand side, we directly obtain the convergence along a suitable subsequence. As the functions w ϕ k i are L 2 ϕ k (Ω; R d )-normalized eigenfunctions, we obtain from (3.1) the following representation of the third product on the righthand side of (4.15): As the sum takes only the indices i = 1, . . . , j − 1 into account, we can again use the induction hypothesis to obtain along a suitable subsequence. Hence, (4.9) directly yields that the third product on the right-hand side of (4.15) converges to zero. The second product can be handled similarly, and we can also show that it tends to zero as k → ∞. In summary, we get This proves (4.13). The claim (4.14) can easily be verified using the induction hypothesis.
In particular, since v = 0, we obtain that v k = 0 for all k ∈ N sufficiently large. For such k ∈ N, we obtain the estimate .
Using (4.13) and (4.14), we conclude from (4.10) that, along a non-relabeled subsequence, In particular, this implies that the subsequence (λ ϕ k j ) k∈N is bounded. Now, for k ∈ N, let u k j denote a L 2 ϕ k (Ω; R d )-normalized eigenfunction to the eigenvector λ ϕ k j . Consequently, due to the Courant-Fischer characterization in Theorem 3.2(b), the sequence (u k j ) k∈N ⊂ H 1 D (Ω; R d ) is bounded. Applying the Banach-Alaoglu theorem, we can thus extract a subsequence such that However, it is a priori not necessarily an eigenfunction to the eigenvalue λ ϕ k j , as the convergence of the corresponding eigenvalues is still unknown.
Proceeding as in Subsection 4.1., we want to show that To verify (4.18), we first observe that due to the orthogonality of eigenfunctions corresponding to different eigenvalues u k j , w ϕ k m ρ(ϕ k ) = 0, for all k ∈ N and for all m = 1, . . . , j * , where j * < j − 1 is the maximal index such that λ ϕ j * < λ j−1 . Recalling the assumptions on (ϕ k ) k∈N , we can use (4.17) and (4.6) to infer that w, u m ρ(ϕ) = 0, for all m = 1, . . . , j * . As j * is chosen maximally, we know from the orthogonality (4.7) that This leads to the representation As the series converges in H 1 D (Ω; R d ), we can use (3.1) to obtain This, however, can be proven completely analogously as in the proof of Lemma 4.1.
In summary, we obtain the convergence along a non-relabeled subsequence. Since the limit does not depend on any subsequence extraction, this convergence holds true for the whole sequence. As in Corollary 4.3, we conclude that u j := w ∈ H 1 D (Ω; R d ) is an eigenfunction to the eigenvalue λ ϕ j . In view of (4.17), this completes the proof.

Local Lipschitz continuity of the eigenvalues
The following lemma shows that all eigenvalues are locally Lipschitz continuous with respect to ϕ.

This means that the mapping
is locally Lipschitz continuous.
where the last inequality holds due to the local Lipschitz continuity of C and ρ, and the boundedness of λ ϕ+h i which follows from Theorem 4.4. Note that the constant C ϕ i may depend on λ ϕ i but not on the eigenfunctions we have chosen, as they were assumed to be normalized.

Suppose now that there exists a zero sequence (h
as k → ∞. For the corresponding sequence of eigenfunctions (w ϕ+h k i ) k∈N ⊂ H 1 D (Ω; R d ) for the eigenvalues (λ ϕ+h k i ) k∈N , we know from Theorem 4.4 that we find a L 2 ϕ (Ω; R d )-normalized eigenfunction w to the eigenvalue λ ϕ i such that as k → ∞, up to subsequence extraction. In particular, for k sufficiently large, we know that the members of this subsequence satisfy and thus, which is an obvious contradiction. This proves the claim.

A sign convention for the eigenfunctions
In the previous analysis there was no need to assume that the eigenspaces are one-dimensional. However, in Section 5, we want to show that the eigenvalues are Fréchet differentiable with respect to the phase-field. Therefore, it will be necessary to assume that for fixed ϕ ∈ H 1 (Ω; R N ) ∩ L ∞ (Ω; R N ) the eigenspace corresponding to the considered eigenvalue λ ϕ i is one-dimensional. In this case the eigenvalue is called simple.
Simplicity of λ ϕ i allows us to choose a corresponding eigenfunction w ϕ i ∈ H 1 D (Ω; R d ) that is normalized with respect to the scalar product on L 2 ϕ (Ω; R d ) and unique up to multiplication by ±1. We call such an eigenfunction a representative corresponding to λ ϕ i . In general, any eigenspace could be higher dimensional. For numerically motivated examples showing that even the simplicity of the first eigenvalue of a scalar elliptic regular PDE is no longer fulfilled in the vector valued case, see [14]. However, in concrete applications, there are physical and numerical justifications for assuming simple eigenvalues. This is due to the fact that nature as well as numerical simulations on computers lead to perturbations of the non-generic case of equal eigenvalues.
As a classical two dimensional example to illustrate this behavior, an eigenvalue problem associated with the Laplacian subject to Dirichlet boundary conditions can be considered. If the domain is a perfect circle, eigenvalues with higher multiplicity will occur. However, as soon as the perfect circular shape of the domain is perturbated by small imperfections, these eigenvalues will become different and simple. For more details see [30].
In the following lemma, we will introduce a condition to fix a sequence of representatives whose elements w is sufficiently close to ϕ. In particular, we see that it is possible to deduce simplicity of the eigenvalues λ ϕ k i in a suitable neighborhood of λ ϕ i . Lemma 4.6. Let i ∈ N and (ϕ k ) k∈N ⊂ H 1 (Ω; R N ) ∩ L ∞ (Ω; R N ) be a sequence such that for k → ∞. Moreover, we assume that λ ϕ i is a simple eigenvalue of (3.1) and let w ϕ i be a corresponding L 2 ϕ (Ω; R d )-normalized eigenfunction. Then for any ε ∈ (0, 1), we can find a K ε i > 0 such that for any k > K ε i , there exists a unique In particular, the eigenvalues λ Proof. Note that we did not make any assumptions on the simplicity of the eigenspaces corresponding to λ ϕ k i for k ∈ N. However, this can be established if ϕ k is close to ϕ by invoking the simplicity of the eigenspace corresponding to λ ϕ i and using the continuity properties known from Theorem 4.4.
In the following, we will assume, without loss of generality, that k is large enough to ensure that all eigenspaces to the eigenvalues λ ϕ k i are simple. If we are now able to find a sequence of representatives w ϕ k i that fulfills (4.19) for a suitable K ε i ∈ N, then the uniqueness assertion is clear since the eigenfunctions are normalized and their sign is fixed by (4.19).
To prove the existence of such a sequence, we argue once more by contradiction. Let ε ∈ (0, 1) be arbitrary and let us assume that there is no K ε i ∈ N such that (4.19) is fulfilled. Hence, after possibly swapping some of the signs, we can extract a subsequence such that Using Theorem 4.4 we obtain a weak limit w of a non-relabeled subsequence of w ϕ k i k∈N and infer from the simplicity of λ ϕ i that w = ±w ϕ i . Hence, using (4.20), we obtain which is obviously a contradiction.
Eventually, this means that condition (4.19) allows us to pick a unique representative w ϕ k i for every k ∈ N sufficiently large such that the obtained sequence fulfills as k → ∞, up to subsequence extraction.
The following corollary is a direct consequence of Lemma 4.6.

Continuity of the eigenfunctions
In view of the sign convention from Corollary 4.7, we can now prove the following continuity result.
. This means that the mapping

be any arbitrary sequence satisfying
we can apply Theorem 4.4 to conclude that along a non-relabeled subsequence. However, as the limit does not depend on the extracted subsequence, this convergence even holds true for the whole sequence. Note that for this reasoning it is essential that all members of the sequence are fixed by the sign convention (4.21). As the sequence (h k ) k∈N was arbitrary, we further infer that Here, the first summand converges to zero because of (4.24). The second summand converges to zero due to the local Lipschitz continuity of ρ and the last summand converges to zero as a consequence of Theorem 4.4. This verifies (4.25) and thus, the proof is complete.

A formal consideration
First of all, we want to discuss the desired differentiability results formally. To obtain the Fréchet derivative of the functional for i ∈ N, we formally differentiate the state equation in the Gâteaux sense. If w is an eigenfunction to the eigenvalue λ ϕ i , we have Computing the first variation of (5.1) with respect to ϕ in the direction h ∈ H 1 (Ω; R N ) ∩ L ∞ (Ω; R N ), and choosing w and η as the L 2 ϕ (Ω; R d )-normalized eigenfunction w ϕ i afterwards, we get Moreover, firstly plugging w = w ϕ i into (5.1), and then computing the first variation with respect to ϕ in the direction h ∈ H 1 (Ω; R N ) ∩ L ∞ (Ω; R N ) reveals that the formal derivative w ϕ i ′ h of w ϕ i has to fulfill the equation In the following, we intend to verify these results rigorously. We already see that formula (5.2) is a priori not well-defined if there are at least two orthogonal eigenfunctions to the eigenvalue λ ϕ i , i.e., if λ ϕ i is not simple. In the following approach, we will see that the simplicity of eigenvalues will play a crucial role in our analysis.

Semi-differentiability of the first eigenvalue
In [26,Sect. 4.2], the concept of semi-differentiability is introduced and applied to the first eigenvalue of an abstract problem discussed there. Semi-differentiability is a concept similar to Gateâux-differentiability but the limit does not need to fulfill any linearity or continuity assumptions, and the variation is only performed along a fixed positive direction. The advantage becomes clear by the following example presented in [21,Sect. 2.5]. We consider the matrix-valued function whose first eigenvalue is not simple at t = 0. Of course, λ 1 is not classically differentiable in t = 0, but we still obtain a well defined limit lim t→0 t>0 This means we can still compute some sort of derivative in a fixed positive direction.
We now give a precise definition of semi-differentiability which can be found, e.g., in [26,Def. 4.6].

Definition 5.1 (Definition of semi-differentiability). Let X, Y be Banach spaces and let D ⊆ X be an open subset. Then the map
In this case we write T ′ (x)h = y(x, h) to denote the semi-derivative of T at the point x with respect to the direction h.
We want to show that in our problem the first eigenvalue also fulfills this weaker notion of differentiability, which will be enough to deduce first-order necessary optimality conditions for the first eigenvalue as we only want to derive convex combinations where t > 0.
The advantage of this approach is that we do not have to assume simplicity of the first eigenvalue in order to obtain semi-differentiability, whereas as illustrated in Section 5.1., we need such simplicity assumptions to obtain classical differentiability.
The semi-differentiability of the first eigenvalue is established by the following lemma.
Theorem 5.2 (Semi-differentiability of the first eigenvalue). Let ϕ, h ∈ H 1 (Ω; R N ) ∩ L ∞ (Ω; R N ) be arbitrary and let us define and thus, the eigenvalue λ ϕ 1 is semi-differentiable with respect to ϕ.
Proof. We first prove that the infimum in (5.4) is actually attained by a minimizer. To this end, let u ∈ F ad be arbitrary, where the feasible set is given as By the differentiability and the local Lipschitz continuity of C, we infer that there exists a constant c ϕ > 0 and t 0 > 0 such that for all t < t 0 , Using (3.1), we conclude that Moreover, due to (2.2) and (2.3), there exists a constant c * ϕ > 0 such that Eventually, combining the above estimates, we conclude that for all u ∈ F ad . This directly implies that the infimum (λ ϕ 1 ) ′ h exists. Hence, we can find a minimizing sequence (u n ) n∈N ⊂ F ad such that Due to (5.6), there exists u * ∈ H 1 D (Ω; R d ) such that as n → ∞, up to a subsequence extraction. In particular, this implies that u * ∈ F ad which leads to This implies that u n → u * even strongly in H 1 D (Ω; R d ). In particular we obtain Hence, the infimum is attained at it suffices to show that there exist functions f, g : R → R with f, g ∈ o(t) as t → 0 such that for all t > 0, By the construction of u * , we first observe that since u * is an L 2 ϕ (Ω; R d )-normalized eigenfunction to the eigenvalue λ ϕ 1 . We compute Now the first three summands on the right-hand side are clearly in o(t) as t → 0 since the eigenvalues converge, and the functions ρ and C are differentiable and locally Lipschitz continuous. For the remaining summands we can use the Courant-Fischer representation for the first eigenvalue which yields .
This implies and thus, (5.9) is established.
To prove (5.10), we argue by contradiction and assume that (5.10) does not hold. Then, there exists ε > 0 and a sequence (t k ) k∈N ⊂ (0, 1] with t k → 0 as k → ∞ such that for all k ∈ N, Then, according to Theorem 4.4, there exists a L 2 ϕ (Ω; R d )-normalized eigenfunction u to the eigenvalue λ ϕ 1 , as well as a sequence (u ϕ+t k h ) k∈N consisting of L 2 ϕ+t k h (Ω; R d )-normalized eigenfunctions to the eigenvalues (λ ϕ+t k h 1 ) k∈N such that as k → ∞, along a non-relabeled subsequence. Recalling the definition of λ ϕ Recalling the convergence property (5.11), that ρ and C are of class C 1 loc , that ϕ → λ ϕ 1 is locally Lipschitz continuous according to Lemma 4.5, and that we conclude that the right-hand side belongs to o(t k ) as k → ∞.
On the other hand we assumed which is obviously a contradiction as the inequality cannot hold for k sufficiently large. This proves (5.10). Now, (5.8) directly follows from (5.9) and (5.10) and thus, the proof is complete.

Fréchet differentiability of eigenvalues and their corresponding eigenfunctions
If the considered eigenvalue is simple, we can even obtain stronger differentiability results in the Fréchet sense. To be precise, if for i ∈ N and ϕ ∈ H 1 (Ω; R N ) ∩ L ∞ (Ω; R N ), the eigenvalue λ ϕ i associated with ϕ is simple, then λ ϕ i and any fixed L 2 ϕ (Ω; R d )-normalized eigenfunction w ϕ i are even Fréchet-differentiable with respect to ϕ. This is established by the following theorem: is well-defined and continuously Fréchet differentiable. Here, w ϑ i denotes the unique L 2 ϑ (Ω; R d )normalized eigenfunction to the eigenvalue λ ϑ i satisfying the sign condition (4.21) written for h = ϑ − ϕ. (5.12) and the Fréchet derivative w ϕ

Moreover, for any
To prove of this theorem, we intend to apply the implicit function theorem (see, e.g., [34,Theorem 4.B]). Therefore, it is essential to show a bijectivity condition. In our setting this condition will be fulfilled if a certain PDE resulting from the eigenvalue equations has a unique solution. To show this existence and uniqueness we need to apply the Fredholm alternative established by Lemma 5.4.
In the following, we use the space along with the canonical embedding In particular, the duality pairing is given by ·, · H −1 ,H 1 = (·, ·) ρ(ϕ) . Lemma 5.4 (Fredholm alternative for the eigenvalue problem). Let ϕ ∈ H 1 (Ω; R N ) ∩ L ∞ (Ω; R N ) be arbitrary and suppose that for i ∈ N, the eigenvalue λ ϕ i is simple. We further fix a L 2 ϕ (Ω; R d )-normalized eigenfunction w ϕ i to the eigenvalue λ ϕ i . Then there exists a solution u ∈ H 1 D (Ω; R d ) of the equation is a continuous, coercive bilinear form, we are able to define a continuous and compact operator that maps any right-hand side g ∈ H −1 (Ω; R d ) onto its unique solution u g ∈ H 1 D (Ω; R d ) of In the following we understand H −1 (Ω; R d ) as a Hilbert space endowed with the scalar product Indeed, this defines a scalar product as L −1 is injective. Note that due to (5.17) and the the fact that E (·) , E (·) C(ϕ) is a scalar product on H 1 D (Ω; R d ), the norm induced by (·, ·) L −1 is equivalent to the canonical operator norm on H −1 (Ω; R d ), which guarantees completeness of this space with respect to this new scalar product.
In the following, we write R(·) and N (·) denote the range and the null space of a linear operator, respectively. It is easy to see that L −1 is self-adjoint with respect to this scalar product. Furthermore the following equivalences are follow by a straightforward computation: Since L −1 is compact, we have that is a Fredholm operator. In particular, we thus know that is closed and Since L −1 is self-adjoint, we infer that where the last equivalence holds since λ ϕ i was assumed to be simple and therefore, the corresponding eigenspace is one-dimensional. In view of (5.19), this means that We further know that Hence, since L −1 f is a solution of (5.17), we have This proves the first assertion.
Let us now assume that f , w ϕ i H −1 ,H 1 = 0 and let denote the orthogonal projection onto the linear subspace w ϕ i span with respect to the scalar product on L 2 ϕ (Ω; R d ). For any solution u ∈ H 1 D (Ω; R d ) of (5.16) we obtain from the decomposition u = u − P ϕ i (u) + P ϕ i (u) that Hence, it also holds that fulfills equality (5.16). Uniqueness of the solution u ⊥ follows from the simplicity of λ ϕ i and the linearity of equation (5.16). In particular, any solution u ∈ H 1 D (Ω; R d ) can be expressed as for some α ∈ R. This completes the proof.
We can now use the Fredholm alternative to prove Theorem 5.3.
Proof of Theorem 5.3. As mentioned above, we intend to apply the implicit function theorem to prove the assertion. To this end, we define the operator Here we canonically understand the first component of the right-hand side as an element of H −1 (Ω; R d ), i.e., To apply the implicit function theorem, we need to show that F is of class C 1 on a suitable neighbourhood of ϕ, w ϕ i , λ ϕ i . For this purpose, we show that all partial Fréchet derivatives are continuous at any point in the domain of definition of F . Formally computing the partial derivatives at a point where the first two identities are to be understood in a weak sense. We can rigorously prove that the above expressions are actually the partial Fréchet derivatives. Here, we will present a detailed proof only for (5.22a) as all other derivatives can be verified analogously. We first notice that for any fixed (ϑ, w, λ) ∈ H 1 (Ω; Indeed, the linearity of the above mapping is clear and the assumptions on C, ρ along with Hölder's inequality imply the existence of a constant C > 0 such that . It further holds that Hence, Hölder's inequality yields which proves that ∂ ϑ F 1 (ϑ, w, λ) is indeed the partial derivative of F 1 with respect to ϑ in the Fréchet sense.
We next need to show that the partial derivative As now all requirements are verified, the implicit function theorem can be applied to the equation F ϕ, w ϕ i , λ ϕ i = 0. It implies that there exists a radius r 0 > 0 such that the mapping , is well-defined, continuously Fréchet differentiable and satisfies F φ, S ϕ i (φ) = 0, for allφ ∈ B r 0 (ϕ). This means that is satisfied for all ϑ ∈ (G m ∩ U c ) and all j ∈ {1, . . . , l}.
Proof. Since G m ∩ U c is convex, it holds that ϕ + t(ϑ − ϕ) ∈ G m ∩ U c for all ϑ ∈ G m ∩ U c and all t ∈ [0, 1]. As the objective functional J ε l is Fréchet differentiable by chain rule, we know that Using (5.12) in Theorem 5.3 it is now straightforward to check that J ′ ϕ (ϑ−ϕ) is identical with the right-hand side of the variational inequality. This completes the proof.

Remark 6.3.
If the first eigenvalue λ ϕ 1 is not simple but only λ ϕ 1 and further simple eigenvalues appear in the objective functional, we can still derive a variational inequality by means of the semi-differentiability established in Theorem 5.2. This is because in the above proof only variations ϕ + t(ϑ − ϕ) with positive t are considered.
If now λ ϕ 1 appears in (P ε l ) but none of the eigenvalues λ ϕ 2 = ... = λ ϕ M does, the term in the variational inequality has to be replaced by Of course, if λ ϕ 1 is simple (i.e., M = 1) both terms coincide. In the following we will only discuss the case of simple eigenvalues, but keep the fact in mind that it is not necessary to require simplicity of the first eigenvalue.

Combination of compliance and eigenvalue optimization
We now want to analyze the optimization problem (K ε l ) (that was introduced in Subsection 2.7.) by establishing results similar to those in Section 6. To this end, we will use the control-to-state operator that was introduced in [7] and maps any ϕ ∈ H 1 (Ω; R N )∩L ∞ (Ω; R N ) onto its corresponding solution u = u(ϕ) of the state equation (2.16). This allows us to consider the reduced optimization problem (with α, β ≥ 0, γ, ε > 0 and m ∈ (0, 1) N ∩ Σ N ), which is obviously equivalent to the original problem (K ε l ). The following theorem ensures the existence of a minimizer to (K ε * l ) or (K ε l ), respectively.