Claire Cury
email: claire.cury.pro@gmail.com

Stanley Durrleman

David M Cash

Marco Lorenzi

Jennifer M Nicholas

Martina Bocchetta

John C Van Swieten

Barbara Borroni

Daniela Galimberti

Mario Masellis

Maria Carmela Tartaglia

James B Rowe

Caroline Graff

Fabrizio Tagliavini

Giovanni B Frisoni

Robert Laforce

Elizabeth Finger

Alexandre De Mendonça

Sandro Sorbi

Sebastien Ourselin

Jonathan D Rohrer

Marc Modat

Spatiotemporal analysis for detection of pre-symptomatic shape changes in neurodegenerative diseases: Initial application to the GENFI cohort

Keywords: Shape analysis, clustering, Computational anatomy, thalamus, spatiotemporal geodesic regression, parallel transport

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

Introduction

Neurodegenerative diseases such as frontotemporal dementia (FTD) present progressive symptoms of behavioural and cognitive dysfunction. These changes follow many years of a clinically silent phase in the disease, where abnormal proteins slowly accumulates within the brain, leading to neurodegenerative processes that ultimately result in loss of function. Reliably identifying presymptomatic changes in individuals could lead to intervention with therapies that could slow, or even halt, the onset of these diseases. However, finding a cohort of presymptomatic individuals guaranteed to develop a form of dementia can be challenging. One common strategy is to investigate people who are at-risk for rare autosomal dominant forms of dementia. Half of these individuals are carriers of the mutation, allowing for comparisons between carriers and non-carriers at various stages within the disease process. In the case of genetic FTD, roughly one third of all cases are caused by autosomal dominant mutations, primarily in three genes: chromosome 9 open reading frame 72 (C9orf72 ), progranulin (GRN ), and microtubule associated protein tau (MAPT ) [START_REF] Rohrer | Phenotypic signatures of genetic frontotemporal dementia[END_REF]. As the name would suggest, in all mutations, there is early involvement of both the frontal and temporal lobes, as well as the insula where differences can be observed as early as ten years before estimated age of expected symptom onset, as shown in Rohrer et al. [START_REF] Rohrer | Presymptomatic cognitive and neuroanatomical changes in genetic frontotemporal dementia in the Genetic Frontotemporal dementia Initiative(GENFI) study: a cross-sectional analysis[END_REF]. However, there are additional structures, such as the thalamus, which also appear to be implicated to some degree early on in the disease process [START_REF] Bocchetta | Thalamic atrophy in frontotemporal dementia -Not just a C9orf72 problem[END_REF]. In many forms of FTD, clinical presentations suggest a left/right asymmetry in terms of which hemisphere is more affected, and this is often supported by evidence of increased atrophy within the affected hemisphere [START_REF] Boccardi | Clinical characteristics of frontotemporal patients with symmetric brain atrophy[END_REF]. However, the affected side is not consistent across all cases, and in some cases, there is no evidence of an asymmetry.

As this asymmetry is likely to start early in the disease process, it must be taken into account when looking to detect early changes with any sensitivity.

One biomarker that shows promise during the presymptomatic phase is measurement of atrophy derived from structural magnetic resonance imaging (MRI) [START_REF] Benzinger | Regional variability of imaging biomarkers in autosomal dominant Alzheimer's disease[END_REF][START_REF] Rohrer | Presymptomatic cognitive and neuroanatomical changes in genetic frontotemporal dementia in the Genetic Frontotemporal dementia Initiative(GENFI) study: a cross-sectional analysis[END_REF][START_REF] Schott | Increased brain atrophy rates in cognitively normal older adults with low cerebrospinal fluid Aβ1-42[END_REF]. Volumes summarizing change within a region of interest (ROI) tend to be more sensitive to early change than voxelwise approaches, but they do not provide any spatial localisation as to where the atrophy is occurring within the ROI. Conversely, voxelwise analysis can provide better spatial localisation, but the mass univariate nature of the analysis requires correction for multiple comparisons to control for false positive findings, which often results in reduced sensitivity. As loss of brain volume will imply a change in the shape of the structure, a third option is to perform the shape analysis over time for a structure of interest. This could provide more spatial information than a single summary measure of volume alone, but does not require the same level of multiple comparisons as a voxelwise analyses. Given the decades long nature of the disease process, it is not yet feasible to measure the complete time course within one individual. Therefore, the pattern of atrophy over the course of the disease must be estimated through spatiotemporal regression models based on large populations of either cross-sectional data or through longitudinal data that covers a smaller segment (i.e. a few years) of the disease process within each individual.

There have been numerous approaches to spatiotemporally model trajectories for ageing and dementia. Some methods model this evolution using dense 4D deformation fields to measure change between timepoints. Lorenzi et al. [START_REF] Lorenzi | 4D registration of serial brain's MR images: a robust measure of changes applied to Alzheimer's disease[END_REF] modelled the 4D deformation fields within a population to obtain subject-specific measurements of atrophy. An extension of this work discriminated spatiotemporal patterns that could be attributed to natural ageing versus those that were related to disease [START_REF] Lorenzi | Disentangling normal aging from Alzheimer's disease in structural magnetic resonance images[END_REF]. Other groups establish point correspondences between subjects on a surface representation, and then apply mixed effects models at those points [START_REF] Datar | Mixed-Effects Shape Models for Estimating Longitudinal Changes in Anatomy[END_REF][START_REF] Muralidharan | Diffeomorphic Shape Trajectories for Improved Longitudinal Segmentation and Statistics[END_REF][START_REF] Younes | Inferring changepoint times of medial temporal lobe morphometric change in preclinical alzheimer's disease[END_REF], providing fixed effects that represent the change across the overall population while allowing individual longitudinal trajectories as random effects. More complex representations of surfaces can be used, as in Durrleman et al. [START_REF] Durrleman | Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape datas[END_REF], they proposed a spatiotemporal regression approach to estimate continuous subject-specific trajectories of longitudinal data.

In our previous work [START_REF] Cury | Spatio-Temporal Shape Analysis of Cross-Sectional Data for Detection of Early Changes in Neurodegenerative Disease[END_REF], we defined the shape of the structure of interest as its 3D outline that is rotation and translation invariant. Differences between shapes were quantified using the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework [START_REF] Trouvé | Diffeomorphisms groups and pattern matching in image analysis[END_REF][START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Glaunès | Large Deformation Diffeomorphic Metric Curve Mapping[END_REF], producing a smooth and invertible continuum between all possible shapes within the population. The smooth representation of these deformations also acted as low-pass filter, reducing the effects of irregularities and errors in the surface boundaries. Overall, our approach consisted of three main steps. First, using all available data, we compute an average shape spatiotemporal trajectory. Second, for every individual shape we evaluate its distance from the mean trajectory. Last, after spatially normalising all the subject-specific distances to the mean, we run a statistical analysis on the subject-specific residuals to assess when a shape starts diverging from normality. This previous work presented a global spatio-temporal analysis, on one side of the brain, without considering a potential left/right asymmetry of the disease. In this paper, we build on the aforementioned framework, which we altered in two main ways.

First, we take into consideration the potential asymmetry of FTD by considering the left and right structures using a common shape representation. Second, we modified our feature extraction method using a clustering approach to ensure we can attribute the recovered differences to substructure of the shape under study, and made a novel local analysis, based on clustering of deformations, which takes better advantage of the LDDMM framework.

We apply this approach to data from the Genetic FTD Initiative (GENFI), an international study of autosomal dominant forms of FTD aimed at collecting multimodal neuroimaging, alongside other biomarkers with the objective of obtaining an improved understanding of the changes that are occurring during the presymptomatic phase of the disease. In general, the expected age of onset of clinical symptoms is estimated by using the average age of onset in the family of the subject, allowing to align the different subjects onto a single time axis. We applied our method to a subcortical structure, the thalamus, which has been shown to present volumetric differences before onset in Rohrer et al. [START_REF] Rohrer | Presymptomatic cognitive and neuroanatomical changes in genetic frontotemporal dementia in the Genetic Frontotemporal dementia Initiative(GENFI) study: a cross-sectional analysis[END_REF]. We used the expected age to onset to characterise the time progression. In the next section, we will present the different steps of the proposed framework before then further describing the experiment and associated results.

Method

We indicate with {(S i , t i )} i∈{0;...;N -1} a set of N shapes associated with a corresponding time point t i . With analogy to classical random-effect-modelling approaches, we assume that each shape is a random realisation of a common underlying spatiotemporal process φ(t):

S i = ρ i (φ(B 0 , t i )) + i ,
where B 0 is a common reference frame, and ρ i is a subject-specific "residual" deformation accounting for individual deviation from the mean shape. We characterise this residual through the diffeomorphism linking the shape S i to the corresponding sample of the common spatiotemporal trajectory at time point t i . We also assume that i is Gaussian randomly distributed noise. In order to identify group-wise differences between the given populations, we rely on the analysis of the subjects-specific residuals deformations

ρ i .
This is a challenging problem, since all ρ i are defined at different time points along the common spatiotemporal trajectory, and therefore cannot be directly compared in a common anatomical framework. Moreover, the optimisation of the functional for the simultaneous estimation of the group-wise trajectory and random effects is not trivial, and would ultimately result in expensive and thus impractical numerical schemes. For these reasons, we propose a serial optimisation of the problem by introducing an efficient numerical framework composed of three steps illustrated in Figure 1.

(i) First, we assume that the residuals deformations ρ i are fixed, and we estimate the common trajectory φ(t). (ii) Second, given the modelled trajectory φ, we estimate the residuals deformations ρ i through non-linear registration between the trajectory point φ(B 0 , t i ) and S i . (iii) Third, we spatially normalise the residual deformations in the common initial reference space B 0 using parallel transport.

The proposed framework relies on the mathematical setting of the Large Diffeomorphic Deformation Metric Mapping (LDDMM) framework and the varifold representation They have to be transported to a common space (i.e. B 0 ) along the geodesic φ, so they can be analysed (step 3).

of shapes (section 2.1). This choice allows a mathematically consistent definition of all steps (section 2.2), namely: (i) the spatiotemporal regression, (ii) the ρ i deformations estimation, and (iii) the normalisation of the initial momentum of ρ i through parallel transport.

Large diffeomorphic deformation metric mapping and varifold representation

The LDDMM framework [START_REF] Trouvé | Diffeomorphisms groups and pattern matching in image analysis[END_REF][START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF] is a mathematical and algorithmic framework based on flows of diffeomorphisms, which allows comparing anatomical shapes as well as performing statistics. The framework used in this paper is a discrete parametrisation of the LDDMM framework, as proposed by Durrleman et al. [START_REF] Durrleman | Optimal Data-Driven Sparse Parameterization of Diffeomorphisms for Population Analysis[END_REF], based on a finite set of N B0 control points overlaid on the 3D space enclosing the initial shape B 0 . The control points number and position are independent from the shapes being deformed as they do not require to be aligned with the shapes' vertices. They are used to define a potentially infinite-dimensional basis for the parametrization of the deformation. Momentum vectors are associated with the control points and are used as weights for the decomposition of a given deformation onto this basis.

Deformation maps ϕ v : R 3 → R 3 are built by integrating time-varying vector fields (v t ) 0≤t≤1 , such that each v(•, t) belongs to a Reproducing Kernel Hilbert Space (RKHS) V with kernel K V . We use a Gaussian kernel for all control points x, y:

k V (x, y) = exp -|x -y| 2 λ 2 Id,
with Id the identity matrix, and λ a scale factor which determines the size of the kernel and therefore the degree of smoothness of the deformations. We define ϕ v (x) = φ v (x, 1) as the diffeomorphism induced by v(x, t) where φ v (x, 1) is the unique solution of the differential equation:

dφ v dt (x, t) = v(φ v (x, t), t), ∀t ∈ [0, 1] with φ v (x, 0) = x, ∀x ∈ R 3 .
Velocity fields v(•, t) are controlled via an energy functional

1 0 v(•, t) 2 V dt
, where • V is a Hilbert norm defined on vector fields of R 3 , which is used as a regularity term in the matching functional to penalise non-regularity. In the LDDMM framework, matching two shapes S and T requires estimating an optimal deformation map φ : R 3 → R 3 such that φ(S) is close to T . This is achieved by optimising

d([ϕ v (S)], [T ]) 2 + γ 1 0 v(•, t) 2 V dt,
where γ balances the regularity of φ v against the spatial proximity d, a similarity measure between the varifold representation of ϕ v (S) and T noted respectively [ϕ v (S)] and [T ].

In a discrete setting, the vector fields v(x, t) corresponding to optimal maps are expressed as combinations of spline parametrised fields that involve the reproducing kernel K V of the space V :

v(x, t) = N B 0 p=1 K V (x, x p (t))α p (t),
where x p (t) = φ v (x p , t) are the trajectories of control points x p . The control points are regularly spaced on a 3D grid overlaid on the space that contains the mesh of the subject S. The control point spacing is defined by the size of the kernel K V . The time-dependent vectors α p (t) ∈ R 3 are referred to as momentum vectors attached to x p .

The full deformation can be encoded by the set of initial momentum vectors α(0) = {α p (0)} 1≤p≤n located at the points {x p } 1≤p≤n . This allows the analysis of a set of deformation maps from a given template to the observed shapes by performing statistics on the initial momentum vectors defined on control points located around the template shape. The process of generating back any deformation map from initial conditions (x p (0), α p (0)), i.e. integrating the geodesic equations, is called geodesic shooting or exponential map and is noted exp xp(0) (α p (0)).

As previously stated, varifolds are used to represent shapes [START_REF] Charon | The Varifold Representation of Nonoriented Shapes for Diffeomorphic Registration[END_REF]. They are nonoriented versions of the representation with currents [START_REF] Vaillant | Surface matching via currents[END_REF], which are used to efficiently model a large range of shapes. To represent a shape S as a varifold, the shape space is embedded into the dual space of a RKHS W, noted W * , and encoded using a set of non-oriented unit normals attached on each vertex of the shape. This kernel-based embedding allows to define a distance between different embedded shapes. Varifolds are robust to varying topologies, do not require point to point correspondences, and embed the shapes in a vector space, which facilitate the interpretation of results. The varifold representation of a discretised mesh composed by M triangles S is noted [S] and writes:

[S](ω) = M k=1 ω(c k )τ (c k ) 2 / τ (c k )
with ω a vector field in W, c k the centre of the triangle k, and τ (c k ) the tangent of the surface S at point c k .

Residual extraction framework

Due to the asymmetry of the disease, the proposed framework has been designed so that it is unbiased to the affected side. For each subject, rather than considering the left or right structure, we build a mean shape by averaging both sides. First, the structure of interest is segmented using the method proposed by Cardoso et al. [START_REF] Cardoso | STEPS: Similarity and truth estimation for propagated segmentations and its application to hippocampal segmentation and brain parcelation[END_REF]. Second, we flip all input T1w brain images and segmentation masks, in order to have all structures, left and right, on the same side. Third, we affinely align the T1w brain images (the originals and the flipped ones) to a subject-specific mid-space [START_REF] Modat | Global image registration using a symmetric block-matching approach[END_REF]. The MNI52 atlas was used to define the mid-space and ensure that all subjects have a similar total intracranial volume (TIV). TIV varies from subject to subject due to normal variability in the population. Alignment to a common mid-space enables to discard this inter-subject variability through normalisation. The obtained affine transformations are then applied to the corresponding segmentation masks. Fourth, we compute a mask covering the area of the structure of interest and its surroundings for all T1w MRI, to estimate a rigid refinement focused on the area of interest. This is achieved by a 6-voxel dilation of the union of all propagated masks to ensure that for each subject, the structure of interest and its surrounding are considered. The rigid refinement step is done using the T1w MRIs rather than the segmented shape. Finally, we extract the meshes of the left (flipped, L i ) and right structures (R i ), and compute the mean shape, by estimating the diffeomorphisms χ (i) v for each subject i, such as χ

(i) v = argmin 1 2 ( [χ vi (L i )] -[S i ] 2 W * + [χ vi (R i )] -[S i ] 2 W * ) + γ 1 0 v i (•, t) 2
V dt with W * the space of varifolds and S i , the obtained subject-specific average shape of the structure of interest, is associated with a temporal information t i , the number of years to the expected onset (EYO) of the subject i.

The computation of the spatiotemporal regression [START_REF] Durrleman | Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape datas[END_REF] requires an initial shape B 0 = {x p } p=1,...,N B 0 as reference. To avoid any bias towards a subject selected as the initial shape, we estimate the initial shape from the 10 subjects who are the furthest away from expected symptom onset, who are all approximately 40 years before their expected onset of clinical symptoms according to EYO. We estimate the centroid of those 10 subjects using the diffeomorphic Iterative Centroid method [START_REF] Cury | Diffeomorphic iterative centroid methods for template estimation on large datasets[END_REF], which estimates a centre of a given population in a reasonable computation time [START_REF] Cury | the Imagen Consortium, Statistical Shape Analysis of Large Datasets Based on Diffeomorphic Iterative Centroids[END_REF].

The spatiotemporal regression of the set of shapes {(S i , t i )} i∈{0;...;N -1} is implemented in the Deformetrica software [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF][START_REF] Routier | Evaluation of morphometric descriptors of deep brain structures for the automatic classification of patients with Alzheimer's disease, mild cognitive impairment and elderly controls[END_REF] 2 . The EYO values are discretised into T time points. Starting from B 0 at time t = 0, a geodesic moving through the positions φ(B 0 , t), ∀t ∈ {0; ...; T } is computed by minimising the discrepancy between the model at time t (i.e. φ(B 0 , t)) and the observed shapes S i :

E(φ v ) = ti d([φ v (B 0 , t i )], [S i ]) 2 + γ v 2 V φ ,
with v the time-varying velocity vector field that belongs to the RKHS V determined by the Gaussian Kernel K V . The initial momentum vectors β 0 (0) = {β 0 p (0)} 1≤p≤N B 0 are defined on the control points grid overlaid on the baseline shape B 0 and fully encode the geodesic regression parametrised by {B 0 ; β 0 (0)}.

We then compute the residuals diffeomorphic deformations ρ i between every observation and the spatio-temporal average shape by estimating a geodesic between φ(B 0 , t i ) and {S i , t i }. This yields a set of trajectories parametrised by {φ(B 0 , t i ); α i (0)} i∈{0;...;N -1} that encode the deformations ρ i from the spatio-temporal regression to all subjects, with α i (0) the initial momentum vectors, where the varying parameter is the step of the deformation. This parameter should not be confused with the time variable corresponding to the EYO and to the time varying deformation of the main spatio-temporal trajectory.

In order to be able to compare this set of momenta, we gather them in the same Euclidean space. This is achieved by transporting all momenta into the initial space of B 0 = φ(B 0 , 0), using a parallel transport method based on Jacobi fields as introduced in [START_REF] Younes | Jacobi fields in groups of diffeomorphisms and applications[END_REF]. Parallel transporting a vector along a curve (the computed trajectory parametrised by (B 0 ; β 0 (0))) consists in translating it across the tangent spaces along the curve by preserving its parallelism, according to a given connection. The Levi-Civita connection is used in the LDDMM framework. The vector is parallel transported along the curve if the connection is null for all steps along the curve [START_REF] Lorenzi | Efficient Parallel Transport of Deformations in Time Series of Images: From Schild's to Pole Ladder[END_REF]. We use Jacobi field instead of the Schild's Ladder method [START_REF] Kheyfets | Schild's Ladder Parallel Transport Procedure for an Arbitrary Connection[END_REF], to avoid the cumulative errors and the excessive computation time due to the computation of Riemannian Logarithms in the LDDMM framework, required for the Schild's Ladder. The cumulative errors would have differed from subject to subject and thus introduce a bias. Indeed, their distances from the baseline shape vary, as they all are at different points along the temporal axis. The Jacobi field, used to transport a vector α i (0) from a time t to the time t 0 = 0 along the geodesic γ, is defined as:

J γ(t) (0, -β 0 (t), α i (0)) = ∂ ∂ exp γ(t) (1/T (-β 0 (t) + α i (0))).
The transported initial momentum vector α i (0) is noted θ i (0). After parallel transporting all residuals, all initial momentum vectors are defined in B 0 .

Feature extraction for statistical analysis

Each transported initial momentum vectors θ i (0) is of size 3 × N B0 , where N B0 is the number of control point used to parametrise the geodesics.

Jacobian determinants are a geometric measure derived from the full deformation tensor that is commonly used to study shrinkage or growth of the surface. In this work

we propose an analysis framework where we decouple the amplitude and the orientation of the deformation. Such an approach will still analyse growth and shrinkage, but also other geometric aspects, such as rotation and torsion, that are not captured by the surface Jacobian. Furthermore, the changes being analysed are residual deformations, which are defined using a purely geometrical spatio-temporal regression. As such, the shape differences that we aim to detect are not necessarily limited to shrinkage or growth, but can be induced by more complex effects.

To analyse direct measures from deformation and to avoid losing statistical power from doing a large number of comparisons, we propose an original clustering approach by grouping the parametrisation (B 0 ; β 0 (0)) of the spatio-temporal regression φ into clusters.

To do so, we defined a similarity measure derived from the positions of the control points x p , the pairwise angles and the magnitudes of the initial momentum vectors {β 0 p (0)} 1≤p≤N B 0 attached to the control point x p . The difference between two control points x p and x q ∀p, q ∈ {1; ...; N B0 } is defined by the euclidean distance, the angle between two vectors is defined by the cosine. The similarity between p and q is defined by

s(p, q) = -5 x p -x q 2 + 2(cos(β 0 p , β 0 q ) + 1) -| β 0 p 2 -β 0 q 2 |.
Parameters are chosen to balance between vector similarity and control point positions and depend on the distance in mm between two points. The distance is determined by the kernel K V so that clusters encompass control points and their momentum vectors within the same area and look alike. To estimate those clusters, we used a spectral clustering method [START_REF] Luxburg | A Tutorial on Spectral Clustering[END_REF] using the discretisation approach presented in [START_REF] Jianbo | Multiclass Spectral Clustering[END_REF] for initialisation, as it has been shown to be more stable than other approaches such as k-means for initialisation. 3000 different initialisations are generated and we select the best one in term of inertia for spectral clustering. We chose 10 clusters as thought this would be a good balance between reducing the number of multiple comparisons while maintaining some spatial specificity in the analyses and equitable clusters. A mean vector is then computed from the parallel transported residuals defined on the control points of the cluster. This is done for each cluster and for each subject. We then obtain N vectors {ν i,k } per cluster k, and 10 vectors per subject i.

For the statistical analysis, we will use two uncorrelated descriptors for the vectors {ν i,k }: the amplitude and the orientation. The orientation of the vectors {ν i,k } is originally represented by 3 angles, one per axis. The angles are then projected via a Principal Component Analysis onto the first eigenvector, therefore the orientation of {ν i,k } considered here is represented by one continuous scalar, leading to the set of responsive variables {O i,k }.

Data and application

As previously mentioned, we applied the proposed framework to the GENFI study and used the thalamus as structure of interest.

Dataset description

All participants included in this study come from the data freeze 1 of the GENFI cohort described in detail in [START_REF] Rohrer | Presymptomatic cognitive and neuroanatomical changes in genetic frontotemporal dementia in the Genetic Frontotemporal dementia Initiative(GENFI) study: a cross-sectional analysis[END_REF]. Initial results from this cohort [START_REF] Rohrer | Presymptomatic cognitive and neuroanatomical changes in genetic frontotemporal dementia in the Genetic Frontotemporal dementia Initiative(GENFI) study: a cross-sectional analysis[END_REF] show volumetric differences in the thalamus at least 5 years before expected age of onset with an effect in all genetic subtypes, and so we chose this well-defined subcortical structure for further analysis. In this paper we used 211 participants, 113 mutation carriers (MAPT=26, GRN=53, C9ORF=34) and 98 non-carriers. All participants have a T1-weighted (T1w) MRI available and an associated expected years to symptom onset (EYO). The EYO, ranging from -40 years to +20 years, is calculated as the difference between the age of the participant at the time of the T1w acquisition and the mean age at onset of affected family members, as in [START_REF] Rohrer | Presymptomatic cognitive and neuroanatomical changes in genetic frontotemporal dementia in the Genetic Frontotemporal dementia Initiative(GENFI) study: a cross-sectional analysis[END_REF]. The median of the age at onset of all subjects is 59.7 years with inter-quartile range IQR = 60.5 -55. Table 1 shows the demographics of the study participants used in this analysis.

Application to the thalamus

As previously mentioned in section 2.2, all T1w MRIs and associated segmentations of the structure of interest, the thalamus, are first aligned to a common space. This enables to normalise for intra-cranial volume differences across subjects. We then extracted the meshes corresponding to the thalamus, composed by around 2, 300 vertices. This resulted in 211 thalamus meshes, representing the mean left and the right shape. Each were associated with the EYO of the corresponding subject as well as mutation status:

non-carrier and mutation carrier (MC). For the spatio-temporal regression, we used 30 time points, which corresponds approximately to one time point every two years. The space of deformations V was defined using a 11mm width kernel, approximately half of the length of the thalamus, which leads to a set of 288 control points. For the space of varifolds we used a 5mm width kernel which covers the size of 2 voxels. This parameter was fixed and thought to be a good compromise between the capture of high frequency changes and the robustness of the approach to noise in the shape segmentation.

Similarly to the volumetric analysis performed by Rohrer et al. [START_REF] Rohrer | Presymptomatic cognitive and neuroanatomical changes in genetic frontotemporal dementia in the Genetic Frontotemporal dementia Initiative(GENFI) study: a cross-sectional analysis[END_REF], we used a mixed effect model to study the shape difference between the non-carriers and mutation carriers.

Amplitude {|ν i,k |} and orientation {O i,k } were used as responsive variables and the fixed effects predictors of interest were mutation carrier status, EYO, interaction between mutation carrier status and EYO, sex and the site in which the subject has been scanned.

A random intercept for family allows values of the marker to be correlated between family members.Correcting for age of subjects is not relevant here, since there is a strong correlation (r = 0.9) between EYO and age. We performed a Wald test for every model, assessing the difference between the mutation carrier group and the non-carrier group, and the evolution of differences across time. For each analysis with statistically significant differences between both groups, further Wald tests were conducted every 5 years as in the volumetric analysis [START_REF] Rohrer | Presymptomatic cognitive and neuroanatomical changes in genetic frontotemporal dementia in the Genetic Frontotemporal dementia Initiative(GENFI) study: a cross-sectional analysis[END_REF] to assess how long before the expected onset we could detect changes between mutation carriers and controls.

Results

Results for the amplitude and the orientation of the residual momentum vectors are presented in Table 2. We found significant differences, after correction for multiple comparisons, in cluster 1 and cluster 4, for both tests; T1:differences between MC and controls and T2: differences over time between MC and controls. Those differences are significant after Bonferroni correction for multiple comparisons (20 tests). Cluster 1 shows differences in the orientation, and no differences in the amplitude, whereas shows significant differences between the mutation carriers and controls, 5 years before EYO (p = 0.03), the uncorrected for this cluster is p = 2e-3, to keep a head to head comparison with the previous studies on this dataset [START_REF] Rohrer | Presymptomatic cognitive and neuroanatomical changes in genetic frontotemporal dementia in the Genetic Frontotemporal dementia Initiative(GENFI) study: a cross-sectional analysis[END_REF][START_REF] Cury | Spatio-Temporal Shape Analysis of Cross-Sectional Data for Detection of Early Changes in Neurodegenerative Disease[END_REF] in which the p-values at -5 EYO was significant but higher than here. The uncorrected p-values show significant differences at 10 years before EYO, with p=0.048 for the orientation of cluster 1. The amplitude between the two groups doesn't differ significantly for the cluster 4 before EYO for corrected p-values, and differs 5 years before onset without correction (p=0.05). Since the number of clusters used (i.e. 10), is an arbitrary choice, we tried to reproduce the results with different number of clusters. We performed the analysis for 2, 4, 6, 8, 10, 12, 14 and 16 clusters which results can be found in supplementary material (https://zenodo.org/record/1324234). For 6 clusters and 16 clusters, there were differences in orientation for one of the clusters which deformation corresponds to the one of cluster 1 (see Figure 3). From 8 clusters to 14 clusters, we found a cluster with strong differences 5 years before the expected onset (p < 0.01) in orientation whose deformation corresponds again to the one of the cluster 1 (p = 0.003). The change in orientation for the deformation recovered within cluster 1 (see Figure 3) appears to be stable for different clustering of the parametrisation of the global spatiotemporal trajectory (https://zenodo.org/record/1324234, Figure 1).

Discussion and conclusion

We applied a novel method of statistical shape analysis to a cohort of individuals with genetic FTD in order to localise any presymptomatic differences present in the shape of the thalamus. From the analysis, we conclude that differences are observed five years before expected symptom onset. While volumetric analysis [START_REF] Rohrer | Presymptomatic cognitive and neuroanatomical changes in genetic frontotemporal dementia in the Genetic Frontotemporal dementia Initiative(GENFI) study: a cross-sectional analysis[END_REF] and our initial shape analysis [START_REF] Cury | Spatio-Temporal Shape Analysis of Cross-Sectional Data for Detection of Early Changes in Neurodegenerative Disease[END_REF] also found these changes, this method showed significance that survived correction for multiple comparisons. The change in shape is primarily attributable to differences in orientation of the deformation rather than changes in amplitude of the deformation, which would imply a simple scaling effect of the region. This result confirms our previous shape analysis in this cohort [START_REF] Cury | Spatio-Temporal Shape Analysis of Cross-Sectional Data for Detection of Early Changes in Neurodegenerative Disease[END_REF] that was performed at a global level through a kernel principal component analysis. The first mode of variation which detected significant shape differences around the same point with respect to EYO did not capture volume differences but only changes in the orientation of the deformation. The results of those studies seem to indicate that shape changes occur before volume changes.

The regions of the thalamus most affected in the analysis are anterior, overlapping with the anterior nuclei group. The main connections of these nuclei are to the prefrontal cortices, an area universally affected in all genetic forms of FTD. To illustrate this purpose, we used the Oxford thalamic connectivity atlas, a thalamic atlas based on its anatomical connectivity to the cerebral cortex [START_REF] Behrens | Non-invasive mapping of connections between human thalamus and cortex using diffusion imaging[END_REF], and displayed at Figure 4 the atlas next to the clusters 1 and 4. Whilst differences are seen in cortical involvement within the different genetic forms of FTD [START_REF] Cash | Patterns of gray matter atrophy in genetic frontotemporal dementia: results from the GENFI study[END_REF], it may well be that this joint analysis of GRN, C9orf72 and MAPT mutations is only identifying thalamic regions jointly affected.

This approach could also be used to explore other regions known to be implicated in FTD, such as the insular cortex, which is located in the lateral sulci and is connected to the limbic system, and to the thalamus. In fact, it would be interesting to analyse the insula and thalamus together, and the insula only, so we could investigate if shape changes in both structures are linked. The small numbers in each group precluded any analysis of the individual genetic types, but it will be important to investigate future data freezes from the GENFI study with larger numbers, particularly the C9orf72 group who have been shown to have early thalamic involvement [START_REF] Cash | Patterns of gray matter atrophy in genetic frontotemporal dementia: results from the GENFI study[END_REF].

Future studies should also evaluate the initial momentum vectors of individual geodesic evolution of shapes from each subject, through longitudinal data. Those individual evolutions would provide information on the differences of evolutions of shape between the mutation carriers and the controls.

(PR/ylr/18575) and Alzheimers Society UK (AS-PG-15-025). We would like to thank the participants and their families for taking part in the GENFI study.

Figure 1 :

 1 Figure1: Overview of the proposed regression approach. The temporal axis indicates the time variable attached to the data, in this case the estimated years to expected symptom onset. The residual deformations (step 2) ρ i parametrised by (φ(B 0 , t i ); α i (0)) computed from the common trajectory (step 1) φ parametrised by (B 0 ; β 0 ), cannot be analysed because they are defined on different spaces i.e. φ(B 0 , t i ). They have to be transported to a common space (i.e. B 0 ) along the geodesic φ, so they can be analysed (step 3).

cluster 4 Figure 2 ,

 42 Figure 2, the p-values and confidence intervals are corrected for multiple comparison across time using Bonferroni correction. The orientation of the cluster 1 deformation

Figure 2 :Figure 3

 23 Figure 2: cluster 1 (orientation component) and cluster 4 (amplitude component) estimates in mutation carriers and controls, by estimated time from expected symptoms onset (EYO). p-values and confident interval are Bonferroni corrected. * : p < 0.05, ** : p < 0.01, *** : p < 0.001

Figure 3 :

 3 Figure 3: Deformation obtained by the momentum vectors (displayed here and coloured by amplitude) of Cluster 1 and Cluster 4. The colour map is in millimetres and indicates the displacement due to the corresponding deformation (blue meshes). The scale for Cluster 1 range from 0 mm to 9 mm, and from 0 mm to 2.8 mm for Cluster 4.

Figure 4 :

 4 Figure 4: Thalamic connectivity atlas, and deformations clusters 1 and 4. The orientation of cluster 1 leads to significant differences between MC and controls 5 years before EYO.

Table 1 :

 1 Data demographics, in absolute values.

		Non-carriers	Mutation carriers
		n=98	n=113
	Males	59	56
	Asymptomatic	98	76
	Age in years (med (IQR))	50.2 (62.1 -36.6=25.5)	52.7 (62.7 -41.1=21.6)
	Years from expected onset:
	≤ -20 years	30	21
	-20 ≤ years ≤ -10	16	21
	-10 ≤ years < 0	23	22
	0 ≤ years	29	49

Table 2 :

 2 p-values with the corresponding χ 2 value, resulting from the Wald tests testing the mutation carrier (MC) differences (test T1), and the evolution of those differences along time (test T2), for the amplitude of the initial momentum vector and its orientation, for the clusters showing at least one significant test. Bold p-values: ≤ 0.05, and starred (*) p-values indicate the corrected threshold for multiple comparisons: ≤ 2.5e-3.

			C1	C2	C4	C6	C7
	Ampl.	T1 T2	p = 0.48 χ 2 df =2 = 1.43 p = 0.24 χ 2 df =1 = 1.37	p = 0.51 χ 2 df =2 = 1.35 χ 2 p =1.5e-3 (*) p = 0.08 df =2 = 12.94 χ 2 df =2 = 5.10 χ 2 p = 0.76 df =2 = 0.55 p = 0.26 p =1.5e-3 (*) p =0.04 p = 0.68 χ 2 df =1 = 1.28 χ 2 df =1 = 10.08 χ 2 df =1 = 4.20 χ 2 df =1 = 0.17
	Orient.	T1 T2	p =2e-4 (*) χ 2 df =2 = 16.60 χ 2 p = 0.12 df =2 = 4.17 χ 2 p = 0.85 df =2 = 0.33 p =9e-4 (*) p =0.05 p = 0.62 χ 2 df =1 = 11.01 χ 2 df =1 = 3.85 χ 2 df =1 = 0.25	p = 0.63 χ 2 df =2 = 0.92 χ 2 p = 0.08 df =2 = 5.06 p = 0.34 p =0.04 χ 2 df =1 = 0.91 χ 2 df =1 = 4.29

http://www.deformetrica.org/

Acknowledgements

Claire Cury is supported by the EU-FP7 project VPH-DARE@IT (FP7-ICT-2011-9-601055). Stanley Durrleman has received funding from the program Investissements davenir ANR-10-IAIHU-06 and the European Unions Horizon 2020 research and innovation programme EuroPOND under grant agreement No 666992, and is funded by the European Research Council (ERC) under grant agreement No 678304. Marco Lorenzi received funding from the EPSRC (EP/J020990/1). Jennifer Nicholas is supported by UK Medical Research Council (grant MR/M023664/1). David Cash is supported by grants from the Alzheimer Society(AS-PG-15-025), Alzheimers Research UK (ARUK-PG2014-1946) and Medical Research Council UK (MR/M023664/1). JBR is supported by the Wellcome Trust (103838). Jonathan D. Rohrer is an MRC Clinician Scientist and has received funding from the NIHR Rare Diseases Translational Research Collaboration. Sebastien Ourselin receives funding from the EPSRC (EP/H046410/1, EP/K005278), the MRC (MR/J01107X/1), the NIHR Biomedical Research Unit (Dementia) at UCL and the National Institute for Health Research University College London Hospitals Biomedical Research Centre (NIHR BRC UCLH/UCL High Impact Initiative-BW.mn.BRC10269). Marc Modat is supported by the UCL Leonard Wolfson Experimental Neurology Centre

List of other GENFI consortium members

Christin Andersson -Department of Clinical Neuroscience, Karolinska Institutet, Stockholm, Sweden


RAPPORT INTERNE LPSC 07-47

Développement de codes de simulations de FFAG spiral et de programmes d'automatisation de balayage de paramètres dans le cadre du projet RACCAM.

Joris Fourrier pour le projet RACCAM

Table des matières.

INTRODUCTION.

I. Développement de procédures de simulations pour FFAG spiral.

1) Modélisation du champ magnétique.

2) Tombées de champ.

3) Calcul des dérivées du champ magnétique. 4) Conclusion.

II. Automatisation des recherches de paramètres dynamiques.

Introduction.

RACCAM [START_REF] Méot | The FFAG R&D and Medical Application Project RACCAM[END_REF]2] (Recherche en Accélérateurs et Applications Médicales) est un projet financé par l'agence nationale de la recherche pour trois ans regroupant trois organismes : un laboratoire de recherche publique : le LPSC (Laboratoire de Physique Subatomique et de Cosmologie), entité de l'IN2P3 / CNRS, un industriel : ΣΦ (Conception et Réalisation d'aimants pour accélérateurs de particules) et le centre hospitalier universitaire de Grenoble. Son but est de développer l'activité de recherche et développement dans le domaine des accélérateurs à champ fixe (FFAG) en France, de participer aux projets internationaux de recherche sur ces mêmes FFAG, tels l'usine à neutrinos [START_REF]A feasibility study of a neutrino factory in Japan[END_REF] ou la construction d'EMMA [START_REF]Introduction to the Non-scaling Electron Model FFAG EMMA, Rob Edgecock and the EMMA Collaboration[END_REF] (modèle en électrons de FFAG) et d'étudier leur utilisation dans le domaine de la radiothérapie des cancers par faisceau de protons [START_REF] Misu | Study of Compact Medical FFAG Accelerators[END_REF].

Dans ce cadre, la mise au point et la réalisation d'un prototype d'aimant FFAG est en cours. La technologie particulière d'un aimant à face spiralée a été choisie en raison de ses avantages en terme de dynamique faisceau et de compacité de la machine permettant ainsi d'en limiter le coût global.

La phase d'étude théorique précédant la construction de l'aimant est en cours d'achèvement, elle a essentiellement consistée dans la mise au point d'outils de simulation de dynamique de particules dans de tels aimants et à l'automatisation des simulations afin de pouvoir étudier de nombreuses configurations de paramètres et leurs conséquences sur la qualité de la dynamique.

Pour ce faire, nous avons eu recours au code de calcul matriciel BeamOptics [START_REF] Beamoptics | A Program for Analytical Beam Optics[END_REF] utilisant le logiciel de calcul formel Mathematica pour les études de paramètres du premier ordre, à savoir les tunes, les fonctions betatroniques et de dispersion, ainsi qu'au logiciel de simulation Zgoubi [START_REF]Zgoubi user's guide[END_REF] afin d'effectuer les études de dynamique complète qui constituent la majeure partie du travail effectué.

Ce logiciel repose sur une technique d'intégration numérique des trajectoires des particules dans différents éléments d'optique modélisés tels que des dipôles, des quadripôles, des aimants de FFAG à secteur radial [START_REF] Zgoubi | Méot for the RACCAM Collaboration[END_REF][START_REF] Lemuet | Developments in the ray-tracing code Zgoubi for multiturn tracking in FFAG rings[END_REF]. Il permet également la simulation de trajectoires dans des cartes de champs magnétiques créées par des codes de simulation 2D ou 3D comme TOSCA, qui donnent des distributions magnétiques issues de calculs éléments finis.

I. Développement de procédures de simulations pour

FFAG spiral.

1) Modélisation du champ magnétique.

Afin de simuler les trajectoires de particules dans un aimant FFAG spiral, une procédure permettant de calculer la composante verticale du champ magnétique ) , ( θ r B Z et ses dérivées dans le plan médian a été mise au point. Cette méthode a été basée sur une procédure existante concernant les FFAG à secteur radial [START_REF] Lemuet | Developments in the ray-tracing code Zgoubi for multiturn tracking in FFAG rings[END_REF].

En se basant sur les modèles de champ magnétique d'aimant FFAG spiral de l'institut MURA [START_REF]Spiral sector magnets[END_REF] : • f défini le flutter magnétique :

⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = θ θ N w R r f R r B r B k Z Z 0 0 0 ln sin 1 ) , ( (1
( )

2 2 Z Z Z B B B f - = (1.2)
• w est une variable dépendant de l'angle spiral. w = ) tan( 1 ξ N ξ étant l'angle spiral défini comme l'angle entre la tangente à la face magnétique et un rayon passant par le centre de la machine (voir Fig. 1).

•

Le facteur

⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + θ N w R r f 0 ln
) ( ) ( ) , ( 0 0 s s e e k Z Z d F d F R r B r B ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = θ (1.3) où : • 0 R , 0 Z B , et k R r ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ 0
ont les mêmes définitions qu'à la page précédente.

•

) ( ), ( ξ étant l'angle spiral défini comme l'angle entre la tangente à la face magnétique et un rayon passant par le centre de la machine (voir Fig. 1).

La forme particulière de ces faces est prise en compte dans le calcul des facteurs de 2) Tombées de champs.

Les facteurs de tombée de champ à l'une des faces magnétiques s'écrivent :

5 5 2 2 1 0 ... ) ( , )] ( exp[ 1 1 ) ( ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + = + = g d C g d C g d C C d P d P d F i i i i i i i (1.

C C

sont les coefficients de champ fuite (coefficients de Enge [START_REF] Deflecting Magnets | Focusing of charged particles[END_REF]).

La distance i d est calculée numériquement en utilisant une méthode permettant de résoudre l'équation en θ [START_REF] Martinache | Optics and field map for a FFAG spiral[END_REF]: Le gap est rendu dépendant de r afin de pouvoir étudier différentes configurations (voir Fig. 6) :

0 )) sin( ) cos( ))( sin( ( )) sin( ) cos( ))( cos( ( = + + + + - + + -θ θ θ R be R e R be Y R e R be R
• gap shaping : k r R g r g ≈ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = κ κ , ) ( 0 0 (1.7) assurant en général la loi k Z r B ≈ • gap à faces parallèles: 0 , ) ( 0 = = κ g r g (1.8) • gap proportionnel à r : 1 , ) ( 0 0 - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = κ R r g r g (1.9)
assurant que les nombres d'onde verticaux soient indépendants de l'énergie.

La focalisation de coin peut être décrite par le formalisme matriciel suivant [START_REF] Brown | A First and Second Order Matrix Theory for the Design of Beam Transport Systems and Charged Particle Spectrometers[END_REF] :

ρ ψ ε ) tan( 43 - = R ; ξ π ε ± - = N pf ) 1 (
(1.10) (-à la face d'entrée , + à la face de sortie) avec ε l'angle de coin (angle entre la face magnétique et la trajectoire de la particule), ρ le rayon de courbure de la trajectoire dans l'aimant et ψ la correction apportée pour les effets de champ de fuite. Cette correction s'écrit [START_REF] Brown | A First and Second Order Matrix Theory for the Design of Beam Transport Systems and Charged Particle Spectrometers[END_REF]:

( )

) cos( ) ( sin 1 2 1 ε ε ρ ψ + = I g (1.11)
avec g le gap magnétique, ρ le rayon de courbure de la trajectoire dans l'aimant , ε l'angle de coin et I 1 l'intégrale suivante :

[ ] ds gB s B B s B I Z Z ∫ ∞ + ∞ - - = 2 0 0 1 ) ( ) (
(1.12) s étant la coordonnée curviligne le long de la trajectoire.

Fig. 6 : Formes possibles du gap : gap shaping (1), gap à faces parallèles(2), gap proportionnel à r(3).

Comme

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = N pf N R π π ρ 2 sin . 2 sin
, le fait de rendre g proportionnel à R rend la correction de champ de fuite indépendante de R et donc de l'énergie. Cela impose donc les nombres d'onde verticaux constants avec l'énergie.

Cette dépendance du gap avec le rayon permet de jouer sur la dépendance du nombre d'onde vertical avec l'énergie. • Fabrication d'une carte de champ 2D dans le plan médian.

) 1 ( ) 2 ( ) 3 ( ) (m R ) (m R ) (m R ) (cm Z ) (cm Z ) (cm Z
Une procédure Fortran dénommée GeneMap permet de calculer le champ magnétique en tout point dans l'aimant et à l'extérieur de l'aimant. Cette carte de champ peut ensuite être directement utilisée dans le code de tracé de trajectoire Zgoubi [START_REF]Zgoubi user's guide[END_REF] afin d'effectuer le tracking pas à pas de particules dans ce modèle d'aimant FFAG spiral.

Cette procédure GeneMap donne les valeurs de champ magnétique dans tout l'aimant grâce à la méthode précédente. Dans ce cas, l'interpolation ne se fait pas avec un maillage « flottant » mais à partir du maillage de la carte.

4) Conclusion.

La 1 ère méthode (FFAG-SPI) autorise la modification facile et automatique des paramètres (voir II.4) et donne une excellente précision.

La 2 ème méthode (carte de champ GeneMap) autorise également la modification rapide et automatique des paramètres mais demande le calcul complet de la carte de champ entre les rayons d'injection et d'extraction même si l'étude se restreint à une énergie précise. Cela implique un supplément de temps de calcul non négligeable si de nombreuses simulations sont nécessaires. Elle offre également une bonne précision numérique.

Une procédure permettant de simuler les trajectoires de particules dans un FFAG spiral a été efficacement mise au point. Elle nous autorise à effectuer du tracking multitour dans de telles structures magnétiques afin d'en étudier les comportements et paramètres dynamiques tels que les orbites fermées, les tunes, les limites de stabilité.

Ces études de dynamique étant itératives, elles demandent un grand nombre de simulations Zgoubi. C'est pour cela qu'une méthode d'automatisation des recherches de paramètres dynamiques a été développée.

II. Automatisation des recherches de paramètres

dynamiques.

1) Recherche d'orbites fermées.

A une énergie donnée, la connaissance de l'orbite fermée est essentielle. En effet, à partir de la connaissance de cette orbite fermée les paramètres dynamiques peuvent être déterminés comme les nombres d'onde, les fonctions beta, de dispersion.

Elle est également à la base des études de dynamique complète comme les recherches d'acceptances horizontale et verticale, les études d'accélération de particules ou de paquets de particules. Nous allons décrire le cheminement mis en oeuvre afin de déterminer ces orbites fermées. Le programme de simulation Zgoubi requiert la connaissance des coordonnées à l'entrée de la cellule composée d'une demi section droite, de l'aimant FFAG puis d'une seconde demi section droite. L'opération est la suivante, la trajectoire de la particule aux coordonnées théoriques est suivie pendant une trentaine de passages dans la cellule FFAG (correspondant à environ 4 tours dans l'anneau représenté en Fig. 1). Ce nombre de passages peut être bien évidemment modifié en fonction de la précision voulue et du temps de calcul rencontré. La trajectoire dans l'espace des phases horizontal (X, X') est généralement une ellipse en mouvement paraxial si les nombres d'onde sont suffisamment éloignés de

... 3 , 2 , 1 , 1 , = ≈ n n Q Q Z X
Le centre de l'ellipse est l'orbite fermée recherchée.

Le programme de recherche d'orbite fermée va donc faire les opérations suivantes :

• Tracking multitour de la particule aux coordonnées théoriques.

• Relevé des coordonnées du centre de l'ellipse.

• Test d'écart entre la particule lancée et le centre de l'ellipse à la précision voulue près. • Ecriture des coordonnées du centre de l'ellipse dans le fichier de simulation. • Tracking multitour de la particule aux nouvelles coordonnées.

Le résultat de cette recherche est le couple de coordonnées

2) Calcul des nombres d'onde.

A partir de l'orbite fermée, il est nécessaire de connaître la valeur des nombres d'onde horizontaux et verticaux Cette méthode consiste à calculer les coefficients de la matrice de transfert de la cellule optique traversée par la particule de référence, c'est-à-dire sur l'orbite fermée, par une méthode d'interpolation entrée / sortie. Elle permet de calculer les termes de 1 er ordre ij R ou de 2 nd ordre ijk T , au-delà la précision de l'interpolation peut être discutable.

Pour le calcul des tunes, on ne s'intéresse qu'aux termes de 1 er ordre, celui-ci s'effectue par la génération de 11 trajectoires autour de l'orbite fermée grâce à une option de Zgoubi.

La simulation de ces 11 trajectoires se fait sur une cellule et les nombres d'onde sont relevés dans le fichier résultat créé par Zgoubi. Les avantages de cette méthode sont la précision des nombres d'onde pour un calcul rapide. L'inconvénient est que cette méthode n'est valable que pour l'optique paraxiale et ne peut être utilisée que pour les calculs de nombres d'onde au voisinage des orbites fermées. Une étude sur des particules éloignées de l'orbite fermée demandera donc une autre méthode basée sur un développement en série de Fourier.

b) Développement en série de Fourier.

Le programme d'analyse de résultat Zpop, associé au logiciel Zgoubi, est capable de calculer les tunes par analyse de Fourier après un tracking multitour (une centaine de passages dans la cellule FFAG sont en général suffisants pour notre étude) 1 . La nécessité d'utiliser Zpop après chaque simulation est longue et fastidieuse, la routine d'analyse de Fourier a donc été reprise pour créer un programme indépendant de Zpop pour l'analyse des tunes.

Dans le fichier résultat généré par Zgoubi, ce programme va lire les coordonnées de la particule tout au long du tracking, effectue un matching d'ellipse dans l'espace des phases ) , ( ' X X , puis le calcul des nombres d'onde associés au mouvement. Les avantages sont la précision et la possibilité de calculer des tunes pour une optique non paraxiale (particules très éloignées de l'orbite fermée, pour une étude de detuning d'amplitude par exemple).

Malheureusement, cette méthode demande un tracking sur une centaine de cellules d'où un temps de calcul augmenté.

De ce fait, la méthode matricielle est utilisée pour le calcul des tunes au voisinage des orbites fermées, tandis que la méthode de développement en série de Fourier sert dans le cas de particules de grande excursion. 3) Recherche des limites de stabilité.

A partir de la connaissance de l'orbite fermée, il est nécessaire de poursuivre les études de dynamique pour déterminer les acceptances dynamiques transverses de la machine, c'est-à-dire les dimensions maximales des faisceaux qui peuvent circuler dans l'anneau. Pour cela, on cherche à déterminer les limites de stabilité horizontale et verticale qui sont données par les trajectoires stables les plus éloignées de l'orbite fermée.

a) Limite de stabilité horizontale.

A partir de l'orbite fermée, on impose un mouvement vertical infinitésimal (afin de laisser apparaître d'éventuels phénomènes de couplage entre les plans horizontaux et verticaux) puis le programme augmente la valeur de X 0 et enfin lance un tracking sur plusieurs centaines de cellules (environ 300 dans notre étude). Si le mouvement est stable, il augmente à nouveau X 0 jusqu'à ce que la dernière position stable soit déterminée. Celle-ci est alors définie comme la limite de stabilité horizontale.

Précédemment, la procédure de recherche de nombres d'onde par développement en série de Fourier a été évoquée. Lors de ce calcul, la trajectoire de la particule dans l'espace des phases ) , ( ' X X est interpolée par une ellipse par une procédure de matching. Cette procédure permet le calcul des paramètres α et β de l'ellipse ainsi que son émittance (surface) ε selon l'équation :

ε β α γ = + + 2 2 ' X ' X X 2 X (2.1)
La même opération est effectuée pour la recherche de limite de stabilité, nous pouvons ainsi connaître, non seulement les nombres d'onde associés à la particule sur la limite de stabilité mais également déterminer l'émittance de celle-ci si la trajectoire est une ellipse. Dans le cas où cette trajectoire n'est pas elliptique comme le montre la Fig. 14, le matching nous donne une indication sur l'émittance, cependant une procédure plus précise devra être développée. Cette valeur nous indique l'acceptance transverse horizontale de la machine.

Grâce à cette détermination automatique de la limite de stabilité horizontale, le temps de calcul est grandement amélioré par rapport à une même étude faite en modifiant à la main les paramètres, la détermination de l'acceptance horizontale est rapide, efficace.

Par ailleurs, on peut lancer plusieurs particules dont les coordonnées sont comprises entre l'orbite fermée et la limite de stabilité afin de faire des études de déformation d'ellipse ou de detuning d'amplitude. Enfin, les éventuels phénomènes de couplage ainsi que la proximité de résonances peuvent être observés grâce à ces études automatiques. De la même manière que précédemment, à partir de l'orbite fermée, le programme augmente la coordonnée verticale Z 0 de la particule lancée jusqu'à ce que le mouvement ne soit plus stable. La dernière trajectoire stable définit la limite de stabilité à laquelle on associe l'acceptance verticale calculée par le matching d'ellipse décrit auparavant.

La prochaine étape sera de marier les deux procédures en simultané. Cela reviendra à chercher la limite de stabilité verticale pour différentes valeurs de X comprises entre l'orbite fermée et la limite de stabilité horizontale. On pourra ainsi voir l'empreinte maximale du faisceau dans le plan transverse ) , ( Z X

. A ce stade, on procédera peut être en lançant un paquet 4D sur-dimensionné et en regardant quelles trajectoires restent stables après plusieurs centaines de tours. On voit que l'accord est bon entre les 2 codes pour des valeurs de k et ξ peu élevées, les régions de stabilité se recoupent et les tunes sont très proches.

Les différences entre les tunes augmentent avec k et ξ , cela reste à être compris mais les hypothèses sont :

• La mimite de validité de l'hypothèse que le rayon de l'orbite est un cercle dans l'approche matricelle. • Les effets des champs de fuite plus forts à fort angle spiral : la particule voyage plus longtemps dans les régions de champ de fuite. L'ouverture dynamique s'effondre à proximité de la résonance

3 2 2 Q X + = .
Grâce à ces données, on peut tracer le diagramme de tunes en incluant les ouvertures dynamiques horizontale et verticale par un dégradé de nuances de gris. Comme on peut le voir ci-dessous, les points les plus clairs sont ceux dont l'ouverture dynamique est la plus faible, ceux les plus foncés ont les ouvertures dynamiques les plus grandes. On voit très bien que la résonance On peut également déterminer l'évolution de l'ouverture dynamique avec l'énergie ainsi que l'acceptance à l'injection, là où le faisceau est le plus gros.

1 ) 3 )

 13 Recherche d'orbites fermées. a) Détermination d'une trajectoire théorique. b) Tracking multitour. 2) Calcul des nombres d'onde. a) Méthode matricielle. b) Transformée de Fourier. Recherche des limites de stabilité. a) Limite de stabilité horizontale. b) Limite de stabilité verticale.

Fig. 1 :

 1 Fig. 1 : Schéma de principe d'un anneau FFAG à secteur spiral à 8 cellules.

R

  est le rayon de référence, pris comme l'excursion radiale de la particule au centre de l'aimant pour l'énergie d'extraction. • 0 Z B est le champ magnétique pris en 0 R au centre de l'aimant, il correspond au champ magnétique maximal vu par la particule ayant l'énergie d'extraction. radiale du champ magnétique, dépendance en puissance de k caractéristique d'un aimant FFAG scaling. • N est le nombre d'aimants dans l'anneau.

  de champ magnétique aux faces d'entrée et de sortie de l'aimant.Les faces de l'aimant ont une forme spirale répondant à l

Fig. 4 :Fig. 5 :

 45 Fig. 4 : Loi de tombée de champ aux faces magnétiques du FFAG spiral.

ω

  étant le positionnement de la face magnétique dans le référentiel choisi.Cela revient à dire que la normale à la face magnétique au point )

Fig. 7 :

 7 Fig. 7 : Champ magnétique dans un aimant FFAG spiral.

Fig. 8 :

 8 Fig. 8 : Interpolation numérique pour la procédure FFAG-SPI.

Fig. 9 :

 9 Fig. 9 : Interpolation numérique pour l'utilisation de cartes de champ.

Fig. 10 :

 10 Fig. 10 : Schéma de l'aimant RACCAM et paramètres pour le calcul des orbites théoriques.

  a) Détermination d'une trajectoire théorique. Une trajectoire théorique est tout d'abord calculée de manière géométrique grâce aux paramètres suivant : le rayon d'extraction des particules (pris comme référence pour la machine) 0 R , le nombre d'aimants N et le packing factor pf .

Fig. 11 :

 11 Fig. 11 : Fichier data de Zgoubi pour la procédure FFAG-SPI.

Fig. 12 :

 12 Fig. 12 : Schéma de principe de la recherche d'orbite fermée.

  qui est ensuite utilisé pour le calcul des nombres d'onde et les recherches de limites de stabilité.

Fig. 13 :

 13 Fig. 13 : Fichier data de Zgoubi pour le formalisme matriciel.

  points de fonctionnement de la future machine. La connaissance des tunes est également intéressante pour les particules éloignées de l'orbite fermée car elle fournit le variation des nombres d'onde avec l'amplitude (detuning d'amplitude). Pour cela, le logiciel Zgoubi propose deux méthodes : le formalisme matriciel et l'utilisation d'un développement en série de Fourier. a) Méthode matricielle.

Fig. 14 :

 14 Fig.14: Schéma de principe de la recherche de limite stabilité horizontale.

Fig. 15 :

 15 Fig. 15 : Schéma de principe de la recherche de limite stabilité verticale.

  Le but de RACCAM est la mise au point d'un prototype d'aimant FFAG à secteur spiral pour une machine de traitement de cancer par bombardement de protons. Avant de lancer la construction d'un tel aimant, les paramètres magnétiques, géométriques et dynamiques doivent être fixés. En particulier il s'agit de déterminer le point de fonctionnement de la future machine.On désire se placer dans une région du diagramme de nombres d'onde) de dynamique est la meilleure possible, c'est-à-dire des régions éloignées de toute résonances nuisibles et pour lesquelles l'ouverture dynamique est grande, régions appropriées à l'injection et à l'extraction du faisceau.De nombreux paramètres géométriques et magnétiques ont une influence sur l'évolution des nombres d'onde) rayon maximal de la machine, le nombre de cellules N, le packing factor pf, l'indice de champ k ou l'angle spiral ξ .Pour cela, on se propose de fixer le rayon maximal, le nombre de cellules et le packing factor pour des considérations techniques : injection, extraction, cavité radiofréquence, diagnostic, encombrement de la machine. A partir de là, on fait varier indépendamment ou simultanément les paramètres magnétique k (indice de champ magnétique) et géométrique ξ (angle spiral donnant la forme des faces magnétiques) dans les procédures de simulation FFAG-SPI ou GeneMap.Cette procédure automatique de balayage de paramètres en ) , ( ξ k permet d'effectuer les opérations automatiques précédentes (recherche d'orbite fermée, de tunes, de limites de stabilité) pour différents couples ) , ( ξ k et donc pour différents points de fonctionnement. A partir de cette étude seront choisis quelques points de fonctionnement possibles pour des investigations plus approfondies, et pour point de départ des études de conception de l'aimant.

Fig. 16 :

 16 Fig. 16 : Exemple de région de stabilité obtenu avec BeamOptics.

Fig. 18 :

 18 Fig. 18 : Comparaison des régions de stabilité et des tunes associés pour BeamOptics (rouge) et Zgoubi (bleu).

:

  longueur de champ de fuite dans l'Eq. 1.5.Les études de balayage sont faites à l'énergie d'extraction, soit 180 MeV.

Fig. 20 :

 20 Fig. 20 : Balayage des tunes autour du point de fonctionnement : ) 78 . 1 Q , 82 . 2 Q ( Z X = =

Fig. 21 :

 21 Fig. 21 : Traversée de la résonance sextupolaire 3 2 2 Q X + = .

  l'acceptance horizontale s'effondre comme la figure précédente le montrait.On peut donc choisir les points de fonctionnement dont les ouvertures dynamiques sont les plus grandes comme candidats possibles pour la conception du prototype d'aimant.

Fig. 22 :

 22 Fig. 22 : Limites de stabilité horizontale (haut) et verticale (bas) dans le diagramme de résonances. Bleu : sextupolaire, vert : octupolaire, violet : dodecapolaire

  Les figures représentant l'anneau et les allures des champs magnétiques ont été faites pour le jeu de paramètres suivant correspondant à une machine pour la protonthérapie à l'étude dans le cadre du projet RACCAM. Certains paramètres seront explicités ci-après.

	Nombre de cellules N											8
	Indice de champ k											4.415
	Angle spiral ξ											50.36°
	Facteur de remplissage pf										0.38
	Rayons mini / max R 1 , R 0										2.8m / 3.55m
	Energie mini / max des protons						17MeV / 180MeV
	Champ magnétique à R 0											1.5T
	Entrefer à R 0											40mm
	Coefficient de Enge C0 … C3							0.15 ; 2.27 ; -0.64 ; 1.16
	forme axiale	e F	(	d	e	),	s F	(	d	s	)	que nous allons expliciter ci-après.
	Le champ magnétique obtenu est illustré à la Fig. 2.

Fig. 2 : Champ magnétique vu par des particules de différentes énergies dans l'aimant FFAG spiral. Tab. 1 : Paramètres du FFAG spiral à l'étude dans ce rapport.

θ

Fig. 3 : Loi de champ magnétique radiale du FFAG spiral.

  De la même manière que précédemment, le programme d'automatisation va effectuer une recherche d'orbite fermée et de nombres d'onde pour chaque . Si le mouvement est stable, alors l'orbite fermée est trouvée et les tunes déterminés par la méthode matricielle de Zgoubi vue précédemment peuvent être calculés.On peut ainsi comparer les résultats du code matriciel Mathematica / BeamOptics et du code de tracé de trajectoire Zgoubi.Pour cette comparaison, les paramètres suivants sont choisis :

	couple	, ( ξ	)

k

Il existe des techniques très rapides, demandant seulement quelques tours, ce qui serait avantageux. Leur application est envisageable.[START_REF]Application de l'analyse en fréquence à l'étude de la dynamique des sources de lumière[END_REF]

Tab. 2 : Paramètres pour la recherche de régions de stabilité.

b) Recherche des régions de stabilité, comparaison avec BeamOptics.

Dans le formalisme matriciel, la stabilité du mouvement périodique horizontal (respectivement vertical) est déterminée au moyen de la demi-trace de la matrice de transfert de la cellule FFAG dans le plan horizontal (respectivement vertical) :

( ξ k peut être estimé de façon approximative par les relations suivantes :

Conclusion.

Le développement d'une procédure de simulation de trajectoires de particules dans un aimant FFAG spiral a été effectué avec succès. Cette procédure permet maintenant d'étudier le comportement dynamique d'une cellule FFAG spiral en déterminant les orbites fermées, les tunes et les limites de stabilité. Les simulations se font par interpolation numérique avec un maillage « flottant » ou par l'utilisation de cartes de champ 2D modélisant le champ magnétique dans le plan médian.

A partir de cette procédure de simulation, l'étude de la dynamique pour différents points de fonctionnement a été rendue possible en élaborant un programme de balayage des paramètres ) , ( ξ k . Cela nous autorise à présent à tracer les régions de stabilité linéaires et les diagrammes de tunes associés afin de choisir des points de fonctionnement éloignés de toute résonance destructrice. Ce programme de balayage permet aussi de déterminer les ouvertures dynamiques pour chaque ) , ( ξ k afin de voir les effets de tels résonances sur l'acceptance de la machine.

Ces outils sont maintenant utilisés quotidiennement dans le cadre du projet RACCAM pour rechercher les points de fonctionnement du futur prototype d'aimant que nous devons construire.

Dans le futur, ces outils seront comparés aux résultats de tracking que nous effectuerons à l'aide de cartes de champ magnétique 3D calculées avec le logiciel TOSCA. Une évolution de ces outils est envisagée afin de leur permettre d'effectuer des études d'évolution des nombres d'onde avec l'amplitude des mouvements transverses des particules (amplitude detuning), des études d'accélération de paquets de particules et des études de défauts.


Joan Manel Ramírez
email: joan-manel.ramirez@u-psud.fr

Qiankun Liu

Vladyslav Vakarin

Xavier Le Roux

Jacopo Frigerio

Andrea Ballabio

Carlos Alonso-Ramos

Enrico Talamas Simola

Laurent Vivien

Giovanni Isella

Joan Manel Ramirez

+ Vladyslav Vakarin

+ Xavier Le Roux

Delphine Marris-Morini

Broadband integrated racetrack ring resonators for long-wave infrared photonics

teaching and research institutions in France or abroad, or from public or private research centers.

Mid infrared (mid-IR) photonics is becoming an interesting research field due to the wide plethora of applications foreseen in different arenas of our society, including high-resolution biosensing, early medical diagnosis, security, free-space optical communications, defense, 3D imaging, and many more [1][2][3]. Among them, the identification of molecules via vibrational spectroscopy in the longwave infrared (LWIR) range using photonic integrated circuits (PICs) has been a subject of intense research due to the expected superior performance compared to other existing approaches [4]. In this regard, silicon (Si) photonics is expected to have a high impact in mid-IR photonics. Indeed, Si photonics provides a reliable and high volume platform that leverages from a mature technology that has been largely employed by now to develop solutions in the near IR wavelength range, notably for telecom/datacom applications [5]. However, the extension of the operation wavelength range of conventional siliconon-insulator (SOI) waveguides towards the mid-IR is not straightforward due to the optical absorption of the buried SiO2 layer for wavelengths beyond λ = 3.6 µm. Different approaches have been proposed to extend Si photonics towards the mid-IR, including the implementation of waveguide engineering strategies to minimize the mode overlap with the bottom cladding [6], the use of pedestal waveguides [7] or the development of suspended membranes [8,9]. While these approaches can be successful below 8 µm wavelength, according to the Si mid-IR transparency window, alternative solutions are required to exploit larger mid-IR wavelengths. For this purpose, germanium (Ge) stands as a promising material owing to its extended transparency up to λ = 15 µm, its CMOS technology maturity, its high third-order nonlinearity and its large refractive index that makes it suitable to implement compact integrated devices with a small footprint. Therefore, germanium-on-silicon has been often chosen as a platform to implement LWIR PICs, as can operate over a wide wavelength range within the molecular fingerprint region (5 µm < λ < 12 µm) with low/moderate propagation loss [10][11][12][13]. Similarly, graded silicon-germanium (SiGe) has also been employed as an enabling integrated platform, as it permits a gradual accommodation of the lattice mismatch between Si and Ge, hence obtaining a Ge-rich SiGe epilayer with good crystalline quality and low defect density [14]. Thus, graded-index SiGe alloys benefit from the advantages of Ge while providing high versatility in terms of mode confinement, as their optical properties can be effectively controlled by careful engineering of the vertical refractive index profile [15]. Recently, we have demonstrated several broadband LWIR photonic components and building blocks using the graded-index SiGe platform, including lowloss and polarization-insensitive waveguides [16,17], Mach-Zehnder interferometers [18], dual-polarization Fourier Transform spectrometers [19] or Bragg grating Fabry-Perot cavities [20]. Furthermore, a preliminary LWIR sensing proof-of-concept has also been realized between 5.2 µm and 6.6 µm wavelengths [21]. In this context, the implementation of micro-ring resonators appears interesting, as it has many utilities in integrated photonics, such as its implementation as i) a notch filter to finely select and reject narrow frequency bands; ii) a compact sensor that effectively increases the interaction length between the propagating mode and the analyte; iii) an elementary unit for optical delay lines; iv) an efficient modulator or v) a nonlinear device, which provides a strong-light matter interaction with the nonlinear material to generate frequency combs around the pumping wavelength [22]. Current development of LWIR ring resonators is still in its infancy, with demonstrations that cover a wavelength span of λ = 2.75 -5.6 µm [23][24][25][26], a range that remains insufficient for the envisioned applications of LWIR integrated platforms [27]. It is worth to note that extending the operation range of ring resonators well into the LWIR wavelength range is not trivial as the use of Ge-based integrated platforms with a relatively low index contrast could compromise the compactness of the ring resonator and ultimately the maximum free-spectral range attainable; In the following, we will show that by a careful design of the resonators LWIR integrated racetrack ring resonators can be implemented.

In this letter, we demonstrate for the first time an integrated LWIR racetrack ring resonator implemented on a Ge-rich SiGe material platform that works at long-wave infrared wavelengths around λ ≈ 8 µm, shows a maximum loaded quality factor (Q) of 3200, offers a high LWIR bandwidth (Δλ) of 900 nm and a maximum extinction ratio (ER) of 10 dB. These results demonstrate that integrated long-wave resonant structures are feasible using Ge-rich graded-index SiGe layers, opening up the way towards advanced LWIR PICs.

Racetrack ring resonators were developed in a 6 µm thick gradedindex SiGe layer grown by Low Energy Plasma Enhanced Chemical Vapor Deposition (LEPECVD) [28]. A linear increase of the Ge concentration along the growth direction from Si up to pure Ge is used to confine light on the upper part of the waveguide. Waveguides are then designed using a rib-like geometry with an etching depth of 3 µm and a width of 4 µm. Further details of this platform can be found elsewhere [29]. Racetrack resonators were designed to operate at λ ≈ 8 µm taking into account several requirements. Firstly, the directional couplers have to provide a power coupling factor (k 2 ) able to compensate for one round-trip loss in the resonator. The gap distance (g) and coupling length (L) have then been chosen as g=1 µm and L=200 µm to ease fabrication process keeping reasonable aspect ratio for the gap definition, while maintaining a compact device footprint and a free-spectral-range (FSR) above the minimum wavelength resolution of the experimental set-up, which is around 1 cm -1 . Figure 1 shows an example of simulation of the Ey component of the even (Fig. 1. a) and odd (Fig. 1. b.) quasi-TM supermodes in the directional coupler. The power coupling factor was deduced for a coupling length of L = 200 µm, as reported in Fig 1 .c. Interestingly, the achieved value, k 2 ≈ 0.2, is compatible with the propagation loss in the resonator cavity of a few dB/cm. A remarkable flat profile of k 2 from 7.5 to 8.5 µm is observed, with a moderate linear increase from k 2 = 0.1 (λ = 7.5 µm) to k 2 = 0.35 (λ = 8.5 µm), denoting the ability of these directional couplers to operate over a broadband at LWIR wavelengths. This phenomenon is attributed to the particular combination of a low material dispersion at these wavelengths and the propagating mode re-adaptation induced by the vertical graded-index profile [15]. Similarly, the bending radius of racetrack ring resonators was chosen considering the trade-off between the reduction of the round-trip loss (by increasing the bending radius) and the minimum FSR value imposed by our experimental setup. Based on these considerations, and taking into account previously reported propagation losses of 2-3 dB/cm [29], 5 different bending radius between R = 50 µm and R = 250 µm were selected.

Racetrack ring resonators were fabricated using e-beam lithography followed by a selective ICP etching. Fig. 2. a. presents a cross-section scanning electron microscopy (SEM) image of the coupling region of a test resonant structure. A top SEM view of a racetrack ring resonator with R = 100 µm is shown in Fig. 2. b. After fabrication, waveguide facets were diced to facilitate butt-coupling of the free-space laser beam. Devices were characterized using an external cavity quantum cascade laser operating in quasi-TM polarization and in pulsed regime (duty cycle of 5 % and repetition rate of 100 kHz), with a linewidth specified to be below 1 cm -1 . A wavelength range from λ = 7.5 µm (1333.3 cm -1 ) to 8.6 µm (1162.8 cm -1 ) was investigated.

The input/output signals were coupled in and out of the chip by means of ZnSe aspheric lenses, obtaining a typical coupling loss of 4 dB/facet. The collected signal was sent through either a high-resolution camera to visualize the spatial profile of the propagating mode or to an HgCdTe (MCT) detector to scan the transmission spectrum. Figure 3 summarizes the normalized quasi-TM transmission of the racetrack ring resonators with increasing radius of curvature in a bottom-up arrangement. The characteristic interferometric pattern of ring resonators is observed, presenting periodic transmission dips that are sometimes partially clipped due to the limited resolution of our setup. Noticeably, the lowest radius of curvature, i. e. R = 50 µm, shows a poor performance with resonance patterns of 4 dB in average and a multimode-like spectrum that deforms the transmission curve. This fact can be attributed to a partial hybridization of the modes, obtaining a non-negligible conversion of the fundamental quasi-TM mode into quasi-TE polarization (in agreement with FEM simulations). Furthermore, a remarkably higher bending loss is to be expected compared to other structures, hence shifting the critical coupling condition to gap values lower than 1 µm. Indeed, we are still able to observe asymmetric resonances for the resonator having an R = 100 µm, which are equally attributed to a mode hybridization effect, although much weaker. Moreover, it is worth noting the superimposed modulation of the interferometric pattern, causing oscillations of the resonance depth, which are typically to the oscillatory behavior of the transmitted signal in the directional couplers [30].

In order to verify the fundamental mode operation of the devices for R > 50 µm, the FSR of the resonances was measured as a function of the wavelength and compared with the theoretical model assuming the following equation:

L n FSR g 2   (1)
Where ng is the group refractive index of the propagating mode, λ is the operating wavelength and L is the round-trip length. Figure 4 presents the dispersion of the FSR for each bending radius, with the theoretical model of equation 1 superimposed assuming an ng = 3.8. As seen, bigger FSR values were measured for lower bending radius, as expected. Moreover, good agreement between the experimental results and the modeling was found, providing evidence that the observed resonances correspond to a single propagating mode inside racetrack ring resonators.

Concerning the transmittance response, resonators with R > 50 µm display deep resonances with an ER up to 10 dB. Figure 3(b) shows a zoomed-in view of a resonance corresponding to the biggest racetrack resonator (R = 250 µm). A lorentzian fit is used to extract the full width at half maximum (FWHM), which will ultimately be used to extract the loaded quality factor:

FWHM Q central   (2)
A maximum loaded quality factor of Q ≈ 3200 is measured at λ ≈ 8.044 µm, which correspond to a FWHM of 0.38 cm -1 . When compared to the laser linewidth used in the experiment (< 1cm -1 ) we can deduce that the measured resonance spectral shape and Q-factor are most probably limited by the experimental set-up. Interestingly, the critical coupling wavelength (around of 8.1-8.2 µm) is compatible with waveguide propagation loss of 3 dB/cm and negligible bend loss for R>150 µm [29]. In conclusion, we have shown the first demonstration of broadband integrated racetrack ring resonators for long-wave infrared photonics that operate around 8 µm wavelength. The resonant structures have been implemented using a 6 µm-thick Ge-rich graded-index SiGe platform that contains a linear increase of the Ge concentration in the vertical direction. Devices show a bandwidth of Δλ ≈ 900 nm, a maximum quality factor of Q = 3200 and an extinction ratio of 10 dB around the critical coupling condition. This work paves the way towards the implementation of multi-target molecular spectroscopic sensors in integrated photonic platforms.

Fig. 1 .

 1 Fig. 1. Even (a) and odd (b) quasi-TM guided supermodes of the directional coupler containing graded-index SiGe waveguides. Waveguide width is 4 µm, the etching depth is 3 µm and the gap between the waveguides is 1 µm. (c) Simulated power coupling factor, k 2 , for a coupling length of L = 200 µm from λ = 7.5 µm to 8.5 µm.

Fig. 2 .

 2 Fig. 2. (a) Cross-section SEM view of a test resonator device in the coupling area. (b) Top SEM view of a racetrack ring resonator.

Fig. 3 .

 3 Fig. 3. (a) Experimental transmission spectra of racetrack ring resonators operating in quasi-TM polarization and arranged according to increasing bending radius. From bottom-up, R = 50 µm (orange curve), R = 100 µm (light blue), R = 150 µm (green), R = 200 µm (violet) and R = 250 µm (grey curve). (b) Zoomed-in view of a single resonance (grey squares), corresponding to a resonator with R = 250 µm, showing a lorentzian fit (red curve) and a calculated quality factor of Q = 3200.

Fig. 4 .

 4 Fig. 4. Measured FSR values as a function of wavelength and the bending radius, being the gray squares R = 250 µm, the light-pink dots R = 200 µm, the green triangles R = 150 µm and the light-blue lozenges R = 100 µm. Straight black lines correspond to the expected spectral response according to equation 1 for an ng = 3.8. Funding. European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (N°639107-INsPIRE). A. Ballabio acknowledges financial support by the TEINVEIN project funded by POR FESR 2014-2020 (ID: 242092). J. Frigerio acknowledges financial support by the European Union Horizon 2020 FET project microSPIRE (ID: 766955).

Acknowledgment. The fabrication of the device was performed at the Plateforme de Micro-Nano-Technologie/C2N, which is partially funded by the "Conseil Général de l'Essonne". This work was partly supported by the French RENATECH network.


Tao Lé

Gérard Bidan

Pascal Gentile

Florence Billon

Catherine Debiemme-Chouvy

Hubert Perrot

Ozlem Sel
email: ozlem.sel@upmc.fr

David Aradilla
email: david.aradilla@uni-goettingen.de

Understanding the energy storage mechanisms of poly(3,4ethylenedioxythiophene)-coated silicon nanowires by electrochemical quartz crystal microbalance

Keywords: SiNWs, PEDOT, EQCM, (pseudo)-supercapacitors

The capacitive properties of electrodes elaborated from the electrochemical deposition of poly(3,4-ethylenedioxythiophene) (PEDOT) coatings onto chemical vapor deposition (CVD)grown silicon nanowires (SiNWs) were investigated for supercapacitor applications. A high areal capacitance value of 17 mF cm -2 (50 F g -1 ) was obtained at a scan rate of 100 mV s -1 in a 3-electrode cell configuration. Furthermore, this 3D hybrid nanostructure was analyzed by electrochemical quartz crystal microbalance (EQCM) to correlate the interfacial ionic exchange mechanisms to their electrochemical performance. It was demonstrated that both anions (BF4 -) and cations (TBA + ) were simultaneously involved in the charge compensation during the oxidation-reduction scans of cyclic voltammetry measurements.

Introduction

During the last decade, doped SiNWs have awakened an enormous interest as promising electrodes in the field of (micro)-supercapacitors (SCs) due to their intrinsic properties in terms of quasi-ideal capacitive behaviour (rectangular shape cyclic voltammograms), wide operating cell voltage (up to 4V), large working temperature range (-40 to 60°C), good technical compatibility with the microelectronics industry (e.g. solder reflow process), high maximal power density (225 mW cm -2 ), excellent coulombic efficiency (~99%) as well as an extraordinary cycling stability (up to millions of galvanostatic charge-discharge cycles) [START_REF] Aradilla | Beyond conventional supercapacitors: Hierarchically conducting polymer-coated 3D nanostructures for integrated on-chip micro-supercapacitors employing ionic liquid electrolytes[END_REF].

Despite of such outstanding electrochemical performances the comprehension of energy storage mechanisms concerning SiNW-based SCs has not yet been explored, which might play a crucial role in unveiling the chemical/physical processes involved during the charge-discharge cycles and thus, in the design of high performance next generation SCs. Within this context, insitu and advanced simulation methods based on cutting-edge spectroscopy, scattering, imaging and electrochemical techniques have recently attracted a great deal of attention [START_REF] Lin | Advanced analytical techniques to characterize materials for electrochemical capacitors[END_REF][START_REF] Salanne | Efficient storage mechanisms for building better supercapacitors[END_REF][START_REF] Forse | New perspectives on the charging mechanisms of supercapacitors[END_REF]. Among them, the electrochemical quartz crystal microbalance (EQCM) is a powerful tool to investigate the ionic exchange at the electrode-electrolyte interface, providing useful information about the capacitive properties of the electrode by analyzing the capacitive current and the resulting mass changes [START_REF] Levi | Quartz crystal microbalance with dissipation monitoring (EQCM-D) for in-situ studies of electrodes for supercapacitors and batteries: A mini-review[END_REF]. To date, in the area of SCs, EQCM has been widely employed in a great variety of carbon (nano)-structures (e.g. carbon nanotubes, activated carbon, carbide-derived carbon, or micro-porous carbon electrodes) or pseudocapacitive material films (e.g. MnO2) [START_REF] Arias | New insights into pseudocapacitive charge-storage mechanisms in Li-birnessite type MnO2 monitored by fast quartz crystal microbalance methods[END_REF]. More recently, we have reported a pioneer study concerning the potential of EQCM in the understanding of the supercapacitive behaviour of vertically-oriented graphene nanosheet (VOGN) electrodes, which had not been reported before in the literature [START_REF] Lé | Unveiling the ionic Exchange mechanisms in vertically-oriented graphene nanosheet supercapacitor electrodes with electrochemical quartz crystal microbalance and acelectrogravimetry[END_REF]. Following these previous developments, in this study, a new approach dealing with PEDOT-coated SiNW electrodes will be investigated, for the first time, by EQCM revealing its excellent electrochemical performance as pseudocapacitor electrode [START_REF] Aradilla | Beyond conventional supercapacitors: Hierarchically conducting polymer-coated 3D nanostructures for integrated on-chip micro-supercapacitors employing ionic liquid electrolytes[END_REF].

Experimental

Growth of SiNWs. Highly n-doped SiNWs were grown in a CVD reactor (EasyTube3000 First Nano, a Division of CVD Equipment Corporation) via gold catalysis using an optimal vapor-liquid-solid deposition procedure [START_REF] Gentile | Effect of HCl on the doping and shape control of silicon nanowires[END_REF]. Briefly, SiNWs were grown from a 12 nm evaporated thin gold film on GaPO4 substrates (specific temperature resistant resonators for the EQCM experiments) using a reaction temperature of 600°C, under 6 Torr total pressure, with Electropolymerization of EDOT on SiNWs. PEDOT films were electrochemically deposited from a propylene carbonate solution containing 10 mM EDOT as the monomer and 0.1M tetrabutylammonium tetrafluoroborate (TBA BF4) as electrolyte. The electropolymerization was conducted in a home-made 3-electrode cell configuration using a multichannel VMP3 potentiostat/galvanostat with Ec-Lab software (BioLogic, France). SiNWs was employed as the working electrode, a Pt grid as counter electrode and a non-aqueous Ag/Ag + reference electrode. The electrochemical deposition of PEDOT was carried out by chronoamperometry using a constant potential of 0.9 V (vs Ag/Ag + ) during 1h. All chemicals were purchased from Sigma Aldrich. where A is the electrode surface (0.2 cm 2 ) and 𝑘 𝑠 is the sensitivity factor (Theoretical value: 7.92 • 10 7 𝐻𝑧 𝑔 -1 𝑐𝑚 2 ) [START_REF] Lé | Unveiling the ionic Exchange mechanisms in vertically-oriented graphene nanosheet supercapacitor electrodes with electrochemical quartz crystal microbalance and acelectrogravimetry[END_REF]. A solution of acetonitrile containing 0.5 M TBABF4 was used as electrolyte for the measurements, with a platinum grid as counter electrode and an acetonitrilebased Ag/Ag + reference electrode. This EQCM setup was employed to perform cyclic voltammetry (CV) combined with simultaneous mass variation measurements in a potential window from -1.5 V to 1V vs Ag/Ag + at a scan rate of 100 mV s -1 . More quantitative information from the EQCM data is obtained with further analysis, i.e. by means of the mass per mole of electron (MPE) estimation. This analysis provides the mass/charge ratio values as a function of applied potential during a CV scan, indicating the apparent mass of the ions involved in the charge compensation. Thus, MPE can be expressed according to the following equation [START_REF] Escobar-Teran | Gravimetric and dynamic deconvolution of global EQCM response of carbon nanotube based electrodes by acelectrogravimetry[END_REF][START_REF] Goubaa | Dynamic resolution of ion transfer in electrochemically reduced graphene oxides revealed by electrogravimetric impedance[END_REF]:

Electrogravimetric characterization with EQCM and evaluation of capacitive

𝑀𝑃𝐸 = 𝐹 ∆𝑚 ∆𝑞 = 𝐹 ∆𝑚/∆𝑡 ∆𝑞/∆𝑡
where F is the Faraday constant and ∆q are the charge variations.

Areal capacitance (AC) was calculated according to the following formula, AC=Q/(ΔVA), where Q is the average voltammetric charge, which is determined by integrating either the oxidative and reduction scans of the corresponding CV curve and ΔV is the potential range.

Results and discussion

Figure 1a and 1b show the morphology of SiNWs grown by CVD and electrochemically PEDOT-coated SiNWs deposited on GaPO4 substrates. As can be seen in Figure 1a, a high density of SiNWs (~10 9 NW cm -2 ) was deposited on the substrate, leading to a total mass of 220 µg cm -2 (determined by the QCM measurements before and after the deposition). The morphological characteristics of SiNWs were analyzed based on the SEM images (Figure 1) and revealed NW length of 10 µm and diameters ranging from 50 to 500 nm. Regarding its surface functionalization by electroactive conducting polymers, a nanometric PEDOT coating (thickness of ~100 nm) with a globular morphology was deposited along all nanowires, which was found to be similar to our previous work dealing with the electropolymerization of PEDOT in presence of ionic liquid electrolytes (e.g. PYR13TFSI) [START_REF] Aradilla | Novel hybrid micro-supercapacitor based on conducting polymer coated silicon nanowires[END_REF]. Thus, according to Figure 1b,a uniform, adherent and conformal polymer coating, corresponding to a deposited PEDOT mass of 125 µg cm -2 , was achieved. The potential of these nanostructures, as capacitive electrodes, was evaluated by CV curves at a scan rate of 100 mV s -1 as displayed in Figure 2. Important differences concerning the working principle between both nanostructures was evidenced. Thus, SiNW electrodes exhibited an excellent electrochemical double layer capacitive behaviour associated to the electroadsorption and de-electroadsorption of ions during the charge-discharge processes (Figure 2a), whereas PEDOT-coated SiNWs reflected a characteristic electrochemical response associated to oxidation-reduction reactions of PEDOT (Faradaic reactions). Thereby, a clear peak at -0.25V (oxidation) vs Ag/Ag + was recorded at a scan rate of 100 mVs -1 (Figure 2b), leading to an AC value of 13 mF cm -2 .This value was clearly found to be higher than that of bare SiNWs (AC: 38 µF cm -2 , Figure 2a), which was expected according to the pseudocapacitive nature of PEDOT.

Concerning the electrogravimetric response, EQCM data could not be obtained in the particular case of bare SiNWs because of the low mass variations. The noise observed in Figure 2a was ascribed to the hydrodynamic contribution and the viscoelastic properties of SiNWs.

However, the SiNWs' electrogravimetric response was modified when covered with PEDOT and consequently significant mass variations of around 2 µg were measured upon cycling with points of zero mass change at 0.1V (oxidation) and -0.25V (reduction), indicating a major anion contribution above 0.1V and cation contribution below -0.25V (Figure 2b), respectively. Above 0.1V the value of the MPE is around 90 g mol -1 , which was found to be in good agreement with the molar mass of BF4 -anions (86.8 g mol -1 ). During a positive scan, these measurements indicate that TBA + cations are expelled from the PEDOT coating during the oxidation peak at -0.25V and over 0.1V, a capacitive process appears, which corresponds to the electroadsorption of BF4 -anions indicated by MPE values.

Conclusions

In this work, for the first time, we have unveiled the mechanisms of energy storage in PEDOT@SiNW-based supercapacitive electrodes by EQCM. The obtained results confirmed the important role of anions (BF4 -) and cations (TBA + ) in the electrochemical performances during the charge-discharge processes. Consequently, EQCM has been demonstrated to be an excellent technique to bring valuable insights into the ionic exchange mechanisms in vertically oriented nanostructured capacitive electrodes for battery and supercapacitor applications. In this direction, further experimental studies will be necessary to adapt and overcome the technical challenges of EQCM for bare SiNWs.  Specific capacitance of 50 F g -1 at a scan rate of 100 mV s -1 .

 Comprehension of mechanisms by electrochemical quartz crystal microbalance.

 Anions and cations play a key role on the electrochemical performance.

  properties. A lab-made QCM device based on a Miller oscillatory circuit was used to measure the resonant frequency change ∆f of the GaPO4 crystal (around 6 MHz), which is then converted into mass changes ∆m using the Sauerbrey equation: ∆𝑓 = -𝑘 𝑠 • ∆𝑚 𝐴

Figure 1 .

 1 Figure 1. SEM images of a) SiNWs grown on a GaPO4/Ti (10 nm)/Pt (400 nm)/Au (12nm) substrate and b) PEDOT-coated SiNWs. c) Schematic illustration of a GaPO4 substrate. Micrographs were recorded by using a Zeiss Ultra 55 microscopy at an accelerating voltage of 4 kV with a tilt angle of 45°.

Figure 2 .

 2 Figure 2. CV coupled with the mass measurement of a) SiNWs and b) PEDOT-coated SiNWs at a scan rate of 100 mV s -1 .

Acknowledgements

This work was supported by the ANR under reference ANR-11-IDEX-0004-02 and by the Cluster of Excellence LABEX MATISSE led by Sorbonne Université.