Crouzeix Michel
email: michel.crouzeix@univ-rennes.fr

Greenbaum Anne

Li Kenan

Michel Crouzeix

Anne Greenbaum
email: greenbau@uw.edu

Kenan Li

Numerical bounds on the

Keywords: 2000 Mathematical subject classifications : 15A60, 15A45, 47A25, 47A30 numerical range, spectral set

come

Introduction

We consider n by n matrices A n with 1's in the strict upper triangle and 0's elsewhere. For n = 3, 4, 5, 6, we have numerically determined upper and lower bounds on the value ψ(A n ) := sup{ p(A n ) : p a polynomial with |p| ≤ 1 in W (A n )}, where W (A n ) denotes the numerical range of A n and • is the spectral norm in C n,n . (We refer to this quantity as the Crouzeix ratio although sometimes this term denotes the reciprocal, max z∈W (An) |p(z)|/ p(A n ) , for a given polynomial p [START_REF] Overton | Local minimizers of the Crouzeix ratio: a nonsmooth optimization case study[END_REF].) Crouzeix's conjecture is that for all square matrices A, ψ(A) ≤ 2. A way to determine ψ(A) is to introduce a Riemann mapping g from the interior of W (A) onto the open unit disk D and to consider the matrix M := g(A). In this case, we can write

ψ(A) = ψ D (M ) := max{ f (M ) : f holomorphic in D with |f | ≤ 1 in D}.
We know that the maximum is realized by a Blaschke product of order n -1 and each choice of such a Blaschke product b provides a lower bound: ψ(A) = ψ D (M ) ≥ b(M ) . From the von Neumann inequality, we easily deduce

ψ D (M ) ≤ ψ cb,D (M ) := min{cond(H) : H ∈ C n,n , H -1 M H ≤ 1}.
Thus, the exhibition of a matrix H satisfying H -1 M H ≤ 1 leads to an upper bound ψ(A) ≤ cond(H) := H • H -1 . This approach has been used to prove ψ(A) ≤ 2 for some classes of matrices [START_REF] Crouzeix | Bounds for analytic functions of matrices[END_REF][START_REF] Choi | A proof of Crouzeix's conjecture for a class of matrices[END_REF][START_REF] Choi | Roots of matrices in the study of GMRES convergence and Crouzeix's conjecture[END_REF][START_REF] Glader | Crouzeix's conjecture holds for tridiagonal 3 × 3 matrices with elliptic numerical range centered at an eigenvalue[END_REF], but this assumes that one knows the matrix M with sufficient accuracy. This is the case when W (A) is an ellipse since there is an analytic formula for g and also for [START_REF] Choi | A proof of Crouzeix's conjecture for a class of matrices[END_REF][START_REF] Choi | Roots of matrices in the study of GMRES convergence and Crouzeix's conjecture[END_REF] where g(A) = cA for a certain constant c. In general, however, there is no simple expression for the boundary of W (A) and for the Riemann mapping g. There are numerical methods for computing g and thus M with high precision, but to guarantee the accuracy would require a complete analysis of all discretization and rounding errors.

In section 2, we consider the matrix

A 3 =   0 1 1 0 0 1 0 0 0   .
Our numerical and analytical work suggests that 1.9956978 < ψ(A 3 ) < 1.9956979, and we are confident in this range, although it does not provide a proof that ψ(A 3 ) ≤ 2, since it relies on numerical computation of g(A).

Another approach is to identify a rational function f 1 such that the image Ω := f 1 (D) of the unit disk is a subset of (and close to) W (A 3 ). Let g 1 be the inverse of f 1 . Then with M 1 = g 1 (A 3 ), we can write

ψ Ω (A 3 ) := sup{ h(A 3 ) : |h| ≤ 1 in Ω} = ψ D (M 1 ) := sup{ f (M 1 ) : |f | ≤ 1 in D}.
Since Ω ⊂ W (A 3 ), we clearly have ψ(A 3 ) ≤ ψ Ω (A 3 ) ≤ ψ cb,D (M 1 ). We are now able to compute M 1 analytically and exhibit a matrix H 1 such that H -1

1 M 1 H 1 = 1 and cond(H 1 ) ≈ 1.9999514. However, we must verify numerically that Ω ⊂ W (A 3 ).

In section 3, we consider more generally the matrices A n for n > 3. They belong to the class of KMS matrices [START_REF] Kac | On the eigenvalues of certain Hermitian forms[END_REF] and they are the matrices in this class for which the boundary of the numerical range contains a line segment [START_REF] Gau | Numerical ranges of KMS matrices[END_REF]. We derive a simple description of their numerical ranges and determine numerically the following bounds:

1.993800 ≤ ψ(A 4 ) ≤ 1.993801, 1.992921 ≤ ψ(A 5 ) ≤ 1.992922, 1.992444 ≤ ψ(A 6 ) ≤ 1.992445.

In section 4, we explain the numerical method used to compute the conformal mapping and we provide the Matlab code used for the computation of M = g(A).

2 Numerical estimates for the matrix A 3

Here we consider the matrix A 3 =   0 1 1 0 0 1 0 0 0   . We will see in the next section that the boundary of its numerical range is the union of a part of an algebraic curve 2 e iθ +e 2iθ 3

: -2π 3 ≤ θ ≤ 2π 3 and of the vertical straight line [-1 2 -i √ 3 6 , -1 2 + i √ 3
6 ]. The algebraic curve is a cardioid ; its Cartesian equation is 27(x 2 +y 2 ) 2 -18(x 2 +y 2 ) -8x -1 = 0.

We denote by g the Riemann mapping from W (A) onto the unit disk D that satisfies g(0) = 0 and g (0) > 0. Then

M = g(A) =   0 a b 0 0 a 0 0 0   , with a = g (0), b = g (0) + 1 2 g (0).
From the numerical computations it appears that a 1.360374515, b 0.710915425 with an accuracy that we empirically estimate better than 10 -8 . Using the Blaschke product f (z) = , we obtain f (M ) = 1.9956978, which (numerically) shows that ψ(A) ≥ 1.9956978. We now choose the matrix

H =   a b/2a -b 2 /8a 3 0 1 -b/2a 2 0 0 1/a   , then H -1 M H =   0 1 0 0 0 1 0 0 0   . Therefore ψ cb,D (M ) ≤ H H -1 = cond(H) 1.995697855.
With the numerical values obtained previously, we believe that we have the two sided estimates 1.9956978 < ψ M (a, b) < 1.9956979. Therefore, it appears from the numerical simulation that W (A) is a (complete) 1.9956979-spectral set for A, that the complete bound [START_REF] Paulsen | Completely bounded maps and operator algebras[END_REF] is the same as the ordinary bound, and that a function which realizes ψ(A) is a Blaschke product of order 2 with 2 real roots. But we have only an empiric estimate of the accuracy that we justify as follows. If we let a(n) be the numerical value of a computed with our program (described further) using n points on the boundary, we have verified that 33 ≤ a(1447)-a(n) n -4 ≤ 65 for values of n between 23 and 1205. This suggests that our method is of order n -4 and suggests that |a(1205) -a| ≤ 2 • 10 -11 . Similarly, it appears that our computation of b is of order n

-4 , 130 ≤ b(n)-b(1447) n -4
≤ 350 and that |b(1205) -b| ≤ 10 -10 .

We turn now to the second attempt which is to consider the image Ω = f 1 (D) of the unit disk by the rational function

f 1 (z) = (c 1 z + c 2 z 2 + • • • + c 7 z 7 )/(1 + d 1 z + • • • + d 7 z 7 )
, with the values c = (0.734, 0.49736, 0.07268, -0.00521, 0.00013, 0.00061, -0.00251), d = (0.32564, -0.03291, 0.01, -0.004, 0.00084, -0.00242, 0.00028).

We let g 1 be the inverse function of f 1 and we set M 1 = g 1 (A). We will see that if Ω is included in W (A), then We are now able to compute a 1 = g 1 (0), b 1 = g 1 (0) + 1 2 g 1 (0), and thus M 1 with an accuracy better than 10 -14 . We choose

ψ(A) ≤ ψ Ω (A) ≤ ψ cb,D (M 1 ) ≤ cond(H 1 ), if H -1 1 M 1 H 1 ≤ 1.
H 1 =   a 1 b 1 /2a 1 -b 2 1 /8a 3 1 0 1 -b 1 /2a 2 1 0 0 1/a 1   , then H -1 1 M 1 H 1 =   0 1 0 0 0 1 0 0 0   .
This gives an estimate ψ(A) ≤ ψ cb,D (M 1 ) ≤ 1.9996222 with an accuracy better than 10 -12 which ensures that W (A) is a 2-spectral set for A.

It remains to show that W (A) contains Ω. For that, we first remark that the set {z : p(z) < 0} with p(z) := 27|z| 4 -18|z| 2 -8 Re z -1 < 0 is the interior of the cardioid and that the rectangle {z :

-1 2 ≤ Re z ≤ 0 and | Im z| ≤ √ 3/6} is contained in W (A).
Taking into account the symmetry with respect to the real axis, in order to show that Ω is contained in W (A), it suffices to show that the set {z = f 1 (e iθ ) : 0 ≤ θ ≤ 3π 4 } is interior to the cardioid and that the set {z = f 1 (e iθ ) : Figure 5:

Curve d dθ p(f 1 (e iθ )), 0 ≤ θ ≤ 3π 4 .
θ ≤ π, an estimate of max d dθ p(f 1 (e iθ )) ≤ 0.015 and thus min Re(f 1 (e iθ )) ≥ -0.4998968 -0.015π/2000 > -0.499921. Note also that this part of the curve clearly satisfies Re(f 1 (e iθ )) ≤ 0 and Im f 1 (e iθ )) ≤ √ 3/6. This shows that the set {z = f 1 (e iθ ) : 3π 4 ≤ θ ≤ π} is interior to W (A). 3 Estimates for the class of matrices A k

Recall [START_REF] Kippenhahn | On the numerical range of a matrix, Translated from the German by Paul F. Zachlin and Michiel E. Hochstenbach[END_REF] that the boundary of the numerical range W (A) is the convex hull of the algebraic curve with tangential equation T (u, v, w) := det(uB+vC+wI) = 0 where we have written A = B+i C, with B and C self-adjoint. For the matrix A k , the corresponding tangential equation T k (u, v, w) = 0, can be obtained from the recursion

T 1 (u, v, w) = w, T k+1 (u, v, w) = w T k (u, v, w) + (w -u+iv 2 ) k -(w -u-iv 2 ) k i v u 2 + v 2 4 .
For instance,

T 3 (u, v, w) = w 3 - 3 4 w(u 2 +v 2 ) + 1 4 (u 2 +v 2 )u, T 4 (u, v, w) = w 4 - 3 2 w 2 (u 2 +v 2 ) + w(u 2 +v 2 )u - 1 16 (u 2 +v 2 )(3u 2 -v 2 ).
We can see, by recursion, that

T k (cos ϕ, sin ϕ, w) = (-1) k sin ϕ Im e -iϕ ( 1 2 e iϕ -w) k .
We now remark that, if ϕ = kθ 2 and w = -1 2 sin((k-1)θ/2) sin (θ/2) , then e -iϕ/k ( 1 2 e iϕ -w) = 1 2 sin(kθ/2) sin θ/2 ∈ R, whence T k (cos ϕ, sin ϕ, w) = 0. Now, we consider the algebraic curve {f k (e iθ ) : |θ| ≤ π} with

f k (z) = ( k-1 j=1 jz k-j )/k = z k+1 -z-k(z 2 -z) k(z-1) 2
. We remark that

k e iθ f k (e iθ ) = k-1 j=1 j(k-j)e i(k-j)θ = e ikθ/2 2 k-1 j=1 j(k-j) cos (k-2j)θ 2 ;
hence, the vector e ikθ/2 is the unit normal at the point f k (e iθ ) and the equation of the tangent at this point is Re (e -ikθ/2 (x+iy-f (e iθ ) = 0, i.e. ux+vy+w = 0, with

u = cos kθ 2 , v = sin kθ 2 , w = -Re e -ikθ/2 f (e iθ ) = -1 2 k-1 j=1 cos( k-2j 2 θ) = -1 2 sin((k-1)θ/2) sin(θ/2)
. Thus, this shows that the algebraic curve which satisfies the tangential equation T k (u, v, w) = 0 is the set:

1 {z : z = 1 k k-1
j=1 j e i(k-j)θ , -π ≤ θ ≤ π}. The points with horizontal tangent are given by θ j = (2j-1)π k , j = 1, . . . , n and are cuspid for j = 2, . . . , k-1; the k-1 points of the algebraic curve on the flat part are the points -1 2 + i 2 tan(jπ/k) , j = 1, . . . , k-1. Note that the matrix A k +A * k +I, whose entries are all 1's, has the simple eigenvalue k (eigenvector (1, 1, • • • , 1) T ), and the eigenvalue 0 of multiplicity k-1 (with orthonormal eigenvectors 1 √ k (e iθ j , e 2iθ j , . . . , e ikθ j ) T , θ j = 2jπ k , j = 1 . . . , k-1). This implies that the point k-1 2 = f k (1) is the extremal right point of W (A k ), and that the boundary has a flat part on the line Re z = -1 2 . The boundary of the numerical range is the union of a part of the algebraic curve and of a straight part:

∂W (A k ) = {z : z = 1 k k-1 j=1 j e i(k-j)θ , |θ| ≤ 2π k } ∪ {z : z = -1 2 (1+iy), -cot( π k ) ≤ y ≤ cot( π k )}.
We turn now to some estimates for the constant corresponding to the matrices A k , k ≤ 6. For We have obtained a lower bound ψ(A 4 ) ≥ 1.9938003 with the Blaschke product corresponding to the 3 coefficients -0.4560323±0.3891911i, 0.2474013. We have also an upper bound: with the values x = b/a 1.5 , y = az-bx/a+c/a 1.5 , z = -0.0735033, t = -0.0231366 and the matrix

H =     a 1.5 xa y t 0 a 0.5 0 z 0 0 a -0.5 -xa 0 0 0 a -1.5     , it holds H -1 g(A 4 )H =     0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0     .
Then, cond(H) 1.9938002 and H -1 g(A 4 )H = 1. We can consider that the estimate 1.993800 ≤ ψ(A 4 ) ≤ 1.993801 is correct.

For A 5 , we have computed

g(A 5 ) =       0 a b c d 0 0 a b c 0 0 0 a b 0 0 0 0 a 0 0 0 0 0       , with a = 1.1170233, b = 0.2325756, c = 0.2187502, d = 0.1895824.
We have obtained a lower bound ψ(A 5 ) ≥ 1.9929216 with the Blaschke product corresponding to the 4 coefficients -0.2583004±0.60451151i, -0.6247827, 0.3295365. We have also an upper bound: with the values u = -0.0194597, w = -0.0384976, g = -0.1091772, h = -0.2503045,

f = ah + ba -2 , v = ag+bh + ca -2 , z = af +ba -1 , t = aw+bg+ch+da -2 , y = av+bf +ca -1 , x = az+b and the matrix H =       a 2 x y t u 0 a z v w 0 0 1 f g 0 0 0 a -1 h 0 0 0 0 a -2       , it holds H -1 g(A 5 )H =       0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0       .
Then, cond(H) 1.9929216 and H -1 g(A 5 )H = 1.

We can consider that the estimate 1.992921 ≤ ψ(A 5 ) ≤ 1.992922 is correct.

For A 6 , we have computed We have obtained a lower bound ψ(A 6 ) ≥ 1.9924447 with the Blaschke product corresponding to the 5 coefficients -0.5859775±0.3199164i, -0.0604565±0.70030221i, 0.3972632, 0.3295365. We have also a upper bound: with the values y = (-0.0163999, -0.0248879, -0.0578414, -0.1294105, -0.243031), x 9 = a y 5 + b a -2.5 , x 8 = a y 4 + b y 5 + c a -2.5 , x 7 = a x 9 + b a -1.5 , x 6 = a y 3 + b y 4 + c y 5 + d a -2.5 , x 5 = a x 8 + b x 9 + c a -1.5 , x 4 = a x 7 + b a -0.5 , x 3 = a y 2 + b y 3 + c y 4 + d y 5 + e a -2.5 , x 2 = a x 6 + b x 8 + c x 9 + d a -1.5 , x 1 = a x 5 + b x 7 + c a -0.5 ; x 0 = a x 4 + b a 0.5 , and the matrix

g(A 6 ) =         0
H =         a 2.5 x 0 x 1 x 2 x 3 y 1 0 a 1.5 x 4 x 5 x 6 y 2 0 0 a .5 x 7 x 8 y 3 0 0 0 a -.5 x 9 y 4 0 0 0 0 a -1.5 y 5 0 0 0 0 0 a -2.5         , it holds H -1 g(A 6 )H =        
0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0

        .
Then, cond(H) 1.9924445 and H -1 g(A 6 )H = 1. We can consider that the estimate 1.992444 ≤ ψ(A 6 ) ≤ 1.992445 is correct. Remark. In each of these cases, the matrix M = g(A k ) satisfies the relation ψ D (M ) = ψ cb,D (M ). This property holds for all d × d matrices M if d ≤ 2, but may fail [START_REF] Crouzeix | Polynomial bounds for small matrices[END_REF] if d ≥ 3. Also here, for the matrix H which realizes ψ cb,D (M ), H -1 M H was a Jordan block, which is not generally the case.

Table 1

A n ) = ψ cb,D (M n ).
4 About the computation of the conformal mapping g

We may write g(z) = z exp(u+iv), with u(z) and v(z) harmonic real-valued functions. Note that u(z) = -log |z| on ∂W (A), which determines u in W (A) in a unique way.

Let us consider a representation ∂W (A) = {σ(θ) ; θ ∈ [0, 2π]} of the boundary and choose λ > 0. If q is a 2π-periodic real-valued function such that,

2π 0 q(θ) log σ(θ)-σ(ϕ) λ dθ = -log |σ(ϕ)|, for all ϕ ∈ [0, 2π[, then it holds u(z) = 2π 0 q(θ) log σ(θ)-z λ
dθ since this integral is clearly harmonic and is equal to -log |z| on ∂W (A). It is known that such a q exists if and only if λ is different from the logarithmic capacity of W (A). Generally we will use this equation with λ = 1 and then rewrite it as

2π 0 q(θ) log σ(θ)-σ(ϕ) e iθ -e iϕ dθ + 2π 0 q(θ) log |e iθ -e iϕ | dθ = -log |σ(ϕ)|, ∀ ϕ ∈ [0, 2π[.
We discretized this equation using a representation σ(θ) of ∂W (A) and approximating q(•) by a trigonometric polynomial q n (•) of degree n, and employing a collocation method at the points θ j , j = 0, 1, . . . , 2n (it is known that an odd number of collocation points is necessary for such a method). So, we get an approximation q j = q n (θ j ) by solving the system 2π 2n + 1

2n j=0 q j log σ(θ j ) -σ(θ i ) e iθ j -e iθ i + 2π 0 q n (θ) log |e iθ -e iθ i | dθ = -log |σ(θ i )|,
for i = 0, 1, . . . , 2n.

We have approximated the first integral by the trapezoidal formula; of course, if j = i, we have to replace log σ(θ j )-σ(θ i ) e iθ j -e iθ i by log |σ (θ i )|. Recall that, for the remaining integral, there holds

2π 0 q n (θ) log |e iθ -e iθ i | dθ = - 2π 2n + 1 2n j=0 c(j -i) q j , with c(k) = c(-k) = n j=1 cos jθ k j .
Then, we obtain the approximation of u from u(z) 2π 2n + 1 2n j=0 q j log |σ(θ j ) -z| and the approximation of the derivatives of g at 0 (note that here v q j σ(θ j ) k ) .

(x) = 0 if x ∈ R) g (0) = exp(u(0)), g (0) = 2 g (0)u (0), g (3) (0) = 3 g (0)(u (0) 2 +u (0)), g (4) (0) = 4 g (0)(u (0) 3 +3u (0)u (0)+u (3) (0)), g (5) (0) = 5 g (0)(u (0) 4 +6u (0) 2 u (0)+4u (0)u (3) (0)+3u (0) 2 +u (4) (0)),
Remark. In order to get q, we have to solve a linear system of the form M q = b. But it could appear that the matrix M is not invertible, or badly conditioned, if the logarithmic capacity of W (A) is close to 1. In this case, we can replace this system by (M -E)q = b-e where E (resp. e) is a matrix (resp. a vector) with all entries equal to 1.

This method is very efficient for analytic boundary (exponential convergence with respect to n, see for instance [START_REF] Cheng | The delta-trigonometric method using the single layer potential representation[END_REF]), but here we have singularities at the transition points between the straight line and the algebraic part, which reduces the order of convergence to O(n -4 ) which is still good. With ϕ k (z) = 1 k k-1 j=1 j z i(k-j) , we have

∂W (A k ) = {z : z = ϕ k (θ), -2π k ≤ θ ≤ 2π k } ∪ {z : z = -1 2 (1+iy), -cot( π k ) ≤ y ≤ cot( π k )}.
We have used 2n+1 points on the algebraic part of ∂W (A k ): z j = ϕ k ( 2πj kn ) for j = -n, . . . , n and 2n 2 equidistant points on the straight part z n+j = z n -jh i, j = 1, . . . , 2n 2 where h = 2 cot( π k )/(2n 2 +1) and n 2 is chosen such that h is as close as possible to |z n -z n-1 |.

5 Program in Matlab for the computation of g(A)

function [gofA,gderivs,nn]= Akstudy(k,n); % For 3 <= k <= 6, forms the kxk matrix A with ones in the % strict upper triangle and zeros elsewhere, and computes % its image gofA under the Riemann mapping from W(A) to the % unit disk with g(0) = 0, g'(0) > 0. gderivs(j), j=1,...,5,

Figure 1 :

 1 Figure 1: The boundary of the numerical range in black, the remaining part of the algebraic curve in dashed blue

Figure 2 :

 2 Figure 2: The boundary of Ω in red, of W (A) in black.

Figure 3 :

 3 Figure 3: Zoom close to the straight line.

3π 4 ≤Figure 4 :

 44 Figure 4: Curve p(f 1 (e iθ )), 0 ≤ θ ≤ 3π 4 .

Figure 6 :

 6 Figure 6: Curve Re(f 1 (e iθ )), 3π 4 ≤ θ ≤ π.

Figure 7 :

 7 Figure 7: Curve d dθ Re(f 1 (e iθ )), 3π 4 ≤ θ ≤ π.

Figure 8 :

 8 Figure 8: W (A 4 ).Figure 9: W (A 5 ).Figure 10: W (A 100 ).

Figure 9 :

 9 Figure 8: W (A 4 ).Figure 9: W (A 5 ).Figure 10: W (A 100 ).

Figure 10 :

 10 Figure 8: W (A 4 ).Figure 9: W (A 5 ).Figure 10: W (A 100 ).

A 4 ,

 4 we have computed the values g (0) = 1.1888506, g (0) = -1.6292742, g (0) = 4.7085601, which gives g(A 4 ) = a = 1.1888506, b = 0.3742134, c = 0.3443362.

,

  with a = 1.0798634, b = 0.1590093, c = 0.1519169, d = 0.1359021 et e = 0.1161184.

q

  j log |σ(θ j )| , u (k) (0) -(k-1

Table 1 :

 1 summarizes our results. Upper and Lower Bounds on ψ(

	n lower bound upper bound difference
	3	1.9956978	1.9956979	10 -7
	4	1.993800	1.993801	10 -6
	5	1.992921	1.992922	10 -6
	6	1.992444	1.992445	10 -6

Thanks to Bernd Beckermann for a useful remark.

Acknowledgment:

The first author conducted this work within the France 2030 framework programme, the Centre Henri Lebesgue ANR-11-LABX-0020-01

% contains the value of the jth derivative of g at 0. % W(A) consists of a cardioid and a vertical line segment, % and 2n+1 points are used to represent the cardioid portion. % Output argument nn is then the total number of discretization % points used to represent the boundary of W(A). Note also % that n may be modified (to make things come out even). M(j,j) = log(zzprim(j)) -c0; for k=j+1:nn, M(k,j) = log(abs((zz(k)-zz(j))/(ee(k)-ee(j)))) -c(k-j); M(j,k) = M(k,j); end; end; % If t < 0, M is badly conditioned, so we translate M. t = 10∧4-cond(M); if t < 0, M = M -ones(M); t, pause, end; % Compute q. b0 = log(abs(zz)); q = -M\b0; % Take account of the translation. if t < 0, b0 = b0 -ones(b0); end; % Compute derivatives of g at 0. gp = exp(sum(q.*b0)); b = -real(sum(q./zz)); c = -real(sum(q./zz.∧2)); d = -2*real(sum(q./zz.∧3)); ed = -6*real(sum(q./zz.∧4));


Valentin Héluain

Grégoire Prévot

Bastien Cabarrou

Olivier Calvayrac

Estelle Taranchon-Clermont

Alain Didier

Sandrine Tabourier-Gouin

Julie Milia

Julien Mazières

Estelle Taranchon- Clermont

Julien Mazieres
email: mazieres.j@chu-toulouse.fr

Clinical and molecular analysis of lung cancers associated with fibrosing interstitial lung disease

Keywords: lung cancer, interstitial lung disease, pulmonary fibrosis, next generation sequencing LC, Lung Cancer, ILD, Interstitial Lung Disease, AEILD, acute exacerbation of ILD, WHO, World Health Organization, CT, Computer Tomography, TNM, Tumors Nodes Metastases, PFS, Progression Free Survival, OS, Overall Survival, DLCO, diffusing capacity for carbon monoxide, IPF, Idiopathic Pulmonary Fibrosis, CTD-ILD, ILD associated with Connective Tissue Disease, SSC, Systemic Sclerosis, UIP, Usual Interstitial Pneumonia, NSIP, Non-Specific Interstitial Pneumonia, CI, Confidence Interval, NR, Not Reached, TMB, Tumour Mutational Burden, Protein Kinases, ERBB2, Erythroblastic Oncogene B2, PI3K/AKT, Phosphoinositide NOTCH1, Notch Homologue

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction:

The link between LC and chronic ILD has been known for many years [START_REF] Meyer | Relationship of interstitial pneumonia honeycombing and atypical epithelial proliferation to cancer of the lung[END_REF]. It is now well established that the existence of ILD poses an independent risk factor for LC despite the frequency of confounding factors such as smoking and age [START_REF] Naccache | Lung cancer and interstitial lung disease: a literature review[END_REF]. This relative risk has been assessed as ranging from 3.5 to 7.5 [START_REF] Naccache | Lung cancer and interstitial lung disease: a literature review[END_REF] in various studies with an estimated incidence of 10% to 20% and a specific mortality of 15% [START_REF] Sato | Long-term results and predictors of survival after surgical resection of patients with lung cancer and interstitial lung diseases[END_REF][START_REF] Sawabata | Japanese lung cancer registry study of 11,663 surgical cases in 2004: demographic and prognosis changes over decade[END_REF]. The issues involved in LC management in such situations are complex and marked by the specific risk of acute exacerbation of ILD (aeILD) and more common adverse effect to anti-cancer treatment [START_REF] Naccache | Lung cancer and interstitial lung disease: a literature review[END_REF]. Current data suggest that the clinical link observed between fibrosing ILD and the onset of LC is supported by the presence of common changes in cell functions involving epithelial-mesenchymal transition and abnormal re-epithelialisation [START_REF] Tzouvelekis | Lung cancer in patients with Idiopathic Pulmonary Fibrosis. A retrospective multicenter study in Greece[END_REF]. The purpose of this work is to both describe a cohort of patients with lung cancer associated with fibrosing ILD and to characterize their molecular profile.

Patients and methods

Our study is an observational, single centre, ambispective cohort study conducted at Toulouse University Hospital. Consecutive cases of LC associated with fibrosing ILD were identified between 2000 and 2018 from clinical databases in the Pulmonology Department, pathological databases in the Pathology Department and from the Regional Cancer Network. The patients enrolled in the study have been diagnosed with pulmonary fibrosis confirmed by thoracic imaging (honeycombing and/or reticulation with peripheral bronchiectasis or bronchiolectasis) regardless of the aetiology, and primary lung cancer diagnosed on the basis of a tumour biopsy sample. Diffuse interstitial pulmonary fibrosis must have been diagnosed before or concomitantly with LC. All primary lung histological subtypes according to the latest WHO Classification of Lung Tumours were accepted. The exclusion criteria were as follows:

unavailable data (no CT imaging or tumour tissue); no evidence of fibrosis on imaging carried out at the time as cancer diagnosis or pulmonary fibrosis being developed after the introduction of first-line treatment; pulmonary metastases of extra-bronchial tumours. The main purpose of this study was the clinical, molecular (main oncogenes, microsatellite instability, tumour mutational burden) and immunohistochemical (expression of PD-L1 and CD8) characterisation of primary lung cancer being developed in the context of pulmonary fibrosis.

Aetiological diagnoses and images of interstitial lung diseases were confirmed following expert appraisals or multidisciplinary discussions in our Institution [START_REF] Lynch | Diagnostic criteria for idiopathic pulmonary fibrosis: a Fleischner Society White Paper[END_REF]. Tumour stage was assessed based on the 8th edition of the TNM Classification of Lung Cancers. Formalin-fixed and paraffin-embedded tumour specimens were extracted from the archives of the Pathology Department of Toulouse University Hospital.

Molecular sequencing

The genetic analysis of tumour tissues was carried out by the Foundation Medicine Inc.

Laboratory using FoundationOne®CDx sequencing (synthesis sequencing, Illumina® HiSeq 4000).

Immunohistochemical analyses

PD-L1 analysis was performed on an automated Bond Polymer Refine Detection System (Leica Biosystems®) after deparaffinisation (clone QR-1, 1-PR292-01, Diagomics®). CD8 analysis was performed with the SP57 RabbitMono clone (Roche®) on the automated Ventana Benchmark Ultra System (Roche®). Results were expressed as a percentage of labelled cells and classified into three categories: 0%, 1-49% and ≥50% respectively.

Statistical analyses

Clinical and biological characteristics were described using standard statistics. Quantitative variables were described by the median and range, and qualitative variables by frequency and percentage. Progression-free survival is defined as the period of time from the start date of the treatment until the disease progression date or date of death from any cause. Living and progression-free patients were censored on the date of the latest progress report or new line of treatment start date. Overall survival is defined as the period of time from the start date of the treatment to the date of death from any cause. Living patients were censored at the date of the latest progress report. Survival data were summarised using the Kaplan-Meier method.

Statistical analyses were performed using STATA v16 software (StataCorp, Texas, USA). The genomic dataset obtained was analysed using the Reactome database (version 72), organised into signalling pathways.

Results

A total of 49 patients diagnosed between 2000 and 2018, were retrospectively analysed between November 2017 and November 2019. Sufficient tumour material was archived to allow both immunohistochemical analysis of 23 tumours and molecular analysis in 14 patients.

The patient selection process is summarised in Figure 1.

Clinical characteristics:

Clinical characteristics are listed in Table 1. The population mainly comprised men (69.4%) with significant smoking exposure in 83.3% of patients. Impaired diffusing capacity for carbon monoxide (DLco) predominated (49%). We mainly found adenocarcinomas (65.3%) and mostly localised stages, I-II (51%).

Concerning ILD, 30.6% of the diagnoses performed were IPF and 22.4% were CTD-ILD, including SSc in 9/11 cases (81.8%). In 20.4% of cases, patients were not investigated for ILD (no pathological diagnosis). CT-Scan mainly revealed usual interstitial pneumonia (UIP) in 53.1% of cases (definite and probable UIP in 28.6% and 24.5% of cases, respectively). Nonspecific interstitial pneumonia (NSIP) and indeterminate UIP accounted for 14.3% and 22.4% of cases, respectively. The primary tumour lesion was located within fibrotic lesions in 77.6% of cases and emphysema lesions were present in 22.4% of cases. In 26% of cases, ILD was diagnosed concomitantly with LC on the basis of a CT-scan. In the remaining cases, a median 3-year follow-up period preceded LC onset.

In the sub-group of non-smoking patients (n=8) the most common diagnosis of ILD was connective tissue disease in 62.5% of cases [SSc, 50% (n=4); rheumatoid arthritis, 12.5% (n=1)], followed by IPF [25% (n=2)].

Therapeutic and progressive characteristics

Forty-one patients (83.7%) received at least one line of treatment for LC. First-line treatment comprised chemotherapy in 58.5% of cases. Radiotherapy was administered to 17.1% of patients and in combination with chemotherapy for one patient (2.4%). Five patients (i.e., 71% of patients undergoing radiotherapy) received stereotactic radiotherapy exclusively. Finally, surgery was proposed in 41.5% of cases depending on individual cases. No patient was offered immunotherapy or targeted therapy as first-line treatment. The response rate to first-line chemotherapy was 38.1% and the disease control rate was 81.0% respectively. A total of five patients (10.2% of patients analysed) developed aeILD during follow-up with a specific mortality rate of 60%. Four occurred immediately after treatment: two cases after surgery, two were iatrogenic causes secondary to chemotherapy and one was idiopathic.

Finally, two cases of grade > II radiation pneumonitis were documented in our cohort, one of which occurred after stereotactic irradiation, with a specific complication rate of 33%. 2b).

The cause of death was known in 50% of the patients who died at the time of the last progress report. 60% of deaths were related to a progressive complication of cancer. 15% were due to a complication of applied cancer treatments (including 2 acute exacerbations of ILD and 1 chronically progressive radiation pneumonitis). Finally, 20% of deaths were due to other causes (infection and cardiovascular problems in particular). Overall, 15% of these deaths were directly related to ILD.

Immunohistochemical characteristics

Twenty-three patients underwent immunohistochemical analysis. 65.2% of the specimens analysed showed no PD-L1 protein expression on tumour cells. Moderate expression (1 to 49%) was documented in 30.4% of patients. Tumour lymphocyte infiltration, as assessed by CD8 labelling, was deemed moderate (labelling between 1 and 49% of cancer cells) in all patients. The results are summarised in Table 2.

Molecular characteristics

Assessed in 14 patients, the tumour mutational burden (TMB) was considered low (0 -5 Mut/Mb) or intermediate (6 -19 Mut/Mb) in 69.2% and 30.8% of the analysed tumours, with a median of 3.78 Mut/Mb, respectively. All assessed patients showed a stable microsatellite status.

At least, one mutation was detected in 112 out of the 324 analysed genes. The most common mutations are shown in Figure 3, and are compared with the TCGA molecular database for lung cancer. TP53 was the most frequently altered gene (71.4%). These mutations were involved in 949 signalling pathways after analysis into the Reactome database. The 25 most significantly impacted signalling pathways are presented in Table 3.

Discussion

We analysed the pathological and clinical characteristics of LC associated with fibrosing ILD.

Adenocarcinoma was the most frequent histological subtype in the assessed population , ahead of squamous cell carcinoma, with percentages comparable to the published data(7)(8).

In the specific IPF subgroup, the squamous subtype was often described as predominant, especially in the basal and peripheral areas where fibrous mutilation is most significant [START_REF] Antoniou | Idiopathic pulmonary fibrosis and lung cancer: a clinical and pathogenesis update[END_REF][START_REF] Tomassetti | The Impact of Lung Cancer on Survival of Idiopathic Pulmonary Fibrosis[END_REF][START_REF] Brown | Idiopathic Pulmonary Fibrosis and Lung Cancer. A Systematic Review and Meta-analysis[END_REF].

In our study, we found a higher proportion (40%) of squamous cell carcinomas in the IPF subgroup with 77% of primary lesions being located in the fibrotic zone.

As previously described(2), it seems that IPF and ILD associated with SSc are were particularly related to LC onset. Only two patients included in our cohort were diagnosed with rheumatoid arthritis with fibrotic ILD and LC. In 2018, in a retrospective cohort of 24 patients presenting CTD-ILD and LC, Watanabe et al. found 54% cases of rheumatoid arthritis compared to only 29% cases of SSc. However, selection bias is extremely likely. Notably, one in five patients (20.4%) presented typical ILD with fibrosis, which was confirmed by imagery. No aetiological diagnosis was given due to insufficient investigations. This may be attributed to the frequency of simultaneous ILD and LC diagnoses in this population (60%) and to the priority given to cancer treatment or a deterioration in the clinical condition requiring palliative management only (in 40% of those patients who underwent no further investigation). UIP (definite or probable UIP) was mainly highlighted on CT-scans. These patients also had a particularly poor prognosis due to treatment-induced aeILD (2).

Overall survival was lower than the values generally reported in the overall population, regardless of the tumour stage [START_REF] Goldstraw | The IASLC Lung Cancer Staging Project: Proposals for Revision of the TNM Stage Groupings in the Forthcoming (Eighth) Edition of the TNM Classification for Lung Cancer[END_REF]. We nevertheless note major differences between tumour stages. Half of LC cases associated with ILD were diagnosed at a localised stage (stages I and II). This difference can be explained by the regular scan follow-up implemented in ILD management with fibrosis. In this population, the survival rate for surgically corrected stages 1 and 2 (34.5 months) is satisfactory and encourages early diagnosis and management for these patients. Regarding metastatic stages, the objective tumour response rate and the rate of disease control under chemotherapy were comparable to those observed in the general population based on previously published data [START_REF] Gandhi | Pembrolizumab plus Chemotherapy in Metastatic Non-Small-Cell Lung Cancer[END_REF][START_REF] Paz-Ares | Pembrolizumab plus Chemotherapy for Squamous Non-Small-Cell Lung Cancer[END_REF]. One notable feature is the overlap of the progression-free survival curves during first-line treatment in stage III and IV LC patients.

This may be attributed to the under-representation of radiotherapy in the management of stage III LC associated with ILD. This is understandable given the much higher risk of aeILD following radiation pneumonitis in this population. In accordance with published data(2), tumour progression was the primary cause of death in our cohort, thus, it highlights the interest in specific treatments for LC associated with fibrosing ILD.

The data observed on PD-L1 expression on tumour cells, tumour lymphocyte infiltrate, tumour mutational burden or microsatellite instability in our cohort are consistent and indicative of anticipated low immunogenicity in these tumours. No prospective data are currently available to assess the efficacy and safety of immunotherapies in the context of fibrosing ILD. In a retrospective analysis of 16 patients with LC related to fibrosing ILD, Guyard et al. found that PD-L1 expression was absent (<1%) or moderate (between 1 and 49%) in 37.5% and 56% of patients respectively [START_REF] Guyard | Morphologic and molecular study of lung cancers associated with idiopathic pulmonary fibrosis and other pulmonary fibroses[END_REF]. These observations are consistent with a different distribution in terms of PD-L1 expression presenting a lower proportion of tumours strongly expressing PD-L1 (>50%) in this population [START_REF] Guibert | Targeting the immune system to treat lung cancer: rationale and clinical experience[END_REF]. Contrastingly, squamous cell carcinomas associated with fibrosing ILD appeared to have significantly higher PD-L1 expression with a lower T-CD8 infiltrate [START_REF] Ueda | Immunosuppressive tumor microenvironment of usual interstitial pneumonia-associated squamous cell carcinoma of the lung[END_REF]. It should, however, be noted that there is currently no single biomarker predictive of immune checkpoint inhibitor activity, and these observations alone do not allow the efficacy of immunotherapy to be determined in this context. In the absence of a biomarker, their use is still not recommended, especially since recent observations show a high rate of severe pulmonary toxicity herein (52).

No molecular abnormality responsible for classic oncogene addiction (EGFR, ALK, ROS1, BRAF for example) was highlighted during the molecular analyses although some patients displayed gene mutation involving KRAS or ERBB2, for example -potential targets for therapies under clinical evaluation [START_REF] Sankar | Molecular therapeutic targets in non-small cell lung cancer[END_REF].

Gene TP53 was the most commonly mutated one within our analysis (in 71.4% of patients). This is comparable to published LC data [START_REF] Zhu | PPM10 Correlation and hiearchical clustering order analysis of top 25 co-occurring genomic alterations in non-small cell lung cancer[END_REF]. TP53 mutations also seem to be linked, in particular, to smoking-related LC [START_REF] Guyard | Morphologic and molecular study of lung cancers associated with idiopathic pulmonary fibrosis and other pulmonary fibroses[END_REF]. It is interesting to note that changes in TP53 and KRAS characterise patients presenting IPF, in tissues devoid of tumour proliferation [START_REF] Tzouvelekis | Lung cancer in patients with Idiopathic Pulmonary Fibrosis. A retrospective multicenter study in Greece[END_REF][START_REF] Otsubo | Paired genetic analysis by nextgeneration sequencing of lung cancer and associated idiopathic pulmonary fibrosis[END_REF].

No mutations in EGFR gene were found, although it is frequently involved in about 10% of BPCs. This observation can probably be explained by the high proportion of male smokers in our cohort but has already been described in patients with IPF-related LC where the frequency of EGFR mutations is significantly lower [START_REF] Kanaji | Impact of idiopathic pulmonary fibrosis on advanced non-small cell lung cancer survival[END_REF][START_REF] Kojima | The pathological features of idiopathic interstitial pneumonia-associated pulmonary adenocarcinomas[END_REF]. Similarly, no ALK-EML4 translocation was detected. The frequency of this mutation does not appear to differ in patients presenting idiopathic ILD according to already published data [START_REF] Kojima | The pathological features of idiopathic interstitial pneumonia-associated pulmonary adenocarcinomas[END_REF], but no firm conclusion can be drawn given that it is relatively rare (approximately 4% of ADC in the general/overall population).

Three isolated ALK mutations with unknown prognosis have nevertheless been identified.

Signalling pathway analysis provides quantitative results, highlights the level and extent of changes within the affected pathways. PI3K/AKT pathway was thus altered to considerable extent. It is known to be involved in fibroproliferative diseases, with a role in the proliferation of human pulmonary fibroblasts, and in cancer by ensuring cell survival [START_REF] Conte | PI3K p110 γ overexpression in idiopathic pulmonary fibrosis lung tissue and fibroblast cells: in vitro effects of its inhibition[END_REF]. TGFβ pathway (Diseases of signal transduction by growth factor receptors and second messengers, Table 3), a major pro-fibrotic signalling pathway and promoter of epithelial-mesenchymal transition [START_REF] Tzouvelekis | Lung cancer in patients with Idiopathic Pulmonary Fibrosis. A retrospective multicenter study in Greece[END_REF], was greatly impacted. This pathway also plays an important role in tumour proliferation (promitotic effect) and progression (crossing the basal membrane and metastasis) and constitutes a research target for the development of NSCLC treatments [START_REF] Eser | TGFβ pathway inhibition in the treatment of non-small cell lung cancer[END_REF]. These observations, which are consistent with literature, suggest that the development of LC in a pulmonary fibrosis setting may be supported by the involvement of common biological mechanisms. An initial randomised trial to assess the merits of Nintedanib and Nab-paclitaxel combination therapy versus paclitaxel alone is on-going [START_REF] Otsubo | Paired genetic analysis by nextgeneration sequencing of lung cancer and associated idiopathic pulmonary fibrosis[END_REF]. This is the first randomised trial conducted specifically in this population. There may also be a prophylactic effect with a documented decrease in the incidence of LC in IPF treated with Pirfenidone [START_REF] Miura | Reduced incidence of lung cancer in patients with idiopathic pulmonary fibrosis treated with pirfenidone[END_REF]. A pre-clinical study demonstrated that blocking PI3K/AKT pathway prevented bleomycin-induced pulmonary fibrosis in rats [START_REF] Wei | A phosphoinositide 3-kinase-γ inhibitor, AS605240 prevents bleomycin-induced pulmonary fibrosis in rats[END_REF].Omipalisib, a specific inhibitor of this pathway is currently being assessed in the treatment of IPF with a recently published phase 1 trial demonstrating acceptable drug tolerability in 17 patients and confirming pharmacodynamic activity on the basis of reduced tomographic fixation of pulmonary fibrosis zones [START_REF] Lukey | A randomised, placebocontrolled study of omipalisib (PI3K/mTOR) in idiopathic pulmonary fibrosis[END_REF]. Finally, Wnt/β-catenin pathway (signal transduction, Table 3), featured in our analyses, is also known to be involved in the development of fibrosis in numerous organs. Specific blocking of Wnt/β-catenin pathway by a monoclonal antibody (OMP-131R10) is the subject of research with interesting pre-clinical results in pulmonary fibrosis models [START_REF] Zhang | Targeting the Wnt signaling pathway through R-spondin 3 identifies an anti-fibrosis treatment strategy for multiple organs[END_REF].

The main interest of this work is the clinical, immuno-histochemical and molecular description of LCs which was developed in the fibrotic ILD context on broad patient selection criteria facilitating analysis of real-life clinical practices in this area. To our knowledge, this is the first study combining immunohistochemical and molecular analyses with such a large cohort. However, several biases exist and unequivocal conclusions cannot, therefore, be drawn. There is a selection bias and, in particular, a "centre effect" one given that this is a single-centre study. There is also information bias in the data collection due to retrospective methodology, which is however more appropriate in this rare disease context. Finally, the absence of a control group, especially for molecular analyses does not allow any conclusion to be drawn with regard to the specific nature of the anomalies detected. Furthermore, the poor number of patients does lead to confusion bias.

Conclusion

These observations suggest a specific phenotype of LC associated with fibrosing ILD. We found a clear majority of patients with idiopathic pulmonary fibrosis or interstitial lung disease associated with SSc, two conditions with a major risk of pulmonary neoplastic complications.

These patients have a poor prognosis regardless of the tumour pathology stage. Nevertheless, our results demonstrate the merits of early detection and management of these neoplastic complications given the high, specific mortality rate associated with cancer. Overall, the predictive markers used to analyse the efficacy of immunotherapies do not suggest that such treatments are beneficial for the study population. The genomic study highlights an over-representation in the tumour specimens of changes in signalling pathways known to be involved in pulmonary fibrosis. This observation suggests a continuum in the development of pulmonary fibrosis towards cancer, and may provide a basis for upcoming specific therapeutic strategies.

Declaration of Competing Interest

All authors have completed the disclosure form. Julien Mazières benefits from funding Roche FMI tests. The other authors have no conflicts of interest to declare.

4 .

 4 The signalling pathways shared by at least three patients are shown in Figure Abnormalities in proliferation signal transduction by mutation of growth factor receptors or second messengers mainly concerned TGF-β (transforming growth factor beta) pathway as well as PDGFR (platelet derived growth factor receptor) signalling pathways, FGFR (fibroblast growth factor receptor), EGFR (epidermal growth factor receptor, or ErbB1), Wnt, MAPK (mitogen-activated protein kinases), ErbB2 (erythroblastic oncogene B2, or HER2), PI3K/AKT (phosphoinositide 3-kinase/protein kinase B), and NOTCH1 (Notch homologue 1).

Acknowledgements

Part of this work was supported by ROCHE SAS by funding FMI tests. ROCHE SAS had no other involvement in the design and the conduct of the research nor the preparation of the article.


Antoine Sabrié
email: antoine.sabrie@ifpen.fralexandre.battiston@ifpen.fr

Alexandre Battiston

Fabien Vidal-Naquet
email: fabien.vidal-naquet@ifpen.fr

Jean-Yves Gauthier
email: jean-yves.gauthier@insa-lyon.frxuefang.shi@insa-lyon.fr

Xuefang Lin-Shi

Analytical Derivation of Output Distortions Harmonics Caused by Voltage Source Inverter Non-Linearities in PMSM Drives (1,2)

Index Terms-Voltage Source Inverter (VSI), Harmonics, Dead-time, Parasitic elements, Wide bandgap Abstract-This paper presents a new direct analytical derivation of low frequency harmonics introduced by voltage source inverters non-linearities in PMSM drives. The proposed voltage drops model is robust and requires fewer calculations than time-based models, and has proven to accurately estimate the experimental voltage drops for SiC based inverter output waveforms. Furthermore, it has been extensively validated experimentally and offers an analytical tool to precisely estimate the voltage drops introduced by the inverter non-linearities.

I. INTRODUCTION

Voltage source inverter fed Permanent Magnet Synchronous Machine (PMSM) coupled with an adequate Pulse-width modulation (PWM) has proven to be a very efficient solution. This is due to its high efficiency and passive filtering of high order harmonics inducing low Total Harmonoc Distorsion (THD) in the output waveforms. However, the inverter is not ideal and can introduce a certain number of voltage drops causing undisered non-linear behavior. Such as the voltage drops of the switches and diodes, the required dead-times, or the turn-on and turn-off delays of the semiconductors. These added non-linearities introduce low frequency harmonics which are difficult to filter out and are at the origin of low frequency current ripple in the machine if they are not compensated by the controller. This induces in turns undesired torque ripple and impacts efficiency.

Numerous work in the literature aims at compensating these undesired effects in the PMSM drives [START_REF] Ji | Control strategies of mitigating dead-time effect on power converters: An overview[END_REF]. The compensation techniques can be based on instantaneous average voltage compensation, as proposed in [START_REF] Lee | A simple and direct dead-time effect compensation scheme in pwm-vsi[END_REF], [START_REF] Zhang | Dead-time compensation of inverters considering snubber and parasitic capacitance[END_REF]. One of the main drawback is that the compensation scheme is dependent on the polarity of the current. Thus very difficult to accurately compensate in nearzero current operation. In [START_REF] Seyyedzadeh | Accurate modeling of the nonlinear characteristic of a voltage source inverter for better performance in near zero currents[END_REF] the authors proposed an analysis and compensation scheme based on look-up tables to improve performances in zero crossings regions. Furthermore, a model-based compensation will never truly mitigate all voltage drops due to parameters identification. An on-line estimator based on the voltage drops is proposed in [START_REF] Urasaki | An adaptive dead-time compensation strategy for voltage source inverter fed motor drives[END_REF] but suffers from the same drawbacks as the average value compensation. Advanced algorithms based on observers aims at compensating the given 6th order harmonics in the dq rotary reference frame [START_REF] Park | Parameter-independent online compensation scheme for dead time and inverter nonlinearity in ipmsm drive through waveform analysis[END_REF], [START_REF] Hao | An smc-eso-based distortion voltage compensation strategy for pwm vsi of pmsm[END_REF]. These methods achieve good compensation in steady state but exhibit poor performances in dynamic responses. As proposed in [START_REF] Wang | A precise voltage distortion compensation strategy for voltage source inverters[END_REF], in which a model-based is coupled with a current observer to compensate voltage drops in both steady state and dynamic responses. But the study does not considerate the parasitic capacitance, which can have a huge impact on voltage drops [START_REF] Wang | Influences on output distortion in voltage source inverter caused by power devices' parasitic capacitance[END_REF]. Methods based on models alone cannot truly compensate the voltage drops. In the other hand, advanced disturbance observer, or repetitive control which are not based on models can accurately compensate in the steady-state but not in dynamic responses.

In this paper a simple analytical model of the voltage drop is proposed based on a frequency analyisis of the voltage drops and undisred non-linear behavior. To the author's knowledge, direct harmonics calculation is complex and requires high computational resources, as can be described for the switching frequency current har-monics derived in [START_REF] Liang | Analytical modeling of current harmonic components in pmsm drive with voltage-source inverter by svpwm technique[END_REF]. In order to improve the dynamic response of such complex compensation algorithm, this paper aims at proposing a simple analytical derivation of the low frequency harmonics of the voltage drops. Furthermore, in the litterture only partial validation of the models and compensation scheme are done on specific operating points. In this paper, a thorough validation in all accessible operating point is done to remove any bias.

The paper is structured as follows : In section II a more complete voltage drop model based on various previous works is presented, then in section III the analytical derivation of the low-order harmonics will be done. Finally, simulations and experiments to validate the robustness of the proposed model are presented in section IV and, section V respectively.

II. NON-LINEARITIES MODELLING IN VSI

Thi section focuses on the development of the voltage drop average model over a swiyching period T sw . It will be done considering an equivalent half bridge as presented in figure Fig. 1. The pole voltage V ao can be expressed by the switched drain-source voltage across the low side MOSFET of the half bridge. The ideal perfect pole voltage V ao in this condition is defined by equation [START_REF] Ji | Control strategies of mitigating dead-time effect on power converters: An overview[END_REF]. Where V DC is the DC bus voltage and C a (t) is the gating signals of the MOSFET.

V ideal ao = V DC • C a (t) (1)
The Fig. 2 summarizes all the considered the voltage drops sources introduced referenced to the ideal voltage. In comparison, the Fig. 3 presents the voltage ∆V ao = V ao -V ideal ao obtained at the output of the half bridge when all presented non-linearities are taken into account.

A. Voltage drops due to dead-times and switching delays

Dead-times are added between the gating signals of the two MOSFETs of the half bridge to prevent any short-circuit of the DC bus voltage. The error introduced in the phase voltage is presented in Fig. 2 in the fourth sub-graph. The dead-times will induce a delay of the establishment of the pole voltage. Since the ideal voltage is the switched voltage corresponding to the DC bus voltage, the error introduced has an amplitude of V DC during the dead time. Due to the way of generating the dead-times, the error will change sign depending on the current polarity. The voltage drop introduced by deadtimes is thus expressed as [START_REF] Lee | A simple and direct dead-time effect compensation scheme in pwm-vsi[END_REF]. Another non-linearity introduced is due to the switching delays at turn-on and turn-off, as presented in Fig. 2 on the third sub-graph.Normally it should vary with different parameters such as the peak phase current I peak , V DC or the temperature, but it would complicate the model and has a lesser impact than the dead-times. It is possible to calculate the error induced by turn-on and turn-off delays as expressed in [START_REF] Zhang | Dead-time compensation of inverters considering snubber and parasitic capacitance[END_REF].

∆V DT ao = T dead T sw • V DC • sign(i a ) (2)
∆V sw ao = t on -t of f T sw • V DC • sign(i a ) (3)
Due to the very close similarity with the dead-time, it is easier to express all of them in the same equation as proposed in [START_REF] Seyyedzadeh | Accurate modeling of the nonlinear characteristic of a voltage source inverter for better performance in near zero currents[END_REF].

∆V DT +sw ao = T dead + t on -t of f T sw • V DC • sign(i a ) (4)

B. Voltage drops due to conduction

The figure Fig. 1 presents the active component in which the current is flowing. As it can be seen it depends on the polarity of the current, the gating signal C a (t), and if takes place during the dead-time. Thus, defining antiparallel diode V DS and V d expressed in [START_REF] Urasaki | An adaptive dead-time compensation strategy for voltage source inverter fed motor drives[END_REF].

V DS = R on DS • i a V d = V do + r d • i a (5)
Where R on DS and r d are the ON-state resistances. One proposed model in [START_REF] Zhao | An accurate approach of nonlinearity compensation for vsi inverter output voltage[END_REF] is done by differentiating the actual time when the MOSFET is ON or OFF from the ideal one. Let us define C a (t) the duration during which the high side MOSFET of phase A is ON during the switching period T sw By definition, the duty cycle is given by [START_REF] Park | Parameter-independent online compensation scheme for dead time and inverter nonlinearity in ipmsm drive through waveform analysis[END_REF].

C a (θ) = 1 T sw Tsw 0 C a (t)dt (6)
Where θ = ωt. If C * a (θ) is the ideal duty cycle that should be applied, then the difference between the real time and the ideal is due to dead-times and the switching delays and can be expressed as [START_REF] Hao | An smc-eso-based distortion voltage compensation strategy for pwm vsi of pmsm[END_REF].

C a (θ) = C a * (θ) -(T dead + t on -t of f ) • sign(i a ) (7)
The voltage drop from the ideal applied voltage can be extensively described in each conduction case by ( 8)

∆V on ao = (C a (θ)R on DS i a + (1 -C a (θ))(V do + r d i a )) • ( 1 + sign(i a ) 2 ) + ((1 -C a (θ))R on DS i a -C a (θ)(V do -r d i a )) • ( 1 -sign(i a ) 2 ) (8)

C. Voltage drops due to Parasitic capacitance

The effect of the parasitic capacitance on output voltage distortions has been studied in previous work solely. The parasitic capacitance plays a role in the switching waveforms during dead-time and switching transitions. The model presented was derived from [START_REF] Wang | Influences on output distortion in voltage source inverter caused by power devices' parasitic capacitance[END_REF]. During a turn-on sequence, for example, the voltage V DS is initially equal to V DC and the parasitic capacitance is charged. If the magnitude of the phase current at the given instant is high enough, the capacitance is very quickly discharged and no effect can be noticed. However, when the magnitude is close enough to zero, the discharge of the parasitic capacitance is much slower and holds its voltage value V DC for a certain amount of time until the end of the dead time when the voltage is clamped to zero.The time needed to discharge the parasitic capacitance is stated in [START_REF] Wang | Influences on output distortion in voltage source inverter caused by power devices' parasitic capacitance[END_REF], where the factor 2 appears because the current needs to discharge the two parasitic capacitance in series from the high and low side MOSFET. And with the hypothesis that the current is fairly distributed among the two. The two last subgraphs of Fig. 2 describes the voltage drop associated.

∆t = 2C oss • V DC i a (9)
From the defined ∆t, we can differentiate if the capacitor has enough time to discharge during the deadtime. This condition is expressed if the magnitude of the phase current i a is enough to discharge the capacitor quickly enough. Thus, drawing the expression of the threshold current given in equation ( 10)

I thr = 2C oss • V DC T dead (10)
The voltage error created is actually opposite to the voltage drops due to the dead-time. Meaning, in one switching cycle, A positive voltage error can be created from the dead-time and a negative one from the parasitic capacitance. The actual voltage error can be derived from the difference of the two. In fact, the simple model of the dead-time and switching delays tend to overestimate and by far the true level of voltage drop. The total voltage drop due to the dead-time taking into account the parasitic capacitance are expressed in [START_REF] Zhao | An accurate approach of nonlinearity compensation for vsi inverter output voltage[END_REF] and ( 12) depending on the magnitude of the current.

When i a ≥ I thr

∆V DT +Coss ao = 1 T sw (T dead + t on -t of f ) • V DC • sign(i a ) - C oss V 2 DC i a (11)
When i a ≤ I thr

∆V DT +Coss ao = i a • (T dead + t on -t of f ) 2 4C oss T sw (12)
The complete averaged model of the voltage loss considering all non-linearites sources is defined in (13) :

∆V ao = ∆V on ao + ∆V DT +Coss ao (13)

III. DIRECT ANALYTICAL DERIVATION OF THE LOW FREQUENCY HARMONICS

The complete model above describes the total voltage drops due to the inverter non-linearities. In this section a frequency analysis of the time-based averaged model will be done in steady state. Thus an anlytical way of calculating the harmonics introduced by the voltage drops will be done. Since the harmonics to compensate can be known in advance a quick analytical derivation of such harmonics could be done.

A. Voltage drops due to conduction

For this purpose, the voltage error ∆V ao needs to be described in the frequency domain based on the knowledge known on the time-based one. The objective is to define from the complete time varying model equation as ∆V ao stated in (13). The proposed method in this paper is to replace the instantaneous values by their Fourier coefficient regarding the desired harmonic. For simplicity and due to the experimental setup, most of the work presented in the following are for the fundamental harmonic. Nevertheless, the procedure is the same for low voltage harmonics like the 5 th , the 7 th , 11 th or 13 th . The harmonic of higher ranks will not be considered, since their amplitude will be very low and their impact negligible. The switching harmonics will also not be derived and would require a lot more complexity in the model and since there are naturally filtered by the load, they have a lesser impact on the torque ripple.

The ideal phase currents references are supposed perfectly sinusoidal and are described in (14).

i a (θ, ϕ) = I peak • sin (θ -ϕ) (14)
The polarity of the current is a term that appears also in the model. Since it is a square wave of amplitude one, its Fourier series is given in (15). F(sign(i a (θ, ϕ))) = 4 π k=1,3,5,7,9,11,13

1 k sin (kθ) (15)
The duty cycle in the case of a Space Vector Pulse Width Modulation (SVPWM) for half the fundamental period is described as follows (16)

C a (θ) =    1 2 (1 + m cos(θ -π 6 )), θϵ[0; π 3 ] 1 2 (1 + m √ 3 cos(θ)), θϵ[ π 3 ; 2π 3 ] 1 2 (1 + m cos(θ + π 6 )), θϵ[ 2π 3 ; π] (16)
The Fourier series of this signal is given in (17).

C a (θ) = 1 2 + m √ 3 sin(θ) + m π 2 ∞ k=3,5,7,.. 1 k 2 sin(3kθ) (17)
From ( 8) a quick analysis of the developpement is made. Then retaining only the terms that have an impact on the fundamental harmonic in [START_REF] Wang | A precise voltage distortion compensation strategy for voltage source inverters[END_REF], the fundamental harmonics of the voltage drop due to conduction can be expression as (18):

∆V on ao (θ) = R on DS + r d 2 i a + R on DS + r d 2 i a sign(i a ) + V d 1 + sign(i a ) 2 + (R on DS -r d )C a (θ) • i a • sign(i a ) -C a (θ) • V do (18)
Considering all the derived harmonic contribution of each terms, the fundamental voltage drop harmonics can be derived. The cross product C a (θ) • i a • sign(i a ) can be derived using the classical description of the Fourier series given in ( 19)-(21).

f (θ) = a 0 + n a n cos( nπθ T ) + n b n sin( nπθ T ) (19)
a n = 1 T T -T f (θ)cos( nπθ T )dθ (20) b n = 1 T T -T f (θ)sin( nπθ T )dθ (21)
Finally, considering the phase of each term, neglecting the dead-time effect on conduction which is in the considered case negligible, a simple analytical description of the fundamental harmonic of the voltage drop introduced by the conduction of semiconductors is given by (22).

F 1 (∆V on ao ) = R on DS + r d 2 I peak + 4 π V d 2 - m √ 3 • sin(ϕ) • V do + (R on DS -r d ) 2 3π I 2 peak m 2 (3 sin 2 ϕ + 1) (22)

B. Voltage drops due to dead-times considering parasitic capacitance and switching delays

Regarding the voltage drops due to the dead time considering the parasitic capacitance effect, the instantaneous function given in [START_REF] Zhao | An accurate approach of nonlinearity compensation for vsi inverter output voltage[END_REF] and (12) contains conditions. Thus, a symbolic calculation of the Fourier coefficient on the different piecewise description has been done. The function being uneven, only the b n coefficient requires calculation and is done from (21). The result of the symbolic integration of the fundamental harmonic is given in (23). Thus, an equation for the direct analytical derivation of the fundamental component is proposed. IV. SIMULATION RESULTS A simulation of a PMSM fed by a three-phase inverter has been carried out in MATLAB Simulink and the SymPowerSystem library. The machine is a traction machine for drives application. The parameters are the one used experimentally and are given in Tab. I. In order to validate the robustness of the approach, the analysis has been done not only on the optimal point or trajectory usually used but on the entire accessible operating point of the machine, even non-optimal ones. For this purpose, a mapping of all couple Direct and Quadratic current i d and i q have been done. The machine is connected in a back to back configuration and sees its mechanical shaft speed imposed at Ω = 1000rpm. The implemented control scheme for the simulation is presented in Fig. 4. The two references of i * d and i * q are thus given by (24).

F b1 (∆V DT +Coss ao ) = I peak T dead 2 4C oss T sw + 3V DC I peak 2 T dead 2 -4C
i * d = I peak sin(δ) i * q = I peak cos(δ) (24)
In order to validate the proposed analytical and direct derivation of the harmonics, a comparison is carried out between the simulation estimated voltage drop referenced to V ideal ao and the on-line calculation of the fundamental component. The time-based voltage drops estimated by the simulation needs to be post processed to obtain the first harmonic. The difference in voltage drop ∆V ao first harmonic between the voltage drop obtained in simulation and the simplified under hypothesis proposed first harmonic calculation is presented in Fig. 5. It can be seen that the derived harmonic calculation is very close to the harmonic obtained with the FFT of the voltage drop obtained in simulation, with less than at most 0.4V of error in the estimation of ∆V ao under worst case scenario in near zero currents and abrupt change in polarity. Thus, simplified direct fundamental analytical derivation under hypothesis engenders small errors. This proves that a simpler derivation of the fundamental harmonics of the voltage drop can be accurately estimated, knowing only basic datasheet parameters of the inverter and simple measurements of the current on any operating point.

V. EXPERIMENTAL RESULTS

In order to validate the results experimentally, the protocol carried was as follows. The output of the PI controllers, meaning the references voltages, V * d and V * q were measured. Furthermore, the three-phase currents I abc and well as all line-to-line voltages U ab , U bc and U ca were measured with a wattmeter giving their respective RMS and mean values. Since the real V d and V q as well as the real i d and i q applied to the machine are needed for the controlled, were also measured and saved. Thus, a database has been created with all needed mean values applied experimentally. The same protocol described in the section IV was carried out for all operating points. Meaning, a mapping of i dq have been done experimentally. Since the experimental result have been done in closed-loop control, the reference values of V * d and V * q will be higher than the real applied V d V q to attain the given reference steady-point. This means that the control will mitigate at least in mean values the voltage drops. In steady-state operation, it is possible to define the real applied V d V q and the references V * d and V * q . Thus, the difference of the control signals references V * d and V * q and the real V d and V q can be used to estimate experimentally the voltage drops due to the inverter nonlinearities as espressed in (25). Finally, to estimate the voltage drop, the norm ∆V 2 d + ∆V 2 q is calculated.

V * d -V d = ∆V d V * q -V q = ∆V q (25)
The obtained mapping, in the i dq plane, of the estimated voltage drops ∆V 2 d + ∆V 2 q experimentally, the simulation and the analytical derivation are presented in Fig. 6 and shows the absolute level of voltage drop estimated. In order to validate the proposed analytical derivation, the contour of error between the experimental voltage drop and the proposed fundamental derivation is plotted in Fig. 7. It can be seen that the derived direct harmonic calculation presents good performances of estimation of the actual experimental voltage drop with less than 0.4V in difference for the majority of operating points in the i d i q map except in near zero currents region. This is due to the loss of precision of the current polarity detection as well as the zero current clamping phenomenon due to the dead time. Furthermore, the intrinsic parameters extracted from the datasheet can not be perfectly estimated and presents variation in real application. The temperature dependency for instance is partially taken into account. This is one limit of the model. However, although not shown here, the added parasitic capacitance effect make the estimation far more precise compared to models that do not take it into account. Furthermore, the error between the models and the experimental could be explained by other factors such as the parametric error that could be also compensated by the control. Meaning that the method to derive the voltage loss experimentally can be improved. But the presented protocol can be achieved with any motor easily without added measurements. Nevertheless, the experimental results show that the analytical derivation of the fundamental component of the voltage drop can approximate the real voltage drop with a difference less than a volt.

One of the advantages of such analytical voltage drops model is that it can be used to study the impact of each parameters on the voltage drops. For demonstration purposed,the purple map plotted in Fig. 6 represents the same voltage drop model without considering the intrinsic capacitance, as is done in most dead-time averaged models. This illustrates the importance of taking it into account, otherwise the compensated voltage drop will be largely overestimated. It also concord with the observation made in [START_REF] Wang | Influences on output distortion in voltage source inverter caused by power devices' parasitic capacitance[END_REF], that the intrinsic capacitance tends to naturally attenuates the voltage drop introduce by dead-times. This also shows that the model can be easily used as a tool to evaluate the individual impacts of each parameters.

VI. CONCLUSION

In this paper, an extensive instantaneous model timebased models based on previous literature work have been proposed. Based on this model, this paper proposes an analytical derivation of the low frequency harmonics that are the harder to filter out and impact the most in the load with current and torque ripple. This paper has made the statement that only the amplitude of the harmonics of the voltage drops are truly necessary to describe the non-linearities introduced by the VSI instead of time based version previously proposed in the literature. The analytical direct derivation of the fundamental harmonic show very good concordance with the time-based model harmonics as well as the one obtained by simulation model of the inverter non-linearities, thus showing the robustness of the approach. Then the analytical harmonics are used to estimate the voltage drops of introduce by a VSI experimentally on all accessible operating points at given speed and DC bus voltage. The results show that the model offers good performance and can accurately estimate the true voltage experimental voltage drop. Errors in the approach are nonetheless engendered, but are due to parametric estimation of the intrinsic parameters that are still needed, as well as the wellknown difficulties of precise current polarity detection and measurement. The proposed work offers a complete voltage drops model that is quite accurate with reduced computation time. Furthermore, the analytical computations can offer a pre-study of the voltage drops for any given inverter. Since the model has been extensively validated both by simulation and experimentally. In order to make the most of this model, it could be paired with advanced disturbance observers to improve dynamic performance of such proposed advance compensation techniques.

Fig. 1 :

 1 Fig. 1: Four principle conduction cases depending on current polarity and gating signals: (a) -for positive current polarity, (b) -for negative current polarity

Fig. 2 :

 2 Fig. 2: Contribution of each considered non-linearity, referenced to the ideal phase voltage (a) -for positive current polarity, (b) -for negative current polarity

Fig. 3 :

 3 Fig. 3: Total voltage drop considering every non-linearity of the VSI : (a) -for positive current polarity, (b) -for negative current polarity

Fig. 4 :Fig. 5 :

 45 Fig. 4: Control scheme setup for both simulation validation

Fig. 6 :Fig. 7 :

 67 Fig.6: Map in the i dq plane of the estimated voltage drops (V). Comparison between the experimental points, the simulation and the proposed direct analytical derivation

TABLE I :

 I Simulated and experimental setup parameters

	Parameter	Symbol	Value
	DC bus Voltage	VDC	350 V
	Switching frequency	fsw	10kHz
	Dead time	T dead	700 ns
	Turn-on delay	ton	120 ns
	Turn-off delay	t of f	100 ns
	MOSFET On resistance	R on DS	3.2 mΩ
	Diode forward voltage	V d 0	0.8
	Diode on resistance	r d	2.3 mΩ
	Intrinsic parasitic capacitance	Coss	25 nF
	Number of pole pairs	p	4
	Motor Inductances	L d	460µH
		Lq	430 µH
	Motor Resistance	Rs	39mΩ
	Flux linkage	ϕ f	0.131 Wb
	Speed reference	Ω *	1000 rpm
	Maximal load Torque	Γ ch	410 Nm