Fatiha Alabau-Boussouira

Matthieu Léautaud
email: leautaud@ann.jussieu.fr

Indirect controllability of locally coupled systems under geometric conditions

We consider systems of two wave/heat/Schrödinger-type equations coupled by a zero order term, only one of them being controlled. We prove an internal and a boundary null-controllability result in any space dimension, provided that both the coupling and the control regions satisfy the Geometric Control Condition. This includes several examples in which these two regions have an empty intersection.

e iθ u 1 -∆ c u 1 + au 1 + δpu 2 = b f dans (0, T ) × Ω, e iθ u 2 -∆ c u 2 + au 2 + pu 1 = 0 dans (0, T ) × Ω, u 1 = u 2 = 0 sur (0, T ) × ∂Ω, (u 1 , u 2 )| t=0 = (u 0 1 , u 0 2 ) dans Ω,

avec θ = 0, ont été étudiées intensivement. Ici, a, b, c, p sont des fonctions réelles régulières de x ∈ Ω, avec b ≥ 0 et p ≥ 0, δ > 0 est un paramètre, -∆ c est un opérateur autoadjoint uniformément elliptique sur Ω, et f est le contrôle. Le résultat général concernant ces systèmes, prouvé par différentes méthodes dans [START_REF] De Teresa | Insensitizing controls for a semilinear heat equation[END_REF][START_REF] Ammar-Khodja | Null controllability of some reaction-diffusion systems with one control force[END_REF][START_REF] González-Burgos | Controllability results for some nonlinear coupled parabolic systems by one control force[END_REF][START_REF] Léautaud | Spectral inequalities for non-selfadjoint elliptic operators and application to the null-controllability of parabolic systems[END_REF] est un théorème de contrôlabilité à zéro dès qu'on suppose {p > 0} ∩ {b > 0} ∅. Qu'en est-il du cas {p > 0} ∩ {b > 0} = ∅? Le second problème auquel on s'intéresse ici est le problème du contrôle frontière associé à (1) (cf.( 4)). S'il semble résolu en dimension 1 d'espace pour a et p constants (voir [START_REF] Fernández-Cara | Boundary controllability of parabolic coupled equations[END_REF]), il reste complètement ouvert en dimension supérieure ou pour des coefficients variables. Pour ces deux problèmes, il semblerait que la théorie des équations paraboliques et les outils utilisés se heurtent pour le moment à de sérieuses difficultés. D'autre part, on sait depuis [START_REF] Russell | A unified boundary controllability theory for hyperbolic and parabolic partial differential equations[END_REF] que les propriétés de contrôlabilité des équations hyperboliques se transmettent aux équations paraboliques. Poursuivant l'étude initiée dans [START_REF] Alabau-Boussouira | A two-level energy method for indirect boundary observability and controllability of weakly coupled hyperbolic systems[END_REF], nous répondons aux deux questions ci-dessus pour le système hyperbolique (3) (consistant à remplacer e iθ u j par u j , pour j = 1, 2, dans (1)) sous certaines hypothèses. Nous en déduisons ensuite une réponse partielle à ces questions pour [START_REF] Alabau-Boussouira | Indirect boundary stabilization of weakly coupled systems[END_REF].

Dans le cadre du contrôle frontière, on renvoie aux systèmes (4) et [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF] ci-dessous, pour lesquelles le contrôle agit par la condition au bord

u 1 | ∂Ω = b ∂ f , avec b ∂ ∈ C ∞ c (∂Ω), b ∂ ≥ 0.
On appellera GCC (resp. GCC ∂ ) la condition de contrôle géométrique interne (resp. frontière) de [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF], que nous rappelons dans la section 2.

Pour formuler nos résultats, on utilisera les hypothèses suivantes: (i) L'Opérateur -∆ c + a est uniformément coercif sur Ω.

(u 0 1 , u 0 2 , u 1 1 , u 1 2 ) ∈ H 1 0 (Ω) × H 2 ∩ H 1 0 (Ω) × L 2 (Ω) × H 1 0 (Ω) (resp. (u 0 1 , u 0 2 , u 1 1 , u 1 2 ) ∈ L 2 (Ω) × H 1 0 (Ω) × H -1 (Ω) × L 2 (Ω)), il existe un contrôle f ∈ L 2 ((0, T ) × Ω) (resp. f ∈ L 2 ((0, T ) × ∂Ω)) tel que la solution de (3) (resp. (5)) vérifie (u 1 , u 2 , u 1 , u 2 )| t=T = 0.
On remarque que les espaces dans lesquels u 1 et u 2 sont contrôlés ne sont pas les mêmes, ce qui est naturel. Ce résultat est montré dans un cadre abstrait (voir Section 3), incluant aussi des systèmes de plaques couplées. On en déduit, grâce aux méthodes de transmutation de [START_REF] Russell | A unified boundary controllability theory for hyperbolic and parabolic partial differential equations[END_REF][START_REF] Miller | Controllability cost of conservative systems: resolvent condition and transmutation[END_REF][START_REF] Phung | Observability and control of Schrödinger equations[END_REF][START_REF] Miller | The control transmutation method and the cost of fast controls[END_REF] les résultats suivants pour deux équations de diffusion ou de Schrödinger couplées.

Theorem 0.2 (Systèmes d'équations de diffusion). On suppose que (i) est satisfaite, que ω p satisfait GCC, et que ω b (resp. Γ b ) satisfait GCC (resp. GCC ∂ ). Il existe alors une constante δ * > 0 telle que pour tout (δ, p + ) satisfaisant

√ δp + < δ * , pour tout T > 0, θ ∈ (-π/2, π/2), pour tous p, b (resp. b ∂ ) satisfaisant (ii) et (iii), et toutes données initiales (u 0 1 , u 0 2 ) ∈ L 2 (Ω) 2 (resp. (u 0 1 , u 0 2 ) ∈ H -1 (Ω) 2 ), il existe un contrôle f ∈ L 2 ((0, T ) × Ω) (resp. f ∈ L 2 ((0, T ) × ∂Ω)) tel que la solution de (1) (resp. (4)) vrifie (u 1 , u 2 )| t=T = 0.
Ce résultat répond donc partiellement aux deux questions posées, fournissant de nombreux exemples pour lesquels {p > 0}∩{b > 0} = ∅ ainsi qu'un résultat de contrôle frontière en toute dimension d'espace pour un couplage variable. Cependant, on notera que les hypothèses géométriques ne sont pas naturelles pour des équations paraboliques. Par conséquent, le théoreme 0.2 n'est qu'un premier pas pour cette étude.

Theorem 0.3 (Systèmes d'équations de Schrödinger). Sous les hypothèse du théorème 0.2, le même résultat de contrôlabilité est valide pour le système (1) (resp. (4)) avec θ = ±π/2, pour des données initiales

(u 0 1 , u 0 2 ) ∈ L 2 (Ω) × H 1 0 (Ω) (resp. (u 0 1 , u 0 2 ) ∈ H -1 (Ω) × L 2 (Ω)).

Introduction

During the last decade, the controllability properties of coupled parabolic equations like

               e iθ u 1 -∆ c u 1 + au 1 + δpu 2 = b f in (0, T ) × Ω, e iθ u 2 -∆ c u 2 + au 2 + pu 1 = 0 in (0, T ) × Ω, u 1 = u 2 = 0 on (0, T ) × ∂Ω, (u 1 , u 2 )| t=0 = (u 0 1 , u 0 2 ) in Ω, (2)
with θ = 0, have been intensively studied. Here, ∆ c is a selfadjoint elliptic operator, and all the parameters are precisely defined in Section 2. The null-controllability problem under view is the following: given a time T > 0 and initial data, is it possible to find a control function f so that the state has been driven to zero in time T ? It has been proved in [START_REF] De Teresa | Insensitizing controls for a semilinear heat equation[END_REF][START_REF] Ammar-Khodja | Null controllability of some reaction-diffusion systems with one control force[END_REF][START_REF] González-Burgos | Controllability results for some nonlinear coupled parabolic systems by one control force[END_REF][START_REF] Léautaud | Spectral inequalities for non-selfadjoint elliptic operators and application to the null-controllability of parabolic systems[END_REF] with different methods that System (2) is null-controllable as soon as {p > 0} ∩ {b > 0} ∅. In these works, the case {p > 0} ∩ {b > 0} = ∅ has been left as an open problem. However, Kavian and de Teresa [START_REF] Kavian | Unique continuation principle for systems of parabolic equations[END_REF] have proved for a cascade system (i.e. taking δ = 0 in (2)) that approximate controllability holds. The natural question is then whether or not nullcontrollability (which is a stronger property) still holds in the case {p > 0} ∩ {b > 0} = ∅?

The second Problem under interest here is the boundary controllability of systems like (2) (or more precisely System (4) below). The recent work [START_REF] Fernández-Cara | Boundary controllability of parabolic coupled equations[END_REF] studies slightly more general systems in one space dimension and with constant coupling coefficients. The cases of higher space dimensions and varying coupling coefficients (and in particular when the coefficients vanish in a neighborhood of the boundary) are to our knowledge completely open.

Concerning these two open problems, it seems that the parabolic theory and associated tools encounter for the moment some essential difficulties.

On the other hand, it is known from [START_REF] Russell | A unified boundary controllability theory for hyperbolic and parabolic partial differential equations[END_REF] that controllability properties can be transferred from hyperbolic equations to parabolic ones. And it seems, at least for boundary control problems, that the theory for coupled hyperbolic equations of the type

               u 1 -∆ c u 1 + au 1 + δpu 2 = b f in (0, T ) × Ω, u 2 -∆ c u 2 + au 2 + pu 1 = 0 in (0, T ) × Ω, u 1 = u 2 = 0 on (0, T ) × ∂Ω, (u 1 , u 2 , u 1 , u 2 )| t=0 = (u 0 1 , u 0 2 , u 1 1 , u 1 2 ) in Ω, (3)
is better understood (see [START_REF] Alabau-Boussouira | A two-level energy method for indirect boundary observability and controllability of weakly coupled hyperbolic systems[END_REF]), even less studied. In the case of varying coefficients and several space-dimensions, the associated stabilization problem is addressed in [START_REF] Alabau-Boussouira | Indirect boundary stabilization of weakly coupled systems[END_REF][START_REF] Alabau-Boussouira | Indirect stabilization of locally coupled wave-type systems[END_REF].

In the present work, we answer these two questions for hyperbolic problems improving the results of [START_REF] Alabau-Boussouira | A two-level energy method for indirect boundary observability and controllability of weakly coupled hyperbolic systems[END_REF][START_REF] Alabau-Boussouira | Indirect stabilization of locally coupled wave-type systems[END_REF], and then deduce a (partial) solution to the two open questions raised above. Indeed, we prove that Systems (2)-( 3) are null-controllable (in appropriate spaces) as soon as {p > 0} and {b > 0} both satisfy the Geometric Control Condition (recalled below) and √ δ p L ∞ (Ω) satisfies a smallness assumption. This contains several examples with {p > 0} ∩ {b > 0} = ∅ in any space-dimension, and partially answers to the first question. We prove as well that the same controllability result holds for boundary control, which partially answers to the second question. Of course, the geometric conditions needed here are essential (and even sharp) for coupled waves, but inappropriate for parabolic equation. However, this is a first step towards a better understanding of these types of systems. In one space dimension in particular, the geometric conditions are reduced to a non-emptiness condition and are hence optimal for parabolic systems as well. Similar results have been obtained simultaneously in [START_REF] Rosier | Exact controllability of a cascade system of conservative equations[END_REF] with different methods in one space dimension for cascade systems.

In this note, we first state our results for wave/heat/Shrödinger-type Systems. Then, we introduced an abstract setting adapted to these problems and give some elements of the proof in this context.

Acknowledgements. The authors want to thank B. Dehman for discussions on the article [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF], S. Ervedoza for having pointed out the papers [START_REF] Miller | Controllability cost of conservative systems: resolvent condition and transmutation[END_REF][START_REF] Miller | The control transmutation method and the cost of fast controls[END_REF], and L. Miller for discussions on these two articles. The first author would like to thank the Fondation des Sciences Mathématiques de Paris and the organizers of the IHP trimester on control of PDE's for their support. The second author wishes to thank O. Glass and J. Le Rousseau for very fruitful discussions and encouragements. Both authors were partially supported by l'Agence Nationale de la Recherche under grant ANR-07-JCJC-0139-01 and the GDRE CONEDP (CNRS/INDAM/UP).

Main results

Let Ω be a bounded domain in R n with smooth (say C ∞ ) boundary (or a smooth connected compact riemannian manifold with or without boundary) and ∆ c = div(c∇) a (negative) Laplace operator (or Laplace Beltrami operator with respect to the riemannian metric) on Ω. Here, c denotes a smooth (say C ∞ ) positive symmetric matrix i.e. in particular

C -1 0 |ξ| 2 ≤ c(x)ξ • ξ ≤ C 0 |ξ| 2 for some C 0 > 1, for all x ∈ Ω, ξ ∈ R n . We consider the control problems (2) with θ ∈ [-π/2, π/2],
including Schrödinger-type systems for θ = ±π/2 and diffusion-type systems for θ ∈ (-π/2, π/2), and (3) consisting in a wave-type system, with only one control force. In these systems, a = a(x), p = p(x) and b = b(x) are smooth real-valued functions on Ω, δ > 0 is a constant parameter and f is the control function, that can act on the system.

We shall also consider the same systems controlled from the boundary through the (smooth) real-valued function b ∂ :

               e iθ u 1 -∆ c u 1 + au 1 + δpu 2 = 0 in (0, T ) × Ω, e iθ u 2 -∆ c u 2 + au 2 + pu 1 = 0 in (0, T ) × Ω, u 1 = b ∂ f , u 2 = 0 on (0, T ) × ∂Ω, (u 1 , u 2 )| t=0 = (u 0 1 , u 0 2 ) in Ω, (4)
               u 1 -∆ c u 1 + au 1 + δpu 2 = 0 in (0, T ) × Ω, u 2 -∆ c u 2 + au 2 + pu 1 = 0 in (0, T ) × Ω, u 1 = b ∂ f , u 2 = 0 on (0, T ) × ∂Ω, (u 1 , u 2 , u 1 , u 2 )| t=0 = (u 0 1 , u 0 2 , u 1 1 , u 1 2 ) in Ω. (5)
We first notice that, on the space L 2 (Ω) 2 endowed with the inner product (u, v 2 , in the sense of semigroup theory. Then, taking f ∈ L 2 ((0, T ) × ∂Ω) the initial-boundary value problem (4), resp. ( 5), is well-posed in 2 , in the sense of transposition solution (see [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF]).

) δ = (u 1 , v 1 ) L 2 (Ω) + δ(u 2 , v 2 ) L 2 (Ω) , the operator A δ = -∆ c + a δp p -∆ c + a , with domain D(A δ ) = H 2 (Ω) ∩ H 1 0 (Ω) 2 , is selfadjoint. As a consequence, for f ∈ L 2 ((0, T ) × Ω), the Cauchy problem (2), resp. (3), is well-posed in L 2 (Ω) 2 , resp. H 1 0 (Ω) 2 × L 2 (Ω)
H -1 (Ω) 2 , resp. L 2 (Ω) 2 × H -1 (Ω)
An important remark to make before addressing the controllability problem is concerned with the regularity of solutions of (3)- [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF]. If one takes for system (3) (resp. ( 5)) an initial data (u 0 1 , u 0

2 , u 1 1 , u 1 2 ) ∈ H 1 0 (Ω) × H 2 ∩ H 1 0 (Ω) × L 2 (Ω) × H 1 0 (Ω) (resp. (u 0 1 , u 0 2 , u 1 1 , u 1 2 ) ∈ L 2 (Ω) × H 1 0 (Ω) × H -1 (Ω) × L 2 (Ω)), and a control f ∈ L 2 ((0, T ) × Ω) (resp. f ∈ L 2 ((0, T ) × ∂Ω)), then the state (u 1 , u 2 , u 1 , u 2 ) remains in the space H 1 0 (Ω) × H 2 ∩ H 1 0 (Ω) × L 2 (Ω) × H 1 0 (Ω) (resp. L 2 (Ω) × H 1 0 (Ω) × H -1 (Ω) × L 2 (Ω)
) for all time. Recall that for Systems (3) and ( 5), the null-controllability is equivalent to the exact controllability. As a consequence, it is not possible, taking for instance zero as initial data to reach any target state in 2 ). The controllability question for (3)-( 5) hence becomes: Starting from rest at time t = 0, is it possible to reach all

H 1 0 (Ω) 2 × L 2 (Ω) 2 (resp. L 2 (Ω) 2 × H -1 (Ω)
H 1 0 (Ω) × H 2 ∩ H 1 0 (Ω) × L 2 (Ω) × H 1 0 (Ω) (resp. L 2 (Ω) × H 1 0 (Ω) × H -1 (Ω) × L 2 (Ω)) in time t = T sufficiently large?
The strategy we adopt here is to prove some controllability results for the hyperbolic systems (3) and ( 5), extending the two-levels energy method introduced in [START_REF] Alabau-Boussouira | A two-level energy method for indirect boundary observability and controllability of weakly coupled hyperbolic systems[END_REF]. Then, using transmutation techniques, we deduce controllability properties of (2) and (4).

To state our results, we recall the classical Geometric Control Conditions GCC (resp. GCC ∂ ), which, according to [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF], is a necessary and sufficient condition for the internal (resp. boundary) observability and controllability of one wave equation. We say that ω ⊂ Ω satisfies GCC (resp. Γ ⊂ ∂Ω satisfies GCC ∂ ) if every generalized geodesic traveling at speed one in Ω meets ω (resp. meets Γ on a non-diffractive point) in finite time.

We shall make the following key assumptions: (i) We have (-∆ c + a)u, u L 2 (Ω) ≥ λ 0 u 2 L 2 (Ω) , for some λ 0 > 0, for all u ∈ L 2 (Ω). In the case where c = Id and a = 0, the best constant λ 0 is 1/C 2 P , where C P is the Poincaré's constant of Ω. We shall also require that the operator A δ satisfies, for

C > 0, A δ (v 1 , v 2 ), (v 1 , v 2 ) δ ≥ C v 1 2 H 1 0 (Ω) + δ v 2 2 H 1 0 (Ω) for all (v 1 , v 2 ) ∈ D(A δ )
. This is the case when assuming

√ δp + < λ 0 .
Theorem 2.1 (Wave-type systems). Suppose that (i) holds, that ω p satisfies GCC and that ω b (resp. Γ b ) satisfies GCC (resp. GCC ∂ ). Then, there exists a constant δ * > 0 such that for all (δ, p + ) satisfying √ δp + < δ * , there exists a time T * > 0 such that for all T > T * , all p, b (resp. b ∂ ) satisfying (ii) and (iii), and all initial data

(u 0 1 , u 0 2 , u 1 1 , u 1 2 ) ∈ H 1 0 (Ω) × H 2 ∩ H 1 0 (Ω) × L 2 (Ω) × H 1 0 (Ω) (resp. (u 0 1 , u 0 2 , u 1 1 , u 1 2 ) ∈ L 2 (Ω) × H 1 0 (Ω) × H -1 (Ω) × L 2 (Ω)), there exists a control f ∈ L 2 ((0, T ) × Ω) (resp. f ∈ L 2 ((0, T ) × ∂Ω)) such that the solution of (3) (resp. (5)) satisfies (u 1 , u 2 , u 1 , u 2 )| t=T = 0.
Another way to formulate this result is to say that, under the assumptions of Theorem 2.1, the reachable set at time T > T * with zero initial data is exactly

H 1 0 (Ω) × H 2 ∩ H 1 0 (Ω) × L 2 (Ω) × H 1 0 (Ω) in the case of L 2 internal control and L 2 (Ω) × H 1 0 (Ω) × H -1 (Ω) × L 2 (Ω)
in the case of L 2 boundary control. Some comments should be made about this result. First this is a generalization of the work [START_REF] Alabau-Boussouira | A two-level energy method for indirect boundary observability and controllability of weakly coupled hyperbolic systems[END_REF] where the coupling coefficients considered have to satisfy p(x) ≥ C > 0 for all x ∈ Ω. The geometric situations covered by Theorem 2.1 are richer, and include in particular several examples of coupling and control regions that do not intersect. Second, the coercivity assumption (i) for -∆ c + a together with the smallness assumption on √ δp + seem to be only technical and inherent to the method we use here. Note by the way that this smallness assumption contains the coercivity assumption for A δ , and allows to consider large p + or large δ (provided that the other is small enough). Moreover, the control time T * we obtain depends on all the parameters of the system, and not only the sets ω p , ω b and Γ b (as it is the case for a single wave equation). Finally, the fact that we consider twice the same elliptic operator ∆ c is a key point in our proof and it is likely that this result does not hold for waves with different speeds (see [START_REF] Alabau-Boussouira | A two-level energy method for indirect boundary observability and controllability of weakly coupled hyperbolic systems[END_REF] for results with different speeds and different operators).

As a consequence of Theorem 2.1 and using transmutation techniques (due to [START_REF] Russell | A unified boundary controllability theory for hyperbolic and parabolic partial differential equations[END_REF][START_REF] Miller | Controllability cost of conservative systems: resolvent condition and transmutation[END_REF] for heat-type equations and to [START_REF] Phung | Observability and control of Schrödinger equations[END_REF][START_REF] Miller | The control transmutation method and the cost of fast controls[END_REF] for Schrödinger-type equations), we can now state the associated results for Systems (2) and (4). Corollary 2.2 (Heat-type systems). Suppose that (i) holds, that ω p satisfies GCC and that ω b (resp. Γ b ) satisfies GCC (resp. GCC ∂ ). Then, there exists a constant δ * > 0 such that for all (δ, p + ) satisfying √ δp + < δ * , for all T > 0, θ ∈ (-π/2, π/2), for all p, b (resp. b ∂ ) satisfying (ii) and (iii), and all initial data

(u 0 1 , u 0 2 ) ∈ L 2 (Ω) 2 (resp. (u 0 1 , u 0 2 ) ∈ H -1 (Ω) 2 ), there exists a control f ∈ L 2 ((0, T ) × Ω) (resp. f ∈ L 2 ((0, T ) × ∂Ω)) such that the solution of (2) (resp. (4)) satisfies (u 1 , u 2 )| t=T = 0.
Corollary 2.3 (Schrödinger-type systems). Under the assumptions of Corollary 2.2, the same null-controllability result holds for System (2) (resp. (4)) with θ = ±π/2, taking initial data [START_REF] Miller | The control transmutation method and the cost of fast controls[END_REF]Theorem 3.4] and the smoothing effect of the heat equation. Corollary 2.3 is a direct consequence of Theorem 2.1, combined with [START_REF] Miller | Controllability cost of conservative systems: resolvent condition and transmutation[END_REF]Theorem 3.1]. Since there is no smoothing effect in this case, we still obtain a controllability result in asymmetric spaces here: the uncontrolled variable u 2 has to be more regular than the other one. This shows that the attainable set from zero for a L 2 internal control (resp. 2 ) is reachable for (2) (resp. ( 4)) with θ = ±π/2 remains open.

(u 0 1 , u 0 2 ) ∈ L 2 (Ω)× H 1 0 (Ω) (resp. (u 0 1 , u 0 2 ) ∈ H -1 (Ω)× L 2 (Ω)). Corollary 2.2 is a direct consequence of Theorem 2.1, combined with
L 2 boundary control) contains L 2 (Ω) × H 1 0 (Ω) (resp. H -1 (Ω) × L 2 (Ω)). Whether or not a general target in L 2 (Ω) 2 (resp. H -1 (Ω)

Abstract setting and ingredients of proof

In this section, we describe the abstract setting (already used in [START_REF] Alabau-Boussouira | Indirect stabilization of locally coupled wave-type systems[END_REF]) in which we prove Theorem 2.1 for Systems (3)-( 5), and define the appropriate spaces and operators. Let H be a Hilbert space and (A, D(A)) a selfadjoint positive operator on H with compact resolvent. We denote by (•, •) H the inner product on H and • H the associated norm. For k ∈ N, we set

H k = D(A k 2 ) endowed with the inner product (•, •) H k = (A k 2 •, A k 2 •) H and associated norm • H k = A k 2 • H .
We define H -k as the dual space of H k with respect to the pivot space H = H 0 , and

• H -k = A -k 2 • H is the norm on H -k .
The operator A can be extended to an isomorphism from H k to H k-2 for any k ≤ 1, still denoted by A. We denote by λ 0 > 0 the largest constant satisfying v 2 H 1 ≥ λ 0 v 2 H for all v ∈ H 1 , that is, the smallest eigenvalue of the selfadjoint positive operator A. We consider that the coupling operator P is bounded on H and denote by P * is its adjoint, p + := P L(H) = P * L(H) . In the following, as in [START_REF] Alabau-Boussouira | A two-level energy method for indirect boundary observability and controllability of weakly coupled hyperbolic systems[END_REF], we shall make use of the different energy levels e k (ϕ(

t)) = 1 2 ϕ(t) 2 H k + ϕ (t) 2 H k-1 , k ∈ Z which are all preserved through time if ϕ is a solution of ϕ + Aϕ = 0.
Before addressing the control problem, let us introduce the adjoint system

         v 1 + Av 1 + δPv 2 = 0, v 2 + Av 2 + P * v 1 = 0, (v 1 , v 2 , v 1 , v 2 )| t=0 = (v 0 1 , v 0 2 , v 1 1 , v 1 2 ) (6)
which shall stand for our observation system. This system can be recast as a first order differential equation V = A δ V, V(0) = V 0 , where

A δ = 0 Id -A δ 0 , A δ = A δP P * A , V = (v 1 , v 2 ), V = (V, V ) = (v 1 , v 2 , v 1 , v 2 ).
Note that the operator A δ is selfadjoint on the space H × H endowed with the weighted inner product V,

Ṽ δ = (v 1 , ṽ1 ) H + δ(v 2 , ṽ2 ) H . Since we have A δ V, V δ = (Av 1 , v 1 ) H +δ(Av 2 , v 2 ) H +2δ(Pv 2 , v 1 ) H ≥ 1-p + √ δ λ 0 v 1 2 H 1 +δ v 2 2
H 1 , we shall suppose that p + √ δ < λ 0 , so that A δ is coercive. Under this assumption, A

1 2 δ V, A 1 2 δ
Ṽ δ defines an inner product on (H 1 ) 2 , equivalent to the natural one. Assuming that P, P * ∈ L(H k ) and writing

H k = (H k ) 2 ×H 2 k-1 , k ∈ Z, the operator A δ is an isomorphism from H k to H k-1 and is skewadjoint on H k , equipped with the inner product (U, V), ( Ũ, Ṽ) H k = (A k 2 δ U, A k 2 δ Ũ) δ + (A k-1 2 δ V, A k-1 2
δ Ṽ) δ . Note that this is an inner product according to the coercivity assumption for A δ , which is equivalent to the natural inner product of H k . Hence, A δ generates a group e tA δ on H k , and the homogeneous problem ( 6) is well-posed in these spaces. An important feature of solutions V(t) of System ( 6) is that all energies E k (V(t)) = 1/2 V(t) 2 H k are positive and preserved through time.

For this system, now studied in H 1 , we shall observe only the state of the first component, i.e. (u 1 , u 1 ), and hence define an observation operator B * ∈ L(H 2 × H, Y), where Y is a Hilbert space, standing for our observation space. This definition is sufficiently general to take into account both the boundary observation problem (taking B * ∈ L(H 2 , Y)) and the internal observation problem (taking B * ∈ L(H, Y)). We assume that B * is an admissible observation for one equation: For all T > 0 there exists a constant C > 0, such that all the solutions ϕ of ϕ + Aϕ = f satisfy ∫

T 0 B * (ϕ, ϕ ) 2 Y dt ≤ C e 1 (ϕ(0)) + e 1 (ϕ(T )) + ∫ T 0 e 1 (ϕ(t))dt + ∫ T 0 f 2 H dt . (A1)
Under this assumption, we have the following lemma.

Lemma 3.1. The operator B * is an admissible observation for (6). More precisely, for all T > 0, there exists a constant C > 0, such that all the solutions of (6) satisfy

T ∫ 0 B * (v 1 , v 1 )(t) 2 Y dt ≤ C e 1 (v 1 (0)) + e 0 (v 2 (0)) . (7)
Note that only the e 0 energy level of the second component v 2 is necessary in this admissibility estimate. Hence, we cannot hope to observe the whole H 1 energy of V and the best observability we can expect only involves e 0 (v 2 ). Our aim is now to prove this inverse inequality of [START_REF] González-Burgos | Controllability results for some nonlinear coupled parabolic systems by one control force[END_REF]. For this, we have to suppose some additional assumptions on the operators P and B * . Let us first precise Assumption (A2), related with the operator P:

We have Pv 2 H ≤ p + (Pv, v) H and there exists an operator Π P ∈ L(H), Π P L(H) = 1, and a number p

-> 0 such that (Pv, v) H ≥ p -Π P v 2 H ∀v ∈ H. (A2)
Note that p -≤ p + = P L(H) and that (A2) implies that the operators P and P * are non-negative. In the applications to coupled wave equations, P is the multiplication by the function p and the operator Π P is the multiplication by the characteristic function 1 ω p . Next, we shall suppose that a single equation is observable both by B * and by Π P in sufficiently large time:

∃T 0 > 0 such that for all T > T 0 there exists a constant C > 0, such that all solutions ϕ of ϕ + Aϕ = 0 satisfy both e 1 (ϕ(0

)) ≤ C ∫ T 0 B * (ϕ, ϕ ) 2 Y and e 1 (ϕ(0)) ≤ C ∫ T 0 Π P ϕ 2 H dt (A3)
In the context of Theorem 2.1, these observability assumptions are satisfied as soon as ω p and ω b satisfy GCC (resp. Γ b satisfies GCC ∂ ). We can now state (without proof) the main result of this note.

Theorem 3.2. Suppose that Assumptions (A1)-(A3) hold. Then there exists a constant δ * such that for all (δ, p + ) satisfying √ δp + < δ * , there exists a time T * such that for all T > T * there exists C > 0, such that for all V 0 ∈ H 1 , the solution V(t) = e tA δ V 0 of (6) satisfies e 1 (v 1 (0)) + e 0 (v 2 (0))

≤ C T ∫ 0 B * (v 1 , v 1 )(t) 2 Y dt. (8)
Applying the Hilbert Uniqueness Method (HUM) of [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF], we deduce now controllability results for the adjoint system. In this context, we have to define more precisely the observation operator. We shall treat two cases: First, B * (v 1 , v 1 ) = B * v 1 with B * ∈ L(H, Y), corresponding to internal observability (with Y = L 2 (Ω)), and second B * (v 1 , v 1 ) = B * v 1 with B * ∈ L(H 2 , Y), corresponding to boundary observability (with Y = L 2 (∂Ω)). In both cases, we define the control operator B as the adjoint of B * , and the control problem reads, for a control function f taking its values in Y,

         u 1 + Au 1 + δPu 2 = B f u 2 + Au 2 + P * u 1 = 0 (u 1 , u 2 , u 1 , u 2 )| t=0 = (u 0 1 , u 0 2 , u 1 1 , u 1 2 ). (9)
This is an abstract version of (3)- [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF]. Note that under this form, System (9) not only contains (3)-( 5), but also locally coupled systems of plate equations, with a distributed or a boundary control.

First case: B * (v 1 , v 1 ) = B * v 1 with B * ∈ L(H, Y). In this case, B ∈ L(Y, H) and the control problem ( 9) is well-posed in H 1 for f ∈ L 2 (0, T ; Y). In this setting, we first deduce from (8) the following other observability estimate for solutions W of ( 6) in H 0 : e 0 (w 1 (0)) + e -1 (w 2 (0)) ≤ C ∫ T 0 B * w 1 (t) 2 Y dt. The internal control result of Theorem 2.1 is then a direct consequence of the HUM since Assumptions (A1)-(A3) are satisfied in this application.

Second case: B * (v 1 , v 1 ) = B * v 1 with B * ∈ L(H 2 , Y). As a consequence of the admissibility inequality [START_REF] González-Burgos | Controllability results for some nonlinear coupled parabolic systems by one control force[END_REF], System (9) is well-posed in H 0 in the sense of transposition solutions. In this setting, the boundary control result of Theorem 2.1 is a direct consequence of the HUM and Theorem 3.2 since Assumptions (A1)-(A3) are satisfied in this application.

  (ii) On a {p > 0} ⊃ ω p pour un ouvert ω p ⊂ Ω et on pose p + := p L ∞ (Ω) . (iii) On a {b > 0} ⊃ ω b (resp. {b ∂ > 0} ⊃ Γ b ) pour un ouvert ω b ⊂ Ω (resp. Γ b ⊂ ∂Ω).

Theorem 0. 1 (

 1 Systèmes d'équations d'ondes). On suppose que (i) est satisfaite, que ω p satisfait GCC, et que ω b (resp. Γ b ) satisfait GCC (resp. GCC ∂ ). Il existe alors une constante δ * > 0 telle que pour tout (δ, p + ) satisfaisant √ δp + < δ * , il existe un temps T * > 0 tel que pour tout T > T * , tous p, b (resp. b ∂ ) satisfaisant (ii) et (iii), et toutes données initiales

  (ii) We have p ≥ 0 on Ω, {p > 0} ⊃ ω p for some open subset ω p ⊂ Ω and set p + := p L ∞ (Ω) ,. (iii) We have b ≥ 0 on Ω, {b > 0} ⊃ ω b (resp. b ∂ ≥ 0 on ∂Ω and {b ∂ > 0} ⊃ Γ b ) for some open subset ω b ⊂ Ω (resp. Γ b ⊂ ∂Ω).


Jordi Vilà-Valls

Laurent Ros

Jean-Marc Brossier

Joint oversampled carrier and time-delay synchronization in digital communications with large excess bandwidth $

Keywords: Carrier synchronization, time-delay synchronization, hybrid Cramér-Rao bound, nonlinear filtering, extended Kalman filter, expectation-maximization, oversampling, GNSS, BOC

This paper deals with the joint estimation of the pair dynamical carrier phase/Doppler shift and the time-delay, in a digital receiver. We consider a Binary Offset Carrier shaping function as used in satellite positioning, which is a time-limited pulse with a large excess bandwidth, and a Data Aided synchronization scenario, where we have a constant time-delay and a Brownian phase evolution with a linear drift. The proposed study is relative to the use of an oversampled signal model after matched filtering, leading to a colored reception noise and a non-stationary power signal. The contribution of this paper is twofold. First, we derive the Hybrid Cramér-Rao Bound for the joint phase/Doppler estimation problem. Then, we propose a method for the joint time-delay/carrier synchronization, which couples an Extended Kalman Filter and an Expectation-Maximization type algorithm. Our numerical results show the potential gain of using the oversampled signal for carrier synchronization, obtaining better performances than using a classical synchronizer, and good time-delay estimation.

Introduction

Synchronization is a fundamental part of global navigation satellite systems (GNSS). In the synchronization step, we estimate some parameters, such as carrier frequency, carrier phase and time-delay, between each visible satellite and the receiver to estimate the corresponding pseudorange. The synchronizer is coupled with a triangulation algorithm to obtain the receiver's position. In this paper, we focus our attention on the joint carrier phase/frequency offset and time-delay estimation problem in a GNSS-type receiver [START_REF] Parkinson | Global Positioning System: Theory and Applications[END_REF], where we consider a Binary Offset Carrier (BOC) shaping function (as used in new civil GNSS Galileo and the modernized GPS) [START_REF]Signal In Space Interface Control Document[END_REF] within a Data-Aided (DA) synchronization scenario.

Many estimation methods for joint carrier and timedelay synchronization have been proposed over the past decades. The time-delay, phase and frequency offset estimation problem is usually solved using Maximum Likelihood (ML) methods.

$ Part of this work was presented in conference [START_REF] Vilà-Valls | On-line Hybrid Cramér-Rao bound for oversampled dynamical phase and frequency offset estimation[END_REF].

Email address: jordi.vila-valls@upc.edu (Jordi Vilà-Valls) Historically, the delay and the frequency shift were assumed to be deterministic. In this context, assuming a known transmitted sequence, the optimal way to estimate these two parameters is to search for the maximum of the ambiguity function (delay/Doppler correlation method) [START_REF] Woodward | Probability and Information Theory, with Application to Radar[END_REF][START_REF] Van Trees | Detection, Estimation, modulation Theory, Part III: Radar-Sonar Signal processing and Gaussian Signals in noise[END_REF].

For time varying parameters, the previous method is still useful to provide an initialization to some tracking procedure and the synchronization is then performed as a two-step procedure: coarse and fine estimation, referring to acquisition and tracking, respectively. The acquisition system provides a first estimate of the time-delay and the Doppler shift and the tracking stage performs a local search for a fine estimation [START_REF]Understanding GPS. Principles and Applications[END_REF]. The tracking is usually based on Phase Locked Loops (PLL)/Delay Locked Loops (DLL) architectures [START_REF] Meyr | Digital Communication Receivers: Synchronization, Channel Estimation and Signal Processing[END_REF][START_REF] Mengali | Synchronization Techniques for Digital Receivers[END_REF]. These methods perform correctly with slowly varying phase errors.

For carrier synchronization, an alternative to the classical synchronization methods [START_REF] Mengali | Synchronization Techniques for Digital Receivers[END_REF][START_REF] Fitz | Planar filtered techniques for burst mode carrier synchronization[END_REF] is to reformulate the problem with a dynamic state-space model and to apply Bayesian estimation methods. Several contributions show the use of Kalman-type solutions [START_REF] Anderson | Optimal Filtering[END_REF] for carrier synchronization [START_REF] Kootsookos | An extended kalman filter for demodulation of polynomial phase signals[END_REF][START_REF] La Scala | Design of an extended kalman filter frequency tracker[END_REF][START_REF] Aghamohammadi | Adaptive synchronization and channel parameter estimation using an extended kalman filter[END_REF][START_REF] Cheng | GPS Carrier phase measurement method based on data transition detection in high dynamic circumstance[END_REF]. For the joint time-delay and carrier estimation, a solution to solve the ML problem within the state-space formulation is to use a recursive Expectation-Maximization (EM) type algorithm [START_REF] Dempster | Maximum likelihood from incompleter data via EM algorithm[END_REF]. Recently this approach has been applied to several problems such as channel estimation in OFDM systems [START_REF] Al-Naffouri | An EM-based forward-backward kalman filter for the estimation of time-variant channels in OFDM[END_REF][START_REF] Al-Naffouri | An forward-backward kalman filter-based STBC MIMO OFDM receiver[END_REF], speech recognition [START_REF] Digalakis | ML estimation of a stochastic linear system with the EM algorithm and its application to speech recognition[END_REF] and several parameter estimation and learning problems [START_REF] Khan | An expectation-maximization algorithm based Kalman smoother approach for event-related desynchronization (ERD) estimation from EEG[END_REF][START_REF] Chen | Direct and EM-based MAP sequence estimation with unknown time-varying channels[END_REF][START_REF] De Freitas | Nonlinear state space learning with EM and neural networks[END_REF].

When having an estimation problem, we need lower bounds on the estimation performance to be used as a benchmark. The family of Cramér-Rao Bounds (CRBs) has been shown to give accurate estimation lower bounds in many scenarios [START_REF] Van Trees | Detection, Estimation and Modulation Theory[END_REF]. For time-varying parameter estimation, an analytical expression of a general on-line recursive Bayesian CRB (BCRB) is given by Tichavský et al. [START_REF] Tichavský | Posterior Cramér-Rao bounds for discrete-time nonlinear filtering[END_REF] and the general formulation for the Hybrid CRB (HCRB), which applies when having both random and deterministic parameters, is derived by Bay et al. [START_REF] Bay | On the Hybrid Cramér Rao Bound and Its Application to Dynamical Phase Estimation[END_REF].

In the literature, most of the lower bounds and the corresponding algorithms assume a white observation noise and a stationary signal.

This contribution, extends and completes the work presented in [START_REF] Vilà-Valls | Oversampled phase tracking in digital communications with large excess bandwidth[END_REF]. We assume a time-varying phase offset modeled as a Brownian motion with a Doppler shift [START_REF] Amblard | Phase tracking: what do we gain from optimality?Particle filtering versus phase-locked loops[END_REF] and a constant time-delay over the observation window. We consider an oversampled (regarding the symbol time interval) signal model after receiver matched filtering. This implies dealing with a colored reception noise and taking into account the non-stationarity of the digital signal power (cyclostationarity when transmitting a random sequence).

Although this scenario is standard in satellite radiolocalization based on a Binary Offset Carrier (BOC) time-limited shaping pulse modulation, there is no theoretical study concerning the performance of the oversampled dynamical phase and frequency offset estimation, and joint time-delay and carrier estimation. In [START_REF] Vilà-Valls | Oversampled phase tracking in digital communications with large excess bandwidth[END_REF], we presented the derivation of a Bayesian CRB for the dynamical phase offset and the EKF that approaches this bound, both presented in a scenario similar to the one treated in this paper.

First, we derive a closed-form expression of the online HCRB for the dynamical phase and frequency offset estimation in the Data Aided (DA) scenario, assuming a Brownian phase evolution with a linear drift (Doppler shift). Secondly, we investigate the use of an EKF based algorithm which can approach this bound (more sophisticated methods such as particle filters or sigma-point Kalman filters are not necessary in this context [START_REF] Vilà-Valls | Oversampled phase tracking in digital communications with large excess bandwidth[END_REF]). We have thus to jointly estimate the colored noise, the dynamical phase and the frequency off-set. And finally, we propose an iterative block method to jointly estimate the carrier phase and the time-delay, coupling an EKF-based algorithm and an EM-type solution.

The study allows to measure the potential gain for carrier synchronization provided by the use of the fractionally-spaced processing after matched filtering, instead of the symbol time-spaced signal and the good performance obtained with the complete solution.

This paper is organized as follows. Section II sets the signal model. Section III, sets the estimation problem. Section IV, first recalls the HCRB expressions and derives the HCRB for this estimation problem, and then presents the EKF and derives the expressions of the filter in the oversampled phase and frequency offset estimation scenario. Section V presents the proposed method for joint time-delay and carrier synchronization. Finally, in section VI, the numerical results are presented and interpreted. The conclusion is given in Section VII.

Notations: the (k, l) th entry of a matrix A is denoted [A] k,l . E x denotes the expectation over x. ∇ θ and ∆ θ ψ represent the first and second-order partial derivatives operator i.e.,

∇ θ = ∂ ∂θ 1 • • • ∂ ∂θ K T and ∆ θ ψ = ∇ ψ ∇ T θ .

Signal Model

We propose the signal model for the transmission of a known sequence {a m } m∈Z over an Additive White Gaussian Noise (AWGN) channel affected by a dynamical carrier phase offset θ(t) (including the Doppler shift) and a time-delay τ(t). For an exhaustive derivation see [START_REF] Vilà-Valls | Oversampled phase tracking in digital communications with large excess bandwidth[END_REF].

Oversampled Signal Model

Discrete-time general formulation

The received complex baseband signal after matched filtering is

y(t) =        T c m a m p(t -mT c -τ(t))e iθ(t) + n(t)        * p * (-t),
(1) where T c , p(t) and n(t) stands for the symbol period, shaping pulse and circular complex Gaussian noise with a known two-sided power spectral density (psd) N 0 .

We assume a shaping pulse p(t) with support in [0, T c ], a constant time-delay within the observation window τ(t) = τ and a slowly varying phase evolution during a period T c , which is a usual assumption in satellite communications because the phase variation (due to oscillators phase noise, Doppler shifts, . . . ) within one 2 symbol period is small. In this case, the received signal can be written as [START_REF] Vilà-Valls | Oversampled phase tracking in digital communications with large excess bandwidth[END_REF]

y(t) = m a m g (t -mT -τ) e iθ(t) + b(t), (2)
where b(t) stands for the filtered colored noise, and g(t) = p(t) * p * (-t). If we consider a full digital synchronization architecture where the received signal is fractionally-spaced at fixed instants t k = k T S , where S is an integer oversampling factor, we can write the received oversampled signal as

y k = A k (τ)e iθ k + b k , (3)
where k refers to t k instants, so y k = y(t k ), θ k = θ(t k ) and b k = b(t k ), and

A k (τ) = m a m g k T S -τ -mT . (4)
Note that the noise b k is colored with variance σ 2 n , where

σ 2 n = N 0 × g(0)
T c is the variance of the AWGN n(t) measured in the noise equivalent bandwidth of the receiver filter p * (-t). We can define the symbol index q = k S , or equivalently, k = qS + s with s the subsymbol index (i.e. the position inside the symbol interval) and s = 0, • • • , S -1. {A k (τ)} k∈Z is a non-stationary power sequence for S > 1, even if {a m } m∈Z is a stationary power symbol sequence (a 2 m = 1).

Discrete-time re-formulation for the noise

The T S -spaced sequence of noise, b k k∈Z , is defined in the previous section from an analog noise n(t). Our motivation now is to replace this time series by another {b k } k∈Z with the same statistical properties, but which can be obtained entirely by a discrete-time formulation. This will be useful for the final state-space model formulation. As we did in [START_REF] Vilà-Valls | Oversampled phase tracking in digital communications with large excess bandwidth[END_REF], we can write that the noise samples b k have the same statistical properties than samples b k , which are obtained by a T S -spaced filtering of the time series n k :

b k = S -1 j=0 Π j n k-j-1 , (5)
where Π j are the coefficients of the filter which depend on the shaping pulse p(t).

Parameter Evolution Model

We consider a constant time-delay τ in the observation window. Concerning the carrier phase, in practice, we have a frequency shift between transmitter's and receiver's carrier oscillator and a Doppler shift due to the relative motion between the satellite and the receiver, so the phase offset is linear with time. We also must consider jitters introduced by oscillators imperfections which can be modeled as a random phase. To take it into account we suppose a Brownian phase offset evolution with a linear drift [START_REF] Amblard | Phase tracking: what do we gain from optimality?Particle filtering versus phase-locked loops[END_REF]:

θ k = θ k-1 + δ k + w k , (6)
δ k = δ k-1 , (7)
where k ≥ 2, w k is an i.i.d. zero-mean Gaussian noise sequences with known variance

σ 2 w
S and δ k is the unknown constant drift. Here σ 2 w stands for the variance growth of the phase noise in one symbol interval. We note that the variance of the Gaussian noise is directly related with the rapidity of evolution of the parameter. We note Σ the N × N covariance matrix of the phase offset evolution,

θ = [θ 1 • • • θ N ] T .

State-Space Model

When using an optimal filtering approach a statespace model formulation is needed. As we want to take into account that the observation noise on the output of the matched filter can be colored, we must include it into the state evolution

The state to be considered includes the phase offset, the drift and the colored noise:

x k = θ k δ k b k ν k • • • ν k-S +1 T , (8)
where

ν k ν k-1 • • • ν k-S +1
T is a sliding vector over an i.i.d noise sequence n k .

We define the state evolution matrix as

M K =                                     1 1 0 0 • • • 0 0 1 0 0 • • • 0 0 0 0 Π 0 Π S -1 0 0 0 0 0 • • • 0 . . . . . . . . . 1 . . . 0 0 0 • • • 0 1 0                                     , (9)
and the state noise as

w k = w k 0 0 n k 0 • • • 0 T .
From this we have that the state evolution and the observation equation of the state-space model are

x k = M K x k-1 + w k , (10)
y k = A k (τ) exp (iθ k ) + b k . (11)
We note that the state equation is linear and the observation equation depends non-linearly on the state. With this formulation we have no observation noise because we have included it in the state.

Estimation problem

In this section, we state the estimation problem and we introduce the proposed solution that we will develop in the following sections.

In general, the objective is to jointly estimate the constant time-delay τ and the states (including the carrier phase θ, which has a dynamical evolution, and the constant linear drift δ, which is hidden in the phase evolution) using the received signal y. We use the state-space model proposed in section 2 (eqs. [START_REF] Anderson | Optimal Filtering[END_REF][START_REF] Kootsookos | An extended kalman filter for demodulation of polynomial phase signals[END_REF]) and we consider that the transmitted symbol sequence is known at the receiver (DA synchronization scenario).

In a positioning context, we are interested in the timedelay τ to obtain the pseudorange estimation between each visible satellite and the receiver. In this case, we have to estimate the carrier phase and Doppler shift to obtain a correct estimation of the time-delay. If the received signal is not perturbed by a carrier phase error and a Doppler shift, the time-delay estimation problem can be solved with a simple correlation method. So what complicates the problem is the presence of these parameters.

If τ is known, the states can be inferred using a Kalman filter. Due to the presence of unobserved data (carrier phase and Doppler shift), the ML method to obtain a time-delay estimate can't be used because the computation of the likelihood function in a closed-form and its maximization w.r.t τ seems to be an intractable problem. To solve this problem we have to resort to iterative methods and the natural solution is to use an EM-type solution.

In the following, we first propose a method for phase and frequency shift estimation considering a known delay, and then we use this solution to propose a method for the joint carrier and time-delay estimation.

Carrier phase and frequency shift estimation

We consider in this section that we have a good timedelay synchronization (known time-delay), and we focus our attention on the joint carrier phase and frequency shift estimation. First, we compute the HCRB to be used as a benchmark on the estimation error, and then we propose a solution based on a Kalman-type algorithm.

Hybrid Cramér-Rao Bound

When dealing with an estimation problem we aim to know the ultimate accuracy that can be achieved by the estimator. The Cramér-Rao Bounds (CRB) provide a lower bound on the Mean Square Error (MSE) achievable by any unbiased estimator. Depending on the nature of the parameters to be estimated we use different bounds of the CRB family. If the vector of parameters is assumed to be deterministic we use the standard CRB and if the vector of parameters is random and an a priori information is available we use the so-called Bayesian CRB [START_REF] Bay | Analytic and Asymptotic Analysis of Bayesian Cramér-Rao Bound for Dynamical Phase Offset Estimation[END_REF]. When dealing with both random and deterministic parameters an Hybrid CRB (HCRB) is used [START_REF] Bay | On the Hybrid Cramér Rao Bound and Its Application to Dynamical Phase Estimation[END_REF]. The CRB suited to our problem is the HCRB as we want to estimate the phase offset evolution vector θ which is a random vector with an a priori probability density function (pdf) p(θ) and the linear drift δ which is a deterministic parameter.

In the on-line synchronization mode, at time k the receiver updates the observation vector y

= [y 1 • • • y k-1 ] T including the new observation y k to obtain the updated vector y = [y 1 • • • y k ] T in order to estimate θ k .
In this section we recall the expression of the Hybrid CRB and we present the closed-form expression of the HCRB for an oversampled dynamical phase and frequency offset estimation problem in a Data Aided scenario.

HCRB: background

We have a set of measurements y and we want to estimate an N-dimensional vector of parameters µ = µ T r µ T d T

. We consider the case where the random (µ r ) and the deterministic (µ d ) parts of the vector of parameters can be statistically dependent. We note µ * d the true value of µ d . The joint probability density of the pair (y, µ) is p y,µ (y, µ) and the a priori pdf of the random part of µ is p(µ r | µ * d ) p(µ r ). If μ(y) is our estimate of µ, the HCRB satisfies the following inequality on the MSE:

E y,µ|µ * d [ μ(y) -µ][ μ(y) -µ] T | µ * d ≥ H -1 (µ * d ), ( 12
)
where H(µ * d ) is the so-called Hybrid Information Matrix (HIM) defined as [START_REF] Bay | On the Hybrid Cramér Rao Bound and Its Application to Dynamical Phase Estimation[END_REF]

H(µ * d ) = E y,µ r |µ * d -∆ µ µ log p(y, µ r | µ d ) | µ * d . ( 13
)
Expanding the log-likelihood the HIM can be rewritten as

H(µ * d ) = E µ r |µ * d F(µ r , µ * d ) +E µ r |µ * d -∆ µ µ log p(µ r | µ d ) | µ * d ,
where F(µ r , µ * d ) is the Fisher Information Matrix (FIM) defined as

F(µ r , µ * d ) = E y|µ r ,µ * d -∆ µ µ log p(y | µ r , µ d ) | µ * d . ( 14
)
We can see that

H(µ * d ) = H D (µ * d ) + H P (µ * d )
, where the first term represents the average information about µ brought by the observations y and the second term represents the information available from the prior knowledge on µ, i.e., p(µ r |µ d ).

The N × N HCRB matrix can be written as

HCRB = H(µ * d ) -1 = H D (µ * d ) + H P (µ * d ) -1 , ( 15
)
where the k th element of the diagonal, [HCRB] k,k represents the lower bound on the estimation of

[µ] k from the observations block y = [y 1 • • • y N ].

HCRB: application to dynamical phase and frequency offset estimation

In this paragraph, a closed-form expression for the HCRB for an on-line fractionally-spaced phase-offset and linear drift estimation problem is presented. On the following we drop the dependence of the different matrices on µ * d = δ * for easier notation. As we consider a constant drift, for the derivation of the HCRB, we note δ k = δ.

We use the model presented in section 2 (eqs. ( 6),( 11)):

θ k = θ k-1 + δ + w k , y k = A k (τ) exp (iθ k ) + b k ,
where, as stated before, b k is a non-white noise with covariance matrix Γ. The index k refers to t k instants and A k (τ) are the coefficients specified in eq. ( 4) with τ = τ, so we can write that A k (τ) = A k .

Comparing this state-space model to the general model presented on the last paragraph, and supposing that we have N available measurements, we identify

µ r = θ = [θ 1 • • • θ N ] T and µ d = δ.
From this the HIM can be rewritten into a (N + 1) × (N + 1) block matrix as [START_REF] Bay | On the Hybrid Cramér Rao Bound and Its Application to Dynamical Phase Estimation[END_REF]

H = H 11 h 12 h 21 H 22 , (16)
where

H 11 = E y,θ|δ * -∆ θ θ log p(y | θ, δ) | δ * +E θ|δ * -∆ θ θ log p(θ | δ * ) , h 12 = h T 21 = E y,θ|δ * -∆ δ θ log p(y | θ, δ) | δ * +E θ|δ * -∆ δ θ log p(θ | δ * ) , H 22 = E y,θ|δ * -∆ δ δ log p(y | θ, δ) | δ * +E θ|δ * -∆ δ δ log p(θ | δ * ) .
So to compute the HIM we need the likelihood function and the a priori pdf. From the model we can write the log-likelihood as

log p(y | θ, δ * ) = log 1 π N |det(Γ)| -y -m H Γ -1 y -m , ( 17
)
where y is the N-dimensional received signal array and m is the mean vector of y, where the

k th component is [m] k = A k e iθ k . The a priori pdf is log p(θ | δ * ) = log p(θ 1 ) + (N -1) log 1 √ 2πσ w -N k=2 (θ k -θ k-1 -δ * ) 2 2σ 2 w . (18)
• Expression of H 11 : we can write that

H 11 = H D 11 + H P 11 , (19)
where

H D 11 = E y,θ|δ * -∆ θ θ log p(y | θ, δ) | δ * , H P 11 = E θ|δ * -∆ θ θ log p(θ | δ * ) .
The first term can be computed from eq.( 17). We note Λ(θ) = log p(y | θ, δ). The first derivative of Λ(θ) with respect to the l th phase parameter is

∂Λ(θ) ∂θ l = ∂ ∂θ l -[y -m] H Γ -1 [y -m] = ∂m H ∂θ l Γ -1 [y -m] + [y -m] H Γ -1 ∂m ∂θ l = 2 ∂m H ∂θ l Γ -1 [y -m] . (20
) If we compute now the derivative with respect to the k th phase parameter we have that

∂ 2 Λ(θ) ∂θ k ∂θ l = 2 ∂ 2 m H ∂θ k ∂θ l Γ -1 [y -m] - ∂m H ∂θ l Γ -1 ∂m ∂θ k . (21)
The (k, l) th element of the matrix H D 11 is

H D 11 k,l = E θ|δ * E y|θ,δ * -∂ 2 Λ(θ) ∂θ k ∂θ l = E θ|δ * 2 ∂m H ∂θ l Γ -1 ∂m ∂θ k .
We note that

∂m H ∂θ l = 0, • • • , 0, -iA * l e -iθ l , 0, • • • , 0 , (22)
∂m ∂θ k = 0, • • • , 0, iA k e iθ k , 0, • • • , 0 T , (23)
with the non-null values on the l th and k th position respectively, and so the coefficients can be written as

H D 11 k,l = E θ|δ * 2 A * l A k • Γ -1 k,l e j(θ k -θ l ) = 2 A * l A k • Γ -1 k,l E θ|δ * e j(θ k -θ l ) .
We can write that

E θ e i(θ k -θ l ) = E θ e i(u T kl θ) = φ (u kl ) , (24)
where u T kl = [0, • • • , 0, (+1), 0, • • • , 0, (-1), 0, • • • , 0], +1 in the k th position and -1 in the l th position of the array, φ(•) is the characteristic function of a Gaussian random variable θ:

φ (u kl ) = exp -1 2 u T kl Σ -1 u kl = exp -1 2 Σ -1 k,k + Σ -1 l,l -2 Σ -1 k,l , (25)
with Σ the covariance matrix of the phase evolution θ. Finally

H D 11 k,l = 2 A * l A k Γ -1 k,l e Ψ , (26)
where

Ψ = - 1 2 Σ -1 k,k + Σ -1 l,l -2 Σ -1 k,l . (27)
We note that the elements Σ -1 k,k are proportional to S /σ 2 w , so for small values of σ 2 w (σ 2 w < 0.1) we have that e Ψ ∼ 0 except when k = l where e Ψ = 1. As an example representing the worst of the cases, when we set S = 1, σ 2 w = 0.1, k = 1 and l = N, e Ψ = 4.54 • 10 -5 .

As we assume that the phase variation is small over the symbol interval (hypothesis done in section 2), we can consider that H D is a diagonal matrix with

H D 11 k,k = 2 | A k | 2 Γ -1 k,k . ( 28
)
On the following we compute the second term of eq. ( 19). From the state evolution eq. ( 6) and assuming that the initial phase θ 1 does not depend on δ, we have that

p(θ | δ * ) = p(θ 1 )Π N k=2 p(θ k | θ k-1 , δ * ), (29)
and due to this expansion we can rewrite the expression as

∆ θ θ ln p(θ, δ * ) = ∆ θ θ ln p(θ 1 ) + N k=1 ∆ θ θ ln p(θ k | θ k-1 , δ * ). ( 30
)
The first term in (eq.30) is a matrix with only one non-zero element, namely, the entry (1,1) which is equal to

∆ θ θ ln p(θ 1 ) 1,1 = ∂ 2 ln p(θ 1 ) ∂θ 2 1 . (31)
The other terms are matrices with only four nonzero elements, namely, the entries (k -1, k -1), (k -1, k), (k, k -1) and (k, k). Due to the Gaussian nature of the noise, one finds

∆ θ θ ln p(θ k | θ k-1 , δ * ) k,k = -S σ 2 w , (32)
∆ θ θ ln p(θ k | θ k-1 , δ * ) k,k-1 = S σ 2 w . (33)
The values for (k -1, k -1) and (k -1, k) are, respectively, the same that for (k, k) and (k -1, k).

Assuming that E θ 1 ∆ θ θ ln p(θ 1 ) = 0 that corresponds to the case of non-informative prior about θ 1 , we obtain that

H P 11 = E θ|δ * -∆ θ θ log p(θ | δ * ) = 1 σ 2 w /S                             1 -1 0 • • • 0 -1 2 -1 . . . . . . 0 . . . . . . . . . 0 . . . -1 2 -1 0 • • • 0 -1 1                             . ( 34
)
• Expression of h 12 : the log-likelihood (eq.17) does not depend on δ so the first term of h 12 is null, so

h 12 = E θ|δ * -∆ δ θ log p(θ | δ * ) .
From the state model we have that

h 12 = S σ 2 w 0 1×N-2 - S σ 2 w T .
• Expression of H 22 : as the log-likelihood does not depend on δ and using eq.( 18) we have that

H 22 = E θ|δ * -∆ δ δ log p(θ | δ * ) = E θ|δ * -∂ 2 ∂δ 2 log p(θ | δ * ) = S (N-1) σ 2 w .
Remarks: as we analyze the estimation problem in a DA scenario the bound depends on the transmitted sequence a. In this paper, we suppose the transmission of a known sequence to analyze the performance of the proposed algorithm and the bound. We note that the HCRB is computed for a specific known sequence. The bound depends on the sequence, the oversampling factor S and the position s inside the symbol interval of the current transmitted symbol (index M):

HCRB θ (a, S , s) = H -1 (a) N,N , (35)
HCRB δ (a, S , s) = H -1 (a) N+1,N+1 , (36)
with N = (M -1) * S + 1 + s.

Extended Kalman Filter for carrier phase and frequency shift estimation

In the sequel, we recall the notation and the basics of the Kalman filter and we derive the EKF [START_REF] Anderson | Optimal Filtering[END_REF] for the oversampled carrier phase and frequency offset estimation.

Kalman background

We consider that we have a system described by the following state-space equations pair

x k+1 = f k (x k ) + w k , y k = g k (x k ) + v k , (37)
where x k is the state vector, w k is a zero-mean white noise with covariance matrix Q k , y k is the observation vector at time k which is a partial and noisy observation of the state x k and v k is the observation noise with covariance matrix R k . Noises w k and v k are supposed to be uncorrelated. The functions f k (•) and g k (•) can be non-linear in a general case.

We note x k|m , the estimation of x k from the observations up to time m, x k|m = x kx k|m , the estimation error and P k|m = E x k|m x T k|m , the covariance matrix of the estimation error. For Gaussian, linear state models, the KF gives the best Mean Square Error (MSE) estimation of the state x k from observations up to time k. For non linear problems, the EKF gives a sub-optimal estimator x k|k in a recursive way: the main idea is to linearize the state-space equations at each iteration in order to transform the filtering problem into a usual Kalman one.

EKF: the algorithm

To derive the EKF, we need to compute ∂ f k (x k ) /∂x k and ∂g k (x k ) /∂x k .

In the state-space model for oversampled phase estimation presented in Section II (eqs. ( 10),( 11)), the state equation is linear, hence ∂ f k (x k ) /∂x k = M K . The state noise covariance Q is independent from k and has only two non-zero elements :

[Q] 1,1 = σ 2 w /S and [Q] 3,3 = σ 2
n . Because we introduced the colored noise b k into the state, there is no observation noise and the covariance matrix R is null. The observation equation is non-linear versus the state, we have to apply a linearization:

g = ∂g k x k|k-1 ∂x k = iA k (τ)e i θ k|k-1 0 1 0 • • • 0 T .
(38) Finally, the EKF expressions for the oversampled algorithm are:

                                           P k|k-1 = M K P k-1|k-1 M H K + Q x k|k-1 = M K x k-1|k-1 K k = P k|k-1 g H gP k|k-1 g H -1 P k|k = I -K k g P k|k-1 x k|k = x k|k-1 + K k y k -A k (τ)e i θ k|k-1 -b k|k-1 (
39) where I is the identity matrix with appropriate dimension.

Joint time-delay and carrier synchronization method

In this section, we propose an iterative block method for joint time-delay and carrier synchronization. The method is inspired by the EM algorithm [START_REF] Dempster | Maximum likelihood from incompleter data via EM algorithm[END_REF], which is an iterative method to find the ML estimate of given desired parameter in the presence of unobserved data or nuisance parameters. The idea behind the algorithm is to augment the observed data with latent data, which can be either missing data or parameter values.

Our method is then based on an iterative optimization of a cost function to find the time-delay (desired parameter in the EM formulation), coupled with the Kalmantype solution proposed in section 4 for the estimation of the carrier phase and the Doppler shift (nuisance parameters in the EM solution).

The proposed method

The main goal is to write a function that only depends on the time-delay, L(τ), that we will optimize iteratively. The starting point is the joint pdf,

p(y, θ; τ) = p(y|θ; τ)p(θ; τ), (40)
From the state-space model, the a priori density p (θ) is

p (θ) = p(θ 1 ) N k=2 p(θ k |θ k-1 ), ( 41
)
where p(θ k |θ k-1 ) are Gaussian densities with mean θ k-1 and variance σ 2 w , so we can write where m(θ, τ) is the mean vector with [m(θ, τ)] k = A k (τ)e iθ k and Γ is the covariance matrix of the observation noise. The product of these densities is

p (θ) = 1 √ 2πσ w N-1 exp        - 1 2σ 2 w N k=2 (θ k -θ k-1 ) 2        , (
p (y|θ; τ) p (θ) = 1 √ 2πσ w N-1 1 π N |det(Γ)| (43) × exp        -y -m(θ, τ) H Γ -1 y -m(θ, τ) - 1 2σ 2 w N k=2 (θ k -θ k-1 ) 2        , and so ln p (y, θ; τ) = C σ 2 w , Γ -y -m(θ, τ) H Γ -1 y -m(θ, τ) - 1 2σ 2 w N k=2 (θ k -θ k-1 ) 2 . ( 44
)
We assume that at the j th iteration of the method, the time-delay estimated at the previous step τ( j-1) , is available. As we want to obtain a function that only depends on τ, and optimize it with respect to this parameter, all the terms in the log-density that do not depend on the time-delay can be omitted. So we consider the function

F (θ, τ) = -y -m(θ, τ) H Γ -1 y -m(θ, τ) = -y H Γ -1 y + y H Γ -1 m(θ, τ) +m(θ, τ) H Γ -1 y -m(θ, τ) H Γ -1 m(θ, τ). ( 45
)
At the j th iteration, we obtain a function that only depends on the time-delay τ by taking the (conditional) expectation of F (θ, τ) with respect to the carrier phase θ (given the knowledge of y, a, τ( j-1) ),

L ( j) (τ) = E θ|y,a,τ ( j-1) [F (θ, τ)] = -y H Γ -1 y + y H Γ -1 E θ|y,a,τ ( j-1) [m(θ, τ)] +E θ|y,a,τ ( j-1) m(θ, τ) H Γ -1 y -E θ|y,a,τ ( j-1) m(θ, τ) H Γ -1 m(θ, τ) , (46)
where the terms E θ|y,a,τ ( j-1) [m(θ, τ)] and E θ|y,a,τ ( j-1) m(θ, τ) H Γ -1 m(θ, τ) only depend on the carrier phase. If we consider that we have an estimate of these two terms, we can estimate the time-delay at the current iteration from the previous time-delay estimate τ( j-1) , by maximizing the cost function

L (|) (τ) w.r.t τ τ( j) = arg max τ L(τ) ( j) . ( 47
)
The method iterates the optimization until the convergence of the sequence (..., τ( j) , τ( j+1) , ...).

Using the characteristic function of a Gaussian density, we can write that

E θ|y,a,τ ( j-1) [m(θ, τ)] k = A k (τ)E θ|y,a,τ ( j-1) e iθ k = A k (τ)e i θ k e -1 2 σ 2 θk
where the estimate of the mean θ ( j) k and the variance σ 2 θk ( j) can be obtained from a Kalman-type solution as presented in Section 4. The remaining term can be written as

E θ|y,a,τ ( j-1) 1 σ 2 n m(θ, τ) H m(θ, τ) = 1 σ 2 n k |A k (τ)| 2
So the cost function to be maximized is

L ( j) (τ) = -y H y -y H diag(d(τ))t( θ ( j) ) -t H ( θ ( j) )diag(d(τ) * )y + k |A k (τ)| 2        = -y -diag(d(τ))t( θ ( j) ) 2 , ( 48
)
where

[d(τ)] k = A k (τ) and [t( θ)] k = e i θ ( j) k e -1 2 σ 2 ( j) θk
. Finally, the proposed method is based on the following optimization

τ( j) = arg max τ        2 y H diag(d(τ))t( θ ( j) ) - k |A k (τ)| 2        . (49)
The iterative block method is sketched in Algorithm 1 and a scheme is given in figure 1. for j=p to p+N-1 do

) ( ˆj τ τ τ τ ) 1 ( ˆ-j τ τ τ τ ) ( ˆj θ ) ( 2 ˆj k θ θ θ θ σ σ σ σ z -1

4:

Estimation of θ j and δ using the Extended Kalman filter in eq.(

Time-delay estimation τ(i) from eq.(49). 7: end for To summarize: we have an iterative approximation of the ML solution for the estimation of the time-delay τ (block method) where we embed a Kalman-type solution to obtain the carrier phase and Doppler shift estimates (operating sample-wise).

This method resembles a classic DA ML method for the estimation of the time-delay after previous estimation and compensation of the carrier phase and Doppler shift. Indeed the cost function in (48) can be regarded as kind of distance between the observations and the expected noiseless signal after a carrier phase correction.

The advantage of our method is that it works iteratively to estimate the carrier phase/Doppler shift and the time-delay. Each operation can then be used to enhance the estimation performance of the other one giving rise to an iterative technique. Actually this intuitive idea is the basis of all the turbo methods for joint parameters estimation. The EM-type formulation used here is a way to arrive to the iterative technique more rigorously [START_REF] Al-Naffouri | An EM-based forward-backward kalman filter for the estimation of time-variant channels in OFDM[END_REF].

Computational complexity

The purpose of this section is to determine the implementation complexity in terms of the number of multiplications needed for our algorithm. The iterative method proposed is composed of two stages: carrier phase and Doppler shift estimation using an extended Kalman filter and time-delay estimation with the maximization of a cost function.

Concerning the first stage, the complexity of the time update step of the Kalman filter (state prediction and covariance of the prediction error) is (n 2

x + 2n 3 x ) and the complexity of the measurement update step (Kalman gain, state update and covariance of the measurement update error) is (2n x + 5n 2

x ), where n x = dim(x k ) is the state dimension (n x = 4, 5 and 7 for S = 1, 2 and 4, respectively, see eq.( 8)). So the Kalman stage complexity is O(Nn 3

x ), where N is the block size. The computational complexity of the optimization stage (eq.( 49)) is (3N), therefore, the bottleneck of the algorithm is the first stage and the overall computational complexity of our algorithm is at each iteration of the EM algorithm O(Nn 3

x ).

Computer Simulations

In this section we show the behavior of the proposed method by considering different scenarios. To assess the method's performance we assume the transmission over an AWGN channel of a M-sequence of length 511 bits generated using a Linear Feedback Shift Register (LFSR) with characteristic polynomial [1021] 8 (octal representation). We consider three oversampling factors (S = 1, 2 and 4) and a BOC shaping pulse (see figure 2). We fix the drift to δ = 0.2.

First we present the performance obtained with the Extended Kalman Filter for joint carrier phase and frequency shift estimation considering two scenarios:

• perfect time-delay synchronization,

• estimated time-delay.

In this case, we plot the Root Mean Square Error (RMSE), obtained over 250 independent Monte Carlo runs, versus the Signal to Noise Ratio (SNR). The SNR corresponds to the Carrier to Noise Ratio ( C N ) at the input of the receiver. In our case, as shaping pulse and . The performances obtained are compared with the HCRB. We compute the bound and the RMSE for the T -spaced symbol reference points for S = 1, 2, 4 (see [START_REF] Vilà-Valls | Oversampled phase tracking in digital communications with large excess bandwidth[END_REF] for the comparison of the estimation using T -spaced symbol reference points and T -spaced symbol mid-points and more scenarios). We note that we compared our method with the standard carrier synchronization Fitz's method [START_REF] Fitz | Planar filtered techniques for burst mode carrier synchronization[END_REF] but this completely fails in our scenarios, which is clear because this method assumes a deterministic phase distortion.

In figure 3 we analyze the HCRB and the EKF behavior for a fixed SNR versus phase-noise variance. We present a scenario with a really low SNR (as used in satellite based positioning), SNR = -20 dB. We can see that using 1 point per symbol (S = 1) the performances on the estimation are far from the theoretical bound. This is because the CRB does not give a good lower bound in the large error regions. In counterpart, with S = 2 and S = 4 we can measure the gain given by the oversampling and the good performances of the algorithm. The gain obtained with the oversampling is greater at small σ 2 w . We also plot on the same figure the performances obtained with an error on the timedelay estimation (τ -τ = T c /8), We can see that the performances are worse than those obtained with a perfect time-delay synchronization but still acceptable with S = 2 and S = 4.

Figure 4 superimposes, versus the SNR, the on-line HCRB and the RMSE obtained with the EKF. We have a phase with a moderate variation, σ 2 w = 0.001 rad 2 . We assume three different scenarios: a perfect time- w . For the estimated (error) time-delay scenarios we note that the performances obtained are a little bit worse than with the perfect synchronization but are still acceptable.

Figure 5 superimposes versus the SNR, the RMSE obtained with the EKF and the on-line HCRB for the Doppler shift estimation. We consider a phase noise variance σ 2 w = 0.001, three oversampling factors S = 1, 2, 4 and a fixed block size N = 511 (we note that the performance on the estimation of the drift depends on the block size N because the parameter to be estimated has a constant value). We note that the performances increase when increasing S and the convergence to the lower bound depends on the error on the time-delay esitmation and the oversampling factor. But in all the cases we obtain a good estimation of the Doppler shift.

After the analysis of the carrier phase and Doppler 10 The convergence and performance of the algorithm mainly depend on the process and measurement noises affecting the system and on the estimation of the carrier phase provided by the first stage. We first consider the joint estimation in a scenario with a SNR = 0 dB where the algorithm almost always converges to the true value (see the table at the end of the section). In figure 6(a), we plot the function to be optimized L(τ) for different iterations of a single realization of the algorithm (with a really fast varying phase offset with variance σ 2 w = 0.1) and in figure 6(b) we show the convergence of the algorithm for different carrier phase evolutions (σ 2 w = 0.001, 0.01 and 0.1). In both cases the true time-delay is τ = T c /8 and the algorithm is initialized at τ init = 3T c /8.

We note that the speed of convergence depends on the phase and observation noises, and that stronger the noises lower the speed of convergence. Indirectly we can see that the good performance of the algorithm directly depends on the carrier estimation stage performance. We also note that in this scenario (where the measurement noise is not really strong) we obtain a good time-delay estimation with few iteration, so the performances obtained with the proposed method are encouraging. But we have to analyze the behavior of the algorithm for lower values of SNR (down to -20dB) and the error associated on the estimation procedure.

On the following we plot the empirical probability density function (e-pdf) of the error on the time-delay estimation for different scenarios, which is computed over 1000 independent Monte Carlo runs. As a reference we fit a Gaussian distribution to the same sam- ple data sets. In figure 7, we plot the e-pdf and the fitted Gaussian distribution for a phase evolution with σ 2 w = 0.001 and a SNR=-10 dB. In figure 8, we consider the same scenario but with a really low SNR=-20 dB. In both cases we use a block size N = 511, the true time-delay is τ = T c /4 and the algorithm is initialized at τ init = T c /2.

We can see that the Gaussian distribution fits correctly with the e-pdf of the time-delay estimation error. On the following table we show the empirical mean and standard deviation computed from the error sample set considering different signal to noise ratios and block sizes: We note that the mean is always equal to the true timedelay but the standard deviation is three times greater in the case SNR= -20dB compared with the SNR=-10 dB case. Even in these conditions (really low SNR scenarios) the performances obtained are still acceptable. The error on the other two cases, SNR=0 dB and SNR= 10 dB, are really weak and the algorithm converges to the true value almost always. We also give the sample mean and standard deviation with different block sizes (for SNR=0 dB). We can see that smaller the block size greater the variance of the estimation error.

Conclusion

In usual transmission systems, the roll-off is between 0 and 100%, however, in the context of satellite positioning systems, like GPS and GALILEO, time limited shaping pulse are used and the Nyquist-Shannon sampling theorem does not apply. These special conditions let us hope a significant receiver synchronization performance improvement when the received signal is oversampled (using more than one sample per symbol).

In this paper, we study the gain due to an oversampling of the received signal for the problem of dynamical carrier phase tracking, and we propose a method for joint carrier and time-delay synchronization. Assuming that the data are known at the receiver, we derive the Hybrid Cramér-Rao Bound for carrier estimation, and we couple a Kalman-based DA algorithm and an EMtype method for joint carrier and time-delay estimation in such an oversampled scenario.

This study shows several improvements when a fractionally-spaced method for phase and frequency offset estimation is used. The estimation MSE decreases as the oversampling factor S increases and the interest of oversampling is more important at low SNR. For S = 1 or 2 samples per symbol, the results obtained with the EKF are close to the theoretical bound for slow and moderate phase evolutions. For S = 4, the HCRB is lower than for S = 2 but the EKF performance does not show the same improvement. We also note the limitations of the algorithm when having an extremely rapidly varying phase evolution with respect to the symbol interval.

We have shown the good performance of the iterative block method proposed for joint time-delay and carrier estimation for high and moderate SNR, and an acceptable estimation performance in low SNR scenarios. This method is based on a two-step iterative method, including the oversampled carrier synchronization solution based on the EKF. The convergence of the method directly depends on the phase and Doppler distortion, and on the estimation of the carrier phase provided by the EKF stage. For a slowly varying carrier phase the algorithm convergence is really fast, and as the phase offset gets stronger the algorithm needs more iterations to converge to the good time-delay value.

  42)In a DA context, the likelihood function w.r.t τ is p(y|a, θ; τ) = N (y; m(θ, τ), Γ) ,

Figure 1 :

 1 Figure 1: Block scheme of the joint estimation method

Figure 2 :

 2 Figure 2: BOC shaping function Π(t) and its autocorrelation g(t)

Figure 3 : 1 σ 2 n

 312 Figure 3: HCRB and RMSE obtained with the EKF, for the carrier phase estimation, versus the phase noise variance for three different oversampling factors S = 1, 2 and 4, and a fixed low SNR = -20 dB. We consider two scenarios: perfecte time-delay synchronization and an error on the delay estimation τ -τ = T c /8.

Figure 4 :

 4 Figure 4: HCRB and RMSE obtained with the EKF, for the carrier phase estimation, versus the SNR for three different oversampling factors S = 1, 2 and 4, and a fixed phase noise variance σ 2 w = 0.001. We consider three scenarios: perfecte time-delay synchronization, an error on the delay estimation τ -τ = T c /8 and an error τ -τ = T c /4.

Figure 5 :

 5 Figure 5: HCRB and RMSE obtained with the EKF, for the estimation of the Doppler shift, versus the SNR for three different oversampling factors S = 1, 2 and 4, and a fixed phase noise variance σ 2 w = 0.001. We consider three scenarios: perfecte time-delay synchronization, an error on the delay estimation τ -τ = T c /8 and an error τ -τ = T c /4.

Figure 6 :

 6 Figure 6: SNR= 0 dB, τ = T c /8 and τ init = 3T c /8. (a) Values of the function to be optimized L(τ) in different iterations for a fast varying phase evolution, σ 2 w = 0.1. (b) Convergence of the algorithm for three phase noise variances: σ 2 w = 0.001, σ 2 w = 0.01 and σ 2 w = 0.1.

Figure 7 :

 7 Figure 7: e-pdf and Gaussian fitted distribution (dotted line) for the time-delay estimation error. SNR=-10 dB and σ 2 w = 0.001.

Figure 8 :

 8 Figure 8: e-pdf and Gaussian fitted distribution (dotted line) for the time-delay estimation error. SNR=-20 dB and σ 2 w = 0.001.

Acknowledgement This work has been partially supported by project 2009 SGR 1236 of the Catalan Administration (AGAUR), by LURGA project of the French ANR (Agence Nationale de la Recherche), by the Spanish Science and Technology Commission: TEC 2010-19171 (MOSAIC) and by the European Commission in the framework of the COST Action IC0803 (RFCSET).


Fouad Demami

Régis Rogel

Anne-Claire Salaün

Laurent Pichon

A-C Salaün

Electrical properties of polysilicon nanowires for devices applications

Keywords: Polysilicon nanowires, spacer, in-situ doping, resistor, TFT)

Introduction

Owing to their physical and electrical properties, SiNWs (silicon nanowires) represent a promising material with strong potential for a large variety of applications. They are currently attracting much attention as components for future nanoelectronic devices such as nanowire field effect transistors [1], photonic and optoelectronic devices [2], and as chemical or biological sensors [3,4]. SiNWs present large advantages with respect to bulk silicon: high surface to volume ratio, possibility of surface functionalization, synthesis compatible with large area technology, leading to the development of innovative electronic devices. SiNWs development could significantly impact areas of electronics, genomics, biomedical diagnosis, drug discovery.... They can be prepared by the top-down approach, using various advanced methods such as e-beam [5], AFM [6] or deep UV [7] lithography. The main disadvantage of these advanced lithographic tools with nanometer size resolution rests on their high cost. The bottom-up approach usually consists in a metal catalytic growth [8]. Thus, VLS (Vapour-Liquid-Solid) growth technique uses metallic nanoparticles as aluminium, nickel, gold… This approach suffers from both metal contamination and the difficulty in precisely positioning the device location.

In this work, SiNWs are synthesized using a low cost technique commonly used in microelectronic industry: the sidewall spacer formation technique. Previous works re-ported the feasibility of undoped silicon NWs and their potential applications for fabrication of field effect transistors [9] as biological sensors [10]. In this case, SiNWs present a polycrystalline structure, and their electrical behaviour is expected to be different from crystalline SiNWs due to higher defect density and carriers confinement at nanometric scale. In addition, thanks to their higher surface to volume ratio compared to thin films, these polysilicon NWs are promising for charged (bio)chemical species detection.

Polysilicon NWs based devices performances are strongly correlated with NWs electrical properties. However, polysilicon NWs electrical behaviour has not been yet studied. Thus, prior to develop polysilicon NWs based devices for high efficient (bio)chemical sensing applications, a well understanding of electronic polysilicon NWs properties is required. Here, we present a study of the polysilicon NWs electrical behaviour as a function of the N-type doping level, and their potential applications in devices (resistors, TFTs) fabrication are demonstrated.

Experimental details

The key nanowire fabrication steps are illustrated in figure 1 (a). At first, a dielectric film A is deposited and patterned into islands by conventional UV lithography. Then, a polysilicon layer is deposited by LPCVD (Low Pressure Chemical Vapour Deposition) technique. Accurate control of the polysilicon layer Reactive Ion Etching (RIE) rate leads to the formation of nanometric size sidewall spacers that can be used as nanowires. The feasibility of this technological step was previously demonstrated [11].

Phosphorus doped polysilicon films used for nanowires are deposited by thermal decomposition of a mixture of pure silane (SiH 4 ) and phosphine (PH 3 ) diluted in helium. The in-situ doping level is controlled by adjusting the PH 3 /SiH 4 mole ratio varying from 0 (undoped film) to 1.2×10 -3 (heavily doped film). The corresponding incorporated phosphorus atoms concentration C p , previously determined from SIMS (Secondary Ions Mass Spectroscopy) analysis, varies from 2×10 16 cm -3 to 2×10 20 cm -3 [12]. Silicon films are deposited in an amorphous state at 550°C and 90 Pa, and then crystallized by a thermal annealing in vacuum at 600°C during 12 hours.

The in-situ doped polysilicon NWs are integrated into resistors (fig. 1 (b)) for electrical characterization. The nanowires are capped by a SiO 2 layer deposited by Atmospheric Pressure Chemical Vapour Deposition (APCVD) technique at 390°C and contact openings are wet etched.

Aluminium is then thermally evaporated and electrodes defined by wet etching.

Static electrical characteristics I(V) are collected at room temperature using a HP 4155 B semiconductor parameter analyzer. For current-temperature measurements, in the range 200K to 530K, samples are placed in a cryostat in vacuum (10 -3 -10 -4 Pa) and dark current is measured using a Keithley 617 electrometer.

Results and discussion

Electrical properties of 10 µm length and 100 nm curvature radius polysilicon NWs are investigated in function of the in-situ doping level. Such polysilicon NWs were previously investigated by scanning electronic microscopy [11].

I(V) characteristics of the polysilicon NWs based resistors are plotted in the figure 2 for 2×10 16 Increasing the doping level from about 1× 10 18 cm -3 results in an abrupt resistivity drop of about 6 orders of magnitude for only a factor of 10 further increase in doping concentration. Beyond that range resistivity decreases following the same doping effect as in crystalline silicon. The shape of the curve, quite similar to previously observed results for in-situ doped polysilicon layers [12][13], indicates that polysilicon NWs doping effect is in accordance with the Seto's theory [14]. In this model, the conduction is thermally activated. Figure 4 (a) presents the temperature dependence of the electrical conductivity, σ, for different C p values. The activation energy, E a , is deduced from the slope of the decreasing linear curve at high temperatures. The variations of E a with C p plotted in figure 4 (b) highlight an evolution in accordance with the Seto's theory. Indeed, for low doping levels, it means below a critical concentration, N*, E a slightly decreases with C p . In this case, generation of trapped carriers dominates. So, the corresponding activation energy is related to the energy barrier needed for emission of carriers into the conduction band from deep levels located into the bad gap associated with residual and/or structural defects. For high doping levels (>N*), free carriers diffusion over the intergranular barrier E b (induced by trapped carriers at grain boundaries) dominates. In this case, E a = E b and it decreases with the doping level. In this model, N* is related to the quality of the polysilicon layer both through the defects density at the grain boundaries and the lateral grain size. In our case N* ~ 10 18 cm -3 , corresponding to the value obtained for previous studies [12][13][14]. These results show that this electrical conduction model can be applied for polysilicon NWs at high temperatures.

For lower temperatures, the slight dependence of σ with 1000/T is commonly described by the conduction in localized states [15]. These states are associated with a high defect density uniformly distributed along the NWs. Two types of defects can be involved either located i) at the surface of the wire and/or ii) in the core of the polysilicon layer. Indeed, in this last case, because the amorphous silicon crystallization process begins at the Si 3 N 4 /amorphous silicon interface, the defects density (including grain boundaries) is higher in the inferior part of the polysilicon layer (few nanometers thick) constituting the nanowire. For the lowest temperatures range σ does not depend on temperature, indicating tunnel conduction through these defects.

Our results highlight the good control of the doping level on the 100 mn long curvature radius polysilicon NWs over a wide range, with a doping effect similar to polysilicon layer [12][13][14]. This study shows that despite the high defects density within these polysilicon NWs, there is no significant doping species segregation at grain boundaries affecting doping efficiency. In addition, thanks to their high surface to volume ratio, such doped NWs are promising candidates for optimized devices devoted to charged (biological and/or chemical) species detection.

Thin Film Transistors (TFTs) are fabricated (see figure 5 (a)) using undoped polysilicon NWs following the synthesis method described in section II. In our case, the substrate is heavily N-type doped crystalline silicon wafer acting as gate electrode. The film A is a 100 nm thick Si 3 N 4 layer deposited by LPCVD technique at 600°C. Source and drain regions are made of heavily in-situ N-type doped polysilicon layer deposited by LPCVD technique. This layer is patterned by reactive ion etching with SF 6 plasma before thermal evaporation of aluminium used for source and drain electrodes. sistor) [START_REF] Sze | Physics of semenoductor Devices[END_REF]. µ is deduced from the maximum slope of the I DS (V GS ) curve (transconductance) plotted in the saturation mode and V T is determined by the intercept of the I DS 1/2 (V GS ) curve with the gate axis voltage. The low field effect mobility (~ 3 cm 2 /Vs) and the high threshold voltage (V T ≈ 8 V) values are related to the poor gate insulator (Si 3 N 4 )/active layer (polysilicon NWs) interface quality. The electrical properties of our polysilicon NWs based TFTs, as well as technology, have to be improved for high efficient chemical and/or biological sensing applications.

Conclusion

The sidewall spacer formation technique allows the feasibility of arrays of parallel horizontal polysilicon NWs using conventional UV lithography. This low cost and compatible with silicon thin film technology patterning method represents an attractive process for polysilicon nanowires fabrication. In-situ doping control of these 100 nm curvature radius polysilicon NWs is demonstrated, and thus it makes them good candidates for the fabrication of electrically controlled thin film devices, particularly for integrated chemical and/or biological sensors.

The next step of this study will consist in the fabrication of lower curvature radius polysilicon NWs (< 50 nm) and to study their specific electrical properties at nanometric scale and their ability for fabrication of optimized high sensitive and low detection level chemical and/or biological sensors. V DS (V)

Figure 1

 1 Figure 1 (a) Fabrication of polysilicon NWs by the sidewall spacer formation technique, (b) polysilicon NWs based resistors.

  Figure 2 I(V) characteristics of the polysilicon NWs based resistors at different phosphorus atoms concentration levels (10 µm length NWs).

Figure 3

 3 Figure 3 Variations of the resistivity of the polysilicon NWs versus the phosphorus atoms concentration.

Figure 4

 4 Figure 4 (a) Arrhenius diagrams of polysilicon NWs conductivity for various phosphorus atoms concentrations. (b) Variations of the corresponding activation energy of the conductivity versus the phosphorus atoms concentration.

  Field effect behaviour displayed on the figure 5 (b) shows a switching ratio ~ 10 4 and good ohmic contacts (see inset of fig. 5 (b)). Threshold voltage, V T , and optimum field effect mobility, µ, are determined according to the classical conduction electrical model used for the MOSFET (Metal Oxide Semiconductor Field Effect Tran

Figure 5

 5 Figure 5 (a) Schematic illustration of the polysilicon NWs based TFT. (b) Transfer characteristics (W/L = 0.1µm/10µm, V DS = 4 V). Inset: output characteristics (V DS = 1 V to 4V, step 1 V).

  Ω.cm for the lightly doped polysilicon NWs.

	1	reaches 10 5
	2
	3
	4
	5
	6
	7
	8
	9
	10
	11
	12
	13
	14
	15
	16
	17
	18
	19
	20
	21
	22
	23
	24
	25
	26
	27
	28
	29
	30
	31
	32
	33
	34
	35
	36
	37
	38
	39
	40
	41
	42
	43
	44
	45
	46
	47
	48
	49
	50
	51
	52
	53
	54
	55
	56
	57

Acknowledgements This work was supported by the Agence Nationale de la Recherche with project Sena (ANR-09-JCJC-0072-01).