Xavier Hadoux

Nathalie Gorretta

Jean-Michel Roger

Ryad Bendoula

Gilles Rabatel

Comparison of the efficacy of spectral pre-treatments for wheat and weed discrimination in outdoor conditions

Keywords: Pre-treatment, Discrimination, Supervised Classification, Reflectance Spectrum, Chemometrics, Hyperspectral

The optimal processing of spectral data often requires specific pre-treatments. In the context of spectral discrimination, results can be greatly improved using the relevant pre-treatment. Most importantly, the pre-treatment must be suited to the nuisance variability that has to be removed. This study focuses on discrimination of weed and wheat using spectra acquired in outdoor conditions with uncontrolled lighting and leaf orientations. Both vegetation spectra are highly similar, and due to the context of acquisition, a lot of spectral variability is present. This nuisance variability is modeled using an additive, multiplicative and noise term, each of which affects the measured spectra. Several pre-treatments were therefore evaluated according to their potential to deal with this variability and their effects were described in the feature space. Finally, results obtained with these pre-treatments combined with two discrimination methods (PLS-LDA and Gaussian SVM) are compared and discussed. Results showed that, thanks to their ability to remove nuisance variability, most pre-treatments are effective in terms of classification accuracy. Gaussian SVM classification results are less influenced by pre-treatments than those of PLS-LDA, since the former compensates the pre-treatment effect by using a different non-linear kernel. For this data-set, the best discrimination result was obtained using the combination logarithm and PLS-LDA. Logarithm actually transforms the multiplicative effect into and additive one, which is then effectively dealt with by PLS-LDA.

Introduction

In the field of spectral data discrimination, vectors to be classified lie in a high-dimensional space which leads to poor classification results when using standard statistical classifiers. In fact, [START_REF] Hughes | On the mean accuracy of statistical pattern recognizers[END_REF] has even demonstrated that as space dimensionality increases there is a corresponding drop in statistical mean accuracy. The problem is even more complex because spectral variables are also highly correlated, which implies that covariance matrices of spectra are rank deficient and are thus not invertible. As a result, usual discrimination methods such as Fisher's Linear Discriminant Analysis (LDA) cannot be applied to spectral data.

To handle these problems, in the chemometrics field, Partial Least Squares (PLS) has been introduced by Herman Wold [START_REF] Wold | Estimation of Principal Components and Related Models by Iterative Least squares[END_REF] in order to solve an over-determined regression problem through dimension reduction. In the case of discrimination, the same idea of applying dimension reduction as a first step has been used in Partial Least Squares Linear Discriminant Analysis (PLS-LDA) [START_REF] Barker | Partial least squares for discrimination[END_REF]. PLS-LDA thus tackles both problems, i.e. high-dimensional vectors and highly correlated vectors, and has led to successful discrimination results in spectroscopy [START_REF] Liu | PLS and dimension reduction for classification[END_REF].

In the pattern recognition field, discrimination problems generally relate to high-dimensional data. For discrimination purposes, Support Vector Machine (SVM) has been introduced in [START_REF] Vapnik | Statistical Learning Theory[END_REF] and has proven its efficiency in finding robust class separators even with highly correlated spectral variables, e.g., in remote sensing applications [START_REF] Plaza | Recent advances in techniques for hyperspectral image processing[END_REF][START_REF] Fauvel | A spatial -spectral kernel-based approach for the classification of remote-sensing images[END_REF][START_REF] Rajan | An Active Learning Approach to Hyperspectral Data Classification[END_REF]. The success of SVM also comes from its ability to deal with non-linearly separable classes by using a kernel that maps the data into a higher dimensional space [START_REF] Aizerman | Theoretical Foundations of the Potential Function Method in Pattern Recognition Learning[END_REF][START_REF] Karimi | Application of support vector machine technology for weed and nitrogen stress detection in corn[END_REF].

The use of spectral properties to discriminate plant species has been reviewed in [START_REF] Zwiggelaar | A review of spectral properties of plants and their potential use for crop/weed discrimination in row-crops[END_REF], since then most studies have revealed that spectral indices only lead to an accurate discrimination when weeds and crop have enough spectral differences [START_REF] Piron | Selection of the most efficient wavelength bands for discriminating weeds from crop[END_REF][START_REF] Pedersen | Weed density estimation from digital images in spring barley[END_REF][START_REF] Thorp | A Review on Remote Sensing of Weeds in Agriculture[END_REF][START_REF] Sui | Ground-based sensing system for weed mapping in cotton[END_REF]. In other cases, finding weeds within crop rows is a much more difficult problem to address, and thus requires a more advanced use of the spectral information [START_REF] Brown | Site-specific weed management: sensing requirements -What do we need to see?[END_REF][START_REF] Slaughter | Autonomous robotic weed control systems: A review[END_REF][START_REF] Hadoux | Weeds-wheat discrimination using hyperspectral imagery[END_REF].

The discrimination methods previously described could be directly applied to the spectra. However, these spectra are often affected by nuisance variability due to the environment which has to be dealt with while at the same time keeping the spectral variability that comes from the objects' chemical composition. For example, accurate spectral discrimination has been obtained under controlled leaf orientations in [START_REF] Franz | The use of local spectral properties of leaves as an aid for identifying weed seedlings in digital images[END_REF], but authors have shown that a major drop in discrimination accuracy occurs as soon as leaf orientation changes. In [START_REF] Vigneau | Potential of field hyperspectral imaging as a non destructive method to assess leaf nitrogen content in Wheat[END_REF], authors have explored the effect of leaf orientation in order to build a prediction model of nitrogen content using hyperspectral imaging. They modeled the effect of outdoor lighting and orientation using a linear model and overcame this effect using standard normal variate (SNV) pre-treatment. Recent studies have evaluated the utility of several pre-treatments primarily for regression purposes. For example, Buddenbaum and Steffens (2012a), Buddenbaum and Steffens (2012b) and [START_REF] Stevens | Measuring soil organic carbon in croplands at regional scale using airborne imaging spectroscopy[END_REF] have reviewed different pre-treatments for soil component estimation. Recent studies [START_REF] Jones | Preprocessing strategies to improve MCR analyses of hyperspectral images[END_REF][START_REF] Vidal | Pre-processing of hyperspectral images. Essential steps before image analysis[END_REF] have also evaluated the efficacy of several pre-treatments, under controlled lighting conditions, for chemical resolution applications.

To the best of our knowledge, the effects of pre-treatments have not been studied in a spectral discrimination context, and especially while taking into account the physical meaning of the nuisance variability that occurs in outdoor conditions with uncontrolled lighting and leaf orientations.

In this paper, we aim to answer the following question: based on the acquired spectral data, which combination of pre-treatments and discrimination method leads to the most accurate and reliable results? To do so, in section 2 we model the nuisance variability due to the outdoor environment and optical measurement variations. In this section, we also present several pre-treatments that are adapted to deal with this variability, their scope and their effect in the feature space. In section 3 our experimental setup and experiment design are detailed. Finally, in section 4, our experimental results obtained with PLS-LDA and Gaussian SVM are presented and discussed. Results obtained with PLS-LDA are also detailed in terms of both the class separability and the bias occurring between training and test set.

Theoretical aspects

Notations

In this paper, the following symbolic conventions are used. Boldface capital characters are used for matrices, e.g. X ; small bold characters are used for column vectors, e.g. x; row vectors are denoted by the transpose notation, e.g. x T ; non bold characters are used for scalars, e.g. matrix elements x ij , vector element x i or indices i. The identity matrix of R n is denoted by I n .

When dealing with spectral data, two representations can be used:

• the individual space (in R): a spectrum is represented as a curve depending on the wavelength: r(λ) . This space should be only used for visualization purposes.

• the feature space (in R p ) in which a spectrum corresponds to a vector r in a p-dimensional space where p is the number of wavelengths. Since multi-variate discrimination methods act in this space, any modification on the spectrum should be understood in this space.

In our study we therefore describe the pre-treatments effect in the feature space and we represent their effects in the individual space (for practical representation reasons).

Optical model

Let us consider a reflectance hyperspectral image which was obtained by using a reference surface for lighting measurement and correction (see [START_REF] Vigneau | Potential of field hyperspectral imaging as a non destructive method to assess leaf nitrogen content in Wheat[END_REF]). Due to the uncontrolled lighting conditions and leaf orientations when acquiring the hyperspectral image, the measured reflectance spectrum for a given leaf can vary a lot. This measured reflectance spectrum r(λ) meas has been represented in [START_REF] Vigneau | Potential of field hyperspectral imaging as a non destructive method to assess leaf nitrogen content in Wheat[END_REF] as a linear function of the leaf reflectance spectrum ρ(λ):

r meas (λ) = m • ρ(λ) (1)
where m is a multiplicative term that models the effect due to the uncontrolled leaf orientation toward the incident light.

In [START_REF] Vigneau | Potential of field hyperspectral imaging as a non destructive method to assess leaf nitrogen content in Wheat[END_REF], the additive term only corresponds to a specular reflection from the sun directly toward the sensor. In our study, another additive term that corresponds to spectrometers common defects is also added to equation 1. It is modeled as a linear trend [START_REF] Barnes | Standard normal variate transformation and detrending of near-infrared diffuse reflectance spectra[END_REF] which is scaled by the wavelength: a 0 + a 1 • λ. Finally, a noise term n(λ) is also added to the model which corresponds to statistical fluctuations in the photoelectron number (Poisson distribution). Note that the noise variance depends on the wavelength and is directly related to the silicon sensitivity of the CCD sensor. The final model of our study is given in equation 2.

r meas (λ) = m • ρ(λ) + a 0 + a 1 • λ + n(λ)
(2)

Spectral pre-treatments

This study is focused on spectral pre-treatments that are able to reduce or remove the nuisance variability sources given in the equation 2, i.e. the multiplicative term m, the additive term a 0 and a 1 and the noise term n(λ) from the measured spectrum. Using this notation, the term ρ(λ) is referred to as useful spectral information while the others are referred to as a nuisance variability (for discrimination purposes).

The pre-treatments evaluated in this study that can cope with some parts of the nuisance variability sources are listed in Table 1. For clarity, a nomenclature to describe these pre-treatments according to their effects and scopes have been used. The term global refers to a pre-treatment that has a unique effect on a given spectrum without differentiating between wavelengths. The term piecewise describes a pre-treatment that uses its neighboring wavelengths. The term local refers to a pre-treatment that treats each wavelength independently. We also differentiate the actual pre-treatment effect using terms: reduce, transform or remove. For example, smoothing reduces the noise piecewise using a numerical filter of a given width, logarithm transforms the multiplicative effect (into additive) locally and centering removes the additive effect globally.

In the following, the main pre-treatments are detailed and their effect in feature spaces are described.

Dealing with additive effects

The centering pre-treatment (also known as constant detrending or detrending of order 0) and linear detrending (or detrending of order 1) have been chosen to deal with additive term a 0 and the linear trend a 0 + λ • a 1 respectively.

Their expressions in the feature space can be expressed using the same notation. Let D k be the matrix

D k = p 0 , p 1 , • • • , p k , where p i = [1 i , 2 i , • • • , p i ] T
Table 1: Pre-treatments evaluated in this study, their effect on spectrum: remove (rm), transform (tr) or decrease (dc) and their scope : global (G), piecewise (P), local (L) as described in Section 2.3.

Pre-treatment Abreviation Additive Multiplicative Noise Constant Linear None R Center D 0 rm -G Detrend linear D 1 rm -G Normalize std N std rm -G Normalize Norm 1 N 1 rm -G Normalize Norm 2 N 2 rm -G Center -Normalize std D 0 -N std rm -G rm -G Center -Normalize Norm 1 D 0 -N 1 rm -G rm -G Center -Normalize Norm 2 D 0 -N 2 rm -G rm -G Logarithm Log tr -L Logarithm -Center Log -D 0 rm -L Smoothing SG 0 dc -P 1st Derivative SG 1 rm -P dc -P 2nd Derivative SG 2 rm -P rm -P dc -P Smoothing -Center SG 0 -D 0 rm -G dc -P Smoothing -Detrend linear SG 0 -D 1 rm -G dc -P Smoothing -Normalize std SG 0 -N std rm -G dc -P Smoothing -Normalize Norm 1 SG 0 -N 1 rm -G dc -P Smoothing -Normalize Norm 2 SG 0 -N 2 rm -G dc -P Smoothing -Center -Normalize std SG 0 -D 0 -N std dc -P Smoothing -Center -Normalize Norm 1 SG 0 -D 0 -N 1 dc -P Smoothing -Center -Normalize Norm 2 SG 0 -D 0 -N 2 dc -P Smoothing -Logarithm SG 0 -Log tr -L dc -P Smoothing -Logarithm -Center SG 0 -Log -D 0 rm -L dc -P
and p corresponds to the number of wavelengths. Detrending of any order k corresponds to projecting the spectrum orthogonal to the matrix D k [START_REF] Boulet | Pretreatments by means of orthogonal projections[END_REF]:

r corr = I P -D k D T k D k -1 D T k r meas (3) If k = 0, p 0 = [1, • • • , 1]
T , the spectrum is projected on an hyperplane orthogonal to the vector composed only of ones. Similarly, linear detrending corresponds in the feature space to the projection of the spectrum orthogonal to the matrix D 1 . Note that in practice, as k increases, special care should be given for numerical inversion of the matrix D T k D k because of large differences between its eigenvalues [START_REF] Meyer | Matrix analysis and applied linear algebra[END_REF].

Expressed in R p , it is clear that the spectral dimension is reduced from p to p -1 and from p to p -2 for the detrend 0 and for the detrend 1 respectively. Indeed, the spectrum is projected orthogonal to a matrix of rank 1 and 2 respectively. Removing these dimensions when dealing with data in such a high-dimensional space might not appear problematic, especially since the actual dimension of visible and near infra red (V-NIR) spectra is usually quite low due to the high correlation between wavelengths [START_REF] Naes | A user friendly guide to multivariate calibration[END_REF]. However, when these pre-treatments are applied, the global intensity information is also definitively removed. Class separability therefore decreases if discriminatory information is colinear to these (removed) directions in the feature space.

Dealing with multiplicative effects

In this study, we considered two types of pre-treatment to deal with the multiplicative term, i.e., normalizing and logarithm transformation.

With the normalizing approach, a single factor computed from all wavelengths is used to normalize the spectrum. Several normalizing factors can be used and lead to different results that depend on the context, e.g. number of wavelengths of the spectral data, presence of outliers, intensity of the reflectance values, etc. Usual normalizing factors are norm 1, norm 2 or standard deviation (std) of the spectrum:

r corr = 1 N (r meas ) r meas (4)
where

N (•) = { • 1 ,
• 2 , std(•)} corresponds to the normalizing factor. Note that by using variation information around the mean, std is less sensitive to centered noise. In the feature space, normalizing using norm 1 corresponds to the projection of the spectrum onto the space orthogonal to the vector [1, • • • , 1] T in R p . When using norm 2, the spectrum is projected onto the unit sphere of dimension p. Finally, normalizing using standard deviation corresponds to the projection of the spectrum onto a (p -1)-dimensional hypersphere of unit radius (centered on its mean). Spectra with different angles (in R p ) before the pre-treatment remain different after projection. Therefore, if discriminatory information between spectra is in a space spanned by a difference in intensity (length of the vector in R p ), the information is lost after projection.

Another approach to deal with multiplicative term is the logarithm transformation which is described as follows:

r corr = log (r meas ) (5)
, where the log(•) operator is applied to every element of r meas individually.

In the feature space, this transformation corresponds to using a logarithmic scale instead of a linear one. Doing so, no spectral dimension is lost and the transformation is reversible. However, logarithm transformation itself does not remove the multiplicative term, it only transforms the multiplicative term into an additive term.

Dealing with noise: smoothing and derivative

Savitzky-Golay smoothing [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF]) performs a piecewise polynomial fitting with specified width and order to the spectrum. This smoothing method has been successfully adopted by the spectroscopy community because it is able to reduce noise while keeping the meaningful variation that occurs at these wavelengths.

In the individual space, it increases the correlation between neighboring wavelengths and reduces random noise leading to smoother curves. In the feature space, it increases the correlation between variables, therefore increasing the matrix column correlations and decreasing the real rank of the spectral data.

Savitzky and Golay have also provided a smoothed version of first and second numerical derivatives with a high signal to noise ratio. Some theoretical advantages and practical properties of Savitzky-Golay derivatives have been studied [START_REF] Luo | Properties of SavitzkyGolay digital differentiators[END_REF]. In our study, first derivatives are used to remove baselines and second derivatives to removes slopes.

Combination of pre-treatments

The previously described pre-treatments taken individually are designed to correct only for a part of the nuisance variability of the model given in equation 2. We must therefore combine them optimally to attempt to remove more of the nuisance variability. Table 1 lists every combination of pre-treatments evaluated in our study.

Because the order of applying different pre-treatments is important, a consistent notation ordered from left to right is used in this paper. For example, Smoothing -Center -Normalize Norm 1, corresponds to apply first a smoothing, then a centering and finally a normalization (with norm 1).

Among these combinations, the most widely used in spectral data analysis is called Standard Normal Variate (SNV) (e.g. see [START_REF] Sun | The vegetation classification in coal mine overburden dump using canopy spectral reflectance[END_REF]) and it 8

Author-produced version of the article published in Computers and Electronics in Agriculture, 2014, N°108(), p.242-249 The original publication is available at http://www.sciencedirect.com Doi: 10.1016/j.compag.2014.08.010 combines centering and normalizing using std:

r SN V = 1 std (r meas ) I P -p 0 p T 0 p 0 -1 p T 0 r meas . (6)
In the feature space, a spectrum is first mapped on the hyperplane of dimension (p -1) orthogonal to [1, • • • , 1] T . The resulting vector is then projected onto a (p -1)-dimensional hypersphere of unit radius. Geometric details of this transformation are well described in [START_REF] Fearn | On the geometry of SNV and MSC[END_REF].

Material and methods

Experimental setup

Short-range in-field hyperspectral images were acquired using a camera (Hyspex V-NIR 1600, Norks Elektro Optikk, Norway) with a spatial resolution of 0.2 mm/pixel at a distance of 1 meter above the ground. The acquired spectra were composed of 160 spectral bands ranging from 415 nm to 1000 nm. Each image was composed of 1600 pixels per row and the number of rows, that depended on the scanning length, was in average 3000. A typical RGB image reconstructed from the hyperspectral image is represented in Figure 1. In this figure, weeds are also highlighted using the Normalized Difference Vegetation Index (NDVI) [START_REF] Rouse | Monitoring vegetation systems in the Great Plains with ERTS[END_REF].

Using a calibrated reference surface in each image (see on the right-hand side of Figure 1), radiance values were transformed into reflectance values. In doing so, the obtained reflectance images were almost independent of the atmospheric conditions and illumination. Image acquisitions were carried out over two consecutive years, at a growth stage that corresponds to the weeding period, and under various lighting conditions, i.e. sunny and cloudy weather, shaded areas, etc. In each image, weeds (dicotyledons), wheat, soil and the reference surface were in the camera field of view as represented in Figure 1. In total, 9 images, in which each weed position was recorded in the field, were used for the study and were randomly separated into two groups: 4 images for the training and 5 images for the test. A ground truth map of weed and wheat was then manually created for each image. Finally, around 1600 spectra for training set and 1600 spectra for the test set were randomly selected and labeled as either wheat or weed according to the ground truth.

We applied every pre-treatment listed in Table 1 to the training and test sets. For Savitzky-Golay smoothing and derivatives, the parameters that gave the best results using a cross-validation [START_REF] Esbensen | Principles of Proper Validation: use and abuse of re-sampling for validation[END_REF] on the training set were a third order polynomial and a kernel width of 17.

Concerning the two-class discrimination problem (wheat and weed), we used Partial Least Square Linear Discriminant Analysis (PLS-LDA) and Support Vector Machine with a Radial Basis Function Kernel (Gaussian SVM). Every model (pre-treatment and classifier) was trained and tuned using a cross-validation on the training set and its performance was then assessed on the independent test set [START_REF] Esbensen | Principles of Proper Validation: use and abuse of re-sampling for validation[END_REF]. Both classifiers, i.e. PLS-LDA and Gaussian SVM were compared in terms of the prediction error (on the test set).

For Gaussian SVM, the model complexity was also reported as a way to understand how complex the discrimination task as a function of the pretreatment.

For PLS-LDA, bias and class separability of the scores were also measured (Λ W ilks ) [START_REF] Everitt | The Cambridge dictionary of statistics[END_REF] on the PLS-LDA scores. In Figure 2 the bias and separability are illustrated for our two class discrimination problem.

Experimental results

Pre-treatment effects on the training spectra are displayed in Figure 3 with green for wheat and red for weed. In each class, the solid line corresponds to the average spectrum and the plain color area represents its variability (standard deviation) within the training set. Raw spectra (without any pre-treatment) are shown in Figure 3a. Additive corrections are represented in Figures 3c (detrend of order 0) and Figure 3d (linear detrend). Figures 3b,3e, 3f and 3g correspond to multiplicative corrections using logarithm, normalizing factors norm 1, norm 2 and std respectively. Smoothing, first and second derivative are shown in Figures 3h, Figures 3i andFigures 3j.

Discrimination results obtained with every combination of pre-treatments (detailed in Table 1) are displayed in Figure 4 using a radar chart. Each edge corresponds to a specific pre-treatment and each color to a discrimination method, i.e. blue for Gaussian SVM and red for PLS-LDA. Gaussian SVM model complexity is defined by the number of support vectors used to create the model divided by the available number of training samples and is displayed (in %) with a dashed blue line. For PLS-LDA, only the rate 11

Author-produced version of the article published in Computers and Electronics in Agriculture, 2014, N°108(), p.242-249 The original publication is available at http://www.sciencedirect.com Doi: 10.1016/j.compag.2014.08.010 Wheat spectra are in green and weed spectra in red. Filled areas correspond to spectral variability within each class. Abbreviations are detailed in Table 1.

(a) R (b) Log (c) D 0 (d) D 1 (e) N 1 (f) N 2 (g) N std (h) SG 0 (i) SG 1 (j) SG 2
Figure 4: Radar chart displaying the rate of good discrimination in (%) for both discrimination methods : PLS-LDA (red), and Gaussian SVM (blue). SVM model complexity (dotted blue) corresponds to the percentage of support vectors used to create the model (in the available training set). Gray area corresponds to result where at least a smoothing has been applied. Abbreviations are detailed in Table 1.

of good discrimination was considered since the same model complexity (6 latent variables) was obtained with every pre-treatment.

The separability obtained with the training (blue) and test set (red) are represented in Figure 5a. The bias obtained from the training set to the test set is represented in Figure 5b. In both radar charts results are expressed in percentages and relative to the results obtained with raw spectra. Note that on the radar chart, a result lying outside the main area corresponds to an unfavorable pre-treatment and a result close to the center corresponds to a useful pre-treatment.

Discussion

Effects on spectra

Raw vegetation spectra in Figure 3a clearly indicate a plateau in the near infra red (NIR) region (750 nm -1000 nm). The increasing slope between 670 nm and 750 nm, called red-edge, is a robust indicator to discriminate vegetation from the background [START_REF] Brown | Site-specific weed management: sensing requirements -What do we need to see?[END_REF], [START_REF] Slaughter | Autonomous robotic weed control systems: A review[END_REF]. The bump between 500 nm and 600 nm corresponds to the green color of vegetation and the strong absorption at around 670 nm indicates the presence of chlorophyll in the leaves. The peak at 760 nm indicates that oxygen absorption was not accurately corrected using the reference surface (due to low signals creating an large amount of quantification noise). Above 900 nm, spectra are noisier because of a lower sensitivity of the CCD sensor in the NIR.

As for the discrimination between these two classes, the NIR plateau is on average higher and the red-edge steeper for wheat spectra than for weed spectra. A slight difference in greenness between both vegetation spectra is also noticeable between 500 nm and 600 nm.

Within each class a large amplitude variation at every wavelength is clearly visible. Comparing Figures 3c andFigures 3d (additive corrections) with Figures 3e, 3f and 3g (multiplicative corrections) indicates that this spectral variation is mostly due to a multiplicative effect.

After smoothing, spectra clearly indicate less random noise in NIR and the 760 nm peak is reduced. However, main spectral features remained, e.g. red edge, plateau and greenness differences between weed and wheat.

With derivatives, visual interpretation of the curves is easier since derivative increases small variations. Boundary effects are also visible at both ends of the spectra due to numerical filtering computational issues. 1.

Discrimination

Curves presented in the previous section demonstrate the difficulty of the discrimination task between wheat and weeds despite any sort of pretreatment. There is in fact very little difference between weed and wheat spectra as well as the variability within each class being large.

However, most pre-treatments evaluated in this study were effective in dealing with the discrimination problem. Both PLS-LDA and Gaussian SVM obtained a discrimination accuracy of more than 80% which demonstrates their potential in spectral analysis. With raw data (R), high discrimination rates are observed for both methods, i.e. 87% for SVM and 80% for PLS-LDA. The SVM model used 60% of the available training spectra to create its model. Around a thousand points in a 160-dimensional vector space is still low and confirms the high correlation of the spectral data (the space is almost empty). PLS-LDA appears more sensitive to the pre-treatment than Gaussian SVM, the latter compensates for the loss of information by creating a more complex non-linear model.

The best pre-treatment is clearly the logarithm transformation, especially when used with PLS-LDA (96%). Logarithm does not suppress multiplicative effect, it transforms this effect into an additive one which is then entirely handled by PLS-LDA.

On the other hand, normalization forced spectra to lie in a lower dimensional subspace by forcing low and high intensity values to be on the same scale. Therefore, some discriminatory features such as NIR plateau value and greenness were diminished. Gaussian SVM model complexity indicates that problems occurred when using normalization of any sort. The obtained model was therefore more complex and only based on subtle features that were still present in the corrected spectra. For PLS-LDA, however, normalization gave better results than with raw spectra. Since PLS-LDA cannot deal with a multiplicative effect, it was better to lose some information than to be perturbed by the multiplicative effect.

For both classifiers, smoothing appears to be an important first pretreatment step by reinforcing the correlation between neighboring wavelengths and by reducing noise. It thus helps the discrimination method to build simpler models that are less sensitive to an error when using a specific wavelength.

Removing additive effects (by centering or removing linear trend) did not improve PLS-LDA discrimination results since it can handle additive effects on its own if they are present within the training set. On the contrary, Gaussian SVM cannot deal with additive effects even if they are present in the training set, i.e. it learns the discrimination boundaries knowing only some examples of translated vectors. By removing additive effects, SVM model complexity decreased therefore helping the model to find a better class separation.

Using the first and second derivative does not significantly improve discrimination results for either classifier. It provides slightly lower Gaussian SVM model complexity while keeping the discrimination result constant. PLS-LDA is not greatly affected by this pre-treatment step since it can handle the same effects on its own.

PLS-LDA separability and bias

In this study, spectra came from images acquired at different times and parts of the field, creating between-image variability. Since spectra from the training set were taken from images different from those of the test set this variability could not be learned by the model. PLS-LDA first step is to reduce the dimension of the space. The class separability and the bias were thus computed in this lower dimensional space avoiding statistical estimation problems that occur in high-dimensional spaces [START_REF] Hughes | On the mean accuracy of statistical pattern recognizers[END_REF]).

In the radar chart of Figure 5a, similar patterns are observed with the training set and the test set. The separability was always worse with the test set than with the training set (except when using the logarithm). This difference indicates a slight over-fitting of the model to the training set even using a proper validation method. None of the pre-treatments was designed to handle the bias occurring between training set and test set, however, from the chart in Figures 5b, we observe they all had a positive influence on the bias.

Generally, removing multiplicative effects gave more separability than with raw data. Logarithm transformation provided the best results in terms of both separability and bias even without smoothing and/or centering. Logarithm kept all the data information by only transforming the multiplicative effect. The discrimination model was therefore based on the main spectral differences and was more robust to small changes occurring between the training set and the test set. With normalizing, results of each factor were nearly equivalent with only a slightly poorer result for standard deviation. Standard deviation, as explained in section 2.3.2, makes more assumptions about the data distributions.

Removing additive effects by detrending or derivative considerably reduced the class separability and gave worse results than with raw data. Indeed, when using PLS-LDA, if the additive effect is present in the training set, there is no need to remove it through any pre-treatment. In terms of bias, almost no improvement is noticeable by removing additive effect.

Smoothing had an large positive effect in terms of bias but had a negative effect in terms of separability. The bad effect on separability was likely caused by the filter being too strong, thus resulting in the removal of some useful discrimination information. On the other hand, smoothing removed small variations that might have occurred between images thus reducing the bias.

Conclusion

This study examines the influence of pre-treatments on the discrimination of wheat and weed spectra in outdoor conditions. We demonstrated that to deal with the modeled nuisance variability the most effective combination of pre-treatment and discrimination methods was the logarithm transformation with PLS-LDA. Logarithm was the only pre-treatment that dealt with multiplicative effect without explicitly projecting the data on a lower dimensional subspace in which useful discriminatory information was lost.

Experimental results also validate the model of measured reflectance for outdoor application with uncontrolled light, and leaf orientations. This approach will however require further studies to increase the understanding of the bias occurring between images. This would help to create a more detailed model and hopefully some treatments that can handle this bias.

Figure 1 :

 1 Figure 1: Illustration of a typical image from the data set. (top) RGB image reconstructed from the hyperspectral image. Note the reference surface used for reflectance correction on the right-hand side. (bottom) NDVI images zoomed on the weeds.

  . The bias corresponded to a shift occurring between the best separator in the training set and the separator that would be the best in the test set. The class separability was measured by computing the Wilk's lamdba 10 Author-produced version of the article published in Computers and Electronics in Agriculture, 2014, N°108(), p.242-249 The original publication is available at http://www.sciencedirect.com Doi: 10.1016/j.compag.2014.08.010

Figure 2 :

 2 Figure 2: Illustration of the bias and separability for PLS-LDA performance assessment. Green and red represent two different classes. Dotted lines represent a training set and plain lines represent a test set.

Figure 3 :

 3 Figure 3: Pre-treatment effects on real reflectance spectra taken from the training set.Wheat spectra are in green and weed spectra in red. Filled areas correspond to spectral variability within each class. Abbreviations are detailed in Table1.

  raining → T est) Bias(T raining → T est)(R)

Figure 5 :

 5 Figure 5: Radar charts obtained with PLS-LDA. (a) the class separability for the training (red) and the test set (blue). The value has been inverted for visualization purpose. (b) the bias between training and test. Results are expressed in percent and relative to the results obtained with raw spectra. Gray area corresponds to results where at least a smoothing has been applied. Abbreviations are detailed in Table1.

Author-produced version of the article published in Computers and Electronics in Agriculture, 2014, N°108(), p.242-249 The original publication is available at http://www.sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.compag.2014.08.010 .08.010

Author-produced version of the article published in Computers and Electronics in Agriculture, 2014, N°108(), p.242-249 The original publication is available at http://www.sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.compag.2014.08.010 .08.010

Author-produced version of the article published in Computers and Electronics in Agriculture, 2014, N°108(), p.242-249 The original publication is available at http://www. sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.compag.2014.08.010 .08.010


R Kohlhaas

A Bertoldi

E Cantin

A Aspect

A Landragin

P Bouyer

Phase locking a clock oscillator to a coherent atomic ensemble

Keywords: numbers: 95.55.Sh, 42.62.Eh, 67.85.-d, 37.25.+k

The sensitivity of an atomic interferometer increases when the phase evolution of its quantum superposition state is measured over a longer interrogation interval. In practice, a limit is set by the measurement process, which returns not the phase, but its projection in terms of population difference on two energetic levels. The phase interval over which the relation can be inverted is thus limited to the interval [-π/2, π/2]; going beyond it introduces an ambiguity in the read out, hence a sensitivity loss. Here, we extend the unambiguous interval to probe the phase evolution of an atomic ensemble using coherence preserving measurements and phase corrections, and demonstrate the phase lock of the clock oscillator to an atomic superposition state. We propose a protocol based on the phase lock to improve atomic clocks under local oscillator noise, and foresee the application to other atomic interferometers such as inertial sensors.

From the first observations of Huygens on the coordinated motion of coupled nonlinear oscillators [START_REF] Huygens | OEuvres Complètes de Christiaan Huygens[END_REF], phase synchronization has evolved to an indispensable tool for time and frequency metrology and starts to be investigated for quantum systems [START_REF] Mari | Measures of quantum synchronization in continuous variable systems[END_REF][START_REF] Xu | Quantum synchronization of two ensembles of atoms[END_REF][START_REF] Cox | Phase diagram for injection locking a superradiant laser[END_REF]. Phase lock loops (PLLs) [START_REF] Guan-Chyun | Phase-locked loop techniques. a survey[END_REF], where a local oscillator (LO) is phase locked to a reference signal, are widely used for the generation of atomic time scales [START_REF] Jones | Carrierenvelope phase control of femtosecond mode-locked lasers and direct optical frequency synthesis[END_REF][START_REF] Diddams | Direct link between microwave and optical frequencies with a 300 thz femtosecond laser comb[END_REF], the synchronization in telecommunication [START_REF] Slavík | All-optical phase and amplitude regenerator for next-generation telecommunications systems[END_REF] or in radio navigation [START_REF] Braasch | GPS receiver architectures and measurements[END_REF]. On the contrary, nowadays in atomic frequency standards, only the frequency of the local oscillator is locked on the atomic resonance. This is directly linked to the quantum nature of the reference system, i.e., the quantum superposition of two internal states of an ensemble of atoms, molecules or ions, which is destroyed by the detection at the end of the interrogation process. A similar limitation exists for any measurement of a quantum system as in magnetometers or inertial sensors [START_REF] Arndt | De Broglie's meter stick: making measurements with matter waves[END_REF]. Adopting PLLs in such systems will improve their long term stability, as the phase is the integral of the frequency, and reduce the constraints on the stability of the LO or of the measured signal. Here, we demonstrate the direct phase lock of a LO to an atomic ensemble, based on repeated coherence preserving measurements of the atomic ensemble. We also study how this technology could improve atomic clocks subject to local oscillator noise.

Phase locking a classical system to a quantum system would give a direct link in metrology to the fundamental oscillations of quantum particles and could lead to enhanced sensitivities and new applications in precision measurements. In an atomic clock, the frequency of a LO is repeatedly referenced to an atomic transition frequency by comparing their respective phase evolutions in an interrogation time T and applying a feedback correction. During the interrogation, the atoms are in a superposition state, and the projection of the relative phase between the LO and the atomic ensemble is measured as a population imbalance of the two clock levels. This leads to a sinusoidal signal, and the phase drift can therefore only be unambiguously determined if it stays within the [-π/2, π/2] interval, hereafter called inversion region. Hence, for a given LO noise, the interrogation time of the atomic transition must be kept short enough such that phase drifts beyond the inversion region are avoided. Currently, LO noise limits the interrogation time in ion [START_REF] Chou | Frequency comparison of two high-accuracy Al+ optical clocks[END_REF] and optical lattice clocks [START_REF] Hinkley | An atomic clock with 10 -18 instability[END_REF][START_REF] Bloom | An optical lattice clock with accuracy and stability at the 10 -18 level[END_REF][START_REF] Ushijima | Cryogenic optical lattice clocks with a relative frequency difference of 1×10 -18[END_REF] and is expected to become a limit for microwave clocks with the recently discovered spin-self rephasing effect [START_REF] Büning | Extended coherence time on the clock transition of optically trapped rubidium[END_REF]. The standard approach to tackle this issue consists in improving the quality of local oscillators [START_REF] Jiang | Making optical atomic clocks more stable with 10 -16 -level laser stabilization[END_REF][START_REF] Kessler | A sub-40-mHz-linewidth laser based on a silicon singlecrystal optical cavity[END_REF][START_REF] Thorpe | Frequency stabilization to 6×10 -16 via spectral-hole burning[END_REF][START_REF] Cole | Tenfold reduction of Brownian noise in highreflectivity optical coatings[END_REF][START_REF] Amairi | Reducing the effect of thermal noise in optical cavities[END_REF]. As an alternative, it has been recently proposed to track and stabilize the LO phase evolution using several atomic ensembles probed with increasing interrogation time [START_REF] Rosenband | Exponential scaling of clock stability with atom number[END_REF][START_REF] Kessler | Heisenberg-limited atom clocks based on entangled qubits[END_REF], or by enhancing the Ramsey interrogation interval by stabilizing the LO either via cascaded frequency corrections [START_REF] Borregaard | Efficient atomic clocks operated with several atomic ensembles[END_REF] or by coherence preserving measurements on the same atomic ensemble and feedback [START_REF] Shiga | Locking the local oscillator phase to the atomic phase via weak measurement[END_REF][START_REF] Shiga | Accelerating the averaging rate of atomic ensemble clock stability using atomic phase lock[END_REF]. We follow the latter proposal, and show how the LO limitations can be bypassed by implementing a phase lock between the LO and the atomic system. We begin with a minimally destructive measurement of the LO phase drift when it is within the inversion region; the measurement readout is then used to correct the LO phase so as to reduce the phase drift. The cycle is repeated using the residual atomic coherence. This results in successive, phase related measurements of the relative phase evolution, and the feedback keeps the phase in the inversion region, leading to an effectively longer interrogation time. We demonstrate this approach with a trapped ensemble of neutral atoms probed on a microwave transition.

The experimental scheme shown in Fig. 1 has been described in [START_REF] Vanderbruggen | Feedback control of trapped coherent atomic ensembles[END_REF]. A cloud of cold 87 Rb atoms is trapped FIG. 1: Experimental scheme. The evolution of the LO phase ϕLO is compared to the phase ϕat of an atomic ensemble in a superposition state using a coherence preserving measurement in a Ramsey spectroscopy sequence. The relative phase is obtained from the read out of the population difference, and is used to implement the phase lock between the two oscillators by applying feedback on a phase actuator on the LO output. The light shift induced by the optical trap and by the probe has been engineered to have a homogeneous measurement on the atomic ensemble [START_REF] Vanderbruggen | Feedback control of trapped coherent atomic ensembles[END_REF].

in an optical potential (see Appendix A), prepared with a π/2 pulse of a resonant microwave field in a balanced superposition state of two hyperfine levels

| ↓ ≡ | F = 1, m F = 0 and | ↑ ≡ | F = 2, m F = 0 of
the electronic ground state, and probed using a nondestructive detection. Typically, 5 × 10 5 atoms at a temperature of 10 µK are used in the measurements reported here. The phase ϕ at of the superposition state oscillates at a frequency of 6.835 GHz corresponding to the energy difference between the | ↓ and | ↑ atomic states, which is the fundamental reference if atoms are protected from perturbations. A microwave LO has a frequency close to the atomic frequency difference, so that the relative phase ϕ = ϕ LO -ϕ at between the two oscillators drifts slowly because of the LO noise. ϕ can be measured using the Ramsey spectroscopy method (see Fig. 2): a second π/2 microwave pulse (projection pulse) maps it onto a population difference, which we read out with a weak optical probe perturbing the atomic quantum state only negligibly and preserving the ensemble coherence [START_REF] Vanderbruggen | Feedback control of trapped coherent atomic ensembles[END_REF][START_REF] Smith | Efficient quantum-state estimation by continuous weak measurement and dynamical control[END_REF][START_REF] Lloyd | Quantum feedback with weak measurements[END_REF]. Unlike for destructive measurements, the interrogation of ϕ can continue in a correlated way, once the action of the projection pulse is inverted using an opposite π/2 microwave pulse (reintroduction pulse), which brings the atomic state back to the previous coherent superposition. Moreover, after each measurement and reintroduction pulse, the phase read out can be used to correct the LO phase. The evolution and manipulation of the atomic ensemble can be illustrated using the Bloch sphere representation (Fig. 2): the collective state of N at two-level atoms in the same pure single particle state (also called coherent spin state (CSS)) forms a pseudo-spin with length J = N at /2, where J z denotes ). The measurement of the relative phase between the CSS and the LO starts when a π/2 microwave pulse around the y axis brings the CSS into a balanced superposition of the | ↓ and | ↑ states, depicted as a vector on the equatorial plane of the Bloch sphere (step 2.). The x axis is chosen to represent the phase of the local oscillator, and ϕ the relative phase between the LO and the atomic superposition that evolves because of the LO noise (step 3.). After an interrogation time T, the projection of ϕ is mapped onto a population difference by a projection π/2 pulse around the x axis and read out with the coherence preserving detection (step 4.). The CSS is rotated back to the equatorial plane by a reintroduction π/2 pulse around the x axis (step 5.), and feedback is applied on the phase of the LO (step 6.). The PLL between the LO and the atomic ensemble consists in the repetition of the steps from 3.) to 6.), potentially till the atomic ensemble shows a residual coherence.

the population difference and ϕ = arcsin (J y /J x ) is the phase difference between the phase of the LO and that of the superposition state. A resonant microwave pulse determines the rotation of J around an axis in the equatorial plane of the Bloch sphere, and the axis direction is set by the phase of the microwave signal. The repetition of the manipulation, measurement and feedback cycle implements the phase lock of the LO on the atomic superposition state, as shown in steps 2.-6. of Fig. 2.

We first show that we can reconstruct the time evolution of the relative phase between the LO and the CSS by monitoring the population difference and without applying feedback. For this purpose, we frequency offset the LO by 100 Hz from the nominal resonance, and periodically measure the projection of the relative phase (sin(ϕ)) with only a small reduction of the atomic ensemble coherence (Fig. 3). Every 1 ms ϕ is mapped to a population difference via a projection π/2 microwave pulse around the x-axis, a weak measurement of J z is performed, and the collective spin is rotated back to the equatorial plane of the Bloch sphere via a reinstroduction π/2 pulse around the x-axis. The π/2 microwave pulses are derived from an amplified version of the LO at 6.835 GHz leading to a pulse length of τ π/2 =47 µs, and the rotation axis is controlled with a quadrature phase shifter. The coherence preserving, dispersive measurement relies on frequency modulation spectroscopy (see Appendix B). Fig. 3 shows in a single experimental run how the spin state evolves around the equator of the Bloch sphere, with the phase ϕ mapped on the normalized population difference J z /J. The S/N of the weak measurements is 20 for a full state coherence and each readout of the relative phase drift reduces the state coherence by 2%. The destructivity from the probe is the main decoherence source till 10 ms, then the inhomogeneous light shift of the dipole trap on the clock states becomes the dominant decoherence source. We next introduce feedback and demonstrate that we can phase lock the LO on the atomic superposition state, and increase the Ramsey interrogation time beyond the limit set by the inversion region between J z and the relative phase. We apply on the local oscillator two types of signals, first a frequency offset, and second periodic phase jumps, and use the output of the coherence preserving measurements to actively minimize ϕ. The phase lock is obtained by controlling the phase of the local oscillator by means of a digital phase shifter (see Appendix D). The feedback is performed after the atomic spin is rotated back to the equatorial plane of the Bloch sphere. When the disturbance applied on the local oscillator con-sists of a frequency offset, there is a linear phase drift between the LO and the atomic phase (Fig 4 ,red). The phase evolution in open loop is reconstructed from the data of Fig. 3 by taking into account the damping on the sinusoidal signal due to the decoherence sources, and knowing that a constant frequency offset is applied on the LO. We remark that sudden sign inversions of the applied frequency offset when ϕ = ±π/2 would produce exactly the same evolution of the population difference, illustrating the need to keep ϕ in the inversion region. In closed loop, the phase drift due to the 100 Hz frequency offset on the LO is periodically reset to zero, with a precision set by the π/32 step size of the digital phase shifter and the uncertainty of the coherence preserving measurements. This results in a saw-tooth-like signal for J z /J (Fig. 4, blue signal). Without phase lock, the phase drift leaves the inversion region after 2.5 ms and rotates several times around the Bloch sphere, whereas with phase lock it stays in the inversion region for all the 22 ms interval shown in the image. When the feedback is active the total phase drift results as the phase measured at the end of the Ramsey interferometer, added to the correction phase shifts on the LO via the feedback controller. We next apply periodic phase jumps of π/3 back and forth on the LO using a second phase shifter. The signal obtained in open loop is shown at the top of Fig. 5. When the feedback controller is active, the jumps detected on the relative phase are corrected to zero, with a precision set by the resolution of the phase shifter and the uncertainty of the weak measurements (Fig. 5, bottom). The solid lines in Fig. 4 and 5 are drawn from the known timing for the applied phase signal and the feedback on the phase. For a combination of a phase drift and phase jumps, the relative phase can leave the inversion region while the noise action cannot be predicted from previous measurements. This is the commonly encountered situation in atomic clocks, and highlights the requirement of feedback on the LO phase to keep track of the relative phase drifts. Without feedback, the Ramsey interrogation time should be kept sufficiently short to avoid ambiguities for the measured phase shift.

We propose now a protocol to efficiently use the phase lock to improve an atomic clock. In a conventional atomic clock, the phase drift ϕ is destructively read out after a single interrogation and feedback is performed on the LO by the addition of a frequency ω FB = -ϕ/T considering unity gain. Our protocol of using the PLL between the LO and the atomic superposition state in an atomic clock is based on the reconstruction of the phase drift experienced by the LO over the extended interrogation time T tot =N×T (N is the number of phase coherent interrogations) by combining the phase shifts applied by feedback and the final phase readout (Fig. 6). The known phase corrections applied to the LO phase serve the dual purpose of keeping the relative phase in the inversion region and giving a coarse estimate of the phase drift during T tot . The final phase measurement ϕ f , together with the phase corrections ϕ FIG. 4: Phase lock between the LO and the atomic superposition state. The evolution of the LO-atom relative phase is reconstructed from the Jz signal in Fig. 3 (red, solid line) and when feedback is applied on the phase of the LO after each measurement (blue, solid line). In the open loop case, the points are reported at the values given by the fit in the previous figure. In the closed loop case, the relative phase is always smaller than π/2, and does not enter in the region represented in grey where it cannot be univocally determined from the measurement. FIG. 5: Correction of phase jumps between the LO and the atomic ensemble. Evolution of the LO-CSS phase when periodic phase jumps of π/3 are applied back and forth on the LO: the points are obtained from the measured population difference, whereas the solid line represents the LO phase measured after the correction phase actuator. The servo controller, implemented using coherence preserving measurements and feedback on a phase actuator, maintains the relative phase close to zero. Above the system is operated in open loop, below in closed loop. In the latter case, the phase corrections are applied 150 µs after each measurement. estimate of the phase drift during T tot . The total phase drift can be computed as

ϕ tot = ϕ f - N i=1 ϕ (i)
FB , where the sum is over all the correction phase shifts applied by the feedback controller and stored in the microcontroller during the sequence. The feedback on the frequency is set accordingly to be ω FB = -ϕ tot /T tot . After the final phase read out, the phase shifter is reset to its initial position to avoid any impact of the intermediate phase shifts on the long term clock stability. The measurement S/N of the total phase drift depends only on the final phase measurement, since it can make use of all the residual ensemble coherence, whereas the S/N of the intermediate measurements has to be only sufficient to keep the atomic state in the inversion region. To further improve the protocol, the clock operation could be optimized by adapting the destructivity of the intermediate measurements to the varying coherence of the atomic sample, or implementing an adaptive measurement protocol for the last measurement [START_REF] Borregaard | Near-Heisenberglimited atomic clocks in the presence of decoherence[END_REF].

For a proof-of-concept demonstration, we run an atomic clock that exploits the PLL between the LO and the atomic superposition state, and the phase reconstruction protocol. The LO signal was deteriorated so as to have an increased phase drift over the interrogation interval (see Appendix C). During the clock operation, the dead time for the preparation of the new ensemble is T d = 1.9 s and the interrogation time is set to T = 1 ms, much shorter than the measured atomic coherence lifetime. As a benchmark, we first run an atomic clock adopting a standard Ramsey interrogation sequence, and without any feedback on the phase. The phase measurement is performed with the coherence preserving detection adopted for the phase lock. The two-sample Allan frequency standard deviation was calculated from the sum of the LO noise and the correction signal applied by the feedback controller. The clock instability reaches a τ -1/2 scaling after a few clock cycles (see Fig. 7, red) at a level consistent with an initial S/N = 20 and considering the coherence decay in the optical trap. The instability is far higher than with state-of-the art atomic clocks since the experimental setup was not explicitly designed for the operation of an atomic clock, and the Allan frequency standard deviation is 1.5 × 10 -9 at 1 s. We then operate an atomic clock making use of a PLL sequence with N=9 successive interrogations, and again with T=1 ms. For simplicity, the intermediate and the final phase readouts are set to have the same measurement strength. The Allan frequency deviation shows a τ -1/2 scaling as expected for atomic clocks, which demonstrates that the phase reconstruction protocol is working properly. In the opposite case, frequency offsets would be corrected only at a short time with the phase actuator, and the stability of the clock would then diverge because of the non-zero dead time interval. The comparison of the instability shows that the clock adopting the phase lock and reconstruction method is at a lower level by a factor (4.76±0.25) with respect to the clock implementing the standard Ramsey interrogation, as shown in Fig. 7 (blue). This is clearly above the factor 3 expected for a succession of 9 uncorrelated phase drift measurements at each clock cycle [START_REF] Lodewyck | Nondestructive measurement of the transition probability in a Sr optical lattice clock[END_REF]. Nevertheless, several detrimental effects prevent us to reach the factor 9 expected for perfectly correlated measurements of the phase drift. The main effects are the reduced S/N due to the cumulated destructivity from the probe and the decay from the optical dipole trap, and the finite phase shifter accuracy, Operation of the clock: at each cycle of duration TC, the relative phase is repeatedly measured in a coherence preserving way during the phase lock interval (above: shaded grey areas). Each interrogation, represented in the inset by a light red peak between the manipulation π/2 pulses, is followed by a phase correction ϕ

(i)
FB on the LO, represented in the inset by a light blue peak in the right after the reintroduction of the spin on the equatorial plane of the Bloch sphere. The final phase readout ϕ f (dark red peak in the inset), whose S/N is set by the residual coherence, together with the previously applied phase shifts on the LO, provides the total phase drift ϕ experienced during the extended interrogation interval Ttot=N×T. The interrogation sequence ends with the application of a frequency correction on the LO (dark blue peak in the inset), then a new atomic ensemble is prepared in the dead time interval TD for the next cycle.

equal to 2.2 • .
The phase lock can be performed as long as the coherence of the state is maintained. For integration times longer than the coherence lifetime of the trapped ensemble, the atomic phase is lost and our locking scheme becomes again a frequency lock, like when the quantum superposition is destroyed by the dectection. In our experiment, the coherence lifetime is limited to 20 ms by the dephasing in the optical dipole trap. Nevertheless, trapped induced dephasing of the atomic state can be suppressed for 87 Rb as reported in [START_REF] Büning | Extended coherence time on the clock transition of optically trapped rubidium[END_REF][START_REF] Radnaev | A quantum memory with telecom-wavelength conversion[END_REF], whereas in an optical lattice it is strongly reduced with the choice of light at the magic wavelength [START_REF] Katori | Ultrastable optical clock with neutral atoms in an engineered light shift trap[END_REF]. In the original proposal to lock the local oscillator phase on the atomic phase [START_REF] Shiga | Locking the local oscillator phase to the atomic phase via weak measurement[END_REF], frequency feedback on the local oscillator after each weak J z measurement is performed. This leads to a longer effective interrogation time, but to a S/N given by the weak measurements, which is lower than that of a measurement at the quantum projection noise and beyond. Our protocol can overcome this limit, and reach projection limited readout while keeping an extended interrogation time, thanks to the feedback on the LO phase.

In the phase lock sequence, several effects must be considered to maximize the S/N of the last measurement while maintaining a high accuracy on the total phase drift over the increased interrogation time: the rotations FIG. 7: Atomic clock implementing a PLL. (Top) Allan frequency standard deviation σ1 for a normal Ramsey clock with interrogation time T=1 ms (red) and σ9 for a clock implementing the phase lock between the LO and the atomic superposition state for 9 successive, correlated interrogations on the same atomic ensemble, for a total interrogation time of 9×T=9 ms (blue). The red and blue lines are fits to the data, with a slope set to τ -1/2 . The dashed black line lies a factor √ 9 below the red curve, and represents the best achievable level if 9 uncorrelated measurements are implemented at each clock cycle. The continuous black line lies a factor 9 below the red curve, and represents the best achievable level if 9 correlated measurements are implemented at each clock cycle. (Bottom) The grey triangles represent the ratio of the Allan deviation for the two clocks; its value is higher than 3 (black dashed line), which is the maximum achievable by realizing at each clock cycle 9 uncorrelated measurements of the phase drift. The black solid line represent the maximum ratio achievable for 9 correlated measurements. operated on the Bloch sphere must be fast, the measurement induced decoherence limited, and the phase shifter used for the correction accurate. The decoherence related to the repeated interrogations of the relative phase can be strongly reduced by the use of an optical cavity to enhance the probe interaction with the atomic ensemble [START_REF] Bohnet | Reduced spin measurement back-action for a phase sensitivity ten times beyond the standard quantum limit[END_REF][START_REF] Lee | Many-atom-cavity QED system with homogeneous atom-cavity coupling[END_REF][START_REF] Leroux | Implementation of cavity squeezing of a collective atomic spin[END_REF]. An optimized clock configuration would consist in a two atomic ensembles using the same LO: the first ensemble provides the information to implement the phase feedback algorithm on the LO; the resulting corrected phase for the LO stays in the inversion region for a much longer period, and this prestabilized LO is used to interrogate the master ensemble with the standard Ramsey sequence. This scheme avoids the requirement of a trade-off between the number of intermediate measurements and the S/N of the final measurement by separating the two problems. The solution promises the same benefits foreseen for the phase reconstruction schemes proposed in [START_REF] Rosenband | Exponential scaling of clock stability with atom number[END_REF][START_REF] Kessler | Heisenberg-limited atom clocks based on entangled qubits[END_REF], but using only a single additional ensemble.

In an atomic clock the phase lock between the LO and the atomic superposition state can reduce the Dick effect [START_REF] Dick | Local oscillator induced degradation of medium-term stability in passive atomic frequency standards[END_REF], i.e. the aliasing of the clock oscillator, thanks to the longer interrogation time. However, the most important advantage of the scheme is the reduction of the decoher-ence related to the local oscillator, which translates to a lower white noise frequency for a fixed detection noise. This can serve to lower the LO stability requirements to the benefit of other parameters, like portability of the experimental setup, or viceversa, to remove the limitation set by the LO to reach ultimate performances. In the latter case, other effects will limit the clock interrogation time. Practically, a first limit is set by the vacuum quality, which can reduce the ensemble coherence via background collisions, but extremely long trapping lifetimes have been already reached for trapped atoms [START_REF] Willems | Creating long-lived neutral-atom traps in a cryogenic environment[END_REF]. Ultimately, with the combination of existing techniques and the method presented in this paper, the 1-10 mHz excited states lifetime expected for alkaline earthlike atoms could be reached, which motivates the search for transitions with a lower linewidth [START_REF] Campbell | Single-ion nuclear clock for metrology at the 19th decimal place[END_REF][START_REF] Derevianko | Highly charged ions as a basis of optical atomic clock-work of exceptional accuracy[END_REF][START_REF] Kozlov | Prospects of building optical atomic clocks using ErI and ErIII[END_REF].

Phase locking the LO to the atomic state preserves classical correlations in time against the decoherence by the local oscillator. The technique and the related enhancement factor could thus be combined to spin squeezing, which improves the clock sensitivity by introducing quantum correlations between the particles to go below the standard quantum limit [START_REF] Leroux | Implementation of cavity squeezing of a collective atomic spin[END_REF][START_REF] Louchet-Chauvet | Entanglement-assisted atomic clock beyond the projection noise limit[END_REF]. More generally, increasing the interrogation time using minimally destructive measurements and feedback on the phase could be applied to other atomic interferometers, such as atomic inertial sensors, where for example the phase of the interrogation lasers could be locked to the phase evolution of matter waves.

In conclusion, using a coherence preserving detection we tracked the phase evolution of an atomic collective superposition, and we reproduced it on a classical replica by introducing feedback to implement a PLL. Phase locking the LO to the atomic superposition state can be a key technology to improve the sensitivity of systems where the atomic phase evolution is compared to that of a classical LO, such as in atomic clocks and inertial sensors, whenever the main decoherence source is determined by the classical subsystem. This development may open new directions and possibilities in several technological fields and as well in basic science.

based on a Spectra Dynamics DLR-100 system as a frequency reference. The DLR-100 relies on an ultra-low noise 100 MHz quartz, locked at low frequency to the 10 th harmonic of a frequency doubled 5 MHz quartz to further improve the phase noise. The 100 MHz signal is multiplied to 7 GHz and then mixed with a tunable synthesizer at 165 MHz to obtain the signal resonant with the transition between the | ↓ and | ↑ state. Frequency noise is added to the LO signal using a frequency modulation port on the synthesizer, with a conversion factor set to 200 Hz/V rms . The noise signal for the demonstration of the clock using the PLL sequence is generated with a signal generator, which produces white frequency noise with a spectral density of 2.7×10 -2 Hz 2 /Hz; this signal is low pass filtered at 1.85 kHz before being added to the LO. This results in a rms phase drift of 430 mrad over 10 ms for the LO.

Appendix D: Feedback controller

The atomic populations on the | ↓ and | ↑ states determine a differential phase shift of the probe sidebands. This results in an amplitude modulation of the beam, which is detected by a photodiode (1591NF, New Focus), amplified (two HMC716LP3E, Hittite) and demodulated (ZX05-73C-C+, Minicircuits). The electronic integration of this signal during the 60 ns interrogation pulse produces a voltage proportional to the average atomic population difference during the probing. Such a voltage is digitized by a 14 bit resolution analog-to-digital converter embedded in the microprocessor (MCU) used to control the feedback loop (ADuC841, Analog Devices). In the experiments implementing the phase lock between the LO and the atomic superposition the MCU controls the phase actuator, a 6 bits step phase shifter with a range of 2π (RFPSHT0204N6, RF-Lambda). The total delay for the feedback is approximately 150 µs depending on the calculation time of the MCU. When running the clock based on the PLL technique, the MCU acts as well on a frequency actuator, which is the frequency modulation input on the 165 MHz synthesizer.

The feedback controller we propose for the clock exploiting the phase lock between the LO and the atomic superposition state consists of a cycle divided in three main steps: in the first one, successive correlated interrogations with probe interval T are realized on the same coherent atomic ensemble, and feedback on the LO phase is applied. The control law is

ϕ (i) LO = ϕ (i-1) LO + ϕ (i) FB = ϕ (i-1) LO + g ϕ ϕ (i)
where ϕ (i) is the estimated phase difference between the LO and the atoms at the i-th cycle and ϕ (i) LO is the phase of the LO only. The values of ϕ (i) FB = g ϕ ϕ (i) are saved in the feedback controller and typically a gain g ϕ = -1 is chosen.

The second step consists in the final phase readout: at the end of the N-th interrogation interval of duration T, a destructive measurement ϕ f of the phase is performed, with the highest possible precision. The saved phase shifts on the LO and ϕ f are used to reconstruct the full phase drift between the LO and the atoms in the total interrogation time, equal to T tot =N×T.

In the third step, one then performs feedback on the frequency as in a conventional atomic clock The cycle then repeats. The important feature of the feedback controller is that the feedback actions on the LO oscillator phase during the interrogation time are saved. They are then used with the output of the final precise measurement to determine the total phase drift. There is no drawback from the uncertainty of the weak measurements, since any feedback errors are detected with the precise final measurement at the end. As already remarked, in the experimental demonstration of the clock operation, we adopted the same probe for the intermediate and the final measurements.

FIG. 2 :

 2 FIG. 2:Bloch's sphere representation of the phase lock between the LO and the atomic superposition state. The phase lock between the classical oscillator and the atomic spin is obtained using repeated, time correlated Ramsey interrogations and feedback. The sequence begins by preparing the atomic CSS in the | ↓ state via optical pumping (step 1.). The measurement of the relative phase between the CSS and the LO starts when a π/2 microwave pulse around the y axis brings the CSS into a balanced superposition of the | ↓ and | ↑ states, depicted as a vector on the equatorial plane of the Bloch sphere (step 2.). The x axis is chosen to represent the phase of the local oscillator, and ϕ the relative phase between the LO and the atomic superposition that evolves because of the LO noise (step 3.). After an interrogation time T, the projection of ϕ is mapped onto a population difference by a projection π/2 pulse around the x axis and read out with the coherence preserving detection (step 4.). The CSS is rotated back to the equatorial plane by a reintroduction π/2 pulse around the x axis (step 5.), and feedback is applied on the phase of the LO (step 6.). The PLL between the LO and the atomic ensemble consists in the repetition of the steps from 3.) to 6.), potentially till the atomic ensemble shows a residual coherence.

FIG. 3 :

 3 FIG.3: Coherence preserving measurement of the relative phase between the LO and the atomic superposition. Real-time measurement of the vertical spin projection to which the relative phase ϕ between the local oscillator and the atomic superposition state is periodically mapped via microwave rotation pulses. The phase precession is induced by setting the LO frequency 100 Hz off the nominal atomic transition frequency. The experimental points are fitted with a sinusoidal evolution, damped by the effect of the Jz measurement and by the residual differential light shift on the clock transition.

FIG. 6 :

 6 FIG.6: Protocol for atomic clock with PLL. Operation of the clock: at each cycle of duration TC, the relative phase is repeatedly measured in a coherence preserving way during the phase lock interval (above: shaded grey areas). Each interrogation, represented in the inset by a light red peak between the manipulation π/2 pulses, is followed by a phase correction ϕ

  cycle g ϕ ϕ (i)and n is the clock cycle. In addition, the phase shift set by the feedback controller on the LO is reset to zero ϕ

Acknowledgments

We acknowledge funding from DGA, CNES, EMRP (QESOCAS), the European Union (EU) (iSENSE), ANR (MINIATOM), LAPHIA (APLL-CLOCK, within ANR-10-IDEX-03-02) and ESF Euroquam. The EMRP is jointly funded by the EMRP participating countries within EURAMET and the European Union. LCFIO and SYRTE are members of the Institut Francilien de Recherche sur les Atomes Froids (IFRAF). We thank T. Vanderbruggen for his contribution during the initial phase of the experiment, as well as G. Santarelli and T. Udem for useful discussions. P.B. acknowledges support from a Chair of Excellence of Région Aquitaine, E.C. from Quantel. Any mention of commercial products does not constitute an endorsement.

Appendix A: Atomic sample preparation 87 Rb atoms laser cooled with a magneto-optical trap are transferred to an optical dipole trap at 1560 nm, that uses a 4 mirror optical resonator to enhance the laser intensity [START_REF] Bernon | Heterodyne nondemolition measurements on cold atomic samples: towards the preparation of non-classical states for atom interferometry[END_REF]. The atoms are trapped at the crossing of two cavity arms with a waist of 100 µm. The ensemble is evaporatively cooled by decreasing the intensity of the dipole trap till a temperature of 10 µK is reached for 5×10 5 atoms in a cloud with 1/e 2 radius of 50 µm. In the last operation before starting the Ramsey interrogation the atoms are optically pumped in the | ↓ state. The sequence to prepare the ensemble in the initial state, which corresponds to the dead time T d in the atomic clock sequence, lasts 1.9 s.

Appendix B: Non-destructive dispersive probe

The measurement of J z is based on the dispersion caused by the trapped atoms on a far off-resonance optical probe. The probe beam has a waist of 47 µm matched to the size of the atomic cloud. It is phase modulated at a frequency of 3.853 GHz and frequency referenced at 3.377 GHz on the red of the the F=1→F'=2 transition; these conditions produce a symmetric mixing of the | ↓ and | ↑ states because of probe induced spontaneous emission, thus avoiding a vertical offset when the Bloch sphere contracts. In this way, each sideband mainly probes the population of one of the two levels, with the same magnitude and opposite sign for the couplings. We cancel the probe induced light shift and the related decoherence by precisely compensating the effect of the carrier with that of the sidebands; this is obtained by setting a modulation depth of 14.8% for the phase modulation. The differential light-shift on the D2 line from the dipole trap at 1560 nm was compensated with light blue detuned from the 5 2 P 3/2 → 4 2 D 5/2,3/2 transitions at 1529 nm. When the total power of the probe is set to 480 µW, it causes the decay of the atomic coherence with a lifetime of 2.85 µs; in the experiment here reported the interrogation pulses, obtained using an amplitude electro-optic modulator, have been set to last 60 ns. This determines for each pulse a 2% destructivity of the ensemble coherence, and a S/N of 20 for the J z measurement on the initial sample of 5×10 5 atoms. The population imbalance read-outs have been normalized to the signal when all atoms were repumped to the state F=2.

Appendix C: Frequency chain The 6.835 GHz frequency used to coherently manipulate the atomic spin is generated by a frequency chain


G Latu

M Mehrenberger

M Ottaviani

E Sonnendrücker

Aligned interpolation and application to drift kinetic semi-Lagrangian simulations with oblique magnetic field in cylindrical geometry

We introduce field aligned interpolation for Semi-Lagrangian schemes, adapting a method developed by to the semi-Lagrangian context. This approach is validated on the constant oblique advection equation and on a 4D drift kinetic model with oblique magnetic field in cylindrical geometry. The strength of this method is that one can reduce the number of points in the longitudinal direction. More precisely, we observe that we gain a factor |n| |n+mι| (where ι is the inverse safety factor), with respect to the classical approach, for the typical function sin(mθ + nϕ).

Introduction

In gyrokinetic simulations, it is observed that solution structures follow the field lines of the (strong) magnetic field and numerical methods have to be adapted to benefit from this fact. Different strategies exist for dealing with field alignement in gyrokinetic codes (see [START_REF] Sánchez | Linear and nonlinear simulations using the EUTERPE gyro kinetic code[END_REF], [START_REF] Görler | The global version of the gyrokinetic turbulence code GENE[END_REF] for example). We explore here an idea developed recently in [START_REF] Hariri | A flux-coordinate independent field-aligned approach to plasma turbulence simulations[END_REF] and adapt it in the context of a semi-Lagrangian code.

Our example of validation will be the following 4D drift-kinetic equation in cylindrical geometry, with oblique magnetic field. We look for f = f (t, r, θ, z, v) satisfying

∂ t f + [φ, f ] + v∇ f -∇ φ∂ v f = 0, with [φ, f ] = - ∂ θ φ rB 0 ∂ r f + ∂ r φ rB 0 ∂ θ f, ∇ = b • ∇, so that ∂ t f - ∂ θ φ rB 0 ∂ r f + ∂ r φ rB 0 + v b θ r ∂ θ f + vb z ∂ z f -b θ ∂ θ φ r + b z ∂ z φ ∂ v f = 0, (1)
for (r, θ, z, v)

∈ [r min , r max ] × [0, 2π] × [0, 2πR] × [-v max , v max ]
. The self-consistent potential φ = φ(r, θ, z) solves the quasi-neutral equation without zonal flow

-∂ 2 r φ + 1 r + ∂ r n 0 n 0 ∂ r φ + 1 r 2 ∂ 2 θ φ + 1 T e φ = 1 n 0 f -f eq dv .
Here the oblique magnetic field B whose norm is B (which can depend on r) writes

B = Bb, b = b z ẑ + b θ θ, b θ = c √ 1 + c 2 , b z = 1 √ 1 + c 2 , c = ιr R ,
and is parametrized by B 0 := Bb z and the rotational transform iota which satisfies

ι = b θ /r b z /R = 1 q ,
where q is called the safety factor. When ι = 0, we get the classical drift kinetic model given in [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Crouseilles | A new fully two-dimensional conservative semi-Lagrangian method: applications on polar grids, from diocotron instability to ITG turbulence[END_REF] for example. A similar model has been simulated in [START_REF] Sánchez | Linear and nonlinear simulations using the EUTERPE gyro kinetic code[END_REF], with ι = 0.8 as an example, using a Particle in Cell method. Equation (1) can be derived from the gyrokinetic equations. Note that some terms are dropped (w.r.t. [START_REF] Sánchez | Linear and nonlinear simulations using the EUTERPE gyro kinetic code[END_REF] for example), which permits to retain the oblique feature while sticking to a same structure of equation as in the case ι = 0.

The strategy that we adopt is to solve the constant oblique advection equation

∂ t f + vb • ∇f = 0, (2)
that enters in equation ( 1) using an interpolation that is aligned along the direction of the magnetic field b. Such strategy permits to reduce the number of points in the z-direction. Note that another possible strategy would be to adapt the grid to magnetic field lines (see [START_REF] Braeunig | Guiding-center simulations on curvilinear meshes[END_REF] for example); one strength of the present approach is that the grid points do not need to be changed, which permits an easier implementation, as for example the Poisson solver (here for the quasi-neutral equation) does not need to be changed.

The numerical scheme is developed in Section 2, and numerical results are shown in Section 3: first, results on the 2D oblique constant advection equation, where numerical errors can be performed and then results on the 4D drift kinetic model with oblique magnetic field. In Section 4, we give the conclusion and the perspectives of this work. In Appendix A, we detail the derivation of (1), from the gyrokinetic equations and give the dispersion relation in Appendix B.

Numerical scheme 2.1 Constant oblique advection

Writing ϕ = z R , we have to solve for g := g(t, θ, ϕ) = f (t, r, θ, Rϕ, v), the constant oblique advection equation [START_REF] Braeunig | Guiding-center simulations on curvilinear meshes[END_REF]. We have

∂ t f + v∇ f = ∂ t f + v b θ r ∂ θ f + vb z ∂ z f = 0, which leads to ∂ t g + v b θ r ∂ θ g + v b z R ∂ ϕ g = ∂ t g + ṽ (ι∂ θ g + ∂ ϕ g) = 0, (3)
with ṽ = vbz R , and g is 2π periodic in θ and ϕ. Let ∆t ∈ R + , and t j = ∆t, ∈ N. We have the relation

g(t + ∆t, θ, ϕ) = g(t , θ -ιṽ∆t, ϕ -ṽ∆t).
Let N θ , N ϕ ∈ N * , and

θ i = 2πi N θ , ϕ j = 2πj
Nϕ , which can be defined for i, j ∈ R. We suppose to know values g ,i,j g(t , θ i , ϕ j ), for i = 0, . . . , N θ -1, j = 0, . . . , N ϕ -1. By periodicity, we can suppose that i, j ∈ Z. We fix two integers r ≤ 0 ≤ s. For j = 0, . . . , N ϕ -1, there exists j 0 ∈ Z and 0 ≤ β < 1 such that ϕ j -ṽ∆t = ϕ j0+β .

We then define ϕ j -ṽ∆t k = ϕ j0+k , k = r, . . . , s. For each k = r, . . . , s, from the values

g ,i,j0+k g(t , θ i , ϕ j -ṽ∆t k ), i = 0, . . . , N θ -1, we compute gi,k g(t , θ i -ιṽ∆t k , ϕ j -ṽ∆t k ) = g(t + ∆t k , θ i , ϕ j ), i = 0, . . . , N θ -1,
by using an interpolation in θ.

For each i = 0, . . . , N θ -1, from the values gi,k g(t , θ i -ιṽ∆t k , ϕ j -ṽ∆t k ), k = r, . . . , s we finally compute g +1,i,j g(t , θ i -ιṽ∆t, ϕ j -ṽ∆t), using an interpolation along the parallel direction: we reconstruct a value gi,β , from the values gi,k , k = r, . . . , s and take g +1,i,j = gi,β .

In the following, we will use Lagrange of degree 2d + 1 LAG(2d+1), for the interpolation in the parallel direction and take r = -d, s = d + 1.

Drift kinetic model with oblique magnetic field

We use a classical backward semi-lagrangian (BSL) scheme as in the case where ι = 0 [START_REF] Crouseilles | A new fully two-dimensional conservative semi-Lagrangian method: applications on polar grids, from diocotron instability to ITG turbulence[END_REF]. The model is implemented in SELALIB [10] and uses a parallelization in r. Advection in z is replaced by an advection along the parallel direction [START_REF] Braeunig | Guiding-center simulations on curvilinear meshes[END_REF]. The term b • ∇φ is computed along the parallel direction, using a finite difference formula.

3 Numerical results

Constant oblique advection

In the case of constant oblique advection, we have to solve [START_REF] Crouseilles | A new fully two-dimensional conservative semi-Lagrangian method: applications on polar grids, from diocotron instability to ITG turbulence[END_REF].

We consider an initial function with a well defined helicity g = g 0 (mθ + nϕ), so that

b • ∇f = m b θ r + n b z R g 0 (mθ + nϕ) = k g 0 (mθ + nϕ),
where

k := b z R (n + ιm) = b z qR (m + qn) ,
In order to have b • ∇f bounded, we look for situations where

|m + qn| ≤ 1,
as in real tokamaks, it is assumed that k will typically be in the range of [-1 q R , 1 q R ]. We will use in the sequel g 0 = sin. The displacement due to advection equation has a main parameter ∆t. A extended set of various ∆t will be investigated because the error of one numerical scheme is much dependant on it. We choose a value for safety factor which is a non-rational surface case, q = √ 2. We look at four configurations for the initial functions

A : (n = 5, m = -7) k ≈ 0.07 1 q R B : (n = 24, m = -34) k ≈ -0.06 1 q R C : (n = 5, m = -6) k ≈ 1.07 1 q R D : (n = 24, m = -33) k ≈ 0.94 1 q R (4)
In the Figures 2 and3, the abcissa of the plots are N ϕ , and the ordinate are the L ∞ norm which is the maximum of the difference of the function computed after one time step versus the analytical function which is here

f (θ, ϕ, ∆t) = sin(m (θ -∆t) + n (ϕ -q ∆t)) .
The maximum is over the grid points and over several values of ∆t. We show the results of the aligned versus the standard (non-aligned) scheme. Some parameters are fixed for this study:

N θ = 400, q = √ 2.
Note that several time steps ∆t are evaluated to approximate well the maximal error one can reach for the set of parameter used (the error is varying much along with the time step).

For left-hand side plots on these two Figures, the k is near 0 and the aligned method is very accurate, even f one takes a low value for N ϕ . If one considers now the right-hand side plots, the k is close to 1 qR . Even if the aligned method gives lower error than standard method, the error is bigger at small N ϕ compared to low k .

The Fig. 2 considers lower frequency than Fig. 3. The behaviour in the two cases are similar except that there is a shift of the curves along the ϕ direction. This is expected, and what is interesting is that the aligned method behaves good (low error) even with low values for N ϕ .

In Fig. 4, the discretisation along θ has been refined (N θ = 2000) compared to Fig. 3. The asymptotic error (at large N ϕ ) is lowered. The aligned method for k close to zero (left-hand side plot) is very accurate with an L ∞ error lower than 10 -6 . A more detailed study is given on Figures 5,6,7. We first use the same configurations, but make vary the degree of interpolation in the Lagrange reconstruction: LAG3, on Figure 5 top), LAG5 on Figure 5 bottom and LAG17 on Figure 6. For the standard method, we use Lagrange interpolation in ϕ direction and stick to cubic splines in the θ direction, as in the case of oblique interpolation. We give the L ∞ error w.r.t. N φ /n for the standard method and w.r.t. N φ /| k | for the aligned method, where we have defined k = n + mι.

Note that k and k only differ by the constant multiplicative factor bz R , which enters as factor in the advection. The error does not here really depend on this factor, as we look here for maximal error over one time step of different size. In the following, we have taken ∆t ∈ {0.1/s}, s = 1, . . . , 100, and solved the equation on [0, 2π] 2 with A 1 = 1 and A 2 = q. So r, R, v can be chosen so that ṽ = q, that is

∂ t g + A 1 ∂ θ g + A 2 ∂ ϕ g,
vb z qR = 1, b z = 1 + r qR 2 .
In this way of presentation, we remark that we have three zones: a first zone, when the error does not decrease, a second zone where the error decreases according to the order of interpolation and then a third zone when the error does no more decrease which depends on N θ /|m|. Note that for a fixed degree and a given N θ /|m|, the error corresponding to the standard or aligned method lies on the same curve. But for a given N φ , there is shift in abscissa which corresponds to the distance between N ϕ /|n| and N ϕ /| k |. So when the error of the standard method is in the left of the curve, the error of the new method lie on the right of it. So, for example, in order to have an error around the error of discretization in θ (when the error begins to saturate: beginning of the third zone), as we are on the same curve, we typically need -1) and the aligned method becomes more and more efficient. For k 0, we have the best result (this corresponds to the previous curve: with n = 24). Increasing again n, the results are still better but less and less (n = 30 and k 6). For n = 100 (not shown), we approach again the curve corresponding to n = 12. We have considered here LAG9, N θ = 200, except for n = 23, where N θ = 400. This to see, that the saturation error diminishes with N θ (on the previous plots, we had seen that the error increases with m).

On Figures 5,6,7, we have considered the following values for N ϕ :

N ϕ ∈ A, A = {2,

Drift kinetic model with oblique magnetic field

The initial function is given by

f (t = 0, r, θ, z, v) = f eq (r, v) 1 + exp - (r -r p ) 2 δr cos mθ + n R z ,
where the equilibrium function is

f eq (r, v) = n 0 (r) exp(-v 2 2Ti(r) ) (2πT i (r)) 1/2 .
The radial profiles {T i , T e , n 0 } have the analytical expressions

P(r) = C P exp -κ P δr P tanh( r -r p δr P ) , P ∈ {T i , T e , n 0 },
where the constants are

C Ti = C Te = 1, C n0 = r max -r min rmax rmin exp(-κ n0 δr n0 tanh( r-rp δrn 0
))dr .

Finally, we consider the parameters of [START_REF] Coulette | Numerical comparisons of gyrokinetic multi-water-bag models[END_REF] (MEDIUM case) We take B 0 = -1. We consider the case ι = 0, n = 1, m = 15, which leads to k = 1 R . We then consider a case ι = 0.8, n = -11, m = 15, which leads to k = bz R (-11 + 0.8

r
• 15) = bz R . Note that b z = 1 √ 1+c 2 , with 0 ≤ c ≤ ι rmax R ≤ 0.05, so that |b z -1| ≤ 1.25 • 10 -3
, and thus the dispersion relation which depend on m and k and not directly on n (see Appendix A) will give almost the same result, which means that both simulations should lead to similar results in the poloidal plane, at least in the linear phase (as it is also observed in [START_REF] Sánchez | Linear and nonlinear simulations using the EUTERPE gyro kinetic code[END_REF]). We take LAG5 for the interpolation along the parallel direction and cubic splines for the interpolation along θ. Finite differences of order 6 are used for the derivative computation along the parallel direction and cubic splines are used otherwise. When ι = 0, we use the classical method with cubic splines for the interpolation along z direction. Such behavior is observed on Figure 8 and Figure 9. We see that the poloidal cut f (t, r, θ, z = 0, v = 0) are similar in the linear phase (Figure 8 left) and the corresponding excited modes are clearly visible in the θ -z cut f (t, r p , θ, z, 0) (Figure 9 right). After a while, we see the effect of the non linear phase. Note that we have not excited the most unstable mode (which is here m = 10, k = 3 R ).We still observe an alignment of the structures (Figure 9 bottom middle/right), and we see that the poloidal cuts are very similar considering N z = 32 or N z = 64. Note that in these figures we use raw data for the visualisation. Since the number of points in N z is purposely low, the corresponding plots in the θ -z plane (Figure 8 middle/right and Figure 9 bottom/right) are necessarily coarse. This is not an indication of numerical problems. Indeed a better visualisation in this plane can be achieved by reconstructing the distribution function on a finer mesh using the field aligned interpolation.

Conclusion and perspectives

We have given some first numerical evidence that the strategy of Hariri-Ottaviani [START_REF] Hariri | A flux-coordinate independent field-aligned approach to plasma turbulence simulations[END_REF] works in the context of semi-Lagrangian gyrokinetic simulations. Validation has been performed on the analytical test of the constant oblique advection and then on a drift kinetic equation with oblique magnetic field. Extension to tokamak configuration in toroidal geometry is the next step of this study. Just before submitting the first preprint, we got knowledge of the paper [START_REF] Kwon | Development of semi-Lagrangian gyrokinetic code for full-f turbulence simulation in general tokamak geometry[END_REF]; there, first results in this direction are already given.

Appendix B: Dispersion equation

We make the following expansions:

f = f 0 + εf 1 + O(ε 2 ), φ = φ 0 + εφ 1 + O(ε 2 ) with f 0 (r, v) = f eq (r, v) = n 0 (r) exp -v 2 2Ti(r) (2πT i (r)) 1/2
, φ 0 = 0.

We obtain

∂ t f 1 - ∂ θ φ 1 rB 0 ∂ r f 0 + vb z ∂ z f 1 + v b θ r ∂ θ f 1 -b θ ∂ θ φ 1 r + b z ∂ z φ 1 ∂ v f 0 = O(ε).
and

-∂ 2 r φ 1 + 1 r + ∂ r n 0 n 0 ∂ r φ 1 + 1 r 2 ∂ 2 θ φ 1 + 1 T e φ 1 = 1 n 0 f 1 dv + O(ε).
We assume that the solutions have the form :

f 1 = f m,n,ω (r, v)e i(mθ+kz-ωt) , φ 1 = φ m,n,ω (r)e i(mθ+kz-ωt)
with k = n R . Then, we obtain

(-ω + kvb z + v mb θ r )f m,n,ω = m rB 0 ∂ r f 0 + b θ m r + b z k ∂ v f 0 φ m,n,ω and
-∂ 2 r φ m,n,ω + 1 r + ∂ r n 0 n 0 ∂ r φ m,n,ω - m 2 r 2 φ m,n,ω + 1 T e φ m,n,ω = 1 n 0 f m,n,ω dv, We get, as k = b θ m r + b z k -∂ 2 r φ m,n,ω + 1 r + ∂ r n 0 n 0 ∂ r φ m,n,ω - m 2 r 2 φ m,n,ω + 1 T e φ m,n,ω = 1 n 0 φ m,n,ω m rB0 ∂ r f 0 + k ∂ v f 0 vk -ω dv
By using the expression of f 0 , we have

I = -v Ti + m k rB0 ( ∂rn0 n0 -∂rTi 2Ti + v 2 ∂rTi 2T 2 i ) v -ω k f 0 dv.
We introduce for ∈ N :

I (u) = 1 n 0 v f 0 v -u f 0 dv, so that I n 0 = - 1 T i I 1 ω k + m k rB 0 ∂ r n 0 n 0 - ∂ r T i 2T i I 0 ω k + ∂ r T i 2T 2 i I 2 ω k .
We use the relations :

I 0 = 1 (2T i ) 1/2 Z u (2T i ) 1/2 , I 1 = 1 + uI 0 , I 2 = u (1 + uI 0 ) , with Z(u) = 1 √ π exp(-x 2 ) x -u dx = i √ π exp(-u 2 )(1 -erf(-iu)), erf(x) = 2 √ π x 0 exp(-t 2 )dt.
The dispersion relation is, putting φ = φ m,n,ω , for convenience,

A = -∂ 2 r φ - 1 r + ∂ r n 0 n 0 ∂ r φ + m 2 r 2 φ + 1 T e φ = - 1 T i (1 + zZ(z)) + m k * rB 0 Z(z) ∂ r n 0 n 0 - ∂ r T i 2T i + z(1 + zZ(z)) ∂ r T i T i φ,

Figure 1 :

 1 Figure 1: Schematic view of the oblique interpolation. Values at green points at time t +1 are updated through values at red points at time t . These values are obtained by computing first values at blue points, using an interpolation in θ (vertical direction) from values at black points (constant advection), at time t . Then values at red points are obtained from values at blue points using Lagrange interpolation along the oblique parallel direction (here LAG5).

Figure 2 :Figure 3 :

 23 Figure 2: Error in L ∞ -norm compared to the analytical solution for advection N θ = 400, q = √ 2, n = 5 and m = -7 (left), m = -6 (right)

Figure 4 :

 4 Figure 4: Error in L ∞ -norm compared to the analytical solution for advection N θ = 2000, q = √ 2, n = 24 and m = -34 (left), m = -33 (right)

Nϕ|

  k | points for the aligned method instead of Nϕ |n| points. The factor of gain is thus |n| | k | = |n| |n + mι| On Figure 7, we fix m = -34 and change the values of n. Now we consider as abscissa only N ϕ /|m|. We then see the same effect. We can remark that when n is smaller than | k | (here n = 5 and k -19), the classical method is more accurate. For n and k similar, both methods have the same accuracy (n = 12 and k -12). Then increasing n, | k | is diminished (n = 23 and k

  min = 0.1, r max = 14.5, κ n0 = 0.055, κ Ti = κ Te = 0.27586, δr Ti = δr Te = δr n0 2 = 1.45, = 10 -6 , R = 239.8081535, r p = r min + r max 2 , δr = 4δr n0 δr Ti .

  with k * = (2T i ) 1/2 k , and z = ω k * , and recalling that k = b θ m r + b z k . Note that the dispersion relation depends on m and k and not directly on n. This means that taking different values of ι and n but with same m and k will lead to the same dispersion relation.

Figure 5 :

 5 Figure 5: L ∞ error vs ϕ /n the standard method (old) and vs N ϕ /|k | for the aligned method (new). Parameters are q = √ 2, N θ = 400. Top: LAG3; bottom: LAG5.

Figure 6 :

 6 Figure 6: L ∞ error vs N ϕ /n for the standard method (old) and vs N ϕ /|k | for the aligned method (new). Parameters are q = √ 2, N θ = 400 and LAG17.

Figure 7 :

 7 Figure 7: L ∞ error vs N ϕ /n using N θ = 200 (or N θ = 200, when n = 23), with parameters q = √ 2, m = 34, LAG9 and n = 5, 12, 23, 30 which corresponds to k -19, -12, -1, 6

Figure 8 :

 8 Figure 8: Poloidal cut (left) and θ -z cut (right) of distribution function at time t = 4000. Top: ι = 0, n = 1, m = 15 on 256 × 512 × 32 × 128 grid. Middle: ι = 0.8, n = -11, m = 15 on 256 × 512 × 32 × 128 grid. Bottom: ι = 0.8, n = -11, m = 15 on 256 × 512 × 64 × 128 grid.

Figure 9 :

 9 Figure 9: Poloidal cut (left) and θ -z cut (right) of distribution function at time t = 6000. Top: ι = 0.8, n = -11, m = 15 on 256 × 512 × 32 × 128 grid. Bottom: ι = 0.8, n = -11, m = 15 on 256 × 512 × 64 × 128 grid.

where φ is the gyro-average operator applied to the electrostatic potential φ.

We then have, supposing that ι = ι(r),

We consider for the moment µ = 0, so that φ = φ. We then get

We then consider the following intermediate model

which we further simplify into

As we have the relation B = B 0 (1 + c 2 ) 1/2 , we get the following model

which corresponds to (1), by writing v instead of v, in order to use a shorter notation.