Da Costa Ramos

Di Meglio

Figueira Da

Reduced Order Model of Laminar Premixed Inverted Conical Flames

An Inverted Conical Flame, anchored at a central bluff-body, in an unconfined burner configuration, is a flame that may be used to further understand more complex flame configurations, such as in aeronautical engines combustion chambers. This article involves the modeling of an unconfined laminar premixed ICF, anchored at a central cylindrical rod. The computational fluid dynamic modeling of this flame configuration involves the solution of transport equations of species mass, momentum and energy, which is computationally expensive. A skeletal methane/air chemical kinetic mechanism is used to capture a plethora of time and length scales. Reduced Order Models (ROM) have been shown to increase the computational efficiency of dynamical systems modeling. Accordingly, in this work a ROM of the steady ICF is developed, using the volume flow rate of the combustible mixture as input variable of the ROM, and the velocity components and the temperature fields as the output. The aim is to optimize the flame modeling computational time, thus paving the way to the reduced order modeling of acoustically excited, unstable, laminar premixed inverted conical flames. A model of the ICF has been developed with Fluent 2019 R2, and the analysis of a diversity of properties and species included in the combustion process has been made. Since a skeletal chemical kinetic model is used, the characterization of the flame is marked by its different property scales, allowing the flame front recognition over the field of CH 2 . However, this CFD model does not present a static convergence behavior, but oscillates over a pseudo-steady state point. Thus an approximation of the statistical steady state has only been achieved by ensemble averaging the results. Therefore, the ROM of the ICF has been developed over a set of averaged data generated with Fluent, and its results shows agreement with CFD results, presenting an average overall error smaller than 3 %.

I. Introduction

C

and acoustic coupling is a major topic of interest for several practical systems, such as domestic boilers and gas turbines combustors, because the flame interactions with flow perturbations may lead to unstable behaviors. The design of such practical systems heavily relies on computational fluid dynamics simulation results. In cases where instabilities may arise, three-dimensional, unsteady simulations which involve dozens of chemical species are often necessary to provide an adequate system characterization [START_REF] Shinjo | Numerical simulation of flame behavior in a lean premixed gas turbine combustor[END_REF]. Even with an ever increasing computational power, performing a large number of such detailed simulations remains unfeasible in the design phase of novel combustors. As a consequence, instability identification could occur late in the design phase, thus leading to expensive reviews.

Lean premixed flames in low emission systems, have been shown to exhibit unstable behavior [START_REF] Birbaud | Dynamics of Confined Premixed Flames Submited to Upstream Acoustic Modulations[END_REF]. Accordingly, studies have been performed aiming to understand the physics of combustion response to incoming velocity/pressure perturbations both in laboratory flames and in combustors [START_REF] Noiray | A Unified Framework for Nonlinear Combustion Instability Analysis Based on the Flame Describing Function[END_REF]. Indeed, investigations on flame instabilities have been carried in turbulent combustors which involve inherent geometrical complexity [START_REF] Lieuwen | Combustion Dynamics in Lean-Premixed Prevaporized (LPP) Gas Turbines[END_REF]. These studies show that mechanisms involving fuel/air ratio fluctuations triggered by pressure oscillations at the injector exhaust, vortical/flame roll up, and interactions with hydrodynamic instability modes are driving sources of unsteady motion [START_REF] Milosavljevic | Pertubation in Combustor Near-Field Aerodynamics as a Main Source of Thermoacoustic Instabilities in Modern Industrial Dry Low NO x Gas Turbine Combustion Systems[END_REF].

However, laminar inverted conical premixed flames feature instability phenomena, representative of more complex situations, but in a simpler flow field framework. Such instabilities have been studied on these flames anchored on a central bluff body [START_REF] Birbaud | Dynamics of Confined Premixed Flames Submited to Upstream Acoustic Modulations[END_REF][START_REF] Durox | Combustion Dynamics of Inverted Conical Flames[END_REF]. In particular, it has been shown that these instabilities are due to velocity fluctuations on the flow, that couples with the flame and generate thermoacoustic instabilities, and that these phenomena can be controlled by forcing the system with specific frequency and velocity (u ). Thereby, studies of inverted conical flames allow to further understand the combustion coupling with acoustic in more complex engineering systems and practical cases, such as gas turbine combustors.

Studies of standard laminar conical flame dynamics instabilities have been performed [START_REF] Schuller | Modeling Tools for the Prediction of Premixed Flame Transfer Functions[END_REF]. However, to the best of the authors knowledge, fewer modeling studies exist concerning inverted conical flames (ICF) anchored on a bluff-body and, also, on reduced order models application to these flames, which is the aim of this paper.

Reduced order model (ROM) is a machine learning process which offers the possibility of significantly decreasing the cost of parametric exploration for complex systems. Such models are built over a few well chosen results in the design space, accurately describing the system of interest, even when large parametric excursions arise [START_REF] Luboz | Personalized Modeling for Real-Time Pressure Ulcer Prevention in Sitting Posture[END_REF].

The use of ROM potentially increases the computation efficiency of the simulation of dynamical systems and it is prevalent in many physical fields [START_REF] Xiao | Non-intrusive Reduced Order Modelling of the Navier-Stokes Equations[END_REF]. A large number of techniques has been proposed for the modeling of different physical cases, such as online and offline models. Online models, which has their learning developed while measurements are being taken, have been previously studied, as for example on Neural Network applied to image recognition [START_REF] Lecun | Deep Learning[END_REF] or nonlinear Autoregressive Exogenous (NARX) Neural Network applied to solar radiation [START_REF] Boussaada | A Nonlinear Autoregressive Exogenous (NARX) Neural Network Model for the Prediction of the Daily Direct Solar Radiation[END_REF], and other statistic models. The offline models, which are developed over measurements that have already been obtained, first decompose the available data, using compressing techniques, such as proper order decomposition (POD) or singular value decomposition (SVD). Then a machine learning approach is applied to learn the studied system. Offline models coupled with POD have been studied on several physical systems, such as on the modeling of pressure ulcer prevention [START_REF] Luboz | Personalized Modeling for Real-Time Pressure Ulcer Prevention in Sitting Posture[END_REF] and thermal behavior of fluid-cooled battery [START_REF] Asgari | Application of POD plus LTI ROM to Battery Thermal Modeling: SISO Case[END_REF].

This paper is devoted to build and to characterize a reduced order model of axisymmetrical, laminar premixed flames with variable dilution by air. Willing to achieve an approximation of the unsteady ICF that has pressure fluctuations varying the flow rate, this reduced order model created has the fuel/air mixture inlet velocity as input variable for the ROM, and the flame velocity components and temperature field as output. This study is thus a first step towards the description of acoustically perturbed flames.

This manuscript is organized as follows; first the inverted conical flame is described and its scales are given. Second, the physical of the ICF and a kinetic chemical model are presented. Then, the computational model described, showing the geometry, mesh, boundary conditions and dynamic mesh adaptation, followed by an explanation of the ROM. Finally, the results from the CFD and from the ROM presented and discussed for the modeling of the static laminar premixed inverted conical flame.

II. Numerical Methodology

The studied inverted conical target flame is presented (section II.A) first, then the experimental setup is given, followed by the flame scales. Then the two used models are outlined, i.e., the computational fluid dynamics model (section II.C) and the reduced order model (section II.D).

A. Inverted Conical Flame

In this paper an axisymmetric flame, depicted at Fig. 1 is studied. In this figure the flame position is visualized by the field of OH mass fraction, which allows to verify that the laminar inverted conical flame (ICF) is anchored by a central bluff-body. Classically, the premixed flame front separates the fuel/air mixture from the combustion products. This flame has been shown to be susceptible to thermo-acoustic instabilities, mainly driven by flame convective vorticity mode [START_REF] Durox | Combustion Dynamics of Inverted Conical Flames[END_REF].

Experimental Setup

Figure 2 shows the dimensions and boundaries of the computational domain. The ICF modeled is based on the experiments of Durox et al. [START_REF] Durox | Combustion Dynamics of Inverted Conical Flames[END_REF], developed on a burner with outlet diameter of D = 22 mm, represented in Fig. 2, fed with a lean methane-air mixture. The premixed fuel enters in the mixture-inlet (Fig. 2), whereas air enters in the air-inlet with a small velocity of v air = 0.3 m/s, so as the air flow in the domain moves downstream. The flame is anchored by a central rod with a diameter of d = 6 mm that is 2 mm higher than the burner outlet, as shown in Fig. 2. In this paper the methane/air flow has a variable volume flow rate ( ∀ d ) and thus a variable mean flow velocity, v d , and constant ambient pressure and temperature of 1 atm and 300 K, respectively. It is important to note that the computational domain (Fig. 2) is divided on two parts; an internal one, that includes the fuel/air inlet tube and the the combustion region, and an external one, which accounts for the air flow and the dilution of the burned gases with air. This division is used to decrease the computational burden of modeling a reactive flow on such sensible flame as the ICF. The inside part of the domain demands high cost, since it is necessary to account for the boundary layers along the bluff-body and the burner walls and, also, the reactions that arise at the flame. The dimension of this inside region have been chosen based on a preliminary trial-and-error procedure, and the boundary and mesh setup of both parts are discussed later on this paper.

Flame Scales

To accurately solve a premixed laminar flame, it is necessary to determine beforehand the associated scales, such as the laminar flame front speed and the thickness for the specific modeled flame. Accordingly, CHEMKIN-Pro has been used to solve an one dimensional steady premixed methane-air flame with an skeletal kinetic chemical model for the combustion (DRM19 kinetic model [START_REF]Reduced Reaction Sets based on GRI-Mech 1.2[END_REF]).

The main features of the obtained solution are as follows: (1) fluid domain spanning from -2 to 5 cm, (2) an adaptive grid control based on the solution curvature and gradient threshold of 0.7 and 0.2 respectively, and (3) the laminar flame front speed S L converged to S L = 40 cm/s. The minimum mesh size is 5.6 µm. This simple modeling with CHEMKIN-Pro using the DRM19 kinetic model, allows the analyses of diverse flame properties and its different scales. However, for the sake of brevity, three of them are discussed here only; temperature, OH and CH 2 mass fraction. This choice has been done based on the large difference between the scale of these properties both on value and lenght, which is seen at Fig. 3, and also due to the characteristics of each of the species shown. Indeed, the OH is often used for experimental flame measurements, and CH 2 is a good indicator of the flame front position and thickness. In order to examine these properties, the fields of the mass fraction of CH 2 , OH and the temperature are given in Figs. 3 and4.

One may observe that, in Fig. 3, CH 2 and OH are absent upstream to 0 cm, where the temperature is 300 K. The temperature then increases up to approximately 800 K and gives rise to the combustion, and significant amounts of CH 2 and OH may be seen. The position of the flame front is highlighted by the increase in CH 2 and OH mass fraction. Further downstream the CH 2 then is consumed and its mass fraction decreases to zero denoting the tendency to thermodynamic equilibrium. The OH, presents a maximum at the flame front, then gradually tends to an equilibrium concentration as x → 5 cm. The temperature increases along the flame until it reaches the adiabatic flame temperature, T ad = 2, 200 K. The mass fraction profiles given in Fig. 4 indicates that the flame reaction zone thickness is around 120 µm. Also, may be inferred that the smaller mesh size required to obtain a well resolved gradient and curvature of the laminar premixed flame properties is of the order 5µm, which is 24 times smaller than the flame reaction thickness.

B. Physical Chemical Model Summary

The studied flames are modeled using the transport equations of species mass, momentum and energy in an axisymmetric reference frame, since axial symmetry is believed to accurately describe the experimental results. The species mass transport equation is solved accounting for the multi-species diffusion coefficients, the Fick's law of diffusion, the Soret thermodiffusion effects, but neglecting the barodiffusion effects since the flame is isobaric. The solved energy transport equation accounts for the non unity Lewis number effects, but thermal radiation is neglected since the modeled flame is a lean premixed flame, which does not presents soot formation. Therefore, the radiative medium emission relies on CO 2 and H 2 O at the burnt gases only, which is significantly smaller than the convective heat transfer contribution. As a consequence, the flame may also be considered transparent regarding the absorption process. The momentum transport equation considers the buoyancy effect, and Newtonian fluid behavior is assumed.

The chemical kinetic process for the lean methane/air mixture is described by a skeletal model. The model used is the DRM19 [START_REF]Reduced Reaction Sets based on GRI-Mech 1.2[END_REF], which is a reduced GRI-Mech 1.2 kinetic model. The GRI-Mech 1.2 model is composed of 32 species and 177 reactions, which can demand a large computational cost. Because of this, a reduced model is used in this work. The DRM19 skeletal kinetic model is composed of 19 (plus N 2 ) species and 84 reactions [START_REF]Reduced Reaction Sets based on GRI-Mech 1.2[END_REF]. Previous assessment of skeletal model accuracy supports this choice [START_REF] Celis | Computational assessment of methane-air reduced chemical kinetic mechanisms for soot production studies[END_REF].

C. Computational Fluid Dynamics Model (CFD)

The numerical model of the inverted conical flame has been developed using ANSYS Fluent 2019 R2. The process of modeling a combustion problem involves the following preliminary steps: the geometry definition and the mesh generation. In this section, the most important steps to the creation of the CFD model are presented.

Since the flame is suposed to be axisymmetric, only a slice of the actual physical domain is modeled, effectively reducing the system to a 2D model, as represented at Fig. 2. The domain is composed of two surfaces with different mesh sizes, which intend to reduce the amount of computational nodes in areas where chemical reaction is absent.

The axisymmetric system has eleven boundaries, as shown at Fig. 2, where two are the rod wall and top (W 4 and W 3 respectively), two are burner internal walls (W 1 and W 2 ) and two are external walls (W 6 and W 5 ), and one represents a slip boundary (Slipboundar y). The remaining boundaries represent the inlets (Airinlet and Mixtureinlet), outlet and axis of symmetry.

The boundary conditions are applied to each boundary of the model, which is represented at Fig. 2, as follows: (1) the inlet boundaries specify the velocity, temperature and mixture molar fraction, specifically the airinlet velocity (v air ) is 0.3 m/s and the mixtureinlet mean flow velocity (v d ) is varied from one simulation to the next, (2) the outlet boundaries specify pressure (P atm = 1 bar) and temperature (T amb = 300 K), whereas the velocity is free to vary, and (3) the wall boundaries conditions for W 1 , W 2 , W 3 and W 4 are no-slip, adiabatic and non catalytic, whereas W 5 , W 6 and Slipboundar y are slip, adiabatic non catalytic walls, since du/dn = 0 and df /dn = 0, f = {T,Y i }. One may note that the airinlet velocity (v air = 0.3 m/s) is chosen based on the Froude number (Fr = v air / √ gL, where g is the acceleration due to gravity and L is the flame length), which might be less or equal one (Fr ≤ 1) to allow for neglecting outer convective effects, such nonzero air velocity value is used to help the simulation convergence, since using zero air velocity led to flashback flow at the outlet boundary.

A characteristic mesh size of 100 µm is prescribed at the combustion region of the domain, i.e., inside and above the burner, and this value have been chosen smaller then the flame front thickness to assure the existence of mesh node inside of the reactive region. However, the analysis of section II.A.2 indicates that a finer mesh is required for the reactive region, indeed using a size of 5 µm throughout the flame front should be necessary, but using this mesh at the whole domain would highly increase the number of mesh nodes and, therefore, the computational burden of the simulation. For this reason, a mesh adaptation tool is applied to the model, allowing the resizing of the mesh at the regions of interest.

The mesh adaption used here is a refining and coarsening procedure used to change the mesh size based on the local temperature gradient modulus. More precisely, at the regions of the mesh where the temperature gradient is higher than 10 K/m, the mesh is refined, and in the regions of the mesh where the temperature gradient is lower than 300 K/m, the mesh is coarsened. This choice of temperature as the mesh adaptation variable is based on a trial-and-error procedure. These thresholds have been chosen by running ICF preliminary computations and analyzing the corresponding results for different meshes. Accordingly, a steady state simulation has been used with mesh adaptation at each 200, 500 and 1, 000 iterations, until the convergence is achieved.

Figure 5 shows a representation of the mesh changes due to the adaption using the temperature gradient. At the left side, the figure shows the base mesh, with a size of 100 µm whereas, at the center, a representation of a mesh already adapted after the flame is developed may be seen. At the right, a zoom on the already adapted mesh is represented. A comparison between Fig. 1 and5, shows that the mesh is refined at the flame front region, where the temperature gradient is high and the combustion process occurs. The initial regular mesh, at the left of Fig. 5, is composed of 509, 473 cells, whereas the mesh, at the last instant of the simulation, presents 2, 427, 310 cells.

D. Reduced Order Model

A reduced order model (ROM) is a simplification of a high-fidelity dynamical model, used for reducing the computational burden or the storage capacity required for complex models, but preserving the essential behavior and dominant effects of the system. The ROM process is a machine learning process which has three phases; learning, validation, consumption. Here, a ROM is computed describing the relation between the mean inlet velocity (or, equivalently, volume flow rate) and the full flame velocity and temperature fields at steady state, using, as input, a dataset generated by CFD models.

To develop a ROM, it is necessary to first choose the snapshots to compose the training set, which is used to model the system, and then divide it in two sets; i.e., the learning and the validation snapshot sets. One may note that the snapshot set is the ensemble of data available for the ROM construction, and each ROM created has its corresponding snapshot set, for instance, snapshots of temperature fields for creating a temperature ROM. Regarding the division procedure, the learning set might consist of 60 % of the available snapshots, distributed through the whole set, thus assuring that all of the ROM input variable range is included at the learning of the system, so as to avoid validation by means of extrapolation. The validation snapshots consists of the remaining 40 % of the original set.

The second step is the application of a singular value decomposition (SVD) algorithm at the learning dataset and, thus, compressing by retaining its higher modes to describe the full field only. Each of the learning snapshots corresponds to a small set of modes coefficients when projected to the reduced SVD base [START_REF] Luboz | Personalized Modeling for Real-Time Pressure Ulcer Prevention in Sitting Posture[END_REF]. Then the learning process of the ROM, is done by using a Genetic Aggregation Response Surface (GARS) to interpolate the previously obtained modes coefficients with respect to the calculation parameterization [START_REF] Acar | Various Approaches for Constructing an Ensemble of Metamodels Using Local Measures[END_REF][START_REF] Viana | Multiple Surrogates: How Cross-Validation Errors Can Help Us to Obtain the Best Predictor[END_REF]. In sequence, the validation computes the Euclidean Norm error between the given training snapshots and the ROM results. Finally, the consumption of the developed ROM to predict the behavior of the system for different values of the input variable, inside of the parametric domain.

Here, four snapshots sets are available for the development of the ICF ROM; the first pair consists of the average velocity component fields, each composed of six or seven snapshots, respectively, whereas the other pair consists of the temperature profiles, also with six or seven snapshots respectively. All of the available snapshots have been developed with steady Fluent simulations of the ICF, with combustible mixture volume flow rates ranging from ∀ 0 = 153.27 cm 3 /s to ∀ 6 = 207.37 cm 3 /s, evenly spaced in the parametric space. Once built, the ROM of the average velocity components and the temperature ROMs should give an accurate approximation of the solution that can be obtained with the CFD solver, but for a different set of volume flow rates.

Upon the ROM process application, the resulting velocity components and temperature fields accuracy depend on the number of snapshots in the training set, the number of modes used in the field compression of the learning dataset, the statistical reduction targeted precision for the flame fields and the interpolation error.

III. Results and Discussion

In this section, results for the modeling of the inverted conical flame with Fluent and for the ROM model are presented. First, the CFD result is given and analyzed, and all of the models that have been developed in this work are listed. Then the StaticROM results for the modeling of the averaged ICF, which are developed using the CFD results as training dataset, are shown and discussed.

A. Laminar Flame Structure

In order to discuss the ICF structure, results have been obtained from the steady axisymmetric laminar premixed inverted conical flame model, with a volume flow rate of ∀ 3 = 180.3 cm 3 /s and a methane/air equivalence ratio of φ = 0.92. Under these reactive mixture flow conditions, the experimental results indicate that a classical premixed laminar flame is anchored at the central rod tip. This flame exhibits an inverted conical shape an gradually consumes the fuel-air mixture. Upon encountering the external air combustion is depleted. At Fig. 6, the model results are presented as a variety of the flame properties. The first three figures (Figs. 6a, 6b and6c) exhibit the temperature and velocity components fields, whereas the others (Figs. 6d, 6e, 6f, 6g, 6h and 6i) show the mass fraction fields of selected chemical species described by the DRM19 kinetic model. These species have been shown in order to highlight the overall flame structure and, in particular, the multiple length scales that arise within the studied premixed flame. Concerning now the mass fraction of CH 4 given, at Fig. 6d, since there is ambient air flowing in the external region of the burner, at the side of the flame that is farther away from the symmetry axis, the fuel air mixture is progressively diluted by the ambient air, whereas it is consumed by the combustion reactions when it is closer to the symmetry axis. Figure 6e gives the field of mass fraction of N 2 , which is present in the air and in the air-fuel mixture, with respective mass fractions of Y air N 2 = 0.7899 and Y mix N 2 = 0.7203. Also, farther away from the symmetry axis of the burner, N 2 mass fraction is the same of the air stream (Y air N 2 = 0.7899), and as it gets closer to the flame front, the mass fraction decreases due to the dilution by CH 4 . In the burner gases region, the air-fuel mixture has a mass fraction of N 2 equal to Y mix N 2 = 0.7203. Downstream the flame front, the N 2 present in the fuel/air mixture is diluted by the burned gases. Figure 6f exhibits the CO mass fraction, which underscores the shape and the position of the flame front, and possibly gives an approximation of the flame thickness. Indeed, this substance is produced as the reactive region of the flame develops, and decreases after the flame front, gradually tending to the equilibrium mass fraction of Y CO = 0.0038, at it is converted to CO 2 . One may note that the field of CO approximate the flame front thickness, since this substance is not completely produced and consumed the combustion region, mixing and reacting with other burned gases downstream the flame, which could give a thickness measurement larger than other more representations.

Figure 6g gives the field of mass fraction of OH, which is produced in the reaction region and could also be an identifier for the position and shape of the flame front. Figure 6h shows the field of mass fraction of CH 2 , which gives a sharp representation of the flame thickness, since it is produced and consumed within the combustion region. The model results indicates a flame thickness of approximately 130 µm, measure obtained with Fluent scaled rule, which agrees with the size of 120 µm presented at section II.A.2, and presenting an error of 8.34 %. In the field of CH 2 it is possible to see that the flame front has a size of approximately 120 µm. Finally, Fig. 6i shows the field of mass fraction of H 2 , which is produced in the combustion and part of it is consumed in the reactive region.

Figure 6 illustrates the various length scales that may result of choosing a kinetic model that accounts for different species and the importance of the adaptive mesh. Included, the above discussed chemical species exhibit different physical scales, which have been accounted for by the chemical reactions represented by the DRM19 chemical kinetic model. Besides, the well resolved flame results is possible thanks to the adaptive mesh procedure, which allows to refine the mesh from a size of 100 µm to a minimum size at the flame region (5 µm), as shown in Fig. 5. This mesh size is able to capture the flame scales since it is a mesh 24 times smaller then the flame thickness.

B. Learning and Validation dataset generation

A total of seven cases of the steady inverted conical flames have been simulated with Fluent 2019 R2. Each simulation took about 15 to 20 days to complete on a Linux-64 SuSE, at ANSYS cluster, Intel(R) Xeon(R) E5-2660 v3, using 108 nodes on 9 machines,. For each case, the nominal mean inlet velocity,v i , of the premixed combustible gases is varied according to the values presented in Table 1. Such a variation of the fuel/air mixture inlet velocity is performed to generate enough ICF data for subsequent use on the ROM development. These values are chosen based on the nominal mean velocity of v d = 2.05 m/s = v 3 , which has been studied in the experimental work [START_REF] Durox | Combustion Dynamics of Inverted Conical Flames[END_REF], by varying the mean inlet velocity of the fuel/air mixture in 5 % multiples of the nominal velocity, v d . On what follows the nominal mixture volume flow rate and the mean velocity values have been denoted as ∀ d and v d , respectively. It should be emphasized that the analysis of the steady simulations shows that model does not converge strictly, but oscillates around a pseudo-steady state point. Indeed, the nine properties presented at Fig. 6 exhibit wrinkles at the flame position, and thismay possibly represent the ICF inherent instabilities, but could also be computational error related to the model. Willing to obtain an approximation of the ICF in its statistical steady state, the ensemble average of the results has been computed, in belief that this average might have a convergent behavior. Figure 7 gives the comparison between a steady snapshot and the ensemble averaged results, for the axial and radial velocity components, and temperature, for the cases with volume flow rates of ∀ 3 = 180.32 cm 3 /s and ∀ 6 = 207.37 cm 3 /s. The purpose here is to illustrate the effect of ensemble averaging on the results. Each of these contours plots exhibit a Fluent steady snapshot result on the left, and the ensemble averaged resultover 1000 iterationson the right side. For the sake of brevity, the results of cases 3 and 6 are compared simultaneously. In Figs. 7a and7d, the axial velocity component (left side) and its ensemble average (right side) are presented. These figures allow to observe that the average flame front tip, i.e., where the fuel/air mixes with the ambient air, is smeared when compared to the steady snapshot. Indeed, in the latter case the wrinkles associated to this region are smaller and smoother than those of the steady snapshot, shown at the left side of Figs. 7a and7d. In Figs. 7b and7e a similar behavior may be seen. Indeed, the radial velocity component presents a rather abrupt change of value on steady snapshot, whereas the averaging has, again, an smoothing effect, rendering the wrinkles weaker and smaller. The temperature fields, shown at Figs. 7c and7f, underscore that these wrinkles at the flame front tip computed by ensemble averaging the results become smaller than those encountered on steady results. The temperature field also highlights that, even though these wrinkles at the flame tip are smeared by the averaging procedure, the flame front still preserves the overall characteristic behavior at the reactive region, such as its anchor point at the bluff-body, the flame front position, shape and length. It should also be noted that for case 6, which corresponds to the highest volume flow rate, ∀ 6 = 207.37 cm 3 /s, the averaging effect is even more remarkable at the flame front tip than for case 3. In particular, for case 6, the ensemble averaged flame front tip averaged velocity is rather blurred and no wrinkles are apparent.

The results shown in Fig. 7 indicate that the adopted averaging procedure allows a smoother approximation of the ICF equilibrium state when is performed a qualitative image comparison with the experimental results [START_REF] Durox | Combustion Dynamics of Inverted Conical Flames[END_REF] presented in Fig. 8. It should be stressed that the flame photograph shown in Fig. 8 exhibits the natural flame luminosity, which could be associated to CH * radicals. Such a chemical species is not presented at the DRM19 chemical mechanism, therefore, a qualitative comparison is possible only. Such qualitative comparison is not as successful as when the steady snapshot is used. These examples illustrate the importance and capacity of computing the average of CFD steady results, ensuring that the obtained model result represents the long term behavior of a system.

C. ROM Analysis

In order to predict the flame behavior, reduced order models have been created using some of the flame properties obtained from the ICF Fluent model results. A total of seven points are available to be used as dataset for the ROM, but the number of data points used on the learning process of this system is varied between four or five points, to show the impact of the dataset size on this process. Here, the ROM input parameter is the fuel/air mixture mean inlet velocity, v i , which values are given, for each of the data points, in Table 1. The ROM then predicts two properties: the flame averaged velocity and averaged temperature fields.

In this study, two models are created for each of the ICF properties. The learning set of the first ROM (denoted ROM4) contains four data points (cases 1, 3, 4, 6 in Table 1), whereas the second on (denoted ROM5) uses five points (cases 0, 1, 3, 4, 6 in Table 1). For each ROM, cases 2 and 5 of Table 1 are used as the validation set. Figure 9 compares the Fluent averaged model results with the ROM prediction of these results, with either four or five learning points, and four SVD modes, for the velocity components and temperature fields. Table 2 reports the errors resulting from the construction of each ROM. A first qualitative comparison between the ROM and Fluent results allow to verify that the overall flow structure is correctly captured when either four or five training points are used. These: important flow structures are; (1) the flame anchoring point at the top of the rod, (2) the angle of the ICF, which is directly related to the laminar flame speed, (3) the location of the smeared flame wrinkles, that marks the dilution region and (4) the outer plume region. More precisely, at Figs. 9d and 9g, the axial velocity, v y , upstream to the flame front is equal for both cases, with a value of v y = 2.34 m/s and to v y = 4.2 m/s immediately after the flame front. However, at the dilution region, the ROM velocity goes from v y = 2.34 m/s to v y = 4.6 m/s, since at this region the flame front is smoother and thicker. Concerning the radial velocity component, v r , at Figs. 9e and9h, abrupt changes are seen on the reactive region, varying from 0.9 m/s upstream the flame, to v r = -0.9 m/s immediately after the flame front, and then slowly increases, achieving a value of v r = 0 m/s further downstream.

The temperature, given at Figs. 9f and 9i, increases from T amb = 300 K, to T ad = 2200 K. Furthermore upstream and downstream the flame front, similar Fluent and ROM behaviors are observed. However, at the dilution region, the ROM flame front is slightly different from the Fluent model. Indeed, the upper region has a flame tip smoother when compared with Fig. 9c. The analysis of the ICF fields predicted by ROM underscores that the temperature is the available predicted property which allows to better identify the flame front position. Indeed, the interface between T amb and T ad is a well defined surface, both for the averaged Fluent and the reduced order models. Therefore, the temperature contour of ROM5, which uses 5 learning data points, for ∀ 5 = 198.35 cm 3 /s, presented on Fig. 10, is used to further understand and compare the flame front position and inclination angle. Figure 10 depicts the temperature contour for the Fluent averaged result in red, and the ROM prediction in yellow. One may note that both curves are anchored at the same position i.e., at the top of the bluff-body. Nevertheless, this figure clearly shows that the Static ROM model, created using five learning data points, is not sufficient to predict the flame front position without error. Indeed, further along the flame, at the dilution region by external air, the red and yellow lines stop overlapping, and exhibit different angles. The measure of these angles in relation to the vertical axis has been approximated by a protractor and are; 28.13°for the Fluent (red) result, and 30.38°for the StaticROM (yellow curves). This value discrepancy corresponds to a total error of 8 % for the ROM result of the flame front angle.

Tables 2 summarize the quantitative prediction capabilities of the ROM developed in this work. For each property field, the corresponding average RMS error of the SVD decomposition of the training set, and Euclidean norm error of the ROM prediction is computed. In Table 2, the prediction of the velocity components shows that the reduction error is zero for the all the properties predicted by ROM4 (ensemble average velocity v and temperature T), which is evident, since the SVD is done over four learning points, and the number of modes used for the ROM is the maximum available, four also. When increasing the number of points, this table allows to verify that the reduction error increases, since the SVD is computed over five learning data points, but only the four largest modes are used on the ROM. Examining now the errors given for the velocity components, in Table 2, the average ROM Euclidean norm error decreases from 2.560 % to 2.470 % when adding only one point to the learning set, which underscores the flame similarities that has been discussed previously for the comparison between ROM4 (Figs. 9d and9e) and ROM5 (Figs. 9g and9h). The temperature prediction average ROM, however, slightly increases from 1.028 % to 1.262 % when changing from four to five learning data points, which also agrees with the high similarity encountered on Figs. 9f and 9i, since the error is less then 2 % for both cases, for the temperature prediction on ROM4 and ROM5.

It should be also emphasized that each of the reduced order models results presented above has been computed in real time, taking around 30 s for the prediction of each created of the ROM. Nevertheless, as underscored in section III.A, a CFD simulation takes about 15 days to compute. This represents a very significant advance on the modeling of flames in an accessible, fast and less expensive way. Furthermore, all the obtained ROM exhibit a rather small error, i.e., less than 3 %, even though a set of 5 data points has been used only.

IV. Conclusions and Perspectives

In this section the conclusions about the steady CFD and reduced order modeling of an inverted conical flame will be presented, and also the next steps towards the dynamic modeling ICF will be listed.

The steady modeling of the inverted conical flame presented on this work, developed with Fluent 2019 R2, shows that it is possible to compute for such complex reactive system, allowing the study of a diverse quantity of species by modeling the chemical of the system with skeletal kinetic models. One may note that the ICF Fluent modeling allows to delve in to the mechanics of flames, enabling a thorough analysis of different physical properties, such as miscellaneous reactive substances (i.e., CH 2 , OH and H 2 ), which can not be achieved experimentally. Also, it has been shown that this flame setup is highly sensible, since the model does not converge to steady point exactly, but oscillates around a pseudo-steady state point. Indeed, the ICF model presents a better convergence behavior when a ensemble average of the steady results is computed. These analysis represents significant advance on the study of combustion, since it yields a characterization of these processes. However, this kind of models have a large inherent cost, that is expressed on the large number of machines necessary for computations and the amount of time required to solve complex physical system.

From the reduced order model analyses, it has been seen that its application on the modeling of the inverted conical flame is not only possible, but also beneficial. Using this type of modeling provides a substantial reduction on the computational burden linked to the modeling of combustion. The ROM developed in this work lead to less than 3 % Euclidean norm ROM error. One may also note that these models are developed using a small amount of data points for learning and validation of the ROM. Since the ICF is a complex and sensible system, the increase of the number of points on the data set could provide even better results.

The successful modeling of the steady inverted conical flame open doors to the development of dynamic study of this system, which is the perspective of this work. Since the creation and validation of the ICF model, either on Fluent and ROM, has now been complete, the modeling of the dynamics of the flame will be developed using ANSYS Software. The next step of this study is first to characterize the structure of the ICF, by means of unsteady Fluent, and use the generated data to train a dynamic ROM.

Fig. 1

 1 Fig. 1 Representative view of the OH mass fraction for the steady inverted conical flame (ICF). Volume flow rate of ∀ d = 180.3 cm 3 /s and fuel/air equivalence ratio φ = 0.92.

Fig. 2

 2 Fig. 2 Computational domain with corresponding dimensions and boundary conditions.

Fig. 3 FieldsFig. 4

 34 Fig. 3 Fields of mass fraction of CH 2 and OH, and temperature (K) obtained for the 1D freely propagating laminar premixed methane/air flame, with equivalence ratio of φ = 0.92.

Fig. 5

 5 Fig. 5 Representation of the mesh adaptation due to temperature gradient (obtained using Fluent). Left: base mesh with a size of 100 µm; center: mesh already adapted after the flame is developed; right: zoom on the adapted mesh.

Fig. 6

 6 Fig. 6 Structure of the laminar premixed inverted conical flame, with a volume flow rate of ∀ d = 180.3 cm 3 /s and a methane/air equivalence ratio of φ = 0.92.

( a )

 a Axial velocity component [m/s] for ∀ 3 = 180.32 cm 3 /s. (b) Radial velocity component [m/s] for ∀ 3 = 180.32 cm 3 /s. (c) Temperature [K] for ∀ 3 = 180.32 cm 3 /s. (d) Axial velocity component [m/s] for ∀ 6 = 207.37 cm 3 /s. (e) Radial velocity component[m/s] for ∀ 6 = 207.37 cm 3 /s. (f) Temperature [K] for ∀ 6 = 207.37 cm 3 /s.

Fig. 7

 7 Fig. 7 Comparison of the flame front structure between Fluent and its average: fields of velocity components and temperature, for two different fuel/air mixture flow rates. Fluent steady snapshot on the left and ensemble average on the right side.

Fig. 8

 8 Fig. 8 Experimental steady inverted conical flame, with equivalent ratio φ = 0.92 and nominal volume flow rate ∀ d = 180.32 cm 3 /s [6].

( a )

 a Averaged axial velocity component obtained with Fluent. v y [-1, 8] m/s. (b) Averaged radial velocity component obtained with Fluent. v r [-2, 5] m/s. (c) Averaged temperature obtained with Fluent. T [295, 2200] K. (d) Axial velocity component obtained with StaticROM with 4 training data points. v y [-1, 8] m/s. (e) Radial velocity component obtained with StaticROM with 4 training data points. v r [-2, 5] m/s. (f) Temperature obtained with Stati-cROM with 4 training data points. T [295, 2200] K. (g) Axial velocity component obtained with StaticROM with 5 training data points. v y [-1, 8] m/s. (h) Radial velocity component obtained with StaticROM with 5 training data points. v r [-2, 5] m/s. (i) Temperature obtained with Stati-cROM with 5 training data points. T [295, 2200] K.

Fig. 9

 9 Fig. 9 Comparison between the Fluent averaged results with ROM prediction for the averaged steady ICF, for a fuel inlet volume flow rate of ∀ 5 = 198.35 cm 3 /s.

Figures

  Figures 9a, 9b and 9c depict the average axial and radial velocity components, as well as temperature of the inverted conical flame with a volume flow rate of ∀ 5 = 198.35 cm 3 /s, modeled with Fluent. Note that, as discussed in Fig.7, a qualitative comparison of the average field yields a smother flame front, with a shape resembling that of the steady experimental results[START_REF] Durox | Combustion Dynamics of Inverted Conical Flames[END_REF]. The Static ROM result of the averaged ICF, for a volume flow rate of ∀ 5 = 198.35 cm 3 /s, is given in Figs. 9d, 9e and 9d for ROM4, which uses 4 training data points, and in Figs. 9g, 9h and 9g shows the results for ROM5, which uses 5 training data points.A first qualitative comparison between the ROM and Fluent results allow to verify that the overall flow structure is correctly captured when either four or five training points are used. These: important flow structures are; (1) the flame anchoring point at the top of the rod, (2) the angle of the ICF, which is directly related to the laminar flame speed, (3) the location of the smeared flame wrinkles, that marks the dilution region and (4) the outer plume region. More precisely, at Figs. 9d and 9g, the axial velocity, v y , upstream to the flame front is equal for both cases, with a value of v y = 2.34 m/s and to v y = 4.2 m/s immediately after the flame front. However, at the dilution region, the ROM velocity

Fig. 10

 10 Fig. 10 Comparison of the contour of the flame temperature for the ensemble averaged steady Fluent result (red) and the StaticROM (yellow) prediction for the ICF case 6, with a volume flow rate of ∀ 5 = 198.35 cm 3 /s.

Table 1 Fuel/air mixture parameters setup for the Fluent modeling of the Inverted Conical Flame.

 1

	Case Number (i)	0	1	2	3 [6]	4	5	6
	∀ i / ∀ d	0.85	0.9	0.95	1	1.05	1.10	1.15
	∀ i [cm 3 /s]	153.27 162.29 171.30 180.32 189.34 198.35 207.37
	v i [m/s]	1.74	1.85	1.95	2.05	2.15	2.26	2.36

Table 2 RMS errors of the SVD and ROM Euclidean norm errors of the learning using 4 modes, for each of the developed models.

 2

	SVD	ROM

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant Agreement nr 766264. This work was performed while L.F. Figueira da Silva was on leave from the Institut Prime (CNRS, France). The authors also gratefully acknowledge the support provided by Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq, under the Research Grant No. 403904/2016-1.


V Arounassalame

M Bouschet

R Alchaar

A Ramiandrasoa

S Bernhardt

C Bataillon

J-P Perez

P Christol

I Ribet-Mohamed

Electro-optical characterizations to study minority carrier transport in Ga-free InAs/InAsSb T2SL XBn midwave infrared photodetector

Keywords: InAs/InAsSb superlattice, XBn, carrier transport, lifetime, infrared detector, bariode, diffusion length, mobility, Quantum efficiency, dark current, midwave infrared

In this communication, we report on electrical and electro-optical characterizations of InAs/InAsSb Type-II superlattice (T2SL) MWIR photodetector, showing a cut-off wavelength at 5 µm. The device, made of a barrier structure in XBn configuration, was grown by molecular beam epitaxy (MBE) on GaSb substrate. At 150K, dark current measurements shows a device in the Shockley-Read-Hall (SRH) regime but with an absolute value comparable to the state-of-the-art. A quantum efficiency of 50% at the wavelength of 3 µm for a 3 µm thick absorption layer is found in simple pass configuration and front-side illumination. Combined with lifetime measurements performed on dedicated samples through time resolved photoluminescence (TRPL) technique, mobility is extracted from these measurements by using a theoretical calculation of the quantum efficiency thanks to Hovel's equations. Such an approach helps us to better understand the hole minority carrier transport in Ga-free T2SL MWIR XBn detector and therefore to improve its performance.

INTRODUCTION

In the past decade, Type-II superlattice (T2SL) infrared (IR) photodetectors based on a barrier design (called XBn) have seen significant improvements of their performance. Several years ago, it was observed that Ga-free InAs/InAsSb T2SL's carrier lifetime was significantly longer 1 than Ga-containing InAs/GaSb's one. This finding triggered intensive developments [START_REF] Ting | Mid-wavelength high operating temperature barrier infrared detector and focal plane array[END_REF][START_REF] Soibel | Mid-wavelength infrared inassb/inas nbn detectors and fpas with very low dark current density[END_REF][START_REF] Ariyawansa | Inas/inassb strained-layer superlattice midwavelength infrared detector for high-temperature operation[END_REF][START_REF] Wu | Mid-wavelength infrared high operating temperature pbn photodetectors based on type-ii inas/inassb superlattice[END_REF][START_REF] Deng | High operating temperature pbn barrier mid-wavelength infrared photodetectors and focal plane array based on inas/inassb strained layer superlattices[END_REF][START_REF] Zavala-Moran | Structural, optical and electrical characterizations of midwave infrared ga-free type-ii inas/inassb superlattice barrier photodetector[END_REF] for IR detectors based on Ga-free T2SL. However, such T2SL system has also its own challenges which need to be investigated, in particular concerning the transport of minority carriers and its relation to design and growth parameters. Transport through the superlattice (SL) in the growth direction is a key element in the design of IR detectors, since the carrier conduction in this direction determines the collection efficiency in the device. Different studies have been realized so far using different theoretical and experimental techniques. [START_REF] Olson | Vertical hole transport and carrier localization in inas/inas 1-x sb x type-ii superlattice heterojunction bipolar transistors[END_REF][START_REF] Casias | Vertical carrier transport in strain-balanced inas/inassb type-ii superlattice material[END_REF][START_REF] Tsai | Study of vertical hole transport in inas/inassb type-ii superlattices by steady-state and time-resolved photoluminescence spectroscopy[END_REF] In this paper, we report the extraction of different transport parameters of an InAs/InAsSb T2SL XBn detector. The XBn device design consists of a unipolar conduction band barrier layer (BL) inserted between n-type absorber and contact layers (AL and CL respectively). The BL with a large bandgap material blocks the flow of majority carriers without impeding the transport of minority carriers (holes here). Dark current measurements show that the device is in the Shockley-Read-Hall (SRH) regime but with an absolute value comparable to the state-of-the-art. A quantum efficiency of 50% at the wavelength of 3 µm for a 3 µm thick absorption layer is found in simple pass configuration and front-side illumination. Combined with lifetime measurements performed on dedicated samples through time resolved photoluminescence (TRPL) technique, mobility is extracted from these measurements by using a theoretical calculation of the quantum efficiency thanks to Hovel's equations. [START_REF] Hovel | Semiconductors and semimetals[END_REF] The mobility values show that miniband transport conduction does not occur and hopping transport is more likely to be the mechanism used by the carriers.

EXPERIMENTAL SET-UP

MATERIAL

In this work, we studied a mid-wave infrared (MWIR) T2SL heterostructures as illustrated in Fig. 1. Figure .1.a, shows the stacking of the MWIR Ga-free T2SL barrier photodetector. The T2SL detector studied in this work was epitaxially grown on 2-inch n-type (Te-doped) GaSb (100) substrate by molecular beam epitaxy (MBE) using a 412 RIBER Machine. From bottom to top, the XBn structure consists of a 400 nm Te-doped (n-type) GaSb buffer layer, which is followed by a 100 nm thick n-type doped InAs/InAs 0.65 Sb 0.35 T2SL and by nonintentionally doped (nid) 3 µm thick AL (536 periods) made of the same T2SL structure. A BL is then made from 100 nm nid AlAs 0.09 Sb 0.91 and finally, the CL of the structure is composed of an 80 nm thick nid T2SL. The AL, CL, and BL are undoped, and the residual doping is expected to be n-type for the two first layers and p-type for the third one, at 3.10 15 cm -3 and 2.10 16 cm -3 , respectively. [START_REF] Zavala-Moran | Structural, optical and electrical characterizations of midwave infrared ga-free type-ii inas/inassb superlattice barrier photodetector[END_REF] More details on the structure can be found in. From epitaxial T2SL structures, circular mesa photodetectors with diameters ranging from 60 µm to 310 µm were fabricated using standard photolithography techniques, they were etched to the bottom of the absorber. Mesa photodetectors were realized by isotropic wet etching using a citric acid solution. No active passivation process has been applied on mesa sidewalls. A polymerized photoresist was spun after the wet etching to protect the mesa surface from ambient air. A metal coating was applied on the n-GaSb substrate and the n-type T2SL cap layer. Such devices were used to perform dark current voltage photoresponse measurements and time resolved photoluminescence (TRPL). In addition to the device samples, dedicated T2SL structures were grown to study material properties (Fig. 1.b). Such samples, consisting of a 3 µm thick InAs/InAsSb T2SL layers sandwiched between two GaSb layers providing carrier barrier confinement, were used to investigate absorption measurements.

EXPERIMENTAL BENCH

To extract material properties as a function of temperature, absorption and minority carrier lifetime measurements were performed on samples (Fig. 1.b) placed in a LN2-cooled JANIS cryostat showed in Fig. 2.b. Figure 2.a and b show the principle used for the absorption measurement. We determined the absorption coefficient through IR transmission measurements using a VERTEX 70v Fourier Transform Infrared (FTIR) spectrometer.

The absorption was found from the ratio of the transmissions measured by the FTIR Spectrometer for two samples, the first one is the substrate only and the second one is the dedicated sample (Fig. 3.5µm]. These measurements were performed in front-side illumination configuration, in which light enters the device through the CL and then the photons not absorbed by the AL are completely absorbed by the substrate. In this configuration, the quantum efficiency is single-pass only. Spectral measurements were performed beforehand using a FTIR Spectrometer. The dark current measurements were realized with the same bench, but without blackbody cavity and by blocking incident flux on the sample.

EXPERIMENTAL RESULTS

Absorption and minority carrier lifetime measurements

Figure 3.a shows the measured absorption coefficient in the temperature range 90K-220K of the InAs/InAsSb absorber region. The absorption was found from the ratio of the transmissions measured by the FTIR Spectrometer for two samples with the same substrate. The coefficient is then deduced using α(λ) = -ln(1 -A(λ))/d. The first sample is the dedicated sample (Fig. 1.b). The second sample is the substrate alone with no epilayers deposited on it. At λ = 3.4µm, the absorption coefficient reaches α = 4800cm -1 at 150K, which is a good absorption value when compared to reported values of InAs/InAsSb with nearby period's thickness [START_REF] Soibel | Temperature dependence of diffusion length and mobility in mid-wavelength inas/inassb superlattice infrared detectors[END_REF][START_REF] Rhiger | Infrared absorption near the bandgap in the inas/inassb superlattice[END_REF][START_REF] Webster | Absorption properties of type-ii inas/inassb superlattices measured by spectroscopic ellipsometry[END_REF] and when compared to InAs/GaSb which reaches an approximate value of 3300cm -1 . [START_REF] Klipstein | Modeling inas/gasb and inas/inassb superlattice infrared detectors[END_REF] The temperature study of the absorption shows only a red-shift of the absorption edge associated with a superlattice bandgap decrease with increasing temperature, the overall absolute value is relatively constant. Lifetime measurements are realized on structure samples will the full XBn design but that were not etched. The lifetime is studied with temperatures (see Fig. 3.b) through time resolved photoluminescence technique. [START_REF] Donetsky | Reduction of interfacial recombination in gainassb/gasb double heterostructures[END_REF][START_REF] Connelly | Direct minority carrier lifetime measurements and recombination mechanisms in long-wave infrared type ii superlattices using time-resolved photoluminescence[END_REF][START_REF] Steenbergen | Significantly improved minority carrier lifetime observed in a long-wavelength infrared iii-v type-ii superlattice comprised of inas/inassb[END_REF] We can observe that the lifetime remains constant in the temperature range 90K-200K with an approximate value of 850 ns and then starts dropping. This temperature independence indicates that minority carrier lifetime is limited by SRH processes as reported in, [START_REF] Olson | Identification of dominant recombination mechanisms in narrow-bandgap inas/inassb type-ii superlattices and inassb alloys[END_REF] in which a similar result was obtained. Improvement in material quality is required to increase the absolute lifetime value.

Dark current and photoresponse measurements

A Keithley 6430 sourcemeter is used to apply bias and measure the current delivered by the photodetectors. Figure 4.a shows the measured dark current density as a function of applied bias for temperatures ranging from 90K to 325K. Under 170K, the shapes of the J-V curves show that the dark current is no longer in the diffusion limited regime. Figure 4.b shows an Arrhenius plot of the temperature dependence of dark current density under -0.4V bias. A fit over the temperature range 170K-250K to the expression Aexp(-E a /k B T) yields an activation energy of E a = 290meV which is close to the bandgap (250 meV at 150K). This slight difference can be explained by the presence of a hole-blocking energy barrier between the absorber layer and the barrier layer, which is discussed below. In the temperature range 110K-150K, the GR regime is identified which indicates that a depletion region is formed at the absorber-barrier interface and extends into the absorber, since it is the less targeted doped region between the two. [START_REF] Bouschet | Influence of pixel etching on electrical and electro-optical performances of a ga-free inas/inassb t2sl barrier photodetector for mid-wave infrared imaging[END_REF] The -0.4V dark current density at 150K is 4.5 × 10 5 A/cm 2 which is approximately 10 times the Rule07. [START_REF] Tennant | rule 07" revisited: Still a good heuristic predictor of p/n hgcdte photodiode performance?[END_REF] Figure 5 shows the experimental QE of the device in temperature range [90K-220K] and under two applied biases V = -0.2V & -0.4V . The detector cutoff wavelength is 5 µm at 150K which is in good agreement with the targeted cutoff wavelength with this design. We can observe that the QE is not fully reached at low temperatures for V = -0.2V , which is probably due to the presence of a valence band offset. At V = -0.4V , the barrier seems to be overcome, the QE increases in temperature and saturates at 200K. That is why, measurements with bias voltage are performed at different temperatures and presented in Fig. 6. High values of QE are obtained compared to the maximum external quantum efficiency that can be reached which is around 55% since no AR coating was applied. Knowing the good absorption properties measured, it is also a sign of good transport properties. The QE increases with applied bias until a certain value is reached and then it starts to "saturate". The behaviour confirms the presence of an unintended barrier in the valence band or valence band offset (VBO). Indeed, as the bias increases the VBO decreases and minority carrier holes can pass the VBO to reach the CL. The VBO is relatively large since the bias, called optimal bias V op , at which saturation occurs is important(V op = -0.3V at 150K for example), which means an important bias has to be applied to completely remove the VBO. Moreover, a peculiar feature is observed : the QE at 0V which is at 0% for low temperatures increases to reach non negligible values (8% at 220K) as the detector's temperature is increased. The thermal energy is not sufficient to help carriers pass the VBO, we suspect a band bending from the depletion region's evolution with temperature. Further investigation is needed to understand this observation. It can be observed that at low temperatures, the QE does not saturate with bias voltage while it is the case at 200K. In order to properly observe this phenomenon, we select the value of QE at λ = 3.4µm and plot them as a function of applied bias in Fig. 6 in the temperature range 77K-220K. The QE increases from 32% at T = 77K and V op = -0.4V to reach 54% at T = 200K-220K and V op = -0.3V . This evolution with temperature is essentially from an increase of the hole minority carrier diffusion length, since absorption measurements showed no significant evolution of the absorption coefficient with temperatures (Fig. 3). The QE increase with bias after V op has been reached is due to the increase of the depletion region width in the absorber. The optimal biases are high compared to other reported values. [START_REF] Soibel | Temperature dependence of diffusion length and mobility in mid-wavelength inas/inassb superlattice infrared detectors[END_REF][START_REF] Rhiger | Carrier transport in the valence band of nbn iii-v superlattice infrared detectors[END_REF] However, the optimal bias can be reduced by adjustments in the AlAsSb composition or doping. At low temperatures the slope of the QE is important after the optimal bias is reached. It means that the diffusion length achieved is not enough and the depletion region helps to collect the minority carriers. Indeed, in the XBn barrier device, photons are absorbed in the absorber region, where electron-hole pairs are generated. Then hole minority carriers should have a sufficient diffusion length compared to the absorber region's length to be able to get past the barrier and be collected by the CL. If the absorber region's length is longer than the diffusion length, some of the photoexcited carriers recombine before reaching the contacts. However, if a depletion region is present, photoexcited carriers are swept into the contact by the electrical field, which increases the QE.

𝑉 𝑏𝑖 = -0,2𝑉 𝑉 𝑏𝑖 = -0,4𝑉

DISCUSSION

In order to extract the diffusion length achieved by the photodetector and confirm our analysis, we decided to simulate the QE of this design with a fully analytic model, using Hovel's expressions. 11 Some assumptions for the calculation have been made which include : Low-level injection, uniform material : spatially invariant parameters, zero excess minority carriers at the depletion edge, uniform absorption coefficient, negligible recombination in the depletion region and negligible majority carrier current. The continuity equation of the current density can then be solved, and the analytical solution for quantum efficiency in the different regions of the bariode can be obtained by using a set of boundary conditions. Figure 7 shows the cross-section of the photodiode used in this model. The n-region consists of the neutral region in the AL, the depletion region is formed in the AL and the p-region is neglected since it forms inside the BL which has very high bandgap compared to AL's bandgap.

Light is incident on the CL layer side (front-side illumination) of the bariode, the continuity equation for the minority holes of the n-type material is :

D p ∂ 2 (∆p) ∂x 2 + µ p E ∂(∆p) ∂x -G p + ∆p τ p = 0 (1)
where D p is the hole diffusion coefficient, µ p is the hole mobility, E the external electric field, G p is the optical generation rate, and τ p the hole recombination lifetime.

The optical generation rate is given by :

-G p (λ) = α(λ)F (λ)(1 -R(λ)) exp(-α(λ)x) (2)
where α is the absorption coefficient, F is the total number of incident photons per unit area per second per unit bandwidth, and R the reflectivity of the air/T2SL interface (no antireflection coating).

Under steady-state conditions, with zero external electric field, Eq. ( 1) can be solved using two boundary conditions : at the junction edge, there is no accumulation of carriers since they are drifted by the electric field inside the depletion region ∆p[x = x cl1 + x bl + x n ] = 0; on the other side, we consider that no surface recombination takes place ∆p[x = x cl1 + x bl + x n + X] = 0, this assumption does not change the shape of the curve. The total QE in the device can be written :

QE = 1 qF (J n + J dr ) = QE n + QE dr (3)
QE dr = (1 -R) × [exp(-αx cl1 ) -exp(α(x cl1 + x n ))] (4)
QE n = (1 -R)αL p (αL p ) 2 -1 × exp(-α(x cl1 + x n )) × αL p - αL p exp(-αX) + sinh(X/L p ) cosh(X/L p ) (5)
with L p the minority carrier hole diffusion length, x cl1 , x bl , x al , x sub are the thicknesses of the contact layer, barrier layer, absorber layer and substrate respectively, X is the thickness of the neutral region in the absorber layer. The Hovel calculated QE is compared to experimental ones at three different temperatures (90 K, 150 K, 200 K) and their corresponding V op in Fig. 8. A very good agreement is found, which shows the accuracy of the selected model.

The hole diffusion length L p is linked to hole mobility µ p and lifetime τ p through the relationship :

L p = k B T q µ p τ p (6)
From Eq.6 we can then deduce the mobility value. Table 1 gives a summary of the values of the transport properties extracted. As it can be observed at 90K, the diffusion length is inferior to the AL thickness, which confirms the important slope of the QE versus the bias in Fig. 6 at low temperature, due to the depletion region extending into the AL. At higher temperatures like 150K, 200K, the diffusion length is superior to the AL thickness, minoirty hole carriers can be collected efficiently. As for the mobility, the values are very small, but they are at state-of-the art values. [START_REF] Casias | Vertical carrier transport in strain-balanced inas/inassb type-ii superlattice material[END_REF][START_REF] Soibel | Temperature dependence of diffusion length and mobility in mid-wavelength inas/inassb superlattice infrared detectors[END_REF] Such low values are reached by the mobility because the hole transport is not through miniband conduction but trap-assisted transport and inter-well hopping transport. [START_REF] Ting | Mid-wavelength high operating temperature barrier infrared detector and focal plane array[END_REF] Mobility was calculated from the experimental values of lifetime and diffusion length. [START_REF] Soibel | Mid-wavelength infrared inassb/inas nbn detectors and fpas with very low dark current density[END_REF] Lifetime was obtained via time-resolved photoluminescence (TRPL) measurement.

CONCLUSION

In summary, we report the performance of a T2SL detector based on InAs/InAsSb T2SL XBn bariode through experimental measurements and modeling. The quantum efficiency model, based on Hovel expressions was presented. We measured the material properties absorption and lifetime and used its results to fit the model to the experimental quantum efficiency data. It allowed us to extract the diffusion length of the minority hole carriers in the absorber layer. We found that the diffusion length is shorter than the absorber layer thickness at low temperatures and comparable or longer for temperatures higher than 150K. We then extracted the mobility values, which are very low due to the transport mechanism not being miniband conduction but inter-well hopping transport.

Figure 1 :

 1 Figure 1: Schematic layout of the epilayers stacking for the structures studied. a) the complete structure for dark current and quantum efficiency studies, b) the dedicated sample for lifetime and absorption studies.

Figure 2

 2 Figure 2: a) Principle of the FTIR transmission measurement in order to get the absorption of the SL absorber layer. b) Photo of the bench used for FTIR transmission measurement. c) Schematic of the illumination of the sample (Fig.1.a) with the blackbody cavity. d) Photo of the experimental bench for quantum efficiency measurement, consisting of the cryostat facing the blackbody cavity.

  T det ↗ (a) Measured absorption coefficient of the InAs/InAsSb dedicated sample (Fig.1.b) for different temperatures. Evolution of the measured minority carrier lifetime with temperature. Error bars are plotted to illustrate the measurement precision with the set-up

Figure 3 :

 3 Figure 3: Results of material properties measurements

  (a) Dark current density (J d ) as a function of applied bias V , measured at temperatures ranging from 90K to 325K. 𝐸𝑎 𝑑𝑖𝑓𝑓 ~290 mev 𝐸𝑎 𝐺𝑅 ~140 mev (b) Arrhenius plot of the dark current density. Diffusion and Generation-Recombination regimes are identified.

Figure 4 :

 4 Figure 4: Temperature dependence of the dark current density (J d ) at V = -0.4V.

Figure 5 :

 5 Figure 5: Evolution of the QE with the temperatures range 77K-220K at fixed voltage. QE is for a frontside illuminated photodetector without antireflection coating.

Figure 6 :

 6 Figure 6: Evolution of the QE at 3.4 µm with the bias voltage in the temperature range 77K-220K. QE is for a frontside illuminated photodetector without antireflection coating.

Figure 7 :

 7 Figure 7: Schematic of the device cross-section for XBn diode. The depletion region is materialized with the region between the dashed lines.

Figure 8 :

 8 Figure 8: Comparison of experimental QE and the Hovel calculated QE presented in this paper for three temperature : 90 K, 150 K, 200 K, and their respective optimal bias V op : -0.4 V, -0.3 V and -0.25 V

Table 1 :

 1 Summary of different extracted transport properties from the analytical QE model. Temperature [K] Diffusion Length 1 [µm] Mobility 2 [cm 2 .V -1 .s -1 ] Lifetime 3 [ns] Diffusion length was extracted from the fit of Hovel equations with experimental QE.

	90	1.55	3.64	850
	150	3.64	12.1	850
	200	7.44	37.7	850
	1

ACKNOWLEDGMENTS

This work was partially supported by the French DGA and by the French National Research Agency (ANR) in the framework of the project HOTMWIR [Contract number N°ANR-18-CE24-0019-02].


Eric Maisonnave

NEMO performance on pre-Exascale processors : NEC SX-Aurora TSUBASA and Fujitsu PRIMEHPC FX700

This note complements a preliminary work , in which the performance of the NEMO oceanonly model was measured on NEC SX-Aurora TSUBASA (20B) and favorably compared to the current scalar platforms. In this document, we resume the study on the same platform and extend it to the sea-ice (SI3) and bio-geo-chemistry (TOP-PISCES) modules of the framework. In a second step, NEMO is also ported and tested on Fujitsu PRIMEHPC FX700, a scalar machine with SIMD vector capability, built from cell phone processors (arm). The NEMO porting on these two new platforms appeared to be easy, but getting optimal performance required more efforts. Interesting performance at low development cost seems to be more affordable on the NEC vector machine. Our result does not allow to determine the level of efficiency we could reach on Fujitsu. The hardware dependency of the comprehensive rewriting needed is another concern on this platform. Beside the never ending race to electricity wasting, we think that there is room for energy efficient and scientifically justified ocean simulations that must be performed on the right equipment.

Table of Contents

In our previous study [START_REF] Maisonnave | NEMO performance optimisation on NEC SX-Aurora TSUBASA , Working Note[END_REF], good vector performance of the NEMO ocean-only model were achieved with a light porting and optimisation effort. The ocean model [START_REF] Madec | NEMO ocean engine[END_REF], initially built in the 90's for vector hardware, was still efficient despite 30 years of modifications. The main limitation found, but less sensible at high resolution, was the vector length bounding to the inner loop length (i.e. number of grid points in longitude) that reduces the possibility of 2D spatial decomposition of the grid in MPI subdomains, which leads to the increase of useless land-only sub-domains. At the opposite, with a limited amount of resources, the better efficiency of our code on the vector platform was obvious.

We propose here to extend this first result by porting and check the vectorisation of the two other modules of the NEMO framework (sea-ice and bio-geo-chemistry). These two modules, more recently added, may not follow the vector compliant coding rules that drove the ocean part building.

In this preliminary work, we also explored the possibility to exploit the full computing node (vector engine for computing + x86 host for disk access) by using the host for disk writing, reserving the vector engine part for what it is made for (computations). This time, we propose to use the official XIOS output server to be able to deliver at the end a standard version of the NEMO code.

Finally, to have a larger picture of what the new technologies are currently proposing, the performance of our ocean model is also tested on a scalar A64FX arm v8.2-A processor. This processor equips the top 500 leader Fugaku Fujitsu PRIMEHPC supercomputer 1 and seems to allow geophysics code to be used with reasonably good performance [START_REF] Dongarra | Report on the Fujitsu Fugaku system[END_REF].

Again, we repeat that this study is not a comprehensive machine benchmarking with NEMO, but a simple attempt to evaluate how the vector potential of our code has evolved.

1-NEC SX-Aurora TSUBASA 1.1 -SI3 and TOP-PISCES performance

For this study, we used one computing node (8 Vector engines -VE-, 64 vector cores) of the vector partition of machine (nesh) owned by the Christian Albrecht University at Kiel 2 . The same NEMO 4.0.4 version that was previously used for ocean only testing is taken as reference to start our study. Its include the modifications already implemented to compile the ocean code on NEC SX-Aurora TSUBASA. This allows to setup (compiling and running) two ocean/sea-ice and ocean/bio-geo-chemistry configurations without any further modification. Our ocean code also includes the enhancements necessary to be vectorised at its best. The realistic BENCH-1 [START_REF] Irrmann | Improving ocean modeling software NEMO 4.0 benchmarking and communication efficiency[END_REF] configuration is preferred, to simplify the input file setup and to avoid any input/output perturbation during the measurements. We first measure the vector performance of the ocean only3 on one node of our machine. Its speed is equal to 380 SYPD. The whole optimised code spends 99% of its time in vector sectors and shows an average vector length (AVL) of 175. It approximately fits the number obtained by B. Lodej in our previous study (315 SYPD, 95% and AVL=177). To get these numbers, the code is compiled as previously with a simple -O2 optimisation. At the opposite of our first study, because this option leads to out of bound errors at runtime on the CAU machine, the key_vectopt_loop pre-processing option is switched off in our configuration. For other details about measurement strategy and details see [START_REF] Maisonnave | NEMO performance optimisation on NEC SX-Aurora TSUBASA , Working Note[END_REF].

The out-of-the-box results are poor, which makes mandatory the checking of the vectorisation, as already done for the ocean code. NEMO is compiled and linked with the -report-all option. On diagnostics files, vectorisation is detailed for each line of the code and can be optimised through directives. Notice that the total amount of lines that we checked in the code is limited by the namelist parameters chosen (the parameters of the official release). A more comprehensive work would be mandatory to ensure the NEC vector compatibility in all configurations.

For sea-ice, three problems are preventing the full vectorisation of the code:

1. The recurrent use of the SIGN FORTRAN instruction. In several time consuming routines (dyn_rhg_evp, dyn_adv_pra, limthd, icethd_dh), the instruction is replaced by the corresponding IF condition, 2. The partial definition of an array (IF without ELSE) in dyn_rhg_evp, 3. The dependency, in a computation of a n index variable, to the n+i index variable ( icethd_zdf_bl99 ). This remains unsolved (and takes approximately 14% of the whole simulation) but the algorithm seems to be also problematic on scalar machine. A general rewriting is expected before proposing a fully vector suitable solution.

In total, we increase the AVL (ocean+sea-ice) from 155 to 160, the time spent on vector sectors from 88 to 97% and the simulation speed from 59 to 161 SYPD.

In TOP-PISCES, in addition to the SIGN issue already mentioned (traadv_mus, trc_sink2), we found inner loops not done on the longitude index but rather on tracers (p4zpoc) or vertical levels (trc_sink2). On p4zpoc, it was possible to re-organise the different loops and ensure long enough vector length by using ji index in the inner loop. The same rewriting was impossible for trc_sink2, which stays for that reason at 96% of vectorisation and AVL=82. We found that various less time consuming routines would benefit for an index re-ordering. Again, this problem also needs to be addressed for routines not called in our parametrisation. This index re-ordering issue seems more problematic in the BGC code than in the sea-ice part.

In total, we increase the AVL (ocean+BGC) from 81 to 170, the time spent on vector sectors from 77 to 96% and the simulation speed from 4 to 45 SYPD.

Considering these gains, for both sea-ice and BGC parts, it seems clear that the vectorisation of much of the NEMO routines is mandatory to take the most of the machine. This implies that (i) the whole NEMO code has to be tested with all parametrisations before being able to certify the full model compliance with NEC SX-Aurora TSUBASA machine and (ii) a reordering of loop indexes must be done to increase the AVL of some routines (particularly in TOP-PISCES) that are currently bottlenecks for performance.

-XIOS porting

Despite our efforts, it was not possible to include XIOS in our coupled configuration. The main reason of this limitation was the unavailability of the external libraries (with appropriate version) required by XIOS (boost, netCDF, hdf5) on both VE and VH processors. Without support of the CAU help desk, it was impossible to install the full set of required libraries. This strongly emphasises the lack of portability of the XIOS IO suite and prevents further attempts to install the full ECHAM or OpenIFS / NEMO coupled model on the Kiel vector machine. However, a new try will be given at DWD, where more libraries seems to be available.

Fujitsu PRIMEHPC FX700

A short allocation is provided by a French computing center for public research (TGCC) on a special partition of their top leading machine 4 . This partition includes 80 nodes of arm FX64 processors 5 . This scalar architecture implement the Scalable Vector Extension (SVE), an SIMD extension that allows vectorisation on up to 8 double precision real. The bandwidth per core is rather large (21Gb/s per core) but to get the most of it, the particular memory access technology supposes a fine tuning of the threads mapping, that NEMO cannot provide due to the lack of OpenMP parallelism. However, we propose to test the model in its standard MPI configuration, assuming that a further implementation of OpenMP could give even better performance (on this topic, see [START_REF] Maisonnave | NEMO 4.0 performance: how to identify and reduce unnecessary communications[END_REF]) . The same NEMO 4.0.4 version than on NEC SX-Aurora TSUBASA is used, without any IO or additional module (no sea-ice, no BGC). This is not the first time that NEMO performance was evaluated on an arm processor [START_REF] Ast | arm CAVIUM THUNDERX2 benchmarks[END_REF] but this previous attempt, on a Cavium ThunderX2, was led with the Intel compiler instead of the native arm compiler, not available at that time. Performance just similar to the Intel Broadwell reference was observed. Similarly, a recent test on Fujitsu A64FX preferred the GNU compiler [START_REF] Banchelli | Cluster of emerging technology: evaluation of a production HPC system based on A64FX[END_REF] and gave deceptive results.

In this study, we propose to use the native Fujitsu compiler 6 with basic optimisation options 7 .

No code modification was necessary for porting, except the adding of the key_nosignedzero pre-processing option.

Figure 1: BENCH-1 scalability and communication time on Fujitsu A64FX armv8.2-A and Intel Ice Lake Xeon Gold 6140 @2.30GHz as reference

We present in Fig. 1 the total time to solution (red line) to complete 100 iterations of the main time loop, for several MPI decompositions in a single node. One resource (core), allocated on a single node, is attributed to each subdomain. We call waiting time (orange line) the time spent by NEMO waiting for the horizontal boundary conditions (halos), communicated from their neighbors subdomains. For comparison, the performance of the same model on one Intel Ice Lake node is plotted in blue. The vector -xCORE-AVX512 option is activated on this machine. Both time to solution and waiting time (mainly load imbalance and time spend in MPI communications) are the same, which does not show a significant improvement but validates the good behavior of our model on the Fujitsu-arm processor, without many porting effort.

To go further, we try to use more than one computing node. Unfortunately, it was not possible to get reasonable performance or even to start simulation on more than 2 nodes of the machine. This can be explained by the relatively unconventional parallel strategy we deploy (full MPI), known to be inefficient on this kind of hardware. In addition, the interconnect network is not the native "Tofu" Fujitsu but the Infiniband EDR (100Go) rather similar to the Infiniband HDR already in use on the TGCC supercomputer. This could also contribute to significantly downgrade our performance. However, we plotted in Fig 2 . the performance measured on two Fujitsu nodes, and compared it with the performance measured on the new Intel Ice Lake platform installed at CERFACS, which again are in the same range.

Perspectives

The NEMO porting on these two new platforms appeared to be easy, but getting optimal performance required more efforts. On NEC vector engine, in addition to the already mentioned vector length related limitation for an optimal MPI decomposition, our tests showed additional bottlenecks in sea-ice and biogeochemistry modules. On Fujitsu, an optimum use of the memory hierarchy still need to be done, but requires a specific NEMO implementation like tiling or OpenMP parallelism.

Interesting performance at low development cost seems to be more affordable on the NEC vector machine. Our result does not allow to determine the level of efficiency we could reach on Fujitsu. The hardware dependency of the comprehensive rewriting needed is another concern on this platform. On a slightly similar architecture [START_REF] Ye | swNEMO_v4.0: an ocean model NEMO for the next generation Sunway supercomputer[END_REF], a very different strategy (parallelism on the vertical levels, asynchronous communications) was recently deployed, with a totally different impact on the code. To conclude, our knowledge were insufficient to be able to establish the power consumption of our model on these platforms. Unlike on AMD or Intel [START_REF] Maisonnave | Energy consumption of the NEMO ocean model measured with the Energy Scope tool[END_REF], we lack here of the appropriate measurement software, that would also tell us how vectorisation [START_REF] Guermouche | Thermal design power and vectorized instructions behavior[END_REF] helps to save energy on NEC machines.

Despite these limitations, we keep thinking that the NEC machine is one of our best current choice for quick NEMO computations, particularly if the configuration resolution allows the use of small size platforms, where interconnection speed is not a bottleneck but benefit can be taken from the high processor throughput. Beside the never ending race to electricity wasting, we think that there is room for energy efficient and scientifically justified ocean simulations that must be performed on the right equipment.

1 -

 1 NEC SX-Aurora TSUBASA.................................................................................................................5 1.1 -SI3 and TOP-PISCES performance..........................................................................................5 1.2 -XIOS porting...............................................................................................................................7 2. Fujitsu PRIMEHPC FX700...................................................................................................................7 3. Perspectives.........................................................................................................................................9 References..............................................................................................................................................11

Figure 2 :

 2 Figure 2: Same than Fig. 1. but on more than one node

https://www.top500.org/system/179807/

nesh is an hybrid computing system, based on a scalar NEC HPC Linux Cluster including GPUs and a NEC SX-Aurora TSUBASA vector system (8 NEC SX vector nodes, 10B first generation technology)

ORCA1 global grid, defined on 362x332x75 points

https://www.top500.org/system/179700/

A64FX armv8.2-A SVE @1.8Ghz, 48 cores per node

frt 4.6.1

-Free -Kfast -CcdRR8 -zcfc=target_sve, with -Kfast equivalent to -O3 -Komitfp,mfunc,eval,fp_relaxed,fz,fp_contract,ilfunc,simd_packed_promotion

Acknowledgments

Thanks to Wonsun Park and Mojib Latif (GEOMAR) for providing the access to the Kiel Christian-Albrecht University (CAU) supercomputer, to Jens-Olaf Beismann (NEC-Germany) and Isabelle d'Ast (CERFACS) for their respective support on the CAU NEC/CERFACS Intel Ice Lake machines, to Vera Maurer (DWD) for providing the NEC compatible libraries for XIOS. This work was granted access to the HPC resources (Fujitsu) of Saclay TGCC super-computing center under the allocation AP010113157 made by GENCI and to the HPC resources (NEC) of the Computing Center of the CAU. The author wishes to acknowledge Thomas Williams and Colin Kelley for the development of the Gnuplot program, which analysis and graphics are displayed in this report, in addition to graphics from Matplotlib, a Sponsored Project of NumFOCUS, a 501(c)(3) non profit charity in the United States. This project did not received funding from the European Union's Horizon 2020 research and innovation program.


Douzal Chouakria

Nacer Hammami

Catherine Garbay
email: catherine.garbay@imag.fr

Local Factorial Analysis of Multivariate Time Series

Introduction

It's quite often that we are faced to datasets where a priori relationship structure is defined on the statistical units. Such structure is generally represented by a neighborhood graph, where vertices are statistical units and the edges connect neighboring units. A pioneering work to include a priori structure in factorial analysis was proposed by [START_REF] Lebart | L'analyse statistique de la contiguïté[END_REF] and [START_REF] Banet | Local and partial principal component analysis and correspondence analysis[END_REF]) with the local and contiguity analysis which mainly consists in decomposing the total variance through the neighboring and non neighboring units. Many other studies follow, for instance, Le [START_REF] Foll | Pondération des distances en analyse factorielle[END_REF] generalizes the local analysis to weighting neighborhood graphs. [START_REF] Mom | Multivariate analysis of spatial patterns: a unified approach to local and global structures[END_REF] proposes a new operator generalizing the discriminant factorial analysis to account for a priori neighboring structure. Escofier and [START_REF] Benali | Smooth factorial analysis and factorial analysis of local differences[END_REF] propose an interesting work on factorial smooth analysis and factorial analysis of local differences. Thioulouse et al. (1995) propose another decomposition of the total variance into local variance and global covariance where units are weighted and centered according to the neighboring relationship structure. This paper focuses on the local analysis for the exploratory of multivariate time series partitioned into a priori classes. The studied structure includes two kinds of relationship: the temporal relationship between observations of a same time series, and the membership relationship between time series belonging to a same class. According to that particular structure, we explore two local analysis approaches based on two distinct neighborhood matrices.

Local Analysis

Notation

Let X = [x ij ] be an (n × p) matrix describing n individuals by p variables. A neighboring structure is defined on the set of n individuals, it is represented by a neighborhood graph, where vertices identify individuals and edges connect neighboring individuals. Let M = [m ij ] be a symmetric (n × n) matrix describing the neighboring relationship between the n individuals: m ij = 1 if i and j are neighbors; 0 otherwise. m i,i is assumed equal to 0. Let N = diag(n 1 , ..., n n ) be a diagonal matrix where n i = n j=1 m ij defines the number of i's neighbors. We consider m = i,j m ij = n i=1 n i as the total number of couple of neighbors.

Local variance/covariance matrix

Let V = [v jj ′ ] be the classical empirical covariance matrix of X. To include the neighboring structure Lebart proposes a decomposition of the covariance matrix into two covariance matrices: the first one, called local covariance, is defined on the neighboring individuals and the second one is defined on the non neighboring individuals. The local covariance matrix

V L = [v L jj ′ ] expression is: V L = 1 2m X ′ (N -M )X with v L jj ′ = 1 2m i,i ′ neighbors (x ij -x i ′ j )(x ij ′ -x i ′ j ′ )
Note that V L is identical to V in the case of a complete structure where any couple (i, i ′ ) are neighbors.

Local analysis consists of the diagonalization of V L to look for the first principal components maximizing the local variance. On the contrary, the diagonalization of V L V to look for the axis minimizing the local variance while maximizing the variance due to the non neighboring individuals defines a extension of the discriminant factorial analysis to include a priori structure. It becomes the classical discriminant factorial analysis in the case of a neighborhood graph composed of disjoint cliques where edges connect individuals belonging to a same class.

Exploratory of multivariate time series

Let's consider a matrix X describing n multivariate time series S 1 , ...S n through p variables X 1 , ..., X p observed at the instants of time t 1 , ..., t k . We assume the set of time series partitioned into C a priori and C(S i ) denoting the S i 's class.

Neighboring structure

According to the above dataset X, we distinguish two a priori relationships: the temporal relationship between the observations of a same time series, and the membership relationship between the time series belonging to a same class. This structure can be represented by a neighborhood hyper-graph, where an hyper-node, representing a time series, is described by a neighborhood temporal graph connecting the observations of a same time series, and an hyper-edge, connecting two neighboring time series, is described by a set of edges connecting the observations of the two time series in question. Let's define the general form of the corresponding neighborhood block matrix M = [M ij ]:

S 1 . . . S n M = S 1 . . . S n    M 11 . . . M 1n . . . . . . . . . M n1 . . . M nn   
where M ij = [m ls ] (i = j) is the (k × k) symmetric matrix defining the membership relationship between the time series S i and S j , and M ii = [m ls ] the symmetric matrix defining the temporal relationship between the observations of S i .

Exploratory analysis purposes

In this paper we propose two exploratory analysis of the set of multivariate time series. The first exploratory analysis purpose is to look for the first axis best separating the time series belonging to different classes; while making closer the time series of a same class. These axis should give back of the central position of the time series not of their temporal behavior. The second exploratory analysis, similarly to the first analysis, looks for the first axis best separating the time series belonging to different classes, while making closer the time series of the same class, but with preserving the temporal behavior characterizing each time series.

With respect to the first exploratory analysis purpose, we propose the following neighborhood block matrices M = [M ij ] definition:

M ij = zero matrix if C(S i ) = C(S j ) unit matrix if C(S i ) = C(S j )
Let V L be the local covariance matrix based on the neighborhood matrix M . The first analysis, consists in diagonalizing V L V to look for the first principal components maximizing the local variance. This local analysis will reveal the main directions separating well the time series belonging to different classes, while bringing closer the series of a same class and the observations of a same time series. As M ii is a zero matrix, factorial projections will reveal some central position of each time series instead of their behaviors. According to the second exploratory analysis purpose, let's consider the following neighborhood block matrices

M 1 = [M 1 ij ] and M 2 = [M 2 ij ] definitions: M 1 ij =    zero matrix if C(S i ) = C(S j ) and (i = j) m ls = 1 if i = j and s = l ± 1; 0 otherwise unit matrix if C(S i ) = C(S j ) M 2 ij =    Identity matrix if C(S i ) = C(S j ) and (i = j) zero matrix if i = j zero matrix if C(S i ) = C(S j ) Let V 1 L , V 2
L be the local covariance matrices based on the neighborhood matrices M 1 and M 2 respectively. To bring out as well as possible the temporal behavior of each time series, we propose to maximize the local covariance V 1 L defined for one on the neighboring temporal observations, and to minimize the local variance V 2 L defined by the time series of a same class, based on their values observed at the same time. Thus, the second exploratory analysis, consists in diagonalizing

V 1 L V 1 L +V 2 L to look for the first principal components maximizing the local variance V 1 L while minimizing the local variance V 2 L .

Application and results

To illustrate the above two variants of the discriminant factorial analysis, we consider the Auslan dataset (http://kdd.ics.uci.edu/databases/auslan/auslan.data.html) describing a set of hand signs language of the deaf Auslan community. Instances of the hand signs were collected using an instrumented glove. Position information is calculated on the basis of ultrasound emissions from emitters on the glove. Each sign is measured through a time series described by 8 variables: x (left/right), y (up/down), z (backward/forward), and roll (the palm pointing up or down) giving the spatial position of the glove, and the four remaining variables thumb, fore middle and ring finger indicating if a finger is fully bent or not. For clarity reasons, we have limited our dataset to the 5 signs "Eat", "Boy", "You", "Wild" and "Pen" of the Adam's session 2, with 6 instances per sign. The 30 collected time series were resampled to 57 points (the maximal observed length). We have processed the two above exploratory analysis on the 30 time series. Figure 1 illustrates, on the left the coordinates on the first factorial plan of the 30 time series according to the first exploratory analysis, and on the right the factorial coordinates obtained through the second exploratory analysis. We use a three-dimensional visualization to show the temporal behavior of the projected time series. We can see that both exploratory analysis succeed to separate well the 5 time series classes, with projecting very closely the time series of each class. Let's focus on the main differences between the two produced factorial plans. Based on the temporal neighboring relationship, the local variance of a projected time series is maximum if the behavior of the projected time series is as close as possible to the initial time series behavior (i.e. behavior in the initial descriptive space). On the contrary, the local variance, of a projected time series is minimum if the projection sub-space is some what orthogonal to the initial time series behavior, where all the

Conclusion

This paper focuses on the exploratory analysis of a set of multivariate time series partitioned into a priori classes. We propose two exploratory analysis. The both exploratory analysis look for a sub-space separating well the time series of different classes; while making closer those belonging to a same class. The main difference, is that the first exploratory analysis looks for a subspace working out the main central position of each time series regardless to their behaviors, whereas on the contrary, the second analysis looks for a subspace projecting as well as possible the behavior characterizing the time series. Based on the obtained discriminant spaces, future work will focus on the classification approaches and on the choice of the appropriate proximity measure according to the two above analysis.

Figure 1 :

 1 Figure 1: First principal components of the 30 time series