Didier Chauveau

Nicolas P A Saby

Thomas G Orton

Blandine Lemercier

Christian Walter

Dominique Arrouays

Large-scale simultaneous hypothesis testing in monitoring carbon content from French soil database -A semi-parametric mixture approach

Keywords: Carbon content, EM algorithms, FDR, finite mixture, France, semi-parametric mixtures, soil monitoring

Investigating the information of the French National Soil Tests database for soil monitoring produces multiple hypothesis testing problems with hundreds or thousands of test responses to consider simultaneously. A largely used concept of error control in such multiple testing is the expected proportion of falsely rejected hypotheses, or False Discovery Rate (FDR). A related notion of local FDR (ℓFDR) can be appropriately represented by considering that the observed p-values come from a twocomponents mixture model where the component corresponding to the null hypothesis is known. In this work, we explore different solutions for FDR estimation. In particular, we introduce a specific version of a semi-parametric Expectation-Maximization (EM) algorithm for ℓFDR estimation, and compare it to classical ℓFDR estimation using parametric mixtures, and conventional FDR approaches. The performances of the different models for estimating the FDR and related criteria are first illustrated on the results of simulated multiple comparison tests. These approaches are then applied to soil carbon content monitoring on our database. The results show that not taking into account the FDR estimation can lead to over-estimation of the number of cantons (locations) subject to a significant change. However, we have detected large numbers of significant changes in the database that occured during the time period of this study. Globally, losses in organic 1 carbon are observed in Northern France, along the Atlantic coastal regions, and to a less extend for the data collected over the North-Eastern regions. The OC increases are more scattered over the territory. We also use the data to estimate the minimum number of samples needed at each period to detect a given change.

Introduction

Political awareness that soil is threatened by increasing pressures has been rising for several years (European Commission, 2002). Indeed, the demand for soil information is increasing continuously (Richer de Forges and [START_REF] Richer De Forges | Analysis of requests for information and data from a national soil data centre in france[END_REF]. Although rates of soil degradation are often slow and only detectable over long timescales, they are often irreversible. Therefore, monitoring soil quality is essential in order to detect adverse changes in their status at an early stage [START_REF] Morvan | Soil monitoring in Europe: A review of existing systems and requirements for harmonisation[END_REF]. There are many generic issues that must addressed when establishing and operating such monitoring systems. Of particular importance is the requirement to detect change in soil over relevant spatial and temporal scales with adequate precision and statistical power [START_REF] Arrouays | Environmental Assessment of Soil for Monitoring: Volume IIa Inventory & Monitoring[END_REF][START_REF] Desaules | Advances in long-term soil-pollution monitoring of Switzerland[END_REF][START_REF] Van Wesemael | How can soil monitoring networks be used to improve predictions of organic carbon pool dynamics and CO(2) fluxes in agricultural soils?[END_REF][START_REF] Arrouays | Generic issues on broad-scale soil monitoring schemes: A review[END_REF]. Analysing the results from existing soil measurement exercises, such as operational soil testing by farmers or fertiliser suppliers is one potential option for detecting large temporal trends in soil characteristics. Therefore, legacy soil data testing results have often been used to assess temporal changes at national and regional levels (eg., for phosphorus by [START_REF] Skinner | Twenty-five years of monitoring pH and nutrient status of soils in England and Wales[END_REF] in England and Wales, [START_REF] Cahoon | Spatial and temporal variability in excessive soil phosphorus levels in eastern North Carolina[END_REF] in eastern North Carolina (USA), [START_REF] Wheeler | Trends in some soil test data over a 14-year period in New Zealand[END_REF] in New Zealand, [START_REF] Lemercier | Soil phosphorus monitoring at the regional level by means of a soil test database[END_REF] in France, [START_REF] Reijneveld | Changes in soil phosphorus status of agricultural land in the Netherlands during the 20th century[END_REF] in The Netherlands; and for carbon, [START_REF] Saby | Changes in soil organic carbon in a mountainous French region, 1990-2004[END_REF] in France; [START_REF] Reijneveld | Soil organic carbon contents of agricultural land in the Netherlands between 1984 and 2004[END_REF] in The Netherlands). However, the conclusions drawn using these kinds of data may be subject to several sources of bias that are inherent in a non-controlled sampling strategy.

Since the exact locations of the carbon content observations are not known, the only geographical information available is the name of the town of the sampled location. A simple statistical framework consists then in gathering the observations per group of commune (called canton) and time period to perform inference tests. In this framework, the null hypothesis could be either that the means or the medians of the two groups are equal. We assume here that the locations are distant enough, so that measurements can be considered as independent (i.e. we do not consider a spatial correlation in the present work). Hence, investigating the information of the French National Soil Tests database for soil monitoring produces multiple hypothesis testing problems with hundreds or thousands of test responses to consider simultaneously. In this multiple inference setup, the unguarded use of single-inference procedures results in an increased false positive rate among the simultaneous tests. The consideration of a global type I error level (say α = 5%) when several tests are performed simultaneously, has been considered first in the well-known context of analysis of variance (Anova), where several levels of a factor are compared two-by-two, leading to n simultaneous tests. It is easy to check that, if n independent tests with level of significance α are applied simultaneously, the achieved FamilyWise Error Rate (FWER), that is the probability of observing at least one false rejection among the n tests, is 1 -(1α) n which quickly increases with n and is already ≈ 99% for n = 100. The historical approach since the early 1950s, called the Bonferroni approximation procedure (see e.g. [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] and references therein) consists in applying each test at a level α/n, resulting in a FWER lesser than α. However, Bonferronitype procedures appear to be too conservative when n gets large because α/n gets too small, leading to too few rejections (as will be illustrated in Section 4). Since then, the point of view on the problem has changed, focusing in the number (or ratio) of erroneous rejections instead of the question wether any error was made, as for the FWER. In this vein, the most popular and largely used concept of error control in such multiple testing inference is nowaday the expected proportion of falsely rejected hypotheses among the simultaneous tests, or False Discovery Rate (FDR, [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF].

Several statistical algorithms have been proposed in the literature for estimating the FDR, the recent and unified procedure based on a nonparametric approach from [START_REF] Strimmer | A unified approach to false discovery rate estimation[END_REF] appearing to be one of the current standard for practitioners. An alternative approach to FDR estimation consists in viewing the problem as the statistical estimation of the parameters of a finite mixture model. Parametric mixtures for FDR estimation have been considered in, e.g., [START_REF] Mclachlan | A simple implementation of a normal mixture approach to differential gene expression in multiclass microarrays[END_REF], and semiparametric versions have also been considered recently, such as in [START_REF] Robin | A semiparametric approach for mixture models: Application to local false discovery rate estimation[END_REF]. The alternative approach that we propose here belongs to that last framework, but takes advantages of very recent developments of statistical algorithms for estimating semi-or nonparametric mixtures. All these standard and new models will be detailed in Section 1.2. In particular, we implement and compare new ideas pertaining to recent semi-and non-parametric extensions of the well-known EM (Expectation-Maximization) algorithm commonly used for estimating parametric mixtures. Some specific ingredients coming from the multiple testing framework can be specified within the EM algorithm machinery.

When monitoring soil quality, the estimation and control of the FDR is of particular importance as false negative or false positive may lead to unwarranted recommendations or to irrelevant policies (see, e.g., [START_REF] Lark | Using expert knowledge with control of false discovery rate to select regressors for prediction of soil properties[END_REF]. Taking the example of organic carbon, such false positive or negative results may lead to errors in national carbon accounting, or to wrong identification of areas at risks where protection policies should be implemented. The basic motivations we have in this paper are thus twofold: (i) To improve analysis of soil carbon content monitoring, in particular by applying recent and new multiple testing procedures and FDR estimation and control, on datasets issued from the French National Soil Tests database. This application field is different than micro-array techniques for which the FDR is more commonly applied. For instance, the proportion of null hypothesis in the population is about 50%, smaller than what is usually observed in microarray setup (e.g., 80% in the gene experiment example studied in [START_REF] Robin | A semiparametric approach for mixture models: Application to local false discovery rate estimation[END_REF]). Our second objective is (ii) to compare recent methods for FDR estimation, such as the challenging nonparametric approach available in the fdrtool package [START_REF] Strimmer | A unified approach to false discovery rate estimation[END_REF], to approaches taking full advantage of an underlying mixture model assumed for the vector of p-values. Note that the alternative FDR control procedures that we develop here are implemented in the software package mixtools [START_REF] Benaglia | mixtools: An R package for analyzing finite mixture models[END_REF] for the R statistical software (R Development Core Team, 2010), and are available in the last public version.

The French National Soil Tests Database

The results of the soil analyses compiled in the French National Soil Test Database are a consequence of requests from farmers and landowners for help to improve the management of their crops and pastures. Samples were taken from topsoil horizons of cultivated or pasture fields, but the specific reasons individual farmers requested soil analysis are not known and therefore the sampling strategy could not be controlled. In each sampled field, 10-15 subsamples of the ploughed layer (or the 0 to 30 cm topsoil layer in the case of pasture) were collected with a hand auger and the extracted soil mixed to provide a composite sample. Samples were sent to laboratories certified by the French Ministry of Agriculture and the results from standardised analytical procedures were stored in the national database.

In the study area, 1,220,039 soil samples were analysed for OC content between 1990 and 2004. All samples were analysed by the wet chromic acid digestion method, utilising an excess amount of potassium dichromate with concentrated sulphuric acid at boiling-point. After the addition of concentrated phosphoric acid, the excess dichromate was back-titrated with ferrous ammonium sulphate. This method is known not to measure all the organic carbon present but only that which is easily oxidisable, and recovery factors have been recently studied by De [START_REF] De Vos | Walkley-black analysis of forest soil organic carbon: recovery, limitations and uncertainty[END_REF] who recommend the use of dry combustion. This was not possible because the study used historical data from samples which were not archived. The OC that is not recovered by wet chromic acid digestion is mainly composed of highly resistant compounds and thus our estimates of change are likely to be little affected by the use of this wet digestion. Moreover, in previous studies we showed that wet chromic acid digestion recovered on average 95% of the OC estimated by dry combustion on various French soils [START_REF] Jolivet | Comparison between analytical methods for organic carbon and organic matter determination in sandy spodosols of france[END_REF][START_REF] Arrouays | The carbon content of topsoil and its geographical distribution in france[END_REF]. The raw data are spatially aggregated at the cantonal level (in France, a canton is an administrative unit combining several municipalities). They are also aggregated per time period, 1990-1997 and 1998-2005. In previous work, [START_REF] Saby | Changes in soil organic carbon in a mountainous French region, 1990-2004[END_REF] applied a different approach regarding splitting the data over time by using 3 time periods. This was adapted to the smallest region under study, where more data were avaible than for the present Nationwide study. In this work, we then used only 2 periods to optimise the number of data per canton. We refer to these data as Carbon data in the following of this paper.

Mixture models in multiple testing FDR evaluation

In multiple testing, we observe a sample of p-values p = (p 1 , . . . , p n ), where each individual observation p i corresponds to the critical probability of the ith test, for which either the null hypothesis H 0 is true (not significant or not interesting), or the alternative hypothesis H 1 is true (interesting i.e. rejection of H 0 ). Since it is not observed whether each hypothesis is true or false, we are in the general framework of statistical inference from missing data: the critical probabilities can be viewed as a sample of observations in [0, 1] coming from a finite mixture model.

Finite mixtures give a flexible way to model a wide variety of random observations (see, e.g., [START_REF] Mclachlan | Finite mixture models[END_REF]. In such models, we assume that n independent measurements X 1 , . . . , X n are observed such that each X i comes from one of m possible sub-populations, called component distributions. Importantly, the component index, 1 though m, is not observed along with X i . It is common to associate to X i the unobserved indicator vector Z i = (Z i1 , . . . , Z im ) where Z ij = 1 if observation i comes from component j (hence Z i have (m -1) zero's and one 1). It is assumed that, unconditional on X i , each P(Z ij = 1) = λ j , the proportion of individuals from component j in the population, and that the probability density function (pdf) of an observation given that it belongs to component j is some density f j . Thus, the pdf of an observation from the mixture is g(x) = m j=1 λ j f j (x), where the λ j 's are strictly positive and sum to unity. The vector Z = (Z 1 , . . . , Z n ) represents the missing data associated to the n observations.

Until recently, only parametric mixture models have been considered and investigated, where parametric means that the f j 's come from a (generally single) parametric family indexed by some parameter ξ, i.e. f j (•) ≡ f ( • ; ξ j ), so that θ = (λ, ξ) = (λ 1 , . . . , λ m , ξ 1 , . . . , ξ m ) is the (finite-dimensional) model parameter to be estimated from the data. Normal mixtures, where f ( • ; ξ j ) is the density of the Gaussian N (µ j , σ 2 j ) with parameters ξ j = (µ j , σ 2 j ), are the most commonly used in this category. In multiple testing, we thus associate to the sample of p-values p the unobserved indicator variables Z where Z i will be the unobserved indicator of the true hypothesis associated with the ith test. For ease of notation we will index the components by {0, 1}, according to the hypotheses, i.e. f j is the pdf of the p-value conditional of H j being true,

Z i = (1, 0) if H 0 is true, and Z i = (0, 1) if H 1 is true, i.e. if H 0 is rejected (interesting ith test).
Hence, p i typically comes from a two-components mixture model with pdf

g θ (p) = λ 0 f 0 (p) + (1 -λ 0 )f 1 (p), (1)
where θ = (λ 0 , f 0 , f 1 ) is the model parameter. Theoretically, given that H 0 is true, the p-value is uniformly distributed over [0, 1],

(p i |"not interesting") = (p i |H 0 true) = (p i |Z i0 = 1) ∼ U [0,1] ,
whichever test is actually performed, so that f 0 (p) = I [0,1] (p) is known and the model parameter reduces to θ = (λ 0 , f 1 ), which makes an important difference between usual mixture model estimation and the present multiple tests setting. The density of (p i |Z i1 = 1) ∼ f 1 may be parametrically specified, as in usual mixture modelling. In the multiple testing setup for instance, Beta distributions have been considered, leading to the Beta-Uniform mixture model (see, e.g. [START_REF] Robin | A semiparametric approach for mixture models: Application to local false discovery rate estimation[END_REF]. However, since often in practice very little is known about the distribution under the alternative hypothesis, we propose here to relax the parametric assumption on f 1 .

Mixing parametric and nonparametric estimates in the FDR context is not new (see Robin et al., 2007, and references therein). In other contexts, extensions of mixture models where the component densities (f j 's) are semior non-parametrically specified have been considered quite recently in the statistical literature [START_REF] Hall | Nonparametric estimation of component distributions in a multivariate mixture[END_REF][START_REF] Hall | Nonparametric inference in multivariate mixtures[END_REF]. Semi-and non-parametric models and associated EM-like estimation algorithms have recently been proposed, in [START_REF] Bordes | A stochastic EM algorithm for a semiparametric mixture model[END_REF] for the univariate case, Benaglia et al. (2009a) and [START_REF] Levine | Maximum smoothed likelihood for multivariate maximum[END_REF] for multivariate and fully nonparametric situations. We extend here these algorithmic ideas to the multiple testing estimation problem. For such a mixture with one component known, the theoretical question of statistical identifiability has been adressed in [START_REF] Bordes | Semiparametric estimation of a two-component mixture model where one component is known[END_REF], where it is proved that identifiability is satisfied at least if f 1 is semi-parametrically specified, as will be detailed later on.

The two-components mixture model given by equation ( 1) is assumed, e.g., in [START_REF] Strimmer | A unified approach to false discovery rate estimation[END_REF] and [START_REF] Robin | A semiparametric approach for mixture models: Application to local false discovery rate estimation[END_REF]. However, one difference with our approach is that the fdrtool estimation method in [START_REF] Strimmer | A unified approach to false discovery rate estimation[END_REF] is grounded on an empirical (nonparametric) estimate of the model cumulative density function (cdf). It does not rely on the missing data aspect induced by the mixture, unlike any EM-based strategy does. In particular, an EM-like algorithm delivers, together with estimates of the mixture parameters (in this case (λ 0 , f 1 )), estimates of the posterior probabilities that each p-value comes from each component. In the present FDR context we can estimate directly the so-called local FDR (ℓFDR, [START_REF] Efron | Empirical bayes analysis of a microarray experiment[END_REF] from these posteriors ℓFDR(p i ) = P("not interesting"|p i ) = P(Z i0 = 1|p i ), for p i 's leading to rejection of H 0 . The FDR can then be computed as the average of the ℓFDR(p i )'s over all the rejected p i 's. [START_REF] Robin | A semiparametric approach for mixture models: Application to local false discovery rate estimation[END_REF] approach is very close to the semi-parametric EM approach we present here, except that they estimate λ 0 separately, and then estimate f 1 using a weighted kernel density estimate (wkde). It appears that this wkde itself is much in the spirit of the density estimation step of our semi-parametric EM-like strategy, following Benaglia et al. (2009a).

Another difference is that these recent FDR control approaches have been compared, e.g. on [START_REF] Robin | A semiparametric approach for mixture models: Application to local false discovery rate estimation[END_REF] and [START_REF] Strimmer | A unified approach to false discovery rate estimation[END_REF], on synthetic p-values coming from possibly unrealistic mixtures, typically Beta-Uniform, Uniform-Uniform, or Uniform-Truncated exponential distributions. The advantage being that in these cases, the pdf of p under H 1 is parametrically known in a simulation context. In an actual situation like in our setting considering the carbon data, the p-values come from two-sample tests that are of different types due to the underlying assumed distribution for the populations that have to be compared (soil carbon levels). Typically, when these can be assumed Gaussian, usual Student t-tests are performed, but when this Gaussian assumption is inappropriate, Mann-Whitney (MW) nonparametric (unpaired) tests are required. In addition, when rejection occurs, the distribution of the critical value under the alternative hypothesis depends on the actual difference (e.g., shift) between the two samples distributions, which is varying and unknown. Our approach relies on the assumption that the complexity of the pdf (p|H 1 true) ∼ f 1 is modelled by the flexibility brought by our nonparametric assumption for it, even though the pdf of (p|H 0 true) remains

U [0,1] .
Finally, note that it is common to consider the probit transform of the p-values since the known component pdf f 0 simply becomes then the density of the standard normal N (0, 1). This is also a good choice in our case for technical reasons since one of our algorithms require kernel density estimates that behave better for observations on the real line instead of [0, 1]. In the present paper we will in practice consider the observed data x = (x 1 , . . . , x n ) where x i = probit(p i ) for i = 1, . . . , n.

EM algorithms

The EM algorithm, as defined in the seminal paper of [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], is more properly understood to be a class of algorithms, a number of which predate even the [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] paper in the literature. These algorithms are designed for maximum likelihood estimation in missing data problems, of which finite mixture problems are canonical examples because the unobserved Z ij 's give an easy interpretation of missing data. Intuitively, EM algorithms replace the unfeasible maximization of the likelihood of the observed data by the maximization of a pseudo likelihood that resembles the likelihood of the complete data, which is itself easy to maximize. For a comprehensive and recent account of EM algorithms, we refer to [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]; here we only briefly describe the finite mixture case.

We consider the complete data ((x 1 , Z 1 ), . . . , (x n , Z n )) = (x, Z) associated with the actually observed sample denoted generically by x (which, for us, represents the vector of probit-transformed p-values from the n statistical tests). In parametric mixture models, the complete data log-likelihood function is easy to maximize over θ. The EM algorithm takes advantage of this by iteratively maximizing instead (since Z is not observed), its expectation conditionnally on x, under the assumption that the parameter driving the random behavior of Z at iteration t is the current value θ (t) . Computing this expectation usually denoted Q(θ|θ (t) ), reduces to the computation of the already introduced posterior probabilities

P (t) ij = P θ (t) (Z ij = 1 | X i = x i ) = λ (t) j f (x i ; ξ (t) j ) m j ′ =1 λ (t) j ′ f (x i ; ξ (t) j ′ ) , (2)
for i = 1, . . . , n and j = 1, . . . m. Each P (t) ij is the posterior probability of individual i coming from component j, given the current θ (t) and the observed data x i . The EM iteration θ (t) → θ (t+1) is defined by

1. E-step: compute Q(θ|θ (t) ) i.e. compute the P (t) ij 's by (2) 2. M-step: set θ (t+1) = argmax θ Q(θ|θ (t) ).
Conveniently, the M-step for finite mixture models always looks partly the same: whatever the assumption on the f j 's, the updates to the mixing proportions are given by λ

(t+1) j = 1 n n i=1 P (t) ij , j = 1, . . . , m. (3)
The updates for the f j or the ξ j parameters depend on the particular form of the component densities and are easy for, e.g., Gaussian mixtures. These will be precised below for our models.

2 Two mixture models for FDR estimation

Gaussian mixture with one component known

This first model that we consider is a 2-component parametric Gaussian mixture with component f 0 forced to N (0, 1), which is a is a slight difference with the Gaussian EM fit from [START_REF] Mclachlan | A simple implementation of a normal mixture approach to differential gene expression in multiclass microarrays[END_REF]. The mixture pdf is then

g θ (x) = λ 0 f 0 (x) + (1 -λ 0 )φ(x; µ 1 , σ 2 1 ), (4)
where φ( • ; µ 1 , σ 2 1 ) denotes here the pdf of the Gaussian N (µ 1 , σ 2 1 ), and the (parametric) model parameter is θ = (λ 0 , µ 1 , σ 2 1 ).

The standard EM algorithm for Gaussian mixtures is straightforward to implement (see, e.g., [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF] and we do not give it here for brevity. The specific version of this algorithm for model (4) turns out to be a very simple case of more general EM-type algorithms for mixtures with constrained parameters, investigated recently by [START_REF] Chauveau | ECM and MM algorithm for mixtures with constrained parameters[END_REF]. Of course constraining here µ 0 = 0 and σ 2 0 = 1 is simple enough that it does not really require these far more general extensions. We only mention it because these new algorithms are already coded in the current version of the mixtools package [START_REF] Benaglia | mixtools: An R package for analyzing finite mixture models[END_REF], and we use it here to run model (4) as a "parametric benchmark" in sections 3 and 4.

Semiparametric mixture with one component known

The purpose of this model is, while keeping f 0 known to be N (0, 1) as before, to relax the parametric assumption on f 1 as explained in Section 1.2. We rename f 1 simply f here, and assume that it is an even (f (x) = f (-x)) but otherwise unspecified density. With such partially nonparametric models, the question of identifiability arises. Basically, this statistical property insures that different values of the parameter must generate different probability distributions of the observable variables. It turns out that identifiability requires here an additional location-shift scalar parameter µ [START_REF] Bordes | Semiparametric estimation of a two-component mixture model where one component is known[END_REF]. The model pdf becomes then

g θ (x) = λ 0 f 0 (x) + (1 -λ 0 )f (x -µ), (5)
with semi-parametric model parameter θ = (λ 0 , µ, f ). This identifiable mixture model may be viewed as a special case of the model for which Benaglia et al. (2009a, section 4.3) have proposed a semi-parametric EM-like algorithm (spEM), that was purely methodological and was called "EM-like" because it does not satisfy a provable convergence property like a genuine EM. We do not get here into details about this methodology, but this original algorithm happened to be very similar to a recent modified version that relies on more solid theoretical grounds, see [START_REF] Chauveau | Estimation for conditional independence multivariate finite mixture models[END_REF]. The algorithm we propose here for our semi-parametric FDR framework is in addition specific for two reasons. First, it requires constraining f 0 to the normal density, and second, to adapt the M-step for the parametric part (λ 0 , µ) and the nonparametric part f for which a weighted kernel density estimate (wkde) is combined with a symmetrization step:

spEM algorithm with one component known Let K be a kernel density function, and h a bandwidth chosen by the user or data-driven (see Section 3). For a current value θ (t) = (λ (t) 0 , µ (t) , f (t) ) of the parameter, the iteration θ (t) → θ (t+1) is as follows: E-step: compute, for i = 1, . . . , n

P (t) i0 = λ (t) 0 λ (t) 0 f 0 (x i ) + (1 -λ (t) 0 )f (x i -µ (t) ) , P (t) i1 = 1 -P (t) i0 M-step 1: parametric part λ (t+1) 0 = 1 n n i=1 P (t) i0 , µ (t+1) = n i=1 P (t) i1 x i n i=1 P (t) i1 (6)
M-step 2: nonparametric density estimation

f (u) = 1 nh n i=1 P (t) i1 K u -x i + µ (t+1) h (7)
and symmetrize, i.e. set f t+1 (u) = ( f (u) + f (-u))/2.

A simulation study

We illustrate the application of the two mixture models presented in Section 2, fitted by these specific EM and spEM algorithms, together with the fdrtool method [START_REF] Strimmer | A unified approach to false discovery rate estimation[END_REF]. We compare the approaches on synthetic data simulated from known models, intended to represent a similar set-up to that presented by the actual Carbon data.

Methodology

We have simulated a multiple testing procedure with n two-samples comparison tests. Let us denote for a single case these two samples by

(Y k 1 , . . . , Y k n k ) i.i.d.∼ F k , for k = 1, 2.
For not interesting cases (H 0 true), F 1 = F 2 , whereas for interesting cases we have simulated F 2 as the distribution F 1 shifted by some δ > 0. For 36% of the test cases, parametric Student t-tests on Gaussian populations (F 1 normal) were simulated. For the remaining 64% cases, F 1 was set to a Student t(5) heavy-tailed distribution, and unpaired nonparametric MW tests were performed. These 36% and 64% proportions correspond approximately to the observed proportions in the carbon data. Each comparison has been based on two small samples (n 1 = 5 and n 2 = 7) where, for interesting cases, the distribution of the second sample was shifted by δ = 3 so that the tests were not too powerful, i.e. rejection of H 0 was not too obvious. We also tried random shifted values δ.

We ran R = 500 replications of n = 2000 such tests, with the proportion of true H 0 set to λ 0 = 0.4 and λ 0 = 0.6, i.e. much larger than in usual microarray experiments, but corresponding roughly to the possible observed proportion in the carbon data. Both EM and semi-parametric EM need an initial parameter value for their first E-step. It is a common practice to provide initial centers from which a binning of the data is performed using a K-means clustering algorithm. Since the probit transform results in normally-distributed x values for not interesting cases, and negative values for the small p-values corresponding to interesting cases, we set these initial centers to (0, -3) for both algorithms, where the value -3 comes from the leftmost peak of the typical histogram of x we obtained.

For comparing the methods on replications, we focused on two criterion motivated by the end user concern: proportion of interesting/not interesting cases and control of the FDR.

Error on the proportion λ 0 : The estimation of the proportion λ 0 of not interesting cases (H 0 true) estimated by each method on the rth replication is denoted λ(r) 0 . We can then evaluate the empirical Mean Square Error for each method by

M SE(λ 0 ) = 1 R R r=1 ( λ(r) 0 -λ 0 ) 2 .
Error on the target FDR level The error on the target level α which is supposed to be achieved by an FDR control procedure may also be evaluated in this simulation context. At replication r, let p (r) = (p

(r) 1 ≤ • • • ≤ p (r)
n ) be the n ordered p-values, and F DR(p) be the FDR estimated by one of the competing methods for p ∈ [0, 1]. In practice, this estimate is used to define how many of the smallest observed p-value have to be rejected, in order to achieve an estimated error level smaller than α. Using the ordered p-values this number of rejections can be defined by the index

d(r) α = dα := max i ∈ {1, . . . , n} : F DR p (r) i ≤ α , (8)
where we will drop the superscript r below to avoid redundancy in some notations. This index corresponds to the largest ordered p-value before which the estimated FDR crosses the level α for the "last time" (note that when this estimated F DR(p) is an increasing function of p, (8) is just the index of the p-value after which the estimated FDR becomes greater than α). Now in our simulation setup, since we observe the complete data (p, Z) (r) , we can compute the true FDR at each replication,

F DR(p (r) i ) = 1 i i ℓ=1 I Z (r) ℓ0 =1 ,
from which we can evaluate the actual error F DR(p (r) dα

) achieved when the dα smallest p-values are rejected at replication r, so that

∆(α) = 1 R R r=1 F DR p (r) dα
α 2 can be viewed as a MSE on the target level α over all the replications.

Results

An illustrative output from a single replication result of our Monte-Carlo experiment is provided in Fig. 1, which displays the behavior of the three FDR estimates we compare, together with the true FDR, BH and rawp-value approaches. The number of rejections for each method based on definition (8) is given in the legend of Fig. 1, together with the value from the true FDR. We can see that for this particular replication, fdrtool and spEM suggest similar numbers of rejections dα = 793 and dα = 798 for the target level α = 10%, both slightly conservative since the true FDR available on these simulated data leads to 817 rejections. This figure also shows graphically how these quantities dα and the error on α depend only on the p-values p, the estimated FDR and the level α. The Benjamini-Hochberg (BH) procedure suggests dα = 680 rejects, which results in an actual error F DR(p 680 ) = 3.82%, that is too conservative. The crude approach based on the uncorrected p-values rejects the dα = 860 smallest p-values, that is non-conservative (as expected) with an actual error F DR(p 860 ) = 11.74%, larger than the desired α = 10%. Global results of the experiment for λ 0 = 0.6 are displayed in Fig. 2 and Table 1. The semi-parametric EM (spEM) with component 1 constrained to the normal N (0, 1) returns better estimates than the Gaussian EM with the same constraint (which is not surprising), but also more precise estimates in terms of bias and variance, for both λ 0 and the target FDR of 10%, than fdrtool. We can see in Table 1 and the right-hand pane of Fig. 2 that spEM is on average less conservative for the target level α than fdrtool, and considerably better than BH or the Gaussian EM. We observed similar results in our second experiment, based on λ 0 = 0.4 (results not presented for conciseness).

We also tried different experiments allowing the distribution F 2 for interesting cases to be shifted by a random variable uniformly distributed. This results in a more complicated component density for interesting cases, rendering it more difficult to recover. When the shift is not too small, e.g., δ ∼ U [3,5] , results are similar to those in Fig. 2. However, this aspect warrants further investigation since for other random shift settings such as δ ∼ U [1,5] we cannot tell which method is preferable between spEM and fdrtool.

FDR estimation for the carbon data

In this section, we detail the numerical results on the carbon data. We applied our two mixture models together with the fdrtool method and the Bonferroni and BH [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] historical procedures that are expected to be too conservative. Numerically, the procedures are identical to the one used in Section 3.1, with respect to the initialization of the EM and spEM algorithms from an empirical binning of the data.

Figure 3 displays the estimated F DR(•)'s together with their associated number of rejected smallest p-values dα as in Fig. 1, except that here we do not know the truth and of course cannot compute the error with respect to the level α. The range of the numbers of rejected hypotheses for a ten percent level α (i.e. less than ten percent of false positives within rejected H 0 among all the n = 2714 tests performed) is somewhat important. Indeed, the Bonferroni procedure suggests 317 rejections among the 2714 tests performed while the fdrtool method suggests 1339 rejections. Moreover, we should reject 1312 null hypothesis on the raw p-values, i.e. without controlling the FDR, 1041 with BH and 1144 with the spEM algorithm.

It is somewhat surprising that fdrtool and spEM seem to behave differently here than in the Monte-Carlo experiment of Section 3, where fdrtool appeared slightly more conservative than spEM. Moreover, the (normally non-conservative) raw p-values approach is here slightly more conservative than fdrtool with 1312 rejects. This comparison with the raw p-values suggests that the estimation given by fdrtool should be taken with caution for these particular data, that may be due to a more complex density of p under the alternative hypothesis.

These results show that not taking into account the FDR estimation can lead to over-estimation of the number of cantons subjected to a significant change (for both increase or decrease). For instance, the difference between the spEM method and more classical one [START_REF] Saby | Changes in soil organic carbon in a mountainous French region, 1990-2004[END_REF][START_REF] Lemercier | Soil phosphorus monitoring at the regional level by means of a soil test database[END_REF] concerns 168 cantons, that is 12.8 % of the area. This over-estimation represents about 16 % for the canton showing an OC increase, and 10% for an OC decrease. Given the need to accurately quantify SOC change for national carbon accounting, this difference is of great importance. Fig 4 maps the results of three of the six procedures tested: raw pvalues, fdrtool and spEM procedures. The Bonferroni and BH [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] historical procedures and gaussEM are not plotted. The spatial structures of the area where the tests suggest significant changes in the carbon content are quite similar on the three maps. This suggests that the spatially isolated cantons where significant changes are detected in the database are more expected to show a false positive test. Globally, losses in organic carbon are observed in Northern France, along the Atlantic coastal regions, and to a lesser extent in the North-Eastern regions. The OC increases are more scattered over the territory.

We also tried to estimate the minimum sample size required at each campaign in order to have a reasonable probability of detecting a given change in the carbon level of a canton between two time periods. For each canton i = 1, . . . , n here, denote by D i the relative difference between the medians of the two samples (between the two time periods), and by S i the smallest size of the corresponding two samples. We tried to fit an Analysis of Covariance model for the predictor D, response variable S, and a categorical factor with two levels, "evolution" or "no evolution" based on the decision rule at a 10% level of the spEM estimate of the FDR. For the numerical predictor, a linear regression model of a log transform, log(S i ) = α log(D i )+ β seemed appropriate. Fig. 5 gives (on the linear scale) the upper prediction intervals for typical levels 50%, . . . , 90%, for the level "evolution" of the categorical predictor. As a matter of example, if we apply the mean yearly decline of 0.6% found by [START_REF] Bellamy | Carbon losses from all soils across England and Wales 1978-2003[END_REF] in England and wales over a 25 yr period , the cumulative relative difference between 1978 and 2003 is about 14% and therefore the 50% or 90% upper prediction interval of the sample sizes needed to detect the changes are 190 and 470, respectively.

Discussion

We applied a large scale Monte-Carlo experiment to compare the recent and standard fdrtool method with mixture models, and in particular the new spEM approach. We have shown that, at least in some cases, spEM performs better in terms of FDR precision and number of rejected null hypotheses. This suggests that in real applications, it would be sensible to at least test both methods to see whether or not they agree. Indeed, for the carbon data, spEM and fdrtool did not agree. The somewhat bizarre results given by fdrtool (less conservative than the raw p-values themselves) can be viewed as an argument to claim that our spEM approach may be more robust in such a case based on a complicated dataset (this complexity may be induced here by p-values coming from two kinds of tests).

Some technical computational details deserves further work. In particular the bandwidth choice in the kernel density estimate involved within the spEM algorithm has been done using Silverman's rule (R default) on the leftmost part of the sample of the probit transform of the p-values, since the values y < c for some c < 0 are mostly associated with the interesting cases (rejections). The tuning of c itself, empirically done here, may be improved.

In this paper, we have detected large areas in France where significant changes of carbon contents occurred during the fifteen years period of this study. Differences in soil carbon content were observed between our two sampling periods. However, the sampling design used does not allow to draw statistical conclusions Nationwide. Indeed, soil carbon changes have already been shown by various studies using numerous paired samples taken exactly at the same location, which is not the case in our study (e.g., [START_REF] Bellamy | Carbon losses from all soils across England and Wales 1978-2003[END_REF]. Therefore, some effective changes might not have been detected using our dataset. Our conclusions might be affected by a bias inherent in the unsupervised sampling strategy associated with the French soil test database: farmers' intentions towards C testing could have changed over time. Moreover, sampling practices may have changed over time (e.g., sampling depth, numbers of subsamples, etc). The spatial and temporal trends described from the national soil test database should be compared with data from two other French soil survey programmes: first, the French Soil Monitoring Network (RMQS) [START_REF] Arrouays | A new initiative in France: a multi-institutional soil quality monitoring network[END_REF], to confirm these trends and assess biases, and second, the inventory programme (IGCS) to take into account the pedological context of soil test values. Indeed, the RMQS grid is more appropriate to a statistical assessment of changes over the national French territory, however it will take at least 10 years before the second sampling campaing is achieved. The IGCS data might be a useful tool since it includes older data (from the 1950's up to now) but it needs to address spatio-temporal correlations effects as the sampling date is related to soil mapping programmes progresses and therefore, the age of data from a given area may be very different to another one.

Figure 1 :Figure 2 :

 12 Figure 1: Example of F DR(p) estimates from five methods, for n = 2000 tests as detailed in Section 3.1. The number of rejections dα defined in (8)for the target level α = 0.1, and actual error achieved by each method are given in the legend. Graphical illustration of computation of dα and error on α are given for the Benjamini-Hochberg (BH) conservative procedure (black dashed) and the non conservative approach using the raw p-values (grey dashed). For each method, dα corresponds to the index (vertical line) located where F DR(•) crosses the level α. Then the level of the horizontal line drawn from each true F DR(p dα ) indicates the actual FDR achieved.

Figure 3 :

 3 Figure3: FDR estimates from six approaches for the carbon data. Bonferroni and BH[START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] historical procedures are also displayed for comparison.

Figure 4 :

 4 Figure 4: Maps of the results of the multiple tests procedures at the cantonal level for a target level α = 0.1. The upper left map corresponds to the raw p-values. The upper right map corresponds to the number of samples used in the procedure. The lower left map uses the FDR estimated by spEM. The lower right uses the FDR estimated by fdrtool. Und corresponds to the cantons where H0 cannot be rejected. ND corresponds to the cantons where there are not enough data to compute the test.

Figure 5 :

 5 Figure5: Upper prediction intervals of levels 50%,. . . ,90%, for the minimim sample sizes needed at each campaign to detect a given change expressed in median relative difference between two time intervals, at a 10% FDR level.

Table 1 :

 1 Mean and MSE's for λ 0 and α from 500 replications of n = 2000 tests, true λ 0 = 0.6.

		spEM	Gaussian EM	fdrtool
	λ0	0.603	0.64	0.609
	M SE(λ 0 ) 0.000149	0.00341	0.000307
	ᾱ	0.0928	0.0424	0.0868
	∆(α)	0.000124	0.00353	0.000350


Radoin Belaouar
email: belaouar@cmap.polytechnique.fr

AND Anne De Bouard
email: debouard@cmap.polytechnique.fr

Arnaud Debussche
email: arnaud.debussche@bretagne.ens-cachan.fr

NUMERICAL ANALYSIS OF THE NONLINEAR SCHR ÖDINGER EQUATION WITH WHITE NOISE DISPERSION

Keywords: Stochastic partial differential equations, numerical analysis, white noise dispersion

This article is devoted to the numerical study of a nonlinear Schrödinger equation in which the coefficient in front of the group velocity dispersion is multiplied by a real valued Gaussian white noise. We first perform the numerical analysis of a semi-discrete Crank-Nicolson scheme in the case when the continuous equation possesses a unique global solution. We prove that the strong order of convergence in probability is equal to one in this case. In a second step, we numerically investigate, in space dimension one, the behavior of the solutions of the equation for different power nonlinearities, corresponding to subcritical, critical or supercritical nonlinearities in the deterministic case. Numerical evidence of a change in the critical power due to the presence of the noise is pointed out.

Introduction

The cubic nonlinear Schrödinger (NLS) equation is a model describing the propagation of light waves in optical fibers. In some circumstances -this is typically the case in dispersion-managed fibers -the group velocity dispersion may vary along the fiber (see [START_REF] Debussche | Nonlinear fiber optics[END_REF][START_REF]Applications of nonlinear fiber optics[END_REF]) and this gives rise in the mathematical model to a NLS equation with a timevarying coefficient in front of the dispersion (Laplace operator). It has been proved in [START_REF] Marty | On a splitting scheme for the nonlinear Schrödinger equation in a random medium[END_REF][START_REF] De Bouard | Debussche The nonlinear Schrödinger equation with white noise dispersion[END_REF] (see also [START_REF] Garnier | Stabilization of dispersion managed solitons in random optical fibers by strong dispersion management[END_REF]) that when this coefficient is random, with mean zero, and under some specific scaling, the limit equation (as the correlation of the noise tends to zero) is a NLS equation with white noise dispersion (i.e. the preceding coefficient becomes a real valued time white noise).

We consider in this paper the more general stochastic nonlinear Schrödinger equation

(1.1) idu + ∆u • dβ + |u| 2σ u dt = 0, u(0) = u 0 ,
where the unknown u is a complex valued random process on a probability space (Ω, F, P) depending on t ≥ 0 and x ∈ R d , and the nonlinear term |u| 2σ u is a power term, σ being a positive real number. The noise term involves a real valued Brownian motion β associated to a filtration (F t ) t≥0 . The product • is a Stratonovich product, and equation (1.1) is formally equivalent to its Itô form:

(1.2) idu + i 2 ∆ 2 u dt + ∆u dβ + |u| 2σ u dt = 0, u(0) = u 0 .

Equation (1.1) has been studied in [START_REF] De Bouard | Debussche The nonlinear Schrödinger equation with white noise dispersion[END_REF] in the case σ < 2/d and in [START_REF] Debussche | 1-D quintic nonlinear Schrödinger equation with white noise dispersion[END_REF] in the case d = 1 and σ = 2, and the existence of a unique global square integrale solution was proved. However, nothing is known about the qualitative behavior of the solutions of (1.1)note that the Hamiltonian structure of the deterministic NLS equation is destroyed by the presence of a time varying dispersion coefficient. While the case of a deterministic time-periodic coefficient has been investigated from the numerical as well as theoretical point of view (see e.g. [START_REF] Kunze | The singular perturbation limit of a variational problem from nonlinear fiber optics[END_REF][START_REF] Turitsyn | Path-averaged chirped optical soliton in dispersion-managed fiber communication lines[END_REF][START_REF] Zharnitsky | Stabilizing effects of dispersion management[END_REF]), few investigations have been performed in the stochastic case (see however [START_REF] Garnier | Stabilization of dispersion managed solitons in random optical fibers by strong dispersion management[END_REF] and [START_REF] Marty | On a splitting scheme for the nonlinear Schrödinger equation in a random medium[END_REF] for numerical simulations in the cubic case), and no theoretical result on the qualitative behavior is known.

These considerations motivated us for a more systematic numerical study of equation (1.1) for general values of the parameter σ. We will be in particular interested in the blow up phenomena. In [START_REF] Debussche | DiMenza Numerical simulation of focusing stochastic nonlinear Schrödinger equations[END_REF] and [START_REF] Barton-Smith | DiMenza Numerical study of two dimensional stochastic NLS equations[END_REF], one of the authors conducted such investigations for a different stochastic equation, in which the noise, which was there a space-time noise, white in time and possibly correlated in space, acted as a linear potential in the equation. It was then observed that in the supercritical case σ > 2/d, as the space correlation of the noise tends to zero, the blow up mechanism could be stopped. Note that it was proved in [START_REF] De Bouard | Debussche Blow-up for the supercritical stochastic nonlinear Schrödinger equation with multiplicative noise[END_REF] that if on the opposite, the noise is spatially sufficiently correlated, then blow up occurs for any initial data with a positive probability, still in the case σ > d/2. We will here do the same investigation in the, somewhat simpler, case of equation (1.1), where the noise only depends on time (see Section 3).

In a first step, we study the order of convergence for equation (1.1), of a semi-discrete Crank-Nicolson scheme, which is one of the schemes used in the numerical simulations of Section 3. As is natural for a Stratonovich noise, we use a mid-point discretization of the noise. The scheme then presents the advantage of preserving the L 2 -norm, which is a conserved quantity of the continuous equation. The same kind of schemes has been considered in [START_REF] De Bouard | A semi-discrete scheme for the stochastic nonlinear Schrödinger equation[END_REF] where the convergence in probability to the solution of the continuous equation is proved, in the case where the noise is a space-time potential and is white in time. The order of convergence of a Crank-Nicolson semi-discrete scheme as also been studied in [START_REF] Gazeau | Strong order of convergence of a semidiscrete scheme for the stochastic Manakov equation[END_REF], for the Manakov-PMD system, a system of NLS equations with a threedimensional time white noise appearing in the equation as a coefficient of the gradient of the solution.This model also appears in nonlinear fiber optics, and the (strong) order of convergence in time was then found to be at least 1/2. The order of convergence of a splitting scheme was studied in [START_REF] Marty | On a splitting scheme for the nonlinear Schrödinger equation in a random medium[END_REF] for equation (1.1), except that the nonlinear term was replaced by a Lipschitz function. Asymptotic preserving properties of the scheme have also been proved ( [START_REF]Duboscq Analyse et simulations numriques d'quations de Schrdinger dterministes et stochastiques[END_REF]). The order of convergence was also found to be at least 1/2, although it is expected to be equal to one in this case, since the noise is one-dimensional. We indeed prove in Section 2 that the strong order in probability for the time semi-discrete Crank-Nicolson scheme is at least 1.

Note that the spatial discretization is very easy to analyse. Indeed, from our argument it is clear that it suffices to study (1.1) with a Lipschitz nonlinearity. Then, classical arguments show that the spatial discretization error is the same as in the deterministic case.

We now give some notations. We will denote by F x the spatial Fourier transform in R d and û denotes also the (space) Fourier transform of the (complex-valued) function u. For any integer (or even real) number s, we denote by H s (R d ) the Sobolev space of complex-valued distributions u on R d satisfying (I -∆) s/2 u ∈ L 2 (R d ) and by • H s its norm. This is a real Hilbert space, which thanks to the Plancherel Theorem is endowed with the inner product

(u, v) s ≡ Re R d (1 + |ξ| 2 ) s û(ξ) v(ξ) dξ.
Note that here we identify C with R 2 , and we will say that a function f defined on

H s (R d ) is in C 2 b (H s ) if f is of class C 2 on H s (R d ; R) × H s (R d ; R
) and has bounded derivatives, i.e. there is a constant C > 0 such that for any u, v, w

∈ H s (R d ) Df (u).v H s ≤ C v H s and D 2 f (u)(v, w) H s ≤ C v H s w H s ,
where D 2 f (u) denotes the second order derivative of f at u, and is a bilinear operator from

H s (R d ) × H s (R d ) into H s (R d ).
With obvious modifications in the definition, we will also consider functions g which belong to C 2 b (H s × H s ). We will in that case denote D 1 g (resp. D 2 g) the partial differential with respect to the first (resp. second) argument of g.

2.

Numerical analysis of the Crank-Nicolson scheme 2.1. Preliminaries and main result. We consider in this section equation (1.1) in the cases where σ is an integer, which makes the power term |u| 2σ u regular in u and allows us to differentiate it as many times as we need. Moreover, we will assume that the existence and uniqueness of a global solution is known, which is the case e.g. if σ = 1 or σ = 2, and the space dimension d is equal to one (see Proposition 2.1 below). We recall (see [START_REF] De Bouard | Debussche The nonlinear Schrödinger equation with white noise dispersion[END_REF][START_REF] Marty | On a splitting scheme for the nonlinear Schrödinger equation in a random medium[END_REF]) that equation (1.1) may also be written in the mild form

(2.1) u(t) = S(t, 0)u 0 + i t 0 S(t, s)(|u| 2σ u(s))ds,
where the random propagator S(t, s) assigns to any F s -measurable L 2 -valued random variable u s the solution at time t > s of the linear equation

(2.2) idu + i 2 ∆ 2 u dt + ∆u dβ = 0, u(s) = u s .
Note that S(t, s) has an explicit expression in Fourier variables given by

F(S(t, s)u s )(ξ) = e -i|ξ| 2 (β(t)-β(s)) ûs (ξ), t ≥ s, ξ ∈ R d ,
(see [START_REF] De Bouard | Debussche The nonlinear Schrödinger equation with white noise dispersion[END_REF]) and this allows to generalize the definition of S(t, s)u to non necessarily F smeasurable random functions u.

Our assumption on global existence and uniqueness result is known to be true in some cases: it has been proved for s = 1, in [START_REF] De Bouard | Debussche The nonlinear Schrödinger equation with white noise dispersion[END_REF] in the case σ = 1 and in [START_REF] Debussche | 1-D quintic nonlinear Schrödinger equation with white noise dispersion[END_REF] for the case σ = 2. The extension of the result to a general Sobolev space H s (R d ) with s ≥ 1, i.e. the fact that the solution u has trajectories a.s. in

C(R + ; H s (R)) if u 0 ∈ H s (R) is an immediate consequence of the embedding H 1 ⊂ L ∞ (R)
and of the isometry of S(t, s) in H s (R) for any s. We recall these results for completeness. We expect that global existence and uniqueness can be proved for σ < 4/d -see the numerical simulations of Section 3 -so that our result would also hold true with σ = 1 for d = 2, 3.

Our aim in this section is to study the following time discretization of equation (1.1), where δt > 0 denotes the time step :

(2.3) i u n+1 -u n δt + χ n √ δt ∆u n+1/2 + g(u n , u n+1 ) = 0,
where the function g is an approximation of the nonlinear term |u| 2σ u given by (2.4)

g(u, v) = 1 0 θ|u| 2 + (1 -θ)|v| 2 σ dθ u + v 2 = 1 σ + 1 |u| 2σ+2 -|v| 2σ+2 |u| 2 -|v| 2 u + v 2 .
Here u n+1/2 = 1 2 (u n + u n+1 ). Moreover, the random variables χ n are defined by

χ n = β((n+1)δt)-β(nδt) √ δt
. Hence the family (χ n ) n ∈ N is an independent family of normal random variables. With these definitions, u n is then an approximation of u(t n ), where t n = nδt and u is the solution of equation (1.1).

The mid-point discretization of the noise term is natural to approximate a Stratonovich product. This scheme is a generalization to equation (1.1) of the classical Crank-Nicolson scheme for the deterministic NLS equation, which is known to be conservative (in the deterministic case) for the L 2 -norm and the Hamiltonian (see e.g. [START_REF] Delfour | Payre Finite difference solutions of a nonlinear Schrödinger equation[END_REF] or [START_REF] Sanz-Serna | Methods for the numerical solution of the nonlinear Schrödinger equation[END_REF]). For the stochastic NLS equation (1.1) with a power nonlinearity, the scheme is conservative for the L 2 -norm (see Section 3).

Another possibility would have been to consider an explicit discretization of the noise term in the Itô formulation (1.2) of the equation. However, such a discretization would not lead to a conservative scheme. On the other hand, the scheme (2.3) may be written as

i u n+1 -u n δt + χ n √ δt ∆u n + 1 2 χ n √ δt ∆(u n+1 -u n ) + g(u n , u n+1 ) = 0
from which we deduce (assuming that each term above is bounded and that g is sufficiently regular and bounded)

∆(u n+1 -u n ) √ δt = iχ n ∆ 2 u n + O( √ δt);
plugging this expression into the preceding equation, we obtain up to the O( √ δt) term the new formulation for the scheme (2.3)

i u n+1 -u n δt + χ n √ δt ∆u n + i 2 χ 2 n ∆ 2 u n + g(u n , u n+1 ) = 0,
which is exactly the Milstein scheme for the Itô equation (1.2) (see [START_REF] Milstein | Approximate integration of stochastic differential equations[END_REF]). Hence we may expect that the scheme has an order of convergence in time equal to one. This is indeed the main result of this section (see Theorem 2.3 below). We now introduce some notations. For a fixed n ∈ N, we denote by S n,δt the linear operator

S n,δt = (I -i √ δt 2 χ n ∆) -1 (I + i √ δt 2 χ n ∆).
This operator is a Fourier multiplier, it is a bounded linear operator on L 2 (R d ) satisfying

F(S n,δt u)(ξ) = 1 -i √ δt 2 χ n |ξ| 2 1 + i √ δt 2 χ n |ξ| 2 û(ξ)
for any u ∈ L 2 (R d ). We also denote by T n,δt the linear operator T n,δt = (I -i

√ δt 2 χ n ∆) -1
, which is also a linear bounded operator in L 2 (R d ). Then the scheme (2.3) may be written as (2.5) u n+1 = S n,δt u n + i(δt)T n,δt g(u n , u n+1 ).

At that point, we should mention that because g is only a locally Lipschitz function on

C 2 -hence the Nemytskii operator (u, v) → g(u, v) is locally Lipschitz in H s (R d ) × H s (R d ) for any s > d/2 -the existence, for a given u n in H s (R d ), of u n+1 ∈ H s (R d ) satisfying (2.5
) is not obvious ; moreover, u n+1 will not be unique in general. In order to define the solution in a simple way, we thus introduce a cut-off in the nonlinear part of the scheme.

Let θ ∈ C ∞ (R + ) with θ ≥ 0, supp θ ⊂ [0, 2] and θ ≡ 1 on [0, 1] ; we then set θ k (x) = θ( x k ), for x ≥ 0 and k ∈ N. Now, in all what follows, the integer s > d/2 is fixed, and we set for k ∈ N and u, v ∈ H s (R d ) :

(2.6)

g k (u, v) = 1 2σ + 2 θ k (|u| 2 H s )θ k (|v| 2 H s ) |u| 2σ+2 -|v| 2σ+2 |u| 2 -|v| 2 u + v 2 .
It is clear that g k is a Lipschitz function on H s (R d ) and satisfies : for any

u 1 , u 2 , v 1 , v 2 ∈ H s (R d ), |g k (u 1 , v 1 ) -g k (u 2 , v 2 )| H s ≤ C(σ, k)(|u 1 -u 2 | H s + |v 1 -v 2 | H s ).
We then replace the scheme (2.5) by (2.7)

u k n+1 = S n,δt u k n + i(δt)T n,δt g k (u k n , u k n+1
). The next proposition may be proved as Proposition 3.1 in [START_REF] Gazeau | Strong order of convergence of a semidiscrete scheme for the stochastic Manakov equation[END_REF], using the integral formulation of (2.7) and a fixed point argument in the subspace of L ∞ (0, T ; H s (R d )) consisting of functions which are constant on each interval [nδt, (n + 1)δt], 0 ≤ n ≤ N -1, with δt = T /N , provided that this latter is less than 1/C(σ, k), where C(σ, k) the Lipschitz constant of g k , together with the fact that the operators S n,δt and T n,δt are bounded in H s (R d ) with bound less than one. For general δt, one obtains the existence of a discrete solution by using a measurable selection theorem (see e.g. [START_REF] De Bouard | Debussche Weak and Strong Order of Convergence of a Semidiscrete Scheme for the Stochastic Nonlinear Schrödinger Equation[END_REF][START_REF] Gazeau | Strong order of convergence of a semidiscrete scheme for the stochastic Manakov equation[END_REF]), but the solution is not unique in general.

Proposition 2.2. Let T > 0 and δt = T /N > 0 be fixed. Let u 0 ∈ H s (R d ) with s > d/2; then for any k ∈ N * , there is an adapted discrete solution u k δt = (u k n ) n=0,...,N satisfying (2.7) for n ≤ N -1, which belongs a.s. to L ∞ (0, T ; H s (R d )), and the solution u k δt is unique if δt is sufficiently small, depending on k. Moreover, for any n ≤ N , the L 2 norm is preserved that is u k n L 2 = u 0 L 2 . Now, we define the stopping time (2.8) τ k δt = inf{nδt, u k n-1 H s ≥ k, or u k n H s ≥ k}, which is nondecreasing with k. It is clear that (u k n ) n=0,.
..,n 0 -1 satisfies the scheme (2.5) as long as n 0 δt ≤ τ k δt . Moreover the discrete solution given by Proposition 2.2 can be constructed in such a way that u k+1 n = u k n for nδt ≤ τ k δt .. Finally, we set for a fixed δt = T /N > 0:

u δt n = u k n if nδt ≤ τ k δt .
This allows to define the discrete solution of the scheme (2.3) up to the stopping time τ

∞ δt ≡ lim k→∞ τ k δt ∧ T . For t ≥ τ ∞ δt , we just set u δt n = u ∞ , where u ∞ is any arbitrary element in H s (R d ).
We are now in position to state our main result, which says that the scheme (2.3) has strong order one in probability.

Theorem 2.3. Let T > 0; assume u 0 ∈ H s (R d ), s > d/2, is such that equation (2.1) has a unique adapted solution u with paths a.s. in C([0, T ] : H s (R d )). For δt = T /N , with N ∈ N * , let t n = nδt, n = 0, ..., N ; then for any n ≤ N , lim C→∞ P u(t n ) -u δt n H s ≥ Cδt = 0, uniformly with respect to δt.
Remark 2.4. We have assumed for simplicity that the solution of the continuous equation (2.1) was defined on a deterministic time interval, but the result would also be true if the solution was only locally defined in time, on a random time interval which is a stopping time (see [START_REF] De Bouard | Debussche Weak and Strong Order of Convergence of a Semidiscrete Scheme for the Stochastic Nonlinear Schrödinger Equation[END_REF][START_REF] Gazeau | Strong order of convergence of a semidiscrete scheme for the stochastic Manakov equation[END_REF]). In particular the result would be true for

s > d/2 in H s (R d ),
or in H s (T d ), where T = R/2πZ, using in all what follows Fourier series instead of the continuous Fourier transform Remark 2.5. We actually expect that a stronger result holds under the assumptions of Theorem 2.3, namely we expect that

lim C→∞ P sup nδt≤T u(t n ) -u δt n H s ≥ Cδt = 0,
uniformly with respect to δt. This would e.g. follow from the arguments in [START_REF] De Bouard | Debussche Weak and Strong Order of Convergence of a Semidiscrete Scheme for the Stochastic Nonlinear Schrödinger Equation[END_REF] if we were able to prove a mean square order one, for the truncated scheme 2.7, in the

L ∞ (0, T ; H s (R d )) norm (see Proposition 2.6) below.
In turn, this latter result would follow from the arguments in [START_REF] Faure | Simulation du mouvement brownien et des diffusions[END_REF], Chapter 5, if we were able to prove that the scheme 2.7 approximates the continuous equation at order one in L p (Ω) for any p ≥ 1. However, the scheme being implicit in the noise, we cannot use the Itô Formula, and obtaining an order in L p (Ω) seems to require very tedious computations, while our argument to get Theorem 2.3 from the mean square result for the truncated scheme given by Proposition 2.6 below is extremely simple.

In order to prove Theorem 2.3 we will thus first prove that the scheme defined by (2.7) gives an approximation of the corresponding continuous equation of order one in mean square. This will be the object of the next subsection.

2.2.

Mean square order for the truncated scheme. In this subsection, we consider the scheme (2.7), where the integer k ∈ N is fixed, and we prove that its solution u k δt given by Proposition 2.2 converges in mean square at order one to the solution of the continuous equation :

(2.9) u(t) = S(t, 0)u 0 + i t 0 S(t, s)f k (u(s)) ds, where f k (u) = θ k ( u 2 H s )|u| 2σ u. Note that, for a fixed k ∈ N, since s > d/2, f k is a Lipschitz function on H s (R d ), and the existence of a unique global adapted solution of (2.9) in L ∞ (Ω; C([0, T ]; H s (R d ))) is immediate, if u 0 ∈ H s (R d ).
Moreover, using an induction argument and the fact that

S(t, s) is an isometry in H α (R d ) for any α, it is not difficult to see that if in addition u 0 ∈ H s+6 (R d ), then u k has paths a.s. in H s+6 (R d ), and there is a deterministic constant C = C(T, k, u 0 H s ) such that (2.10) u k C([0,T ];H s+6 ) ≤ C, a.s.
Proposition 2.6. Let T > 0 and u 0 ∈ H s+6 (R d ). Let k ∈ N be fixed, u k be the unique adapted solution of (2.9) with paths a.s. in C([0, T ]; H s+6 (R d )), and u k δt = (u k n ) n=0,..,N with δt = T /N be an adapted discrete solution of (2.7) given by Proposition 2.2, where g k is defined by (2.6). Then there is a constant C which may depend on T and sup t∈[0,T ] E( u(t) 4 H s+6 ) such that for any δt = T /N ,

sup nδt≤T E( u k n -u k (t n ) 2 H s ) ≤ C(δt) 2
where t n = nδt.

Proof of Proposition 2.6. In all the proof, since the index k refereeing to the cut-off is fixed, we ignore it and simply denote u k n by u n , u k by u, and so on. Moreover, all the constants appearing in the subsection may depend on k, even though this is not always explicitly stated. We then denote by e n the error e n = u n -u(t n ) where we recall that t n = nδt. In order to prove Proposition 2.6, we first prove that there are positive constants ε, K 1 and K 2 , that possibly depend on sup t∈[0,T ] E( u(t) 4 H s+6 ) such that (2.11)

E( e n+1 2 H s ) ≤ (1 + (δt)K 1 )E( e n 2 H s ) + (δt) 3 K 2 for all δt ≤ ε and all n ∈ N such that nδt ≤ T .
Using the mild form of the scheme (2.7), and the expression (2.12)

u(t n+1 ) = S(t n+1 , t n )u(t n ) + i t n+1 tn S(t n+1 , s)f k (u(s))ds.
which is deduced from (2.9), one may express e n+1 in terms of e n as follows :

(

= i(δt)T n,δt (g k (u n , u n+1 ) -g k (u(t n ), u(t n+1 ))) .
We will estimate separately the contribution of each of the above terms to the left hand side of (2.11). Indeed, we may write using (2.13) :

(2. [START_REF] Milstein | Approximate integration of stochastic differential equations[END_REF])

E( e n+1 2 H s ) = 5 j=1 E( e n,δt,j 2
H s ) + 2 1≤i<j≤5
E(e n,δt,i , e n,δt,j ) s where we recall that (., .) s denotes the inner product in H s (R d ).

First, it is obvious from its expression in Fourier space (we recall that χ n is a real valued Gaussian random variable) that the operator S n,δt is a unitary operator on H s (R d ). It follows (2.20) E( e n,δt,1

2 H s ) = E( e n 2
H s ). The second term, which corresponds to the approximation of the linear equation, is estimated thanks to the next lemma.

Lemma 2.7. There is a constant C such that for any random function v ∈ L 4 (Ω;

H s+6 (R d )), E( S(t n+1 , t n )v -S n,δt v 2 H s ) ≤ C(δt) 3 (E( v 4 H s+6 )) 1/2 . If moreover, v is F tn -measurable, then E( S(t n+1 , t n )v -S n,δt v 2 H s ) ≤ C(δt) 3 E( v 2 H s+6 ).
Applying Lemma 2.7 with v = u(t n ) immediately leads to the estimate

(2.21) E( e n,δt,2 2
H s ) ≤ C(δt) 3 sup t∈[0,T ] E( u(t) 2 H s+6 ),
for any n such that nδt ≤ T .

Proof of Lemma 2.7. Thanks to the expression of the operator S n,δt as a Fourier multiplier, we may write, for any ξ ∈ R d :

F(S n,δt v)(ξ) = 1 -i √ δt 2 χ n |ξ| 2 1 + i √ δt 2 χ n |ξ| 2 v(ξ) = e -2i arctan( √ δt 2 |ξ| 2 χn) v(ξ).
Hence the Parseval Theorem implies that for any random function v ∈ L 2 (Ω;

H s+6 (R d )), E( S(t n+1 , t n )v -S n,δt v 2 H s ) = R d E(|e -i √ δt|ξ| 2 χn -e -2i arctan( √ δt 2 |ξ| 2 χn) | 2 |v(ξ)| 2 )(1 + |ξ| 2 ) s dξ
Now, since |1 -e ix | 2 ≤ C|x| 2 for any real x, and | arctan y -y| ≤ C|y| 3 for any y ∈ R, where the constants C are independent of x and y, the previous equality implies

E( S(t n+1 , t n )v -S n,δt v 2 H s ) ≤ C(δt) 3 E χ 6 n R d (1 + |ξ| 2 ) s |ξ| 12 |v(ξ)| 2 dξ .
If moreover, v is independent of χ n , then the right hand side above is bounded by

C(δt) 3 E(χ 6 n )E( v 2 H s+6
) and we may conclude in the general case using Cauchy-Schwarz inequality.

We now turn to the estimate of E( e n,δt,3 2 L 2 ). For this purpose, we will use a Taylor expansion of the function (in fact the Nemytskii operator) g k . Indeed, we recall that since σ is an integer, there exists a real valued polynomial function p σ (x 1 , x 2 ) such that, for all u, v ∈ H s (R d ) and all x ∈ R d , (2.22)

g k (u, v)(x) = θ k (|u| 2 H s )θ k (|v| 2 H s ) p σ (|u(x)| 2 , |v(x)| 2 ).
It easily follows that g k is differentiable on H s (R d )×H s (R d ), and for any (u, v) ∈ H s ×H s , and any w ∈ H s , the function

D u g k (u, v).w belongs to H s (R d ) and ∀x ∈ R d , (D u g k (u, v).w)(x) = θ k (|u| 2 H s )θ k (|v| 2 H s ) ∂ x 1 p σ (|u(x)| 2 , |v(x)| 2 ) 2Re(ū(x)w(x)) + 2θ ′ k (|u| 2 H s )θ k (|v| 2 H s ) (u, w) s p σ (|u(x)| 2 , |v(x)| 2
). A similar expression holds for D v g k (u, v).w, and we deduce that there is a constant C(k), which only depends on k, such that for any (u, v) ∈ H s × H s , and any (w

1 , w 2 ) ∈ H s × H s , Dg k (u, v).(w 1 , w 2 ) H s ≤ C(k)( w 1 H s + w 2 H s ), in other words, g k ∈ C 1 b (H s × H s ).
In the same way, it may be checked that

g k ∈ C 2 b (H s × H s ), i.e. for any (u, v) ∈ H s × H s , and any (w 1 , w 2 ) ∈ H s × H s , D 2 g k (u, v).(w 1 , w 2 ) 2 H s ≤ C(k)( w 1 2
H s + w 2 2
H s ). Then, by Taylor expansion, we may write, for any s ∈ (t n , t n+1 ), (2.23)

f k (u(s)) = g k (u(s), u(s)) = g k (u(t n ), u(t n+1 )) + 1 0 Dg k (u(t n ) + θE n (s), u(t n+1 ) + θE n+1 (s))(E n (s), E n+1 (s))dθ where E n (s) = u(s) -u(t n ).
H s ) ≤ C(δt) t n+1 tn 1 0 E Dg k (u(t n ) + θE n (s), u(t n+1 ) + θE n+1 (s))(E n (s), E n+1 (s)) 2 H s dθds ≤ C(δt) t n+1 tn E( E n (s) 2 H s + E n+1 (s) 2 H s )ds. 2
In order to estimate the right hand side above, E n (s) is written using again the mild form (2.9) of the continuous equation as :

(2.25) E n (s) = u(s) -u(t n ) = S(s, t n )u(t n ) -u(t n ) + i s tn S(s, σ)f k (u(σ))dσ
and we use the following lemma.

Lemma 2.8. There is a constant C such that for any random function v ∈ L 4 (Ω; H s+2 (R d )), and for any s ∈ (t n , t n+1 ),

E( S(s, t n )v -v 2 H s ) ≤ C(δt)(E( v 4 H s+2 )) 1/2 . If moreover, v is F tn -measurable, then E( S(s, t n )v -v 2 H s ) ≤ C(δt)E( v 2 H s+2 ), and
E( S(s, t n )v -v 4 H s ) ≤ C(δt) 2 E( v 4 H s+2
). Proof of Lemma 2.8. The first inequality requires the same arguments as the proof of Lemma 2.7 : the inequality |e ix -1| 2 ≤ C|x| 2 , which is true for any real x, implies, for any v ∈ L 4 (Ω; H s+2 (R d )), and for any s ∈ (t n , t n+1 ) :

(2.26)

E(|F x (S(s, t n )v -v)(ξ)| 2 ) = E(|e -i(β(s)-β(tn))|ξ| 2 -1| 2 |v(ξ)| 2 ) ≤ CE(|β(s) -β(t n )| 2 |ξ| 4 |v(ξ)| 2 ).
Parseval Theorem and Cauchy-Schwarz inequality then easily lead to the first inequalities of Lemma 2.8. On the other hand, the above argument and Parseval Theorem imply, for any

F tn -measurable v ∈ L 4 (Ω; H s+2 (R d )), E( S(s, t n )v -v 4 H s ) = E R d R d |e i(β(s)-β(tn))|ξ| 2 -1| 2 |e i(β(s)-β(tn))|η| 2 -1| 2 (1 + |ξ| 2 ) s |v(ξ)| 2 (1 + |η| 2 ) s |v(η)| 2 dξdη ≤ C R d R d (1 + |ξ| 2 ) s (1 + |η| 2 ) s |ξ| 4 |η| 4 E(|β(s) -β(t n )| 4 |v(ξ)| 2 |v(η)| 2 )dξdη
which is bounded above by

C(δt) 2 E R d (1 + |ξ| 2 ) s+2 |v(ξ)| 2 dξ 2 = C(δt) 2 E( v 4 H s+2 ),
using again the independence of v and β(s) -β(t n ). The third inequality of Lemma 2.8 follows.

Lemma 2.8, (2.25) and the fact that

f k is Lipschitz in H s (R d ) with f k (0) = 0 imme- diately imply the bound (2.27) E( E n (s) 2 H s ) ≤ C(δt) sup t∈[0,T ] E( u(t) 2 H s+2 ),
for any s ∈ (t n , t n+1 ). Modifying the above arguments in an obvious way, one may also prove that

(2.28)

E( E n+1 (s) 2 H s ) ≤ C(δt) sup t∈[0,T ] E( u(t)
H s ) ≤ 2E( t n+1 tn (I -S(t n+1 , s))g k (u(t n ), u(t n+1 ))ds 2 H s ) + 2(δt) 2 E( (I -T n,δt )g k (u(t n ), u(t n+1 )) 2 H s ds = I + II. 2
The first term on the right hand side above is bounded thanks to the previous arguments : using Cauchy-Schwarz inequality,

I ≤ C(δt)E t n+1 tn (I -S(t n+1 , s))g k (u(t n ), u(t n+1 ))
I ≤ C(δt) 3 sup t∈[0,T ] (E( u(t) 4 H s+2 ) 1/2 .
In order to estimate II, we again use Parseval identity, and noticing that for any ξ ∈ R d ,

1 1 + i √ δt 2 χ n |ξ| 2 -1 2 ≤ δt 4 χ 2 n |ξ| 4 ,
we obtain, again with the use of Cauchy-Schwarz inequality and the fact that

g k is Lipschitz in H s+2 (R d ) × H s+2 (R d ), II ≤ 2(δt) 2 E R d 1 1 + i √ δt 2 χ n |ξ| 2 -1 2 |F x (g(u(t n ), u(t n+1 ))(ξ)| 2 dξ ≤ C(δt) 3 (E(χ 4 n ) 1/2 ) sup t∈[0,T ] [E( u(t) 4 H s+2 )] 1/2 .
From (2.30) and these two estimates, we deduce

(2.32) E( e n,δt,4 2
H s ) ≤ C(δt) 3 1 + sup t∈[0,T ] E( u(t) 4 H s+2 ) .
The estimate of the last term, given by (2.18), simply follows from the boundedness of the operator T n,δt in H s (R d ) and the fact that the function

g k is Lipschitz in H s (R d )×H s (R d ); (2.33) E( e n,δt,5 2
H s ) ≤ C(δt) 2 E( e n 2 H s + e n+1 2
H s ). We now wish to estimate the second term in (2.19), i.e. to consider the H s -inner products of the terms e n,δt,j , for 1 ≤ j ≤ 5. We first note that all the terms of the form E(e n,δt,j , e n,δt,k ) s with 2 ≤ j, k ≤ 4 may be bounded thanks to the Cauchy-Schwarz inequality, since from (2.21), (2.29)

))u(t n ) s = Re E R d (1 + |ξ| 2 ) s ên (ξ)(1 -e -2i arctan( √ δt 2 χn|ξ| 2 )+i √ δtχn|ξ| 2 ) ū(t n , ξ)dξ
and since χ n is independent of û(t n , ξ), the previous term is equal to

Re R d E(1 -e i √ δtχn|ξ| 2 -2i arctan( √ δt 2 χn|ξ| 2 ) )E((1 + |ξ| 2 ) s ên (ξ) ū(t n , ξ))dξ. Since arctan x -x = - x 3 3 + x 5 ε(x)
for some bounded function ε(x) on R, it follows

1 -e i √ δtχn|ξ| 2 -2i arctan( √ δt 2 χn|ξ| 2 ) = 1 -e i 12 (δt) 3/2 χ 3 n |ξ| 6 +i(δt) 5/2 χ 5 n |ξ| 10 ε( √ δt 2 χn|ξ| 2 ) = -i 12 (δt)
) s | ≤ C(δt) 5/2 E( e n 2 H s ) 1/2 E( u(t n ) 2 H s+5 ) 1/2 ≤ C 1 (δt)E( e n 2 H s ) + C 2 (δt) 4 sup t∈[0,T ] E( u(t)
S(t n+1 , s)(g k (u(t n ), u(t n+1 )) -f k (u(s)))ds s = I + II + III with (2.37) I = E e n , i t n+1 tn [f k (u(t n )) -f k (u(s))]ds s , (2.38) II = (δt)E e n , i[g k (u(t n ), u(t n+1 )) -f k (u(t n ))] s and (2.39) III = E e n , i t n+1 tn (I -S * n,δt S(t n+1 , s))[g k (u(t n ), u(t n+1 )) -f k (u(s))]ds s .
The use of a Taylor expansion for f k -which is easily seen to belong to C 2 b (H s (R d ))allows to write (2.40)

I = E e n , i t n+1 tn Df k (u(t n ))E n (s)ds s +E e n , i t n+1 tn 1 0 (1 -θ)D 2 f k (u(t n ) + θE n (s))(E n (s), E n (s))dθds s = I 1 + I 2 .
Thanks to the Cauchy-Schwarz inequality and the fact that

(2.41) E( E n (s) 4 H s ) ≤ C(δt) 2 sup t∈[0,T ] E( u(t) 4 H s+2 ),
which follows easily from Lemma 2.8 and (2.25), the following estimate holds :

(2.42)

I 2 ≤ C(δt) 2 E( e n 2 H s ) 1/2 sup t∈[0,T ] E( u(t) 4 H s+2 ) 1/2 .
In order to estimate the term I 1 , we use again, the expression (2.25) of E n which allows to write

E e n , t n+1 tn Df k (u(t n ))E n (s)ds s = E e n , t n+1 tn Df k (u(t n ))(S(s, t n ) -I)u(t n )ds s +E e n , t n+1 tn Df k (u(t n )).i s tn S(s, σ)f k (u(σ))dσds s .
The second term in the right hand side above is bounded, again thanks to the Cauchy-Schwarz inequality, by

C(δt) 2 E( e n 2 H s ) 1/2 sup t∈[0,T ] E( u(t) 2 H s ) 1/2 .
The first term is equal to

E e n , t n+1 tn Df k (u(t n ))E(S(s, t n ) -I)u(t n )ds s
and the use of Cauchy-Schwarz inequality, together with the fact that Df k (u) is a bounded operator in H s (R d ) uniformly in u, allows to bound this term by

C(δt) 1/2 E e n 2 H s ) 1/2 E( t n+1 tn E(S(s, t n ) -I)u(t n ) 2 H s ds 1/2 .
Then, Parseval Theorem and the fact that, or all s ∈ (t n , t n+1 ),

(2.43)

|E(F x ((S(s, t n ) -I)u(t n ))(ξ) F tn )| = |E(e i(β(s)-β(tn))|ξ| 2 -1)û(t n , ξ)| = |e -(s-tn)|ξ| 4 /2 -1||û(t n , ξ)| ≤ C(δt)|ξ| 4 |û(t n , ξ)|
lead to the estimate :

(2.44)

I 1 ≤ C(δt) 2 E( e n 2 H s ) 1/2 sup t∈[0,T ] E( u(t) 2 H s+4 ) 1/2 .
The term II in (2.36) may be estimated in the same way. In order to bound the third term, we use the following lemma.

Lemma 2.9. There exists a constant C > 0 such that for any random function v ∈ L 4 (Ω; H s+2 (R d )), and for any s ∈ (t n , t n+1 ),

E( (I -S * n,δt S(t n+1 , s))v 2 H s ) ≤ C(δt)(E( v 4 H s+2 )) 1/2 . Applying Lemma 2.9 with v = g k (u(t n ), u(t n+1 )) -f k (u(s))
E( u(t) 4 H s+4 )].
Proof of Lemma 2.9. The proof is based on similar arguments as those of Lemma 2.7 : Parseval theorem and the fact that |1 -e ix | 2 ≤ C|x| 2 for any real x allows to write, for any v ∈ L 4 (Ω;

H s+2 (R d )), E( (I -S * n,δt S(t n+1 , s))v 2 H s ) ≤ CE R d 2 arctan( √ δt 2 χ n |ξ| 2 ) + (β(t n+1 ) -β(s))|ξ| 2 2 (1 + |ξ| 2 ) s |v(ξ)| 2 dξ ≤ C(δt)(E(χ 2 n + |β(t n+1 )-β(s)| 2 t n+1 -s ) 2 ) 1/2 E( v 4 H s+2 )
I = E Im R d (1 + |ξ| 2 ) s ēn (ξ)F x E S * n,δt ((δt)T n,δt - t n+1 tn S(t n+1 , s)ds)f k (u(t n )) dξ. Now, for a fixed ξ ∈ R d , F x E S * n,δt ((δt)T n,δt - t n+1 tn S(t n+1 , s)ds) f k (u(t n )) (ξ) = (m 1,n,δt (ξ) + m 2,n,δt (ξ)) f k (u(t n ))(ξ)
where

m 1,n,δt (ξ) = (δt)E (1 -i √ δt 2 χ n |ξ| 2 ) -1 -1 and m 2,n,δt (ξ) = t n+1 tn E(1 -e -2i arctan( √ δt 2 χn|ξ| 2 )+i(β(t n+1 )-β(s))|ξ| 2 )ds.
On the one hand, we use the fact that 1

1 -ix -1 = ix + x 2 ε(x)
for some bounded function ε on R to obtain

(2.48) |m 1,n,δt (ξ)| ≤ C(δt) 2 |ξ| 4 .
On the other hand, we successively use the equalities

arctan x = x + x 3 η(x)
and 1 -e iy = -iy + y 2 ν(y) for some bounded functions η and ν on R, to get

(2.49) |m 2,n,δt (ξ)| ≤ C(δt) 2 |ξ| 4 .
We deduce from (2.47), (2.48), (2.49) and Cauchy-Schwarz inequality that

I ≤ C(δt) 2 E( e n 2 H s ) 1/2 E( f k (u(t n )) 2 H s+4 ) 1/2 .
Because s is a positive integer, one may easily check that there is a constant

C(k) such that f k (v) H s+4 ≤ C(k) v 2 H s+4 , ∀v ∈ H s+4 (R d ), and I ≤ C(δt) 2 E( e n 2 H s ) 1/2 sup t∈[0,T ] (E( u(t) 2 H s+4 )) 1/2 .
It follows

(2.50) E e n,δt,1 , e n,δt,4 s ≤ C(δt) 2 E( e n 2 H s ) 1/2 [ sup t∈[0,T ] E( u(t) 4 H s+4 ) + 1].
The inner product of the first and fifth terms in (2. [START_REF] Milstein | Approximate integration of stochastic differential equations[END_REF]) is easily estimated thanks to the Cauchy-Schwarz inequality, the fact that g k is a Lipschitz function in H s × H s , and the unitarity of the operators S * n,δt and T n,δt in H s (R d ) :

(

E e n,δt,1 , e n,δt,5 s

= (δt)E e n , iS * n,δt T n,δt [g(u n , u n+1 ) -g(u(t n ), u(t n+1 ))] s ≤ C(δt)E( e n 2 H s ) 1/2 E( e n 2 H s + e n+1 2 H s ) 1/2 ≤ C(δt)E( e n 2 H s + e n+1 2 H s ).
Finally, all terms of the form E e n,δt,j , e n,δt,5 with 2 ≤ j ≤ 5 are bounded using Cauchy Schwarz inequality and (2.21), (2.29), (2.32) 2.52), we deduce the existence of constants Kj , for j = 1, 2, 3, depending on k and sup t∈[0,T ] E( u(t) 4 H s+6 ) (which also depends on k) such that for any n with nδt ≤ T ,

E( e n+1 2 H s ) ≤ (1 + K1 δt)E( e n 2 H s ) + K2 (δt) 3 + K3 (δt)E( e n+1 2 H s ).
We deduce that (2.11) holds provided T /N = δt ≤ 1/2 K3 with K 1 = 3 K1 + 2 K3 and K 2 = 2 K2 . The discrete Gronwall Lemma concludes the proof of Proposition 2.6 for δt = T /N ≤ 1/2 K3 . The conclusion is obviously still true for all N < 2T K3 , since the continuous solution u has paths a.s. in C([0, T ]; H s (R d )) and the discrete solution u δt given by Proposition 2.2 has paths a.s. in L ∞ (0, T ; H s (R d )).

2.3.

Proof of Theorem 2.3. We recall that the discrete solution is now defined by u δt n = u k n if nδt ≤ τ k δt , and τ k δt is defined in (2.8), for any k ∈ N. On the other hand, by Proposition (2.6), we have for any k ∈ N, recalling that δt = T /N :

E sup n≤N u k (t n ) -u k n 2 H s ≤ n≤N E( u k (t n ) -u k n 2 H s ) ≤ C(k)T (δt)
and this show that sup n≤N u(t n ) -u δt n H s converges to zero in probability as δt tends to zero : indeed, if ε > 0 is fixed, with ε ≤ 1, then for k ∈ N,

P sup n≤N u(t n ) -u δt n H s ≥ ε ≤ P sup n≤N u(t n ) H s ≥ k + P sup n≤N u k+1 (t n ) -u k+1 n H s ≥ ε ≤ P sup n≤N u(t n ) H s ≥ k + T C(k)(δt)ε -2 .
Hence, choosing first k sufficiently large so that

P(sup n≤N u(t n ) H s ≥ k) ≤ ε/2 - this is possible since u ∈ C([0, T ]; H s (R d )) a.s.
-then δt sufficiently small so that T C(k)(δt)ε -2 ≤ ε/2, the convergence in probability is proved. We deduce from this convergence in probability, taking again δt = T /N for some

N ∈ N, that (2.53) lim K→∞ P sup n≤N u δt n H s ≥ K = 0
uniformly with respect to δt = T /N . Now, we can prove that P( u(t n ) -u δt n H s ≥ Cδt) tends to zero as C tends to infinity, uniformly with respect to δt : indeed, thanks to Markov inequality 2 .

P( u(t n ) -u δt n H s ≥ Cδt) ≤ P sup n≤N u(t n ) H s ≥ k + P sup n≤N u δt n H s ≥ k + E( u k (t n ) -u k n 2 H s ) C 2 (δt)
We conclude choosing first k sufficiently large so that the two first terms are small, thanks to (2.53), then for k fixed, choosing C sufficiently large, so that the last term, which is bounded by C(k)/C 2 thanks to Proposition 2.6, is in turn small.

Numerical simulations in dimension one

In this section, we consider equation (1.1) in space dimension one, i.e.

(3.1) idφ + ∂ 2 x φ • dβ + |φ| 2σ φ = 0, for different values of σ. We numerically illustrate the order of convergence given in Section 2, and we investigate the behavior of the solutions in terms of the power σ of the nonlinear term.

3.1. Numerical schemes. Three different schemes have been used for the numerical simulations. These are classical schemes for the deterministic nonlinear Schrödinger equation, that we have adapted to the simulation of the equation with Stratonovich noise (3.1). For each scheme, we work on the spatial domain [-L x , L x ] and we use a regular mesh in space. We of course need the use of boundary conditions in space, even though our aim is to approximate the solution of equation (3.1) set in R. We will work with either Dirichlet or periodic boundary conditions, and take L x sufficient large so that the solution stays confined in the interval [-L x , L x ] up to the final time of the computations. The spatial mesh step size is δx = 2Lx

Nx where N x = 2M x is an even integer, the time step being δt > 0, and we denote the grid points in space and time by x j = -L x + jδx, j = 0, 1, . . . , N x , and t k = kδt, k = 0, 1, 2 . . . with x 0 = -L x , x Nx = L x . Furthermore, φ k i is the approximations of φ(t k , x i ). Time-splitting spectral discretization. We consider the equation with periodic boundary conditions, and use pseudo-spectral discretization in space. The relations between φ k i and its discrete Fourier coefficients in space φk ℓ are given by (3.2)

φ k i = Mx-1 ℓ=-Mx φk ℓ e iξ ℓ x i , for i = 1, • • • N x , and
(3.3) φ k ℓ = 1 N x Nx i=1 φ k i e -iξ ℓ x i , ℓ = -M x , .., M x -1,
where ξ ℓ = πℓ Lx , for ℓ = -M x , . . . , 0, . . . , M x -1, gives the regular mesh grid for the frequencies.

Equation (3.1) is solved with the classical Strang splitting scheme. One first solves (3.4) idφ + ∂ 2 x φ • dβ = 0, on a time interval of length δt/2. In the second step, one solves (3.5) i∂ t φ + |φ| 2σ φ = 0

φ 0 = φ *
where φ * is the solution of (3.4) at time t = t k + δt/2. The last step consists in solving again (3.4) over a half time step with initial data given by the solution of the second step. Equation (3.4) is discretized in space by the Fourier spectral method and integrated exactly in time. For t ∈ [t k , t k+1 ], multiplying (3.5) by φ, the conjugate of φ and taking the imaginary part, we obtain for x ∈ [-L x , L x ] and t ≥ 0 :

d dt |φ(t, x)| 2 = 0.
With this conservation law, solving (3.5), amounts to solve :

(3.6) i∂ t φ + |φ * | 2σ φ = 0.
Hence the algorithm may be summarized as follows.

For the first step of the splitting, we compute

φ k l = 1 N x Nx i=1 φ k i e -iξ l x i , l = -M x , .., M x -1, φ * l = φ k l exp -i δt 2 ξ 2 l χ k φ * i = Mx-1 l=-Mx φ * l e iξ l x i , for i = 1, • • • N x ,
where δt 2 χ k = β(t k + δt 2 ) -β(t k ), so that (χ k ) k=0,1,... is a family of independent normal random variables, and ( U l ) denote the Fourier coefficients of a periodic vector U = (U i ).

For the second step, we compute

φ j k+1 = φ * j exp iδt|φ * | 2σ . (3.7)
Finally for the last step :

(3.8) φ k+1 i = Mx-1 l=-Mx φ l k+1 exp i ξ l x i - δt 2 ξ 2 l χ k ) , for i = 1, • • • N x .
Crank-Nicolson scheme. We consider the Crank-Nicolson scheme which was analyzed in Section 2, together with finite differences centered discretization in space. We use here Dirichlet boundary conditions which are easily discretized with finite differences. Therefore, ∂ 2

x is approximated by the classical centered finite difference operator D + D -:

(D + D -φ) i = φ i+1 -2φ i + φ i-1 (δx) 2 .
The scheme reads:

(3.9) i φ k+1 -φ k δt + χ k √ δt D + D - φ k+1 + φ k 2 + g(φ k+1 , φ k ) = 0
where here √ δtχ k = β(t k + δt) -β(t k ), and g is defined by (2.4). Note that the computation of φ k+1 requires the resolution of a nonlinear problem for which a fixed point procedure is used.

Crank-Nicolson scheme with relaxation. For the last scheme, we use a fractional-step Crank-Nicolson scheme with relaxation introduced by C. Besse for the deterministic nonlinear Schrödinger equation [START_REF] Besse | Schéma de relaxation pour l'équation de Schrödinger non linéaire et les systèmes de Davey et Stewartson[END_REF]. The advantage is that the scheme is linearly implicit and requires only the resolution of a linear system. It reads in our case :

(3.10) i φ k+1 -φ k δt + χ k √ δt D + D - φ k+1 + φ k 2 + R k+ 1 2 φ k+1 + φ k 2 = 0,
where again √ δtχ k = β(t k + δt) -β(t k ), and the auxiliary function R is given by

R k+ 1 2 + R k-1 2 2 = |φ k | 2σ .
3.2. Numerical results. We here report the results of several numerical experiments. In order to avoid numerical reflections due to the boundaries, we take a domain [-L x , L x ] large enough.

In order to satisfy the required independence properties of the family of Gaussian random variables χ k , these are simulated with the Netlib Library (http://netlib.org) which uses a random generator routine with period 10 32 -1. This has to be compared with the maximum size of the random vector we use which is of order 10 9 .

3.2.1. Strong order of convergence. In the first experiment, we compute the strong numerical order of the three numerical schemes described above and in particular, we recover numerically the result obtained in Theorem 2.3. In the deterministic case, the schemes are known to be consistent and of order 2 in time (see [START_REF] Besse | Schéma de relaxation pour l'équation de Schrödinger non linéaire et les systèmes de Davey et Stewartson[END_REF], [START_REF] Besse | Order estimates in time of splitting methods for the nonlinear Schrödinger equation[END_REF], [START_REF] Delfour | Payre Finite difference solutions of a nonlinear Schrödinger equation[END_REF] ) which is shown in figure 1.

We choose the computational domain [-L x , L x ] = [-30, 30] and an initial data with a Gaussian profile φ(0, x) = exp(-3x 2 ). As our aim is to identify only the numerical order in time, we take a fine spatial mesh with δx = 10 -3 .

Since there is no explicit solution for the stochastic equation (3.1), we first compute a reference solution φ ref on a fine mesh in time δt ref = 10 -5 , that we compare to approximations of the same equation on 6 coarser grids in time by using a time step δt = 2 p δt ref with p = 1, ..., 6. The Brownian path (one trajectory) is kept fixed for each approximation. In figure 1, we plot four convergence curves corresponding to the logarithm of the relative errors |ψ ref -ψ 2 p δt | 2 L 2 at time T = 0.5 for different schemes, in terms of the logarithm of the number of points (in time) used in the discretization. The slopes of the three curves corresponding to the stochastic schemes are close to -1 and we recover that the strong order of the Crank Nicolson scheme is 1 in time. Thus, the result agrees with the theoretical analysis of the previous section.

The same experiment has been done with many other trajectories and the same order of convergence has always been observed.

3.2.2.

Influence of the power nonlinearities σ. We now investigate the influence of the power nonlinearity σ on the behavior of the numerical solution of (3.1). For these experiments, we only use the Crank-Nicolson scheme for the following values of σ : σ = 1, σ = 2, σ = 3, σ = 3.9 and σ = 4, which we conjecture to be the critical case. Indeed, due to the scaling invariance λ -1/2 β(λt) of the Brownian motion β, one may easily see that equation (3.1) with σ = 4 is invariant under the scale change u λ (t, x) = √ λu(λ 4 t, λx) for any λ > 0, ans so is the L 2 -norm in dimension one. This suggests that equation (3.1) with σ = 4 should be critical in L 2 (R). This is different from the deterministic case for which the L 2 -critical power is σ = 2. For each simulation (σ = 1, σ = 2, σ = 3 and σ = 3.9), we take a Gaussian initial state given by φ(0, x) = exp(-3x 2 ) on the computational domain [-L x , L x ] = [-10π, 10π] with N x = 2 12 points in space and a time step δt = 0.4δx 2 with a final time T which is chosen in such a way that the numerical solution stays confined in [-L x , L x ] on the whole interval [0, T ]. In order to understand the noise effects, all the experiments are done with one realization of a Brownian path. Note that all the different simulations that we have performed have given the same qualitative behavior.

Figure 2 shows the profile of the space-time evolution of the solution in the case σ = 1. We can see clearly the effect of the noise on the dispersion. It seems that dispersion and nonlinearity stay well balanced for a while, resulting in alternating contractions and expansions of the solution. Here the solutions stays confined until T = 5, then boundary effects occur (which would suggest that dispersion finally wins) and we had to stop the computation. The modulus of the Fourier modes of the solution are plot in figure 6. In the cases σ = 2, σ = 3 and σ = 3.9 (see figure 3, figure 4 and figure 5), the behavior of the solution in time is similar, although the final computation time may be different from one case to the other. The final amplitude seems much smaller than the initial amplitude although the modulus of the Fourier modes of the solution do not seem to move with time (see figure 6 and figure 7 ), which indicates that after a while there is no energy transfer between the different Fourier modes.

For the case σ = 4, which is the expected critical case (see the above comments), it is much more difficult to perform numerical simulations since the numerical method has to be consistent with the small spatial scales of the blow up structures. Here we do not use spatial local refinement as in [START_REF] Debussche | DiMenza Numerical simulation of focusing stochastic nonlinear Schrödinger equations[END_REF] . Since we are in dimension 1 in space, and it is not too expensive, we only refine globally the time-space discretization. We take here δx = 10 -4 and δt = 0.1δx 2 . We take a Gaussian inital data φ(0, x) = 1.4 exp(-3x 2 ). In figure 8 we plot the modulus of the initial state and of the solution at time t = 0.035. Figure 9 shows the modulus of the spatial spectrum of the solution at three different times; there seems to be a collapse of this spatial spectrum at time t = 0.035. Finally, figure 10 represents the evolution with time of the discrete H 1 -norm, which indicates that the solution reaches very high amplitudes.

We finally plot in figure 11 the evolution with time of the expected amplitude of the solution for σ = 1. We use the notation < > for the empirical average which is the numerical approximation of the the average amplitude E(max x |φ(t, x)|) where N is number of computed trajectories. For the numerical simulation we take a Gaussian inital data given by φ(0, x) = exp(-3x 2 ) on the computational domain [-L x , L x ] = [-20π, 20π] with N x = 2 12 points in space and a time step δt = 0.4δx 2 with a final time T = 1. The average is taken over 1000 trajectories. There seems to be a rather strong diffusion of the expected amplitude as time evolves, but no clear time scale for the diffusion phenomenon can be extracted from the numerical simulations.

Proposition 2 . 1 .

 21 Let d = 1, σ = 1 or 2, and u 0 ∈ H s (R) with s ≥ 1. Then there exists a unique adapted solution of equation (2.1) with paths a.s. in C(R + ; H s (R)).

Figure 1 .

 1 Figure 1. Strong order of convergence

Figure 2 .

 2 Figure 2. σ = 1 : space-time evolution of the solution (left) and contour plot of the solution (right).

Figure 3 .

 3 Figure 3. σ = 2 : space-time evolution of the solution (left) and contour plot of the solution (right).

Figure 4 .

 4 Figure 4. σ = 3 : space-time evolution of the solution (left) and contour plot of the solution (right).

Figure 5 .

 5 Figure 5. σ = 3.9 : space-time evolution of the solution (left) and the contour plot of the solution (right).

Figure 6 .

 6 Figure 6. Space-time evolution of the modulus of the spatial spectrum of the solution; left : σ = 1, right : σ = 2.

Figure 7 .Figure 8 .

 78 Figure 7. Space-time evolution of the modulus of the spatial spectrum of the solution; left : σ = 3, right : σ = 3.9.

Figure 9 .

 9 Figure 9. σ = 4. Evolution of the spatial spectrum of the solution at times t = 0, t = 0.00012 and t = 0.035.

Figure 10 . σ = 4 :

 104 Figure 10. σ = 4 : Evolution of the discrete H 1 norm.

Figure 11 .

 11 Figure 11. σ = 1 : Time evolution of < max x |φ(t, x)| > .

  Let us estimate the expectation of the product of the first and second terms in(2.19). Using the expressions of those terms given in (2.14) and (2.15), together with Plancherel equality and the Fourier transform of the operators S n,δt and S(t n+1 , t n ), we obtain E(e n,δt,1 , e n,δt,2 ) s = E S n,δt e n , (S n,δt -S(t n+1 , t n ))u(t n ) s = E e n , S * n,δt (S n,δt -S(t n+1 , t n

			and (2.32), for any j with 2 ≤ j ≤ 4,
	(2.34)	E( e n,δt,j	2 H

s ) ≤ C(δt) 3 1 + sup t∈[0,T ] E( u(t) 4 H s+6 ) .

  1/2 by Cauchy-Schwarz inequality. In order to estimate the term E e n,δt,1 , e n,δt,4 s , we use (2.14), (2.17) and a Taylor expansion to decompose it as E e n,δt,1 , e n,δt,4 s= I + II + III n+1 , s)ds D v g k (u(t n ), u(t n ))E n (t n+1 ) s (u(t n ), u(t n ))(E n (t n+1 )) 2 dθ swhere E n is defined in (2.25). The two last terms are bounded as above, using Cauchy-Schwarz inequality and the fact that (δt)T n,δt -n+1 , s)ds is of order (δt) 3/2 while E n (t n+1 ) is of order (δt) 1/2 (see(2.27) and (2.41)). The first term is bounded thanks to arguments similar to those we have used to obtain (2.35). Indeed, Plancherel Theorem implies(2.47)

			t n+1
	tn S(t and
	III = E e n , iS * n,δt (δt)T n,δt -	t n+1 tn	S(t n+1 , s)ds	1 vv g k t n+1 0 (1-θ)D 2 tn S(t

with

I = E e n , iS

* n,δt (δt)T n,δt -t n+1 tn S(t n+1 , s)ds f k (u(t n )) s , II = E e n , iS * n,δt (δt)T n,δt -

Acknowledgement : This work was supported by the ANR project STOSYMAP (ANR-2011-BS01-015-03).