Oleg Szehr
email: oleg.szehr@posteo.de

Rachid Zarouf
email: rachid.zarouf@univ-amu.fr

Optimal estimates for condition numbers and norms of resolvents in terms of the spectrum

Keywords: Condition number, Resolvent, Toeplitz matrix, Model matrix, Blaschke product 2010 Mathematics Subject Classification: Primary: 15A60; Secondary:30D55 Contents

est destinée au dépôt

I. INTRODUCTION

Let M n be the set of complex n × n matrices and let ||T || denote the spectral norm of T ∈ M n . We denote by σ = σ(T ) the spectrum of T and by m T its minimal polynomial. We denote by |m T | the degree of m T . In this article we are interested in estimates of the type

||R(ζ, T )|| ≤ Φ(|m T |, σ, ζ) (1)
where R(ζ, T ) = (ζ -T ) -1 denotes the resolvent of T at point ζ ∈ C -σ and Φ is a function of |m T |, σ and ζ. We present a general method that yields optimal estimates of the above type and allows us to pinpoint the "worst" matrices for this problem. As it turns out for any ζ ∈ σ among these matrices are so-called analytic Toeplitz matrices. Generally, Toeplitz matrices are characterized by the existence of a sequence of complex numbers a = (a k ) k=n-1 k=-n+1 with

T =       a 0 a -1 . . a -n+1 a 1 . . . . . . . . . . . . . a -1 a n-1 . . a 1 a 0       .
We call T a an analytic Toeplitz matrix if a k = 0 for all k < 0 and we will denote by T n ⊂ M n the set of Toeplitz matrices and by T a n ⊂ T n the analytic Toeplitz matrices. We will think of an estimate of the above type as an extremal problem. That is, under certain constraints on the set of admissible matrices T and fixed ζ, we are looking to maximize ||R(ζ, T )|| over this set. To guarantee a finite bound we at least need a constraint on the norm and the spectrum of possible T . In our discussion we will generally impose that ||T || ≤ 1 and call such T a contraction. Note that this condition can always be achieved by normalization T / ||T ||. C n ⊂ M n denotes the set of contractions. To ensure that the resolvent is finite, ζ must be separated from σ. Depending on the situation it will be convenient to measure this separation in Euclidean distance (2)

The upper bound can be seen using polar decomposition [START_REF] Nikolski | Condition numbers of large matrices and analytic capacities[END_REF]p. 137]. To prove the reverse inequality, one can use A. Horn's theorem (see [AH] and also [START_REF] Marshall | Inequalities: theory of majorization and its applications[END_REF]Chap. 9,Sect. E]) that gives a converse to the famous Weyl inequalities [START_REF] Nikolski | Condition numbers of large matrices and analytic capacities[END_REF]p. 137].

The case |ζ| = 1 was studied by Davies and Simon and applied to characterize the localization of zeros of random orthogonal polynomials on the unit circle [DS]. The corresponding estimate for |ζ| = 1 is

sup {d(ζ, σ) ||R(ζ, T )|| : T ∈ C n } = cot( π 4n
).

(3)

Note that |ζ| = 1 implies that we have p(ζ, σ) = 1 for any σ such that we could add this condition to the optimization problem. In this paper we provide an approach, which 1) allows us to treat the problems (2) and (3) simultaneously and generalize them to any |ζ| ≤ 1,

2) strengthens the result (3) for |ζ| = 1, by considering the sup over T ∈ C n with fixed minimal polynomial and

3) provides explicit analytic Toeplitz matrices that achieve the estimates (1), ( 2) and (3).

Already for ζ = 0 it is a non-trivial question which matrices achieve the estimate (2). A theoretic classification of such matrices was given in [START_REF] Nikolski | Condition numbers of large matrices and analytic capacities[END_REF]. However, it seems rather difficult to compute explicit matrix representations from this classification. It was suspected in [RZ] that the estimate (2) is not reached by any Toeplitz matrix; this claim is refuted here.

II.

RESOLVENT BOUNDS

Before we state our main results we need to introduce some further notation. We denote by D = {z ∈ C : |z| < 1} the open unit disc in the complex plane, D is its closure and ∂D its boundary. For any finite set σ ⊂ D we define

s(ζ, σ) := max 1 -|λ| 2 |1 -λζ| : λ ∈ σ, p(ζ, λ) = p(ζ, σ) .
The function s is related to the Stolz angle [START_REF] Conway | Functions of one complex variable II[END_REF]p. 23]

between ζ ∈ D -σ and λ ∈ σ. For T ∈ C n the function ζ → s(ζ, σ(T )) is bounded on D by 1 + ρ(T ) ≤ 2, where ρ(T ) := max{|λ| : λ ∈ σ(T )}
is the spectral radius of T . Of particular importance to our discussion will be the analytic Toeplitz matrix X r, β given entry-wise by

(X r, β ) ij = 0 if i < j r n-1 if i = j βr n-(i-j+1) if i > j (4) with r ∈ [0, 1] and β ∈ [0, 2].
The spectral norm of X r, β is computed below in Proposition II.6.

T * =          λ 0 . . . . . . 0 1 -λ 2 λ . . . . . . . -λ(1 -λ 2 ) 1 -λ 2 λ . . . . . . . . . . . . . . . . . . 0 (-λ) n-2 (1 -λ 2 ) . . . -λ(1 -λ 2 ) 1 -λ 2 λ         
, is a contraction of minimal polynomial m = (z -λ) n that achieves (5).

Remarks:

1. Setting ζ = 0 in Theorem II.1 we recover Kronecker's result for condition numbers. Moreover, an analytic Toeplitz matrix is extremal for the Kronecker maximization problem (2). This answers the questions raised in [RZ].

2. Considering Theorem II.1 for ζ ∈ ∂D, and bounding

X 1, s(ζ,σ(T )) ≤ X 1, 1+ρ(T ) ≤ ||X 1, 2 || ,
we recover the resolvent estimate of [DS] as ||X 1, 2 || = cot( π 4n ) (cf. Proposition II.6 below). However, if e.g. ρ(T ) is known we can improve on [DS].

3. Equality is achieved by a so-called model matrix corresponding to a fully degenerate spectrum. This matrix arises as a representation of the compressed backward shift operator of H 2 with respect to the Malmquist-Walsh basis, see Section III below for details. The model matrix is a natural generalization of a simple Jordan block of size n to the situation when the spectrum is σ = {a} ⊂ D.

The next proposition treats the case σ(T ) ⊂ ∂D.

Proposition II.2. Let T ∈ C n with spectrum σ ⊂ ∂D. Then for any ζ ∈ C -σ ||R(ζ, T )|| ≤ 1 d(ζ, σ) . ( 6
)
It follows that if T is a contraction such that σ(T ) ⊂ ∂D then it satisfies a Linear Resolvent Growth with constant 1, see [BN]. Clearly, any diagonal matrix in C n achieves the inequality (6).

Theorem II.1 has the disadvantage that the upper bound depends on the slightly unyielding quantity ||X r,β ||. Taking together Theorem II.4 with Proposition II.6 (both below) leads directly to the following.

||R(ζ, T )|| ≤ 1 d(1, σζ) 1 r |m| (1 -r|ζ|) .
The estimate is sharp in the limit of large |m|.

B. Constrained maximization of the resolvent

For any given ζ ∈ D and r ∈ [0, 1] we study the quantity

R(ζ, r) := sup {d(1, σζ) ||R(ζ, T )|| : T ∈ C n , p(ζ, σ) = r} , (7)
where σ = { λ : λ ∈ σ}. R also depends on n but we keep the notation simple and do not write this explicitly.

Theorem II.4. For any ζ ∈ D and r ∈ (0, 1) we have that

R(ζ, r) = 1 r n ||X r,βmax || ,
where β max = 1-r 2 1-r|ζ| . The estimate is achieved by an analytic Toeplitz matrix and in the limit of large matrices we find

R(ζ, r) ∼ 1 r n 1 1 -r |ζ| as n → ∞.
Recall that Kronecker's result (2) is

R(0, r) = 1 r n .
Davies and Simon (3) have shown that

R(ζ, 1) = cot( π 4n ), ζ ∈ ∂D -{1},
which can be seen in the above theorem by continuously moving ζ to ∂D. Note, however, the fundamental difference in the large n behavior of the resolvent for ζ ∈ D and ζ ∈ ∂D. This is a consequence of the fact that for ζ inside the unit disk ||X r,β || is bounded from above, see Proposition II.6.

Finally for ζ ∈ ∂D and a polynomial P σ with zero set σ so that p(ζ, σ) = 1 we can compute the constant

K(ζ, P σ ) = sup {d(ζ, σ) ||R(ζ, T )|| : T ∈ C n , m T = P σ } , (8)
which is bounded by R(ζ, r).

Proposition II.5. Let P σ be a polynomial with zero set σ = σ 1 ∪ σ 2 , σ 1 ⊂ D and σ 2 ⊂ ∂D so that P σ =P σ 1 P σ 2 . For any ζ ∈ ∂D -σ, we have

K(ζ, P σ ) ≤ max(1, X 1, s(ζ, σ 1 ) ) ≤ ||X 1, 2 || = cot( π 4|P σ 1 | ), (9)
where

X 1, s(ζ, σ) , X 1, 2 ∈ M |Pσ 1 | (C). Moreover, fixing ζ ∈ ∂D-σ, n 1 , n 2 ≥ 1 and 0 ≤ ρ 1 < 1 and taking the supremum of K(ζ, P σ ) over all sets σ = σ 1 ∪ σ 2 with n 1 = |P σ 1 | and n 2 = |P σ 2 | such that max λ∈σ 1 |λ| ≤ ρ 1 we get: sup K(ζ, P σ ) = ||X 1, 1+ρ 1 || (10)
where X 1, 1+ρ 1 ∈ M n 1 (C). In particular, the supremum of the above quantity taken over

0 ≤ ρ 1 < 1 is equal to cot( π 4n 1 ). C. Computing ||X r, β ||
In this subsection we study the quantity ||X r, β || for r ∈ (0, 1] and β ∈ [0, 2]. We provide an implicit formula for ||X r, β ||.

Proposition II.6. Let X r, β with r ∈ (0, 1] and β ∈ [0, 2] be the n × n matrix introduced in (4). The following assertions hold.

i) If β = 1 -r 2 then ||X r, β || = 1.
ii) At fixed n, ||X r, β || grows monotonically with r and β.

iii) At fixed r and 1 -r 2 ≤ β, ||X r, β || grows monotonically with n. iv) Let U n denote the n-th degree Chebychev polynomial of second kind, i.e. U n (cos(θ)) = sin ((n+1)θ) sin (θ)

. If β = 1 -r 2 then ||X r, β || = |λ * |
, where λ * is the largest in magnitude solution of the polynomial equation

µ n+1 rU n γ 2µ + µ n (λ 2 -r 2n-2 (β -1) 2 )U n-1 γ 2µ = 0 (11) with µ = λ 2 + r 2n-2 (β -1) and γ = -λ 2 (r + 1/r) + r 2n-2 (r + 1 r (β -1) 2 ).
v) If r ∈ (0, 1) and 1 -r 2 ≤ β we have the limit

||X r, β || → β 1 -r 2 , n → ∞.
vi) For any n, r ∈ (0, 1), β max ∈ (0, 2] and 1 -r 2 ≤ β ≤ β max we have

||X r, β || ≤ β max 1 -r 2 .
The equations in iv) (see ( 11)) are discussed in Lemma II.7 below. For instance we can recover known values of ||X 1, β || from Proposition II.6, iv). In this case

||X 1, β || = 1 2 (β -2) 2 + β 2 cot(θ * /2) ,
where θ * is the unique solution of ( 11) in [ 2n-1 2n π, π) and it follows, see [OS], that

||X 1, 0 || = 1, ||X 1, 1 || = 1 2 sin( π 4n+2 ) , ||X 1, 2 || = cot( π 4n
).

Lemma II.7. Let r ∈ (0, 1] and β ∈ (0, 2] such that β -1 + r 2 = 0. We consider the equation

cot(nθ) + (2 -β)r r 2 + (β -1) 1 sin(θ) + r 2 -(β -1) r 2 + (β -1) cot(θ) = 0. ( 12
)
1. If r = 1, then (12) has exactly 2n distinct solutions in [-π, π).

If

(a) if 1 < β < 1 + r, then (12) has exactly 2n -2 distinct solutions in [-π, π) if and only if n > β-1+r 2 (1-r)(r+β-1) . If n ≤ β-1+r 2
(1-r)(r+β-1) , then (12) has exactly 2n distinct solutions in [-π, π).

(b) If β > 1+r, and n ≥ 2, then (12) has exactly 2n-4 distinct solutions in [-π, π) if and only if n > β-1+r 2 (1-r)(r+β-1) . If β-1+r 2 (1+r)(β-1-r) < n ≤ β-1+r 2 (1-r)(r+β-1) , then (12) has exactly 2n -2 distinct solutions in [-π, π). Finally, if n ≤ β-1+r 2 (1+r)(β-1-r) , then (12) has exactly 2n distinct solutions in [-π, π).

If β < 1, and r /

∈ {1, 1 -β}, we distinguish two subcases:

(a) if 1 -r 2 < β < 1, then (12) has exactly 2n -2 distinct solutions in [-π, π) if and only if n > β-1+r 2 (1-r)(r+β-1) , (b) if 0 < β < 1 -r 2 , we distinguish two subcases: i. if 0 < β < 1 -r, then (12) has exactly 2n -4 distinct solutions in [-π, π) if and only if n > β-1+r 2 (1-r)(r+β-1) . If n ≤ β-1+r 2 (1-r)(r+β-1) , then (12) has exactly 2n distinct solutions in [-π, π). ii. If 1 -r < β < 1 -r 2 , then (12) has exactly 2n -2 distinct solutions on [-π, π) for any n ≥ 2. 5. If β = 1 then (12) has exactly 2n -2 distinct solutions in [-π, π) if and only if n > r 1-r . If n ≤ r 1-r , then (12) has exactly 2n distinct solutions in [-π, π).

III. MODEL SPACES AND MODEL OPERATORS

This section lays down the theoretical framework on which our results are footed. We refer to [START_REF] Nikolski | Treatise on the shift operator[END_REF] for an in-depth study of the topic. To the minimal polynomial m = m T of T ∈ M n we associate a Blaschke product

B(z) := i z -λ i 1 -λi z .
The product is taken over all i such that the corresponding linear factor z -λ i occurs in the minimal polynomial m. Thus, the numerator of B as defined here is exactly the associated minimal polynomial. The space of holomorphic functions on D is denoted by Hol(D). The Hardy spaces considered here are

H 2 := f ∈ Hol(D) : ||f || 2 H 2 := sup 0≤r<1 1 2π 2π 0 |f (re iφ )| 2 dφ < ∞ , and
H ∞ := f ∈ Hol(D) : ||f || H∞ := sup z∈D |f (z)| < ∞ .
Let B be the Blaschke product associated to the minimal polynomial of T ∈ M n . We define the |m|-dimensional model space

K B := H 2 ⊖ BH 2 := H 2 ∩ (BH 2 ) ⊥ ,
where we employ the usual scalar product from L 2 (∂D). The model operator M B acts on K B as

M B : K B → K B f → M B (f ) := P B (zf ),
where P B denotes the orthogonal projection on K B . In other words, M B is the compression of the multiplication operation by z to the model space K B . Multiplication by z has operator norm 1 so that ||M B || ≤ 1. Moreover, it is not hard to show that the eigenvalues of M B are exactly the zeros of the corresponding Blaschke product, see [START_REF] Nikolski | Operators, functions, and systems: an easy reading[END_REF].

For our subsequent discussion we will assume that T can be diagonalized. This does not result in any difficulties since the upper bounds obtained in the special case extend by continuity to bounds for non-diagonalizable T , see [OS]. One natural orthonormal basis for K B (for diagonalizable T ) is the Malmquist-Walsh basis {e k } k=1,...,|m| with ([NN3, p. 137])

e k (z) := (1 -|λ k | 2 ) 1/2 1 -λk z k-1 i=1 z -λ i 1 -λi z with e 1 (z) = (1 -|λ 1 | 2 ) 1/2 1 -λ1 z .
IV. PROOFS A. Proofs of Theorem II.1 and Proposition II.2

We recall [START_REF] Szehr | Eigenvalue estimates for the resolvent of a non-normal matrix[END_REF]Theorem III.2] that given a finite Blaschke product B = N i=1 z-ν i 1-ν i z , of degree N ≥ 1, with zeros ν i ∈ D, the components of the resolvent of the model operator M B at any point ζ ∈ C -{ν 1 , ..., ν N } with respect to the Malmquist-Walsh basis are given by

ζ -M B -1 1≤i, j≤N = 1 B(ζ)      0 if i < j 1 1-ν i ζ s =i ζ-νs 1-νsζ if i = j (1-|ν i | 2 ) 1/2 1-ν i ζ (1-|ν j | 2 ) 1/2 1-ν j ζ s / ∈{j, ...i} ζ-νs 1-νsζ if i > j . (13)
We begin by proving Proposition II.2.

Proof of Proposition II.2. We consider T ∈ C n and assume for a moment that σ(T ) ⊂ D.

It follows from [START_REF] Szehr | Eigenvalue estimates for the resolvent of a non-normal matrix[END_REF]Theorem III.2] that for any ζ ∈ C -σ(T ) the resolvent of T is bounded by

||R(ζ, T )|| ≤ ||R(ζ, M B )|| ,
where B is the Blaschke product associated with m T . We conclude from (13) that

ζ -M B -1 1≤i, j≤|m| =      0 if i < j 1 ζ-λ i if i = j (1-|λ i | 2 ) 1/2 1-λ i ζ (1-|λ j | 2 ) 1/2 1-λ j ζ i k=j 1-λk ζ ζ-λ k if i > j .
Taking the limit |λ i | → 1 we have that for ζ ∈ C -σ(T ) the norm of the resolvent is bounded by the norm of the diagonal matrix

   1 ζ-λ 1 . . . 1 ζ-λn    , which is 1 d(ζ, σ) .
Proof of Theorem II.1. We begin by showing that (5) holds. As before we suppose σ(T ) ⊂ D. The more general case σ(T ) ⊂ D follows immediately by continuously moving the eigenvalue towards ∂D. As before we have that

||R(ζ, T )|| ≤ ||R(ζ, M B )|| and that ζ -M B -1 1≤i, j≤|m| =        0 if i < j 1 1-λ i ζ 1 b λ i if i = j (1-|λ i | 2 ) 1/2 1-λ i ζ (1-|λ j | 2 ) 1/2 1-λ j ζ 1 i k=j b λ k if i > j , where we abbreviated b λ i = ζ-λ i 1-λ i ζ . Recall that for any n × n matrices A = (a ij ) and A = (a ′ ij ), the condition |a ij | ≤ a ′ ij implies that ||A|| ≤ ||A ′ ||. Hence, we estimate ζ -M B -1 i, j ≤        0 if i < j 1 |1-λ i ζ| 1 |b λ i | if i = j (1-|b λ i | 2 ) 1 2 (1-|b λ j | 2 ) 1 2 1-|ζ| 2 1 i k=j |b λ k | if i > j ,
where we made use of the well-known [JG] fact that

1 -|λ| 2 |1 -λζ| 2 = 1 -|b λ (ζ)| 2 1 -|ζ| 2 , λ = ζ.
By assumption min i |b λ i | = p(ζ, σ) = r and for any l with |b λ l | = r we have

ζ -M B -1 i, j ≤        0 if i < j 1 |1-λ i ζ| 1 r if i = j 1-|b λ l (ζ)| 2 1-|ζ| 2 1 r i-j+1 if i > j ,
and since

1 -|b λ l (ζ)| 2 1 -|ζ| 2 = 1 -|λ l | 2 |1 -λ l ζ| 2 ≤ 1 d(1, σζ) 1 -|λ l | 2 |1 -λ l ζ| ≤ β d(1, σζ)
,

with β = s(ζ, σ) we conclude that ζ -M B -1 i, j ≤ 1 d(1, σζ) 1 r |m|    0 if i < j r |m|-1 if i = j βr |m|-(i-j+1) if i > j .
Next we discuss the sharpness of (5). We consider ζ ∈ [-1, 1] and the model operator

M B with B = b n λ and λ ∈ (-1, 1). In this case d(1, σζ) = |1 -λζ| and p(ζ, σ) = |b λ (ζ)| and s(ζ, σ) = 1-|λ| 2 |1-λζ| . We find that ζ -M B -1 i, j = 1 d(1, σζ) 1 B(ζ)      0 if i < j (b λ (ζ)) n-1 if i = j 1-|λ| 2 |1-λζ| (b λ (ζ)) n-(i-j+1) if i > j . If b λ (ζ) > 0 then the equality ||R(ζ, M B )|| = 1 d(1, σζ) 1 r n ||X r, β || is clear. If b λ (ζ) <
R n (ζ, r) = R n (|ζ| , r). Proof. Indeed, writing ζ = |ζ| exp(iθ) we have 1 -λζ = 1 -exp(iθ) λ |ζ| and ||R(ζ, T )|| = (|ζ| -exp(-iθ)T ) -1 , which entails that d(1, σ(T )ζ) ||R(ζ, T )|| = d(1, σ(T ′ ) |ζ|) R(|ζ| , T ′ ) with T ′ = exp(-iθ)T . Clearly T ∈ C n ⇔ T ′ ∈ C n and p(ζ, σ(T )) = p(|ζ| , σ(T ′ )). The same argument implies T ∈ T a n ⇔ T ′ ∈ T a n
The lemma says that for our optimization task it is sufficient to consider real ζ. Proof of Proposition II.5. We first prove (9). We consider T ∈ C n with spectrum σ and assume for a moment that σ 2 ⊂ D. It follows from [START_REF] Szehr | Eigenvalue estimates for the resolvent of a non-normal matrix[END_REF]Theorem III.2] that for any ζ ∈ C -σ the resolvent of T is bounded by

||R(ζ, T )|| ≤ ||R(ζ, M B )|| ,
where B is the Blaschke product associated with P σ . We conclude from (13) that the matrix of the resolvent of the model operator M B at any point ζ ∈ C -σ (with respect to the Malmquist-Walsh basis) is given by a block-diagonal matrix

A 1 0 0 A 2 ,
where

A 1 ∈ M |Pσ 1 | (C) and A 2 ∈ M |Pσ 2 | (C)
are respectly defined by (taking now the limit as |µ| → 1 for µ ∈ σ 2 ):

(A 1 ) 1≤i, j≤|Pσ 1 | =      0 if i < j 1 ζ-λ i if i = j (1-|λ i | 2 ) 1/2 1-λ i ζ (1-|λ j | 2 ) 1/2 1-λ j ζ i k=j 1-λk ζ ζ-λ k if i > j , and
A 2 =     1 ζ-µ 1 . . . 1 ζ-µn 2     ,
where µ i are the unimodular eigenvalues of M B . To conclude, it remains to remark that the spectral norm of R(ζ, T ) is nothing but

max(||A 1 || , ||A 2 ||), that ||A 2 || = 1 d(ζ, σ 2 )
, while the inequality

||A 1 || ≤ 1 d(ζ, σ 1 ) X 1, s(ζ,σ 1 ) ,
follows directly from Theorem II.1. This proves the first inequality in (9). The second one is due to the fact that s(ζ, σ 1 ) is bounded from above by 2 and the equality refers to the resolvent estimate in [DS] or more generally to our Proposition II.6. Now, let us prove (10). We consider

T ∈ C n with spectrum σ = σ 1 ∪ σ 2 , n 1 = |P σ 1 | and n 2 = |P σ 2 | such that max λ∈σ 1 |λ| ≤ ρ 1 . We fix ζ ∈ ∂D -σ.
The last inequality, combined with the trivial observation that s(ζ, σ 1 ) ≤ ρ 1 , proves that

||A 1 || ≤ ||X 1, ρ 1 || .
To complete the proof of (10), we first remark that by rotation invariance (as in the proof of Lemma IV.1) the sup defined in the statement of (10) does not depend on ζ and coincides with its value at point ζ = 1. We consider the spectrum σ = σ 1 ∪ σ 2 defined by setting P σ 1 = (z -ρ 1 ) n 1 and P σ 2 = (z + 1) n 2 . Let B be the Blaschke product associated with P σ . We repeat the proof of (9) to the model operator M B . The matrix of its resolvent at point 1 (with respect to the Malmquist-Walsh basis) is given by the block-diagonal matrix

A 1 0 0 A 2 ,
where

(A 1 ) 1≤i, j≤n 1 = 1 d(1, σ 1 ) 1 B(1)      0 if i < j (b ρ 1 (1)) n-1 if i = j 1-|ρ 1 | 2 |1-ρ 1 | (b ρ 1 (1)) n-(i-j+1) if i > j = X 1, 1+ρ 1 |1 -ρ 1 | .
and A 2 = 1 2 I n 2 , where I n 2 is the identity matrix of M n 2 (C). Thus,

||R(1, M B )|| = max( 1 |1 -ρ 1 | ||X 1, 1+ρ 1 || , 1 2 ), but d(1, σ) = d(1, σ 1 ) = |1 -ρ 1 |
, and as a consequence:

d(1, σ) ||R(1, M B )|| = max(||X 1, 1+ρ 1 || , |1 -ρ 1 | 2 ).
The above maximum is equal to ||X 1, 1+ρ 1 || : indeed, according to Proposition II.6,

||X 1, 1+ρ 1 || ≥ ||X 1, 1 || = 1 2 sin( π 4n 1 +2 )
, and the latter quantity

1 2 sin( π 4n 1 +2
) is always greater than 1 2 (increasing with n 1 ) since n 1 ≥ 1. The last part of the statement relies again on Proposition II.6 and the fact that

ρ → ||X 1, 1+ρ || is increasing on [0, 1].
Proof of Proposition II.6. Our approach is guided by the methods developed in [EE, DS, OS]. Often it will be convenient to transform the Toeplitz matrix X r, β to an associated Hankel matrix X r, β , see [DS, OS]. We set X r, β = X r, β J, where

J =    1 . . . 1    .
X is a Hermitian Hankel matrix that satisfies || X|| = ||X r, β ||. Hence, in our discussion we can focus on the largest in magnitude eigenvalue of X. (We omit indices in X for simpler notation.) i) In the case β = 1 -r 2 we have

1 r n X r, β = 1 r n    0 if i < j r n-1 if i = j (1 -r 2 )r n-(i-j+1) if i > j = M -1 b n -r 1≤i, j≤n ,
where the second equality is (13). It is well known [NN1, Theorem 3.12, p.147], [START_REF] Szehr | Eigenvalue estimates for the resolvent of a non-normal matrix[END_REF]Lemma III.5] that (see the cited articles for a discussion)

M -1 B = ||1/z|| H∞/BH∞ = 1 |B(0)| = λ∈σ 1 |λ| .
In particular, we find

M -1 b n -r = 1 r n and ||X r, β || = 1. ii)
X is an entry-wise non-negative matrix. A standard result in Perron-Froebenius theory [START_REF] Minc | Nonnegative matrices, Wiley-Interscience series in discrete mathematics and optimization[END_REF]Theorem 4.2] asserts that the spectral radius of such matrices is a non-negative eigenvalue and the corresponding eigenvector is entry-wise non-negative1 . Furthermore the spectral radius is monotonically increasing in the matrix-entries [HM, Corollary 2.1], i.e. X ≤ B entry-wise implies ρ( X) ≤ ρ(B).

iii) Here, we show a simple proof for β > 1. The more general discussion under the condition 1 -r 2 ≤ β can be found below. Let Xn ∈ M n and Xn+1 ∈ M iv) We turn to the case r 2 +β -1 = 0. The discussion is guided by [EE]. The eigenvalues of X2 are the eigenvalues of X squared. Hence, we are looking for the largest λ 2 such that

0 = det( X 2 -λ 2 ) = det( X -λ) det( X + λ).
The matrix X 2 -λ 2 can be computed explicitly. We display this matrix for n = 5 for simplicity, see Appendix V for the general matrix. If n = 5 we have

X2 -λ 2 =          -λ 2 + β 2 1-r 2•4 1-r 2 + r 8 β 2 r 1-r 2•3 1-r 2 + βr 7 β 2 r 2 1-r 2•2 1-r 2 + βr 6 β 2 r 3 1-r 2•1 1-r 2 + βr 5 βr 4 β 2 r 1-r 2•3 1-r 2 + βr 7 -λ 2 + β 2 r 2 1-r 2•3 1-r 2 + r 8 β 2 r 3 1-r 2•2 1-r 2 + βr 7 β 2 r 4 1-r 2•1 1-r 2 + βr 6 βr 5 β 2 r 2 1-r 2•2 1-r 2 + βr 6 β 2 r 3 1-r 2•2 1-r 2 + βr 7 -λ 2 + β 2 r 4 1-r 2•2 1-r 2 + r 8 β 2 r 5 1-r 2•1 1-r 2 + βr 7 βr 6 β 2 r 3 1-r 2•1 1-r 2 + βr 5 β 2 r 4 1-r 2•1 1-r 2 + βr 6 β 2 r 5 1-r 2.1 1-r 2 + βr 7 -λ 2 + β 2 r 6 1-r 2•1 1-r 2 + r 8 βr 7 βr 4 βr 5 βr 6 βr 7 -λ 2 + r 8          .
The determinant of the above matrix is transformed by the following steps.

1. We divide the k th column and the k th row of ( X 2 -λ 2 ) by r k-1 .

2. In the determinant resulting from (1), we subtract the k th column from the (k -1) th and leave the n th unchanged.

3. In the determinant resulting from (2), we subtract the (k -1) th row from the k th and leave the n th unchanged. We obtain the determinant of a tri-diagonal n × n matrix.

4. In the determinant resulting from (3), we multiply the k th row by r 2k-1 .

We obtain that det( X2 -λ 2 ) = 0 iff

0 = det          γ α ′ 0 • • • 0 α γ α ′ . . . . . . 0 α . . . . . . 0 . . . . . . . . . γ α ′ 0 • • • 0 α y          = det          γ µ 0 • • • 0 µ γ µ . . . . . . 0 µ . . . . . . 0 . . . . . . . . . γ µ 0 • • • 0 µ y         
,

where α = λ 2 r + r 2n-1 (β -1), α ′ = α r 2 , y = -rλ 2 + r 2n-1 and γ = -λ 2 (r + 1/r)+ r 2n-2 (r + 1 r (β -1) 2
) and we set µ = α/r in the last step. We have that

det          γ µ 0 • • • 0 µ γ µ . . . . . . 0 µ . . . . . . 0 . . . . . . . . . γ µ 0 • • • 0 µ y          = det          γ µ 0 • • • 0 µ γ µ . . . . . . 0 µ . . . . . . 0 . . . . . . . . . γ µ 0 • • • 0 µ γ          + det          γ µ 0 • • • 0 µ γ µ . . . . . . 0 µ . . . . . . 0 . . . . . . . . . γ 0 0 • • • 0 µ y -γ          = det          γ µ 0 • • • 0 µ γ µ . . . . . . 0 µ . . . . . . 0 . . . . . . . . . γ µ 0 • • • 0 µ γ          + (y -γ) det          γ µ 0 • • • 0 µ γ µ . . . . . . 0 µ . . . . . . 0 . . . . . . . . . γ µ 0 • • • 0 µ γ         
, where in the last sum the first determinant is of order n, while the second is of order n -1.

The assumption µ = 0 is equivalent to λ 2 = r 2n-2 (1 -β) and the determinant equation turns into γy = 0. In this case we obtain a solution if β = 0 or r 2 + β -1 = 0 but these choices of β are treated differently. (Note that λ = 1 implies µ = 0.) For the determinant of a tri-diagonal Toeplitz matrix we have the well known formula [MV, EP, EE, OS] det

         a 1 0 • • • 0 1 a 1 . . . . . . 0 1 . . . . . . 0 . . . . . . . . . a 1 0 • • • 0 1 a          = U n (a/2),
where U n denotes the n-th degree Chebychev polynomial of the second kind. Hence the eigenvalues of X correspond to the roots of the polynomial equation

µ n+1 U n γ 2µ + µ n (y -γ)U n-1 γ 2µ = 0. v)
We make use of properties of U n to reveal properties of the roots of equation ( 11). Let us consider λ 2 so that there is real θ with γ/2µ = cos(θ), i.e.

λ 2 = r 2(n-1) (β -1) 2 + r 2 -2r(β -1) cos(θ) 1 + r 2 + 2r cos(θ) . ( 14
) Using U n (cos(θ)) = sin ((n+1)θ) sin (θ)
Equation ( 11) turns into ), which leads to [START_REF] Mason | Chebychev polynomials[END_REF]Section 1.4] [START_REF] Mason | Chebychev polynomials[END_REF]Section 1.4] to find

sin((n + 1)θ) sin(θ) + 1 r λ 2 -r 2n-2 (β -1) 2 λ 2 + r 2n-2 (β -1) sin(nθ) sin(θ) = 0, which plugging back λ 2 is cot(nθ) + 1 tan(θ) r 2 -(β -1) r 2 + β -1 + (2 -β)r r 2 + β -1 1 sin(θ) = 0. (15) Similarly, if γ/2µ ∈ [1, ∞) we set γ/2µ = cosh(θ) with θ ∈ [0, ∞
λ 2 = r 2(n-1) (β -1) 2 + r 2 -2r(β -1) cosh(θ) 1 + r 2 + 2r cosh(θ) (16) with sinh ((n + 1)θ) sinh (nθ) = - (2 -β)r + 2(1 -β) cosh(θ) r 2 + β -1 . (17) Finally, in case that γ/2µ ∈ [-1, -∞) we set γ/2µ = -cosh(θ) with θ ∈ [0, ∞) and recall that U n (-x) = (-1) n U n (x)
λ 2 = r 2(n-1) (β -1) 2 + r 2 + 2r(β -1) cosh(θ) 1 + r 2 -2r cosh(θ) (18) with sinh ((n + 1)θ) sinh (nθ) = (2 -β)r -2(1 -β) cosh(θ) r 2 + β -1 . ( 19
)
We consider the first case. Any solution θ k of ( 15) leads to a solution λ2 k of det( X2λ 2 ) = 0 via ( 14). An important point about Equation ( 14) is that all such solutions |λ k | 2 can be bounded by a function that decays exponentially with n. This allows us to study properties of || X|| by counting solutions of equation ( 15). We observe that

1 sin(θ) → +∞ as θ → 0 + and 1 tan(θ) → +∞ as θ → 0 + 1 sin(θ) → +∞ as θ → π - and 1 tan(θ) → -∞ as θ → π -.
It follows that if 1 -r 2 < β < 1 + r then (see Lemma II.7 for details)

1 tan(θ) r 2 -(β -1) r 2 + β -1 + (2 -β)r r 2 + β -1 1 sin(θ) → +∞ as θ → 0 + 1 tan(θ) r 2 -(β -1) r 2 + β -1 + (2 -β)r r 2 + β -1 1 sin(θ) → +∞ as θ → π -.
As -cot(nθ) has n branches in the interval [0, π] and since

-cot(nθ) → -∞ as θ → 0 + -cot(nθ) → +∞ as θ → π -,
we conclude that Equation (15) has n -1 solutions θ k ∈ (0, π) (if n is large enough, see Lemma II.7). The corresponding eigenvalues λ k of X approach 0 exponentially fast, so that for the largest eigenvalue of X we have λ * → Tr( X). On the other hand direct consideration of the trace shows that at fixed β and r we have

Tr( X) → β 1 -r 2 as n → ∞.
If β > 1 + r we have only n -2 real solutions θ k to Equation (15), see Lemma II.7. We consider Equations ( 16), (17) and Equations ( 18), ( 19). As n gets large sinh ((n+1)θ) sinh (nθ) quickly converges to e θ with θ ∈ [0, ∞). By the intermediate value theorem (17) has a solution for β > 1 + r and sufficiently large n. Any λ 2 as in Equation ( 16) decays exponentially with n. Hence, in case β > 1 + r we obtain n -1 eigenvalues that decay exponentially with n, as before. In the asymptotic limit of large n we find that e (θ * ) = 1/r becomes a solution of ( 19) and λ * in (18) does not decay exponentially with n, which identifies the latter as the largest in magnitude solution for large enough n. Hence, we have that ||X|| → β/(1 -r 2 ) as n → ∞.

iii) (second part) We have already seen that ||X r,β || grows monotonically with n if β > 1. Hence, we consider 1 -r 2 < β < 1. As outlined in v) there are for large enough 2 n, see Lemma II.7, n-1 solutions θ k to Equation ( 15) and the corresponding eigenvalues, see ( 14), can be bounded as

r 2(n-1) 1 - β 1 + r 2 ≤ |λ k | 2 ≤ r 2(n-1) 1 - β 1 -r 2 . ( 20
)
We write X(n) and X(n-1) for the matrix X in n and n -1 dimensions. Let ||•|| F R denote the Frobenius norm. Direct computation shows that

X(n) 2 F R -X(n-1) 2 F R = nr 2(n-1) + (n -1)β 2 r 2(n-2) -(n -1)r 2(n-2) .
On the other hand it is well known that || X(n

) || 2 F R = n i=1 |λ (n) i | 2
, where the λ

(n) i
are the eigenvalues of X(n) . Writing λ * ,(n) for the largest in magnitude eigenvalue of X(n) and making use the mentioned properties and the estimates (20) we find that

|λ * ,(n) | 2 -|λ * ,(n-1) | 2 ≥ nr 2(n-1) + (n -1)β 2 r 2(n-2) -(n -1)r 2(n-2) -(n -1)r 2(n-1) 1 - β 1 -r 2 + (n -2)r 2(n-2) 1 - β 1 + r 2 .
The lower bound is a degree-2 polynomial in β, which is convex if r ∈ (0, 1/2). Hence, the lower bound is an increasing function of β ∈ [1 -r 2 , 1]. As β = 1 -r 2 implies |λ * ,(n) | 2 -|λ * ,(n-1) | 2 = 0 this proves monotonicity in n if r ∈ (0, 1/2). The range for r can be extended by using better bounds in (20). Here we note that at the endpoints of the interval for β we already have established monotonicity for any r. Furthermore the norm is a monotonous, convex functions in β. Any of the branches of -cot(nθ) is convex and the terms

r 2 -(β -1) r 2 + β -1 and (2 -β)r r 2 + β -1 in (15) are decreasing in β ∈ [1 -r 2 , 1]. It follows that the solutions θ (n) k of (15) satisfy θ (n) k (β 0 ) -θ (n) k (β 1 ) ≤ θ (n-1) k (β 0 ) -θ (n-1) k (β 1 )
for β 0 ≤ β 1 . In conclusion the quantity

nr 2(n-1) + (n -1)β 2 r 2(n-2) -(n -1)r 2(n-2) + n-2 k=1 |λ (n-1) k | 2 - n-1 k=1 |λ (n) k | 2 is increasing in β. vi) is a direct consequence of iii) and v).
Proof of Lemma II.7. For fixed β ∈ [0, 2], r ∈ (0, 1] we for now focus our attention on real values of θ and write shortly

ϕ r, β (θ) = r 2 -(β -1) r 2 + (β -1) cot(θ) + (2 -β)r β -1 + r 2 1 sin(θ) such that Equation (12) becomes ϕ r, β (θ) = -cot(nθ), θ ∈ [-π, π).
We start with some elementary observations.

1. We restrict to the interval [0, π) as all functions are symmetric with respect to the point (π, 0).

2. cot(nθ) has poles in [0, π) at θ k = k n π, k = 1, ..., n -1.
3. The functions 1 sin(θ) and cot(θ) have poles at 0 and π, cot(θ) ∼ 1 sin(θ) as θ tends to 0 + and cot(θ) ∼ -1 sin(θ) as θ tends to π -.

4. Looking at the polynomials:

P β (r) = r 2 -(2 -β)r -(β -1) = (r -1)(r -(1 -β)),
and

Q β (r) = r 2 + (2 -β)r -(β -1) = (r + 1)(r -(β -1)),
we have that P β (r) < 0 for r ∈ (1 -β, 1) and P β (r) ≥ 0 elsewhere, while Q β (r) < 0 for r ∈ (-1, β -1) and Q β (r) ≥ 0 elsewhere.

5. Taking (c) and (d) into account, we have

ϕ r, β (θ) ∼ θ→0 + r 2 + (2 -β)r -(β -1) β -1 + r 2 1 sin(θ) = Q β (r) β -1 + r 2 1 sin(θ) , while ϕ r, β (θ) ∼ θ→π -- r 2 -(2 -β)r -(β -1) β -1 + r 2 1 sin(θ) = -P β (r) β -1 + r 2 1 sin(θ)
.

We now prove the points in the lemma.

1. The case r = 1. In this case, the structure of ( 12) is discussed in [OS, Section B].

We now suppose that r ∈ (0, 1).

2. The case r = |β -1| / ∈ {0, 2}.

(a) The subcase β < 1. We have that ϕ r, β (θ) = (-2 + β) (cos(θ) + 1) sin(θ) β , and as a consequence, ϕ r, β < 0 on (0, π), lim θ→0 + ϕ r, β (θ) = -∞, while lim θ→π -ϕ r, β (θ) = 0. Thus, the graph of ϕ r, β intersects n -1 times with the graph of θ → -cot(nθ) on [ π 2n , π). Moreover,

ϕ r, β (θ) -cot(nθ) = (2 -β) sin(θ) sin(nθ) β cos(nθ)(1 -cos(θ)) ,
and in particular, for any θ ∈ (0, π 2n ),

ϕ r, β (θ)
-cot(nθ) > 1, which can be seen by computing derivatives with respect to θ: θ → ϕ r, β (θ) -cot(nθ) is increasing on (0, π 2n ), with

lim θ→0 + ϕ r, β (θ) -cot(nθ) = 2n(2 -β) β ≥ 2n > 1.
The graph of ϕ r, β does not intersect the one of θ → -cot(nθ) on [0, π 2n ). We conclude that in this case (12) has 2n -2 real solutions θ k in [-π, π). From now on, we suppose that r = 1 -β, and r = β -1. Let us first notice that in case β -1 + r 2 = 0, we have

lim θ→0 + ϕ r, β (θ) -cot(nθ) = n(1 + r)(β -1 -r) β -1 + r 2 ,
and

lim θ→π - ϕ r, β (θ) -cot(nθ) = n(1 -r)(r + β -1) β -1 + r 2 .
Below, we repeat the reasonings of (2) for the remaining cases (3), ( 4) and ( 5).

3. The case β > 1, and r / ∈ {1, β -1}. With (e), P β (r) < 0, for any r ∈ [0, 1). In particular, for any r ∈ [0, 1),

lim θ→π - ϕ r, β (θ) = +∞, while lim θ→0 + ϕ r, β (θ) = -∞ if r ∈ [0, β -1) +∞ if r ∈ (β-1, 1) , because of the property of Q β (r) in (d).
(a) The subcase 1 < β < 1 + r. In this case, ϕ r, β > 0 on (0, π) and

lim θ→0 + ϕ r, β (θ) = lim θ→π - ϕ r, β (θ) = +∞.
Moreover, ϕ r, β has a global minimum at some point θ 0 ∈ (0, π) such that ϕ r, β is decreasing on (0, θ 0 ] and increasing on [θ 0 , π). Thus, the curve of ϕ r, β intersects the one of θ → -cot(nθ) exactly n -1 times on (0, (2n-1)π 2n

]. On the interval ( (2n-1)π 2n , π), the curve of ϕ r, β does not intersect the one of

θ → -cot(nθ) if and only if lim θ→π - ϕ r, β (θ) -cot(nθ) > 1, (else, the number of real intersections on [0, π) is n). Indeed as before, θ → ϕ r, β (θ)
-cot(nθ) is (positive, continuous and) decreasing on ( (2n-1)π 2n , π). We conclude that in this case (12) has 2n -2 real solutions θ k in [-π, π) if and only if n > β-1+r 2

(1-r)(r+β-1) .

(b) The subcase β > 1 + r. In this case, ϕ r, β is increasing on (0, π), and lim θ→0 + ϕ r, β (θ) = -∞, while lim θ→π -ϕ r, β (θ) = +∞. Thus, the curve of ϕ r, β intersects the one of θ → -cot(nθ) exactly n -2 times on [ π 2n , (2n-1)π 2n

].

Moreover, on (0, π 2n ) ∪ (

(2n-1)π 2n , π), θ → ϕ r, β (θ)
-cot(nθ) is positive and continuous; it is increasing on (0, π 2n ) and decreasing on ( (2n-1)π 2n , π). Thus, the curve of ϕ r, β does not intersect the one of θ → -cot(nθ) on (0, π 2n ) (resp. on ( (2n-1)π 2n , π))

if and only if

lim θ→0 + ϕ r, β (θ) -cot(nθ) > 1 (resp. lim θ→π - ϕ r, β (θ)
-cot(nθ) > 1). We conclude that in this case (12) has 2n -

4 real solutions θ k in [-π, π) if and only if n > β-1+r 2 (1-r)(r+β-1) , since β-1+r 2 (1+r)(β-1-r) < β-1+r 2 (1-r)(r+β-1) .
4. The case β < 1 and r / ∈ {1, 1 -β}. With (e), Q β (r) ≥ 0, for any r ∈ [0, 1). In particular, for any r ∈ [0, 1),

lim θ→0 + ϕ r, β (θ) = +∞ if β -1 + r 2 > 0 -∞ if β -1 + r 2 < 0 , while lim θ→π - ϕ r, β (θ) =    +∞ if β -1 + r 2 > 0 +∞ if β -1 + r 2 < 0 and r ∈ [0, 1 -β) -∞ if β -1 + r 2 < 0 and r ∈ (1 -β, 1)
, because of the property of P β (r) in (d).

(a) The subcase 1 -r 2 < β < 1. In this case, ϕ r, β > 0 on (0, π), and lim θ→0 + ϕ r, β (θ) = lim θ→π -ϕ r, β (θ) = +∞. Moreover, ϕ r, β has a global minimum at some point θ 0 ∈ (0, π) such that ϕ r, β is decreasing on (0, θ 0 ] and increasing on [θ 0 , π). We conclude as in (3)-(a) that in this case, (12) has 2n -

2 real solutions θ k in [-π, π) if and only if n > β-1+r 2 (1-r)(r+β-1) . (b) The subcase 0 < β < 1 -r 2 .
i. If 0 < β < 1 -r, ϕ r, β is increasing on (0, π), lim θ→0 + ϕ r, β (θ) = -∞, and lim θ→π -ϕ r, β (θ) = +∞. We conclude by a similar reasoning as the one in (3)-(b) that in this case also (12) has 2n -

4 real solutions θ k in [-π, π) if and only if n > β-1+r 2 (1-r)(r+β-1) , since β-1+r 2 (1+r)(β-1-r) < β-1+r 2 (1-r)(r+β-1) . ii. If 1 -r < β < 1 -r 2 , ϕ r, β (θ) < 0, ∀θ ∈ (0, π), and lim θ→0 + ϕ r, β (θ) = lim θ→π -ϕ r, β (θ) = -∞.
Moreover, ϕ r, β has a global maximum θ 0 ∈ (0, π) such that ϕ r, β is increasing on (0, θ 0 ] and decreasing on [θ 0 , π). Thus, the curve of ϕ r, β intersects the one of θ → -cot(nθ) exactly n -1 times on [ π 2n , π). in the interval (0, π 2n ), the curve of ϕ r, β does not intersect the one of θ → -cot(nθ) if and only if lim θ→0 + ϕ r, β (θ) -cot(nθ) > 1. The latter condition is always satisfied since n(1+r)(β-1-r) β-1+r 2 = n(1+r)(1-β+r) 1-β-r 2 > n. We conclude that in this case, (12) has 2n -2 real solutions θ k on [-π, π).

5.

The case β = 1. In this case, ϕ r, β (θ) = 1 + rcos(θ) rsin(θ) ,

and thus ϕ r, β (θ) > 0, ∀θ ∈ (0, π), and lim θ→0 + ϕ r, β (θ) = lim θ→π -ϕ r, β (θ) = +∞. Moreover, ϕ r, β has a global minimum at some point θ 0 ∈ (0, π) such that ϕ r, β is decreasing on (0, θ 0 ] and increasing on [θ 0 , π). In particular, the graph of ϕ r, β intersects the one of θ → -cot(nθ) exactly n -1 times on (0, (2n-1)π 2n

]. On the interval ( (2n-1)π 2n , π), the function θ → ϕ r, β (θ)

-cot(nθ) is decreasing and thus the curve of ϕ r, β does not intersect the one of θ → -cot(nθ) if and only if lim θ→π -ϕ r, β (θ) -cot(nθ) > 1. We conclude that in this case, (12) has 2n -2 real solutions θ k on [-π, π) if and only if n > r 1-r (else (12) has 2n real solutions in [-π, π)).

V. APPENDIX

The columns C i , i = 1, . . . , n of the n × n matrix X -λ 2 are:

C 1 =           
-λ 2 + β 2 (1 + r 2 + . . . + r 2.(n-2) ) + r 2.(n-1) β 2 r(1 + r 2 + r 2.(n-3) ) + β r 2.(n-1)-1 β 2 r 2 (1 + r 2 + r 2.(n-4) ) + β r 2.(n-1)-2 . . .

β 2 r n-3 (1 + r 2 ) + β r n+1 β 2 r n-2 + β r n β r n-1            , C 2 =           
β 2 r(1 + r 2 + r 2.(n-3) ) + β r 2.(n-1)-1 -λ 2 + β 2 r 2 (1 + r 2 + . . . + r 2.(n-3) ) + r 2.(n-1) β 2 r 3 (1 + r 2 + r 2.(n-4) ) + β r 2.(n-1)-1 . . .

β 2 r n-2 (1 + r 2 ) + β r n+2 β 2 r n-1 + β r n+1 β r n            , C 3 =              
β 2 r 2 (1 + r 2 + r 2.(n-4) ) + β r 2.(n-1)-2 β 2 r 3 (1 + r 2 + r 2.(n-4) ) + β r 2.(n-1)-1 -λ 2 + β 2 r 4 (1 + r 2 + . . . + r 2.(n-4) ) + r 2.(n-1) β 2 r 5 (1 + r 2 + . . . + r 2.(n-5) ) + β r 2.(n-1)-1 . . .

β 2 r n-1 (1 + r 2 ) + β r n+3 β 2 r n + β r n+2 β r n+1               , C 4 =                
β 2 r 3 (1 + r 2 + r 2.(n-5) ) + β r 2.(n-1)-3 β 2 r 4 (1 + r 2 + r 2.(n-5) ) + β r 2.(n-1)-2 β 2 r 5 (1 + r 2 + r 2.(n-5) ) + β r 2.(n-1)-1 -λ 2 + β 2 r 6 (1 + r 2 + . . . + r 2.(n-5) ) + r 2.(n-1) β 2 r 7 (1 + r 2 + r 2.(n-6) ) + β r 2.(n-1)-1 . . .

β 2 r n (1 + r 2 ) + β r n+4 β 2 r n+1 + β r n+3 β r n+2                 , . . . . . . C n-3 =                 
β 2 r n-4 (1 + r 2 + r 4 ) + β r n+2 β 2 r n-3 (1 + r 2 + r 4 ) + β r n+3 β 2 r n-2 (1 + r 2 + r 4 ) + β r n+4 . . . β 2 r 2n-10 (1 + r 2 + r 4 ) + r 2n-4 β 2 r 2n-9 (1 + r 2 + r 4 ) + r 2n-3 -λ 2 + β 2 r 2n-8 (1 + r 2 + r 4 ) + r 2n-2 β 2 r 2n-7 (1 + r 2 ) + β r 2n-3 β 2 r 2n-6 + β r 2n-4 β r 2n-5

                 , C n-2 =                   
β 2 r n-3 (1 + r 2 ) + β r n+1 β 2 r n-2 (1 + r 2 ) + β r n+2 . . . . . . β 2 r 2n-9 (1 + r 2 ) + β r 2n-5 β 2 r 2n-8 (1 + r 2 ) + β r 2n-4 β 2 r 2n-7 (1 + r 2 ) + β r 2n-3 -λ 2 + β 2 r 2n-6 (1 + r 2 ) + r 2n-2 β 2 r 2n-5 + β r 2n-3 β r 2n-4

                   , C n-1 =                     
β 2 r n-2 + βr n β 2 r n-1 + βr n+1 β 2 r n + βr n+2 . . .

                    
.

  d(z, w) := |z -w| , z, w ∈ C or pseudo-hyperbolic distance p(z, w) := z -w 1 -zw , z, w ∈ C. The case ζ = 0 in (1) corresponds to the study of the condition number CN (T ) = ||T || • T -1 . It was shown by Kronecker in the 19th century that sup T -1 : T ∈ C n , p(0, σ) = r = 1 r n .

A

  . Optimal upper bounds Theorem II.1. Let T ∈ C n with minimal polynomial m and spectrum σ. For any ζ ∈ D-σ, it holds for the resolvent of T that ||R(ζ, T )|| ≤ 1 d(1, σζ) 1 r |m| ||X r, β || , (5) where X r, β is the |m| × |m| analytic Toeplitz matrix defined in (4), r = p(ζ, σ) and β = s(ζ, σ). Moreover, for any ζ ∈ [-1, 1] and λ ∈ (-1, 1), λ = ζ the analytic n × n Toeplitz matrix

  Theorem II.3. Let T ∈ C n with minimal polynomial m and spectrum σ. For any ζ ∈ D-σ let p(ζ, σ) = r denote the pseudo-hyperbolic distance of ζ to σ. Then it holds for the resolvent of T that

  0 equality still holds because the Toeplitz matrices a = (a k ) k≥0 and b = ((-1) k a k ) k≥0 have the same norm. To see this write D = diag(-1, ..., -1) and note that DT a D † = T b . As D is unitary ||T a || = ||T b ||. B. Proofs of Theorem II.4 and Propositions II.5, II.6 Lemma IV.1. For any ζ ∈ D and r ∈ (0, 1)

  Proof of Theorem II.4. We consider a matrix T ∈ C n with p(ζ, σ(T )) = r and the function s in (5). The function λ → 1-|λ| 2 |1-λ|ζ|| is maximized on the pseudo-hyperbolic circle of center |ζ| and radius r at λ max = |ζ|-r 1-r|ζ| and takes the valueβ max = 1-|λmax| 2 1-λmax|ζ| = 1-r 2 1-r|ζ| . Thus using Lemma IV.1 and Theorem II.1 we find that R n (ζ, r) = R n (|ζ|, r) ≤ 1 r n ||X r,βmax ||. According to the proof of Theorem II.1 the model operator corresponding to Blaschke product B = b n λ with λ ∈ (-1, 1), λ = |ζ| satisfies ||R(|ζ|, M B )|| = 1 d(1, σ|ζ|) 1 r n ||X r, β || . We choose λ max = |ζ|-r 1-r|ζ| and have p(|ζ| , σ) = r and s(|ζ| , σ) = β max = 1-r 2 1-r|ζ| . It follows that R n (ζ, r) ≥ 1 r n ||X r, β || and the "worst" matrix for the extremal problem in (7) is the analytic Toeplitz matrix M B with B = b n λmax . The monotonicity properties and the asymptotic behaviour of R n (ζ, r) are consequences of our study of ||X r, βmax ||, which is summarized in Proposition II.6.

  n+1 . We show || Xn || ≤ || Xn+1 ||. Let x * with ||x * || = 1 and Xn x * = || Xn ||x * . Computing entries of vectors and using β ≥ 1 we see that Xn x * ≤ Xn+1 x * 0 .

( b )

 b The subcase β > 1. We have that ϕ r, β (θ) = (-2 + β) (cos(θ) -1) sin(θ) β , and as a consequence, ϕ r, β is increasing on (0, π), ϕ r, β > 0 on (0, π),lim θ→0 + ϕ r, β (θ) = 0, while lim θ→π -ϕ r, β (θ) = +∞. Hence the graph of ϕ r, β intersects n -1 times the one of -cot(nθ) on (0, (2n-1)π 2n ]. Moreover, ϕ r, β (θ) -cot(nθ) = (2 -β) sin(θ) sin(nθ) β cos(nθ)(1 + cos(θ)) ,and in particular, for any θ ∈ ( theorem guarantees that the graph of ϕ r, β does not intersect the one of θ → -cot(nθ) on ( (2n-1)π 2n , π) if and only if 2n(2-β) β) ). We conclude that in this case (12) has 2n -2 real solutions θ k in [-π, π) if and only if n > β 2(2-β) .

  2n-5 + β r 2n-3 -λ 2 + β 2 r 2n-4 + r 2n-

In fact, as can be seen easily by computing powers, X is primitive.

When there are n solutions to (15), monotonicity of ||X r,β || can be readily be seen in (14) and (15).


Amira Essaid

Arnaud Martin

Boutheina Grégory Smits

Ben Yaghlane

Boutheina Ben Yaghlane

A Distance-Based Decision in the Credal Level

Keywords: belief function theory, imprecise decision, distance

come   L'archive ouverte pluridisciplinaire

Introduction

Belief function theory [START_REF] Dempster | Upper and Lower probabilities induced by a multivalued mapping[END_REF][START_REF] Shafer | A mathematical theory of evidence[END_REF] allows us to represent all kinds of ignorance and offers rules for combining several imperfect information provided by different sources in order to get a more coherent one. The combination process helps to make decisions later. Decision making consists in selecting, for a given problem, the most suitable actions to take. Today, we are often confronted with the challenge of making decisions in cases where information is imprecise or even not available. In [START_REF] Smets | The Transferable Belief Model[END_REF], Smets proposed the transferable belief model (TBM) as an interpretation of the theory of belief functions. The TBM emphasizes a distinction between knowledge modeling and decision making. Accordingly, we distinguish the credal level and the pignistic level. In the credal level, knowledge is represented as belief functions and then combined. The pignistic level corresponds to decision making, a stage in which belief functions are transformed into probability functions.

The pignistic probability, the maximum of credibility and the maximum of plausibility are rules that allow a decision on a singleton of the frame of discernment. Sometimes and depending on application domains, it seems to be more convenient to decide on composite hypotheses rather than a simple one. In the literature, there are few works that propose a rule or an approach for making decision on a union of hypotheses [START_REF] Denoeux | Analysis of evidence-theoric decision rules for pattern classification[END_REF][START_REF] Appriou | Approche générique de la gestion de l'incertain dans les processus de fusion multisenseur[END_REF][START_REF] Martin | Decision support with belief functions theory for seabed characterization[END_REF]. Recently, we proposed a decision rule based on a distance measure [START_REF] Essaid | Uncertainty in ontology matching: a decision rule based approach[END_REF]. This rule calculates the distance between a combined mass function and a categorical one. The most likely hypothesis to choose is the hypothesis whose categorical mass function is the nearest to the combined one.

The main topic of this paper is to demonstrate that our proposed decision rule is a particular case of that detailed in [START_REF] Denoeux | Analysis of evidence-theoric decision rules for pattern classification[END_REF] and to extend our rule so that it becomes able to give decisions even with no categorical mass functions. We present also our experiments on mass functions generated randomly as well as on real databases.

The remainder of this paper is organized as follows: in section 2 we recall the basic concepts of belief function theory. Section 3 presents our decision rule based on a distance measure proposed in [START_REF] Essaid | Uncertainty in ontology matching: a decision rule based approach[END_REF]. In section 4, we demonstrate that our proposed rule is a particular case of that proposed in [START_REF] Denoeux | Analysis of evidence-theoric decision rules for pattern classification[END_REF]. Section 5 presents experiments and the main results. Section 6 concludes the paper.

The theory of belief functions

The theory of belief functions [START_REF] Dempster | Upper and Lower probabilities induced by a multivalued mapping[END_REF][START_REF] Shafer | A mathematical theory of evidence[END_REF] is a general mathematical framework for representing beliefs and reasoning under uncertainty. In this section, we recall some concepts of this theory.

The frame of discernment Θ = {θ 1 , θ 2 , . . . , θ n } is a set of n elementary hypotheses related to a given problem. These hypotheses are exhaustive and mutually exclusive. The power set of Θ, denoted by 2 Θ is the set containing singleton hypotheses of Θ, all the disjunctions of these hypotheses as well as the empty set.

The Basic belief assignment (bba), denoted by m is a mass function defined on 2 Θ . It affects a value from [0, 1] to each subset. It is defined as:

A⊆2 Θ m(A) = 1. (1)
A focal element A is an element of 2 Θ such that m(A) > 0. A categorical bba is a bba with a unique focal element such that m(A) = 1. When this focal element is a disjunction of hypotheses then the bba models imprecision.

Based on the basic belief assignment, other belief functions (credibility function ad plausibility function) can be deduced.

-Credibility function bel(A) expresses the total belief that one allocates to A.

It is a mapping from elements of 2 Θ to [0, 1] such that:

bel(A) = B⊆A,B =∅ m(B). (2)
-Plausibility function pl(A) is defined as:

pl(A) = A∩B =∅ m(B). (3)
The plausibility function measures the maximum amount of belief that supports the proposition A by taking into account all the elements that do not contradict. The value pl(A) quantifies the maximum amount of belief that might support a subset A of Θ.

The theory of belief function is a useful tool for data fusion. In fact, for a given problem and for the same frame of discernment, it is possible to get a mass function synthesizing knowledge from separate and independent sources of information through applying a combination rule. Mainly, there exists three modes of combination:

-Conjunctive combination is used when two sources are distinct and fully reliable. In [START_REF] Smets | The Combination of Evidence in the Transferable Belief Model[END_REF], the author proposed the conjunctive combination rule which is defined as:

m 1 ∩ 2 (A) = B∩C=A m 1 (B) × m 2 (C). (4)
The Dempster's rule of combination [START_REF] Dempster | Upper and Lower probabilities induced by a multivalued mapping[END_REF] is a normalized form of the rule described previously and is defined as:

m 1⊕2 (A) =            B∩C=A m 1 (B) × m 2 (C) 1- B∩C=∅ m 1 (B) × m 2 (C) ∀A ⊆ Θ, A = ∅ 0 if A = ∅ (5)
This rule is normalized through 1-B∩C=∅ m 1 (B)×m 2 (C) and it works under the closed world assumption where all the possible hypotheses of the studied problem are supposed to be enumerated on Θ. -Disjunctive combination: In [START_REF] Smets | Belief functions: The disjunctive rule of combination and the generalized Bayesian theorem[END_REF], Smets introduced the disjunctive combination rule which combines mass functions when an unknown source is unreliable. This rule is defined as:

m 1 ∪ 2 (A) = B∪C=A m 1 (B) × m 2 (C) (6)
-Mixed combination: In [START_REF] Dubois | Representation and combination of uncertainty with belief functions and possibility measures[END_REF], the authors proposed a compromise in order to consider the benefits of the two combination modes previously described. This combination is given for every A ∈ 2 Θ by the following formula:

   m DP (A) = m 1 ∩ (A) + B∩C=∅,B∪C=A m 1 (B)m 2 (C) ∀A ∈ 2 Θ , A = ∅ m DP (∅) = 0 (7)

Decision Making in the theory of belief functions

In the transferable belief model, decision is made on the pignistic level where the belief functions are transformed into a probability function, named pignistic probability. This latter, noted as BetP is defined for each X ∈ 2 Θ , X = 0 as:

betP (X) = Y ∈2 Θ ,Y =∅ |X ∩ Y | |Y | m(Y ) 1 -m(∅) (8)
where |Y | represents the cardinality of Y .

Based on the obtained pignistic probability, we select the most suitable hypothesis with the maximum BetP. This decision results from applying tools of decision theory [START_REF] Denoeux | Analysis of evidence-theoric decision rules for pattern classification[END_REF]. In fact, if we consider an entity represented by a feature vector x. A is a finite set of possible actions A = {a 1 , . . . , a N } and Θ a finite set of hypotheses, Θ = {θ 1 , . . . , θ M }. An action a j corresponds to the action of choosing the hypothesis θ j . But, if we select a i as an action whereas the hypothesis to be considered is rather θ j then the loss occurred is λ(a i |θ j ). The expected loss associated with the choice of the action a i is defined as:

R betP (a i |x) = θj ∈Θ λ(a i |θ j )BetP (θ j ). (9)
Then, the decision consists in selecting the action which minimizes the expected loss. In addition to minimizing pignistic expected loss, other risks are presented in [START_REF] Denoeux | Analysis of evidence-theoric decision rules for pattern classification[END_REF]. Decision can be made on composite hypotheses [START_REF] Appriou | Approche générique de la gestion de l'incertain dans les processus de fusion multisenseur[END_REF][START_REF] Martin | Decision support with belief functions theory for seabed characterization[END_REF]. We present in this paper the Appriou's rule [START_REF] Appriou | Approche générique de la gestion de l'incertain dans les processus de fusion multisenseur[END_REF] which helps to choose a solution of a given problem by considering all the elements contained in 2 Θ . This approach weights the decision functions (maximum of credibility, maximum of plausibility and maximum of pignistic probability) by an utility function depending on the cardinality of the elements. A ∈ 2 Θ is chosen if:

A = argmax X∈2 Θ (m d (X)pl(X)) (10)
where m d is a mass defined by:

m d (X) = K d λ X 1 |X| r (11)
The value r is a parameter in [0, 1] helping to choose a decision which varies from a total indecision when r is equal to 0 and a decision based on a singleton when r is equal 1. λ X helps to integrate the lack of knowledge about one of the elements of 2 Θ . K d is a normalization factor and pl(X) is a plausibility function.

In the following, we present our decision rule based on a distance measure.

Decision rule based on a distance measure

In [START_REF] Essaid | Uncertainty in ontology matching: a decision rule based approach[END_REF], we proposed a decision rule based on a distance measure. It is defined as:

A = argmin(d(m comb , m A )) (12)
This rule aims at deciding on a union of singletons. It is based on the use of categorical bba which helps to adjust the degree of imprecision that has to be kept when deciding. Depending on cases, we can decide on unions of two elements or three elements, etc. The rule calculates the distance between a combined bba m comb and a categorical one m A . The minimum distance is kept and the decision corresponds to the categorical bba's element having the lowest distance with the combined bba. The rule is applied as follows:

-We consider the elements of 2 Θ . In some applications, 2 Θ can be of a large cardinality. For this reason, we may choose some elements to work on. For example, we can keep the elements of 2 Θ whose cardinality is less or equal to 2. -For each selected element, we construct its corresponding categorical bba.

-Finally, we apply Jousselme distance [START_REF] Jousselme | A New Distance Between Two Bodies of Evidence[END_REF] to calculate the distance between the combined bba and a categorical bba. The distance with the minimum value is kept. The most likely hypothesis to select is the hypothesis whose categorical bba is the nearest to the combined bba.

Jousselme distance is defined for two bbas m 1 and m 2 as follows:

d(m 1 , m 2 ) = 1 2 (m 1 -m 2 ) t D(m 1 -m 2 ) ( 13
)
where D is a matrix based on Jaccard distance as a similarity measure between focal elements. This matrix is defined as:

D(A, B) = 1 if A=B=∅ |A∩B| |A∪B| ∀A, B ∈ 2 Θ (14)
In this paper, we propose to apply the rule through two different manners:

-Distance type 1 is calculated with categorical bbas (m(A) = 1) for all elements of 2 Θ except Θ to have an imprecise result rather than a total ignorance. -Distance type 2 is calculated with simple bbas such as m(A) = α, m(Θ) = 1α.

In the following, we show that our proposed rule can be seen as a particular case of that proposed in section 3. Jousselme distance can be written as:

d(m 1 , m 2 ) = 1 2 Y ⊆Θ X⊆Θ |X ∩ Y | |X ∪ Y | m(X)m(Y ) (15)
If we consider the expected loss of choosing a i , then it can be written as:

R betP (a i |x) = Y ∈Θ λ(a i |Y )BetP (Y ). R betP (a i |x) = Y ∈Θ λ(a i |Y ) X∈Θ |X ∩ Y | |X| m(X) 1 -m(∅) . R betP (a i |x) = Y ∈Θ X∈Θ λ(a i |Y ) |X ∩ Y | |X| m(X) 1 -m(∅) . ( 16
)
The equation relative to decision is equal to that for the risk for a value of λ that has to be equal to:

λ(a i |Y ) = |X|(1 -m(∅)) |X ∪ Y | m(X) (17)
In this section, we showed that for a particular value of λ, our proposed decision rule can be considered as a particular case of that proposed in [START_REF] Denoeux | Analysis of evidence-theoric decision rules for pattern classification[END_REF]. In the following section, we give experiments and present comparisons between our decision rule based on a distance measure and that presented in [START_REF] Appriou | Approche générique de la gestion de l'incertain dans les processus de fusion multisenseur[END_REF].

Experiments

Experiments on generated mass functions

We tested the proposed rule [START_REF] Essaid | Uncertainty in ontology matching: a decision rule based approach[END_REF] on a set of mass functions generated randomly. To generate the bbas, one needs to specify the cardinality of the frame of discernment, the number of mass functions to be generated as well as the number of focal elements. The generated bbas are then combined. We use the Dempster's rule of combination, the disjunctive rule and the mixed rule. Suppose we have a frame of discernment represented as Θ = {θ 1 , θ 2 , θ 3 } and three different sources for which we generate their corresponding bbas as given in Table 1. We apply combination rules and we get the results illustrated in Table 2. Once the combination is performed, we can make decision. In Table 3, we compare between the results of three decision rules, namely the pignistic probability, the Appriou's rule with r equal to 0.5 as well as our proposed decision rule based on distance measure.

Table 3 shows the decision results obtained after applying some combination rules. We depict from this table that not all the time the rule proposed by Appriou gives a decision on a composite hypotheses. In fact, as shown in Table 3, the application of disjunctive rule as well as the mixed rule lead to a decision on a singleton which is θ 1 . This is completely different from what we obtain when we apply our proposed rule which promotes a decision on union of singletons when combining bbas. The obtained results seems to be convenient especially that the disjunctive and the mixed rules help to get results on unions of singletons.

Experiments on real databases

To test our proposed decision rule, we do some experiments on real databases (IRIS 1 and HaberMan's survival 2 ). Iris is a dataset contaning 150 instances, 4 attributes and 3 classes where each class refers to a type of iris plant. HaberMan is a dataset containing results study conducted at the University of Chicago's Billings Hospital on the survival of patients who had undergone surgery for breast cancer. This dataset contains 306 instances, 3 attributes and 2 classes (1: patient survived 5 years or longer, 2: patient died within 5 years). For the classification, our experiments are handled in two different manners.

-First, we apply the k-NN classifier [START_REF] Denoeux | A k-nearest neighbor classification rule based on Dempster-Shafer Theory[END_REF]. The results are illustrated in a confusion matrix as shown in Table 4 (left side). -Second, we modify the k-NN classifier's algorithm based on the use of Dempster rule of combination, to make it able to combine belief functions through the mixed rule. Then, Appriou's rule and our proposed decision rule are applied to make decision. Results are illustrated in Table 4. The same tests are done for HaberMan's survival dataset. The results of applying k-NN classifier, Appriou's rule and our decision rule are given respectively in Table 5. For the classification of 40 sets chosen randomly from Iris, we remark that with the k-NN classifier, all the sets having θ 1 and θ 3 as corresponding classes are well classified and only two originally belonging to class θ 2 were classified as θ 3 . Appriou's rule gives a good classification for sets originally belonging to classes θ 1 and θ 2 and thus promoting a result on singletons rather than on a union of singletons.

Considering the results obtained when applying our decision rule based on a distance type 1, we note that only 2 sets are not well classified and that 3 have θ 2 ∪ θ 3 as a class. The obtained results are good because our method is based on an imprecise decision which is underlined by the fact of obtaining θ 2 ∪ θ 3 as a class. Considering HaberMan's survival dataset, we note that the k-NN classifier, Appriou's rule as well as our decision rule give the same results where among the sets originally belonging to θ 1 , 34 are well classified and among the 18 belonging to θ 2 , only 6 are well classified. We obtain the same results as the other rules because the HaberMan's survival dataset has only two classes and our method is based on getting imprecise decisions and excluding the ignorance.

All the experiments given previously are based on the use of distance type 1. The results shown below are based on distance type 2. In fact, we consider a simple bba and each time, we assign a value α to an element of 2 Θ . The tested rule on Iris as illustrated in Table 6 (left side) gives better results with an α < 0.8. In addition to that, we obtained decisions on a union of singletons. The tests done on HaberMan's survival as given in Table 6 (right side) shows that with α > 0.5, we obtain a better rate of good classification although we did not obtain a good classification for the class θ 2 and no set belongs to Θ. We aim in the future to make experiments on other datasets because HaberMan's survival, for example, does only have 2 classes, so we do not have enough imprecise elements.

Conclusion

In this paper, we presented a rule based on a distance measure. This decision rule helps to choose the most likely hypothesis based on the calculation of the distance between a combined bba and a categorical bba. The aim of the proposed decision rule is to give results on composite hypotheses. In this paper, we demonstrated that our proposed rule can be seen as a particular case of that proposed in [START_REF] Denoeux | Analysis of evidence-theoric decision rules for pattern classification[END_REF]. We presented also the different experiments handled on generated mass functions as well as on real databases.

Table 4 .

 4 Confusion Matrices for Iriskθ1 θ2 θ1 ∪ θ2 θ3 θ1 ∪ θ3 θ2 ∪ θ3 θ1 11 θ1 θ2 θ1 ∪ θ2 θ3 θ1 ∪ θ3 θ2 ∪ θ3 θ110

Table 1 .

 1 Three sources with their bbas

		S1	S2	S3
	θ1	0.410 0.223 0.034
	θ2	0.006 0.108 0.300
	θ1 ∪ θ2 0.039 0.027 0.057
	θ3	0.026 0.093 0.128
	θ1 ∪ θ3 0.094 0.062 0.04
	θ2 ∪ θ3 0.199 0.153 0.004
	θ1 ∪ θ2 ∪ θ3 0.226 0.334 0.437

Table 2 .

 2 Combination results

		Dempster rule Disjunctive rule Mixed rule
	θ1	0.369	0.003	0.208
	θ2	0.227	0	0.128
	θ1 ∪ θ2	0.025	0.061	0.075
	θ3	0.168	0	0.094
	θ1 ∪ θ3	0.049	0.037	0.064
	θ2 ∪ θ3	0.103	0.035	0.093
	θ1 ∪ θ2 ∪ θ3	0.059	0.864	0.338

Table 3 .

 3 Decision results

		Pignistic Appriou rule Rule based on
		Probability		distance measure
	Dempster rule	θ1	θ1 ∪ θ2	θ1
	Disjunctive rule	θ1	θ1	θ1 ∪ θ2
	Mixed rule	θ1	θ1	θ1 ∪ θ2

Table 5 .

 5 Confusion Matrices for HaberMan's survival

	k-NN classifier	Appriou's rule	Our decision rule
	θ1	θ2	θ1 θ2 Θ	θ1 θ2	Θ
	θ1 34	4	θ1 34 4	0	θ1 34 4	0
	θ2 12	6	θ2 12 6	0	θ2 12 6	0

Table 6 .

 6 Rates of good classification

	α < 0.8 α >= 0.8	α <= 0.2 α ∈ [0.3, 0.5] α > 0.5
	Iris 0.95	0.675	HaberMan's survival 0.786	0.803	0.821

http://archive.ics.uci.edu/ml/datasets/Iris

http://archive.ics.uci.edu/ml/datasets/Haberman%27s+Survival


C Faber

P Boulanger

C Attaccalite

E Cannuccia

I Duchemin

T Deutsch

X Blase

Exploring approximations to the GW self-energy ionic gradients

Keywords: numbers: 71, 15, Qe,71, 38, -k

The accuracy of the many-body perturbation theory GW formalism to calculate electron-phonon coupling matrix elements has been recently demonstrated in the case of a few important systems. However, the related computational costs are high and thus represent strong limitations to its widespread application. In the present study, we explore two less demanding alternatives for the calculation of electron-phonon coupling matrix elements on the many-body perturbation theory level. Namely, we test the accuracy of the static Coulomb-hole plus screened-exchange (COHSEX) approximation and further of the constant screening approach, where variations of the screened Coulomb potential W upon small changes of the atomic positions along the vibrational eigenmodes are neglected. We find this latter approximation to be the most reliable, whereas the static COHSEX ansatz leads to substantial errors. Our conclusions are validated in a few paradigmatic cases: diamond, graphene and the C60 fullerene. These findings open the way for combining the present many-body perturbation approach with efficient linear-response theories.

I. INTRODUCTION

Electron-phonon coupling (EPC) occupies a prominent role in various fields of condensed matter physics, including phonon-mediated superconductivity, photoemission band gap renormalization, current carriers inelastic scattering or the lifetime of hot electrons. Concerning the calculation of EPC matrix elements on the ab initio level, up to now, mainly density functional theory (DFT) 1 and its perturbative linear response extensions (DFPT) 2,3 have been applied, providing important information at the microscopic level.

Recently, several studies questioned the accuracy of the DFT-based EPC matrix elements when obtained with (semi)local functionals such as LDA or PBE. 4 By way of example, the electron-phonon coupling involving specific electronic and phonon modes in graphene 5,6 and graphite, 7 the value of the electron-phonon coupling potential to states at the Femi level in the electron-doped fullerene, 8,9 superconducting bismuthates and transitionmetal chloronitrides, 10 or the renormalization of the photoemission band structure of pentacene 11 and diamond crystals, 12 were shown to be affected by a significant underestimation of the EPC matrix elements when calculated within DFT and (semi)local functionals.

As a cure to such problems, many-body perturbation theory (MBPT) techniques within the so-called GW approximation [13][14][15][16][17][18] showed a clear improvement when compared to available experimental data. 5,7,[9][10][11][12] Unfortunately, the kind of theories that are available within DFT, and in particular the powerful DFPT formalism, 2,3 are not yet available within the framework of MBPT. Existing GW calculations of EPC matrix elements have been therefore based on the frozen-phonon approach, necessitating a stepwise displacement of the atoms along the phonon modes with an explicit evaluation of the electronic structure on the GW level for each step. This demanding approach cannot be reconciled with very large unit cells and is thus only feasible for zone-center or zone-boundary phonon modes. Going beyond the frozenphonon approach, in order to access EPC matrix elements on the GW level for a large number of electron and phonon wave vectors, remains a considerable challenge to the ab initio community. Significant work is still needed to come up with a scheme that may allow to generalize the pioneering above-mentioned studies to a much larger set of systems.

In the present study, we explore the merits of the static coulomb-hole plus screened exchange (COHSEX) 14,15,[19][START_REF] March | [END_REF][START_REF] Bruneval | [END_REF][22][23][24][25][26][27] approximation to the GW selfenergy for the calculation of the electron-phonon coupling matrix elements in the isolated fullerene molecule, in diamond and in graphene. Our definitions of the relevant EPC matrix elements follow previous studies related to the superconducting transition in the fullerides, 8,9 the zero-point motion renormalization of the gap in diamond 12 and the gap opening through electron-phonon coupling in graphene, 5 respectively. For the same set of systems, we also explore the accuracy of the constantscreening approximation, namely the assumption that the screened Coulomb potential can be considered, to linear order, as a constant upon small changes of the atomic positions. Whereas COHSEX leads to non-negligible discrepancies, this latter approximation will be shown to be robust and accurate. The present results offer important perspectives in making the GW formalism amenable to the study of condensed matter phenomena involving electron-phonon coupling.

II. METHODOLOGY

We briefly introduce the many-body perturbation theory Green's function formalism, providing a solid framework for the calculation of quasiparticle energies E. In such an approach, the one-body quasiparticle eigenvalue equation reads:

-∇ 2 2 + V ion (r) φ(r) + V H (r)φ(r) + dr ′ Σ(r, r ′ ; E)φ(r ′ ) = Eφ(r),
where V ion and V H are the ionic and Hartree potential, respectively. The self-energy Σ(r, r ′ ; E) replaces the well-known exchange-correlation potential of density functional theory or the exchange operator in the Hartree-Fock formalism. In general, it is non-local, energy-dependent and non-Hermitian. Derived within Schwinger's functional derivative approach to perturbation theory, 13 the GW approximation to the self-energy leads to:

Σ GW (r, r ′ ; E) = i 2π dωe iδω W (r, r ′ ; ω)G(r, r ′ ; E + ω),
where G and W are the time-ordered one-particle Green's function and the dynamically screened Coulomb potential, respectively, and δ = 0 + a small positive infinitesimal.

The Coulomb-hole plus screened-exchange (COHSEX) representation of the GW self-energy was originally introduced using a time representation of G, W and Σ. 14,[START_REF] March | [END_REF][START_REF] Bruneval | [END_REF] We follow Hybertsen and Louie 15 by using the following spectral representations of G and W :

G(r, r ′ ; E + ω) = n φ n (r)φ * n (r ′ ) E + ω -ε n -iδsgn(ε n -µ) W (r, r ′ ; ω) = v(r, r ′ ) + ∞ 0 dω ′ 2ω ′ B(r, r ′ ; ω ′ ) ω 2 -(ω ′ -iδ) 2 ,
where G has been written in terms of one-body eigenstates φ n and eigenenergies ε n , typically starting DFT Kohn-Sham or Hartree-Fock solutions, and W in terms of its spectral function B. These expressions allow to obtain the pole structure of both G and W . From the residue theorem, one then rapidly obtains an exact decomposition of Σ:

Σ SEX (r, r ′ ; E) = - occ n W (r, r ′ ; E -ε n )φ n (r)φ * n (r ′ ) Σ COH (r, r ′ ; E) = n φ n (r)φ * n (r ′ )P ∞ 0 dω ′ B(r, r ′ ; ω ′ ) E -ε n -ω ′ ,
where P indicates the principal value. Σ SEX , which involves a summation over the occupied states only, originates from the poles of G. It is called the screened exchange interaction in analogy to the bare exchange term that can be obtained by replacing W with the energyindependent bare Coulomb potential v. Σ COH originates from the poles of W and represents the Coulomb-hole contribution, since it can be shown to be related to the interaction of an electron with its related adiabatically built correlation hole.

The static approximation to the exact COHSEX decomposition assumes that (E -ε n ) ≃ 0 for all (n), leading to simplified static screened exchange and Coulomb-hole expressions, namely:

Σ SEX static (r, r ′ ; 0) = - occ n W (r, r ′ ; 0)φ n (r)φ * n (r ′ ) Σ COH static (r, r ′ ; 0) = 1 2 n φ n (r)φ * n (r ′ ) W (r, r ′ ; 0) = 1 2 δ(r -r ′ ) W (r, r ′ ; 0).
Here, W = (W -v) is the difference between the screened and bare Coulomb potential. As a result, besides being a static approximation, this Coulomb-hole term is also local in space. Such a static COHSEX approximation (labeled COHSEX here below) was shown to yield too large gaps in the case of semiconductors. 15 By way of example, in the present case of the C 60 molecule, the COHSEX gap is found to be 5.3 eV, i.e. 0.4 eV larger than the ∼4.9 eV experimental gap. 67 Nevertheless, this has to be compared to the starting 1.6 eV DFT-LDA Kohn-Sham gap which is dramatically too small. While it cannot be claimed that the static COHSEX approach is a good approximation to absolute quasiparticle energies, we emphasize that we are interested in quasiparticle energy differences upon small (infinitesimal) atomic lattice motions. The main assumption on which we rely to calculate the electron-phonon coupling within the COHSEX approximation is that the variations of the dynamical contribution to the self-energy can be neglected. This can be a priori rationalized by emphasizing that dynamical interactions are driven by the plasmons dynamics, collective excitations less sensitive to small atomic displacements than single-particle excitation energies and wave functions. It remains, however, that besides the static approximation, the spatially local character of the static COH term is at odds with the nonlocality of the full GW self-energy.

In order to further justify the following results concerning the constant screening approach, namely the second approximation we explore in this study, it is instructive to consider the GW plus Bethe-Salpeter formalism. [28][29][30] The latter is a many-body perturbation theory approach concerned with describing the linear response of a system with respect to a time-dependent external perturbation and thus the Bethe-Salpeter equation (BSE) is the MBPT analogue to time-dependent DFT. At the heart of the GW /BSE approach lies the variation (∂GW/∂λ), where the "perturbation" λ is the one-body Green's function G. The most common approximation that has been shown to be remarkably accurate [31][32][33] is to replace the GW self-energy by its static COHSEX approximation, and to consider further that (∂W/∂λ) = 0. It is such a simplified scheme we aim to explore in the present study, differing in the fact that the perturbation λ is now induced by a vibrational distortion of the system.

III. TECHNICAL DETAILS

The many-body GW and COHSEX calculations for the isolated C 60 fullerene are performed using the Fiesta code, an implementation of the GW formalism within a Gaussian basis. [34][35][36] We start from DFT-LDA eigenstates calculated with the Siesta package 37 and a triplezeta plus polarization (TZP) basis [START_REF]3 shells of s and p orbitals[END_REF] for the description of the valence orbitals combined with standard normconserving pseudopotentials. [START_REF] Troullier | [END_REF] We use the resolution of the identity technique (RI-SVS) with an even-tempered auxiliary basis set formed of four Gaussians per angular channel (up to d character) with exponents between 0.2 and 3.2 Bohr -2 . With such a basis, our G 0 W 0 @LDA ionization potential and HOMO-LUMO gap are found to be 7.33 eV and 4.39 eV (B3LYP geometry), respectively, starting from a 1.66 eV DFT-LDA gap. This is in close agreement with the 7.45 eV and 4.40 eV plane wave G 0 W 0 @LDA values of Ref. 40 (LDA geometry). The G 0 W 0 @LDA gap remains smaller than the ∼4.9 eV experimental value, an issue that will be addressed by the partially self-consistent scheme used in the present study, as discussed below.

Following Ref. 8, we use the relaxed structure and phonon eigenmodes generated within a DFT-B3LYP approach using a 6-311G(d) basis. This was shown to provide vibrational frequencies in excellent agreement with Raman measurements. 8,41 Concerning the calculation of the EPC matrix elements, our present study follows the approach described in Refs. 8,9, where the electronphonon coupling strength associated with the (electrondoped) lowest unoccupied 3-fold degenerate t 1u molecular orbital (LUMO) has been explored within DFT and GW . The specific choice of the LUMO level is dictated by the physics of the phonon-mediated superconducting transition in the fullerides. 42,43 We adopt the standard definition of the effective electron-phonon coupling potential V ep that enters e.g. the McMillan formula for the superconducting transition temperature through the dimensionless parameter λ = N (E F )V ep , where N (E F ) is the density of states at the Fermi level. Namely, V ep reads in the so-called molecular limit: 44,45

V ep = 1 9M ν 1 ω 2 ν 3 m=1 ∂ǫ m ∂u ν 2 . (1)
Here, M is the mass of the carbon ions, ω ν is the frequency of the vibrational mode with index (ν) and u ν is the vibrational polarization vector. As already discussed before, the electronic states with energy (ǫ m ) are limited to the 3-fold degenerate LUMO level. From group theory analysis it follows that only the H g and A g vibrational modes can couple to these states. A schematic representation of the t 1u level splitting as a function of the deformation amplitude along a H g mode is provided in Fig. 1(a). The EPC matrix elements are consequently related to the slopes (∂ǫ m /∂u ν ), showing a strong dependence on the formalism adopted to calculate the energies ǫ m . These energy derivatives are calculated within the frozen phonon approach, using a symmetric five points finite-difference formula. Such an effective electron-phonon potential was calculated for C 60 in Refs. 8,9, both within DFT using diverse (semi)local and hybrid exchange-correlation functionals and GW . While DFT-LDA values were shown to significantly underestimate the coupling energy, DFT calculations performed with global hybrids containing 20%-30% of exact exchange and the GW approach were found 8,9 to compare favorably to gas phase experimental data. 46,47 As a drawback of DFT calculations with hybrids, it was evidenced that the resulting V ep potential would quickly increase with the amount of exact exchange, yielding the usual question of the proper choice of the functional parameters for a given system.

In the present study, we are concerned with exploring simplified GW schemes, i.e approximations reducing the needed computational effort. As such, quasiparticle energies are evaluated at a single-shot COHSEX and a partially self-consistent evCOHSEX level, where the quasiparticle energies are reinjected self-consistently in the construction of G and W , while keeping the starting Kohn-Sham LDA wave functions unchanged. Our reference point are partially self-consistent evGW calculations with self-consistency on the eigenvalues as presented in Ref. 9 for the calculation of the EPC matrix elements in C 60 . Such a simplified approach to full selfconsistency, justified by the dramatically too small starting point DFT-LDA gap for small isolated molecules, has been shown in the case of gas phase organic systems to produce more accurate ionization potentials and electronic affinities, 34,35,[48][49][50][51] together with improved optical excitation energies at the GW /BSE level, 36,[52][53][54] when LDA or PBE Kohn-Sham eigenstates are chosen as starting points. Consistently, when starting from DFT-LDA, EPC matrix elements in C 60 were also found to be in closer agreement with experiment within evGW as compared to G 0 W 0 . 9

For the periodic diamond and graphene systems, we use the Yambo code which implements many-body perturbation theory within a plane wave formalism. 55 The starting DFT Kohn-Sham eigenstates are generated with the Abinit code. 56 In the case of graphene, the running parameters are identical to those used in a previous GW study of the electron-phonon coupling in this material, 5 namely we use a = 2.46 Å as lattice parameter, 15 a.u. distance between the layers, 57 36 × 36 × 1 k-points for the GW calculations in the Brillouin zone (BZ) corresponding to the primitive (2-atoms) cell and the nearest equivalent k-grid in the BZ corresponding to the supercell needed to describe zone-boundary (q = K) phonons in the frozen-phonon approach. The energy cutoff on the plane wave basis was set to 60 Ry. All conduction bands within an energy range of 45 eV above the Fermi level were used to build G and W , and a cutoff of 4 Ha for the construction of the dielectric function. The Godby-Needs plasmon pole model was employed for the dynamical dielectric constant entering in W .

As schematically described in Fig. 1(b), we focus on the splitting ∆E of the degenerate occupied and unoccupied levels at the Dirac point caused by the coupling to the qν = Γ-E 2g optical mode at zone center together with the zone-boundary qν = K-A ′ 1 phonon, for which the strongest renormalization of the coupling constant have been observed upon replacing the DFT-LDA approach by the G 0 W 0 formalism. 5,6 The EPC matrix elements of interest are thus here related to the derivatives: (∂∆E/∂u qν ), namely (see Ref. 5):

D 2 qν = 1 A qν ∂∆E ∂u νq 2 , (2)
expressed in (eV / Å) 2 , where A qν equals 16 and 8 for qν = Γ-E 2g and qν = K-A ′ 1 , respectively. We performed "single-shot"G 0 W 0 and COHSEX calculations on top of DFT-LDA results, namely the standard perturba-tive (non self-consistent) GW approach to obtain quasiparticle energies.

For diamond, following Antonius and coworkers, 12 we have used a relaxed DFT-LDA geometry (lattice parameter a = 3.591 Å) with a plane wave cutoff energy of 80 Ry and norm-conserving pseudopotentials. All calculations were performed with a 8x4x4 uniform Monkhorst-Pack grid to describe the supercell necessary for the frozenphonon simulation of the (q = X) phonon. We have used 160 bands and a cutoff of 8 Ha for the construction of the dielectric function and the GW self-energy. Following Refs. 12 and 58 where the diamond band-gap renormalization by the zero-point motion was studied, we define the relevant EPC matrix elements as the second derivative of the electronic eigenvalues with respect to the phonon displacements at equilibrium:

∂ǫ mk ∂n νq = h 2M ω νq ∂ 2 ∂u 2 νq ǫ mk , (3)
where n νq is the occupation of the phonon state (ν) of wave vector q. The EPC matrix elements are now related to the temperature dependent renormalization of the electronic bands. We evaluate the second-derivative of the eigenvalues using a symmetric five-point finitedifference formula with displacements of 0.01, 0.02 and 0.03 Å. In Ref. 59, it was already shown that standard DFT and DFPT methods with (semi)local functionals strongly underestimate (by ca. 30%) the zero-point renormalization of the direct gap of diamond compared to experiment, while Ref. 12 demonstrated that G 0 W 0 calculations restored the good agreement between theory and experiment. We focus here below on the renormalization of the lowest conduction band edge at (mk = X 1c ) and (mk = Γ 15c ) by the coupling to the (νq = X 4 ) phonon modes at the zone-boundary (see Fig. 1c). It is for such a phonon wave vector that the EPC matrix elements were found to be large and that also the GW correction was shown to be the most significant. 12 Following Ref. ) phonon mode and are expressed in (eV / Å) 2 . Number in parenthesis are obtained with a plasmon pole energy of 7 eV, close to the graphene π-plasmon energy, instead of the 27 eV originally used in Ref. 5 (see text). For diamond, these are the results of equation ( 3) for the X1c and Γ15c conduction states and for the X4 zone-boundary phonon mode only (values in eV). This corresponds to twice the zero-point renormalization arising from this given phonon mode (see Ref. 12).

averaged over the degenerate electronic manifold considered (namely the 3 degenerate levels at zone-center and the 2 degenerate levels at the X 1c zone boundary).

IV. RESULTS

A. Reference GW calculations

We first reproduce for validation and for reference the previously published GW calculations related to EPC matrix elements in graphene, diamond and C 60 . Our results are compiled in Table I for diamond and graphene and in Table II for C 60 , respectively. Our LDA and G 0 W 0 values for graphene are close to those found in Ref. 5, namely e.g. 201 eV/ Å (G 0 W 0 value, present study) for the largest matrix element with the (K-A ′ 1 ) phonon, to be compared to 193 eV/ Å in the previous study. 57 The difference can be explained by an increased five points, instead of only two, finite-difference formula, and by the increase of the dielectric matrix size energy cutoff, from 2 to 4 Ha cutoff. We also verified that changing the Godby-Needs plasmon model input finite frequency, from the (default) 27 eV value in Ref. 5 to a 7 eV value closer to the π-plasmon resonance in graphene, does not change significantly the calculated EPC strength (see numbers in parenthesis in the Table I). Such differences are negligible with respect to the more than 100% increase as compared to the LDA value.

For diamond, the EPC contribution associated with the mk = Γ 15c conduction states are in good agreement with the LDA and G 0 W 0 results of Antonius and coworkers, 12 i.e. we find a 11 meV difference with respect to their EPC in the case of LDA, while we find a difference of 8 meV for the G 0 W 0 results. 68 These small variations can be ascribed to the different pseudopotentials (Ref. 59 showed that pseudopotentials can lead to errors up to 50 meV), and to the different convergence parameters.

Since Ref. 12 focused on the renormalization of the direct band gap at the zone-center, comparison with the present results for the phonon coupling to the X 1c elec-tronic state cannot be made. Our findings are, however, very consistent with what is observed for the Γ 15c state, namely a dramatic decrease of the coupling strength with the X 4 phonon mode which, again, is the dominant coupling mode. The coupling energy is indeed found to be reduced by 52% and 49% from LDA to GW for the Γ 15c and X 1c states, respectively. Since the X 1c state is closer to the true conduction band minimum, this is a strong indication that zero-point motion renormalization of the indirect gap will certainly also be strongly affected by the GW correction. The full study of such an effect is beyond the scope of the present paper.

For C 60 , the total evGW coupling potential is within 8% of that found by Ref. 9 as a result of the larger TZP basis used here. [START_REF]3 shells of s and p orbitals[END_REF] Our evGW V ep potential is found to increase by 44% with respect to the corresponding LDA value, to be compared with the 48% increase obtained in Ref. 9 with a smaller basis. This indicates the good convergence of the GW correction to the LDA value. We note that while DFT and GW calculations are known to converge very differently with respect to basis size, we compare in the present study approaches which are much more similar, namely "standard and approximated" GW calculations, indicating certainly an even better convergence of the differences we are interested in.

B. The COHSEX approximation

We now explore approximations to the full GW calculations performed in Refs. 5,9,12 and reproduced here above. We start by the case of diamond. As compared to the G 0 W 0 calculations, the static COHSEX approximation is shown to induce errors larger than 100% for the relevant electron-phonon coupling energy to the Γ 15c and X 1c states. In both cases, the static COHSEX approach strongly overestimates the electron-phonon coupling. A representation of the errors as compared to the G 0 W 0 calculations is provided in Fig. 2.

The case of graphene is of specific interest since, in great contrast to diamond and C 60 , the fundamental gap reduces to zero in the equilibrium geometry and opens

Γ 15c state X 1c state LDA G 0 W 0 G 0 W 0 (W) COHSEX COHSEX(W) LDA G 0 W 0 G 0 W 0 (W) COHSEX COHSEX(W)
Diamond FIG. 2: Relative error of the different approximations with respect to the G0W0 results for the "thermal renormalization" EPC energies (Eqn. 3) for the lowest conduction band of diamond.

upon phonon distortions. For the optical phonon mode at the zone-center, the COHSEX approximation leads again to an increase as compared to the G 0 W 0 value, even though smaller than in the diamond case with a 20% error. Such an error becomes larger for the K-A ′ 1 mode with a 29% increase of the EPC-related gap-opening rate. In the present graphene case, the error induced by the COHSEX approach remains smaller than the 34% and 56% reduction observed upon using DFT-LDA instead of G 0 W 0 . However, the close to 30% error observed for the COHSEX calculation of D 2 K is clearly significant (see Fig. 3), questioning again the applicability of the static COHSEX formalism for EPC matrix elements calculations.

We now consider the C 60 case, where we compare evGW and evCOHSEX calculations. A global increase of about 7% from evGW to its static approximation is observed. Besides the H g (4) and A g (1) modes showing very small couplings, the static COHSEX approximation

Γ (E 2g ) mode K (A 1 ') mode LDA G 0 W 0 G 0 W 0 (W) COHSEX COHSEX(W) LDA G 0 W 0 G 0 W 0 (W) COHSEX COHSEX(W)
Graphene FIG. 3: Relative error of the different approximations with respect to the G0W0 results for the electron-phonon coupling in graphene.

is observed again to systematically increase the coupling strength as compared to evGW . 62 This result is consistent with what was observed for diamond and graphene, even though clearly the error appears to vary from one system to another.

C. The constant screening approximation

Following the above-mentioned analogy with the Bethe-Salpeter formalism, we finally test another approximation, that is the constant screening approach. Namely, we explicitly calculate the screened Coulomb potential for the undistorted structure and upon changing the positions of the ions along the (R qν ) vibrational eigenmodes, we keep W frozen to the undistorted value. Formally, this amounts to assuming that: ∂GW/∂u qν ≈ (∂G/∂u qν )W . Such an approach will be labeled GW (W ) or COHSEX(W ) in the following. In the case of the fullerene, where calculations are performed within a Gaussian basis, special care must be taken in the implementation of such a constant-screening approximation. In the auxiliary basis, or resolution of the identity approach, the bare and screened Coulomb potentials are expressed in terms of an atom-centered auxiliary basis with the following relations:

[W ] β,β ′ = drdr ′ β(r)W (r, r ′ )β ′ (r ′ ), (4)
W (r, r ′ ) = β,β ′ β(r) S -1 [W ]S -1 β,β ′ β ′ (r ′ ), (5)
where S is the overlap matrix in the auxiliary basis. Using now the notation: W and β for the screened Coulomb potential and the auxiliary basis for the slightly distorted system, the assumption: W (r, r ′ ) ≃ W (r, r ′ ) leads straightforwardly to the condition:

[W ] β,β ′ ≃ S ββ [W ] β,β ′ S β ′ β ′ , (6)
where S ββ =< β|β > is an overlap matrix between the auxiliary bases for the perturbed and unperturbed systems, respectively. For C 60 again, we test the constant screening approximation at the COHSEX level only, namely comparing evCOHSEX and evCOHSEX(W) calculations of the V ep energies. While the plasmon-pole approach used for diamond and graphene is based on a fixed finite frequency pole value independent of the geometry, the real axis poles contribution to the correlation energy in the contour deformation approach (see Refs. 19,34) implemented in the Fiesta package changes from one structure to another, leading to difficulties when trying to implement the constant-W approach within the frozen-phonon scheme.

The results of the constant-screening approximation are compiled in Tables I and II. For C 60 , an excellent agreement is obtained within evCOHSEX(W ) compared to the corresponding evCOHSEX calculations.

Comparing the total coupling, evCOHSEX(W ) and ev-COHSEX agree within 1.5 %.

We now explore the constant-screening approximation for periodic diamond. Results are again compiled in Table I. For diamond, we observe a 3.4% and 7.5% difference between G 0 W 0 (W ) and G 0 W 0 for the X 1c and Γ 15c states, respectively. At the COHSEX level, the constant screening approximation yields errors of the order of 7.0% and 10.2%, respectively. In summary, while the COHSEX approximation was failing for diamond, the constant-W approximation leads again to very reasonnable results, with errors well within 10% at the GW level. As illustrated in Fig. 2, this error is much smaller than the error induced by the LDA approximation.

We finally address the case of graphene. Concerning the coupling with the zone-center optical mode, the constant-W approach leads to a 5.7% and 2.1% error respectively when comparing G 0 W 0 (W ) to G 0 W 0 and COHSEX(W) to COHSEX. In the case of coupling with the zone-corner K-A ′ 1 phonon mode, the discrepancy is of 3% and 1.5% applying the constant-screening approximation to GW and COHSEX approaches, respectively. Notice that the results for the constant-screening of the K-A ′ 1 phonon are affected by a numerical error of about 4%, therefore the discrepancy could be slightly larger. Anyway this difference is much smaller than the error induced by the standard DFT-LDA calculations, but also much smaller than the error induced by the static COH-SEX formalism as compared to the reference GW calculation.

V. DISCUSSION

Clearly, among the two approximations explored here above, the constant-screening stands as a much better approach than the static COHSEX for the calculation of electron-phonon coupling matrix elements within the GW formalism. Overall, the largest error induced by the static COHSEX approximation is larger than 100% in the diamond case, while it is reduced to 7% in the case of the constant-screening approximation to full GW calculations.

The static COHSEX approximation is known to always overestimate band gaps. As proposed in Ref. 8, a too large band gap should lead to an underscreening of the electron-phonon interaction and consequently to enhanced EPC matrix elements. This is consistent with the fact that, within Hartree-Fock or with increasing exact exchange in hybrid functionals, the EPC matrix elements are found to steeply increase. The Hartree-Fock D 2 Γ and D 2 K coupling constants in graphene were found to be about 5 and 30 times larger, respectively, as compared to their GW analog, 5 while in C 60 increasing the percentage of exact exchange from 20% to 30% in hybrid functionals was found to enhance the V ep energy by about 15%. 8 Such an interpretation matches the observation that in graphene and C 60 , the COHSEX EPC coupling constants are larger than their corresponding GW reference. In the case of diamond, where we are concerned with the second-order derivatives, the effect of assuming a static approximation may be more difficult to interpret.

An important observation further is that the Coulombhole term is a local potential in the static COHSEX approximation, washing out spatial local-field effects, a observation that may potentially allow to understand the failure of the COHSEX approximation to reproduce accurately the evolution of the band structure with phonon deformation. In any case, the large variations of the induced errors from diamond to graphene and C 60 is still to be understood, and we will just stand here on the observation that the static COHSEX approximation cannot be trusted to improve on DFT-LDA calculations.

Concerning the neglect of the gradient of W with respect to the ionic positions, such an approximation was already tested by Ismail-Beigi and Louie in the context of excited state ionic forces within the Bethe-Salpeter formalism, 63 showing good agreement with explicit finite-difference "exact" BSE calculations for small CO and NH 3 molecules. It is commonly assumed that the GW /BSE approach is much more resistant to approximations on W as compared to the GW approach for (charged) excitations. This is due to cancellations of errors between the electron-electron and electron-hole interactions. Namely, any error introduced in W is expected to affect excitonic interactions and quasiparticle gaps in opposite ways. Clearly, the present GW study of the variations of a given quasiparticle energy with respect to ionic positions cannot benefit from such cancellation of errors. Still, the constant-screening approximation turns to be a reliable approach to save on computational cost.

An important consequence of the present findings is that once the screened Coulomb potential W (r, r ′ ; ω) is built for the equilibrium geometry, the calculation of the variations of the quasiparticle energies with respect to the perturbation (λ) only requires the evaluation of the variations of the Green's function G with respect to the perturbation. This can be performed within standard DFPT techniques, at least in the case of non-selfconsistent G 0 W 0 calculations where the Green's function assumes an explicit form in function of the input DFT eigenstates. This may invite, for systems such as C 60 , to use G 0 W 0 calculations starting from DFT eigenstates obtained with hybrid functionals, which have been shown to be a better starting point for finite-size systems. 54,[64][65][66]

VI. CONCLUSIONS

We have explored two approximations for calculating self-energy gradients of interest for electron-phonon coupling, namely the two simplifications commonly used in the GW /Bethe-Salpeter calculations, that is the static COHSEX and the constant screening approximations. We explored these approaches in the case of diamond, graphene and C 60 . Our findings suggest that the COH-SEX approximation cannot be trusted to improve on the DFT-LDA values as clearly illustrated in particular in the case of diamond. On the contrary, the constant screening hypothesis, namely assuming that W remains to first order constant with respect to small ionic displacements, seems to be reliable, with a discrepancy no larger than 10%, even in the difficult case of graphene where the phonon perturbation dramatically affects the Dirac cone and the (semi)metallic nature of the graphene sheet. Even though a deeper understanding of the origin of the specific difficulties uncountered by the static COHSEX approximation would allow to better rationalize the validity and limits of the tested approximations, the present results offer promising perspectives to carry on such many-body evaluations of the electron-phonon coupling gradients with much reduced computer cost on realistic systems, including the study of periodic systems with arbitrary wave vector perturbation. Further studies are however required on a larger set of systems and physical observables in order to better assess the interest, with respect to common DFT calculations, of using the GW approach, and its various approximations, for calculating self-energy gradients with respect to ionic positions.

FIG. 1 :

 1 FIG. 1: Schematic representation of changes in the electronic structure under phonon distortion: a) in C60, the 3-fold LUMO levels (blue) split under a distortion (red) along a Hg phonon mode. b) In graphene, an energy gap ∆E at the Dirac point is opened (red) when distorting the equilibrium structure (blue) along the Γ -A ′ 1 or K -E2g phonon modes. c) In diamond, the lowest lying conduction bands at Γ/X (blue) are renormalised by their coupling to a zone-boundary phonon mode. For clarity, only the states showing the largest shift are plotted (red dot, see text).

FullereneFIG. 4 :

 4 FIG.4: Relative error of the different approximations with respect to the evGW results for the Hg-related V ep of equation(1). For convenience, a comparison to the experimental data is provided.

TABLE II :

 II Mode-resolved V ep ν coupling potential associated with the 3-fold degenerate C60 LUMO and the corresponding V ep totals (in meV). For group symmetry reasons, only the Hg and Ag vibrational modes couple to this electronic state. a Ref. 9, b Ref. 46, Table V, c Ref. 47.

Acknowledgments. P.B. warmly acknowledges G. Antonius for providing the diamond phonon eigenmodes used in this study and for guiding our comparison with Ref. 12. C. A. acknowledges Daniele Varsano for the correction of some bugs in the Yambo code. C.F. is indebted to the French CNRS and CEA for Ph.D funding. P.B. acknowledges a postdoctoral fellowship from the French National Research Agency under Contract No. ANR-2012-BS04 PANELS. Computing time has been provided by the "Curie" national GENCI-IDRIS supercomputing center under contract No. i2012096655 and a PRACE European project under contract No. 2012071258. X.B. acknowledges illuminating discussions with Steven G. Louie and G. Antonius.

ciated with the (single-shot) COHSEX approach, starting from DFT-LDA eigenstates. Clearly, the EPC coupling energies approach from below the evGW value with an error not larger than that induced by the evCOHSEX technique. This stems from a cancellation of errors between the COH-SEX approximation and the lack of self-consistency which leads to overscreening when building W directly from the DFT-LDA spectrum. 63 The band gap is calculated as the difference of the experimental ionization energy and electronic affinity. We adopt the most recent experimental data. 68 Direct comparison with Figure 2 of Ref. 12 is not possible, since the figure presents the total electron-phonon coupling at wavevector q = X, i.e. all phonon modes are considered and not just X4 as in this study. The comparison was made directly with the data used in Ref. 12 via personal correspondence.