Benoît Sijobert
email: benoit.sijobert@inria.fr

Charles Fattal
email: cfattal@lachataigneraie.fr

Anne Daubigney
email: anne.daubigney@gmail.com

Christine Azevedo Coste

Christine Azevedo-Coste

Christine Azevedo
email: christine.azevedo@inria.fr

Participation to the first Cybathlon: an overview of the FREEWHEELS team FES-cycling solution

Keywords: FES cycling, spinal cord injuries, Cybathlon

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.

In 2016 an academic research team (CAMIN, INRIA, Montpellier, France) and a rehabilitation center team (COS DIVIO, Dijon, France) have participated into Cybathlon competition, in the FES-bike discipline under the name of FreeWheels team. [START_REF]FREEwheels official website[END_REF] Despite more than 14 years investigating on Functional Electrical Stimulation (FES) applications and more broadly in the field of Neuroprostheses, [START_REF] Zhang | Evoked electromyography-based closed-loop torque control in functional electrical stimulation[END_REF][START_REF] Guiraud | Implanted functional electrical stimulation: case report of a paraplegic patient with complete SCI after 9 years[END_REF][START_REF] Coste | Continuous gait cycle index estimation for electrical stimulation assisted foot drop correction[END_REF][START_REF] Bó | On the Use of Fixed-Intensity Functional Electrical Stimulation for Attenuating Essential Tremor[END_REF][START_REF] Sijobert | Effects of sensitive electrical stimulation-based somatosensory cueing in Parkinson's disease gait and Freezing of Gait assessment[END_REF] using FES for cycling was a first for the research team. The FES-bike discipline consisted of a series of races between two spinal cord injured (SCI) participants propelling a cycling device by means of their lower limbs muscular contractions elicited by FES. The participants had to cover a 750 meters distance in less than 8 minutes, starting from a 10° descending ramp. Mobile and stationary cycling using FES have been widely investigated since the 1980s through many studies. [START_REF] Glaser | Physiologic aspects of spinal cord injury and functional neuromuscular stimulation[END_REF][START_REF] Fornusek | Cardiovascular and Metabolic Responses During Functional Electric Stimulation Cycling at Different Cadences[END_REF][START_REF] Berkelmans | Fes cycling[END_REF][START_REF] Newham | FES cycling[END_REF][START_REF] Wilder | Functional electrical stimulation cycle ergometer exercise for spinal cord injured patients[END_REF][START_REF] Guimarães | Towards Parameters and Protocols to Recommend FES-Cycling in Cases of Paraplegia: A Preliminary Report[END_REF][START_REF] Coste | FES in Europe and Beyond: Current Translational Research[END_REF] It has been shown that FES-cycling of subjects with SCI results in physiological and psychological positive effects such as cardiovascular training, decrease in pressure sores occurrence and selfesteem improvements. [START_REF] Figoni | Acute hemodynamic responses of spinal cord injured individuals to functional neuromuscular stimulation-induced knee extension exercise[END_REF][START_REF] Pl | Modes, benefits, and risks of voluntary electrically induced exercise in persons with spinal cord injury[END_REF][START_REF] Peri | Can FES-Augmented Active Cycling Training Improve Locomotion in Post-Acute Elderly Stroke Patients?[END_REF] However, the use of this technology has often remained restricted to indoor and stationary ergometers in clinical contexts, partly due to the small amount (10-25 W) of power produced, [START_REF] Eser | Influence of different stimulation frequencies on power output and fatigue during FES-cycling in recently injured SCI people[END_REF][START_REF] Pollack | Aerobic training effects of electrically induced lower extremity exercises in spinal cord injured people[END_REF] the requirement of experimented users and the lack of commonly available affordable FES outdoor bikes. Cybathlon FES-bike challenge motivated considerable efforts in technical developments to allow the coverage of a significant distance at a suitable speed. This event promoted the activity as recreational with potential mobility applications. In the following we describe the different characteristics of our pilot and device.

Material and Methods

This study was conducted in accordance with the principles of the Declaration of Helsinki and approved by the local Ethics Committee (CPP Dijon, ID RCB: 2016-A00279-42), including the participant selection, physical preparation, training and participation into Cybathlon. One subject was included, the so-called pilot.

Inclusion and exclusion criteria

The Cybathlon competition rules constrained the participant to present an ASIA A or B spinal cord injury, with a complete loss of motor function. We have furthermore defined the following inclusion criteria: age ≥ 18 years and < 65 years, complete traumatic lesion > 12 months, neurological level T2 to T12, stable medical status, stimulable sublesional muscles and passive joint movement in the lower limbs. Exclusion criteria were defined as follows: body mass index ≥ 30, pressure ulcer, thrombophlebitis, neurogenic paraosteoarthropathy, cardiovascular disease, muscle disease, DXA T-score <-2.5, hip or knee arthroplasty, epilepsy, hypotension, lower limb fracture within the past 12 months, a pacemaker or other implant, and pregnancy. No betablocker treatment was authorized during the experiment.

Pilot profile

The pilot was selected in respect with the inclusion criteria listed in the previous paragraph. He was 47 years -266 - In order to later facilitate diffusion of the approach we also wanted the device to be affordable, foldable and sufficiently light for transportation. We have selected the ICE Adventure 26 . This product has a mesh seat, designed to provide an optimum back support and comfort. The tension in the back could be adjusted. The cover features is made on a breathable fabrics to maximise airflow. A VICAIR®-type pressure relief cushion was chosen by the team occupational and physiotherapists and then fixed on the seat to avoid sliding. Hase® pedals with calf support were adapted to ensure a sufficient lateral locking of the legs while holding the ankle joint at 90°. The ICE trike was initially equipped with a 26" rear wheel, we changed for a 24" one to tilt the bike and elevate the crank relatively to the pilot's position (Fig. 2). Low rolling resistance tires (Schwalbe® Kojak 24 x 1.35") were also mounted instead of the original tires.

A few weeks before Cybathlon competition we decided to modify the rear wheel transmission, initially a freewheel, to a fixed-gear. This decision was made in consequences of difficulties to pass the "dead points" during cycling phase (0° and 180°, full extension/full flexion). Our pilot could not succeed to complete the pedaling cycle once fatigue appeared and sometimes from the beginning of training. Furthermore, the cocontraction balance between hamstrings and quadriceps muscle activations has a critical influence at these particular dead points and with time the equilibrium initially tuned is becoming less optimal (fatigue, movements on the seat…). The fixed-gear wheel was also supported by Szecsi 2007 et al. [START_REF] Szecsi | Functional output improvement in FES cycling by means of forced smooth pedaling[END_REF] results, where the forces smoothing introduced by the use of a fixed-gear improved endurance. The gear ratio was set to be equal to 1.22 (number of rear wheel spins for one complete crank spin, or ratio between the chainring and the rear sprocket). Calculated with the 24" rear wheel, we chose this ratio in order to best match between the optimal literature cadence, [START_REF] Szecsi | Functional output improvement in FES cycling by means of forced smooth pedaling[END_REF] and the speed needed to reach 750 meters in 8 minutes, i.e. about 47 RPM (5.6km/h). The classical approach of FES-cycling is to predefine muscular activation patterns regarding the crank angle needed for the timing of the stimulation. The Berkelbike Pro® is equipped with such a crank sensor that can be connected to a stimulator commercialized by the company. A similar crank angle encoder (Baumer MDFK08) has been adapted on the crank axis of the FreeWheels trike in order to provide pedaling angle as an input to the stimulator. In addition, we have instrumented the device with different kind of sensors used to monitor the performances. Two speedometers were set, one with the screen oriented towards the pilot and a second one oriented to be visible to a person walking aside the trike. Finally, we added an ANT+ cardiac strap monitor to follow heart rate in response to effort.

Stimulation Parameters

During this study, we used the Berkelbike® FES-Box stimulator, an 8-channels battery powered stimulator. In cycling mode, a pre-programmed pattern based on crank angle, triggers the stimulation on and off cycles. Electrical stimulation was delivered on each leg via three stimulation channels composed each one of two platinum stimulating electrodes, size 4.5x10cm, respectively located along the rectus femori/vastus lateralis muscles, the vastus medialis and the hamstrings (Fig. 3). A musculocutaneous flap was taken from the gluteus maximus as a treatment of a sacral pressure ulcer which prevented us to consider this muscle group as a candidate to electrical stimulation. Stimulation pattern was designed with individualized quadriceps.

On each leg, the first channel dedicated to the rectus femori/vastus lateralis muscles was set to be effective during the pushing phase (full flexion to vertical crank), meanwhile the second channel was meant to later activate the vastus medialis, from the time where the

Training information

A one year (12 mo.) training program was designed, divided into two stages:

1. from months M1 to M6: home based, two to three times a week, stimulation program of the lower limb muscles (quadriceps / hamstrings) during 30 minutes of isometric training using CEFAR® Physio 4 stimulator and the following stimulation parameters: rectangular charge-balanced biphasic pulse stimulation with a pulsewidth of 300µS per phase and a ramp up, sent at 30Hz during 10s followed by 3s of rest. At the start of each month, a mapping determined the intensity of stimulation needed to obtain contraction of the electrostimulated muscle and the intensity needed to reach a 4/5 score on the MRC (Medical Research Council) scale.

from M7 to M11: FES-cycling training program, two

to three times a week, on an instrumented trike. Either set on a home trainer or rolling outdoor and indoor on a flat surface. Stimulations parameters were modulated on the basis of progress.

From the beginning of the training, psychological and functional variables were also continuously monitored (self esteem, life quality, etc… ).

Results

Performances FES-bike race was divided into two main steps, a qualification race and a final race. Our pilot was successfully qualified after a first race (670m in 430s, Fig. 4) and could take part in the final race. The average speed during the final race was of 5.71 km/h where he reached a maximal speed of 6.26km/h. He finished 6th out of 12, accomplishing the whole race in 467 sec. Our cyclist reached his best performance one week before Cybathlon when he succeeded to achieve a distance of 1820m on the ICE trike in stationary mode and 760m rolling outdoor on a flat surface (Fig. 5). Overcoming some technical issues, the functional performances dramatically and continuously increased from M11 (Fig. 5). Speed and endurance goals were successfully reached at M12 and kept evolving at M13 (1080m in 14min, Fig. 6).

Feasibility

The total cost of the device including the tricycle and the stimulator is estimated to 7650€ (4500€ for the stimulator, 3000 euros for the tricycle, and the rest for mechanical adaptations). The solution has the advantage to be easily duplicated. The choices made on its adaptation were driven by the Cybathlon objective: to cover 750m in less than 8min. The system should be adapted and optimized to allow better performances based on the acquired knowledge. The stimulation patterns are designed through an individual manual tuning, with no general rule, making the procedure, at this stage, hardly feasible by an uninitiated user outside of a clinical environment

Discussion

Our participation into this human adventure was a team success. Despite numerous practical issues the initial goal was reached. We have not yet been able to objectivize and quantify all the technical choices and their consequences made to reach the Cybathlon goals.

During the development phase we equipped the tricycle with inertial sensors (one on the crank and a second one on the rear wheel) and pedal power meters (Powertap P1, Fig. 7), using ANT+ wireless communication protocol and a specifically developed software. The original idea was to modulate stimulation parameters and patterns regarding pedaling cycle and produced forces. However, the feature which was supposed to give us a real time access to the stimulator parameters had not been ready on time before Cybathlon competition. Moreover, the pedal power meters were not sensitive enough to ensure an accurate measurement of the weak produced forces and to be usable as input to close loop control. Meanwhile, the different figures monitoring speed during training and races (Fig. 4 and Fig. 6) show the manual control of intensity level regarding cycling speed was really close to the initially automated considered one (consign of average speed set at 5.6 km/h, average recorded speed during final race of 5.71 km/h). The covering surface of the track also represented an important factor which drastically affected the obtained performances all along training and race sessions and could not have been quantified. The pilot training was part of a research protocol in which a longitudinal observation was achieved. The aim was to follow the evolution of performances as well as physiological and psychological parameters. Our pilot demonstrated a global acceptance of the participation into this 1 year protocol and an improvement along the experience of psychological variables such as life quality and self-esteem could be observed. Taking part in a competition as the final aim of the training program demonstrated a real benefit on training constraints tolerance. [START_REF] Coste | Comparison of strategies and performance of functional electrical stimulation cycling in spinal cord injury pilots for competition in the first ever Cybathlon[END_REF][START_REF] Berkelmans | Strategies and performances of Functional Electrical Stimulation Cycling using the BerkelBike with Spinal Cord Injury in a competition context (CYBATHLON)[END_REF][START_REF] Laubacher | The road to Cybathlon 2016 -Functional electrical stimulation cycling Team IRPT/SPZ[END_REF][START_REF] Guimarães | FES Bike Race preparation to Cybathlon 2016 by EMA team: a short case report[END_REF][START_REF] Metani | Functional electrical stimulation cycling strategies tested during preparation for the First Cybathlon Competition -a practical report from team ENS de Lyon[END_REF][START_REF] Mcdaniel | Cycle Training Using Implanted Neural Prostheses: Team Cleveland[END_REF][START_REF] Arnin | Development of practical functional electrical stimulation cycling systems based on an electromyography study of the Cybathlon 2016[END_REF][START_REF] Leung | The Effectiveness of Functional Electrical Stimulation (FES) in On-Off Mode for enhancing the cycling performance of Team Phoenix at 2016 Cybathlon[END_REF] This experience and the interactions with other teams breaks new grounds to innovative and generic FES-cycling based training programs for an important range of SCI subjects. However, the constraints and the reality behind people with SCI daily lives are far from the Cybathlon FES bike context and require numerous improvements such as further investigations to enable a safe recreational use in complete autonomy of this activity.

Fig 1 .Fig 2 .Fig 3 .

 123 Fig 1. The two trikes used (left) during the initial training program, the Berkelbike Pro® and (right) for the final training program and at the competition, the FreeWheels trike old (1m78, 74kg) at the time of inclusion and presented a spinal cord lesion T3 (AIS A) since 1995.Tricycle descriptionTwo distinct cycling devices were used during the study (Fig.1). A the Berkelbike Pro®, a commercial tricycle designed for persons with physical limitations. B the FreeWheels trike, an instrumented recumbent cycle adapted from a commercial device, the ICE

Fig 4 .

 4 Fig 4. Average speed (km/h) computed on each lap (from L1 to L10) during the Cybathlon qualification (top) and final (bottom) races

Fig 5 .Fig 6 .

 56 Fig 5. Maximum covered distance during the training sessions, from M9 to M13 (from week 11th to week 29th of the FES-cycling program)

Fig 7 .

 7 Fig 7. Power is computed via force pedals, cadence and speed are monitored via two inertial sensors (on the crank and rear wheel)

Acknowledgments

This research was supported by a Neuroglia foundation for neurosciences grant and an INRIA internal financial support. Part of this work has been done in collaboration with EMA team (University of Brasilia) with the help of CACAO associate team (INRIA)

List of acronyms

Author's contributions

BS is the principle author and responsible for the technical integration and software design of the FreeWheels project. CF and AD were responsible for the clinical protocol, pilot training and follow up. CA was the leader of the team. All the authors read and approved the final version of the manuscript.

Conflict of Interest

The authors declare no potential conflict of interests.


Olivier Collier

Arnak S Dalalyan

Estimating linear functionals of a sparse family of Poisson means

Keywords: MSC 2010 subject classifications: Primary 62J05; secondary 62G05 Nonasymptotic minimax estimation, linear functional, group-sparsity, thresholding, Poisson processes

Assume that we observe a sample of size n composed of pdimensional signals, each signal having independent entries drawn from a scaled Poisson distribution with an unknown intensity. We are interested in estimating the sum of the n unknown intensity vectors, under the assumption that most of them coincide with a given "background" signal. The number s of p-dimensional signals different from the background signal plays the role of sparsity and the goal is to leverage this sparsity assumption in order to improve the quality of estimation as compared to the naive estimator that computes the sum of the observed signals.

We first introduce the group hard thresholding estimator and analyze its mean squared error measured by the squared Euclidean norm. We establish a nonasymptotic upper bound showing that the risk is at most of the order of σ 2 (sp + s 2 √ p) log 3/2 (np). We then establish lower bounds on the minimax risk over a properly defined class of collections of s-sparse signals. These lower bounds match with the upper bound, up to logarithmic terms, when the dimension p is fixed or of larger order than s 2 . In the case where the dimension p increases but remains of smaller order than s 2 , our results show a gap between the lower and the upper bounds, which can be up to order √ p.

INTRODUCTION AND PROBLEM FORMULATION

Let µ = (µ 1 , . . . , µ p ) be any vector from R p with positive entries. In what follows, we denote by P p (µ) the distribution of a random vector X with independent entries X j drawn from a Poisson distribution P(µ j ), for j = 1, . . . , p. We consider that the available data is composed Modal'X, UPL, Univ Paris Nanterre, F92000 Nanterre France. 3 avenue Pierre Larousse, 92245 Malakoff, France.

of n independent vectors X 1 , . . . , X n randomly drawn from R p such that

X i iid ∼ σ 2 P p (σ -2 µ i ), i = 1, . . . , n,
where σ > 0 is a known parameter that can be interpreted as the noise level. The use of this term is justified by the fact that when σ goes to zero, we have σ 2 P p (σ -2 µ) ≈ N p µ, σ 2 diag(µ) . We will work under the assumption that most intensities are known and equal to a given vector µ 0 , which can be thought of as a background signal. We denote by X = [X 1 , . . . , X n ] the p × n matrix obtained by concatenating the vectors X i for i = 1, . . . , n.

Assumption (S): for a given vector µ 0 ∈ (0, ∞) p , the set S = {i : µ i = µ 0 } is a very small subset of [n] := {1, . . . , n}.

The unknown parameter in this problem is the p × n matrix M = [µ 1 , . . . , µ n ] obtained by concatenating the signal vectors µ i . However, we will not aim at estimating the matrix M. Instead, we aim at estimating a linear functional of the intensities:

L(M) = i∈[n] (µ i -µ 0 ) = i∈S (µ i -µ 0 ) = M1 n -nµ 0 .
This functional has a clear meaning: it is the superposition of all the signals contained in X i 's beyond the background signal µ 0 . Notice that if we knew the set S, it would be possible to use the oracle

L S = i∈S (X i -µ 0 ), (1)
which would lead to the risk

E L S -L(M) 2 2 = i∈S j∈[p] Var(X ij ) = σ 2 i∈S µ i 1 p ≤ µ ∞ σ 2 sp,
where µ ∞ is an upper bound on max i µ i ∞ and s is the cardinality of S. On the other hand, one can use the naive estimator

L [n] = i∈[n] (X i -µ 0 ) = L(X)
which does not use at all the sparsity assumption (S). It has a quadratic risk given by

E L [n] -L(M) 2 2 = i∈[n] j∈[p] Var(X ij ) = σ 2 i∈[n] µ i 1 p ≤ µ ∞ σ 2 np.
One can remark that the rate σ 2 sp we have got for the oracle is independent of n, which obviously can not continue to be true when S is unknown and should somehow be inferred from the data. An important statistical question is thus the following.

Question (Q1): What kind of quantitative improvement can we get upon the rate σ 2 np of the naive estimator by leveraging the sparsity assumption?

In the framework considered in the present work, for answering Question (Q1), we are allowed to use any estimator of L(M) which requires only the knowledge of σ and µ 0 . However, in practice, it is common to give advantage to estimators having small computational complexity. Therefore, another question to which we will give a partial answer in this work is:

Question (Q2): What kind of quantitative improvement can we get upon the rate σ 2 np of the naive estimator by leveraging the sparsity assumption and by constraining ourselves to estimators that are computable in polynomial time?

In the statement of the last question, the term "polynomial time" refers to the running time as a function of s, n and p.

Relation to previous work

Various aspects of the problem of estimation of a linear functional of an unknown highdimensional or even infinite-dimensional parameter were studied in the literature, mostly focusing on the case of a functional taking real values (as opposed to the vector valued functional considered in the present work). Early results for smooth functionals were obtained by [START_REF] Koshevnik | On a non-parametric analogue of the information matrix[END_REF]. Minimax estimation of linear functionals over various classes and models were thoroughly analyzed by [START_REF] Butucea | Adaptive estimation of linear functionals in the convolution model and applications[END_REF]; Cai andLow (2004, 2005); [START_REF] Donoho | On Minimax Estimation of Linear Functionals[END_REF]; [START_REF] Efromovich | Adaptive estimates of linear functionals[END_REF]; [START_REF] Golubev | An oracle approach to adaptive estimation of linear functionals in a gaussian model[END_REF]; [START_REF] Juditsky | Nonparametric estimation by convex programming[END_REF]; [START_REF] Klemela | Sharp adaptive estimation of linear functionals[END_REF]; [START_REF] Laurent | Adaptive estimation of linear functionals by model selection[END_REF]. There is also a vast literature on studying the problem of estimating quadratic functionals [START_REF] Bickel | Estimating integrated squared density derivatives: Sharp best order of convergence estimates[END_REF][START_REF] Cai | Optimal adaptive estimation of a quadratic functional[END_REF][START_REF] Donoho | Minimax quadratic estimation of a quadratic functional[END_REF][START_REF] Laurent | Adaptive estimation of a quadratic functional by model selection[END_REF]. Since the estimators of (quadratic) functionals can be often used as test statistics, the problem of estimating functionals has close relations with the problem of testing that were successfully exploited in [START_REF] Collier | Curve registration by nonparametric goodness-of-fit testing[END_REF]Comminges andDalalyan, 2012, 2013;[START_REF] Lepski | On estimation of the Lr norm of a regression function[END_REF]. The problem of estimation of nonsmooth functionals was also tackled in the literature, see [START_REF] Cai | Testing composite hypotheses, hermite polynomials and optimal estimation of a nonsmooth functional[END_REF].

Most papers cited above deal with the Gaussian white noise model, Gaussian sequence model, or the model in which the observations are iid with a common (unknown) density function. Surprisingly, the results on estimation of functionals of the intensity of a Poisson process are very scarce. (We are only aware of the paper [START_REF] Kutoyants | Estimation of linear functionals of poisson processes[END_REF], which in an asymptotic set-up established the asymptotic normality and efficiency of a natural plug-in estimator of a real valued smooth linear functional.) This is surprising, given the relevance of the Poisson distribution for modeling many real data sets arising, for instance, in image processing and biology. Furthermore, there are by now several comprehensive studies of nonparametric estimation and testing for Poisson processes, see [START_REF] Birgé | Model selection for Poisson processes[END_REF]; [START_REF] Ingster | Nonparametric hypothesis testing for intensity of the poisson process[END_REF]; [START_REF] Kutoyants | Statistical Inference for Spatial Poisson Processes[END_REF]; Reynaud-Bouret and Rivoirard (2010) and the references therein. In addition, several statistical experiments related to Poisson processes were proved to be asymptotically equivalent to the Gaussian white noise or other important statistical experiments [START_REF] Brown | Equivalence theory for density estimation, poisson processes and gaussian white noise with drift[END_REF][START_REF] Grama | Asymptotic equivalence for nonparametric generalized linear models[END_REF]. However, the models related to the Poisson distribution demonstrate some notable differences as compared to the Gaussian models; the most simple differences are the heteroscedasticity of the Poisson distribution1

and the fact that it has heavier tails than the Gaussian distribution. The impact of these differences on the statistical problems is amplified in high dimensional settings, as it can be seen in [START_REF] Ivanoff | Adaptive lasso and group-lasso for functional poisson regression[END_REF][START_REF] Jia | The lasso under poisson-like heteroscedasticity[END_REF].

To complete this section, let us note that the investigation of the statistical problems related to functionals of high-dimensional parameters under various types of sparsity constraints was recently carried out in several papers. The case of real valued linear and quadratic functionals was studied by [START_REF] Collier | Minimax estimation of linear and quadratic functionals on sparsity classes[END_REF] and [START_REF] Collier | Optimal adaptive estimation of linear functionals under sparsity[END_REF], focusing on the Gaussian sequence model. [START_REF] Verzelen | Adaptive estimation of High-Dimensional Signal-to-Noise Ratios[END_REF] analyzed the problem of the signal-to-noise ratio estimation in the linear regression model under various assumptions on the design. In a companion paper of the present submission, [START_REF] Collier | Rate-optimal estimation of p-dimensional linear functionals in gaussian models[END_REF] considered the problem of a vector valued linear functional estimation in the Gaussian sequence model. The relations with this latter work will be discussed in more details in Section 4 below.

Agenda

The rest of the paper is organized as follows. In Section 2, we introduce an estimator of the functional L(M) based on thresholding the columns of X with small Euclidean norm. The same section contains the statement and the proof of a nonasymptotic upper bound on the expected risk of the aforementioned estimator. Lower bounds on the risk are presented in Section 3. Some avenues for further research and a summary of the contributions of the present work are given in Section 4. The proofs of the technical lemmas are gathered in Section 5.

Notation

We use uppercase boldface letters for matrices and uppercase italic boldface letters for vectors.

For any vector v and for every q ∈ (0, ∞), we use the standard definitions of the q -norms v q q = j |v j | q , along with v ∞ = max j |v j |. For every n ∈ N, 1 n (resp. 0 n ) is the vector in R n all the entries of which are equal to 1 (resp. 0). For a matrix M, we denote by M the spectral norm of M (equal to its largest singular value), and by M F its Frobenius norm (equal to the Euclidean norm of the vector composed of its singular values). For any matrix M, we denote by L(M) the linear functional defined as the sum of the columns of M.

GROUP HARD THRESHOLDING ESTIMATOR

This section describes the group hard thresholding estimator, that is computationally tractable and achieves the optimal convergence rate in a large variety of cases. The underlying idea is rather simple: the Euclidean distance between each signal X i and the background signal µ 0 is compared to a threshold; the sum of all the signals for which the distance exceeds the threshold is the retained estimator.

Thus, the group hard thresholding method analysed in this section can be seen as a two-step estimator, that first estimates the set S and then substitutes the estimator of S in the oracle L S considered in (1). To perform the first step, we choose a threshold λ > 0 and define

S(λ) = i ∈ [n] : X i -µ 0 2 ≥ σ( X i 1 + λ) 1/2 .
(2)

We will refer to the plug-in estimator L S(λ) as the group hard thresholding estimator, and will denote it by

L GHT = i∈ S(λ) (X i -µ 0 ). ( 3
)
It is clear that the choice of the threshold λ has a strong impact on the statistical accuracy of the estimator L GHT . One can already note that the threshold to which the distance X i -µ 0 2 is compared depends on the signal X i ; this is due to the heteroscedasticity of the models based on observations drawn from the Poisson distribution.

Prior to describing in a quantitative way the impact of the threshold on the expected risk, let us emphasize that the computation of L GHT requires O(np) operations; so it can be done in polynomial time. Indeed, the computation of each of the n norms X i -µ 0 2 requires O(p) operations, so S(λ) can be computed using O(np) operations. The same is true for the sum in the right hand side of (3).

The following theorem establishes the performance of the previously defined thresholding estimator, for a suitably chosen threshold λ.

Theorem 1. Suppose σ 3 µ -3/2 ∞ ≤ 40 log 3/2 (2np) ≤ 0.9σ -3 µ 3/2 ∞ . The group hard thresholding estimator (2)-(3) with λ = 40 µ 0 ∞ √ p log 3/2 (2np), (4)
has a mean squared error bounded as follows:

sup M∈M(p,n,s,µ∞) E L GHT -L(M) 2 2 ≤ 170µ ∞ σ 2 s 2 √ p log 3/2 (2np) + 6µ ∞ σ 2 sp.
Proof. Let us introduce the auxiliary notation ξ i = σ -1 (X i -µ i ) and write the decomposition

L GHT -L(M) = i∈S (X i -µ i ) - i∈S\ S(λ) (X i -µ 0 ) + i∈ S(λ)\S (X i -µ 0 ) = i∈S∩ S(λ) σξ i - i∈S\ S(λ) (µ i -µ 0 ) + i∈ S(λ)\S σξ i .
We have

L GHT -L(M) 2 ≤ σ i∈S∩ S(λ) ξ i 2 T 1 + i∈S\ S(λ) (µ i -µ 0 ) 2 T 2 +σ i∈ S(λ)\S ξ i 2 T 3 . ( 5
)
Let us use the notation Ξ S = [ξ i : i ∈ S] for the matrix obtained by concatenating the vectors ξ i with i ∈ S, and denote by ξ j S its j-th row. We can write

T 2 1 = Ξ S 1 S(λ) 2 2 ≤ s Ξ S 2 = s j∈[p] (ξ j S ) ξ j S ,
where 1 S(λ) is the vector in {0, 1} S having zeros at positions i ∈ S(λ) and zeros elsewhere.

One can observe that the s × s matrices (ξ j S ) ξ j S are independent and each of them is in expectation equal to the diagonal matrix diag(µ j S ). Since the operator norm of a diagonal matrix is equal to its largest (in absolute value) entry, we have

E T 2 1 ≤ µ ∞ sp + sE j∈[p] (ξ j S ) ξ j S -diag(µ j S ) ≤ µ ∞ sp + s E j∈[p] (ξ j S ) ξ j S -diag(µ j S ) 2 1/2
.

Upper bounding the spectral norm by the Frobenius norm, we arrive at

E j∈[p] (ξ j S ) ξ j S -diag(µ j S ) 2 ≤ E j∈[p] (ξ j S ) ξ j S -diag(µ j S ) 2 F = j∈[p] E (ξ j S ) ξ j S -diag(µ j S ) 2 F .
Using the fact that the fourth order central moment of the Poisson distribution P(a) is a + 3a 2 and that

ξ ij ∼ σ(P(σ -2 µ 0j ) -σ -2 µ 0j ), one can check that E[ξ 4 ij ] = σ 2 µ ij + 3µ 2 ij ≤ (2µ ∞ ) 2 , (6)
where the last inequality holds true since σ 2 ≤ (0.9) 1/3 µ ∞ .

In view of ( 6) and the Cauchy-Schwarz inequality, it holds that E[ (ξ j S )

ξ j S -diag(µ j S ) 2 F ] ≤ 2.5sµ 2 ∞ + (s 2 -s)µ 2 ∞ ≤ (s + 1) 2 µ 2 ∞ .
Putting these estimates together, we get

E T 2 1 ≤ µ ∞ sp + µ ∞ s(s + 1) √ p ≤ µ ∞ sp + 2µ ∞ s 2 √ p. (7)
We would like to stress here that this upper bound is quite rough, because of the step in which we upper bounded the spectral norm by the Frobenius one. Using more elaborated arguments on the expected spectral norm of a random matrix, such as the Rudelson inequality [START_REF] Rudelson | Random vectors in the isotropic position[END_REF] and (Boucheron et al., 2013, Cor. 13.9) combined with the symmetrisation argument (Koltchinskii, 2011, Theorem 2.1), it is possible to replace the rate s 2 √ p by s( √ s + √ p). However, we do not give more details on this computations since the rate s 2 √ p appears anyway as an upper bound on E[T 2 ].

To upper bound the term T 2 in (5), we define

τ i = σ 2 ( X i 1 + λ) and u = 40 log 3/2 (2np), so that τ i = σ 2 ( X i 1 + µ 0 ∞ √ p u). Note that i ∈ S \ S(λ) implies µ i -µ 0 2 2 < τ i -σ 2 ξ i 2 -2σ(µ i -µ 0 ) ξ i < τ i -σ 2 ξ i 2 + 2σ|(µ i -µ 0 ) ξ i | < τ i -σ 2 ξ i 2 + 0.5 µ i -µ 0 2 2 + 2σ 2 |(µ i -µ 0 ) ξ i | 2 µ i -µ 0 2 2 .
In other terms, setting θ i = (µ i -µ 0 )/ µ i -µ 0 2 , we have

S \ S(λ) ⊂ µ i -µ 0 2 2 < 2τ i -2σ 2 ξ i 2 + 4σ 2 |θ i ξ i | 2 .
As a consequence,

T 2 2 ≤ s i∈S\ S(λ) µ i -µ 0 2 2 ≤ 2µ ∞ s 2 σ 2 √ p u + 2σ 2 s i∈S ( X i 1 -ξ i 2 2 ) + + 4σ 2 s i∈S (θ i ξ i ) 2 .
Taking the expectation of both sides, in view of the independence of the components of ξ i and those of X i and the condition σ 2 ≤ µ ∞ , we arrive at

E T 2 2 ≤ 2µ ∞ σ 2 s 2 √ p u + 2µ ∞ σ 2 s 2 6p + 4µ ∞ σ 2 s 2 ≤ 89µ ∞ σ 2 s 2 √ p log 3/2 (2np). ( 8
)
The third term in (5) corresponds to the second kind error in the problem of support estimation and its expectation can be upper bounded using the Cauchy-Schwarz inequality as follows

E[T 2 3 ] ≤ n i ∈S E[ ξ i 2 2 1 ξ i 2 2 ≥τ i /σ 2 ] = n i ∈S j∈[p] E[ξ 2 ij 1 ξ i 2 2 ≥τ i /σ 2 ] ≤ n i ∈S j∈[p] E[ξ 4 ij ]P( ξ i 2 2 ≥ τ i /σ 2 ) 1/2 .
In view of (6), this yields

E[T 2 3 ] ≤ 2µ ∞ np i ∈S P( ξ i 2 2 ≥ τ i /σ 2 ) 1/2 . ( 9
)
For the last probability, we have

P( ξ i 2 2 ≥ τ i /σ 2 ) = P( ξ i 2 2 -X i 1 ≥ µ 0 ∞ √ p u).
To upper bound this probability, we apply Lemma 1 below 2 to η = σ -2 X j and ν = σ -2 µ j .

Lemma 1. Let η 1 , . . . , η p be independent random variables drawn from the Poisson distributions P(ν 1 ), . . . , P(ν p ), respectively. We set η = (η 1 , . . . , η p ) and ν = (ν 1 , . . . , ν p ). Then, for every u satisfying ν

-3/2 ∞ ≤ u ≤ 0.9ν 3/2 ∞ , we have P η -ν 2 2 -η 1 ≥ ν ∞ √ p u ≤ (2p + 1)e -cu 2/3 .
where c ≥ 0.38 is a universal constant.

2 The proof of Lemma 1 is postponed to Section 5

This yields that for

σ 3 µ -3/2 ∞ ≤ u ≤ 0.9σ -3 µ 3/2
∞ , we have

P( ξ i 2 2 ≥ τ i /σ 2 ) ≤ (2p + 1)e -cu 2/3 ,
where c ≥ 0.38 is an absolute constant. In view of (9), this implies the following upper bound on the expectation of T 3 :

E[T 2 3 ] ≤ 2µ ∞ n 2 p(2p + 1) 1/2 e -0.5cu 2/3 ≤ 4µ ∞ n 2 p 2 e -0.5cu 2/3 .
The choice u = 40 log 3/2 (2np) ensures that

E[T 2 3 ] ≤ 4µ ∞ n 2 p 2 e -2 log(2np) ≤ µ ∞ , (10)
which is negligible with respect to the upper bound on E[T 2 1 ] obtained above.

The result follows from upper bounds ( 5), ( 7), ( 8) and ( 10), as well as the choice of u. Indeed, by the Minkowskii inequality, we get

E[ L GHT -L(M) 2 2 ] 1/2 ≤ 11.85σµ 1/2 ∞ sp 1/4 log 3/4 (2np) + σ(µ ∞ sp) 1/2 .
Taking the square of both sides and using the inequality (a + b) 2 ≤ ( 6 /5)a 2 + 6b 2 , we get the claim of the theorem.

LOWER BOUNDS ON THE MINIMAX RISK OVER THE SPARSITY SET

Let us denote by M(n, p, s, µ 0 , µ ∞ ) the set of all p × n matrices with positive entries bounded by µ ∞ and having at most s columns different from µ 0 . This implicitly means, of course, that µ ∞ ≥ µ ∞ , otherwise the set is empty. The parameter set M(n, p, s, µ 0 , µ ∞ ) is the sparsity class over which the minimax risk will be considered. The result of the previous section implies that the minimax risk over M(n, p, s, µ 0 , µ ∞ ) is at most of order µ ∞ σ 2 s(p + s √ p), up to a logarithmic factor. In the present section, we will establish two lower bounds on the minimax risk over M(n, p, s, µ 0 , µ ∞ ).

The first lower bound shows that the minimax risk is at least of order s 2 , up to log factor. This implies that the group hard thresholding estimator defined in previous section is rate-optimal when p is fixed and s goes to infinity. The proof of this lower bound relies on the comparison of the Kullback-Leibler divergence between the probability distribution of X corresponding to M 0 = [µ 0 , . . . , µ 0 ], and the one corresponding to a matrix M, the n columns of which are randomly drawn from the set {µ 0 , µ 1 } with proportion (n -s) : s. When µ 1 is well chosen, the two aforementioned distributions are close, whereas the corresponding linear functional L(M 0 ) is sufficiently away from L(M).

Theorem 2. Assume that s ≥ 128, then

inf L sup M∈M(n,p,s,µ 0 ) E M L -L(M) 2 ≥ µ 0 ∞ 513 σ 2 s 2 log 1 + n 2s 2 ,
where the infimum is taken over all estimators of the linear functional.

Proof. We will make use of some lemmas, the statement of which requires additional notation.

As usual in lower bounds, we will make use of the Bayes risk: for every measure π on the set of p × n real matrices, we define the probability measure

P π ( • ) = R p×n P M ( • ) π(dM).
Let F : R p×n + → R m be any functional, not necessarily linear, defined on the set of matrices.

Lemma 2. Assume that π 0 and π 1 are two probability measures on the set of p × n matrices, R p×n + , and set M = max i Var π i (F ). Furthermore, suppose that KL(P π 1 , P π 0 ) ≤ 1 /8 and, for some v > 0,

E π 0 (F ) -E π 1 (F ) 2 > 2 v + 2 √ M .
Then, for any estimator F of F (M), we have

sup M∈M E M F -F (M) 2 2 ≥ 1 /8 -max i π i (M ) v 2 .
Lemma 3. Let µ 0 > 0 and > 0 be two real numbers. Define π 0 = δ ⊗n µ 0 as the Dirac mass concentrated in (µ 0 , . . . , µ 0 ) ∈ R n + and π 1 = λ ⊗n as the product measure of the mixture

λ = 1 - s 2n δ µ 0 + s 2n δ µ 0 +σ 2 ε .
Then

KL(P π 1 P π 0 ) ≤ s 2 4n e σ 2 ε 2 /µ 0 -1 .
We are now in a position to establish the lower bound of Theorem 2. For simplicity and without loss of generality, let us assume that the first entry of µ 0 , µ 01 , is the largest one. We choose two probability measures π p 0 and π p 1 on the set of p × n (intensity) matrices with positive entries, R p×n + , in the following way. The measure π p 0 is simply the Dirac mass at M 0 = (µ 0 , . . . , µ 0 , µ ∞ ). As for π p 1 , its marginal distribution of the first row equals the measure π 1 defined in Lemma 3, whereas all the other rows are deterministically equal to the corresponding rows of the matrix M 0 . So, for two matrices M and M randomly drawn from π p 0 and π p 1 , all the rows except the first are equal (and equal to the corresponding row of M 0 ). This implies that choosing the real ε by the formula

ε 2 = (µ 01 /σ 2 ) log(1 + n /2s 2 ),
we get (in view of Lemma 3) KL(P π p 1 P π p 0 ) ≤ 1/8. Furthermore, we have

E π 0 [L(M)] = nµ 0 , E π 1 [L(M)] = nµ 0 + 1 2 sσ 2 ε 0 p-1 ,
and

Var π 0 [L(M)] = 0, Var π 1 [L(M)] = 1 2 1 - s 2n sσ 4 ε 2 .
This means that, after upper bounding s by s 2 /128, the constant M of Lemma 2 can be chosen as M = 2 -8 s 2 σ 4 ε 2 . Hence, for v = ( 1 /8)σ 2 sε, we have

E π 0 [L(M)] -E π 1 [L(M)] 2 = sσ 2 ε 2 = 2v + 4 √ M .
This implies that the two conditions of Lemma 2 are fulfilled with

v = σ 2 sε 8 = σs √ µ 01 8 log 1/2 (1 + n /2s 2 ).
Therefore, taking into account the fact that π p 0 (M(n, p, s, µ 0 , µ ∞ )) = 1, we have sup M∈M(n,p,s,µ 0 ,µ∞)

E M L -L(M) 2 2 ≥ 1 /8 -π p 1 (M(n, p, s, µ 0 , µ ∞ ) ) v 2 .
To evaluate the probability π p 1 (M(n, p, s, µ 0 , µ ∞ ) ) , we remark that it is equal to P(ζ > s), where ζ is a random variable drawn from the binomial distribution B(n, ( s /2n)). Using the Bernstein inequality (Shorack and Wellner, 1986, p. 440), one can check that P(ζ > s) ≤ e -3s/5 ≤ e -76 , provided that s ≥ 128. This completes the proof.

It follows from Theorem 2 that the factor s 2 present in the upper bound of the minimax risk is unavoidable. The next result establishes another lower bound showing the the term sp present in the upper bound is unavoidable as well.

Theorem 3. For all 1 ≤ s ≤ n, if p ≥ 16, then inf L sup M∈M(n,p,s,µ 0 ,µ∞)

E M L -L(M) 2 ≥ 2 -14 µ ∞ σ 2 sp.
Proof. We set ε 2 = 2 -7 ( σ 2 /sµ∞) and, for any T ⊂ [p], define the matrix M T by

M T i,j =      µ ∞ (1 -ε), if (i, j) ∈ T × [s] µ ∞ , if (i, j) ∈ T × [s], µ 0,i , j ∈ [s].
We also set M 0 = M ∅ . It is clear that all these matrices M T belong to the sparsity class M(n, p, s, µ 0 , µ ∞ ). In addition, for any T , we have

KL P M T P M 0 = sµ ∞ σ 2 j∈T (1 -ε) log(1 -ε) + ε ≤ spε 2 µ ∞ σ 2 .
According to the Varshamov-Gilbert bound (Tsybakov, 2009, Lemma 2.9), there exist m = [e p/8 ] + 1 sets T 0 , . . . , T m ⊂ [p] such that T 0 = ∅ and With this choice of sets, we have for all i = j,

3 i = j =⇒ |T i T j | ≥ p 8 .
L M T i -L M T j 2 ≥ s 2 |T i T j |ε 2 µ 2 ∞ ≥ s 2 pε 2 µ 2 ∞ 2 3 .
Finally, using Theorem 2.5 in [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF] with α = 1/16, we get inf L sup M∈M(n,p,s,µ 0 ,µ∞)

E M L -L(M) 2 ≥ (1 -e -p/16 ) × 1 - 1 8 - 1 √ p s 2 pε 2 µ 2 ∞ 2 5 ≥ 2 -14 σ 2 spµ ∞ .
and the result follows.

CONCLUDING REMARKS AND OUTLOOK

Combining the results of the previous sections one can draw the optimality diagram, see Fig. 1.

In particular, Theorem 1 and Theorem 2 imply that the group hard thresholding estimator (3), (2) with the adaptively chosen threshold (4) is rate optimal when the term sp dominates s 2 √ p log 3/2 (2np). This corresponds to the "dimension dominates squared sparsity" regime p ≥ s 2 log 3 (2np). According to Theorem 1 and Theorem 3, the estimator defined by (3), ( 2) is rate optimal in the low dimensional setting p = O(1) as well. The range of values of p in which the determination of the minimax rate remains open corresponds to the regime p → ∞ so that p/s 2 = O(1).

In a companion paper [START_REF] Collier | Rate-optimal estimation of p-dimensional linear functionals in gaussian models[END_REF], the problem of estimating the linear functional L(M) is considered, when the observations are Gaussian vectors with unknown means µ i . It turns out that in the Gaussian setting the optimal scaling of the minimax risk is ps + s 2 , up to a logarithmic factor. However, this rate is achieved by an estimator that has a computational complexity that scales exponentially with n. Furthermore, in the Gaussian model, the upper bound relies heavily on th fact that the Gaussian distribution is sub-Gaussian.

The fact that the Poisson distribution is not sub-Gaussian (even though it is sub-exponential), does not allow us to apply the methodology developed in the Gaussian model.

It is worth mentioning here that although the risk bounds (both lower and upper) depend on the sparsity level s and the largest intensity µ ∞ , the group hard thresholding estimator is completely adaptive with respect to these parameters.

The problem we have considered in this work assumes that the noise variance σ 2 and the background signal µ 0 are known. It is not very difficult to check that the result of Theorem 1 can be extended to the case where these parameters are unknown, but an additional data set X1 , . . . , Xm is available, such that Xi 's are iid vectors drawn from the distribution σ 2 P p (µ 0 /σ 2 ). One can then estimate the background signal µ 0 by the empirical mean of Xi 's and σ 2 by the method of moments based on the second order moment. Next, these estimates can be plugged in the group hard thresholding estimator. If m is large enough, the error estimating µ 0 and σ 2 is negligible with respect to the error of estimating L(M). Determining more precise risk bounds clarifying the impact of m on the risk is an interesting avenue for future research.

PROOFS OF THE LEMMAS

Proof of Lemma 1. Let us introduce an auxiliary real number γ > 1.5 2 and denote

Z j = (η j -ν j -0.5) 2 , Z j,γ = Z j ∧ γ, Zj,γ = (Z j -γ) +
as well as

ν j,γ = E[Z j,γ ], νj,γ = E[ Zj,γ ].
It is clear that ν j,γ + νj,γ = ν j + 0.25 and we have

η -ν 2 2 -η 1 = j∈[p] (Z j,γ -ν j,γ ) + j∈[p] ( Zj,γ -νj,γ ) = Z γ 1 -ν γ 1 + Zγ 1 -νγ 1 . (11)
Since the random variables Z j,γ are independent and bounded, with P(Z j,γ ∈ [0, γ]) = 1, the Hoeffding inequality yields

P Z γ 1 -ν γ 1 > t ≤ e -2t 2 /(pγ 2 ) . (12)
On the other hand, it is clear that

P j∈[p] ( Zj,γ -νj,γ ) > 0 ≤ P max j∈[p] Z j > γ ≤ p max j∈[p] P |η j -ν j -0.5| > √ γ ≤ p max j∈[p] P |η j -ν j | > √ γ -0.5 ≤ p max j∈[p] P |η j -ν j | > (2/3) √ γ .
Using the Bennet inequality for the Poisson random variables-see, for instance, [START_REF] Pollard | A few good inequalities[END_REF], Problem 15)-we get 4 P |η j -ν j | > (2/3) √ γ ≤ P η j > ν j + (2/3) √ γ + P η j < ν j -(2/3) √ γ ≤ e -1.54γ/(9ν j ) + e -2γ/(9ν j ) ≤ 2e -γ/(6ν j ) , ∀γ ≤ (9/4)ν 2 ∞ . (13) Combining ( 11), ( 12) and ( 13), we get

P η -ν 2 2 -η 1 ≥ t ≤ inf γ∈[(3/2) 2 ,(9/4)ν 2 ∞ ]
e -2t 2 /(pγ 2 ) + 2pe -γ/(6ν∞) .

Let us set t = ν ∞ √ p u and γ 3 = 12 u 2 ν 3 ∞ . For u such that ν

-3/2 ∞ ≤ u ≤ 0.9ν 3/2
∞ , the aforementioned value of γ belongs to the range over which the infimum of the last display is computed. This yields P ην 2 2 -η 1 ≥ ν ∞ √ p u ≤ (2p + 1)e -( 12) 1/3 u 2/3 /6 , and the claim of the lemma follows.

Proof of Lemma 2. Using Markov's inequality, we get

sup M∈M v -2 E M F -F (M) 2 2 ≥ sup M∈M P M F -F (M) 2 ≥ v ≥ max i=0,1 P π i F -F (M) 2 ≥ v -max i=0,1 π i (M ).
Now, using the test statistic ψ = arg min i∈{0,1} F -E π i (F ) 2 and the notation M = max i Var π i (F ) = max i E π i F -E π i [F ] 2 2 , we have

P π i (ψ = i) = P π i F -E π 1-i (F ) 2 ≤ F -E π i (F ) 2 ≤ P π i F -E π i (F ) 2 ≥ v + M/α .
Finally, using the triangle inequality again,

P π i (ψ = i) ≤ P π i F -F (M) 2 ≥ v + P π i F (M) -E π i (F ) 2 > √ α -1 M ≤ α + P π i F -F (M) 2 ≥ v , so that sup M∈M P M F -F (M) 2 ≥ v ≥ max i=0,1 P π i (ψ = i) -α -max i=0,1 π i (M ).
The result is now a direct application of Theorem 2.2 in [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF] using the assumption on the Kullback-Leibler divergence.

Proof of Lemma 3. Let us denote by Q µ , for every vector µ ∈ (0, ∞), the Poisson distribution P(σ -2 µ), and by Q λ = (0,∞) Q µ λ(dµ), the λ-mixture of Poisson distributions. We have, for every X ∈ N,

d Q λ d Q µ 0 (X) = 1 + s 2n d Q µ 0 +σ 2 d Q µ 0 (X) -1 = 1 + s 2n 1 + σ 2 ε µ 0 X e -ε -1 .
Using the inequality log(1 + t) ≤ t, we get

KL(Q λ Q µ 0 ) ≤ s 2n N 1 + σ 2 ε µ 0 X e -ε Q λ (dX) - s 2n .
On the other hand, for every µ ∈ (0, ∞) p , we have

1 + σ 2 ε µ 0 X Q µ (dX) = k∈N 1 + σ 2 ε µ 0 k µ k σ 2k k! e -µ/σ 2 = exp µε /µ 0 , so that KL(Q λ Q µ 0 ) ≤ s 2n 1 - s 2n + s 2
4n 2 e σ 2 2 /µ 0 -

s 2n = s 2
4n 2 e σ 2 ε 2 /µ 0 -1 .

The result follows since, by independence, KL(P π 1 P π 0 ) = nKL(Q λ Q µ 0 ).

Fig 1 .

 1 Fig 1. The rate optimality diagram in the minimax sense for the group hard thresholding estimator.

By heteroscedasticity we understand here the fact that if the mean is unknown, then the variance is unknown as well and these two parameters are interrelated

In what follows, stands for the symmetric difference of two sets.

We use here the fact that the Bennet ψ * satisfies ψ * (x) ≥ 0.772 when x ∈ [0, 1].

ACKNOWLEDGEMENTS O. Collier's research has been conducted as part of the project Labex MME-DII (ANR11-LBX-0023-01). The work of A. Dalalyan was partially supported by the grant Investissements d'Avenir (ANR-11IDEX-0003/Labex Ecodec/ANR-11-LABX-0047).


Christoph Weiss

Niklas Teichmann

Generation of mesoscopic superpositions of a binary Bose-Einstein condensate in a slightly asymmetric double well

Keywords: numbers: 03.75.Gg, 03.75.Lm, 03.75.Mn Keywords: Mesoscopic entanglement, Bose-Einstein condensation, double-well potential

A previous publication [Europhysics Letters 78, 10009 (2007)] suggested to coherently generate mesoscopic superpositions of a two-component Bose-Einstein condensate in a double well under perfectly symmetric conditions. However, already tiny asymmetries can destroy the entanglement properties of the ground state. Nevertheless, even under more realistic conditions, the scheme is demonstrated numerically to generate mesoscopic superpositions.

I. INTRODUCTION

Generation of mesoscopic entanglement in a Bose-Einstein condensate (BEC) [1,2,3,4,5,6,7,[START_REF] Creffield | Quantum control and entanglement using periodic driving fields[END_REF] is still a challenge of fundamental research. A possible system for which entanglement could be achieved is a BEC in a double-well potential [6,[START_REF] Gati | [END_REF]. While attractive condensates can become unstable [10], for repulsive single species condensates the ground state is not the desired superposition (see, e.g., Ref [3]). Nevertheless, by phase engineering fidelities (i.e. the probability to be in an ideal superposition) of about 60% have been obtained numerically [6].

In Ref. [11] we suggested to use a binary condensate to produce Bell-like mesoscopic superpositions with fidelities above 95% by using a method of tunneling control which was very recently realized experimentally in an optical lattice [START_REF] Lignier | Dynamical control of matter-wave tunneling in periodic potentials[END_REF]. However, to achieve the entanglement generation, the ground state of a two-component condensate of exactly the same number of particles of each kind was followed adiabatically in a perfectly symmetric double well potential [START_REF]As the coherent control was performed via time-periodic potential differences[END_REF]. Given the fact that even very small asymmetries can destroy entanglement properties of eigen-functions [START_REF] Dounas-Frazer | Tunneling resonances and entanglement dynamics of cold bosons in the double well[END_REF], this rises the question if the generation of mesoscopic entanglement proposed in Ref. [11] would also be possible under more realistic conditions, i.e., slightly different particle numbers and a small tilt.

The focus of this manuscript is to investigate the scheme of entanglement generation including asymmetries in both the potential and the particle numbers. In Sec. II the model used to describe a two-component BEC in a double well is introduced. While Sec. III explains how the mesoscopic superpositions can be characterized, Sec. IV outlines the noticeable differences between the ideal case considered in Ref. [11] and the more realistic * Electronic address: christoph.weiss@lkb.ens.fr (slightly asymmetric) situation considered here for which the entanglement generation no longer is entirely adiabatic.

II. MODEL, CONTROLLING THE ENTANGLEMENT GENERATION

To describe a binary Bose-Einstein condensate in a double well, we use the Hamiltonian in two-mode approximation (Ref. [START_REF] Milburn | [END_REF] cf. Ref. [16]):

Ĥ = - Ω 2 â1 â † 2 + â † 1 â2 + κ A â † 1 â † 1 â1 â1 + â † 2 â † 2 â2 â2 - Ω 2 b1 b † 2 + b † 1 b2 + κ B b † 1 b † 1 b1 b1 + b † 2 b † 2 b2 b2 +2 κ AB â † 1 â1 b † 1 b1 + â † 2 â2 b † 2 b2 + (ε + µf (t)) â † 2 â2 -â † 1 â1 + b † 2 b2 -b † 1 b1 , (1)
where â( †) j ( b( †) j ) creates, respectively annihilates, a boson of type A (B) in well j, Ω is the tunneling frequency, κ A (κ B ) is the interaction parameter for atoms of type A (B) and κ AB models the interaction between distinguishable particles. Such Hamiltonians have also been used to describe entanglement generation beyond BECphysics [17]. The (unwanted) small tilt which might not be avoidable experimentally is represented by ε.

The time-dependent potential difference ∝ f (t) can be used to perform tunneling control via resonant tunneling [7,18] (cf. Ref. [4]). References [START_REF] Creffield | Quantum control and entanglement using periodic driving fields[END_REF]11] generate multi-particle entanglement by periodic driving,

f (t) = sin(ωt) , (2)
to effectively renormalize the tunneling frequency Ω by the J 0 -Bessel function in the high frequency limit to

Ω eff = ΩJ 0 (2µ/ω) . (3)
This mechanism for tunneling control (cf. [19]), which on the single particle level has been investigated theoretically already for some 16 years [20,21,22,23], has very recently been observed experimentally for a Bose-Einstein condensate in an optical lattice [START_REF] Lignier | Dynamical control of matter-wave tunneling in periodic potentials[END_REF]. Another possibility to achieve the entanglement generation of Ref. [11] is to decrease the tunneling frequency by increasing the distance d between the two wells. Approximating the wave-functions of particles of mass m in the two wells by the ground states of harmonic oscillators of oscillation frequency ω and oscillator length ℓ osc = /(2m ω), one obtains [3,[START_REF] Milburn | [END_REF]]

Ω ≃ ω d ℓ osc 2 π exp - d 2 ℓ 2 osc (4)
for the tunneling frequency and

κ = a 4 √ πmℓ 3 osc , ( 5
)
where a is the s-wave scattering length, for the interaction parameters of Eq. ( 1) as a function of Ω. A possible realization of such a model could be two hyperfine-states of, e.g, 87 Rb. Feshbach resonances for interaction between different hyper-fine states have been observed experimentally [24], thus providing a way to slightly increase the interaction between distinguishable particles as compared to identical particles.

To coherently generate entanglement, Ref. [11] suggests to decrease the (effective) tunneling frequency for such a situation, thus adiabatically following the ground state towards highly entangled mesoscopic states. As decoherence-times tend to rapidly decrease with increasing particle numbers, not too large particle numbers should be used to investigate entanglement generation. Experimentally, Bose-Einstein condensates of the order of 100 particles have been realized [25]; in Ref. [11] Inset: The same situation for an untilted double well. Even without the tilt, the ratio κAB/κA has to be larger than in Ref. [11] in order to ensure that the Bell-state |ψ+ still is the ground state for a symmetric potential.

III. IDENTIFYING MESOSCOPIC SUPERPOSITIONS

For the symmetric conditions of Ref. [11] the fidelity obtained by projecting the wave-function |ψ on the Bell state

|ψ ± = 1 √ 2 |0, N A A |N B , 0 B ± |N A , 0 A |0, N B B (6)
provides an excellent measure of entanglement generation,

p ± = | Ψ|ψ ± | 2 , (7)
if p ± > 0.5. In Fig. 2 the tunneling frequency in an asymmetric double well containing a binary BEC with N A = N B is reduced linearly during half the time and during the other half Ω held at zero:

Ω = 1 -τ τramp : 0 ≤ τ ≤ τ ramp 0 : τ > τ ramp , τ = tΩ . ( 8
)
For a slightly tilted double well, the ground state cannot be one of the two superpositions of Eq. ( 6): In the limit Ω → 0 the two wave-functions used in the superpositions (6) do have different energies. Thus, even if a perfect superposition was reached, the phase would oscillate similar to Fig. 2. The fact that the fidelities stay below 1 indicates that more Fock states then just those included in Eq. ( 6). Increasing the tilt by, say, a factor of ten would further increase the oscillation frequency as the energy difference between both wells increases while increasing it by a much larger factor would destroy the generation of mesoscopic superpositions. Periodic driving introduces an additional frequency, thus leading to even larger oscillations frequencies in the fidelities.

For the purpose of identifying mesoscopic superpositions independent of these oscillations, we define the modified fidelity:

p fid = 1 2 ( √ p 1 + √ p 2 ) 2 , (9)
where

p 1 = ( Ψ|0, N A A |N B , 0 B ) 2 (10)
and

p 2 = ( Ψ|N A , 0 A |0, N B B ) 2 . ( 11
)
As for the usual fidelity, this measure is 1 for a perfect superposition and already p fid ≈ 0.5 might signify that there is no superposition at all. The modified fidelity ( 9) can be derived by projecting the wavefunction on the superposition

|ψ ϕ = 1 √ 2 |0, N A |N B , 0 + e iϕ |N A , 0 |0, N B (12)
and maximizing the resulting projection

p ϕ = | Ψ|ψ ϕ | 2 (13)
as a function of ϕ.

IV. ENTANGLEMENT GENERATION ONLY PARTIALLY ADIABATIC

In order to understand why generation of mesoscopic superpositions is still possible, Fig. 3 shows the modified fidelity [START_REF] Gati | [END_REF] for different ramping times τ ramp (see Eq. ( 8)). If the ramping is done on too short timescales, the oscillations characteristic for too fast, i.e. nonadiabatic, ramping occur and high fidelities cannot be reached. If, on the other hand, the ramping is done too slowly, the system does not reach a mesoscopic superpositions either as it adiabatically follows the ground state for too long and thus becomes rather asymmetric. For moderate ramping, the system first adiabatically follows the ground state thus leading to a wave-function with a bimodal probability distribution. Then, the symmetry of the wave-function is conserved (instead of tunneling through energetically higher states with half the particles of at least one kind in both of the wells) to finally reach a mesoscopic superposition.

If the tunneling control is done via slowly variing the amplitude µ of periodic driving as in Ref. [11] (cf. Eqs. ( 2) and ( 3)), i.e., Ω = const., ω ≫ Ω and A ramping time of τramp = 50 (dotted line) is much too short as can be seen in the oscillations indicating non-adiabaticity. On the other hand τramp = 50000 (dashed line) is too slow; it adiabatically follows the ground state up to values where it becomes rather asymmetric. This leads to a final state with a modified fidelity well below 90%. Inset: For a ramping time of τramp = 1000, the final modified fidelity is shown as a function of the particle number difference ∆ with NA = 50 -∆/2 and NB = 50 + ∆/2. Up to ∆/50 = 20%, the value used for most of the figures, a transfer to highly entangled states is achieved. For higher ∆, the ramping-time and/or the ratio of the interaction parameters would have to be adjusted.

µ ≃ 2.4048ω 2 τ τ ramp , 0 ≤ τ ≤ τ ramp , (14)
where 2.4048 . . . is the first zero of the J 0 -Bessel function (see Eq. ( 3)), the generation of mesoscopic entanglement also occurs. For the parameters in Fig. 2 with a ramping time of τ ramp = 1000, Fig. 4 shows the emergence of an entangled wave-function from an initial state for which the probability distribution is not bimodal.

V. CONCLUSIONS

While already two independent Bose-Einstein condensates provide interesting effects (see, e.g. Ref. [26]), twocomponent BECs provide an even wider variety of effects [27,28]. Here it has been shown that for a binary condensate in a double-well potential, the generation of mesoscopic entanglement suggested in Ref. [11] works also under less ideal conditions, i.e., if the particle numbers are different and if the double well is not perfectly symmetric.

For a given experimental situation, the linear ramping used here in the spirit of coherent control theory can be optimized both toward short ramping times and large fidelities, e.g., by methods of optimal control theory [29].

FIG. 1 :

 1 FIG. 1: (color online) Sketch of one of the target states for entanglement generation in a double well: a Bell-like superposition of two wave-functions for a binary Bose-Einstein condensate.

FIG. 2 :

 2 FIG.2: Fidelity (7) as a function of dimensionless time τ = Ωt if the tunneling frequency is reduced linearly from Ω(τ = 0) = 1 to Ω(τ = 500) = 0 and then held at this value. Solid line: p+, dashed line p-. The double well is slightly asymmetric (ε = 0.001) with NA = 45, NB = 55, κA = κB = 0.04; κAB = 0.05. Inset: The same situation for an untilted double well. Even without the tilt, the ratio κAB/κA has to be larger than in Ref.[11] in order to ensure that the Bell-state |ψ+ still is the ground state for a symmetric potential.

FIG. 3 :

 3 FIG.3: Modified fidelity (9) as a function of time t. For ramping time τramp = 1000 (solid line), p fid reaches 95.61% by first adiabatically following the ground state and then nonadiabatically conserving the symmetry of the wave-function. The other parameters are chosen as in Fig 2.A ramping time of τramp = 50 (dotted line) is much too short as can be seen in the oscillations indicating non-adiabaticity. On the other hand τramp = 50000 (dashed line) is too slow; it adiabatically follows the ground state up to values where it becomes rather asymmetric. This leads to a final state with a modified fidelity well below 90%. Inset: For a ramping time of τramp = 1000, the final modified fidelity is shown as a function of the particle number difference ∆ with NA = 50 -∆/2 and NB = 50 + ∆/2. Up to ∆/50 = 20%, the value used for most of the figures, a transfer to highly entangled states is achieved. For higher ∆, the ramping-time and/or the ratio of the interaction parameters would have to be adjusted.

FIG. 4 :

 4 FIG. 4: (color online)The projection of the wave-function on Fock-states |i, NA -i |j, NB -j for a binary condensate for which the tunneling control is performed via periodic driving as suggested in Ref.[11] (see Eqs. (14) and (3)). While the initial state at τ = 0 (top) is no bimodal distribution at all, a mesoscopic superposition is generated already at τ = 400 (bottom). The parameters where chosen as in Fig3with a ramping time of τramp = 1000. The driving frequency ω = 40Ω; at τ = τramp the final fidelity is again well above 90%.

  FIG. 4: (color online)The projection of the wave-function on Fock-states |i, NA -i |j, NB -j for a binary condensate for which the tunneling control is performed via periodic driving as suggested in Ref.[11] (see Eqs. (14) and (3)). While the initial state at τ = 0 (top) is no bimodal distribution at all, a mesoscopic superposition is generated already at τ = 400 (bottom). The parameters where chosen as in Fig3with a ramping time of τramp = 1000. The driving frequency ω = 40Ω; at τ = τramp the final fidelity is again well above 90%.

Acknowledgments We would like to thank Y. Castin for his support. N. T. thanks the Institut Henri Poincare-Centre Emile Borel for hospitality and support. Funding by the European Union (C. W. through contract MEIF-CT-2006-038407 and N. T. through contract MEST-CT-2005-019755) is gratefully acknowledged.