Sara Teixeira

Marie-Anne Cambon-Bonavita

Ester A Serrão

Daniel Desbruyeres

Sophie Arnaud-Haond

Daniel Desbruyéres

Recent population expansion and connectivity in the hydrothermal shrimp Rimicaris exoculata along the Mid Atlantic Ridge

Keywords: Bottleneck, deep sea biogeography, dispersal, genetic diversity, hydrothermal vents, marine biogeography, Mid-Atlantic Ridge, Rimicaris exoculata, spatio-temporal distribution

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Dispersal plays a central role in both the dynamics and evolution of species. Effective dispersal (i.e., gene flow) attenuates genetic differences among populations and increases the genetic variability within populations [START_REF] Billiard | Evolution of migration under kin selection and local adaptation[END_REF]. Dispersal is a risky strategy as immigrants may never reach a suitable habitat, or face higher mortality and lower reproductive success in the receiver population [START_REF] Ims | Spatial synchronization of vole population dynamics by predatory birds[END_REF][START_REF] Garant | The multifarious effects of dispersal and gene flow on contemporary adaptation[END_REF]. Yet, dispersal can also offer clear advantages, particularly in fragmented ephemeral habitats, as it allows areas to be explored for more favourable habitats and extinction to be avoided by escaping locally deteriorating physical or biotic conditions, such as competition [START_REF] Hamilton | Dispersal in stable habitats[END_REF][START_REF] Olivieri | Metapopulation genetics and the evolution of dispersal[END_REF]. By enabling populations to randomly spread their progeny more evenly among different sites dispersal may indeed represent a favourable bet-hedging strategy in temporally variable or ephemeral ecosystems [START_REF] Ronce | Life history evolution in metapopulations. Ecology, genetics, and evolution of metapopulations[END_REF].

Deep-sea hydrothermal vents are an extreme example of a spatially and temporally unstable environment. Environments of this type are scattered worldwide along the global mid-ocean ridge system, and are sometimes thousands of kilometres apart [START_REF] Tunnicliffe | Influence of sea-floor spreading on the global hydrothermal vent fauna[END_REF]Van Dover, 2000). They are extreme and unique habitats, with communities composed of specialized animals, dependent on chemoautotrophic microorganisms that exploit the reduced compounds (mainly H 2 S and CH 4 ) abundant in vent fluids [START_REF] Desbruyères | A review of the distribution of hydrothermal vent communities along the northern Mid-Atlantic Ridge: dispersal vs. environmental controls[END_REF]. Each hydrothermal vent has a characteristic profile that may differ in depth, temperature and/or chemical composition (Van Dover, 2000), adding abiotic heterogeneity to the spatial component of habitat fragmentation that potentially influences the patchy community distribution and limited effective dispersal. Besides the distance separating successive hydrothermal vents, communities may also become genetically isolated across transform faults, which act as potential barriers between populations and limit effective dispersal among sites [START_REF] Plouviez | Comparative phylogeography among hydrothermal vent species along the East Pacific Rise reveals vicariant processes and population expansion in the South[END_REF]. This hypothesis is supported by the high levels of endemism found across biogeographic provinces (Van Dover et al., 2002;[START_REF] Fabri | Combining the Biocean database for deep-sea benthic data with the online Ocean Biogeographic Information System[END_REF][START_REF] Bachraty | Biogeographic relationships among deep-sea hydrothermal vent faunas at global scale[END_REF]Moalic et al. submitted) as well as by the rather high genetic differentiation found across some neighbouring sites (e.g. [START_REF] Jollivet | Hydrothermal-vent alvinellid polychaete dispersal in the eastern Pacific. 2. A metapopulation model based on habitat shifts[END_REF]Won et al., 2003;[START_REF] Hurtado | Distinct patterns of genetic differentiation among annelids of eastern Pacific hydrothermal vents[END_REF][START_REF] Matabos | Geographic clines and stepping-stone patterns detected along the East Pacific Rise in the vetigastropod Lepetodrilus elevatus reflect species crypticism[END_REF]Young et al., 2008;[START_REF] Plouviez | Comparative phylogeography among hydrothermal vent species along the East Pacific Rise reveals vicariant processes and population expansion in the South[END_REF], which indeed suggest a limited effective dispersal over large distances. On the other hand, depending on the sea floor geology and magmatic and tectonic events, individual hydrothermal vents have a relatively short lifespan of 1 to 100 years for fast spreading ridges such as the East Pacific Rise [START_REF] Hessler | Temporal change in megafauna at the rose garden hydrothermal vent (Galapagos Rift; eastern tropical Pacific)[END_REF][START_REF] Lalou | New age data for Mid-Atlantic Ridge hydrothermal sites: TAG and Snakepit chronology revisited[END_REF], and up to 300 years in more stable, slow spreading ridges such as the Mid-Atlantic Ridge [START_REF] Brazelton | Archaea and bacteria with surprising microdiversity show shifts in dominance over 1,000-year time scales in hydrothermal chimneys[END_REF]. These common temporal fluctuations in vent activity and changes in the geochemical profile of their habitat prevent long-term local population persistence [START_REF] Tunnicliffe | The biology of hydrothermal vents: ecology and evolution[END_REF][START_REF] Jollivet | Specific and genetic diversity at deep-sea hydrothermal vents: an overview[END_REF].

Those species that escape local extinction due to habitat disappearance and which are capable of (re)colonizing newly active sites may also have the capacity for long-distance effective migration.

Species associated with hydrothermal vents should therefore also have good colonization abilities with high rates of dispersal, growth and fecundity (Vrijenhoek, 1997;Tyler & Young, 1999). Indeed, in 1991 the formation of new hydrothermal vents in the East Pacific Rise was studied and, astonishingly, microbe-derived material was found in active areas of the vent immediately following their formation. Within the first year, mobile vent and non-vent fauna such as amphipods, copepods, crabs and non-vent shrimps then proliferated (Shank et al., 1998a). These vastly differing sources of influence lead to contradictory expectations about strategy, which have yet to be tested. On one hand, low effective dispersal and high self recruitment may be favoured across limited distances due to the patchy nature of vent distribution and abiotic characteristics. On the other hand, larger scale dispersal potential may be selected for in the long term, allowing species persistence despite instability of vent activity. Finally, temporal instability [START_REF] Desbruyères | Biology and ecology of the "Pompeii worm" (Alvinella pompejana Desbruyères and Laubier), a normal dweller of an extreme deep-sea environment: A synthesis of current knowledge and recent developments[END_REF], although apparently less pronounced in the Mid-Atlantic Ridge (MAR) than in the East Pacific Rise (EPR) [START_REF] Copley | Assessment of decadal-scale ecological change at a deep Mid-Atlantic hydrothermal vent and reproductive time-series in the shrimp Rimicaris exoculata[END_REF] on one site), may also result in short-term changes in genetic composition.

Along the Mid-Atlantic Ridge (MAR), vent communities are dominated by shrimps, mussels, anemones and crabs [START_REF] Tunnicliffe | The biology of hydrothermal vents: ecology and evolution[END_REF]Van Dover, 1995). The first described and the most abundant species is the shrimp Rimicaris exoculata Williams & Rona, 1986. This species is often present in high-density swarms (> 1500 individuals/m 2 ) around high-temperature sulphide chimneys at deeper hydrothermal vent sites [START_REF] Segonzac | Enigma of the trophic adaptation of the shrimp Alvinocarididae in hydrothermal areas along the Mid-Atlantic Ridge[END_REF][START_REF] Copley | Spatial and interannual variation in the faunal distribution at Broken Spur vent field (29ºN, Mid-Atlantic Ridge)[END_REF]. The species has highly modified eyes that form a single dorsal organ, which has been suggested to detect black-body radiation associated with high temperature vents (Van Dover et al., 1989;Van Dover, 1995). This shrimp does not seem to be a scavenger or a predator because its gut is filled with minerals, mainly iron oxides and sulphurs, and its nutrition is hypothesized to be derived from several epibiotic bacterial sources, mainly those located in the gill chambers and gut (see Durand et al., 2010 and references therein). Although hydrothermal vents at MAR have been known since 1985 [START_REF] Rona | Black smokers, massive sulphides and vent biota at the Mid-Atlantic Ridge[END_REF], little is known about population genetic structure of its most abundant shrimp species. Genetic studies conducted to date have focused on clarifying the systematic status of the juvenile morphotype of R. exoculata [START_REF] Creasey | Genetic comparison of two populations of the deepsea vent shrimp Rimicaris exoculata (Decapoda: Bresiliidae) from the Mid-Atlantic Ridge[END_REF]Shank et al., 1998b). The allozyme analyses conducted in these studies showed very low genetic diversity within sites and essentially no genetic differentiation between populations, suggesting that there is high gene flow.

The main goal of this study was to infer the spatial and temporal distribution of genetic diversity of R. exoculata from different hydrothermal vents along the Mid-Atlantic Ridge, in order to assess the extent of effective population sizes, rates of dispersal and recent demographic histories at these sites.

We used mitochondrial cytochrome c oxidase subunit I (COI) as a molecular marker to estimate the level of genetic variability in natural populations, to test for genetic subdivision and to obtain insights into population history, dynamics and demography at the ridge level. In order to assess the effect of potential temporal instability of hydrothermal activity and population dynamics of associated species, two of the five populations studied were analysed at time intervals of one or ten years.

Materials and Methods

Sampling and DNA extraction

Samples were collected from five locations corresponding to different hydrothermal vent fields (Rainbow, Logatchev, TAG, Ashadzé and South MAR) along the Mid-Atlantic Ridge (Fig. 1, Table 1). For two locations (Rainbow and Logatchev), samples from different years were used to assess temporal variability; a time interval of 10 years (Rainbow 1997 and2007) and 3 consecutive one year intervals (Logatchev 2005(Logatchev , 2006(Logatchev , 2007) ) were used. Samples were collected using a slurp-gun from the ROV (Remotely Operated Vehicle) or the Nautile. Prior to each dive, the bowls used for collecting the shrimps were aseptically washed with ethanol (96 %) before being filled with sterile seawater. Once on board, live specimens were either entirely frozen or immediately dissected into body components under sterile conditions and frozen or stored in 70 % alcohol. DNA extraction was performed using the CTAB (cetyl trimethyl ammonium bromide) method [START_REF] Doyle | Isolation of plant DNA from fresh tissue[END_REF]) on muscle tissue. Sampling size could not be pre-defined to a constant value across all sampling sites and times due to the overall difficulty in obtaining samples from deep hydrothermal vents; in some sites sampling size was thus limited to the shrimp samples available from Ifremer research cruises for each particular location and time, and (for the 4 ºS location) some samples kindly provided by Max Plank Institute (N. Dubilier).

Polymerase chain reaction and sequencing

Mitochondrial COI was amplified using the universal primers LCO1490 and HCO2198 described by [START_REF] Folmer | Conserved primers for PCR amplification of mitochondrial DNA from different invertebrate phyla[END_REF]. All polymerase chain reaction (PCR) amplifications were carried out in 50 µL volumes containing 50 ng DNA, 1× reaction buffer (GoTaq, Promega), 0.2 mM of each dNTP, 2.5 mM MgCl 2 , 0.4 U Taq DNA polymerase (GoTaq, Promega, Madison, WI, USA) and 0.6 µM of each primer. The PCR amplification was conducted in a Perkin-Elmer Gene Amp System 7200 (Waltham, MA, USA) with the following program: 2 min. at 95 ºC; 35 cycles composed of 1 min denaturation at 95 ºC, 1 min at 52 ºC and 1.5 min elongation at 72 ºC, followed by a final elongation step of 7 min at 72 ºC. PCR products were sent to be purified and sequenced commercially at Macrogen, Inc.

(Seoul, Korea). Sequences were proofread and aligned using Geneious software (Drummond et al., 2009). The 18S rDNA (~1.7 kb) was amplified for 4 individuals of all sampled locations and years (except for the Ashadzé population, where only the three individuals available were amplified) using 18S-82F/18S-1498R primers [START_REF] López-García | Autochthonous eukaryotic diversity in hydrothermal sediment and experimental microcolonizers at the Mid-Atlantic Ridge[END_REF] and sequenced in an ABI 3130 XL automatic sequencer (Applied Biosystems, Foster City, CA, USA), to provide hints into nuclear variation. The nuclear regions ITS 1 and 2 (~931 bp) were also amplified for four individuals of three populations (Rainbow, Logatchev and South MAR) using the primer pair ITS 1 and 4 [START_REF] Kumar | Use of PCR targeting of internal transcribed spacer regions and single-stranded conformation polymorphism analysis of sequence variation in different regions of rRNA genes in fungi for rapid diagnosis of mycotic keratitis[END_REF]. The obtained sequences were cloned (PGem-T easy vector, Promega) following the supplier"s protocol; four clones per individual were sequenced using the ABI 3130 XL automatic sequencer (Applied Biosystems).

No diversity was detected across the 4 individuals from the 5 sites sequenced for 1200 bp of 18S rDNA. The amplified regions ITS1 and 2 proved to be too variable to be used: all the sequences obtained were different, including sequences from different clones of the same individual.

Data analysis

ARLEQUIN version 3 [START_REF] Schneider | ): Software for population genetics data analysis[END_REF] was used to estimate gene diversity and conduct statistical tests. For each population, the following statistics were computed: haplotype (h) [START_REF] Nei | Molecular evolutionary genetics[END_REF] and nucleotide diversities (π 2 ) [START_REF] Nei | Molecular evolutionary genetics[END_REF], and mean number of pairwise differences (π 1 ) [START_REF] Tajima | Evolutionary relationships of DNA sequences in finite populations[END_REF]. For populations with significantly different sample sizes, a re-sampling method was applied using R to generate files for 100 random re-samplings of the population haplotype frequencies to a minimum common size. The 100 generated files were then analysed with the batch file function of ARLEQUIN software to re-estimate gene diversity in these sub-samples. To assess population differentiation, pairwise F ST values were calculated following the method of [START_REF] Hudson | Estimation of levels of gene flow from DNAsequence data[END_REF] and exact tests of differentiation were conducted following the method of [START_REF] Raymond | Genepop (version 3.2): population genetics software for exact tests and ecumenicism[END_REF].

In order to analyse and illustrate the phylogenetic relationships among haplotypes, i) a medianjoining network was constructed using Network v. 4.1.0.9 [START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF] to infer the most parsimonious branch connections between the sampled haplotypes, and ii) a neighbour-joining phylogram was constructed using the mutation model HKY (Geneious software; Drummond et al., 2009) and the bresiliid sister species Alvinocaris sp. as an outgroup.

We tested whether the sampled populations fitted mutation drift equilibrium expectations, and assessed which population demographic model was most appropriate, i.e. the size-stationary model or the expansion model. We determined Fu"s F s [START_REF] Fu | New statistical tests of neutrality for DNA samples from a population[END_REF]) and Tajima"s D [START_REF] Tajima | Statistical-method for testing the neutral mutation hypothesis by DNA polymorphism[END_REF], which are sensitive to departures from both selective neutrality and from population size equilibrium caused by expansions or bottlenecks [START_REF] Tajima | The amount of DNA polymorphism maintained in a finite population when the neutral mutation rate varies among sites[END_REF][START_REF] Fu | Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection[END_REF]. Both statistics result in negative values after a population expansion [START_REF] Ray | Intra-deme molecular diversity in spatially expanding populations[END_REF] or a selective sweep, whereas positive values are expected under balancing selection of recent bottlenecks. ARLEQUIN was used to perform tests and estimate the associated P-values. Mismatch curves were drawn for each population and observed values were compared to the expected curves fitted for the constant population size model and for the population growth/decline model, using the DnaSP 4.0 program [START_REF] Rozas | DnaSP, DNA polymorphism analyses by the coalescent and other methods[END_REF].

Results

Genetic diversity and population differentiation

A total of 47 haplotypes were recovered from the 152 individuals analysed, 30 (63.8 %) of the haplotypes were "private", i.e. unique to a single locality (Table 2; Fig. 1; accession numbers: HM125910 to HM125956). The most common haplotype was present in all population samples and was the most abundant haplotype across all populations, identified in 49 % of all the sampled individuals across the entire study region. This main haplotype was also central to all other haplotypes, most of these represented by a single individual and divergent by a single point mutation, leading to this haplotype being the core of a star-like cluster topography (Fig. 2). The remaining haplotypes (36.2 %), although not unique, were mainly (27.7 %) represented by only one individual in each population. Accordingly, the neighbour-joining phylogram revealed extremely unresolved clusters supported by poor bootstrap values (see Appendix S1 in Supporting Information).

At a temporal scale of a 10-year sampling interval (Rainbow 1997 and2007), the unique haplotypes (present in one year and not on the other) represented 56 % of the total haplotypes found, whereas on the temporal scale of a one year sampling interval (Logatchev 2005(Logatchev , 2006(Logatchev , 2007) ) the proportion of unique haplotypes was higher (88.5 %).

Haplotype diversity (h) was high, ranging from 0 to 1 in populations with small sample sizes (Ashadzé: 3 individuals and 1 haplotype; South MAR: 4 individuals and 4 haplotypes) and 0.69 to 0.82 in populations with larger sample sizes (>11 individuals). Considering these rather high values of haplotype diversity, nucleotide diversity (π 2 ) was low, ranging from 0 to 0.0033 (Table 2). Highest diversities were recorded in the South MAR population, although the sample size was very low (4 individuals), followed by the Rainbow 2007 population. Populations sampled in more than one year presented different diversity indices across years. The Rainbow population sampled in 2007 showed more polymorphic sites, and higher haplotype and nucleotide diversity than when it was sampled in 1997, even though the number of samples was lower in 2007 (26 individual) than in 1997 (32 individual). A random re-sampling method was therefore not needed for Rainbow, as the sample with the fewest individuals already presented the higher diversity indices. Random resampling was, however, performed for Logatchev, where samples had more contrasting sizes (12 to 31 individuals) with no obvious diversity trend through time. In order to dissociate a putative biological effect from a sampling size effect, the random re-sampling (100 times) was therefore performed to a minimum common size of 12 individuals per sampling time, and no significant differences were found in any of the diversity indices over time (Table 3).

No significant population differentiation was detected between any of the sampling sites or years using pairwise F ST estimates (Table 4).

Historical demography

Demographic analyses suggested the occurrence of population expansions, as highly significant negative values for Tajima"s D and Fu"s F S were obtained for all populations (Table 2). Mismatch distribution curves were consistent with models of population expansion, displaying a unimodal distribution (Fig. 3). When comparing the expected and observed mismatch curves under constant population size versus the growth/decline model (Fig. 4), there is an excess of rare variants in the observed curves in relation to that expected under constant population size. Furthermore, the observed curves seem to better follow the expected curve for the growth/decline model.

Discussion

The analysis of mtDNA haplotypes from a total of 152 individuals from five discrete vent localities along the Mid-Atlantic Ridge showed an overall high genetic variability within the studied populations. Both the number of haplotypes per population (1 to 16) and gene diversity (0.69 to 0.82) are within the range of values found for other vent organisms studied with the same marker, including tubeworms (e.g. Riftia pachyptila, Alvinella pompejana), limpets (e.g. Lepetodrillus elevatus) and bivalves (e.g. Bathymodiolus azoricus, Calyptogena magnifica) (e.g. [START_REF] Hurtado | Distinct patterns of genetic differentiation among annelids of eastern Pacific hydrothermal vents[END_REF][START_REF] Johnson | Migration, isolation, and speciation of hydrothermal vent limpets (Gastropoda; Lepetodrilidae) across the Blanco transform fault[END_REF][START_REF] Plouviez | Comparative phylogeography among hydrothermal vent species along the East Pacific Rise reveals vicariant processes and population expansion in the South[END_REF]. However, these results contradict those obtained for Rimicaris exoculata with allozymes, where genetic variability was amongst the lowest reported for vent organisms [START_REF] Creasey | Genetic comparison of two populations of the deepsea vent shrimp Rimicaris exoculata (Decapoda: Bresiliidae) from the Mid-Atlantic Ridge[END_REF]. This may be due to the different properties of the distinct markers because allozymes may exhibit lower polymorphism than mtDNA, possibly associated with protein constraints. Mitochondrial DNA polymorphism is classically considered as less influenced by possible selective processes [START_REF] Ballard | The incomplete natural history of mitochondria[END_REF], although it was recently demonstrated that this may not be a systematic rule [START_REF] Galtier | Mitochondrial DNA as a marker of molecular diversity: a reappraisal[END_REF]. In any case, DNA sequences retain more polymorphisms than proteins, as polymorphic DNA sequences include synonymous mutations; higher variability can, therefore, be detected by DNA sequences than by allozymes. On the other hand, a differential imprint of demographic events on mitochondrial versus nuclear DNA may be evoked, although lower diversity is usually expected for mtDNA in relation to nuclear DNA, due to the difference in effective population size. Such discrepancies may also be affected by the more restricted sampling that was available for the allozyme study, which included only two sampling locations despite the higher overall number of samples analysed (448 individuals).

The rather unresolved phylogram, with low bootstrap values, shows a generally low divergence without a geographical imprint. Together with the non-significant F ST values obtained in this study, these results indicate a lack of geographic partitioning of the numerous but slightly divergent haplotypes among the sampled populations. This lack of genetic structure may be due to several nonmutually exclusive factors: high levels of effective dispersal (gene flow) between all populations of the Mid-Atlantic ridge, i.e., panmixia at the ridge level; large effective population sizes (N e ); and/or recent common population history. The life cycle of Rimicaris is poorly understood. Some authors have suggested that R. exoculata might feed in the upper water column during larval dispersal [START_REF] Dixon | Results of DNA analyses conducted on vent shrimp postlarvae collected above the Broken Spur vent field during the CD95 cruise, August 1995[END_REF], and others that the juvenile stages of these shrimp differ physiologically from the adults [START_REF] Pond | Wax-ester reserves facilitate dispersal of hydrothermal vent shrimps[END_REF][START_REF] Allen | Lipid profiles of hydrothermal vent shrimps[END_REF] in traits such as wax-ester reserves and possibly photosynthetically derived lipids. These findings led the authors to hypothesize that the species had a long planktotrophic phase with long-distance migration achieved via the surface, so that some of the barriers usually considered, (depth, seafloor topology and transform faults) would not constrain this species. Yet, such a hypothesis implies either a huge risk factor in the case of passive dispersal, or an extremely sophisticated mechanism of active and deterministic migration (that has not been identified so far) allowing a return to hydrothermal habitat at the end of larval development. In any case, whether dispersal occurs via the surface or not, the little information available to date on R. exoculata suggests a significant dispersal potential either at larval or at adult stage (although no observations of adults away from hydrothermal sites have been reported so far).

The low degree of genetic differentiation among populations may also be attributed to a relatively recent common history together with large effective population sizes, which limit genetic drift [START_REF] Crow | Basic concepts in population, quantitative, and evolutionary genetics[END_REF]. The scenario of (1) high levels of gene flow and that of (2) low gene flow since a recent divergence and large N e may both leave similar signatures of low genetic structure of populations. The hypothesis of large effective population sizes limiting drift is concordant with the population biology of this species, as R. exoculata is reported to have remarkably high population densities of over 1500 individuals / m 2 [START_REF] Copley | Spatial and interannual variation in the faunal distribution at Broken Spur vent field (29ºN, Mid-Atlantic Ridge)[END_REF]. Such large population sizes may, in the absence of temporal fluctuation or significant variance in reproductive success, result in large effective population sizes, limiting the differentiating effects of drift. Yet, the haplotype phylogram and network, together with the analyses of population history, indicate a recently very low effective population size as they show: (1) a star-like haplotype network where many new haplotypes are derived from one dominant haplotype via single mutation steps, indicating that they are the product of recent mutation events; (2) negative and highly significant values of Tajima"s D and Fu"s F S tests, as expected under population bottleneck/expansion scenarios; and (3) the expected and observed mismatch curves that follow a growth/decline model. This profile supports a scenario of relatively young populations sharing a recent common history of bottleneck or founder events, and expansion.

Recent study of epibionts of the R. exoculata gill chamber [START_REF] Petersen | Dual symbiosis of the vent shrimp Rimicaris exoculata with filamentous gamma-and epsilonproteobacteria at four Mid-Atlantic Ridge hydrothermal vent fields[END_REF] showed spatial clustering of bacterial sequences across the different sites (Rainbow, Logatchev and South MAR).

Such geographical segregation of the symbionts compared to host homogeneity may be due either to the influence of environment on symbionts, or to the very recent common history of R. exoculata populations with recent divergence only detectable on short generation symbionts but still undetectable on the genetics of their hosts.

The analysis of temporal variability revealed a lack of genetic differentiation among samples separated by time steps of one (Logatchev) and ten years (Rainbow), although within population diversity increased within a decade in Rainbow. In contrast to Logatchev, where no significant differences were observed across three consecutive years, Rainbow exhibited a significantly higher genetic diversity in 2007 compared with ten years earlier, despite the lower sample size and equal number of haplotypes (but more evenly distributed) in 2007. Genetic diversity, quantified both in terms of allelic richness and heterozygosity, increased significantly over the ten year interval in Rainbow, in agreement with a possible ongoing effect of population expansion (Table 2). Such trends were not detected in the diversity indices for the one year interval samples, possibly indicating that temporal variations in diversity are not sufficiently large to be detected over the shorter time scale of one year intervals. The pattern of variability across temporal samples therefore suggests a contemporary expansion, over one decade, at least in Rainbow. However, population growth tests are known not to perform well when population expansion is very recent [START_REF] Ramos-Onsins | Statistical properties of new neutrality tests against population growth[END_REF], and any expansion episodes detected using the growth decline model are, therefore, likely to be older than a decade. Altogether, with the star-like cluster of haplotypes revealing the accumulation of mutations over a recent evolutionary scale, these results imply either successive episodes of expansion or one long-term expansion episode affecting the pattern of haplotype diversity at both ecological and evolutionary time scales.

Considering the demographic history of these apparently rather young genetic lineages, which have probably undergone rather recent recolonization and/or expansion events, the lack of genetic differentiation among sites does not allow us to make inferences about the ongoing occurrence of extensive or restricted gene flow among the hydrothermal sites distributed along the Mid-Atlantic

Ridge. Yet, in the absence of any selective hypothesis explaining the star-like pattern of COI polymorphism, the hypothesis of such a recent and global (re)colonization also implies a high dispersal potential per se in the rather recent past. Recent bottlenecks and/or expansion have also been observed across populations of several hydrothermal invertebrates along other oceanic ridges (Young et al., 2008;[START_REF] Faure | Speciation in the deep sea: multi-locus analysis of divergence and gene flow between two hybridizing species of hydrothermal vent mussels[END_REF][START_REF] Plouviez | Comparative phylogeography among hydrothermal vent species along the East Pacific Rise reveals vicariant processes and population expansion in the South[END_REF] and such profiles of young lineages may be a rule rather than an exception for hydrothermal vent-associated species, due to the dynamic nature of vents subject to recurrent extinction events. In other ridges, such as the East Pacific Rise, most species studied to date also display a clear North/South genetic differentiation and in some cases even highly significant F ST values within the north and south populations [START_REF] Hurtado | Distinct patterns of genetic differentiation among annelids of eastern Pacific hydrothermal vents[END_REF][START_REF] Matabos | Geographic clines and stepping-stone patterns detected along the East Pacific Rise in the vetigastropod Lepetodrilus elevatus reflect species crypticism[END_REF][START_REF] Plouviez | Comparative phylogeography among hydrothermal vent species along the East Pacific Rise reveals vicariant processes and population expansion in the South[END_REF]. This genetic differentiation reported among populations of different vent species may be due to a lower dispersal potential and/or effective population size of these organisms, or to a higher occurrence or strength of physical and abiotic barriers to gene flow, such as transform faults, in the East Pacific Rise than in the Mid-Atlantic Ridge (but see O" [START_REF] O"mullan | A hybrid zone between hydrothermal vent mussels (Bivalvia: Mytilidae) from the Mid-Atlantic Ridge[END_REF]. Other species, however, do not exhibit such cleavage (Won et al., 2003;[START_REF] Hurtado | Distinct patterns of genetic differentiation among annelids of eastern Pacific hydrothermal vents[END_REF][START_REF] Plouviez | Comparative phylogeography among hydrothermal vent species along the East Pacific Rise reveals vicariant processes and population expansion in the South[END_REF], which may indicate a species-specific characteristic (such as dispersal potential or population size), or might simply reflect different sampling schemes. Large gaps in the sampling were indeed observed for most species where a north/south cleavage was shown and such a sampling scheme was demonstrated by [START_REF] Audzijonyte | When gaps really are gaps: statistical phylogeography of hydrothermal vent invertebrates[END_REF] to have the potential to result in artificially inflated F ST values. Yet, in the study of the sister mussel species Bathymodiolus azoricus and B. puteoserpentis and their hybrids [START_REF] O"mullan | A hybrid zone between hydrothermal vent mussels (Bivalvia: Mytilidae) from the Mid-Atlantic Ridge[END_REF], the clear distinction in species distribution between north (Rainbow, Lucky strike and Menez Gwen vents) and south (Snake Pit and Logatchev vents) along the MAR, suggests the existence of an ancient vicariance event that apparently does not affect Rimicaris exoculata, possibly due to a more recent history, higher effective population sizes, or higher dispersal potential in this shrimp species.

As few species have been studied so far in the MAR compared with the more extensively studied East Pacific Rise (EPR), a generalized comparison is not yet possible, but a species-specific hypothesis may be put forward to explain the lack of structure observed for R. exoculata. High dispersal potential, including long-distance migration as a bet-hedging strategy enabling escape from extinction, is thought to be a factor of risk in a fragmented environment due to the possible loss of progeny in non suitable habitats. In R. exoculata, this risk may be compensated by a high fecundity, as suggested for other species [START_REF] Llodra | Reproductive biology of three caridean shrimp, Rimicaris exoculata, Chorocaris chacei and Mirocaris fortunata (Caridea: Decapoda), from hydrothermal vents[END_REF], and/or by an environmentally driven delay in larval development, suggested by [START_REF] O"connor | Temperature control of larval dispersal and the implications for marine ecology, evolution, and conservation[END_REF] to be a possible adaptation associated with large scale dispersal. Moreover, it has been suggested that R. exoculata may detect black-body radiation associated with high-temperature vents (Van Dover et al., 1989[START_REF] Herring | Are vent shrimps blinded by science?[END_REF]; this may be a specific adaptation of this species to locate its exceptional habitat during large scale dispersal, allowing the benefit of the large scale dispersal to be kept while minimizing the associated risk of loss at the time of settlement. However, empirical data have shown that the detection of temperature in seawater may only be possible within a small area around the chimneys where the seawater is actually warmer [START_REF] Segonzac | Enigma of the trophic adaptation of the shrimp Alvinocarididae in hydrothermal areas along the Mid-Atlantic Ridge[END_REF]. The putative mechanism that may allow R. exoculata larvae to settle in their habitat after this hypothetical long-distance dispersal is therefore not elucidated, it may be temperature-associated or based on chemical signals linked to the presence of adults, as has been demonstrated for other marine invertebrates such as the oyster Ostrea puelchana, where preferential settlement of spat occurs close to the adults due to a chemical stimulus [START_REF] Pascual | Evidence of a chemically mediated adult-larval interaction triggering settlement in Ostrea puelchana: applications in hatchery production[END_REF].

This study makes a novel contribution to our understanding of the biogeography of the Mid-Atlantic Ridge by providing new insight into the population biology and biogeographic structure of a key species of shrimp from the unique, extreme and remote marine habitats of hydrothermal vents. Our results indicate that, across this very remote biogeographic region, all vent populations of this shrimp have a recent common history derived from a bottleneck or founder event, followed by a demographic expansion toward apparently large effective population sizes. Associated with high dispersal potential, such history possibly results in the absence of significant drift effects.

Interestingly, these results may indicate a possible common source or centre of origin of a (re)colonization event at the origin of the lineages sampled from the MAR in our study. The identification of such a potential source population is, unfortunately, rendered particularly difficult in these ecosystems because the sampling of hydrothermal vents is necessarily as fragmented in space and time as the ecosystems themselves are. Finally, more species should be studied on the MAR to compare results with the more extensively studied EPR (East Pacific Rise) in order to test for the possible influence of the shape of ridges and transform faults on the pattern of evolution and genetic differentiation of species associated with hydrothermal vents.

FiguresFigure 1 .

 1 Figures

Figure 2 .

 2 Figure 2. Haplotype network of the obtained mtDNA haplotypes of Rimicaris exoculata sampled

Figure 3 .

 3 Figure 3. Mismatch distribution of nucleotide pairwise differences for all the mitochondrial COI

Figure 4 .

 4 Figure 4. Observed and simulated mismatch curves of nucleotide pairwise differences of

Table 2 .

 2 Genetic diversity indices calculated for each sampled population and year, of the deep-sea shrimp Rimicaris exoculata sampled along the Mid-Atlantic Ridge. (n) sample size; (Nh) number of haplotypes; (Nph) number of private haplotypes, i.e. the number of haplotypes that only appeared in one population; (k) number of polymorphic sites; (h) haplotype diversity; (π 1 ) mean number of pairwise differences; (π 2 ) nucleotide diversity. Neutrality and population expansion tests: D= Tajima"s D-test; F S = Fu"s F S test. All values obtained for the neutrality tests were significant at the 5 % level.

	Population Year n k Nh	h	π 1	π 2	D	F S	Nph
	Rainbow	2007 26 18 14 0.82±0.07 1.67 0.0033 -2.292(0.002) -11.01(0)	3
		1997 32 14 14 0.69±0.09 1.05 0.0021 -2.298(0.001) -13.87(0)	7
	Logatchev 2007 31 20 16 0.77±0.08 1.41 0.0028 -2.487(0)	-15.38(0)	8
		2006 12 9	7	0.77±0.13 1.50 0.0030 -2.016(0.004) -3.52(0.003) 3
		2005 13 8	8	0.81±0.11 1.36 0.0027 -1.831(0.012) -5.25(0)	3
	TAG	2005 31 16 16 0.77±0.08 1.09 0.0022 -2.454(0)	-18.32(0)	6
	Ashadzé	2007 3	0	1	0	0	0	0	0	0
	SouthMAR 2006 4	3	4	1.00±0.18 1.50 0.0030 -0.754(0.23)	-2.37(0.008) 0

Table 3 .

 3 Genetic diversity indices for the temporal samples of the population of Rimicaris exoculata from Logatchev. Diversity indices for the sampling years of 2005 and 2007 are given as averages

	(±SD) due to the random re-sampling (100 times) to a minimum common size of 12 individuals per
	population. (Nh) number of haplotypes; (k) number of polymorphic sites; (h) haplotype diversity;
	(π 2 ) nucleotide diversity.
	Year	Nh	k	h	π 2
	2005	5.26±1.1	5.1±1.46	0.74±0.12	0.0024
	2006	7	9	0.77±0.13	0.0030
	2007	5.9±1.51	6.78±2.25	0.73±0.15	0.0026

Table 4 .

 4 Pairwise F ST values based on haplotype frequencies obtained for the species Rimicaris exoculata sampled along the Mid-Atlantic Ridge. None of the F ST values were significant at the 5 % level.

	Pops	Rainbow	Rainbow	Logatchev	Logatchev	Logatchev	TAG Ashadzé South
		07	97	07	06	05		MAR
	Rainbow07
	Rainbow97 0.00245
	Logatchev07 -0.0115	-0.01243
	Logatchev06 -0.0184	-0.0221	-0.0298
	Logatchev05 -0.0335	-0.0144	-0.0227	-0.0360
	TAG	-0.0147	-0.0138	-0.0122	-0.0262	-0.0260
	Ashadzé	0.0480	-0.0230	0.0140	0.0145	0.0331	0.014
							0
	South MAR -0.0589	0.0452	0.0031	-0.0147	-0.0495	0.003	0.250
							1

Acknowledgements

The authors wish to thank Alex Canton, Rui Candeias, Olivier Mouchel, Valérie Cueff-Gauchard and Stéphane L"Harridon, for technical help in the lab; Filipe Alberto for help with R; Didier Jollivet and Stéphane Hourdez for interesting discussions; Jillian Petersen and Nicole Dubilier for kindly providing samples from South MAR; the captain and crews of the vessels R/V Atalante, Pourquoi Pas? and Meteor, and the crew of the submersible Nautile and ROV Victor6000; and the cruise chiefs D. Desbruyers, A. Godfroy, Y. Fouquet, F. Gaill and J. Dyment for the cruises MARVEL 1997, EXOMAR 2005, SERPENTINE 2007, and MoMARDREAM-Naut 2007, respectively. We also thank Robert Vrijenhoek and an anonymous referee for the helpful comments and suggestions on this manuscript and Helen McCombie for the professional editing of this manuscript. S.T. was supported by a postdoctoral fellowship from FCT, Portuguese Science Foundation and ESF (European Social Fund), S.T, M-A C-B. and S.A-H were supported by Ifremer funds (DeepOases, ANR project).

Supporting information

Additional Supporting Information may be found in the online version of this article: Appendix S1 Neighbour-joining phylogram constructed with 152 COI mitochondrial sequences of

Rimicaris exoculata

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer-reviewed and may be re-organized for online delivery, but are not copy-edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Biosketches

The researchers involved in this study have a long-standing interest in Biogeography, Ecology and Table 1. Geographic coordinates, depth and cruise details for the populations of organisms sampled along the Mid-Atlantic Ridge.


Laurent Ros

Eric-Pierre Simon

SECOND-ORDER MODELING FOR RAYLEIGH FLAT FADING CHANNEL ESTIMATION WITH KALMAN FILTER

Keywords: Rayleigh Channel estimation, Jakes' Doppler spectrum, Flat fading, Kalman Filter, steady-state

This paper deals with channel estimation over flat fading Rayleigh channel with Jakes' Doppler spectrum. Many estimation algorithms exploit the time-domain correlation of the channel by employing Kalman Filter based on an approximation of the time-varying channel. A common method used in the literature is based on a first-order auto-regressive (AR1) model for the channel approximation, combined with a Correlation Matching (CM) criterion to fix the value of the AR1-parameter. In this paper, we propose first to replace the AR1 model by a specific second-order (Or2) model (as the one used for the phase estimation in presence of frequency offset), which is more appropriate for slow fading variations. Secondly, we propose a criterion based on the Minimization of the Asymptotic Variance (MAV) of the Kalman estimator to fix the parameter of the Or2-model. Closed-form expressions of the optimum Or2-parameter and of the corresponding Mean Square Error (MSE) are derived for a given channel state (Doppler spread, SNR). MSE theoretical analysis and simulation results prove a significant improvement with the Or2-MAV approach compared to more conventional AR1-CM approach (or even AR1-MAV optimized approach) in terms of MSE, especially for slow fading variations.

INTRODUCTION

The use of Kalman Filter (KF) for the channel estimation problem has received great attention in recent years in the wireless communication literature. It concerns a very large and various range of systems, as for example Multi-Input-Multi-Output systems [START_REF] Komninakis | Multi-Input Multi-Output Fading Channel Tracking and Equalization Using Kalman Estimation[END_REF] [START_REF] Liu | Space-time coding and Kalman Filtering for time selective fading channels[END_REF], or Orthogonal Frequency Division Multiplexing systems [START_REF] Chen | Kalman filter Channel estimator for OFDM systems in time and frequency-selective fading environment[END_REF] [START_REF] Al-Naffouri | An EM-based Forward-Backward Kalman Filter for the Estimation of Time-Variant Channels in OFDM[END_REF][7] [START_REF] Ros | Paths Complex Gain Tracking Algorithms for OFDM Receiver in Slowly-Varying Channels[END_REF], or again GSM systems [START_REF] Abeida | An EM algorithm for path delay and complex gain estimation of a slowly varying fading channel for CPM signals[END_REF]. All the aforementioned works based their Kalman Filter on the autoregressive (AR) approximation of the widely accepted Rayleigh fading channel with Jakes' Doppler spectrum [START_REF] Jakes | Microwave Mobile Communications[END_REF], as developed in [START_REF] Baddour | Autoregressive modeling for fading channel simulation[END_REF]. Moreover, the tuning of the AR model parameters with respect to a given normalized Doppler frequency ( f d T ) is always based on a same criterion, called Correlation Matching (CM) criterion in this paper (see section 3.1.2). Note also that this AR-CM method to approach the correlated fading channels is also used today for the computer simulation of the "Rayleigh fading channel model with Jakes' Doppler spectrum", as for example in the actual version of Matlab (function "rayleighchan"). The authors of [START_REF] Baddour | Autoregressive modeling for fading channel simulation[END_REF] insisted that low orders are appropriate for narrowband Doppler fading processes, and [START_REF] Wang | On verifying the First Order Markovian Assumption for a Rayleigh Fading Channel Model[END_REF] showed that a first order assumption is enough to capture most of the channel tap dynamics. Finally, in the KF based channel tracking problem, the combination "AR1 model" (or Gauss-Markov assumption) and "CM criterion" is often retained [START_REF] Komninakis | Multi-Input Multi-Output Fading Channel Tracking and Equalization Using Kalman Estimation[END_REF][5] [START_REF] Al-Naffouri | An EM-based Forward-Backward Kalman Filter for the Estimation of Time-Variant Channels in OFDM[END_REF][7] [START_REF] Abeida | An EM algorithm for path delay and complex gain estimation of a slowly varying fading channel for CPM signals[END_REF]... But for the case of slow fading scenario, the MSE results with the so-called AR1(CM)-Kalman estimator were somewhat disappointing (compared to Bayesian Cramer Rao Bound (BCRB), or to the MSE of lower complexity algorithms), as we pointed out recently in [START_REF] Ros | Paths Complex Gain Tracking Algorithms for OFDM Receiver in Slowly-Varying Channels[END_REF]-Fig. 2 (see also [START_REF] Abeida | An EM algorithm for path delay and complex gain estimation of a slowly varying fading channel for CPM signals[END_REF]-Fig. 3). It should be noted that the slow channel variations assumption ( f d T ≤ 10 -2 ) corresponds to the most common transmission scenario (e.g. [3][4][5] [START_REF] Al-Naffouri | An EM-based Forward-Backward Kalman Filter for the Estimation of Time-Variant Channels in OFDM[END_REF]) where the channel variation within one symbol duration can be neglected (and then also the symbol distortion in a mono-carrier scenario, or the Inter-Carrier-Interference in a multi-carrier system [START_REF] Hijazi | Joint Data QR-Detection and Kalman Estimation for OFDM Time-varying Rayleigh Channel Complex Gains[END_REF]). Note also that some authors propose to modify slightly the CM constraint to better approximate the original process in some sense (for example [START_REF] Baddour | Autoregressive modeling for fading channel simulation[END_REF] propose to add a positive ε to the zero-th autocorrelation lag but without explicit analysis to fix the value of ε). The authors of [START_REF] Ghandour | Estimation de canal vanouissement lent pour les communications radio-mobiles[END_REF][START_REF] Ghandour-Haidar | On the Use of first-order Autoregressive Modeling for Rayleigh Flat Fading Channel Estimation with Kalman filter[END_REF] explained analytically for the flat fading case the previous disappointing results and propose to compute the AR1-parameter under a Minimum Asymptotic Variance (MAV) criterion without imposing the CM constraint. The resulting MSE performance was then improved, but remains far from the BCRB.

In this paper, we propose a second-order model based KF designed for the flat fading Rayleigh channel estimation under the slow channel variations scenario. The 2nd-order model, called Or2-model in this paper, generalizes (for a complex parameter) the one used for the phase estimation scalar problem in presence of frequency offset [13][14]. This model is based on a Brownian evolution of the slope of the parameter to be tracked. The main motivation for this model is that for the slow fading Rayleigh channel ( f d T ≤ 10 -2 assumed), the channel Complex Gains (CGs) exhibit strong trend behavior, i.e. they continue in some direction during several symbols, as shown in Fig. 1-left. So it is satisfactory that the Or2-model takes into account a drift (or slope) into the CG variation contrary to the AR1-model, with a (slow) possible evolution of this slope. An important question that arises is about the choice (or tuning) of the Or2 model-parameter for a given channel state (Doppler frequency f d T of the Jakes' spectrum, Signal to Noise Ratio SNR). Considering the tracking mode where the KF can be regarded as a linear time-invariant filter, we give a MSE analysis that permits to express the asymptotic (or steady state) tracking error variance in terms of static and dynamic error variance. Then the Or2-parameter that permits to minimize the asymptotic global variance of the KF estimator (i.e. under MAV criterion) is derived. Simulation results and comparison to literature AR1-Kalman filters and to BCRB validate the proposed algorithm and the theoretical analysis.

This paper is organized as follows: section 2 describes the model and objectives, section 3 recalls Kalman Filters equations for AR1-KF and Or2-KF. Section 4 gives MSE analysis of the Or2-KF and the proposed optimization. Finally, the different results are discussed in Section 5.

MODEL AND ESTIMATION OBJECTIVES

We consider the estimation of a flat Rayleigh fading channel. The discrete-time observation is 1 :

y (n) = α (n) + N (n) ( 1
)
where n is the symbol time index, N (n) is a zero-mean additive white circular complex Gaussian noise with variance σ2 N , and α (n) is a zero-mean circular Gaussian channel CG 1 Model (1) assumes that symbols are normalized and known (or decided), additionally to flat fading assumption. Although this model is admittedly simplistic, it can be applied to different (more involved) contexts, such as pilot-aided multicarrier systems in frequency-selective wireless channels.

with variance σ 2 α = 1. The normalized Doppler frequency of this channel is f d T , where T is the symbol period. A Jakes' Doppler spectrum is assumed for this channel:

Γ α ( f ) = ⎧ ⎨ ⎩ σ 2 α π f d 1-f f d 2 if | f | < f d 0 i f | f | > f d (2)
The autocorrelation coefficient R α [m] of the stationary CG α is then defined for lag m by:

R α [m] = E{α (n) .α (n-m) * } = σ 2 α J 0 (2π f d T.m) (3)
where J 0 is the zeroth-order Bessel function of the first kind.

Given the observation model ( 1) and the Doppler spectrum statistical constraint (2) for the dynamic evolution of the CG, we look for an on-line unbiased estimation α(n) of α (n) based on KF. The variance σ 2

ε def = E |ε (n) | 2 of the estimation error ε (n) def = α (n) -α(n) will be investigated.

KALMAN FILTERS

Review of the AR1-Kalman Filter

Dynamical model for approaching the CG variation

The time-varying CG α can be approached by a first order autoregressive (AR1) model α AR1 :

α AR1 (n) = γ.α AR1 (n-1) + e (n) (4)
where γ is a positive real such that γ ∈ [0; 1[ and e (n) is a white circular complex Gaussian noise with variance such that σ 2 e = (1 -γ 2 ).σ 2 α . According to (4), the AR1-coefficient γ verifies

γ = R α AR1 [1] R α AR1 [0] (5)
The specific equations of the AR1-KF can be found in [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF].

Correlation Matching criterion [1][3]

The CM criterion consists in imposing that the autocorrelation coefficients R α AR1 [m] of the approximated AR1 process α AR1 perfectly match the sampled autocorrelation function R α [m] (which is a Bessel function) of the true CG α, for lags 2 m ∈ {-1, 0, 1} . Using the hermitian symmetry of a stationary correlation matrix, it results that

R α AR1 [0] = R α [0],
and

R α AR1 [1] = R α [1]. Regarding equations (5)&(3) the AR1- coefficient γ noted γ CM results in γ CM = J 0 (2π f d T ) (6)
It is shown in [START_REF] Ghandour | Estimation de canal vanouissement lent pour les communications radio-mobiles[END_REF][START_REF] Ghandour-Haidar | On the Use of first-order Autoregressive Modeling for Rayleigh Flat Fading Channel Estimation with Kalman filter[END_REF] that for low Doppler and medium SNR, the asymptotic variance of the AR1(CM)-Kalman Filter is:

σ 2 ε (AR1(CM)) ≈ σ 2 N + π √ 2 f d T σ α σ N (7)
This standard choice of J 0 (2π f d T ) for γ results from the CM constraint. However, note that imposing the matching of three taps (for p = 1) of the approximate and original processes auto-correlation functions does not ensure the best similitude (in terms of euclidean distance for example) between the two functions. This is especially true for low f d T << 1 (see [START_REF] Komninakis | Multi-Input Multi-Output Fading Channel Tracking and Equalization Using Kalman Estimation[END_REF], Fig. 1) where the three taps are very close to the value 1 (since J 0 (2π

f d T ) ≈ 1 -1 4 .(2π f d T ) 2 ).

Minimum Asymptotic

Variance criterion [START_REF] Ghandour | Estimation de canal vanouissement lent pour les communications radio-mobiles[END_REF][START_REF] Ghandour-Haidar | On the Use of first-order Autoregressive Modeling for Rayleigh Flat Fading Channel Estimation with Kalman filter[END_REF] For low Doppler and low SNR, [START_REF] Ghandour | Estimation de canal vanouissement lent pour les communications radio-mobiles[END_REF][START_REF] Ghandour-Haidar | On the Use of first-order Autoregressive Modeling for Rayleigh Flat Fading Channel Estimation with Kalman filter[END_REF] showed that the AR1coefficient γ noted γ MAV that allows to minimize the asymptotic estimation error variance σ 2 ε of the AR1-KF is

γ MAV = 1 -4 (π f d T ) 4 σ 2 N /σ 2 α 1 3 (8)
and the closed-form expression of the corresponding asymptotic estimation error variance is

σ 2 ε (AR1(MAV )) ≈ 3 2 .(σ 2 α ) 1 3 . π f d T σ 2 N 2 3 (9)

Or2-Kalman Filter

Dynamical model for approaching the CG variation

In case of slow fading, the CG variations look linear during a few symbols. It could then be more appropriate to consider 2nd-order model including a (slightly mobile) linear drift δ , where the Jakes'process α (n) is approached by α(n) as

α(n) = α(n-1) + δ (n-1) (10) δ (n) = δ (n-1) + u (n) (11)
where u (n) is zero mean Gaussian complex circular with a variance σ 2 u . And the parameter of this model, σ 2 u , has to be calibrated such that the slope (or drift) of the CG variation, δ (n) , changes slowly with time n, according to the actual value of f d T . The 2nd-order model of the CG evolution can be re-formulated in a state-space model. The state vector to be considered includes the CG and the drift,

a (n) = [ α(n) , δ (n) ] T .
The state evolution matrix is M = 1 1 0 1 and the state-noise

vector is u (n) = [0, u (n) ]
T . The observation matrix with size

1 × 2 is S (n) = s T def = [1, 0]
. The state evolution ( 10)&( 11) and the observation (1) become in a state-space formulation:

a (n) = Ma (n-1) + u (n) ( 12
)
y (n) = s T a (n) + N (n) (13)

Equations of the Or2-Kalman Filter

Regarding the state-space formulation ( 12)&( 13), the two stages of the so called Or2-Kalman Filter are given by [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]:

Time Update Equations:

â(n|n-1) = Mâ (n-1|n-1) (14)
P (n|n-1) = MP (n-1|n-1) M H + U (15)
Measurement Update Equations:

K (n) = P (n|n-1) s s T P (n|n-1) s + σ 2 N ( 16
)
â(n|n) = â(n|n-1) + K (n) y (n) -s T â(n|n-1) (17)
P (n|n) = I 2 -K (n) s T P (n|n-1) (18)
where

K (n) = k 1(n) k 2(n)
is the Kalman gain vector, U = 0 0 0 σ 2 u , and P (n|n) and P (n|n-1) are respectively the 2 × 2 a posteriori and predicted error covariance matrices.

ANALYSIS AND OPTIMIZATION OF THE OR2-KF

Steady-state Or2-KF equations

Since the linear system ( 12)&( 13) is observable and controllable, an asymptotic regime is quickly reached ( [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]), for which the Kalman gain and covariance matrices become constant. For our second-order model, we can write for sufficiently large n that

K (n) = K (n+1) = K ∞ def = k 1 k 2 ,
α(n|n) = α(n-1|n-1) + δ(n-1|n-1) + k 1 .v ε(n) (19) δ(n|n) = δ(n-1|n-1) + k 2 .v ε(n) (20) with v ε(n) def = y (n) -( α(n-1|n-1) + δ(n-1|n-1) ) (21)
The filter is equivalent to a second-order tracking loop parametrized by 2 constants k 1 and k 2 , that play respectively the role of the Proportional (P) and Integral (I) gains of the loop ( [START_REF] Patapoutian | On phase-locked loops and Kalman filters[END_REF][14] [START_REF] Ros | Paths Complex Gain Tracking Algorithms for OFDM Receiver in Slowly-Varying Channels[END_REF]). But here the values of these 2 parameters can not directly be chosen: they depend on the channel state ( f d T , SNR = σ 2 α /σ 2 N ) and on the model noise variance σ 2 u , which is the only parameter that the user can tune. In the following, we aim to express (k 1 , k 2 ) and then the asymptotic error variance σ 2 ε with respect to σ 2 u , in the perspective of an optimization.

Steady-state Kalman gain

To find the steady-state Kalman gain vector (k 1 , k 2 ) we need first to find the matrix P' ∞ , and then to resolve the so-called Riccati equations ( 16)&( 15)&(18) rewritten for n → ∞ as: We can easily deduce that k 2 P 11 = k 1 P 21 , k 1 P 12 = P 22 and k 2 P 12 = σ 2 u and thus (using the hermitian symmetry of the covariance matrices and ( 22)) that k 1 , k 2 , P' ∞ , and P ∞ are real-valued objects. It results that P 21 = P 12 = σ u P 11 + σ 2 N , and equation ( 22) provides an exact expression of the KF gain with respect to P 11 and σ 2 N , where P 11 can be found after some manipulations of the Riccati equations as the solution of ( [START_REF] Qian | Steady-State Performance of Kalman Filter for DPLL[END_REF]):

k
P 4 11 = σ 2 u (P 11 + σ 2 N )(P 11 + 2σ 2 N ) 2 (25)
Now, in order to obtain more tractable formulation, we make the assumption P 11 << σ 2 N , which means that we have a low Kalman gain k 1 << 1 (see ( 22)). This assumption a priori agrees with a low Doppler ( f d T ≤ 10 -2 ) and low SNR (σ N ∈ [0.01, 1]) scenario. We get then from (25) the approximation

P 11 ≈ σ 2 N 2σ u /σ N , which is injected in (22)
to finally obtain:

k 1 ≈ P 11 σ 2 N ≈ 2σ u σ N (26) k 2 ≈ P 21 σ 2 N ≈ σ u σ N ( 27
)
Note that k 2 , k 1 are then linked by k 1 ≈ √ 2k 2 . Note also that our assumption implies that the model noise is chosen weak compared to the observation noise, i.e. σ 2 u << σ 2 N (see ( 26)).

Steady-state MSE and optimization

Steady-state transfer function of the KF and MSE

Using ( 20)&( 21), the equation ( 19) of the steady-state Kalman Filter can be expressed in Z-Transform domain by:

α(z).[1 -z -1 ] = [k 1 + k 2 .z -1 1 -z -1 ].v ε (z) (28)
The error signal (21) can be rewritten versus the estimation error from ( 21)&( 19)&(1) as

v ε(n) = 1 1 -k 1 .{α (n) -α(n|n) } + 1 1 -k 1 .N (n) (29)
Combining ( 29) and (28) leads in the Z-transform domain to α(z) = L(z).(α(z) + N(z)) and then to the estimation error

ε(z) = α(z) -α(z) = (1 -L(z)).α(z) -L(z).N(z) (30)
where L(z) is the transfer function of the steady-state KF:

L(z) = k 1 -k 2 1-k 1 .(1 -z -1 ) + k 2 1-k 1 (1 -z -1 ) 2 + k 1 -k 2 1-k 1 .(1 -z -1 ) + k 2 1-k 1 (31)
In using the previous assumptions, we have that 0 < k 2 << k 1 ≈ √ 2k 2 << 1, and then L(z) is a low-pass filter which can be approximated in frequency domain (using z = e j2π f T and then 1z -1 ≈ j2π f T for f T << 1) by

L(e j2π f T ) ≈ √ 2k 2 .( j2π f T ) + k 2 ( j2π f T ) 2 + √ 2k 2 .( j2π f T ) + k 2 (32)
This corresponds to a second-order low-pass filter with normalized natural pulsation

ω n T ≈ √ k 2 ≈ k 1 √ 2 ,
and with imposed damping factor ζ = √ 2 2 . Using (30), the asymptotic MSE is divided into two parts:

σ 2 ε = σ 2 εα + σ 2 εN ( 33
)
• the dynamic error variance σ 2 εα is due to the variations of the CG α (n) filtered by the high pass filter 1 -L(z):

σ 2 εα def = + 1 2T -1 2T |1 -L(e j2π f T )| 2 .Γ α ( f )d f ≈ + f d -f d (2π f T ) 4 k 2 2 .Γ α ( f )d f = 3 8 .σ 2 α . ( 2π
f d T ) 4 k 2 2 (34)
The approximation in (34) is made in using the previous approximation of L(.) in (32), and in assuming low Doppler f d T << 1, and then a Kalman gain (or a natural pulsation) such that 2π f d T << √ k 2 << 1. The exact integration is then computed from Γ α ( f ) defined in (2) by the way of the variable change cosθ = ( f / f d ).

• the static error variance σ 2 εN is due to the additive noise N (n) filtered by the low pass filter -L(z):

σ 2 εN def = σ 2 N × T + 1 2T -1 2T |L(e j2π f T )| 2 d f (35) = σ 2 N × 2k 2 1 -3k 1 k 2 + 2k 2 k 1 .[4 -(2k 1 + k 2 )] (36) ≈ σ 2 N × 3 √ 2 4
. k 2 (37)

The exact integration (36) is obtained from exact L(z) in (31) in using method of the book [START_REF] Jury | Theory and Application of the Z-Transform Method[END_REF]. The approximation in (37) is obtained in using 2k 2 ≈ k 2 1 and k 2 1 << k 1 << 1. (38)

And the corresponding optimal MSE results in

σ 2 ε (Or2(MAV )) = 15 8 ( √ 2π) 4 5 . σ 2 α 1 5 • σ 2 N • f d T 4/5 (39)

SIMULATION RESULTS AND DISCUSSION

Fig. 2 and Fig. 3 compare the MSE obtained with the AR1-KF (with CM and MAV criteria) and with the Or2-KF (with MAV criterion) by the way of Monte-Carlo simulations. We also plot the on-line BCRB ( [START_REF] Hijazi | Bayesian Cramer-Rao Bounds for OFDM Slowly Varying Multi-path Rayleigh Channel Estimation[END_REF]) as reference. We verify by these figures that:

• With the CM criterion, the MSE of the AR1-KF is approximately constant with respect to the Doppler frequency for f d T ≤ 10 -2 (which agrees with (7) for medium SNR), and is far from the BCRB for lower f d T , as was observed in [9] [START_REF] Ros | Paths Complex Gain Tracking Algorithms for OFDM Receiver in Slowly-Varying Channels[END_REF]. Furthermore, the use of the MAV versus CM criterion improves the AR1-KF strongly, especially for lower Doppler frequencies and low SNRs, which corroborates the study of [START_REF] Ghandour | Estimation de canal vanouissement lent pour les communications radio-mobiles[END_REF][START_REF] Ghandour-Haidar | On the Use of first-order Autoregressive Modeling for Rayleigh Flat Fading Channel Estimation with Kalman filter[END_REF].

• For the Or2-KF, the MSE computed by Monte-Carlo simulation is very close to σ 2 ε (Or2(MAV )) obtained by the closed-form expression (39), so we validate our theoretical analysis and our approximations.

• The most important point is that with the proposed Or2-KF and its optimization, we are closer to the BCRB than with the AR1-KF (under CM or MAV criteria), especially for lower normalized Doppler frequencies f d T ≤ 10 -2 . This point reveals the advantage of using an appropriate 2nd-order modeling versus a first-order one in slow fading scenario, as we had already pointed out in [START_REF] Ros | Paths Complex Gain Tracking Algorithms for OFDM Receiver in Slowly-Varying Channels[END_REF]. The present study shows that this 2nd-order model allows a KF MSE that is proportional to the ( 4 5 ) power of f d T , versus the ( 23 ) power for the 1st-order AR model, in agreement with (39) and (9).

CONCLUSION

The modeling for the Kalman Filter based Rayleigh channel estimation was investigated, together with the tuning of the parameters (γ or σ 2 u ) of the models. The simple case of a flat fading channel was considered, but the results can be applied or generalized to more complex systems. In the literature, the 1st-order Auto-Regressive (AR1) model is often used, with a Bessel AR1-coefficient resulting from the standard Correlation Matching criterion. This choice seems to be convenient for high Doppler frequencies and high SNR, but not for lower ones. For most conventional Doppler speeds ( f d T ≤ 10 -2 ), we have demonstrated the interest of using a specific 2ndorder model (based on a Brownian evolution of the Complex Gain slope) versus 1st-order model. We have given closed form expressions of the optimal model-parameter that minimizes the asymptotic estimation error variance assuming Rayleigh-Jakes channel. It is demonstrated that the MSE of the 2nd-order KF changes proportionally to the ( 4 5 ) power of the product ( f d T × σ 2 N ), where σ 2 N is the observation noise variance. Simulation results proved an amelioration in terms of MSE, in the choice we proposed.

Fig. 1 .

 1 Fig. 1. Example (one realization) of the temporal evolution of a (Rayleigh-Jakes) Complex Gain (real part) over a window of 6 symbol durations, for a slow fading case (left, with f d T = 10 -3 ), and a very fast fading case (right, with f d T = 0.3)

Fig. 2 .

 2 Fig. 2. MSE of the Kalman Filters for the proposed Or2(MAV)-KF compared to literature AR1(CM)-KF and AR1(MAV)-KF versus SNR for f d T = 10 -4 (a), f d T = 10 -3 (b), and f d T = 10 -2 (c)

Fig. 3 .

 3 Fig. 3. MSE in function of f d T , for SNR = 20 dB.

  11 P 12 P 21 P 22 = P 11 + P 12 + P 21 + P 22 P 12 + P 22 P 21 + P 22 P 22 + σ 2

	1 =	P 11 P 11 + σ 2 N	and k 2 =	P 21 P 11 + σ 2 N	(22)
						(23)
					u
	P 11 P 12 P 21 P 22	=	(1 -k 1 )P 11 -k 2 P 11 + P 21 -k 2 P 12 + P 22 (1 -k1)P 12	(24)

P

For a more general AR-p model (with then more degrees of freedom than the AR1-model), the coefficients of the AR-p process are calculated[START_REF] Baddour | Autoregressive modeling for fading channel simulation[END_REF] by still imposing the correlation matching constraint for lags m ∈ {-p, ...., -1, 0, 1, ...., p}, but followed by the resolution of the Yule-Walker equations in order to minimize the prediction error (model noise).


Pierre-Yves Angrand
email: angrand@enst.fr

Jacques Sakarovitch
email: sakarovitch@enst.fr.

On the enumerating series of an abstract numeration system

It is known that any rational abstract numeration system is faithfully, and effectively, represented by an N-rational series. A simple proof of this result is given which yields a representation of this series which in turn allows a simple computation of the value of words in this system and easy constructions for the recognition of recognisable sets of numbers.

It is also shown that conversely it is decidable whether an N-rational series corresponds to a rational abstract numeration system.

Introduction

In order to state our result, we have first to recall the definition -due to Lecomte and Rigo [START_REF] Lecomte | Numeration systems on a regular language[END_REF] -of an abstract numeration system and, in order to motivate it, the more common one of numeration systems.

Numbers do exist independently of the way we represent them, and operations on numbers are defined independently of the way they are computed. The role of a numeration system is to set a framework in which numbers are represented by words (over a suitable alphabet) allowing to describe operations on numbers as algorithms on the representations, that is, on words.

The most common numeration system -in our modern times -is the k-ary system where numbers are given their representation in base k, that is, written as words over the alphabet A k = {0, 1, . . . , k -1} and which do not start with 0 (but for the representation of 0 itself). The sequence of the representations of the integers in the binary system is: {0, 1, 10, 11, 100, 101, 110, . . .} .

While keeping the notion of position numeration system, the k-ary systems can be generalised by replacing the sequence (k n ) n 0 with some increasing sequence U = (U n ) n 0 of integers such that U 0 = 1 . Using a greedy algorithm, every integer n is then given a representation in the 'base' U , called its U -representation and denoted by n U . A well-known example is the Fibonacci numeration system based on the sequence F = (F n ) n 0 of Fibonacci numbers starting with F 0 = 1 and F 1 = 2 . In this system, every positive integer is given a canonical representation which is computed by the greedy algorithm and which is characterised by the fact it does not contain 11 as a factor. The sequence of the representations of the integers in the Fibonacci system is: {0, 1, 10, 100, 101, 1000, 1001, 1010, . . .} .

It is possible to look at these two numeration systems, the 2-ary system and the Fibonacci system, independently from the sequences (2 n ) n 0 and (F n ) n 0 and the greedy algorithm, and by just considering the set of words that represent the integers: 1{0, 1} * ∪{0} in the first case, 1{0, 1} * \{0, 1} * 11{0, 1} * ∪{0} in the second case and by enumerating the element of this set in the radix order. 1 In both cases, every integer will be given the same representation without reference to the way this representation is computed. It is the language of all representations that matters and this naturally leads to the definition of abstract numeration systems.

Definition 1 ([14]

). An abstract numeration system (or ANS for short) is a triple S = (L, A, <) where A is an alphabet equipped with a total order < and L is an infinite language of A * .

The system S allows to define a one-to-one correspondence between N and L by associating every integer n with the (n + 1)-th word of L in the radix order defined on A * by <. This representation of n is denoted by n S and conversely the corresponding value of a word w of L is denoted by π S (w). Of course, the following holds:

π S (w) S = w and π S ( n S ) = n .
In most cases, the alphabet A and the order < on A are fixed and understood and we speak of the ANS defined by the language L and we use the simpler notations n L and π L (w).

If L is a rational language of A * , we say that the ANS is rational. Beyond the irrepressible appeal to generalisation and abstraction, a true motivation that supports the definition of ANS is to understand which properties of a numeration system depend upon the whole language of the representations only, and which are more directly related to the way the representation of every number is computed. For instance, we have shown in a previous paper [START_REF] Angrand | Radix enumeration of rational languages[END_REF] that the successor function in a rational ANS is a piecewise cosequential function, whereas the characterisation of those systems for which this successor function is co-sequential is known in the case of β-numeration systems (cf. [START_REF] Frougny | On the sequentiality of the successor function[END_REF]) but seems to be out of reach for arbitrary rational ANS so far.

The purpose of this paper is to set up even tighter bonds between rational abstract numeration systems and classical automata theory. We reach this goal via the definition of the enumerating series of a numeration system and with the use of its representation in the case it is rational. Definition 2. Let S = (L, A, <) be an abstract numeration system. The enumerating series of S is the N-series over A * denoted by E S and defined by:

E S = w∈L (π S (w) + 1) w .
As above, the notation can be simplified as E L = w∈L (π L (w) + 1) w .

Remark 1. The above definition has been taken so that the language L is entirely determined by E L . Indeed, L = supp (E L ) .

One certainly could have taken E L = w∈L (π L (w)) w as a definition for E L .

All the results we are going to describe would have been valid and it may have looked more natural. But we would have lost the information on the first word of L, that is, the representation of 0.

The starting point of our work is a direct proof of the following result (the definition of N-rational series will be recalled below).

Theorem 1 ([7]

). The enumerating series of a rational abstract numeration system is an N-rational series.

In [START_REF] Choffrut | Rational transductions and complexity of counting problems[END_REF], Theorem 1 was a corollary of constructions set up for establishing the rationality or algebraicity of a family of counting problems by means of rational transductions. Theorem 1 was also given another and specialised proof in [START_REF] Rigo | Numeration systems on a regular language : arithmetic operations, recognizability and formal power series[END_REF]. Even if both this and the original proofs are effective, the one we give below in Section 3 amounts to compute directly a representation of E L from a representation of (the characteristic series of) L and also to give an even more compact algorithm for calculating the coefficient of a word w in E L , that is, the value of w in the system L increased by 1. We then deduce from this latter algorithm the construction of the automaton that recognises the set of representations in the system L of a recognisable set of numbers (Section 4). It is to be noted that the same last construction was also given in [START_REF] Krieger | Decimation of languages and state complexity[END_REF] (cf. Remark 6).

The next result plays the role of a converse of Theorem 1: of course not every N-rational series is the enumerating series of a rational number system, but one can at least know when it is the case. Theorem 2. It is decidable whether an N-rational series is the enumerating series of a (rational) abstract numeration system or not.

We end the paper with some problems that are directly inspired by Theorem 1.

Preliminary and notation

This paper makes use of several notions of automata theory such as unambiguous, or deterministic, automata, rational series and languages, with which the reader is supposed to be familiar. Definitions that are not given here are to be found in reference books such as [START_REF] Eilenberg | Automata, Languages and Machines[END_REF][START_REF] Berstel | Rational Series and their Languages[END_REF][START_REF] Sakarovitch | Eléments de théorie des automates[END_REF]. Our notation are mainly those used in [START_REF] Sakarovitch | Eléments de théorie des automates[END_REF].

In the sequel, A is a finite alphabet, A * the free monoid generated by A, 1 A * the empty word, identity of A * . The length of a word w in A * is denoted by |w|. Let A be totally ordered by <. The radix order ≺ on A * is defined by:

2 u ≺ v if either |u| < |v| , or |u| = |v| , u = w a u ′ , v = w b v ′ and a < b .
The radix order is a well order, that is, every non empty subset of A * has a smallest element for ≺ and can thus be used to enumerate any subset of A * . Let K be a semiring; for instance, N, the semiring of non negative integers. A (K-)series s (over A * ) is a map from A * to K, and the image of a word w by s is called the coefficient of w in s and is denoted by <s, w>. The set of series over A * with coefficients in K is denoted by K A * . The support of a series s is the language, denoted by supp s, which contains those words whose coefficient in s is different from 0 K . Conversely, the characteristic series of a language L of A * is the N-series, denoted by L, defined by <L, w> = 1 if w is in L and <L, w> = 0 otherwise.

We call (K-)representation, of dimension n, a triple (λ, µ, ν) where µ is a morphism µ : A * → K n×n from A * to the n × n-matrices with entries in K, and λ and ν are two vectors of dimension n with entries in K, λ a row vector and ν a column vector. A series s in K A * is (K-)recognisable if there exists a representation (λ, µ, ν) such that, for every w in A * , <s, w> = λ • µ (w) • ν .

A series is (K-)rational if it is the behaviour of a finite (K-)automaton, that is, an automaton with multiplicity in K (the behaviour of an automaton A is the series where the coefficient of a word w is the sum of the multiplicities of all computations in A with label w). Finite K-automata whose transitions are labelled by letters and K-representations are two ways to describe the same concept. 3 The illustration given with Example 1 suffices for the definition. As every K-automaton is equivalent to one which is labelled by letters, the families of K-rational and Krecognisable series coincide.

Example 2. The language L 1 of words with an even number of b's is recognised by the automaton A 1 drawn at Figure 1. The representation (λ 1 , µ 1 , ν 1 ) associated with A 1 is The proof of Theorem 1, as given in [START_REF] Choffrut | Rational transductions and complexity of counting problems[END_REF] where Theorem 1 is (part of) Corollary 8, is based on the construction of an unambiguous rational transduction that associates to every word u all words v that are greater than u in the radix order. From this, it is easy to derive that the image of the characteristic series of a rational language L is a recognisable series, and equal to E L -up to the intersection (or Hadamard product) with L. The advantage of this construction is that it can be applied to unambiguous context-free languages and to various other counting functions as well. The inconvenient is that it does not provide directly the representation of E L , although it is very similar to the one we develop below.

λ 1 = 1 0 , µ 1 (a) = 1 0 0 1 , µ 1 (b) = 0 1 1 0 , ν 1 = 1 0 . b b a a
In [START_REF] Rigo | Numeration systems on a regular language : arithmetic operations, recognizability and formal power series[END_REF], Theorem 1 is Proposition 29; its proof is more direct than in [START_REF] Choffrut | Rational transductions and complexity of counting problems[END_REF] in the sense it does not rely on the rational transduction machinery but makes use instead of the characterisation of recognisable series as those which belong to a finitely generated stable submodule of K A * . But this proof yields neither the representation of E L nor a simple mean to compute it.

Preparation

If a is a letter of A, let us denote by A a the set of letters of A smaller than a:

A a = {b ∈ A | b < a} .
If u be a word of A * , let us denote by P (u) the set of words of A * (strictly) smaller than u in the radix order:

P (u) = {v ∈ A * | v ≺ u} .
This set P (u) can be defined by induction on the length of u by the following remark. Any word smaller than u followed by any letter is smaller than u a, and so is u b for any letter b smaller than a, and the empty word is also smaller than u a. These three sets are pairwise disjoint and any word smaller than u a falls in one of them. Altogether, we have proved the following lemma:

4 Lemma 3. ∀u ∈ A * , ∀a ∈ A P (u a) = 1 A * ∪ uA a ∪ P (u)A .
Let L be a rational language of A * and (λ, µ, ν) the N-representation which corresponds to an unambiguous finite automaton which recognises L:

∀w ∈ A * λ • µ (w) • ν = 1 ⇐⇒ w ∈ L .
We use the following notation: if K is a (finite) subset of A * , then µ (K) = w∈K µ (w) . As (λ, µ, ν) corresponds to an unambiguous automaton, we have:

∀K ⊆ A * λ • µ (K) • ν = w∈K λ • µ (w) • ν = card (K ∩ L) . (1)

Proof of Theorem 1

Let S = (L, A, <) be a rational ANS, A an unambiguous automaton that recognises L, (λ, µ, ν) the corresponding N-representation, and k its dimension. From (1) follows:

∀u ∈ A * λ • µ (P (u)) • ν = card ({v ∈ A * | v ∈ L and v u}) .
and thus:

∀w ∈ L λ • µ (P (w)) • ν = π L (w) .
From Lemma 3 follows:

5 ∀u ∈ A * , ∀a ∈ A λ • µ (P (u a)) • ν = λ • µ (1 A * ) • ν + λ • µ (u) • µ (A a ) • ν + λ • µ (P (u)) • µ (A) • ν . ( 2
)
Let σ = µ (A) and, for every a in A, σ a = µ (A a ). Thus (2) is rewritten as:

∀u ∈ A * , ∀a ∈ A λ • µ (P (u a)) • ν = λ • ν + λ • µ (u) • σ a • ν + λ • µ (P (u)) • σ • ν . (3)
Let (η, κ, ζ) be the representation of dimension 2k + 1 described by the following (1, k, k)-block decomposition:

η = 1 λ 0 , ∀a ∈ A κ (a) =   1 0 λ 0 µ (a) σ a 0 0 σ   , ζ =   0 0 ν   .
It is routine to verify, by induction on the length of u, and based on Lemma 3, that λ • µ (P (u))

• ν = η • κ (u) • ζ for every u in A * . Let now ξ =   1 0 ν
  and let s be the series realised by (η, κ, ξ):

∀u ∈ A * <s, u> = 1 + card ({v ∈ A * | v ∈ L and v ≺ u}) .
5 cf. Remark 6 below.

In order to get the enumerating series of L, we must retain the words that belong to L only, that is, to make the Hadamard product with the characteristic series L of L:

E S = s ⊙ L ,
and E S is N-rational as the Hadamard product of two N-rational series (this is often referred to as (another) Schützenberger Theorem6 ).

Remark 2. The construction underlying the proof yields for E L an N-representation of dimension 2 k 2 + k.

Computation of the value of a word

The description of a N-rational series s by a N-representation gives a way to compute the coefficient of any word w in s. Theorem 1 thus solves ipso facto the problem of computing the value π L (w) of a word w in a rational abstract number system L (which occupies the whole Sect. 2 in [START_REF] Lecomte | Abstract numeration systems[END_REF]).

If s has a representation (χ, ω, φ) of dimension n, and if w is of length ℓ, the general algorithm consists in computing χ

• ω(w i+1 ) = (χ • ω(w i )) • ω(a i+1 ) for i = 0 to i = ℓ -1
, where a i is the i-th letter of w and w i its prefix of length i. Every step costs 2 n 2 operations, thus in total, roughly 2 ℓ n 2 operations.

It would be not such a good idea, however, to apply this general algorithm to the representation of dimension 2k 2 +k we have obtained in the proof of Theorem 1 above. Its particular form allows, in fact, to compute with vectors and matrix of dimension k only.

Given as above the unambiguous automaton A of dimension k which recognises L and the corresponding N-representation (λ, µ, ν), we associate a pair α(w), γ(w) with every w in A * , where α(w) and γ(w) are two (row) vectors of dimension k, α(w) with entries in {0, 1}, γ(w) with entries in N. The pair α(w), γ(w) is computed by induction on the length of w in the following way. Let ℓ be the length of w, let α(1 A * ) = λ , β(1 A * ) = λ , and γ(1 A * ) = 0 , and, for every 0 i < ℓ, let

α(w i+1 ) = α(w i ) • µ (a i+1 ) , β(w i+1 ) = α(w i ) • σ a i+1 , and γ(w i+1 ) = λ + β(w i+1 ) + γ(w i ) • σ .
All α(w) have entries in {0, 1} since A is unambiguous. As a simple reformulation of the preceding subsection, we have π L (w) = γ(w)•ν if α(w)•ν = 1 , that is, if w is recognised by A and thus in L, π L (w) undefined otherwise. This algorithm, that is, the computation of α(w), γ(w) , costs roughly 6 ℓ k 2 operations.

Example 3. Let us consider again the language L 1 and the DFA A 1 of Figure 1. We have thus:

σ a = 0 0 0 0 , σ b = 1 0 0 1 , σ = 1 1 1 1 .
The computation of π L 1 (bbabb) for instance takes the following steps.

i a i α i β i γ i 0 (1, 0) (1, 0) (0, 0) 1 b (0, 1) (1, 0) (2, 0) 2 b (1, 0) (0, 1) (3, 3) i a i α i β i γ i 3 a (1, 0) (0, 0) (7, 6) 4 b (0, 1) (1, 0) (15, 13) 5 b (1, 0) (0, 1) (29, 29) And finally π L 1 (bbabb) = (29, 29) • ν 1 = 29.
Remark 3. The computation of α(w), γ(w) is very similar to the construction called product of an automaton by a skew action in [START_REF] Sakarovitch | On the decomposition of k-valued rational relations[END_REF][START_REF] De Souza | Etude structurelle des transducteurs de norme bornée[END_REF].

Representation of recognisable subsets of numbers

If s is an N-rational series, that is, a map s : A * → N , it is well known that for any recognisable set of numbers X, s -1 (X) is a rational set of A * (see [3, Corol. III.2.4], [8, Th. VI.10.1] or [17, Corol. III,4,21], for instance). Theorem 1 thus directly implies the following statement, which has also been proved without reference to it in [START_REF] Lecomte | Numeration systems on a regular language[END_REF] and in [START_REF] Krieger | Decimation of languages and state complexity[END_REF].

Corollary 4 ([14]

). A recognisable set of numbers is L-recognisable in any rational abstract numeration system L.

If Corollary 4 requires formally no proof after the characterisation of rational abstract numeration systems given by Theorem 1, it is interesting to further investigate the construction which, given L and a recognisable set of numbers X computes an automaton which recognises the set X L . The computation method used in the preceding section (which is not the mere application of the general result that yields Corollary 4) allows to establish easily the following statement. Proposition 5. Let L be a rational language over A * recognised by a deterministic automaton of dimension k. For any integers p and r < p, let X p,r = pN + r be the set of integers congruent to r modulo p. Then the language X p,r L of representations of numbers in X p,r is recognised by a deterministic automaton with at most k p k states.

Proof. Let A be an automaton, with set of states Q of cardinal k, which recognises L and (λ, µ, ν) the corresponding N-representation. If A is deterministic, then λ and µ are row monomial and so are all α(w), for w in A * , which are thus in 1-1 correspondence with the elements of Q.

Let C be the automaton whose set of states is

R = { α(w), δ(w) | w ∈ A * }
where δ(w) = γ(w) mod p .

Thus, R ⊆ Q × (Z/pZ) k . The transitions of C are defined by, for every a in A:

∀w ∈ A * , ∀a ∈ A α(w), δ(w) a --→ C α(wa), δ(wa) .
The initial state of C is (λ, 0) and its final states are those α(w), δ(w) where α(w) is final in A and δ(w) • ν = r mod p . It then follows that the language accepted by C is X p,r L .

Remark 4. If we start from an unambiguous automaton A of dimension k, the same method yields a deterministic automaton C with at most 2 k p k states.

Example 4. The automaton built in this way from A 1 and for the recognisable set of numbers 3 N + 1 is the automaton C 1 shown at Figure 2 (this automaton is not minimal; its minimal quotient has only 8 states). Remark 5. In [START_REF] Lecomte | Abstract numeration systems[END_REF], another construction has been given for the same purpose.

The automaton D built with this other method and which recognises X p,r L is not deterministic, but codeterministic and has, roughly, k p k+1 states. Since D is codeterministic, its determinisation yields the minimal automaton of X p,r L and thus, thanks to Proposition 5, does not produce an exponential blow-up. We do not know of a direct proof of this fact.

Remark 6. After the submitted version was written (and sent), we have learned of the reference [START_REF] Krieger | Decimation of languages and state complexity[END_REF]. Not only Corollary 4 is established there, but with a method of proof which is very similar to ours. Our Lemma 3 is Lemma 1 in [START_REF] Krieger | Decimation of languages and state complexity[END_REF]. The term representation is not used there but the matrices µ(a), σ a and σ are defined (under other notation) and used to give the same proof of Equation (2) (Lemma 2 in [START_REF] Krieger | Decimation of languages and state complexity[END_REF]). Afterwards, [START_REF] Krieger | Decimation of languages and state complexity[END_REF] develops in another direction than this paper: it proves lower bounds for the state complexity of X p,r L and shows that the property corresponding to Corollary 4 does not hold for context-free languages.

Remark 7. If the numeration system considered is a positional numeration system (and still a rational one), and under some supplementary hypotheses, then the exact number of states for the minimal automaton of X p,r L can be computed (cf. [START_REF] Charlier | State complexity of testing divisibility[END_REF]).

Proof of Theorem 2

The image of a rational language by a rational relation (or transduction) is a rational language; this classical result, due to Nivat and called Evaluation Theorem in [START_REF] Eilenberg | Automata, Languages and Machines[END_REF], extends to rational series, as we state now (cf. [START_REF] Sakarovitch | Eléments de théorie des automates[END_REF]). if it is defined, is a K-rational series over B * .

It is this result that was used in [START_REF] Choffrut | Rational transductions and complexity of counting problems[END_REF] for the proof of Theorem 1.

Proof of Theorem 2. Let s in NRat A * and L = supp s in Rat A * . The set L is totally ordered by the radix order L = {w 0 < w 1 < w 2 < . . . < w n < . . .} and Succ L is the function from A * into itself whose domain is L and which maps every w i to w i+1 . It is well-known that Succ L is a rational function ( [START_REF] Berthé | On the cost and complexity of the successor function[END_REF][START_REF] Frougny | Number representation and finite automata[END_REF]) and hence unambiguous ( [START_REF] Eilenberg | Automata, Languages and Machines[END_REF][START_REF] Sakarovitch | Eléments de théorie des automates[END_REF]). It then follows that the series Succ L (s) = i=∞ i=0 <s, w i > w i+1 is an N-rational series and t = s -Succ L (s) is a Z-rational series. Now, s is the enumerating series E L of the abstract numeration system L if, and only if, for every positive integer i, <t, w i > = 1 , that is, if, and only if, t -L \ {w 0 } = 0 , a condition which is known to be decidable as Z is a sub(semi)ring of a field (cf. [START_REF] Eilenberg | Automata, Languages and Machines[END_REF][START_REF] Sakarovitch | Eléments de théorie des automates[END_REF]).

Problems and future work

Looking at abstract number systems as N-rational series naturally leads to two families of questions. The first family consists in questions on N-rational series which ask to which extent the series is related to abstract number systems; the second in questions which generalise to N-rational series questions that are usually considered for (abstract) numeration systems.

An example of questions in the first family is to ask if it is decidable whether a given N-rational series is the enumeration (in a radix ordering) of its (rational) support in a certain, and unknown, abstract numeration system. This seems to be rather a difficult problem. An obvious necessary condition for a series to be a positive instance of this problem is itself a non trivial problem that can be formulated in the following way.

Example 1 .

 1 Let A = {a, b} , with a < b and let L 1 be the language of words with an even number of b's: L 1 = {w ∈ {a, b} * | |w| b ≡ 0 mod 2} . The sequence of the representations of the integers is: {ε, a, aa, bb, aaa, abb, bab, . . .} and, for instance, 18 L 1 = a a b a b and π L 1 (b b a b b) = 29 .

Figure 1 :

 1 Figure 1: A DFA accepting words with an even number of b's

Figure 2 :

 2 Figure 2: A DFA recognising the set 3N + 1 in the ANS L 1 .

Proposition 6 .

 6 Let ϕ : A * → B * be an unambiguous rational relation and s a K-rational series over A * . Then the seriesϕ(s) = w∈A * <s, w> ϕ(w) = u∈B * <s, ϕ -1 (u)> u ,(4)

The definition of radix order will be given below.

Notice that ≺ is not reflexive and is not the order but the strict part of the radix order.

This is true only because A * is a free monoid.

cf. Remark 6 below.

cf. [3, Th. I.5.3], [8, Prop. VI.7.1] or [17, Cor. III.3.9].

Acknowledgements

The authors are grateful to Sylvain Lombardy for the discussions with him on this problem. They also have greatly benefited from having access to the chapter of Lecomte and Rigo on abstract numeration systems in the book Combinatorics, Automata and Number Theory to be published by Cambridge University Press.

The authors are particularly indebted to Michel Rigo who, after seeing a first version of this paper, gave them the information that Theorem 1 was already known, the pertinent references, and the encouragements to nevertheless write down a complete version.

Conjecture 7. It is decidable whether an N-rational series is a monotone increasing function (for a given order of letters).

A result due to Honkala [START_REF] Honkala | A decision method for the recognizability of sets defined by number systems[END_REF] provides a kind of converse of Corollary 4 in the case of p-ary numeration systems and states that it decidable whether a precognisable set of numbers is recognisable. The generalisation of this result to larger class of numeration systems has been recently studied in [START_REF] Bell | A decision problem for ultimately periodic sets in non-standard numeration systems[END_REF][START_REF] Charlier | Abstract numeration systems: recognizability, decidability, multidimensional S-automatic sequences, and real numbers[END_REF]. Its generalisation to abstract number systems has been stated as a problem in [START_REF] Honkala | A note on decidability questions related to abstract numeration systems[END_REF]. It is also a typical example of a question in the second family.

Conjecture 8. It is decidable whether the set of coefficients of an N-rational series is a recognisable set of numbers.

Summary

In this short paper, we have presented a new idea for the study of abstract number systems, which brings to the subject the whole power of weighted automata theory. In return, the subject of abstract number systems naturally opens new questions for the theory of N-rational series.