Thierry Gallay

Denis Serre

The numerical measure of a complex matrix

We introduce a natural probability measure over the numerical range of a complex matrix A ∈ M n (C). This numerical measure µ A can be defined as the law of the random variable AX, X ∈ C when the vector X ∈ C n is uniformly distributed on the unit sphere. If the matrix A is normal, we show that µ A has a piecewise polynomial density f A , which can be identified with a multivariate B-spline. In the general (nonnormal) case, we relate the Radon transform of µ A to the spectrum of a family of Hermitian matrices, and we deduce an explicit representation formula for the numerical density which is appropriate for theoretical and computational purposes. As an application, we show that the density f A is polynomial in some regions of the complex plane which can be characterized geometrically, and we recover some known results about lacunas of symmetric hyperbolic systems in 2 + 1 dimensions. Finally, we prove under general assumptions that the numerical measure of a matrix A ∈ M n (C) concentrates to a Dirac mass as the size n goes to infinity.

Introduction

If A ∈ M n (C) is a complex square matrix of size n ∈ N * , the numerical range of A is the compact subset of the complex plane defined by

W (A) = Ax, x ∈ C x ∈ C n , x = 1 ,
where x, y = y * x is the usual scalar product in C n and x = x, x 1/2 . It is quite obvious that W (A) ⊃ σ(A), where σ(A) (the spectrum of A) is the collection of all eigenvalues of A, and that W (A) = W (U * AU ) for any unitary matrix U ∈ U n (C). Moreover, a celebrated result due to Toeplitz [START_REF] Toeplitz | Das algebraische Analogon zu einem Satze von Fejér[END_REF] and Hausdorff [START_REF] Hausdorff | Der Wertevorrat einer Bilinearform[END_REF] asserts that W (A) is always a convex subset of the complex plane. In particular, W (A) contains the convex hull of σ(A), and it is easy to verify that W (A) = conv(σ(A)) if the matrix A is normal, namely AA * = A * A. The interested reader is referred to Chapter 1 of [START_REF] Horn | Topics in matrix analysis[END_REF] for a detailed discussion of the various properties of the numerical range, including complete proofs.

Let ∂B n = {x ∈ C n | x = 1} be the unit sphere in C n , considered as a real manifold of dimension 2n -1. By definition, the numerical range W (A) is the image of the numerical map Φ A : ∂B n → C defined by Φ A (x) = Ax, x , x ∈ ∂B n .

The algebraic and geometric properties of the map Φ A have been extensively studied, see [START_REF] Kippenhahn | Über den Wertevorrat einer Matrix[END_REF][START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF][START_REF] Bazer | Lacunas of the Riemann matrix of symmetric-hyperbolic systems in two space variables[END_REF][START_REF] Fiedler | Geometry of the numerical range of matrices[END_REF][START_REF] Jonckheere | Differential topology of numerical range[END_REF][START_REF] Joswig | On the numerical range map[END_REF]. In particular, the set of all critical values of Φ A , which we denote by Σ A ⊂ C, has received 1 much attention, because this is an interesting object which contains a lot of information on the matrix A. For instance, it is known that ∂W (A) ⊂ Σ A and W (A) = conv(Σ A ). As is shown in [START_REF] Joswig | On the numerical range map[END_REF], there exists a real algebraic curve C A ⊂ C ≃ R 2 with the property that Σ A = C A ∪ C ′ A , where C ′ A denotes the set of all line segments joining pairs of points of C A at which C A has the same tangent line [START_REF] Joswig | On the numerical range map[END_REF]. Under generic assumptions on A, the bitangent set C ′ A is empty, and the critical set Σ A is therefore the union of a finite number of closed curves, one of which is the boundary of the numerical range W (A). This distinguished curve is smooth, and encloses all the other ones in its interior. We refer to Section 5 below for more details on the geometry of the singular set, and to Section 7 for a few concrete examples.

Our purpose in this paper is to introduce another mathematical quantity which is naturally related to the numerical map Φ A . Given A ∈ M n (C), the numerical measure of A is the probability measure µ A on C defined by the formula

C φ(z) dµ A (z) = ∂B n φ( Ax, x ) dσ(x) , (1)
for all continuous functions φ : C → C. Here σ denotes the Euclidian measure on the unit sphere ∂B n , normalized as a probability measure. In words, the numerical measure is thus the image under the numerical map of the normalized Euclidean measure on the unit sphere. Equivalently, if X is a random variable that is uniformly distributed on ∂B n , the numerical measure µ A is just the distribution of the random variable AX, X ∈ C. This probabilistic interpretation will be useful later, especially in Section 8. Our first goal is to establish a few general properties of the numerical measure µ A . It is clear by construction that µ A is invariant under unitary conjugations of A, namely µ U * AU = µ A for all U ∈ U n (C). This is precisely the reason why we used the Euclidean measure on ∂B n in the definition [START_REF] Atiyah | Hyperbolic differential equations and algebraic geometry (after Petrowsky)[END_REF]. It is also easy to verify that the support of µ A is exactly the numerical range W (A), see Section 2 below. Less obvious, perhaps, is the fact that µ A is absolutely continuous with respect to the Lebesgue measure λ on W (A), so that we can define the numerical density f A as the Radon-Nikodym derivative of µ A with respect to λ (in the particular situation where W (A) reduces to a line segment Γ, we understand λ as the one-dimensional Lebesgue measure on Γ, see Section 2.) We also prove that the numerical density f A is strictly positive in the interior of A, a property that can be interpreted as a strong version of Hausdorff's theorem [START_REF] Hausdorff | Der Wertevorrat einer Bilinearform[END_REF]. Finally, we shall see that the singular support of µ A is contained in the critical set Σ A , which means that the numerical density f A is smooth outside Σ A . In fact, we conjecture that sing supp(µ A ) = Σ A for all A ∈ M n (C), but this has not been proved yet.

After these general properties have been established, our next goal is to give a more precise description of the numerical density f A . For this purpose, it is convenient to distinguish between various cases:

1. (The scalar case) If W (A) is reduced to a single point {z}, then A = zI n (where I n ∈ M n (C) denotes the identity matrix) and µ A = δ z . In this trivial situation, of course, there is no need to introduce a numerical density.

2. (The Hermitian case) Assume that n ≥ 2 and that W (A) ⊂ C is a line segment. Then we can find z ∈ C, θ ∈ [0, 2π], and a Hermitian matrix H such that A = zI n + e iθ H. If λ 1 ≤ λ 2 ≤ • • • ≤ λ n are the eigenvalues of H, we shall see in Section 3 that the numerical measure µ H is absolutely continuous with respect to Lebesgue's measure on R, and that the corresponding density f H is exactly the normalized B-spline of degree n -2 with knots λ 1 , . . . , λ n [START_REF] De Boor | A Practical Guide to Splines[END_REF]. In particular, f H is polynomial of degree n -2 on each interval [λ i , λ i+1 ], vanishes identically outside [λ 1 , λ n ], and is continuous at each point λ i together with its derivatives up to order d i = n -2m i , where m i ≥ 1 is the multiplicity of λ i as an eigenvalue of H (if d i < 0, then f H is discontinuous at λ i .) This gives an explicit representation of the numerical measure µ H , and the measure µ A is the image of µ H under the affine isometry w → z + e iθ w.

3. (The normal case) Suppose now that n ≥ 3 and that A ∈ M n (C) is a normal matrix whose spectrum σ(A) is not contained in a line segment. Then W (A) = conv(σ(A)) is a convex polygon with nonempty interior, and it turns out that the numerical density f A is the bivariate B-spline of degree n -3 whose knots are the eigenvalues of A. Here we refer to the work of W. Dahmen [START_REF] Dahmen | On multivariate B-splines[END_REF] for the definition and the main properties of multivariate B-splines. In this particular case, the critical set Σ A is thus the collection of all line segments joining pairs of eigenvalues of A, and the density f A is polynomial of degree n -3 in each connected component of C \ Σ A . In the generic situation where no straight line contains more than two eigenvalues of A, one can show that f A is continuous together with its derivatives up to order n -4 (and is discontinuous on ∂W (A) if n = 3.) 4. (The nonnormal case) Finally, we consider the most interesting situation where the matrix A ∈ M n (C) is not normal. In that case, there is no explicit formula for the numerical density, but the problem can be reduced in some sense to the Hermitian case by the following simple observation. For any θ ∈ S 1 = R/(2πZ), let H(θ) be the Hermitian matrix defined by

H(θ) = 1 2 (e -iθ A + e iθ A * ) = A 1 cos(θ) + A 2 sin(θ) , (2)
where A 1 = (A + A * )/2 and A 2 = (A -A * )/(2i). Then Re (e -iθ Ax, x ) = H(θ)x, x for all x ∈ ∂B n . Now, if the random variable X is uniformly distributed on ∂B n , the distribution of H(θ)X, X is by definition the numerical measure µ H(θ) , whereas the distribution of Re (e -iθ AX, X ) is easily identified as the two-dimensional Radon transform of the numerical measure µ A , with angular parameter θ ∈ [0, 2π]. We thus have

Rµ A (θ) = µ H(θ) , θ ∈ S 1 , (3)
where R denotes the two-dimensional Radon transformation. Since the numerical density of H(θ) is known to be the B-spline based on the eigenvalues λ 1 (θ), . . . , λ n (θ) of H(θ), we can reconstruct the numerical measure µ A by inverting the Radon transformation in (3), using the well-known backprojection method which plays an important role in tomography [START_REF] Helgason | The Radon transform, 2nd edition[END_REF]. This provides a useful representation formula for the numerical density, as well as an efficient algorithm for numerical calculations, see Section 4 for more details.

It is worth mentioning here that the critical set Σ A can be conveniently characterized using the family of Hermitian matrices H(θ) associated with A. Indeed, if we define the eigenvalues λ j (θ) in such a way that they depend analytically on θ, one can shown that the algebraic curve C A which generates Σ A is given by

C A = e iθ (λ j (θ) + iλ ′ j (θ)) j ∈ {1, . . . , n} , θ ∈ [0, π] ,
see [START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF][START_REF] Jonckheere | Differential topology of numerical range[END_REF][START_REF] Joswig | On the numerical range map[END_REF] and Section 5 below. In the generic case where λ 1 (θ) < λ 2 (θ) < • • • < λ n (θ) for all θ ∈ [0, π], the bitangent set C ′ A is empty and Σ A = C A . As was already mentioned, the numerical density f A is smooth (in fact, real-analytic) on each connected component of C \ Σ A . The regularity across Σ A is more difficult to study, but we shall show in Section 6.2 that f A is everywhere of class C n-3 if n ≥ 3 and A ∈ M n (C) satisfies some generic hypotheses, which exclude in particular the case of normal matrices. In addition, for an arbitrary matrix of size n, we shall prove that all derivatives of f A of order n -2 vanish identically in some The numerical density f A is represented for a typical matrix A ∈ M 3 (R), given by (55) below. In the contour plot (left), the exterior ovate curve is the boundary of the numerical range W (A), and the other component of the critical set Σ A is the interior cuspidal triangle. The three-dimensional plot (right) confirms that the numerical density is continuous, positive inside the numerical range, and constant over the cuspidal triangle, in agreement with the results of Section 6.

distinguished regions, which have the following geometric characterization. For any z ∈ C \ Σ A , let N (z) be the number of straight lines containing z which are tangent to the curve C A , see Eq. ( 39) below for a precise definition where possible multiplicities are taken into account. It is easy to verify that N (z) is constant in each connected component of C \ Σ A , and that N (z) ≤ n [START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF]. The distinguished regions where f A is polynomial of degree n -3 (if n ≥ 3) or f A ≡ 0 (if n = 2) are exactly those connected components of C \ Σ A on which N (z) takes its maximal value n. This remarkable property of the numerical density, which is one of our main results, will be established in Section 6.1. The geometric condition N (z) = n is always satisfied in the complement of the numerical range, where f A vanishes identically, but for many matrices of size n ≥ 3 is it also met in some regions inside W (A). For instance, in the three-dimensional example represented in Fig. 1, it is easy to verify that N (z) = 3 if z is outside W (A) or inside the cuspidal triangle, and N (z) = 1 in the intermediate region where the numerical density is not constant.

At this point, it is necessary to make a connection with the theory of lacunas of symmetric hyperbolic systems of partial differential equations in 2 + 1 variables. Given A ∈ M n (C), we consider the following system of linear PDE's in R t × R 2

x :

∂ t u + A 1 ∂ x 1 u + A 2 ∂ x 2 u = 0 , (4)
where A 1 , A 2 are as in ( 2) and u = (u 1 , . . . , u n ) ⊤ : R × R 2 → R n . The fundamental solution of ( 4) is the unique (matrix-valued) distribution E supported in the half-space R + × R 2 which satisfies

∂ t E + A 1 ∂ x 1 E + A 2 ∂ x 2 E = I n δ t=0 ⊗ δ x=0 . (5)
One can show that E(t, x) is homogeneous of degree -2 in t and x. It is thus sufficient to consider the time-one trace E * = E(1, •), which is a distribution on R 2 . Due to the finite speed of propagation, it is well-known that E * is zero outside a compact set of R 2 , but it may also happen that E * vanishes identically in some regions inside the domain of influence of the origin. Such regions are called lacunas of the hyperbolic system (4). The properties of the fundamental solution of symmetric hyperbolic systems have been studied by many authors, see e.g. [START_REF] Petrowsky | On the diffusion of waves and lacunas for hyperbolic equations[END_REF][START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF][START_REF] Bazer | Lacunas of the Riemann matrix of symmetric-hyperbolic systems in two space variables[END_REF][START_REF] Atiyah | Lacunas for hyperbolic differential operators with constant coefficients[END_REF][START_REF] Atiyah | Lacunas for hyperbolic differential operators with constant coefficients[END_REF]. In the particular case of system (4), J. Bazer and D. Yen have shown that, if one identifies C with R 2 , the singular support of the distribution E * is contained in the critical set Σ A , and the (stable) lacunas of system (4) are exactly the regions described above where the numerical density f A is polynomial of degree n -3. This remarkable coincidence is of course not fortuitous. In Section 6.3, we explain it by showing that the fundamental solution E can be expressed as a linear combination of derivatives of order n -1 of a homogeneous extension of the numerical density f A . This connection allows us to recover some of the main results of [START_REF] Bazer | Lacunas of the Riemann matrix of symmetric-hyperbolic systems in two space variables[END_REF], and therefore confirms that the numerical measure is a natural quantity attached to the matrix A. One might even argue that µ A contains more information than E * , since for instance supp(µ A ) = W (A) while supp(E * ) is in general strictly smaller and not necessarily convex, see Section 7. Similarly, we believe that sing supp(µ A ) always coincide with Σ A , while sing supp(E * ) is usually smaller.

A final question that is worth investigating is the behavior of the numerical measure µ A when the size of the matrix A goes to infinity. Here of course, specific assumptions have to be made in order to obtain convergence results. Suppose for instance that {A n } n≥1 is a sequence of complex matrices with A n ∈ M n (C), Tr(A n ) = 0, and A n ≤ C for all n ≥ 1. If, for each n ≥ 1, X n is a random variable that is uniformly distributed on ∂B n , we show in Section 8 that the complex variable A n X n , X n converges almost surely to zero as n → ∞. This is reminiscent of the strong law of large numbers in probability theory. Under slightly stronger assumptions, we also establish the analog of the central limit theorem in this context. Our convergence results mean that µ A is very close to δ z when dim(A) is large, where z is the barycenter of σ(A). This explains why plotting Ax, x for randomly chosen points x ∈ ∂B n is a very unefficient algorithm for determining the numerical range

W (A) if A is a large matrix!
The rest of the paper is organized as follows. In Section 2, we establish some general properties of the numerical measure. Section 3 is devoted to the particular situations where the matrix A is Hermitian or normal. The nonnormal case is treated in Sections 4-6, which constitute the core of the paper. In Section 4, we derive a representation formula for the numerical density using the inversion of the Radon transformation. Section 5 collects a few results on the geometry of the critical set Σ A , which are mainly borrowed from [START_REF] Kippenhahn | Über den Wertevorrat einer Matrix[END_REF][START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF][START_REF] Joswig | On the numerical range map[END_REF]. These informations are used in Section 6 to derive an explicit formula for the derivatives of order n -2 of the numerical density, which allows us to obtain generic regularity results and to express the fundamental solution of the hyperbolic system (4) in terms of derivatives of the numerical density. To illustrate our results, a few explicit examples are treated in Section 7. Finally, we investigate in Section 8 the concentration properties of the numerical density for large matrices, and we discuss in Section 9 a possible extension of our results to hyperbolic polynomials with an arbitrary number of variables.

Acknowledgements. This work has benefited of stimulating discussions with several of our colleagues, including Y. Colin de Verdière, F. Faure, and A. Joye.

General properties of the numerical measure

In this section, we establish a few general properties of the numerical measure of a complex matrix. In particular, we show that µ A is absolutely continuous with respect to the Lebesgue measure on W (A), and we prove a direct sum formula which will be useful later.

Support and regularity properties

We first show that the support of the numerical measure always coincides with the numerical range of the matrix.

Lemma 2.1 For any

A ∈ M n (C), one has supp(µ A ) = W (A). Proof. If V = C \ W (A), then Φ -1 A (V ) = ∅, hence µ A (V ) = σ(Φ -1 A (V )) = 0. This shows that supp(µ A ) ⊂ W (A). Conversely, if V ⊂ C is any open set such that V ∩ W (A) = ∅, then Φ -1 A (V ) is a nonempty open subset of ∂B n , hence µ A (V ) = σ(Φ -1 A (V )) > 0. Thus W (A) ⊂ supp(µ A ). Our next goal is to locate the singular support of µ A . We recall that x ∈ ∂B n is a regular point of Φ A if the differential map d x Φ A : T * x ∂B n → C is onto.
Otherwise, we say that x is a critical point. The following characterization will be useful: Lemma 2.2 [START_REF] Jonckheere | Differential topology of numerical range[END_REF][START_REF] Joswig | On the numerical range map[END_REF] Let A ∈ M n (C). 1) A point x ∈ ∂B n is a critical point of the numerical map Φ A if and only if x is an eigenvector of the Hermitian matrix H(θ) defined in (2) for some θ ∈ [0, π].

2) The differential of Φ A vanishes at x ∈ ∂B n if and only if x is an eigenvector of both A and A * .

In other words, the range of the differential d x Φ A has (real) dimension 1 if and only if x is an eigenvector of H(θ) for a unique θ ∈ [0, π), and is reduced to {0} if and only if x is an eigenvector of H(θ) for all θ ∈ [0, π]. The proof is neither new nor difficult, but we shall repeat it here in order to introduce some notation that will be needed later on.

Proof. Since Φ A (e iθ x) = Φ A (x) for all θ ∈ [0, 2π], we can consider the numerical map as acting on the quotient space ∂B n /S 1 ≃ CP n-1 [START_REF] Jonckheere | Differential topology of numerical range[END_REF]. Thus, to detect the critical points of Φ A , we study the reduced map ΦA :

∂B n /S 1 → C defined by ΦA ([x]) = Φ A (x) = A 1 x, x + i A 2 x, x , x ∈ ∂B n ,
where [x] = {e iθ x | θ ∈ S 1 } and A 1 , A 2 are the Hermitian matrices introduced in (2). If x ∈ ∂B n , the tangent space to ∂B n /S 1 at [x] is just the (2n -2)-dimensional affine subspace {[x + y] | y ∈ C n , x, y = 0}. Thus, using the definition above of ΦA , it is straightforward to verify that, for all y ∈ C n with x, y = 0, one has

1 2 d [x] ΦA (y) = (A 1 x|y) + i(A 2 x|y) = (v 1 (x)|y) + i(v 2 (x)|y) , (6)
where (x|y) = Re x, y = (Re x) ⊤ (Re y) + (Im x) ⊤ (Im y) denotes the real scalar product in C n ≃ R 2n , and

v 1 (x) = A 1 x -A 1 x, x x , v 2 (x) = A 2 x -A 2 x, x x .
Of course, replacing A 1 x, A 2 x with v 1 (x), v 2 (x) has no effect in (6) since x, y = 0, but after this substitution we can let y run over the whole of C n without increasing the range. So our task is reduced to computing the rank of the R-linear map y → (v 1 (x)|y) + i(v 2 (x)|y), which is just the rank of the 2 × 2 matrix

D(x) = (v 1 (x)|v 1 (x)) (v 1 (x)|v 2 (x)) (v 2 (x)|v 1 (x)) (v 2 (x)|v 2 (x)) . (7)
By the Cauchy-Schwarz inequality, the positive matrix D(x) is singular if and only if there exists θ ∈ [0, π] such that v 1 (x) cos θ + v 2 (x) sin θ = 0, which exactly means that x is an eigenvector of H(θ). Moreover, D(x) = 0 if and only if v 1 (x) = v 2 (x) = 0, which is equivalent to saying that x is an eigenvector of both A 1 and A 2 , hence of both A and A * .

Let Σ A ⊂ C denote the set of all critical values of Φ A , namely Σ A = Φ A (Γ(A)) where Γ(A) ⊂ ∂B n is the set of all critical points of Φ A . Our next result is:

Lemma 2.3 If A ∈ M n (C), then sing supp(µ A ) ⊂ Σ A .
Proof. If the numerical range W (A) is reduced to a line segment or to a single point, then Σ A = W (A), hence sing supp(µ A ) ⊂ supp(µ A ) = Σ A by Lemma 2.1. Thus, we assume from now on that W (A) has nonempty interior. By Sard's lemma, the critical set Σ A is then a compact subset of W (A) with zero Lebesgue measure. We have to show that there exists a smooth density function

f A ≥ 0 such that dµ A (z) = f A (z) dz on C \ Σ A . Clearly, we must have f A = 0 on C \ W (A). If z ∈ W (A) \ Σ A , then N z := Φ -1 A (z)
is a compact submanifold of ∂B n of codimension 2, which depends smoothly on z. Using classical arguments, involving a partition of unity and the Implicit Function Theorem, it is not difficult to verify that, for any continuous function φ with supp(φ) ⊂ W (A) \ Σ A , one has

∂B n φ( Ax, x ) dσ(x) = 1 ω n C φ(z) Nz dν(x) 2∆(x) 1/2 dz ,
where ω n = 2π n /((n-1)!) is the total measure of ∂B n , ν is the (2n-3)-dimensional Euclidean measure on the submanifold N z , and ∆(x) = det D(x) where D(x) is the 2 × 2 matrix defined in [START_REF] De Boor | A Practical Guide to Splines[END_REF]. Remark that 2∆(x) 1/2 = λ1 (x) λ2 (x), where λ1 (x), λ2 (x) are the singular values of the differential map d x Φ A . In view of (1), we conclude that dµ A (z) = f A (z) dz on C \ Σ A , where

f A (z) = 1 ω n Nz dν(x) 2∆(x) 1/2 , z ∈ C \ Σ A . (8)
It is easily verified that the density f A is smooth and strictly positive on W (A) \ Σ A .

The results obtained so far are summarized in the following proposition, which also asserts that the numerical measure is absolutely continuous with respect to Lebesgue's measure on W (A).

Proposition 2.4

Let A ∈ M n (C). 1) If the numerical range W (A) has nonempty interior, the numerical measure µ A is absolutely continuous with respect to the (two-dimensional) Lebesgue measure on C ≃ R 2 . The numerical density f A = dµ A /dz is smooth outside the critical set Σ A .

2) If A is a nonscalar Hermitian matrix, then W (A) ⊂ R and the numerical measure is absolutely continuous with respect to the (one-dimensional) Lebesgue measure on R. The numerical density f A = dµ A /dx is smooth outside the spectrum σ(A). Remark 2.5 As is explained in the introduction, Proposition 2.4 covers all interesting cases. Indeed, if the numerical range has empty interior, then either W (A) is reduced to a single point, in which case A is a scalar matrix and µ A is just a Dirac mass, or W (A) is a line segment of nonzero length, in which case A can be reduced to a nonscalar Hermitian matrix by a simple affine transformation.

Proof. Using the same notations as in Lemmas 2.2 and 2.3, we observe that Γ(A) = {x ∈ ∂B n | ∆(x) = 0}, where ∆(x) = det D(x) is a polynomial in the 2n variables Re x i , Im x i (i = 1, . . . , n). Thus one of the following two situations must occur: 1) Γ(A) is an algebraic submanifold of ∂B n of codimension at least 1. By Sard's lemma, this is the case if and only if W (A) has nonempty interior. In that situation, since we already know that µ A has a smooth density outside the critical set Σ A , we only need to show that µ A (Σ A ) = 0. Given ǫ > 0, let Γ ǫ (A) = {x ∈ ∂B n | dist(x, Γ(A)) ≤ ǫ}, where "dist" denotes here the geodesic distance on the unit sphere. We decompose

Φ -1 A (Σ A ) = Φ -1 A (Σ A ) ∩ Γ ǫ (A) ∪ Φ -1 A (Σ A ) ∩ Γ ǫ (A) c = E 1 (ǫ) ∪ E 2 (ǫ) , where Γ ǫ (A) c = ∂B n \ Γ ǫ (A). Since σ(Γ(A)) = 0, we have σ(E 1 (ǫ)) ≤ σ(Γ ǫ (A)) → 0 as ǫ → 0.
Moreover, the proof of Lemma 2.3 shows E 2 (ǫ) is a codimension two submanifold of ∂B n , so that σ(E 2 (ǫ)) = 0 for any ǫ > 0. Using the definition of the numerical measure, we conclude that

µ A (Σ A ) = σ(Φ -1 A (Σ A )) = 0. 2) Γ(A) = ∂B n .
This is the case if and only if W (A) has empty interior, and without loss of generality we can then assume that the matrix A is Hermitian. Since W (A) ⊂ R, it is more natural here to consider Φ A as a map from ∂B n into R. If we do that, then repeating the proofs of Lemmas 2.2 and 2.3 we easily find that the critical points of Φ A are exactly the eigenvectors of A. Moreover, the numerical measure has a smooth density on R \ σ(A), and is absolutely continuous with respect to the Lebesgue measure dx if A is not a scalar matrix. We skip the details here, because the Hermitian case will be treated in full details in Section 3 below.

To conclude this section, we show that the numerical density is strictly positive in the interior of W (A). A little care is needed in the formulation of that result, because as we shall see in Section 4.2 the numerical density need not be a continuous function.

Proposition 2.6 If z 0 ∈ C is an interior point of W (A), then lim inf ǫ→0 1 ǫ 2 µ A ({z ∈ C | |z -z 0 | ≤ ǫ}) > 0 . (9)
Proof. If z 0 is an interior point of W (A), it is shown in [START_REF] Joswig | On the numerical range map[END_REF]Proposition 2.11] that the preimage Φ -1 A (z 0 ) contains at least one regular point x 0 . Let V be an open geodesic ball centered at x 0 ∈ ∂B n whose closure does not intersect Γ(A). Proceeding as in the proof of Lemma 2.3, we find

µ A ({z ∈ C | |z -z 0 | ≤ ǫ}) ≥ 1 ω n |z-z 0 |≤ǫ Nz∩V dν(x) 2∆(x) 1/2 dz .
If ǫ > 0 is sufficiently small, the integral inside the curly brackets is a smooth and positive function of z, and (9) follows.

The direct sum formula

Let p, q ∈ N * and n = p + q. Given A ∈ M p (C) and B ∈ M q (C), the direct orthogonal sum of A and B is the matrix A ⊕ B ∈ M n (C) defined by

A ⊕ B = A 0 0 B .
In this situation, we have a formula for the numerical measure µ A⊕B in terms of µ A and µ B .

Proposition 2.7 For any φ ∈ C 0 (C), we have [START_REF] Gardner | The Brunn-Minkoski inequality[END_REF] where B(p, q) is Euler's beta function

C φ(z) dµ A⊕B (z) = 1 B(p, q) C C 1 0 φ(tz ′ + (1 -t)z ′′ ) t p-1 (1 -t) q-1 dt dµ A (z ′ ) dµ B (z ′′ ) ,
B(p, q) = 1 0 t p-1 (1 -t) q-1 dt = Γ(p)Γ(q) Γ(p + q) .
Proof. Any unit vector x ∈ ∂B n can be written as

x = √ t u √ 1 -t v ,
where u ∈ ∂B p , v ∈ ∂B q , and t ∈ [0, 1]. Up to negligible sets, the map (u, v, t) → x defines a diffeomorphism from ∂B p × ∂B q × [0, 1] onto ∂B n . With this parametrization it is not difficult to verify that the Euclidean measure on ∂B n has the following expression

dσ n (x) = 1 2 t p-1 (1 -t) q-1 dt dσ p (u) dσ q (v) . Equivalently, since σ n (∂B n ) = ω n = 2π n /Γ(n), the normalized Euclidean measure satisfies dσ n (x) = 1 B(p, q) t p-1 (1 -t) q-1 dt dσ p (u) dσ q (v) .
Thus, using definition (1) and the fact that (A ⊕ B)x, x = t Au, u + (1t) Bv, v , we easily obtain

C φ(z) dµ A⊕B (z) = 1 B(p, q) ∂B p ∂B q 1 0 φ t Au, u + (1 -t) Bv, v t p-1 (1 -t) q-1 dt dσ q (v) dσ p (u) = 1 B(p, q) C C 1 0 φ(tz ′ + (1 -t)z ′′ ) t p-1 (1 -t) q-1 dt dµ A (z ′ ) dµ B (z ′′ ) ,
which is the desired result.

As an application, if we choose φ(z) = e -iξ•z in Proposition 2.7, we obtain the following relation between the Fourier transforms of the measures µ A , µ B and µ A⊕B :

μA⊕B (ξ) = 1 B(p, q) 1 0 μA (tξ) μB ((1 -t)ξ) t p-1 (1 -t) q-1 dt , ξ ∈ R 2 . ( 11
)
The formula given in Proposition 2.7 can be generalized in a straightforward way to a direct sum with an arbitrary number of terms. Assume that

A = A 1 ⊕ • • • ⊕ A k where A j ∈ M p j (C), so that A ∈ M n (C) with n = p 1 + • • • + p k . Setting p = (p 1 , . . . , p k ), we denote B(p) = D k-1 t p 1 -1 1 • • • t p k -1 k dt 1 . . . dt k-1 = Γ(p 1 ) • • • Γ(p k ) Γ(p 1 + • • • + p k ) ,
where

t k = 1 -(t 1 + • • • + t k-1
) and D k-1 denotes the (k-1)-dimensional simplex

D k-1 = (t 1 , . . . , t k-1 ) ∈ R k-1 + t 1 + • • • + t k-1 ≤ 1 . (12)
Using [START_REF] Gardner | The Brunn-Minkoski inequality[END_REF] and proceeding by induction over k, we easily obtain the general formula

C φ(z) dµ A (z) = 1 B(p) C k D k-1 φ(t 1 z 1 + • • • + t k z k ) t p 1 -1 1 • • • t p k -1 k dt 1 . . . dt k-1 dµ A 1 (z 1 ) . . . dµ A k (z k ) , ( 13
)
where it is understood again that

t k = 1 -(t 1 + • • • + t k-1 ).
3 The numerical density of a normal matrix

If A ∈ M n (C
) is a normal matrix, the numerical measure µ A is entirely determined by the spectrum σ(A) = {λ 1 , . . . , λ n }. Indeed, we know that A is unitarily equivalent to the diagonal matrix diag(λ 1 , . . . , λ n ), and that a unitary conjugation does not affect the numerical measure. Using this observation and the direct sum formula of Section 2.2, we shall prove that the numerical density of A is a piecewise polynomial function, which can be characterized as a multivariate B-spline whose knots are the eigenvalues of A. We begin with the important particular case where all eigenvalues of A are colinear.

The Hermitian case

If A ∈ M n (C) is a Hermitian matrix, then W (A) ⊂ R and the numerical measure µ A is therefore supported on the real axis. Assuming that A is not a multiple of the identity matrix, we show in this section that µ A is absolutely continuous with respect to Lebesgue's measure on R, and we give a simple characterization of the numerical density f A = dµ A /dx. The result is:

Proposition 3.1 If A ∈ M n (C
) is a nonscalar Hermitian matrix, the numerical density f A : R → R + is the normalized B-spline of degree n -2 whose knots are the eigenvalues of A.

To make the statement clear, we briefly recall the definition and some elementary properties of the classical B-splines [START_REF] De Boor | A Practical Guide to Splines[END_REF]. If λ 1 , . . . , λ n ∈ R are pairwise distinct, the (n-1) th divided difference of a continuous function g at the points λ 1 , . . . , λ n is the quantity

δ n-1 [λ 1 , . . . , λ n ]g = n j=1 g(λ j ) k =j (λ j -λ k ) . ( 14
)
This is the leading coefficient of the unique polynomial of degree at most n-1 which agrees with g at the points λ 1 , . . . , λ n . It is easy to verify that the right-hand side of ( 14) is a completely symmetric function of the variables λ j . If g ∈ C n-1 (R), the divided difference can be extended by continuity to arbitrary (not necessarily distinct) values of λ 1 , . . . , λ n , and we have the integral formula:

δ n-1 [λ 1 , . . . , λ n ]g = D n-1 g (n-1) (t 1 λ 1 + • • • + t n λ n ) dt 1 . . . dt n-1 , (15)
where D n-1 is the (n-1)-dimensional simplex defined in [START_REF] Helgason | The Radon transform, 2nd edition[END_REF] and

t n = 1 -(t 1 + • • • + t n-1 ).
In what follows, we shall always assume that the set S = {λ 1 , . . . , λ n } is not reduced to a single point, so that the (n-1) th divided difference is well-defined as soon as g is of class C n-2 in a neighborhood of S. With these notations, the normalized B-spline of degree n -2 with knots λ 1 , . . . , λ n is the function B : R → R defined by

B(x) ≡ B[λ 1 , . . . , λ n ](x) = (n -1) δ n-1 [λ 1 , . . . , λ n ](• -x) n-2 + , x ∈ R , (16)
where

(• -x) n-2 + denotes the map y → max(0, y -x) n-2 . If λ 1 ≤ • • • ≤ λ n , it is not difficult to show that B(x) vanishes identically outside [λ 1 , λ n ],
and coincides with a polynomial of degree at most n -2 on each nonempty interval (λ j , λ j+1 ). Moreover, if m j denotes the multiplicity of λ j in S, one can verify that B(x) is continuous at x = λ j together with its derivatives up to order

d j = n -2 -m j , provided d j ≥ 0. If d j = -1, then B(x) is discontinuous at λ j . Finally, we shall see below that B(x) is positive on (λ 1 , λ n ) and that R B(x) dx = 1.
Proof of Proposition 3.1. Let A ∈ M n (C) be a normal matrix with eigenvalues λ 1 , . . . , λ n . To compute the numerical density, we can assume without loss of generality that A is diagonal, namely

A = A 1 ⊕ • • • ⊕ A n with A j = λ j ∈ M 1 (C).
Thus we can use the direct sum formula [START_REF] Hausdorff | Der Wertevorrat einer Bilinearform[END_REF] with k = n and p = (1, . . . , 1). Since B(p) -1 = (n -1)! and µ A j = δ λ j for j = 1, . . . , n, we obtain the relation

C φ(z) dµ A (z) = (n -1)! D n-1 φ(t 1 λ 1 + • • • + t n λ n ) dt 1 . . . dt n-1 , (17)
for any continuous function φ :

C → C. Assume now that A is Hermitian, so that λ 1 , . . . , λ n ∈ R, and that σ(A) is not reduced to a single point. If B = B[λ 1 , . . . λ n ] is the normalized B-spline defined by (16), we claim that R φ(x)B[λ 1 , . . . , λ n ](x) dx = (n -1)! D n-1 φ(t 1 λ 1 + • • • + t n λ n ) dt 1 . . . dt n-1 . (18)
Indeed, it is clearly sufficient to prove [START_REF] Kato | Perturbation theory for linear operators[END_REF] for compactly supported functions φ ∈ C 0 (R). Moreover, since both members of (18) depend continuously on λ 1 , . . . , λ n , we can also assume that the eigenvalues of A are all distinct. In that case, it follows immediately from ( 14), ( 16) that

R φ(x)B[λ 1 , . . . , λ n ](x) dx = (n -1) δ n-1 [λ 1 , . . . , λ n ]Φ , where Φ(y) = R φ(x)(y -x) n-2 + dx. Since Φ (n-1) = (n -2)! φ,
we can use the integral formula (15) to evaluate the divided difference in the right-hand side, and we obtain [START_REF] Kato | Perturbation theory for linear operators[END_REF]. Now, comparing ( 17) and ( 18), we conclude that µ A is absolutely continuous with respect to Lebesgue's measure on R, and that the numerical density

f A = dµ A /dx is precisely the normalized B-spline B[λ 1 , . . . , λ n ].
Incidentally, the argument above shows that B is positive on its support (see Proposition 2.6) and that R B(x) dx = 1.

The quasi-Hermitian case

We say that a matrix A ∈ M n (C) is quasi-Hermitian if the numerical range W (A) ⊂ C has empty interior, i.e. W (A) is a single point or a line segment. In such a case, there exist z ∈ C, θ ∈ S 1 , and a Hermitian matrix H such that A = zI n + e iθ H (in particular, A is normal). Indeed, if we choose z, θ such that W (A) ⊂ z + e iθ R, the matrix H = e -iθ (A -zI n ) satisfies W (H) ⊂ R and is therefore Hermitian. Since Ax, x = z + e iθ Hx, x for all x ∈ ∂B n , it is clear from (1) that the numerical measure µ A is just the image of µ H under the affine isometry w → z + e iθ w. Combining this remark with Proposition 3.1, we thus obtain a precise characterization of the numerical measure of any quasi-Hermitian matrix.

The normal case

Finally, we consider the case of a normal matrix A ∈ M n (C) whose numerical range W (A) has nonempty interior. This, of course, is possible only if n ≥ 3. By Proposition 2.4, the numerical measure µ A is absolutely continuous with respect to Lebesgue's measure, and we have the following characterization of the numerical density f A = dµ A /dz:

Proposition 3.2 If A ∈ M n (C) is a normal matrix whose numerical range W (A) ⊂ C has nonempty interior, the numerical density f A : C → R + is the bivariate B-spline of degree n -3 whose knots are the eigenvalues of A.
The reader is referred here to the work of W. Dahmen [START_REF] Dahmen | On multivariate B-splines[END_REF], where multivariate B-splines are defined and studied in detail. To make the connection with the numerical density of a normal matrix, we use the relation [START_REF] Joswig | On the numerical range map[END_REF], which corresponds to formula (2.2) in [START_REF] Dahmen | On multivariate B-splines[END_REF]. In the rest of this section, we assume that A ∈ M n (C) is a normal matrix whose eigenvalues λ 1 , . . . , λ n are not colinear, and we often identify the complex plane C with R 2 .

Proof of Proposition 3.2. Assume first that n = 3. Then W (A) ⊂ C ≃ R 2 is the 2-simplex with vertices λ 1 , λ 2 , λ 3 , and using the change of variables

w = t 1 λ 1 + t 2 λ 2 + (1 -t 1 -t 2 )λ 3 in (17) we easily obtain C φ(z) dµ A (z) = 1 |W (A)| W (A) φ(w) dw .
This shows that the numerical measure µ A is uniformly distributed on W (A). The numerical density is thus a multiple of the characteristic function of W (A), which (by definition) is the bivariate B-spline of degree zero with knots λ 1 , λ 2 , λ 3 .

We now assume that n ≥ 4. Then we can choose n vectors v 1 , . . . , v n ∈ R n-1 such that 1) The (n-1)-simplex S ⊂ R n-1 with vertices v 1 , . . . , v n has unit volume; 2) P v i = λ i for i = 1, . . . , n, where P : R n-1 → R 2 is defined by P (x 1 , . . . , x n-1 ) = (x 1 , x 2 ). This (elementary) claim is proved in [6, Section 2]. The simplex S ⊂ R n-1 is not uniquely defined, but any choice satisfies P S = W (A) = conv(σ(A)). Returning to [START_REF] Joswig | On the numerical range map[END_REF], we have

C φ(z) dµ A (z) = (n -1)! D n-1 φ(P (t 1 v 1 + • • • + t n v n )) dt 1 . . . dt n-1 = S φ(P w) dw ,
where the second equality is obtained by applying the change of variables (t 1 , . . . , t n-1 )

→ w = t 1 v 1 + • • • + t n v n , with t n = 1 -(t 1 + • • • + t n-1
). This shows that the numerical measure is the image under the projection P of the Lebesgue measure on the simplex S ⊂ R n-1 . Given z ∈ C ≃ R 2 , the numerical density f A (z) is thus the (n-3)-dimensional measure of the simplex P -1 z ∩ S. This is precisely the definition given in [START_REF] Dahmen | On multivariate B-splines[END_REF] of a bivariate B-spline of degree n -3 with knots λ 1 , . . . , λ n .

We conclude by listing a few properties of the numerical density f A , which follow from [6, Theorem 4.1]. Let Σ A ⊂ C be the union of all line segments joining pairs of eigenvalues of A. Using Lemma 2.2, it is straightforward to verify that Σ A is exactly the set of critical values of the numerical map Φ A . Then f A is a polynomial of total degree n -3 in each connected component of C \ Σ A . Moreover, if ℓ ⊂ C is a straight line passing through a pair of eigenvalues of A, the numerical density is continuous across ℓ ∩ W (A) together with its derivatives up to order n -2m ≥ 0, where m ≥ 2 is the number of eigenvalues of A (counted with multiplicities) which belong to ℓ. If m = n -1, then f A is discontinuous on ℓ ∩ W (A). In particular, in the generic case where no straight line contains more than two eigenvalues of A, the numerical density

f A is of class C n-4 if n ≥ 4.
Remark. In the Hermitian and normal cases, the fact that the numerical density f A is a projection of the characteristic function of a convex set, which is a log-concave function, implies that the density is itself a log-concave function. We thus have the following inequality

f A (λz + (1 -λ)z ′ ) ≥ f A (z) λ f A (z ′ ) 1-λ , (19)
for all z, z ′ ∈ C and all λ ∈ (0, 1). This follows from the Prékopa-Leindler inequality, see e.g. [10, Section 9]. We warn the reader that this property does not extend to nonnormal matrices, as we can see already from the two-dimensional case considered in Section 4.1.

The Radon transform of the numerical measure

Our purpose in this section is to derive a representation formula for the numerical density of a nonnormal matrix A ∈ M n (C). Our approach is based on a natural expression of the Radon transform of the numerical measure µ A in terms of the Hermitian matrices H(θ) defined in [START_REF] Atiyah | Lacunas for hyperbolic differential operators with constant coefficients[END_REF]. By definition, the Radon transform of µ A is the family Rµ A = {Rµ A (θ) | θ ∈ S 1 }, where Rµ A (θ) denotes the Borel measure on R defined by

(Rµ A (θ))(I) = µ A ({z ∈ C | Re (e -iθ z) ∈ I}) ,
for any open set I ⊂ R. In other words, Rµ A (θ) is the image of the measure µ A under the orthogonal projection in C ≃ R 2 onto the line e iθ R. The fundamental observation is:

Proposition 4.1 For any θ ∈ S 1 = R/(2πZ), one has Rµ A (θ) = µ H(θ) .
Proof. The definition (2) implies that Re (e -iθ Ax, x ) = H(θ)x, x for any x ∈ ∂B n . Thus, for any open set I ⊂ R, we have

(Rµ A (θ))(I) = σ({x ∈ ∂B n | Re (e -iθ Ax, x ) ∈ I}) = σ({x ∈ ∂B n | H(θ)x, x ∈ I}) = µ H(θ) (I) ,
which proves the claim.

Proposition 4.1 shows that the numerical measure of an arbitrary matrix A ∈ M n (C) is entirely determined by the one-dimensional measures associated with the Hermitian matrices {H(θ) | θ ∈ S 1 }. If we assume that W (A) has nonempty interior, which is always the case if A is nonnormal, the matrix H(θ) is nonscalar for every θ ∈ S 1 and Proposition 3.1 shows that its numerical density is the Bspline B[λ 1 (θ), . . . , λ n (θ)], where λ 1 (θ), . . . , λ n (θ) are the eigenvalues of H(θ). In that case, the result of Proposition 4.1 can be stated in the following equivalent form

R f A (e iθ (x + iy)) dy = B[λ 1 (θ), . . . , λ n (θ)](x) , (20)
where equality holds for all θ ∈ S 1 and almost all x ∈ R since the numerical density f A belongs to L 1 (C).

Our goal is to invert the Radon transform [START_REF] Lax | Differential equations, difference equations and matrix theory[END_REF] to obtain a representation formula for the numerical density f A . The general results established in [START_REF] Helgason | The Radon transform, 2nd edition[END_REF] show that

(Jf A )(x + iy) = 1 4π S 1 B[λ 1 (θ), . . . , λ n (θ)](x cos θ + y sin θ) dθ , (21)
where

J = (-∆) -1/2 is the Riesz potential defined by (Jf )(z) = 1 2π C |z -z ′ | -1 f (z ′ ) dz ′ .
The idea is thus to apply the nonlocal operator (-∆) 1/2 to both sides of ( 21), but since we are not dealing with smooth functions we have to differentiate in the sense of distributions. As a preliminary remark, if f (x, y) = g(x cos θ + y sin θ) for some test function g : R → R and some fixed θ ∈ S 1 , a direct calculation shows that (∆f )(x, y) = g ′′ (x cos θ + y sin θ), and a standard interpolation argument allows us to conclude that (-∆) 1/2 f (x, y) = Hg ′ (x cos θ + y sin θ), where Hg ′ denotes the Hilbert transform of the derivative g ′ . Using this observation, we easily obtain the representation formula

f A (x + iy) = 1 4π S 1 HB ′ [λ 1 (θ), . . . , λ n (θ)](x cos θ + y sin θ) dθ , (22)
where both sides define integrable functions of z = x + iy ∈ C, and equality holds almost everywhere. So we have shown:

Proposition 4.2 If the numerical range of a matrix A ∈ M n (C) has nonempty interior, the numerical density of A can be represented as in [START_REF] Petrowsky | On the diffusion of waves and lacunas for hyperbolic equations[END_REF] for almost all (x, y) ∈ R 2 .

In the rest of this section, we shall apply Proposition 4.2 to compute the numerical density of a two-dimensional nonnormal matrix. We shall also consider the interesting particular situation where the numerical density is radially symmetric, in which case the representation formula takes a simpler form. Proposition 4.2 will be used again in Section 6 to derive some important properties of the numerical measure in the general case.

The two-dimensional case

Let A ∈ M 2 (C), and assume that the numerical range W (A) has nonempty interior. As is well-known [START_REF] Horn | Topics in matrix analysis[END_REF], W (A) is then a filled ellipse, and without loss of generality we can assume that this ellipse is centered at the origin and that its major axis is aligned with the real axis of the complex plane. In that case, up to a unitary conjugation, the matrix A has the following form

A = -c 2b 0 c , W (A) = x + iy ∈ C x 2 a 2 + y 2 b 2 ≤ 1 , (23)
where b > 0, c ≥ 0, and

a = (b 2 + c 2 ) 1/2 . For θ ∈ S 1 = R/(2πZ), let H(θ) = 1 2 (e -iθ A + e iθ A * ) = -c cos θ b e -iθ b e iθ c cos θ .
The eigenvalues of H(θ) are ±λ(θ), where λ

(θ) = (b 2 + c 2 cos 2 θ) 1/2 . It follows that B[-λ(θ), λ(θ)](s) = 1 2λ(θ) 1 [-λ(θ),λ(θ)] (s) .
According to [START_REF] Petrowsky | On the diffusion of waves and lacunas for hyperbolic equations[END_REF], we have to differentiate this expression with respect to s, which yields a linear combination of Dirac masses, and to apply the Hilbert transformation, which is the convolution with the distribution p.v. 1 πs . We obtain

HB ′ [-λ(θ), λ(θ)](s) = 1 2πλ(θ) p.v. 1 s + λ(θ) -p.v. 1 s -λ(θ) = 1 π p.v. 1 λ(θ) 2 -s 2 , hence f A (x + iy) = 1 4π 2 p.v. S 1 1 b 2 + c 2 cos 2 θ -(x cos θ + y sin θ) 2 dθ . ( 24
)
It remains to compute the right-hand side of [START_REF] Teicher | Almost certain convergence in double arrays[END_REF], which is a simple exercise in complex analysis. Setting z = x + iy and w = e 2iθ , we first observe that

b 2 + c 2 cos 2 θ -(x cos θ + y sin θ) 2 = 1 4w (c 2 -z 2 ) + 2w(a 2 + b 2 -|z| 2 ) + w 2 (c 2 -z2 ) , hence f A (z) = 1 iπ 2 p.v. |w|=1 1 (c 2 -z 2 ) + 2w(a 2 + b 2 -|z| 2 ) + w 2 (c 2 -z2 ) dw . (25)
In [START_REF] Toeplitz | Das algebraische Analogon zu einem Satze von Fejér[END_REF], the roots of the denominator are

w ± = 1 c 2 -z2 -a 2 -b 2 + |z| 2 ∓ (a 2 + b 2 -|z| 2 ) 2 -(c 2 -z 2 )(c 2 -z2 ) = 1 c 2 -z2 -a 2 -b 2 + |z| 2 ∓ 4(a 2 b 2 -b 2 x 2 -a 2 y 2 ) . (26)
We can therefore distinguish between two cases:

1. The point z = x + iy belongs to the interior of W (A). Then the expression under the square root is positive, and it is easy to verify that |w + | > 1, |w -| < 1 (in the limiting case where z = ±c is a focus of the ellipse, one can set w -= 0 and w + = ∞.) Thus the principal value in [START_REF] Toeplitz | Das algebraische Analogon zu einem Satze von Fejér[END_REF] is not needed, and the residue theorem shows that

f A (z) = 1 π 1 (a 2 + b 2 -|z| 2 ) + (c 2 -z2 )w - = 1 2π 1 a 2 b 2 -b 2 x 2 -a 2 y 2 . ( 27
)
2. The point z = x + iy lies outside W (A). Then the expression under the square root in [START_REF] Unser | On the Asymptotic Convergence of B-Spline Wavelets to Gabor Functions[END_REF] is negative, and one verifies that |w ± | = 1. Thus the integrand in ( 25) is holomorphic outside the unit circle {|w| = 1} and decreases like 1/|w| 2 at infinity. It follows that

I z (r) := 1 iπ 2 |w|=r 1 (c 2 -z 2 ) + 2w(a 2 + b 2 -|z| 2 ) + w 2 (c 2 -z2 ) dw = 0 , for all r = 1, hence f A (z) = 1 2 (I z (1+) + I z (1-)) = 0
. Summarizing, we have shown that the numerical density of the matrix [START_REF] Serre | Matrices. Theory and applications[END_REF] is the function f A ∈ L 1 (C) defined by [START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF] inside W (A), and vanishing identically outside W (A). Note that sing supp(µ A ) = ∂W (A) = Σ A , in agreement with Lemma 2.3, and that f A (z) blows up when z converges to the boundary of W (A) from inside. In particular f A is not log-concave, in contrast to what happens when A is normal.

The radially symmetric case

It sometimes happens that the numerical range of a matrix A ∈ M n (C) is a disk in the complex plane and that the numerical density is radially symmetric about the center. In such a case, it is possible to obtain a representation formula which is simpler than [START_REF] Petrowsky | On the diffusion of waves and lacunas for hyperbolic equations[END_REF]. Indeed, assume that the disk W (A) is centered at the origin, and let R > 0 denote the numerical radius of A (or any larger positive number). We set

f A (z) = F A (R 2 -|z| 2 ) , z ∈ C , |z| ≤ R , (28)
where F A : [0, R 2 ] → R + has to be determined. Since the numerical measure is invariant under rotations about the origin, it follows from Proposition 4.1 that the projected measure µ H(θ) does not depend on θ. In analogy with (28), if f H denotes the numerical density of H(θ) for any θ ∈ S 1 , we set

f H (x) = F H (R 2 -x 2 ) , x ∈ R , |x| ≤ R . ( 29
)
We then have the following result:

Proposition 4.3
If the numerical density of a nonscalar matrix A ∈ M n (C) is radially symmetric, the functions F A , F H defined in (28), (29) satisfy the relations

F H (t) = t 0 F A (t -s) 1 √ s ds , F A (s) = 1 π d ds s 0 F H (s -t) 1 √ t dt , ( 30
)
for t, s ∈ [0, R 2 ].
Proof. By the definition of the Radon transformation, we have f H (x) = R f A (x+iy) dy for all x ∈ R. Using (28), (29) and the support property, we thus find

F H (R 2 -x 2 ) = y 2 ≤R 2 -x 2 F A (R 2 -x 2 -y 2 ) dy = R 2 -x 2 0 F A (R 2 -x 2 -s) 1 √ s ds , for all x ∈ [-R, R]. Setting t = R 2 -x 2
, we obtain the first relation in (30). So far, we have shown that

F H = π 1/2 IF A , where I is the Riesz potential (If )(t) = 1 √ π t 0 f (t -s) 1 √ s ds , t > 0 .
Now, it is well known that (I 2 f )(t) = t 0 f (s) ds, see e.g. [12, Chapter V.5], thus the second relation in (30) follows from the first one.

Examples. Let A = (a ij ) be a complex matrix, and assume that there exists a nonzero integer k such that a ij = 0 whenever ji = k. Then W (A) is a disk centered at the origin, and the numerical density f A is radially symmetric. Indeed, given any θ ∈ S 1 , the map

x =    x 1 . . . x n    → x θ =    e iθ x 1 . . . e inθ x n   
is a measure-preserving isomorphism of the unit sphere ∂B n , and our assumption on A implies that Ax θ , x θ = e ikθ Ax, x for all x ∈ ∂B n . This proves that a rotation of angle kθ about the origin does not affect the numerical measure µ A . Since θ is arbitrary, the numerical density f A is necessarily radially symmetric.

The simplest example in this category is the 2 × 2 Jordan block

A 2 = 0 2 0 0 . ( 31
)
Here W (A 2 ) = {z ∈ C | |z| ≤ 1}, and applying [START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF] with a = b = 1 we find

f A 2 (z) = 1 2π 1 1 -|z| 2 1 {|z|<1} . (32)
Alternatively, since the eigenvalues of

H 2 = 1 2 (A 2 + A * 2 ) are ±1, we have f H 2 = 1 2 1 [-1,1]
and applying Proposition 4.3 we easily obtain (32).

As a more interesting application, consider the 3 × 3 matrix

A 3 =   0 a 0 0 0 b 0 0 0   , (33)
where a, b ∈ C and |a| + |b| > 0. Multiplying A 3 with a positive constant, we can assume that

|a| 2 + |b| 2 = 4. Then W (A 3 ) = {z ∈ C | |z| ≤ 1}, and
f A 3 (z) = 1 π log 1 + 1 -|z| 2 |z| 1 {0<|z|<1} . (34)
Indeed, the eigenvalues of H 3 = 1 2 (A 3 + A * 3 ) are 0 and ±1, hence f H 3 (x) = B[-1, 0, 1](x) = (1-|x|) + . Applying Proposition 4.3, we obtain (34) by a straightforward calculation. Alternatively, if we take (a, b) = (0, 2), we can use the direct sum formula [START_REF] Gardner | The Brunn-Minkoski inequality[END_REF] to deduce (34) from (32). Here again, we have sing supp(µ A 3 ) = ∂W (A 3 ) ∪ {0} = Σ A 3 , in agreement with Lemma 2.3. Remark that f A 3 (z) = O((1-|z|) 1/2 ) as |z| → 1-, so that the singularity of the numerical measure at the boundary is weaker than it was for A 2 . This reflects the fact that f H 3 (x) = (1 -|x|) + is Lipschitz continuous, whereas f H 2 (x) = 1 [-1,+1] (x) had jump discontinuities. However, we observe that the (logarithmic) singularity of f A 3 at the origin is much stronger than the (square root) singularity at the boundary of W (A 3 ). As we shall see in Section 5, this is a nongeneric concentration phenomenon due to the fact that the component of the critical set Σ A 3 associated with the eigenvalue 0 of H(θ) is reduced to a single point. We have here the rare instance of an unbounded numerical density for a matrix of size n ≥ 3.

Another interesting conclusion that can be drawn from this example is that the numerical density µ A does not determine the matrix

A ∈ M n (C) up to unitary conjugations if n ≥ 3. Indeed, if we set (a, b) = (0, 2) and (a, b) = ( √ 2, √ 2
) in the definition (33) of A 3 , the resulting matrices are not even similar, yet they have the same numerical density, given by (34). not claim much originality here: the material of this section is essentially borrowed from [START_REF] Kippenhahn | Über den Wertevorrat einer Matrix[END_REF][START_REF] Bazer | Lacunas of the Riemann matrix of symmetric-hyperbolic systems in two space variables[END_REF][START_REF] Jonckheere | Differential topology of numerical range[END_REF][START_REF] Joswig | On the numerical range map[END_REF], and is reproduced below for the reader's convenience.

As was shown by Kippenhahn [START_REF] Kippenhahn | Über den Wertevorrat einer Matrix[END_REF], the singular set Σ A has a natural description in terms of the eigenvalues λ 1 (θ), . . . , λ n (θ) of the Hermitian matrices H(θ) defined in [START_REF] Atiyah | Lacunas for hyperbolic differential operators with constant coefficients[END_REF]. To see that, we first recall that these eigenvalues can be numbered in such a way that they are real-analytic functions of θ ∈ R, see [START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF]. By analyticity, any two eigenvalues either coincide for all θ ∈ R or cross at most a finite number of times on each compact interval. As a consequence, there exists an integer m ≤ n such that H(θ) has exactly m distinct eigenvalues for all θ ∈ [0, π) \ Θ, where Θ ⊂ [0, π) is a finite set. Since H(θ + π) = -H(θ), it follows that we can number the eigenvalues in such a way that

σ(H(θ)) = {λ 1 (θ), . . . , λ m (θ)} , for all θ ∈ R ,
where λ 1 (θ), . . . , λ m (θ) are pairwise distinct and have constant multiplicities outside the crossing set Θ + πZ. Let τ : {1, . . . , m} → {1, . . . , m} be the permutation defined by

{λ τ (1) (θ), . . . , λ τ (m) (θ)} = {-λ 1 (θ + π), . . . , -λ m (θ + π)} , (35)
for any θ ∈ [0, π) \ Θ (hence for all θ ∈ R). If we decompose τ into disjoint cycles C 1 , . . . , C k , we can associate with each cycle C J its length ℓ J and its multiplicity m J , the latter being defined as the multiplicity of λ j (θ) as an eigenvalue of H(θ) for any j ∈ C J . By construction, we have

ℓ 1 + • • • + ℓ k = m , and
ℓ 1 m 1 + • • • + ℓ k m k = n .
Moreover, if j ∈ C J , then λ j (θ + ℓ J π) = (-1) ℓ J λ j (θ), hence λ j (θ) is periodic with period ℓ J π if ℓ J is even and 2ℓ J π if ℓ J is odd. Note however that these periods are not necessarily minimal. Now, we associate with each cycle C J of the permutation τ a closed curve C J ⊂ C defined by

C J = e iθ (λ j (θ) + iλ ′ j (θ)) θ ∈ R , (36)
where j ∈ {1, . . . , m} is any element of the cycle C J . Since λ j (θ + 2ℓ J π) = λ j (θ), it is clear that C J is indeed a closed curve, and the definition of the permutation τ shows that the right-hand side of (36) does not depend on the choice of j ∈ C J . Equivalently, we can define C J as the union over all j ∈ C J of the curve segments {e iθ (λ

j (θ) + iλ ′ j (θ)) | θ ∈ [0, π]}. Let C A = C 1 ∪ . . . ∪ C k ⊂ C ,
and let C ′ A ⊂ C be the bitangent set of C A , namely the set of all line segments joining pairs of points of C A at which C A has the same tangent line. With these definitions, we have the following useful characterization of the singular set Σ A :

Proposition 5.1 [17, Theorem 3.5] Σ A = C A ∪ C ′ A .
Proof. According to Lemma 2.2, Σ A is the set of all complex numbers of the form Ax, x where x ∈ ∂B n is a normalized eigenvector of the Hermitian matrix H(θ) for some θ ∈ [0, π]. To describe that set, fix j ∈ {1, . . . , m}, θ 0 ∈ R, and for θ in a neighborhood of θ 0 let

z j (θ) = Ax j (θ), x j (θ) ,
where x j (θ) is a normalized eigenvector of H(θ) associated with the eigenvalue λ j (θ) and depending smoothly on θ. General results in perturbation theory imply that such an eigenvector indeed exists [START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF]. Using the definition (2) of H(θ), we find

Re (e -iθ z j (θ)) = Re e -iθ Ax j (θ), x j (θ) = H(θ)x j (θ), x j (θ) = λ j (θ) ,

Im (e -iθ z j (θ)) = Im e -iθ Ax j (θ),

x j (θ) = H ′ (θ)x j (θ), x j (θ) = λ ′ j (θ) ,
since H(θ)x ′ j (θ), x j (θ) + H(θ)x j (θ), x ′ j (θ) = 2λ j (θ) Re x j (θ), x ′ j (θ) = 0 due to the normalization condition. Thus

z j (θ) = e iθ (λ j (θ) + iλ ′ j (θ)) , hence z ′ j (θ) = ie iθ (λ j (θ) + λ ′′ j (θ)) . ( 37
)
These relations show that the curve

C J = {z j (θ) | θ ∈ R} is tangent, at each point z j (θ), to the straight line L j (θ) = {e iθ (λ j (θ) + iα) | α ∈ R} = {z ∈ C | Re (ze -iθ ) = λ j (θ)} . (38)
In other words C J is the envelope of the family of straight lines L j (θ), for any j ∈ C J . Since C J ⊂ Σ A by construction, we have shown that Σ A contains the curve

C A = C 1 ∪ . . . ∪ C k .
However, it is important to realize that Σ A can be larger than C A if the crossing set Θ ⊂ [0, π) defined above is nonempty. Indeed, assume that λ j (θ 0 ) = λ p (θ 0 ) for some θ 0 ∈ Θ and some j, p ∈ {1, . . . , m} with j = p. For θ in a neighborhood of θ 0 , let x j (θ), x p (θ) be smooth, normalized eigenvectors of H(θ) corresponding to λ j (θ), λ p (θ) respectively. Using the same notations as above, we have z j (θ 0 ) = z p (θ 0 ) in general, because λ ′ j (θ 0 ) = λ ′ p (θ 0 ). Now, since x j (θ), x p (θ) = 0 whenever λ j (θ) = λ p (θ), we also have x j (θ 0 ), x p (θ 0 ) = 0 by continuity. In particular, if α, β ∈ C satisfy |α| 2 + |β| 2 = 1, then x = αx j (θ 0 ) + βx p (θ 0 ) is a normalized eigenvector of H(θ 0 ), and a direct calculation yields

Im e -iθ 0 Ax, x = H ′ (θ 0 )x, x = |α| 2 λ ′ j (θ 0 ) + |β| 2 λ ′ p (θ 0 ) ,
whereas Re e -iθ 0 Ax, x = H(θ 0 )x, x = λ j (θ 0 ) = λ p (θ 0 ). This shows that Σ A contains the line segment [z j (θ 0 ), z p (θ 0 )], which by construction is tangent to the curve C J at z j (θ 0 ) and to the curve C P at z p (θ 0 ). Repeating the same argument for all eigenvalue crossings, we conclude that Σ A contains the whole bitangent set C ′ A . Finally, it is clear from Lemma 2.2 that all points of Σ A either belong to C A or to C ′ A .

Examples:

1. (The generic case) As is well known, within the space of all Hermitian matrices of size n ≥ 2, the set of matrices having a multiple eigenvalue is a finite union of submanifolds of codimension at least three. This implies that, for a generic matrix A ∈ M n (C), the Hermitian matrices H(θ) defined by (2) will have simple eigenvalues for all θ ∈ S 1 [START_REF] Jonckheere | Differential topology of numerical range[END_REF]. In that situation, we denote by λ 1 (θ) < λ 2 (θ) < . . . < λ n (θ) the eigenvalues of H(θ). Using the notations introduced above, we have m = n and the permutation τ defined in (35) is simply

τ = 1 2 . . . n n n -1 . . . 1 .
If n = 2k is even, τ has k cycles with ℓ 1 = • • • = ℓ k = 2; if n = 2k -1 is odd, τ has one fixed point and k -1 cycles of length 2. In all cases, the multiplicities m 1 , . . . , m k are all equal to 1, and the bitangent set C ′ A is empty by assumption. Thus Σ A = C A is the union of k = [ n+1 2 ] closed curves. It is not difficult to prove that the curve C 1 associated with the cycle C 1 = (1 n) is smooth, strictly convex, and contains all the other curves C 2 , . . . , C k in its interior [START_REF] Kippenhahn | Über den Wertevorrat einer Matrix[END_REF]. In particular, C 1 = ∂W (A). Consider now the curve C J = {z j (θ) | θ ∈ R} for some j = 1, n. If δ j (θ) = λ j (θ) + λ ′′ j (θ) is not identically zero, the formulas (37) show that the curvature of C J at any regular point z j (θ) is strictly positive:

κ j (θ) = 1 |z ′ j (θ)| = 1 |λ j (θ) + λ ′′ j (θ)| > 0 .
This means that the tangent vector z ′ j (θ) always rotates counterclockwise when θ is increased. Nevertheless, the whole curve C J is not convex in general, because it may have a finite number of singular points corresponding to zeros of δ j (θ). As is easily verified, simple zeros of δ j (θ) correspond to cusp points of the curve C J , see e.g. Fig. 1 where a generic example with n = 3 is represented. On the other hand, if δ j (θ) vanishes identically, then z ′ j (θ) ≡ 0 and the curve C J reduces to a single point. Under our generic assumptions, this can happen only if n is odd and j = (n + 1)/2. As an example, the singular set of the matrix A 3 defined in (33) includes the isolated point {0}.

(The normal case)

In contrast with the previous example, we now consider the particular case where the matrix A is normal. If λ 1 , . . . , λ m denote the distinct eigenvalues of A, it is straightforward to verify that the eigenvalues of the Hermitian matrices H(θ) defined in (2) are simply λ j (θ) = Re (λ j e -iθ ), j = 1, . . . , m. In view of (37), this means that z j (θ) = λ j for all θ ∈ R, hence the curve C j is reduced to the single point {λ j } for all j = 1, . . . , m. Needless to say, the permutation τ defined by ( 35) is the identity here. It follows that C A = σ(A) = {λ 1 , . . . , λ m }, and proceeding as in the proof of Proposition 5.1 we easily see that C ′ A is the set of all line segments joigning pairs of eigenvalues of A. This is in full agreement with the conclusions of Section 3.3.

Remark:

We have seen that the curve C A which generates the singular set Σ A is the envelope of the family of straight lines L j (θ) defined in (38). This geometric construction can be formulated in an equivalent way [START_REF] Atiyah | Hyperbolic differential equations and algebraic geometry (after Petrowsky)[END_REF], which is more conceptual and worth mentioning here. Assume for simplicity that the matrix A ∈ M n (C) is "generic" in the sense of Example 1 above, and let a(ξ) be the homogeneous polynomial of degree n defined by

a(ξ) = det ξ 0 I n + ξ 1 A 1 + ξ 2 A 2 , ξ = (ξ 0 , ξ 1 , ξ 2 ) ∈ R 3 ,
where A 1 , A 2 are as in [START_REF] Atiyah | Lacunas for hyperbolic differential operators with constant coefficients[END_REF]. Since a(ξ) has real coefficients, the equation a(ξ) = 0 defines, in homogeneous coordinates, an algebraic curve Γ in the projective plane RP 2 . Moreover, our genericity assumption on A implies that this curve is nonsingular: for each ξ ∈ Γ, the tangent line to Γ at ξ is uniquely defined and satisfies, in homogeneous coordinates, an equation of the form x 0 ξ 0 + x 1 ξ 1 + x 2 ξ 2 = 0. The set of all x = (x 0 , x 1 , x 2 ) obtained in this way is again a curve Γ ′ in RP 2 (the dual curve of Γ) given by the equation b(x) = 0 for some homogeneous polynomial b. As is shown in [START_REF] Fiedler | Geometry of the numerical range of matrices[END_REF], the degree of b does not exceed n(n -1) if n ≥ 2. Now, it is rather straightforward to verify that the curve C A defined as the envelope of the family of straight lines (38) is nothing but the restriction of the projective curve Γ ′ to the subspace

x 0 = 1, namely z = x + iy ∈ C A if and only if b(1, x, y) = 0. Thus C A is a real algebraic curve in C ≃ R 2 of degree at most n(n -1) if n ≥ 2.
In the language of partial differential equations, the algebraic variety Γ is the characteristic variety of the symmetric hyperbolic system (4), and we shall see in Section 6 that the dual variety Γ ′ is related to the singular support of the fundamental solution of (4). Note that our genericity assumption on A precisely means that system (4) is strictly hyperbolic.

We conclude this section with a brief discussion of the number of tangent lines to the algebraic curve C A which can be drawn from a given point. We recall that C A = C 1 ∪ . . . ∪ C k , where each C J is a closed curve associated with the cycle C J of the permutation (35). For all J ∈ {1, . . . , k} and all z ∈ C \ C J , we denote by N J (z) the number of straight lines that are tangent to the curve C J and contain the point z. Note that, since C J was itself defined as the envelope of a family of straight lines, the tangent line to C J is well defined even at singular points. In the degenerate case where C J reduces to a single point {z J }, the set of tangent lines should be understood as the pencil of all straight lines through z J . Now, if z ∈ C \ C A , we denote by N (z) the total number of tangents to the curve C A = C 1 ∪ . . . ∪ C k that can be drawn from the point z, with multiplicities taken into account:

N (z) = m 1 N 1 (z) + • • • + m k N k (z) , z ∈ C \ C A . (39)
The following elementary properties of N J (z) and N (z) will be useful.

Proposition 5.2 For each J ∈ {1, . . . , k}, the number

N J (z) is constant in each connected component of C \ C J . Moreover N (z) ≤ n for all z ∈ C \ C A .
Proof. Fix J ∈ {1, . . . , k} and pick j ∈ C J . For any z ∈ C \ C J , N J (z) is the number of zeros of the function

f j (θ, z) = λ j (θ) -Re (ze -iθ )
for θ in the interval [0, ℓ J π). When z is varied, this number can only change if f j (θ, z) has a double zero for some θ, but this would mean that z = z j (θ) = e iθ (λ j (θ) + iλ ′ j (θ)) ∈ C J , thus contradicting our assumption. Therefore N J (z) is necessarily constant in each connected component of C \ C J . On the other hand, we have the identity

det H(θ) -Re (ze -iθ )I n = k J=1 j∈C J f j (θ, z) m J . ( 40
)
Fix z ∈ C \ C A , and consider both sides of (40) as functions of θ. If multiplicities are taken into account, the number of zeros of the right-hand side on the interval [0, π) is precisely N (z). But the left-hand side, being a trigonometric polynomial of degree at most n, cannot have more than n zeros in [0, π). This proves the claim.

Remark 5.3 As was mentioned in the proof of Proposition 5.1, if z ∈ C \ C A it is possible that f j (θ, z) = f p (θ, z) = 0 for some θ ∈ [0, π) and some j = p. This is the case, in particular, whenever

z ∈ C ′ A . However, if z ∈ C \ Σ A , we have ∂ θ f j (θ, z) • ∂ θ f p (θ, z) > 0 whenever f j (θ, z) = f p (θ, z) = 0 .
Indeed, replacing A with A -zI n , we can assume without loss of generality that z = 0. If λ j (θ) = λ p (θ) = 0 for some θ ∈ [0, π) and some j = p, then λ ′ j (θ) and λ ′ p (θ) have necessarily the same sign, otherwise (37) would imply that the origin belongs to the line segment [z j (θ), z p (θ)], thus contradicting our assumption that z / ∈ C ′ A .

Proposition 5.2 asserts that the integer N J (z) can only change if z crosses the curve C J . In fact, if the crossing occurs at a regular point z ∈ C J , it is not difficult to verify that N J (z) is decreased by two units if z crosses C J from the convex to the concave side (i.e., in the direction of the local center of curvature), and increased by two units if z crosses C J in the converse direction, see [4, Section 4.1] or Section 6.2 below. These simple rules give an efficient algorithmic way to compute N J (z), and hence N (z), in concrete examples. Consider for instance Fig. 2, where the singular set Σ A of a generic matrix

A ∈ M 5 (R) is represented. Here Σ A = C A = C 1 ∪ C 2 ∪ C 3 , where C 1 is the boundary of W (A), C 2
is the closed curve with two swallowtails, and C 3 is the triangle with three cusps. The set C \ Σ A has 11 connected components, on which N (z) is equal to 5, 3, or 1. Using the crossing rules above, and the fact that m 1 = m 2 = m 3 in the present case, it is easy to verify that N (z) = 5 in six different regions: inside both swallowtails, inside the three tips of the triangle, and outside W (A). Such a result is definitely more cumbersome to obtain by counting directly the number of tangents to C A from a given point.

Qualitative properties of the numerical density

Equipped with the results of Sections 4 and 5, we now derive some of the main properties of the numerical density f A of an arbitrary matrix A ∈ M n (C). We first establish an explicit formula for the derivatives of order n -2, which allows us to prove that f A is polynomial in some distinguished regions of the complex plane which can be characterized geometrically. We next show that, for a generic matrix

A ∈ M n (C), the density f A is of class C n-3 if n ≥ 3.
Finally, as announced in the introduction, we prove that the fundamental solution of the linear hyperbolic system (4) can be represented in terms of derivatives of the numerical density f A . In particular, the lacunas of system (4) are precisely the regions where f A is polynomial of degree n -3.

Polynomial regions

Most of what we know about the numerical density of nonnormal matrices is based on the following result:

Proposition 6.1 Let A ∈ M n (C)
, where n ≥ 2, and let P = P (∂ x , ∂ y ) be a homogeneous differential operator of degree n -2. Then, for all z = x + iy ∈ C \ Σ A , we have

(Pf A )(z) = - (n-1)! 4π 2 f.p. S 1 P (cos θ, sin θ) ∆(θ, z) dθ , (41)
where

∆(θ, z) = det H(θ) -Re (ze -iθ )I n , θ ∈ S 1 . (42)
In (41) the symbol f.p. denotes the finite part of the integral in the sense of Hadamard, but in many situations it is sufficient to take simply the Cauchy principal value, as in [START_REF] Teicher | Almost certain convergence in double arrays[END_REF]. This is the case in particular for generic matrices in the sense of Section 5 (Example 1), because all multiplicities m 1 , . . . , m k are all equal to one and formula (40) then shows that the map θ → ∆(θ, z) has only simple zeros on

S 1 if z ∈ C \ Σ A .
In fact, using the analyticity of the integrand, it is possible to rewrite (41) in a slightly different form which is appropriate for further analysis. As in Section 4.1, we set w = e 2iθ and we observe that

P (cos θ, sin θ) = w 1-n 2 P (w) , ∆(θ, z) = w -n 2 ∆(w, z) ,
where (if z = x + iy)

P (w) = P w + 1 2 , w -1 2i , ∆(w, z) = det A 2 + wA * 2 -x w + 1 2 + y w -1 2i I n .
Thus ( 41) is equivalent to

(Pf A )(z) = i(n-1)! 4π 2 f.p. |w|=1 P (w) ∆(w, z) dw ≡ i(n-1)! 8π 2 |w|=1-ǫ P (w) ∆(w, z) dw + i(n-1)! 8π 2 |w|=1+ǫ P (w) ∆(w, z) dw , ( 43
)
where ǫ is any positive number that is sufficiently small, depending on z. In the particular case where P = 1 and A ∈ M 2 (C) is given by ( 23), we recover [START_REF] Toeplitz | Das algebraische Analogon zu einem Satze von Fejér[END_REF].

Proof of Proposition 6.1. We first consider the generic situation where the Hermitian matrices H(θ) have simple eigenvalues λ 1 (θ) < λ 2 (θ) < • • • < λ n (θ) for all θ ∈ S 1 . In that case, the B-spline representing the numerical density of µ H (θ) can be given an explicit expression using the divided difference formula ( 14), ( 16):

B[λ 1 (θ), . . . , λ n (θ)](s) = (n -1) n j=1 (λ j (θ) -s) n-2 + k =j (λ j (θ) -λ k (θ))
.

In particular, differentiating (n -1) times with respect to s, we find

B (n-1) [λ 1 (θ), . . . , λ n (θ)](s) = (n -1)! n j=1 δ(s -λ j (θ)) k =j (λ k (θ) -λ j (θ)) . ( 44
)
Let P = P (∂ x , ∂ y ) be a homogeneous differential operator of degree n -2. From the representation formula [START_REF] Petrowsky | On the diffusion of waves and lacunas for hyperbolic equations[END_REF], we deduce at least formally

(Pf A )(x + iy) = 1 4π S 1 P (cos θ, sin θ) HB (n-1) [λ 1 (θ), . . . , λ n (θ)](x cos θ + y sin θ) dθ . ( 45
)
To evaluate the integrand in the right-hand side, we start from (44) and recall that the Hilbert transform (with respect to the variable s) of the Dirac measure δ(sλ) is the distribution 1 π p.v. 1 s-λ . We also use the identity

n j=1 1 µ j k =j (µ k -µ j ) = n j=1 1 µ j ,
which holds for any collection of pairwise distinct nonzero complex numbers µ 1 , . . . , µ n . We thus find

HB (n-1) [λ 1 (θ), . . . , λ n (θ)](s) = - (n-1)! π p.v. n j=1 1 λ j (θ) -s = - (n-1)! π p.v. 1 det(H(θ) -sI n )
.

Setting s = x cos θ + y sin θ and inserting this expression into (45), we obtain (41).

The calculations so far are formal, but they can be justified if we assume that z = x + iy ∈ C \ Σ A . In that case, we know from the proof of Proposition 5.2 that the map θ → ∆(θ, z) has only simple zeros on S 1 , because this is the case for each of the factors f j (θ, z) in (40) and, by assumption, the eigenvalues λ j (θ) are all distinct for θ ∈ S 1 . It follows that the integral in ( 41) is well-defined in the sense of Cauchy's principal value, and depends smoothly on z ∈ C \ Σ A . This in turn implies that the density f A given by ( 22) is smooth on C \ Σ A , as it should be, and that the calculations above are correct.

To conclude the proof of Proposition 6.1, it remains to verify that (41) or (43) holds for an arbitrary matrix A ∈ M n (C). To do that, we first observe that the singular set Σ A ⊂ C is an upper-semicontinuous function of A in the sense that

δ(Σ B , Σ A ) ≡ sup z∈Σ B dist(z, Σ A ) ---→ B→A 0 .
Moreover the numerical density f A , together with its derivatives, depends continuously on A in C\Σ A . These rather classical facts can be established using, for instance, the representation formula (8) (we omit the details). On the other hand, it is not difficult to verify that the right-hand side of (43) depends continuously on A for each z ∈ C \ Σ A . The crucial point here is that the polynomial map w → ∆(w, z) keeps the same number of zeros on the unit circle (counted with multiplicities) if the matrix A is slightly varied; in particular, we can choose the same ǫ > 0 in (43) for all matrices in a neighborhood of A. This property can be established using the factorization (40), Remark 5.3, and general results for perturbations of eigenvalues of Hermitian matrices, see [START_REF] Kato | Perturbation theory for linear operators[END_REF][START_REF] Serre | Matrices. Theory and applications[END_REF]. Now, given A ∈ M n (C), there exists a sequence {A ℓ } ℓ∈N of generic matrices converging to A as ℓ → ∞. If z ∈ C \ Σ A , we know that equation (43) holds for A ℓ if ℓ is sufficiently large, hence taking the limit ℓ → ∞ and using the continuity properties mentioned above we obtain the desired equality.

As a consequence of Proposition 6.1, we now establish an important property of the numerical density in the regions of the complex plane where the number N (z) defined in (39) takes its maximal value n.

Corollary 6.2 Given

A ∈ M n (C), let Π A = z ∈ C \ Σ A N (z) = n . ( 46
)
If n ≥ 3, the numerical density f A is polynomial of degree at most n -3 on each connected component of

Π A . If n = 2, then f A = 0 on Π A .
Proof. We shall prove that Pf A vanishes identically on Π A for any homogeneous differential operator of order n -2. Indeed, in view of ( 40), (42), the assumption z ∈ Π A implies that the map θ → ∆(θ, z) has exactly n zeros (counting multiplicities) on [0, π). Equivalently, the polynomial ∆(w, z) has exactly n zeros on the unit circle {|w| = 1}. But since ∆(w, z) has degree n, this polynomial has no zeros outside the unit circle, and using Cauchy's theorem we conclude that the first integral in the last member of (43) vanishes. The second integral is also zero, because the numerator is a polynomial of degree at most n -2, while the denominator has degree exactly n, hence the integrand P (w)/ ∆(w, z) decays at least like |w| -2 as |w| → ∞. Thus Pf A ≡ 0 on Π A , and the conclusion follows.

Remark 6.3 If z = x + iy ∈ σ(A), then the polynomial ∆(w, z) has degree strictly less than n, and it follows from the above proof that z /

∈ Π A . Thus σ(A) ∩ Π A = ∅.
The set Π A is never empty, because it always contains the complement of the numerical range W (A), where the density f A vanishes identically [START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF][START_REF] Bazer | Lacunas of the Riemann matrix of symmetric-hyperbolic systems in two space variables[END_REF]. Moreover, in many situations, one or several components of Π A are contained in W (A), in which case Corollary 6.2 gives nontrivial informations on the numerical density. For instance, if A ∈ M n (C) is a normal matrix whose numerical range has nonempty interior, then N (z) = n for all z ∈ C \ σ(A), and it follows from Corollary 6.2 that f A is polynomial of degree at most n -3 in each connected component of C \ Σ A , in agreement with Proposition 3.2. In the same spirit, if

A = A 1 ⊕ A 2 with A 1 ∈ M n 1 (C) and A 2 ∈ M n 2 (C), it is easily verified that N (z) = n = n 1 + n 2 if z / ∈ W (A 1 ) ∪ W (A 2 ), hence f A is piecewise polynomial outside W (A 1 ) ∪ W (A 2 )
. Finally, Fig. 1 shows a typical example of a matrix A ∈ M 3 (C) for which Π A has a component inside W (A), on which the density f A is identically constant by Corollary 6.2.

Generic regularity results

Our purpose here is to establish regularity results for the numerical density f A in the whole complex plane, and not only outside the singular set Σ A . For simplicity, we assume henceforth that our matrix A ∈ M n (C) enjoys the following (generic) properties: H1: The eigenvalues of (2) satisfy λ 1 (θ) < . . . < λ n (θ) for all θ ∈ S 1 .

H2: For all j ∈ {1, . . . , n}, the function λ j (θ) + λ ′′ j (θ) is not identically zero. The second assumption guarantees that the curve (36) associated with λ j does not degenerate to a single point. Since the two-dimensional case n = 2 is completely treated in Section 4.1, we can assume that n ≥ 3 and we have the following result:

Proposition 6.4 Assume that A ∈ M n (C) satisfies H1, H2 above. If n ≥ 3, the numerical density f A : C → R + is of class C n-3 .
We have already seen that both hypotheses H1, H2 are necessary, in general, for the conclusion of Proposition 6.4 to hold. For instance, if A ∈ M 3 (C) is a normal matrix whose numerical range W (A) has nonempty interior, the numerical density f A is proportional to the characteristic function of W (A) and is therefore discontinuous on ∂W (A). A more subtle example is provided by the matrix A 3 defined in (33): if |a| 2 + |b| 2 = 4, we have here λ 1 (θ) = -1, λ 2 (θ) = 0, λ 3 (θ) = 1 for all θ ∈ S 1 . Thus H1 is satisfied, but obviously not H2, and the explicit formula (34) shows that the numerical density of A 3 is discontinuous at the origin.

Proof of Proposition 6.4. Assume that n ≥ 3 and that A ∈ M n (C) satisfies H1, H2. If Q = Q(∂ x , ∂ y ) is a homogeneous differential operator of degree n -3, we have as in (45):

(Qf A )(x + iy) = 1 4π S 1 Q(cos θ, sin θ) HB (n-2) [λ 1 (θ), . . . , λ n (θ)](x cos θ + y sin θ) dθ , (47)
where

B (n-2) [λ 1 (θ), . . . , λ n (θ)](s) = -(n -1)! n j=1 H(λ j (θ) -s) k =j (λ k (θ) -λ j (θ)) . (48)
Here H : R → [0, 1] denotes the Heaviside function. Taking the Hilbert transform of (48) with respect to s and using (47), we thus find

(Qf A )(z) = (n-1)! 4π 2 n j=1 S 1 Q(cos θ, sin θ) k =j (λ k (θ) -λ j (θ)) log |f j (θ, z)| dθ ≡ n j=1 S 1 h j (θ) log |f j (θ, z)| dθ , (49)
where f j (θ, z) = λ j (θ) -Re(ze -iθ ) and h j : S 1 → R is a smooth function. Note that the integrand in (49) is 2π-periodic, because by assumption H1 this is the case for all eigenvalues λ 1 (θ), . . . , λ n (θ). It remains to show that each integral in the right-hand side of (49) defines a continuous function of z ∈ C. Replacing A by A -zI n (an operation which does not affect the properties H1, H2), we see that it is sufficient to prove continuity at z = 0. This in turn is obvious if λ j (θ) = f j (θ, 0) does not vanish, so from now on we focus on the case where λ j (θ) has (isolated) zeros on S 1 . Using a partition of unity, we can treat each zero separately, so it is sufficient to consider the case where λ j (θ) has a single zero of order q ≥ 1 at θ = 0, and h j (θ) is localized near the origin. By analyticity, for z close to zero we have the factorization

f j (θ, z) = g(θ, z) q p=1 (θ -µ p (z)) ,
where the (possibly complex) roots µ p (z) depend continuously on z, with µ p (0) = 0, and g(•, z) does not vanish in a neighborhood of zero. The quantity we have to study is therefore

q p=1 R h j (θ) log |θ -µ p (z)| dθ + R h j (θ) log |g(θ, z)| dθ .
The last term is clearly a continuous function of z. In the integral involving µ p , we make the change of variables θ = t + Re µ p (z) and observe that log |t| ≤ log |ti Im µ p (z)| ≤ 0 in a neighborhood of (t, z) = (0, 0), hence continuity with respect to z follows from Lebesgue's dominated convergence theorem.

Once continuity of the derivatives of order n -3 has been established, we can obtain further regularity results for the numerical density by using the representation formula (41) or (43). In view of Lemma 2.3, it is sufficient to study the density in a neighborhood of a point z ∈ Σ A . Since Σ A = C A = C 1 ∪ . . . ∪ C k by assumption H1, there exists j ∈ {1, . . . , k} such that z ∈ C j , and for simplicity we assume that z / ∈ C p for all p = j and that the curve C j has no self-intersection at z. This means that the function θ → f j (θ, z) has a zero of order m ≥ 2 at some point θ ∈ S 1 , and only simple zeros for θ = θ mod π; moreover, if p = j, θ → f p (θ, z) has only simple zeros and does not vanish for θ = θ. Without loss of generality, we assume from now on that z = 0 and θ = 0, and we first consider the simplest case where m = 2. This means that λ j (θ) = α 2 θ 2 + O(θ 3 ) near θ = 0, for some α = 0. If z = x+iy is sufficiently small, it follows that the analytic function θ → f j (θ) = λ j (θ)-(x cos θ+y sin θ) has exactly two zeros θ ± (z) in a neighborhood of the origin, which satisfy

θ ± (z) = 1 α y ± sign(α) y 2 + 2αx + O(|x| + |y| 2 ) .
The critical curve C j , which is the set of all points z for which θ + (z) = θ -(z), is therefore given by the equation x = -1 2α y 2 + O(y 3 ) in a small ball B around the origin. Moreover B \ C j = B r ∪ B c , where B r is the set of all z ∈ B for which the roots θ ± (z) are real and distinct, whereas z ∈ B c when θ ± (z) are complex conjugate with nonzero imaginary part. As is easily verified, the local center of curvature of C j is located on the side of B c . Now, let P = P (∂ x , ∂ y ) be a homogeneous differential operator of degree n -2, and consider the expression of (Pf A )(z) given by (43). Assume that, in the right-hand side, the parameter ǫ > 0 is chosen in such a way that the slit annulus z) for all z ∈ B c , but that A 2ǫ does not contain any other root of the determinant ∆(w, z) for z ∈ B (these conditions are easily achieved by choosing first ǫ and then B sufficiently small). Under these assumptions, the right-hand side of (43) defines a smooth function of z ∈ B, which coincides with (Pf A )(z) if z ∈ B r but not if z ∈ B c . Indeed, in the latter case, we have to consider in addition the roots w ± (z) of ∆(w, z) which are not taken into account by the fixed integration contours in (43) since |1 -|w ± (z)|| < ǫ. Using Cauchy's theorem, we easily obtain

A ǫ = {w ∈ C | 0 < |1 -|w|| < ǫ} contains the points w ± (z) = e 2iθ ± (
(Pf A )(z) = (Pf A ) reg (z) + (n-1)! 4π Q(w -(z), z) + Q(w + (z), z) w -(z) -w + (z) , z ∈ B c ,
where (Pf A ) reg (z) denotes the regular part of (Pf A )(z), given by the right-hand side of (43) with fixed ǫ, and

Q(w, z) = P (w)(w -w + (z))(w -w -(z)) ∆(w, z) .
In particular, we have for all z ∈ B c

(Pf A )(z) -(Pf A ) reg (z) ≤ C |w + (z) -w -(z)| ≤ C |θ + (z) -θ -(z)| ≤ C dist(z, C j ) 1/2 ,
and this estimate is sharp if P (1) = P (1, 0) = 0, because in that case Q(w, z) does not vanish near w = 1 if z ∈ B. Summarizing, we have reached the following important conclusion: Near a regular point z ∈ C j of the critical set Σ A of a generic matrix A ∈ M n (C), the derivatives of order n -2 of the numerical density f A are smooth on the convex side of the curve C j , and blow up like dist(z, C j ) -1/2 on the concave side. If n = 2, this is in full agreement with the explicit formula [START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF] obtained in Section 4.1. If n ≥ 3, we deduce after integration that the numerical density f A is of class C n-5/2 in a neighborhood of such a point z.

Using the same techniques, it is also possible to study the regularity of the numerical density near more singular points z ∈ Σ A . On the typical example represented in Fig. 2, we see that the following two cases have to be analyzed: i) Crossings, which occur when two critical curves C j and C p intersect transversally at z, or when a curve C j has a self-intersection at z. This situation can be treated exactly as before, except that one has to consider four distinct regions near z, instead of two.

ii) Cusps, which arise whenever one of the functions λ j (θ) + λ ′′ j (θ) has a simple zero. Here we can repeat the analysis above, assuming that λ j (θ) = α 3 θ 3 + O(θ 4 ) near θ = 0. In a neighborhood of z = 0, we find that f j (θ, z) has either three real roots (for z inside a cuspidal domain with tip at z), or one real and two complex conjugate roots (outside the cusp). Using the same argument as before, we conclude that Pf A is smooth inside the cusp, but blows up on the other side of the critical curve C j . Altogether, the numerical density f A is of class C n-8/3 near the cusp, see [START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF]Section 4.3] for a similar analysis of the singularities of the fundamental solution of (4).

Under generic assumptions on the matrix A, all intersections are transversal and all cusps are non degenerate, so that the singular set C A is a generic curve in the sense of real algebraic geometry, and the singularities of the numerical density f A can be completely analyzed using the techniques described bove. In particular, the derivatives of order n -2 of the numerical density f A are locally integrable, and since we know from Proposition 6.4 that f A ∈ C n-3 we conclude that the relation (41) holds everywhere (in the sense of distributions), and not only on the complement of Σ A .

Connexion with the fundamental solution

As was explained in the introduction, the numerical measure of a matrix A ∈ M n (C) is related to the fundamental solution of the hyperbolic system (4). This connexion can be established rigorously by comparing the expression (41) for the derivatives of the numerical density f A with the representation formulas for the fundamental solution E(t, x) of ( 4), which can be found e.g. in [START_REF] Bazer | Lacunas of the Riemann matrix of symmetric-hyperbolic systems in two space variables[END_REF][START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF].

In what follows, we identify C with R 2 , and we denote by x = (x 1 , x 2 ) the points of the Euclidean plane. With a slight abuse of notation, we write f A (x) instead of f A (z), and we consider as subsets of R 2 the various regions associated with A, such as W (A) or Σ A . To derive a representation formula for the fundamental solution E(t, x), we take the Radon transform of (5) with respect to x ∈ R 2 , and obtain the one-dimensional hyperbolic system

∂ t Ẽ(t, s, θ) + H(θ)∂ s Ẽ(t, s, θ) = I n δ t=0 ⊗ δ s=0 ,
where H(θ) is given by (2) and Ẽ(t, s, θ) denotes the Radon transform of E(t, x). Using the method of characteristics, we easily find Ẽ(t, s, θ) = n j=1 δ(stλ j (θ))P j (θ) , t ≥ 0, where λ 1 (θ), . . . , λ n (θ) are the eigenvalues of H(θ) and P 1 (θ), . . . , P n (θ) the corresponding spectral projections. Now, if we invert the Radon transform as in [START_REF] Petrowsky | On the diffusion of waves and lacunas for hyperbolic equations[END_REF] and use the identity

n j=1 1 s -tλ j (θ) P j (θ) = (sI n -tH(θ)) -1 ,
we arrive at the representation formula

E(t, x) = - 1 4π 2 f.p. S 1 (x 1 cos θ + x 2 sin θ)I n -tH(θ) -2 dθ , t ≥ 0 , (50)
which coincides with Eq. (4.4a) in [START_REF] Bazer | Lacunas of the Riemann matrix of symmetric-hyperbolic systems in two space variables[END_REF]. Arguing as in Section 6.1, one can show that equation ( 50) is rigorously satisfied for all (t, x) ∈ (0, ∞) × R 2 with x t ∈ R 2 \ Σ A , and that E(t, x) is smooth in that region of space-time.

To compare the numerical density with the fundamental solution, it is natural to extend f A to a homogeneous function of space and time by setting

F A (t, x) = t n-3 f A x t , t > 0 , x ∈ R 2 . ( 51
)
We then have the following result: Proposition 6.5 For any A ∈ M n (C), there exists a matrix-valued homogeneous polynomial Q of degree n -1 such that

E(t, x) = Q(∂ t , ∂ x 1 , ∂ x 2 )F A (t, x) , (52)
for all (t, x) ∈ (0, ∞) × R 2 with x t ∈ R 2 \ Σ A .
In particular, if Ω is a connected component of the region Π A ⊂ R 2 defined by (46), we know from Corollary 6.2 that f A (x) is polynomial of degree at most n -3 in Ω, and (51) then shows that F A (t, x) is also polynomial of degree at most n -3 in the half-cone C + (Ω) = {(t, x) ∈ (0, ∞) × R 2 | x t ∈ Ω}. By Proposition 6.5, we conclude that E(t, x) = 0 in C + (Ω). We thus have Corollary 6.6 The fundamental solution E(t, x) of a matrix A ∈ M n (C) vanishes for all (t, x) ∈ (0, ∞) × R 2 such that x t belongs to the region Π A ⊂ R 2 defined by (46).

In the language of partial differential equations, the domain of influence of the origin for system (4) is the half-cone

C + (D) = {(t, x) ∈ (0, ∞) × R 2 | x t ∈ D}, where D ⊂ R 2 is the complement of the largest connected open region on which E * = E(1, •) vanishes. If in addition E * = 0 in
some open region L ⊂ D, we say that L is a lacuna of the hyperbolic system [START_REF] Bazer | Lacunas of the Riemann matrix of symmetric-hyperbolic systems in two space variables[END_REF]. With this terminology, Corollary 6.6 asserts that each connected component of Π A either lies outside the domain of influence of the origin, or is a lacuna. This important geometric characterization of lacunas is originally due to Petrowsky [START_REF] Petrowsky | On the diffusion of waves and lacunas for hyperbolic equations[END_REF], and was thoroughly discussed in [START_REF] Bazer | Lacunas of the Riemann matrix of symmetric-hyperbolic systems in two space variables[END_REF][START_REF] Yen | The Riemann matrix of (2 + 1)-dimensional symmetric-hyperbolic systems[END_REF] and from a more algebraic point of view in [START_REF] Atiyah | Lacunas for hyperbolic differential operators with constant coefficients[END_REF][START_REF] Atiyah | Lacunas for hyperbolic differential operators with constant coefficients[END_REF]. Strictly speaking, it gives only a sufficient condition for the occurence of lacunas, but further work allows to show that all stable lacunas satisfy this criterion. In other words, in generic situations, the converse of Corollary 6.6 also holds: each open region on which E * vanishes (in particular, any lacuna) belongs to Π A .

As a typical example, consider the matrix A ∈ M 3 (C) defined by (55), whose numerical density is represented in Fig. 1. The polynomial region Π A has just two components here: the exterior of D = W (A), where both f A and E * vanish, and the interior of the cuspidal triangle, where E * = 0 and f A is identically constant, which is therefore a lacuna. A more complicated situation is depicted in Fig. 2, where the polynomial region Π A has six connected components, among which five correspond to lacunas. In a different spirit, it is also interesting to consider the nongeneric matrix (57) whose numerical density is studied in Section 7. The domain of influence of the origin is the union of the closed unit disk and a single point {a}, so the numerical range W (A) = conv(D) is substantially larger than D if |a| > 1, see Fig. 5. In that case, the region Π A has again two components, none of which is a lacuna: the exterior of W (A), and the interior of the triangular region W (A)\D where f A is identically constant. We see on these examples that the numerical density allows to distinguish between various regions where the fundamental solution vanishes identically, and which are nevertheless of rather different nature.

Before proving Proposition 6.5, we briefly verify its conclusion on a simple example. If A ∈ M 2 (C) is defined by (31), system (4) reduces to

∂ t u 1 + ∂ x 1 u 2 -i∂ x 2 u 2 = 0 , ∂ t u 2 + ∂ x 1 u 1 + i∂ x 2 u 1 = 0 .
Combining both equations, one verifies that ∂ 2 t u j = ∆u j for j = 1, 2, hence the fundamental solution E(t, x) can easily be computed using Poisson's formula for the solution of the wave equation in two dimensions. In agreement with Proposition 6.5, the result is:

E(t, x) = ∂ t -∂ x 1 + i∂ x 2 -∂ x 1 -i∂ x 2 ∂ t F A (t, x) , |x| < t ,
where according to (32), (51)

F A (t, x) = 1 t f A x t = 1 2π 1 t 2 -|x| 2 1 {|x|<t} .
Proof of Proposition 6.5. If n = 1 the conclusion is trivial, because both E * and f A vanish identically outside Σ A (which is reduced to a single point), so we assume henceforth that n ≥ 2.

If P = P (∂ t , ∂ x 1 , ∂ x 2
) is a homogeneous differential operator of degree n -2, then using (51) it is straightforward to verify that

(PF A )(t, x) = 1 t P - x t • ∇ ξ , ∂ ξ 1 , ∂ ξ 2 f A ξ= x t , (53)
whenever x t ∈ R 2 \ Σ A .
Here ξ ∈ R 2 denotes the argument of the function f A , which has to be replaced by x t after differentiation. We warn the reader that equality (53) holds only if P is of degree n -2. Applying Proposition 6.1, we deduce that

P (∂ t , ∂ x 1 , ∂ x 2 )F A (t, x) = - (n-1)! 4π 2 t f.p. S 1 P (-x t • e θ , cos θ, sin θ) det(H(θ) -x t • e θ I n ) dθ , (54)
where e θ = (cos θ, sin θ). On the other hand, starting from (50), we observe that

- x t • e θ I n -H(θ) -2 = ∂ ∂s (s + x t • e θ )I n -H(θ) -1 s=0
.

By Cramer's rule, the inverse of the matrix SI n -H(θ) = SI n -A 1 cos θ -A 2 sin θ, with S = s + x t • e θ , has the following form

SI n -H(θ) -1 = 1 ∆(θ, S)
S P 0 (S, cos θ, sin θ) + cos θ P 1 (S, cos θ, sin θ) + sin θ P 2 (S, cos θ, sin θ) , where ∆(θ, S) = det(H(θ)-SI n ) and P 0 , P 1 , P 2 are matrix-valued homogeneous polynomials of degree n -2. The idea is now to insert this expansion into the right-hand side of (50) and to use (54) to express the result as a derivative of order n -1 of the numerical density F A . We begin with the term involving P 1 , and remark that cos θ ∂ ∂s

P 1 (s + x t • e θ , cos θ, sin θ) ∆(θ, s + x t • e θ ) s=0 = t ∂ ∂x 1 P 1 ( x t • e θ , cos θ, sin θ) ∆(θ, x t • e θ )
.

The corresponding contribution to (50) is thus

E 1 (t, x) = 1 4π 2 t ∂ ∂x 1 f.p. S 1 P 1 ( x t • e θ , cos θ, sin θ) ∆(θ, x t • e θ ) dθ = - 1 (n-1)! ∂ ∂x 1 P 1 (-∂ t , ∂ x 1 , ∂ x 2 )F A (t, x) .
Similarly, the term involving P 2 gives the contribution

E 2 (t, x) = - 1 (n-1)! ∂ ∂x 2 P 2 (-∂ t , ∂ x 1 , ∂ x 2 )F A (t, x) .
Finally, to treat the expression containing P 0 , we observe that ∂ ∂s

(s + x t • e θ )P 0 (s + x t • e θ , cos θ, sin θ) ∆(θ, s + x t • e θ ) s=0 = (1 + x • ∇ x ) P 0 ( x t • e θ , cos θ, sin θ) ∆(θ, x t • e θ )
, and using ( 53), ( 54) we obtain the contribution

E 0 (t, x) = - 1 (n-1)! t (1 + x • ∇ x )P 0 (-∂ t , ∂ x 1 , ∂ x 2 )F A (t, x) = 1 (n-1)! ∂ ∂t P 0 (-∂ t , ∂ x 1 , ∂ x 2 )F A (t, x) .
Recalling that E(t, x) = E 0 (t, x) + E 1 (t, x) + E 2 (t, x), we arrive at (52).

Three-dimensional examples

To illustrate the results of the previous sections, we consider here four concrete examples which, according to [19, Section I.7], give a complete picture of what can happen for three-dimensional matrices. The two-dimensional case, which is much simpler, was already treated in Section 4.1, and the references [START_REF] Atiyah | Lacunas for hyperbolic differential operators with constant coefficients[END_REF][START_REF] Bazer | Lacunas of the Riemann matrix of symmetric-hyperbolic systems in two space variables[END_REF][START_REF] Jonckheere | Differential topology of numerical range[END_REF] include a detailed study of the singular set Σ A for a few higher-dimensional examples.

Example 1. We first consider the 3 × 3 matrix

A =   -1.5 1 0 -1 1 1 0 -1 0.5   , (55)
which is generic in the sense that hypotheses H1, H2 in Section 6.2 are fulfilled. The eigenvalues of the associated Hermitian matrix H(θ) satisfy λ 1 (θ) < λ 2 (θ) < λ 3 (θ) for all θ ∈ [0, π], as can be seen from Fig. 3 for which the critical set Σ A can be computed exactly. Indeed, if H(θ) is the Hermitian matrix (2), it is easy to verify that det(λI

(right). The critical set Σ

A = C A = C 1 ∪ C 2 is
A =   0 1 1 0 0 1 0 0 0   , (56)
3 -H(θ)) = λ 3 -3 4 λ -1 4 cos θ, hence λ 1 (θ) = cos θ + 2π 3 , λ 2 (θ) = cos θ -2π 3 , λ 3 (θ) = cos θ 3 .
Thus the permutation (35) is just a cycle τ = (1 2 3), and applying (37) we easily find that the critical curve C A is the cardioid defined by

C A = { 1 3 (2e iφ + e 2iφ ) | φ ∈ S 1 }. Since λ 1 (0) = λ 2 (0), the bitangent set C ′
A is not empty and consists of the line segment joining the points -1/2 ± i/(2 √ 3), see Fig.

(left). Altogether we have Σ

A = C A ∪ C ′
A , and we observe that Σ A encloses a convex region of the complex plane which is of course the numerical range W (A). The index N (z) defined by (39) is equal to 3 outside C A and to 1 inside. The numerical density f A vanishes outside W (A) and is equal to a nonzero constant inside the cuspidal region, in agreement with Propositions 2.6 and 6.1. In particular, f A is discontinuous along the line segment C ′ A .

- Example 3. The matrices considered so far were unitarily irreducible. In contrast, the matrix

A =   0 2 0 0 0 0 0 0 a   , where a ∈ C , (57)
is the direct orthogonal sum of the two-dimensional Jordan block (31) and the one-dimensional matrix (a). The numerical density f A can therefore be computed using Proposition 2.7, and without loss of generality we can assume that a ≥ 0. However, we have to distinguish between three cases: i) If 0 ≤ a < 1, the numerical range W (A) is the closed unit disk, and the numerical density has the following expression:

f A (z) = 1 π √ 1 -a 2 argch 1 -az 1 (1 -az 1 ) 2 -(1 -|z| 2 )(1 -a 2 ) , |z| < 1 ,
which reduces to (34) when a = 0. In particular f A vanishes on the unit circle, has a logarithmic singularity at the point {a}, and is otherwise smooth. The singular set Σ A is the union of the unit circle and the point {a}.

ii) In the limiting case a = 1, the numerical range is still the closed unit disk, but the formula

f A (z) = 1 π 1 -|z| 2 1 -z 1 , |z| < 1 ,
shows that the numerical density has now an algebraic singularity at the boundary point z = 1.

iii) When a > 1, the numerical range W (A) is the convex hull of the union of the unit disk and the exterior point {a}. Within this region, the numerical density satisfies

f A (z) = 1 π √ a 2 -1 arccos 1 -az 1 (1 -az 1 ) 2 + (1 -|z| 2 )(a 2 -1) , when |z| < 1 , and f A (z) = 1/ √ a 2 -1 when |z| > 1.
As is easily verified, the eigenvalues of the Hermitian matrix H(θ) are λ 1 (θ) = -1, λ 2 (θ) = 1, and λ 3 (θ) = a cos θ, see Fig. 5 (right). The algebraic curve C A consists of the unit circle (associated with λ 1 , λ 2 ) and the point {a} (corresponding to λ 3 ), but the bitangent set C ′ A is not empty and consists of two line segments, see Fig. 5 (left). Example 4. As a final example, we consider the case of a normal matrix A ∈ M 3 (C) whose eigenvalues λ 1 , λ 2 , λ 3 are not colinear. Then the numerical range is the triangle with vertices {λ 1 , λ 2 , λ 3 } and the numerical density is a multiple of the characteristic function of W (A). In that situation Σ A = C A ∪C ′ A , where C A is the set of all vertices and C ′ A the set of all edges of the triangle.

Corollary 8.2 Assume that A n ∈ M n (C) and that A n 2 /n → 0 as n → ∞. Then the measure µ Anδ µ An converges weakly to zero as n → ∞. We recall that the numerical radius of a matrix A ∈ M n (C) is defined by

R(A) = sup{|z| | z ∈ W (A)} = sup{| Ax, x | | x ∈ ∂B n } ,
and satisfies R(A) ≤ A ≤ 2R(A) [START_REF] Horn | Topics in matrix analysis[END_REF][START_REF] Serre | Matrices. Theory and applications[END_REF]. Thus, a sequence of matrices A n ∈ M n (C) is uniformly bounded (in the operator norm) if and only if the numerical ranges W (A n ) are all contained in a bounded region of the complex plane. Under this assumption, Proposition 8.1 shows that the variance of µ An is O(1/n) as n → ∞, so that the numerical measure is asymptotically concentrated in a disk of radius O(1/ √ n) around the mean µ An .

In the introduction, we have observed that the numerical measure µ An is the distribution of the random variable A n X n , X n ∈ C when the vector X n is uniformly distributed on the unit sphere ∂B n . With this interpretation, Corollary 8.2 is reminiscent of the weak law of large numbers in probability theory. Under slightly stronger assumptions, it is also possible to obtain a pointwise convergence result in the spirit of the strong law of large numbers. Without loss of generality, we assume from now on that Tr(A n ) = 0 for all n ≥ 1, so that the measure µ An is centered at the origin.

Proposition 8.3 Assume that the matrices

A n ∈ M n (C) satisfy Tr(A n ) = 0 for all n ≥ 1 and sup n≥1 (log n) A n n 1/2 < ∞ . (59)
If for each n ≥ 1 the random variable X n is uniformy distributed on the unit sphere ∂B n , then A n X n , X n converges almost surely to zero as n → ∞.

Proof. It is clearly sufficient to prove the result for Hermitian matrices A n , because the general case then follows by considering the real and imaginary parts of A n X n , X n . We thus assume that A n = A * n for all n ≥ 1, and we denote by λ n,1 , . . . , λ n,n the eigenvalues of A n . For each n ≥ 1, let Y n,1 , . . . , Y n,n be independent and identically distributed complex random variables with density function f Y (z) = π -1 e -|z| 2 , z ∈ C. In particular, we have E(|Y n,m | 2k ) = k! for each k ∈ N. Since the Euclidean measure on ∂B n is the projection on the unit sphere of the standard Gaussian measure in C n , we obtain a uniformly distributed random variable on ∂B n by setting

X n = U n Y n / Y n , where Y n = (Y n,1 , . . . , Y n,n ) ⊤ and U n ∈ U n (C) is a unitary matrix such that U * n A n U n = diag(λ n,1 , . . . , λ n,n ). Thus A n X n , X n = λ n,1 |Y n,1 | 2 + • • • + λ n,n |Y n,n | 2 |Y n,1 | 2 + • • • + |Y n,n | 2 = P n Q n ,
where

P n = 1 n n m=1 λ n,m |Y n,m | 2 = 1 n n m=1 λ n,m (|Y n,m | 2 -1) , Q n = 1 n n m=1 |Y n,m | 2 .
By the strong law of large numbers, the denominator Q n converges almost surely to 1 as n → ∞, hence it remains to show that the numerator P n converges almost surely to zero. But this follows from classical theorems on the limiting behavior of weighted sums of independent random variables, see [START_REF] Chow | Limiting behavior of weighted sums of independent random variables[END_REF][START_REF] Teicher | Almost certain convergence in double arrays[END_REF]. Since for each n ≥ 1 the random variables

Central limit theorems

The results established so far show that, for a sequence of traceless matrices A n ∈ M n (C), the numerical measure µ An tends to concentrate on the origin as n → ∞. Under stronger assumptions, we now prove that the rescaled measure µ √ nAn converges to a Gaussian distribution, as in the classical central limit theorem. We first consider the Hermitian case, which is somewhat simpler. Proposition 8.4 Let {A n } n≥1 be a sequence of Hermitian matrices such that A n ∈ M n (C) and Tr(A n ) = 0. We assume that

(log n) A n n 1/2 ---→ n→∞ 0 , and
1 n Tr(A 2 n ) ---→ n→∞ σ 2 > 0 . ( 60
)
Then the rescaled numerical measure µ √ nAn converges weakly to the normal distribution N (0, σ 2 ) as n → ∞.

Proof. We use the same notations as in the proof of Proposition 8.3. For each n ≥ 1, the numerical measure µ An is the distribution of the random variable

Z n = λ n,1 |Y n,1 | 2 + • • • + λ n,n |Y n,n | 2 |Y n,1 | 2 + • • • + |Y n,n | 2 = P n Q n ,
where λ n,1 , . . . , λ n,n denote the eigenvalues of A n and Y n,1 , . . . , Y n,n are independent complex random variables with density function f Y (z) = π -1 e -|z| 2 . Since Q n converges almost surely to 1 as n → ∞, we have to show that n 1/2 P n converges in law to N (0, σ 2 ).

To do that, we use the Lindeberg-Feller theorem for triangular arrays of random variables [8, Section 2. Invoking the Lindeberg-Feller theorem, we conclude that n 1/2 P n = X n,1 + • • • + X n,n converges in law to N (0, σ 2 ), which is the desired result.

Remark. Since the numerical density of a Hermitian matrix is a B-spline, Proposition 8.4 shows under very general assumptions that B-splines of degree n satisfy a central limit theorem in the limit n → ∞. In the particular case of uniform B-splines, this result was obtained by Unser et. al. in [START_REF] Unser | On the Asymptotic Convergence of B-Spline Wavelets to Gabor Functions[END_REF].

Before considering more general matrices, we would like to mention an alternative proof of Proposition 8.4 which has its own interest. The starting point is a very nice formula for the moments of the numerical measure of a Hermitian matrix A ∈ M n (C). Fix k ∈ N and let λ 1 , . . . , λ n ∈ R denote the eigenvalues of A. Using [START_REF] Joswig | On the numerical range map[END_REF] or [START_REF] Kato | Perturbation theory for linear operators[END_REF] with φ(x) = x k , we find

R x k dµ A (x) = (n -1)! D n-1 (t 1 λ 1 + • • • + t n λ n ) k dt 1 . . . dt n-1 ,
where D n-1 is the (n-1)-dimensional simplex defined in [START_REF] Helgason | The Radon transform, 2nd edition[END_REF] and t n = 1-(t 1 +• • •+t n-1 ). To evaluate the right-hand side, we apply the multinomial formula

(X 1 + • • • + X n ) k = |α|=k k! α! X α ,
which shows that the k-th moment of the numerical measure µ A is the complete symmetric homogeneous polynomial of degree k in the variables λ 1 , . . . , λ n , divided by the combinatorial factor n+k-1 k which is just the number of terms in the sum. Using the Newton identities, the right-hand side of (61) can be decomposed into as a sum of products of elementary symmetric polynomials of the form p ℓ = λ ℓ 1 + • • • + λ ℓ n , see [START_REF] Macdonald | Symmetric functions and Hall polynomials[END_REF]Eq. (2.14')]. If we assume that Tr(A) = 0, then p 1 = 0 and the number of nonzero terms in the sum is considerably reduced. Using these remarks, it is not difficult to show that, under the assumptions of Proposition 8.3, the k-th moment of the rescaled numerical measure µ √ nAn satisfies

R x k dµ √ nAn (x) ---→ n→∞ 0 if k is odd , 2 -k/2 σ k k! (k/2)! if k is even .
Since Tr(A * n A n ) = n -1 and Tr(A 2 n ) = 0, the assumptions of Proposition 8.5 are satisfied with a = 1 and b = 0. Thus the numerical measure of √ nA n converges weakly to the normal distribution π -1 e -|z| 2 dz as n → ∞. In this example, the measure µ An has in fact a radially symmetric density for all n ≥ 1, see Section 4.2.

Perspectives

As a conclusion, we briefly mention a natural extension of our work, which is left for future investigation. Recall that a homogeneous polynomial P ∈ R[X 0 , X 1 , . . . , X d ] of total degree n is hyperbolic in the direction e 0 = (1, 0, . . . , 0) if, on the one hand, it has partial degree n with respect to the first variable X 0 , and on the other hand the n roots of the univariate polynomial t → P (t, y) are real for every vector y ∈ R d . Hyperbolic polynomials arise as principal symbols of hyperbolic differential operators of order n in d space variables, see [START_REF] Gårding | Linear hyperbolic partial differential equations with constant coefficients[END_REF]. As an example, if A 1 , . . . , A d ∈ M n (C) are Hermitian matrices, the polynomial P (X 0 , X 1 , . . . , X d ) = det(X 0 I n -

X 1 A 1 -• • • -X d A d ) (64)
is hyperbolic. In the particular case where d = 2, it has been conjectured in [START_REF] Lax | Differential equations, difference equations and matrix theory[END_REF], and proved in [START_REF] Helton | Linear matrix inequality representation of sets[END_REF], that all monic hyperbolic polynomials are of the form (64). This is no longer true if d ≥ 3.

It might be argued that a large part of our work is not really about matrices, but rather concerns hyperbolic polynomials in 2 + 1 variables. Indeed, if A ∈ M n (C) and A 1 , A 2 are as in (2), the eigenvalues λ 1 (θ), . . . , λ n (θ) of the Hermitian matrix H(θ) are the solutions of the equation P A (λ, cos θ, sin θ) = 0 , θ ∈ S 1 , where P A (X 0 , X 1 , X 2 ) = det(X 0 I n -X 1 A 1 -X 2 A 2 ) is the hyperbolic polynomial associated with A 1 , A 2 . As was shown in Section 4, the numerical measure µ A is entirely determined by the eigenvalues λ j (θ), hence by the polynomial P A . This in turn suggests a natural way to associate with any hyperbolic polynomial P of degree n in d + 1 variables a probability measure µ P on R d , which coincides with the numerical measure µ A when d = 2 and P = P A . Given a unit vector ω ∈ S d-1 ⊂ R d , let λ 1 (ω), . . . , λ n (ω) ∈ R be the roots of the polynomial equation P (λ, ω) = 0, and let B ω (s) = B[λ 1 (ω), . . . , λ n (ω)](s) be the normalized B-spline with knots λ 1 (ω), . . . , λ n (ω). The "numerical measure" of P is then the unique probability measure µ P on R d whose Radon transform satisfies In this generalized setting, the counterpart of normal matrices is the case where the polynomial P split into linear factors

P (X) = n k=1 (X 0 -v k • (X 1 , . . . , X d )) , with v k ∈ R d .
When n ≥ d + 1 and the vectors v k span the affine space R d , the density of µ P with respect to the Lebesgue measure is the multivariate B-spline in d variables, whose nodes are the v k 's. It is piecewise polynomial of degree nd -1. In the generic situation where any (d + 1)-uplet of vectors v k is an affine basis, it is of class C n-d-2 , see [START_REF] Dahmen | On multivariate B-splines[END_REF]. Again the density is log-concave in this case.

As in the two-dimensional case, the measure µ P can be expressed in terms of B ω using the backprojection method. The example above suggests that, as the space dimension increases, the measure µ P becomes more singular. Then, the inversion formula has to be understood in the sense of distributions. In the three-dimensional case d = 3, we arrive at the simple expression

µ P = - 1 8π 2 ∆ x S 2 B ω (x • ω) dσ(ω) , (65)
where dσ denotes the Euclidean measure on the unit sphere S 2 . As an example, in the particular situation where P = X 2 0 -X 2 1 -X 2 2 -X 2 3 , which corresponds to the differential operator ∂ 2 t -∆ of the wave equation, we obtain µ P = 1 4π dσ(x). This shows that, when d ≥ 3, the support of µ P does not need to be convex. Because the wave equation satisfies the Huyghens Principle, this example suggests that the link between polynomial regions of the density and lacunas of differential operators persists in higher dimensions.

Figure 1 :

 1 Figure 1: The numerical density f A is represented for a typical matrix A ∈ M 3 (R), given by (55) below. In

Figure 2 :

 2 Figure 2: The singular set Σ A (left) and the eigenvalues λ 1 (θ), . . . , λ 5 (θ) (right) are represented for a typical matrix A ∈ M 5 (R). The singular set consists of three closed curves, one of which (the boundary of W (A)) is smooth, and the other two have cusps. In this example, the set Π A defined in (46) has six connected components.

Figure 3 :

 3 Figure 3: The critical set Σ A (left) and the eigenvalues λ j (θ) (right) are represented for the matrix (55).

Figure 4 :

 4 Figure 4: The critical set Σ A (left) and the eigenvalues λ j (θ) (right) are represented for the matrix (56).

Figure 5 :

 5 Figure 5: The critical set Σ A (left) and the eigenvalues λ j (θ) (right) are represented for the matrix (57) with a = 2.

  X n,m = |Y n,m | 2 -1 (1 ≤ m ≤ n) are independent, have zero mean and finite second order moment, and since by (59) the coefficients a n,m = n -1 λ n,m satisfy max 1≤m≤n |a n,m | ≤ C n 1/2 log n , the results of [24, Section 3] imply that P n = n m=1 a n,m X n,m converges almost surely to zero as n → ∞.

  4.b]. Let X n,m = 1 √ n λ n,m (|Y n,m | 2 -1) , 1 ≤ m ≤ n , so that n 1/2 P n = X n,1 + • • • + X n,n .For each fixed n ≥ 1, the random variables X n,m are independent and satisfy E(X n,m ) = 0 for m = 1, . . . , n. Moreover, any ǫ > 0, we haveE(|X n,m | 2 ; |X n,m | ≥ ǫ) = λ 2 n,mπnDn,m,ǫ (|z| 2 -1) 2 e -|z| 2 dz , where D n,m,ǫ = {z ∈ C | λ 2 n,m (|z| 2 -1) 2 ≥ nǫ 2 }. Thus, using the first assumption in (60), we obtain by a direct calculation n m=1 E(|X n,m | 2 ; |X n,m | ≥ ǫ) ≤ n sup 1≤m≤n E(|X n,m | 2 ; |X n,m | ≥ ǫ) ---→ n→∞ 0 .

  where the sum runs over all multi-indices α ∈ N n of order |α| = α 1 + • • • + α n = k. Here we use the standard notationsX α = X α 1 1 • • • X αn n and α! = (α 1 !) • • • (α n !). Now, it is not difficult to verify thatD n-1 t α dt 1 . . . dt n-1 = α! (n + k -1)! ,for any α ∈ N n with |α| = k. We thus obtain the following identityR x k dµ A (x) = k! (n -1)! (n + k -1)! |α|=k λ α ,

(

  Rµ P )(ω, ds) = B ω (s) ds ,where, by definition, (Rµ P )(ω,I) = µ P ({x ∈ R d | x • ω ∈ I}) for any interval I ⊂ R.

  an algebraic curve of degree 6 consisting of a smooth ovate curve C 1 enclosing a cuspidal triangle C 2 (Fig.3, left). The component C 1 corresponds to the eigenvalues λ 1 (θ), λ 3 (θ) while C 2 is associated with λ 2 (θ). All multiplicities are equal to one, and the number N (z) defined by (39) is equal to 3 outside C 1 and inside C 2 , and to 1 in the intermediate region. The numerical density f A is continuous, identically constant inside C 2 , and vanishes outside C 1 . Moreover f A is Hölder continuous with exponent 1/2 across C 1 and C 2 , except at the cusps. The level lines of f A are represented in Fig.1.

	Example 2. We next consider a nongeneric matrix

The geometry of the singular setWe know from Lemma 2.3 that the numerical density of a matrix A ∈ M n (C) is smooth outside the set Σ A of all critical values of the numerical map Φ A . In this section, we describe a few geometric properties of the singular set Σ A which will be needed to formulate our main results in Section

We do

Statistical properties of the numerical measure

In this section, we study the numerical measure from a statistical point of view, and we establish various convergence results which show that the measure µ A of a large matrix A ∈ M n (C) is concentrated in a neighborhood of the barycenter of the spectrum σ(A).

Concentration phenomena

Proposition 8.1 For any A ∈ M n (C), the first moment of the probability measure µ A is the normalized trace of A:

and the variance of µ A is given by

Proof. Let a jk denote the entries of A. Applying definition (1) with φ(z) = z, we have to compute the average of jk a jk x j x k over the unit sphere. By symmetry, the average of x j x k is equal to zero if j = k and to 1/n if j = k. Thus

To compute the second moment, we take φ(z) = |z| 2 and proceed in exactly the same way. By symmetry, the average of x j x k x ℓ x m is zero unless j = k and ℓ = m, or j = m and ℓ = k. Moreover, it is easy to verify that the average of

This gives the desired result, since Var(µ

Now we consider a sequence of matrices

where

As a consequence, if we suppose that A n 2 = o(n) as n → ∞, it follows from (58) that the variance of µ An converges to zero as n → ∞. This gives: Since the moments in the right-hand side are those of the normal law, we conclude that µ √ nAn converges weakly to N (0, σ 2 ) as n → ∞.

We now consider general matrices A n ∈ M n (C), and obtain a central limit theorem by applying Proposition 8.4 to the Radon transform of A n . Proposition 8.5 Let {A n } n≥1 be a sequence of matrices satisfying A n ∈ M n (C), Tr(A n ) = 0, and (log n) A n /n 1/2 → 0 as n → ∞. We assume that

Then the rescaled numerical measure µ √ nAn converges weakly to the Gaussian measure f ∞ (z) dz as n → ∞, where

Proof. For each θ ∈ S 1 , the Hermitian matrices H n (θ) = 1 2 (e -iθ A n + e iθ A * n ) satisfy the assumptions of Proposition 8.4, with

Let us denote by µ n and µ n,θ the numerical measures of √ n A n and √ n H n (θ), respectively. Since µ n,θ is the Radon transform of µ n , the two-dimensional Fourier transform μn (ξ) for ξ = re iθ is precisely the one-dimensional Fourier transform of µ n,θ evaluated at r [START_REF] Helgason | The Radon transform, 2nd edition[END_REF]. Thus, applying Proposition 8.4, we find

for any r ≥ 0, θ ∈ S 1 . This shows that µ n converges weakly as n → ∞ to the measure µ ∞ on C defined by μ∞ (ξ) = e -1 2 |ξ| 2 σ(θ) 2 = e -a+b 4 Re(ξ) 2 e -a-b 4 Im(ξ) 2 , ξ ∈ C . Inverting the Fourier transform, this gives dµ ∞ = f ∞ (z) dz with f ∞ as in (63).

Remarks.

1. If we assume for simplicity that A n is uniformly bounded, we can suppose (up to extracting a subsequence) that 1 n Tr(A * n A n ) converges as n → ∞ to some a ≥ 0. However, we have to assume in (62) that a > 0 in order to get a universal Gaussian limit.

2.

Similarly, the first assumption in (62) implies that 1 n Tr(A 2 n ) converges, after extracting a subsequence, to some b ∈ C with |b| ≤ a. Multiplying A n by a unit complex number, we can assume that 0 ≤ b ≤ a, but in the borderline case where b = a the limiting measure f ∞ (z) dz should be replaced (2πa) -1/2 e -x 2 /(2a) dx ⊗ δ y=0 .

Example. Let us consider the Jordan block of size n:

.