François Durand
email: francois.durand@inria.fr.

Benoît Kloeckner
email: benoit.kloeckner@ujf-grenoble.fr.

Fabien Mathieu
email: fabien.mathieu@alcatel-lucent.com

Ludovic Noirie
email: ludovic.noirie@alcatel-lucent.com.

GEOMETRY ON THE UTILITY SPACE *

We examine the geometrical properties of the space of expected utilities over a finite set of options, which is commonly used to model the preferences of an agent. We focus on the case where options are assumed to be symmetrical a priori. Specifically, we prove that the only Riemannian metric that respects the geometrical properties and the natural symmetries of the utility space is the round metric.

I Introduction

I.A Motivation

In order to represent an agent's preferences over options in a context of uncertainty, a simple and elegant model was formally defined by Von Neumann et Morgenstern: expected utilities. For any agent, a numerical utility is associated to each option; the utility of a lottery over the options is computed as an expected value. See Von Neumann and [START_REF] Neumann | Theory of games and economic behavior[END_REF]; [START_REF] Fishburn | Utility theory for decision making[END_REF]; [START_REF] Kreps | A Course in Microeconomic Theory[END_REF]; [START_REF] Mas-Colell | Microeconomic theory[END_REF].

For a great deal of applications, options are financial rewards or quantities of one or several goods. As a consequence, there is a natural structure over the space of options: the structure of the real line or R m , at least as a topological, or even metrical, vector space. In such cases, it is natural to study, for example, possible derivatives of the utility with respect to the reward: monotony, concavity or convexity (risk attitude), etc.

However, for some other applications, such as voting systems, options can be, for example, candidates for a position or possible decisions about a given 1 question. Without additional information, we may consider that options are symmetrical a priori. In that case, to the best of our knowledge, little attention has been paid to the geometrical properties of the utility space.

Under the assumptions of Von Neumann-Morgenstern theorem (that will be recalled), an agent's utility vector is defined up to two constants; choosing a specific normalization is arbitrary. In particular, utilities of two distinct agents are essentially incomparable without an additional normalization assumption. Geometrically, the utility space has the properties of a projective space.

A priori, this space has no natural metric. It is a pity, because it is generally convenient to work in a space with as much structure as possible. More specifically, with a metric, it is quite easy to endow the space with various probability laws, in order to model the "culture" of a population. In this paper, our goal is to give a better geometrical understanding of the utility space and to show the limited ways of endowing it with a natural metric.

I.B Contributions

Firstly, we investigate various geometrical properties of the utility space:

• We show that the utility space has a projective structure and may be seen as a quotient of the dual of the space of pairs of lotteries over the candidates;

• We naturally define an inversion operation, that corresponds to reversing preferences while keeping their intensities, and a summation operation, that is proved to preserve unanimous preferences;

• We remark that the utility space has a pathological topology when keeping the indifference point but a spherical topology when removing it.

Secondly, since the utility space is a variety, it is a natural wish to endow it with a Riemannian metric. We show that the only Riemannian representation that preserves the natural projective properties and the a priori symmetry between the candidates is the round metric.

I.C Plan

The rest of the paper is organized as follows. Von Neumann-Morgenstern theorem is recalled in section II. Several properties of this model are studied: duality with the space of pairs of lotteries in section III, inversion and summation operators in section IV, topology in section V. In section VI, it is shown what Riemannian metrics respect the natural symmetries of this space. Finally, in section VII, formulas are given that allow to use this spherical representation for applications.

In the appendix, the reader will find tables recalling the main notations of this paper.

II Von Neumann-Morgenstern model

In this section, we define some notations and we recall Von Neumann and Morgenstern theorem, which allows to represent an agent's preferences over probabilistic options by the means of an expected utility form; see Von Neumann and [START_REF] Neumann | Theory of games and economic behavior[END_REF].

We first formalize the available options and lotteries over them.

Notations II.1 (candidates and lotteries). For any k ∈ N, we denote 1, k the integer interval from 1 to k.

Let m ∈ N \ {0} and C m = 1, m . Elements in C m represent mutually exclusive options called candidates. Let L m = {(L 1 , . . . , L m ) ∈ R + m s.t. m j=1 L j = 1}. It is the set of lotteries over the candidates in C m .
For any two lotteries L = (L 1 , . . . , L m ) and M = (M 1 , . . . , M m ), for any α ∈ [0, 1], we naturally define their reduced compound lottery as a barycenter:

αL + (1 -α)M = αL 1 + (1 -α)M 1 , . . . , αL m + (1 -α)M m .
The preferences of an agent over the space of lotteries are represented by a binary relation over L m . Let ≺ be the strict relation associated to , defined as: L ≺ M ⇔ L M and not M L.

We also need to introduce a few additional notations.

Notations II.2 (linear algebra). We denote -→ 1 the vector whose m coordinates are 1 and J the m × m matrix whose m 2 coefficients are 1.

When E is a part of a vector space, we denote vect(E) the linear span of E.

The canonical inner product of -→ u and

- → v is denoted - → u | - → v . The canonical Euclidean norm of - → u is denoted - → u .
Definition II.3 (completeness). We say that relation is complete if and only

if ∀(L, M ) ∈ L m 2 : L M or M L.
Definition II.4 (transitivity). We say that relation is transitive if and only if ∀(L, M, N ) ∈ L m 3 :

L M and M N ⇒ L N.

Definition II.5 (archimedeanness). We say that relation is archimedean if and only if ∀(L, M, N ) ∈ L m 3 :

L ≺ M and M ≺ N ⇒ ∃ε ∈ ]0, 1[ s.t. (1 -ε)L + εN ≺ M ≺ εL + (1 -ε)N.
Definition II.6 (independence of irrelevant alternatives). We say that is independent of irrelevant alternatives if and only if

∀(L, M, N ) ∈ L m 3 , ∀α ∈ ]0, 1]: L ≺ M ⇒ αL + (1 -α)N ≺ αM + (1 -α)N.
Definition II.7 (utility vector). Let be a binary relation over

L m . Let - → u = (u 1 , . . . , u m ) ∈ R m .
We say that -→ u is a utility vector representing if and only, for any two lotteries L = (L 1 , . . . , L m ) and M = (M 1 , . . . , M m ):

L M ⇔ m j=1 L j u j ≤ m j=1 M j u j .
Theorem II.8 (Von Neumann and Morgenstern). Let be a binary relation over L m .

The following conditions are equivalent.

1. The binary relation is complete, transitive, archimedean and independent of irrelevant alternatives.

2. There exists -→ u ∈ R m a utility vector representing .

When they are met, -→ u is defined up to an additive constant and a positive

multiplicative constant: if - → u ∈ R m is a utility vector representing , then - → v ∈ R m is a utility vector representing if and only if ∃a ∈ ]0, +∞[, ∃b ∈ R s.t. - → v = a - → u + b - → 1 .
For a proof, see Von [START_REF] Von Neumann | Theory of Games and Economic Behavior[END_REF], [START_REF] Mas-Colell | Microeconomic theory[END_REF] or [START_REF] Kreps | A Course in Microeconomic Theory[END_REF].

Discussing whether completeness, transitivity and archimedeanness assumptions are experimentally valid is out of the scope of this paper: see [START_REF] Fishburn | Nonlinear preference and utility theory[END_REF]; [START_REF] Mas-Colell | Microeconomic theory[END_REF]. Let us note that, while keeping a sheer simplicity, this model is sophisticated enough to deal with relative intensities of the agent's preferences. This gives more information than just a total preorder over the candidates.

Definition II.9 (utility space). Let ≈ be the equivalence relation defined by

∀( - → u , - → v ) ∈ (R m ) 2 : - → u ≈ - → v ⇔ ∃a ∈ ]0, +∞[, ∃b ∈ R s.t. - → v = a - → u + b - → 1 .
We call utility space over m candidates, and we denote U m , the quotient space R m / ≈.

We call canonical projection from R m to U m the function:

≈ π : R m → U m - → u → { - → v ∈ R m s.t. - → v ≈ - → u }. For any - → u ∈ R m and ≈ u = ≈ π( - → u ),

III Duality

In this section, we remark that the utility space is a dual of the space of pairs of lotteries. This will be especially helpful to demonstrate theorem IV.9, characterizing the summation operator that will be defined in section IV.

In the example represented in figure I, we consider m = 3 candidates and -→ u = ( 53 ; -1 3 ; -4 3 ). The great triangle is the space of lotteries L m , that is, the simplex defined by L i = 1 and ∀i, L i ≥ 0. Hatchings are the agent's indifference lines: she is indifferent between any pair of lotteries on the same line (see [START_REF] Mas-Colell | Microeconomic theory[END_REF], section 6.B). Here, we have drawn a utility vector -→ u that is in the plane of the great triangle and orthogonal to the indifference lines. This is not mandatory, since -→ u can be arbitrarily chosen in ≈ u, but it is a quite natural choice, since the component of -→ u in the direction -→ 1 , orthogonal to the simplex, has no meaning in terms of preferences. The utility vector -→ u may be seen as a (uniform) gradient of preference: at each point, it reveals in what directions one can find lotteries that are preferred by the agent. However, only the direction of -→ u is important, whereas its norm has no specific meaning; as a consequence, the utility space is not exactly a dual space, but rather a quotient of a dual space, as we will see more formally.

L 1 L 2 L 3 1 1 1 L m - → u
Definition III.1 (bipoint). For any two lotteries

L = (L 1 , . . . , L m ) and M = (M 1 , . . . , M m ), we call bipoint from L to M the vector in R m : --→ LM = (M 1 -L 1 , . . . , M m -L m ).
Definition III.2 (tangent polytope and hyperplane). We call tangent polytope of L m the set T of bipoints of L m :

T = { --→ LM , (L, M ) ∈ L m 2 }.
We call tangent hyperplane of L m :

H = {(∆ 1 , . . . , ∆ m ) ∈ R m s.t. m j=1 ∆ j = 0}.
Remark III.3 (geometrical interpretation). In figure I, T is the set of the bipoints of the great triangle, seen as a part of a vector space (whereas L m is a part of an affine space). Hyperplane H is the whole vector plan containing T .

Lemma III.4 (T generates H). ∀ - → ∆ ∈ H, ∃ --→ LM ∈ T , ∃λ ∈ ]0, +∞[ s.t. - → ∆ = λ --→ LM .
Proof. Polytope T contains a neighborhood of the origin in vector space H. Definition III.5 (linear forms on H). Let H be the dual space of H, that is, the set of linear forms on H.

Let - → u = (u 1 , . . . , u m ) ∈ R m .
We call linear form associated to -→ u the following element of H :

- → u | : H → R - → ∆ = (∆ 1 , . . . , ∆ m ) → - → u | - → ∆ = m j=1 u j ∆ j . We call positive half-hyperplane associated to - → u : - → u +H = { - → ∆ ∈ H s.t. - → u | - → ∆ ≥ 0}.
Proposition III.6 (preferences expressed in terms of linear forms). Let (L, M ) ∈ L m 2 and -→ u ∈ R m . Then:

L ≈ u M ⇔ - → u | --→ LM ≥ 0 ⇔ --→ LM ∈ - → u +H .
Proof. This follows easily from the definitions.

Proposition III.7 (utility space and linear forms). For any (f, g) ∈ (H ) 2 , we denote f ≈ g if and only ∃a ∈ ]0, +∞[ s.t. g = af . We denote

≈ π(f ) = {g ∈ H s.t. f ≈ g}.
For any ( -→ u , -→ v ) ∈ (R m ) 2 , we have:

- → u ≈ - → v ⇔ - → u | ≈ - → v |.
The following application is a bijection:

Θ : U m → H / ≈ ≈ π( - → u ) → ≈ π( - → u |). Proof. - → u ≈ - → v ⇔ ∃a ∈ ]0, +∞[, ∃b ∈ R s.t. - → v -a - → u = b - → 1 ⇔ ∃a ∈ ]0, +∞[ s.t. - → v -a - → u is orthogonal to H ⇔ ∃a ∈ ]0, +∞[ s.t. ∀ - → ∆ ∈ H, - → v | - → ∆ = a - → u | - → ∆ ⇔ - → u | ≈ - → v |.
The implication ⇒ proves that Θ is correctly defined: indeed, if

≈ π( - → u ) = ≈ π( - → v ), then ≈ π( - → u |) = ≈ π( - → v |). The implication ⇐ ensures that Θ is injective.
Finally, Θ is obviously surjective.

Remark III.8 (interpretation). Hence, the utility space may be seen as a quotient of the dual H of the tangent space H of the lotteries L m . A utility vector may be seen, up to a positive constant, as a linear form over H that reveals, for any point in the space of lotteries, in what directions the agent can find lotteries that she prefers.

IV Inversion and summation operators

As a quotient of R m , the utility space inherits natural operations from R m : inversion and summation. We will see that both these quotient operators have an intuitive meaning regarding preferences.

Definition IV.1 (inversion). In U m , the inversion operator is defined as:

-: U m → U m ≈ π( - → u ) → ≈ π(-- → u ).
Remark IV.2 (the inversion operator is a bijection). The inversion operator is correctly defined and it is a bijection. Indeed,

≈ π( - → u ) = ≈ π( - → v ) if and only if ≈ π(-- → u ) = ≈ π(-- → v ).
Remark IV.3 (interpretation in terms of preferences). Considering the additive inverse amounts to reverting the agent's preferences, without modifying their relative intensities. Now, we want to push the summation operator from R m to the quotient U m . The idea is quite simple and general: considering Converting back to the quotient, we can get for example

≈ π(2 - → u + - → v ) and ≈ π( - → u + 3 - → v
), which are generally distinct. As a consequence, the set of outputs is not the utility space U m , but rather the set P(U m ) of its subsets. This example shows how we could define the sum of two elements ≈ u and ≈ v. In order to be more general, we will define the sum of any number of elements of U m . Hence we also take P(U m ) as the set of inputs.

Definition IV.4 (summation). On the subsets of U m , we define the summation operator as:

: P(U m ) → P(U m ) A → ≈ π ( n i=1 - → u i ) , n ∈ N, ( -→ u 1 , . . . , -→ u n ) ∈ ≈ π -1 (A) n .
Example IV.5 (sum of two utility vectors

- → u - → v ≈ u, ≈ v} -→ h 1 -→ h 2 -→ h 3
Figure II: Sum of two utility vectors in the utility space.

Remark IV.6 (geometrical interpretation). Geometrically, the sum is the quotient of the convex hull of the inputs. Note that this convex hull is actually a convex cone. A bit later, we will see its interpretation in terms of preferences.

Remark IV.7 (closed cone). Our convention is to consider the closed cone: for example, ≈ u fits in our definition. That would not be the case if we took only

≈ π(a - → u + a - → v + b)
, where a > 0. In the same spirit, ≈ π( -→ 0 ) fits in our definition. Generally, that would not be the case if we took only

≈ π(a - → u + a - → v + b)
, where a > 0 and a > 0. The stated purpose of this convention is to have a concise wording for theorem IV.9 that follows.

Proposition IV.8 (basic properties of the sum). Let A ∈ P(U m ). Then:

1. ≈ π( - → 0 ) ∈ A; 2. A ⊂ A.
Proof.

1. Consider n = 0 in the definition.

For

≈ u ∈ A, if it sufficient to consider n = 1 and -→ u 1 ∈ ≈ π -1 ( ≈ u) in order to prove that ≈ u ∈ A.
We now prove that, if A is the set of the utility vectors of a population, then A is the set of utility vectors that respect the unanimous preferences of A.

Theorem IV.9 (characterization of the sum). Let A ∈ P(U m ) and

≈ v ∈ U m . The following conditions are equivalent. 1. ≈ v ∈ A. 2. ∀(L, M ) ∈ L m 2 : ∀ ≈ u ∈ A, L ≈ u M ⇒ L ≈ v M. Proof. Let - → v ∈ ≈ π -1 ( ≈ v).
We have the following equivalences.

• ∀(L, M ) ∈ L m 2 , ∀ ≈ u ∈ A, L ≈ u M ⇒ L ≈ v M , • ∀ --→ LM ∈ T , ∀ - → u ∈ ≈ π -1 (A), - → u | --→ LM ≥ 0 ⇒ - → v | --→ LM ≥ 0, • ∀ - → ∆ ∈ H, ∀ - → u ∈ ≈ π -1 (A), - → u | - → ∆ ≥ 0 ⇒ - → v | - → ∆ ≥ 0 (cf. lemma III.4), • - → u ∈ ≈ π -1 (A) - → u +H ⊂ - → v +H , • - → v is in the convex cone of ≈ π -1 (A), • ≈ v ∈ A.
Example IV.10 (a utility vector and its inverse). Let us consider a non null However, we have just seen that the sum is the subset of utility vectors preserving the unanimous preferences over lotteries. Intuitively, we could think that, since ≈ u and -≈ u seem to always disagree, any utility vector ≈ v respects the empty set of their common preferences; so, their sum should be the whole space.

But this is a false assumption. For example, let us consider -→ u = (1, 0, . . . , 0).

About any two lotteries L and M , inverse opinions ≈ u and -≈ u agree if and only if L 1 = M 1 : in that case, both of them are indifferent between L and M . Hence, ≈ u and -≈ u do have a characteristic common feature: they pay attention to candidate 1 only, even if their opinions about her are diverging.

V Topology

In this section, we remark that the utility space has a pathological topology when keeping the indifference point, but a spherical topology when removing it. This will be used in section VI to demonstrate theorem VI.5, characterizing the suitable Riemannian metrics. Definition V.1 (T 0 space). We say that a topological space is a T 0 space or a Kolmogorov space if and only if, for any two distinct points, at least one of them has an open neighborhood not containing the other. See [START_REF] Guénard | Compléments d'analyse[END_REF]. Definition V.2 (T 1 space). We say that a topological space is a T 1 space if and only if, for any two distinct points, each of them has an open neighborhood not containing the other. See [START_REF] Guénard | Compléments d'analyse[END_REF].

Proposition V.3 (topology of the utility space). We assume m ≥ 2. The utility space U m = R m / ≈ is endowed with the quotient topology.

U m \{ ≈ 0} has the same topology as the sphere of dimension m-2. In particular it is a T 1 space.

U m is a T 0 space but not a T 1 space.

Proof. Let S m-2 = { - → u = (u 1 , . . . , u m ) ∈ R m s.t. m j=1 u j = 0, m j=1 u 2 j = 1}. Then: Θ : S m-2 → U m \ { ≈ 0} - → u → ≈ π( - → u )

VI Riemannian representation

Since the utility space is a variety, it is a natural desire to endow it with a Riemannian metric. In this section, we prove that there is a limited choice of metrics that are coherent with the natural properties of the space and with the a priori symmetry between the candidates. Definition VI.1 (line). We assume that m ≥ 3. The image of a 2-dimensional subspace of

H in U m \ { ≈ 0} by the canonical projection ≈ π is called a line in U m \ { ≈ 0}.
Remark VI.2 (connection with the sum operator). This notion is deeply connected to the summation operator defined in section IV: indeed, the sum of two non antipodal points in U m \ { ≈ 0} is a segment of the line joining them, the shortest one in the sense of the round metric. We will show that this property holds true only for this metric. Remark VI.4 (geometrical interpretation). In order to get an intuitive vision of this metric, one can represent any position ≈ u by a vector -→ u that verify:

u i = 0, u 2 i = 1.
We obtain a (m -2)-dimensional sphere and we consider the metric induced by the canonical Euclidean metric of R m . That is, distances are measured on the surface of the sphere, using the restriction of the canonical inner product on each tangent space.

Let us remark that any inner product on R m induces a Riemannian metric on U m \{ ≈ 0} by the same process and that any two of these metrics are isometric; however, they are not equal, except if these inner products are multiple of each other. The difference between isometry and equality is essentially the same as between two loops of equal length: there exists an application from one to the other that preserves distances along the curve, but they are not generally equal.

We will prove that for any m ≥ 4, the spherical representation is the only one that is coherent with the natural properties of the space and that respects the a priori symmetry between candidates. Proof. Since the implication 2 ⇒ 1 is obvious, we now prove 1 ⇒ 2. The deep result behind this is a classical theorem of Beltrami, which dates back to the middle of the nineteenth century: see Beltrami (1866) and Beltrami (1869).

Beltrami's theorem precisely tells us that condition 1a implies that U m \ { ≈ 0} has constant curvature. Note that this result is in fact more subtle in dimension 2 (that is, for m = 4) than in higher dimensions; see Spivak (1979a), Theorem 1.18 andSpivak (1979b), Theorem 7.2 for proofs.

Remark VII.4 (utilities have no absolute meaning). Hence, for any point ≈ u in the utility space, we use a canonical representative belonging to H and whose Euclidean norm is equal to 1. For example, we might have the following representatives:

-→ u = (0.71, -0.71, 0) , -→ v = (0.82, -0.41, -0.41) .

Since v 1 > u 1 , does it mean that an agent with preferences ≈ v "appreciates" candidate 1 more than an agent with preferences ≈ u? We don't say anything like this: it is well known that utilities belonging to two agents are essentially incomparable: see Von Neumann and [START_REF] Neumann | Theory of games and economic behavior[END_REF]. Taking canonical representatives on the (m -2)-dimensional sphere is only used to compute distances between two points in the utility space.

Regarding the other question considered in this section, how to draw vectors according to a uniform probability law over U m \ { ≈ 0}, it is sufficient to use a uniform law on the unit sphere in H.

To conclude about this spherical representation, here is an example of a law defined by its density. Given a vector -→ u 0 in the unit sphere of H and κ a nonnegative real number, the Von Mises-Fisher distribution of pole -→ u 0 and concentration κ is defined by the following density with respect to the uniform law on the unit sphere in H:

p( - → u ) = C κ e κ - → u | -→ u0 ,
where C κ is a normalization constant.

It can be proved that, given the mean resultant vector of a distribution over the sphere, Von Mises-Fisher distribution maximizes the entropy, in the same way that, in the Euclidean space, Gaussian distribution maximizes the entropy among laws with given mean and standard deviation. Hence, without additional information, it is the "natural" distribution that should be used. In order to draw data according to this law for simulation purposes, one may use Ulrich algorithm modified by Wood: see [START_REF] Ulrich | Computer generation of distributions on the m-sphere[END_REF] and [START_REF] Wood | Simulation of the von Mises Fisher distribution[END_REF].

VIII Conclusion

We have studied the geometrical properties of the classical model of expected utilities, introduced by Von Neumann and Morgenstern, when candidates are considered symmetrical a priori. We have remarked that the utility space may be seen as a dual of the space of lotteries, that inversion and summation operators inherited from R m have a natural interpretation in terms of preferences and that the space has a spherical topology when the indifference point is removed.

We have proved that the only Riemannian representation that respects the projective lines naturally defined by the summation operator and the symmetry between candidates is a round sphere.

All these considerations lay on the principle to add as little information as possible in the system, especially by respecting the a priori symmetry between candidates. This does not imply that the spherical representation of the utility space U m is the most relevant one in order to study a specific situation. Indeed, as soon as one has additional information (for example, a model that places candidates in a political spectrum), it is natural to include it in the model. However, if one wishes, for example, to study a voting system in all generality, without focusing on its application in a specific field, it looks natural to consider a utility space with a metric as neutral as possible, like the one defined in this paper by the spherical representation.

  we denote without ambiguity ≈ u the binary relation over L m represented by -→ u in the sense of theorem II.8 (Von Neumann-Morgenstern).

Figure I :

 I Figure I: Space of lotteries for 3 candidates.

v

  in U m , their antecedents are taken in R m thanks to ≈ π -1 , the sum is computed in R m , then the result is converted back into the quotient space U m , thanks to ≈ π. However, the result is not unique. Indeed, let us take arbitrary representatives -In order to compute the sum, we can think of any representatives. So, possible sums are a -→ u + a -→ v + b + b , where a and a are positive and where b + b is any real number.

≈u

  and let us examine the sum of ≈ u and its additive inverse -≈ u. Applying the definition, we see that the sum consists of ≈

  NotationsVI.3 (round metric). The quotient R m /vect( -→ 1 ) is identified to H and endowed with the inner product inherited from the canonical one of R m . The utility space U m \ { ≈ 0} is identified to the unit sphere of H and endowed with the induced Riemannian structure. Let us call ξ 0 this Riemannian metric on U m \ { ≈ 0}.

  Theorem VI.5 (Riemannian representation of the utility space). We assume that m ≥ 4. Let ξ be a Riemannian metric on U m \ { ≈ 0}. Conditions 1 and 2 are equivalent. 1. (a) The lines of U m \ { ≈ 0} are geodesics of ξ; and (b) for any permutation σ of 1, m , the action Φ σ induced on U m \ { ≈ 0} by (u 1 , . . . , u m ) → (u σ(1) , . . . , u σ(m) ) is an isometry. 2. ∃λ ∈ ]0, +∞[ s.t. ξ = λξ 0 .

  ). Let us consider U 4 . In figure II, for the purpose of visualization, we represent its projection in H, which is permitted by the choice of normalization constants b.

	let (	-→ h 1 ,	-→ h 2 ,	-→ h 3 ) be an orthonormal base of H.	Since H is a 3-dimensional space,
	For two non-trivial utility vectors tiplicators a allow to choose representatives -≈ u and ≈ v, the choice of normalization mul-→ u and -→ v whose Euclidean norm
	is 1.
	In this representation, the sum { u, ≈ to vectors a -→ u + a -→ v , where a and a are nonnegative. Indeed, we took a ≈ v} consists of utilities corresponding
	representation in H, so all normalization constants b vanish. Moreover, a, a or both can be equal to zero because our definition allows to ignore -→ u , -→ v or both.
	Up to taking representatives of unitary norm for non-trivial utility vectors, let
	us note that the sum { u, ∼ point -→ 0 of total indifference. ∼ v} can be represented by the dotted line and the

Since U m \ { ≈ 0} is a topological sphere, this constant curvature must be positive. Up to multiplying ξ by a constant, we can assume that this constant curvature is 1. As a consequence, there is an isometry Ψ : S m-2 → U m \ { ≈ 0}, where S m-2 is the unit sphere of R m-1 , endowed with its usual round metric. The function Ψ obviously maps geodesic to geodesics, and by a standard projective geometry argument there is a linear map Λ : R m-1 → H inducing Ψ, that is such that:

where

Using Ψ to push the canonical inner product of R m-1 , we get that there exists an inner product

with the unit sphere defined in H by φ and restricting φ to it.

The last thing to prove is that φ is the inner product coming from the canonical one on R m . Note that hypothesis 1b is mandatory, since any inner product on H does induce on U m \ { ≈ 0} a Riemannian metric satisfying 1a. Each canonical basis vector -→ e j = (0, . . . , 1, . . . , 0) defines a point in U m \ { ≈ 0} and a half-line j in H. Condition 1b ensures that these half-lines are permuted by isometries of (H, φ). In particular, there are vectors -→ u j ∈ j that have constant pairwise distance (for φ), so that they form a regular simplex and j -→ u j = -→ 0 . By a standard projective geometry argument, -→ u 1 , . . . , ---→ u m-1 is up to multiplication by a scalar the unique basis of H such that -→ u j ∈ j and j<m -→ u j ∈m . Now consider the canonical inner product φ 0 on H that comes from the canonical one on R m . Since permutations of coordinates are isometries, we get that the vectors -→ v j = Π( -→ e j ) (where Π is now the orthogonal projection from R m to H) form a regular simplex for φ 0 , so that j -→ v j = -→ 0 . It follows that -→ u j = λ -→ v j for some λ > 0 and all j. We deduce that the -→ u j form a regular simplex for both φ and φ 0 , which must therefore be multiple from each other.

Remark VI.6 (utility space for three candidates). However, the implication 1 ⇒ 2 in the theorem is not true for m = 3. For each non null utility vector, let us consider its representative verifying min(u i ) = 0 and max(u

is endowed with the metric induced on these edges by the canonical inner product on R 3 . Then conditions 1a and 1b of the theorem are met, but not condition 2.

VII Riemannian representation in practice

The Riemann representation introduced in the previous section may be used for two main purposes.

• One may want to compute distances between two utility vectors.

• Once the space is endowed with a metric, it is endowed with a natural probability measure: the uniform measure in the sense of this metric. Hence, one can use the uniform measure itself or define other measures by their densities with respect to the uniform measure.

In this section, we detail how to deal with these two aspects in practice.

Remark VII.1 (indifference point). Since the point ≈ 0 is a geometrical singularity, it is difficult to include it naturally in such a measure. If one wants to take it into account, the easiest way is to draw it with a given probability and to use some given measure over U m \ { ≈ 0} in the other cases.

Notations VII.2 (orthogonal projection on H). We denote D 0 the matrix of the orthogonal projection on H:

Then, in the sense of metric λξ 0 of theorem VI.5, the distance between

Proof. First of all, we compute the projection

Then we normalize -→ u in order to get it on the unit sphere of H, that is,

Now, we just have to compute the distance between -→ u and -→ v on the surface of the unit sphere and to multiply by λ: we get d(

Vector whose m coordinates are 1. vect(E) Linear span of E, where E is a part of a vector space.


Axel Kröner
email: axel.kroener@oeaw.ac.at

Karl Kunisch
email: karl.kunisch@uni-graz.at

A minimum effort optimal control problem for the wave equation

Keywords: Optimal control, wave equation, semi-smooth Newton methods, finite elements

A minimum effort optimal control problem for the undamped wave equation is considered which involves L ∞ -control costs. Since the problem is non-differentiable a regularized problem is introduced. Uniqueness of the solution of the regularized problem is proven and the convergence of the regularized solutions is analyzed. Further, a semi-smooth Newton method is formulated to solve the regularized problems and its superlinear convergence is shown. Thereby special attention has to be paid to the well-posedness of the Newton iteration. Numerical examples confirm the theoretical results.

Q = J × Ω, and Σ = J × ∂Ω. We consider the following problem

             min (y,u)∈X×U 1 2 C ωo y -z 2 L 2 (Q) + α 2 u 2 L ∞ (Q) , s.t.
Ay = B ωc u in Q, y(0) = y 0 in Ω,

Cy = 0 on Σ (P 1 )
for given state space X, control space U , initial point y 0 , parameter α > 0 and desired state z ∈ L 2 (Q). A denotes the wave operator, B ωc the control operator, C ωo : X → L 2 (Q) an observation operator, where ω o , ω c ⊂ Ω describe the area of observation and control, and C : X → L 2 (Q) denotes a boundary operator. A detailed formulation in a functional analytic setting is given in the following section.

The interpretation of the cost functional in (P 1 ) can be described as minimizing the tracking error by means of a control which is pointwise as small as possible. The appearance of the L ∞ -control costs leads to nondifferentiability. The analytic and efficient numerical treatment of this nonsmooth problem by a semi-smooth Newton method stands in the focus of this work. We prove superlinear convergence of this iterative method and present numerical examples.

Numerical methods for minimum effort problems in the context of ordinary differential equation are developed in publications, see, e.g., Neustadt [START_REF] Neustadt | Minimum effort control systems[END_REF], and Ito and Kunisch [START_REF] Ito | Minimal effort problems and their treatment by semi-smooth Newton methods[END_REF] and the references given there. In the context of partial differential equations there exist only few results, see the publications on elliptic equations by Grund and Rösch [START_REF] Grund | Optimal control of a linear elliptic equation with a supremum norm functional[END_REF] and Clason, Ito, and Kunisch [START_REF] Clason | A minimum effort optimal control problem for elliptic PDEs[END_REF].

The literature for numerical methods for optimal control of the wave equation is significantly less rich than for that for equations of parabolic type. Let us mention some selected contributions for the wave equation. In [START_REF] Gugat | Penalty techniques for state constrained optimal control problems with the wave equation[END_REF] Gugat treats state constrained optimal control problems by penalty techniques. Gugat and Leugering [START_REF] Gugat | L ∞ -norm minimal control of the wave equation: on the weakness of the bang-bang principle[END_REF] analyze bang-bang properties for L ∞ -norm minimal control problems for exact and approximate controllability problems and give numerical results. Time optimal control problems are considered by Kunisch and Wachsmuth [START_REF] Kunisch | On time optimal control of the wave equation and its numerical realization as parametric optimization problem[END_REF] and semi-smooth Newton methods for control constrained optimal control problems with L 2 -control costs in Kröner, Kunisch, and Vexler [START_REF] Kröner | Semi-smooth Newton methods for optimal control of the wave equation with control constraints[END_REF]. Gerdts, Greif, and Pesch in [START_REF] Gerdts | Numerical optimal control of the wave equation: Optimal boundary control of a string to rest in finite time[END_REF] present numerical results driving a string to rest and give further relevant references. A detailed analysis of discretization issues for controllability problems related to the wave equation is contained in the work by Zuazua, see e.g. [START_REF] Zuazua | Propagation, observation, and control of waves approximated by finite difference methods[END_REF] and Ervedoza and Zuazua [START_REF] Ervedoza | On the Numerical Approximations of Exact Controls for Waves[END_REF].

We will present an equivalent formulation of the minimal effort problem (P 1 ) with a state equation having a bilinear structure and controls satisfying pointwise constraints, i.e. we move the difficulty of nondifferentiability of the control costs in the cost functional to additional constraints. To solve the problem we apply a semi-smooth Newton method.

Special attention has to be paid to the well-posedness of the iteration of the semi-smooth Newton scheme. The lack of smoothing properties of the equations (when compared to elliptic equations) requires careful choice of the function space setting to achieve superlinear convergent algorithms; in particular, the problem is not reformulated in the reduced form as in [START_REF] Clason | A minimum effort optimal control problem for elliptic PDEs[END_REF] and the well-definedness of the Newton updates has to be shown.

The restriction to dimensions d ≤ 4 is due to the Sobolev embedding theorem which is needed in Lemma 5.1 to verify Newton differentiability.

The paper is organized as follows. In Section 2 we make some preliminary remarks, in Section 3 we formulate the minimal effort problem, in Section 4 we present a regularized problem, in Section 5 we formulate the semi-smooth Newton method, in Section 6 we discretize the problem, and in Section 7 we present numerical examples.

Preliminaries

In this paper we use the usual notations for Lebesgue and Sobolev spaces. Further, we set

Y 1 = H 1 0 (Ω) × L 2 (Ω), P 1 = L 2 (Ω) × H 1 0 (Ω), Y 0 = L 2 (Ω) × H -1 (Ω), P 0 = H -1 (Ω) × L 2 (Ω).
There holds the following relation between these spaces

(Y 1 ) * = P 0 , (Y 0 ) * = P 1 ,
where * indicates the dual space. For a Banach space Y we set

L 2 (Y ) = L 2 (J, Y ), H 1 (Y ) = H 1 (J, Y ).
Further, we introduce the following operators

A : X = L 2 (Y 1 ) ∩ H 1 (Y 0 ) -→ Y = L 2 (Y 0 ) ⊕ H 1 (P 0 ) * , A * : Y * = L 2 (P 1 ) ∩ H 1 (P 0 ) -→ X * = L 2 (P 0 ) ⊕ H 1 (Y 0 ) * with A = ∂ t -id -∆ ∂ t , A * = -∂ t -∆ -id -∂ t (2.1)
for the Laplacian (-∆) : H 1 0 (Ω) → H -1 (Ω) and identity map id :

L 2 (Q) → L 2 (Q). For (y, p) = (y 1 , y 2 , p 1 , p 2 ) ∈ X × Y * there holds the relation Ay, p Y,Y * + (y(0), p(0)) L 2 (Ω) = y, A * p X,X * + (y(T ), p(T )) L 2 (Ω) , (2.2) since ∂ t y 1 -y 2 ∂ t y 2 -∆y 1 , p 1 p 2 Y,Y * + (y(0), p(0)) L 2 (Ω) = y 1 y 2 , -∂ t p 1 -∆p 2 -∂ t p 2 -p 1 X,X * + (y(T ), p(T )) L 2 (Ω) . (2.3)
We introduce the observation and control operator

C ωo = χ ωo id, 0 : L 2 → L 2 , (2.4)
B ωc = 0, χ ωc id : L 2 → L 2 (2.5)
with the characteristic functions χ ωo and χ ωc of J × ω o and J × ω c for given nonempty open subsets ω o , ω c ⊂ Ω. Here we used the notation

L 2 = L 2 (L 2 (Ω) × L 2 (Ω))
. Further, we define the boundary operators

C = C Ω , B = B Ω .
For the inner product in L 2 (Q) we write

(•, •) = (•, •) L 2 (Q) .
Throughout this paper C > 0 denotes a generic constant.

The minimum effort problem for the wave equation

In this section we present the minimum effort problem (P 1 ) in detail and formulate an equivalent problem in which we move the difficulty of the nondifferentiability of the control costs to additional control constraints. Furthermore we derive the optimality system for the latter problem.

To make the minimum effort problem (P 1 ) precise we choose U = L ∞ (Q), y 0 ∈ Y 1 and the operators A, B ωc C ωo , C and B as defined in the previous section.

Problem (P 1 ) can be formulated equivalently as

                   min (y,u,c)∈X×U ×R + 0 J(y, c) = 1 2 C ωo y -z 2 L 2 (Q) + α 2 c 2 , s.t. Ay = B ωc u in Q, y(0) = y 0 in Ω, Cy = 0 on Σ, u U ≤ c. (P 2 )
Except for the case c = 0 problem (P 2 ) is equivalent to problem (P 3 ) given by

                   min (y,u,c)∈X×U ×R + 0 J(y, c) = 1 2 C ωo y -z 2 L 2 (Q) + α 2 c 2 , s.t. Ay = cB ωc u in Q, y(0) = y 0 in Ω, Cy = 0 on Σ, u U ≤ 1. (P 3 )
By standard arguments the existence of a solution (y * , u * , c * ) ∈ X × U × R + 0 of problem (P 3 ) can be verified. The reformulation of the problem has the advantage that the domain of the control space does not depend on the parameter c. Remark 3.1 For c = 0 any control u with u U ≤ 1 is a minimizer of (P 3 ). To avoid this case we assume that

J(y * , c * ) < 1 2 z 2 L 2 (Q) (3.1) for a solution (y * , u * , c * ). If c = 0 problem (P 3 ) reduces to                    min (y,u)∈X×U J(y) = 1 2 C ωo y -z 2 L 2 (Q) , s.t. Ay = 0 in Q, y(0) = y 0 in Ω, Cy = 0 on Σ, u U ≤ 1.
Thus, y is determined by the equation and u can be chosen arbitrarily as far as the pointwise constraints are satisfied. If (3.1) holds, we have

1 2 C ωo y -z 2 L 2 (Q) + α 2 c 2 < 1 2 z 2 L 2 (Q)
and with c = 0 this leads to the contradiction

1 2 z 2 L 2 (Q) < 1 2 z 2 L 2 (Q) .
By standard techniques the optimality system can be derived.

Lemma 3.1 The optimality system for problem (P 3 ) is given by

           A * p + C * ωo C ωo y -C * ωo z = 0, p(T ) = 0, Bp| Σ = 0, (-B * ωc p, δu -u) ≥ 0 for all δu with δu L ∞ (Q) ≤ 1, αc -(u, B * ωc p) = 0, Ay -cB ωc u = 0, y(0) = y 0 , Cy| Σ = 0 (3.2) with p ∈ Y * .
From the pointwise inspection of the second relation in the optimality system (3.2) we obtain for (t, x)

∈ J × Ω u(t, x) =    1 if B * ωc p(t, x) > 0, -1 if B * ωc p(t, x) < 0, s ∈ [-1, 1] if B * ωc p(t, x) = 0 (3.3) or equivalently u ∈ Sgn(B * ωc p) with Sgn(s) =    1 if s > 0, -1 if s < 0, [-1, 1] if s = 0.
Eliminating the control we obtain the reduced system

     A * p + C * ωo C ωo y -C * ωo z = 0, p(T ) = 0, Bp| Σ = 0, αc -B * ωc p L 1 (Q) = 0, Ay -cB ωc sgn(B * ωc p) ∋ 0, y(0) = y 0 , Cy| Σ = 0. (3.4)
Under certain conditions the solution of (P 3 ) is unique.

Lemma 3.2 For c = 0 and ω c ⊂ ω o the solution of problem (P 3 ) is unique if we set the control to zero on Q \ (J × ω c ).

Remark 3.2

The value of the control on the domain Q \ (J × ω c ) has no influence on the solution of the control problem as far as u L ∞ (Q\(J×ωc)) ≤ u L ∞ (J×ωc) . To obtain uniqueness we set u ≡ 0 on Q \ (J × ω c ).

Proof of Lemma 3.2 Because of the equivalence of (P 3 ) and (P 1 ) for c = 0 it suffices to prove uniqueness for the latter one. Let S : U → L 2 (Q) be the control-to-state mapping for the state equation given in (P 1 ). Further, let (y, u) be a solution of (P 1 ) with the cost given by

F (u) := 1 2 C ωo S(u) -z 2 L 2 (Q) + α 2 u 2 L ∞ (Q) .
For ω o = Ω the map C ωo S is injective and we have strict convexity of

1 2 C ωo S(u) -z Ay 1 = B ωc u 1 in Q, y 1 (0) = y 0 in Ω, Cy 1 = 0 on Σ, Ay 2 = B ωc u 2 in Q, y 2 (0) = y 0 in Ω, Cy 2 = 0 on Σ for (y i , u i ) ∈ X × U , i = 1, 2, with C ωo y 1 = C ωo y 2 . (3.5)
This implies, that y 1 -y 2 = 0 on J × ω o . Hence, A(y 1 -y 2 ) = 0 on J × ω o and thus, u 1 = u 2 on J × ω c . Consequently, we derive y 1 = y 2 on Q. Thus C ωo S is injective and we are in the situation as above.

Remark 3.3 For general ω o ⊂ Ω and ω c ⊂ Ω we cannot expect uniqueness due to the finite speed of propagation. Consider a one dimensional domain Ω = (0, L), L > 0, with ω o = (0, ε), ε > 0 small, and ω c = Ω. Let (y, u) be a corresponding solution of (P 1 ) with y = 0 in an open subset of Q \ (J × ω o ). Then there exists an open set J in Q \ (J × ω o ) in which the adjoint state p does not vanish. Thus we have u = 0 on J. Let ŷ be the solution of

Aŷ = B ωc g in Q, ŷ(0) = 0 in Ω, ŷ = 0 on Σ with g = -sgn(u)η in B δ , 0 else
for δ, η > 0 and B δ ⊂ J and the usual sgn-function. Here, B δ denotes a ball with radius δ with respect to the topology of

Q. Then u + g L ∞ (Q) = u L ∞ (Q)
and C ωo (ŷ + y) = C ωo y for δ, η > 0 sufficiently small. Thus we obtain a second solution (u + g, ŷ).

The regularized minimum effort problem

The optimality system in (3.4) is not (in a generalized sense) differentiable. Therefore we consider a regularized minimum effort problem given by

                         min (y,u,c)∈X×U ×R + 0 J β (y, u, c) = 1 2 C ωo y -z 2 L 2 (Q) + βc 2 u 2 L 2 (Q) + α 2 c 2 , s.t. Ay = cB ωc u in Q, y(0) = y 0 in Ω, Cy = 0 on Σ, u L ∞ (Q) ≤ 1 (P reg ) for y 0 ∈ Y 1 , parameters α, β > 0, and z ∈ L 2 (Q).
The existence of a solution follows by standard arguments which we denote by (y β , u β , c β ).

Remark 4.1 The regularization term scales linearly with the parameter c. Alternative regularizations, where the penalty term is constant or quadratic in c, are discussed in [START_REF] Clason | A minimum effort optimal control problem for elliptic PDEs[END_REF]. Remark 4.2 Let (y β , u β , c β ) be a solution of (P reg ). To exclude the case c β = 0 for β sufficiently small we assume z ≡ 0 and

J(y * , c * ) < 1 2 z 2 L 2 (Q) . (4.1) If c β = 0 we have 1 2 z 2 L 2 (Q) = 1 2 C ωo y β -z 2 L 2 (Q) + βc β 2 u β 2 L 2 (Q) + α 2 c 2 β ≤ 1 2 C ωo y * -z 2 L 2 (Q) + βc * 2 u * 2 L 2 (Q) + α 2 (c * ) 2 ≤ J(y * , c * ) + β(c * ) 2 2 meas(Q) ≤ 1 2 z 2 L 2 (Q) + J(y * , c * ) -z 2 L 2 (Q) + β(c * ) 2 2 meas(Q)
which is a contradiction to (4.1) for all β > 0 sufficiently small.

The optimality system for the regularized problem is given by

             (βu β -B * ωc p β , u -u β ) ≥ 0 for all u with u L ∞ (Q) ≤ 1, A * p β + C * ωo (C ωo y β -z) = 0, p β (T ) = 0, Bp β | Σ = 0, αc β + β 2 u β 2 L 2 (Q) -(u β , B * ωc p β ) = 0, Ay β -c β B ωc u β = 0, y β (0) = y 0 , Cy β | Σ = 0 with p β ∈ Y * .
By pointwise inspection of the first relation we have

u β = sgn β (B * ωc p β ) =    1 B * ωc p β (t, x) > β, -1 B * ωc p β (t, x) < -β, 1 β B * ωc p β (t, x) |B * ωc p(t, x)| ≤ β. (4.2)
Before we prove uniqueness of a solution of (P reg ) we recall the following wellknown property. Let N = { 0 } × L 2 (Q). Then we can introduce the inverse operator

A -1 : N → L 2 , f → y,
where y ∈ X, y(0) = 0, is the unique solution of

Ay, ϕ Y,Y * = (f, ϕ) L2 ∀ϕ ∈ Y * .
By a priori estimates, see, e.g., Lions and Magenes [12, p. 265], we obtain that A -1 is a bounded linear operator. Consequently, there exists a well-defined dual operator (A -1 ) * : L 2 → N satisfying

((A -1 ) * w, v) L2 = (w, A -1 v) L2 (4.3) for w ∈ L 2 and v ∈ N .
Using this property we can guarantee uniqueness of a solution of the regularized problem under certain conditions. The uniqueness is not obvious because of the bilinear structure of the state equation. Lemma 4.1 Let (y β , u β , c β ) be a solution of (P reg ). Then y β and u β are uniquely determined by c β . Conversely, c β and y β are uniquely determined by u β . Further, for α > 0 sufficiently large there exists a unique solution of problem (P reg ).

Proof To prove uniqueness we use a Taylor expansion argument as in [1, Appendix A]. To utilize (4.3) we need to transform (P reg ) into a problem with homogeneous initial condition. For this purpose let ȳ ∈ X be the solution of

     Aȳ = 0 in Q, ȳ(0) = y 0 in Ω, Cȳ = 0 on Σ. (4.4)
We set z = -C ωo ȳ + z and introduce problem (P hom ) given by

                         min (y,u,c)∈X×U ×R + 0 J β (y, u, c) = 1 2 C ωo y -z 2 L 2 (Q) + βc 2 u 2 L 2 (Q) + α 2 c 2 , s.t. Ay = cB ωc u in Q, y(0) = 0 in Ω, Cy = 0 on Σ, u U ≤ 1. (P hom )
The control problems (P reg ) and (P hom ) are equivalent. Thus, without restriction of generality we can assume that the initial state y 0 is zero.

We define the reduced cost

F (u, c) = 1 2 C ωo A -1 (cB ωc u) -z 2 L 2 (Q) + βc 2 u 2 L 2 (Q) + α 2 c 2 .
To shorten notations we set

M = C ωo A -1 B ωc , i.e. M : L 2 (Q) Bω c ----→ N A -1 ----→ L 2 Cω o ----→ L 2 (Q).
Since M is a linear, bounded operator and using (4.3) we derive the optimality conditions

c β (βu β -M * z + c β M * Mu β , u -u β ) ≥ 0 for all u L ∞ (Q) ≤ 1, (4.5
)

αc β + β 2 u β 2 L 2 (Q) -(u β , M * z) + c β Mu β 2 L 2 (Q) = 0. (4.6)
The partial derivatives of F at (u β , c β ) are given by

F uu = c 2 β M * M + βc β id, F cc = Mu β 2 L 2 (Q) + α, F uc = 2c β M * Mu β -M * z + βu β , F uuc = 2c β M * M + β id, F ccu = 2M * Mu β , F uucc = 2M * M.
Let (u, c) be an admissible pair. Then we set

û = u -u β , ĉ = c -c β .
The Taylor expansion is given as follows, we use the fact, that F c (u β , c β ) = 0, see (4.6), and that the derivatives commute

F (u, c) -F (u β , c β ) = c β (βu β -M * z + c β M * Mu β , û) + c 2 β 2 Mû 2 L 2 (Q) + βc β 2 û 2 L 2 (Q) + 1 2 Mu β 2 L 2 (Q) + α ĉ2 + (2c β M * Mu β -M * z + βu β , û)ĉ + 3 6 2c β Mû 2 L 2 (Q) ĉ + β û 2 L 2 (Q) ĉ + (Mu β , Mû)ĉ 2 + 6 24 2 Mû 2 L 2 (Q) ĉ2 .
Using twice (4.5) we further have

F (u, c) -F (u β , c β ) ≥ c 2 β 2 Mû 2 L 2 (Q) + β 2 (c β + ĉ) û 2 L 2 (Q) + + 1 2 Mu β 2 L 2 (Q) + α ĉ2 + c β (Mu β , Mû)ĉ + c β Mû 2 L 2 (Q) ĉ + (Mu β , Mû)ĉ 2 + 1 2 Mû 2 L 2 (Q) ĉ2 . With (Mu β , Mû) = ( √ ηMu β , ( √ η) -1 Mû) ≥ - η 2 Mu β 2 L 2 (Q) - 1 2η Mû 2 L 2 (Q)
we obtain

F (u, c) -F (u β , c β ) ≥ c 2 β 2 (1 - 1 η ) Mû 2 L 2 (Q) + β 2 c û 2 L 2 (Q) + 1 2 (α -η Mu β 2 L 2 (Q) )ĉ 2 + c β Mû 2 L 2 (Q) ĉ. (4.7) Set K := sup { Mu 2 L 2 (Q) | u L ∞ (Q) ≤ 1 } and choose η = 1. For α > K 2 we have F (u, c) -F (u β , c β ) ≥ β 2 c û 2 L 2 (Q) + 1 2 (α -K 2 )ĉ 2 + c β Mû 2 L 2 (Q) ĉ ≥ 0. (4.8)
Let (u β , c β ) and (u β ′ , c β ′ ) be two solutions. Then we obtain from (4.8), that

(u β , c β ) = (u β ′ , c β ′ ).
From (4.7) we see, that for c β = c β ′ also u β = u β ′ for any η ≥ 1 and for u β = u β ′ we have c β = c β ′ for η > 0 sufficiently small. These last two statements do not require any assumption on α.

In the following we analyze convergence of the the solution of (P reg ) for β → 0 and proceed as in [START_REF] Clason | A minimum effort optimal control problem for elliptic PDEs[END_REF]. Lemma 4.2 For β > 0 let (y β , u β , c β ) denote a solution of (P reg ). Further let (P 3 ) have a solution (y * , u * , c * ) which we associate with β = 0 and also denote by (y 0 , u 0 , c 0 ). Then for any 0 ≤ β ≤ β ′ we have

c β ′ u β ′ 2 L 2 (Q) ≤ c β u β 2 L 2 (Q) , (4.9)
J(y β , c β ) ≤ J(y β ′ , c β ′ ), (4.10)
J(y β , c β ) + βc β 2 u β 2 L 2 (Q) ≤ J(y * , c * ) + βc * 2 u * 2 L 2 (Q) . (4.11)
Proof We recall the proof from [START_REF] Clason | A minimum effort optimal control problem for elliptic PDEs[END_REF] and apply it to the time-dependent case. Since (y β , u β , c β ) is a solution of (P reg ) and (y 0 , u 0 , c 0 ) a solution of (P 3 ), respectively, we have for 0

≤ β ≤ β ′ that J(y β , c β ) + βc β 2 u β 2 L 2 (Q) ≤ J(y β ′ , u β ′ ) + βc β ′ 2 u β ′ 2 L 2 (Q) .
Thus, further

J(y β , c β ) + βc β 2 u β 2 L 2 (Q) + (β ′ -β)c β ′ 2 u β ′ 2 L 2 (Q) ≤ J(y β ′ , c β ′ ) + β ′ c β ′ 2 u β ′ 2 L 2 (Q) ≤ J(y β , c β ) + β ′ c β 2 u β 2 L 2 (Q) .
(4.12)

From the outer inequality we have (β

′ -β)(c β ′ u β ′ 2 L 2 (Q) -c β u β 2 L 2 (Q) )
≤ 0 implying the first assertion. From (4.12) we derive

J(y β , c β ) -J(y β ′ , c β ′ ) ≤ β(c β ′ u β ′ 2 L 2 (Q) -c β u β 2 L 2 (Q) ) (4.13)
and the right hand side is smaller than or equal to zero by the previous result and thus (4.10) follows. Assertion (4.11) follows from the last inequality in (4.12) by setting β ′ = β and β = 0.

After this preparation we prove convergence of minimizers of (P reg ) following [START_REF] Clason | A minimum effort optimal control problem for elliptic PDEs[END_REF]. Theorem 4.1 Let (P 3 ) have a unique solution. Then any selection of solutions

{ (y β , u β , c β ) } β>0 of (P reg ) are bounded in X × L ∞ (Q) × R +
0 and converges weak * to the solution of (P 3 ) for β → 0. It converges strongly in

X × L q (Q) × R + 0 for q ∈ [1, ∞).
Proof Let ŷ ∈ X be the solution of (4.4). The point (ŷ, û, ĉ) = (ŷ, 0, 0) is feasible for the constraints. Thus we have

C ωo y β -z 2 L 2 (Q) + βc β u β 2 L 2 (Q) + αc 2 β ≤ C ωo ŷ -z 2 L 2 (Q)
and consequently, the boundedness of c β follows. The controls u β are bounded by the constant 1 in L ∞ (Q) and hence, y β is bounded in X.

Therefore, there exists (ȳ, ū, c)

∈ X × L ∞ (Q) × R + 0 such that for a subse- quence there holds (y β , u β , c β ) ⇀ * (ȳ, ū, c) in X × L ∞ (Q) × R.
By passing to the limit in the equation we obtain that

(ȳ, ū, c) is a solution of      Ay = cB ωc u in Q, y(0) = y 0 in Ω, Cy = 0 on Σ. Since the L ∞ -norm is weak * lower semicontinuous, we have ū L ∞ (Q) ≤ 1.
Further, by the weak lower semicontinuity of

J β : L 2 (Q) × L 2 (Q) × R + 0 → R we derive that (ȳ, ū, c
) is a solution of (P 3 ). Uniqueness of the solution of (P 3 ) implies that (ȳ, ū, c) = (y * , u * , c * ). Thus we have proved weak * convergence.

For strong convergence insert the weak limit (u * , c * ) in (4.9) with β = 0 and obtain for all β ′ > 0 from the lower semicontinuity of the norm that

c β ′ u β ′ 2 L 2 (Q) ≤ c * u * 2 L 2 (Q) ≤ lim inf β ′ →0 c β ′ u β ′ 2 L 2 (Q) .
This implies lim sup

β ′ →0 u β ′ 2 L 2 (Q) ≤ u * 2 L 2 (Q) ≤ lim inf β ′ →0 u β ′ 2 L 2 (Q)
and consequently, strong convergence in L 2 (Q). Using

u β ′ -u * L p (Q) ≤ u β ′ -u * L 2 (Q) u β ′ -u * L ∞ (Q) we obtain strong convergence of u β ′ → u * in every L q (Q), q ∈ [1, ∞).
Furthermore, this implies strong convergence of the corresponding state y β ′ .

Remark 4.3 In case of non-unique solvability of (P 3 ) the assertion of Theorem 4.1 remains true if we consider the convergence of subsequences to a solution of (P 3 ).

From the strong convergence of u β we can derive a convergence rate for the error in the cost functional.

Corollary 4.1 Let (P 3 ) have a unique solution (y * , u * , c * ). Then there holds

J(y β , c β ) -J(y * , c * ) = o(β)
for β → 0.

Proof From (4.13) we have

0 ≤ J(y β , c β ) -J(y * , u * ) ≤ β(c β u β L 2 (Q) -c * u * L 2 (Q) )
which proves the assertion.

From know on we will assume, that problem (P reg ) has a unique solution.

5 Semi-smooth Newton method

In this section we formulate the semi-smooth Newton method and prove its superlinear convergence. To keep notations simple we omit the index β for the solution of the regularized problem. Using

p L 1 β (Q) = Q |p(t, x)| β dxdt, |p(t, x)| β =    p(t, x) -β 2 if p(t, x) > β, -p(t, x) -β 2 if p(t, x) < -β, 1 2β p(t, x) 2 if |p(t, x)| ≤ β
we reformulate the optimality system for the regularized problem. We eliminate the control u using (4.2) and obtain

A * p + C * ωo (C ωo y -z) = 0, p(T ) = 0, Bp| Σ = 0, (5.1) αc -B * ωc p L 1 β (Q) = 0, (5.2)
Ay -cB ωc sgn β (B * ωc p) = 0, y(0) = y 0 Cy| Σ = 0.

(5.3) To write the system equivalently as an operator equation we set

W = X × Y * 0 × R + , Z = X * × R × Y × Y 1 .
For convenience of the reader we recall that

Y 1 = H 1 0 (Ω) × L 2 (Ω), Y 0 = L 2 (Ω) × H -1 (Ω), X = L 2 (Y 1 ) ∩ H 1 (Y 0 ), P 1 = (Y 0 ) * , P 0 = (Y 1 ) * , Y * 0 = { p ∈ L 2 (P 1 ) ∩ H 1 (P 0 ) | p(T ) = 0 }, and Y = L 2 (Y 0 ) ⊕ H 1 (P 0 ) * .
Then, we can define the operator T by

T : W -→ Z, T(x) = T(y, p, c) =     A * p + C * ωo C ωo y -C * ωo z αc -B * ωc p L 1 β (Q) Ay -cB ωc sgn β (B * ωc p) y(0) -y 0     (5.4)
and obtain (5.1)-( 5.3) equivalently as

T(x) = 0 (5.5) for x ∈ W .
To formulate the semi-smooth Newton method we need Newton differentiability of the operator T. Let

W R = X × Y * 0 × R.
Lemma 5.1 The operator T is Newton differentiable, i.e. for all x ∈ W and h ∈ W R there holds

T(x + h) -T(x) -T ′ (x + h)h Z = o(h) for h W R → 0. (5.6)
Proof The operator max : For the min operator an analog Newton derivative can be obtained. Since

L p (Q) → L q (Q), p > q ≥ 1 is Newton differentiable with derivative (D N max(0, v -β)h)(t, x) = h(t, x), v(t, x) > β, 0, v(t, x) ≤ β for v, h ∈ L p (Q), β ∈ R,
sgn β (v) = 1 β (v -max(0, v -β) -min(0, v + β))
we obtain for the operator

sgn β : L p (Q) → L q (Q), p > q ≥ 1,
the Newton derivative

D N sgn β (p)h (t, x) = 0, |p(t, x)| > β, 1 β h(t, x), |p(t, x)| ≤ β for v, h ∈ L p (Q), β ∈ R + ,
• L 1 β (Q) : L p (Q) → R, p ≥ 4
is Newton differentiable with Newton derivative

D N ( v L 1 β (Q) )h = (sgn β (v), h) for v, h ∈ L p (Q),
p ∈ C(H 1 (Ω)) ֒→ L q (Q) for q = 2d d-2 the mappings p → B * ωc p → B * ωc p L 1 β (Q) , Y * → L 4 (Q) → R, p → B * ωc p → cB ωc sgn β (B * ωc p), Y * → L 4 (Q) → L 2 (Q) → X * (5.7)
for d ≤ 4 are Newton differentiable. Consequently, we obtain the assertion.

To formulate the semi-smooth Newton method we consider

T ′ : W → L(W R , Z) (5.8)
with

T ′ (x)(δy, δp, δc) =     A * δp + C * ωo C ωo δy αδc -(sgn β (B * ωc p), B * ωc δp) Aδy -δcB ωc sgn β (B * ωc p) -c β B ωc B * ωc δpχ Ip δy(0)     (5.9)
for x = (y, p, c) ∈ W and (δy, δp, δc) ∈ W R . Here χ Ip denotes the characteristic function for the set

I p = { (t, x) ∈ Q | |B * ωc p(t, x)| ≤ β } . (5.10)
The operator T ′ (x) is invertible on its image as we see in the next lemma. The proof is presented in the appendix. Lemma 5.2 For x ∈ W the operator

T ′ (x) : W R → Im(T ′ (x)) ⊂ Z
is bijective and we can define

T ′ (x) -1 : Im(T ′ (x)) → W R .
Furthermore, there holds the following estimate

T ′ (x) -1 (z) W R ≤ C z Z (5.11) for z ∈ Im(T ′ (x)) ∩ Z 1 uniformly in x ∈ W ,
where

Z 1 = S × R × M × { 0 } ⊂ Z and S = { (χ ωo v, 0) | v ∈ L 2 (L 2 (Ω)) }, M = { (0, v) | v ∈ L 2 (H -1 (Ω)) }.
Directly applying the Newton method to equation (5.5) leads to the iteration

T ′ (x k )(δx) = -T(x k ),
(5.12)

x k+1 = δx + x k , x 0 ∈ W, (5.13)
where in every Newton step the following system

A * δp + C * ωo C ωo δy = -A * p k -C * ωo C ωo y k + C * ωo z, (5.14) δp(T ) = 0, Bδp| Σ = 0, αδc -(sgn β (B * ωc p k ), B * ωc δp) = -αc k + B * ωc p k L 1 β (Q) , Aδy -δcB ωc sgn β (B * ωc p k ) - c k β B ωc B * ωc δpχ I k = -Ay k + c k B ωc sgn β (B * ωc p k ), (5.15
)

δy(0) = -y k (0) + y 0 , (5.16
)

Cδy| Σ = 0
with I k = I p k has to be solved. To simplify the system we reformulate it equivalently as follows

A * p k+1 + C * ωo C ωo y k+1 = C * ωo z, (5.17)
p k+1 (T ) = 0, (5.18)

Bp k+1 | Σ = 0, (5.19) αc k+1 -(sgn β (B * ωc p k ), B * ωc p k+1 ) = -(sgn β (B * ωc p k ), B * ωc p k ) + B * ωc p k L 1 β (Q) , (5.20
)

Ay k+1 -c k+1 B ωc sgn β (B * ωc p k ) - c k β B ωc B * ωc p k+1 χ I k = - c k β B ωc B * ωc p k χ I k , (5.21)
y k+1 (0) = y 0 , (5.22)
Cy| Σ = 0, (5.23)

δy = y k+1 -y k (5.24)
for k ∈ N 0 . The iterates p k+1 are solutions of wave equations with right hand side in L 2 (Q) × { 0 } and the iterates y k+1 are solutions of wave equations with right hand side in { 0 } × L 2 (Q) for all k ≥ 0 and initialization (y 0 , p 0 , c 0 ) ∈ W . Under certain conditions the well-definedness of the Newton iteration can be shown. Lemma 5.3 For x 0 ∈ W the Newton iterates x k satisfy

x k ∈ W, T(x k ) ∈ Im(T ′ (x k )) for k ∈ N 0 if c k > 0.
Remark 5.1 Let x * be the solution of (P reg ). In Theorem 5.1 we will show that for β and x 0 -x * W R sufficiently small the iterates c k remain positive.

Proof of Lemma 5.3 For given iterate x k ∈ W we consider the control problem

                           min (y,u,c)∈X×U ×R + 0 J(c, u, y) = 1 2 C ωo y -z 2 L 2 (Q) + βc k 2 u 2 L 2 (Q) + α 2 |c -z 1 | 2 , s.t. Ay -cB ωc sgn β (B * ωc p k ) -c k B ωc uχ I k = z 2 in Q, y(0) = y 0 in Ω, Cy| Σ = 0 on Σ (5.25) with I k = I p k , z 2 ∈ L 2 ({ 0 } × L 2 (Ω)), z 1 ∈ R, y 0 ∈ Y 1 , x k = (y k , p k , c k
) and α, β, y 0 as in (5.17)- (5.23). This problem has a unique solution (y, u, c).

The optimality system of (5.25) is given by (5.17)-( 5.23) if we choose

z 1 = - 1 α (sgn β (B * ωc p k ), B * ωc p k ) -B * ωc p k L 1 β (Q) , (5.26)
z 2 = - c k β B ωc B * ωc p k χ I k .
(5.27) From (5.17)-( 5.19) we derive that p ∈ Y * 0 . This implies x k ∈ W for all Newton iterates. Since (5.12)-(5.13) and (5.17)-(5.24) are equivalent, the second assertion follows, when setting δy = y -y k .

To apply (5.11) we need T(x k ) ∈ Z 1 . For k ≥ 1 this follows immediately from (5.17)-(5.23). To obtain T(x 0 ) ∈ Z 1 , we choose x 0 = (y 0 , p 0 , c 0 ) ∈ W , such that y 0 (0) = y 0 ,

(5.28)

∂ t y 0 1 = y 0 2 ,
(5.29)

A * p 0 + C * ωo C ωo y 0 -C * ωo z = 0. (5.30)
To prove superlinear convergence of the Newton method we need the following estimate.

Lemma 5.4 Let x * ∈ W be the solution of (P reg ) and let x 0 ∈ W satisfy (5.28)-(5.30). Then the Newton iterates satisfy

x k+1 -x * = -T ′ (x k ) -1 T(x k ) -T(x * ) -T ′ (x k )(x k -x * ) (5.31)
and there holds the following estimate

x k+1 -x * W R ≤ C T(x k ) -T(x * ) -T ′ (x k )(x k -x * ) Z (5.32) for k ∈ N 0 if c k > 0 with x k = (y k , p k , c k ).
Proof There holds T(x * ) = 0 and T ′ (x k )(x k -x * ) ∈ Im(T ′ (x k )), and according to Lemma 5.3 we have T(x k ) ∈ Im(T ′ (x k )). Consequently,

T(x * ) -T(x k ) -T ′ (x k )(x k -x * ) ∈ Im(T ′ (x k ))
for all k ∈ N 0 . Further, we derive from (5.17)-( 5.24) and (5.28)-(5.30) that for k ∈ N 0

T(x * ) -T(x k ) -T ′ (x k )(x k -x * ) ∈ Im(T ′ (x k )) ∩ Z 1 .
Thus, the assertion follows with Lemma 5.2.

The superlinear convergence of the Newton method is shown in the next main theorem.

Theorem 5.1 Let x * = (y * , p * , c * ) be the solution of (P reg ) with z ≡ 0 and β sufficiently small, such that c * > 0 (cf. Remark 4.2). Further let x 0 ∈ W satisfy (5.28)-(5.30) and let x 0 -x * W R be sufficiently small. Then the iterates x k = (y k , p k , c k ) ∈ W of the semi-smooth Newton method (5.12)-(5.13) are well defined and they satisfy To realize the semi-smooth Newton method we introduce the active sets

x k+1 -x * W R ≤ o( x k -x * W R ) (5.33) for x k -x * W R → 0. Proof From
A + k = { (t, x) ∈ Q | B * ωc p k (t, x) > β } , A - k = { (t, x) ∈ Q | B * ωc p k (t, x) < -β } for iterates p k ∈ Y * . With I k = I p k (cf. the definition in (5.10)) we have Q = I k ∪ A + k ∪ A - k .
The Newton method is realized as presented in Algorithm 5.1. Remark 5.2 The solution of system (5.17)-(5.24) in Step 7 of Algorithm 5.1 can be found by solving the control problem (5.25) if we assume that the scalar c is always positive. Algorithm 5.1 Semi-smooth Newton algorithm with path-following 1: Choose n = 0, y 0 = (y 01 , y 02 ) ∈ X satisfying (5.28) and (5.29), c 0 ∈ R + , q ∈ (0, 1), tol, tol β , β 0 ∈ R + , and n ∈ N. 2: For given y 0 solve the adjoint equation (5.1) and obtain p 0 ∈ Y * . 3: repeat 4:

Set k = 0 and (y 0 , p 0 , c 0 ) = (y n , p n , c n ). 5: repeat 6:

Compute the active and inactive sets A + k , A - k , and I k :

A + k = { (t, x) ∈ Q | B * ωc p k (t, x) > β } , I k = { (t, x) ∈ Q | |B * ωc p k (t, x)| ≤ β } , A - k = { (t, x) ∈ Q | B * ωc p k (t, x) < -β } .

7:

Solve for x k = (y k , p k , c k ) system (5.17)-( 5.24) and obtain

x k+1 = (y k+1 , p k+1 , c k+1 ).

8:

Set k = k + 1. 9:

until x k -x k-1 W R < tol. 10:
Set (y n+1 , p n+1 , c n+1 ) = x k .

11:

Compute u k+1 = sgn β (B * ωc (p k+1 )). 12:

Set β n+1 = qβn.

13:

Set n = n + 1. 14: until β n+1 < tol β or n > n.

Discretization

To realize Algorithm 5.1 numerically we present the discretization of problem (5.25) for data given by (5.26)-(5.27).

For the discretization of the state equation we apply a continuous Galerkin method following Kröner, Kunisch, and Vexler [START_REF] Kröner | Semi-smooth Newton methods for optimal control of the wave equation with control constraints[END_REF]. For temporal discretization we apply a Petrov-Galerkin method with continuous piecewise linear ansatz functions and discontinuous (in time) piecewise constant test functions. For the spatial discretization we use conforming linear finite elements. Let

J = {0} ∪ J 1 ∪ • • • ∪ J M be a partition of the time interval J = [0, T ] with subintervals J m = (t m-1 , t m ] of size k m and time points 0 = t 0 < t 1 < • • • < t M -1 < t M = T.
We define the time discretization parameter k as a piecewise constant function by setting k| Jm = k m for m = 1, . . . , M . Further, for

0 = l 0 < l 1 < • • • < l N -1 < l N = L let T h = L 1 ∪ • • • ∪ L N
be a partition of the space interval Ω = (0, L) with subintervals L n = (l n-1 , l n ) of size h n and h = max n=1,...,M h n . We construct on the mesh T h a conforming finite element space V h in a standard way by setting

V h = v ∈ H 1 0 (Ω) v| Ln ∈ P 1 (L n ) .
Then the discrete ansatz and test space are given by

X kh = { v ∈ C( J, L 2 (Ω)) | v| Jm ∈ P 1 (J m , V h ) } , X kh = { v ∈ L 2 (J, H 1 0 (Ω)) | v| Jm ∈ P 0 (J m , V h ) and v(0) ∈ L 2 (Ω) } ,
where P r (J m , V h ) denotes the space of all polynomials of degree lower or equal r = 0, 1 defined on J m with values in V h . For the discretization of the control space we set U kh = X kh .

In the following we present the discrete optimality system for (5.25) assuming that the iterates c k are positive. With the notation

(•, •) Jm := Jm (•, •) L 2 (Ω) dt for the Newton iterates c k ∈ R, u k kh ∈ U kh , y k kh = (y k 1 , y k 2 ) ∈ X kh × X kh , and p k kh = (p k 1 , p k 2 ) ∈ X kh × X kh , k ∈ N, the adjoint equation is given by - M -1 m=0 (ψ 1 (t m ), p k+1 1 (t m+1 ) -p k+1 1 (t m )) L 2 (Ω) + (∇ψ 1 , ∇p k+1 2 ) + (ψ 1 (t M ), p k+1 1 (t M )) L 2 (Ω) = -(ψ 1 , χ ωo (y k+1 1 -z)) ∀ψ 1 ∈ X kh , (6.1)
- M -1 m=0 (ψ 2 (t m ), p k+1 2 (t m+1 ) -p k+1 2 (t m )) L 2 (Ω) -(ψ 2 , p k+1 1 ) + (ψ 2 (t M ), p k 2 (t M )) L 2 (Ω) = 0 ∀ψ 2 ∈ X kh , (6.2)
the optimality conditions by

αc k+1 -(sgn β (χ ωc p k 2 ), χ ωc p k+1 2 ) = -(sgn β (χ ωc p k 2 ), χ ωc p k 2 ) + χ ωc p k 2 L 1 β (Q) , (6.3)
(βu k+1 kh , τ u) = (χ I kh χ ωc p k+1 2 , τ u) ∀τ u ∈ U kh , (6.4)
for I kh = I p kh and the state equation by

M m=1 (∂ t y k+1 1 , ξ 1 ) Jm -(y k+1 2 , ξ 1 )+(y k+1 1 (0)-y 0,1 , ξ 1 (0)) L 2 (Ω) = 0 ∀ξ 1 ∈ X kh , (6.5)
M m=1 (∂ t y k+1 2 , ξ 2 ) Jm + (∇y k+1 1 , ∇ξ 2 ) + (y k+1 2 (0) -y 0,2 , ξ 2 (0)) L 2 (Ω) -c k+1 (sgn β (χ ωc p k 2 ), ξ 2 ) -c k (χ ωc u k+1 kh χ I kh , ξ 2 ) = -c k (sgn β (χ ωc p k 2 ), ξ 2 ) ∀ξ 2 ∈ X kh (6.6)
with y 0 = (y 0,1 , y 0,2 ). When evaluating the time integrals by a trapedoizal rule the time stepping scheme for the state equation results in a Crank Nicolson scheme.

To solve the system (6.1)-(6.6) we introduce the control-to-state operator for the discrete state equation (6.5)-(6.6)

S k kh : U kh × R → L 2 (Q), (u kh , c) → y 1
and the discrete reduced cost functional

j k kh : U kh × R → R + 0 , j k kh (u kh , c) = 1 2 χ ωo S kh (u kh , c) 2 L 2 (Q) + βc 2 u kh 2 U + α 2 |c -z 1 | 2 ,
with

z 1 = - 1 α (sgn β (χ ωc p k 2 ), χ ωc p k 2 ) -χ ωc p k 2 L 1 β (Q) ,
where p k kh = (p k 1 , p k 2 ) results from the previous iterate. Then the solution of the system is given as a solution of the reduced problem min j k kh (u kh , c), (u kh , c) ∈ U kh × R.

The necessary optimality condition is given by

(j k kh ) ′ (u kh , c)(δu, δc) = 0 ∀(δu, δc) ∈ U kh × R.
We solve this reduced problem by a classical Newton method, i.e. the Newton update (τ u, τ c) ∈ U kh × R is given by

(j k kh ) ′′ (u kh , c)(τ u, τ c, δu, δc) = -(j k kh ) ′ (u kh , c)(δu, δc) ∀(δu, δc) ∈ U kh × R. (6.7)
The explicit representations of the derivatives of the reduced cost functional are given in the Appendix 8.2.

Numerical examples

In this section we present numerical examples confirming the theoretical results from above. In the first three examples we consider the convergence behaviour of the Newton iteration in the inner loop of Algorithm 5.1, i.e. we consider the case with path iteration number n = 0. Further, we present an example in which we consider the algorithm with n large and analyze the behaviour for β → 0. The computations are done by using MATLAB, for the plot in Figure 7.3 the optimization library RoDoBo [START_REF] Rodobo | A C++ library for optimization with stationary and nonstationary PDEs with interface to[END_REF] was used.

Example 7.1 Let the data be given as follows z(t, x) = -sin(2πx), α = 10 -2 , β = 10 -3 , y 0 (x) = (x(1 -x), 0), n = 0 for x ∈ Ω = (0, 1) and T = 1. The control and observation area is given by ω o = (0, 1), ω c = (0, 1).

As an initial point for the algorithm we choose y 0 (t, x) = (x(1 -x), 0) , c 0 = 10 satisfying (5.28) and (5.29). We discretize our problem as presented in the previous section on uniform meshes with N = 256 and M = 255. In Table 7.1 we see the errors in the scalar e k c = |c k -c * |, in the state e k y = y k -y * X and in the adjoint state e k p = p k -p * L 2 (P 1 )∩L 2 (P 0 ) in every Newton iteration k. For the exact solution (y * , p * , c * ) we choose the 8th iterate. We do not consider the full norm of Y * 0 for the adjoint state, since we discretize the adjoint state by piecewise constants in time. By am we denote the number of mesh points in set A -and by ap the number of mesh points in set A + . As the stopping criterion for the Newton iteration we choose tol = 10 -9 . If we go beyond this tolerance the residuums in the conjugate gradient method to solve the Newton equation (6.7) reach the machine accuracy. The behaviour of the errors presented in Table 7.1 indicate superlinear convergence.

Example 7.2 In this example we keep the data as above except for ω o = (0, 1), ω c = (0, 1/3), i.e. the control domain is a subset of the domain of observation, cf. Lemma 3.2. In Table 7.2 the behaviour of the errors of the Newton iterates are shown. For the exact solution we choose the 6th Newton iterate and as in the previous superlinear convergence.

We note that in these three examples above am and ap are identified before we stop. In fact not only the cardinality of the sets A -and A + stagnates but the sets themselves are identified. Example 7.4 In this example we apply a simple path-following strategy by choosing in every iteration the new regularization parameter by the rule

β n+1 = qβ n , n ∈ N 0 ,
with some given q ∈ (0, 1) and β 0 > 0.

We choose z = -1, α = 10 -2 , β 0 = 10 -1 , y 0 (x) = (sin(2πx), -4 sin(2πx)), n = 6 (7.1) for x ∈ Ω = (0, 1) and T = 1. Further we set q = 0.2 and ω c = ω o = Ω. We solve the problem on a spatial and temporal mesh with N = 100 and M = 127. For initialization we choose y 0 (t, x) = ((t -1) 4 sin(2πx), 4(t -1) 3 sin(2πx)), c 0 = 1 (7.2) for (t, x) ∈ Q. The results are presented in Table 7.4. For decreasing β the corresponding values of the cost functional and the behaviour of the error e J βn = J(u βn , c βn ) -J(u * , c * ) is shown. For the exact solution (u * , c * ) we take (u β6 , c β6 ). Further the number of Newton steps ns is presented. The values of the cost functional decrease which confirms the theoretical result in (4.10). Further, the behaviour of the errors indicates superlinear convergence for β → 0, which confirms the result of Corollary 4.1. The number of active points in A -is larger than in A + which we expect for the given desired state.

For β smaller than presented in Table 7.4 the number of active and inactive nodes remains constant up to 3 switching nodes, however one looses the superlinear convergence. The number of Newton steps decreases which relies on the fact that the iteration is nested. For a non nested iteration (i.e. (y n , p n , c n ) = (y 0 , p 0 , c 0 ) for all n) the number of Newton steps is increasing for smaller β, see Table 7.5. In Figure 7.1 we compare for time horizon T = 2 the state of the regularized control problem for data given in (7.1), (7.2) and β = 4•10 -3 with the solution of the state equation for u ≡ 0. The plots show the behaviour of the state with If we go beyond the time horizon T = 2 the tracking of the desired state by the optimal state of the regularized problem further improves.

In Figure 7.2 we see the corresponding optimal control of the regularized problem which is nearly of bang-bang type.

Figure 7.3 shows the optimal state for problem (P 1 ) when replacing the L ∞ -by L 2 -control costs with α given as in (7.1). The tracking of the desired state is nearly the same as in case of the regularized problem presented in Figure 7.1. But we see that in some parts the deflection in positive direction is less than for the regularized problem. This reflects our expectation, since the L 2 -control space is larger than the L ∞ -space and thus allows a better approximation of the desired state.

We also tested the case with ω o = [0, 1/3] and ω c = [2/3, 1] and observed similar numerical behaviour for different initializations y 0 and c 0 . (8.1) with x = (y, p, c) and z 0 ≡ 0. Existence of a unique solution follows by considering the reduced functional j(δu, δc) = J(δc, δu, δy(δu, δc)), where δy is the solution to the constraining partial differential equation as a function of (δu, δc). The solution is necessarily zero.

In the second step we prove the estimate (5.11). Let x = (y, p, c) ∈ W , δx = (δy, δp, δc) ∈ W R and z = (z 0 , z 1 , z 2 , 0) ∈ Im(T ′ (x)) ∩ Z 1 . Then the equation T ′ (x) -1 (z) = δx is equivalent to the following system

A * δp + C * ωo C ωo δy = C * ωo z 0 , (8.2

  and (t, x) ∈ Q. The mapping w : R → R, s → |w(s)| β defines a differentiable Nemytskii operator from L p (Q) to L 2 (Q) for p ≥ 4 according to Tröltzsch [15, Chapter 4.3.3]. Thus, the mapping

  the estimates (5.6) and (5.32) and the positiveness of c * > 0 we conclude that c k > 0 for all k ∈ N if x 0 -x * W R is sufficiently small. Thus by Lemma 5.3 all iterates are in W . Estimate (5.33) follows from Lemma 5.1, Lemma 5.4 and [8, Proof of Theorem 8.16].

Fig. 7 . 1 :

 71 Fig. 7.1: State for the controlled and uncontrolled (u ≡ 0) problem

Fig. 7 Fig. 7 . 3 : 8 Appendix 8 . 1 2 C ωo δy -z 0 2 L 2 ( 2 ,

 7738812222 Fig. 7.2: Optimal control

2 )ωo δy 2 L 2 (= -z 2 2 L 2 -χ ωc p k 2 L 1 β 2 ), p k+1 2 ) 2 ), δp k+1 2 )- 1 )

 2222222122221 ) δp(T ) = 0,Bδp| Σ = 0, αδc -(sgn β (B * ωc p), B * ωc δp) = z 1 ,(8.3)Aδy -δcB ωc sgn β (B * ωc p) -c β B ωc B * ωc δpχ Ip = z 2 ,with δy and (8.4) with -δp and adding both equations we obtain C Q) + δc(sgn β (B * ωc p), B ωc δp) , δp Y,Y * + (C * ωo z 0 , δy). (8.5)Here we used (2.2) and that δy(0) = 0 and δp(T ) = 0. By multiplying (8.3) with -1 α (sgn β (B * ωc p), B * ωc δp) and adding it to (8.5) we haveC ωo δy (Q) + 1 α (sgn β (B * ωc p), B ωc δp) 2 ≤ z 2 Y δp Y * + C * ωo z 0 L 2 (Q) C ωo δy L 2 (Q) + 1 α |(sgn β (B * ωc p), B * ωc δp)| |z 1 |.From the priori estimate in[12, p. 265] we haveδp Y * ≤ C C * ωo (C ωo δy -z 0 ) L2 (8.6)then the first and second derivative of j kh at a point (u kh , c) ∈ U kh × R are given by(j k kh ) ′ (u kh , c)(δu, δc) = βc k (u kh , δu) + α(c + 1 α (sgn β (χ ωc p k 2 ), χ ωc p k 2 ))δc (Q) δc -δc(sgn β (χ ωc p k -c k (χ ωc δuχ I kh , pk+1kh ) ′′ (u kh , c)(τ u, τ c, δu, δc) = βc k (τ u, δu)+αδc•τ c-δc(sgn β (χ ωc p k c k (χ ωc τ uχ I kh , δp k+1for δu, τ u ∈ U kh and δc, τ c ∈ R.

  see Clason, Ito, and Kunisch [1]. Further, since B * ωc

Table 7 . 1 :

 71 Error of the Newton iterates

	k	c	e k c	e k c /e	(k-1) c	e k y	e k y /e	(k-1) y	e k p	e k p /e	(k-1) p	am	ap
	0	10.0	8.49		-	9.63		-	0.74		-	42025	17450
	1	1.214853 2.95e-01	3.48e-02	8.30e+00	8.63e-01	2.32e-01	3.12e-01	32298	25740
	2	1.393382 1.17e-01	3.95e-01	7.51e-01	9.05e-02	3.12e-02	1.34e-01	22186	33754
	3	1.489242 2.08e-02	1.78e-01	3.14e-01	4.18e-01	8.52e-03	2.73e-01	26101	31073
	4	1.509587 4.43e-04	2.13e-02	3.28e-02	1.04e-01	1.04e-03	1.22e-01	24874	31490
	5	1.510022 8.19e-06	1.85e-02	1.48e-03	4.52e-02	4.40e-05	4.22e-02	25120	31450
	6	1.510031 9.14e-08	1.12e-02	5.53e-06	3.73e-03	9.35e-08	2.13e-03	25133	31448
	7	1.510031 1.62e-10	1.77e-03	3.70e-10	6.69e-05	6.97e-12	7.45e-05	25133	31448
	8	1.510031	0.00	0.00	0.00	0.00	0.00	0.00	25133	31448

Table 7 . 2 :

 72 Error of the Newton iterates In this example we choose the data as above except forω o = (1/2, 1), ω c = (0, 1/3),i.e. ω c ⊂ ω o . Further, we set tol = 10 -7 for the reason already mentioned in Example 7.1. The behaviour of the errors of the Newton iterates is presented in Table7.3. As the exact solution we take the 4th iterate and again we obtain

	k	c	e k c	e c k /e c (k-1)	e k y	e y k /e y (k-1)	e k p	e p k /e p (k-1)	am	ap
	0	10.0	9.62	-	11.64	-	8.67e-01	-	19923	0
	1	0.346298 2.88e-02	3.00e-03	4.43e-01	3.80e-02	1.71e-02	1.97e-02	15559	2782
	2	0.374087 1.06e-03	3.67e-02	2.81e-02	6.35e-02	1.21e-03	7.10e-02	15892	1699
	3	0.375135 9.41e-06	8.89e-03	2.13e-04	7.59e-03	8.98e-06	7.40e-03	15875	1784
	4	0.375145 4.24e-09	4.51e-04	5.15e-08	2.41e-04	1.64e-09	1.83e-04	15875	1785
	5	0.375145 3.36e-13	7.93e-05	1.31e-12	2.54e-05	4.86e-14	2.96e-05	15875	1785
	6	0.375145	0.00	0.00	0.00	0.00	0.00	0.00	15875	1785
	example the iterates converge superlinearly.
	Example 7.3

Table 7 . 3 :

 73 Error of the Newton iterates

	k	c	e k c	e k c /e k-1 c	e k y	e k y /e k-1 y	e k p	e k p /e k-1 p	am	ap
	0	10.0	9.45	-	10.67	-	0.41	-	0	10075
	1	0.491219	1.28e-02	1.35e-03	3.91e+00	3.67e-01	1.81e-02	4.36e-02	0	11444
	2	0.503892	1.57e-04	1.22e-02	1.07e-02	2.73e-03	2.63e-05	1.45e-03	0	11927
	3	0.504049	1.10e-08	7.01e-05	2.01e-08	1.89e-06	4.69e-10	1.78e-05	0	11927
	4	0.504049	0.00	0.00	0.00	0.00	0.00	0.00	0	11927

Table 7 . 4 :

 74 Error in the cost functional

	n	βn	J(u e J βn /e J β n-1	c βn	am	ap	ns
	0	1.0	0.25171	3.77e-02	-	2.94061	2590	0	3
	1	2.0e-02	0.22346	9.47e-03	0.25	3.70135	6833	0	3
	2	4.0e-03	0.21623	2.24e-03	0.24	3.96424	8720	722	3
	3	8.0e-04	0.21444	4.53e-04	0.20	4.04065	9702	1563	2
	4	1.6e-04	0.21407	8.57e-05	0.18	4.05790	10164	1915	2
	5	3.2e-05	0.21401	1.47e-05	0.17	4.06136	10371	2054	2
	6	6.4e-06	0.21400	0.00	0.00	4.06209	10520	2120	1

βn , c βn ) J(u βn , c βn ) -J(u * , c * )

Table 7 . 5 :

 75 Newton steps for decreasing β

	β	1.0e-00	2.0e-02	4.0e-03	8.0e-04	1.6e-04	3.2e-05	6.4e-06
	ns	3	5	5	5	6	7	12

The second author was supported in part by the Austrian Science Fund (FWF) under grant SFB F32 (SFB "Mathematical Optimization and Applications in Biomedical Sciences").

with Y * = L 2 (P 1 ) ∩ H 1 (P 0 ). Using Young's inequality we further derive

and hence,

This implies

and together with (8.3)

Finally, from (8.4), (8.8), (8.9)

Tangent and additional adjoint equations

Let

2 ) be the solution of the tangent equation

and

) of the additional adjoint


Loic Zussa
email: loic.zussa@mines-stetienne.fr

Amine Dehbaoui
email: amine.dehbaoui@mines-stetienne.fr

Karim Tobich
email: karim.tobich@lirmm.fr

Jean-Max Dutertre
email: jean-max.dutertre@mines-stetienne.fr

Philippe Maurine
email: philippe.maurine@lirmm.fr

Ludovic Guillaume-Sage
email: ludovic.guillaume-sage@lirmm.fr

Jessy Clediere
email: jessy.clediere@cea.fr

Assia Tria
email: assia.tria@cea.fr

Efficiency of a Glitch Detector against Electromagnetic Fault Injection

The use of electromagnetic glitches has recently emerged as an effective fault injection technique for the purpose of conducting physical attacks against integrated circuits. First research works have shown that electromagnetic faults are induced by timing constraint violations and that they are also located in the vicinity of the injection probe. This paper reports the study of the efficiency of a glitch detector against EM injection. This detector was originally designed to detect any attempt of inducing timing violations by means of clock or power glitches. Because electromagnetic disturbances are more local than global, the use of a single detector proved to be inefficient. Our subsequent investigation of the use of several detectors to obtain a full fault detection coverage is reported, it also provides further insights into the properties of electromagnetic injection and into the key role played by the injection probe.

I. INTRODUCTION

Since the early warning of [START_REF] Quisquater | Eddy current for Magnetic Analysis with Active Sensor[END_REF] [START_REF] Quisquater | Eddy current for Magnetic Analysis with Active Sensor[END_REF], the use of electromagnetic (EM) glitches has recently emerged as an effective fault injection technique for the purpose of conducting physical attacks against ICs [START_REF] Dehbaoui | Electromagnetic transient faults injection on a hardware and a software implementations of aes[END_REF], [START_REF] Dehbaoui | Electromagnetic glitch on the aes round counter[END_REF], [START_REF] Bayon | Contactless electromagnetic active attack on ring oscillator based true random number generator[END_REF]. These latter works indicate that the mechanism related to the injection of faults involves timing constraint violations. The timing constraint violation is induced by a transient underpowering of the target created by the EM disturbances. This effect is highly correlated with the quality of the coupling between the supply network of the target (the victim) and the injection probe (the aggressor). Further, the underpowering seems to be more local than global: faults are located in the vicinity of the injection probe.

The novelty of this threat explains that no countermeasure dedicated to cope with EM injection has been yet proposed (to the best of our knowledge).

There is two other common fault injection means related to timing violations: clock and power supply glitches. They both have a global effect (i.e. the disturbance affects the whole chip). A delay-based countermeasure (CM) has been recently proposed and validated by [START_REF] Endo | An efficient countermeasure against fault sensitivity analysis using configurable delay blocks[END_REF] to cope with this kind of timing violation. However, the question of its efficiency against EM glitches was to be raised. Indeed, an EM disturbance located away from the actual implementation of the CM may induce a fault without triggering an alarm.

This paper reports an evaluation of a delay-based CM against EM glitches. Because a single CM was insufficient to detect with a high level of confidence EM induced faults, we have investigated the use of several CMs to attain this purpose. Conducting these experiments also provides many further insights into the properties of EM injection: how local its effect is and how the design of the injection probe may influence the process.

The contributions of this paper to that research field are:

• the disclosure of guidelines to implement a delay-based countermeasure against EM injection,

• a study and evidences of its local effect,

• an illustration of the key role of the injection probe,

• a further assessment of the actual threat related to EM injection. This article is organized as follows. Section II recalls some basics related to timing violations, describes the delay-based countermeasure, and reports its efficiency against clock and power glitches. Section III describes the experimental set-up and protocol, presents the experimental results and provides an interpretation. Finally section IV concludes the paper.

II. PRELIMINARIES

This section reminds the mechanisms involved in fault injection by timing constraint violation. It also describes the principle of the delay-based countermeasure we designed to cope with this injection technique. Its efficiency against clock and power supply glitches is reported likewise.

A. Timing constraint violation

Most digital ICs use one or several synchronous clock signals to synchronize their internal operations. Data are released from launch register banks on a clock rising edge, processed by the logic, and finally latched by capture registers on the next clock rising edge. Thus, in first approximation, the clock period (T clk ) has to be greater than the largest propagation delay of the logic (i.e. its critical time D pM ax ) to ensure correct operation. A precise writing of this timing constraint requires to take into account four other parameters: D clk2q , the time spent by a register to release a data after the clock rising edge; T skew , the clock skew between launch and capture registers; T jitter the clock jitter; and T setup , the setup time of the capture register. The setup time measures the amount of time a register's input data must be stable before the clock's rising edge to ensure reliable operation. This constraint expressed in Eq. 1 is called the setup time constraint (note that we do not describe here the clock pulse width constraint and the hold time constraint for the sake of brevity).

T clk > D clk2q + D pM ax + T setup + T skew + T jitter (1)
The time margin related to eq. 1 is called the time slack (i.e. the difference between required arrival time and actual arrival time at the input data of a register). A positive slack means that eq. 1 is fulfilled, conversely a negative slack is led to its violation.

The violation of this timing constraint is a straightforward means to inject faults into a circuit. A setup time violation arises if the last signal transition is too close to the clock rising edge (The DFF's output undergoes a metastable behavior [START_REF] Horstmann | Metastability behavior of cmos asic flip-flops in theory and test[END_REF]: it may stabilize either at a high or low state. An error may occur or not). Then, increasing the stress applied to the chip, an early latching arises. There is no signal transition during the setup time and an erroneous value is latched. Hereafter, we will refer to constraint timing violation for both cases.It exists many means to obtain a timing constraint violation for the purpose of injecting faults into an IC. The most common are clock and power supply glitches that induce a transient violation of eq. 1. A clock glitch [START_REF] Endo | An on-chip glitchy-clock generator for testing fault injection attacks[END_REF] consists in reducing temporarily the clock period (left handside of eq. 1) to obtain a negative slack, whereas a power glitch [START_REF] Zussa | From physical stresses to timing constraints violation[END_REF] induces a transient increase of the logic propagation times (right handside of eq. 1). The next sections report the principle of a countermeasure based on monitoring the advent of timing violations due to clock or power glitches, and some experimental validations.

B. The delay-based countermeasure

We have designed a delay-based countermeasure (CM) against timing violations from the implementation reported by Endo et al. in [START_REF] Endo | An efficient countermeasure against fault sensitivity analysis using configurable delay blocks[END_REF] (see also [START_REF] Selmane | Security evaluation of application-specific integrated circuits and field programmable gate arrays against setup time violation attacks[END_REF]). The principle of this CM consists in detecting the violation of a guarding delay prior to any timing violation. The clock signal is used as a reference to be able to draw comparisons between the guarding delay and the clock period (T clk ). In normal operation the guarding delay is set greater than the critical time (D pM ax ), but smaller than T clk . By doing so, if T clk is decreased for the purpose of inducing a timing violation, it will have to be shorter than the guarding delay (what we call a guarding delay violation) before inducing a timing violation. Hence, if the CM is able to detect this guarding delay violation, it will also be able to detect any fault injection attempt by clock glitches. Similarly, if the power supply voltage is decreased for the purpose of injecting a fault by increasing the target's critical time, the guarding delay which depends upon the voltage supply in a similar way will be increased as well. As a consequence the violation of the guarding delay will arise and be detected first. The schematic of the CM depicted in Fig. 1-a fits with the above mechanisms. The guarding delay (denoted delay) is implemented with the circuit logic. It is used to obtain a delayed clock (denoted DCK) from the clock signal (denoted CK) where DCK(t) = CK(tdelay). A D flip-flop (DFF) is used as a phase comparator between CK (connected to its data input) and DCK (connected to its clock input). The design is tuned in order to comply with the timing given in Eq. 2.

D pM ax < delay < T clk (2)
In normal operation, as depicted in Fig. 1-b, the DFF output (denoted alarm) is low. Fig. 1-c illustrates the detection of a power supply glitch. As the power supply voltage is decreased the guarding delay is increased, it goes larger than T clk . Thus a high level is latched by the DFF on the next rising edge of DCK: the alarm is triggered at high level indicating a fault injection attempt. For the sake of brevity we have not sketched the detection of a clock glitch which is very similar. This countermeasure was first designed to be embedded in an FPGA running the AES encryption algorithm as described in the next subsection. We also report in subsection II-D experiments of its validity against clock and power glitches. The next section devoted to the study of the effect of EM glitches describes its efficiency against such perturbations. It will be exposed that this CM may be used likewise as a test pattern to study the local effects of EM glitches.

C. Implementation

The experiments reported in this work were carried out on an FPGA (Xilinx Spartan 700). At first, the CM was embedded along with a hardware implementation of the AES algorithm (its 128-bits version). The device was operated at 100MHz at its nominal power supply (1.2V).

The CM primitive was implemented as a hard macro [START_REF]Xilinx. FPGA Editor Guide[END_REF]. It allows to change its location within the logic slices of the FPGA without modifying its actual architecture after a proper tuning of the guarding delay has been found. Furthermore, it ensured that its design was identical after each of the synthesize and routing steps we have ran. The Xilinx software 'FPGA Editor' enables the creation of hard macros, which are manually placed, routed and configured designs, which can be instantiated multiple times in an FPGA design. This feature can be exploited to design an exactly defined delay based alarm and instantiate it several times.

In a second time, while studying more precisely the effect of EM glitches, several CM primitives were embedded as described in subsection III-A. Thanks to the use of hard macros, every instance of the CM was identical (however because of within-die variations [START_REF] Sedcole | Within-die delay variability in 90nm fpgas and beyond[END_REF] the various CMs may have different thresholds of triggering).

D. Clock and power supply glitches

As shown in [START_REF] Zussa | From physical stresses to timing constraints violation[END_REF] clock and power supply glitches lead to timing constraint violations that both affect the whole die (i.e. they have a global effect). Thus, a single countermeasure should be sufficient to detect any fault injection attempt related to these techniques. For the purpose of drawing comparisons while studying the effects of EM glitches, we have first carried out injection experiments by means of clock and power glitches. The experimental setup is described in Fig. 6.

Clock glitch injections (several thousands with random data) were conducted during the 9 th round of the AES by shortening progressively the corresponding clock period from 10ns to 7.2ns. As expected the alarm was always triggered prior to the injection of the first fault. Fig. 2 reports an example of the obtained results. The occurrence rates of the triggering of the alarm (in red) and of the injection of a fault (in blue) are drawn versus the clock period of the attacked round. A one hundred percent of detection was recorded before that the fault occurrence rate started to rise.

Similar experiments were made by means of negative power supply glitches of 10ns centered on the AES 9 th round. Alike, as shown in Fig. 3, were the occurrence rates are drawn versus the amplitude of the power glitches, the triggering of the alarm was always prior to any fault injection.

As a result, it appeared that a single glitch detector is sufficient to protect the whole FPGA.

III. ELECTROMAGNETIC GLITCHES

This new fault injection means takes advantage of the presence of the numerous horizontal and vertical loops, in the routing of the power and ground networks, to inject sudden variations of the current flowing through. These sudden current variations produce voltage drops and ground bounces that alter the propagation of signals through the logic and therefore induce timing faults. Thus, they produce an effect similar to that of voltage glitches. However, because of the locality of the EM coupling, EM spikes can be locally injected, anywhere, in the IC resulting in an injection that could be sufficiently local to get round a single glitch detector usually located in the neighborhood of supply pads or of the on-chip voltage regulator.

The following experiments have several goals. The first one is to characterize the spatial limitation of the aforementioned CM. This also demonstrates the locality of EM glitches, our second goal. Finally, the last goal of these experiments is to demonstrate that handmade probes designed according to some of the recommendations of [START_REF] Omarouayache | Magnetic microprobe design for em fault attackmagnetic microprobe design for em fault attack[END_REF] are able to concentrate intense and powerful magnetic fields on reduced areas.

A. The FPGA implementation EM injection experiments were conducted on an FPGA implementing the AES algorithm and five glitch detectors. Fig. 4-b shows the floorplan, the placement and the routing of the AES and glitch detectors mapped into the FPGA. The location of the detectors is highlighted with red squares. Fig. 4-a presents a front view of the opened device after chemical decapsulation, it allows determining the size and location of the silicon die (underlined in black). This choice was done in order to try to detect any EM glitch penetrating in the IC at any point of the IC surface. At this point, this choice of five detectors as well as their position was done from purely geometric considerations.

B. Electromagnetic Glitch

The EM injection probes were handmade probes. They are made up of a ferrite core end around which a copper wire is enclosed to form 4-5 turns (see Fig. 5). During the experiments two different probes were used. One, denoted 'G', with a diameter of 3000µm and a flat tip end; the second, denoted 'F', with a ferrite core of 500µm with a sharp tip end of diameter of 300µm according to recommendations of [START_REF] Omarouayache | Magnetic microprobe design for em fault attackmagnetic microprobe design for em fault attack[END_REF]. To produce EM glitches a voltage pulse generator (4A and 20V-200V pulse capability) was adopted. It delivers pulses with a duration ranging from 10ns to 100ns to the EM probes. To access the EM glitch susceptibility of the IC as well as the efficiency of the glitch detectors, a motorized stage was used to scan the IC surface and its neighborhood using the EM glitch Fig. 5. The two EM injection probes platform (see Fig. 6). The scanned area was 14.4 × 14.4mm 2 while the silicon die area (white square in cartographies) was 5.6×5mm 2 . A same displacement step of 400µm was chosen along the XY axis. At each location of the injection probe, 20 EM glitches were induced. As a result it takes about 220min. to perform a full cartography.

Using this platform, several EM glitch cartographies (targeting the AES 9 th round) were launched with both probes G and F and pulses of different amplitudes were used to modify both the power and the resolution of the injections. During these cartographies, the output of the glitch detectors (alarms) as well as the ciphertexts were monitored.

C. Spatial limitation of glitch detectors

Fig. 7 displays the cartography showing the triggering rate of the alarm located in the bottom right part of the floorplan while Fig. 8 gives a cartography of the fault injection rate. Both were obtained with probe G and 200V pulse amplitude. 1296 probe locations were investigated to draw the cartography, corresponding to 25360 EM glitches. The alarm was triggered 7308 times. From Fig. 7 one may conclude that one glitch detector is not sufficient to detect all EM glitches produced above the die (the white rectangle) and its close vicinity. From Fig. 8 one may conclude that this alarm, which was detecting absolutely all voltage or clock glitches, is not able to detect all EM glitches producing faults. There is, indeed, several locations in the right bottom part of the investigated area that have one hundred fault injection rates while having a zero detection rate. From the 1935 faults actually injected, 80 were undetected. Moreover, this demonstrates that EM glitch has a local effect contrarily to voltage or clock glitches which are global.

D. Enhancing the robustness of the IC against EM glitches

The above results show that our detector has a wide, but insufficient, area of detection. Indeed all EM glitches attempts are not detected. It is therefore mandatory to integrate several detectors. This explains why five detectors have been integrated. Fig. 9 shows the detection rate of EM glitches obtained by considering the answers of the five detectors rather than the answer of a single one. The detection area was clearly

E. Number of detectors

Considering the above results, one may wonder about the number of glitch detectors a designer has to integrate to secure an IC. If the answer depends of course on the quality of the detectors, it also depends on the resolution of the EM injection probes that the adversaries may use and of the pulse amplitude Indeed Fig. 10 shows the triggering rate of at least one alarm among the five alarms with probe G and F and pulse amplitudes equal to 200V, 100V or 75V while Fig. 11 gives the fault injection rates. Note that these rates are drawn according the same color scale used in previous cartographies.

25530 EM glitches were injected to draw each cartography. With probe G at 200V pulse amplitude, the alarm triggering rate (at least one alarm out of five triggered) was 32%, 1995 faults were induced which were always detected. At 100V amplitude, the alarm triggering rate was 21%, 1052 faults were injected. 8 of these injected faults were undetected. With probe F at 100V pulse amplitude, the alarm triggering rate was 6% and 133 faults were induced from which 12 were not detected. At 75V pulse amplitude with probe F, the alarm triggering rate was 5.1% and 115 faults were induced. Five faults were undetected. The decrease of the alarm triggering rate as the pulse amplitude or the probe size were decreased was expected because a smaller voltage pulse amplitude induces a weaker EM disturbance. However, the design of probe F makes it possible to defeat the detection capability of the five glitch detectors simultaneously with a success rate slightly above 10%.

Therefore the use of five glitch detectors may not be sufficient or their placement not optimal depending on the resolution of the injection probe. However, despite succeeding in injecting a small number of undetected faults, the feasibility of the attack is questionable. Indeed, many fault injection attempts would have been detected first and an adapted answer to that threat should have been taken [START_REF] Endo | An on-chip glitchy-clock generator for testing fault injection attacks[END_REF].

One may also observe that the detection area is slightly narrower when reducing the pulse amplitude and much more narrow when the resolution of the EM probe is enhanced. This latter experimental results confirms that the use of ferrite rods with a sharp end allows to concentrate the magnetic field on reduced surface (as expected from the simulation results of [START_REF] Omarouayache | Magnetic microprobe design for em fault attackmagnetic microprobe design for em fault attack[END_REF]) and therefore:

• to enhance the spatial resolution of EM injection probes, • to get around efficient glitch detectors designed to detect clock and power spikes.

IV. CONCLUSION

In this paper, EM glitch was compared to voltage and clock glitch fault injection techniques. If clock and voltage glitches can be thwarted by the use of a single delay-based glitch detector, we experimentally highlighted that the use of a single detector may not be sufficient to thwart EM glitch based fault attacks. Additionally, we had experimentally demonstrated that the use of ferrite rods with sharp end allows enhancing the spatial resolution of the EM glitches. Finally, the efficiency of glitch detector (based on the monitoring of delay at least) in detecting EM pulses has been pointed out. However, the number of detectors to be used as well as their optimal placement on the IC surface remain an open problems.

Fig. 1 .

 1 Fig. 1. Delay-based countermeasure principle

Fig. 2 .Fig. 3 .

 23 Fig. 2. Detection of clock glitches.

Fig. 4 .

 4 Fig. 4. (a) Test chip after front side decapsulation and (b) AES (blue) and CMs (red squares) implementation within the design floor plan

Fig. 6 .Fig. 7 .

 67 Fig. 6. Experimental setup

Fig. 8 .

 8 Fig. 8. Fault occurrence rate w.r.t. the probe tip end position (Probe G, pulse amplitude equal to 200V).

Fig. 9 .

 9 Fig. 9. Cartography of the alarm triggering rate -5 detectors (Probe G, pulse amplitude equal to 200V).

Fig. 10 .

 10 Fig. 10. Rate of triggering of at least one alarm among the five alarms with Probes G (flat tip end 3000µm) and F (sharp tip end 300µm) with pulse amplitude equal to 200V, 100V or 75V.

Fig. 11 .

 11 Fig. 11. Fault injection rate with Probes G (flat tip end 3000µm) and F (sharp tip end 300µm) with pulse amplitude equal to 200V, 100V or 75V.

TABLE I OVERHEAD

 I

	OF THE GLITCH DETECTORS
		Number	Overhead
		of used slices
	None Protected AES	2741
	AES + 1 Glitch detector	2750	0.3 %
	AES + 5 Glitch detectors	2785	1.6 %
	extended. Table I reports the overhead in terms of the area used
	by the glitch detectors.

ACKNOWLEDGMENT

We kindly acknowledge the french Agence Nationale de la Recherche for funding this research work with the E-MATA HARI project (ANR-12-INSE-0005) and the EMAISeCi project (ANR-10-SEGI-0012).