Sylvie Benzoni-Gavage

Jean-François Coulombel
email: jean-francois.coulombel@univ-nantes.fr

The amplitude equation for weakly nonlinear reversible phase boundaries

This technical note is a complement to an earlier paper [Benzoni-Gavage & Rosini

Introduction

The reader is assumed to be familiar with the weakly nonlinear theory developed in [START_REF] Benzoni-Gavage | Weakly nonlinear surface waves and subsonic phase boundaries[END_REF], which is very much inspired from the seminal work [START_REF] Hunter | Nonlinear surface waves[END_REF]. In particular, we keep the same notation as in [3, Section 2], and consider the liquid-vapor phase transition model described in [START_REF] Benzoni-Gavage | Weakly nonlinear surface waves and subsonic phase boundaries[END_REF]Paragraphs 3.1,3.2]. Our goal is to make explicit the amplitude equation, which corresponds to [START_REF] Benzoni-Gavage | Weakly nonlinear surface waves and subsonic phase boundaries[END_REF]Equation (2.20)] in an abstract framework. We adopt slightly different conventions compared with [START_REF] Benzoni-Gavage | Weakly nonlinear surface waves and subsonic phase boundaries[END_REF]Paragraph 3.3]. Namely, the eigenmodes with nonzero real part for the state ahead of the phase transition are1 :

β - 1 := a ℓ -i u ℓ η 0 c 2 ℓ -u 2 ℓ , β + 1 := -a ℓ -i u ℓ η 0 c 2 ℓ -u 2 ℓ = -β - 1 , a ℓ := -c ℓ (c 2 ℓ -u 2 ℓ ) |η| 2 -η 2 0 . (1)
The eigenmodes with nonzero real part for the state behind the phase transition are:

β - 2 := -a r + i u r η 0 c 2 r -u 2 r , β + 2 := a r + i u r η 0 c 2 r -u 2 r = -β - 2 , a r := c r (c 2 r -u 2 r ) |η| 2 -η 2 0 . (2)
The remaining, purely imaginary, eigenmodes are:

β + 3 = • • • = β + d+1 := i η 0 u ℓ , β - 3 = • • • = β - d+1 := - i η 0 u r .
The corresponding right eigenvectors are:

R ± 1 := r ± 1 0 , r - 1 :=   -i η 0 + u ℓ β - 1 i c 2 ℓ η -a ℓ   , r + 1 :=   i η 0 -u ℓ β + 1 -i c 2 ℓ η -a ℓ   = r - 1 , (3)
R ± 2 := 0 r ± 2 , r - 2 :=   -i η 0 -u r β - 2 i c 2 r η -a r   , r + 2 :=   i η 0 + u r β + 2 -i c 2 r η -a r   = r - 2 , (4)
R + j := r + j 0 , R - j := 0 r - j , j = 3, . . . , d + 1 , r + j :=   0 η 0 ěj-2 u ℓ η • ěj-2   , r - j :=   0 η 0 ěj-2 u r η • ěj-2   , (5)
with ě1 := η and the d -2 vectors ě2 , . . . , ěd-1 ∈ R d-1 span η⊥ . The left eigenvectors are:

L ± 1 := ℓ ± 1 0 , ℓ - 1 := c 2 ℓ -u 2 ℓ 2 a ℓ (u ℓ a ℓ + i c 2 ℓ η 0 )   i η 0 -2 u ℓ β + 1 -i η β + 1   , ℓ + 1 := c 2 ℓ -u 2 ℓ 2 a ℓ (u ℓ a ℓ -i c 2 ℓ η 0 )   -i η 0 + 2 u ℓ β - 1 i η -β - 1   = ℓ - 1 , L ± 2 := 0 ℓ ± 2 , ℓ - 2 := c 2 r -u 2 r 2 a r (u r a r + i c 2 r η 0 )   -i η 0 -2 u r β + 2 i η β + 2   , ℓ + 2 := c 2 r -u 2 r 2 a r (u r a r -i c 2 r η 0 )   i η 0 + 2 u r β - 2 -i η -β - 2   = ℓ - 2 , L + j := ℓ + j 0 , L - j := 0 ℓ - j , j = 3, . . . , d + 1 , ℓ + 3 := 1 η 2 0 + u 2 ℓ |η| 2   u ℓ -η 0 /(u ℓ |η| 2 ) η -1   , ℓ - 3 := 1 η 2 0 + u 2 r |η| 2   -u r η 0 /(u r |η| 2 ) η 1   , ℓ + j := - 1 u ℓ η 0   0 ě′ j-2 0   , ℓ - j := 1 u r η 0   0 ě′ j-2 0   , j = 4, . . . , d + 1 ,
After linearizing the jump conditions we are left with the matrices

H(v) :=     0 0 1 0 0 -1 0 u ℓ I d-1 0 0 -u r I d-1 0 c 2 ℓ -u 2 ℓ 0 2 u ℓ -(c 2 r -u 2 r ) 0 -2 u r (c 2 ℓ -u 2 ℓ ) u ℓ 0 u 2 ℓ + µ -(c 2 r -u 2 r ) u r 0 -u 2 r -µ     , (6)
where

µ := 1 2 u 2 ℓ + g(ρ ℓ ) = 1 2 u 2 r + g(ρ r ) , (7)
and

J(v) η :=     [ρ] η 0 [p] η 0 (µ [ρ] -[p]) η 0     . ( 8
)
Here we have used the relation 1 2

j [u] + [ρ f ] = µ [ρ] -[p] ,
where µ is defined in (7), in order to simplify the last entry of J(v) η.

The Lopatinskii determinant is defined by

∆(η) := det J(v) η H(v) R - 1 . . . H(v) R - d+1 . (9)
Each column vector in the above determinant is computed by using (8), ( 6) and ( 3), (4), (5) for the definition of the right eigenvectors. Then some simple manipulations on the rows and columns of the above determinant yield (see [START_REF] Benzoni-Gavage | Stability of multi-dimensional phase transitions in a van der Waals fluid[END_REF] for similar computations):

∆(η) = (-u r η 0 ) d-2 u r det ě1 . . . ěd-1 det     [ρ] η 0 a ℓ a r |η| 2 0 u ℓ a ℓ + i η 0 c 2 ℓ u r a r + i η 0 c 2 r 2 u r |η| 2 -[p] η 0 i η 0 u ℓ c 2 ℓ i η 0 u r c 2 r u 2 r |η| 2 [p] -i u ℓ c 2 ℓ -i u r c 2 r η 0     .
For future use, we introduce the quantity

Υ := (-u r η 0 ) d-2 u r det ě1 . . . ěd-1 ∈ R \ {0} . ( 10
)
The expression of the Lopatinskii determinant then reduces to

∆(η) = -[ρ] [u] Υ (η 2 0 + u 2 r |η| 2 ) u ℓ u r a ℓ a r + c 2 ℓ c 2 r η 2 0 . (11)
We fix a root η of the Lopatinskii determinant (see [START_REF] Benzoni-Gavage | Stability of multi-dimensional phase transitions in a van der Waals fluid[END_REF] for the properties of such roots, in particular the location of | η 0 |/| η | with respect to various velocities associated with the phase transition). From now on, an underline refers to evaluation at the root η of the Lopatinskii determinant.

We now define a vector σ ∈ C d+2 by computing some of the minors of ∆( η ). More precisely, the vectors H(v) R - 1 , . . . , H(v) R - d+1 are linearly independent and we can thus define a vector σ

∈ C d+2 \ {0} satisfying ∀ X ∈ C d+2 , det X H(v) R - 1 . . . H(v) R - d+1 = σ * X .
This vector σ can be computed explicitly by performing some elementary manipulations on each minor of (9). We do not give the detailed calculations but rather give the expression of σ. We find

σ * = Υ D 1 Ď ηT D d+1 D d+2 , (12)
where Υ denotes the quantity Υ in (10) evaluated at the frequency η, and2

D 1 + µ D d+2 := -( η 2 0 + u 2 r | η| 2 ) [u] c 2 ℓ c 2 r η 0 -i u ℓ u r (c 2 r a ℓ -c 2 ℓ a r ) , Ď := -[u] u r a ℓ (u r a r -i c 2 r η 0 ) + a r (u ℓ a ℓ -i c 2 ℓ η 0 ) , D d+1 := -i ( η 2 0 + u 2 r | η| 2 ) (u r c 2 r a ℓ -u ℓ c 2 ℓ a r ) , D d+2 := [u] η 0 (a ℓ a r + c 2 ℓ c 2 r | η| 2 ) + i η 2 0 (c 2 r a ℓ -c 2 ℓ a r ) -i | η| 2 u r (u ℓ -2 u r ) c 2 r a ℓ + u ℓ u r c 2 ℓ a r . (13)
It remains to compute the coefficients γ 1 , γ 2 satisfying:

J(v) η + γ 1 H(v) R - 1 + γ 2 H(v) R - 2 = 0 .
Observe that our convention differs from that in [3, page 1477]. We get:

γ 1 := [ρ] u r η 0 u r a ℓ -i c 2 ℓ η 0 = -i [ρ] c 2 r η 2 0 a ℓ (u ℓ a r -i c 2 r η 0 ) , γ 2 := -[ρ] u ℓ η 0 u ℓ a r -i c 2 r η 0 = i [ρ] c 2 ℓ η 2 0 a r (u r a ℓ -i c 2 ℓ η 0 ) , ( 14
)
where the equalities follow from the relation u ℓ u r a ℓ a r + c 2 ℓ c 2 r η 2 0 = 0 that is satisfied by the root η of the Lopatinskii determinant. For notational convenience, we also set

γ 3 = • • • = γ d+1 := 0 .
Following [3, Proposition 2.2], the evolution of a weakly nonlinear phase transition is governed by a scalar amplitude w obeying a nonlocal Burgers equation:

a 0 (k) ∂ τ w(τ, k) + R a 1 (k -k ′ , k ′ ) w(τ, k -k ′ ) w(τ, k ′ ) dk ′ = 0 , ( 15
)
where a 0 and a 1 are given by Equations (2.24) and (2.25) in [3, page 1471]. With the present notation, this yields

a 0 (k) = α 0 /(i k) if k > 0, α 0 /(i k) if k < 0,
and α 0 is a complex number whose definition is recalled in Equation ( 17) below. The expression of the kernel a 1 is recalled and made explicit in Section 3 below.

2 Computation of the coefficient α 0 Proposition 1. The coefficient α 0 in the expression of a 0 is given by

α 0 = - [ρ] [u] Υ η 0 ( η 2 0 + u 2 r | η| 2 ) u 2 ℓ u 2 r a 2 ℓ c 2 ℓ + a 2 r c 2 r + 2 c 2 ℓ c 2 r η 2 0 , (16)
and it coincides with the derivative of the Lopatinskii determinant ∆ with respect to η 0 at its root η . In particular, α 0 is a nonzero real number.

Proof. We recall that the expression of α 0 is

α 0 = σ * [ f0 (v)] + i σ * H(v) R + p (L + p ) * γ q R - q β + p -β - q , ( 17
)
where we use Einstein's summation convention over repeated indices. We first observe that the Hermitian product (L + p ) * R - q , q = 1, 2, vanishes as soon as p is larger than 4. In the same way, the products (L + 1 ) * R - 2 , (L + 3 ) * R - 2 and (L + 2 ) * R - 1 vanish so the expression of α 0 reduces to the sum of four terms:

α 0 = σ * [ f0 (v)]+i σ * H(v) R + 2 (L + 2 ) * γ 2 R - 2 β + 2 -β - 2 +i σ * H(v) R + 1 (L + 1 ) * γ 1 R - 1 β + 1 -β - 1 +i σ * H(v) R + 3 (L + 3 ) * γ 1 R - 1 β + 3 -β - 1 . (18)
We now compute each of these four quantities separately. Using

[ f0 (v)] = 1 η 0 J(v) η - 1 η 0     0 [p] η 0 0     ,
and the orthogonality relation σ * J(v) η = 0, we get

σ * [ f0 (v)] = -[ρ] Υ u ℓ u r | η| 2 η 0 Ď = [ρ] [u] Υ u ℓ u 2 r | η| 2 η 0 a ℓ (u r a r -i c 2 r η 0 ) + a r (u ℓ a ℓ -i c 2 ℓ η 0 ) . ( 19
)
We now turn to the second term on the right in (18). There holds

(L + 2 ) * γ 2 R - 2 β + 2 -β - 2 = - c 2 r | η| 2 (u r a r -i c 2 r η 0 ) 2 a 2 r ( η 2 0 + u 2 r | η| 2 ) γ 2 ,
and we also compute

i σ * H(v) R + 2 = i σ * H(v) R - 2 = 2 Υ c 2 r η 0 (D d+1 + u r D d+2 ) -u r | η| 2 Ď = 2 [u] Υ c 2 r a r ( η 2 0 + u 2 r | η| 2 ) (u r a ℓ -i c 2 ℓ η 0 ) , (20)
where the latter relation is obtained by using the expressions (13). Recalling the expression (14) of γ 2 , we obtain

i σ * H(v) R + 2 (L + 2 ) * γ 2 R - 2 β + 2 -β - 2 = [ρ] [u] Υ i u ℓ u r a ℓ a r c 2 r | η| 2 (u r a r -i c 2 r η 0 ) . ( 21
)
We now examine the third term on the right in (18). There holds:

(L + 1 ) * γ 1 R - 1 β + 1 -β - 1 = - c 2 ℓ | η| 2 (u ℓ a ℓ -i c 2 ℓ η 0 ) 2 a 2 ℓ ( η 2 0 + u 2 ℓ | η| 2 ) γ 1 ,
and we also compute

i σ * H(v) R + 1 = i σ * H(v) R - 1 = 2 Υ c 2 ℓ u ℓ | η| 2 Ď -η 0 (D d+1 + u ℓ D d+2 ) = -2 [u] Υ c 2 ℓ a ℓ ( η 2 0 + u 2 r | η| 2 ) (u ℓ a r -i c 2 r η 0 ) . ( 22
)
Using the expression (14) of γ 1 , we obtain

i σ * H(v) R + 1 (L + 1 ) * γ 1 R - 1 β + 1 -β - 1 = [ρ] [u] Υ i u ℓ u r a r a ℓ c 2 ℓ | η| 2 η 2 0 + u 2 r | η| 2 η 2 0 + u 2 ℓ | η| 2 (u ℓ a ℓ -i c 2 ℓ η 0 ) . ( 23
)
It remains to compute the last term on the right in (18) and to add the four expressions. First we compute (L + 3 ) * γ 1 R -

1 β + 3 -β - 1 = - c 2 ℓ η 2 0 + u 2 ℓ | η| 2 γ 1
, and we also compute

σ * H(v) R + 3 = -[u] Υ u ℓ | η| 2 (2 D d+1 + (u ℓ + u r ) D d+2 ) , (24)
where we have used the relation (which amounts to σ * H(v) R - 3 = 0):

D 1 + η 0 Ď + 2 u r D d+1 + (µ + u 2 r ) D d+2 = 0 .
The expression (24) can be factorized by using the definitions (13) of D d+1 , D d+2 , and we obtain

σ * H(v) R + 3 = - [u] 2 Υ u ℓ | η| 2 η 0 ( η 2 0 -u ℓ u r | η| 2 ) a r (u r a ℓ -i c 2 ℓ η 0 ) + a ℓ (u ℓ a r -i c 2 r η 0 ) . (25)
Using ( 14), we derive the expression

i σ * H(v) R + 3 (L + 3 ) * γ 1 R - 1 β + 3 -β - 1 = [ρ] [u] Υ i [u] u ℓ c 2 ℓ | η| 2 η 2 0 -u ℓ u r | η| 2 η 2 0 + u 2 ℓ | η| 2 (u r a r -i c 2 r η 0 ) . ( 26
)
According to the decomposition (18), the coefficient α 0 is the sum of the four quantities in ( 19), ( 21), ( 23) and (26). Factorizing [ρ] [u] Υ in each term, we first observe that the imaginary part of the sum equals zero. We can thus simplify α 0 by retaining only the real part of each term. This leads to the expression

η 0 α 0 [ρ] [u] Υ = -c 2 ℓ c 2 r η 2 0 (u 2 r + u ℓ u r ) | η| 2 + u ℓ u r a ℓ c 2 r η 2 0 a r c 2 r | η| 2 + u ℓ u r a r c 2 ℓ η 2 0 a ℓ c 2 ℓ | η| 2 η 2 0 + u 2 r | η| 2 η 2 0 + u 2 ℓ | η| 2 + (u ℓ u r -u 2 ℓ ) | η| 2 c 2 ℓ c 2 r η 2 0 η 2 0 -u ℓ u r | η| 2 η 2 0 + u 2 ℓ | η| 2 .
At this stage, some elementary manipulations lead to the expression (16) of α 0 .

The link between α 0 and the partial derivative ∂ η 0 ∆( η) comes from the relation

η 0 ∂ ∂η 0 (u ℓ u r a ℓ a r + c 2 ℓ c 2 r η 2 0 ) η = u 2 ℓ u 2 r a 2 ℓ c 2 ℓ + a 2 r c 2 r + 2 c 2 ℓ c 2 r η 2 0 ,
which is obtained by differentiating (11) and the expressions (1), ( 2) with respect to η 0 and then evaluating at η.

Since the coefficient α 0 is real, we obtain a 0 (k) = α 0 /(i k) for all k = 0. In particular, the amplitude equation ( 15) reduces to

∂ τ w(τ, k) + i k α 0 R a 1 (k -k ′ , k ′ ) w(τ, k -k ′ ) w(τ, k ′ ) dk ′ = 0 .
As far as smooth solutions are concerned, the Cauchy problem associated with this kind of nonlocal Burgers equation is known to be locally well-posed under rather simple algebraic conditions (see [START_REF] Benzoni-Gavage | Local well-posedness of nonlocal Burgers equations[END_REF]). These conditions are invariant under multiplication by a nonzero real constant, so that it is sufficient to investigate whether they are satisfied by the slightly simpler kernel 4 π a 1 . This is the purpose of the next section.

Computation of the quadratic kernel

We define the kernel q(k, k ′ ) := 4 π a 1 (k, k ′ ). Following [START_REF] Benzoni-Gavage | Weakly nonlinear surface waves and subsonic phase boundaries[END_REF], we can decompose q as follows:

q(k, k ′ ) = 5 j=1 q j (k, k ′ ) ,
where the kernels q 1 , . . . , q 5 are given by3

q 1 (k, k ′ ) := σ(k + k ′ ) * d-1 j=0 η j d f j (v r ) • ( r + (k, 0) + r + (k ′ , 0)) -d f j (v ℓ ) • ( r -(k, 0) + r -(k ′ , 0)) , (27)
q 2 (k, k ′ ) := -σ(k + k ′ ) * d 2 f d (v r ) • ( r + (k, 0), r + (k ′ , 0)) -d 2 f d (v ℓ ) • ( r -(k, 0), r -(k ′ , 0)) , (28)
q 3 (k, k ′ ) := i (k + k ′ ) +∞ 0 L(k + k ′ , z) dA(v, η) • r(k, z) • r(k ′ , z) dz , (29)
q 4 (k, k ′ ) := +∞ 0 L(k + k ′ , z) ∂ ∂z d Ȃd (v) • r(k, z) • r(k ′ , z) dz , (30)
q 5 (k, k ′ ) := - +∞ 0 L(k + k ′ , z) Ȃ(v, η) ∂ r ∂z (k, z) + ∂ r ∂z (k ′ , z) dz . ( 31
)
We first examine the kernel q 1 and derive its expression for all values of (k, k ′ ).

Lemma 1. Let us define the quantity

Q := 2 [ρ] [u] Υ ( η 2 0 + u 2 r | η| 2 ) ( β - 1 + β - 2 ) i u ℓ u r a ℓ a r u ℓ a r + i c 2 r η 0 u ℓ a r -i c 2 r η 0 . ( 32
)
Then the kernel q 1 in (27) satisfies

q 1 (k, k ′ ) = 0 if k > 0 and k ′ > 0, Q if k > 0, k ′ < 0 and k + k ′ > 0.
Proof. The function g 0 is known to be an entropy for the isothermal Euler equations with corresponding flux (g 1 , . . . , g d ). We thus have the relations

∀ j = 0, . . . , d , dg 0 (v) A j (v) = dg j (v) ,
where we use the convention A 0 (v) = I for all v. Using this relation between the Jacobian matrices, we get

d-1 j=0 η j d f j (v r ) • r - 2 = d-1 j=0 η j A j (v r ) r - 2 dg 0 (v r ) d-1 j=0 η j A j (v r ) r - 2 = i β - 2 A d (v r ) r - 2 dg d (v r ) • r - 2 = -i β - 2 H(v) R - 2 .
Similarly we have

d-1 j=0 η j d f j (v ℓ ) • r - 1 = -i β - 1 H(v) R - 1 .
For k > 0 and k ′ > 0, we thus get

q 1 (k, k ′ ) = 2 σ * d-1 j=0 γ 2 η j d f j (v r ) • r - 2 -2 σ * d-1 j=0 γ 1 η j d f j (v ℓ ) • r - 1 = -2 i γ 2 β - 2 σ * H(v) R - 2 + 2 i γ 1 β - 1 σ * H(v) R - 1 = 0 , because σ is orthogonal to both H(v) R - 1 and H(v) R - 2 .
Let us now consider the case k > 0, k ′ < 0 and k + k ′ > 0. Using the same relations as above for the differentials dg j , we obtain

q 1 (k, k ′ ) = σ * d-1 j=0 η j d f j (v r ) • (γ 2 r - 2 + γ 2 r + 2 ) -σ * d-1 j=0 η j d f j (v ℓ ) • (γ 1 r - 1 + γ 1 r + 1 ) = -i γ 2 β + 2 σ * H(v) R + 2 + i γ 1 β + 1 σ * H(v) R + 1 .
The Hermitian products σ * H(v) R + 2 and σ * H(v) R + 1 have already been computed in the proof of Proposition 1, see (20) and ( 22). We then obtain

q 1 (k, k ′ ) = -2 [u] Υ ( η 2 0 + u 2 r | η| 2 ) c 2 r a r (u r a ℓ -i c 2 ℓ η 0 ) γ 2 β + 2 -2 [u] Υ ( η 2 0 + u 2 r | η| 2 ) c 2 ℓ a ℓ (u ℓ a r -i c 2 r η 0 ) γ 1 β + 1 .
We use the definition (14) to obtain

c 2 r a r (u r a ℓ -i c 2 ℓ η 0 ) γ 2 = c 2 ℓ a ℓ (u ℓ a r -i c 2 r η 0 ) γ 1 = i [ρ] c 2 ℓ c 2 r η 2 0 u ℓ a r -i c 2 r η 0 u ℓ a r + i c 2 r η 0 ,
and the claim follows using the relation

c 2 ℓ c 2 r η 2 0 = -u ℓ u r a ℓ a r .
Deriving the expression of the kernels q 3 , q 4 , q 5 requires the expression of the row vector L(k + k ′ , z), which we derive right now. Lemma 2. For k > 0, there holds

L(k, z) = ω 1 γ 1 exp(-k β + 1 z) ℓ 1 + ω 3 γ 1 exp(-k β + 3 z) ℓ 3 ω 2 γ 2 exp(-k β + 2 z) ℓ 2 ,
where we have set

ℓ 1 := i η 0 -2 u ℓ β + 1 -i ηT β + 1 , ℓ 3 := -u 2 ℓ | η| 2 η 0 ηT u ℓ | η| 2 , ℓ 2 := -i η 0 -2 u r β + 2 i ηT β + 2 ,
and

ω 1 := [ρ] [u] Υ η 2 0 + u 2 r | η| 2 η 2 0 + u 2 ℓ | η| 2 i u ℓ η 0 (u r a r -i c 2 r η 0 ), ω 3 := [ρ] [u] 2 Υ η 2 0 -u ℓ u r | η| 2 η 2 0 + u 2 ℓ | η| 2 (u r a r -i c 2 r η 0 ), ω 2 := [ρ] [u] Υ i u r η 0 (u ℓ a ℓ -i c 2 ℓ η 0 ).
Proof. We first observe that σ is orthogonal to the vectors H(v) R + p for p ≥ 4, so for k > 0 the expression of L(k, z) reduces to

L(k, z) = 3 p=1 σ * H(v) R + p exp(-k β + p z) (L + p ) * .
The expression of the products σ * H(v) R + p , p = 1, 2, 3 can be found in (20), ( 22), (25), and we use the definitions of the left eigenvectors L + p to derive the expressions given in Lemma 2.

We now examine the kernel q 5 , which, as q 1 but unlike q 2 , q 3 , q 4 , does not contain any term in p ′′ (ρ ℓ,r ).

Lemma 3. With Q defined in (32), the kernel q 5 in (31) satisfies

q 5 (k, k ′ ) =    0 if k > 0 and k ′ > 0, Q k ′ k if k > 0, k ′ < 0 and k + k ′ > 0.
Proof. For k > 0 and k ′ > 0, we compute

Ȃ(v, η) ∂ r ∂z (k, z) + ∂ r ∂z (k ′ , z) = i γ 1 ( β - 1 ) 2 (k exp(k β - 1 z) + k ′ exp(k ′ β - 1 z)) A d (v ℓ ) r - 1 i γ 2 ( β - 2 ) 2 (k exp(k β - 2 z) + k ′ exp(k ′ β - 2 z)) A d (v r ) r - 2 .
Using the orthogonality properties

ℓ 1 A d (v ℓ ) r - 1 = ℓ 3 A d (v ℓ ) r - 1 = ℓ 3 A d (v r ) r - 2 = 0, we get q 5 (k, k ′ ) = 0 if k > 0 and k ′ > 0 because the integrand in (31) vanishes.
Let us now consider the case k > 0, k ′ < 0 and k + k ′ > 0. From the previous argument, we still find that the term L(k + k ′ , z) Ȃ(v, η) ∂ z r(k, z) vanishes. We thus get

q 5 (k, k ′ ) = - +∞ 0 L(k + k ′ , z) Ȃ(v, η) ∂ r ∂z (k ′ , z) dz = -i k ′ +∞ 0 L(k + k ′ , z) γ 1 ( β + 1 ) 2 exp(k ′ β + 1 z) A d (v ℓ ) r + 1 γ 2 ( β + 2 ) 2 exp(k ′ β + 2 z) A d (v r ) r + 2 dz .
We now use the expression of L(k + k ′ , z) in Lemma 2. The expression of q 5 (k, k ′ ) is simplified by recalling the orthogonality property ℓ 3 A d (v ℓ ) r + 1 = 0 and we get

q 5 (k, k ′ ) = -i k ′ +∞ 0 γ 1 γ 1 ω 1 ℓ 1 A d (v ℓ ) r + 1 ( β + 1 ) 2 exp(-k β + 1 z) dz -i k ′ +∞ 0 γ 2 γ 2 ω 2 ℓ 2 A d (v r ) r + 2 ( β + 2 ) 2 exp(-k β + 2 z) dz , = -i k ′ k γ 1 γ 1 ω 1 ℓ 1 A d (v ℓ ) r + 1 β + 1 + γ 2 γ 2 ω 2 ℓ 2 A d (v r ) r + 2 β + 2 .
The conclusion of Lemma 3 then follows from the relations

γ 1 γ 1 = - γ 2 γ 2 = - u ℓ a r -i c 2 r η 0 u ℓ a r + i c 2 r η 0 , ω 1 ℓ 1 A d (v ℓ ) r + 1 = -ω 2 ℓ 2 A d (v r ) r + 2 = 2 [ρ] [u] Υ ( η 2 0 + u 2 r | η| 2 ) u ℓ u r a ℓ a r .
We immediately get Corollary 1. With Q defined in (32), the kernels q 1 , q 5 in (27), (31) satisfy

(q 1 + q 5 )(k, k ′ ) =    0 if k > 0 and k ′ > 0, Q 1 + k ′ k if k > 0, k ′ < 0 and k + k ′ > 0.
Our goal now is to derive an explicit expression for the kernels q 2 , q 3 , q 4 in (28), (29), (30). This is more intricate because these kernels are quadratic with respect to the vector r(k, z) and there is more algebra involved to obtain a factorized expression in each region of the (k, k ′ )-plane. We begin with some preliminary computations that will be useful in Propositions 2 and 3 below. Lemma 4. The coordinates (13) of the vector σ satisfy

γ 1 Ď = -[ρ] [u] u r η 0 (u r a r -i c 2 r η 0 ) , γ 2 Ď = [ρ] [u] u r η 0 (u ℓ a ℓ -i c 2 ℓ η 0 ) , γ 1 D d+1 = -[ρ] u r ( η 2 0 + u 2 r | η| 2 ) (u ℓ a r + i c 2 r η 0 ) , γ 2 D d+1 = -[ρ] u ℓ ( η 2 0 + u 2 r | η| 2 ) (u r a ℓ + i c 2 ℓ η 0 ) , γ 1 D d+2 = [ρ] ( η 2 0 + u 2 r | η| 2 ) (u r a r + i c 2 r η 0 ) -[ρ] [u] u r | η| 2 (u r a r -i c 2 r η 0 ) , γ 2 D d+2 = [ρ] ( η 2 0 + u 2 r | η| 2 ) (u ℓ a ℓ + i c 2 ℓ η 0 ) + [ρ] [u] u r | η| 2 (u ℓ a ℓ -i c 2 ℓ η 0 ) .
Proof. From the definition (13), there holds

Ď = [u] u r a r (u r a ℓ -i c 2 ℓ η 0 ) + a ℓ (u ℓ a r -i c 2 r η 0 ) ,
and we then use the expression (14) of γ 1 , γ 2 to compute γ 1 Ď and γ 2 Ď. The expressions of γ 1,2 D d+1 follow from the observation that D d+1 satisfies

η 0 D d+1 = -( η 2 0 + u 2 r | η| 2 ) (u r a ℓ -i c 2 ℓ η 0 ) (u ℓ a r + i c 2 r η 0 ) = ( η 2 0 + u 2 r | η| 2 ) (u r a ℓ + i c 2 ℓ η 0 ) (u ℓ a r -i c 2 r η 0 ) .
We then use again (14).

To compute the product γ 1 D d+2 , we recall the relation (22) which we found in the proof of Proposition 1. It reads

η 0 (D d+1 + u ℓ D d+2 ) -u ℓ | η| 2 Ď = [u] a ℓ ( η 2 0 + u 2 r | η| 2 ) (u ℓ a r -i c 2 r η 0 )
. We multiply the latter relation by γ 1 and use the previous expressions of γ 1 Ď and γ 1 D d+1 to obtain that of γ 1 D d+2 . The expression of γ 2 D d+2 is then easily deduced by using γ 2 /γ 1 = i c 2 ℓ η 2 0 /(u r a r ).

We now compute the kernels q 2 , q 3 , q 4 in the region {k > 0, k ′ > 0}.

Proposition 2. Let us define the quantities

Q ℓ := [ρ] [u] Υ u ℓ u r a r a ℓ ( η 2 0 + u 2 r | η| 2 ) ( η 2 0 + u 2 ℓ | η| 2 ) γ 1 (i η 0 -u ℓ β - 1 ) , Q r := [ρ] [u] Υ u ℓ u r a ℓ a r ( η 2 0 + u 2 r | η| 2 ) 2 γ 2 (i η 0 + u r β - 2 ) , (33)
Q ♯ := 2 [ρ] Υ ( η 2 0 + u 2 r | η| 2 ) ( η 2 0 + u ℓ u r | η| 2 ) i c 2 ℓ c 2 r η 0 c 2 r γ 2 ρ r u r - c 2 ℓ γ 1 ρ ℓ u ℓ .
Then the kernels q 2 , q 3 , q 4 defined in (28), ( 29) and (30) satisfy

(q 2 + q 3 + q 4 )(k, k ′ ) = p ′′ (ρ ℓ ) 2 + c 2 ℓ ρ ℓ Q ℓ + p ′′ (ρ r ) 2 + c 2 r ρ r Q r + Q ♯ ,
for all k > 0 and k ′ > 0.

Proof. We first recall the expressions of the second differentials that are involved in the kernels q 2 , q 3 , q 4 , see [3, page 1480]:

d-1 k=1 ηk d 2 f k (v ℓ,r ) • (v, v) = p ′′ (ρ ℓ,r )   0 ρ 2 η 0   + 2 ρ ℓ,r   0 η • ǰ ǰ η • ǰ (j d -u ℓ,r ρ)   , (34)
d 2 f d (v ℓ,r ) • (v, v) = d 2 f d (v ℓ,r ) • (v, v) d 2 g d (v ℓ,r ) • (v, v) = p ′′ (ρ ℓ,r )     0 0 ρ 2 u ℓ,r ρ 2     + 2 ρ ℓ,r     0 (j d -u ℓ,r ρ) ǰ (j d -u ℓ,r ρ) 2 0     + 1 ρ ℓ,r     0 0 0 3 u ℓ,r (j d -u ℓ,r ρ) 2 -u ℓ,r c 2 ℓ,r ρ 2 + 2 c 2 ℓ,r ρ j d + u ℓ,r ǰ • ǰ    . ( 35
)
The proof of Proposition 2 splits in several steps. We assume from now on that k and k ′ are both positive and we wish to compute the expression of the kernel q 2 + q 3 + q 4 .

• Step 1: computation of the p ′′ (ρ ℓ ) factor. In this first step, we collect all the terms that involve p ′′ (ρ ℓ ) in q 2 + q 3 + q 4 . The contribution of the kernel q 2 equals Υ

(D d+1 + u ℓ D d+2 ) γ 2 1 (i η 0 -u ℓ β - 1 ) 2 = -[ρ] [u] Υ γ 1 (i η 0 -u ℓ β - 1 ) 2 ( η 2 0 + u 2 r | η| 2 ) i c 2 r η 0 + (u r a r -i c 2 r η 0 ) u ℓ u r | η| 2 , ( 36
)
where we have used Lemma 4 to compute the product (D d+1 + u ℓ D d+2 ) γ 1 . The contribution of the kernel q 3 equals i

(k + k ′ ) +∞ 0 ω 1 γ 1 (-i ηT ) e -(k+k ′ ) β + 1 z γ 2 1 (i η 0 -u ℓ β - 1 ) 2 η e (k+k ′ ) β - 1 z dz + i (k + k ′ ) +∞ 0 ω 3 γ 1 ( η 0 ηT ) e -(k+k ′ ) β + 3 z γ 2 1 (i η 0 -u ℓ β - 1 ) 2 η e (k+k ′ ) β - 1 z dz = γ 1 (i η 0 -u ℓ β - 1 ) 2 ω 1 | η| 2 β + 1 -β - 1 + ω 3 | η| 2 i η 0 β + 3 -β - 1 .
Similarly, the contribution of the kernel q 4 reads

- +∞ 0 ω 1 γ 1 β + 1 e -(k+k ′ ) β + 1 z γ 2 1 (i η 0 -u ℓ β - 1 ) 2 (k + k ′ ) β - 1 e (k+k ′ ) β - 1 z dz - +∞ 0 ω 3 γ 1 u ℓ | η| 2 e -(k+k ′ ) β + 3 z γ 2 1 (i η 0 -u ℓ β - 1 ) 2 (k + k ′ ) β - 1 e (k+k ′ ) β - 1 z dz = -γ 1 (i η 0 -u ℓ β - 1 ) 2 ω 1 β + 1 β - 1 β + 1 -β - 1 + ω 3 | η| 2 u ℓ β - 1 β + 3 -β - 1 .
Adding the contributions of q 3 and q 4 gives the term

γ 1 (i η 0 -u ℓ β - 1 ) 2 - 1 2 a ℓ a 2 ℓ c 2 ℓ + c 2 ℓ | η| 2 ω 1 + u ℓ | η| 2 ω 3 . (37)
We now use the definitions of ω 1 and ω 3 , see Lemma 2, and add the contributions in (36) and (37) in order to obtain the p ′′ (ρ ℓ ) term in q 2 + q 3 + q 4 . We first obtain that the sum of the right hand side of (36) and the expression in (37) equals

-[ρ] [u] Υ γ 1 (i η 0 -u ℓ β - 1 ) 2 η 2 0 + u 2 r | η| 2 2 c 2 ℓ ( η 2 0 + u 2 ℓ | η| 2 ) (u ℓ a ℓ + i c 2 ℓ η 0 ) c 2 r η 2 0 + u r a r u ℓ a ℓ u 2 ℓ c 2 ℓ | η| 2 .
This last expression is simplified a little further by observing that we have

(i η 0 -u ℓ β - 1 ) (u ℓ a ℓ + i c 2 ℓ η 0 ) = -1 c 2 ℓ -u 2 ℓ (u 2 ℓ a 2 ℓ + c 4 ℓ η 2 0 ) = -c 2 ℓ ( η 2 0 + u 2 ℓ | η| 2 ) , and
c 2 r η 2 0 + u r a r u ℓ a ℓ u 2 ℓ c 2 ℓ | η| 2 = u ℓ u r a r a ℓ c 2 ℓ | η| 2 - a 2 ℓ c 2 ℓ = u ℓ u r a r a ℓ ( η 2 0 + u 2 ℓ | η| 2 ) .
Eventually, we find that the sum of all the terms that involve p ′′ (ρ ℓ ) factorizes as Q ℓ /2 where Q ℓ is defined in (33).

• Step 2: computation of the p ′′ (ρ r ) factor. We follow the same strategy as in the first step and compute the contribution that involves p ′′ (ρ r ) in each kernel. The contribution of the kernel q 2 equals

-Υ (D d+1 +u r D d+2 ) γ 2 2 (i η 0 +u r β - 2 ) 2 = [ρ] [u] Υ γ 2 (i η 0 +u r β - 2 ) 2 u ℓ a ℓ c 2 r i η 0 u r a r -u 2 r c 2 r | η| 2 , ( 38
)
where we have used Lemma 4 to simplify (D d+1 + u r D d+2 ) γ 2 . The contribution of the kernel q 3 equals i

(k + k ′ ) +∞ 0 ω 2 γ 2 (i ηT ) e -(k+k ′ ) β + 2 z γ 2 2 (i η 0 + u r β - 2 ) 2 η e (k+k ′ ) β - 2 z dz = -γ 2 (i η 0 + u r β - 2 ) 2 ω 2 | η| 2 β + 2 -β - 2 ,
and the contribution of the kernel q 4 reads +∞ 0

ω 2 γ 2 β + 2 e -(k+k ′ ) β + 2 z γ 2 2 (i η 0 + u r β - 2 ) 2 (k + k ′ ) β - 2 e (k+k ′ ) β - 2 z dz = γ 2 (i η 0 + u r β - 2 ) 2 ω 2 β + 2 β - 2 β + 2 -β - 2 .
Adding the contributions of q 3 and q 4 gives the term

-γ 2 (i η 0 + u r β - 2 ) 2 ω 2 2 a r a 2 r c 2 r + c 2 r | η| 2 .
When we add the latter term with the expression in (38), we obtain

-[ρ] [u] Υ γ 2 (i η 0 + u r β - 2 ) 2 ( η 2 0 + u 2 r | η| 2 ) i u r η 0 2 a r (u ℓ a ℓ + i c 2 ℓ η 0 ) ,
and this quantity is further simplified by using the relation

(i η 0 + u r β - 2 ) i η 0 (u ℓ a ℓ + i c 2 ℓ η 0 ) = -u ℓ a ℓ u 2 r a 2 r + c 4 r η 2 0 c 2 r (c 2 r -u 2 r ) = -u ℓ a ℓ ( η 2 0 + u 2 r | η| 2 ) .
Eventually, we find that the sum of all the terms that involve p ′′ (ρ r ) factorizes as Q r /2.

• Step 3: computation of the remaining terms. In order to prove Proposition 2, we can assume from now on, and without loss of generality that p ′′ (ρ ℓ ) = p ′′ (ρ r ) = 0 in (34) and ( 35). With this simplification, we compute

d-1 k=1 ηk d 2 f k (v ℓ ) • (r - 1 , r - 1 ) = - 2 i c 4 ℓ | η| 2 ρ ℓ   0 -i η β - 1   , (39)
d 2 f d (v ℓ ) • (r - 1 , r - 1 ) = 2 c 2 ℓ ρ ℓ     0 -i c 2 ℓ β - 1 η c 2 ℓ | η| 2 + (i η 0 -u ℓ β - 1 ) 2 i c 2 ℓ η 0 β - 1 + u ℓ (i η 0 -u ℓ β - 1 ) 2     , ( 40
) d-1 k=1 ηk d 2 f k (v r ) • (r - 2 , r - 2 ) = 2 i c 4 r | η| 2 ρ r   0 i η β - 2   , (41)
d 2 f d (v r ) • (r - 2 , r - 2 ) = 2 c 2 r ρ r     0 i c 2 r β - 2 η c 2 r | η| 2 + (i η 0 + u r β - 2 ) 2 -i c 2 r η 0 β - 2 + u r (i η 0 + u r β - 2 ) 2     , (42)
where in ( 40) and ( 42), we have used the relations

c 2 ℓ ( β - 1 ) 2 = c 2 ℓ | η| 2 + (i η 0 -u ℓ β - 1 ) 2 , c 2 r ( β - 2 ) 2 = c 2 r | η| 2 + (i η 0 + u r β - 2 ) 2 .
With these expressions, let us look first at the kernel q 2 in (28). Using (40), we compute

σ * d 2 f d (v ℓ ) • (r - 1 , r - 1 ) γ 2 1 = 2 c 2 ℓ ρ ℓ Υ γ 1 i c 2 ℓ β - 1 ( η 0 γ 1 D d+2 -| η| 2 γ 1 Ď) + c 2 ℓ | η| 2 γ 1 D d+1 + (i η 0 -u ℓ β - 1 ) 2 γ 1 (D d+1 + u ℓ D d+2 ) ,
and Lemma 4 turns this expression into

σ * d 2 f d (v ℓ ) • (r - 1 , r - 1 ) γ 2 1 = 2 c 2 ℓ ρ ℓ [ρ] Υ ( η 2 0 + u 2 r | η| 2 ) γ 1 u r a r (u ℓ a ℓ + i c 2 ℓ η 0 ) β - 1 -u r c 2 ℓ | η| 2 (u ℓ a r + i c 2 r η 0 ) - 2 c 2 ℓ ρ ℓ [ρ] [u] Υ γ 1 (i η 0 -u ℓ β - 1 ) 2 ( η 2 0 + u 2 r | η| 2 ) i c 2 r η 0 + u ℓ u r | η| 2 (u r a r -i c 2 r η 0 ) . (43)
Similarly, we derive the relation

σ * d 2 f d (v r ) • (r - 2 , r - 2 ) γ 2 2 = - 2 c 2 r ρ r [ρ] Υ ( η 2 0 + u 2 r | η| 2 ) γ 2 u ℓ a ℓ (u r a r + i c 2 r η 0 ) β - 2 + u ℓ c 2 r | η| 2 (u r a ℓ + i c 2 ℓ η 0 ) + 2 c 2 r ρ r [ρ] [u] Υ γ 2 (i η 0 + u r β - 2 ) 2 ( η 2 0 + u 2 r | η| 2 ) i c 2 ℓ η 0 + u 2 r | η| 2 (u ℓ a ℓ -i c 2 ℓ η 0 ) . ( 44
)
The kernels q 3 , q 4 both read as a sum of two contributions, one from the system ahead of the phase boundary, and one from the system behind the phase boundary. We compute each of these contributions separately in order to combine them with either (43) or (44). Using (39), the 'left' contribution of q 3 equals i

(k + k ′ ) +∞ 0 ω 1 γ 1 ℓ 1 e -(k+k ′ ) β + 1 z γ 2 1 -2 i c 4 ℓ | η| 2 ρ ℓ   0 -i η β - 1   e (k+k ′ ) β - 1 z dz + i (k + k ′ ) +∞ 0 ω 3 γ 1 ℓ 3 e -(k+k ′ ) β + 3 z γ 2 1 -2 i c 4 ℓ | η| 2 ρ ℓ   0 -i η β - 1   e (k+k ′ ) β - 1 z dz = 2 c 4 ℓ ρ ℓ γ 1 ω 1 | η| 2 β + 1 -β - 1 ℓ 1   0 -i η β - 1   + 2 c 4 ℓ ρ ℓ γ 1 ω 3 | η| 2 β + 3 -β - 1 ℓ 3   0 -i η β - 1   , (45)
and, similarly (using now (41) rather than (39)), the 'right' contribution of q 3 equals i

(k + k ′ ) +∞ 0 ω 2 γ 2 ℓ 2 e -(k+k ′ ) β + 2 z γ 2 2 2 i c 4 r | η| 2 ρ r   0 i η β - 2   e (k+k ′ ) β - 2 z dz = - 2 c 4 r ρ r γ 2 ω 2 | η| 2 β + 2 -β - 2 ℓ 2   0 i η β - 2   . (46)
The kernel q 4 is computed by first observing that the vectors d

2 f d (v ℓ )•(r - 1 , r - 1 ), and d 2 f d (v r )•(r - 2 , r -
2 ) are obtained by retaining only the three first coordinates in (40) and (42):

d 2 f d (v ℓ ) • (r - 1 , r - 1 ) = 2 c 4 ℓ β - 1 ρ ℓ   0 -i η β - 1   , d 2 f d (v r ) • (r - 2 , r - 2 ) = 2 c 4 r β - 2 ρ r   0 i η β - 2   .
Consequently, the 'left' contribution of q 4 equals

- +∞ 0 ω 1 γ 1 ℓ 1 e -(k+k ′ ) β + 1 z γ 2 1 2 c 4 ℓ β - 1 ρ ℓ   0 -i η β - 1   (k + k ′ ) β - 1 e (k+k ′ ) β - 1 z dz - +∞ 0 ω 3 γ 1 ℓ 1 e -(k+k ′ ) β + 3 z γ 2 1 2 c 4 ℓ β - 1 ρ ℓ   0 -i η β - 1   (k + k ′ ) β - 1 e (k+k ′ ) β - 1 z dz = - 2 c 4 ℓ ρ ℓ γ 1 ω 1 ( β - 1 ) 2 β + 1 -β - 1 ℓ 1   0 -i η β - 1   - 2 c 4 ℓ ρ ℓ γ 1 ω 3 ( β - 1 ) 2 β + 3 -β - 1 ℓ 3   0 -i η β - 1   , (47)
and the 'right' contribution of q 4 equals

2 c 4 r ρ r γ 2 ω 2 ( β - 2 ) 2 β + 2 -β - 2 ℓ 2   0 i η β - 2   . (48)
The 'left' contribution of q 3 + q 4 is obtained by adding the expressions in ( 45) and ( 47), which gives

- 2 c 4 ℓ ρ ℓ γ 1 ω 1 ( β - 1 ) 2 -| η| 2 β + 1 -β - 1 ℓ 1   0 -i η β - 1   - 2 c 4 ℓ ρ ℓ γ 1 ω 3 ( β - 1 ) 2 -| η| 2 β + 3 -β - 1 ℓ 3   0 -i η β - 1   = - 2 c 2 ℓ ρ ℓ γ 1 (i η 0 -u ℓ β - 1 ) 2    ω 1 ℓ 1   0 -i η β - 1   + ω 3 ℓ 3   0 -i η β - 1      = c 2 ℓ ρ ℓ a ℓ γ 1 (i η 0 -u ℓ β - 1 ) 2 ( η 2 0 + u 2 ℓ | η| 2 -2 c 2 ℓ | η| 2 ) ω 1 + 2 u ℓ a ℓ | η| 2 ω 3 . (49)
The 'right' contribution of q 3 + q 4 is obtained by adding the expressions in ( 46) and (48), which gives

2 c 4 r ρ r γ 2 ω 2 ( β - 2 ) 2 -| η| 2 β + 2 -β - 2 ℓ 2   0 i η β - 2   = 2 c 2 r ρ r γ 2 ω 2 (i η 0 + u r β - 2 ) 2 β + 2 β - 2 -| η| 2 β + 2 -β - 2 = c 2 r ρ r a r γ 2 (i η 0 + u r β - 2 ) 2 ( η 2 0 + u 2 r | η| 2 -2 c 2 r | η| 2 ) ω 2 . ( 50
)
We can now compute the 'left' contribution of the full kernel q 2 + q 3 + q 4 by combining the expression in (43) with the one in (49). We use here the expressions of ω 1 and ω 3 given in Lemma 2. The sum of (43) and (49) reads

A ℓ := 2 c 2 ℓ ρ ℓ [ρ] Υ ( η 2 0 + u 2 r | η| 2 ) γ 1 u r a r (u ℓ a ℓ + i c 2 ℓ η 0 ) β - 1 -u r c 2 ℓ | η| 2 (u ℓ a r + i c 2 r η 0 ) - 2 c 2 ℓ ρ ℓ [ρ] [u] Υ γ 1 (i η 0 -u ℓ β - 1 ) 2 ( η 2 0 + u 2 r | η| 2 ) i c 2 r η 0 + u ℓ u r | η| 2 (u r a r -i c 2 r η 0 ) + c 2 ℓ ρ ℓ a ℓ [ρ] [u] Υ γ 1 (i η 0 -u ℓ β - 1 ) 2 2 u ℓ a ℓ | η| 2 [u] η 2 0 -u ℓ u r | η| 2 η 2 0 + u 2 ℓ | η| 2 (u r a r -i c 2 r η 0 ) +( η 2 0 + u 2 ℓ | η| 2 -2 c 2 ℓ | η| 2 ) i u ℓ η 0 (u r a r -i c 2 r η 0 ) η 2 0 + u 2 r | η| 2 η 2 0 + u 2 ℓ | η| 2 .
The three last rows in the definition of A ℓ are factorized after a little bit of algebra, and we get

A ℓ = 2 c 2 ℓ ρ ℓ [ρ] Υ ( η 2 0 + u 2 r | η| 2 ) γ 1 u r a r (u ℓ a ℓ + i c 2 ℓ η 0 ) β - 1 -u r c 2 ℓ | η| 2 (u ℓ a r + i c 2 r η 0 ) - 1 ρ ℓ [ρ] [u] Υ γ 1 (i η 0 -u ℓ β - 1 ) 2 u ℓ u r a r a ℓ (u ℓ a ℓ + i c 2 ℓ η 0 ) ( η 2 0 + u 2 r | η| 2 ) .
We have thus shown that the 'left' contribution A ℓ of q 2 + q 3 + q 4 reads

A ℓ = 2 c 2 ℓ ρ ℓ [ρ] Υ ( η 2 0 + u 2 r | η| 2 ) γ 1 u r a r (u ℓ a ℓ + i c 2 ℓ η 0 ) β - 1 -u r c 2 ℓ | η| 2 (u ℓ a r + i c 2 r η 0 ) + c 2 ℓ ρ ℓ Q ℓ . (51)
The 'right' contribution of q 2 + q 3 + q 4 by combining the expression in (44) (with a minus sign, recall the definition (28) of the kernel q 2 ) with the one in (50):

A r := 2 c 2 r ρ r [ρ] Υ ( η 2 0 + u 2 r | η| 2 ) γ 2 u ℓ a ℓ (u r a r + i c 2 r η 0 ) β - 2 + u ℓ c 2 r | η| 2 (u r a ℓ + i c 2 ℓ η 0 ) - 2 c 2 r ρ r [ρ] [u] Υ γ 2 (i η 0 + u r β - 2 ) 2 ( η 2 0 + u 2 r | η| 2 ) i c 2 ℓ η 0 + u 2 r | η| 2 (u ℓ a ℓ -i c 2 ℓ η 0 ) + c 2 r ρ r a r γ 2 (i η 0 + u r β - 2 ) 2 ( η 2 0 + u 2 r | η| 2 -2 c 2 r | η| 2 ) ω 2 .
Once again, the last two rows are factorized after a few calculations that we skip, and we get

A r = 2 c 2 r ρ r [ρ] Υ ( η 2 0 + u 2 r | η| 2 ) γ 2 u ℓ a ℓ (u r a r + i c 2 r η 0 ) β - 2 + u ℓ c 2 r | η| 2 (u r a ℓ + i c 2 ℓ η 0 ) - 1 ρ r [ρ] [u] Υ γ 2 (i η 0 + u r β - 2 ) 2 u ℓ u r a ℓ a r (u r a r + i c 2 r η 0 ) ( η 2 0 + u 2 r | η| 2 ) .
When we add the latter term with the expression in (55), we obtain the constant term (independent of (k, k ′ )):

-

[ρ] [u] Υ γ 2 |i η 0 + u r β - 2 | 2 ( η 2 0 + u 2 r | η| 2 ) i u r η 0 2 a r (u ℓ a ℓ + i c 2 ℓ η 0 ) = 1 2 Q r .
The only term that depends on (k, k ′ ) arises in (56) and equals

k ′ k γ 2 |i η 0 + u r β - 2 | 2 ω 2 (i η 0 + u r β + 2 ) 2 2 c 2 r a r (c 2 r -u 2 r ) = k ′ 2 k γ 2 (i η 0 + u r β + 2 ) 2 η 2 0 + u 2 r | η| 2 a r ω 2 = k ′ 2 k Q r .
The sum of all the terms that involve p ′′ (ρ r ) factorizes as claimed in Proposition 3.

• Step 3: computation of the remaining terms. In order to prove Proposition 3, we can assume from now on, and without loss of generality that p ′′ (ρ ℓ ) = p ′′ (ρ r ) = 0 in (34) and (35). With this simplification, we compute

d-1 k=1 ηk d 2 f k (v ℓ ) • (r - 1 , r + 1 ) = -i c 4 ℓ | η| 2 ρ ℓ   0 2 i η -( β + 1 + β - 1 )   , (57)
d 2 f d (v ℓ ) • (r - 1 , r + 1 ) = c 4 ℓ ρ ℓ     0 i ( β + 1 + β - 1 ) η -2 β + 1 β - 1 -2 u ℓ β + 1 β - 1     , (58)
d-1 k=1 ηk d 2 f k (v r ) • (r - 2 , r + 2 ) = -i c 4 r | η| 2 ρ r   0 2 i η β + 2 + β - 2   , (59)
d 2 f d (v r ) • (r - 2 , r + 2 ) = - c 4 r ρ r     0 i ( β + 2 + β - 2 ) η 2 β + 2 β - 2 2 u r β + 2 β - 2     , (60)
As was done in Step 3 of the proof of Proposition 2, we are going to split the computations between the 'left' and 'right' contributions. Let us first concentrate on all 'left' terms. The 'left' contribution of q 2 is computed from (58) and equals

σ * d 2 f d (v ℓ ) • (r - 1 , r + 1 ) |γ 1 | 2 = σ * c 4 ℓ ρ ℓ     0 i ( β + 1 + β - 1 ) η -2 β + 1 β - 1 -2 u ℓ β + 1 β - 1     |γ 1 | 2 = 2 c 4 ℓ ρ ℓ (c 2 ℓ -u 2 ℓ ) Υ γ 1 u ℓ η 0 | η| 2 γ 1 Ď -( η 2 0 -c 2 ℓ | η| 2 ) γ 1 (D d+1 + u ℓ D d+2 ) .
Using Lemma 4, we find that the 'left' contribution of q 2 is given by

2 c 4 ℓ ρ ℓ (c 2 ℓ -u 2 ℓ ) [ρ] [u] Υ γ 1 i c 2 r η 0 ( η 2 0 -c 2 ℓ | η| 2 ) ( η 2 0 + u 2 r | η| 2 ) -u ℓ u r c 2 ℓ | η| 4 (u r a r -i c 2 r η 0 ) . (61)
Similarly, we use (60) to compute the 'right' contribution of q 2 , which yields

-σ * d 2 f d (v ℓ ) • (r - 2 , r + 2 ) |γ 2 | 2 = 2 c 4 r ρ r (c 2 r -u 2 r ) [ρ] [u] Υ γ 2 × (i c 2 ℓ η 3 0 + u ℓ a ℓ u 2 r | η| 2 ) ( η 2 0 -c 2 r | η| 2 ) -η 2 0 u 2 r | η| 2 (u ℓ a ℓ -i c 2 ℓ η 0 ) . (62)
Using (57), the 'left' contribution of q 3 equals i

(k + k ′ ) +∞ 0 ω 1 γ 1 ℓ 1 e -(k+k ′ ) β + 1 z -i c 4 ℓ | η| 2 ρ ℓ   0 2 i η -( β + 1 + β - 1 )   |γ 1 | 2 e k β - 1 z e k ′ β + 1 z dz + i (k + k ′ ) +∞ 0 ω 3 γ 1 ℓ 3 e -(k+k ′ ) β + 3 z -i c 4 ℓ | η| 2 ρ ℓ   0 2 i η -( β + 1 + β - 1 )   |γ 1 | 2 e k β - 1 z e k ′ β + 1 z dz = c 4 ℓ ρ ℓ γ 1 ω 1 ℓ 1   0 2 i η -( β + 1 + β - 1 )   | η| 2 β + 1 -β - 1 1 + k ′ k + c 4 ℓ ρ ℓ γ 1 ω 3 ℓ 1   0 2 i η -( β + 1 + β - 1 )   (k + k ′ ) | η | 2 k ( β + 3 -β - 1 ) + k ′ ( β + 3 -β + 1 ) ⋄ , (63)
and, similarly (using now (59) rather than (57)), the 'right' contribution of q 3 equals i

(k + k ′ ) +∞ 0 ω 2 γ 2 ℓ 2 e -(k+k ′ ) β + 2 z -i c 4 r | η| 2 ρ r   0 2 i η β + 2 + β - 2   |γ 2 | 2 e k β - 2 z e k ′ β + 2 z dz = c 4 r ρ r γ 2 ω 2 ℓ 2   0 2 i η β + 2 + β - 2   | η| 2 β + 2 -β - 2 1 + k ′ k . (64)
The 'left' contribution of q 4 is computed by retaining the three first rows in (58):

- +∞ 0 ω 1 γ 1 ℓ 1 e -(k+k ′ ) β + 1 z c 4 ℓ ρ ℓ   0 i ( β + 1 + β - 1 ) η -2 β + 1 β - 1   |γ 1 | 2 (k β - 1 + k ′ β + 1 ) e k β - 1 z e k ′ β + 1 z dz - +∞ 0 ω 3 γ 1 ℓ 1 e -(k+k ′ ) β + 3 z c 4 ℓ ρ ℓ   0 i ( β + 1 + β - 1 ) η -2 β + 1 β - 1   |γ 1 | 2 (k β - 1 + k ′ β + 1 ) e k β - 1 z e k ′ β + 1 z dz = - c 4 ℓ ρ ℓ γ 1 ω 1 ℓ 1   0 i ( β + 1 + β - 1 ) η -2 β + 1 β - 1   1 β + 1 -β - 1 β - 1 + k ′ k β + 1 - c 4 ℓ ρ ℓ γ 1 ω 3 ℓ 1   0 i ( β + 1 + β - 1 ) η -2 β + 1 β - 1   k β - 1 + k ′ β + 1 k ( β + 3 -β - 1 ) + k ′ ( β + 3 -β + 1 ) ⋄ , (65)
and the 'right' contribution of q 4 equals

- c 4 r ρ r γ 2 ω 2 ℓ 2   0 i ( β + 2 + β - 2 ) η 2 β + 2 β - 2   1 β + 2 -β - 2 β - 2 + k ′ k β + 2 = c 4 r ρ r γ 2 ω 2 ℓ 2   0 i( β + 2 + β - 2 ) η 2 β + 2 β - 2   - c 4 r ρ r γ 2 ω 2 ℓ 2   0 i( β + 2 + β - 2 ) η 2 β + 2 β - 2   β + 2 β + 2 -β - 2 1 + k ′ k . (66)
The 'left' contribution of q 3 + q 4 is computed as follows. We first observe that the difference between the diamond term in (63) and the diamond term in (65) reads

c 2 ℓ ρ ℓ γ 1 ω 3 u ℓ | η| 2 c 2 ℓ | η| 2 -c 2 ℓ β + 1 β - 1 -(i η 0 -u ℓ β + 1 ) (i η 0 -u ℓ β - 1 ) = 2 c 6 ℓ ρ ℓ (c 2 ℓ -u 2 ℓ ) u ℓ | η| 4 γ 1 ω 3 .
Consequently the 'left' contribution of q 3 + q 4 , which corresponds to the sum of (63) and (65), equals

c 4 ℓ γ 1 ω 1 ρ ℓ ( β + 1 -β - 1 )    | η| 2 ℓ 1   0 2 i η -( β + 1 + β - 1 )   -β + 1 ℓ 1   0 i ( β + 1 + β - 1 ) η -2 β + 1 β - 1      1 + k ′ k + c 4 ℓ ρ ℓ γ 1 ω 1 ℓ 1   0 i ( β + 1 + β - 1 ) η -2 β + 1 β - 1   + 2 c 6 ℓ ρ ℓ (c 2 ℓ -u 2 ℓ ) u ℓ | η| 4 γ 1 ω 3 = 2 c 4 ℓ ρ ℓ ( β + 1 -β - 1 ) γ 1 ω 1 ( β + 1 β - 1 -| η| 2 ) (( β + 1 ) 2 -| η| 2 ) 1 + k ′ k - 2 c 4 ℓ ρ ℓ (c 2 ℓ -u 2 ℓ ) γ 1 ω 1 ( η 2 0 -c 2 ℓ | η| 2 ) β + 1 + i u ℓ η 0 | η| 2 + 2 c 6 ℓ ρ ℓ (c 2 ℓ -u 2 ℓ ) u ℓ | η| 4 γ 1 ω 3 .
so that (70), which is the relation we wish to prove, reads B ℓ + B r = 0. Using the relations u ℓ a r + i c 2 r η 0 u ℓ a r -i c 2 r η 0 = -u r a ℓ + i c 2 ℓ η 0 u r a ℓ -i c 2 ℓ η 0 , (u r a ℓ -i c 2 ℓ η 0 ) γ 1 = -ρ ℓ [u] η 0 , we get u ℓ (u r a ℓ -i c 2 ℓ η 0 ) B ℓ = -i η 0 (u ℓ a ℓ + i c 2 ℓ η 0 ) .

Similarly, we compute u r (u ℓ a r -i c 2 r η 0 ) B r = -i η 0 (u r a r + i c 2 r η 0 ) .

Combining the relations (71) and (72), we get u ℓ u r (u r a ℓ -i c 2 ℓ η 0 ) (u ℓ a r -i c 2 r η 0 ) (B ℓ + B r ) = -i η 0 u r (u ℓ a r -i c 2 r η 0 ) (u ℓ a ℓ + i c 2 ℓ η 0 ) -i η 0 u ℓ (u r a ℓ -i c 2 ℓ η 0 ) (u r a r + i c 2 r η 0 ) = -i η 0 (u ℓ + u r ) (u ℓ u r a ℓ a r + c 2 ℓ c 2 r η 2 0 ) = 0 .

This means that (70) holds, and consequently the expression of the kernel q is as claimed in Proposition 4. The verification of Hunter's stability condition q(1, 0 + ) = q(1, 0 -) is then straightforward, and local well-posedness in H 2 (R) follows from the main result in [START_REF] Benzoni-Gavage | Local well-posedness of nonlocal Burgers equations[END_REF].

The main difference here with[START_REF] Benzoni-Gavage | Weakly nonlinear surface waves and subsonic phase boundaries[END_REF] Paragraph 3.3] is our sign convention for a ℓ , the latter being denoted by α ℓ in[START_REF] Benzoni-Gavage | Weakly nonlinear surface waves and subsonic phase boundaries[END_REF] Paragraph 3.3].

Since a ℓ is negative and a r is positive, D d+1 is nonzero and we can thus check that σ is a nonzero vector.

We keep the notation of[START_REF] Benzoni-Gavage | Weakly nonlinear surface waves and subsonic phase boundaries[END_REF] for the functions r±, r and so on.

We have thus shown that the 'left' contribution A ℓ of q 2 + q 3 + q 4 reads

The final step of the proof consists in simplifying the remaining terms in A ℓ and A r . More specifically, the first factor in the expression (51) of A ℓ can be simplied as follows:

Similarly, the first factor in the expression (52) of A r can be simplied by using

Using these two last simplifications, we can add (51) and (52) and obtain

with Q ♯ defined in (33). This completes the proof of Proposition 2.

We now compute the kernel q 2 + q 3 + q 4 in the domain {k > 0, k ′ > 0, k + k ′ > 0}.

Proposition 3. Let Q ℓ and Q r be defined in (33) and let us define

Then the kernels q 2 , q 3 , q 4 defined in (28), ( 29) and (30) satisfy

Proof. We split again the proof in several steps, as was done in the proof of Proposition 2.

• Step 1: computation of the p ′′ (ρ ℓ ) factor. We collect again the contributions that involve p ′′ (ρ ℓ ). The contribution of the kernel q 2 equals Υ

where the product (D d+1 + u ℓ D d+2 ) γ 1 has already been computed when deriving (36).

The contribution of the kernel q 3 equals i

.

Similarly, the contribution of the kernel q 4 reads

.

Adding the contributions of q 3 and q 4 gives the term

We now add the contributions in ( 53) and (54) in order to obtain the p ′′ (ρ ℓ ) term in q 2 + q 3 + q 4 . There is first a constant term that is independent of (k, k ′ ), and this term is entirely similar to the one derived when adding (36) and (37) (see Step 1 in the proof of Proposition 2). Namely, the constant term in the sum of (53) and (54) equals

The last contribution in the p ′′ (ρ ℓ ) term is the one that depends on (k, k ′ ) in (54), that is

We have thus shown that the p ′′ (ρ ℓ ) term in (q 2 + q 3 + q 4 )(k, k ′ ) is as claimed in Proposition 3.

• Step 2: computation of the p ′′ (ρ r ) factor. The contribution of the kernel q 2 equals

After a little bit of algebra, the 'left' contribution of q 3 + q 4 is found to be equal to

When combined with (61), we have thus shown that the 'left' contribution of the kernel q 2 + q 3 + q 4 reads

Let us now compute the 'right' contribution of the kernel q 3 + q 4 . We add the expressions in ( 64) and (66), which gives

The first row in (68) is exactly the opposite of (62), that is of the 'right' contribution of q 2 . In other words, we have found that the 'right' contribution of q 2 + q 3 + q 4 is given by the second row in (68), which reads

The kernel q 2 + q 3 + q 4 is the sum of the latter quantity and the right hand side of (67). Collecting the terms, this completes the proof of Proposition 3.

Conclusion

There remains a final simplification in order to achieve the final form of the kernel q. Our main result reads as follows.

Proposition 4. With Q ℓ , Q r , Q ♯ defined in Propositions 2, let us define

Then the kernel q = 4 π a 1 satisfies

In particular, q satisfies Hunter's stability condition q(1, 0 + ) = q(1, 0 -), and (15) is locally well-posed in H 2 (R).

Proof. In view of the expression of q given in Corollary 2, we only need to show the relation

) u ℓ u r , we are led to showing that the following relation holds:

where we use the notation j := ρ ℓ u ℓ = ρ r u r to denote the mass flux across the phase boundary. The verification of (69) proceeds as follows. We first combine the first and third row in (69) by recalling the definitions of a ℓ , a r , see ( 1) and (2). Thus verifying (69) amounts to showing

Let us now define


Keywords: resilience, fault tolerance, eigenvalue problems, linear algebra, HPC, numerical methods, Arnoldi, IRAM, subspace iteration, Jacobi-Davidson

 for linear systems. For a given numerical scheme, the IR strategies consist in extracting relevant spectral information from available data after a fault. After data extraction, a well selected part of the missing data is regenerated through interpolation strategies to constitute meaningful input to restart the numerical algorithm. A main feature of this numerical remedy that it does not require extra resources, i.e., computational unit or computing time, when no fault occurs. In this paper, we revisit a few state-of-the-art methods for solving large sparse eigenvalue problems namely the Arnoldi methods, subspace iteration methods and the Jacobi-Davidson method, in the light of our IR strategies. For each considered eigensolver, we adapt the IR strategies to regenerate as much spectral information as possible. Through intensive numerical experiments, we illustrate the qualitative behavior of the resulting schemes when the number of faults and the amount of lost data are varied.

Stratégies d'interpolation-restart pour des solveurs aux valeurs propres résilients

Résumé : Le calcul des paires propres (valeurs propres et vecteurs propres) de matrices creuses de grande taille intervient dans de nombreuses applications scientifiques et d'ingénierie, par exemple dans les problèmes d'analyse de stabilité. Pour des simulations de grande dimension en physique des matériaux ou en thermo-acoustique, ces calculs peuvent durer plusieurs heures même sur des supercalculateurs. Pour les futurs systèmes de calcul haute performance à grande échelle, il est attendu que l'intervalle de temps entre deux pannes consécutives diminue de sorte que plusieurs pannes pouraient intervenir au cours d'une même résolution. Il est donc nécessaire de concevoir des solveurs parallèles de problèmes aux valeurs propres qui survivent aux pannes. Dans ce contexte, nous considérons essentiellement des solutions numériques pour en évaluer la pertinence, celles-ci pourraient se combiner avec des techniques de résilience plus classiques pour des mises en oeuvre efficaces qui passent à l'échelle. Afin de se concentrer sur des solutions numériques, nous faisons l'hypothèse qu'un mécanisme distinct assure la détection des pannes et qu'une couche système s'occupe de ramener l'environnement (processus, . . . ) dans un état fonctionnel. Une fois que le système est de nouveau opérationnel, après une panne, notre principal objectif est d'étudier des stratégies robustes de résilience qui permettent de converger sans avoir à redémarrer complètement le calcul. À cette fin, nous étendons les techniques Interpolation-Restart (IR) introduites dans un précédent travail [START_REF] Agullo | Towards resilient parallel linear Krylov solvers: recover-restart strategies[END_REF] pour les systèmes linéaires. Les stratégies IR consistent à regéner des valeurs pour les données perdues en utilisant les données qui sont disponibles; l'ensemble de ces données est utilisé pour redemarrer efficacement les différents solveurs. Un atout de ces approches est qu'elles ne nécessitent pas de ressources supplémentaires (calculatoire ou mémoire) en l'absence de pannes. Dans ce travail, à la lumière des stratégies d'interpolation-restart, nous avons revisité des méthodes itératives de l'état de l'art à savoir les méthodes Arnoldi, les méthodes d'itérations de sous-espace et la méthode de Jacobi-Davidson. Pour chacune de ces méthodes, nous avons adapté les stratégies d'interpolation-restart pour régénérer autant d'information spectrale que possible. Nous illustrons, grâce à de nombreuses expériences, le comportement numérique qualitatif de nos stratégies lorsque le nombre de pannes et la quantité de données perdues varient.

Mots-clés :

résilience; tolérance aux pannes; problèmes aux valeurs propres; algèbre linéaire; HPC; méthodes numériques; Arnoldi; IRAM; itération de sous-espace; Jacobi-Davidson.

Introduction

The computation of eigenpairs (eigenvalues and eigenvectors) of large sparse matrices is involved in many scientific and engineering applications such as when stability analysis is a concern. To name a few, it appears in structural dynamics, thermodynamics, thermo-acoustics, quantum chemistry. With the permanent increase of the computational power of high performance computing (HPC) systems by using a larger and larger number of CPU cores or specialized processing units, HPC applications are increasingly prone to faults. To cope with these unstable situations, numerical simulations need to be equipped with resilient mechanisms that enable them to complete a calculation even though several faults might occur during its realization. We investigate in this paper numerical remedies in the framework of the solution of eigenproblems. We consider the standard eigenproblem of the form:

Au = λu,
where A ∈ C n×n , with u = 0, u ∈ C n and λ ∈ C. The couple (λ, u) is called an eigenpair of A where the vector u is an eigenvector with the associated eigenvalue λ.

In this paper, we extend the interpolation-restart (IR) strategies introduced for the solution of linear systems [START_REF] Agullo | On the resilience of a parallel sparse hybrid solver[END_REF][START_REF] Agullo | Towards resilient parallel linear Krylov solvers: recover-restart strategies[END_REF] to a few state-of-the-art eigensolvers. More precisely, the Arnoldi [START_REF] Arnoldi | The principle of minimized iterations in the solution of the matrix eigenvalue problem[END_REF], Implicitly restarted Arnoldi [START_REF] Lehoucq | Deflation techniques for an implicitly re-started Arnoldi iteration[END_REF], subspace iteration [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems[END_REF], and Jacobi-Davidson algorithms have been revisited to make them resilient in the presence of faults. Most of the considered eigensolvers naturally implement a restart mechanism to cope with memory constraints. We incorporate an additional restart in our scheme after faults. We attempt to regenerate as much as possible relevant spectral information to perform the restart effectively and efficiently.

The remainder of the paper is structured as follows: in Section 2 we describe how the interpolation techniques can be extended to regenerate meaningful spectral information. We briefly present the eigensolvers that we have considered in Section 3 and detail how the recovery ideas can be tuned for each one of them. Section 4 is devoted to the numerical experiments; we discuss the robustness of the various resilient numerical schemes and conclude with some perspectives in Section 5.

Interpolation-restart principles

In this section, we describe how the interpolation strategies can be used to regenerate meaningful spectral information. Contrarily to what has been proposed for the Krylov linear solvers where only a meaningful iterate is computed to serve as a new initial guess for restarting the iterations [START_REF] Agullo | On the resilience of a parallel sparse hybrid solver[END_REF][START_REF] Agullo | Towards resilient parallel linear Krylov solvers: recover-restart strategies[END_REF], more flexibility exists in the framework of eigensolution where similar ideas can be adapted to better exploit the numerical features of the individual eigensolvers. The main reasons are that some of the considered eigensolvers do not rely on a central equality or a sophisticated short term recurrence (such as Conjugate Gradient for linear system solution). Furthermore, we consider also situations where a few nev eigenpairs are sought, which also provides additional freedom. We present in details different variants for selecting and computing the relevant subspaces to perform the restart for each particular considered eigensolver in Section 3.

Context

We introduce the governing ideas that underline the design of the IR strategies. For the sake of exposure, we restrict ourselves to parallel distributed environments, although these strategies can be extended to other HPC contexts, assuming that the faults can be detected and reported by a lower layer of the operating system stack [START_REF] Jaulmes | Exploiting asynchrony from exact forward recovery for DUE in iterative solvers[END_REF].

Assumption 1. In our parallel computational context, all the vectors or matrices of dimension n are distributed by blocks of rows in the memory of the different computing nodes but scalars or low dimensional matrices are replicated.

Inria

According to Assumption 1, the eigenvector u and the matrix A are distributed according to a block-row partition as well as all vectors of dimension n generated during the solution whereas scalars (for example λ) or low dimensional matrices are replicated on all nodes. Let N be the number of partitions, such that each block-row is mapped to a computing node. For all p, p ∈ [1, N ], I p denotes the set of row indices mapped to node p. With respect to this notation, node p stores the block-row A Ip,: and u Ip as well as the entries of all the vectors involved in the solver associated with the corresponding row indices of this block-row. If the block A Ip,Iq contains at least one non zero entry, node p is referred to as neighbor of node q as communication will occur between those two nodes to perform a parallel sparse matrix-vector product.

When a fault occurs on a node, all data in its memory are lost. We consider the formalism proposed in [START_REF] Langou | Recovery Patterns for Iterative Methods in a Parallel Unstable Environment[END_REF] where lost data are classified into three categories: the computational environment, the static data and the dynamic data. The computational environment is all the data needed to perform the computation (code of the program, environment variables, . . . ). The static data are those that are set up during the initialization phase and that remain unchanged during the computation such as the coefficient matrix A. The dynamic data is all data whose value may change during the computation. The iterates λ and u are examples of dynamic data. In Figure 1a, we depict a block-row distribution on four nodes. The data in blue is the static data associated with the eigenproblem (i.e., the matrix) while the data in green is the dynamic data (here, only the eigenpair (λ, u) is shown). If node P 1 fails, the first block-row of A as well as the first entries of u are lost (in black in Figure 1b). The eigenvalue λ is a scalar replicated on all nodes and thus remains available on non failed ones. Figure 1: General interpolation scheme. The matrix is initially distributed with a block row partition, here on four nodes (a). When a fault occurs on node P 1 , the corresponding data is lost (b). Whereas static data can be immediately restored, dynamic data that has been lost cannot and we investigate numerical strategies for interpolating it (c).

We assume that, when a fault occurs, the failed node is replaced and the associated computational environment and static data are restored. In Figure 1c for instance, the first matrix block-row is restored as it is a static data. However, the eigenvector iterate u, being a dynamic data, its entries u I1 are definitely lost and we present strategies for regenerating them through some interpolations that follow. Our strategies do not attempt to interpolate all the dynamic data but for each eigensolver, we study which part of the dynamic data could be interpolated. Because not all data is regenerated and/or that the interpolation computed approximation of the lost data, part of the information remains lost. As a consequence, it is furthermore required to restart the numerical process. We will present in detail propositions for selecting (or producing) relevant data for performing that restart for each particular eigensolver considered in Section 3. We assume in the rest of this paper that a fault occurs during iteration f + 1 and the proposed interpolation strategies are thus based on the numerical values computed by the algorithms at iteration f .

Interpolation methods

The IR strategies consist in interpolating lost data by using non corrupted data. Let u (f ) be an approximated eigenvector when a fault occurs. After the fault, the entries of u (f ) are correct, except those stored on the failed node p. Assuming that in a parallel distributed environment, the current eigenvalue λ f is naturally replicated in the memory of the different computing nodes, we present two strategies to compute a new approximate solution. The first strategy, referred to as linear interpolation and denoted LI, consists in solving a local linear system associated with the submatrices A Ip,Ip of the failed node. The second one relies on the solution of a least squares interpolation and is denoted LSI. Those two alternatives result from considering (λ f , u (f ) ) as an exact eigenpair. We may have a block-row view point, which defines the LI variant. If node p fails, LI computes a new approximation of the eigenvector u (LI) as follows

     u (LI) Iq = u (f ) Iq for q = p, u (LI) Ip = (A Ip,Ip -λI Ip,Ip ) -1 (- q =q A Ip,Iq u (f ) Iq ).
Alternatively, we can have a block-column point of view, which leads to the LSI variant that computes u (LSI) via

     u (LSI) Iq = u (f ) Iq for q = p, u (LSI) Ip = argmin u Ip (A :,Ip -λI :,Ip )u Ip - q =p (A :,Iq -λI :,Iq )u (f ) Iq ) .
Here, I ∈ C n×n is the identity matrix and we furthermore assume that (A Ip,Ip -λI Ip,Ip ) is non singular and that (A :,Ip -λI :,Ip ) has full rank. The matrix involved in the least squares problem, (A :,Ip -λI :,Ip ), is sparse of dimensions |J p | × |I p | where its number of rows |J p | depends on the sparsity structure of A :,Ip . Consequently the LSI strategy may have a higher computational cost. In the rest of this paper, we use IR indifferently to denote LI and/or LSI.

Interpolation-restart strategies for general purpose eigensolvers

In this section, we briefly describe a few of the most popular eigensolvers for general sparse matrices that we have considered in this work to compute nev eigenpairs. For each eigensolver, we describe which data are regenerated after a fault to make them resilient. We also briefly present the numerical approach and associated algorithm. We then describe how the IR strategies can be applied to compute a few of the key numerical quantities of the solver when some data is lost as well as how this data can be used to perform an effective and efficient restart.

Subspace iterations to compute nev eigenpairs

Brief description of the numerical algorithm: The subspace iteration method is a block variant of the power method [START_REF] Stewart | Matrix algorithms -Volume II: Eigensystems[END_REF]. It starts with an initial block of nev + s (with s ≥ 0) linearly independent vectors corresponding to matrix

U (0) = [u (0) 1 , . . . , u (0)
nev+s ] ∈ C n×(nev+s) . Under certain assumptions [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems[END_REF], the sequence of U (k) = A k U (0) generated by the algorithm, converges to the nev + s eigenpairs of A associated with the eigenvalues of the largest magnitude. To guarantee the full column rank in U (k) for large values of k, the Q factor of its QR factorization may be used at each iteration.

Inria

Algorithm 1 Basic subspace iteration

1: Choose U (0) = [u (0) 1 , ..., u (0) nev+s ] ∈ C n×(nev+s) 2: for k = 0 . . . until convergence do 3: Orthonormalize U (k) 4: Compute W (k) = AU (k) 5: Form the Rayleigh quotient C (k) = W (k) H AW (k) 6:
Compute the eigenvectors G (k) = [g 1 , . . . , g nev+s ] of C (k) and eigenvalues σ(C (k) ) = (λ 1 , . . . , λ nev+s )

7:

Update Ritz vectors :

U (k+1) = W (k) G (k) 8: end for
To compute the eigenpairs associated with the smallest eigenvalues in magnitude, or eigenpairs associated with a given set of eigenvalues in a given region of the complex plane, the basic subspace iteration depicted in Algorithm 1 is no longer appropriate. For example, to compute eigenpairs associated the eigenvalues nearest to σ ∈ C, it is possible to combine the subspace iterations with the shift-invert technique [START_REF] Freitag | Shift-invert Arnoldi's method with preconditioned iterative solves[END_REF]. With shift-invert spectral transformation, the subspace iteration method will be applied to the matrix (A -σI) -1 . Alternatively, a polynomial acceleration [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF], can be used as a preconditioning technique to approximate eigenvalues near σ. This polynomial acceleration consists in applying a given polynomial P to the matrix A. In Algorithm 1 line 4 would change and become W (k) = P(A)U (k) . The polynomial should act as a filter to damp eigencomponents in some undesired part of the spectrum.

Interpolation-restart policy: In the subspace iteration method depicted in Algorithms 1, according to Assumption 1, the Ritz vectors U (k) are distributed, whereas the Rayleigh quotient C (k) and Ritz values are replicated. When a fault occurs, we distinguish two cases. During an iteration, a fault may occur before or after the computation of the Rayleigh quotient C (f +1) .

1. When a fault occurs before the computation of the Rayleigh quotient C (f +1) (Algorithm 1, lines 2 to 5) the surviving nodes cannot compute the Rayleigh quotient C (f +1) because entries of W (f +1) are missing. In this case, we consider the available entries of the Ritz vectors U (f ) and its corresponding eigenvalues σ(C (f ) ). We interpolate the m Ritz vectors individually (u (f ) , 1 ≤ ≤ nev + s) using LI or LSI. In the particular case of the computation of polynomial acceleration, all the computation in the filtering step are lost.

2. When a fault occurs after the computation of the Rayleigh quotient C (f +1) (Algorithm 1, line 6 to 7) all surviving nodes can compute the entries of U (f +1) relying on a local replicate of C (f +1) and the local entries of W (f +1) . The missing entries of each Ritz vector (u (f +1) , 1 ≤ ≤ nev + s) can be individually interpolated using LI or LSI relying on the corresponding eigenvalues σ(C (f +1) ).

After the interpolation, the subspace iteration algorithm is restarted with the matrix

U (IR) = [u (IR) 1 , . . . , u (IR)
nev+s ] ∈ C n×(nev+s) until convergence or the next fault.

Arnoldi method to compute one eigenpair

Brief description of the numerical algorithm: The Arnoldi method is an efficient procedure for approximating eigenvalues lying in the periphery of the spectrum of A. For a prescribed dimension m, referred to as the restart, the method starts from an initial guest vector v 1 of norm one and builds an orthonormal basis V m of the search space spanned by span{v 1 , Av 1 , ..., A m-1 v 1 } as well as an upper Hessenberg matrix H m ∈ C m×m . Those matrices satisfy the so called Arnoldi equality:

AV m = V m H m + β m v m+1 e T m , (1)
where e m denotes the last column of the m × m identity matrix. By construction, H m is the Rayleigh quotient matrix associated with V m that is used to compute approximate eigenpairs of A. If the computed Ritz pair corresponding to the targeted eigenvalue is not accurate enough, the procedure is restarted using the best current Ritz vector as the new initial guess. This space expansion and restart sequence is repeated until convergence.

Interpolation-restart policy: According to Assumption 1, we assume that the low dimension Hessenberg matrix H k is replicated on each node. Consequently, regardless of the step in which the fault occurs during the iteration, each surviving node q can solve the eigenproblem H f g = λg redundantly, then compute its entries of the Ritz vector u Iq = V f (I p , :)g. The next step is the computation of the Ritz vector entries allocated on the failed node using LI or LSI. The resulting vector becomes the new initial guess to restart the Arnoldi iterations.

Implicitly restarted Arnoldi method to compute nev eigenpairs

Brief description of the numerical algorithm: Developed by Lehoucq and Sorensen in [START_REF] Lehoucq | Deflation techniques for an implicitly re-started Arnoldi iteration[END_REF], the implicitly restarted Arnoldi method (IRAM) depicted in Algorithm 2 is commonly used for the solution of large eigenvalue problems. IRAM is an extension of the Arnoldi method that starts with an Arnoldi equality of size m. From the Arnoldi equality of dimension m, IRAM performs a contraction of the equality from size m down to size m (nev ≤ m ≤ m). This is achieved by applying a polynomial filter of degree = m -m that reduces the size of the Arnoldi equality down to m (see Algorithm 2, line 12), that is then expanded again to size m before checking the accuracy of the Ritz eigenpairs. The expansion and contraction steps are repeated until convergence. The contraction step acts as a filter to focus the spectral information in a target region of the spectrum from the m dimensional Krylov subspace while maintaining the Arnoldi equality.

Algorithm 2 Implicitly restarted Arnoldi method with restart m

1: Compute Arnoldi equality AV m = V m H m + f m e T m . 2: for k = 0, 1, . . . , until convergence, do 3:
Compute σ(H m ) and select shifts µ 1 , . . . , µ ( = m -m).

4:

Q = I m 5: for i = 1, . . . , do 6: QR Factorize Q i R i = H m -µ i I 7: H m = Q i H H m Q i 8: Q = QQ i 9:
end for 10:

β m = H m ( m + , m) 11: f m = v m+ β m + f m Q(m, m) 12: V m = V m Q(:, 1 : m); H m = H m (1 : m, 1 : m) 13:
Starting with 1. The contraction implicitly implements the shifted QR mechanism. When a fault occurs during the contraction step each surviving node q can compute this step concurrently and redundantly using the replicated Hessenberg matrix H m . The outcome of this local calculation is that all surviving nodes have V m(I q , :) as well as the reduced Hessenberg matrix H m.

AV m = V mH m + f me T m, perform steps of Arnoldi algorithm to get AV m = V m H m + f m e T

Inria

2. When the fault occurs in the expansion phase during step f + 1 with m < f + 1 < m, each surviving node can compute the implicitly shift QR update of H f to compute H m using the shift defined by σ(H f ). Using the result of the shift QR calculations, each node q can also compute V m(I q , :).

From V m(I q , :) and H m, the surviving nodes may then compute eigenvectors G = [g 1 , . . . , g nev ] and eigenvalues (λ 1 , . . . , λ nev ) of H m. Consequently, the entries of the Ritz vectors are computed by

U (f ) (I q , :) = V m(I q , :)G. (2)
The missing entries I p of Ritz vectors may be interpolated using either LI or LSI of the interpolation algorithms. In that situation, the available entries of V m no longer satisfy any Arnoldi equality.

To take into account all the available spectral information, we compute the linear combination of the interpolated eigenvectors, u = m j=1 u

(IR) j

, and restart with the normalized linear combination v 1 = u u as initial vector. The motivation is that Arnoldi converges within the first expansion if started from a vector v that is a linear combination of k ≤ m eigenvectors [START_REF] Stewart | Matrix algorithms -Volume II: Eigensystems[END_REF].

3.4

The Jacobi-Davidson method to compute nev eigenpairs Brief description of the numerical algorithm: The Jacobi-Davidson method, proposed by Sleijpen and van der Vorst in [START_REF] Gerard | A Jacobi-Davidson iteration method for linear eigenvalue problems[END_REF], is a widely used eigensolver, especially for eigenpairs in the interior of the spectrum. The basic ingredient of Jacobi-Davidson is depicted in Algorithm 3 for the computation of one eigenpair whose eigenvalue is close to a given target τ . It starts with a given normalized vector v and constructs a basis V extended using the Jacobi orthogonal correction method. At each iteration the algorithm computes the Ritz pairs associated with V and selects the eigenpair whose eigenvalue is the closest to the target τ . Select Ritz pair (λ k , u k ) such that λ k is the closest to τ 5:

r k = Au k -λ k u k 6:
Perform a few steps of GMRES to solve the correction equation

(I -u k u H k )(A -τ I)(I -u k u H k )v = -r k , so that v ⊥ u k 7:
Compute w by orthonormalizing v against

V k : w = v -V k (V H k v) 8:
Set V k+1 = [V k , w] 9: end for Algorithm 3 converges to one eigenpair. If more than one eigenpair needs to be computed, Algorithm 3 can be accommodated to compute a partial Schur decomposition of A. In that respect, the next iterations are enforced to generate a search space orthogonal to the space spanned by the nconv already converged eigenvectors. This is achieved by representing this space using the corresponding Schur vectors. Let AZ nconv = Z nconv T nconv denote the partial Schur form where the columns of the orthonormal matrix Z nconv span the converged eigenspace and the diagonal of the upper triangular matrix T nconv are the associated converged eigenvalues.

Algorithm 4 corresponds to the Jacobi-Davidson style QR algorithm presented in [START_REF] Diederik | Jacobi-Davidson style QR and QZ algorithms for the partial reduction of matrix pencils[END_REF]. It is conceived to be used by a higher level routine that decides the number of wanted eigenpairs nev, the target point τ , the maximum and the minimum size of the basis V , etc. The inputs for the algorithm are the existing converged Schur vectors Z nconv of A, the current size k of the basis V k , W k = AV k , and the Rayleigh quotient C k . The focal point τ , the maximum dimension m affordable for the space spanned by V , the size of the restarted basis m (1 ≤ m < m) and the maximum number of restarts allowed are also provided. Outputs are µ, z and t, such that

A(Z nconv z) = (Z nconv z) T nconv t 0 µ
is a partial Schur decomposition of one higher dimension. for k = k init , . . . , m do 5:

% Computation of the Schur decomposition CQ = QT % so that the eigenvalues on the diagonal of T % are sorted by increasing distance to τ

[Q, T ] = SortSchur(C(1 : k, 1 : k), τ, k), 6:
Choose µ = T (1, 1) and g = Q(:, 1), the Ritz pair closest to τ 7:

Approximate eigenvector of A: z = V (:, 1 : k)g, and Az: y = W (:, 1 : k)g

8:

Compute the residual r = y -µz, orthogonalize it against Z(:,1:nconv) and compute its norm: rnorm=norm(r)

9:
% Convergence test:

10:
if rnorm is small enough then 11:

nconv = nconv + 1 12:
% Prepare outputs and deflate: %No convergence reached and k < m. Solve the correction equation:

13: t = Z H y; V = V (:, 1 : k)Q(:, 2 : k); W = W (:, 1 : k)Q(:, 2 : k); C = T (2 : k, 2 : k).
(I -zz H )(A ⊥ -τ I)(I -zz H )v = -r 19:
Orthogonalize v against V (:, 1 : k) and Z(:, 1 : nconv) 20:

% Extend the Rayleigh basis and the Rayleigh quotient:

21:

V (:, k + 1) = v, W (:, k + 1) = Av, C(k + 1,1:k)=v H W (:,1:k), C(1:k,k + 1)=V (:,1:k) H W (:,k + 1), C(k + 1,k + 1)=v H W (:,k + 1) 22:
end for 23: end while The higher level routine must furnish the necessary inputs to Algorithm 4. If the process starts from the beginning, there are then two situations. The first one corresponds to the case when the computation starts from a single random vector. Then the higher level routine computes an

Inria Arnoldi decomposition of size m AV m = V mH m + βv m+1 e T m,
and Jacobi-Davidson starts with

U = [ ], V = V m, W = AV and C = H m.
The second case is when the process starts from a given number k of initial vectors. The initial block of vectors is then orthonormalized to obtain V k and the process can start as indicated previously, with

Z = [ ], V = V k , W = AV k and C = V H k AV k = V H k W .
Once a partial Schur form of size nev is available, the eigenpairs (λ , u ) (with = 1, . . . , nev) of A can be computed. The eigenvalue λ is the Ritz value of T nev associated with the Ritz eigenvector g so that u = Z nev g .

Interpolation-restart policy: According to Assumption 1, the Schur vectors Z nconv = [z 1 , . . . , z nconv ], and the basis V f = [v 1 , . . . , v f ] are distributed among the computing units as the matrix T nconv ∈ C nconv×nconv , and the Rayleigh quotient matrix C f ∈ C f ×f are replicated.

The Jacobi-Davidson algorithm enables more possibilities to regenerate a meaningful context for the restart after a fault. There are mainly two reasons. First, the algorithm does not rely on an equality that is incremented at each iteration such as Arnoldi; preserving such an incremental equality after a fault is very challenging. Second, the algorithm can start from a set of vectors and its convergence will be fast if these vectors are rich in the sought spectral information.

When the fault occurs on node p while nconv (nconv > 0) Schur vectors were converged, good approximations of the associated converged eigenvectors can easily be computed as follows. Each non-faulty node q performs: 1. the spectral decomposition of the partial Schur matrix T nconv

T nconv G nconv = G nconv D with G nconv = [g 1 , . . . , g nconv ],
2. and the computation of its entries of the converged eigenvectors u (I q ) = Z nconv (I p , :)g for = 1, . . . , nconv.

The missing entries of the eigenvectors can be computed using IR to build U

(IR) nconv = [u (IR) 1 , . . . , u (IR) nconv ].
In addition to U (IR) nconv , further information can be extracted from the search space V f and the Rayleigh quotient matrix C f . Following the same methodology, spectral information built from C f and V f can be computed to generate additional directions to expand the initial search space (U (IR) nconv ) used to restart the Jacobi-Davidson algorithm. Each non-faulty node q computes: 1. the sorted Schur decomposition of C f , that writes C f Gf = Gf S f so that the eigenvalues on the diagonal of S f are sorted by increasing distance to τ , (Algorithm 4, line 5), 2. and the entries of the Ritz vectors ũ (I q ) = V f (I q , :)g for = 1, . . . , s, where s is the number of Ritz vectors we want to interpolate. Because Gf has been sorted, these vectors may be considered as the s best candidates to expand

U (IR) nconv .
That is, ũ1 is the Ritz vector associated with S(1, 1) which is the Ritz value the closest to the target τ , that is improved by Jacobi-Davidson iterations.

In addition, the missing entries I p of the Ritz vectors ũ can be computed using LI or LSI,

U (IR) = [u (IR) 1 , . . . , u (IR) nconv , ũ(IR) 1 , . . . , ũ (IR) s
]. Once U (IR) has been computed, the vectors in U (IR) are then orthonormalized to obtain V restart . The Jacobi-Davidson algorithm can be restarted with nconv then, the nconv already converged Schur vectors will be retrieved in the initial basis V nconv . In floating point arithmetic, there is no guarantee of retrieving the already converged nconv Schur vectors by restarting with V nconv , although in practice this is likely to happen as we will see in the numerical experiments.

Z = [ ], V = V restart , W restart = AV restart , C = V H restart W.
Although in principle s has only to satisfy 0 ≤ s ≤ f, because of Remark 1, a natural choice for s is s = nev (we interpolate nconv + nev vectors) so that the initial search space after a fault will be at least of dimension nev.

Numerical experiments

In this section, we investigate the numerical behavior of the eigensolvers in the presence of faults when the IR policies are applied. In Section 4.1, we present results for subspace iteration methods and study the robustness of the IR strategies when converging eigenvectors associated with both smallest and largest eigenvalues. We assess the robustness of our resilient Arnoldi procedure in Section 4.2, whereas Section 4.3 analyzes the robustness of the resilient algorithm designed for IRAM. In Section 4.4, we discuss results obtained for Jacobi-Davidson and the impact of different variants on the convergence behavior. We define

η(u (k) , λ k ) = Au (k) -λ k u (k) |λ k | , (3)
the scaled residual associated with the approximate eigenpair (λ k , u (k) ) for nonzero eigenvalue approximation. This scaled residual is a lower bound of the normwise backward error of the eigenpair (u (k) , λ k ). Given a threshold ε, the widely used stopping criterion to detect convergence is defined by

η(u (k) , λ k ) ≤ ε.
For all the experiments reported in this section we set ε = 10 -6 . We have performed extensive numerical experiments and only report here on the qualitative numerical behavior observed on a few examples that are representative of our observations. Although many test matrices have been considered for evaluating the qualitative behavior of the resilient schemes, we kept only one example in this section that comes from thermo-acoustic instabilities calculation in combustion chambers [START_REF] Salas | Physical and numerical aspects of thermoacoustic instabilities in annular combustion chambers[END_REF][START_REF] Salas | Spectral recycling strategies for the solution of nonlinear eigenproblems in thermoacoustics[END_REF]. Indeed this test case exhibits many illustrative features. The matrix is unsymmetric; its spectrum lies in the left plane and it has small eigenvalues close to zero that can be computed without shift invert techniques using the different eigensolvers we have considered. Although its size is rather small (n = 1500), it exhibits numerical difficulties that are encountered on real life large scale problems [START_REF] Salas | Physical and numerical aspects of thermoacoustic instabilities in annular combustion chambers[END_REF][START_REF] Salas | Spectral recycling strategies for the solution of nonlinear eigenproblems in thermoacoustics[END_REF].

In order to study the numerical behavior of the IR strategy and illustrate their possible robustness and weaknesses we consider three additional executions in our numerical experiments. The first one is the non faulty (NF) execution case that provides the best expectation. In order to evaluate the impact of the restart and the quality of the interpolated values, we also report on what is referred to as the Enforced Restart (ER) execution. It consists in enforcing the solver to restart at iteration f using the available quantities at this iterations, that are:

• for the subspace iteration, ER is equivalent to NF because this numerical scheme is a fixed point iteration,

Inria

• for regular Arnoldi, it corresponds to classical Arnoldi with variable restarts,

• for IRAM and Jacobi-Davidson, the ER calculation is identical to the IR strategies except that none entries are interpolated; everything is computed using the values available at iteration f , Finally, to illustrate the benefit and robustness of the interpolation policies, we also report on the numerical behavior of the so called Reset variant, where all the quantities interpolated by the IR strategies are simply replaced by random values. Faults are injected using the Weibull probability distribution [START_REF] Weibull | A statistical distribution function of wide applicability[END_REF] that is supposed to be the most realistic probabilistic model that characterizes the normal behavior of large-scale computational platforms. We vary some of its parameters to increase or decrease the number of faults for a given calculation. Finally, for the sake of fair comparison, for a set of experiments (e.g, Figure 2a), faults are injected at the same iterations (e.g, 485 and 968) and during the same instructions for all cases (except NF, of course).

Resilient subspace iteration methods to compute nev eigenpairs

In this section, we analyze the robustness of the proposed resilient IR subspace iteration methods in the presence of faults. To analyze the robustness of our strategies, we simulate stressful conditions by increasing the fault rate and the volume of lost data. We present results for two variants of subspace iteration methods:

1. The subspace iteration with Chebyshev polynomial acceleration is used for the computation of the five eigenpairs corresponding to the smallest eigenvalues (Figures 2 and3). In practice, a certain amount of information about the spectrum is needed in order to build the ellipse associated with the Chebyshev filter polynomial. The ellipse chosen must be as small as possible, but large enough to enclose the unwanted part of the spectrum. We mention that for thermo-acoustic calculations [START_REF] Salas | Physical and numerical aspects of thermoacoustic instabilities in annular combustion chambers[END_REF][START_REF] Salas | Spectral recycling strategies for the solution of nonlinear eigenproblems in thermoacoustics[END_REF] this prerequisite information is computationally affordable.

2. The classical method for the computation of the five eigenpairs corresponding to the largest magnitude eigenvalues (Figures 4 and5).

For both calculations (the five largest and the five smallest eigenvalues), we report the maximum of the individual scaled residual norms (defined by Equation ( 3)) of the five Ritz pairs at each iteration. The execution ends when the five Ritz pairs satisfy the stopping criterion, i.e., when the maximum of the scaled residual norms is lower than the selected threshold ε.

When converging the Ritz pair associated with the smallest Ritz value, the scaled residual norm increases during the first iterations before it begin to decrease slowly reaching the target threshold in 420 iterations as in the NF case (Figures 2 and3). With the reset strategy, the convergence history is characterized by large jumps after each fault that strongly delay the convergence. The delay is significant even with a very small fault rate (Figure 2a). When the fault rate is relatively large (Figure 2d), this strategy is likely to stagnate completely. It does exhibit large convergence peak even when a little amount of data is lost as it can be observed in Figure 3 when the amount of lost data ranges from 0.2% to 6%. Contrary to Reset, both LI and LSI are extremely robust and resilient. Indeed, regardless the number of faults and the volume of lost data, LI and LSI almost consistently overlap with ER and NF, except in the presence of a very large fault rate (Figure 2d). However, the resilience capability of LI and LSI is preserved because they overlap with ER when varying either the fault rate (Figure 2) or the volume of lost data (Figure 3).

When converging the Ritz pair associated with the largest eigenvalues, NF converges in 485 iterations as depicted in Figures 4 and5. The Reset strategy again exhibits large peaks in the scaled residual norm after each fault, but this time, it can converge when only a few faults occur (Figure 4a) or only a small amount of data is lost (Figure 5). Regarding the robustness of both IR strategies, the convergence histories of LI and LSI almost consistently overlap again the NF curve regardless the fault rate (Figure 4) or the volume of lost data (Figure 5).

Arnoldi method to compute one eigenpair

In this section, we assess the robustness of our resilient Arnoldi for computing the eigenpair associated with the largest eigenvalue in magnitude. Because this is an easy calculation for Arnoldi, we select a rather small restart parameter m = 7 so that the convergence is not too fast in order to leave room to inject a few faults. One iteration consists of building a Krylov basis of size m followed by the approximation of the desired eigenpair. Because this computation requires only a few (outer) iterations we consider one single fault rate and we focus on the impact of proportion of lost data. Because the fault rate is constant the number of faults displayed in the different plots might differ depending on the convergence penalty they induce for the different resilient strategies. We report the convergence histories in Figure 6. NF converges smoothly in 9 (outer) iterations whereas Reset strategy exhibits a large peak in the scaled residual norm after each fault. When the amount of lost data is low (0.2% in Figure 6a), the Reset penalty remains reasonable (2 extra iterations), but it becomes more significant (7 extra iterations) if that amount increases (6% in Figure 6d). Because ER has to restart more often than NF, its convergence history exhibits some whereas the fault rate is constant (7 faults). LI, LSI, ER and NF coincide.

delay compared to NF. On the other hand, both IR strategies are again extremely robust. Indeed, LI and LSI convergence coincide with ER, regardless the proportion of lost data (Figure 6). Note that if the proportion of lost data is very large (Figure 6d), LI and LSI may slightly differ from ER. The fact that LI and LSI convergence coincide with ER, indicates that the spectral information regenerated by the LI and LSI is as good as the one computed by the regular solver.

Implicitly restarted Arnoldi method to compute nev eigenpairs

To investigate the robustness of IR strategies designed for IRAM, we compute the five eigenpairs (nev = 5) that correspond to the largest magnitude eigenvalues. At each iteration, we report the maximum of the scaled residual norms of those five sought eigenpairs. We consider a restart parameter m=10 (see Algorithm 2). One iteration thus consists of building a Krylov subspace of size 10, followed by the computation of the approximate eigenpairs. If the eigenpairs do not satisfy the stopping criterion, the next iteration starts with a contracted Arnoldi equality of size m = 5.

The NF calculation computes the five sought eigenvectors in 11 (outer) iterations (index k in Algorithm 2). The Reset strategy exhibits a large peak in the scaled residual norm after each fault, its scaled residual norm increases further than the initial one. As a consequence, convergence is very much delayed. Furthermore, Reset is also sensitive to the amount of lost data. The larger the volume of lost data, the more its convergence is delayed (Figure 7). On the other hand, both IR strategies are much more robust than Reset. However, they still require a few extra iterations than NF. Because the high coincidence with ER, it can be concluded that this slight penalty is not due to the quality of interpolation but to the necessity of restarting with the information of the five dimension space compressed in one single direction.

Jacobi-Davidson method to compute nev eigenpairs

In this section, we investigate the resilience of the IR strategies designed for Jacobi-Davidson. In all the experiments, we seek for the five (nev = 5) eigenpairs whose eigenvalues are the closest to zero (τ = 0). The correction equation is solved using 30 iterations of GMRES. To facilitate the readability and the analysis of the convergence histories plotted in this section, we use vertical green lines to indicate the convergence of new eigenpairs (such as iterations 95, 130, 165, 200 and 242 in Figure 8a), and vertical red lines to indicate faulty iterations (such as iterations 148 and 228 for the sixth and ninth fault, respectively, in Figure 8a). According to Remark 1, although very likely to happen, there is no guarantee of retrieving all the already converged Schur vectors in the basis used to restart. As a consequence, we indicate the number of Schur vectors retrieved in the basis used to restart in red color under the vertical red line corresponding to the associated fault. For instance, 2 already converged Schur vectors are immediately retrieved at restart, after the fault at iteration 148 in Figure 8a. In the Jacobi-Davidson method there is some flexibility for selecting the number of vectors (i.e., the dimension of the space generated for restarting) that can be interpolated after a fault. For our experiments we choose to interpolate the converged Schur vectors as well as nev of the best candidates for Schur vectors extracted from the search space V f . For the calculation of the five smallest eigenvalues, the NF algorithm converges in 210 iterations while faulty executions have extra iterations. For the sake of comparison, we consider only the first 300 iterations of all the runs so that the graphs have exactly the same scales and range of iteration count.

We first consider all the different restarting strategies. We report in Figure 8 their convergence histories in different subplots for a fixed fault rate to evaluate the quality of the basis used for the restarts. The curves in Figure 8a show the impact of the enforced restarts (35 additional iterations for ER compared to NF) and will serve as reference to evaluate the quality and relevance of the interpolated directions considered for LI, LSI and Reset. The first comment goes to the Reset approach that completely fails and is unable to compute any eigenpair. The LI interpolation behaves much better but only succeeds to compute four out of the five eigenpairs in 300 iterations. For this example, LSI is slightly more robust and computes the five eigenpairs with a few extra iterations compared to ER. For both IR approaches it can be observed that the converged Schur Finally we notice that even if, in Figure 8, experiments are performed with the same fault rate, the number of faults varies depending on the interpolation strategy. Indeed, the less robust the strategy, the more the convergence is delayed, which increases the probability for additional faults to occur. In order to reduce the number of graphs to illustrate the impact of the amount of lost data (see Figure 9) and the influence of the fault rate (see Figure 10), we only consider the LSI approach in the rest of this section. The influence of the volume of lost data is displayed in Figure 9 where its proportion varies from 0.2 % to 6 % while the fault rate remains constant. As one would expect, the general trend is: the larger the amount of lost data, the larger the convergence penalty. A large amount of lost data at each fault prevents LSI from converging the five eigenpairs; only four are eventually computed within 300 iterations for 6 % lost. Notice that allowing for more than the 300 iterations would probably enable the solver to succeed in the calculation of the five eigenpairs. Figure 9: Impact of the amount of lost data on the resilience of LSI using (nev + nconv) regenerated vectors when converging the five eigenpairs associated with the smallest eigenvalues using Jacobi-Davidson. The volume of lost data varies from 0.2% to 6% whereas the fault rate is constant. All converged Schur vectors are found immediately after interpolation followed by restart.

Finally, in Figure 10 we depict the convergence histories of the various restarting strategies when the fault rate varies leading to a number of faults that varies from 3 to 24; as expected, the larger the number of faults, the larger the convergence delay. However, the IR policy is rather robust and succeeds in converging the five eigenpairs in less than 300 iterations. Despite the robustness of the resilient schemes, peaks of the residual norm associated with the current Schur vector are often observed after each fault. A possible remedy consists of using a hybrid approach where we interpolate the nconv Schur vector while reusing the best candidate Schur vector available in V f (as if we had checkpointed this single direction) when the fault occurs and interpolate nev-1 directions but the first to recover additional meaningful spectral information from V f . This procedure is beneficial for improving the convergence independently of the size of the space used for the restart. It is best highlighted in the most difficult situations that are large amount of lost data (see Figure 11) and large number of faults (see Figure 12). In these figures, the convergence of the scaled residual norm no longer exhibits peak after the faults when the best Schur candidate vector is used in the set of directions for the restart. The most relevant information on the next Schur vector to converge seems to be concentrated in the current best Schur candidate vector. Intensive experiments show that the interpolation strategies do not succeed in regenerating accurately enough the lost entries of this special direction of V f . Consequently, a reasonable tradeoff may consist in checkpointing only the current Schur vector in order to increase the robustness Inria of the resilient Jacobi-Davidson method.

Concluding remarks

Many scientific and engineering applications require the computation of eigenpairs of large sparse matrices. The objective of the paper has been to study numerical schemes suitable for the design of resilient parallel eigensolvers. For that purpose, we have proposed two interpolation procedures to regenerate meaningful spectral information for restarting the eigensolver after a fault. To evaluate the qualitative behavior of the resilient schemes, we have simulated stressful conditions by increasing the fault rate and the volume of lost data.

We have considered two variants of subspace iteration methods. On the one hand, we have considered the subspace iteration with Chebyshev polynomial acceleration for the computation of eigenpairs corresponding to the smallest eigenvalues. On the other hand, we have considered the classical method for the computation of eigenpairs corresponding to the largest magnitude eigenvalues. For both methods, the Reset strategy strongly penalizes the convergence at each fault, while both LI and LSI are extremely robust and resilient, regardless the number of faults and the volume of lost data. The same numerical behavior is observed for our resilient Arnoldi for computing the eigenpair associated with the largest eigenvalue in magnitude. Our LI/LSI resilient IRAM for the computation of a few eigenpairs are much more robust than Reset. However, they do exhibit a slight penalty, not due to the quality of interpolation, but to the restarting policy that Inria leads to compressing in a single direction the eigenspace under calculation when faults occur.

We have had a stronger emphasis on the Jacobi-Davidson method. The motivation is twofold: the Jacobi-Davidson method is widely used in many real-life applications, and, in addition, it offers some flexibility to select different spectral information to construct an efficient restart mechanism after a fault. We have observed that despite the increase of the amount of recovered data, the peak of the residual norm associated with the current Schur vector persists after a fault. For a possible remedy of these effects, we have designed a hybrid approach that consists in combining interpolation techniques with classical checkpoint for a single vector. This illustrates that numerical resilient strategies can be effectively combined with the state-of-the-art fault tolerant policies to design efficient and robust methods as demonstrated for the Jacobi-Davidson algorithm by the combination of light checkpointing and a numerical resilience.

Finally, the effectiveness of the proposed approaches has been illustrated from a numerical view point. Their effectiveness and scalable implementation in a parallel computing environment deserves to be studied once the fault tolerant supports have been made available in the MPI standard.

m 14 :

 14 end for Interpolation-restart policy: When a fault occurs during iteration f + 1, it may be during the expansion of the Krylov subspace (Algorithm 2, line 13) or during the contraction step (Algorithm 2, line 3 to 12).

Algorithm 3 3 :

 33 [λ, u] = Basic-Jacobi-Davidson(v, τ ) Jacobi-Davidson algorithm to compute the eigenvalue of A closest to a target value τ 1: Set V 1 = [v] 2: for k = 1, 2, . . . ,until convergence do Compute Rayleigh quotient: C k = V H k AV k , and eigenpairs of C k 4:

Algorithm 4 1 :

 41 [µ, z, t]=JDQR(Z(:,1:nconv),V (:,1:k),W (:,1:k),C(1:k,1:k), τ , m,m,maxiter) Jacobi-Davidson style QR algorithm for expansion of partial Schur decomposition Set iter = 0; k init = k; tr = 0 2: while iter<maxiter do 3: iter = iter + 1 4:

14 : return 15 :

 1415 else if k = m then 16: % Restart: V (:, 1 : m) = V (:, 1 : m)Q(:, 1 : m); W (:, 1 : m) = W (:, 1 : m)Q(:, 1 : m); C(1 : m, 1 : m) = T (1 : m, 1 : m); k init = m

Remark 1 .

 1 Let us assume that the partial Schur decomposition has converged in exact arithmetic (AZ nconv = Z nconv T nconv ), and that the nconv eigenpairs are exact solutions (Au = λ u ) still in exact arithmetic. Under this assumption, the eigenvectors (u(IR) ) computed by IR are the same exact eigenvectors as long as A(I p , I p )-λ I Ip,Ip is nonsingular or A :,Ip -λ I :,Ip is full column rank for LI and LSI, respectively. As a consequence, if Jacobi-Davidson is restarted with the initial basis V nconv obtained from the orthonormalization of the vectors of U (IR)

Figure 2 :

 2 Figure2: Impact of the fault rate on the resilience of IR strategies when converging the five eigenpairs associated with the smallest eigenvalues using subspace iteration method accelerated with Chebyshev. A proportion of 0.8 % of data is lost at each fault. Note that LI, LSI, ER and NF almost coincide.

Figure 3 :

 3 Figure3: Impact of the amount of lost data on the resilience of IR strategies when converging the five eigenpairs associated with the smallest eigenvalues using subspace iteration method accelerated with Chebyshev. The volume of lost data varies from 0.2% to 6% whereas the fault rate is constant (7 faults). LI, LSI, ER and NF coincide.

Figure 4 :

 4 Figure 4: Impact of the fault rate on the resilience of IR strategies when converging the five eigenpairs associated with the largest eigenvalues using the basic subspace iteration method. A proportion of 0.8 % of data is lost at each fault. LI, LSI and NF coincide in (a), (b) and (c).

Figure 5 :

 5 Figure5: Impact of the amount of lost data on the resilience of IR strategies when converging the five eigenpairs associated with the largest eigenvalues using basic subspace iteration method. The volume of lost data varies from 0.2% to 6% whereas the fault rate is constant (9 faults). LI, LSI and NF coincide.

  6% of lost data(5 faults)

Figure 6 :

 6 Figure 6: Impact of the amount of lost data on the resilience of IR strategies when converging the eigenpair associated with the largest eigenvalue using Arnoldi method. LI, LSI and NF coincide in (a) (b) and (c).

Figure 7 :

 7 Figure 7: Impact of the amount of lost data on the resilience of IR strategies when converging the five eigenpairs with largest eigenvalues using IRAM. The volume of lost data varies from 0.2% to 6% whereas the fault rate is constant. LI, LSI and ER coincide in (a) and (b).

Figure 8 :

 8 Figure8: Comparison of IR strategies using nev regenerated vectors when converging the five eigenpairs associated with the smallest eigenvalues using Jacobi-Davidson. The fault rate is the same for all strategies and a proportion of 0.8 % data is lost at each fault.

Figure 10 :

 10 Figure10: Impact of the fault rate on the resilience of LSI using (nev + nconv) regenerated vectors when converging the five eigenpairs associated with the smallest eigenvalues using Jacobi-Davidson. The fault rate varies whereas a proportion of 0.2 % of data is lost at each fault. At each fault, all the already converged Schur vectors are retrieved in the basis of restart.

Figure 11 :

 11 Figure 11: Impact of keeping the best Schur vector candidate in the search space after a fault combined with nev + nconv -1 interpolated directions. Top two plots nev interpolated directions, bottom two plots nev -1 interpolated directions plus the best Schur candidate. Calculation of the five eigenpairs associated with the smallest eigenvalues using Jacobi-Davidson. The fault rate is constant over all sub-figures. The proportion of lost data is either 3% in Figures (a) and (c) or 6% in Figures (b) and (d). All converged Schur vectors are found immediately after interpolation followed by restart.

Figure 12 :

 12 Figure12: Impact of keeping the best Schur vector candidate in the search space after a fault combined with nev + nconv -1 interpolated directions. Top two plots nev + nconv interpolated directions, bottom two plots nev + nconv -1 interpolated directions plus the best Schur candidate. Calculation of the five eigenpairs associated with the smallest eigenvalues using Jacobi-Davidson with a proportion of 0.2 % of data is lost at each fault.

RR n°8625

Inria

avenue de la Vielle Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

This work was partially supported by the French research agency ANR in the framework of the RESCUE project (ANR-10-BLANC-0301), in particular the PhD thesis of the fourth author was funded by this project. This research also benefited from the G8-ECS project.