P Fogel
email: paulfogel@hotmail.com

C Geissler
email: cgeissler@advestis.com

P Cotte
email: pcotte@advestis.com

G Luta

Applying separative non-negative matrix factorization to extra-financial data

Keywords: Interpretability, Features Engineering, ESG data, Clustering, Dimension reduction, Machine Learning

We present here an original application of the non-negative matrix factorization (NMF) method, for the case of extra-financial data. These data are subject to high correlations between co-variables, as well as between observations. NMF provides a much more relevant clustering of co-variables and observations than a simple principal component analysis (PCA). In addition, we show that an initial data separation step before applying NMF further improves the quality of the clustering.

Introduction

A question frequently encountered in economics is that of business classification. The financial industry (in particular index providers) has developed several standards for classifying companies based on their dominant industrial activity, such as the Global Industry Classification Sector (GICS R ). These classifications are widely used by investment fund managers as a tool for diversifying economic risk, or alternatively to create thematic portfolios focused on one or a small number of industry sectors.

The purpose of this paper is to outline a method for classifying companies based on their extra-financial profiles. To this end, we will use scores from extra-financial rating agencies. These scores quantify the behaviour of issuers in the three main pillars of Environment (E), Social (S) and Governance (G).

For the data-scientist, these data constitute a rich and complex field of study. The richness is linked to the variety of themes noted, ranging from the treatment of employees and suppliers, to the carbon footprint and the independence of the Board of Directors. These ratings provide information that cannot be derived from the industrial sector alone, and provide a new viewpoint to discriminate between companies.

The complexity stems from the fact that those data are not provided uniformly for all issuers, as some specific scores only make sense for very specific sectors. On the other hand, behavioural scores are both numerous (typically between 50 and 100) and sometimes highly correlated.

We describe here an unsupervised clustering experiment based on a matrix of extra-financial data. The data constituting the matrix are normalized between 0 and 100, and therefore lend themselves well to the use of non-negative matrix factorization (NMF) methods. This factorization provides a grouping of both co-variables and companies. The results are compared with those of a Principal Component Analysis (PCA). We also show, following a path initiated by [START_REF] Fogel | Applications of a novel clustering approach using non-negative matrix factorization to environmental research in public health[END_REF], that a tensor variant of NMF can be used to isolate groups of issuers characterized by either high or low scores on some co-variables, while ordinary NMF would mainly focus on groups of high scores.

Presentation of the data

We use ESG scores from the research of Vigeo-Eiris, whom we thank here for making these data available for this study.

The data consist of 77 scores spread over the three pillars E, S and G. They concern 699 European stocks. The scores are positive numbers between 0 and 100, with 0 representing the worst possible score and 100 the best. These quantities do not have absolute significance but are essentially used to establish a ranking between issuers.

Nomenclature of variables

We will denote by X ESG the matrix formed by the 77 variables or 'features' as columns. This matrix is of type (699, 77). The generic term X ESG i,j is the value of the feature j for stock i.

Vigeo's behavioural scores are grouped into 6 areas, with taxonomy being reflected in the denomination of variables:

• Environment: variables names starting with env .

• Corporate Governance: "cg "

• Social: Business Behaviour: "c s "

• Social: Human rights: "hrts "

• Social: Human Resources: "hr "

• Social: Community involvement: "cin " More information about the scoring methodology can be found at: https://vigeo-eiris.com/about-us-2/ methodology-quality-assurance/.

The interesting feature here is that we have a prior fundamental knowledge about variables that will be used as a benchmark to evaluate the relevance of any NMF-based clustering. Similarly, the fundamental information about the industrial sector the companies belong to will also be used as a prior knowledge.

Missing data

Variables are scores that rate the company's activity according to different criteria. Some of these criteria have no relevance outside a particular industry sector, and are therefore not reported for all issuers. The next two charts give a visual summary of the average fill rates grouped by sector (Figure 2) and by features (Figure 3). The sectors are presented along lines, by decreasing rate, while the variables are presented in columns, also by decreasing rate. As shown on chart 2, the consumer staples sector has the least amount of missing data. In contrast, several environment-related scores in the finance or telecommunications sector have little or no information at all, i.e. a fill rate close to zero. twenty two features out of 77 have a fill rate less than 75 %, the average being 81.1 %. Before performing any statistical analyses such as linear regressions, it is necessary to apply a missing values imputation strategy. A default strategy is to fill in the missing data for a variable with the mean or median, taken over the entire sample, of the non-missing values of that variable.

A first analysis using Principal Components Analysis (PCA)

Principal component analysis of a X data matrix focuses on the uncorrelated directions that explain the highest proportion of variance in the data. The directions are linear combinations of the original variables, which are mutually orthogonal. These principal directions are none other than the eigenvectors of the M = XX T matrix.

In the case of the X ESG matrix, we need to fill in the missing data with non-empty numerical values, because PCA algorithms are not meant to handle missing values. To do this, we apply the most neutral imputation strategy with respect to the objective of explaining variance: for each variable, a missing observation is replaced by the average of the variable, calculated on the set of non-missing data.

A first look at the results of a standard PCA allows the following observations to be made:

• The observation of the cumulated variance explained according to the rank of the component shows that it increases rather slowly. At least twenty two components are needed to explain 90 % of the total variance (Figure 4).

• The first component has coefficients that are all positive except one. All other components, which are mutually orthogonal, have positive and negative coefficients, at least half of which are greater in absolute value than 10 % of the largest coefficient. (Figure 5).

Except for the first component (thick blue line) for which all coefficients but one are positive, all other components have coefficients of both signs, which makes their interpretability difficult. 4 Factoring the data through Non negative Matrix Factorization (NMF)

Principles

Let us assume the available data are represented by an X matrix of type n × f , i.e. n rows and f columns.

We also assume that these elements are positive or null and bounded. With that in mind, we can re-write X (here shown as a 3 × 4 example matrix, with rows and columns indices starting at 0) as

  x 0 0 x 1 0 x 2 0 x 3 0 x 0 1 x 1 1 x 2 1 x 3 1 x 0 2 x 1 2 x 2 2 x 3 2   X =   w 0 0 w 1 0 w 0 1 w 1 1 w 0 2 w 1 2   W h 0 0 h 0 1 h 0 2 h 0 3 h 1 0 h 1 1 h 1 2 h 1 3 H T +   0 0 1 0 2 0 3 0 0 1 1 1 2 1 3 1 0 2 1 2 2 2 3 2   =X-W H T (1)
A non-negative factorization of X (Lee D. [1999]) then solves the following problem: Find two matrices W of type n × c and H of type f × c minimizing the quadratic norm of the last term of the right hand of Equation ( 1)

|| || 2 2 := X -W H T 2 2 (2)
H and W being subject to the constraints of only having non-negative entries.

A few remarks:

• The quadratic norm used here is restricted to numeric entries (i.e. represents the square root of the sum of squared non-missing values). This slight extension of the definition allows to deal with input matrices having some missing values.

• The objective function in 2 may also include additional by L1 and L2 penalization terms applied to the W and H matrices [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF]). The aim of this penalization is to eliminate noise coming from small and generally unstable coefficients, and to concentrate the signal on the most significant contributors. This objective is called sparsity.

• This decomposition is an approximation, not an equality. Any matrices W and H that are solutions to (2) minimize the quadratic error of the approximation.

• It is common in the literature to refer to rows of W as scores, and to those of H as loadings.

• The decomposition is not necessarily unique, even with the positivity constraint. One can indeed use any D matrix having non-negative coefficients, as well as its inverse D -1 , and then obtain identical decompositions of the form:

X = W DD -1 H T + .
Different examples of such a transformation would be:

-A change of unit in the variables, with a diagonal D matrix representing the equivalence coefficient between two measuring systems.

-A permutation in the components, D being then represented by a permutation matrix.

• c ≤ min(n, f ) is an integer representing the number of selected components, also called rank of the decomposition. Like in most clustering methods, this parameter represents a degree of freedom; determining the most adequate number of components requires an additional optimality criterion than the one for error, as the error generally decreases with the rank.

• The n rows of X will be referred to as observations, the f columns of X as variables or features, and the c columns of W as meta-variables or meta-features.

It is customary to consider the rows of X as n observations of objects from the same category (e.g. customers, patients, companies, etc). Each object has f attributes, and the data contained in the matrix correspond to the values taken by the attributes for the corresponding objects.

The main effect of this decomposition is to decrease the information necessary to describe an observation. The original observations of the X matrix can be recovered, at the cost of an approximation, as linear combinations with positive coefficients of the rows of the H matrix. A reduction in size has thus been achieved: part of the information contained in the original matrix, comprising f features per individual, can be approximately summarized by a smaller set of c features per individual.

A first NMF run

We start with a first NMF decomposition on the 500 × 77 scores matrix X. Two prior steps are handling missing values in X, and choosing a rank c.

We use the Python package nmtf developed by one of the authors. This package is public and available at https://pypi.org/project/nmtf/. It can be directly installed with the pip command: pip install nmtf.

The optimisation algorithm used in this package takes advantage of the extended definition of the quadratic norm and can handle missing values in input matrices, as opposed to the NMF version 0.17 of Scikit-Learn, which raises an error in case of missing entries in the matrix.

This ability to handle missing values avoids having to deal with the delicate issue of imputation of these missing values. In the field of extra-financial scores, a missing value is not simply data for which a measurement is missing because it did not take place. It can be data whose measurement simply does not make sense. To give a concrete example from the environmental domain, scores for waste reprocessing or contribution to deforestation are not relevant for some industrial sectors.

Another difference between nmtf and Scikit-Learn lies in the way sparsity is specified. In Scikit-Learn, the objective function is modified by addition of penalty coefficients that have an indirect effect on the sparsity of the W and H matrices. In nmtf, on the contrary, the sparsity objective on H or on W is provided as a parameter to the model, as a positive coefficient less than 1. For example, a parameter value of 0.9 on H will provide, if the resolution is possible, an H matrix with about 90% of zero entries.

The rank is set here at 6. This choice is not totally arbitrary as a first trial, it is exactly the number of topics set by Vigeo-Eiris for the taxonomy of the scores.

The solution X provided by the algorithm satisfies:

X -W H T 2 2 = 5.7 × 10 -2 X 2 2 .
In other words, the squared approximation error for rank c = 6 is 5.7 % of the squared norm of X. Referring to the performances of PCA on chart 4, it is necessary to aggregate the first 32 factors in order to explain a fraction of 100 % -5.7 % = 94.3 % of the total variance. This illustrates one of the main strengths of NMF: relaxing the requirement of factor orthogonality provides in return a much stronger explanatory power.

Clustering observations and features

An interesting induced property of NMF is the ability to cluster variables and observations in a natural way.

Regarding the clustering of features, the principle is to link each of the f variables with the meta-variable (indexed between 1 and c) having the highest coefficient on this variable. Similarly, each observation is assigned the index of the meta-variable with the highest contribution to the observation.

This grouping scheme has the disadvantage of being dependent on the particular system of units retained for the meta-variables, i.e. a transformation by an arbitrary D matrix. In order to get around this disadvantage, the authors of [5] introduce the notion of leverage, which represents "the ability of a row or column to specifically exert an influence on a particular element without a corresponding increase in the influence on other components". Leverage is a dimensionless quantity between 0 and 1. The clustering scheme based on maximum leverage is more stable than that based simply on the maximum coefficient.

The most visible benefit that we obtain from the non-negative matrix is the replacement of the initial 77 features with a reduced subset of 6 meta-features; the initial features can themselves be written as positive linear combinations of these meta-features.

Figure 6 shows the distribution of the coefficients (the 'loadings') of the initial 77 features over the 6 new meta-features, sorted by decreasing order. We observe that as the rank of the meta-features increases, fewer and fewer of the original features have significant coefficients on them.

Introducing Separative NMF

One of the few weaknesses of NMF is its asymmetry of treatment between high and low values. By construction, standard NMF will group observations by looking for the meta-features having highest scores. The attribution of a cluster will be more discriminant for companies characterized by a high score on a meta-feature (e.g. high environmental scores), and less discriminant for companies characterized by a low score. But when it comes to ESG scores, knowing that a company has a low score on a meta-feature is just as important as the information provided by a high score.

In [START_REF] Fogel | Applications of a novel clustering approach using non-negative matrix factorization to environmental research in public health[END_REF], an interesting variant of NMF is proposed in order to go around this obstacle. The method is called PosNegNMF. The idea is to separate the original matrix and to rewrite it as the sum of three components:

X = X + -X -+ X 0 (3)
with the following notations:

• X 0 is the "baseline", or the matrix of columns-wise medians of X.

• X + = M ax(0, X -X 0 ), the matrix of positive high deviations from X 0 .

• X -= M ax(0, X 0 -X), the matrix of positive low deviations from X 0 .

Having done this, PosNegNMF will perform a single decomposition of the matrix containing positive and negative deviations. This is equivalent to two independent decompositions, yielding two H matrices with a common W matrix. More precisely:

X = def [X + |X -] = W [H + |H -] T + (4) X = X + -X -+ X 0 = W (H + -H -) T + X 0 + (5)
. The notation [A|B] represents the r × 2c horizontal concatenation of two r × c matrices A and B. The first equation describes the NMF decomposition of X. The second one reconstructs the original matrix as the sum of baseline, high and low deviations.

The new approach proposed here relies on the more general k(= 3)-dimensional algorithmic version of NMF, known as Non-negative Tensor Factorization (NTF, [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF]), to approximate the original matrix by a three-dimensional n × f × 2 tensor X:

X = def [X + |X -] = W ⊗ H ⊗ Q + (6)
Note that all elements of X are non-negative. Similarly to the NMF-based algorithm described above, this decomposition provides a canonical clustering based on the meta-features provided by H matrix. Let us emphasize that this H matrix is common to the two slices of the tensor X. In other words, what we obtain is not the same as PosNegNMF, which performs two NMF decompositions on X + and X -respectively. The H matrix describes common features having an influence on either higher or lower deviations of the scores from their median.

PosNegNMF yields an approximation X + -X -= W H T with H = H + -H -. Thus, each column of H may have signed entries, which complicates interpretation since the features under study are non-negative scores. In contrast, our NTF-based approximation yields a non-negative H matrix in which each column describes common features that have an impact on either higher or lower deviations of the scores from their median, as determined by the sign of the difference Q[1, column] -Q[2, column]. The non-negativity of the H-loadings makes interpretation as straightforward as it is for the standard NMF. For this reason, we propose to call our approach Separative NMF, or S 2 NMF, since it aims at separating patterns of features with either positive or negative deviations.

Precision and sparsity

We will now evaluate the comparative advantage brought by S 2 NMF, by varying the level of control of parsimony of the variables. More precisely, we observe how the approximation error and the composition of the meta-features evolve as the model is constrained to work with the most significant features. In order to focus on the sole effect of the parsimony, we start from a fully populated data matrix. Missing data are replaced for each feature, by the average of non-missing values of this feature. Doing this prevents from any unwanted numerical side-effect that can appear in the presence of a sparsity constraint applied to a matrix with empty entries.

To perform this experiment, we proceed as follows:

• The number of components is set to 6.

• We run in parallel an ordinary NMF model with its S 2 NMF counterpart.

• Both models are first launched with a sparsity target of 0.9 on the matrix H of meta-features. To achieve such a sparse representation, each vector of H is approximated after each update of the numerical resolution by a close vector (in the sense of quadratic norm) with the required degree of sparsity [START_REF]Block coordinate descent for sparse nmf[END_REF], [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF]). By doing so, some features show zero coefficients on all meta-features in matrix H. This is the sparse version of the model.

• The features showing only null loadings at the previous step are dropped from the matrix X, and the model is run again on the reduced matrix with no sparsity constraint: this is the restricted version of the model.

Figure 7 summarizes the steps of this process. Table 2 presents the results of this experiment. The columns of the tables are defined as follows:

• Dimension is the number c of components, here set at 6 for all models.

• Features refers to the number of features effectively included in the factorization.

• H-sparsity is a tuning parameter of the NMF/S 2 NMF algorithm, representing a target level of sparsity, or a proportion of null values, in matrix H.

• The Error column is the relative squared error of the approximation.

We see that the ordinary NMF model collapses when the sparsity is constrained at the high level of 90%. The mean squared error of the factorization becomes greater than the variance of the original data. In order to recover an acceptable error level, we have to drop 41 features out of 77 that do not contribute to the composition of any meta-feature, leaving us with 36 significant features. By contrast, starting directly with a Separative NMF approach in sparse mode, effectively reduces the error to an acceptable level of 9%. Forty five features are flagged as irrelevant, leaving the dense model with a remaining 32 columns.

Clusterization efficiency

The advantage of S 2 NMF over NMF is more easily seen by looking at the ordered heatmaps presentations of the matrix. Figures 8 and9 show NMF-and S 2 NMF-based heatmaps of the restricted matrix. The restriction takes out the features that have been eliminated in the first step with a target sparsity set at 90 %.

It is visually clear that S 2 NMF has a greater ability than NMF to discern and separate clusters, at both ends of the score range (blue or red aggregates).

If we try to formalise this visual intuition, the mathematical translation is that the distribution of scores within clusters is more homogeneous than in the totality of the data. This homogeneity can be quantified thanks to the notion of entropy.

Let us recall its definition: for any probability distribution p on a finite set:

X = {x 1 , . . . , x n }
the entropy of the distribution is given by

E(p) = i=n i=1 p(x i ) × log 1 p(x i ) (7)
Information theory tells us that log(1/p(x)) is (up to an additive constant) the optimal number of bits required to encode the value x, given its probability of occurrence p(x). From this viewpoint, E(p) appears to be the expected encoding length of a random pick from X. The minimum possible entropy is 0, attained for a distribution concentrated at 100 % on a single value. The maximum entropy is realized for a uniform distribution on X, and its value is then log(n). Let us use a normalized entropy defined by

E * (p) = E(p) log(n) (8)
The efficiency of a NMF/S 2 NMF-based clustering can be assessed by looking into the diagonal blocks of the sorted matrix. We expect to observe less dispersion, hence less entropy, among the values a) taken by the features belonging to a cluster i, and b) restricted to those observations also belonging to cluster i. The rationale here is to group together observations having high scores on a given meta-feature (read in matrix W ), hence having likely high scores on the features having a high loading in this meta-feature (read in matrix H). The advantage of S 2 NMF vs NMF is in the ability to group together not only the high values, but also the low values which can contain important information as well.

"Clustering entropy" shown in Table 3 represents the relative entropy of the clustering induced by each model. It is calculated as follows:

• The normalized entropy of each diagonal cluster (i, i) is calculated. The i-th diagonal cluster is a submatrix restricted to features in cluster, i and to observations having their dominant score on feature i.

• The clustering entropy is the average of relative entropies measured on all diagonal clusters, weighted by their sizes.

• This average entropy on diagonal clusters can be compared with the normalized entropy of the whole matrix. The difference defines the "Entropy delta".

Entropy delta Restricted NMF -0.77% Restricted S 2 NMF -8.03% It is not a surprise to observe lower average entropies in the clusters, compared to the respective original matrices. This is an indirect proof of the efficiency of general NMF methods applied to this particular dataset. Moreover, we see that the loss in entropy is higher (-8.03% vs -0.77%) with S 2 NMF compared to NMF, and this is a numerical confirmation of the visual difference between the respective heatmaps.

Interpreting meta-features

Let us select as final model a 6-components S 2 NMF, after elimination of the features from an initial run with sparsity 0.9. By taking a closer look at the definition of the meta-features, some coherence appears in the most influential features, i.e. those having the highest (positive) coefficients. 3 show the 6 most influential features of each meta-feature. Interestingly, the topics of these influential features are interrelated. We can therefore assign a main theme to each meta-feature. For example, the theme related to feature V1 is mostly about Human Rights (hrts) and Human Resources (hr), while V2 is clearly concentrated on Environment (env). S 2 NMF has traded dimension for interpretability, by drastically lowering the dimension (from 77 to 6) at the price of getting more fuzzy, but still providing meaningful features.

Rank

Chart 10 shows the correspondence, or incidence matrix, between the 6 clusters identified by S 2 NMF and the 11 industrial sectors of the stocks. Each cell shows the proportion of the stocks in a given sector (in rows) classified in a given cluster (in columns). It is visually clear that the two clusterings are not mutually redundant. S 2 NMF provides a data-driven, endogenous way of grouping stocks with respect to their ESG scores. For example, we see that the Utilities sector has about half of its members classified on the cluster 1, whose dominant theme is Human Rights (table 3). Similarly, the cluster 2 (Environment) is only dominant for the Health Care sector, and poorly represented in Utilities and Materials.

Conclusion and future work

We have illustrated here the application of NMF to extra-financial score data. In contrast to principal component analysis, which focuses on risk reduction by creating mutually orthogonal factors, NMF allows the information contained in the data to be summarized faithfully using meta-variables, where the original variables appear as linear combinations with positive coefficients. These meta-variables are not constrained by mutual orthogonality; they induce a natural clustering of the original variables and observations, at the cost of much lower error than would be produced by a PCA of the same dimension. Moreover, the S 2 NMF variant presented here allows, by symmetrizing the role of low and high scores, to further improve the results in terms of error and cluster relevance.

An induced result is the interpretability of the meta-features thus created, as shown by the comparison of their definitions with the prior knowledge of the meaning of the variables, which was not taken into account by the decomposition.

Finally, we have shown that an initial application of NMF with a high level of sparsity allows for the removal of variables that are insufficiently relevant.

Future directions of research are:

• The determination of heuristics independent of any prior knowledge of the domain, to optimize the number of components retained in the factorization. Possible heuristics are combinations of sparsity, robustness, precision and interpretability metrics.

• Generalizing the use of S 2 NMF as a pre-processing of variables, in a features engineering step upstream of a supervised learning process. The aim here is twofold: to reduce computation time and, above all, to reduce the risk of overlearning. Replacing a large number of partially redundant features by a much smaller number of interpretable features containing globally the same amount of information is a good policy in machine learning, regardless of the domain.

• The clustering induced by S 2 NMF, being independent from sectors, can induce a smart strategy for the imputation of missing data.

• S 2 NMF is applied to 2-dimensional data, and involves the use of 3-dimensional tensors to perform the factorization. A seemingly interesting generalization is to transpose respectively PosNegNMF and S 2 NMF to 3-dimensional data (with the addition of a time axis). The same way S 2 NMF uses an additional dimension to factorize the duplicated matrix (as described in equation 3), the 3-dimensional factorization would be achieved with 4-dimensional tensors, and yield PosNegNTF and SNTF.

We thank Vigeo-Eiris for making some datasets available to the authors. We also thank Vincent Margot for useful discussions and remarks.

Figure 1 Figure 2 :

 12 Figure 1 Figure 2: Each cell represents the average fill rate of a given feature within a given sector.

Figure 3 :

 3 Figure 3: The average fill rate by feature.

Figure 4 :Figure 5 :

 45 Figure 4: Cumulative explained variance of the first n components or eigenvectors.

Figure 6 :

 6 Figure 6: Coefficients of the 6 meta-features sorted by decreasing order.

Figure 7 :

 7 Figure 7: A schematic representation of the various factorizations applied to the initial data. Rectangles with rounded corners represent input or output data, while rectangles with slanted edges represent algorithmic treatments.

Figure 8 :Figure 9 :

 89 Figure 8: The original matrix reordered according to NMF clusters, and restricted to the 37 features preserved by the sparsity condition.

Figure 10 :

 10 Figure 10: Incidence matrix between industrial sectors and S 2 NMF clusters.

Table 1

 1 illustrates the fill rate (number of stocks for which the variable is filled in, divided by the number of stocks). The fill rate is calculated as an average per industry sector.

		Average fill rate
	Consumer Staples	88.41 %
	Consumer Discretionary	85.52 %
	Materials	84.36 %
	Health Care	83.89 %
	Utilities	83.54 %
	Industrials	82.37 %
	Energy	79.22 %
	Real Estate	79.02 %
	Communication Services	78.03 %
	Information Technology	77.18 %
	Financials	76.21 %

Table 1 :

 1

Sectors sorted by decreasing average fill rate, the average being taken over all features. The global average fill rate is 81.1 %

Table 2 :

 2 Comparison of NMF and S 2 NMF in sparse and dense versions.

Table 3 :

 3 Relative entropies for NMF and S 2 NMF restricted models.

Table 3 :

 3 Assigning themes to meta-features based on their dominant features

Table


Meysam Mayahi
email: meysam.mayahi@inria.fr

Antonio Costanzo
email: antonio.costanzo@inria.fr

Valeria Loscrí
email: valeria.loscri@inria.fr

Anna Maria Vegni
email: annamaria.vegni@uniroma3.it

An Interference to Noise Ratio Handover mechanism for Mobile Visible Light Communication Networks

Keywords: Visible light communication, Interferencebased handover, Adaptive Modulation Scheme

Visible Light Communication (VLC) is recognized as a complement communication technology that meets many hurdles of Radio Frequency (RF) systems. For example, VLC overcomes the limited RF spectrum resources for mobile users, even though its application demands an effective mechanism to handover among the different devices, in order to keep the connectivity consistent. In this paper, a new handover mechanism is introduced for mobile VLC nodes, based on a real time evaluation of the Interference-to-Noise Ratio (INR). This mechanism is coupled with four adaptive modulation schemes to further enhance the data rate. The communication system is evaluated at 0.15 mps and 0.3 mps velocities for different parameters i.e., handover rate, delivered data per handover and handover delay ratio, through an extensive simulated campaign. An experimental validation carried out through software-defined approach, considering a small-scale scenario and low power Light Emitting Diodes (LEDs). A perfect results alignment between simulation and measurement show how suggested INR-based mechanism overcomes classical SNR-based handover for all performance parameters evaluated in this study.

I. INTRODUCTION

Since a couple of decades, the spectrum crunch issue is encouraging the alternative solutions to meet the lack of sufficient radio resources for growing number of mobile devices and applications [START_REF] Matheus | Visible Light Communication: Concepts, Applications and Challenges[END_REF]. Among the others, Visible light Communication (VLC) is considered as a reliable complement to RF technologies not only because of the exploitation of the existing lighting infrastructures but also because of the potential availability of the huge optical spectrum [START_REF] Haas | What is LiFi?[END_REF]. On the other hand, the presence of the Line of Sight (LoS) between the transmitters (i.e., Light Emitting Diodes -LEDs) and the receivers (i.e., Photo Diodes -PD, or cameras) seems mandatory in most of applications [START_REF] Cui | Line-of-sight visible light communication system design and demonstration[END_REF].

The main challenge for the LoS requirement is the mobility, where the mobile receiver in VLC has to maintain a high-speed connection in uninterrupted manner within the coverage area of a single access point (AP), as well as when it is moving across multiple transmitters [START_REF] Haas | What is LiFi?[END_REF]. The principal solution for switching the communication session from one AP to another is the handover mechanism. The handover procedures are generally pushed by the degradation of the quality of the current link, followed by the selection of a higher quality host link (if there is any), and finally executing the handover through exchanging control packets among the nodes that share the link [START_REF] Dang | Handover Procedure and Algorithm in Vehicle to Infrastructure Visible Light Communication[END_REF]. The link quality evaluation is determined by the error probability ratio, which follows the distribution of the link disturbances such as noise and interference [START_REF] Rahaim | Optical interference analysis in Visible Light Communication networks[END_REF].

Back to the VLC literature, the maximum received power indicators, such as the signal-to-noise ratio (SNR) with Gaussian noise distribution are the well-known link quality metrics, even though overhearing the disturbances from other VLC transmitters offer less accurate performance evaluation. In [START_REF] Wu | Smart Handover for Hybrid LiFi and WiFi Networks[END_REF], the maximum received power has been used to make soft horizontal handover among coordinated multipoints (CoMP), by adjusting the handover power margin and time-to-trigger parameters. The received signal intensity (RSI) extension in [START_REF] Vegni | Handover in VLC systems with cooperating mobile devices[END_REF] was inspired by the RF received signal strength (RSS) concept, and exploited in RSI-based handover technique for mobile VLC systems. Signal to interference and noise ratio (SINR) is used in [START_REF] Demir | CoMP-Based Dynamic Handover for Vehicular VLC Networks[END_REF] to manage vertical handover between Li-Fi and Wi-Fi and to reduce the handover rate through the concept of handover skipping. The SINR evaluates the communication link more precisely, specially when the noise and interference power levels are comparable. Consequently, the interference is modeled as a Gaussian random process too, with finite independent signal aggregation, where no individual signal dominates over the rest, [START_REF] Win | A Mathematical Theory of Network Interference and Its Applications[END_REF]. At the same time, the scenarios where the disturbances are dominated by one or multiple sources are commonly present in VLC systems. When the power of individual interferer or accumulative interference is at least one order of magnitude larger than the noise, the system is interference dominant. In this case, INR ought to properly describe the interference distribution [START_REF] Rahaim | Optical interference analysis in Visible Light Communication networks[END_REF]. Nevertheless, the Quality of Service (QoS) constraints was exploited for handover management in hybrid RF/VLC indoor systems [START_REF] Hammouda | Link Selection in Hybrid RF/VLC Systems Under Statistical Queueing Constraints[END_REF].

In this paper, INR is employed for the first time as a handover decision metric. Link quality analysis based on INR allows to identify the disturbance distribution, which can belong to two cases i.e., • higher INR, which refers to an interference dominated scenario, and •• lower INR, where the noise and interference have similar strength and might be treated in the same way. An important parameter to be considered in interference dominated case is the Interferer-to-Interference Ratio (IIR), in order to verify the dominant interferer among all the interferers. Later on, the dominant interferer will be selected according to the application requirements in terms of reliability and stability of the communication system. Regarding the data rate, an additional improvement is obtained via adaptive modulation technique switching among On-Off Keying (OOK) and Pulse Shift Keying (PSK) with different modulation orders. The evaluation system is performed through extensive simulations based on NS3 simulator tool [START_REF] Riley | The ns-3 Network Simulator[END_REF] and an experimental validation in a small-scale version scenario. The good match of the simulation and experiment results proves the faithful performance evaluation of the proposed handover and transmission mechanism.

The reminder of this paper is organised as follows: Section II presents the system model via algorithms and analytical description. The validation of the system model is presented in Section III by comparing the results of simulations and experimental testbed in different scenarios. Finally, conclusions are drawn at the end of the paper.

II. SYSTEM MODEL

In this section, we first describe the VLC propagation channel and the front-end model, then we will explain in detail the proposed handover mechanism. Consider the VLC system at Fig. 1 consisting of N transmitters and a receiver. At the generic transmitter i, the LED (fixed at position X i ) transforms the electrical signals to the light intensity in order to provide a background illumination along with visible light communication across the coverage area up to the maximum distance (i.e., d max [m]) depending on the half power semiangle (i.e., ϕ max ). Intensity modulation constrains the instantaneous output signal power S i (t) [W] between P t,min and P t,max with average optical signal power P t,i = E[S i (t)]. At the receiver side, for handover mechanism analysis purpose the photo diode (PD) assumes to move along the vector x with constant velocity V P D and instantaneous position, X P D . It generates an electric signal proportional to the average optical power received from LED i, P r,i , estimated as:

P r,i = h Lam,i • P t,i , (1)
where the channel impulse response between LED i and the PD (i.e., h Lam,i ) is restricted to the Lambertian emission of order m i.e.,

h Lam,i = (m + 1)A 2πd 2 i cos m (ϕ i )T s (ψ i )g(ψ i ) cos (ψ i ), (2)
where A [cm 2 ] is the area of the PD which lays at distance d i [m] far from the LED i, T s and g are the gains of the PD optical filter and the optical concentrator, respectively at ϕ i degree of radiance and ψ i degree of acceptance angles.

In addition to the path loss, the VLC links suffers from the Additive White Gaussian Noise (AWGN) components including the thermal σ thermal and the shot noises σ shot , which are expressed through a total contribution σ n as: Given Gaussian noise dominant system, the performance of the link between transmitter i and the receiver is evaluated through the SN R i expression i.e.,

σ 2 n = σ 2 shot + σ 2 thermal , (3)
SN R i = (r • P r,i ) 2 σ 2 n , (4)
where r [A/W] stands for PD responsivity Interference in VLC systems arises when the PD detects the neighbour transmitter signals. A more precise link evaluation should consider both the interference and the background noise at the surrounding. Network interference is commonly modeled as a Gaussian distribution of independent processes with no signal prevailing [START_REF] Win | A Mathematical Theory of Network Interference and Its Applications[END_REF]. In such scenarios, the central limit theorem is respected and the variance of the aggregate interference current (i.e., σ 2

I [A 2 ]
) is used to define signalto-interference ratio (SIR i ) evaluation of the link between LED i and the photo diode i.e.,

SIR i = (r • P r,i ) 2 σ 2 I , (5)
Performance evaluation of PD -LED i link when the disturbances from both noise and interference are in the same order of magnitude is feasible through SIN R i i.e.,

SIN R i = (r • P r,i ) 2 σ 2 I + σ 2 n , (6)
The system is assumed to be interference dominant if the power of the interference is at least ten times larger than the noise power [START_REF] Rahaim | Optical interference analysis in Visible Light Communication networks[END_REF], then it is straightforward to define IN R as a dominant disturbance identifier i.e.,

IN R = N i=1 (r • P r,i ) 2 σ 2 n , (7)
where N interferers are assumed to disturb each other.

Algorithm 1 Interference-based handover Require: α, T c , t s measure IN R while IN R ≥ 10 do measure IIR if IIR ≥ β then measure SN R i if SN R i ≥ α then estimate t c if t c ≥ T c then T c = t c candidate i end if end if end if IN R = IN R -SN R i end while return T c , i if null then backoff t s end if
Gaussian interference assumption facilitates the link performance evaluation by employing the SN R i , SIR i and SIN R i interchangeably. However this assumption is invalid in the presence of a dominant interferer. In such scenarios, the interference is properly modeled by the distribution of the dominant interferer [START_REF] Rahaim | Optical interference analysis in Visible Light Communication networks[END_REF]. In order to specify the dominant interferer, the IIR is calculated as:

IIR = N max i=1 P r,i Σ i P r,i , (8)
Low IIR fits well with the normal distribution assumption, while higher IIR proves less normal distribution and the Gaussian assumption leads to a heavier tail in the distribution.

In extreme case when IIR → 1, the only noticeable interference comes from a single dominant interferer. Finally, once the dominant disturbance of the system has been identified, the interference-based handover mechanisms relies on the appropriate metric in order to evaluate the link quality. Algorithm 1 describes the main steps of INR-based handover mechanism. The mobile node (P D) in Fig. 1 is crossing multiple VLC coverage areas while moving along x -axis. The quality of the link degrades as it moves away from a fixed AP (LED i) due to the path loss and the front-end orientation. In order to maintain the connectivity link with a predefined Quality of Service (QoS) level, the receiver P D may assess the link quality frequently (i.e. every t s seconds) and modify the transmission parameters, if required [START_REF] Mayahi | Link adaptive protocol for V2LC[END_REF]. The INR-based algorithm evaluates the link performance based on the prominent component of the disturbances. To this end, it starts with IN R calculation to determine if the system is noise-dominant or interference-dominant. In order to consider the system interference dominant, the interference power must be 10 times bigger than the noise power. In case that IN R < 10, the noise and the interference power are comparable, and we treat them equally as a noise.

Algorithm 2 Adaptive Modulation Scheme Require: i, α while SN R i ≥ α do adapt the modulation Transmit Receive measure SN R i end while call handover At the interference dominant system, there could be single or multiple dominant interferer. The IIR is able to identify the dominant interferer if IIR value is bigger than a threshold β, which is the minimum amount of IIR the system can detect and it is experimentally selected according to the node density(ρ l ). The interferers with IIR ≤ β are not dominant and are considered as a noise. The dominant interferer i with IIR ≥ β will be nominated for handover if it satisfies the QoS requirements of safety applications, expressed in terms of reliability and stability. A reliable link has SN R i ≥ α [dB], where α is the minimum SNR level required for safety applications at maximum error probability, equal to 10 -3 [START_REF] Song | Performance Analysis of the IEEE 802.11p Multichannel MAC Protocol in Vehicular Ad Hoc Networks[END_REF]. In order to increase the lifetime of a single link, the mobile node should avoid unnecessary handovers [START_REF] Sadiq | An Intelligent Vertical Handover Scheme for Audio and Video Streaming in Heterogeneous Vehicular Networks[END_REF]. Skipping unstable VLC sessions will increase the connection time i.e., t c [s], estimated as

t c = d max + (X P D -X i ) V P D , (9)
Setting the minimum required connection time equal to T c [s], the dominant interferer i will be candidate as AP if IN R i ≥ α and t c ≥ T c . This process will be successively reapplied in case of multiple dominant interferers. Finally, the dominant interferer i with the maximum connection time is returned to adaptive modulation scheme, as explained in Algorithm 2. On top of handover mechanism, an Adaptive Modulation Scheme (ADM) protocol triggers different modulation schemes based on Table I, in order to obtain a proper compromise between the data rate and the communication robustness. The relative bitrates are appropriately selected at the simulator to cope with the testbed limited components. After co-channel interference mitigation among different LEDs using frequency division multiplexing [START_REF] Costanzo | Visible Light Indoor Positioning in a Noise-aware Environment[END_REF], the proper modulation scheme is used for transmission, while the low SN R i amount requires handover decision.

SN R max = N max i=1 (rPr,i) 2 σ 2 n if SN R ≥ α then estimate t c if t c ≥ T c then T c = t c candidate i end if end if return T c , i if null then backoff t s end if

III. VALIDATION

The proposed INR-based protocol for transmission and handover has been implemented in NS3 simulator [START_REF] Aldalbahi | Visible Light Communication Module: An Open Source Extension to the ns3 Network Simulator With Real System Validation[END_REF] as well as a real specific scenario for proper experimental validation. In particular, we compare simulated and experimental protocol performance in terms of • handover rate, namely the number of successful handovers per minute, •• average delivered data per handover, namely the total delivered data over the total number of handover, and • • • handover delay ratio, namely the ratio between time required for performing the handover and total connection time.

In practice, we use an array of commercial low-power LEDs as APs, in order to build a small-scale setup. Each optical transmitter is connected to a dedicated output port of an Arduino 1 in order to control and drive separately each LED. A simple receiver consists of a low-cost PD and other passive elements is used for delivering the data to the signal process unity, mainly composed by a USRP 2920 and a Low Frequency daughterboard. The distance between the ceiling of the APs setup and the plane of the mobile receiver is constant, while the number of LEDs and the distance between them is properly modified to allow performance evaluation for different linear node densities (ρ l ).

As a benchmark, we model the SNR-based handover mechanism with the same parameters, and we compared both mechanisms in simulation and experimental tests. In the SNR-based case, the handover decision is based on the maximum received power and the connection time, while, for the sake of fairness, we applied the ADM on top of both handover mechanisms. Table II shows the main system parameters used in the simulation and experimental tests.

Assuming the testbed scenarios of Fig. 2, different configurations of LEDs represent different node's linear densities while feeding the signal to the USRP at different time interval manipulates different speeds. As the mobile node moves along the coverage area of a given AP, the SNR metric is used to measure the link QoS and allocate the modulation scheme for the next transmission bursts, whereas the interference is mitigated based on the frequency shift keying. When the error probability increases as a result of degradation of the link quality, the INR-based handover nominates an alternative AP with the highest connection time. On the other hand the SNR-based algorithm is suggesting the link with the highest received power that last at least T c seconds, otherwise it backs off the scanning for t s seconds. In Fig. 3, we compare the handover rate (i.e., HO rate [min -1 ]) obtained with the INR-based and the SNRbased handover schemes, for different speeds. It can be shown that the handover rate has a reverse relationship with the average connection time t c [s] as follows:

HO rate = 60 t c , (10)
Regarding the excellent matching between the simulation and the testbed results, it is admitted that, on the same speed as the AP density increases, the SNR-based handover rate increases linearly. According to Eq. ( 9) and ( 10), going from lower to higher ρ l , the distance between node P D (Rx) and node i (Tx) is decreasing. Since d max and V P D are fixed, the average connection time declines accordingly, and as a result the handover rate increases. On the other side, the INRbased handover rate increases in the sparse AP deployment and saturates in the dense AP distribution. At the dense AP installation (X P D -X i ≤ d max ) the longest t c is determined by the farthest AP which stands at most d max meter far from the mobile node. Higher velocity corresponds to the lower connection time, and hence higher handover rate. In Fig. 4, we assess the impact of the node density on the average delivered data per handover at speed 0.3 mps and 0.15 mps. In both medium access schemes, the transmission carried out when the front-end device is well covered by the AP (i.e., SN R ≥ α). Giving t c the available time for transmission, at rapid gait (i.e., 0.3 mps) the connection time is shorter than the one of lower speed (i.e., 0.15 mps) providing bigger portion of time for transmissions in different bitrate according to Table I. As a result, the average delivered data is larger for the higher speed than the lower one. It is also obvious from Fig. 4 that the INR-based mechanism always transmits more data per access point than the SNR-based scenario, that is not because of higher data rate obtained with our proposed scheme (i.e., the transmission scheme is the same at both proposed scheme and the benchmark), but because the INR-based mechanism makes much less handovers for the same AP distribution. As the ρ l increases the average delivered data decreases in both schemes with slight incline for INR-based scheme.

Regarding the handover delay, the time efficiency of a given handover mechanism is defined as the time required for performing the handover(HO delay ) to the whole connection time offered by that handover (t c ) [START_REF] Hien | Handover in outdoor Visible Light Communication system[END_REF] i.e., HO delay ratio = HO delay t c × 100,

HO delay = t HO req + t HO res ,

with t HO req and t HO res defined as the time necessary to • make an handover request and •• provide an handover response, respectively. Assuming the same frame size used for handover request and handover response (i.e., t HO req = t HO res ) included in the MAC frame format specified by IEEE standard for short-range wireless optical communication using VLC 802.15.7 [START_REF]IEEE Standard for Local and Metropolitan Area Networks-Part 15.7: Short-Range Wireless Optical Communication Using Visible Light[END_REF], the reserved bits in the control frame of Fig. 5 is used to specify the velocity of the mobile node, as well as d max of each LED rather than the number of interferers N for the INR-based mechanism. The required time to perform a handover is the handover frame size (i.e., HO f rame size [bits]) over the bitrate at both the transmitter and the receiver (i.e., M ) ). Given the modulation scheme M (M ∈ {OOK, 4P SK, 8P SK, 16P SK}) and R b,tx (M ) = R b,rx (M ) the handover request time is simplified through eq. ( 13),

R b (M ) = 1 R b,tx (M ) + 1 R b,rx 1
t HO req (M ) = 2 HO f rame size R b,tx (M ) (13)
The average handover time is calculated as follows:

t HO req = M (t HO req (M ) • HO rate (M )) HO rate , (14)
where HO rate (M ) is the handover rate on the modulation scheme M . The evolution of the handover delay ratio comparing to AP density for different velocities is depicted in Fig. 6. At higher velocity, the connection time is lower resulting in less efficient handover delay ratio. For a given velocity, at the low density scenarios the SNR-based performs better due to the shorter handover frame size. As the node density increases, the number of SNR-based handover grows linearly. Simultaneously, the dense AP deployment boost the received power at the front-end device which in return applies higher modulation levels and increases the bitrate, although the handover rate progresses faster than the bitrate, leading to HOrate R b ≥ 1. That is why the SNR-based handover delay ratio is monotonically increasing. On the other side, the interference-based handover mechanism substantially outperforms the SNR-based handover by recovering the higher overhead with the the lower HO rate . The ratio HOrate R b in INR-based scheme is less than unity because the HO rate reaches the saturation level at the high density, while the R b continues increasing leading to decreasing delay ratio.

IV. CONCLUSIONS

A novel handover technique for mobile visible light communication nodes based on the INR and IIR is proposed and validated through the simulation and the experimental tests. A fair comparison between INR-based and SNR-based handover mechanisms is presented in terms of handover rate, average delivered data per handover, and handover delay ratio. The results in both simulation and experimental tests confirmed that the proposed algorithm avoids redundant handovers, and consequently, improves the average delivered data per handover, reducing significantly the overhead information of the network. At the same time, the handover delay ratio showed up to 1.5 times enhancement as compared to the SNR-based handover in high density scenarios (i.e., ρ l = 10 m -1 ).

Fig. 1 .

 1 Fig. 1. Proposed VLC network architecture.

10 Fig. 2 .

 102 Fig. 2. Testbed configurations for performance evaluation, considering different ρ l .

Fig. 3 .

 3 Fig. 3. Handover rate versus node density.

Fig. 4 .

 4 Fig. 4. Delivered data per handover versus node density.

Fig. 5 .

 5 Fig. 5. General MAC frame 802.15.7 and the modified control field.

Fig. 6 .

 6 Fig. 6. Handover delay ratio versus node density.

ACKNOWLEDGMENT

This article is based on the work from COST Action NEWFOCUS CA19111, supported by COST (European Cooperation in Science and Technology).


Daniel Prost

Measurement of the Amplitude and Polarization of the Electric Field by Electromagnetic Infrared Thermography

Keywords: Antenna, axial ratio, infrared thermography, near-field measurement, polarization, thin film pattern

The ElectroMagnetic InfraRed (EMIR) method developed and used at ONERA provides a map of the amplitude of the electric field radiated by an antenna: a thin film of low conductivity is placed in front of the EM source, the electric field induces currents which heat the film by Joule effect; this heating is measured by an infrared camera. A low-frequency modulation of the signal associated with synchronous demodulation of the images allows eliminating the parasitic thermal phenomena. By using anisotropic films with conductivity patterns, we propose here for the first time a vector measurement of the electric field. In the case of linear polarization, this technique allows to obtain the field direction, and in the case of circular polarization, to measure the axial ratio. Furthermore, we show, with the help of a patterned array, that we can obtain the spatial structure of the near field, both in amplitude and direction.

I. INTRODUCTION

T HE ElectroMagnetic InfraRed (EMIR) method, which has been developed by the French Aerospace Lab, ONERA, in the 90s, consists in submitting a weakly conductive thin film to an electromagnetic (EM) field. Due to induced currents, the heating of the film, proportional to the absorbed power (and therefore to the square of the field amplitude), can be recorded with an infrared camera. To increase the spatial resolution, and eliminate the effects due to thermal processes other than the local heating by the incident field (convection, conduction), the EM source is modulated at low frequency (a few Hz) while the thermal frames are being demodulated at the same frequency (lock-in thermography) [START_REF] Balageas | EMIR: A photothermal tool for electromagnetic phenomena characterization[END_REF] (see Fig. 1). This technique allows visualizing EM maps, on a large 2-D scale (thin films of 1 m 2 have been used). This method can be used for antenna radiation pattern [START_REF] Prost | Imaging electric and magnetic near field of radiating structures by infrared thermography[END_REF], [START_REF] Prost | EMIR field imaging for diagnosis and characterization of space antennas[END_REF] measurement, or field leakage detection and in general, to rapidly obtain EM field characterization. However, using homogeneous and isotropic thin films only allows getting the (square of) amplitude of the field. This is the reason why anisotropic films have first been tested and patented by ONERA [START_REF] Prost | Dispositif pour révéler des variations spatiales de polarisation d'un rayonnement électromagnétique[END_REF]. Those films are made of dielectric film partially coated with linear strips of a thin conductive layer. Such a grating does not behave like a polarizer since the electrical conductivity is weak enough but it absorbs a part of the energy of the field component parallel to the strips. Therefore, to obtain a thermal measurement, the grating had to be parallel to the field, assuming a constant linear field polarization (plane wave condition typically). Finding the direction of the field implied performing many measurements with different orientations of the grating, as in [START_REF] Levesque | Capteur de champ électromagnétique par thermographie infrarouge[END_REF] and [START_REF] Levesque | Vector characterization of electromagnetic fields by infrared thermography[END_REF].

To obtain the field polarization in a general case (including circular-or elliptic-polarization), we have recently developed new films based on a specific star pattern that can directly give the field orientation (in case of linear polarization) or the axial ratio (in case of circular polarization). Moreover, by using films with such pattern-based networks, a vector mapping of the field can be obtained (especially in the case where the field structure is not that of a plane wave and exhibits a spatial variation). Although we did not measure the phase with this method, there are "plane-to-plane" (PtP) algorithms [START_REF] Bucci | Fast and accurate near-field-far-field transformation by sampling interpolation of planepolar measurements[END_REF], [START_REF] Will | Phase measurements of electromagnetic fields using infrared imaging techniques and microwave holography[END_REF] that can give it from two amplitude measurements at different distances from the source: the thermal measurements are processed by the PtP iterative Fourier technique to determine the phase and then the complete complex far-field pattern of an antenna for example [START_REF] Yaccarino | Phaseless bi-polar near-field measurements using a two-plane squared amplitude interpolation algorithm[END_REF]. However, these algorithms have not been tested in the case of anisotropic films. This interesting question could be addressed in future work.

II. ANISOTROPIC FILMS

A. Shape and Simulations

Beyond the grating structure, we considered that a way to discriminate the direction of the electric field was to use specific patterns of lines with different directions. Simple crosses such as two lines crossing at 90 • do not allow direction discrimination, since the quadratic behavior of the electric field gives two possible directions for the same components on an orthonormal basis. Therefore we focused on "star" patterns with 8 to 16 branches (see Fig. 2). The goal of such shapes is obvious: only the branches parallel to a component of the field will heat, and the heating will be proportional to the square of the amplitude of the in-plane electric field. Resonant behavior is not expected thanks to the very weak conductivity of the pattern (with a resistance per square of almost 500 ). To optimize the shape and size pattern, we did some simulations (using HFSS full wave simulation software). Fig. 2 shows three cases for the surface loss density in the branches.

1) If the branches are too thick, the contrast between parallel (H) and perpendicular branches (V) to the field is not sufficient. 2) Using a too conductive pattern exhibits a resonant behavior that must be avoided (loss and heating, will be locally much too high). 3) Optimizing the surface conductivity and branch thickness provides the expected behavior: weak loss in the branches perpendicular to the field, maximum heating in the branches parallel to the field, and the same value divided by almost 2 on the 45 • branches. These simulations allowed us to converge toward an optimized design: stars with 8 or 12 branches of 200 μm width, a diameter of about 4 cm and a surface resistance of a few hundred Ohm. These rules also guided us for the etching of the samples.

B. Pattern Realization and Characteristics

The patterns are obtained by laser printing of 50-μm wide Kapton samples. The laser beam burns the Kapton locally and the resulting carbonized strips correspond to the desired pattern. Suitable laser beam parameters (power, moving celerity) are used to control the surface impedance; several passes of the beam are necessary to obtain the correct linewidth. Figs. 3 and4 show, for example, a 12-branch star with a diameter of 5 cm. The regularity of the pattern branches is checked using a profile projector: the width is found to vary in the 200-215 μm range, which is within the required accuracy for our application.

III. PATTERN AND NETWORK STRUCTURE CHARACTERIZATION

A. Star Pattern Performance and Management

In order to characterize the performance of the star patterns, we consider the measurement of the near field of a linearly polarized horn antenna. In this application, the star patterned film is located in front of the opening of the horn antenna. This horn antenna [START_REF] Albert | Design, manufacturing and test of a spline-profile square horn for focal array applications[END_REF] has a square aperture and provides linear electric field polarization in the Ku-band (12-18 GHz). We voluntarily place the pattern in a "random" orientation, such that none of the branches is parallel to any side of the horn aperture (actually the orientation of the pattern does not matter, as will be shown later).

The infrared demodulated image is given in Fig. 5. The maximum heating is obtained at the center of the star. The corresponding value is high enough to allow a quantitative measurement of the electric field amplitude. Indeed, the heating is related to the tangential electric field. Through a thermal analysis, we can compute the absorbed surface power density, in (W/m 2 ) [1]

P abs = T (4h 2 + (Cρeω) 2 ) ( 1
)
where h ≈ 10 W/K/m 2 is the heat transfer coefficient; ρ = 1420 kg/m 3 is the film density; C = 1090 J/kg/K is the film heat capacity; e = 50 μm is the film thickness; ω = 0.5 rad/s is the modulation low frequency. The absorbed surface power density is also given (peak value) by

P abs = E 2 /Z s .
(

) 2
Combining the two previous equations gives

E = k √ T . ( 3
)
With our film (with surface impedance Z s = 200 ) we obtain k ∼ 210. N.B. This factor k can alternatively be calibrated using a reference measurement.

Moreover, as it is located in the center area of the horn aperture, we expect the pattern to give the theoretical field with maximal level. We repeated this measurement for different incident powers, from 0.5 to 5 W and compared values with HFSS simulation. Fig. 6 shows that the amplitude is correctly obtained from the heating value.

To obtain the field direction, we integrate the measured heating on each branch of the star. The integration lines are detected by matching the star pattern to the infrared image: the optimized location corresponds to the highest value of the sum of the 12 integrals. Fig. 7 represents the heating of the branches and the corresponding integration lines (the average heating T is given for each branch).

Taking the square roots of those heatings gives 12 components with associated angles of the electric field (Fig. 8).

One way to process those components is to consider three orthogonal bases (corresponding to the three red, blue, and green, colors in Fig. 8). On each basis, there are two possible resulting vectors for the field, since one component only gives the orientation but not the direction. For instance, combining 1 and 4 (or 7 and 10) blue components from Fig. 8 gives V2 direction in Fig. 9, while combining 1 and 10 (or 4 and 7) gives V1. The same process with the green and red basis provides also two possible vectors V3 and V4 (for the green) and V5 and V6 (for the red). Finally, we obtain two possibilities for the field direction, in each of the three basis. Only one orientation is compatible with the three orthonormal bases, and corresponds to the field direction. Here only V1 (from blue basis), V4 (from green), and V5 (from red) are compatible. The mean value V1, V4, V5 gives the field direction (-0.44 • versus the horizontal axis in this case).

To check the accuracy of the detection, we duplicated this measurement with another sample with only eight branches. Moreover, we achieved two more measurements, the second one after having slightly turned the film (by almost 1 • ).

We obtained on Fig. 10 an angle equal to -0.67 • (respectively -0.47 • ) with the eight-branch star pattern (respectively after moving it slightly and rotating it by 1 • ). As we measured -0.44 • with the 12-branch star pattern, it appears that this method gives the field orientation with an accuracy better than 0.3 • .

B. Star-Pattern-Based Network Performance

In this paragraph, we present another kind of film, with an array of several stars. The purpose of this array is to allow a spatial resolution of the field structure, for cases where this structure is not homogeneous (particularly when we are out of the plane wave regime). Each pattern or element of the network will provide a local measurement of the electric field. Having such a network with many stars is a topology problem, since the star must be as large as possible, but the distance between two centers is as small as possible since it determines the spatial resolution. It appears that three-branch stars permit a quite compact network: for example in Fig. 11 the elements are distant of 1 cm from each other while the branches have a length of 7 mm.

The way each three-branch pattern element gives the local field direction is slightly different from the eight and 12-branch patterns presented before. Indeed, we may still define three coordinate systems, but they are not orthogonal anymore. Therefore, the measured heating on a given branch cannot directly give a field component. Let us define e 1 , e 1 , and e 1 , the three directions of the three branches of the star. The heating on branch e 1 (respectively e 2 ) for instance, corresponds to the square of the scalar product

x 1 = E| e 1 (respectively x 2 = E| e 2 ). (4)
However, here we are on a non-orthogonal basis ( e 1 , e 2 ). The Gram-Schmidt orthogonalization process [START_REF] Wong | An application of orthogonalization process to the theory of least squares[END_REF] is required to obtain the components on a (local) orthogonal basis; we obtain

x = x 1 x 2 = G 12 E ( e 1 , e 2 ) with G 12 e 1 | e 1 e 2 | e 1 e 1 | e 2 e 2 | e 2 (5)
where G 12 is the Gram-Schmidt matrix (two branches are 120 • angle-shifted)

G 12 = 1 -1/2 -1/2 1 (6) E ( e 1 , e 2 ) = G -1 12 x. ( 7
)
The last step consists of a rotation to come back in the global ( i , j) reference coordinate system, which corresponds to (using the usual R rotation matrix)

E ( i, j ) = R 12 E ( e 1 , e 2 ) . ( 8
)
This finally gives E x and E y components in the reference system but uncertainty still remains since two field orientations are possible The same process is achieved on the two other bases ( e 2 , e 3 ) and ( e 3 , e 1 ), with corresponding G 23 and G 31 Gram-Schmidt and R 23 and R 31 rotation matrices. Having those three solutions in the three coordinate systems makes it possible to remove the uncertainty, as in Section III-A. Fig. 12(b) illustrates for instance what can be obtained for the horn antenna of Section III-A. In order to improve the spatial resolution, we collected several measurements at different positions of the film (half a step size displacements, in H than in V direction): Fig. 12(c) and(d).

E ( i, j) = E x i + E y jor E = E x i -E y j. ( 9
)
To illustrate how this array can resolve a complex electric field structure, at a much lower λ scale, we also tested it on a high impedance surface (HIS) metamaterial structure. It consists of a microstrip line with a special structure of mushrooms-type resonant elements localized inside the substrate, between the ground plane and the line [START_REF] Sievenpiper | High-impedance electromagnetic surfaces with a forbidden frequency band[END_REF]. This kind of structure is developed in order to reduce the size of radiating elements, in particular in the GNSS context (where Con- trolled Radiation Pattern Antennas-CRPAs-are often used), by limiting the coupling between the elements. To develop and study such structures, the measurement of the near electric field is necessary. Due to its small dimensions, a high spatial resolution is very hard to obtain with a standard probe measurement technique. Infrared thermography is therefore particularly interesting in this case. To this extent, in [START_REF] Prost | Electric field imaging of a high impedance surface for GNSS array decoupling application[END_REF], infrared frames of the electric field amplitude had been obtained using the "traditional" EMIR method. With a star-patterned network film, we can now exhibit both the amplitude and the local direction of the electric field, with accurate details at a very low scale. As seen in Fig. 13, this direction changes with frequency. In particular, resonant and passing frequencies correspond to very different electric field structures.

C. Evaluation of Axial Ratio for Circular Polarization

In a circular polarized field, each branch of the star is expected to heat equivalently. However, in the very common case of imperfect circular polarization, it is of interest to estimate the axial ratio (the ratio between the major and minor axes), usually given by

A R = 20Log(E x /E y ). ( 10
)
Let us assume an elliptical polarization of the electric field. The patterned film is located in the wave plane (perpendicular to the wave vector). In a ( i, j ) reference coordinate system corresponding to the orientation α of the ellipse, the electric field can be written as

E ( i, j) = E x (t) i + E y (t) j ( 11
)
with E x (t) = E x .si n(ωt) and E y (t) = E y .cos(ωt).

In the same coordinate system, the direction vector of a given branch θ i of the star is written as (see Fig. 14)

v θ i = cos(θ i ) i + sin(θ i ) j. ( 12
)
The squared scalar product gives the expected heating of the branch

T θ i = E| -→ v θ i 2 = (E x (t) i + E y (t) j| cos(θ i ) i + sin(θ i ) j) 2 . ( 13
)
Then, after integration over time, we find

T θ i = E 2 x cos 2 (θ i ) + E 2 y sin 2 (θ i ). ( 14
)
The EMIR measurement thereby provides measurements of the heating on each branch T m (θ i ). We thus must compare these measurements with the computed heating figures of ( 14), which depend on the three unknowns (α, E x , E y ). Placed in the actual elliptical polarization coordinate system, the following error function E(α, E x , E y ) must be minimum (n being the number of branches):

E α, E x , E y = n i=1 [T m (θ i )-T(θ i )]. ( 15
)
The n angles θ i are known to be within the angular offset α, and the T(θ i ) depend also on E x and E y .

After having minimized this three-variables function (for instance using the MATLAB "fminsearch" function or equivalent), we obtain ellipse parameters E x , E y , and α, (and consequently the orientation of the ellipse and its axial ratio).

Fig. 15 shows the 12-branch star patterned film located in front of a circularly polarized Ka-band horn antenna. The camera (not visible in the picture) is the film. Note that the horn-to-film distance can be reduced to the far/near-field boundary. This could be useful in the case where we measure the AR as a function of angle: indeed, we only present here boresight measurements of the AR, but for many applications it is important to know how the AR depends on the angle (to obtain the AR Beamwidth). The best solution is, therefore, after having placed the film and the camera precisely in the antenna boresight angle, to rotate the antenna itself (this is much easier than moving the film and/or IR camera). Thus, the thin film and the camera sensor remain strictly parallel, and we avoid any optical distortion of the star pattern. The closer the film is, the wider the range of angles that can be explored.

The measurement at 19 GHz (middle of the band) is given in Fig. 16:

The minimization process applied to EM mapping (obtained from Fig. 16) gives the axial ratio AR = 0.19 dB (Fig. 17).

Another measurement at the very beginning of the frequency range of the antenna ( f = 17.3 GHz) illustrates the degra- dation of the axial ratio that reaches 0.45 dB. In this case, we repeated the measurement with different powers (and a much lower horn-to-film distance). As expected (Fig. 18), the radiated power only increases the E x and E y values but it has obviously no impact on the orientation and axial ratio. This also proves the robustness of the measurement.

To conclude, we compared the axial ratios obtained using this original EMIR method with traditional RF measurements, over the whole antenna band (17.5-20.5 GHz). Compared to those traditional RF techniques, the EMIR technique offers an alternative process: Fig. 19 shows a good consistency between the two methods, although an uncertainty of about 0.2% persists on EMIR measurements. This uncertainty can be attributed to several factors: 1) the limited number of camera pixels corresponding to a branch of the star pattern, which also depends on the relative position of the branch with respect to the camera frame and thus varies between branches (this problem could be solved using a higher camera resolution); 2) the widths of the branches are not strictly equal (see Section II-B), leading to under or overestimation of the heating; and 3) the tolerances of the positioning system, although extremely small (we have taken care on this), are not zero. However, the usual specification that the axial ratio should be less than 0.5 dB in the band can be validated by the EMIR method.

IV. CONCLUSION

Infrared thermography is now a well-proven method to obtain maps of the electric (or magnetic) field amplitude. We have shown here that it can also give the polarization of the field, i.e., its direction (for linear polarization) or its axial ratio (for circular polarization). This is achieved by placing in the wave plane an insulating thin film on which weakly conductive patterns are etched. The advantage of this method is that it is relatively easy to implement, as it requires neither a very rigorous alignment nor a perfectly controlled displacement of the sensor with respect to the radiating source. Moreover, in the case of a specific near-field structure, a network of these patterns can also reveal the spatial distribution of the field, even with details much smaller than the wavelength.

Fig. 1 .

 1 Fig. 1. EMIR method.

Fig. 2 .

 2 Fig. 2. |E| on star pattern and profiles along branches: (a) too wide lines leading to transverse currents, (b) too high conductivity leading to resonant behavior, and (c) optimized design.

Fig. 3 .

 3 Fig. 3. Star pattern, zoomed-in view on 200-μm wide branch using a profile projector.

Fig. 4 .

 4 Fig. 4. Horn antenna and 12-branch patterned Kapton film.

Fig. 5 .

 5 Fig. 5. Thermal image of the 12-branch star located near the horn antenna aperture.

Fig. 6 .

 6 Fig. 6. Field components obtained from the thermal measurement of the 12-branch star.

Fig. 7 .

 7 Fig. 7. E-field from IR measurement and simulation, for P i = 0.5-5 W.

Fig. 8 .

 8 Fig. 8. Integration lines and heating along the optimized 12-branch star.

Fig. 9 .

 9 Fig. 9. Possible directions of the field in the three bases, and resulting horizontal direction detected ( f = 15 GHz).

Fig. 10 .

 10 Fig. 10. (a) E-field direction given by an eight-branch star pattern and (b) second measurement on the right after applying a 1 • offset angle on the film.

Fig. 11 .

 11 Fig. 11. Compact network of three-branch star patterns, thermal image, and three heating corresponding to one of the patterns.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Horn antenna near field: (a) amplitude "traditional" EMIR measurement [2], (b) amplitude and direction (array of three-branch stars), (c) same with two measurements collected for two films positions, and (d) three positions.

Fig. 14 .

 14 Fig. 14. Offset angle α and branch angles θ i on a global basis.

Fig. 15

 15 Fig. 15. 12-branch star-patterned film for circular polarization measurement.

Fig. 16 .

 16 Fig. 16. Thermal (scaled and root squared) image of a 12-branch star pattern.

Fig. 17 .

 17 Fig. 17. Ellipse obtained from Fig. 16 measurement ( f = 19 GHz).

Fig. 18 .

 18 Fig. 18. (a) Ellipses and axial ratios for different powers P i . (b) Amplitudes of E x and E y proportional to √ P i . ( f = 17.3 GHz).

Fig. 19 .

 19 Fig. 19. Axial ratio as a function of frequency: comparison of EMIR and usual RF methods.

ACKNOWLEDGMENT

The author would like to thank Thales Alenia Space (TAS) for making available the horn antennas and providing measurement facilities. He was grateful to the CIRIMAT Lab, Center Inter-universitaire de Recherche et d'Ingénierie des Matériaux, UMR CNRS 5085, for providing the printed films that were used in this work.


A Brioude

Frédéric Lequevre

Jacques Mugnier

Charles Bovier

J Dumas

Jean Claude Plénet

C B O V I E R

J C Plenet

Densification and nanocrystallisation of sol-gel ZrO2 thin films studied by surface plasmon polariton-assisted Raman spectroscopy

Introduction

The Surface Plasmon Polariton (SP p ) excited by evanescent waves on smooth metal surfaces and out-coupled in attenuated total reflection (ATR) configuration improves the limited sensitivity of backscattering Raman investigations when studying ultra thin films. Thus, SP p -assisted Raman spectroscopy has been used to study many systems as molecular overlayers adsorbed on well defined metal surfaces [1][2][3][4]. The enhancement of the Raman intensity [5][6][7] is due to the amplification of the electromagnetic field at the metal-dielectric interface due to the SP p electromagnetic resonance.

This enhancement can be obtained in the so-called Kretschmann geometry [8] which consists of a prism and a thin metal film covered by the dielectric which is to be studied. This set-up allows excitation of the metalthin film interface SP p using laser incident light coming through the prism. This resonant coupling occurs at a specific incident angle (θ AT R ). The exact value of θ AT R depends on a relation between the wavelength of the incident light, the dielectric constants and thicknesses of each layers making up the samples [9]. a e-mail: plenet@dpm.univ-lyon1.fr b CNRS UMR 5620

In the first part, we describe sample preparation constituted by silver layer deposition and ZrO 2 very thin film preparation by the Sol-Gel (SG) process. Densification of ZrO 2 ultra-thin sol-gel films is studied, measuring the evolution of refractive index and thickness of the dielectric film after each annealing step. In the last part, we present Raman scattering spectroscopy results when studying two samples of different thicknesses. The estimated enhancement factor of the Raman scattering intensity (around 85) is due to the SP p excitation and is sufficient to determine the ZrO 2 films phases. Transmission Electron Microscopy (TEM) observations are conducted on the same Zirconia layers in order to check and complete the optical and spectroscopic results respectively.

Samples preparation and heat treatment procedure

Thin silver films were electron-gun deposited onto a precleaned glass ( Pyrex, optical index 1.472 at 514.5n m ) under vacuum of 10 -6 torr. The deposited layer thickness (50 ± 0.5 nm) and the evaporation rate (0.5n m /s) were measured by a quartz balance. This thickness was determined in order to improve the light coupling efficiency. The silver growth rate was determine to procedure a smooth metal surface [10,11]. Very thin ZrO 2 films have been prepared by SG route and deposited by a dip-coating method [START_REF] Brinker | Sol-Gel Science[END_REF][START_REF] Klein | Sol-Gel technology for thin films, fibers preforms, electronics and speciality shapes[END_REF] onto the silver layer. The sol was obtained by mixing zirconium n-propoxide [Zr(OC 3 H 7 ) 4 , 70% solution in propanol (Fluka)] [START_REF] Tohge | [END_REF] with a stabilizer to prevent ZrO 2 precipitation. In this study, zirconium alkoxide was modified with acetylacetone (acacH) according to the exothermic reaction [15]:

Zr(n -OC 3 H 7 ) 4 + acacH → Zr(n -OC 3 H 7 ) 3 +(acacH) + C 3 H 7 OH.
When the hydrolysis and condensation reactions have been performed, 2-propanol (Prolabo) was added to adjust the sol viscosity at the wanted value. This last step during the process was necessary to adjust the ultra thin films thickness in the range around few nanometers [16]. After the dip coating deposition, the films were dried at a temperature of 100 o C. Such films are continuous and crack free as proved by micro-optical observations. In this paper, we present results on three samples labelled A, B and C. As it has been presented above, film thicknesses are adjusted by addition of 2-propanol to the starting SG solution. The sample A and B has been obtained by diluting the sol with a volume ratio 6/1 for 2-propanol and sol respectively. In order to elaborate the thicker sample C, a 4/1 ratio has been selected. The sample A was annealed under Infra-Red (IR) lamp at a temperature around 750 o C under a pur O 2 flow. The total heat treatment duration was 28 minutes by step of 4 min mainly. This annealing time does not take into account the temperature rise time of about 10 seconds and the room temperature decrease of 1 minute. On the one hand, the aim of this annealing process by IR irradiation is to induce the nanocrystal growth into the film without destroying the silver layer and on the other hand to eliminate carbon residues due to SG method. The samples B and C were annealed for 15 minutes under IR lamp at a temperature around 500 o C under continuous pure O 2 flow and then were heat-treated in a tubular furnace under oxygen flow at 700 o C for 5 minutes. The aim of this last annealing treatment is to achieve crystallisation as it has been already shown elsewhere [17].

Thickness and refractive index measurement

A just dried SG film exhibits an amorphous structure [START_REF] Brinker | Sol-Gel Science[END_REF][START_REF] Klein | Sol-Gel technology for thin films, fibers preforms, electronics and speciality shapes[END_REF]. Therefore, the aim of this study is to determine the crystallisation temperature for a well defined annealing procedure. Thicknesses and refractive index evolution of the sample A were measured using an SPR optical set-up already described in a previous paper [16]. These measurements presented in Figure 1 have been performed at a wavelength of 594.5 nm (He-Ne laser). the same time sample A thickness continuously decreases.

A break in the slope of the two curves is observed, suggesting the onset of a phase transformation (see Fig. 1). After 8 min, the optical index of the sample A roughly remains constant, and its thickness just weakly decreases. After a complete annealing treatment as described above the films thicknesses were 12 ± 0.5n ma n d2 6 .6 ± 0.5 nm for the samples B and C respectively. Their refractive indices were 2.06 ± 0.01 and 1.92 ± 0.01 respectively confirming that the ultra-thin film refractive index is strongly dependant on the films thickness. Differences in their measured refractive indices have already been explained [18] and confirmed that the thinner the SG film is, the denser it is.

Transmission Electron Microscopy studies

In order to compare ZrO 2 ultra-thin films behaviour, and effects of annealing procedure, and thickness samples A, B, C have been analysed by electron microscopy. Conventional Transmission Electron Microscopy (CTEM) studies have been performed using a TOPCON EM-002B working at 200 kV to obtain the general picture on the changes in crystalline-phases. In this case ZrO 2 very thin films have been peeled off from their substrates and deposited onto a 3 mm copper grid for direct observation. Sample A CTEM analysis result is presented in Figure 2a. The associated diffraction pattern (Fig. 2b) shows that sample A exhibits a tetragonal phase (space groupe P4 2 /nmc, [19]). It can be observed that the rings are not well-defined and not continuous. This means that the sample contains isolated nanocrystals so that only few crystalline orientations are detected in the diffraction pattern. The large diffuse and intense rings are the signature of the predominant amorphous phase. Due to the not well-resolved nanocrystals shape in the CTEM micrography, it is no possible to determine the mean diameter of each nanocrystal. Each dark zone contains several particles (see Fig. 2a). So an HRTEM study of the same area must be conducted and is presented in Figure 3a. This picture shows several crystals randomly oriented with sizes in the range between 5 to 10 nm. The corresponding Fourier transformation of this image is shown in Figure 3b with indexed points corresponding to (101) atomic planes. Taking into account that the film A is a mixture of an amorphous phase and of nanocristals, the difference between the measured refractive index and bulk one (2.05 and 2.20) [20,21] respectively can be explained.

As films exhibit different refractive indices for the same heat treatment, their structures have been investigated by CTEM. CTEM analysis results are presented in Figures 4a and5a respectively. The diffraction pattern associated with Figure 4a (sample (B)) shows clearly a mixture of characteristic rings for the tetragonal phase and the monoclinic one (space groupe P 2 1 /c). For sample (C) (Fig. 5b) the tetragonal phase is only observed. It is to notice that for thicker films used as planar optical waveguide and exhibiting a thickness in the range around 150 nm, monoclinic phase only appears at temperature annealing greater than 600 o C and for an annealing time of 1 hour [22,23]. For sample B, annealing, temperature and duration required for the crystallisation are lower than for thicker films, thus the thermal behaviour of the SG films depends on their thicknesses. We have already observed this phenomenon for TiO 2 SG films [17].

Micro-Raman spectroscopy experiment

The device used to perform the Raman experiment is presented in Figure 6. The Raman scattering measurements were obtained with an XY Dilor triple spectrometer followed by a Nitrogen cooled CCD (Charge Coupled Device) multichannel detector. The 2016 type Spectra-Physics Argon laser beam (λ = 514.5 nm, operating at 40 mW) is incident on a hemispherical prism at the plasmon resonance angle θ AT R . The samples are positioned on the top of the prism on the glass side with an index liquid (glycerol) of refractive index (n =1.474 at λ = 514.5 nm) very close to one. Details of the experimental device have been presented and discussed elsewhere [17]. An horizontal micrometric shift of a planar mirror (M1) controls vertically the incidence of laser beam on the parabolic mirror (M2), so that it focuses precisely the light at the silver/ZrO 2 interface through the prism (P) and along the optical microscope axis. The light is then collected with the microscope objective of the Raman microspectrometer.

The main advantage [17] of micro-Raman spectroscopy experiment for such thin layer studies is that it is a non destructive analysis method, in comparison to TEM observations. The results of the micro-Raman spectroscopy experiments are shown in Figure 7. For the thinner film (sample B), no Raman bands can be distinguished from the background noise. For the thicker film (sample C), the spectrum exhibits wide bands at around 240, 459 and 640 cm -1 . The observed bands are assigned to the tetragonal ZrO 2 phase [22,24]. The broadness of the bands is due to the very small size of the particles, a few nanometers. Even though the Raman signal enhancement factor is about 85, the thickness (e.g. the low quantity of active and probably the very weak Raman scattering cross section of ZrO 2 phase does not allow to obtain a significantly good signal for sample B. The fact, that the bands assigned to the monoclinic phase have not been observed, may have two reasons: we can assume that Raman scattering cross section of monoclinic phase is lower than tetragonal one. In this case, even though the proportion of monoclinic and tetragonal phase has not been determined by microscopy diffraction pattern study on which rings of each crystalline phase are clearly well-defined without noticeable intensity differences. We can conclude that the number of Raman active sites is lower than the one present in a film of the same thickness but only composed of tetragonal nanocrystallites. Raman scattering excited by surface plasmons seems to be a limited spectroscopic technic and has to be completed by microscopy studies. For materials as ZrO 2 , thicknesses below 20 nanometers seems to be difficult to investigate only by Raman spectroscopy. Secondly, if we assume that the Raman scattering cross section of the monoclinic phase is at least equal to tetragonal one, it means that the proportion of monoclinic nanocrystals is very low. In this case, Raman scattering study can pro-

Conclusion

This paper has illustrated an example of preparation of dense ultra thin films of ZrO 2 using the sol-gel route. The precise control of each annealing process step allows to obtain higher refractive indices than in thicker film case. The layer densification and nanocrystallization is realized by IR lamp heat treatment whereas the whole crystallisation is achieved by classical heating. CTEM and HRTEM studies coupled to micro-Raman spectroscopy enhanced by SP p excitation, is proposed as an interesting complementary structural analysis technique for ultra thin films.

Fig. 1 .

 1 Fig. 1. Measured refractive index of sample A and corresponding thickness versus treatment annealing time.

Fig. 2 .

 2 Fig. 2. (a) Conventional Transmission Electron Microscopy image of sample A. (b) Associated diffraction pattern.

Fig. 3 .

 3 Fig. 3. (a) High Resolution Transmission Electron Microscopy image of sample A. (b) Associated diffraction pattern.

Fig. 4 .Fig. 5 .

 45 Fig. 4. (a) Conventional Transmission Electron Microscopy image of sample B. (b) Associated diffraction pattern.

Fig. 6 .

 6 Fig. 6. Optical device used for the micro-Raman experiment.

Fig. 7 .

 7 Fig. 7. Raman spectra of sample B and C.