Jialun Zhou

Salem Said

Yannick Berthoumieu

Riemannian information gradient methods for the parameter estimation of ECD : Some applications in image processing

Keywords: elliptically-contoured distribution, Riemannian information gradient, large-scale dataset, image re-colorization, texture classification

Elliptically-contoured distributions (ECD) play a significant role, in computer vision, image processing, radar, and biomedical signal processing. Maximum likelihood estimation (MLE) of ECD leads to a system of non-linear equations, most-often addressed using fixed-point (FP) methods. Unfortunately, the computation time required for these methods is unacceptably long, for large-scale or high-dimensional datasets. To overcome this difficulty, the present work introduces a Riemannian optimisation method, the information stochastic gradient (ISG). The ISG is an online (recursive) method, which achieves the same performance as MLE, for large-scale datasets, while requiring modest memory and time resources. To develop the ISG method, the Riemannian information gradient is derived taking into account the product manifold associated to the underlying parameter space of the ECD. From this information gradient definition, we define also, the information deterministic gradient (IDG), an offline (batch) method, which is an alternative, for moderate-sized datasets. The present work formulates these two methods, and demonstrates their performance through numerical simulations. Two

Introduction

The family of Elliptically-contoured distributions (ECD) was originally introduced in [START_REF] Kelker | Distribution theory of spherical distributions and a locationscale parameter generalization[END_REF], and investigated in [START_REF] Fang | Generalized multivariate analysis[END_REF][START_REF] Fang | Symmetric multivariate and related distributions[END_REF]. It contains many widely-used statistical distributions, such as elliptical Gamma, Pearson type II, and elliptical multivariate logistic distributions. In terms of applications, the most popular classes of ECD are multivariate generalized Gaussian distributions (MGGD), and multivariate Student-T distributions [START_REF] Gómez | A multivariate generalization of the power exponential family of distributions[END_REF][START_REF] Sánchez-Manzano | A matrix variate generalization of the power exponential family of distributions[END_REF][START_REF] Kotz | Multivariate t-distributions and their applications[END_REF]. These are location-scale distributions, and are further parameterised by a shape parameter, or a degrees of freedom parameter.

MGGD are used in image processing, as models for wavelet and curvelet coefficients, and as models for three-channel color vectors, in image denoising, context-based image retrieval, image thresholding, texture classification, and image quality assessment [START_REF] Boubchir | Multivariate statistical modeling of images with the curvelet transform[END_REF][START_REF] Cho | Image denoising based on wavelet shrinkage using neighbor and level dependency[END_REF][START_REF] Verdoolaege | Multiscale colour texture retrieval using the geodesic distance between multivariate generalized gaussian models[END_REF][START_REF] Bazi | Image thresholding based on the em algorithm and the generalized gaussian distribution[END_REF][START_REF] Scharcanski | A wavelet-based approach for analyzing industrial stochastic textures with applications[END_REF][START_REF] Gupta | Generalized gaussian scale mixtures: A model for wavelet coefficients of natural images[END_REF]. MGGD are also used in video coding and denoising, radar signal processing, and biomedical signal processing [START_REF] Desai | Robust gaussian and non-gaussian matched subspace detection[END_REF][START_REF] Elguebaly | Bayesian learning of generalized gaussian mixture models on biomedical images[END_REF]. Some applications of Student-T distributions are presented in [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF], involving image denoising. In radar imaging, the Student-T distribution, so-called G 0 model within the family of spherically invariant random vectors (SIRVs), is largely exploited in the context of SAR or PolSAR imaging, for tasks such as despeckling, classification, segmentation or detection [START_REF] Bombrun | Hierarchical segmentation of polarimetric sar images using heterogeneous clutter models[END_REF][START_REF] Fernández-Michelli | Unsupervised polarimetric sar image classification using G 0 p mixture model[END_REF][START_REF] Chen | A novel statistical texture feature for sar building damage assessment in different polarization modes[END_REF].

Because ECD have been successful in real-world image and signal processing applications, much attention has been devoted to developing efficient methods for estimating their parameters. The vast majority of works, dedicated to this estimation problem, are focused on estimating the scatter matrix, considering the other parameters, i.e. location and shape parameters, as known.

In terms of maximum-likelihood estimation, two main classes of algorithms have been studied. Fixed-point (FP) algorithms, and gradient descent algorithms have been proposed, based on the geometric properties of the manifold of positive definite matrices [START_REF] Tyler | A distribution-free m-estimator of multivariate scatter[END_REF][START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF][START_REF] Sra | Geometric optimisation on positive definite matrices for elliptically contoured distributions[END_REF][START_REF] Zhang | Multivariate generalized gaussian distribution: Convexity and graphical models[END_REF][START_REF] Sra | Conic geometric optimization on the manifold of positive definite matrices[END_REF]. For MGGD, when the location parameter is equal to zero and the shape parameter is given, the uniqueness of the maximum-likelihood estimator has been shown, under a restriction on the value of the shape parameter [START_REF] Pascal | Parameter estimation for multivariate generalized gaussian distributions[END_REF]. In this case, a method of moments has also been developed [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized gaussian distributions with an application to multicomponent texture discrimination[END_REF]. For Student-T distributions, with a known degrees of freedom parameter, a fixed-point method for parameter estimation is given in [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF], where the existence and uniqueness, of location and scatter maximum-likelihood estimates, is shown for a fixed degrees of freedom parameter, superior to 1.

There is a shared drawback, in all of the maximum-likelihood estimation methods, just mentioned [START_REF] Sra | Geometric optimisation on positive definite matrices for elliptically contoured distributions[END_REF][START_REF] Zhang | Multivariate generalized gaussian distribution: Convexity and graphical models[END_REF][START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF][START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized gaussian distributions with an application to multicomponent texture discrimination[END_REF][START_REF] Pascal | Parameter estimation for multivariate generalized gaussian distributions[END_REF]. Specifically, these methods work well for datasets of moderate size and dimension, but require excessive resources in memory and time, for large-scale datasets, e.g. with a few millions of samples, an order of magnitude commonly encountered in optical or SAR image processing. This issue can be so severe as to make any of these methods inapplicable.

Mainly, this is due to the fact that all of these methods are off-line, or batch, estimation methods. They require access to the whole dataset, at once, for each iteration, and therefore consume increasing time and memory resources, in order to converge to a useful estimate, as the dataset grows large.

In order to overcome this drawback, the present work builds on the ideas from Riemannian stochastic optimisation, proposed in [START_REF] Bonnabel | Stochastic gradient descent on riemannian manifolds[END_REF][START_REF] Tripuraneni | Averaging stochastic gradient descent on riemannian manifolds[END_REF]. The problem of estimating the parameters of an ECD is viewed as the problem of minimising the Kullback-Leibler divergence, between the true (unknown) and estimate distributions. When this problem is addressed using Riemannian stochastic optimisation, each iteration of a stochastic optimisation method requires access to only one datapoint (one sample), instead of the whole dataset. In this way, the present work proposes a recursive method, for estimating the parameters of ECD (each time a new sample is processed, this sample is used to update the current estimate).

The proposed method will be called the information stochastic gradient (ISG). In its simplest form, it is an improvement of a previous method, used to estimate the scatter matrix, when the location and shape parameters are known [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a riemannian manifold[END_REF]. In this paper, we consider also the general case where the scatter matrix, the location and shape parameters are unknown. The ISG method relies on two main ideas :

• The greatest difficulty, in using recursive methods, is that they may require a careful choice of step-sizes. The standard Riemannian stochastic gradient method (as in [START_REF] Bonnabel | Stochastic gradient descent on riemannian manifolds[END_REF]), is very sensitive to the choice of step-sizes.

However, using the information gradient (also called the natural gradi-ent [START_REF] Amari | Information geometry and its applications[END_REF][START_REF] Amari | Natural gradient works efficiently in learning[END_REF]) leads to an automatic choice of step-sizes, which guarantees optimal performance. The ISG method implements the information gradient, relying on the Fisher information metric (or matrix), of the ECD model.

• The parameter space of an ECD model does not only contain the scatter matrix, but also location and shape parameters. In the case of MGGD or Student-T models, this parameter space is a product space, made up of triplets: (scatter matrix, location parameter, shape/degrees of freedom parameter). Since the geodesic curves of this space do not have a tractable expression [START_REF] Verdoolaege | On the geometry of multivariate generalized gaussian models[END_REF], an intuitive idea is to update each one of the three parameters, in its own turn, in an alternating fashion.

To understand the benefit of combining these two ideas, consider the special case of MGGD. In this case, a method of moments (MM) was used for the joint estimation of all three parameters [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized gaussian distributions with an application to multicomponent texture discrimination[END_REF], while their maximum-likelihood esitmation (MLE) was studied in [START_REF] Pascal | Parameter estimation for multivariate generalized gaussian distributions[END_REF]. It is well known that MLE performs better than MM in most scenarios [START_REF] Fuhrer | Estimating the linear-quadratic inventory model maximum likelihood versus generalized method of moments[END_REF]. However, as mentioned above, MLE cannot be applied to large scale datasets, due to its computational requirements. The ISG method strikes a balance between the low complexity of MM, and the stronger performance of MLE. For example, in the case where the scatter matrix and the location parameter are unknown, the complexity of ISG is comparable to that of the MM, while its performance is similar to that of the MLE, when the number of available samples is sufficiently large.

In other words, the size of the dataset is leveraged as a source of information, rather than as a computational burden.

The two ideas which underly the ISG method (discussed above), are also implemented in an offline (batch) method, called the information deterministic gradient (IDG) method. While its complexity (and therefore computation time) is much higher than the ISG method, the IDG method consistently outperforms other methods, even in the case where all three parameters are unknown.

The main results of the present work are given in Section 4. A detailed comparison of various estimation methods (MM, FP, ISG, IDG), based on computer experiments, is carried out in Section 5. Two image processing applications with real datasets are presented in Section 6.

First, Sections 2 and 3 define the general estimation problem for ECD models, and introduce some necessary geometric concepts.

The estimation problem

The ECD family

ECD is a general family of probability distributions that contains many important sub-families. The name ECD comes from the fact that when an ECD has a probability density function, the contours (level surfaces) of this function are ellipsoids.

The location, or expectation, parameter µ of an ECD determines the centre of these ellipsoids, while the axes of these ellipsoids are proportional to the eigenvalues of the inverse of the scatter matrix Σ. The shape parameter β determines the factor for this proportionality (β is the degrees of freedom parameter, for Student-T distributions).

Let X be a m-dimensional random vector that follows a ECD model. Denote θ = (µ, Σ, β) the parameters of this ECD, and Θ = R m × P m × R + its parametric space, where P m is the set of all symmetric positive definite matrices of size m × m. If X has a probability density function, then this takes on the following form

p(x; θ) = c(β) |Σ| -1/2 g [δ x (µ, Σ), β] (1)
where c(β) is a normalizing factor which depends only on β, and

δ x (µ, Σ) = (x -µ) † Σ -1 (x -µ).
The density generator g depends on the specific subfamily of ECD distributions, for example

g [δ x (µ, Σ), β] = exp - 1 2 δ β x for MGGD g [δ x (µ, Σ), β] = 1 + 1 β δ x -β+m 2 for Student-T 2
Recall the definition of the KL divergence

D(θ) = R p p(x; θ * ) ln p(x; θ * ) p(x; θ) dx = E θ * [ (θ * ; x)] -E θ * [ (θ; x)] (3)
Where (θ; x) = log p(x; θ) is the log-likelihood,

(θ; x) = α(β) - 1 2 log det(Σ) + h(δ x , β) (4)
with α(β) = log c(β) and h = log g. In the following, the KL divergence

(3) will be minimised using Riemannian information gradient descent. Some

Riemannian Information geometry concepts are recalled in the following section.

Necessary geometric concepts

The gradient descent method on Riemannian manifolds is based on the following update rule [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]

θ n+1 = R θn (α n+1 u(θ n )) (5)
Here, the smooth mapping R θ from the tangent space T θ Θ to Θ is required to be a retraction, in the sense that it verifies

R θ (0 θ ) = θ (6a) DR θ (0 θ ) = Id T θ Θ (6b)
where 0 θ denotes the zero element in T θ Θ, and Id T θ Θ denotes the identity mapping on T θ Θ. Each vector u(θ n ) belongs to the tangent space T θn Θ, and provides the direction of descent. In the present work, -u(θ n ) is the Riemannian information gradient, derived using the Fisher information metric.

The positive scalar α n is the step-size. The aim of equation ( 5) is to generate a sequence (θ n ) n≥0 ∈ Θ that converges to a stationary point θ * of the cost function (under some restrictions on u and α).

For our estimation problem, the model has three different parameters, which belong to three different Riemannian manifolds. Precisely, the parameter space is the product manifold Θ = R m × P m × R + . Therefore, tractable expressions of the Fisher information metric, and of the intrinsic geodesic map, on this product manifold, are needed. However, Θ does not support any such expressions [START_REF] Verdoolaege | On the geometry of multivariate generalized gaussian models[END_REF]. As the global Fisher information metric has not a closed form, we propose to use the product metric

u, v θ = u µ , v µ µ + u Σ , v Σ Σ + u β , v β β (7)
where

u = (u µ , u Σ , u β ) and v = (v µ , v Σ , v β ) are tangent vectors at the point θ = (µ, Σ, β). The metrics •, • µ , •, • Σ and •, • β are respectively the in-
trinsic Fisher information metrics of their corresponding sub-spaces. For the location parameter µ in R m , its information metric is expressed in terms of the usual Euclidean metric,

u µ , v µ µ = I µ u † µ Σ -1 v µ (8)
where •, • denotes the scalar product in Euclidean space. The information constant I µ is

I µ = - 4 m E ∂ 2 h(δ x , β) ∂δ 2 x δ x -2E ∂h(δ x , β) ∂δ x (9)
with h = log g. As for Σ ∈ P m , the Fisher information metric •, • Σ for the ECD model is defined by the Riemannian geometry of P m [START_REF] Berkane | Geodesic estimation in elliptical distributions[END_REF].

U Σ , V Σ Σ = I 1 tr(Σ -1 U Σ Σ -1 V Σ ) + I 2 tr(Σ -1 U Σ )tr(Σ -1 V Σ ) (10)
Here the constants I 1 > 0 and I 2 0 depend on the particular model under consideration, as follows

I 1 = 2A m(m + 2) I 2 = A m(m + 2) - 1 4 A = E ∂h(δ x , β) ∂δ x δ x 2 (11)
The shape factor β belongs to R + , so the Fisher information metric is given by

u β , v β β = I β u β v β (12)
with the information constant

I β = -E ∂ 2 α(β) ∂β 2 + ∂ 2 h(δ x , β) ∂β 2 (13)
Now, the information gradient ∇ θ D(θ) with respect to the product metric ( 7) is obtained by solving the following equation,

d D(θ) v = ∇ θ D(θ), v θ (14)
where the scalar product on the right-hand side is given by ( 7), and d is the differential form of D. Precisely, this product information gradient has the following form

∇ θ D(θ) = (∇ µ D(θ), ∇ Σ D(θ), ∇ β D(θ)) (15)
The first component ∇ µ D(θ) is expressed as

∇ µ D(θ) = -I -1 µ Σ E [G µ (θ; x)] (16)
where I µ is given in equation ( 9), and vector G µ (θ; x) is actually the gradient in the classic Euclidean sense

G µ (θ; x) = 2 ∂h(δ x , β) ∂δ x Σ -1 (x -µ) (17)
The second component ∇ Σ D(θ) is a bit more complicated (see Figure 3, for an illustration of the following computations)

∇ Σ D(θ) = -J -1 1 E [G Σ (Σ; x)] ⊥ -J -1 2 E [G Σ (θ; x)] (18) 10
where J 1 = I 1 and J 2 = I 1 + mI 2 , in terms of I 1 and I 2 given in [START_REF] Scharcanski | A wavelet-based approach for analyzing industrial stochastic textures with applications[END_REF], and where ⊥ and denote the following decomposition of G Σ (θ; x),

[G Σ (θ; x)] = 1 m tr Σ -1 G Σ (θ; x) Σ (19) [G Σ (θ; x)] ⊥ = G Σ (θ; x) -[G Σ (θ; x)] ( 20
)
in terms of

G Σ (θ; x) = - 1 2 Σ - ∂h(δ x , β) ∂δ x S x with S x = (x -µ)(x -µ) † (21)
Finally, for the third component,

∇ β D(θ) = -I -1 β E [G β (β; x)] (22)
where I β was given in [START_REF] Desai | Robust gaussian and non-gaussian matched subspace detection[END_REF], and

G β (θ; x) = ∂α(β) ∂β - ∂h(δ x , β) ∂β (23)
With regard to the retraction R θ , it will be defined as the product Riemannian exponential map,

R θ : T θ Θ -→ Θ u =      u µ u Σ u β      -→      Exp µ (u µ ) Exp Σ (u Σ ) Exp β (u β )      ( 24
)
where u is the direction of descent. The exponential map on R p (a Euclidean space) reduces to vector addition

Exp µ (u µ ) = µ + u µ (25)
The exponential map on P m is defined as follows [START_REF] Pennec | A riemannian framework for tensor computing[END_REF]:

Exp Σ (u Σ ) = Σ exp Σ -1 u Σ (26)
As for β, since it belongs to R + , the corresponding exponential map is a 1-dimensional version of ( 26)

Exp β (v β ) = β exp(β -1 v β ) (27)
All these three exponential map Exp verify the properties ( 6), therefore the their direct product (24) also verifies these properties, and is a well-defined retraction. Finally, the Riemannian distance associated to the metric ( 7) is given by

d 2 (θ 1 , θ 2 ) = d 2 R m (µ 1 , µ 2 ) + d 2 Pm (Σ 1 , Σ 2 ) + d 2 R + (β 1 , β 2 ) ( 28
)
For µ, the information distance is proportional to the Euclidean distance in

R p d 2 R m (µ 1 , µ 2 ) = I µ (µ 1 -µ 2 ) † (µ 1 -µ 2 ) (29a)
for µ 1 , µ 2 ∈ R m , with the constant I µ given by [START_REF] Verdoolaege | Multiscale colour texture retrieval using the geodesic distance between multivariate generalized gaussian models[END_REF]. For Σ, the information distance is defined as in [START_REF] Mostajeran | Ordering positive definite matrices[END_REF]

d 2 Pm (Σ 1 , Σ 2 ) = I 1 tr log(Σ -1 1 Σ 2 ) 2 + I 2 tr 2 log(Σ -1 1 Σ 2 ) (29b)
for Σ 1 , Σ 2 ∈ P m , where the constants I 1 and I 2 are given by [START_REF] Scharcanski | A wavelet-based approach for analyzing industrial stochastic textures with applications[END_REF], and the function log denotes the symmetric matrix logarithm. Finally, for β, the information distance is given by

d 2 R + (β 1 , β 2 ) = I β log 2 β -1 1 β 2 (29c)
for β 1 , β 2 ∈ R + , where I β is given in [START_REF] Desai | Robust gaussian and non-gaussian matched subspace detection[END_REF]. With the necessary geometric concepts now in place, the next section will introduce our estimation algorithms.

The IDG and ISG methods

This section will describe the IDG and ISG methods, and discuss their main properties. The IDG method (information deterministic gradient) is a deterministic gradient method, and the ISG method (information stochastic gradient) is a stochastic gradient method.

When the direction of descent is chosen according to [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF], the updated es-

timates θ k+1 = (µ k+1 , Σ k+1 , β k+1 ) rely on the current estimates θ k = (µ k , Σ k , β k ),
through the following alternating optimisation scheme step 1 :

µ k+1 ← (µ k , Σ k , β k ) step 2 : Σ k+1 ← (µ k+1 , Σ k , β k ) step 3 : β k+1 ← (µ k+1 , Σ k+1 , β k ) (30)

Deterministic gradient

The IDG method is a second-order offline method, somewhat similar to a Newton method. In the Newton method, the direction of descent is found by solving the Newton equation [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. In the IDG method, the Hessian in the Newton equation is approximated by the Fisher information metric (or matrix) I(θ).

Since IDG is an offline method, it choses a direction of descent which depends on the complete dataset. The cost function ( 2) is reformulated, by replacing the KL divergence, with the empirical average of -(θ, x n ). This empirical average is denoted by D(θ),

D(θ) = - 1 N N n=1 (θ n ; x n ) (31)
If the current estimate is θ k , the direction of descent is given by its three components

∇ µ D(θ k ) = -I -1 µ Σ 1 N N n=1 G µ (θ k ; x n ) (32a) ∇ Σ D(θ k ) = -J -1 1 1 N N n=1 [G Σ (θ k ; x n )] ⊥ -J -1 2 1 N N n=1 [G Σ (θ k ; x n )] (32b) ∇ β D(θ k ) = -I -1 β 1 N N n=1 G β (θ k ; x n ) (32c)
which are the same as ( 16), ( 18), ( 22), but with expectations replaced by empirical averages. Using the expressions (32a), (32b), (32c), the IDG algorithm can now be stated as follows.

Algorithm 1 Information Deterministic Gradient (IDG) algorithm

Input: A dataset X = (x 1 , • • • , x N ), an initialization θ 0 ∈ S 0 ⊂ Θ;
Output: The estimate θ;

1: for k = 0, 1, 2 • • • , K do 2: µ k+1 ← µ k -α µ ∇ µ D(µ k , Σ k , β k ); 3: Σ k+1 ← Σ k exp -Σ -1 k α Σ ∇ Σ D(µ k+1 , Σ k , β k ) ; 4:
β k+1 ← β k exp -β -1 k α β ∇ β D(µ k+1 , Σ k+1 , β k ) ; 5: end for 6: θ ← (µ K+1 , Σ K+1 , β K+1 )
In this algorithm, α denotes the step-size, which is selected according to the Armijo-Goldstein rule, and S 0 denotes a neighborhood of the true parameter value θ * . The following Proposition 1 states the convergence of Algorithm 1. is approximated by the Fisher information metric. Therefore, one should expect the θ k converge to θ * with a superlinear rate of convergence, just like the Newton method dose. Precisely, if θ = (Σ) or θ = (µ, Σ), with a fixed shape parameter β * , then, under the assumptions of Proposition 1, one should expect Algorithm 1 to generate a sequence (θ k ) k≥0 converging superlinearly to θ * . This is essentially due to Theorem 6.3.2 in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], and will be observed experimentally in Section 5 (see Figure 1).

Stochastic gradient

The ISG method is an online quasi-Newton method. For each update, only one sample or one mini-batch is used. Here, the cost function remains the same as in equation [START_REF] Fang | Symmetric multivariate and related distributions[END_REF]. For the current estimate θ k = (µ k , Σ k , β k ) the stochastic information gradients are

∇ µ (θ; x n ) = I -1 µ Σ G µ (x n ; θ) (33a) ∇ Σ (θ; x n ) =J -1 1 [G Σ (x n ; θ)] ⊥ + J -1 2 [G Σ (x n ; θ)] (33b) ∇ β (θ; x n ) = I -1 β G β (x n ; θ) (33c)
Accordingly, the expected direction of descent is E θ * [∇ θ (θ; x)], which is equal to 0 at the global minimum θ * . As in the classic stochastic gradient descent method, the step-size α n = a n is strictly positive, decreasing, and verifies the usual conditions

α n = ∞ α 2 n < ∞
Using the expressions (33a), (33b), (33c), the ISG algorithm can now be stated as follows.

Algorithm 2 ISG algorithm

Input: A dataset X = (x 1 , • • • , x N ), an initialization θ 0 ∈ S 0 ⊂ Θ, the coefficient a > 0;
Output: The estimate θ;

1: for n = 0, 1, 2 • • • , N do 2:
α n+1 ← a n+1 ;

3:

µ n+1 ← µ k + α n+1 ∇ µ (µ n , Σ n , β n , x n ); 4: Σ n+1 ← Σ n exp (Σ -1 n α n+1 ∇ Σ (µ n+1 , Σ n , β n , x n ));
5:

β n+1 ← β n exp (β -1 n α n+1 ∇ β (µ n+1 , Σ n+1 , β n , x n )); 6: end for 7: θ ← (µ N +1 , Σ N +1 , β N +1 );
Remark that, the descending direction is -∇ θ (θ; x), and the double negative Since D(θ) has an isolated stationary point at θ = θ * , the Hessian at point θ = θ * can be expressed in normal coordinates

H ij = ∂ 2 D ∂θ i ∂θ j θ i =0 (34)
The matrix H = (H ij ) is positive definite [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. With these notations, the rate of convergence is given by the following proposition.

Proposition 3. Under the assumptions of Proposition 2, if a > 1 2λ , where λ > 0 is the smallest eigenvalue of H,

E[d 2 (θ * , θ n )] = O(n -1 ) (35)
Here, d(•, •) stands for the product distance in [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a riemannian manifold[END_REF], and the "big O" notation means that there exist K > 0 and n 0 > 0 such that

∀n n 0 E[d 2 (θ * , θ n )] K n
In terms of the normal coordinates (θ i ), let the direction of descent ∇ θ (θ * ; x)

at the point θ = θ * have components (u i (θ * )). Let G * = (G * ij ), be the matrix G * ij = E θ * u i (θ * )u j (θ * ) (36)
Then, the following proposition gives the asymptotic normality of the ISG algorithm.

Proposition 4 (asymptotic normality). Under the assumptions of Propositions 2 and 3, the distribution of the re-scaled coordinates (n

1 2 θ i ) i∈{1,••• ,d}
converges to a centred d-variate normal distribution, where d is the dimension of Θ, with covariance matrix G given by the following Lyabunov equation

AG + GA = -a 2 G * (37)
Here

, A = (A ij ) with A ij = 1 2 δ ij -aH ij (δ denotes Kronecker's delta).
The proofs of Propositions 4 and 3 are discussed in Appendix C. For the case θ = (Σ) or θ = (µ, Σ), the product metric ( 7) coincides with the information metric of the ECD model. Then, the assumptions of Proposition 5 in [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a riemannian manifold[END_REF] are satisfied, and the following corollary may be obtained.

Corollary 1. For the ECD model, parameterised by θ = (Σ) or θ = (µ, Σ), with a fixed β * , the product metric [START_REF] Boubchir | Multivariate statistical modeling of images with the curvelet transform[END_REF] coincides with the information metric.

1. the rate in equation ( 35) holds, whenever a > 1/2.

2. if a = 1 the distribution of the re-scaled coordinates (n 1/2 θ i ) converges to a centred d-variate normal distribution, with covariance matrix equal to the identity G * = I d , and the recursive estimates θ n are asymptotically efficient.

Note that, Item 2) of Corollary 1 implies that the distribution of nd 2 (θ * , θ n ) converges to a χ 2 -distribution with d degrees of freedom.

nd 2 (θ * , θ n ) ⇒ χ 2 m(m + 1) 2 for θ = (Σ) (38a) nd 2 (θ * , θ n ) ⇒ χ 2 m(m + 1) 2 + m for θ = (µ, Σ) (38b)
This provides a practical means of confirming the asymptotic normality of the estimators θ n . The function d 2 (•, •) denotes the square information distance, here the same as (28).

Global convergence analysis

This section studies the global convergence of the IDG and ISG algorithms, for two specific families of distributions, MGGD and Student-T. The main results are stated in the following two tables. For the cases indicated 1 and2, the cost function D(θ) (or D(θ)) has a unique stationary point, at θ * , which is the global minimizer. This will be obtained from the following development.

First, for the case of θ = (Σ) with known µ * and β * , let

f (δ x , β) = 1 g(δ x , β) (39)
then, for the MGGD model

f (δ x , β) = exp 1 2 δ β x β > 0 (40a)
and for the Student-T model,

f (δ x , β) = 1 + δ x β β+m 2 β > -m (40b)
The following proposition introduces a sufficient condition for the KL divergence D(Σ) and its empirical approximation D(Σ) to be geodesically strictly convex.

Proposition 5. assume that the function f : R + → R + in (39) verifies the following condition : for any ϕ :

R → R + ϕ strictly log-convex ⇒ f • ϕ strictly log-convex (41)
Then, the KL divergence D(Σ) (and its approximation D) is geodesically strictly convex.

In particular, the unique global minimum, and the unique stationary point, of D(Σ) is at the true Σ * . This proposition 5 directly yields the following corollary, for the specific MGGD model and Student-T model, by plugging (40a) and (40b) into (41).

Corollary 2. the KL divergence D(Σ) and its empirical approximation D(θ)

are geodesically strictly convex, with unique global minimum (and unique stationary point), in both of the following cases.

1. X is distributed according to an MGGD model, with scatter matrix Σ * and with shape parameter β > 0.

2. X is distributed according to a Student-T model, with scatter matrix Σ * and degree of freedom β > -m.

Thus, when Σ is unknown and β satisfies the conditions of Corollary 2, this corollary implies the global convergence of Algorithms 1 and 2. Precisely, these algorithms will always converge to the true value θ * of the parameter θ.

For the more complicated situation θ = (µ, Σ), global convergence does not always hold. The cost function D(θ) is not geodesically convex, but may be reformulated, using a new matrix argument [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF]

S =   Σ + µµ † µ µ † 1   (42)
If the new random vector y is given by

y = x † , 1 † (43)
then the cost function can be reformulated as

D(θ) = - 1 2 log det(S) -log f (δ y ) (44)
where

δ y = y † S -1 y = (x -µ) † Σ -1 (x -µ) + 1 (45) and f (δ y ) = exp 1 2 (δ y -1) β for MGGD (46a) Then for Student-T is f (δ y ) = 1 - 1 β + δ y β β+m 2
for Student-T (46b)

In [START_REF] Laus | Multivariate myriad filters based on parameter estimation of student-t distributions[END_REF], the minimization of D(θ) was proven to be equivalent to the minimization of D(θ). Replacing the new function f into (41), the following corollary is obtained. 1. X is distributed according to an MGGD model, with expectation and scatter matrix (µ * , Σ * ) and with fixed shape parameter β > 0.

2. X is distributed according to a Student-T model, with expectation and scatter matrix (µ * , Σ * ) and with the fixed degree of freedom β > 0.

For these two cases, global convergence is then guaranteed.

Finally, for the most complicated case, θ = (µ, Σ, β), the cost function is always non-convex. Moreover, we have verified experimentally that it has multiple stationary points in Θ = R m × P m × R + . Therefore, the correct estimation can only be guaranteed when the initial value θ 0 is close enough to the global minimum θ * .

Computer experiments

This section presents a set of computer experiments, which confirm the theoretical results of Section 4, and provide a detailed comparison of the ISG and IDG estimation methods, with the already existing MM and FP. For every experiment, 1000 Monte Carlo trials were carried out. For each trial, the dataset X = {x 1 , • • • , x N } is independent and identically distributed, according to true parameters (µ * , Σ * , β * ). The dimension m of x n is taken equal to 10. The true µ * is randomly chosen from a multivariate normal distribution. The scatter Σ * is defined as Σ(i, j) = ρ |i-j| for i, j ∈ {1, m}, and ρ ∼ U(0.2, 0.8). The shape parameter β * is uniformly selected from the intervals [0.2, 5] for MGGD and for Student-T.

The first experiment confirms the super-linear rate of convergence of IDG, for a dataset, distributed according to the MGGD model, which contains N = 10 4 samples. The initial value θ 0 is defined as the MM estimate, using 10% of the entire dataset. Figure 1(a) presents the case of θ = (Σ) with known (µ * , β * ). The IDG method converges after only two iterations, and if the same accuracy needs to be achieved, the deterministic gradient method (not using the information gradient) requires at least 88 iterations. For the case of θ = (µ, Σ) with known (β * ), things are similar. and 2(b) (that is to say, when β * is fixed), the step-size coefficient a always equals 1, which satisfies the condition in 1. For the case of unknown β , the step-size coefficient a is taken much larger, in order to meet the conditions of proposition 3. In fact, here, a = 100.

For the case of θ = (Σ) and θ = (µ, Σ), Figures 3(a In the third experiment, we compare the efficiency of the IDG and ISG For MGGD, the MM was given in [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized gaussian distributions with an application to multicomponent texture discrimination[END_REF], and the FP method in [START_REF] Pascal | Parameter estimation for multivariate generalized gaussian distributions[END_REF].

In Figures 4(a However, the experimental results show that the mini-batch has no significant effect on the accuracy of ISG.

As for computational time, information gradient methods have a significant advantage. The computational time of the ISG algorithm is similar to

d 2 (θ * , θ n ) 0 d 2 (θ * , θ n ) → c 0 and ∇ θ D(θ n ) 0 and ∇ θ D(θ n ) 0 θ = (µ, Σ, β) 73% 27%
that of MM, and is significantly less than that of FP. Meanwhile, its accuracy is significantly better than that of MM. In most experiments, the accuracy of ISG is similar to, or even better than, FP. Although the computational time of IDG is greater than that of ISG, it is comparable to that of FP, while, in most cases, IDG can achieve the best accuracy, among the four estimation methods considered.

Application with real dataset

In addition to experimental simulations, we also applied our methods to real datasets.

Color transformation

The first application is to color transformation for image editing with MGGD models, which was investigated in [START_REF] Hristova | Transformation of the multivariate generalized gaussian distribution for image editing[END_REF]. Its goal is to replace the color distribution of the input image by that one of a target image. The main idea is to fit the input and the target distributions, with two different MGGD models. Then, the transformation between these two MGGDs is implemented by a linear Monge-Kantorovich transformation for Σ, and a stochastic transformation for β. Specifically, this conversion can be threedimensional (3D), for RGB images, or five-dimensional (5D), when spatial gradient-field information is included.

Starting with the 3D rgb case, Figure 7 presents the transformed images and some of their details. The detail (a1) clearly shows that the cloud 'drawn' by MM appears too green. Similarly, FP also presents a green appearance, in detail (a2). On the contrary, the two gradient methods, i.e. IDG and ISG methods, show pure white cloud color in (a3) and (a4). Note also the difference in the amount of blue in the shadows on the grass. Too much blue is mixed with the shadow, in MM's output detail (b1). In details (b2),(b3),(b4), the results of MLE methods lead to a more natural appearance.

From the point of view of the present work, the most interesting aspect of this application is in term of computational time. The recursive (online)

ISG method takes about 10 seconds for two images (input and output). In contrast, FP and IDG each require more than two hours. In other words, ISG has a decisive advantage, in terms of time consumption.

Then, gradient-field information was included, so the transformation came to involve 5D, which consist in three color components (of CIELAB) and two components of the image spatial gradient field (dx and dy). For this application, the shape parameter of the MGGD model was supposed to fixed. Figure 8 presents the four different implementations. It can be observed that the output of the three MLE methods is significantly better than that of MM. In the transformed result of MM, the hue is darker and greener. MLE results are better, since the frost on the grass is whiter and appears more natural, and the forest on the mountain in the image also appears darker. The two images in Figure 8 have more than 1.2 × 10 6 pixels (i.e. 1.2 × 10 6 samples).

The FP and IDG need more than 4 hours to run, on the these two images.

The ISG method needs only 21 seconds.

We also considered an application to full HD images. In this case, as demonstrated in Figure 9, the advantages of the ISG algorithm were signif- was not feasible to run FP and IDG, with the entire dataset. Rather, the estimation was done on subsets of the complete dataset. These two subsets have 4 × 10 5 samples, that are randomly taken from the original images. In the autumn leaves obtained using FP and IDG, the yellow color has obviously been smeared. ISG is more natural, in which the yellow color is more uniform, and it is closer to the style of the target image.

Classification

MGGD are also used for texture modeling [START_REF] Kwitt | Testing a multivariate model for wavelet coefficients[END_REF][START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized gaussian distributions with an application to multicomponent texture discrimination[END_REF]. Without going into the details of presently existing classification methods, we attempted to use an MGGD representation, in order to distinguish between different groups of textures. Three groups of textures are selected from the VisTex database [START_REF]Mit vision and modeling group, vision texture[END_REF], We have reason to believe that, in this scenario, ISG has achieved the same performance as FP and simultaneously it used less time.

E u(θ 1 , x), ∇ θ D(θ 1 ) = -∇ Σ D(µ 1 , Σ 0 , β 0 ) 2 E u(θ 2 , x), ∇ θ D(θ 2 ) = -∇ β D(µ 1 , Σ 1 , β 0 ) 2 and so on, for k ≥ 3. This shows that (B.1) is verified.

Appendix C. Proof of propositions 3 and 4

As for Proposition 2, this is an application of Remark 2 in [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a riemannian manifold[END_REF]. According to this remark, in order to obtain the mean-square rate and the asymptotic normality, it is enough to show the mean vector field X(θ) = E θ * [u(θ, x)] has an attractive stationary point at θ = θ * . Since u(θ, x) = ∇ θ (θ, x)

E θ * [u(θ; x)] =      ∇ µ D(θ) ∇ Σ D(θ) ∇ β D(θ)      (C.1)
The covariant derivative of this vector field at the point θ = θ * is equal to the Hessian H(θ * ), which is positive definite. Therefore, the results of Propositions 3 and 4 follow by Remark 2 in [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a riemannian manifold[END_REF].

Plugging p+1 i=1 w 2 i e -q i t + 1 into the reformulated f

f • δ y (t) = exp    1 2 p+1 i=1 w 2 i e -q i t β    (E.6)
This function is proved to be log-convex in equation (D.3). Therefore, condition (41) is verified since β > 0 for MGGD model.

2) For Student-T, plugging (46b) into (41),

log f • ϕ (t) = β + m 2 1 - 1 β + 1 β exp(ψ(t)) (E.7)
Therefore, condition (41) is verified since β > 0.

. 2 .

 2 Problem formulation Parameter estimation will be formulated as the problem of minimising the Kullback-Leibler divergence D(θ * ||θ), denoted D(θ) for short. That is to say, the estimator θ is sought which is the solution of the following minimisation problem θ = arg min θ∈Θ D(θ)

Proposition 1 .

 1 Assume the cost function[START_REF] Verdoolaege | On the geometry of multivariate generalized gaussian models[END_REF] has an isolated stationary point θ = θ * in some neighborhood S 0 ⊂ Θ, and that the estimates (θ k ) k≥0 remain within S 0 . Then, for the sequence (θ k ) k≥0 generated by Algorithm 1,lim k→∞ θ k = θ *Appendix A sketches a proof of this convergence. For the case θ = (Σ) or θ = (µ, Σ), near the true value θ * , the Hessian of the function D(θ)

2 . 2 . 2 +

 222 sign is simplified as positive in the algorithm. The compact and convex set S 0 is a neighborhood of θ * , in which the cost function D(θ) has an isolated stationary point θ = θ * . The following proposition 2 states the convergence of Algorithm Proposition Assume the function D(θ) has an isolated stationary point at θ = θ * in S 0 , and that the estimates (θ n ) n≥0 remain within S 0 . Then, lim θ n = θ * almost surely. The proof of this convergence is discussed in Appendix B. Note that S 0 admits a system of normal coordinates (θ i ; i = 1, • • • , d) with origin at θ * , where d is the dimension of the parameter space Θ, d = m(m+1) m + 1.

Table 2 :

 2 Convergence analysis: Student-T Student θ = (Σ) Globally for β > -m θ = (µ, Σ) Globally for β > 0 in Tables

Corollary 3 .

 3 the KL divergence D(µ, Σ) (and D(µ, Σ)) has a unique global minimum (and unique stationary point) at (µ * , Σ * ), in both of the following cases.

Figure 1 (Figure 1 :

 11 Figure 1: The superlinear convergence rate for IDG

  )

  ) and 3(b) confirm the chi-squared limit distribution in corollary[START_REF] Kelker | Distribution theory of spherical distributions and a locationscale parameter generalization[END_REF]. The samples x n being matrices of size m × m with m = 10, the dashed blue curve is the probability density of a chi-squared distribution with 55 and 65 degrees of freedom, for Figures3(a) and 3(b) respectively. The solid lines are the smoothed histograms of N d 2 (θ * , θN ) where N = 10 5 . These "estimated p.d.f." coincide very closely with the theoretical chi-squared probability density.

( a )Figure 2 :

 a2 Figure 2: Linearly convergence rate for ISG

Figure 3 :

 3 Figure 3: Validation by fitting a Chi-2

3

 3 ), 4(b) and 4(c), the x-axis denotes the size of the dataset, and the y-axis denotes the expectation of the square distance between θ * and the estimated θ. This expectation is approximated by the average of 10 Monte Carlo trials. For the cases θ = (Σ) and θ = (µ, Σ), the IDG and ISG algorithms show a better accuracy. When θ = (µ, Σ, β), the accuracy of the MLE method is still significantly better than MM, and the accuracies of IDG and FP coincide. However, the accuracy of ISG is not as good as as FP or IDG. This phenomenon may be explained theoretically. Indeed, when θ = (µ, Σ, β), the product metric does not coincide with the information metric of the ECD model, and this leads to a less efficient estimation. The fluctuations of the curves in Figure4(c) are quite significant. This means the variance of the final estimate θ is significant.

( a )Figure 4 :

 a4 Figure 4: Efficiency comparison

Figure 5 :

 5 Figure 5: Additional experiments

  icant. The result of MM failed to achieve the color of the autumn leaves in the target image, showing cyan instead of yellow. Since the input image and the target image have more than 4 × 10 6 pixels (that is 4 × 10 6 samples)

Figure 7 Figure 8 :

 78 Figure 7: 3D transformation

11 paintings, 18

 18 fabrics, and 11 terrains. Each texture is considered as an RGB 3-dimensional image, modeled by an MGGD, whose parameters θ = (µ, Σ, β) are estimated by two MLE methods, i.e. FP and ISG. Then, the scatter matrices Σ are normalized by their trace, i.e. M = p tr(Σ) Σ (in order to

Figure 9 :Figure 10 :

 910 Figure 9: 5D transformation of images that have large size

Table 1 :

 1 Convergence analysis: MGGD

Table 3 :

 3 Percentage of 'correct' estimates

	correct estimates incorrect estimates

Acknowledgement

At the end, we thank the support from the ANR MAR-GARITA under Grant ANR-17-ASTR-0015 for our works.

Let (θ k ) k≥0 be an infinite sequence generated by Algorithm 1. Recall the retraction R θ defined in [START_REF] Pascal | Parameter estimation for multivariate generalized gaussian distributions[END_REF]. Consider the sequence of tangent vectors (η k ) k≥0 where η k belongs to T θ k Θ, and

, and so on.

Then, the sequence (θ k ) is given as in Algorithm 1 of [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF],

with step-size t k chosen according to Armijo-Goldstein rule (note that t 0 = α µ , t 1 = α Σ , and t 2 = α β , etc.).

The sequence (θ k ) remains within the neighborhood S 0 of θ * . Without loss of generality, assume this neighborhood is compact. According to Corollary 4.3.2 in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], if the sequence (η k ) k≥0 is gradient-related,

Then, since θ * is the only stationary point of the cost function [START_REF] Verdoolaege | On the geometry of multivariate generalized gaussian models[END_REF] in S 0 , it follows that lim k→∞ θ k = θ * , as required. To show that the sequence (η k ) k≥0 is gradient-related, note that

and so on, for k ≥ 3. In other words, the scalar product between η k and ∇ θ D(θ k ) is always strictly negative. Therefore, the sequence (η k ) k≥0 is gradient-related.

Appendix B. Proof of proposition 2

The proof is a direct application of Remark 2, concerning Proposition 1, in [START_REF] Zhou | Fast, asymptotically efficient, recursive estimation in a riemannian manifold[END_REF]. According to this remark, if u(θ n , x) denotes the direction of descent, and if

then lim θ n = θ * almost surely. Here (compare to the proof of Proposition 1), the direction of descent is given by u(θ

) and so on. Therefore, the expectation in (B.1) is equal to

Appendix D. Proof of proposition 5

For the case of θ = (Σ), the geodesic convexity of the cost function D(θ)

(or of D(θ)) follows by proving -(θ; x) is geodesically strictly convex in θ = (Σ) for any x.

To do this, for any fixed x, denote g(θ) = -(θ; x). Recall that, geodesic curves on P m are of the form [START_REF] Pennec | A riemannian framework for tensor computing[END_REF] γ :R → P m t →A exp(tr)A † (D.1)

where exp denotes the matrix exponential map, A is an invertible matrix, and r is a diagonal matrix, both of same size as Σ. Then, g(θ) is geodesically convex if and only if the composition (g • γ)(t) is always a convex function with respect to t. Moreover, geodesic strict convexity is defined in exactly the same way. The composition (g • γ)(t) can be expressed

where

x has components u i , and r i are the diagonal elements of r. The function ϕ : R → R + is strictly log-convex, because it is the Laplace transform of a positive measure [START_REF] Shiryayev | Probability[END_REF]

where µ = m i=1 u 2 i δ r i , and δ r i is the Dirac measure concentrated at r i . Assume that the function f verifies Condition (41). Then, since ϕ is strictly log-convex, f •ϕ is strictly log-convex. Thus, the term log [(f • ϕ)(t)] of (D.2) is a strictly convex function of the real variable t. Since the term tr(r) t 2 of (D.2) amounts to an affine function of t, it is now clear that (g •γ)(t) is a strictly convex function of the real variable t, for any geodesic curve γ : R → P m . Finally, since x was chosen arbitrarily, -(θ; x) is geodesically strictly convex in θ = (Σ) for each x. Therefore, D(θ) and D(θ) are both geodesically strictly convex.

Appendix E. Proof of corollary 2 and 3

For the case of θ = (Σ), note that ϕ : R → R + is strictly log-convex if and only if ϕ(t) = exp(ψ(t)) where ψ : R → R is strictly convex. 1) plugging (40a) into (41),

Therefore, condition (41) is verified since β > 0.

2) plugging (40b) into (41),

Therefore, condition (41) is verified since β + m > 0.

For the case of θ = (µ, Σ), as mentioned above, the function f is reformulated. Then, the same strategy is applied for this reformulated f . 1) For MGGD, recall the geodesic curve for reformulated matrix S(t),

where exp denotes the matrix exponential map, B is an invertible matrix, and s is a diagonal matrix, both of same size as S.

δ y (t) = y † S -1 y = p+1 i=1 v 2 i e -s i t with v = B -1 y (E.4)

According to equation (45), we have δ y > 1. Therefore, ∃w ∈ R p+1 and ∃q ∈ (0, +∞) p+1 (e.g. w = (u, 0) and q = (r, 1) ) such that p+1 i=1 v 2 i e -s i t = p+1 i=1 w 2 i e -q i t + 1 (E.5)


A Gatera

J Belmans

S Boussa

F Davin

W De Cock

V De Florio

F Doucet

L Parez

F Pompon

A Ponton

D Vandeplassche

E Verhagen

Mol Sck Cen

Belgium F Bouly

E Froidefond

LPSC A Plaçais

M Ben Abdillah

C Joly

L Perrot

C Zhang

MYRRHA-MINERVA INJECTOR STATUS AND COMMISSIONING *

 at SCK CEN, Belgium, aims at coupling a 600 MeV proton accelerator to a subcritical fission core operating at a thermal power of 60 MW. The nominal proton beam for this ADS has an intensity of 4 mA and is delivered in a quasi-CW mode. MYRRHA's linac is designed to be fault tolerant thanks to redundancy implemented in parallel at low energy and serially in the superconducting linac. Phase 1 of the project, named MINERVA, will realise a 100 MeV, 4 mA superconducting linac with the mission of demonstrating the ADS requirements in terms of reliability and of fault tolerance. As part of the reliability optimisation program the integrated prototyping of the MINERVA injector is ongoing at SCK CEN in Louvainla-Neuve, Belgium. The injector test stand aims at testing sequentially all the elements composing the front-end of the injector. This contribution will highlight the beam dynamics choices in MINERVA's injector and their impact on ongoing commissioning activities.

INTRODUCTION

The MYRRHA Accelerator is a 600 MeV proton linear accelerator part of a flexible irradiation facility being created by SCK CEN in Belgium as a test-bed for transmutation and as a fast-spectrum facility for material and fuel developments. This ADS consists of a low energy section (injector or Linac front-end) characterized by low beam velocity and * Part of this work supported by the European Commission Framework Programme H2020, MYRTE project nr. 662186 † angelique.gatera@sckcen.be a highly modular medium-and-high energy section (a socalled main linac) connecting a series of independently controlled accelerating superconducting cavities. Figure 1 depicts the accelerator and its two sections.

When compared with other accelerators of the same class [START_REF] Vandeplassche | Accelerator Driven Systems[END_REF], a major distinctive character of this system is given by its considerably harder reliability requirements: in MYRRHA, beam trips of a duration exceeding 3 seconds are considered as a system failure leading to a plant shutdown. In terms of MTBF, this corresponds to requiring a mean time between consecutive failures of at least 250 hours during the system's operative cycle (90 days). A common strategy to meet reliability requirements is given by fault-tolerance, namely the systematic adoption of provisions to compensate for unexpected deviations of the system state originating inside or outside of the system's boundaries. Given the complexity of our linac system, two fault-tolerance design patterns have been selected in function of the characteristics of the section they shall be applied to:

• In the low-energy section, hot spares shall be used as a form of parallel redundancy: two equivalent 16. [START_REF] Bouly | Commissioning of the MYRRHA Low Energy Beam Transport Line and Space Charge Compensation Experiments[END_REF] MeV injectors with fast switching capabilities [START_REF] Froidefond | Study and Design of a Fast Switching Magnet for the MYRRHA Project[END_REF] shall be put in place.

• In the main linac section, the availability of similarly structured cavities shall be capitalized so as to realize a serial redundancy strategy, with the functions of failing cavities to be compensated by their four nearest neighbors [START_REF] Bouly | Method Development for Cavity Failure Compensation in a Superconducting Linac[END_REF].

The MYRRHA project is scheduled to proceed in stages, with a first accelerator only based on a single injector connected to a superconducting section and an extraction line to a set of target facilities. In the first MYRRHA stage, MIN-ERVA, the main linac shall be based on single-spoke cavities and reach 100 MeV. In a second stage, the accelerator shall be extended with double-spoke and elliptical cavities up to a full-fledged 600 MeV configuration. A final stage will see the accelerator connected to a subcritical nuclear reactor.

THE INJECTOR TEST STAND

The MINERVA injector accelerates the beam up to 16.6 MeV. Its front-end is composed of an ECR proton source, a 2.6 m long LEBT (low energy beam transport line) [START_REF] Bouly | Commissioning of the MYRRHA Low Energy Beam Transport Line and Space Charge Compensation Experiments[END_REF] and a four-rod RFQ accelerating the beam to

RFQ BEAM COMMISSIONING

The first stage of RFQ beam commissioning, including full power commissioning, was carried out in 2020. As presented in Ref. [START_REF] Gatera | Minerva (MYRRHA Phase 1) RFQ Beam Commissioning[END_REF], all the measured performances were according to expectations and no problems were encountered when pushing duty cycle to CW. The RFQ transmission at theoretical optimal power (110 kW or 44 kV) was 95 %, which could be increased to 98 % at 125 kW (48 kV). In Figure 3: RFQ transmission maps with RFQ output isocurrent curves when scanning the solenoids in the LEBT for residual gas (top picture) / injected argon (bottom picture) in the LEBT.

1.5 MeV [START_REF] Zhang | New Reference Design of the European ADS RFQ Accelerator for MYRRHA[END_REF]. After the RFQ comes a small matching section with quarter wave resonators, called MEBT1, followed by a series of normal conducting CH cavities [START_REF] Kümpel | Current Status of MYRRHA Cavities[END_REF].

The front-end and the first seven CH cavities are being installed and tested on a dedicated test stand in Louvain-la-Neuve, Belgium, on the Cyclotron Resources Center (CRC) site. Today we have installed and are currently testing up to MEBT1 (see Fig. 2). We started by recommissioning the LEBT that had already been installed and commissioned at LPSC in Grenoble, France [START_REF] Bouly | Commissioning of the MYRRHA Low Energy Beam Transport Line and Space Charge Compensation Experiments[END_REF]. The two main objectives of the recommissioning have been [START_REF] De Bruyn | The Belgian MYRRHA ADS project: recent developments and future perspectives[END_REF] the study of the beam matching from the LEBT to the RFQ and its transmission, and (2) the study of space charge compensation transients, achieved by observing in detail the effects of enhanced space charge compensation when injecting gas into the LEBT.

For a defined current at the source exit, one can demonstrate by scanning through combinations of the two solenoids, that RFQ beam transmission is better optimised when argon is injected into the LEBT. Transmission maps without and with argon are shown in Fig. 3 where it is clear that transmission at high current is improved. Previously [START_REF] Gatera | Minerva (MYRRHA Phase 1) RFQ Beam Commissioning[END_REF], we also measured that adding argon allowed smaller steady-state emittance to be achieved. But that said, gas addition improves neither the transmission at nominal current nor the current rise-time at the RFQ exit. Therefore, unless there is a need for high beam current (>5 mA), the added value of injecting argon is limited.

the last months, we found out that running the RFQ also improved the transmission performances. We can now reach 98.9 % transmission at 4 mA around our new nominal RF power: 120 kW (46.5 kV).

The second stage of commissioning, carried out from Mid-2021, had as main objective the measurement of beam energy at the RFQ exit by time of flight (ToF).

The Time of Flight Method

The bunch's time of flight between two beam phase pickup systems (BPM or phase probe) can be estimated based on the phase shift measured at both detectors assuming that the number of RF buckets between the two is well known and that the phase offset due to cables and electronics is well calibrated [START_REF] Bozzolan | BPM Time of Flight Measurements for Setting-up the RF Cavities of the CERN Linac4[END_REF]. If the bucket number is not known, a third pick-up is needed to determine it. The setup in Fig. 4 is composed of two BPM (at black cables location) developed for MYRRHA during the MYRTE project [START_REF] Abdillah | Development of a Lowbeta BPM for MYRTE Project[END_REF]. It is installed at the end of MEBT1, after the second quarter wave rebuncher (see Fig. 2). A third BPM placed between the two rebuncher cavities has been used for bucket count verification.

RFQ Energy Measurement

The ToF setup described in the previous section was used for RFQ energy measurement with both bunching cavities switched off and detuned.

The measured output energy shows small fluctuations with respect to the RF power level inside the RFQ as can be seen from Fig. 5. These fluctuations are believed to come from synchrotron oscillations due to input energy mismatch. Below 90 kW, the beam bunching is not complete and the phase measurement shows more and more statistical variations until the system becomes unable of locking the phase.

The measured reference energy at nominal power (120 kW) is 1.494±0.003 MeV. The uncertainty margin covers statistical uncertainties and phase and distance errors from both BPMS. The obtained result is fairly close to the expected 1.5 MeV and is fully in the acceptance range of the downstream cavities. It was previously demonstrated in equivalent installations that the oscillation in function of RF power is sensitive to the RFQ input energy [START_REF] Orduz | Beam Commissioning SPIRAL2[END_REF]. In an effort to minimise the oscillation, we attempted an input energy scan by changing the source platform voltage.

As can be seen from the results in Fig. 6, increasing the platform voltage from 30.0 kV to 30.4 kV improves slightly energy stability. No improvement was observed passed 30.4 kV. Being the energy stability at 30.0 kV already good (only ±4 keV variation), improving it to ±3 keV variation does not seem as an improvement worth pursuing. The final decision on source platform voltage optimisation will have to wait for emittance measurements. We have in the same measurement campaign tried to test RFQ energy stability both in time and for different input beam currents. The energy variation in time was better than 1 keV over an hour of measurement and across different measurements (for the same LEBT and RFQ setpoints). Our next step is to install diagnostics to measure the longitudinal emittance. This should also help to conclude on the source platform voltage choice.

MEBT1

MEBT1, for first medium energy transport line, is the first matching section right after the RFQ. Its main goal is to match the beam into the longitudinal acceptance of the CH section linac. It is composed of two quarter wave rebunchers (QWR1 and QWR2) and one quadrupole triplet as shown in Fig. 7. The full RF conditioning of QWR1 and QWR2 was done in less than a day each. For these cavities, we use the same type of solid state amplifiers as for the RFQ [START_REF] Vandeplassche | Integrated Prototyping in View of the 100 MeV Linac for Myrrha Phase 1[END_REF] and dedicated in-house cavity tuners and LLRF systems on RFSoC FPGA [14].

The tuning of MEBT1 was done in two steps. First, the triplet was fine tuned around its beam dynamics determined setpoints in order to optimize beam transmission through MEBT1. Then, the rebunchers were tuned in phase and amplitude using the same ToF setup as for energy measurement. They were tuned so as to match an energy variation pattern simulated by the Tracewin code [START_REF] Uriot | Status of TraceWin Code[END_REF] with very good agreement as shown in Fig. 8.

Based on the good agreement between measurements and simulations, we are confident that the beam has the right characteristics for entering the first accelerating CH cavity. We should soon receive diagnostics that will allow us to confirm this with transverse emittance and bunch shape measurements.

CONCLUSION

In the last months, we have achieved the measurement of RFQ beam energy and tuning of our first RF cavities using a ToF setup composed of two BPM. A third independent BPM was used for bunch count or energy range confirmation.

The measured RFQ energy is 1.494 MeV which is in the acceptance range. Added to alomst 99 % transmission at nominal current and good behavior in CW mode, the RFQ has so far performed very well through all the milestones.

We are currently preparing the transverse emittance and bunch shape test bench. The diagnostics to be used are expected on site by the end of the year and the measurements are planned for the first quarter of 2022.

Figure 1 :

 1 Figure 1: Aceelerator scheme for the 100 MeV (MINERVA) and 600 MeV (MYRRHA ADS) phases [2].

Figure 2 :

 2 Figure 2: The injector test stand for MINERVA's front-end linac from the ECR source to the 7th accelerating cavity (up to 5.9 MeV).

Figure 4 :

 4 Figure 4: Setup for ToF measurement between two BPM.

Figure 5 :

 5 Figure 5: RFQ output energy with respect to RF power inside the RFQ.

Figure 6 :

 6 Figure 6: Fluctuations due to synchrotron oscillation for different RFQ input energy defined by the source platform voltage.

  61 st ICFA ABDW on High-Intensity and High-Brightness Hadron Beams HB2021, Batavia, IL, USA JACoW Publishing ISBN: 978-3-95450-225-7 ISSN: 2673-5571 doi:10.18429/JACoW-HB2021-WEBC3 WEBC3 Content from this work may be used under the terms of the CC BY 3.0 licence (© 2021). Any distribution of this work must maintain attribution to the author(s), title of the work, publisher, and DOI

Figure 7 :

 7 Figure 7: MEBT1 composed of two quarter wave rebunchers and one quadrupole triplet as currently installed at the injector test stand.

Figure 8 :

 8 Figure 8: Measured relative phase versus energy gain before (tentative) and after tuning a quarter wave rebuncher to match Tracewin simulation.

Content from this work may be used under the terms of the CC BY 3.0 licence (© 2021). Any distribution of this work must maintain attribution to the author(s), title of the work, publisher, and DOI

Commissioning and Operations


D Romanin

M Calandra

A W Chin

Excitonic switching across a Z 2 topological phase transition: from Mott-Wannier to Frenkel excitons in organic materials

Controllable topological phase transitions in low dimensional materials allows for unconventional optical selection rules and exciton series. We show that the broad mixing of single-particle contributions to the exciton states occurring across a Z2 phase transition in 1D polyacene polymers leads to an extreme tunability of optical properties, with exciton dispersions ranging from gigantic bandwidth (≈ 1.5 eV) to practically zero. The acene length transversal to the periodic axis controls the exciton localization, allowing for a cross-over from Mott-Wannier to Frenkel-like excitons and large changes in optical gaps, singlet-triplet splittings and the orbital structure of the real-space exciton states. Our work opens appealing perspectives for the design of novel optolectronic devices.

I. INTRODUCTION

Band topology is now considered to be as fundamental as band structure [START_REF] Vanderbilt | Berry phases in electronic structure theory: electric polarization, orbital magnetization and topological insulators[END_REF][START_REF] Breunig | [END_REF], underlying emergent quasiparticle phenomena such as robust edge states and spin-chain physics in organic systems (e.g. graphene nanoribbons) [3][4][5]. Moreover, the impact of topological band structures on the physics of excitons (i.e. an electron-hole photo-excited state bound by Coulomb interaction) has attracted a lot of attention in the last decade [6][7][8][9][10][11][12][13][14][15][16][17]. As a matter of fact, it has been shown for low dimensional materials that band topology dramatically changes their photophysics, breaking conventional optical selection rules and leading to novel exciton series [18,19]. Thus, understanding and controlling the effects of topological phase transitions on the optical response of materials represents a new powerful tool for the design of novel optoelectronic devices.

Generally speaking, π-conjugated organic semiconductors form a highly promising class of materials for the development of cheap, flexible and non-toxic optoelectronic devices such as photovoltaics, sensors and solid state lighting [START_REF] Köhler | Electronic processes in organic semiconductors: An introduction[END_REF][START_REF] Bronstein | [END_REF]. Experimental evidence for a topological (Z 2 ) phase transition has recently been presented in a series of 1D polyacene polymers in which the monomer size and chain length could be systematically varied [22,23], representing a physical realization of the well known Su-Schrieffer-Heeger (SSH) model [START_REF] Vanderbilt | Berry phases in electronic structure theory: electric polarization, orbital magnetization and topological insulators[END_REF]. This transition is associated with a dramatic closing and reopening of the polymer quasiparticle gap, as the size of the monomer increases. Nevertheless, the effects of a topological phase transition on the photophysics of such low-dimensional organic systems has, to date, remained unexplored and, based on the discussion above, could be relevant for the design of future carbon-based optoelectronic devices.

Figure 1 shows the structure of the polyacene polymers * romanin@insp.jussieu.fr † m.calandrabuonaura@unitn.it ‡ chin@insp.jussieu.fr of Ref. [22], which are composed by acene monomers connected by ethynylene bridges. The number of aromatic rings in each monomer is labeled n. Experimental chain lengths for both anthracene (n = 3) and pentacene (n = 5) vary from a few to a few-hundred monomers in the bridge direction (one dimensional periodicity). The HOMO-LUMO gap for acene monomers decreases smoothly as n increases, tending to a finite limit as n → ∞ [24]. However, it was observed in Refs. [22,23] that the band gap in polyacene polymers on Au(111) almost closes completely when going from n = 3 (anthracene) to n = 5 (pentacene). STM images revealed this to be accompanied by a marked change in the bondlength alternation (BLA, defined as the difference be-tween the single and triple bond lengths in the ethynylene link) of the polymer bridges, as well as a transition from an aromatic to quinoidal configuration for the central carbon rings in the monomers (Fig. 1). Crucially, these observations were rationalized in Ref. [22,23] by relating the electronic structure of the polymers to an effective Su-Schreiffer-Heeger (SSH) model [START_REF] Vanderbilt | Berry phases in electronic structure theory: electric polarization, orbital magnetization and topological insulators[END_REF]. The rationale for this correspondence is summarized in the Supplemental Material [START_REF]model, structural properties, CAMB3LYP vs. evGW bandstructures, antiferromagnetic groundstate and exciton wavefunctions[END_REF] (SM, Sec. S2), but the key point is that this well-known model possesses a Z 2 topological phase transition as the ratio t 1 /t 2 between the intra-(HOMO-LUMO gap, t 1 ) and inter-monomer (bridge-mediated, t 2 ) electron hopping amplitudes is tuned. Changing n, and thus t 1 , allows the polymers to cross the transition. The topological nature of the n = 5 acene polymers was further confirmed by the detection of localized 'edge states' in sufficiently long-but-finite pentacene (n = 5) polymers, as well as band inversion (vide infra) of their conduction and valence states [22,23].

In this paper, by employing manybody perturbation theory, we investigate the excitonic properties across the Z 2 topological transition occurring in polyacene polymers. We show that the topological transition is accompanied by a switching from extended Mott-Wannier excitons to localized Frenkel excitations and exciton dispersions that can be tuned from gigantic bandwidth (≈ 1.5 eV) to practically zero, simply by controlling the acene length n (see Fig. 1). This extreme tunability of optical properties occurs whilst also maintaining the desirable characteristics (large binding energy, strong optical absorption, large singlet-triplet splittings) that make organics attractive for applications [START_REF] Ratner | [END_REF]. Even more strikingly, we also show that a defining feature of non-trivial band topology -band inversion -can be seen directly in the experimentally measurable real-space wave functions of the excitons [27].

In Sec. II we summarize the computational details both for the structural, electronic (Sec. II A) and optical (Sec. II B) properties. We then discuss in the detail in Sec. III A the electronic band structure in vacuum of the various polymers under study (i.e. n = 3, 5, 7) and in Sec. III B we also analyze the effect due to the screening of the environment. After that, in Sec. III B and Sec. III C we analize the effects of the topological phase transition on the exciton wavefunctions and their dispersion respectively. In Sec. III D we study the origin of a negative dispersive triplet close to and after the quantum critical point of the topological phase transition. Finally, in Sec. IV we report our conclusions.

II. COMPUTATIONAL METHODS

A. Structural and electronic properties

Density functional theory calculations on the groundstate geometry (i.e. unit cell and atomic positions) are performed in vacuum with localized gaussian basis, using CAMB3LYP hybrid xc functional [28,29] and CRYSTAL17 [30,31]. The CAMB3LYP hybrid functional includes long-range corrections to the B3LYP, which is particularly suited for carbon-based materials, and has been shown to provide accurate results for charge transfer excitations, excitations to Rydberg states, polarizability of long chains and vibrational properties comparable to those obtained via Quantum Montecarlo methods [32,33] for systems in vacuum. We have employed the triple-ξ-polarized Gaussian type basis set [34] with real space integration tolerances of 10-10-10-25-50 and an energy tolerance of 10 -12 Ha for the total energy convergence. Each polymer has its principal axis along the x direction, with ∼ 500 Å of vacuum along the non-periodic y and z directions in order to have a proper one-dimensional system. The Brillouin zone is sampled with a uniform mesh of 50 × 1 × 1 k-points. Relaxation of the atomic coordinates is performed by imposing the polymers to be flat, since they have been synthesized on a surface. The out-ofplane relaxation is discussed in the SM [START_REF]model, structural properties, CAMB3LYP vs. evGW bandstructures, antiferromagnetic groundstate and exciton wavefunctions[END_REF] (Sec. S3). Finally, we consider here a paramagnetic ground state. The possible occurrence of magnetic states at very low temperature [35,36] is discussed in the SM [START_REF]model, structural properties, CAMB3LYP vs. evGW bandstructures, antiferromagnetic groundstate and exciton wavefunctions[END_REF] (Sec. S5), which however does not hinder the results presented here.

B. Optical properties

Many-body corrections to the single-particle bandstructure are computed within the GW approximation using the Yambo [37,38] code and the plasmonpole approximation [39]. Eigenvalues and wavefunctions are initially computed with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation (xc) functional and Quantum ESPRESSO [40,41]. Each polymer has its principal axis along the x direction, with ∼ 10 Å of vacuum along the non-periodic y and z directions as well as a cut-off on the Coulomb potential [42] in order to have a proper one-dimensional system. In this case we have employed a norm-conserving pseudopotential with the kinetic energy cut-off set to 70 Ry and a threshold of 10 -10 Ry on the total energy.

We performed self-consistent GW calculations on eigenvalues only (evGW ) [43] for both the Green's function G and the screened electron-electron interaction W in the Plasmon-Pole approximation (PPA). Excitonic effects are then evaluated by solving the Bethe-Salpeter [44,45] equation (BSE) in the Tamm-Dancoff approximation on top of the evGW band structure (a summary of the theory is given in the SM [START_REF]model, structural properties, CAMB3LYP vs. evGW bandstructures, antiferromagnetic groundstate and exciton wavefunctions[END_REF], Sec. S1). The number of k-points has been fixed to 100 in order to converge the first two singlet excitonic peaks in the absorption spectrum. Finally, we want to specify that results for the excitons are obtained via a dynamical screening calculated using QP eigenvalues for the BSE equation. Detailed numerical parametrization of the BSE is reported in the SM [START_REF]model, structural properties, CAMB3LYP vs. evGW bandstructures, antiferromagnetic groundstate and exciton wavefunctions[END_REF] (Sec. S1).

III. RESULTS AND DISCUSSION

A. Electronic structure While our main focus will be on excitonic properties, it is the quasi-particle band structure that undergoes the topological transition, and understanding its dependence on n provides important insights for our later discussions. Indeed, our aim is to study theoretically the effects of the topological phase transition on the optical properties of polyacene polymers. It is moreover important to note that in order to keep our focus on the essential, intrinsic physics of the transition and its potential impact on excitonic states, we have carried out all of our simulations in vacuum, avoiding further complications due to a surrounding environment. The possible effects of a substrate on the electronic structure is discussed later. for for polyacenes polymers in vacuum as a function of n. All lengths are given in Å, while all energies are given in eV. Hybrid refers to the CAMB3LYP functional [28,29], evGW to the GW method with self-consistency on the eigenvalues and BS to the Bethe-Salpeter formalism.

n c BLA ∆E ∆E SI EB SI -TI (Hybrid) (Hybrid) (Hybrid) (evGW) (BS) (BS) (BS) 3 6
In Fig. 2 we report the electronic band-structure obtained for n = 3, 5, 7 at the evGW level, together with the HOMO/LUMO wavefunctions both at the center (k = 0) and at the border (k = π/c) of the Brillouin zone. Tab. I reports the result of the single-particle band-gap, obtained via DFT (CAMB3LYP) and DFT+many-body perturbation theory (self-consistent GW on-eigenvaluesonly, i.e. evGW) for n = 3, 5, 7 in vacuum. All polymers show a direct single-particle band-gap at the edge (k = π/c) of the Brillouin zone: the gaps decrease from n = 3 to n = 5, then increase for n = 7, supporting a picture of a topological phase transition (band closure) between poly-pentacene and poly-heptacene.

We also note that while CAMB3LYP and evGW give very similar gaps for n = 3, the band gap for evGW becomes smaller than the CAMB3LYP gap for n = 5, 7, suggesting increasing electronic correlations close to the transition: approximate functionals cannot fully take these into account, so more accurate many-body techniques yield larger corrections. The largest deviations occurs for n = 5 (≈ 39%), which is the closest to the topological transition (smallest QP gap). Comparison (Experimental [22]) (evGW) between the band structures computed via CAMB3LYP and evGW for n = 3 -7 are presented in the SM [START_REF]model, structural properties, CAMB3LYP vs. evGW bandstructures, antiferromagnetic groundstate and exciton wavefunctions[END_REF] (Sec. S4).

B. Screening due to environment

If we compare the evGW estimate of the single-particle band-gap with the experimental values there seems to be a substantial difference (Tab. II). However, our computations are performed in vacuum while from the experimental point of view these polymers are deposited on Au(111) surfaces [22], which acts as a polarization environment, screening long-range interactions and strongly modifying correlation energies [46]. Therefore, we are not taking into account the effects of the environment on the quasiparticle properties. Indeed it has been shown [46,47] in the case of anthracene and pentacene molecules that, when deposited on metallic surfaces, their quasi-particle band-gap was subjected to a renormalization between 50% and 80% with respect to the gas phase. If we apply such renormalization to our evGW values we find a quasiparticle band-gap between ∼ 2.66 eV and ∼ 1.67 eV for poly-anthracene (n = 3) and ∼ 1.22 eV and ∼ 0.77 eV for poly-pentacene (n = 5), which are closer to the experimental results. The final discrepancy with respect to the experimental bandgap might be related to anharmonic effects, which can further reduce the energy difference between the HOMO and LUMO.

This reasoning suggests that the properties of the substrates on which the polymers are grown or deposed could provide another powerful and external control handle. Indeed, the screening from the substrate could substantially reduce the effects of the strong, long-range interactions, allowing for additional modulation of excitonic properties and of the critical acene length n c at which the topological transition occurs: Cirera et al.'s placement of n = 5 in the non-trivial phase is likely related to the metallic substrate inducing a gap reduction.

Notice however that neither in Ref. [22] nor in this work anharmonic quantum fluctuations are taken into account: these can play an important role in geometry optimization as well as in the Helmoltz free energy [32,33,48], especially close to the quantum critical point. in reciprocal space).

It is clear that the v and c states have inverted orbital character when passing from n = 5 to n = 7. This band inversion is the defining feature of the topologically non-trivial phase of the periodic SSH model [START_REF] Vanderbilt | Berry phases in electronic structure theory: electric polarization, orbital magnetization and topological insulators[END_REF], and occurs smoothly across the band structure, i.e. the c (v) orbitals of n = 7 at k = 0 are built from monomer LUMO (HOMO) orbitals and evolve into the inverted order shown in Fig. 3a as k → k = π/c [START_REF]) monomer orbitals due to the π orbitals of the bridge are present in all cases[END_REF]. No band in-version is seen for n = 3 and n = 5. Thus, in vacuum, the Z 2 transition occurs for monomers lengths between n = 5 (trivial) and n = 7 (nontrivial). This is further supported by the behaviour of the BLA of the bridge atoms (see Tab I), which dramatically reduces from n = 5 → n = 7, indicating the switch from alternating single and triple bonds to a nearly cumulenic double-double bond resonance form. Complete structural characterization can be found in the SM [START_REF]model, structural properties, CAMB3LYP vs. evGW bandstructures, antiferromagnetic groundstate and exciton wavefunctions[END_REF] (Sec. S3). As the k-space distribution of e-h transitions in both S I states is confined to the band edges, the exciton wave functions closely resemble the conduction band states at k = π/c, and thus also show band inversion as n = 5 → n = 7. Indeed, the exciton wave function receives contributions from single-particle orbitals taking part to the excitation and, as a consequence of the band inversion at the border of the Brillouin zone, we can observe as well a corresponding modification of the orbital character in the exciton wavefunction in real space.

However, a more dramatic effect of topological band inversion is shown in Fig. 3b, which shows the third bright excition (S III ) for n = 7. Here, we observe a very broad distribution of e-h transitions in the k-space wave function, whose width is larger than the region over which the band inversion occurs. This leads to contributions from single particle transitions in both the inverted and noninverted regions of the c/v bands. As a consequence, we observe a mixture of HOMO/LUMO states in the exciton wavefunction. These potentially measureable orbital textures are always absent in the topologically-trivial n = 3 and n = 5 polymers as no band exchange happens in their electronic structure. Representations of all singlet and triplet exciton wavefunctions for n = 3, 5, 7 is given in the SM [START_REF]model, structural properties, CAMB3LYP vs. evGW bandstructures, antiferromagnetic groundstate and exciton wavefunctions[END_REF] (Sec. S6).

D. Exciton dispersions and singlet fission

Fig. 4 compares the exciton dispersions E λ (q) for n = 5 and n = 7 as a function of the center-of-mass momentum q. Singlet excitons are optically spin-allowed, and have been plotted to take into account their generalized optical strength: is the electric dipole operator (here, aligned along polymer repeat axis) and A λ c,v (k, q) is the amplitude of a given electron-hole transition in the excitonic (λ) wave function in k space.

f λ (q) = c,v,k A λ c,v (k, q) ⟨c, k + q| D |v, k⟩
For n = 5, we find a rather congested spectrum containing both bright and dark singlet excitons, as well as two triplet bands. Notably, the lowest bright singlet exciton (S I ) is extraordinarily dispersive and remains 'optically' active over more than 60% of its bandwidth (1.6 eV, c.f. ≈ 0.1 eV in monomer films [START_REF] Refaely-Abramson | [END_REF]). Note that only states within the optical light cone can actually be excited by light, but the finite optical transition dipole moments of larger-q excitons will still be important for processes such as inter-polymer förster energy transfer and near-field, nano optical experiments on intrachain photodynamics [START_REF] Davoody | Physics, Simulation, and Photonic Engineering of Photovoltaic Devices IV[END_REF][START_REF] Hanewinkel | [END_REF]. Triplet excitons are also highly dispersive, with band widths of ≈ 0.7 eV, c.f. ≈ 0.03 eV in films [START_REF] Refaely-Abramson | [END_REF]. On the other side of the topological transition (n = 7), the exciton dispersions become significantly flatter, matching the underlying band structure (Fig. 2). However, the bright singlets maintain their oscillator strengths over a much larger momentum range due to their Frenkel-like nature (small Bohr radius).

Acene-based films and crystals are also known to undergo the process of singlet fission (SF) [START_REF] Refaely-Abramson | [END_REF]53], which allows a singlet exciton to decay into a pair of spinentangled triplet excitons [54] and could be used to increase the efficiency of organic photovoltaics [55]. Here, we find that the strong dispersion of the n = 5 polymer permits -on energetic grounds -the decay of S I at finite-q into both T I T I and excited T I T II triplet pairs (see Fig. 4) [54]: although the q = 0 S I state lies below the threshold energy for SF (≈ 2 × E T I ), as shown on Fig. 4, its rapid dispersion creates a large, energetically alllowed phase space for SF at finite momentum. Finite-q SF has previously been shown to be an effective mechanism in pentacene (n = 5) thin films [START_REF] Refaely-Abramson | [END_REF], even though the joint phase space is much smaller, due to the narrow dispersion of singlet excitons in films [START_REF] Refaely-Abramson | [END_REF]. Strikingly, the n = 5 polymer also permits ( on energetic grounds) the decay of S I at large-q into an excited T I T II triplet pair (Fig. 4). Recently, Pandya et al. have shown that optical excitation of a (T I T I ) → (T I T n ) transition induces spatial separation of the tightly bound triplet pairs that form the lowest excited state of blue polydiacetylenes, suggesting that the new SF pathway available in the n = 5 polymer could also spontaneously generate 'free' triplets pairs that preserve their germinate spin entanglement across real space [54].

Another critical aspect of SF is the matrix element for the conversion of a singlet Frenkel-like exciton, which -in acenes -is thought to arise via a super-exchange mechanism involving a higher-lying 'charge transfer' singlet exciton. While determination of SF rates is beyond the scope of this work (see Ref. [START_REF] Refaely-Abramson | [END_REF]), we do point out that the matrix elements for SF typically involve (virtual) superexchange coupling via a higher-lying charge-transfer exciton [56]. Indeed, this matrix element takes the generic form:

χ SF = (M Sn S I M T T Sn )(∆E S I Sn ) -1 (2)
where M y x is a Coulombic transition amplitude between states x and y that are separated in energy by ∆E xy . As shown in Sec. S6, the S III state in polypentacene (i.e. n = 5) has strong charge-transfer character, corresponding to a large M T T S III , and strongly mixes with S I near the avoided crossing (∆E S I S III ≈ 0). Consequently, a large super-exchange χ SF could be also be expected above threshold, leading to concomitantly large SF rates.

We also note that, for n = 5, S III at q = 0 is already above threshold and is optically bright, permitting direct SF from this CT state [56].

On the other side of the topological transition (n = 7), the larger relative energy of the T I states reduces the phase space for finite-q SF, although the S III band remains above threshold for all momenta. We also note that the bright singlets maintain their oscillator strengths over a much larger momentum range for n = 7 with respect to n = 5, due to the smaller excitonic Bohr radius. However, the key change of relevance for SF is the effect of band inversion on the SF matrix elements: as SF matrix elements are constructed from the coherent summation of Coulombic scattering amplitudes of electrons and holes [START_REF] Refaely-Abramson | [END_REF], this 'texture' in the electron density could lead to strong interference effects, which -to the best of our knowledge -have not previously been identified or studied. Such topologically induced interference effects would also radically affect the matrix elements for any finite-momentum scattering processes, such as the exciton-phonon interaction and related properties, e.g. resonance Raman spectra.

E. Negatively dispersive excitons

Another interesting feature in Fig. 4 is the existence of negatively dispersing T II excitons. For n = 5 this negative dispersion is confined to a small region around q = 0, while for n = 7 it extends across the entire band, decreasing across the band by ≈ 40 meV (K B T ≈ 25 meV at 298K), so that phonon scattering should spontaneously relax these excitons into stable finite-momenta states. At low temperatures, this could significantly modify the mechanism of triplet transport. A weak negative dispersion of the T II in a pentacene crystal was also observed, though not analysed, in Ref. [START_REF] Refaely-Abramson | [END_REF].

With a view to understanding and controlling these anomalous dispersions, Fig. 5a compares the k-space distributions of e-h transitions in the T II state on the n = 7 quasiparticle band structure. At q = 0, the T II state is built by e-h pairs clustered around the band edge and zone center, and at q = π/c the dominant transitions lie in the centre of the band [START_REF]Concerning q = 0.5 (i.e. q = π/c, in units of 2π/c), in Fig. 5a light blue (dark blue) dots represent electrons excited from band v at states 0[END_REF]. According to the effective Bethe-Salpeter Hamiltonian, the energy of the triplet exciton states contains two contributions: the (non-interacting) kinetic energy of the e-h quasiparticles and the matrix element between the conduction c and valence v band of the Fourier transform of the screened direct Coulomb interaction K X cvk c ′ v ′ k ′ at k-points k and k ′ = k + q. Due to the flatness of the band structure for poly-heptacene, there is essentially no change in the carrier's kinetic energy as T II (q = 0) → T II (q = π/c), whereas this normally provides a strong, positive contribution to the exciton dispersion, as it does for polypentacene.

At the same time, the direct Coulomb interaction ma- trix element that mixes e-h transitions in the excitonic wave function arises from a screened q-dependent interaction. Fig. 5d shows the screened q-dependent direct interaction averaged over triplet states for the n = 7:

K X (q) = cvk c ′ v ′ k ′ (A λ c ′ ,v ′ (k ′ , q)) * A λ c,v (k, q)K X cvk c ′ v ′ k ′ (3)
showing that K X (q) is a continuously decreasing function of q for T II , i.e. the binding energy of the exciton becomes greater at the band edge. Combined with the neutral change in carrier Kinetic energy (K.E.) -which is induced by the flat bands that appear above the topological transition -we then see that the negative dispersion arises from reduced screening at short wavelengths (large q). This interpretation is visualized directly in Fig. 5b&c, where the shrinkage (increased binding energy/decreasing Bohr radius) of the exciton can be seen in the real-space wave function as q → π/c. Thus, it is the flat bands that must appear above the topological transition (the disconnected dimer limit, in the SSH model) that promote anomalous exciton dynamics. We also note that negative dispersions in molecular acene crystals have previously been reported for low-lying singlet excitons, and were shown to arise from mixing (level repulsion) with higher-lying CT excitons [START_REF] Roth | [END_REF]59]. While the congested spectrum for n = 5 may partially promote a similar mechanism in the non-monotonic dispersion of T II and the rather flat S III , the CT bands are well separated for n = 7, and the shrinkage of the wave function at the band edge is inconsistent with growing admixtures of CT states with enhanced e-h separations.

IV. CONCLUSION

In conclusion, we have shown that the topological transition occurring in ethynylene-bridged acene polymers of-fers extreme tunability of optoelectronic properties, and we believe that the rather general mapping between 1D molecule-bridge structures and the SSH model should extend this tunability to many other organic constructs [22]. Additionally, the properties of the substrates on which the polymers are grown or deposed could provide another powerful and external control handle. Indeed, the screening from the substrate could substantially reduce the effects of the strong, long-range vacuum interactions native to organics, allowing for additional modulation of excitonic properties and -most importantly -of the critical acene length n c at which the topological transition occurs. Additionally, based on the findings of previous work on carbynes [32] and finite-sized acene polymers [23], the cross-over from acetylenic to cumulenic bridging at the transition suggests that phonons and anharmonicity could also play an important, temperaturedependent role in closing the charge gap, and this will be discussed in future work.

FIG. 1 .

 1 FIG. 1. Topological phase transition of polyacene polymers: as the acene length n increases, the bridge configuration switches from ethylynic (-C ≡ C-) to cumulenic (= C = C =), while the central ring goes from an aromatic to a quinoidal form. Above a critical length nc, the topologically non-trivial phase (see text) emerges, showing a smooth change of orbital character (red,blue) across the Brillouin zone (band inversion).

kFIG. 3

 3 FIG.2. Self-consistent GW (evGW) quasi-particle bandstucture and HOMOs/LUMOs orbital wavefunctions (with relative phases) at the center (k = 0 and at the border (k = 0.5) of the Brillouin zone for ethynylene-bridged poly-anthracene (n=3), poly-pentacene (n=5) and poly-heptacene (n=7). Grey (white) spheres represent carbon (hydrogen) atoms.

  The excitation and binding energies of the lowest (optically bright) singlet (S I ) and triplet (T I ) excitons are given in Tab I for n = 3, 5, 7 . The excitonic gaps follow the same trend as the band structure, as do the large excitonic binding energies (∼ 0.3 -0.6 eV). Strong 1D coulomb interactions generate large singlet-triplet exchange splittings (∼ 0.3/0.6 eV) in all polymers (Tab. I). Strikingly, the optical gap for n = 5 is dramatically reduced (0.77 eV, c.f. 1.8 eV in molecular crystals), and the strong increase in the S I binding energies from n = 5 → n = 7 (Tab. I) implies a transition from Mott-Wannier to Frenkel excitons. Indeed the Bohr's radius of S I goes from ∼ 76 Å (i.e. ∼ 11 unit cells) to ∼ 28 Å (i.e. ∼ 4 unit cells) when passing from n = 5 to n = 7. The change in localization of the excitonic wave function is clearly seen in both real and reciprocal space in Fig. 3a (right panel).

2 ( 1 )

 21 FIG. 4. Exciton dispersions for n = 5 (a) and n = 7 (b) showing the energy threshold for triplet pair creation (horizontal lines). Singlet exciton dispersions are coloured according to their normalized optical oscillator strength (color bar), with the brightest state having unit strength.

FIG. 5 .

 5 FIG.5. Poly-heptacene (n=7): (a) single-particle transitions responsible for TII at q = 0 (red) and q = π/c (blue) on top of the folded band structure in the irreducible first Brillouin zone; (b-c) Modulus squared of the exciton wave functions for TII at q = 0 (b) and q = π/c (c) showing contraction of the exciton at finite momentum. Holes are localized at circled position; (d) Screened direct Coulomb interaction K X averaged over triplet eigenfunctions as a function of q, obtained with Yambopy.

ACKNOWLEDGMENTS

The authors acknowledge support from the ANR project ACCEPT (Grant No. ANR-19-CE24-0028). This work was granted access to the HPC resources of IDRIS, CINES and TGCC under the allocation 2021-A0100912417 made by GENCI. MC acknowledges PRACE (Grant no. 2021240020) for awarding him access to irene Rome in France. D. R. thanks D. Varsano and F. Paleari for fruitful dicussion on the GW+BSE computations. The authors acknowledge M. Casula for fruitful discussion on the possible insurgence of an antiferromagnetic groundstate. The data that supports the findings of this study are available from the corresponding author on reasonable request.