Romain Schuster
email: romain.schuster@irstea.fr

Dominique Heitz

Etienne Mémin

Motion Estimation Under Location Uncertainty, Application To Large-Scale Characterization Of A Mixing Layer

Keywords: Large-scale PIV, Optical Flow, Stochastic Modeling, Mixing Layer

published or not. The documents may come

L'archive ouverte

Introduction

The increasing necessity to understand and control indoor airflows requires the development of new tools or methods that can describe exhaustively flow features inside buildings. To develop such methods we naturally think about performing particle image velocimetry (PIV) inside buildings. As a matter of fact PIV is a very performant method developed and used in research laboratories to characterize fluid flows. A detailed description on advanced PIV principles can be found in the text book [START_REF] Adrian | Particle-imaging techniques for experimental fluid mechanics[END_REF]. However PIV has been designed and optimized to study flows under control. Its transposition to study free flows inside buildings is not immediate.

Several difficulties appear when we try to perform it in situ in real conditions. Firstly, the dimensions of fields that we want to investigate inside buildings are generally really bigger (≥ possible to increase the number of cameras as already tested by [START_REF] Lee | Large-scale PIV measurements of ventilation flow inside the passenger compartment of a real car[END_REF]. Actually it will make the set-up more complex and we aim at making it simple. Another option is to move backward the camera in order to extend the investigated area. Naturally by doing this the resolution will decrease. Particles within the image will appear smaller and the scattered light will be weakened. To overcome this problem, bigger particles like helium filled soap bubbles (HFSB) can be used. As a matter of fact these particles are brighter and follow the flow with a high fidelity [START_REF] Scarano | On the use of heliumfilled soap bubbles for large-scale tomographic PIV in wind tunnel experiments[END_REF]. This tracer is becoming widely used for large-scale applications. Another option is to use dense smoke seeding. By following smoke density gradients movements it is possible to retrieve quite well the velocity. In both cases the combination of this type of tracer with an expansion of the field of interest leads to a large-scale observation of the flow. Then motion estimation from image sequences has to be adapted.

Usually velocity fields are extracted from image sequences by performing cross-correlation between successive images. It is a robust technique but it only provides large-scale velocity fields due to the size of interrogation windows. Moreover cross-correlation is not suited for smoke images, especially when luminance gradients are weak [START_REF] Heitz | Variational fluid flow measurements from image sequences : synopsis and perspectives[END_REF]. Optical flow is a motion estimation technique which comes from the computer vision community and was first introduced by [START_REF] Horn | Determining optical flow[END_REF]. Initially it was based on the hypothesis of luminance consistency during motion. This hypothesis gives rise to the optical flow constraint equation (1) also called data term that links luminance 𝐼(𝒙, 𝑡) variations to the motion field 𝒖(𝒙, 𝑡) such that

𝜕𝐼 𝜕𝑡 + 𝒖 ⋅ 𝛁𝐼 = 0 (1)
Optical flow technique computes a dense motion field (1 vector/pixel) so the resolution is higher compared to cross-correlation. When the brightness consistency hypothesis is well suited for rigid motions, it is not enough elaborate to cope with fluid motions. Then several attempts to improve fluid motion estimation using an optical flow approach were performed and have been reported in the review of [START_REF] Heitz | Variational fluid flow measurements from image sequences : synopsis and perspectives[END_REF]. Still with the goal of providing a more accurate estimation of fluid motion we have to consider a specific aspect of fluid flows which present a multi-scale behavior. As a matter of fact, due to camera frame rate and size of the field of view the spatio-temporal resolution is often too weak to see the small scales of the flow. This is even more true when smoke is used as a tracer. We cannot see unresolved scales but we have to model it to have a better estimation of coarser ones. A few works has been dedicated to the modeling of subgrid scales for motion estimation. This field was first adressed by [START_REF] Cassisa | Subgrid scale formulation of optical flow for the study of turbulent flow[END_REF] who introduced a data term based on passive scalar transport equation where the action of [START_REF] Bardina | Improved subgrid-scale models for large-eddy simulation[END_REF].

In this paper we propose a new approach to derive scales interactions and improve motion estimation. This approach is based on [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] work who decomposed the Eulerian flow velocity field into a large-scale smooth component and a rough small-scale turbulent component.

Given this new random field (Resseguier et al., 2016) derived a new stochastic transport equation for a passive scalar that we will use as our new optical flow data term.

The rest of this article is organized as follows. In Sect. 2, the optical flow formulation under location uncertainty is derived in detail. Then evaluation on synthetic images is performed. To finish an experimental evaluation on a mixing layer is demonstrated and discussed.

Methodology description

Data term based on stochastic transport

Hereafter we will use the following notations :

-ℝ ! , for the image domain;

-𝑓: ℝ ! → ℝ ! , will be a stochastic process representing the image luminance;

-𝑿 ! 𝑥 : Ω×ℝ ! → Ω, will be a stochastic process representing the Lagrangian trajectory of a particle located at 𝑥 at the time 𝑡 = 0 (𝑿 ! 𝑥 = 𝑥);

The model introduced by (Mémin, 2014) is based on a decomposition of the Eulerian velocity field into two components :

𝒗 = 𝒘 + 𝝈𝑩, (2)
where -𝑩, is a random small-scale component, rapidly decorrelated in time, correlated in space, Gaussian, inhomogeneous and anisotropic, modeled by a Brownian process in time and a white noise process in space.

-
A particle transported in this stochastic Eulerian velocity field will have the following stochastic Lagrangian infinitesimal displacement :

d𝑿 ! = 𝒘 𝑿 ! , 𝑡 d𝑡 + 𝝈 𝑿 ! , 𝑡 d𝑩 ! , (3)
where 𝝈(𝑿 ! , 𝑡) is a diffusion operator that specifies the spatial correlations of the velocity uncertainty.

The conservation of a scalar quantity implies that the material derivative has to be equal to zero.

To express this material derivative (Resseguier et al., 2016) introduced a stochastic transport operator. Assuming that the luminance is proportional to the scalar density [START_REF] Liu | Fluid flow and optical flow[END_REF],

we can use this transport operator to express the luminance transport which is given by :

d ! 𝑓 + 𝒘 * d𝑡 + 𝝈d𝑩 ! ⋅ 𝛁𝑓 -𝛁 ⋅ 1 2 𝒂𝛁𝑓 d𝑡 = 0, (4)
where

𝒘 * = 𝒘 - ! ! 𝛁 ⋅ 𝒂 is a modified large-scale velocity component due to small-scale inhomogeneities.
The first term of the left-hand side of (4) is the temporal increment of 𝑓. The second term is the advection by both modified large-scale velocity and random small-scale velocity. The third term is a diffusion term where 𝒂 = 𝛔𝛔 ! is the diffusion tensor. Notice that the mathematical expectation of (4) is given by :

𝔼 d ! 𝑓 + 𝒘 * d𝑡 ⋅ 𝛁𝑓 -𝛁 ⋅ 1 2 𝒂𝛁𝑓 d𝑡 = 0, (5)
because the mathematical expectation of a Brownian process is null. So we obtain a new optical flow data term from (4). This single equation is not enough to estimate the motion field (a single equation for two unknown components). To solve this ill-posed problem we are looking for the velocity 𝒘 * which minimizes :

𝔼 d ! 𝑓 + 𝒘 * 𝑑𝑡 + 𝝈𝑑𝑩 ! ⋅ 𝜵𝑓 -𝜵 ⋅ 1 2 𝒂𝜵𝑓 𝑑𝑡 ! .
19th International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 16 -19, 2018

We use the local Lucas-Kanade [START_REF] Lucas | An Iterative Image Registration Technique with an Application to Stereo Vision[END_REF] resolution technique where the mathematical expectation can be seen as a Gaussian spatial filter 𝑔 ! ! , 𝜎 ! defining the size of the local interrogation window. Then it is easy to show that the minimization of

𝔼 d ! 𝑓 + 𝒘 * 𝑑𝑡 + 𝝈𝑑𝑩 ! ⋅ 𝜵𝑓 -𝜵 ⋅ 1 2 𝒂𝜵𝑓 𝑑𝑡 ! ,
with respect to 𝒘 * assumed locally constant is equivalent to the minimization of

𝑔 ! ! * d ! 𝑓 + 𝒘 * 𝑑𝑡 + 𝝈𝑑𝑩 ! ⋅ 𝜵𝑓 -𝜵 ⋅ 1 2 𝒂𝜵𝑓 𝑑𝑡 ! .
This latter leads to a linear system of two equations for two unknowns whose solution is given by :

𝑢 * 𝑣 * = 𝑔 ! ! * !" !" d𝑡 ! !" !" !" !" d𝑡 ! !" !" !" !" d𝑡 ! !" !" d𝑡 ! !! × 𝑔 ! ! * - !" !" d𝑡 d ! 𝑓 -∇ ⋅ ! ! 𝒂𝜵𝑓 d𝑡 - !" !" d𝑡 d ! 𝑓 -∇ ⋅ ! ! 𝒂𝜵𝑓 d𝑡 . (6)

Diffusion tensor estimation

The choice of the uncertainty model can be varied. We chose to begin with two simple models :

an isotropic diffusion and an anisotropic one without cross terms. Hence, in the isotropic case the diffusion tensor 𝒂(𝒙) takes the form below :

𝛼(𝒙) 0 0 𝛼(𝒙) ,
with the condition 𝛼(𝒙) ! ≥ 0 because 𝒂(𝒙) has to be symmetric positive definite. Under the assumption of a local homogeneity of 𝒂(𝒙) in a close neighborhood of the investigated pixel, the equation (4) becomes :

d ! 𝑓 + 𝒘 * d𝑡 + 𝝈d𝑩 ! ⋅ 𝛁𝑓 -𝛁 ⋅ 1 2 𝛼Δ𝑓 d𝑡 = 0. ( 7
)
19th International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 16 -19, 2018

Then we estimate 𝛼 locally by minimizing the expectation of the squared left-hand side of (7) with respect to 𝛼, which gives :

α = 2 𝑔 ! ! * Δ𝑓d𝑡 ! ×𝑔 ! ! * d ! 𝑓 + 𝒘 * d𝑡 ⋅ ∇𝑓 ∆𝑓d𝑡 . (8)
In the case of an anisotropic diffusion, tensor 𝒂(𝒙) takes the following form :

𝛼 ! (𝒙) 0 0 𝛼 ! (𝒙) ,
with the condition 𝛼 ! (𝒙)𝛼 ! (𝒙) ≥ 0 because 𝒂 𝒙 has to be symmetric positive definite.With a similar assumption as used for the isotropic case, the equation (4) becomes :

d ! 𝑓 + 𝒘 * d𝑡 + 𝝈d𝑩 ! ⋅ 𝛁𝑓 - 1 2 𝛼 ! 𝜕 ! 𝑓 𝜕𝑥 ! + 𝛼 ! 𝜕 ! 𝑓 𝜕𝑦 ! d𝑡 = 0. (9)
Then we estimate 𝛼 ! , 𝛼 ! locally by minimizing the expectation of the squared left-hand side of (9) with respect to 𝛼 ! , 𝛼 ! , which gives rise to a linear system (10) that we solve thanks to the Cramer rule :

𝛼 ! 𝛼 ! = 𝑔 ! ! * ! ! ! ! ! !" ! d𝑡 ! ! ! ! ! ! !" ! ! ! ! !" ! d𝑡 ! ! ! ! ! ! !" ! ! ! ! !" ! d𝑡 ! ! ! ! ! ! !" ! d𝑡 ! !! × 𝑔 ! ! * ! ! ! !" ! d𝑡 d ! 𝑓 + 𝒘 * d𝑡 ⋅ 𝛁𝑓 ! ! ! !" ! d𝑡 d ! 𝑓 + 𝒘 * d𝑡 ⋅ 𝛁𝑓 . (10)
A more complex form of the diffusion matrix can be considered. It is currently under investigation but we can already say that it gives rise to a complex and time-consuming constrained minimization problem.

Validation with synthetic images

Tutrb2D sequence

In order to evaluate our novel motion estimation method we use synthetic images database « Turb2D » coming from the FLUID project [START_REF] Heitz | Final report on the evaluation of the tasks of the workpackage 2[END_REF]. This database contains 256×256

Evaluation criteria

In order to assess statistically the performances of our motion estimation method we use as validation criteria the averaged magnitude error « Root Mean Square Error (RMSE) », the averaged angular error « Average Angular Error (AAE) » and the « Average Vorticity Error (AVE) ». Here are their mathematical definition :

𝑅𝑀𝑆𝐸 𝑡 = 1 𝑁 𝒖 𝒄 𝒙, 𝑡 -𝒖 𝒆 (𝒙, 𝑡) ! 𝒙∈! , (11)
𝐴𝐴𝐸 𝑡 = 1 𝑁 cos !! 𝒖 𝒄 (𝒙, 𝑡) ⋅ 𝒖 𝒆 (𝒙, 𝑡) 𝒖 𝒄 (𝒙, 𝑡) 𝒖 𝒆 (𝒙, 𝑡) 𝒙∈! , ( 12
) 𝐴𝑉𝐸 𝑡 = 1 𝑁 𝜔 ! 𝒙, 𝑡 -𝜔 ! (𝒙, 𝑡) ! 𝒙∈! , ( 13
)
where 𝑁 stands for the whole image pixel number, and where the subscripts 𝒄 and 𝒆 stand for the groundtruth field and the estimated field, respectively.

Results

(a)

(b) Concerning the particle image sequence, figure 3 Note that these synthetic cases are not necessarily relevant to simulate twin experiments with large-scale motion since the DNS provides the smallest fluid motion scales. However, the scalar and particles distributions could introduce some large scale and anisotropic effects. For the scalar case, LKMU gave better estimation by simulating the isotropic molecular diffusivity of the scalar (Sc = 0.7) via the diffusion tensor, which also explained no big differences between isotropic and anisotropic estimations after t=50 when the scalar was more homogeneously distributed. For the particle case, there was no real large-scale assumption (the size of the interrogation windows was small and the particle density was high), thus only slight differences were observed between LK and LKMU estimations and mainly explained by anisotropic effects.

(a) In this part we evaluate the performances of our method on real images of an isothermal mixing layer. This experiment was carried out in the wind tunnel HABV of Irstea-Rennes. The dimensions of the field of interest were 1m×1m. To enlighten this area we developed in collaboration with the french company EFFILUX an array of LEDs coupled with a linear Fresnel lens. The system emits a lightsheet of 1 cm thick at a wavelength around 465 nm. It was synchronized with a camera photron at 125 fps. A pulse duration of 1 ms was selected. The camera recorded 16 bits images of size 1024×1024 pixels. We tested two kinds of seeding: smoke and HFSB, see figure 6. The HFSB generator system was designed by LaVision. For the smoke we used a fog generator and designed two seeding ramps to contain the smoke in the lightsheet axis. The overall set-up is shown in figure 4. In figure 5 we display typical images we obtained for smoke and HFSB.

Results and discussion

Two sequences of velocity fields were computed with our new motion estimator : one for smoke seeding and one for HFSB seeding. From the mean velocity field exhibited in figure 7 we can derive some features of the mixing layer.

𝑈 𝑥, 𝑦 = 1 2 1 -𝑒𝑟𝑓 𝜋 𝛿 ! 𝑥 (𝑦 -𝑦 ! (𝑥)) 𝑈 ! -𝑈 ! + 𝑈 ! (14)
By fitting the law of (14) with the mean velocity field displayed in figure 7 and figure 8 we obtain characteristic parameters of the mixing layer. First we can retrieve the high velocity U ! and low velocity U ! which are very close to the wind tunnel set values reported in table 1. The two main first order parameters are the vorticity thickness δ ! (x) and the mixing layer axis y ! (x). r is the ratio

! ! ! ! and λ = !!! !!! .
The good agreement between measurements and the similitude law means that we are in the auto-similarity zone. By analysing more deeply table 1 we see that smoke seeding measurements are very close to the set values. No differences can be noticed between the istotropic and anisotropic versions of the motion estimation method. Concerning the HFSB seeding we can notice that the velocities are slightly over-estimated. This may come from the expulsion flow of the bubble generator array of nozzles which gives an initial velocity to the bubbles and which could perturbate the wind tunnel flow. For these reasons we plane to study the influence of bubble generator nozzles on the wind tunnel flow. Tab. 1 Mixing layer parameters.

To validate the estimation of the vorticity thickness we compared our results to hot-wire anemometry measurements from [START_REF] Heitz | Etude expérimental du sillage d'un barreau cylindrique se développant dans une couche de mélange plane turbulente[END_REF] and [START_REF] Sodjavi | Etude expérimentale de la turbulence dans une couche de mélange anisotherme[END_REF] x-x ν (mm) Smoke, anisotropic LKMU Smoke, isotropic LKMU Smoke, classic LK HFSB, anisotropic LKMU HWA [START_REF] Heitz | Etude expérimental du sillage d'un barreau cylindrique se développant dans une couche de mélange plane turbulente[END_REF] Pitot [START_REF] Heitz | Etude expérimental du sillage d'un barreau cylindrique se développant dans une couche de mélange plane turbulente[END_REF] 19th International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 16 -19, 2018 the Irstea wind tunnel and perturbated by the seeding ramps, and the other formerly investigated with HWA in one of the wind tunnel of university of Poitiers [START_REF] Heitz | Etude expérimental du sillage d'un barreau cylindrique se développant dans une couche de mélange plane turbulente[END_REF].

Nevertheless the global trend of the curves was consistent.

! ! ! ∆! ! , ! ! ! ∆! ! and ! ! ! ! ∆! ! (from top to bottom).
Comparisons between hot-wire anemometry measurements from Heitz (1999) for 𝜆 = 0.2 and smoke image velocimetry with anisotropic LKMU for 𝜆 = 0.18.

To test the accuracy limit of the method we investigated the kurtosis and skewness. As observed in figure 11 the proposed method is not able to properly retrieve the Kurtosis and the Skewness of the streamwise velocity component. We then performed a large-scale experiment on a mixing layer to evaluate the potential of our method with real images. Specific devices (illumination system, seeding system) were developed to carry out this large-scale experiment. First large-scale results are very encouraging, the main features of the mixing layer could be retrieved. It particularly shows that smoke image velocimetry measurements can be performed with an alternative motion estimation method.

However no significant effect of the stochastic data model could be observed with these real

particleFig 1

 1 Fig 1 Synthetic images of scalar (left) and particles (right) extracted from turb2D sequence of the FLUID project (Heitz et al., 2007).

Fig. 2

 2 Fig. 2 RMSE (a), AAE (b) and AVE (c) versus time for the three methods applied to scalar images. The classic LK method is the brightness consistency based Lucas-Kanade resolution (black line). LKMU is the stochastic transport equation based Lucas-Kanade resolution (blue and red lines).

Fig. 3

 3 Fig. 3 RMSE (a), AAE (b) and AVE (c) versus time for the three methods applied to particle images.

Fig. 4

 4 Fig. 4 Spectrum analysis of the turbulent velocity for DNS scalar image sequence at time t=100. Spectra of the error for the same data are shown in inset.

Fig. 5 Fig. 6

 56 Fig. 5 Experimental set-up.

Fig. 7

 7 Fig. 7 Mean velocity field computed from instantaneous fields obtained with a smoke seeding (note that mean field for HFSB is quite similar). Black dashed lines delimit the vorticity thickness of the mixing layer. The red one represents the mixing layer axis.

Fig. 8

 8 Fig. 8 Dimensionless stacked mean velocity profiles.

Fig. 9

 9 Fig. 9 Vorticity thickness evolution.

Fig. 10

 10 Fig. 10 Turbulence intensities

Fig. 11

 11 Fig. 11 Kurtosis, ! ! ! ! ! ! ! (a) and Skewness, ! ! ! ! ! ! ! ! (b).

  on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 16 -19, 20185. ConclusionIn this paper we presented a new formulation of optical flow for large-scale turbulent flows analysis. This formulation is based on a data term derived from the stochastic transport equation of a passive scalar. This model takes into account the interactions between unresolved scales and resolved ones. Our new motion estimator improved the estimation results for the scalar and particles synthetic images.

  experiments. We will continue to investigate the potential effects of the stochastic data term on a synthetic 3D flow and on a new measurements campaign in our wind tunnel. There will be a deeper validation by performing comparison with HWA measurements carried out in the same flow (ie in the same wind tunnel) and by investigating possible perturbations introduced by the intrusive seeding systems. This study will let us precisely characterize large-scale turbulent flow estimations extracted from both smoke and HFSB image sequences.

  19th International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY16 -19, 2018 unresolved scales is modeled by a trubulent diffusion computed a priori. If the estimation is improved, the main drawback of this approach comes from the necessity to know the power spectrum of the flow to compute the turbulent diffusion which is not necessarily known. In the continuation of this work[START_REF] Zillé | Observation models based on scale interactions for optical flow estimation[END_REF] proposed to model subgrid scales action during the multi-resolution process. Two options are exposed. The first one is a two-step procedure where the velocity field is computed once without subgrid modeling. Thanks to this motion field the subgrid scales are assessed and taken into account during the second step. The other option uses a subgrid scales model based on Smagorinsky coefficient. Despite the improvement they bring these strategies are strongly dissipative. Finally (Chen et al., 2015) improved the model of (Cassisa et al., 2011) using mix models introduced in LES by

  𝒘, is a smooth large-scale component, 19th International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 16 -19, 2018

  Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 16 -19, 2018

	19th International
		0.5
									anisotropic LKMU
									isotropic LKMU
		0.48							classic LK
		0.46
		0.44
	RMSE (pixel)	0.4 0.42
		0.38
		0.36
		0.34
		0.32
		20	30	40	50	60	70	80	90	100
						t
		14
									anisotropic LKMU
									isotropic LKMU
									classic LK
		13
		12
	AAE (°)	11
		10
		9
		8
		20	30	40	50	60	70	80	90	100
						t

  Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics・LISBON | PORTUGAL ・JULY 16 -19, 2018

	19th International
		0.1
									anisotropic LKMU
									isotropic LKMU
									classic LK
		0.095
		0.09
	AVE	0.085
		0.08
		0.075
		0.07
		20	30	40	50	60	70	80	90	100
						t

also shows the RMSE, AAE and AVE time evolution (frames 20 to 100). As we can see LKMU slightly improves the estimation, but not significantly (around 2%-3% of error reduction). However in this case anisotropic LKMU is always better than isotropic LKMU. In both cases (scalar and particles) the best improvement is observed for the vorticity error. The proposed model can better catch the vortex structures.

  where mixing layers with similar 𝜆 were studied. The characteristic parameters are reported in table 2. First-order characteristic parameters for different mixing layer with similar 𝜆.

	Parameters	Smoke	HFSB	(Sodjavi, 2013)	(Heitz, 1999)
	𝜆	0.18	0.17	0.2	0.2
	d𝛿 ! d𝑥	0.03	0,029	0.033	0.0336
	𝜎	59,08	61,12	53,71	52.75
	𝜎 !	10,63	10,39	10,74	10.55
	Tab. 2


Pharmacodynamics of regulatory T cells in mice and humans treated with low-dose IL-2

To the Editor: IL-2 is the major cytokine controlling homeostasis and function of regulatory T (Treg) cells. [START_REF] Malek | Interleukin-2 receptor signaling: at the iterface between tolerance and immunity[END_REF] At low dose, IL-2 specifically activates Treg cells and improves autoimmunity in mice and humans. [START_REF] Klatzmann | The promise of low-dose interleukin-2 therapy for autoimmune and inflammatory diseases[END_REF][START_REF] Yu | Selective IL-2 responsiveness of regulatory T cells through multiple intrinsic mechanisms supports the use of low-dose IL-2 therapy in type 1 diabetes[END_REF][START_REF] Saadoun | Regulatory T-cell responses to low-dose interleukin-2 in HCV-induced vasculitis[END_REF] As autoimmune diseases (AIDs) intrinsically reveal Treg-cell insufficiency, [START_REF] Klatzmann | The promise of low-dose interleukin-2 therapy for autoimmune and inflammatory diseases[END_REF] IL-2 has broad potential for their treatment. Because Treg cells control all immune responses and inflammation, low-dose IL-2 (ld-IL2) also has therapeutic potential for inflammatory diseases (IDs) [START_REF] Klatzmann | The promise of low-dose interleukin-2 therapy for autoimmune and inflammatory diseases[END_REF] and allergy [START_REF] Bonnet | Low-dose IL-2 induces regulatory T cell-mediated control of experimental food allergy[END_REF] and in the field of transplantation. [START_REF] Klatzmann | The promise of low-dose interleukin-2 therapy for autoimmune and inflammatory diseases[END_REF][START_REF] Koreth | Interleukin-2 and regulatory T cells in graft-versus-host disease[END_REF] However, given the pleiotropic effects of IL-2 on other immune cell types that also respond to IL-2 at higher doses, such as CD4 1 and CD8 1 effector T (Teff) cells, natural killer cells, and group 2 innate lymphoid cells, [START_REF] Klatzmann | The promise of low-dose interleukin-2 therapy for autoimmune and inflammatory diseases[END_REF] and given its short half-life, [START_REF] Konrad | Pharmacokinetics of recombinant interleukin 2 in humans[END_REF] finding a dose and schedule of administration that maintain a proper balance of Treg/Teff cells over time is key to the therapeutic use of ld-IL2. We aimed to rationally define IL-2 doses and schemes of administration that will lead to long-term improvement of the Teff-/Treg-cell balance.

We first investigated the dynamics of immune cells in C57BL/6 mice receiving an ''induction'' dosage of IL-2 aimed at a brisk stimulation of Treg cells, consisting of a 5-day course of daily intraperitoneal injection of 25,000, 50,000, or 100,000 IU (Fig 1,A) from day 0 to day 4. The percentage of blood Treg cells (defined as CD25 1 FOXP3 1 cells among CD4 1 T cells [START_REF] Koreth | Interleukin-2 and regulatory T cells in graft-versus-host disease[END_REF][START_REF] Rosenzwajg | Low-dose interleukin-2 fosters a dose-dependent regulatory T cell tuned milieu in T1D patients[END_REF] ) was assessed 24 hours after each administration and results are expressed as relative variation compared with the PBS-treated group. The maximum increase in Treg cells was observed at day 5 and was dependent on the cumulative dose received. We also analyzed the effects of IL-2 on blood Treg cells at 2 hours after IL-2 injections in the same mice. Blood Treg-cell proportions were decreased at 2 hours (see Fig E1 in this article's Online Repository at www.jacionline.org). This suggested a recirculation of Treg cells immediately after IL-2 injection and prompted us to study the dynamics of Treg cells in secondary lymphoid organs. In spleen, Treg-cell proportions were increased at 2 hours (see Fig E1). Collectively, these results indicate that there is an early recirculation or margination of blood Treg cells upon ld-IL2 administration. We then assessed what would be the effects of a ''maintenance'' dosage aimed at sustaining a long-term effect on Treg cells, administered as repeated single IL-2 administrations given after an induction course (see Fig E1). The 50,000 IU IL-2 induction course resulted in a 1.71 6 0.04-fold increase in blood %Treg cells (P < .0001) that was down to a 1.16 6 0.03-fold increase at day 8 (P 5 .02). A single reinjection of IL-2 at day 8 increased the Treg-cell proportion 24 hours thereafter 1.17 6 0.02-fold (P 5 .004). Further injections every 3 days sustained elevated levels of Treg cells over baseline values throughout the treatment period (from 1.16-to 1.33-fold increase). Spacing the injections every 7 or 15 days led to reduced Treg-cell increase during maintenance treatment.

We used these in vivo data to model the dose/scheme effects of IL-2 on Treg cells in mice during induction and maintenance phases. Data of Treg-cell expansion versus time were modeled and simulated using the NLME R package (R version 3.1.3, R development core team). A single-compartment model with first-order input and output was chosen. The evaluation of goodness-of-fit graphs showed that this model fitted the data well (Fig 1,B). We then used this model to simulate in silico reinjections every 3 days (Fig 1,C). The prediction fitted very well the actual in vivo data obtained using this treatment scheme (Fig 1 ,C).

We then aimed to adapt this model validated in mice to modeling human data. We used data from our previous study of ld-IL2 in type 1 diabetes (NCT01353833) that assessed the effects of a 5-day course of daily IL-2 injections. [START_REF] Rosenzwajg | Low-dose interleukin-2 fosters a dose-dependent regulatory T cell tuned milieu in T1D patients[END_REF][START_REF] Hartemann | Low-dose interleukin 2 in patients with type 1 diabetes: a phase 1/2 randomised, double-blind, placebo-controlled trial[END_REF] The expansion of blood Treg cells for each dose is presented in Fig 2,A (see Fig E2 ). For example, with 1 million international unit (MIU) IL-2 ''virtual'' injections every 3 days, Treg cells were expected to remain stable with a delta variation of %Treg cells ranging from 3 to 4 (ie, a 60%-80% increase); with 3 MIU every 3 days, %Treg cells were predicted to increase continuously; with reinjections every 15 days, Treg cells were expected to (1) return to their baseline approximately at day 20 with 0.33 MIU; (2) be stable with a delta variation of %Treg cells ranging from 0.5 to 1.5 (ie, 10%-30% increase) with 1 MIU; and (3) to be stable with a delta variation of % of Treg cells ranging from 2.5 to 3.5 (ie, 50%-70% increase) using 3 MIU.

On the basis of safety data from our previous clinical trials and our modeling, we chose the dose of 1 MIU/injection, with once-a-fortnight injections during the maintenance course to perform a novel clinical trial of ld-IL2 in patients with mild forms of various AIDs (rheumatoid arthritis, ankylosing spondylitis, systemic lupus erythematosus, psoriasis, Behcet disease, Wegener granulomatosis, Takayasu disease, Crohn disease, ulcerative colitis, autoimmune hepatitis, sclerosing cholangitis, and Gougerot-Sj€ ogren syndrome; Transreg clinical trial: NCT01988506). Treg-cell dynamics over time in the first 40 treated patients is presented in Fig 2 ,B. Note that Treg-cell evaluations were always made just before an IL-2 injection, thus 14 days away from the last injection. These values therefore underestimate the overall effect during the 14-day period. As predicted from above, %Treg cells was increased during the maintenance phase in treated patients (delta of 1.12% 6 0.25%, P < .0001 at day 31, and 1.325% 6 0.33%, P < .0001 at day 91), representing a minimum 20% increase for a prediction that was between 10% and 30%.

Thus, we have defined and validated a model for predicting the dynamic Treg-cell response to IL-2 that should help navigate its use in numerous therapeutic targets, from various AIDs to transplantation and allergy. Our pharmacodynamic model can predict the effects of any combination of any of the treatment parameters: dose, length of the induction phase, timing of maintenance injections. Also, our results show that the response to IL-2 is predictable across multiple AIDs and that we have validated a therapeutic scheme that provides long-term improvement in the Teff-/Treg-cell balance.

  , and was modeled using the model developed for mouse studies (see Fig E2 in this article's Online Repository at www.jacionline.org). Results are expressed as delta variation of Treg-cell frequency in CD4 1 T cells; that is, % of Treg cells in CD4 1 of the day minus % of Treg cells in CD4 1 at day 0 and vertical dotted lines represent IL-2 injection. We then predicted in silico the potential effects of different doses and the periodicities of IL-2 reinjections given after an induction course

FIG 1 .

 1 FIG 1. Dynamics and modeling of Treg cells in mice treated with IL2. A, Dynamics of Treg cells in blood during a 5-day course of IL-2 (pool of 3 different experiments; n 5 10). Observed and in silico modeling of Treg-cell dynamics in blood (B) after a 5-day course of daily injections of 50,000 IU of IL-2, and (C) followed by reinjections of 50,000 IU of IL-2 every 3 days.

FIG 2 .

 2 FIG 2. Dynamics and modeling of CD4 1 Treg cells in humans after ld-IL2. A, Observed Treg-cell dynamics in humans after a 5-day course of ld-IL2 (doses ranging from 0.33 MIU/d to 3 MIU/d) (n 5 6 patients/treatment). B, Comparison of in silico modeling and real Treg-cell dynamics (n 5 40) after reinjections of 1 MIU of IL-2 every 15 days after a 5-day course of IL-2.

FIG E1 .

 E1 FIG E1. Dynamics and modeling of Treg cells in mice treated with IL2. Results are expressed as relative variation compared with the PBS-treated group. These data are from a pool of 3 different experiments (n 5 10) using C57BL/6 mice. Treg-cell dynamics in blood at 2 and 24 hours postinjection (A) and Treg-cell dynamics in blood and spleen at 2 hours postinjection (B), during a 5-day course of 50,000 IU/d of IL-2. C, Dynamics of Treg cells in blood during a 5-day course of 50,000 IU of IL-2 and then reinjection every 3, 7, or 15 days with 50,000 IU of IL-2. Comparison between IL-2 and PBS groups was done using the Mann-Whitney test (*P < _ .05; **P < _ .01; ***P < _ .0001).

FIG E2 .

 E2 FIG E2. Dynamics and modeling of CD4 1 Treg cells in humans after ld-IL2. A and B, Vertical dotted lines represent IL-2 injection. Blood sampling was performed just before IL-2 injection. Results are expressed as delta variation of Treg-cell frequency in CD4 1 T cells, that is, % of Treg cells in CD4 1 of the day minus % of Treg cells in CD4 1 at day 0. Fig E2, A, Modeling of the results obtained in humans after a 5-day course of ld-IL2 (doses ranging from 0.33 MIU/d to 3 MIU/d). Fig E2, B, In silico modeling of Treg-cell dynamics after injections of ld-IL2 every 3 or 15 days after a 5-day course of ld-IL2, using 0.33, 1, or 3 MIU/injection.

We are grateful to Alice Barrateau for excellent technical assistance. We thank the animal care team members Christelle Enond, Bocar Kane, Olivier Bregerie, and Serban Morosan from the Centre d'Exploration Fonctionnelle, UPMC Paris 06. We thank the entire DF-IL2 and Transreg study groups for the performance of the clinical trials that contributed to this study (groups are listed in this article's Online Repository at www.jacionline.org).

Guillaume Churlaud, PharmD, PhD a,b * Chadi Abbara, PharmD, PhD c * Pierre-Axel Vinot, MSc a,b Gwladys Fourcade, PhD b Paul-Gydeon Ritvo, PhD b Roberta Lorenzon, MD a,b Michelle Rosenzwajg, MD, PhD a,b Bertrand Diquet, PharmD, PhD c David Klatzmann, MD, PhD a,b From a AP-HP, Hôpital Piti e-Salpêtri ere, Biotherapy (CIC-BTi) and Inflammation-Immunopathology-Biotherapy Department (i2B), Paris, b Sorbonne Universit e, UPMC Univ Paris 06, Immunology-Immunopathology-Immunotherapy (I3), Paris, and c Service de Pharmacologie, Toxicologie et centre de pharmacovigilance, CHU et facult e de sant e, Universit e d'Angers, Angers, France. E-mail: david. klatzmann@sorbonne-universite.fr. *These authors contributed equally to this work. The project was funded by the LabEx Transimmunom (grant no. ANR-11-IDEX-0004-02) and European Research Council Advanced Grant TRiPoD (322856) to D.K. G.C. was supported by the Acad emie Nationale de M edecine. Disclosure of potential conflict of interest: M. Rosenzwajg and D. Klatzmann are inventors in a patent application related to the therapeutic use of ld-IL2 for treating autoimmune-related or inflammatory disorders, which belongs to their academic institutions and has been licensed to ILTOO pharma. G. Churlaud, M. Rosenzwajg, and D. Klatzmann hold shares in ILTOO pharma. The rest of the authors declare that they have no relevant conflicts of interest.


Romain Schuster
email: romain.schuster@irstea.fr

Dominique Heitz

Philippe Georgeault

Etienne Mémin

ON-SITE AIRFLOW MEASUREMENT OF A LABORATORY FUME HOOD USING CUSTOMIZED LARGE-SCALE IMAGE-BASED VELOCIMETRY

Keywords: fume hood, large-scale PIV, optical flow

come

On-site airflow measurement of a laboratory fume hood using customized large-scale image-based velocimetry

INTRODUCTION

Laboratory fume hoods are protection devices generally used in places where harmful substances are handled. They are aimed at limiting the exposure of humans to contaminants by capturing and removing hazardous fumes or dust generated inside the enclosure thanks to a suction flow. To insure their good performance through time, fume hoods are subject to periodical controls. These controls usually consist of a three-part evaluation process composed by a face flow rate evaluation, a smoke visualization and a tracer gas test. This protocol does not give any information on the presence or the location of potential leakages. To detect potential pollutants leakages it is necessary to understand the aerodynamics features of the fume hood. It can be done by flow visualization as shown by some studies [START_REF] Tseng | Correlation Between Airflow Patterns and Performance of a Laboratory Fume Hood[END_REF][START_REF] Chen | Installation of a Flow Control Device in an Inclined Air-Curtain Fume Hood to Control Wake-Induced Exposure[END_REF]. However in these studies flow visualization is only qualitative. In fact it is possible to extract quantitative flow motion estimation from flow visualization [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]. This is what we propose in this paper.

In this study we propose to perform image-based air velocities measurements on a laboratory fume hood in operating conditions. Traditionally particle image velocimetry (PIV) is carried out to estimate fluid motion. This method is widely used to study turbulent flows in fluid mechanics community. PIV principles and advancements have been carefully explained and detailed in the book Particle Image Velocimetry [START_REF] Mémin | Particle Image Velocimetry[END_REF]. Image-based velocimetry techniques like PIV are nonintrusive and provide instantaneous velocity fields. PIV was originally designed for laboratories and well-controlled flows, and its transposition to an on-site measurement tool is not immediate. Some ventilation or air conditioning systems inside cars, planes or rooms were investigated by PIV [START_REF] Cao | Particle image velocimetry measurement of indoor airflow field: A review of the technologies and applications[END_REF] but experiments were mostly carried out in research laboratories on specifically designed models and rarely in real conditions. Thus performing PIV in real conditions needs some adaptations.

As fields of investigation are getting bigger, classical particles like oil droplets cannot be used as tracer, they are too small. Two options are then possible: using bigger particles like Helium Filled Soap Bubbles (HFSB) or using smoke. The advantage of HFSB is that they are big enough (diameter around 300 µm) to strongly scatter the light and being visible from far. Secondly thanks to the helium, bubbles have the neutral buoyancy needed to follow the flow with high fidelity [START_REF] Scarano | On the use of helium-filled soap bubbles for large-scale tomographic PIV in wind tunnel experiments[END_REF]. The interest of employing a smoke generator is it simplicity and availability. Smoke is compound of a high density of oil droplet haze scattering the light and providing smooth images, difficult to analyse with correlation techniques but from which fluid motion can be estimated via optical flow approaches [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF].

The following study is aimed at demonstrating the feasibility of an on-site PIV-like measurement of a fume hood airflow in operating conditions. We will compare the two kinds of tracers mentioned above: HFSB and smoke, in order to determine what tracer is the best suited for large-scale on-site flow characterization. The rest of this paper is organized as follows. In Section 2, the image-based method is described in detail. Then results and performances of the method are presented and commented in Section 3. Finally, conclusion is discussed in Section 4.

METHODS

The method consists of bringing to light the invisible airflow generated by a laboratory fume hood. To make it possible, the flow is seeded with a tracer (particles or smoke). In the same time a light sheet is emitted by an illumination system in order to enlighten the tracers in the zone of interest. The light scattered by the tracer is recorded by a camera at a given frame rate. By computing tracer's displacements between successive frames thanks to a motion estimation method it is possible to obtain an estimation of the velocity field in the illuminated slice. Naturally as we only illuminate a plane and use only one camera we can only retrieve two components of the velocity in two directions (2D2C). However 2D2C PIV is relevant for strongly bidimensional flows, as it is the case for fume hood.

Experimental setup

Usually, to perform a PIV measurement a Nd:YAG pulsed laser at 532 nm combined with a cylindrical lens is employed in order to generate a light sheet. The light delivered by the laser is powerful enough for large-scale applications but the whole light system is not easy to move everywhere and is very expensive. Therefore we have chosen to develop our own illumination system more transportable, easy handling and less expensive. Our device consists of an array of LEDs coupled with a linear Fresnel lens that produces a collimated light sheet. LEDs were chosen to emit at 432 nm because of their better efficiency at this wavelength. The LEDs array generates a 1 m high and about 10 mm thick light sheet. The illumination system can be run both in trigger mode and in continuous mode. The advantage of trigger mode is the possibility of scaling up the light power for a short duration. In our configuration the LEDs array was placed vertically facing the fume hood front opening.

The HFSB generator that we used produces small bubbles down to 300 µm diameter size with a high seeding density (up to 2 million bubbles per second).

Image recordings were acquired with a CCD camera equipped with a 35 mm focal lens and storing 1024 x 1024 pixels images at 60, 125 and 250 fps. The camera was positioned perpendicular to the illuminated plane at a distance of 1,5 m from it.

The laboratory fume hood used for this study was located at the Irstea center of Rennes. We were constrained by the geometry of the room. Experiments were conducted at sash window heights of 500 mm. For this study the fume hood was unoccupied. We think that performing our method with a manikin should be relevant to see the influence of human body on the flow, but this will be carried out in further studies.

Motion estimation from image sequences method

As mentioned above, we seeded the region of interest with a tracer in order to make the fluid flow visible. With a well-chosen tracer that can closely follow the fluid motion, we can assume that the tracer motion equals the fluid motion. In this section we expose our strategy to estimate the fluid motion from observations of tracer displacements during an image sequence. This issue is frequently addressed, in PIV community, by a correlation based processing between two images. Here, to compute the motion from two successive snapshots we used an optical flow method. Optical flow is a method that comes from the computer vision community [START_REF] Horn | Fluid flow dynamics under location uncertainty[END_REF]. Precisely what is called optical flow is the apparent motion between two successive images. Optical flow is based on an observation model that links the spatio-temporal tracer luminance variations to the apparent fluid motion field. One of the benefits of using optical flow is that it provides a dense motion field (1 vector/pixel), much better than correlation resolution limited to 0.13 vector/pixel in the best case. Moreover, compared to correlation approaches optical flow provides a framework to use physicalbased models [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF], which improves the estimation. Several observation models have been investigated to cope with complex fluid motion. In our case we developed a largescale observation model that takes into account the effects of unresolved scales and errors, on the scale of the measured velocity. In this model introduced by Mémin (2014), non-observed scales are expressed thanks to stochastic formalism. Based on the passive scalar transport equation where subpixel scales are assumed random, this model is well suited for large-scale observations encountered in real world. This new stochastic multiscale data term has already shown good performances for the motion estimations of turbulent flows [START_REF] Cai | Motion estimation under location uncertainty for turbulent fluid flows[END_REF]. Concerning the optimization strategy we chose a local resolution scheme exposed by [START_REF] Lucas | An Iterative Image Registration Technique with an Application to Stereo Vision[END_REF], which is well-suited for parallelization (because each pixel is processed independently), this is interesting for real-time in situ applications.

RESULTS AND DISCUSSION

Mean flow analysis

Figure 1 shows two mean velocity fields obtained from HFSB seeding and smoke seeding. As we can see, the two presented fields are not exactly the same although they are supposed to represent the same flow. In the central region of the sash opening the two mean fields are quite similar, but at the bottom and at the top the two vector fields differ. It comes from the fact that the smoke was not well distributed everywhere during the image acquisition process. As a matter of fact the smoke generator was driven manually and it was not easy to constantly cover the entire zone of interest. Moreover images need to have a lot of luminance gradients, so the smoke generation should be pulsed in order to create smoke structures and not a diffuse haze. For this reasons some regions like bottom or top of the fume hood were not always well seeded leading to underestimated vectors. In opposition the HFSB generation was run automatically so it was easier to insure a dense homogeneous and continuous seeding in the whole field of observation.

Figure 2 shows vertical profiles of the mean velocity magnitude taken between the fume hood workbench (height = 0 mm) and the glass inferior edge (height = 500 mm). The goal is to compare our method to a typical measurement that is commonly done during periodical controls namely spot hotwire anemometry. As we can see, the HFSB seeding gave results close enough to the HWA (except at extremities where the airspeed was too weak to be measured by the hot-wire anemometer) but with a clearly better spatial resolution. As expected, the smoke gave the global trend but slightly underestimated the velocity for the reasons mentioned above.

Turbulent flow features analysis

Figure 3a is a map of the flow turbulence intensity defined as <u' 2 >/U out 2 where <u' 2 > stands for the quadratic time average of instantaneous fluctuations defined as u'(t)=u(t)-U mean . We can clearly distinguish areas where the turbulence is stronger. The obvious one is located just behind the edge of the bench where the flow is highly turbulent due to the corner wake and bench boundary layer. Another area, but with lower intensity, corresponds to the wake of the bottom edge of the sash. Finally a third big area at the top inside the fume hood behind the glass seems to present a high level of turbulence. As we will see next, turbulence activity is a sign of flow unsteadiness, which can give rise to eddies and in the worst case cause leakages. The main interest of an image-based method compared to HWA is the possibility of visualizing the space and time flow evolution and thus following dynamic flow features, which is difficult with HWA. Figure 3b and3c shows instantaneous streamlines and vorticity fields taken from the HFSB sequence. As we can see we can retrieve flow features that were already observed through visualization by [START_REF] Tseng | Correlation Between Airflow Patterns and Performance of a Laboratory Fume Hood[END_REF]

CONCLUSIONS

In this study we have demonstrated the possibility of performing on-site large-scale PIV measurements on an operating fume hood. We tested two seeding systems: HFSB and smoke. Results were better with bubbles than with smoke. However the smoke generator can be improved and when

Recirculation bubble

Vortex shedding in bottom sash edge wake

Coherent flow structures in bench corner wake and bench boundary layer the smoke distribution is well structured the vector fields can be quite good [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]. Thanks to our method we could first retrieve the mean velocity field over an extended area. There are clear advantages over HWA: we got a large distribution of two-components velocities on a whole plane when for HWA we only got one-component sparse spot measurements. We also could quantify the turbulence intensity and some aerodynamic features of the flow like recirculation bubbles or coherent flow structures, which is very important for leakages detection.

To complete this study, the next step could be to perform our method on several planes (vertical and horizontal) in order to have a 3D map of the flow features. We also would like to test the method with a manikin standing in front the opening to see the effect of the presence of a human body on the flow.

In conclusion, this study let us believe that the exposed method could be widely applied for periodical control of laboratory fume hoods or other ventilation system inside buildings. To our minds it delivers a much more complete flow characterization than actual periodical controls based on hot-wire anemometry.

Figure 1 .

 1 Figure 1. Mean velocity fields obtained for an opening height of 500 mm with a continuous light. Results for two kinds of seeding are shown: HFSB (left) and smoke (right).

Figure 2 .

 2 Figure 2. Vertical profiles of the U x mean velocity component taken at the glass door position on the x-axis. Results from HFSB seeding (red line and orange dashed-line) and Smoke seeding (blue line)are compared to sparse hot-wire anemometry measurements (green errorbars).

Figure 3 .

 3 Figure3ais a map of the flow turbulence intensity defined as <u' 2 >/U out 2 where <u' 2 > stands for the quadratic time average of instantaneous fluctuations defined as u'(t)=u(t)-U mean . We can clearly distinguish areas where the turbulence is stronger. The obvious one is located just behind the edge of the bench where the flow is highly turbulent due to the corner wake and bench boundary layer. Another area, but with lower intensity, corresponds to the wake of the bottom edge of the sash. Finally a third big area at the top inside the fume hood behind the glass seems to present a high level of turbulence. As we will see next, turbulence activity is a sign of flow unsteadiness, which can give rise to eddies and in the worst case cause leakages. The main interest of an image-based method compared to HWA is the possibility of visualizing the space and time flow evolution and thus following dynamic flow features, which is difficult with HWA. Figure3b and 3cshows instantaneous streamlines and vorticity fields taken from the HFSB sequence. As we can see we can retrieve flow features that were already observed through visualization by[START_REF] Tseng | Correlation Between Airflow Patterns and Performance of a Laboratory Fume Hood[END_REF] and by[START_REF] Chen | Installation of a Flow Control Device in an Inclined Air-Curtain Fume Hood to Control Wake-Induced Exposure[END_REF], like recirculation area behind the bench corner or vortex shedding in the bottom edge of sash wake. The intensity and persistence in time of such flow structures give an indication about the location of potential leakages.


T Chattopadhyay

- D Fraix-Burnet

- S Mondal

Typeset using L A T E X default style in AASTeX62 Unsupervised Classification of Galaxies. I. ICA feature selection

Keywords: galaxies, general -galaxies, evolution -methods, statistical

Subjective classification of galaxies can mislead us in the quest of the origin regarding formation and evolution of galaxies since this is necessarily limited to a few features. The human mind is not able to apprehend the complex correlations in a manyfold parameter space, and multivariate analyses are the best tools to understand the differences among various kinds of objects. In this series of papers, an objective classification of 362,923 galaxies from the Value Added Galaxy Catalogue (VAGC) is carried out with the help of two methods of multivariate analysis. First, Independent Component Analysis (ICA) is used to determine a set of derived independent components that are linear combinations of 47 observed features (viz. ionized lines, Lick indices, photometric and morphological properties, star formation rates etc.) of the galaxies. Subsequently, a K-means cluster analysis is applied on the nine independent components to obtain ten distinct and homogeneous groups. In this first paper, we describe the methods and the main results. It appears that the nine Independent Components represent a complete physical description of galaxies (velocity dispersion, ionisation, metallicity, surface brightness and structure). We find that our ten groups can be essentially placed into traditional and empirical classes (from colour-magnitude and emission-line diagnostic diagrams, early-vs late-types) despite the classical corresponding features (colour, line ratios and morphology) being not significantly correlated with the nine Independent Components. More detailed physical interpretation of the groups will be performed in subsequent papers.

INTRODUCTION

Investigating the formation and evolution of galaxies is becoming a complicated process with the increased availability of huge databases as a result of instrumental improvements. A good understanding of the underlying physical processes requires synthetic numerical simulations of the formation and evolution of galaxies in an evolving Universe. According to various studies (e.g. see reviews in [START_REF] Silk | [END_REF]Genel et al. 2014), classical formation of galaxies have been proposed to follow five major models: (i) the monolithic collapse model, (ii) the major merger model, (iii) the multiphase dissipational collapse model, (iv) accretion and (v) in situ hierarchical merging. But no model uniquely explains the formation of all galaxies.

The historical and still most common approach to the classification of galaxies, is based on physical criteria, such as apparent features (i.e. morphology, emission line properties etc.) or more or less understood processes (starbursts, Active Galactic Nuclei (AGN) etc.). Hubble's subjective classification based on galaxy morphology ignores many significant observable features such as kinematics, chemical composition etc. Physically based classifications also rely on only very few attributes among the numerous ones available: colour, size, environment, star formation rate, emission lines. All these classifications are certainly useful to study a specific property, but they cannot embrace the whole diversity and complexity of processes that shaped the galaxies.

With the advent of multi-wavelength and multivariate databases, the goal of many studies has been to find the properties that best characterize the established classes. One such example is given by the Baldwin, Phillips & Tarlevich (BPT) diagrams (Baldwin et al. 1981;Veilleux & Osterbrock 1987) built from a few ratios of emission lines such as [OIII]/Hα, [OIII]/Hβ, [NII]/Hβ, [SII]/Hα and [OI]/Hα. Delimitation lines have been defined empirically to separate different kinds of ionizing sources (AGNs, Low-Ionization Nuclear Emission-Line Regions (LINERs), Seyfert galaxies, star forming regions). But a multivariate analysis using [OIII]/Hα, [NII]/Hβ and EW(Hα) have shown that the empirical classes are not all statistically supported (de Souza et al. 2017). One probable explanation is that several kinds of ionizing sources can be present in a single galaxy, so that a classification built from a two-feature diagram cannot reasonably translate the real complexity of even a single galaxy.

These classification approaches can be called object-based, to be opposed to data-driven science in which some characteristics are revealed by the accumulation of data and attributes. The so-called bimodality of galaxies (Strateva et al. 2001;Eales et al. 2018) is an example of the latter: the observations of millions of galaxies have revealed that for several properties, galaxies show two distribution peaks that do not easily match pre-established physically motivated classes (e.g. Eales et al. 2018;Evans et al. 2018). By chance, this bimodal structure can be seen with the eye on 2D diagrams.

Extragalactic studies have now entered the statistical era challenging our approaches to understand the physics of galaxies and their evolution. In most cases, there is no simple mathematical relationship between a given physical process and the relevant observable attributes. In addition, these attributes are generally multiple so that classifications based on only one or two features are not adequate to represent the physics of galaxies. Despite this caveat, many studies have embarked into complex statistical techniques (Bayesian tools, Deep Learning) to tackle the avalanche of data that astronomy is about to generate with the goal to obtain an automatic classification (e.g. Brescia et al. 2015). However, these approaches, most often multivariate, use supervised learning techniques, i.e. the algorithm is trained on labelled data. The problem here is that the labels come from pre-established classifications based on very few features, in other words the algorithm is trained to find the results of human subjectivity (e.g. Cavuoti et al. 2014).

In this context, one is tempted to apply statistical unsupervised classification, such as a multivariate partitioning analysis to find homogeneous groups, not focusing on only one aspect of the physics of galaxies, but by exploring the cluster structure of the dataset (Fraix-Burnet et al. 2015). The goal is to use data-driven science to study the formation and evolutionary history of galaxies. One basic tool is Principal Component Analysis (PCA) and it has been used by many authors (e.g. Whitmore 1984;Watanabe et al. 1985;Cabanac et al. 2002;Chattopadhyay & Chattopadhyay 2006;Peth et al. 2015). However this is not a classification tool. It may be useful to represent pre-established classes using some of the Principal Components, but its main interest is a reduction of the dimensionality. When the number of attributes is large, distances (e.g. euclidean) become less discriminant so that cluster structures cannot be detected (this is one aspect of the curse of dimensionality). Feature selection, for instance through an objective dimensionality reduction, is then compulsory. But the use of Principal Components to perform a clustering (unsupervised classification) is not recommended since the components with the largest eigenvalues define the axes of maximum variance in the dataset, and these axes are generally not the most discriminative ones to reveal the cluster structure (Chang 1983).

Regarding unsupervised classification, some attempts have been made by K-means cluster analysis (Ellis et al. 2005;Chattopadhyay & Chattopadhyay 2007;Chattopadhyay et al. 2007Chattopadhyay et al. , 2008;;Chattopadhyay et al. 2009a,b;Sánchez Almeida et al. 2010;Fraix-Burnet et al. 2010[START_REF] Silk | [END_REF]De et al. 2016). Though sophisticated statistical techniques are being developed steadily, unsupervised learning approaches are not widely used across the astronomical community (Fraix-Burnet et al. 2015;de Souza et al. 2017) probably because of the difficulty to relate the multivariate results with physical models. This issue can however be overcome, for instance by comparing statistically the multivariate results with the outcomes of numerical simulations of galaxy and cosmic evolution (e.g. Fraix-Burnet et al. 2012;Genel et al. 2014).

A partitioning of objects into robust groups is possible when the features are independent. Observationally the information is usually summarized into broad-band fluxes (magnitudes), slopes (colours), medium-band and line fluxes (Lick indexes) etc, because they can be easily measured and can generally be yielded by models and numerical simulations. However, their relationship to intrinsic physical processes and their mutual influences are most often quite complex. As a consequence, these direct observables cannot suffice to compare two galaxies and the multivariate aspect of the physics of galaxies must be included as well. This is exactly the purpose of the Independent Component Analysis (ICA) that we will use in this paper.

In this study, we have taken a large dataset from Sloan Digital Sky Survey (SDSS) data archive including various observables regarding morphology, chemical composition and kinematics and used multivariate statistical techniques to explore and explain the underlying diversities. This work presents several novelties for unsupervised classification in astrophysics:

• a large dataset of galaxies

• a large number of features

• the Independent Component Analysis (ICA)

• the adequation of the method (ICA) to the fact that our data are non-Gaussian Some other studies have used the SDSS to get a large dataset (e.g. Cavuoti et al. 2014) or use many features (up to 4520!, D'Isanto et al. 2018) but they perform supervised classification which, in the case of continuous values, is equivalent to a regression problem. The ICA technique has already been used in astrophysics, mainly for source separation (Pires et al. 2006;Pike et al. 2017;Martins-Filho et al. 2018;[START_REF] Sheldon | American Astronomical Society Meeting Abstracts #[END_REF] and dimensionality reduction (Richardson et al. 2016;Sarro et al. 2018) but rarely for unsupervised classification (Mu 2007;Das et al. 2015).

This paper is organized as follows. A brief description of the dataset is given in Section 2. The methods are described in Section 3. The results and discussion are included in Sections 4 and 5, respectively. Finally, conclusions are traced in Section 6.

DATASET

The NYU Value-Added Galaxy Catalogue (VAGC Blanton et al. 2005;Padmanabhan et al. 2008;Abazajian et al. 2009) is a cross-matched collection of galaxy catalogues maintained for the study of galaxy formation and evolution1 . It is based on the Sloan Digital Sky Survey Data Release 7 (SDSS-DR72 ).

In the raw table, 2,506,754 objects are available. We have selected only galaxies, by disregarding QSOs and stars using the SDSS flag and subsequently by removing some obviously wrong classification (such as entries with aberrant redshifts or magnitudes). Non-galaxy objects probably remains, but they should show up in the multivariate classification. We ended up with 865,333 entries. We have then restricted the sample to z < 0.2 to avoid too much shift between the wavelength ranges in which the magnitudes are obtained. Finally, only data with a good signal to noise ratio (median S/N per pixel of the whole spectrum > 10) were kept. This leaves us with 362,923 galaxies.

We had to limit the number of attributes to keep the computation tractable. For this, we eliminated redundant properties, i.e. features that bear the same information (such as the Sersic profile in different bands), and selected only a few photometric bands and colours. We were careful to keep most of the physical information so that it does not impact much our analysis based on dimensionality reduction through the ICA. The final set used in our analysis consists of 49 attributes which cover photometry, spectroscopy, morphology, chemical composition and kinematics. Star formation rates and specific star formation rates are also included but not used in the ICA analysis itself since they are not observable features. All these attributes are described in Table A1 and details are given on the source website3 . Please note that the equivalent width is negative if the line is in emission, EW(Hα) and EW(NII 6584 ) are not corrected for possible blending of these two lines, and lines which are absent in spectra are given the value 0.

STATISTICAL ANALYSES

Shapiro-Wilk test

The non-Gaussian nature of the dataset has been explored by the statistical Shapiro-Wilk test (Shapiro & Wilk 1965) in which the test statistics is defined by W

= n i=1 a i x 2 i / n i=1 (x i -x) 2
, where n is the number of observations,

x i 's are ordered sample values and a i 's are the constants generated from the order statistics of a sample from normal distribution. In this work a multivariate extension (Alva & Estrada 2009) has been used. The p value of the test is 2.17×10 -13 , which is sufficiently small to confidently reject the null hypothesis. Therefore, the dataset is found to be non-Gaussian in nature.

Independent Component Analysis

We have already mentioned that Principal Component Analysis (PCA) has been applied by many authors (Brosche 1973;Whitmore 1984;Murtagh & Heck 1987, etc) but it is not appropriate for clustering and classification (Chang 1983). Furthermore it is applicable to Gaussian data which is not the present case. On the other hand, Independent Component Analysis (ICA) is also a dimension reduction technique, i.e. it reduces the number of observed features p to a number m (m << p) of new variables, the components, but assumes non-Gaussian data (Pfister et al. 2018;Hyvärinen 1998Hyvärinen , 1999a,b),b). The second important difference with PCA is that in addition of being uncorrelated, the components are also mutually independent (for details regarding comparison between PCA and ICA see Section 3 of Chattopadhyay et al. 2013, and references therein).

Mathematically speaking, let X 1 , X 2 , X 3 , ..., X p be p random vectors (here p features, p = 47) and n (here 362,923) be the number of observations of each X i , (i = 1, 2, 3, ..., p).

Let X = AS, where S = [S 1 , S 2 , S 3 , ..., S p ] ′ is a random vector of hidden components S i , (i = 1, 2, 3, ..., p) such that S i 's are mutually independent and A is a non singular matrix. Then the objective of ICA is to find S by inverting A, i.e., S = A -1 X = WX. Since ICA is able to separate independent components (sources) present in a signal, W is called the unmixing matrix (for more details see Comon 1994;Chattopadhyay et al. 2013, and references therein). Independence is obtained by maximizing the non-Gaussianity using negentropy.

Presently there is no good method available for the determination of the optimum number of ICs. In this work, the optimum number of ICs have been chosen by the optimum number of Principal Components (PCs) (Albazzaz & Wang 2004;Chattopadhyay et al. 2013;Eloyan & Ghosh 2013), to find m (m << p) (Chattopadhyay & Chattopadhyay 2007;Babu et al. 2009;Fraix-Burnet et al. 2010;Chattopadhyay et al. 2010Chattopadhyay et al. , 2013)). We have first performed PCA to find the significant number of ICs. In PCA the maximum variation with significantly high eigenvalue (viz. λ ∼ 1) was found to be almost 90% for nine PCs. Hence, we have chosen nine ICs (i.e. here m = 9) for cluster analysis (CA).

Here it is worth mentioning that the choice of the ICs is a difficult question (see Kairov et al. 2017) but we have checked that the set of most correlated original features (Section 4.1 and Table 1) is stable through multiple runs of the algorithm. Note that for our dataset PCA takes 132 seconds whereas ICA takes less than 1 minute (Intel(R) Core TM i3-2100 CPU @ 3.10GHz, RAM: 2.00 GB, Windows 10, 64 bit operating system) and we have used the library functions publicly available for ICA and PCA in R.

K-means cluster analysis

K-means cluster analysis (CA) is a multivariate technique for finding coherent groups in a dataset giving information of the underlying structure. In this method:

• Each object must belong to a single cluster.

• Each cluster must contain at least one object.

• All the objects are distributed among K clusters by satisfying the first two properties. (2010,2012,2013); Das et al. (2015).

The number K of groups is an input to the algorithm. The optimum value of K is found as follows. For a particular choice of K initially we compute a distortion measure d K , given by,

d K = (1/p)min x E[(x K -c K ) ′ (x K -c K )]
which is the distance of x K vector (data point) from the centroid c K of the corresponding group. Then we compute the value

of jump J K = (d -p/2 K -d -p/2 K-1
) where p is the total number of variables (here total number of independent components under consideration i.e. p = 9). The above d K and J K values are computed for several choices of K starting from K = 1, 2, 3, etc. The maximum value of the curve obtained by plotting J K versus K gives the optimum value of K (Sugar & James 2003).

In this study, we have performed a K-means CA with respect to the ICs and have found the optimum number of groups to be K = 10. We name the groups K1 to K10.

RESULTS

Properties of the ICs

Performing Independent Component Analysis for the VAGC dataset, we have taken nine significant ICs. Then we have computed the correlation coefficients of each component with the attributes to understand their physical meaning. We list only the most influencial attributes in Table 1.

From Table 1, it is clear that the nine ICs represent five kinds of properties: 1) velocity dispersion (ICA 2, ICA 3, ICA 7), 2) ionisation (ICA 4, ICA 9), 3) metallicity (ICA 1, ICA 5), 4) surface brightness (ICA 8) and 5) structural properties (ICA 6). This shows that the ICA successfully retrieve the main physical ingredients of galaxies, while also taking into account their interactions through the linear combinations of the observed features.

This result is indeed very significant. Usually, a given physical process (star formation, nuclear activity, kinematics, radiation characteristics, etc) is investigated through an obvious observable feature, leading by the way to the traditional classifications mentioned in Section 1. Here, each of the ICs is dominated by some features specific to a particular physical property of galaxies, without any a priori selection. In addition, each IC is not limited to the dominant features, it also includes the other ones with some lower weights thus taking into account the complex interplay between observed attributes.

The absence of some traditional features such as colour, line ratios or morphology among the dominant features is striking and will be discussed in Sections 4.3, 4.4 and 4.5 respectively.

These nine independent components are used instead of the initial 47 attributes for Cluster Analysis, very substantially reducing the dimensionality of the dataset while keeping all the physical information.

Properties of the galaxies in the ten groups

The cluster analysis divided the galaxies into ten groups, K1-K10. The distribution of galaxies within these groups is given in Table 2. It appears that the four groups K4, K5, K3, K8, in decreasing importance, already gather 82% of the objects (52% with K4 and K5 only), the K1 group being very small.

Figure A1 shows the distribution of the ten groups in the nine ICs. IC1 and IC3 do not seem very discriminant, but all the others explain the specificities found by the CA algorithm. Many dimension reduction techniques creates new variables (here the IC components) which are linear or non-linear combinations of the original features. They are a mathematical representation of the data that is not intended to have a physical meaning. This is a mathematical space in which the data appear structured in groups. The physical interpretation should aim at understanding these groups in the physical space, not in the IC space. It should be clear that traditional approaches fail since each group is characterized by a complex multivariate physics. It is certainly simpler to summarize galaxies as blue, star forming, or spiral, but this is obviously very limited. The solution is thus to have more complex descriptions of galaxy classes, and use models and numerical simulations to interpret the results of multivariate analyses.

The boxplots shown in Figure A2 summarize the statistics for the ten groups of the 47 attributes used for the clustering analysis plus SFR, specSFR, Sersic r50 R and the redshift. The order and the names of the groups given in Table 2 is arbitrary and has been chosen to smoothen the variation of the boxplots from group K1 to K10 for most observed features. This nomenclature thus bears absolutely no physical or temporal evolutionary relationships between these groups.

It is interesting to note that the dispersions are nearly always relatively small, indicating that the groups found by the cluster analysis are quite homogeneous. This is particularly striking for the biggest groups, K4 and K5. There are often large overlaps, but also clearly separated properties distributions between groups in many instances. For some properties, such as Hα abs , J-K and H-K, there is very little variation from group to group.

The small group K1 stands out in many features, especially in EW(OIII 5007 ). This group is however close to the groups K2 and K3 for many properties, such as Ca K abs , Sersic n R or Lick G4300. With K4 they share a low velocity dispersion with a clearcut split with respect to the other groups K5 to K10. This split is also visible for Sersic n R and also somehow on the redshift diagram.

The specific SFR shows a remarkably regular decrease from K1 to K10, while D4000 n has the opposite behaviour. Groups K6 and K7 gather galaxies which are far much larger than all the other ones. These two groups are not so different in the other observable features, except may be that K7 has a higher average redshift, pointing to possible measurement errors. It is important to recall that most attributes are determined automatically, so aberrant values may appear. However, these two groups are not very small, making such errors (hopefully) improbable. Also, such large sizes are not uncommon (e.g. Guo et al. 2009). Further analyses of the galaxies of these groups, for instance using images, should be done before any physical interpretation. This is beyond the scope of this first paper, we simply emphasize that the CA analysis is able to separate these objects into distinct groups.

Colour-magnitude diagram

The well known bimodality of galaxies is seen on the colour-magnitude diagrams in Fig. 1 with a crescent shape of the distribution of the whole sample, with the so-called red and blue branches (or sequences).

The groups K5, K8, K9 and K10 clearly belong to the red branch, while K2 to K3 are essentially on the blue one. The group K4 spans both branches, but peaks mainly on the red branch and includes the region in between often called the green valley. Interestingly, the groups K6 and K7 that have the largest galaxies in our sample, are both part of this green valley. The very small group K1 appears peculiar, with a very blue part, and 25% of its galaxies belonging to the red branch.

The correspondence between our groups and the rough usual division in red, green or blue regions of the plot is quite good despite the fact that no colour nor magnitude are involved in the independent components used for the classification. This is very likely due to a certain degree of correlation among some properties of galaxies, which multivariate clustering, by determining more subtle and objective categories, is able to stress. Some groups extend to several regions, such as the big group K4. However, the bivariate colour-magnitude diagram alone cannot reasonably encompass all the physics of galaxies. In this respect, it is interesting to compare our groups 2014) on which some properties, like SFR or EW(Hδ), are plotted. Clearly these properties do not peak and are not confined to only one of the blue, red or green regions, similarly to our groups.

Emission-line diagnostic diagrams

Two emission line ratios were recommended by Baldwin et al. (1981, hereafter BPT) and are often used to discriminate between star forming and AGN-dominated galaxies. This diagram is also the basis for more refined empirical classifications based on theoretical population synthesis and photoionization models (Veilleux & Osterbrock 1987;Kauffmann et al. 2003;Kewley et al. 2006;Kewley et al. 2013). Two other similar diagrams, BPT-SII and BPT-OI, were proposed by Veilleux & Osterbrock (1987). They will be presented in subsequent papers. 2006) (right to left respectively). The bottom left zone is for star forming regions, the upper zone is for AGNs and Seyfert galaxies, and the bottom right zone is for LINERs.

The standard classification scheme on this diagram is based on equations of curves separating the different classes. These cuts are sharp, somewhat arbitrary. Its main goal is to distinguish between different ionization source (basically thermal or non-thermal) while the properties used for our clustering analysis are not limited to this peculiar aspect of galaxies. It also appears that the classification based on these diagrams is not clearcut. For instance, "the current LINER classification scheme encompasses two or more types of galaxies, or galaxies at different stages in evolution" (Kewley et al. 2006). Cid Fernandes et al. (2010) proposed a new cut between Seyfert and LINERs to resolve this ambiguous class. They also present an interesting and critical discussion of the classifications based on the BPT diagrams. An important reminder in their discussion is that the star-forming region delimitation is rather arbitrary in the lower part of the diagram where most of the galaxies lie. Also, the upper right part of the diagram is not composed of pure AGNs. Finally, different cuts are proposed by different authors (e.g. Kewley et al. 2001;Lamareille 2010). Interestingly, all these works tend to suggest that more parameters are probably required to fully understand different kinds of galaxies and ionization processes (Richardson et al. 2016).

Indeed, an objective classification with soft frontiers performed with the diagnostics line ratios indicates a somewhat different picture with only four categories, some including for instance sub-divisions like strong or weak AGNs (de Souza et al. 2017). However while being multivariate, this study is only three-dimensional.

An additional limitation of the classifications based on single kinds of ionization processes is that the galaxies are indeed very probably a mixture of different regions (Belfiore et al. 2016).

These limitations on the significance of the cuts should be kept in mind when comparing our multivariate clustering results and diagnostics diagrams (Fig. 2). We see that K2 and K3 have no AGN, hence are pure star forming galaxies, whereas the very small group K1 is quite peculiar with two peaks, one in the pure star forming region, the other one in the pure AGN zone. Groups K4 and K5 are intermediate between star forming galaxies and AGNs. All the other groups (K6, K7, K8, K9 and K10) are LINERs (Fig. 2).

Our result is in the line of the discussion found in the literature as presented above in the sense that our groups can be clearly placed in the different zones of the BPT diagram with some fuzziness. This agreement is remarkable because no line ratios are involved in the ICs used for our classification, only H α , [OIII] and [NII] but not Hβ.

Morphology

In Fig. 3, we plot the density distribution of the concentration index C = Sersic r90 R / Sersic r50 R for all the ten groups. Late-type morphologies are characterized by a low C (< 2.6, Strateva et al. 2001). Note the strong dichotomy with two peaks at C≃ 2.5 and 7.5.

The groups K6 and K7 are nearly entirely made of early-type galaxies while all the other groups are mixtures of both categories, with a higher fraction of late-type galaxies in K2 and K3 and of early-type ones in K5, K10, K8, K9.

The dichotomy between early-and late-type galaxies is thus roughly recovered despite the fact that no morphological indicator such as the concentration index and Sersic n R are present in the ICs. This absence probably explains why some of our groups have a mixture of both types, but it is also an indication that the morphology is not very discriminative when the whole multivariate physics of galaxies is considered.

DISCUSSION

In this first paper, we intend to focus mainly on the technique used for the multivariate classification of our large sample, and to provide a rather general description of the physical properties of the ten groups we have found to demonstrate that the classification obtained from the IC components is physically meaningful. Naturally, the multivariate grouping seems more complicated than traditional classifications, so that more detailed studies on some specific aspects, such as nucleus activity or SFR, will be made in subsequent papers. For now, from the results described in the previous section, we can distinguish four categories of groups that mainly match usual classes of galaxies.

Groups K1, K2, K3

The three groups K1, K2, K3 belong to the so-called blue sequence and are active sites of star formation as supported by high Hβ abs , Hδ abs , Hγ abs , specSFR, low metallicity and high EW(OIII 5007 ).

They are also the less massive objects of the sample, especially K1. They have higher values of Hδ abs thus containing young stellar populations but the widths of its distribution are not similar in these groups. There is maximum scatter in K1 and minimum scatter in K2. They all have a high SFR, but the very small group K1 has a particularly high specSFR. The large scatter in EW(OIII 5007 ) in the latter group indicates that the galaxies belonging to this K1 group contain gas at high temperature to form new generation of stars. Hence K1 has a stellar population which is a mixture of various ages as a result of star bursts of recent origin.

Though K2 and K3 are dominated by late type galaxies, they also have small fractions of spheroidals. The galaxies in these groups have a low Sersic n R lying between 0 and 2 which indicates that the bulges in the spheroidals are not also very pronounced, i.e., these galaxies may be in a formation stage. The two K4 and K5 groups gather about half of our sample. From the colour-magnitude diagram (viz. Fig. 1), it is clear that they are mainly on the red branch with some extensions towards the green and blue regions, especially for K4. However, their Sersic n R indexes lie between 2 -4 which show that some of the galaxies are young and in the formative stages of their bulges and some of them have well developed bulges. The SFR and metallicities in these galaxies are intermediate between highest and lowest values. The two groups differ largely in their light element abundances: K5 has larger value than K4. Regarding Hδ abs the values lie between 0.5 and 2.5 with pronounced scatter. Groups K4 and K5 have globally intermediate values of SFR and more particularly of specSFR with respect to the other groups. Most of the properties like metallicity, SFR, Hδ abs , specSFR, D4000 n and C have large scatters which indicate that K4 and K5 have a certain mixture of old and young populations, as well as of active and passive galaxies with well defined bulges for K5 and more pseudo bulges for K4.

Groups K4, K5

Galaxies in the K5 group have higher velocity dispersion hence are more massive than those belonging the group K4. The group K5 is more abundant in helium enriched population, which is a signature of second generation stars. So K5 has slightly older population than K4. This is also evident from the D4000 n peaks: K5 has a larger peak at higher D4000 n whereas the opposite is true for K4. Also C values show that the two groups have populations of late type spirals, but in much higher proportion in K4.

As a conclusion, despite the fact that K4 and K5 galaxies mainly belong to the red sequence, and would consequently be considered as quenched, they have somewhat average values for most features, and possess some star forming objects.

This result may seem inconsistent, but this is the outcome of a multivariate analysis that gathers galaxies sharing not merely a few properties, but many. However, we cannot exclude that this apparent inconsistency might be due to two possible limitations of our study. The first one is that the clustering algorithm used (K-means) may not be the most suitable technique for the sample, and in particular for the subsamples made with K4 and K5 galaxies. The second explanation might be that even with the information contained in the 47 attributes used in our study and summarized in the 9 ICs used for the classification, it is not possible to distinguish sub-classes within these two big groups. This may be due to several reasons: lack of more distinctive information, uncertainties that smear out differences among variable values of different groups, or the fact that the observables used here are integrated values over an entire galaxy probably mixing up several different regions within these big and complex systems. Further investigations will be made in another paper.

Groups K6, K7

Groups K6 and K7 are globally similar with mostly average values. They are often close to groups K4 and K5, but they have remarkable high C together with very high Sersic r50 R and Sersic r90 R. Hence they are made up of very big and large galaxies, and thus spheroidals. They occupy the same region in the colour-magnitude diagram, essentially corresponding to upper part of the so-called green valley. They have a relatively old population as based on D4000 n , but with a small fraction of younger stellar populations and some star formation ongoing.

These galaxies thus appear to be in a transitional phase of being quenched but they are really big and massive galaxies which have no equivalent in the red branch members in our sample (red branch that is usually considered as the end fate of all galaxies). The groups K6 and K7 may thus represent the dead end of big galaxies.

Groups K8, K9, K10

The three groups K8, K9 and K10 are rather specific: they have the lowest specSFR, high velocity dispersion, high metallicity, high Sersic indexes and C >> 2.6 for most of their members. All this indicates that they are early type spheroidal galaxies and quenched. They are more concentrated at the center and massive in nature. They have a high abundance of oxygen (viz. Fig. A2, EW(OIII 4363 ), especially in K10) which might be due to the explosion of massive supernovae (Pop II objects). They have low Hδ abs values with small scatters and high D4000 n values. They have well developed bulges as seen from their high Sersic n R values (∼ 4, Fig. A2). The forbidden line are pronounced in K10, indicating that these galaxies have enough neutral gas.

Independent components

We have seen that the nine ICs used to perform the cluster analysis represent a complete physical description of galaxies: velocity dispersion, ionisation, metallicity, surface brightness and structure. It is important to realize that they have not been selected on this basis. It should also be clear that using the nine ICs describing five properties is not equivalent to selecting five observables representing these properties, or even nine to match the number of ICs, since the latter are linear combinations of the observable features that reveal some latent multivariate relationships.

The most remarkable outcome is that the groups built from the ICs very satisfactorily match the traditional categorization of galaxies, despite the fact that none of the traditional observables is involved in the nine components. We might see this correspondence from two points of view: either the multivariate analysis is unnecessary, or, as we tend to think, the traditional approach is only an approximation of the complexity of galaxies and is thus incomplete. We obtain more groups than the usual classification, but this is because we are not able to see into a high-dimension space.

Among the observable features involved in the ICs, it is striking to see that the equivalent widths of several lines are important. Also the velocity dispersion and the line widths take their share. The size (radius) and mass (with the J magnitude) seem to matter. On the contrary, the absence of any colour, of the concentration index, related to the global morphology, or even the index D4000 n , is remarkable.

SUMMARY AND CONCLUSION

In this study, we have classified a large dataset of galaxies with a large number (47) of morphological, photometric and spectroscopic parameters compiled from VAGC/SDSS data archive. We have used two sophisticated statistical methods. Firstly we have performed a dimension reduction using the ICA analysis which is appropriate for a non-Gaussian dataset such as this one. Secondly, k-means cluster analysis has been used to find structures in the parameter space defined by nine Independent Components to obtain ten coherent groups.

The most correlated observable features involved in the nine ICs represent a complete physical description of galaxies (velocity dispersion, ionisation, metallicity, surface brightness and structure). Despite neither morphology nor colour nor emission-line ratios are among the most influential features in the ICs, our ten groups can be essentially identified as early-or late-type dominated and placed in the colour-magnitude and the line diagnostic (BPT) diagrams in agreement with the literature.

However, our classification is not empirical and is based on a rather complete physical description of the galaxies. The ten groups then allow for a more refined interpretation of galaxy diversity.

Recovering known divisions of galaxies gives us confidence in our unsupervised multivariate clustering analysis, and subsequent papers will be devoted to the astrophysical interpretation and implications of each of the groups we have identified. In particular, the mixture of traditional populations in some of our groups must be understood in the multivariate frame, compared to models and numerical simulations of galaxies that only can provide a multivariate physical understanding, without forgetting the context of evolution. Concentration Index: ratio between Sersic r90 R and Sersic r50 R uz Magnitude u minus magnitude z (u-z) J -H H -K a See http://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/SDSS line.html for more details. b The flux f in nanomaggies is defined in Blanton et al. (2005) as magnitude= 22.5 -5 log 10(f ).

  Details of algorithm and applications are found in MacQueen (1967); Chattopadhyay et al. (2009b); Chattopadhyay et al.

Figure 1 .

 1 Figure 1. Colour-magnitude (u-z vs. U) diagrams of the whole sample and for each of the groups.

Figure 2 .

 2 Figure 2. The distribution of galaxies in the groups K1 to K10 in the BPT diagram: [OIII]/Hα vs [NII]/Hβ. Curves are from Kewley et al. (2001); Kauffmann et al. (2003); Stasińska et al. (2006) (right to left respectively). The bottom left zone is for star forming regions, the upper zone is for AGNs and Seyfert galaxies, and the bottom right zone is for LINERs.

Figure 3 .

 3 Figure 3. Density distribution of the concentration index C of galaxies in the groups for K5 to K7.

Figure A2 .

 A2 Figure A1. Boxplots for the IC components.

Table 1 .

 1 Observed attributes with highest correlation coefficients (in parentheses) for each Independent Component.

	IC	Most influencial attributes
	ICA1 EW(OII3729)(0.83), Lick Fe5015 (0.53)
	ICA2 v disp (0.97), J (-0.7)
	ICA3 σ balmer (-0.60)
	ICA4 EW(OIII5007) (0.98)
	ICA5 Lick Fe5015 (0.83)
	ICA6 Sersic r90 R(-0.99)
	ICA7 σ forb (-0.60), σ balmer (0.42)
	ICA8 Sersic amp R (-0.99)
	ICA9 EW(NII6584)(0.79), EW(Hα)(0.75),
		EW(SII6731)(0.72)

Table 2 .

 2 Distribution of galaxies in the ten groups found with K-means.

	Group Number of galaxies Percentage
	K1	1,375	0.38
	K2	16,325	4.50
	K3	68,050	18.75
	K4	109,188	30.08
	K5	79,576	21.93
	K6	17,336	4.78
	K7	7,846	2.16
	K8	40,045	11.03
	K9	10,552	2.91
	K10	12,630	3.48

Table A1 .

 A1 All the parameters of the dataset used for the analysis. a

	Parameter v disp	Description Estimated velocity dispersion from spectrum
	u J Sersic amp R Sersic n R sigma balmer sigma forb	u absolute magnitude (log of intensity) J absolute magnitude The best fit to the variable "A" in band R (nanomaggies b /arcsec 2 ): describes the radial distribution of light The best fit to the Sersic index "n" in band R Velocity dispersion (σ not FWHM) measured simultaneously in all the Balmer lines in km/s Velocity dispersion (σ not FWHM) measured simultaneously in all
	oii 3729 seqw neiii 3869 seqw oiii 4363 seqw (EW(OIII 4363 ))	the forbidden lines in km/s The equivalent width of the continuum-subtracted emission line with the other emission lines subtracted off " "
	oiii 5007 seqw (EW(OIII 5007 )) hei 5876 seqw oi 6300 seqw h alpha seqw (EW(Hα)) nii 6584 seqw sii 6731 seqw Hδ abs Hγ abs Hβ abs Hα abs Ca K abs Ca h abs Na d abs Lick CN2 Lick Ca4227 Lick G4300 Lick Fe4383 Ca4455 Lick Fe4531 Lick c4668 Lick Hb Lick Fe5015 Lick Mgb Lick Fe5270 Lick Fe5335 Lick Fe5406 Lick Fe5709 Lick Fe5782 Lick NaD Lick Hd A D4000n SFR specSFR Sersic r0 R	" " " " " " The best fit to the variable r 0 in R band (arcsec) Equivalent width for absorption lines " " " " " " Stellar absorption line (Lick) index " " " " " " " " " " " " " " " " Å The break in the spectrum at 4000 Star Formation Rate Specific Star Formation Rate
	Sersic r90 R C	90% light radius of best fit model in R band (arcsec)

http://sdss.physics.nyu.edu/vagc/

http://classic.sdss.org/dr7/

http://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/raw data.html

ACKNOWLEDGEMENTS

We thank Fabrice Lamareille and Emmanuel Davoust for having compiled for us the sample used in this analysis. We warmly thank Malgorzata Siudek for invaluable discussions. TK and DFB are thankful to the Indo-French Centre for Applied Mathematics (Bangalore, India) for providing partial support. We thank the two referees for their detailed and very constructive reports.

Funding for the Sloan Digital Sky Survey (SDSS) has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington.

APPENDIX

A. COMPLEMENTARY TABLES AND FIGURES