Bernhard Beckermann

Daniel Kressner
email: daniel.kressner@epfl.ch

AND Marcel Schweitzer
email: marcel.schweitzer@epfl.ch.

LOW-RANK UPDATES OF MATRIX FUNCTIONS

Keywords: matrix function, low-rank update, Krylov subspace method, tensorized Krylov subspace, matrix exponential, Markov function, graph communicability AMS subject classifications. 15A16, 65D30, 65F30, 65F60

We consider the task of updating a matrix function f (A) when the matrix A ∈ C n×n is subject to a low-rank modification. In other words, we aim at approximating f (A + D) -f (A) for a matrix D of rank k n. The approach proposed in this paper attains efficiency by projecting onto tensorized Krylov subspaces produced by matrix-vector multiplications with A and A * . We prove the approximations obtained from m steps of the proposed methods are exact if f is a polynomial of degree at most m and use this as a basis for proving a variety of convergence results, in particular for the matrix exponential and for Markov functions. We illustrate the performance of our method by considering various examples from network analysis, where our approach can be used to cheaply update centrality and communicability measures.

1. Introduction. This work is concerned with the problem of updating a matrix function f (A) when A ∈ C n×n is subject to a low-rank modification. More specifically, given D ∈ C n×n of rank r n we aim at computing the difference

f (A + D) -f (A), (1.1)
at a cost significantly smaller than the cost of computing f (A + D) from scratch. Throughout this work, we assume that f : Ω → C is analytic on some domain Ω containing Λ(A) and Λ(A + D), the spectra of A and A + D. The generality of (1.1) includes a wide scope of applications. We will particularly focus on measuring network properties, such as centrality and communicability, via applying the matrix exponential to the adjacency matrix of an undirected graph [START_REF] Estrada | Network properties revealed through matrix functions[END_REF]. Removing/inserting individual edges or nodes in the graph corresponds to rank-2 updates of A. Being able to solve (1.1) efficiently therefore allows to quickly update these measures. In fact, this task is explicitly needed when measuring the betweenness of a node [START_REF] Estrada | Network properties revealed through matrix functions[END_REF]Sec. 2.3].

For f (z) = z -1 , the well-known Sherman-Morrison-Woodbury formula implies that the difference (1.1) has rank r and can be directly obtained from applying f (A) = A -1 to the low-rank factors of D. In particular, when D has rank 1 and can be written as D = bc * with vectors b and c, we have

(A + bc * ) -1 -A -1 = - 1 1 + c * A -1 b A -1 bc * A -1 . (1.2)
Unless A is a scalar multiple of the identity matrix [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]Theorem 1.35], there is no such simple relation between f (A) and f (A + bc * ) for a general analytic function f . In the case of a rational function f ≡ r of degree δ, Bernstein and Van Loan [START_REF] Bernstein | Rational matrix functions and rank-1 updates[END_REF] show that r(A + bc * ) -r(A) has rank at most δ and provide an explicit formula for the rank-δ correction. This construction is based on explicit Krylov bases and potentially prone to numerical instability for larger degrees. The Cauchy integral representation

f (A) = 1 2πi Γ f (z)(zI -A) -1 dz (1.3)
provides a link between matrix functions and (shifted) inverses. In [START_REF] Stange | Beiträge zu effizienten Algorithmen basierend auf Rang-1 Aufdatierungen[END_REF], a combination of numerical quadrature with (1.2) has been explored for approximating exp(A + bc * ) -exp(A). However, such an approach suffers from a number of drawbacks. Most importantly, the choice of a contour Γ that is equally good for A and A + bc * is highly nontrivial, in particular in the non-Hermitian case. Also, the evaluation of the approximation amounts to the solution of a differently shifted linear system for each quadrature point. Although not used in our algorithms, integral representations will play a role in their convergence analysis, see Section 5.

The approach proposed in this work avoids the need for choosing a contour or an explicit rational approximation of f . Inspired by existing Krylov subspace methods for matrix equations [START_REF] Simoncini | Computational methods for linear matrix equations[END_REF] and linear systems with tensor product structure [START_REF] Kressner | Krylov subspace methods for linear systems with tensor product structure[END_REF], we make use of tensorized subspaces. More specifically, given orthonormal bases U m , V m for m-dimensional Krylov subspaces involving the matrices A, A * and the starting vectors b, c, we construct an approximation of the form

f (A + bc * ) -f (A) ≈ U m X m (f )V * m ,
with a well-chosen small matrix X m (f ) ∈ C m×m . In turn, when m n, the computational cost is dominated by m matrix-vector products with A and A * . As we will demonstrate for a variety of examples, this dramatically reduces the computational effort compared to computing f (A +bc * ) from scratch, even when only selected quantities -such as the diagonal or trace of f (A + bc * ) -are required.

The rest of this paper is organized as follows. In section 2, we introduce Krylov subspace algorithms for approximating the update (1.1) when D is a matrix of rank one, both for the Hermitian and the non-Hermitian case, and discuss the extension to higher rank. In section 3, we prove that the approximation obtained from m steps of our Krylov subspace algorithm is exact when f is a polynomial of degree at most m. This result forms the basis of the convergence analysis presented in sections 4 and 5. Precisely, the convergence of our Krylov approximations is linked to certain polynomial approximation problems in section 4, while section 5 contains results that are based on exploiting an integral representation of f (A), e.g., for general analytic functions and for Markov functions. Applications for our methods are described in Section 6, with special focus on up-/down-dating of communicability measures in network analysis, and the efficiency of our approach is illustrated by numerical experiments from this area. Concluding remarks are given in Section 7.

We use • to denote the Euclidean norm of a vector or the induced spectral norm of a matrix.

Krylov projection algorithms.

In the following, we assume that f (A) or at least the part relevant to the application (e.g., the diagonal of f (A)) have already been computed. For the moment, we continue to assume a rank-1 modification, that is, we consider the approximation of f (A+bc * )-f (A) for vectors b, c. In Remark 2.3 below, we will comment on the extension to higher ranks. The first step of our algorithm consists of constructing an orthonormal basis for a Krylov subspace of A with starting vector b. Given m ≤ n, such a Krylov subspace takes the form

U m := K m (A, b) = span{b, Ab, A 2 b, . . . , A m-1 b}.
For simplicity, we suppose that U m has dimension m, which is generically satisfied. Applying m steps of the Lanczos method [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear differential and integral operators[END_REF] (possibly with reorthogonalization) results in the Lanczos relation

AU m = U m G m + β m+1 u m+1 e * m , (2.1)
where the columns of U m ∈ C n×m form an orthonormal basis of U m , the matrix

G m = U * m AU m ∈ C m×m is tridiagonal
, and e m denotes the mth unit vector of length m. The first column of U m is a scalar multiple of b and, without loss of generality, we may assume that U * m b = b e 1 with the first unit vector e 1 ∈ C m . The second step of our algorithm then chooses an approximation in the tensorized subspace U m ⊗ U m , that is,

f (A + bb * ) -f (A) ≈ U m X m (f )U * m (2.2)
for some matrix X m (f ) ∈ C m×m . In the spirit of Krylov subspace methods for matrix equations [START_REF] Simoncini | Computational methods for linear matrix equations[END_REF], we choose X m (f ) as the solution of the compressed problem:

X m (f ) = f U * m (A + bb * )U m -f U * m AU m = f G m + b 2 e 1 e * 1 -f (G m ), (2.3)
where we assume that f is defined on the spectra of G m and G m + b 2 e 1 e * 1 . Below, in Theorem 3.2, we will see that this choice of X m (f ) leads to an approximation that is exact when f is a polynomial of degree at most m.

The resulting method is summarized in Algorithm 1. The tridiagonal matrix G m from (2.1) is built from the orthogonalization coefficients α j , β j as

G m =        α 1 β 2 β 2 α 2 β 3 . . . . . . . . . β m-1 α m-1 β m β m α m        .
A trivial modification of this algorithm can be used to approximate

f (A -bb * ) - f (A).
Algorithm 1 represents the most basic form of the Lanczos process; in particular, it does not employ any reorthogonalization. It is well known that such short recurrences may suffer from severe loss of orthogonality in the presence of round-off errors, so that it can be advisable to use reorthogonalization strategies [START_REF] Paige | The computation of eigenvalues and eigenvectors of very large sparse matrices[END_REF][START_REF] Simon | Analysis of the symmetric Lanczos algorithm with reorthogonalization methods[END_REF]. The most straightforward to retain numerical orthogonality amongst the basis vectors is to store all basis vectors and perform full reorthogonalization in each step of the method.

An alternative is to perform no reorthogonalization at all. In contrast to, e.g., the CG method for linear systems, there are no short recurrences for the iterates U m X m (f )U * m available. In turn, this requires to use a two-pass strategy for forming

Algorithm 1 Krylov subspace approximation of f (A + bb * ) -f (A) for Hermitian A 1: u 0 = 0. 2: u 1 = (1/ b )b. 3: β 1 = 0.
4: for j = 1, . . . , m do 5:

w j = Au j -β j u j-1 .

6:

α j = u * j w j .

7:

w j = w j -α j u j .

8:

β j+1 = w j .

9: U m X m (f )U * m if one wants to avoid storing the full basis U m . In a first pass, the tridiagonal matrix G m is assembled (while storing only three basis vectors at a time and discarding the older ones) and X m (f ) is formed. In a second pass, the basis vectors are computed anew (again only storing three at a time), and the parts of interest of U m X m (f )U * m (e.g., its diagonal) are gradually computed. This approach doubles the number of matrix-vector products but reduces the storage requirement from O(mn) to O(m 2 + n).

u j+1 = (1/β j+1 )w j . 10: end for 11: Compute matrix function X m (f ) = f G m + b 2 e 1 e * 1 -f (G m ) 12: Return U m X m (f )U * m .
Example 2.1.

We illustrate the (possible) impact of reorthogonalization on convergence. Consider the diagonal matrices A 1 , A 2 ∈ C 100×100 , where the eigenvalues of A 1 are equidistantly spaced in the interval [10 -3 , 10 3 ] and the eigenvalues of A 2 are logarithmically spaced in [10 -3 , 10 3 ]. We approximate f (A i + bb * ) -f (A i ), i = 1, 2 where f (z) = exp(-z) and b is a random vector of unit norm using Lanczos without reorthogonalization and with full reorthogonalization, as explained above. The resulting convergence curves are shown in Figure 2.1. For A 1 , both methods are observed to behave identically, which is well explained by the fact that the eigenvalue distribution of A 1 does not favor the convergence of Ritz values, which has a close link to loss of orthogonality [START_REF] Meurant | The Lanczos and conjugate gradient algorithms in finite precision arithmetic[END_REF][START_REF] Paige | The computation of eigenvalues and eigenvectors of very large sparse matrices[END_REF]. For A 2 , the eigenvalue distribution favors the convergence of Ritz values to larger eigenvalues and, in turn, convergence degrades after some time when using no reorthogonalization. Still, the method eventually converges to the same accuracy, it just needs more iterations. This phenomenon is also known from the conjugate gradient method in finite precision arithmetic [START_REF] Meurant | The Lanczos and conjugate gradient algorithms in finite precision arithmetic[END_REF] and when approximating f (A)b [START_REF] Druskin | Error bounds in the simple Lanczos procedure for computing functions of symmetric matrices and eigenvalues[END_REF]. To avoid that orthogonality issues distort our findings, we use the Lanczos method with full reorthogonalization for all numerical experiments in the rest of this paper.

The computational cost of m steps of Algorithm 1 is as follows:

• m matrix-vector products with A (2m, when the two-pass approach is used), requiring O(m • nnz(A)) operations for a sparse matrix A with nnz(A) nonzeros; • O(m 2 n) operations for the orthogonalization procedure when using full reorthogonalization (and only O(mn) when the two-pass approach is used, because in each step orthogonalization against only two previous basis vectors is performed) • the computation of X m (f ), which requires the evaluation of two functions of m × m matrices, which, depending on the function f , typically needs at most O(m 3 ) operations. One should avoid forming U m X m (f )U * m explicitly and we expect that this is actually not needed in most applications involving a large sparse matrix A. For example, the computation of communicability measures discussed in Section 6.1 only requires the diagonal entries of U m X m (f )U * m , which can be computed directly from U m , X m (f ) with O(m 2 n) operations.

2.2. Non-Hermitian rank-1 updates. In the general non-Hermitian case, we construct orthonormal bases for the two (polynomial) Krylov subspaces U m := K m (A, b) and V m := K m (A * , c). Applying the Arnoldi method with reorthogonalization results in the Arnoldi relations (2.1) and, additionally,

A * V m = V m H m + h m+1,m v m+1 e * m , (2.4)
where the columns of V m ∈ C n×m form an orthonormal basis of V m . Note that both G m = U * m AU m and H m = V * m A * V m are now m × m upper Hessenberg matrices. The approximation is chosen in the tensorized subspace U m ⊗ V m , that is,

f (A + bc * ) -f (A) ≈ U m X m (f )V * m (2.5)
for some matrix X m (f ) ∈ C m×m . Concerning the choice of X m (f ), it turns out that the non-Hermitian analogue of (2.3) does not have favorable theoretical properties. For example, the polynomial exactness property mentioned above for (2.3) and proven in section 3 does not hold for such a choice. We will use a different choice, motivated by the following simple result.

Lemma 2.2. Let A ∈ C n×n , and b, c ∈ C n . Define the block matrix

A := A bc * 0 A + bc * .
(2.6)

Then

f (A) = f (A) f (A + bc * ) -f (A) 0 f (A + bc * ) .
Proof. Letting Γ denote a contour that encloses both Λ(A) and Λ(A + bc * ), the Cauchy integral formula (1.3) applied to f (A) yields

f (A) = 1 2πi Γ f (z)(zI -A) -1 dz = 1 2πi Γ f (z) (zI -A) -1 -(zI -A) -1 bc * (zI -A -bc * ) -1 0 (zI -A -bc * ) -1 dz = f (A) -1 2πi Γ f (z)(zI -A) -1 bc * (zI -A -bc * ) -1 dz 0 f (A + bc * ) . (2.7)
On the other hand, combining the Cauchy integral formula for f (A + bc * ) and f (A) with the second resolvent identity yields

f (A + bc * ) -f (A) = 1 2πi Γ f (z) (zI -A -bc * ) -1 -(zI -A) -1 dz = - 1 2πi Γ f (z)(zI -A -bc * ) -1 bc * (zI -A) -1 dz = - 1 2πi Γ f (z)(zI -A) -1 bc * (zI -A -bc * ) -1 dz, (2.8)
which matches the (1,2) block of (2.7) and thus completes the proof. The result of Lemma 2.2 shows that the desired update is contained in the (1,2) block of f (A). This motivates us to choose X m (f ) as the (1,2) block of f applied to the compression of

A onto U m ⊕ V m : f U m 0 0 V m * A bc * 0 A + bc * U m 0 0 V m = f G m b c e 1 e * 1 0 H * m + c V * m be * 1 = f (G m ) X m (f ) 0 f (H * m + c V * m be * 1 ) , (2.9)
where we again assume that f is defined on the spectra of G m and H * m + c V * m be * 1 . Using Lemma 2.2, it is straightforward to see that this choice of X m (f ) coincides in the Hermitian case with the one from section 2.1.

Algorithm 2 summarizes the proposed procedure, where we omit the algorithmic details for the Arnoldi method for the sake of brevity; see, e.g., [START_REF] Golub | Matrix Computations[END_REF].

Algorithm 2 Krylov subspace approximation of f (A + bc * ) -f (A) 1: Perform m steps of the Arnoldi method to compute an orthonormal basis U m of K m (A, b) and G m = U * m AU m . 2: Perform m steps of the Arnoldi method to compute an orthonormal basis V m of K m (A * , c) and H m = V * m A * V m . 3: Compute matrix function F m = f G m b c e 1 e * 1 0 H * m + c V * m be * 1 . 4: Set X m (f ) = F m (1 : m, m + 1 : 2m). 5: Return U m X m (f )V * m .
The computational effort of Algorithm 2 compares to that of Algorithm 1 (with full reorthogonalization) as follows. In contrast to the Hermitian case, we now need to build two Krylov spaces instead of one, meaning that the number of matrix vector products necessary for m steps of the method increases from m to 2m. The orthogonalization cost is the same as in the Hermitian case (with full reorthogonalization). As the number of operations needed for approximating a (dense) matrix function often grows cubically in the matrix size it is typically about four times as expensive to compute one matrix function of size 2m × 2m than computing two matrix functions of size m×m. The fact that the matrix for which we need to evaluate f in Algorithm 2 is non-Hermitian and not tridiagonal may increase the cost further. However, as long as m n, the main cost of Algorithm 2 consists of performing matrix-vector products, so that we can expect it to take roughly twice the computation time of Algorithm 1 for a problem of the same size and structure. i , e.g., by computing a singular value decomposition of D and then apply k times Algorithm 2 in order to subsequently incorporate each of the k rank-one modifications. Note that the ith step of this procedure requires working with the Krylov subspaces

K m (A + i-1 j=1 b j c * j , b i ) and K m (A * + i-1 j=1 c j b * j , c i ) which in general do not coincide with K m (A, b i ) and K m (A * , c i ).
Some care is required in order to make use of Algorithm 1 for a Hermitian rankk matrix D. After a preprocessing step (see, e.g., [6, Section 2.3]) one can write

D = k i=1 b i b * i - k i= k+1 b i b * i .
First, Algorithm 1 is applied to incorporate the first k terms followed by a slight modification of this algorithm to incorporate the last k -k terms 2. Block Krylov subspaces [START_REF] Gutknecht | Block Krylov space methods for linear systems with multiple right-hand sides: An introduction[END_REF][START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF][START_REF] Simoncini | An iterative method for nonsymmetric systems with multiple right-hand sides[END_REF] offer a conceptually different way of dealing with a rank-k modification. Instead of the Arnoldi method, a block Arnoldi method is used for computing orthonormal bases of the block Krylov subspaces K m (A, B) = range[B, AB, . . . , A m-1 B] and K m (A * , C). The approximation of f (A + D) -f (A) is computed by projecting onto the tensorized block Krylov subspace K m (A, B) ⊗ K m (A * , C), leading to a straightforward extension of Algorithm 2. Block Krylov subspace methods are more complicated to implement, see, e.g., [START_REF] Gutknecht | Block Krylov space methods for linear systems with multiple right-hand sides: An introduction[END_REF] for some of the issues one has to take into account, but they offer (at least) one major advantage. Even though the product of A with an n × k matrix requires the same number of operations as k individual matrix vector products, it often performs much faster on a computer, benefitting from a more "cache-friendly" memory access pattern, see, e.g., [START_REF] Baker | On improving linear solver performance: A block variant of GMRES[END_REF]. Again some care is required in the symmetric case, to derive a block variant of Algorithm 1.

Stopping criterion. A simple stopping criterion for Algorithm 1 or Algorithm 2 is based on the error estimate

f (A + bc * ) -U m X m (f )V * m ≈ U m+d X m+d (f )V * m+d -U m X m (f )V * m (2.10)
for some small integer d ≥ 1. This error estimate is similar in spirit to error estimates for linear systems proposed in [START_REF] Golub | Matrices, moments and quadrature II; How to compute the norm of the error in iterative methods[END_REF][START_REF]The computation of bounds for the norm of the error in the conjugate gradient algorithm[END_REF]. Evaluating the right-hand side of (2.10) only requires forming the coefficient matrices X m (f ), X m+d (f ) defined in (2.2) or (2.5), because

U m+d X m+d (f )V * m+d -U m X m (f )V * m = U m+d X m+d (f )V * m+d -U m+d X m (f ) 0 0 0 V * m+d = X m+d (f ) - X m (f ) 0 0 0 .
We illustrate the behavior of the error estimator (2.10) by means of a simple numerical experiment.

Example 2.4. We choose A ∈ R 400×400 as the standard finite difference discretization of the two-dimensional Laplace operator on the unit grid, and b, c as vectors with random, normal distributed entries. We use our proposed Krylov subspace algorithm to approximate (A + bc * ) -1/2 -A -1/2 and show the convergence history as well as the proposed error estimator for d = 1, 2, 3 in Figure 2.2. We observe that the difference-based error estimators follow the norm of the exact error very closely and thus give very accurate results already for such small values of d.

We remark that the error estimator (2.10) is clearly heuristic and can be overly optimistic, especially in situations where the iteration (almost) stagnates, although it works very well in most practical situations. Of course, when already m + d iterations of the method have been performed, one will typically use the approximation U m+d X m+d V * m+d from the (m + d)th step instead of the approximation from the mth step (for which the error is estimated), as it can be expected to be a more accurate approximation. In fact, in Example 2.4, for d = 3 the error of the iterate X m+d lies below the error estimate (2.10) for all but four iterations.

3. Exactness properties of Krylov subspace approximations for the update. In this section, we show that the approximation (2.5) is exact when f is a polynomial of degree at most m. This serves two purposes: On the one hand, this justifies the particular choice of approximation we made. On the other hand, this provides the fundamental basis for performing our convergence analyses in section 4.

As a tool for proving polynomial exactness, we utilize the following matrix identity, which does not seem to be widely known.

Proposition 3.1. Let M, N ∈ C n×n . Then, for any j ∈ N,

M j -N j = j-1 k=0 N j-1-k (M -N )M k .
Proof. The proof is by induction on j. For j = 0, the assertion trivially holds. Suppose now that the assertion holds for j -1:

M j-1 -N j-1 = j-2 k=0 N j-2-k (M -N )M k . (3.1)
Multiplying both sides of (3.1) by N from the left and adding M j , we find

M j + N M j-1 -N j = M j + j-2 k=0 N j-1-k (M -N )M k .
Finally, subtracting M j-1 N on both sides and noting that

M j -N M j-1 = N 0 (M - N )M j-1 completes the proof. Theorem 3.2. Let A ∈ C n×n , b, c ∈ C n .
Then the Krylov subspace approximation returned by Algorithm 2 is exact for all p ∈ Π m , where Π m denotes the space of all polynomials of degree at most m, i.e.,

p(A + bc * ) -p(A) = U m X m (p)V * m .
Proof. By linearity, it suffices to show that the result holds for every monomial p j (z) = z j , j = 0, . . . , m. We recall that X m (p j ) is the (1, 2)-block of the matrix

G m b c e 1 e * 1 0 (H * m + c V * m be * 1 ) j .
In particular, X m (p 1 ) = b c e 1 e * 1 . By (2.9), the result of the theorem is shown if we can establish

(A + bc * ) j -A j = U m X m (p j )V * m . (3.2) Considering G m b c e 1 e * 1 0 (H * m + c V * m be * 1 ) j = G j-1 m X m (p j-1 ) 0 (H * m + c V * m be * 1 ) j-1 G m X m (p 1 ) 0 (H * m + c V * m be * 1 )
, we find the relation

X m (p j ) = G j-1 m X m (p 1 ) + X m (p j-1 )(H * m + c V * m be * 1 ).
Resolving this recursion gives

X m (p j ) = j-1 k=0 G j-1-k m ( b c e 1 e * 1 )(H * m + c V * m be * 1 ) k . (3.3)
Recalling that G m is the compression of A onto K m (A, b), a well-known polynomial exactness property (see, e.g., [START_REF] Th | Computing functions of matrices using Krylov subspace methods[END_REF][START_REF] Saad | Analysis of some Krylov subspace approximations to the matrix exponential operator[END_REF]) of Krylov subspace approximations yields

b U m G m e 1 = A b for all = 0, . . . , m -1. (3.4)
Similarly, by noting that

H m + c e 1 bV * m is the compression of (A + bc * ) * onto K m (A * , c) = K m ((A + bc * ) * , c), we obtain c e * 1 (H * m + c V * m be * 1 ) V * m = c * (A +
U m X m (p j )V * m = j-1 k=0 U m G j-1-k m ( b c e 1 e * 1 )(H * m + c V * m be * 1 ) k V * m which by (3.4)-(3.5) gives for j ≤ m the relation U m X m (p j )V * m = j-1 k=0 A j-1-k bc * (A + bc * ) k = (A + bc * ) j -A j ,
where we used Proposition 3.1 in the second equality. This establishes (3.2) and thus completes the proof.

4. Convergence analysis based on polynomial approximation problems. In this section, we give bounds for the error of the approximations (2.2) and (2.5). The results are based on the polynomial exactness property from Theorem 3.2 and connect the approximation quality of (2.2) and (2.5) to certain polynomial approximation problems.

The convergence results in this section rely on a theorem by Crouzeix [START_REF] Crouzeix | Numerical range and functional calculus in Hilbert space[END_REF][START_REF] Crouzeix | The numerical range as a spectral set[END_REF], for which we first recall some basic concepts. The field of values (or numerical range) A ∈ C n×n is defined as the set

W (A) := {x * Ax : x = 1},
which is a convex and compact subset of C containing all eigenvalues of A. Further, we make use of the supremum norm

f E := sup z∈E |f (z)|
on a subset E ⊂ C for which f is defined. Using these definitions, Crouzeix's theorem states that

f (A) ≤ C f E (4.1)
with a constant C ≤ 1 + √ 2 for any function f which is analytic in a neighborhood of E ⊇ W (A). Notice that if A is Hermitian then we can compute f (A) in terms of the eigenvalues of A, and thus obviously (4.1) holds with C = 1.

We now use (4.1) together with the polynomial exactness property proven in Theorem 3.2 to obtain bounds on the error at the mth step of our method:

E m (f ) := f (A + bc * ) -f (A) -U m X m (f )V * m . (4.2)
Theorem 4.1. Let A be Hermitian and let f be defined on a compact set E containing W (A) ∪ W (A + bb * ). Then, the error

(4.2) with U m = V m and X m (f ) from (2.3) satisfies E m (f ) ≤ 4 min p∈Πm f -p E .
Proof. By the polynomial exactness property of U m X m (f )U * m stated in Theorem 3.2, we have

E m (f ) = E m (f ) -p(A + bb * ) + p(A) + U m X m (p)U * m = E m (f -p) for any polynomial p ∈ Π m . For arbitrary p ∈ Π m we thus have E m (f ) = E m (f -p) = (f -p)(A + bb * ) -(f -p)(A) -U m X m (f -p)U * m ≤ (f -p)(A + bb * ) -(f -p)(A) -U m X m (f -p)U * m ≤ 2 f -p E + U m X m (f -p)U * m , (4.3)
where the last inequality follows from (4.1). By (2.2) and (2.3), we have

U X m (f -p)U * m = U m (f -p) G m + b 2 e 1 e * 1 -(f -p)(G m ) U * m . In turn, U m X m (f -p)U * m = X m (f -p) ≤ 2 f -p E , (4.4)
where, using W (G m ) ⊆ W (A) and W (G m + b e 1 e * 1 ) ⊆ W (A + bb * ), we again applied (4.1). Inserting (4.4) into (4.3) and taking the minimum over all p ∈ Π m completes the proof.

Theorem 4.1 allows to derive convergence bounds from known polynomial approximation results for analytic functions. The obtained convergence rates will typically be exponential for functions which are analytic in a neighborhood of E and superlinear for entire functions like the exponential.

To extend Theorem 4.1 to the non-Hermitian case, we have to assume that f is analytic on the field of values of the block matrix A from (2.6).

Theorem 4.2. Let A :=

A bc * 0 A + bc * and let f be analytic in a neighborhood of a compact set E containing W (A). Then, the error (4.2), with U m , X m (f ), and V m computed by Algorithm 2, satisfies

E m (f ) ≤ 2C min p∈Πm f -p E with a constant C ≤ 1 + √ 2. Proof. By Lemma 2.2, the update f (A + bc * ) -f (A) is the (1,2) block of f (A), which we denote by f (A) 1,2 . Letting W m = U m 0 0 V m , we note that the columns of W m are orthonormal and U m X m (f )V * m = [W m f (W * m AW m )W * m ] 1,2 holds by the definition of X m (f ).
As in the proof of Theorem 4.1 we use the fact that E m (f ) = E m (f -p) for any p ∈ Π m . Thus, we obtain for arbitrary p ∈ Π m that

E m (f ) = E m (f -p) = (f -p)(A + bc * ) -(f -p)(A) -U m X m (f -p)V * m ≤ [(f -p)(A)] 1,2 + [W m (f -p)(W * m AW m )W * m ] 1,2 ≤ (f -p)(A) + (f -p)(W * m AW m ) ≤ 2C f -p E
where the last inequality follows from Crouzeix's theorem, using W (W * m AW m ) ⊆ W (A) ⊆ E. Taking the minimum over all p ∈ Π m gives the desired result.

Note that the matrix A from Theorem 4.2 can be easily block-diagonalized: Unfortunately, this does not seem to admit any meaningful conclusion about the numerical range of A. In fact, we are not aware of any tight relationship between W (A) and the numerical ranges of A, A + bc * . Writing

T -1 AT = A A + bc * ,
A = A 0 0 A + 0 bc * 0 bc * = A 0 0 A + b b [0 * , c * ]
we obtain from [29, Section 1.0.1 & Property 1.2.10] the inclusion

W (A) ⊆ W (A) + W (uv * ), u = b b , v = 0 c
where + refers to the Minkowski sum of sets. In general, the field of values of a rank-one matrix uv * is an ellipse with focal points 0 and v * u and minor semi-axis

1 2 ( u 2 v 2 -|v * u| 2 ) 1/2
, which in our special case amounts to focal points 0 and c * b and minor semi-axis 1 2 (2 b 2 c 2 -|c * b| 2 ) 1/2 . In Section 6.2, we will illustrate for an example of practical relevance that W (A) can have rather undesirable properties, to the extent that Theorem 4.2 is of little help in understanding the convergence of our algorithms. In the next section, we therefore derive convergence bounds by an approach that avoids the dependence on W (A) and only depends on W (A) and W (A + bc * ). While the second set may still be larger than W (A), it is at least always smaller than W (A).

Convergence results based on integral representations.

We begin by considering results based on the Cauchy integral formula in Section 5.1 and then focus on the special case of Markov functions in Section 5.2.

Convergence results based on the Cauchy integral formula.

Let f be analytic on a domain E containing the eigenvalues of A and A+bc * . We recall (2.8):

f (A + bc * ) -f (A) = - 1 2πi Γ f (z)(zI -A) -1 bc * (zI -A -bc * ) -1 dz, (5.1)
with Γ = ∂E. The integrand in (5.1) involves solutions of shifted linear systems. Letting

(zI -A)x (z) = b and (zI -A -bc * ) * y (z) = c, (5.2)
we can compactly write the right-hand side of (5.1) as

f (A + bc * ) -f (A) = 1 2πi Γ f (z)x (z)y (z) * dz. (5.3)
Now, consider the FOM [START_REF] Saad | Krylov subspace methods for solving large unsymmetric linear systems[END_REF] approximations for (5.2), given by

x m (z) := b U m (zI -G m ) -1 e 1 , y m (z) := c V m (zI -H m -c e 1 b * V m ) -1 e 1 ,
where we used the Arnoldi decompositions (2.1) and (2.4). Recalling that X m (f ) is defined as the (1,2) block of

f G m b c e 1 e * 1 0 H * m + c V * m be * 1
and using contour integration, as in the proof of Lemma 2.2, we find that

U m X m (f )V * m = - 1 2πi Γ f (z)x m (z)y m (z) * dz (5.4)
In other words, the approximation to the update matrix can be interpreted as the integral over outer products of FOM approximations for the shifted linear systems (5.2).

Inserting (5.3) and (5.4) into the definition (4.2) of the error gives

E m (f ) = - 1 2πi Γ f (z) x (z)y (z) * -x m (z)y m (z) * dz = - 1 2πi Γ f (z) x (z)(y (z) -y m (z)) * + (x (z) -x m (z))y m (z) * dz. (5.5)
Convergence estimates can now be obtained by taking norms in the contour integral (5.5) and bounding all occurring terms. To do so, we first introduce some notation. In the following, E is a compact, convex set containing both W (A) and W (A + bc * ). The closed unit disk is denoted by D. Denoting by C := C ∪ {∞} the extended complex plane, let ψ be the conformal mapping from C \ D onto C \ E, normalized at infinity such that ψ(∞) = ∞, ψ (∞) > 0. Note that ψ (w) exists for almost all w on ∂D. We are now in the position to state the following norm bounds on the quantities involved in (5.5). Lemma 5.1. For all m ≥ 2 and z = ψ(u) ∈ C \ E we have

max x (z) b , x m (z) b , y (z) c , y m (z) c ≤ 1 dist(z, E) ≤ |u|/ψ (∞) (|u| -1) 2 , (5.6) max x (z) -x m (z) b , y (z) -y m (z) c ≤ 4|u| 1-m |ψ (u)|(|u| 2 -1) ≤ 4|u| -m dist(z, E) , (5.7)
where dist(z, E) denotes the distance of z from E.

Proof. The first inequality in (5.6) follows for x (z) immediately from the fact that

x (z) = (zI -A) -1 b ≤ (zI -A) -1 b ≤ b dist(z, E) ,
and analogously for x m (z), y (z), and y m (z). For the second inequality in (5.6), we recall a result by Kühnau (see Theorem 3.1 in [START_REF] Toh | The Kreiss matrix theorem on a general complex domain[END_REF] and its proof):

|u| -1 1 + 1/|u| ≤ dist(ψ(u), E) |ψ (u)| ≤ |u| 2 -1 |u| , (5.8)
which, in fact, does not require E to be convex. By convexity of E, we also have the inequality

ψ (u) ψ (∞) -1 ≤ 1 |u| 2
(5.9) due to Grötzsch and Golusin; see [START_REF] Kővari | On Faber polynomials and Faber expansions[END_REF]Section 2]. Combining these two inequalities leads to

1 dist(z, E) ≤ 1 + 1/|u| |ψ (u)|(|u| -1) ≤ |u| ψ (∞)(|u| -1) 2 ,
showing (5.6).

We now turn to the second set (5.7) of inequalities. We will make use of properties of the Faber transform F of a function analytic in the open unit disk and continuous on the closed unit disk; see [START_REF] Beckermann | Error estimation and evaluation of matrix functions via the Faber transform[END_REF][START_REF] Ellacott | On the Faber transform and efficient numerical rational approximation[END_REF] for more details. Letting

G(v) = 1 ψ (u)(u-v) we obtain g(ζ) := F(G)(ζ) = 1 2πi ∂E G(ψ -1 ( ζ)) d ζ ζ -ζ (5.10) = 1 2πi |v[=1 G(v) ψ (v) dv ψ(v) -ζ = 1 ψ(u) -ζ ,
where the last equality follows from the residue theorem; see also [5, Eqns (2.5), (2.7)].

Let P be defined by the formula

P (v) = u ψ (u)(|u| 2 -1) v u m-1 + 1 ψ (u)(v -u) v u m-1 -1 ,
depending on the parameter u. Then it is straightforward to verify that P is a polynomial of degree at most m -1, and that

G(v) -P (v) = 1 ψ (u)(|u| 2 -1) vu -1 u -v v u m-1
vanishes at 0 because m ≥ 2. Thus, with the notation of [5, §2], and p = F(P ) ∈ Π m-1 we get that F(G -P ) = F + (G -P ) = g -p, using (5.10). Now Theorem 2.1 in [START_REF] Beckermann | Error estimation and evaluation of matrix functions via the Faber transform[END_REF], being related to the fundamental work of Crouzeix [START_REF] Crouzeix | Numerical range and functional calculus in Hilbert space[END_REF][START_REF] Crouzeix | The numerical range as a spectral set[END_REF], implies

(zI -A) -1 -p(A) = (g -p)(A) ≤ 2 max |v|=1 |G(v) -P (v)| = 2|u| 1-m |ψ (u)|(|u| 2 -1)
.

Since W (G m ) ⊂ W (A), the same upper bound is obtained for (zI -G m ) -1 -p(G m ) .

Because of p ∈ Π m-1 and the exactness property (3.4), we have that p(A)b = b V m p(G m )e 1 and thus

x (z) -x m (z) b = x (z) -p(A)b -(x m (z) -b V m p(G m )e 1 ) b ≤ (zI -A) -1 -p(A) + (zI -G m ) -1 -p(G m ) ≤ 4|u| 1-m |ψ (u)|(|u| 2 -1)
,

that is, we have shown the first inequality of (5.7). The second inequality follows from a combination with the second inequality in (5.8). The inequalities in (5.7) involving y (z) instead of x (z) are proven in a completely analogous fashion.

In order to state our first explicit convergence bound, consider the (compact) level sets E r , which are defined via their complements as

E r = C \ E c r , where E c r := {z ∈ C \ E : |ψ -1 (z)| > r}. Theorem 5.2. Suppose that f is analytic in E R . Then, for 1 < r < R E m (f ) ≤ R -m-1 16/ψ (∞) (1 -r/R)(1 -1/r) 3 max z∈Γ R |f (z)| b c , where Γ R = ∂E R .
Proof. Let Γ = ∂E r . Then we have, by taking norms in (5.5),

E m (f ) ≤ 1 2π Γ |f (z)| x (z) y (z)-y m (z) + y m (z) x (z)-x m (z) | dz|. (5.11)
Because of the polynomial exactness of the Krylov approximation (2.5) proven in Theorem 3.2, the error of our Krylov approximation is the same for f and f -p, for any p ∈ Π m . This allows us to conclude from (5.11) that

E m (f ) ≤ f -p Γ 2π Γ x (z) y (z) -y m (z) + y m (z) x (z) -x m (z) | dz|.
(5.12) Applying the substitution z = ψ(u) with |u| = r, we obtain from Lemma 5.1 the upper bound

|u|=r x (ψ(u)) y (ψ(u)) -y m (ψ(u)) + y m (ψ(u)) x (ψ(u)) -x m (ψ(u)) |ψ (u) du| 2π ≤ b c |u|=r 8r 2-m ψ (∞)(r 2 -1)(r -1) 2 | du| 2π ≤ 8r -m-1 ψ (∞)(1 -1/r) 3 b c (5.13)
for the integral in (5.12). Further, we can use a partial Faber sum (cf. [START_REF] Beckermann | Error estimation and evaluation of matrix functions via the Faber transform[END_REF]Remark 3.3] with ρ = R/r) to find the Bernstein-type estimate for the best polynomial approximation on the convex set E r

min p∈Πm f -p Γ ≤ r R m 2r R -r max z∈Γ R |f (z)|.
(5.14)

Inserting (5.13) and (5.14) into (5.12) then yields the desired result. Let us illustrate Theorem 5.2 for some particular sets E and the particular case of the exponential function f (z) = exp(z). We suppose in the following that E is convex and symmetric with respect to the real axis. Then ψ(R) ∈ R is the element in E R with the largest real part, and hence

R -m-1 max z∈Γ R |f (z)| = e ψ(R) /R m+1 .
Our aim will be to choose R > 1 to make the above right-hand side small. This task has been (implicitly) accomplished by Hochbruck and Lubich [28, Section 3] for various families of sets like real and purely imaginary intervals, disks, and wedgeshaped sets with a corner at ψ(1), see also the analysis of [START_REF] Beckermann | Error estimation and evaluation of matrix functions via the Faber transform[END_REF]Section 4] where also general convex domains with corners have been considered.

Our upper bounds will be stated in terms of ψ(1), the element of E of largest real part, and in terms of ρ = ψ (∞), the logarithmic capacity of E (which is increasing as the set becomes larger). For the four families of sets mentioned above, explicit formulas are known for ψ(R) in terms of ψ(1), ρ and R.

We start with some negative result. By convexity, the function r → rψ (r) is known to be increasing. Provided that ψ (1) ≥ m + 1, elementary calculus implies that R → e ψ(R) /R m+1 is increasing for R ∈ [1, ∞), and thus Theorem 5.2 is not useful in this case. However, for m + 1 2ρ ≥ ψ (1), we get the following result of superlinear convergence.

Corollary 5.3. Let E be convex and symmetric with respect to the real axis. Then, under the conditions of Theorem 5. Proof. From the Grötzsch and Golusin inequality (5.9), we find for 1 ≤ r ≤ R ψ (r) -ρ 1 + 1 r 2 ≤ 0, where ρ = ψ (∞).

(5.16)

Integrating (5.16) from 1 to R then yields

ψ(R) ≤ ψ(1) + ρ R - 1 R .
By arguing as in the proof of [START_REF] Hochbruck | On Krylov subspace approximations to the matrix exponential operator[END_REF]Theorem 4] we find for R = m+1 ρ that e ψ(R) /R m+1 ≤ e ψ(1) ( eρ m+1 ) m+1 . In addition, since (1

-1/ √ R) -4 ≤ (1 -e -1/2 ) -4 ≤ 42, we conclude that 16/ρ (1 -1/ √ R) 4 e ψ(R) R m+1 ≤ 672 ρ e ψ(1) ρe m + 1 m+1 ,
and our claim follows from Theorem 5.2 with r = √ R. Example 5.4.

We illustrate the result of Corollary 5.3 by a simple numerical experiment. We choose A ∈ C 100×100 as a diagonal matrix with eigenvalues equidistantly spaced in [-20, 0] and b as a random vector of unit norm, resulting in Λ(A -bb * ) ⊆ [-20.2, 0] =: E. We apply Algorithm 1 (modified to account for the update -bb * ) for approximating exp(A -bb * ) -exp(A) and report the resulting convergence curve, together with the error bound (5.15) in Figure 5.1. We note that the convergence rate is predicted very accurately, but that the magnitude of the error is severly overestimated due to the large constant in (5.15), which we expect to not be optimal. We thus expect that in general the convergence slope of our bound will be quite accurate, while there is a (large) constant distance between the actual convergence curve and the error bound, something which happens quite frequently for bounds based on the field of values. For practical purposes, we thus suggest to ignore the constant in (5.15).

We remark that by combining the result of Theorem 4.1 with the result of [5, Corollary 4.1] on polynomial approximation of the exponential function, we find a bound that predicts essentially the same convergence rate as the bound (5.15), albeit with a much smaller constant. The real use of the integral representation is thus in the non-Hermitian case, where it allows to circumvent the reliance on W (A).

In view of known results for disks (see [28, Example after Theorem 5]), we do not expect E m (f ) to be small in general for m + 1 < ψ (1). However, in the case of a corner at ψ(1), the right-most element of E, we have a different regime of convergence for m + 1

ρ. Here we consider only the wedge-like set E = E(α, ρ) with

ψ(w) = ψ(1) + ρw 1 - 1 w α , 1 < α ≤ 2, (5.17)
having an outer angle of απ at ψ(1), by slightly improving [28, Theorem 6] (though following [5, Corollay 4.2] one could include more general E having such an angle).

Notice that

E(2, ρ) = [ψ(1) -4ρ, ψ (1)
] is a real interval of capacity ρ. Corollary 5.5. Consider the wedge-like set E = E(α, ρ) for α ∈ (1, 2] and ρ > 0. Then, under the conditions of Theorem 5.2 and f (z) = exp(z), we have for

m + 1 -4/α ∈ [αρ 1/α , αρ] E m (f ) ≤ (4ρ 1/α ) 4 ρ exp ψ(1) -(α -1) m + 1 -4 α αρ 1/α α α-1 b c . (5.18)
Proof. We consider for R > 1 the strictly increasing function

u = u(R) = ψ(R) -ψ(1) ρ 1/α = R 1/α -R -1+1/α , and notice that u ≤ 1 implies that R 1/α = u + R -1+1/α ≤ u + 1 ≤ 2. Hence, provided that u ∈ (0, 1), e ψ(R)-ψ(1) R m+1 (1 -1/ √ R) 4 ≤ 16e ψ(R)-ψ(1) R m+1 (1 -1/R) 4 ≤ 16 u 4 exp ρu α + m + 1 - 4 α log 1 R ≤ 16 u 4 exp ρu α + m + 1 - 4 α 1 R -1 = 16 u 4 exp ρu α -m + 1 - 4 α u R 1/α ≤ 16 u 4 exp ρu α -m + 1 - 4 α u 2 .
The argument of the exponential function on the right takes its minimum at by our assumption on m. Hence

u α-1 = m + 1 -4/α 2αρ ∈ [ρ -α-1 α , 1]
e ψ(R)-ψ(1) R m+1 (1 -1/ √ R) 4 ≤ 16 u 4 exp -(α -1) m + 1 - 4 α u 2α ≤ 16ρ 4/α exp -(α -1) m + 1 -4 α 2αρ 1/α α α-1
, and our claim follows from Theorem 5.2. Example 5.6.

We illustrate the result of Corollary 5.5 for a diagonal matrix A ∈ C 1000×1000 with eigenvalues in the wedge-like set (5.17) with ψ(1) = 0, α = 3/2 and ρ = 100. The vector b is again chosen as a random vector of unit norm, which results in W (A -bb * ) ⊆ E where E is a wedge-like set corresponding to ψ(1) = 0, α = 1.5 and ρ = 101. We apply Algorithm 2 (modified to account for the update -bb * ) for approximating exp(A -bb * ) -exp(A) and report the resulting convergence curve together with the error bound (5.18) in Figure 5.2. We again observe that the error norm is overestimated by a few orders of magnitude, while the convergence slope is predicted quite well (although not as accurately as in the Hermitian case in Example 5.4).

In the particular case of symmetric data, A = A * and c = b, the smallest set containing both W (A) and

W (A + bb * ) is the interval E = [λ min (A), λ max (A + bb * )],
where λ min (•) and λ max (•) denote the smallest and largest eigenvalues of a Hermitian matrix. For the exponential function a refined analysis would be possible, using the fact that the coefficients of the Chebyshev series of f (z) = exp(z) are explicitly known in terms of Bessel functions. We believe, however, that we do not gain much qualitatively compared to the results of the two preceding corollaries with the choices

ψ(1) = λ max (A + bb * ), ρ = λ max (A + bb * ) -λ min (A) 4 , α = 2.
5.2. Convergence results for Markov functions. The approach outlined above extends to other integral representations of f . In particular, a Markov function can be written as

f (x) = β α dµ(z) x -z , ( 5
f (z) = z -γ = sin(γπ) π 0 -∞ (-x) -γ dx z -x
for γ ∈ (0, 1) or the logarithm

f (z) = 1 z log(1 + z) = -1 -∞ (-1/x) dx z -x ,
see, e.g., [START_REF] Ch | Potential Theory on Locally Compact Abelian Groups[END_REF][START_REF] Henrici | Applied and Computational Complex Analysis[END_REF] for more details and other examples of Markov functions.

Combining the result of Lemma 5.1 with the integral representation of the error (5.5) now allows to obtain convergence estimates.

Theorem 5.7. Let E be symmetric with respect to the real axis and let ω ∈ E ∩ R denote the element of E with smallest real part. For a Markov function f with α < β < ω in (5. [START_REF] Estrada | Subgraph centrality in complex networks[END_REF]), it holds that

E m (f ) ≤ 8 |f (ω)| |φ(β)| m b c ,
where φ is mapping conformally from C \ E onto C \ D.

Proof. Let us first show that we may take Γ = [α, β] in equation (5.5). Take r > 1 such that f is analytic in E r , that is, [α, β] ∩ E r = ∅. Since all expressions x (z), y (z) * , x m (z), y m (z) * are analytic outside E and decay like 1/z at ∞, we get by exchanging integration and using the Cauchy residual theorem

E m (f ) = 1 2πi ∂Er β α dµ(t) t -z x (z)(y (z) -y m (z)) * + (x (z) -x m (z))y m (z) * dz = β α x (t)(y (t) -y m (t)) * + (x (t) -x m (t))y m (t) * dµ(t).
By taking norms, we obtain

E m (f ) ≤ β α ( x (t) y (t) -y m (t) + y m (t) x (t) -x m (t)) ) dµ(t) (5.20)
Applying the first inequality in (5.6) and the second inequality in (5.7) to the individual terms in (5.20) then yields

E m (f ) ≤ 8 b c β α |u| -m dµ(t) dist(t, E) 2 .
(5.21)

Using the fact that t = ψ(u), i.e., u = φ(t) and the fact that dist(t, E) = ω -t for t ∈ [α, β], the inequality (5.21) yields

E m (f ) ≤ 8 b c β α 1 |φ(t)| m dµ(t) (ω -t) 2 .
Since the function 1/|φ(t)| is monotonically increasing on [α, β] (see, e.g., the proof of Theorem 6.1 in [START_REF] Beckermann | Error estimation and evaluation of matrix functions via the Faber transform[END_REF]), we further have

E m (f ) ≤ 8 b c |φ(β)| m β α dµ(t) (ω -t) 2 .
(5.22)

The result of the theorem now follows by noting that the integral on the right-hand side of (5.22) is exactly |f (ω)|, see, e.g., [START_REF] Alzer | Some classes of completely monotonic functions[END_REF]. We now show how the bound from Theorem 5.7 simplifies when E is an ellipse on the right of [α, β] (or, as a special case, an interval). This gives a more explicit idea of the convergence behavior one can expect.

Corollary 5.8. Let f be a Markov function (5.19) with -∞ ≤ α < β < ∞, let E be an ellipse

E = {z ∈ C : |z -σ + τ | + |z -σ -τ | ≤ τ (ρ + ρ -1 )}, (5.23)
where σ ∈ R, τ > 0, ρ ≥ 1.

Let

ω := σ - τ 2 (ρ + ρ -1
).

(5.24)

If β < ω and E is symmetric to the real axis we have

E m (f ) ≤ 8 |f (ω)| b c ρ • √ κ -1 √ κ + 1 m , κ = |β -σ| + τ |β -σ| + τ . (5.25)
Proof. Obviously, ω defined in (5.24) is the element of smallest real part in E ∩ R. The conformal mapping φ for E is given by the inverse Joukowski mapping

φ(z) = ρ -1 ζ + ζ 2 -1 , ζ = z -σ τ .
By straight-forward algebraic manipulations, we find that

1 |φ(β)| = ρ √ κ -1 √ κ + 1
with κ as defined in (5.25). Applying Theorem (5.7) then gives the desired result. A further simplification is possible in the Hermitian positive definite case, where the ellipse E is an interval. We assume β ≤ 0 and bound 1/|φ(β)| ≤ 1/|φ(0)| in the following result, as this gives a result which closely resembles the classical convergence bound for the conjugate gradient method.

Corollary 5.9. Let f be a Markov function (5.19) with -∞ ≤ α < β ≤ 0, let A be Hermitian positive definite and let b = c * . Further, let

κ * = λ max (A + bb * ) λ min (A)
. We then have Proof. In the Hermitian positive definite case, W (A) = [λ min (A), λ max (A)] and

E m (f ) ≤ 8 |f (0)| b 2 √ κ * -1 √ κ * + 1 m . ( 5
W (A + bb * ) = [λ min (A + bb * ), λ max (A + bb * )]. As λ min (A) ≤ λ min (A + bb * ) and λ max (A) ≤ λ max (A + bb * ), we can thus take E = [λ min (A), λ max (A + bb * )].
This is a special case of an ellipse (5.23) with

ρ = 1, σ = τ = λ max (A + bb * ) -λ min (A) 2 .
(5.27)

Noting that λ min (A) is the smallest element in E ∩ R = E and inserting (5.27) into (5.25) gives the desired result after some simple calculations.

Example 5.10. To illustrate the result of Corollary 5.9, consider a diagonal matrix A ∈ R 100×100 with equidistantly spaced eigenvalues in the interval [0.1, 10]. As in the previous examples, we choose b as a random vector of unit norm, which results in Λ(A + bb * ) ⊆ [0.1, 10.1] =: E. We approximate (A + bb * ) -1/2 -A -1/2 and report the resulting convergence curve, together with the error bound (5.26) in Figure 5.3. We can observe the typical shortcoming of bounds of the form (5.26) for Markov functions; they only predict linear convergence, while superlinear convergence due to spectral adaption can be observed in this case, see, e.g., [START_REF] Beckermann | Superlinear convergence of the rational Arnoldi method for the approximation of matrix functions[END_REF]. In addition, as we also observed for the exponential function, the large constant in the bound (5.26) leads to a overestimation of the order of magnitude of the error.

Remark 5.11. We make two remarks concerning Corollary 5.9.

Applications and numerical experiments.

In this section, we present possible applications for the developed methods, with a special focus on the computation of communicability measures in network analysis. All experiments are performed in MATLAB R2016b on a Linux machine with Intel Core i7-6700 CPU and 32 GB main memory.

6.1. Updating network communicability measures. Matrix functions, especially the matrix exponential, play an important role in network analysis, see [START_REF] Estrada | Network properties revealed through matrix functions[END_REF] and the references therein.

Given an undirected graph G = (V, E) with V = {1, . . . , n} and E ⊆ V × V , we let A ∈ R n×n be the adjacency matrix of G defined by a ij = a ji = 1 if (i, j) ∈ E and a ij = 0 otherwise. Note that A is symmetric. Introduced in [START_REF] Estrada | Subgraph centrality in complex networks[END_REF], the subgraph centrality of node i is given by [exp(A)] ii trace(exp(A))

.

This quantity represents a weighted average of the number of closed walks which connect node i to itself, see, e.g., [START_REF] Estrada | Network properties revealed through matrix functions[END_REF] for details. An interesting question is how the subgraph centralities in a graph change when an edge is added to or removed from the graph; see for example [START_REF] Arrigo | Updating and downdating techniques for optimizing network communicability[END_REF], where the effect of adding/removing edges on the socalled total communicability is studied. The addition of an edge (i, j), with i = j, corresponds to the rank-two modification

A + BC * , where B = [e i , e j ], C = [e j , e i ], (6.2)
of the adjacency matrix. Analogously, the removal of an edge (i, j) corresponds to a rank-two modification. In turn, these tasks fit perfectly into the framework considered in this paper.

In our experiments:

• Following [START_REF] Benzi | Quadrature rule-based bounds for functions of adjacency matrices[END_REF][START_REF] Benzi | Ranking hubs and authorities using matrix functions[END_REF], the diagonal entries of exp(A), needed for evaluating (6.1) for the original matrix A, are estimated by combining the Lanczos method with quadrature [START_REF] Benzi | Bounds for the entries of matrix functions with applications to preconditioning[END_REF][START_REF]Matrices, Moments and Quadrature with Applications[END_REF]. Specifically, we use the implementation of these methods from the mmq toolbox [START_REF] Meurant | MMQ toolbox[END_REF]. • To handle the rank-two modification (6.2) we perform two consecutive rankone updates with Algorithm 1, as outlined in Remark 2.3. Note that while B = C, the resulting matrix BC * is Hermitian. Thus, we can perform a preprocessing step in order to obtain two Hermitian rank-one updates and then apply Algorithm 1. As said before, it is not necessary to explicitly form the matrix U m X m (exp)V * m , we only evaluate its diagonal entries. We now suppose that the diagonal of the matrix exponential has been computed beforehand (in an expensive offline calculation). We then add or remove edges from the graph and compare the computation time needed for updating them using Algorithm 1 to the time needed for recomputing them from scratch using the mmq toolbox. When using the mmq toolbox, we perform five Lanczos iterations per diagonal entry (this value is, e.g., also used in [START_REF] Benzi | Quadrature rule-based bounds for functions of adjacency matrices[END_REF]) and can be expected to give a rather rough approximation of the exact value. When using Algorithm 1, we aim for an accuracy of 10 -6 according to the stopping criterion from Section 2.3, evaluated with d = 2. We observed that typically between 10 and 30 iterations of the algorithm are necessary to fulfill this criterion. We use five different networks available from the SuiteSparse Matrix Collection [START_REF] Davis | The SuiteSparse matrix collection[END_REF], and for each network we randomly choose ten edges that are added to or removed from the network (i.e., in order to compute the update, we have to call Algorithm 1 a total of 20 times). Table 6.1 summarizes the obtained results. For all networks under consideration we observe that our algorithm clearly outperforms recomputation.

To also illustrate how the convergence profile of the method looks like and to judge the quality of the error estimate (2.10) in a practical situation, we show the convergence curve of the first update performed for the network Pajek/USpowerGrid together with the error estimate in Figure 6.1. We can observe that the method converges very smoothly and that the error estimate is almost identical to the exact error norm for all iterations. 6.2. One-dimensional convection diffusion equation. In the following, we illustrate the convergence of Algorithm 2 for a non-Hermitian example. Consider the one-dimensional convection diffusion equation

u -c • u = f in [0, 1],
u(0) = u(1) = 1, Therefore, our convergence estimates cannot be expected to be sharp in this case. We also report the shape of W (A) in Figure 6.2 in order to illustrate the superiority of the approach from Section 5 over the result of Theorem 4.2, which depends on the field of values of A. Comparing W (A + bc * ) and W (A), we observe that the latter set is substantially larger, especially for large values of c. An especially unfavorable case appears for c = 60, where it turns out that W (A) is not a subset of the left complex half-plane anymore (and contains 0), unlike W (A) and W (A + bc * ). For an entire function like the exponential function, this only leads to worse convergence bounds. For other functions, like the inverse square root having a singularity at 0, this can be more problematic, to the extent that Theorem 4.2 is not applicable. Figure 6.2 also shows the difference-based error estimate (2.10) for d = 2. While this estimate closely follows the exact error curve in later iterations, it severely underestimates the error in the first few iterations. This is due to the fact that the method almost stagnates in these iterations, and is a typical shortcoming of difference-based error estimates.

Conclusions.

We have proposed Krylov subspace methods for approximating f (A+D)-f (A) for a low-rank matrix D. We proved that the resulting approximation is exact for a polynomial of a certain degree and used this to derive a variety of convergence results, which either link the convergence of our method to polynomial approximation problems or exploit results on the error in the full orthogonalization method in conjunction with an integral representation of f .

In numerical experiments-in particular on applications from network analysiswe have illustrated that our approach can dramatically reduce the cost of computing f (A + D) or portions thereof. We expect that our algorithms will prove useful in other application areas and we will explore this in future work. Another interesting topic for future research is the use of extended and rational Krylov subspaces in the proposed method and an analysis of the resulting convergence behavior.

2. 1 .

 1 Hermitian rank-1 updates. Because it is conceptually simpler, we first discuss the Hermitian case: A = A * and b = c.

Fig. 2 . 1 .

 21 Fig. 2.1. Convergence curves of simple Lanczos and Lanczos with full reorthogonalization for diagonal matrices with equidistantly spaced (left) and logarithmically spaced (right) eigenvalues in [10 -3 , 10 3 ].

Remark 2 . 3 .

 23 There are two different ways of extending our approach from a rank-one modification to a general rank-k modification f (A + D) with D = BC * and B, C ∈ C n×k : 1. By letting b 1 , . . . , b k and c 1 , . . . , c k denote the columns of B and C, respectively, we can write D as a sum of k rank-one matrices, D = k i=1 b i c *

3 Fig. 2 . 2 .

 322 Fig. 2.2. Comparison of exact error and difference-based error estimates for the update of the inverse square root described in Example 2.4.

  bc * ) for all = 0, . . . , m -1. (3.5) Multiplying (3.3) by U m from the left and by V * m from the right then gives

  T having the very modest 2-norm condition number κ(T ) = 1+

Fig. 5 . 1 .

 51 Fig. 5.1. Exact error norm and error bound (5.15) for the update of the matrix exponential described in Example 5.4.

  2 and f (z) = exp(z), we have for m+1 ≥ eρ E m (f )

Fig. 5 . 2 .

 52 Fig. 5.2. Exact error norm and error bound (5.18) for the update of the matrix exponential described in Example 5.6.

Fig. 5 . 3 .

 53 Fig. 5.3. Exact error norm and error bound (5.26) for the update of the matrix inverse square root described in Example 5.10.

Fig. 6 . 1 .

 61 Fig. 6.1. Convergence curve of Algorithm 2 for computing a rank-one update of the matrix Pajek/USpowerGrid together with the error estimate (2.10).

  .19) where µ is a positive measure with support in the interval [α, β] with -∞ ≤ α < β < ∞. Any such Markov function is analytic in C \ [α, β]. Examples of Markov functions include inverse fractional powers

  1. Similar to what we pointed out already for the exponential function, if we use known approximation results for Markov functions, like, e.g., [5, Theorem 6.1 & Remark 6.3] together with Theorem 4.1, we obtain the same convergence factor (with a different constant) as in Corollary 5.9 in the Hermitian case. 2. It is possible to obtain a slightly sharper version of (5.26) in which the constant 8 is replaced by 4, and κ * is replaced by the maximum of the Euclidean norm condition numbers of A and A+bb * by explicitly exploiting the connection to the conjugate gradient method [26] and integrating over the classical CG convergence bounds, see, e.g., [20, Lemma 4.1 & Theorem 4.3] for a similar technique. We refrain from giving the details of this, as the improvement over the bound (5.26) is quite marginal in most cases.

Table 6 . 1

 61 Number n of nodes and number k of edges of the networks used in our experiments, together with the computation time for updating the subgraph centrality of all nodes, when modifying 10 edges, using Algorithm 1 vs. recomputing them from scratch using the mmq toolbox.

	Network			n	k Algorithm 1	mmq
	Gleich/Minnesota	2,642	6,606	0.15 s	1.80 s
	Pajek/Erdos992	6,100	15,030	0.62 s	6.22 s
	Pajek/USpowerGrid	4,941	13,188	0.48 s	4.96 s
	SNAP/ca-HepTh	9,877	51,971	1.25 s	14.58 s
	SNAP/email-Enron	36,692 367,662	1.96 s 147.47 s
		10 2
					Exact error norm
		10 0			Error estimate, d = 2
		10 -2
	Error norm	10 -6 10 -4
		10 -8
	10 -10	0	5	10	15
					m

discretized by centered finite differences combined with a second-order scheme for the convection term. We choose the convection coefficient to be c = 10 and discretize the equation on 256 interior grid points. We then consider the rank-one modifications A + bc * of the discretized operator that change the convection coefficient c in the middle, at the 128th grid point, to (a) c = 20, (b) c = 40, (c) c = 60. On the lefthand side of Figure 6.2, the field of values of A, A + bc * and the matrix A from (2.6) is shown for these modifications, and on the right-hand side the convergence history of Algorithm 2 for approximating exp(A + bc * ) -exp(A) is given.

As shown in Figure 6.2, the field of values of A lies in a wedge-like set in left halfplane, as considered in Corollary 5.5, with a rather small inner angle at the origin. The angle increases when considering A + bc * and this increase becomes more pronounced as c grows. In light of Corollary 5.5, we would thus expect slower convergence for larger values of c. However, as the convergence curves on the right-hand side of Figure 6.2 show, the convergence behavior is virtually identical for all three test cases.


Cagatay Haratoka
email: cagatay.haratoka01@univ-poitiers.fr

Karl Joulain

Younès Ezzahri
email: younes.ezzahri@univ-poitiers.fr

Thermal rectification of a vanadium dioxide based conductive thermal diode under transient Dirichlet conditions

Keywords: Thermal rectification, thermal diode, phase change materials (PCM), vanadium dioxide, thermal conductivity, specific heat

We focus in this work on analysing the behaviour of a thermal diode and consider its operation in the conductive heat transfer regime. The thermal rectification of a thermal diode made of VO2 and SiO2, is investigated in detail by changing different intrinsic and extrinsic parameters under the Dirichlet boundary conditions. The behaviour of the thermal diode is numerically examined according to each parameter individually. We found that the thermal rectification factor can rise to very high values in the first moments of the transient process before reaching steady-state conditions. This feature makes the study of the thermal diode behaviour under transient conditions very promising to pave the way for a better control of practical thermal management systems.

Introduction

Thermal management is becoming one of the most important engineering endeavours, especially after encountering problems such as global warming, energy crisis and system safety [START_REF] Li | Transforming heat transfer with thermal metamaterials and devices[END_REF]. Thermal management is essential to keep the temperature within a certain level for various applications for both performance, longevity of operation and safety such as electronic packages [START_REF] Li | Transforming heat transfer with thermal metamaterials and devices[END_REF], photovoltaic panels [START_REF] Browne | Phase change materials for photovoltaic thermal management[END_REF], batteries [START_REF] Rao | A review of power battery thermal energy management[END_REF] and thermal storage systems [START_REF] Pasupathy | Phase change material-based building architecture for thermal management in residential and commercial establishments[END_REF]. The thermal diode, which is analogous to the electronic diode, is an attractive and encouraging device to generate thermal rectification for different thermal management applications. It is a two-terminal system that rectifies the heat flux depending on the direction of the temperature gradient between the boundaries of the two terminals. An ideal thermal diode allows heat transfer in one direction according to the difference of temperatures at the terminals, and it blocks it in the opposite direction when the temperatures are exchanged [START_REF] Wong | A review of state of the art thermal diodes and their potential applications[END_REF][START_REF] Wehmeyer | Thermal diodes, regulators, and switches: Physical mechanisms and potential applications[END_REF][START_REF] Roberts | A review of thermal rectification observations and models in solid materials[END_REF].

Thermal rectification is a process that occurs as soon as the device thermal properties depend on the temperature so that the heat flux depends on the direction of the temperature gradient [START_REF] Roberts | A review of thermal rectification observations and models in solid materials[END_REF]. The topic was firstly studied in 1936 by Starr [START_REF] Starr | The copper oxide rectifier[END_REF]. Starr examined thermal rectification of Copper/Cuprous Oxide interfaces experimentally. Among thermal rectification devices, researches on thermal diodes have dramatically increased in the past decade [START_REF] Li | Transforming heat transfer with thermal metamaterials and devices[END_REF][START_REF] Wong | A review of state of the art thermal diodes and their potential applications[END_REF][START_REF] Wehmeyer | Thermal diodes, regulators, and switches: Physical mechanisms and potential applications[END_REF]. There are thermal diodes based on the principles of the three modes of heat transfer: conduction, convection and radiation [START_REF] Wong | A review of state of the art thermal diodes and their potential applications[END_REF][START_REF] Wehmeyer | Thermal diodes, regulators, and switches: Physical mechanisms and potential applications[END_REF][START_REF] Roberts | A review of thermal rectification observations and models in solid materials[END_REF]. Graphene based heterostructures, in which thermal rectification occurs due to asymmetrical phonon propagation at the material interface, are widely under investigation in the literature as an example for the principle of conduction [START_REF] Ding | Interfacial thermal conductance in graphene/mos2 heterostructures[END_REF][START_REF] Li | Interface thermal transport of graphene-based intralayer heterostructures[END_REF][START_REF] Chen | Highly efficient thermal rectification in carbon/boron nitride heteronanotubes[END_REF][START_REF] Zhang | Thermal transport in graphene/h-bn lateral heterostructures with interface compositional diffusion[END_REF][START_REF] Liang | Multilayer in-plane graphene/hexagonal boron nitride heterostructures: Insights into the interfacial thermal transport properties[END_REF]. In addition, there are also studies on functionally graded materials to obtain thermal rectification [START_REF] Go | On the Condition for Thermal Rectification Using Bulk Materials[END_REF][START_REF] Gupta | Recent development in modeling and analysis of functionally graded materials and structures[END_REF][START_REF] Naebe | Functionally graded materials: A review of fabrication and properties[END_REF][START_REF] Dettori | Thermal rectification in silicon by a graded distribution of defects[END_REF][START_REF] Yousefi | Non-equilibrium molecular dynamics study on radial thermal conductivity and thermal rectification of graphene[END_REF]. These are the materials whose composition along an axis changes and causes thermal rectification due to different material properties along the axis. Although functionally graded materials have a potential for mass production due to being constructed by one material, they need more experimental studies [START_REF] Wong | A review of state of the art thermal diodes and their potential applications[END_REF]. Moreover, thermal rectification can also be obtained using superconductor materials, but there are limited studies about thermal diodes based on superconductors [START_REF] Wong | A review of state of the art thermal diodes and their potential applications[END_REF].

Phase change materials (PCM) are another alternative to generate thermal rectification. It is worth mentioning that PCM have different contexts in different communities. They include solid-liquid transition as well as crystal structure transformation (crystalline-amorphous transition and dielectric-metallic transition). Consequently, they manifest different thermal, electrical and optical properties according to their critical temperature of transition [START_REF] Raoux | Phase change materials[END_REF]. Therefore, thermal diodes based on PCM constitute a beneficial solution for thermal management due to their altering characteristics. Among PCM, vanadium dioxide (VO2), which is one of the most popular PCM, has recently drawn great attention thanks to the closeness of its critical transition temperature to room temperature. It is an ideal candidate to manipulate conductive heat flux and in which thermal rectification can occur. VO2 shows a metal-insulator transition (MIT) around a critical temperature, which is To = 341 K [START_REF] Morin | Oxides which show a metal-to-insulator transition at the neel temperature[END_REF]. The crystalline structure of VO2 shows monoclinic and tetragonal-rutile crystalline structures below and above To, respectively [START_REF] Andersson | Studies on vanadium oxides. ii. the crystal structure of vanadium dioxide[END_REF]. The metallic phase results in a high thermal conductivity above To, whereas the dielectric phase results in a low thermal conductivity below it [START_REF] Oh | Thermal conductivity and dynamic heat capacity across the metal-insulator transition in thin 27 film vo2[END_REF].

In the literature, available studies mostly focused on exploring conductive thermal rectification under steady-state conditions. Yang et al. [START_REF] Yang | Optimal thermal rectification of heterojunctions under fourier law[END_REF] modelled a thermal diode that consists of two materials with linear temperature-dependent thermal conductivities based on Fourier's law of heat conduction and carried out a study to optimize the thermal rectification. They obtained a maximum rectification ratio of 0.5 [START_REF] Yang | Optimal thermal rectification of heterojunctions under fourier law[END_REF]. However, reaching this high rectification ratio is difficult and costly, because there is a need for a high temperature difference larger than 100 K to obtain it [START_REF] Yang | Optimal thermal rectification of heterojunctions under fourier law[END_REF]. As a practical advantage of PCM, this drawback can be overcome thanks to their thermal conductivity alteration to some extent. Thus, there are detailed investigations to use VO2 in thermal diodes as a single PCM or in combination with another PCM. [START_REF] Kasali | Conductive thermal diode based on two phase-change materials[END_REF][START_REF] Kang | Thermal rectification via heterojunctions of solid-state phase-change materials[END_REF][START_REF] Cottrill | Dual phase change thermal diodes for enhanced rectification ratios: Theory and experiment[END_REF][START_REF] Ordonez-Miranda | Conductive thermal diode based on the thermal hysteresis of vo2 and nitinol[END_REF][START_REF] Kasali | Spherical and cylindrical conductive thermal diodes based on two phase-change materials[END_REF]. Ordonez-Miranda et al. [START_REF] Ordonez-Miranda | Conductive thermal diode based on the thermal hysteresis of vo2 and nitinol[END_REF] examined a thermal diode with a single PCM, (VO2 or nitinol) under steady-state conditions. They reached 19.7% of thermal rectification and concluded that a higher thermal conductivity contrast, smaller thermal hysteresis and faster phase transition can generate higher thermal rectification. Moreover, Kasali et al. [START_REF] Kasali | Conductive thermal diode based on two phase-change materials[END_REF] theoretically investigated the thermal diode that was made by two PCM as the combination of VO2 and Polyethylene. In their study, they reached a thermal rectification of 60%. Later, they analysed the influence of the diode geometry. They obtained a thermal rectification of 63.5% and 63.2% for spherical and cylindrical thermal diodes, respectively, that were made by VO2

and Polyethylene. These values of thermal rectification, which are very close, are higher than the one obtained in the plane geometry [START_REF] Kasali | Spherical and cylindrical conductive thermal diodes based on two phase-change materials[END_REF].

Nevertheless, not much effort has been devoted to the study of conductive thermal rectification under transient conditions. The latter needs to be well understood to be able to utilize thermal components practically and efficiently.

Herrera et al. [START_REF] Herrera | Thermal Rectification Under Transient Conditions: The Role of Thermal Capacitance and Thermal Conductivity[END_REF] published the first study related to thermal diodes under transient conditions. They calculated the thermal rectification of thermal diodes whose thermal conductivity varies linearly with temperature by applying finite difference method. Also, they introduced time-dependent thermal rectification factor and observed a higher rectification factor under transient conditions [START_REF] Herrera | Thermal Rectification Under Transient Conditions: The Role of Thermal Capacitance and Thermal Conductivity[END_REF].

Klinar et al. [START_REF] Klinar | Toward a solidstate thermal diode for room-temperature magnetocaloric energy conversion[END_REF] carried out a numerical study to investigate the influence of a thermal diode in a magnetocaloric refrigeration device. They analysed various temperature dependent thermal conductivity profiles and found out that the thermal rectification reaches higher values under transient conditions than steady-state conditions.

There are many studies that examine the thermal rectification of thermal diodes in steady-state conditions and the interest for transient analysis is increasing. However, we have not encountered any prior work that studied thermal diodes while considering properties of VO2 in full aspects under transient conditions. Therefore, the objective as well as the novelty of this study is to investigate a thermal diode based on a single PCM, under transient conditions.

The thermal diode is made of VO2 and SiO2 and we consider the Dirichlet boundary conditions to be applied. The thermal conductivity and the specific heat of VO2 depend strongly on temperature, mainly in the vicinity of the critical transition temperature [START_REF] Morin | Oxides which show a metal-to-insulator transition at the neel temperature[END_REF]. On the other hand, the changes of the thermal conductivity and the specific heat of SiO2 with temperature are negligible. The thermal rectification of the thermal diode is investigated in detail by changing different intrinsic and extrinsic parameters.

Theoretical modelling

We consider a thermal diode in a plane geometry that is made of a PCM and non-PCM in this study as stated before. We assume this thermal diode to be placed between two thermal baths in a vacuum chamber so that we can neglect any heat exchange due to convection. We also assume heat exchange due to radiation to be negligible and consider conduction to be the only heat transfer mode between the two thermal baths. Therefore, there is a conductive heat exchange between the two thermal baths due to their temperature difference Th -Tc, as illustrated in

Nondimensionalization of the transient process

The nonlinear heat diffusion equation for VO2 is given by:

𝜌 𝑐 (𝑇) = 𝑘(𝑇) (1)
The thermal conductivity, k(T ), and specific heat capacity, cp(T ), of VO2 vary 120 with temperature. The thermal conductivity of VO2 can be described using the following equation [START_REF] Ordonez-Miranda | Conductive thermal diode based on the thermal hysteresis of vo2 and nitinol[END_REF].

𝑘(𝑇) = 𝑘 + ( ) (2)
where kd and km represent the thermal conductivity of VO2 in the dielectric and metallic phases, respectively. To is the phase critical transition temperature and β is the phase transition slope. The values of these parameters are as 125 follows: 3.6 W.m 1 .K 1 for kd, 6 W.m 1 .K 1 for km and 1.6 K 1 for β. The values of kd and km were measured experimentally by Oh et al. [START_REF] Oh | Thermal conductivity and dynamic heat capacity across the metal-insulator transition in thin 27 film vo2[END_REF], while the value of β corresponds to the best fit of the experimental data using Eq.2 [START_REF] Ordonez-Miranda | Conductive thermal diode based on the thermal hysteresis of vo2 and nitinol[END_REF].

In addition, T0 takes the values Toh = 342.9 K and Toc = 334.9 K during the heating and cooling processes respectively due to 8 K hysteresis effect of VO2 between these processes [START_REF] Ordonez-Miranda | Conductive thermal diode based on the thermal hysteresis of vo2 and nitinol[END_REF]. Therefore, Toh and Toc represent the transition temperatures of the heating and cooling processes, respectively.

As reported in Ref [START_REF] Ordonez-Miranda | Modeling of the electrical conductivity, thermal conductivity, and specific heat capacity of vo2[END_REF], the expression of the specific heat capacity of VO2 is given by:

𝑐 (𝑇) = 𝑐 (𝑇) + 𝑇 + 𝐿 (3)
where cd, ρv, δ, γ and Ln represent the specific heat capacity predicted by Debye model, density of VO2, coefficient of thermal expansion, isothermal compressibility and latent heat, respectively. Since the variation of cd within the range of the considered boundary conditions (BC) is negligible, we assume it to be constant and we take the average of the values of cd at the BC temperatures as 0.704 kJ.kg 1 .K 1 . Moreover, ρv of VO2 changes according to its crystalline structure. However, for simplicity sake throughout the study, we assume ρv to be also constant and we take the average of the values of the dielectric monoclinic and metallic tetragonal-rutile crystalline structures as 4612 kg.m 3 . Ln equals 51.45 kJ.kg 1 for the heating process of VO2 [START_REF] Berglund | Electronic properties of vo2 near the semiconductor-metal transition[END_REF] and we assume the same value for the cooling process. As explained in the study of Ordonez-Miranda et al. [START_REF] Ordonez-Miranda | Modeling of the electrical conductivity, thermal conductivity, and specific heat capacity of vo2[END_REF], the second term in the right side of Eq.3 can be neglected, because the values of this term are much smaller than the values of cd within the phase transition temperatures interval. Thus, its contribution is neglected throughout the study. The temperature behaviours of the thermal conductivity and specific heat capacity of VO2 are reported in Figs. 2(a) and 2(b), respectively. Eq.1 can be nondimensionalized easily as

1 + ∆ 𝜃 * = * [1 + (𝑘 -1)𝑓(𝜃)] * (4)
where the dimensionless temperature is defined by θ = (T -To)/∆To, the dimensionless length by x * = x/l and the dimensionless time by t * = 𝛼 t/ l 2 . ∆To is the difference between the transition temperatures of heating and cooling processes: ∆To = Toh -Toc. 𝛼 is the average thermal diffusivity between the one of the dielectric structure of VO2 and SiO2, which is 𝛼 = (αd + αs)/2. αd and αs represent the thermal diffusivity of the dielectric structure of VO2 and SiO2, respectively. They are given by αd = kd/(ρvcd) and αs = ks/(ρscs). ρv and ρs are the densities of VO2 and SiO2 in order. The thermal conductivity, ks, specific heat capacity, cs, and density, ρs, of SiO2 are 1.5 W.m 1 .K

Eq.(4-5) are numerically solved by using implicit finite difference method and discretized in time and space for many nodes. Finally, the normalized heat flux density (ratio of the heat flux density to the quantity ∆To/l) can be evaluated by using Eq.6 and Eq.7 for forward and backward configurations.

𝛷 * = -𝑘 * (6) 𝛷 * = -𝑘 [1 + (𝑘 -1)𝑓(𝜃)] * (7)
where + and -signs represent forward and backward fluxes, respectively. Even though the quantities expressed in Eq.6 and Eq.7 have the dimension of a thermal conductivity, we will refer to them as nondimensional heat flux densities henceforth in this paper.

The interface thermal resistance describing the temperature discontinuity at the interface between the two slabs of VO2 and SiO2 is also modelled according to the study of Yuan et al. [START_REF] Yuan | Heat transfer analysis in multi-layered materials with interfacial thermal resistance[END_REF]. The authors developed two one-dimensional models, that represent the interface thermal resistance with a virtual layer or an artificial layer, to analyse the heat transfer in multi layered materials including the interface thermal resistance. In our study, we create an artificial layer with a width of ∆x between the two materials. Half of the length of the artificial layer is overlaid in the left material and the other half is in the right material, which are VO2 and SiO2, respectively. The artificial layer has the modified thermal properties according to the features of VO2 and SiO2 to represent their effects at the interface of these two materials as explained in the study of Yuan et al. [START_REF] Yuan | Heat transfer analysis in multi-layered materials with interfacial thermal resistance[END_REF]. We define an effective thermal conductivity, ke, that takes into account the effects of the thermal conductivities of the artificial layer, VO2 and SiO2 as [START_REF] Yuan | Heat transfer analysis in multi-layered materials with interfacial thermal resistance[END_REF] = + + [START_REF] Starr | The copper oxide rectifier[END_REF] ki is the thermal conductivity of the artificial layer and can be defined as ki = ∆x/r, where r is the interface thermal resistance. Details of this model are explained in the reference [START_REF] Yuan | Heat transfer analysis in multi-layered materials with interfacial thermal resistance[END_REF].

Nondimensionalization of the steady-state process

The heat conduction equation is solved under steady-state conditions semianalytically in order to validate the results of the steady-state regime that were numerically obtained during the transient regime analysis. The following Eqs. [START_REF] Ding | Interfacial thermal conductance in graphene/mos2 heterostructures[END_REF][START_REF] Li | Interface thermal transport of graphene-based intralayer heterostructures[END_REF][START_REF] Chen | Highly efficient thermal rectification in carbon/boron nitride heteronanotubes[END_REF][START_REF] Zhang | Thermal transport in graphene/h-bn lateral heterostructures with interface compositional diffusion[END_REF] are obtained after Eq.1 was nondimensionalized under steady-state conditions:

𝜉 𝜃 , + 𝑘 + ( ) 𝜃 , = 𝜉(𝜃 ) + 𝑘 𝜃 +

𝜉 𝜃 , + 𝑘 + ( ) 𝜃 , = 𝜉(𝜃 ) + 𝑘 𝜃 +

where 𝜉(𝜃) = (𝑘 -1)𝑙𝑛 1 + 𝑒 ∆ / [𝛽∆𝑇 ]. 𝜃 , and 𝜃 , are the nondi- mensional temperatures of the hot and cold sides at the interface of two materials in the forward and backward configurations, respectively, as illustrated in Fig. 1 above. θh and θc are the hot and cold nondimensional BC temperatures, respectively. In addition, η is defined by r/l.

Eqs. 9 and 10 allow determining 𝜃 , and 𝜃 , as a function of θh, θc and the rest of the system parameters. Once 𝜃 , and 𝜃 , are determined, we can obtain straightforwardly 𝜃 , and 𝜃 , via the following Eqs. 11 and 12:

𝜃 , = , (11)
𝜃 , = , (12)

Thermal rectification analysis

In the absence of an internal heat source, the heat flux remains constant ss and backward heat fluxes to the maximum heat flux between the two configurations. However, this is not valid for transient conditions as emphasized by Herrera et al. [START_REF] Herrera | Thermal Rectification Under Transient Conditions: The Role of Thermal Capacitance and Thermal Conductivity[END_REF], because heat flux can vary throughout the system under these conditions. Thus, we analysed the thermal rectification on the basis of the heat flux density at the cold end by considering the fact that heat spontaneously moves from the hot side to the cold side according to the second law of thermodynamics. Consequently, we define the thermal rectification as

𝑅(𝑡 * ) = | * * | * , * (13)
The thermal rectification is evaluated using three different metrics as proposed by Herrera et al.

where 𝑡 * and 𝑡 * represent the times for which Rss and the first peak in the thermal rectification profile are reached, respectively. It is also worth mentioning that the first peak in the thermal rectification profile as a function of time coincides very often with Rmax. As a general rule, we considered t * to be the time at which the full numerical solution of the thermal rectification in the transient regime analysis equals 99% of the semi-analytical solution obtained under steady-state conditions. In the case of the steady-state analysis, we can easily show that the thermal rectification factor is given by:

𝑅 = ( ) , , , , , (15)

Results and discussion

In our analysis of the time dependent thermal rectification behaviour of the conditions, we explore the consequences of two initial conditions: T (t = 0) = Tc leading to a heating process of the full system and T (t = 0) = Th which leads to a cooling process. We investigate the impact of varying six different parameters; BC change, β, Ln, kr, ∆To and r. During the study, only one parameter is changed in each analysis and the rest is kept unchanged. We firstly consider four different BC cases according to the real features of VO2 and SiO2. Then, we deeply analyse the remaining parameters with five different cases individually by considering the BC to be Th = 354 K and Tc = 324 K. The cases for each parameter are shown in Table 1.

In the first analysis of the four different BC, we consider Th and Tc for the heating and cooling conditions, separately, according to the real characteristic values of VO2 and SiO2. We fix these values with respect to To, which takes the value of 342.9 K for Toh and the value of 334.9 K for Toc, as stated before. In order to check the effect of the phase transition of VO2 on the thermal rectification in detail, we specify the values of temperatures for the BC in various intervals surrounding the phase transition area of VO2, which was previously shown in

(i) (ii) (iii) (iv) (v) β [K 1 ] 0.5 1 1.6 3 5 Case 3 -Ln (i) (ii) (iii) (iv) (v)
Ln [kJ.kg 1 ] 66.89 (+30%) 59.17 (+15%) 51.45 43.73 (-15%) 36.02 (-30%) Case 4 -kr

(i) (ii) (iii) (iv) (v) kr 1.2 1.5 1.667 3 5 Case 5 -∆T (i) (ii) (iii) (iv) (v) ∆To [K] 2 5 8 13 17

Heating analysis

Toh [K] 336.9 339.9 342.9 347.9 351.9

Cooling analysis

Toc [K] 340.9 337.9 334.9 329.9 325.9

Case 6 -r The specific heat capacity cp has a direct effect on heat transfer under transient conditions. In the case of VO2 and as we discussed it above, the latent heat Ln plays a key role in the behaviour of cp. We report in Figs. 4(a Indeed, it is interesting to note that in the first four cases of the cooling process, 𝑅 increases by increasing ∆To. Nevertheless, the fifth case shows a sudden decrease. This sudden decrease can also be seen in the fifth case of the heating process. Besides, in the fourth and fifth cases, the amplitude of the first peak in the thermal rectification time profile in the heating process, is lower than the value at the steady-state regime. A behaviour that we encountered also while varying the ratio kr as shown in Fig. 5(a). Moreover, Rss decreases by increasing ∆To for both heating and cooling processes. The analysis until this point was carried out by assuming there is no interface thermal resistance between the two slabs of VO2 and SiO2 composing the conductive thermal diode. This constitutes the sixth and last parameter we considered in our parametric study and it is modelled by using the artificial layer approach of Yuan et al. [START_REF] Yuan | Heat transfer analysis in multi-layered materials with interfacial thermal resistance[END_REF] as mentioned in the theory section above.

(i) (ii) (iii) (iv) (v) r [m 2 .K.W 1 ] 10 1 10
The general trend is similar for the different values of r for heating and cooling processes as reported in Figs.7(a) and 7(b), respectively. Although this effect is not obvious in the cooling process, we get higher Rmax and lower Rss in the heating process when the value of r increases. Furthermore, especially in the cooling process, R tends to saturate at its steady-state value faster by increasing r. Nevertheless, this effect seems to be almost insignificant except between the first two cases. Similarly, it is only between these two cases, where the timeintegrated thermal rectification factor 𝑅 shows a dependence on r in the cooling process. For the rest of the cases, 𝑅 remains almost constant as for Rmax and Rss. We summarize all the calculated values of Rmax, Rss and 𝑅 for the different cases of each of the six parameters investigated in this study in Table 2 below.

i) (ii) (iii) (iv) (v) (i) (ii) (iii) (iv) (v) (
i) (ii) (iii) (iv) (v) (i) (ii) (iii) (iv) (
Rmax 0.23 0.21 0.20 0.18 0.16 0.75 0.74 0.73 0.71 0.70 Rss 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.

(i) (ii) (iii) (iv) (v) (i) (ii) (iii) (iv) (v)
Rmax 0.23 0.20 0.20 0.20 0.20 0.72 0.72 0.72 0.72 0.72 Rss 0.10 0.11 0.12 0.12 0.12 0.09 0.11 0.12 0.12 0.12 R̃ 0.09 0.09 0.09 0.09 0.09 0.59 0.56 0.56 0.56 0.56

Conclusion

In this study, we have theoretically analysed the thermal rectification factor of a conductive thermal diode in a plane geometry, that is made of VO2 and

Fig. 1 .

 1 The heat flux flows from the PCM to the non-PCM in the forward configuration. Conversely, it flows in the opposite direction in the backward configuration. Heat fluxes in forward and backward configurations are driven by the thermal conductivities of the PCM and non-PCM. The total length of the model is 2l and both materials have equal lengths, which mean half of the total length, as shown in Fig.1.

Figure 1 :

 1 Figure 1: Schematic representation of the studied thermal diode in the forward and backward configurations.

Figure 2 :

 2 Figure 2: Computed behaviours of the (a) thermal conductivity and (b) specific heat capacity of VO2 as functions of temperature.

  [START_REF] Herrera | Thermal Rectification Under Transient Conditions: The Role of Thermal Capacitance and Thermal Conductivity[END_REF]: (a) The first metric is Rmax which represents the maximum thermal rectification factor, (b) the second metric is Rss which is the thermal rectification factor at steady-state conditions, and (c) the third metric is the time-integrated thermal rectification factor, 𝑅 , which allows a better quantification of the capacitive effects, 𝑅 = ∫ 𝑅(𝑡 * )𝑑𝑡 * * *

Fig. 2 .

 2 Fig.2. We report in Figs.3, the computed time behaviours of the nondimensional heat flux density at the cold boundary [Figs.3(a) and 3(c)] in the forward and backward configurations and the thermal rectification factor [Figs.3(b) and 3(d)] for both the heating [Figs.3(a) and 3(b)] and cooling [Figs.3(c) and 3(d)] processes of the full thermal diode. As our modelling is based on Fourier's law, we only account for times beyond the diffusion time in order to avoid the drawback of the parabolic nature of the heat equation. The diffusion time corresponds to the time required for the temperature front to diffuse from the hot boundary to the cold one and can roughly be estimated as l 2 /𝛼 . The general trend is the same for the different values of BC in each process and each configuration. In the heating process, the nondimensional heat flux density increases with time to saturate at the steady-state values for both forward and backward configurations, while there is a decreasing behaviour in the cooling process. As expected the value of nondimensional heat flux density increases by increasing the tem-

2 10 3 10 4 10 5 14 heat

 14 flux density seems to saturate faster (slower) by increasing ∆T for the forward (backward) configuration in the heating process, the opposite occurs in the cooling process. The behaviour of the thermal rectification factor R shows interesting nonlinear features, which are consequences of the behaviour of the 275 nondimensional heat flux densities in the forward and backward configurations for each process. In both the heating and cooling processes, R manifests a maximum Rmax, goes to zero, then increases again to reach the steady-state

Figure 3 :

 3 Figure 3: Computed time behaviours of the nondimensional heat flux density and thermal rectification factor in the heating [(a) and (b)] and cooling [(c) and (d)] processes, respectively, for different values of BCs.

) and 4 (

 4 b), the computed time behaviours of the thermal rectification factor R in the heating and cooling processes, respectively, for different values of Ln. The main features in the time profile of R, remain similar to the ones shown in Figs.3(b) and 3(d). While varying Ln affects more Rmax in the heating process than in the cooling process similarly to the effect of β, it influences the time at which R vanishes in both processes. Additionally, Rmax and 𝑅 decrease by decreasing Ln for both processes. In the steady-state regime, the effect of Ln disappears and all the curves collapse on each other as expected. cp has no impact on the heat transfer under steady-state conditions.

Figure 4 :

 4 Figure 4: Computed time behaviours of the thermal rectification factor in the heating (a) and cooling (b) processes, for different values of the parameter Ln.

Figure 5 :

 5 Figure 5: Computed time behaviours of the thermal rectification factor in the heating (a) and cooling (b) processes, for different values of the parameter kr .

Figure 6 :

 6 Figure 6: Computed time behaviours of the thermal rectification factor in the heating (a) and cooling (b) processes, for different values of the parameter ∆To.

Figure 7 :

 7 Figure 7: Computed time behaviours of the thermal rectification factor in the heating (a) and cooling (b) processes, for different values of the parameter r.

  SiO2 considering the Dirichlet boundary conditions under transient conditions.The heat flux density behaviour at the cold side of the thermal diode and its consequence on the thermal rectification factor of the thermal diode have been investigated deeply by considering various intrinsic and extrinsic parameters. It has been shown that the thermal conductivity and the specific heat capacity have a major impact on the nonlinear time behaviour of the thermal rectification under transient conditions. Furthermore, we obtained higher thermal rectification under transient conditions than steady-state conditions for almost all parameters values. There is an effect of the interface thermal resistance on the thermal rectification for both steady-state and transient conditions, but this effect seems to be minor generally in comparison with the other intrinsic parameters of the conductive thermal diode. Considering the reasonable outcomes of this theoretical study, further research should be conducted on this type of conductive thermal diodes especially under different transient boundary conditions in order to demonstrate the full benefits of the transient analysis for a better control of practical thermal management systems. Pulsed and dynamical boundary conditions constitute a further step to understand and explore to get more insights onto the behaviour of the conductive thermal diode based on a phase change material and the phase change effects on the thermal rectification under these conditions. Acknowledgement This article pertains to the France 2030 program "Investissements d'Avenir" (LABEX INTERACTIFS, reference ANR-11-LABX-0017-01). This article pertains to the French government program "Investissements d'Avenir" (EUR IN-TREE, reference ANR-18-EURE-0010).

Table 1 :

 1 Investigated parameters and cases.

		Case 1 -BC
	(i)	(ii)	(iii)	(iv)
		Heating analysis

Table 2 :

 2 The calculated values for the maximum thermal rectification factor, thermal recti-

	fication factor at steady-state conditions and time-integrated thermal rectification factor for
	all investigated cases.
			Heating analysis		Cooling analysis
				Case 1 -BC
		(i)	(ii)	(iii) (iv)	(i)	(ii)	(iii) (iv)
	Rmax 0.55 0.40 0.33 0.29	0.83 0.78 0.74 0.72
	Rss	0.11 0.12 0.12 0.12	0.11 0.12 0.12 0.12
	R̃	1.84 0.86 0.35 0.23	4.58 2.54 1.88 1.54
				Case 2 -β

throughout the system under steady-state conditions. Therefore, thermal rectification can be calculated using the ratio of the difference between forward

the different parameters, but they all impact the whole time profile of the thermal rectification, ranging from the early transient regime to the steady-state regime. Rmax, Rss and 𝑅 all change by varying each of the parameters. 𝑅 which quantifies the capacitive effect of the conductive thermal diode, gives a central information about the duration and magnitude of the thermal rectifica-395 tion dynamics span.

suggested VO2 based conductive thermal diode under the transient Dirichlet

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.