François Dupuy

Abdelkader Ameraoui

Kamal Boukhetala

Jean-François Dupuy

Bayesian estimation of the tail index of a heavy tailed distribution under random censoring

Keywords: Extreme value theory, asymptotic properties, maximum a posteriori estimator, mean posterior estimator, simulations

come

Bayesian estimation of the tail index of a heavy tailed distribution under random censoring

Abdelkader Ameraoui, Kamal Boukhetala, Jean-

Introduction

Tail index estimation is one of the most important issues in extreme value theory. The tail index measures the thickness of the tail of a probability distribution function and thus plays a crucial role for evaluating the risk of occurrence of extremes events. In particular, estimation of the tail index constitutes usually a first step in an extreme value analysis. A vast literature has been dedicated to this topic. Recent overviews can be found in the monographs [START_REF] Beirlant | Statistics of Extremes[END_REF] and [START_REF] De Haan | Extreme Value Theory[END_REF].

Let F be the cumulative distribution function (cdf) of some non-negative random variable Y . We assume that F is heavy-tailed, that is, there exists a constant α > 0 such that

1 -F (x) = x -α (x), (1)
where is a slowly varying function at infinity:

lim x→∞ (tx) (x)
= 1 for all t > 0.

If (1) holds, we have:

lim x→∞ 1 -F (tx) 1 -F (x) = t -α
for all t > 0 and we say that F is regularly varying at infinity with tail index α, which we denote by F ∈ R -α . The positive number γ := α -1 is called the extreme value index (EVI) of F . The conditions above amount to assuming that the distribution function F is in the max-domain of attraction of a Fréchet distribution. Such distribution functions are useful in practice for investigating phenomena where exceptional values have a significant occurrence frequency. Examples include the number of claims in insurance [START_REF] Embrechts | Modelling Extremal Events[END_REF], transmission times in telecommunications [START_REF] Resnick | Heavy-Tail Phenomena[END_REF], log-returns of price speculation [START_REF] Embrechts | Modelling Extremal Events[END_REF]. Several estimators have been proposed for the tail index α, or equivalently, for the EVI γ, such as Pickands estimator [START_REF] Pickands | Statistical Inference Using Extreme Order Statistics[END_REF], Dekkers et al. (or moment) estimator [START_REF] Dekkers | A moment estimator for the index of an extremevalue distribution[END_REF] and Hill's estimator [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF], which is the most popular estimator of γ in model [START_REF] Andersen | Statistical Models Based on Counting Processes[END_REF]. Let X 1 , X 2 , . . . , X n be independent and identically distributed (iid thereafter) random variables with common cdf F . Let k ∈ {2, . . . , n} and X n,1 ≤ X n,2 ≤ . . . ≤ X n,n be the ordered statistics associated to the sample X 1 , X 2 , . . . , X n . Hill's estimator is defined as

H(k) := 1 k k i=1 log(X n,n-i+1 ) -log(X n,n-k ). (2)
Consistency of Hill's estimator was proved in [START_REF] Mason | Laws of large numbers for sums of extreme values[END_REF] under the regular variation condition [START_REF] Andersen | Statistical Models Based on Counting Processes[END_REF]. Its asymptotic normality was further established under an additional condition known as the second-order regular variation condition (see [START_REF] De Haan | Extreme Value Theory[END_REF], page 117).

In this paper, we address estimation of the tail index α when observations X 1 , X 2 , . . . , X n are randomly right-censored. Censoring commonly occurs in the analysis of event time data. For example, X may represent the duration until the occurrence of some event of interest, such as death of a patient. If a patient is still alive or has dropped out of the study for some reason when the data are collected, the variable of interest X is not available. An appropriate way to model this situation is to introduce a random variable Y (called a censoring random variable) such that the observations consist of pairs (Z i , δ i ), 1 ≤ i ≤ n where Z i = min(X i , Y i ), δ i = 1 {Xi≤Yi} and 1 is the indicator function. Estimation of the EVI with censored data was considered in [START_REF] Beirlant | Estimation of the extreme value index and extreme quantiles under random censoring[END_REF], [START_REF] Brahimi | On the asymptotic normality of Hill's estimator of the tail index under random censoring[END_REF], [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF], [START_REF] Gomes | Estimation of the extreme value index for randomly censored data[END_REF] and [START_REF] Gomes | Censoring estimators of a positive tail index[END_REF]. For example, [START_REF] Beirlant | Estimation of the extreme value index and extreme quantiles under random censoring[END_REF] [START_REF] Beirlant | Estimation of the extreme value index and extreme quantiles under random censoring[END_REF] proposed to estimate γ by the following modified version of Hill's estimator:

H C (k) := k i=1 log(Z n,n-i+1 ) -log(Z n,n-k ) k i=1 δ [n-i+1] , (3)
where Z n,1 ≤ Z n,2 ≤ . . . ≤ Z n,n are the ordered statistics associated to the censored sample Z 1 , Z 2 , . . . , Z n and δ [n-i+1] is the concomitant value of δ associated with Z n,n-i+1 . More generally, Einmahl et al. (2008) [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF] and Gomes and Neves (2011) [START_REF] Gomes | Estimation of the extreme value index for randomly censored data[END_REF] proposed to estimate γ by γ(C)

Z (k) = γZ (k) p , (4)
where γZ (k) is any of the classical EVI estimators calculated on the censored observations Z 1 , . . . , Z n and p = 1

k k i=1 δ [n-i+1]
is the proportion of uncensored values in the k largest observations of Z. Obviously, the estimator (4) coincides with the adapted Hill's estimator (3) when γZ (k) is Hill's estimator [START_REF] Beirlant | Statistics of Extremes[END_REF]. In this paper, we adopt a completely different approach and investigate Bayesian estimation of the tail index α := γ -1 .

Bayesian estimation provides an alternative to frequentist methods. In the context of extreme values analysis without censoring, Bayesian estimators have been investigated in [START_REF] Cabras | A Bayesian approach for estimating extreme quantiles under a semiparametric mixture model[END_REF], [START_REF] Coles | Bayesian methods in extreme value modelling: a review and new developments[END_REF], [START_REF] Diebolt | Quasi-conjugate bayes estimates for gpd parameters and application to heavy tails modelling[END_REF], [START_REF] Do Nascimento | A semiparametric Bayesian approach to extreme value estimation[END_REF], [START_REF] Stephenson | Bayesian inference for extremes: accounting for the three extremal types[END_REF], [START_REF] De | Parameter estimation of the generalized Pareto distribution-Part II[END_REF]. See also [START_REF] Beirlant | Statistics of Extremes[END_REF] (chapter 11). Applications include operational risk modeling [START_REF] Ergashev | A Bayesian approach to extreme value estimation in operational risk modeling[END_REF], hydrology [START_REF] Liang | Application of Bayesian approach to hydrological frequency analysis[END_REF] and market indices modeling [START_REF] So | Bayesian analysis of tail asymmetry based on a threshold extreme value model[END_REF]. But to the best of our knowledge, no Bayesian estimator of the tail index α has been proposed when censoring is present. The present work intends to fill this gap.

We construct several Bayesian estimators of α in model [START_REF] Andersen | Statistical Models Based on Counting Processes[END_REF]. Bayesian estimation requires specifying a prior distribution for the unknown parameter. We investigate various priors (Jeffreys, maximal data information and conjugate Gamma priors), leading to several maximum a posteriori and mean posterior estimators of α. We establish consistency of these estimators and asymptotic normality of their posterior distribution. Their finite-sample performance are assessed via simulations. Tail index estimation requires choosing an appropriate threshold for defining excesses. We propose a randomization method for tackling this issue. This procedure is evaluated by simulations. Finally, we illustrate our methodology on a set of AIDS survival data.

The paper is organized as follows. In Section 2, we construct our estimators and we establish their asymptotic properties. All proofs are deferred to an appendix. Section 3 reports the results of a comprehensive simulation study. The randomization procedure of the threshold is proposed and assessed by simulations in Section 4. Application to AIDS data is carried out in Section 5. A discussion and some perspectives are given in Section 6.

Bayesian estimation of the tail index

We first introduce some notations. Let X 1 , X 2 , . . . , X n be n iid copies of a non-negative random variable X with cdf F . The probability density function of X is denoted by f . We assume that F is heavy-tailed with tail index α, i.e. F ∈ R -α . Let Y 1 , Y 2 , . . . , Y n be an iid sample of the non-negative censoring random variable Y with cdf G and density function g. We assume that G is also heavy-tailed, with positive tail index β, i.e. G ∈ R -β . We assume that X and Y are independent and that we observe the n independent pairs

(Z i , δ i ), 1 ≤ i ≤ n,
where Z i = min(X i , Y i ) and δ i = 1 {Xi≤Yi} . We let Z n,1 ≤ Z n,2 ≤ . . . ≤ Z n,n be the ordered statistics associated to the sample (Z 1 , Z 2 , . . . , Z n ) and δ [n-i+1] be the concomitant value of δ associated with

Z n,n-i+1 , that is, δ j = δ [n-i+1] if Z j = Z n,n-i+1
. Let H be the cdf of Z. Note that by independence of X and Y , H is also heavy-tailed and H ∈ R -(α+β) . Now, we define the random variables E j,u = Zj u given Z j > u, which represent relative excesses over the threshold u. Let N u denote the number of such excesses. Then we can obtain the partial likelihood of α based on the sample (E, ∆) = (e j,u , δ [n-Nu+j] ) j=1,...,Nu (see [START_REF] Andersen | Statistical Models Based on Counting Processes[END_REF], [START_REF] Beirlant | Estimation of the extreme value index and extreme quantiles under random censoring[END_REF]):

L(E, ∆|α) = Nu j=1 αe j,u -(α+1) δ [n-Nu +j] e j,u -α 1-δ [n-Nu +j] = α Nu j=1 δ [n-Nu +j]   Nu j=1 e j,u   -α   Nu j=1 e j,u -δ [n-Nu+j]   . (5)
In the Bayesian framework, unknown parameters are considered as random variables whose probability density function (the so-called prior density) represent the initial set of beliefs about the parameters. The prior distribution is updated by using the information contained in the data, yielding the so-called posterior density of the parameters. In our setting, we provide the unknown tail index α with a prior density π(α). Then, using Bayes' theorem, we obtain the posterior density π(α|E, ∆) of α, which writes as:

π(α|E, ∆) = L(E, ∆|α)π(α) Λ L(E, ∆|α)π(α)dα
, where Λ is the support of the distribution of α and the denominator Λ L(E, ∆|α)π(α)dα is a normalization constant independent of α. The posterior distribution of α is proportional to the product of the partial likelihood (5) and the prior, namely:

π(α|E, ∆) ∝ L(E, ∆|α)π(α).
Choosing the prior density is a central issue in Bayesian estimation. When the information available for prior elicitation is minimal, one can use objective (or non-informative) priors, such as Jeffreys prior [START_REF] Jeffreys | Theory of Probability[END_REF] and the maximal data information (MDI) prior introduced by Zellner (1971) (see [START_REF] Zellner | An Introduction to Bayesian Inference in Econometrics[END_REF] and [START_REF] Zellner | Models, prior information, and Bayesian analysis[END_REF]). On the other hand, when prior choice can be based on experts opinion, so-called subjective prior distributions can be used, such as conjugate priors. Conjugate priors are such that the posterior distribution of the parameter is in the same family as the prior. In Subsections 2.1 and 2.2, we investigate MDI and Jeffreys priors for the estimation of the tail index α under random censoring. In Subsection 2.3, we investiagte a conjugate Gamma prior.

Maximal data information (MDI) prior

The MDI prior was defined so as to maximize the average information in the data density relative to that in the prior (see [START_REF] Zellner | An Introduction to Bayesian Inference in Econometrics[END_REF]). If L(E, ∆|α) is the likelihood of a single observation (E, ∆), then π(α) ∝ exp [E(log L(E, ∆|α)]. In our setting, the MDI prior for α is:

π(α) ∝ exp log(α) - 1 α . (6)
Using ( 5) and ( 6), the posterior density of α based on MDI prior is given by: (α|E, ∆)

π(α|E, ∆) ∝ α Nu j=1 δ [n-Nu+j] exp   -α
∝ (1 + b u ) log(α) -a u α - 1 α .
From this, we can define the maximum a posteriori (or MAP) tail index estimator as:

α(MDI) M AP = arg max α (α|E, ∆) := B + √ B 2 + 4A 2A , (7)
where A = a u and B = 1 + b u . Letting u = Z n,n-k in [START_REF] Deheuvels | Almost sure convergence of the Hill estimator[END_REF], we obtain the following formal definition of our MAP estimator of α under MDI prior:

Definition 2.1. Let k := k n be a sequence such that k → ∞ and k n → 0 as n → ∞.
Then the MAP estimator of the tail index α under MDI prior is given by: α(MDI)

M AP := B + √ B 2 + 4A 2A , (8)
where

A = k i=1 log Z n,n-k+i Z n,n-k and B = 1 + k i=1 δ [n-k+i] .
Remark 2.1. In order to simplify notations, the dependency of α(MDI) M AP on k and n is not made explicit. This convention will also be adopted for the other estimators defined below.

Jeffreys prior

A second kind of non-informative prior is Jeffreys prior [START_REF] Jeffreys | Theory of Probability[END_REF], which has an interesting invariant reparametrization property. This prior is proportional to the square root of Fisher's information. If L(E, ∆|α) is the likelihood of a single observation (E, ∆), Jeffreys prior can be written as:

π(α) ∝ -E ∂ 2 ∂α 2 log L(E, ∆|α) 1 2 ∝ 1 α . (9)
Using ( 5) and ( 9), the posterior density of α based on Jeffreys prior is given by:

π(α|E, ∆) ∝ α ( Nu j=1 δ [n-Nu +j] )-1 exp   -α Nu j=1 log e j,u   ,
which coincides (up to some normalizing constants) with the probability density function of the Gamma distribution Γ(

Nu j=1 δ [n-Nu+j] ,
Nu j=1 log e j,u ). Based on this posterior, we construct two classical Bayesian estimators of α, namely the mean posterior estimator (MPE) and the maximum a posteriori (MAP) estimator.

• mean posterior estimator (MPE): this estimator is defined as the mean of the posterior distribution of α:

α(J) M P E := Nu j=1 δ [n-Nu+j] Nu j=1 log e j,u . (10)
Letting u = Z n,n-k in [START_REF] Diebolt | Quasi-conjugate bayes estimates for gpd parameters and application to heavy tails modelling[END_REF], we obtain the following formal definition of our MP estimator of the tail index α:

Definition 2.2. Let k := k n be a sequence such that k → ∞ and k n → 0 as n → ∞.
Then the MP estimator of the tail index α under Jeffreys prior is given by: α(J)

M P E := k i=1 δ [n-i+1] k i=1 log Z n,n-i+1 Z n,n-k . ( 11
)
We note that ( α(J) M P E ) -1 coincides with the estimator of the EVI γ proposed by Einmahl [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF] when censoring occurs.

• maximum a posteriori estimator (MAP): this estimator is defined as the mode of the posterior distribution of α:

α(J) M AP := arg max α π(α|E, ∆) = Nu j=1 δ [n-Nu+j] -1 Nu j=1 log e j,u . (12)
Letting u = Z n,n-k in [START_REF] Embrechts | Modelling Extremal Events[END_REF], we obtain the following formal definition of our MAP estimator of the tail index α under Jeffreys prior:

Definition 2.3. Let k := k n be a sequence such that k → ∞ and k n → 0 as n → ∞.
Then the MAP estimator of the tail index α under Jeffreys prior is given by: α(J)

M AP := k i=1 δ [n-i+1] -1 k i=1 log Z n,n-i+1 Z n,n-k . ( 13
)

Conjugate prior

In Subsections 2.1 and 2.2, prior elicitation is essentially data-driven since Jeffreys and MDI priors only take account of the information provided by the observations for estimating the tail index α. On the other hand, subjective priors (such as conjugate priors) attempt to bring prior knowledge about the phenomenon under study into the problem. Conjugate priors are particularly useful due to their mathematical convenience. Here, we propose a conjugate prior from the Gamma family Γ(•, •). If α ∼ Γ(a, b) (with a, b > 0), we can easily obtain the posterior distribution of α given the observations (E, ∆), as:

α|(E, ∆) ∼ Γ   Nu j=1 δ [n-Nu+j] + a, Nu j=1 log e j,u + b   .
Based on this, we propose two estimators of α under a conjugate Gamma prior. Definition 2.4. Let k := k n be a sequence such that k → ∞ and k n → 0 as n → ∞. Let the prior density π(α) be a Γ(a, b) distribution. Then:

1. the MP estimator of the tail index α is given by: α(c)

M P E := k i=1 δ [n-i+1] + a k i=1 log Z n,n-i+1 Z n,n-k + b , (14)
2. the MAP estimator of the tail index α is given by:

α(c) M AP := k i=1 δ [n-i+1] + a -1 k i=1 log Z n,n-i+1 Z n,n-k + b . (15)

Asymptotic properties of the Bayesian estimators

In this section, we investigate asymptotic properties of the proposed estimators. We first establish their strong consistency. Proof is given in Appendix A.

Theorem 2.1. Let k := k n be a sequence such that k → ∞ and k n → 0 as n → ∞. Then the estimators α(MDI) M AP , α(J) M P E , α(J) M AP , α (c)
M P E and α(c) M AP given by ( 8), ( 11), ( 13), ( 14) and ( 15) respectively converge almost surely to α as n → ∞. Now, we discuss the limiting behavior of the posterior distribution of the proposed estimators. We need some further notations. Let

L k (α) = log α k i=1 δ [n-k+i] -α k i=1 log e [n-k+i] - k i=1 δ [n-k+i] log e [n-k+i] be the log-likelihood L(E, ∆|α) calculated with the threshold u = Z n,n-k . Let I k (α) = -1 k ∂ 2 ∂α 2 L k (α), I(α) = lim k→∞ I k (α)
and α 0 denote the "true" tail index value. We assume that α 0 belongs to an open interval A = (a 1 , a 2 ) of R, where 0 < a 1 < a 2 < ∞. Finally, let p := P(X < Y ) be the probability that X is uncensored and d -→ denote convergence in distribution. Then the following holds:

Theorem 2.2. Let k := k n be a sequence such that k → ∞ and k n → 0 as n → ∞.
Let π(α) be a prior density with finite expectation and such that π(•) is continuous and positive at α 0 . Let α• be any of the estimators α(MDI) M AP , α(J)

M P E , α(J) M AP , α(c) M P E , α(c) M AP . Then as n → ∞, √ k (α 0 -α• ) d -→ N 0, 1 pI(α 0 )
.

Proof is given in Appendix A.

Simulation study

In this section, we assess the finite-sample performance of the proposed estimators ( 8), ( 11), ( 13), ( 14) and ( 15) of the tail index of a heavy-tailed distribution under random censoring. More precisely, we assess influence of the censoring proportion and threshold choice on the behaviour of our estimators.

Study design. The simulation design is as follows. Let X and Y be independent random variables with cdf F and G respectively. X is the variable of interest and Y is the censoring random variable. We assume that F and G are heavy-tailed with tail index α > 0 and β > 0 respectively. Let p := P(X < Y ) be the probability that X is uncensored. Some straightforward calculations yield β ≈ α 1-p p (e.g., [START_REF] Gomes | Estimation of the extreme value index for randomly censored data[END_REF]). Our simulation process is as follows:

1. We simulate a sample of n independent copies of (Z, δ), where Z = min(X, Y ), δ = 1 {X≤Y } and the censoring proportion is 1 -p. We consider various values for p, namely p = 0.9, 0.75, 0.5, 0.25 (that is, we allow the percentage of censoring in the right tail of X to be 10%, 25%, 50% and 75% respectively).

2.

For each p, we compute the various proposed estimates, by incrementing the threshold (or fraction level) k n from kmin n to kmax n . For notational simplicity, we use the following notations:

• α1 (k) is the MPE estimator (11) in Definition 2.2, • α2 (k) is the MAP estimator (13) in Definition 2.3, • α3 (k) is the MAP estimator (8) in Definition 2.1, • α4 (k) is the MPE estimator (14) in Definition 2.4,
• α5 (k) is the MAP estimator (15) in Definition 2.4.

3.

Steps 1-2 are repeated m times, so that we obtain m realisations of each α (k), for each p.

4.

For each p, we calculate the empirical bias and RMSE (Root Mean Square Error) of each α (k) (for k = k min , . . . , k max and = 1, 2, 3, 4, 5) over the m estimates.

We consider the following two simulation settings:

• X and Y are both distributed as Fréchet random variables, with cdf F (x) = exp(-x -α ) and

G(x) = exp(-x -β ) respectively,
• X and Y are both distributed as Generalized Pareto Distribution (GPD) with cdf

F (x) = 1 -(1 + 1 α x) -α and G(x) = 1 -(1 + 1 β x) -β respectively.
In either case, we take α = 0.5 and α = 0.75 and for each p = 0.9, 0.75, 0.5, 0.25, we take β = α 1-p p . The simulations are conducted using the statistical software R [START_REF]R: A language and environment for statistical computing[END_REF]. Results are provided for a sample size n = 1000 and m = 1000 simulated samples, with k min = 20 and k max = 250.

Results for Fréchet model. For each = 1, . . . , 5, we plot the empirical bias and RMSE of α as functions of k for p = 0.9, 0.75, 0.5, 0.25 (Figure 1: α = 0.5 and Figure 2 : α = 0.75). From these figures, it appears as expected that the bias and RMSE of the proposed estimators decrease when the censoring proportion decreases. The MAP estimator under MDI prior α3 (see [START_REF] Dekkers | A moment estimator for the index of an extremevalue distribution[END_REF], Definition 2.1) might be regarded as the best among all estimators α . In particular, this estimator outperforms all others in terms of RMSE for almost every k.

In Tables 1 and2, we report the averaged (over the m simulated samples) value, empirical bias and RMSE of each α at the optimal fraction level k opt = arg min k RM SE[α (k)] (Table 1: α = 0.5 and Table 2: α = 0.75). These tables confirm the superiority of the MAP estimator under MDI prior α3 . In particular, at the optimal fraction level k opt , α3 appears to be less biased than all other estimators, for every p.

Results for the GPD model. Results for the GPD are given in Appendix B. The conclusions are similar to conclusions for Fréchet model.

Randomizing the fraction level

In the previous section, the (n -k)-th order statistic Z n,n-k was used as threshold in the excesses e j,u and k (or equivalently, the fraction level k/n) was chosen so as to minimize the RMSE. Obviously, this procedure cannot be applied in a practical data analysis since the true tail index α used to calculate the RMSE is unknown. Thus, in this section, we develop an alternative approach for tackling this threshold issue.

We propose to consider the fraction level t := k/n as a random quantity. As k/n lies in ]0, 1[, we suggest to use a Beta distribution for t, with probability density function:

π u (t) = 1 B(a, b) (t -t min ) a-1 (t max -t) b-1 (t max -t min ) a+b-1 1 {tmin t tmax} , (16)
where B(a, b) :=

1 0 s a-1 (1 -s) b-1 ds (a > 0, b > 0)
is the Beta function and t min and t max are suitable bounds for the fraction level. Then we propose the following Monte Carlo approach for infering on α from a sample of observations (Z i , δ i ) i=1,...,n .

First, we simulate x ∼ Beta(a, b) and we calculate t x = (t max -t min )x + t min . Then we obtain k x = [nt x ] (where [ • ] denotes the integer part) and we calculate the estimate α with k = k x . This procedure is repeated N times, which yields N Monte Carlo realisations α(1) , . . . , α(N) of α . Finally, α can be estimated from these N realisations (by taking their empirical mean, median or mode, for example).

We illustrate this procedure in a short simulation study. A sample of size n = 1000 of (Z, δ) is obtained by simulating n observations of X and Y from Pareto distributions with tail index α = 2 and β = α 1-p p respectively, where p = 0.25, 0.5, 0.75 and 0.9 (that is, the censoring fraction is 75%, 50%, 25% and 10% respectively). Based on results of Section 3, we restrict our attention to objective priors. Hence, the proposed randomization procedure is applied to the MP estimator of α under Jeffreys prior (i.e., α1 ) and to the MAP estimator under Jeffreys and MDI priors ( α2 and α3 respectively). We take t min = 0.01 and t max = 0.4. Note that under these values, k is allowed to range from 10 to 400, with a higher concentration of values around k = 300. This is a reasonable choice in view of the sample size. We also take a = 4 and b = 1.5 (a procedure for choosing a and b will be discussed in Section 5). Finally, we let N = 10000.

For each p ∈ {0.25, 0.5, 0.75, 0.9} and l ∈ {1, 2, 3}, we obtain the empirical mean, median and mode of the N Monte Carlo estimates α(j) , j = 1, . . . , N (see Table 3) and we plot their histogram and estimated density (see Figure 3). From Figure 3, histograms are concentrated around the true tail index value α = 2 (represented by the vertical dotted line) in every case. Moreover, the empirical mean, median and mode of the α(j) , j = 1, . . . , N appear to provide satisfactory approximations of α, for every p and = 1, 2, 3.

Overall, combining the proposed Bayesian estimators with the randomization tool described above seems to provide a relevant approach for estimating α. This approach is now illustrated on a real data set.

Jeffreys prior (MPE)

An application to AIDS survival data

In this section, we illustrate our methodology on a set of AIDS survival data. The dataset contains n = 2843 patients diagnosed with AIDS in Australia before 1 July 1991. The source of the data is Dr P.J. Solomon and the Australian National Centre in HIV Epidemiology and Clinical Research (see [START_REF] Venables | Modern Applied Statistics with S[END_REF]). The data are available in the R package "MASS" [START_REF]R: A language and environment for statistical computing[END_REF]. Information on each patient includes gender, date of diagnosis, age at diagnosis, date of death or end of observation and an indicator which equals 1 if the patient died and 0 otherwise. 1761 patients died. The other survival times are right-censored. Our objective is to estimate the tail index α of the cdf F of the survival time. We calculate our estimators on the whole dataset of patients, without distinction of sex or age.

Considering only objective priors for α, we calculate the MP and MAP estimators under Jeffreys prior ( α1 and α2 respectively) and the MAP estimator under MDI prior ( α3 ). To choose the fraction level k/n, we first rely on a method proposed by [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF] [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF]. Then we use the randomizing approach described in Section 4. [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF] [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF] propose to plot the proportion p(k) = Three distinct phases can be distinguished on Figure [START_REF] Brahimi | On the asymptotic normality of Hill's estimator of the tail index under random censoring[END_REF]. In phase (A), the behaviour of k → p(k) is somewhat erratic. Then, there is stable part (phase (B), where k ranges approximately from 70 to 170 and the average p(k) is approximately 0.261) and an increasing part (phase (C)). Einmahl et al. (2008) [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF] suggest to estimate α based on the k-values of phase (B). Indeed, the estimate of α should be quite stable within this range of k. On Figure [START_REF] Cabras | A Bayesian approach for estimating extreme quantiles under a semiparametric mixture model[END_REF], we plot the α1 (k), α2 (k) and α3 (k) as functions of k. We clearly distinguish the same three phases (A), (B) and (C) as in Figure [START_REF] Brahimi | On the asymptotic normality of Hill's estimator of the tail index under random censoring[END_REF]. All three estimates appear to be close to each other and quite stable in phase (B). The average value of α (k) on phase (B) is approximately 1.26 (for = 1, 2, 3). Based on Einmahl et al. [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF] empirical methodology, one may thus retain 1.26 as an estimate of α.

Australian Aids Survival Data

Australian Aids Survival Data

We now apply our randomization procedure. We take t min = 0.01 and t max = 0.4, so that k ranges from 28 to 1137 (this range includes and goes far beyond the stable phase (B)). We consider several values for the parameters a, b > 0 of the Gamma distribution used to randomize the fraction level, namely (a, b) ∈ S := {1, 1.5, 2, 2.5, 3, 3.5, 4} × {1, 1.5, 2, 2.5, 3, 3.5, 4}. For every (a, b) ∈ S, we apply the proposed randomization procedure with N = 10000 and we the standard deviation (Std.) of the resulting N Monte Carlo estimates α(j) , j = 1, . . . , N . Then, we retain the value (a * , b * ) such that the standard deviation is minimum. Finally, we propose to infer on α from the α(j) , j = 1, . . . , N obtained with (a, b) = (a * , b * ).

Table 4 reports standard deviations for α1 . The minimum is achieved when (a, b) = (4, 3) (Std. = 0.117). Standard deviations for the pairs (4, 2) (Std. = 0.120) and (4, 2.5) (Std. = 0.118) are close to the minimum, thus we also retain these values for making inference on α. The minimum standard deviations for α2 and α3 are achieved for the same values of (a, b). The corresponding tables are omitted for conciseness. We can now infer on α from the α(j) , j = 1, . . . , N obtained with (a * , b * ) (by taking their empirical mean, median or mode, for example). We report some summary statistics for α1 , α2 and α3 when (a * , b * ) = (4, 2), (4, 2.5), (4, 3) (see Table 5). We also plot the corresponding histograms of the α(j) , j = 1, . . . , N (see Figure 6). From Figure 6, histograms are quite symmetric and show limited spread. Thus, it seems reasonable to estimate α by any of the usual measures of central tendency. Moreover, from Table 5, the empirical means and medians of the α ( = 1, 2, 3) are close to each other, for every (a, b) = (4, 2), (4, 2.5), [START_REF] Brahimi | On the asymptotic normality of Hill's estimator of the tail index under random censoring[END_REF][START_REF] Beirlant | Estimation of the extreme value index and extreme quantiles under random censoring[END_REF]. From these results, one may retain 1.4 as a reasonable and consensual (over the α , = 1, 2, 3) approximation of α.

Std

Conclusion and perspectives

In this paper, we address estimation of the tail index of a heavy-tailed distribution when data are randomly right-censored. We constructed two kinds of Bayesian estimators, namely mean posterior and maximum a posteriori estimators, for various objective and subjective prior distributions. Our simulation results indicate better performance of estimators based on objective MDI and Jeffreys priors. We also proposed an original automatic procedure for selecting the threshold (or fraction level k/n) required to estimate the tail index. The whole procedure was applied to a set of real data and provided coherent results. Now, several issues still deserve attention. First, we mentioned in introduction that estimation of the tail index constitutes usually a first step in an extreme value analysis. Extreme quantile estimation often comes as a second step. Extreme quantiles of the distribution of X are quantities of the form

F ← (1 -p) = inf{y : F (y) ≥ 1 -p},
where p is so small that this quantile falls beyond the range of the observed data X 1 , . . . , X n . Estimation of extreme quantiles has become a crucial issue in many domains (e.g., determination of the value at risk in financial risk management, determination of the return period of extreme precipitation in meteorology. . . ). We are currently working on Bayesian estimation of extreme quantiles from censored data.

In practice, it often arises that some covariate information W is available to the investigator and the distribution of X depends on W . In this case, interest turn to estimation of the conditional tail index α(w). Adapting our Bayesian approach to this setting also constitutes a topic for our future research.

Appendix A. Proofs.

In this appendix, we outline the proofs of Theorem We first establish asymptotic normality of α(J) M P E . Proof proceeds along the same lines as proof of Theorem 4.2 in [START_REF] Ghosh | An Introduction to Bayesian Analysis[END_REF], which establishes asymptotic normality of the posterior distribution of a strongly consistent solution of a likelihood equation (note that α(J) M P E coincides with the solution of the likelihood equation ∂L k (α) ∂α = 0 and is strongly consistent by Theorem 2.1). Some regularity conditions are needed. First, we note that the support of the density of (Z, δ) is the same for all α ∈ A. Moreover, L 1 (α) is thrice differentiable with respect to α and if U = (α 0 -δ, α 0 + δ) is a neighborhood of α 0 , we have

sup α∈U ∂ 3 L1(α) ∂α 3
< M (Z, δ), where EM (Z, δ) < ∞. Finally, for any δ > 0, there exists an > 0 such that for all sufficiently large k, we have

sup |α-α(J) M P E |>δ 1 k L k (α) -L k (α (J) M P E ) < -. (17)
To see this, consider the Taylor expansion

L k (α) = L k (α (J) M P E ) + (α - α(J) M P E ) ∂L k (α (J) M P E ) ∂α + (α - α(J) M P E ) 2 2 ∂ 2 L k (α) ∂α 2 ,
where α belongs to the line segment between α and α(J)

M P E . Then 1 k L k (α) -L k (α (J) M P E ) = - (α - α(J) M P E ) 2 2α 2 1 k k i=1 δ [n-k+i] .
Let ξ > 0 be arbitrary small. For k sufficiently large, 1 k k i=1 δ [n-k+i] < p + ξ and thus, for all α such that |α -α(J)

M P E | > δ, 1 k L k (α) -L k (α (J) M P E ) < - δ 2 2a 2 2 (p + ξ).
Finally, [START_REF] De Haan | Extreme Value Theory[END_REF] follows by letting := δ 2 2a 2

2

(p + ξ). We now establish asymptotic normality of

W k := √ k(α 0 - α(J) M P E ).
The posterior density of W k can be written as

π W k (w|E, ∆) := C -1 k π(α (J) M P E + k -1/2 w) exp L k (α (J) M P E + k -1/2 w) -L k (α (J) M P E ) , where C k = R π(α (J) M P E +k -1/2 w) exp[L k (α (J) M P E +k -1/2 w)-L k (α (J)
M P E )]dw. We show that π W k (w|E, ∆) converges in L 1 to the probability density function of the random variable N (0, 1 pI(α0) ), that is, we show that

I k := R π W k (w|E, ∆) - pI(α 0 ) 2π exp - w 2 2 pI(α 0 ) dw -→ 0 as k → ∞.
Convergence in L 1 -norm will imply that W k converges in distribution to N (0, 1 pI(α0) ). Let

g k (w) = π(α (J) M P E + k -1/2 w) exp L k (α (J) M P E + k -1/2 w) -L k ( α(J) M P E ) -π(α 0 ) exp - w 2 2 pI(α 0 ) .
Then we have:

I k ≤ C -1 k R |g k (w)| dw + R C -1 k π(α 0 ) - pI(α 0 ) 2π exp - w 2 2 pI(α 0 ) dw. (18)
We first prove that

R |g k (w)| dw → 0 as n → ∞ (and thus, as k → ∞). Let R = R 1 ∪ R 2 , where R 1 = {w : |w| > k 1/2 δ} and R 2 = {w : |w| < k 1/2 δ} for some δ > 0. We show that Rj |g k (w)| dw → 0 as n → ∞, for j = 1, 2. Note first that R1 |g k (w)| dw ≤ R1 π( α(J) M P E + k -1/2 w) exp L k (α (J) M P E + k -1/2 w) -L k (α (J) M P E ) dw + R1 π(α 0 ) exp - w 2 2 pI(α 0 ) dw ≤ exp(-k ) + π(α 0 ) R1 exp - w 2 2 pI(α 0 ) dw,
where the second inequality follows from [START_REF] De Haan | Extreme Value Theory[END_REF]. The term R1 exp -w 2 2 pI(α 0 ) dw is proportional to Ψ(-k 1/2 δ), where Ψ is the cdf of the normal distribution N (0, 1 pI(α0) ). Therefore, this term converges to 0 as n → ∞ and finally, R1 |g k (w)| dw → 0 as n → ∞. The proof that R2 |g k (w)| dw → 0 as n → ∞ is similar and is therefore omitted. Now, we consider the second term on the right-hand side of [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF]. Note first that

C k -π(α 0 ) 2π pI(α 0 ) ≤ R |g k (w)| dw.
This implies that C -1 k π(α 0 ) → pI(α0)

2π
as n → ∞ and thus, R |C -1 k π(α 0 )-pI(α0) 2π | exp[-w 2 2 pI(α 0 )] dw → 0 as n → ∞. Finally, I k → 0 as n → ∞ and thus, W k := √ k(α -α(J) M P E ) converges in distribution to N (0, 1 pI(α0) ).

We now establish asymptotic normality of α(J) M AP , α(MDI) M AP , α(c) M P E and α(c) M AP . Let α• be any of these estimators. We first show that a condition similar to [START_REF] De Haan | Extreme Value Theory[END_REF] converges a.s. to 0 as n → ∞. Morever, α(J) M P E maximizes L k (α), thus for k sufficiently large, we have:

1 k L k ( α(J) M P E ) -L k (α • ) < 2 ,
where is as in [START_REF] De Haan | Extreme Value Theory[END_REF]. Now,

1 k [L k (α) -L k ( α• )] = 1 k L k (α) -L k ( α(J) M P E ) + 1 k L k (α (J) M P E ) -L k (α • ) .
For k sufficiently large, |α -α(J) M P E | > δ and |α -α• | > δ hold simultaneously since α(J) M P E -α• converges to 0 a.s.. Thus, for k sufficiently large, we have:

sup |α-α•|>δ 1 k [L k (α) -L k (α • )] < -+ 2 = - 2 .
Weak convergence of √ k(α 0 -α• ) can now be proved along the same lines as for α(J) M P E . Details are omitted.

Appendix B. Results for the GPD model.

Results are provided for a sample size n = 1000 and m = 1000 simulated samples. For each = 1, 2, 3, 4, 5, we plot the empirical bias and RMSE of α as a function of k for p = 0.9, 0.75, 0.5, 0.25 (Figure 7: α = 0.5 and Figure 8: α = 0.75).

In Table 6, we provide the averaged value (over the m simulated samples), empirical bias and RMSE of each α at the optimal fraction level k opt = arg min k RM SE[α (k)] (for α = 0.5). Results for α = 0.75 yield similar conclusions and are therefore omitted.

k

  letting a u := Nu j=1 log e j,u and b u := Nu j=1 δ [n-Nu+j] , the log-posterior density (•|E, ∆) of α satisfies:

Figure 1 :

 1 Figure 1: Empirical bias (left) and RMSE (right) of α (l = 1, . . . , 5) for α = 0.5 and p = 0.25 (1st line:), p = 0.5 (2nd line), p = 0.75 (3rd line) and p = 0.9 (4th line).

Figure 2 :

 2 Figure 2: Empirical bias (left) and RMSE (right) of α (l = 1, . . . , 5) for α = 0.75 and p = 0.25 (1st line:), p = 0.5 (2nd line), p = 0.75 (3rd line) and p = 0.9 (4th line).

  α1

Figure 3 :

 3 Figure 3: Randomization of the fraction level: histograms of the Monte Carlo estimates α(j) , j = 1, . . . , N of the tail index α = 2 for p = 0.25 (1st line), p = 0.5 (2nd line), p = 0.75 (3rd line), p = 0.9 (4th line).

δ i 1

 1 {Zi>Z n,n-k } of non-censored observations in the k largest Z i 's as a function of k (see Figure (4)).

Figure 4 :

 4 Figure 4: Proportion of uncensored Z i 's among the k largest observations.

Figure 5 :

 5 Figure 5: Tail index estimation for Australian AIDS survival data.

Table 5 :

 5 Summary statistics for the Monte Carlo estimates α(j) , j = 1, . . . , N obtained with (a, b) = (4, 2) (see ( * )), (a, b) = (4, 2.5) (see ( * * )) and (a, b) = (4, 3) (see ( * * * )).

Figure 6 :

 6 Figure 6: Randomization of the fraction level: histograms of the α ( = 1, 2, 3) with (a * , b * ) = (4, 2) (1st line), (a * , b * ) = (4, 2.5) (2nd line), (a * , b * ) = (4, 3) (3rd line).

  holds for α• . Since α• and α(J) M P E are strongly consistent,1 k L k (α (J) M P E ) -L k (α • ) = log α(J) M P E -log α• 1 k k i=1 δ [n-k+i] -α(J) M P E -α• 1 k k i=1 log e [n-k+i]

Figure 7 :

 7 Figure 7: Empirical bias (left) and RMSE (right) of α (l = 1, . . . , 5) for α = 0.5 and p = 0.25 (1st line:), p = 0.5 (2nd line), p = 0.75 (3rd line) and p = 0.9 (4th line).

Figure 8 :

 8 Figure 8: Empirical bias (left) and RMSE (right) of α (l = 1, . . . , 5) for α = 0.75 and p = 0.25 (1st line:), p = 0.5 (2nd line), p = 0.75 (3rd line) and p = 0.9 (4th line).

Table 1 :

 1 Optimal results for Fréchet model with n = 1000, m = 1000 and α = 0.5.

			k opt Estimator	RMSE	Bias
		α1 133 0.4701402 0.04985720 -0.02985977
	p=0.9	α2 133 0.4663334 0.05199783 -0.03366664 α3 137 0.4804375 0.04398932 -0.01956245 α4 133 0.4703423 0.04950047 -0.02965768
		α5 133 0.4665644 0.05162741 -0.03343560
	p=0.75	α1 131 0.4610741 0.05884770 -0.03892594 α2 131 0.4565036 0.06172861 -0.04349636 α3 131 0.4752784 0.05071124 -0.02472156 α4 130 0.4616714 0.05834371 -0.03832862
		α5 130 0.4571051 0.06120118 -0.04289495
		α1 113 0.4428444 0.08150309 -0.05715557
	p=0.5	α2 113 0.4349585 0.08687673 -0.06504148 α3 118 0.4639250 0.06592663 -0.03607495 α4 111 0.4480398 0.08020518 -0.05196024
		α5 113 0.4359118 0.08560104 -0.06408818
	p=0.25	α1 α2 144 0.3704152 0.1468569 -0.12958479 89 0.3979008 0.1343862 -0.10209924 α3 86 0.4656990 0.0900378 -0.03430104 α4 89 0.4017657 0.1292944 -0.09823434
		α5 144 0.3734453 0.1434458 -0.12655470
			k opt Estimator	RMSE	Bias
		α1 137 0.7102997 0.07354345 -0.03970025
	p=0.9	α2 137 0.7048372 0.07626555 -0.04516285 α3 137 0.7233166 0.06780636 -0.02668343 α4 137 0.7079724 0.07410573 -0.04202759
		α5 120 0.7076827 0.07689530 -0.04231728
	p=0.75	α1 132 0.6930671 0.08929557 -0.05693289 α2 132 0.6863281 0.09336546 -0.06367190 α3 132 0.7093273 0.08017017 -0.04067273 α4 132 0.6904320 0.09015946 -0.05956802
		α5 132 0.6837832 0.09433285 -0.06621679
		α1	74	0.6909873 0.1176917 -0.05901272
	p=0.5	α2 105 0.6599687 0.1265239 -0.09003127 α3 115 0.6940075 0.1005289 -0.05599250 α4 74 0.6839079 0.1178655 -0.06609215
		α5 105 0.6558685 0.1277539 -0.09413154
	p=0.25	α1 α2 177 0.5438101 0.2246926 -0.20618988 89 0.5910295 0.2071481 -0.15897046 α3 89 0.6687678 0.1513115 -0.08123217 α4 89 0.5854156 0.2064865 -0.16458441
		α5 177 0.5424875 0.2248853 -0.20751253

Table 2 :

 2 Optimal results for Fréchet model with n = 1000, m = 1000 and α = 0.75.

Table 3 :

 3 Randomization of the fraction level: empirical mean, median and mode of the Monte Carlo estimates α(j) , j = 1, . . . , N of the tail index α = 2.

Table 4 :

 4 Standard deviations of the Monte Carlo estimates α(j)

	.353 0.381	0.410	0.443	0.456 0.463 0.469
	1.5 0.243 0.282	0.315	0.342	0.372 0.395 0.409
	2.0 0.194 0.202	0.233	0.263	0.286 0.312 0.341
	2.5 0.156 0.166	0.175	0.194	0.121 0.244 0.263
	3.0 0.134 0.141	0.144	0.152	0.167 0.187 0.205
	3.5 0.127 0.130	0.127	0.128	0.134 0.141 0.159
	4.0 0.122 0.123 0.120 0.118 0.117 0.124 0.127

1 , j = 1, . . . , N for various values of (a, b). The minimum standard deviations are indicated in bold.

  . Min. 1st Qu. Mode Median Mean 3rd Qu. Max.

		α1 0.120 1.112	1.272	1.293	1.390	1.381	1.470 2.602
	( * )	α2 0.119 1.111 α3 0.121 1.116	1.268 1.277	1.287 1.298	1.386 1.398	1.378 1.388	1.466 2.570 1.478 2.645
		α1 0.118 1.119	1.306	1.327	1.398	1.402	1.487	2.614
	( * * )	α2 0.117 1.118 α3 0.119 1.122	1.303 1.314	1.436 1.307	1.395 1.405	1.398 1.410	1.483 2.583 1.495 2.657
	( * * * )	α1 0.117 1.124 α2 0.117 1.123 α3 0.119 1.128	1.323 1.317 1.335	1.361 1.458 1.471	1.425 1.422 1.435	1.419 1.414 1.427	1.506 1.502 2.557 2.589 1.511 2.633

  2.1 and Theorem 2.2.

	Proof of Theorem 2.1.
	We prove strong consistency of	α(MDI) M AP . First, it follows from [9] (Theorem 2.3) and [11] that 1 k	k i=1 δ [n-i+1]
	converges almost surely (a.s.) to α α+β as n → ∞. From [7], 1 k	k i=1 log	Zn,n-i+1 Z n,n-k	converges a.s. to the
	EVI 1 α+β of Z, as n → ∞. By the continuous mapping theorem (see [28], for example),	α(MDI) M AP converges
	a.s. to α α+β • (α + β) = α as n → ∞. Proofs for	α(J) M P E ,	α(J) M AP ,	α(c) M P E and	α(c) M AP are similar and are thus
	omitted.
	Proof of Theorem 2.2.

Corresponding author : aameraoui@usthb.dz

Acknowledgements

The authors acknowledge financial support from the CNEPRU project "Modélisation -Simulation Actuarielle de Risques (MSAR)" (project number B00220120012, funded by Agence Nationale de Développement de la Recherche Universitaire, Algeria) and from the project SESA (project number EDC25895, funded by Centre National de la Recherche Scientifique -CNRS, France -and Direction de la Post-Graduation et de la Recherche-Formation -DPGRF, Algeria). Part of this work was made while 20 A. Ameraoui was visiting INSA of Rennes and while J.-F. Dupuy was visiting the Université des Sciences et de la Technologie Houari Boumediene, Algiers.

Bias p=0.9 α1 241 0.4871029 0.03316159 -0.012897068 α2 241 0.4848721 0.03397750 -0.015127864 α3 242 0.4937372 0.03125152 -0.006262834 α4 241 0.4871525 0.03301891 -0.012847470 α5 241 0.4849317 0.03383208 -0.015068322 p=0.75 α1 197 0.4813512 0.04283156 -0.01864877 α2 174 0.4817153 0.04407021 -0.01828467 α3 229 0.4820538 0.03934719 -0.01794616 α4 175 0.4853335 0.04255187 -0.01466651 α5 174 0.4818332 0.04375507 -0.01816682 p=0.5 α1 149 0.4650962 0.06667390 -0.03490379 α2 148 0.4593110 0.06983460 -0.04068903 α3 196 0.4655255 0.05698297 -0.03447449 α4 149 0.4654849 0.06588210 -0.03451512 α5 148 0.4597745 0.06900169 -0.04022548 p=0.25 α1 114 0.4453977 0.09310355 -0.05460229 α2 150 0.4228396 0.10027361 -0.07716037 α3 144 0.4660381 0.07223843 -0.03396190 α4 114 0.4467343 0.09076946 -0.05326572 α5 150 0.4242566 0.09843514 -0.07574341 Table 6: Optimal results for GPD model with n = 1000, m = 1000 and α = 0.5.


Oleksandr Stepanenko

Emmanuel Quillier

Hervé Tronche

Pascal Baldi

Marc De Micheli

Highly confining planar proton exchanged waveguides on z-cut LiNbO 3 without degrading its nonlinear coefficients: the HISoPE process

Highly confining planar proton exchanged waveguides on z-cut LiNbO 3 without degrading its nonlinear coefficients: the HISoPE process.

Oleksandr Stepanenko, Emmanuel Quillier, Hervé Tronche, Pascal Baldi, and Marc De Micheli

Abstract-We propose a new proton exchange process to reproducibly realize on lithium niobate highly confining (δne = 0.1) planar waveguides without degrading the nonlinear coefficient of the crystal. Waveguides were fabricated with different acidic baths and the crystallographic changes, the index profiles and the local SHG response were measured. We show that for certain fabrications conditions, characteristic of the so-called High Index Soft Proton Exchange process (HISoPE), the high index modification can be obtained without degrading the crystal nonlinear coefficient. This result should open the way to the fabrication of very efficient nonlinear devices.

Index Terms-waveguide, lithium niobate, integrated optics, proton exchange, nonlinear optics.

I. INTRODUCTION

L ITHIUM Niobate (LN) is known and widely used in integrated optics for its high nonlinear and electro-optic coefficients. Commercial devices are based on two fabrication techniques: the Ti-indiffusion [START_REF] Schmidt | Metal-diffused optical waveguides in LiNbO 3[END_REF], [START_REF] Neyer | High-speed cutoff modulator using a Tidiffused LiNbO 3 channel waveguide[END_REF], which allows obtaining very low losses for telecom modulators at 1.5 µm, and the Annealed Proton Exchange (APE) process [START_REF] Bortz | Annealed proton-exchanged LiNbO 3 waveguides[END_REF], which allows fabricating waveguides with a high polarization rejection and a better resistance to the optical damage, for applications using a shorter wavelength. In the labs, to improve the nonlinear efficiency, people have used Soft Proton Exchange (SPE) and Reverse Proton Exchange waveguides, which allow improving the confinement by respectively showing a slightly higher index increase or a symmetric index profile. All these techniques induce no or very small modifications of the crystallographic properties and are not affecting the electro-optics and nonlinear properties of the crystal but lead to gradient index profile waveguides with low index increase δn e (between 0.01 and 0.03).

Nowadays no LN device is using the most confining waveguides presenting a step index profile with the highest δn e =0.1 reported so far [START_REF] Baldi | Proton exchanged waveguides in LiNbO 3 and LiTaO 3 for integrated lasers and nonlinear frequency converters[END_REF]. Indeed such waveguides are obtained when the proton exchange is performed in pure acid (generally benzoic acid) heated between 200 and 300 • C and in this case the nonlinear coefficients are totally degraded [START_REF] Laurell | Loss of optical nonlinearity in proton-exchanged LiNbO 3 waveguides[END_REF]. Such waveguides are called PE waveguides in this paper.

Up to now, Vapor-phase Proton Exchange (VPE) is the only technique that proposed the possibility to combine a strong confinement (δn e =0.09) with preserved nonlinearities [START_REF] Rams | Preparation of proton-exchange LiNbO 3 waveguides in benzoic acid vapor[END_REF]. However, the control of the fabrication conditions of VPE The authors are with LPMC, CNRS UMR 7336, University Nice Sophia Antipolis, 06108 Nice, France (e-mail: oleksandr.stepanenko@unice.fr).

waveguides is rather complicated, which results in its very poor reproducibility [START_REF] Chanvillard | Interactions parametriques guidees de grande efficacite : utilisation de l'échange protonique doux sur niobate de lithium inverse periodiquement[END_REF] and explains the fact that no functional device realized by VPE has been reported yet.

In this paper we present a new waveguide fabrication technique on LN, which allows combining a high optical confinement with δn e ≥ 0.09, undegraded nonlinear coefficient and a good reproducibility. We call this process HISoPE for High Index Soft Proton Exchange.

II. FABRICATION OF THE WAVEGUIDES

In order to compare the properties of the HISoPE waveguides with those of the other members of the PE family, APE, SPE, etc., a set of planar waveguides was fabricated on z-cut LN samples, cut from 3" Gooch & Housego Palo Alto, (formerly Crystal Technology) wafers using a diamond wire saw. As proton sources, we used benzoic acid (BA) melts buffered with lithium benzoate (LB) with ρ LB = m(LB)/[m(LB) + m(BA)] ranging from 0 to 3%. In order to limit the number of varying parameters, all the samples were processed at 300 • C following the sealed ampules process [START_REF] Li | Réalisation sur LiNbO 3 de guides dondes présentant une forte variation d'indice et de trés faibles pertes[END_REF] using an inert metallic container that stands the benzoic acid vapor pressure much better than the Pyrex ampules.

The container is filled with the BA+LB mixture in its bottom part and the samples to process are fixed in the upper one, evacuated down to a pressure of 3 mbar, which allows controlling the atmosphere in the container, and finally hermetically closed. Then it is placed into a furnace at 300 • C and after 30 minutes, when both the melt and the samples have reached the exchange temperature, it is turned upside down, dipping the crystals into the melt and allowing the exchange to start. The exchange time can vary from several hours to several days depending on the aimed waveguide depth and the melt composition. When the exchange is finished the container is turned upside down again, allowed to cool and then open to get out the samples. In this work, we concentrated on samples produced with amounts of LB such that the working point is near but inferior to the threshold concentration ρ th = 2.7% above which we obtain SPE waveguides at 300 • C [START_REF] Baldi | Proton exchanged waveguides in LiNbO 3 and LiTaO 3 for integrated lasers and nonlinear frequency converters[END_REF].

The samples were fabricated using the following values of ρ LB : 0% to have a reference as a PE waveguide, 2.9% to have a reference as a SPE waveguide, and from 2% to 2.4% to have an indication of the range in which it is possible to produce HISoPE waveguides. Values below but near 2% have been shown to produce poor quality waveguides [START_REF] Baldi | Proton exchanged waveguides in LiNbO 3 and LiTaO 3 for integrated lasers and nonlinear frequency converters[END_REF] and were not used. Values between 2.4% and 2.7% were also not used as they are too close to ρ th % for SPE and may lead to waveguides which are also of poor quality.

III. CRYSTALLOGRAPHIC CHARACTERIZATION BY X-RAY DIFFRACTION

The crystalline structure of the planar waveguides has been studied using X-Ray diffraction method, which allows measuring the modifications of the lattice parameters of the exchanged layer. We used the diffraction from the (00.12) planes parallel to the surface of the z-cut LN samples, which, in the substrate, is observed at the Bragg angle θ = 40.85 • . This value is used as origin in the reported graphs (Fig. 1), and in order to identify the obtained crystallographic phases, the diffraction spectra are compared to the ones reported in [START_REF] Korkishko | Relationships between structural and optical properties of proton-exchanged waveguides on z-cut lithium niobate[END_REF].

Structural characterizations show a clear difference between the PE, HISoPE and SPE waveguides. SPE waveguides are associated to a graded variation of the cell parameter, which is represented by a gradient extending down to an angle of -150 indicating that all the waveguides is in the α phase (see Fig. 1).

-1000 -800 -600 -400 -200 0 200 PE waveguiding layer is formed of two crystallographic phases, β 1 represented by a broad peak around -580 and β 2 with a well defined peak at -840 . One has to note that the β 1 sublayer is not of the same crystalline quality as the β 2 sublayer. This may impact the optical properties. Note also that in this case, the lattice mismatch is very important and therefore important strains and stresses are present in this kind of waveguides, and produce unwanted effects [START_REF] Chen | Loss mechanisms and hybrid modes in high-δne proton-exchanged planar waveguides[END_REF] in the devices fabricated using this recipe.

HISoPE waveguides exhibit one sharp intense peak with the same FWHM as the substrate peak indicating that the crystalline quality has not been affected by the process (see Fig. 1). The position of this peak corresponds to the κ 2 phase. The substrate peak is no longer symmetric, and as for the SPE waveguides, a certain intensity can be seen down to an angle of -150 , which indicates that a small amount of the modified α phase is present in the waveguide. We can then expect a complex index profile with different sections associated to the different crystallographic layers. We can also note that the spectrum of the HISoPE waveguides are very similar to the X-ray rocking curves of the VPE waveguides [START_REF] Tsou | Structural and optical characterization for vapor-phase proton exchanged lithium niobate waveguides[END_REF] although the fabrication methods are completely different. At this stage of the characterization, similar nonlinear properties can be expected.

In this work we have let ρ LB vary from 2 to 2.4%. The Xray spectra (Fig. 2) show that the cell modification increases when ρ LB is reduced without any other dramatic change, which indicate that all the waveguides produced in this range of concentration are of the HISoPE type.

IV. INDEX PROFILE RECONSTRUCTION USING M-LINES

MEASUREMENTS

The effective indices of the modes of the different waveguides have been measured using a two-prisms coupling technique [START_REF] Tien | Theory of prism-film coupler and thin-film light guides[END_REF] in order to evaluate their propagation losses and to reconstruct their index profiles. Light scattering along the propagation allows estimating the propagation losses. A typical value is < 1 dB/cm.

The index reconstruction is a two steps process. The first step consists in using the IWKB method [START_REF] White | Optical waveguide refractive index profiles determined from measurement of mode indices: a simple analysis[END_REF], in order to evaluate the shape of the profile and the value of the depth and the index increase δn e of the different layers. This step cannot give a precise estimation when the index profile contains a step index part but as soon as they are thick enough to support a mode at the characterization wavelength, the different layers are clearly identified.

The second step consists in calculating the effective indices of a waveguide with the obtained complex index profile shape and looking for the best fit with the measured values. This step is used to optimize the profile given by the IWKB method, without adding new layers.

PE waveguides are formed of two different crystallographic layers (β 1 and β 2 phases) but the index profile is a simple step. SPE waveguides (α phase) present a graded variation of the cell parameter and a graded index profile (Fig. 4). For the HISoPE waveguides, we obtained a more complex index profile. The crystallographic characterizations indicate that there are two layers to take into account (Fig. 2), one showing a well defined cell parameters modification (κ 2 ) and an other one very similar to the α phase layer observed for SPE waveguides. In the index profile, both layers are clearly identified and we took them into consideration to show that the κ 2 layer is responsible for the step part of the index profile with a δn e = 0.1 at λ = 633 nm, while the graded layer is well fitted with an exponential function as for SPE waveguides (Fig. 3). It is worth noting that the values of the index increase involved in the profiles do not depend on ρ LB .

An important parameter for a waveguide fabrication technique being the reproducibility we test it by realizing samples with the same fabrication parameters and performed index profile characterization. The results are reported in Fig. 3, and show that the process is pretty well controlled.

0

V. SURFACE SHG RESPONSE

The local value of the d 33 nonlinear coefficient was probed by surface SHG microscopy using the experimental set-up described in our previous works [START_REF] Hadi | Quasi-phase-matched parametric interactions in proton-exchanged lithium niobate waveguides[END_REF]. This setup, which has been improved and is now known as scanning SHG microscopy, is widely used to visualize domains in ferroelectric crystals [START_REF] Denev | Probing ferroelectrics using optical second harmonic generation[END_REF]. A laser beam at 1.55 µm is focused on the polished end facets of the waveguides and scanned through the different layers: substrate, exchanged layers, air. This characterization technique permits to compare the value of the nonlinear coefficient in the waveguiding layers with its value in the substrate. The results are given in Figs. 45. The intensities of the fundamental and the second harmonic are normalized by their maxima and plotted as a function of the waveguide depth. The index profiles are also plotted to see the correlations between the modifications of the index and of the nonlinear coefficient.

For the PE waveguide we obtain a SHG response that is strongly reduced in the exchanged layer (Fig. 4). In other work, people have reported no signal in the exchanged layer, but it is worth to note that the waveguide tested here is fabricated at 300 • C in a hermetically closed container while the other samples were fabricated at lower temperature in open containers. As we know from previous works on x-cut LN substrates, increasing the exchange temperature improves the crystalline quality of the exchanged layer [START_REF] Chen | Loss mechanisms and hybrid modes in high-δne proton-exchanged planar waveguides[END_REF], [START_REF] Micheli | Crystalline and optical quality of proton exchanged waveguides[END_REF]. Some differences in the nonlinear response may be due to this difference in the fabrication conditions despite the index profiles are similar.

For the SPE waveguide (Fig. 4), we see that the SH signal follows exactly the fundamental signal and as soon as the beam is completely focused in the material the SHG signal is constant showing no variation between the exchanged area and the substrate. This indicates that the α layer presents the same nonlinear coefficient as the substrate, which has been verified in many other experiments using these nonlinear waveguides [START_REF] Korkishko | The SHG-response of different phases in proton exchanged lithium niobate waveguides[END_REF]. It is worth to note that in these surface SHG experiments, one always observe a maximum of the SHG signal when the spot crosses an interface with a discontinuity of the d 33 along the propagation direction [START_REF] Kaneshiro | Threedimensional observations of polar domain structures using a confocal second-harmonic generation interference microscope[END_REF]. This is related to the phase profile in the focused beam and the sign of the phase mismatch in the material. It is therefore very important to verify that during the scan the focus point is always at the same distance of the polished end face. This is proved by the fact that the SH signal is constant when the scan explore the substrate.

With the HISoPE waveguides (Fig. 5), the SHG response shows a more or less sharp peak after the air-crystal interface but no variation at the waveguide-substrate interface. The value of this maximum is probably affected by the index modification in the crystal and the associated guiding effect, despite in the experiments nothing is done to optimize the coupling into the waveguide. Therefore, without a careful modeling of the reflected SHG signal it is impossible to extract information about the exact value of the nonlinear coefficient near the surface but it is possible to say that d 33 is not reduced as in the PE waveguides.

In [START_REF] Hadi | Quasi-phase-matched parametric interactions in proton-exchanged lithium niobate waveguides[END_REF], we also observed such a SHG peak with APE waveguides. We explained this peak by an extra backscattering coming from the bad quality of the produced waveguides. In this case the phenomena was also associated with a strong degradation of the retro-reflected fundamental beam, which was not observed with any of the samples we tested here. At the waveguide substrate interface, the absence of variation of the SH signal, allows saying that there is no variation of the d 33 coefficient on the same area. Therefore, despite the presence of the SHG peak is difficult to explain, these results show that HISoPE process allows fabricating low loss planar waveguides in z-cut LN exhibiting a δn e = 0.1 and no reduction of the nonlinear coefficient. This result is very important as it opens the way to the fabrication of very efficient nonlinear devices.

VI. CONCLUSION

In this paper we have presented a preliminary work aiming at fabricating highly confining waveguides on LN without destroying the nonlinear properties of the crystal. By measuring the crystallographic properties of the produced waveguides as well as their linear and nonlinear properties, we have shown that adding an appropriate amount of lithium benzoate in benzoic acid allows obtaining this result. This study confirms that very different layers can exhibit the same index profile, which can explain some "contradictions" reported in the literature, and show that there is a strong correlation between crystallographic and nonlinear properties. Finally, we have demonstrated that the HISoPE process, which produces highly confining planar waveguides with preserved nonlinearities, is quite reproducible what makes of him a good candidate to realize highly efficient nonlinear devices.

Fig. 1 .

 1 Fig.1. Typical X-ray rocking curves of SPE, PE and HISoPE planar waveguides. The temperature of fabrication is the same for the three samples and is 300 • C. Lithium benzoate concentration was ρ LB = 2.9%, 2.4% and 0% for SPE, HISoPE and PE waveguides, respectively.

Fig. 3 .

 3 Fig. 3. Index profiles (at λ = 633 nm) of 4 planar HISoPE waveguides realized at the same conditions: 300 • C in ρ LB =2% melt during 17 hours. All the sample were fabricated at different time during the year.

Fig. 4 .

 4 Fig. 4. Surface SHG profiling of PE (a) and SPE (b) waveguide and the corresponding index profiles (at λ = 633 nm). The samples were fabricated with ρ LB = 0% at 360 • C during 30 minutes for PE, and with ρ LB = 2.9% at 300 • C during 3 days for SPE.

Fig. 5 .

 5 Fig. 5. Surface SHG profiling of HISoPE waveguides and their index profiles at λ = 633 nm. The samples were fabricated using ρ LB : (a) 2.4% during 35 hours; (b) 2.2% during 72 hours; (c) 2% during 40 hours. Fabrication temperature is 300 • C.

ACKNOWLEDGMENT

The authors would like to thank Maud Nemoz (CNRS-CRHEA, rue Bernard Grégory, 06560 VALBONNE, France) for her valuable contribution to this research providing a useful help in XRD characterizations.


Parametric downconversion in submicronic periodically poled (PP) quadratic materials allows the generation of signal and idler waves in opposite directions. The distributed feedback mechanism enables a backward mirrorless optical parametric oscillator (BMOPO), which has been experimentally realized using a periodically-poled KTP crystal with 800 nm periodicity. A remarkable spectral property of the BMOPO is that the bandwidth of the forward co-propagating wave is comparable to the pump bandwidth, whereas that of the backward wave is typically several orders of magnitude narrower. In the backward idler configuration realized in PPKTP the backward idler exhibits a coherence gain of two orders of magnitude. In this paper we will consider the backward signal configuration, where the effect of coherence enhancement of the backward signal wave is especially pronounced for exact group velocity matching of the co-propagating pump and idler waves. This QPM scheme requires a 335 nm scale poling periodicity which may be achieved in GaN waveguides, combining e-beam lithography and epitaxy. However, during the lithography step necessary to obtain the PP material, stitching faults may occur. The BMOPO device is then formed by a sequence of 200 microns long PP elements stuck together and we assume that the stitching faults do not exceed half of the poling period, in order to prevent the parametric reverse process. The original three-wave QPM model with the explicit periodicity of the nonlinear parametric coupling is necessary for taking into account the periodicity faults in the junctions. Thanks to the High Performance Computing facility available at Nice University, we analyse the parametric downconversion efficiency and the coherence gain of the backward signal versus the stochastic pump bandwidth in a large pump bandwidth range. We show that if the stitching faults do not exceed a quarter of the poling period, i.e. less than 80 nm, which is technically feasible, the efficiency is comparable to the BMOPO in the perfect PP material. In any case, for an incoherent pump of broad spectral bandwidth we obtain a coherence gain of almost three orders of magnitude. c 2014 Optical Society of America OCIS codes: 030. 1640, 190.4410, 190.4970, 190.5970, 230.4320

Introduction

Optical parametric oscillators based on periodically poled (PP) nonlinear quadratic materials, with submicronic poling period allow the generation of quasi-phase-matched (QPM) signal and idler waves in opposite directions [START_REF] Harris | Proposed backward wave oscillation in the infrared[END_REF][START_REF] Ding | Backward optical parametric oscillators and amplifiers[END_REF]. The distributed feedback mechanism enables backward mirrorless optical parametric oscillation (BMOPO). This has been experimentally realized in the backward idler configuration using a periodically-poled KTiOPO 4 (PPKTP) crystal with Λ G = 800 nm periodicity as a QPM nonlinear medium [START_REF] Canalias | Mirrorless optical parametric oscillator[END_REF]. A remarkable spectral property of the BMOPO is that the co-propagating forward wave has a bandwidth comparable to the pump bandwidth, whereas that of the backward wave is typically several orders of magnitude narrower. Two orders of magnitude have been obtained in this backward idler configuration, recently demonstrated both experimentally and numerically for a pump wave whose bandwidth is stretched by deterministic phase modulation in a PPKTP BMOPO for relatively narrowband (1.2 THz) linearly-chirped pump pulses [START_REF] Strömqvist | Coherent phasemodulation transfer in counterpropagating parametric down-conversion[END_REF] and relatively broadband (4 THz) pump pulses too [START_REF] Strömqvist | Temporal coherence in mirrorless optical parametric oscillators[END_REF].

In this paper we will consider the backward signal configuration, where the effect of coherence enhancement of the backward signal wave is especially pronounced owing to the possibility of exact group velocity matching of the forward co-propagating pump and idler waves. But the short nanometric periods required for QPM have been only experimentally achieved, up to now in GaN, combining e-beam lithography and epitaxy [START_REF] Chowdhury | Second-harmonic generation in periodically poled GaN[END_REF]- [START_REF] Pezzagna | Refractive indices and elasto-optic coefficients of GaN studied by optical waveguiding[END_REF]. Asymmetric GaN waveguides grown on (0001) sapphire by molecular beam epitaxy or by metal-organic chemical vapor deposition are promising structures for applications in nonlinear optics. Periodically poled GaN (PPGaN) was used for second harmonic generation [START_REF] Chowdhury | Second-harmonic generation in periodically poled GaN[END_REF] and the possibility to grow submicronic PP GaN films [START_REF] Pezzagna | Submicron periodic poling and chemical patterning of GaN[END_REF] opens the door for parametric counterpropagating oscillation thanks to reasonably large second order susceptibility (χ (2) 33 = 10.6 pm/V measured in [START_REF] Pezzagna | Refractive indices and elasto-optic coefficients of GaN studied by optical waveguiding[END_REF]), comparable to that of KTP.

However, during the fabrication process of the PPGaN waveguide, stitching errors cannot be avoided for centimeter long interaction lengths. Indeed, during the lithography step necessary to obtain the PP material, the realization of centimeter long devices presenting a submicronic periodicity is not possible without moving the sample under the writing beam. This displacement introduces the stitching faults. We have studied the dynamics in a PP-GaN BMOPO formed by a sequence of 200 microns long PP elements (PPGaN bricks) stuck together, so that the stitching faults do not exceed half of the poling period Λ G in order to avoid the switch from parametric downconversion to the reverse process of sum frequency generation, which combines signal and idler to create pump [START_REF] Gallo | All-optical diode based on second-harmonic generation in an asymmetric waveguide[END_REF][START_REF] Subramani | Analysis of a phase reversal quasi-phase matching device for the dual peak second harmonic response[END_REF]. The defects will also affect the parametric downconversion efficiency and the phenomenon of coherence enhancement which decrease as the amplitude of the perturbations increases. In order to take into account the periodicity perturbations leading to phase shifts in the junctions between the perfect QPM PP segments, we must consider the original QPM model [START_REF] Armstrong | Interaction between light waves in a nonlinear dielectric[END_REF][START_REF] Myers | Quasi-phase-matched optical parametric oscillators in bulk periodically poled LiNbO3[END_REF] with the complete explicit dependece of the nonlinear coupling parameter, without substituting the grating periodicity by an effective quadratic nonlinearity, i.e. without introducing the usual exact phase-matching approximation [START_REF] Montes | Backward optical parametric efficiency in quasi-phase-matched GaN waveguide presenting stitching faults[END_REF].

Our aim here is to study the coherence transfer for a broad bandwidth stochastic pump.

The generation of a coherent wave in the process of parametric down-conversion of an incoherent wave by the convection-induced phase-locking mechanism has been the object of theoretical and numerical studies [START_REF] Picozzi | Parametric Three-Wave Soliton Generated from Incoherent Light[END_REF]- [START_REF] Montes | Coherent signal from incoherently cw-pumped singly resonant Ti:LiNbO3 integrated optical parametric oscillators[END_REF]. In accordance with this mechanism, the phase modulation in the pump is coherently transferred to the forward parametric wave, while the backward wave retains a narrow bandwidth and high coherence. Some recent studies show the parametric amplification of the signal wave by conversion of a broadband incoherent pump to the group velocity matched idler wave in the forward three-wave interaction [START_REF] Pyragaite | Parametric amplification in the field of incoherent ligth[END_REF]- [START_REF] Piskarskas | Generation of coherent waves by frequency up-conversion and down-conversion of incoherent light3[END_REF] and using the generalized nonlinear envelope equation [START_REF] Wabnitz | Control of signal coherence in parametric frequency mixing with incoherent pumps: narrowband mid-infrared light generation by downconversion of broadband amplified spontaneous emission source at 1550 nm[END_REF]. Thanks to the 1200 multiprocessor computer facility (HPC) of the "Centre de Calcul Interactif" hosted by the "Université Nice Sophia Antipolis", the coherence transfer mechanism is studied here for broad incoherent pump pulses in a large spectral bandwidth range [0.08 GHz -6.0 THz]. We analyse the downconversion efficiency versus the stochastic pump bandwidth for fragmented PPGaN BMOPOs whose stitched junctions have a perturbation length lower than Λ G /2, and we compare them to the perfect PPGaN BMOPO. For a stochastic pump of broad spectral bandwidth and maximum perturbation length of 80 nm, which is technically feasible, we obtain a coherence gain of almost three orders of magnitude.

BMOPO model for a QPM fragmented waveguide

In a BMOPO, the pump wave is downconverted into forward-propagating and backwardpropagating parametric waves, in accordance with energy conservation, ω p = ω s + ω i and counterpropagating quasi-phase matching, k p + k sk i = K G , where j = p, s, i denote the pump, the backward signal and the forward idler waves, respectively. K G = 2π/Λ G is the grating wave vector given by Λ G , the modulation period of the second-order nonlinear coefficient. The coherence enhancement effect of the backward wave is especially pronounced when the pump wavelength and the index profile of the waveguide are chosen so that the group velocities of the pump and the forward idler wave are exactly matched [START_REF] Picozzi | Coherence properties of the parametric threewave interaction driven from an incoherent pump[END_REF]. This configuration may be achieved in a GaN waveguide with a group velocity dispersion modified by a particular profile of the refractive index [START_REF] Picozzi | Influence of dispersion on the resonant interaction between three incoherent waves[END_REF], a technique widely used in fibers to compensate the material dispersion with the waveguide dispersion [START_REF] Monerie | Propagation in doubly clad single-mode fibers[END_REF]. For z-polarized waves in PPGaN, matching group velocities can be achieved by designing the index profile so that the pump and the forward wave are on different sides of the maximum on the group-velocity curve shown in Fig. 1. This shape of group-velocity curve may be obtained in the guided PPGaN structure composed of a ridge GaN waveguide of 2 μm in the width and 1.8 μm in the depth on a AlInN buffer which has been realized on sapphire substrate by epitaxy. The Sellmeier relation given in [START_REF] Pezzagna | Refractive indices and elasto-optic coefficients of GaN studied by optical waveguiding[END_REF] is used for material dispersion, and finite waveguide dispersion corrections [START_REF] Picozzi | Influence of dispersion on the resonant interaction between three incoherent waves[END_REF] yield the group-velocity curve of Fig. 1 This may be realized in the GaN structure combining e-beam lithography and epitaxy [START_REF] Pezzagna | Submicron periodic poling and chemical patterning of GaN[END_REF][START_REF] Pezzagna | Refractive indices and elasto-optic coefficients of GaN studied by optical waveguiding[END_REF]. To fulfil QPM condition with the set of wavelengths chosen here we need a periodicity of the grating of Λ G = 335 nm. Fig. 1. Group velocity for z-polarized waves in GaN as function of the wavelength, calculated from the Sellmeier expansion [START_REF] Pezzagna | Refractive indices and elasto-optic coefficients of GaN studied by optical waveguiding[END_REF] with a particular profile of the refractive index of the guide [START_REF] Picozzi | Influence of dispersion on the resonant interaction between three incoherent waves[END_REF]. The triangles mark the configuration of perfectely matched pump and idler group velocities (|v pv i | = 0), and are labeled P for the pump, S for the backward signal, and I for the idler.

The pump intensity required to reach threshold for BMOPO operation, I p,th , can be estimated from the plane-wave monochromatic model for the undepleted pump parametric approximation [START_REF] Ding | Backward optical parametric oscillators and amplifiers[END_REF] as

I p,th = cε 0 n p n f n b λ f λ b /2L 2 d 2 ef f
, where c is the speed of light in vacuum, ε 0 is the permittivity of free space, n j denotes the refractive indices, λ j denotes the backward and forward wavelengths, L is the nonlinear interaction length and d ef f is the effective quadratic nonlinearity.

Numerical simulations of the nonlinear three-wave interactions in a BMOPO solve the coupled three wave equations in the slowly-varying-envelope approximation. The parameters used in the model are those of a periodically poled GaN waveguide (PPGaN) grown by epitaxy on sapphire substrates [START_REF] Chowdhury | Second-harmonic generation in periodically poled GaN[END_REF] and presenting submicronic periodicity of the χ (2) coefficient [START_REF] Pezzagna | Submicron periodic poling and chemical patterning of GaN[END_REF]. The effective nonlinear coefficient is then given by d ef f = (2/π)χ 33 , d ef f varies from 7 pm/V [START_REF] Pezzagna | Refractive indices and elasto-optic coefficients of GaN studied by optical waveguiding[END_REF] to 3.8 pm/V [START_REF] Abe | Accurate measurement of quadratic nonlinear-optical coefficients of gallium nitrid[END_REF] (±45% of accuracy). For the mechanism studied here we will take the value of d ef f = 5 pm/V. In order to account for stitching faults in the joints, or periodicity perturbations, simulated through phase shifts between the perfect QPM segments, we call on the original QPM model [START_REF] Armstrong | Interaction between light waves in a nonlinear dielectric[END_REF][START_REF] Myers | Quasi-phase-matched optical parametric oscillators in bulk periodically poled LiNbO3[END_REF], without substituting the grating periodicity of the nonlinear coupling parameter by an effective quadratic nonlinearity. Indeed, we have shown -within the exact phase-matching model -, that the parametric downconversion efficiency and the coherence enhancement is as high as for the perfect single piece PP waveguide whatever be the number of segments separated by uniformly polarized sections [START_REF] Montes | Backward optical parametric efficiency in quasi-phase-matched GaN waveguide presenting stitching faults[END_REF]. This would correspond to an ideal case where the separation of the perfect PP segments has the length of an entire number of grating periods.

The field amplitudes, A j (j = p, s, i), of the pump and the parametric backward signal wave and forward idler wave, evolve with the coupled equations:

(∂ t + v p ∂ x + γ p + iβ p ∂ tt ) A p = -σ p cos[K G x + δφ m (x)]A s A i exp(iΔKx) (∂ t -v s ∂ x + γ s + iβ s ∂ tt ) A s = σ s cos[K G x + δφ m (x)]A p A * i exp(-iΔKx) (1)
(∂ t + v i ∂ x + γ i + iβ i ∂ tt ) A i = σ i cos[K G x + δφ m (x)]A p A * s exp(-iΔKx),
where γ j and β j ≡ v j β 2,j /2 are the damping and dispersion coefficients, and σ j = 2πd ef f v j /λ j n j are the nonlinear coupling coefficients in the PP elements, with the spatial period Λ G = 2π/K G of the periodic poled χ

(2)

33 domains and ΔK = k p + k sk i . By expanding the second-order susceptibility in Fourier series and retaining the first order we obtain χ [START_REF] Ding | Backward optical parametric oscillators and amplifiers[END_REF]. Inside the PP elements σ j cos(K G x) is periodically perfect, so that QPM is not affected by fluctuating nonlinearity [START_REF] Kh | Evolution of spatial solitons in quadratic media with fluctuating nonlinearity[END_REF] or random effects [START_REF] Helmfrid | Influence of randomly varying domain lengths and nonuniform effective index on second-harmonic generation in quasi-phase-matching waveguides[END_REF][START_REF] Conti | Parametric self-trapping in the presence of randomized quasi phase matching[END_REF], which would lead to effective nonlinear losses. The stitching faults introduced through finite phase shift perturbations δφ m (x) (m=1,N-1) between the N stuck PP segments (cf. Fig. 2) are only located at each junction m. We take them to be randomly distributed in the interval (0, pπ)), with p < 1 to avoid the parametric reverse process. In order to preserve efficient paramet-ric amplification and coherence enhancement we will limit ourselves to phase perturbations δφ m (x) ≤ π/2, which correspond to perturbation lengths ≤ Λ G /4. By expanding the secondorder susceptibility of the domains in Fourier series and retaining the first order grating we obtain χ [START_REF] Ding | Backward optical parametric oscillators and amplifiers[END_REF]. The exact phase-matched model equations [START_REF] Montes | Backward optical parametric efficiency in quasi-phase-matched GaN waveguide presenting stitching faults[END_REF] are obtained from equations (1) by imposing the matching condition exp(iK G x) = exp(iΔKx) but differ by non-phase-matched terms. As known [START_REF] Fejer | Quasi-phase-matched second harmonic generation: tuning and tolerances[END_REF], the parametric growth rate governed by ( 1) is lower than that for the exact QPM model.

(2) F χ (2) ef f cos(K G x)
(2) F χ (2) ef f cos(K G x)
The input parameters in the model are the properties of the nonlinear medium and the pump amplitude at the input face, A p (x = 0, t), generating outputs of A p (x = L, t), A i (x = L, t) and A s (x = 0, t). The nonlinear counterpropagation dynamics in the quadratic BMOPO in the presence of group-velocity dispersion, which introduces second-order time derivatives, makes use of the numerical scheme which combines the trajectories method and 4th-order Runge Kutta algorithm for the treatment of the nonlinear counterpropagating three-wave interaction, extensively used for the treatment of stimulated Brillouin backscattering problems [START_REF] Montes | Dissipative three-wave structures in stimulated backscattering. I. A subluminous solitary attractor[END_REF][START_REF] Montes | Dissipative three-wave structures in stimulated backscattering. II. Superluminous and subluminous solitons[END_REF], with fast Fourier transformation (FFT) to account for the groupvelocity dispersion effects in the spectral domain [START_REF] Montes | Model for Coherence Transfer in a Backward Optical Parametric Oscillator[END_REF].

The nonlinear medium is a GaN structured ensemble composed by a train of N periodicallypoled GaN elements of length disposed in a straight row as shown in Fig. 2. Simulations for this segmented structure of length L= N , where a random phase-shift δφ m (x) (m=1,N-1) takes place at each jonction m, representing the stitching fault between the N stuck PP segments, associated to a physical stitching length δL m (x) = Λ G (δφ m (x)/(2π), will be compared to the case of a one piece perfect PP medium of length L without phase-shifts. For L = 64 = 1.25 cm, each PP element has a length = 195 μm. The waveguide has a section S eff = 3.6 μm 2 .

We choose the three QPM waves marked on the Fig. 1, and we model BMOPO operation with a broadband stochastic pump pulse with a FWHM temporal length of 150 ps, and where the spectral FWHM width Δν p (0) centered at λ p = 1.06μm varies in a large spectral bandwidth range [0.08 GHz -6.0 THz], with n p = 2.30856, v p /c = 0.41337, β 2,p = 0.286 ps 2 /m; and with λ s = 1.56242μm, n s = 2.26580, v s /c = 0.41872, β 2,s = 0.0647 ps 2 /m; λ i = 3.29791μm, n i = 2.12353, v i /c = 0.41337, β 2,i = -0.418 ps 2 /m. The QPM period as low as Λ G = 335 nm may be realized in GaN by combining e-beam lithography and epitaxy [START_REF] Pezzagna | Refractive indices and elasto-optic coefficients of GaN studied by optical waveguiding[END_REF].

The pump pulse is issued from a Gaussian stochastic process translationally invariant with zero mean A p (x = 0, t) and exponential autocorrelation function

A p (x = 0, t + t)A * p (x = 0, t) = A 0 p 2 exp (-|t|/τ c )
, where τ c = 1/πΔν p is the correlation time. Note that, while the pump phase fluctuations are absorbed by the idler wave thanks to the phase-locking mechanism, the generation of the signal still remains affected by the intensity fluctuations of the pump. In the presence of small group velocity difference between the pump and the signal, which is the case for the forward OPO, pump intensity fluctuations are partly transferred to the signal component [START_REF] Montes | Ultra-coherent output from an incoherent cw-pumped singly resonant optical parametric oscillator[END_REF][START_REF] Montes | Coherent signal from incoherently cw-pumped singly resonant Ti:LiNbO3 integrated optical parametric oscillators[END_REF]. Conversely, in the presence of the strong group velocity difference of the backward MOPO, intensity fluctuations of the pump are averaged by the strong convection effect, which thus leads to an enhancement of the coherence of the backward wave [START_REF] Strömqvist | Temporal coherence in mirrorless optical parametric oscillators[END_REF] [START_REF] Montes | Backward optical parametric efficiency in quasi-phase-matched GaN waveguide presenting stitching faults[END_REF]. The higher the incoherence of the pump, the larger the amount of coherence transfer as studied in section V.B of Ref. [START_REF] Picozzi | Coherence properties of the parametric threewave interaction driven from an incoherent pump[END_REF].

Stochastic BMOPO operation in both continuous and segmented structures are simulated with a pump pulse of intensity I p = 4 GW/cm 2 , the BMOPO threshold intensity being I p,th = 2 GW/cm 2 . The peak power yields 140 W, a pulsed energy of 20 nJ, and a fluence of 0.6 J/cm 2 which remains below the damage threshold for GaN [START_REF] Wang | Femtosecond pulsed laserinduced periodic surface structures on GaN/sapphire[END_REF][START_REF] Ozono | High-speed ablation etching of GaN semiconductor using femtosecond laser[END_REF].

The resonant three-wave interaction inside the PPGaN BMOPO composed of the 64 PP segments is illustrated in Fig 3 . We show in Fig 3(a), at time t = 254 ns near the maximum of the nonlinear interaction, the forward stochastic pump amplitude, the induced stochastic forward idler amplitude and the resulting high coherent backward amplitude. The phaselocking behaviour inside the entire BMOPO length L is illustrated in Fig 3(b) and zoomed in the the first L/10 length in Fig 3(c). We remark that the forward idler phase is a replica (phase conjugate) of the stochastic pump phase, whereas the backward signal phase remains almost constant throughout the medium. As expected from the convection-induced phaselocking mechanism, even for the fragmented BMOPO, the stochastic phase modulation in the pump is essentially transferred to the forward idler, while the phase of the backward signal is almost constant.

To avoid other nonlinearities, such as stimulated Raman scattering, which is the most harmful nonlinearity competing with the parametric downconversion [START_REF] Canalias | Mirrorless optical parametric oscillator[END_REF], we can further decrease the pump intensity and increase the interaction length. With I p = 2 GW/cm 2 and a fluence of 0.3 J/cm 2 , and for a pump depletion of 15% we attain a coherence gain Δν p (L)/Δν s (0) = 173.

Efficiency and coherence gain within the quasi-phase-matching model

In order to account for stitching faults which are able to arise in the joints, simulated through phase shifts or finite perturbations between the perfect QPM segments we use the original QPM model of Eqs.(1), without substituting the grating periodicity of the nonlinear coupling parameter by an effective quadratic nonlinearity. Within this QPM model, the pump depletion or parametric downconversion efficiency versus initial pump bandwidth for the 64 times segmented BMOPO is compared to the one piece perfect PP waveguide shown in ,N-1) and the direct comparison to the upper curve for the perfect PP waveguide yields a relative efficiency of 63%. We obtain that the downconversion efficiency and the coherence enhancement are limited by the stitching defects as much as their length increase approaching Λ G /2. For Λ G /2 the downconversion efficiency vanishes. It is interesting to recall that for the ideal case where the stitching length is an entire number of Λ G , the exact phase-matching model provides a parametric downconversion efficiency and a coherence signal enhancement as high as for the single piece perfect PP waveguide, provided that the sum of the PP segment lengths be equal to the length of the uniform PP waveguide [START_REF] Montes | Backward optical parametric efficiency in quasi-phase-matched GaN waveguide presenting stitching faults[END_REF]. Each point of the curves plotted in The apparent random variation of the pump depletion or downconversion efficiency around the mean linear fitting is due to stochastic resonances in the nonlinear medium associated to the initial phase. Its amplitude lowers when the mean value is performed over more initial phases.

The pump and the parametric spectra for the 64 times segmented BMOPO are illustrated in Fig 7 to Fig 9 for increasing pump bandwidth input and increasing phase-shift amplitude perturbation. As we can see, the higher the incoherence of the pump (larger spectral bandwidth), the larger the amount of coherence enhancement of the signal wave. By integrating the spectra, it is found that the conversion into parametric waves for the five chosen cases shown in Table 1, decrease with increasing pump bandwidth. For an incoherent pump spectrum of FWHM Δν p (L) = 4.54 THz, phase-shift amplitude perturbations in the joints ≤ π/2, and a pump depletion of 0.0256, Figure 9 shows a coherence gain Δν p /Δν s = 467, which is comparable to the perfect PPGaN BMOPO. This maximum stitching of π/2, i.e. 80 nm length, is technically attainable. Table 1. Parametric downconversion efficiency and coherence gain in the 64 times segmented GaN BMOPO with stitching faults of maximum π/4, of maximum π/3 (marked with * ), and of maximum π/2 (marked with * * ). The coherence gain attains Δν p /Δν s = 347, for a pump depletion or efficiency of 0.0759 (cf. Table 1). The coherence gain attains Δν p /Δν s = 467, for a pump depletion or efficiency of 0.0256 (cf. Table 1).

Δν p (0) Δν p (L) Δν s (0) Δν i (L) 1 -Ip(L) Ip(0) Is(0) Ip(0) I i (L) Ip(0) Δνp(L) Δνs(0) [THz] [THz] [GHz] [

Summary

We have numerically simulated the dynamics of downconversion efficiency and coherence signal enhancement in fragmented BMOPOs formed of a realistic PPGaN waveguide by using stochastic pump pulses from 0.08 GHz to 6.0 THz input bandwidth. We have performed the study with the original quasi-phase-matching model of Eqs. [START_REF] Harris | Proposed backward wave oscillation in the infrared[END_REF], which allow to consider finite phase perturbations in the junctions between the PPGaN segments in order to account for possible stitching faults. This realistic model assess the limits of the phenomenon of coherence enhancement when the phase perturbations approach half of the poling period because the parametric downconversion is replaced by sum frequency, its reverse process. For stitching errors less than the quarter of the grating period (i.e. less than 80 nm), which is technically acceptable, the efficiency remains comparable to the perfect PP BMOPO, and we may obtain a coherence gain of almost three orders of magnitude.

Fig. 2 .

 2 Fig. 2. The quasi-phase-matching (QPM) is obtained by a phase-reversal sequence of the ferroelectric domains with a grating periodicity Λ G = Λ QP M = 335 nm. The segmented backward MOPO is composed of a train of periodically poled GaN elements of length = 195 μm. Stitching faults between the N PP elements stuck together, are simulated by a phase-shift perturbation δφ m (x) (m=1,N-1) at each joint m, whose amplitude is randomly distributed in the interval (0,pπ), with p < 1 to avoid the parametric reverse process. It is associated to a stitching length Λ G δφ m (x)/(2π).

  Fig 4(a), and repeated in the upper curves of Fig 5 and Fig 6 for direct comparison. In Fig 4(b), it is shown for random phase perturbations in the junctions 0 < δφ m (x) < π/4 (m=1,N-1) and presents a relative efficiency of 93% with respect to the perfect PP waveguide shown in Fig 4(a). In the lower curve of Fig 5 it is shown for random phase perturbations in the junctions 0 < δφ m (x) < π/3 (m=1,N-1) and the direct comparison to the upper curve for the perfect PP waveguide yields a relative efficiency of 87%. In the lower curve of Fig 6 it is shown for random phase perturbations in the junctions 0 < δφ m (x) < π/2 (m=1

  Fig 4(a)(b), Fig 5, and Fig 6, is a mean value over three different initial phases and requires 21 hours CPU numerical calculation due to the combined Runge-Kutta-FFT procedure [30]. Each curve contains 284 points and the total CPU time for one processor would be 21 × 284 = 5964 CPU hours for each curve. The use of High Performance Computing facility (a 1200 multiprocessor computer) has reduced this time by more than two orders of magnitude.

Fig. 3 .

 3 Fig. 3. (a) Amplitude distributions inside the 64 times segmented GaN BMOPO for the forward stochastic pump, the forward induced stochastic idler and the resulting coherent backward signal; (b) phase-locking inside the entire BMOPO length L, and (c) zoomed in the first L/10 length, near the maximum of the three-wave interaction (at time t = 254 ps).

Fig. 4 .

 4 Fig. 4. Parametric downconversion efficiency [pump depletion 1-I p (L)/I p (0)] versus input pump spectral width : (a) for the one piece perfect PPGaN BMOPO, and (b) for the 64 times segmented PPGaN BMOPO with the random phase perturbation in the junctions of amplitude 0 < δΦ m (x) ≤ π/4. The linear fitting presents a decreasing slope of -0.0375 and relative efficiency of 93% with respect to the one piece perfect PPGaN BMOPO of Fig 4(a).

Fig. 5 .

 5 Fig. 5. Parametric downconversion efficiency [pump depletion 1-I p (L)/I p (0)] versus input pump spectral width, for the 64 times segmented PPGaN BMOPO with the random phase perturbation in the junctions of amplitude 0 < δΦ m (x) ≤ π/3. The linear fitting presents a decreasing slope of -0.0369 and relative efficiency of 87% with respect to the one piece perfect PPGaN BMOPO (upper curve).

Fig. 6 .

 6 Fig.6. Parametric downconversion efficiency [pump depletion 1-I p (L)/I p (0)] versus input pump spectral width, for the 64 times segmented PPGaN BMOPO with the random phase perturbation in the junctions of amplitude 0 < δΦ m (x) ≤ π/2. The linear fitting presents a decreasing slope of -0.0320 and relative efficiency of 63% with respect to the one piece perfect PPGaN BMOPO (upper curve).

Fig. 7 .

 7 Fig. 7. Stitched BMOPO (segmented in 64 elements with maximum perturbation δφ m (x) = π/4) for a stochastic pump input of FWHM Δν p (0) = 4.36 THz, (a) Output incoherent pump spectrum fitted by a Gaussian envelope of FWHM Δν p (L) = 4.25 THz, (b) backward signal spectrum with Δν s (0) = 16.30 GHz, and (c) forward idler spectrum with a Gaussian fitting of FWHM Δν i (L) = 4.23 THz. The coherence gain attains Δν p /Δν s = 312 for a pump depletion or efficiency of 0.0944 (cf.Table 1).

Fig. 8 .

 8 Fig. 8. Stitched BMOPO (segmented in 64 elements with maximum perturbation δφ m (x) = π/3) for a stochastic pump input of FWHM Δν p (0) = 4.36 THz : (a) Output incoherent pump spectrum fitted by a Gaussian envelope of FWHM Δν p (L) = 4.51 THz, (b) backward signal spectrum with Δν s (0) = 12.97 GHz, and (c) forward idler spectrum with a Gaussian fitting of FWHM Δν i (L) = 4.62 THz.The coherence gain attains Δν p /Δν s = 347, for a pump depletion or efficiency of 0.0759 (cf. Table1).

Fig. 9 .

 9 Fig. 9. Stitched BMOPO (segmented in 64 elements with maximum perturbation δφ m (x) = π/2) for a stochastic pump input of FWHM Δν p (0) = 4.36 THz : (a) Output incoherent pump spectrum fitted by a Gaussian envelope of FWHM Δν p (L) = 4.54 THz, (b) backward signal spectrum with Δν s (0) = 9.72 GHz, and (c) forward idler spectrum with a Gaussian fitting of FWHM Δν i (L) = 4.62 THz.The coherence gain attains Δν p /Δν s = 467, for a pump depletion or efficiency of 0.0256 (cf. Table1).

Table 1

 1

).

Acknowledgements

The authors thank Pascal Baldi and Krishna Thyagarajan for stimulating discussions.

This work was granted access to the HPC and visualization resources of "Centre de Calcul Interactif" hosted by "Université Nice Sophia Antipolis".