Medina Hadjem

Osman Salem

Farid Naït-Abdesselam

An ECG Monitoring System For Prediction Of Cardiac Anomalies Using WBAN

Keywords: WBAN, ECG, CVD, Bayesian Network Classifier

Cardiovascular diseases (CVD) are known to be the most widespread causes to death. Therefore, detecting earlier signs of cardiac anomalies is of prominent importance to ease the treatment of any cardiac complication or take appropriate actions. Electrocardiogram (ECG) is used by doctors as an important diagnosis tool and in most cases, it's recorded and analyzed at hospital after the appearance of first symptoms or recorded by patients using a device named holter ECG and analyzed afterward by doctors. In fact, there is a lack of systems able to capture ECG and analyze it remotely before the onset of severe symptoms. With the development of wearable sensor devices having wireless transmission capabilities, there is a need to develop real time systems able to accurately analyze ECG and detect cardiac abnormalities. In this paper, we propose a new CVD detection system using Wireless Body Area Networks (WBAN) technology. This system processes the captured ECG using filtering and Undecimated Wavelet Transform (UWT) techniques to remove noises and extract nine main ECG diagnosis parameters, then the system uses a Bayesian Network Classifier model to classify ECG based on its parameters into four different classes: Normal, Premature Atrial Contraction (PAC), Premature Ventricular Contraction (PVC) and Myocardial Infarction (MI). The experimental results on ECGs from real patients databases show that the average detection rate (TPR) is 96.1% for an average false alarm rate (FPR) of 1.3%.

I. INTRODUCTION

Cardiovascular diseases (CVD) provoke serious anomalies that involves the heart or the blood vessels. According to the World Health Organization [START_REF]World health organisation[END_REF], these diseases are the leading causes of death worldwide. In 2008, an estimated 17.3 million people died from CVD and almost 23.3 million people will die annually from CVD by 2030. Prevention and early detection are key factors in reducing the number of victims and the associated healthcare costs.

Recent technological advances in wireless networks, sensors conception and miniaturization can revolutionize the way how healthcare services are deployed and delivered. In fact, several manufacturers propose today sensors with wireless transmission capabilities able to collect various physiological parameters in real time such as Heart Rate, Temperature, Blood Pressure, Oxygen Saturation etc. These sensors are non-invasive and allow remote in-home, outdoor or ambulatory monitoring. Also, they can be deployed in a Wireless Body Area Network (WBAN) where several sensors are attached to the human body and transmit collected data to a gateway device (Smartphone, Tablet, PDA, etc.) which processes and transmits gathered information in real time to remote healthcare databases so professionals can access and analyze these data and make appropriate medical decisions. An example of a WBAN architecture for ECG monitoring is shown in figure 1.

ECG

Remote Monitoring Center

Fig. 1. A WBAN Architecture For ECG Monitoring WBAN systems is a new promising technology and various applications can be considered but it's still very challenging and several issues are open. One of the major one is the automatic and early detection of medical anomalies symptoms in a WBAN environment, especially in the case of CVD, this issue is of utmost importance for patients and health professionals and is highly supported by governments and research.

One of the most frequent CVD is arrhythmia which is a trouble of the cardiac rhythm that can take many forms, the most common are Premature Ventricular Contractions (PVC) and Premature Atrial Contractions (PAC). These two anomalies are due to the premature discharge of an electrical impulse in the atrium of the heart in the case of PAC [START_REF] Nainggolan | Risk factors identified for premature atrial contractions[END_REF] and the ventricles of the heart in the case of PVC [START_REF] James | Premature ventricular contraction[END_REF]. This results in an extra heartbeat that occurs earlier than the next regular beat should have occurred. Another serious heart anomaly is Myocardial Infarction (MI) commonly known as heart attack. MI is an acute ischemic heart disease that happens when the arteries that supply blood to heart muscle become hardened and narrowed and can't provide correctly the blood and oxygen to the heart muscle which can suffer from permanent damages. Many other CVD exist and are detailed in [START_REF] Morris | ABC of clinical electrocardiography[END_REF].

Symptoms of CVD can be observed and confirmed by various medical examinations made by doctors, the principal one is the Electrocardiogram (ECG). The ECG represents the contraction and relaxation of cardiac muscle resulting from the depolarization and re-polarization of myocardial cells. These electrical changes are recorded via electrodes placed on the limbs and chest wall and are transcribed on a graph paper.

The most widely used ECG recording technique is the 12lead ECG. This technique uses 10 recording sites from which are recorded 12 signals called leads among which 6 Precordial Leads (V1, V2, V3, V4, V5, V6), 3 Limb leads (I, II, III) and 3 augmented Limb leads (aVR, aVL, aVF). Each lead views the heart from a different angle and is composed of 5 waves that refer to a depolarization or re-polarization of some region in the heart, these waves, are symbolized P , Q, R, S, and T and are separated by intervals and segments. Figure 2 illustrates a typical one-cycle ECG waveform with the 5 waves and the main intervals and segments. ECG should be acquired and correctly interpreted as soon as possible. In practice, traditional techniques record ECG for a relatively short time (minutes in hospitals to 24h or 48h hours with a Holter ECG) and some episodic symptoms my not occur during the monitoring. Furthermore, ECG is generally recorded after first symptoms so heart damages may already have occurred. Therefore, It is needed to have a sufficiently prolonged ECG recording systems with real time monitoring and analysis capabilities.

The aim of this paper is to introduce an efficient WBAN ECG monitoring system with the ability to distinguish between normal ECGs from three cardiac anomalies discussed above: PAC, PVC and MI. The rest of this paper is organized as follows: Section II surveys related work. Section III presents the proposed approach and Section IV discusses experimental results. Finally, Section V concludes this paper.

II. RELATED WORK

During the last decade, a lot of out-hospital wireless ECG monitoring systems have been proposed by researchers and industries. From our study, we found that the majority of which capture, process and send ECG data to specialists in order to remotely follow their patients, but they generally don't detect cardiac anomalies within the collected data or achieves a basic ECG analysis like heart rate variability. This is the main limitation of these systems like CodeBlue [START_REF] Malan | Codeblue: An ad hoc sensor network infrastructure for emergency medical care[END_REF], LiveNet [START_REF] Sung | Livenet: Health and lifestyle networking through distributed mobile devices[END_REF], MyHeart [START_REF] Luprano | Combination of body sensor networks and on-body signal processing algorithms: the practical case of myheart project[END_REF], MEDISN [START_REF] Ko | Medisn: Medical emergency detection in sensor networks[END_REF] or industrial products like AliveCor [START_REF] Saxon | Ubiquitous wireless ecg recording: a powerful tool physicians should embrace[END_REF], PhysioMem [START_REF]Physiomem pm 1000[END_REF] or LifeMonitor [START_REF]The eq02 lifemonitor[END_REF].

More recently, some works start to address the problem of ECG anomalies detection in a wireless or mobile environment, such as in [START_REF] Sankari | Heartsaver: A mobile cardiac monitoring system for auto-detection of atrial fibrillation, myocardial infarction, and atrio-ventricular block[END_REF], a system named HeartSaver is proposed, this system analyzes ECG in real time and includes automated detection of three cardiac pathologies (Atrial fibrillation, Myocardial infarction, and Atrio-ventricular block) but this system was tested using Matlab simulated pathological ECG instead of real ECG data that is affected by various environment parameters and much more complex to analyze. In [START_REF] Ashrafuzzaman | Heart attack detection using smart phone[END_REF], a heart attack detection application using Smartphone is proposed, this system is based on placing index fingers on mobile camera, the application will calculate the heart rate based on detection of peaks of blood but this technique is limited and simplistic because detection of heart attacks is complex and needs analyzing specific features in ECG. In [START_REF] Zhou | A real-time continuous cardiac arrhythmias detection system: Recad[END_REF], the RECAD system is designed, it is a real-time wireless monitoring ECG system for Arrhythmia detection. The system in [START_REF] Oresko | A wearable smartphone-based platform for real-time cardiovascular disease detection via electrocardiogram processing[END_REF] is a wearable Smartphone based platform for real-time heart beat anomaly detection as normal, PVC, right bundle branch block beat (RBBB), paced beat (PACE), and PFUS beats. Some industry devices for arrhythmia detection are available on the market like CardioNet [START_REF]Mcot cardionet[END_REF], Smartheart [START_REF]Smartheart tm[END_REF] and BodyGuardian [START_REF]Bodyguardian remote monitoring system[END_REF].

In this study, we also review works that propose algorithms and models for ECG analysis and anomalies detection. In most cases, classification algorithms are used to classify ECG into normal or abnormal : In [START_REF] Jovic | Electrocardiogram analysis using a combination of statistical, geometric, and nonlinear heart rate variability features[END_REF], Jovic et al. proposed a classification method based on analysis of a combination of Heart Rate Variability and approximate entropy features in order to classify ECG in four classes (normal heart rhythm, arrhythmia, supra-ventricular arrhythmia, and congestive heart failure). ECG records from online databases were analyzed by seven clustering and classification algorithms. Results show that the top 3 accurate classification methods are Random Forest (RF) with 99.6%, Bayesian Networks with 99.4% and Support Vector Machines (SVM) with 98.4%. In [START_REF] Jadhav | Artificial neural network models based cardiac arrhythmia disease diagnosis from ecg signal data[END_REF], Jadhav et al. proposed an arrhythmia disease diagnosis method using Artificial Neural Network (ANN) classifier of 12-lead ECG recordings data. Their model classifies ECG as normal or abnormal, they used UCI Arrhythmia data set to train and test three different Multilayer perceptron ANN models. They obtained an accuracy and sensitivity of 86.67% and 93.75% respectively with the two first models and a specificity of 93.1% for the third model. In [START_REF] Giri | Automated diagnosis of coronary artery disease affected patients using lda, pca, ica and discrete wavelet transform[END_REF] a methodology is proposed for the automatic detection of normal and Coronary Artery Disease conditions bu testing different combinations of data reduction techniques (PCA, LDA, ICA) and four different classifiers (SVM, GMM, PNN, KNN). Their results showed that the combination of ICA with GMM resulted in highest accuracy of 96.8% compared to other combinations.

These models achieve a good accuracy but they are generally not suitable for a WBAN environment due to their high computational complexity. In this contribution, a methodology has been developed to combine the WBAN technology with statistical classification technique to design a new ECG monitoring system with cardiac anomaly detection capabilities that:

-Monitors real time ECG in a WBAN environment.

-Predicts 3 cardiac anomalies : PVC, PAC and MI.

-Uses numerically efficient algorithms.

III. PROPOSED APPROACH

We consider a real deployment scenario where wireless sensors send ECG measurements to a Smartphone in real time. ). -A U serLayer for patients and remote healthcare professionals who will receive alerts and patients data in real time and take appropriate decisions about the patient condition. In this paper, the proposed approach focuses only on the application layer of our system. This layer is described in Figure 3 and achieve 3 main steps detailed bellow. -The power line interference around 50Hz or 60Hz.

-The baseline wandering with a frequency bellow 0.5Hz.

-Electrodes motion artifacts ranging from 1Hz to 10Hz.

-Electromyographic (EMG) noise from 25Hz to 100Hz. Power line noise comes from electronic circuits of ECG sensors, it is a narrow-band signal centered at a frequency of 50Hz or 60Hz with a bandwidth of less than 1Hz, it has a form of a sinusoid with a significant amplitude compared to the magnitude of the digitized signal. Generally this noise is filtered by ECG signal acquisition hardware.

The baseline wandering usually comes from patient movements and respiration at frequencies between 0.15Hz and 0.5Hz and causes changes in iso-electric line position of the measured ECG. It can be suppressed without loss in the original signal. It is possible to remove it by a high-pass digital filter or by the use of a standard Wavelet Transform (WT). We choose to use a 0.5Hz FIR (Finite Impulse Response) high-pass digital filter.

Other noises are more difficult to remove because they may be a complex stochastic processes within a wideband and traditional digital filters can't remove them. Instead, Discrete Wavelet Transform (DWT) technique is widely used or more recently Undecimated Wavelet Transform (UWT) which has a better balance between smoothness and accuracy than the DWT and ensures no loss of signal sharpest features. More details on its use for ECG denoising can be found in [START_REF] Raj | Ecg signal denoising using undecimated wavelet transform[END_REF]. In our system, we use UWT for both pre-processing and features extraction stage where this technique will be further detailed.

Once the signal cleaned and denoised, the system extracts the main ECG features. In this step, we first detect all QRS complex peaks and then detect other waves P and T . As previously mentioned, the features extraction algorithm uses UWT, a variant of Wavelet Transform (WT) which is one of the most popular methods for time frequency transformations. It can provide the frequencies of a signal and the time associated to those frequencies, making it very convenient for its application in numerous fields. WT should allow only changes in time extension but not shape. This is performed by choosing suitable basis functions ψ. Generally, the WT can be expressed with the following equation:

[W ψ f ] (a, b) = 1 |a| ∞ -∞ ψ x -b a f (x)dx (1)
The wavelet coefficients c jk are then given by:

c jk = [W ψ f ] 2 -j , k2 -j (2)
a = 2 -j is called the binary or dyadic dilation and b = k2 -j is called the binary or dyadic position. The basis function ψ can be regarded as an impulse response of a system with which the function x(t) has been filtered. Two implementations of WT are commonly known : Discrete Wavelet Transform (DWT) and Continuous Wavelet Transform (CWT). Detailed information about WT and its applications can be found in [START_REF] Akansu | Emerging applications of wavelets: A review[END_REF].

The UWT is an improved design of DWT, it is implemented in our system based on the method proposed in [START_REF] Bhyri | Estimation of ecg features using labview[END_REF]. The method first decomposes the ECG signal by 8 level Daubechies wavelets and then reconstructs the signal from approximation coefficients and detail coefficients of all frequency bands. After, peaks detection step is performed to detect QRS complex, P and T peaks by properly setting a threshold that specifies the value which the algorithm uses to reject peaks of a particular amplitude, and a windows width that specifies the current number of samples of the signal for the detection of onsets and offsets. The steps performed by each QRS, P and T peak detector algorithm are described bellow : -Perform UWT with the Daubechies wavelet on input ECG.

-Detect the zero crossing points in the detail coefficients.

-Set the zero crossing points at the largest scale as the coarse estimation of the real peaks. -Searches the fine scale for the corresponding nearest zero crossing point for each detected point. -Repeats step 4 until it reaches the finest scale, which is the first level. Once all peaks detected, the onsets and offsets are calculated by selecting an optimum window width and using the peaks as references inside the windows in order to scan the waves on both the sides of the peaks to get the zero crossing points for obtaining onsets and offsets. If the signal is not crossing zero line, then the minimum values in that window are considered as onsets and offsets. At the end of these two steps, the 6 features bellow are extracted from each one-cycle ECG (heart-beat): -P onset: Onset of the P wave of ECG. -P of f set: Offset of the P wave of ECG.

-QRSonset: Onset of the QRS complex of ECG.

-QRSof f set: Offset of the QRS complex of ECG.

-T onset: Onset of the T wave of ECG. -T of f set: Offset of the T wave of ECG. From these features, amplitudes and durations of each ECG wave, segment and interval (see Figure 2) can be calculated easily. We choose to calculate the 9 most relevant parameters used by cardiologists [START_REF] Morris | ABC of clinical electrocardiography[END_REF] A Bayesian network is a statistical model that represents a set of random variables and their conditional dependencies via a Directed Acyclic Graph (DAG), this is why it can be defined as a probabilistic directed acyclic graphical model. In a medical context, a Bayesian Network could represent the probabilistic relationships between diseases and symptoms. Given symptoms, the network can be used to compute the probabilities of the presence of various diseases. In a Bayesian Network, the nodes of the graphs represent random variables (they may be observable quantities, latent variables, unknown parameters or hypotheses) and the edges represent conditional dependencies. Nodes which are not connected represent variables which are conditionally independent of each other. Each node is associated with a probability function that takes as input a particular set of values for the node's parent variables and gives the probability of the variable represented by the node. Bayesian Networks that model sequences of variables are called Dynamic Bayesian Networks.

Let G = (H, E) a Directed Acyclic Graph (or DAG) where H are the nodes and E the edges and let X = (X h ) h∈H a set of variables indexed by H.

X is a Bayesian Network if its joint probability density function can be written as a product of the individual density functions, conditional on their parent variables:

p(x) = h∈H p x h |x pa(h) (3)
where pa(h) is the set of parents of h.

For any set of random variables, the probability of any member of a joint distribution can be calculated from conditional probabilities using the chain rule as follows:

P(X 1 = x 1 , . . . , X n = x n ) = n h=1 P (X h = x h | X h+1 = x h+1 , . . . , X n = x n ) (4)
In our model, the 9 parameters extracted represents the nodes of the Bayesian Network. Details about application of Bayesian network approach can be found in [START_REF] Barbini | Bayesian Approach in Medicine and Health Management,Current Topics in Public Health[END_REF].

Algorithm 1 summarizes all steps achieved by the application layer of our system.

IV. EXPERIMENTAL RESULTS

In order to evaluate our proposed approach we used 3 medical databases from the Physionet [START_REF] Goldberger | PhysioBank, Phys-ioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals[END_REF]. The European ST-T Database (EDB) used for subjects with Myocardial Infarction, the St. Petersburg Institute of Cardiological Technics 12-lead Arrhythmia Database (INCARTDB) for subjects with PVC and PAC and the MIT-BIH Arrhythmia Database (MITDB) for normal subjects.

Figure 4 shows various ECG signals including normal ECG from MITDB, MI ECG from EDB which is characterized by an elevation of the ST Segment as shown on the figure, PVC and PAC ECG from INCARTDB.

Figure 5 shows the previous signals processed by the 0.5Hz FIR high pass digital filter and UWT. We can see that the resulting ECG signals contain little baseline wandering but keeps the main characteristics of the original ECG signals, we can see also that the wideband noises are suppressed while almost all the information of the ECG signal are maintained. 6 shows the extracted features P onset, P of f et, Rpeak, QRSonset, QRSof f et, T onset and T of f et for two previous ECGs (normal and MI). As described in Section III, these values were obtained by applying UWT to first extract QRS complex peaks and then deduct P and T peaks and their corresponding onsets and offsets.

Figure 7 shows the calculated amplitude parameters P amp, QRSamp, T amp, ST amp for the previous normal, MI, PVC and PAC ECGs during 15 minutes, representing about 1000 ECG cycles or heart beats. As we can see in the figure, the amplitudes are stable enough for normal ECG but vary significantly from one anomaly to another. For example, ECG with MI has ST amp parameter greater than 0.5mV while it is negative or around 0mV for PVC and PAC.

Figure 8 shows the calculated duration parameters P dur, QRSdur, T dur, P Rdur, QT dur for the previous ECGs during the same period (1000 ECG cycles). We can also see differences between normal and abnormal ECG, especially for QT dur and T dur which are almost greater than 0.45s and 0.3s for ECG with MI while are generally less than 0.4s and 0.25s in other cases.

To evaluate the performance of our BN prediction model, we first pre-processed and extracted parameters of all records of the 3 databases EDB, MITDB and INCARTDB. Then, we selected a data set of 4451 ECG consisting of 896 ECG with MI, 1207 ECG with PVC, 1238 ECG with PAC and 1110 normal ECG. Each ECG is characterized by the 9 parameters mentioned above. 5% of the data set was used as a training set (TS) and the remaining 95% as a test set. We are interested for each class to measure the accuracy parameters bellow :

True Positive Rate (Sensibility) = T P T P + F N

False Positive Rate (False Alarm Rate) = F P F P + T N [START_REF] Sung | Livenet: Health and lifestyle networking through distributed mobile devices[END_REF] True Negative Rate (Specificity) = T N F P + T N [START_REF] Luprano | Combination of body sensor networks and on-body signal processing algorithms: the practical case of myheart project[END_REF] Positive Predictive Value (Precision) = T P T P + F P

Where T P , F P , T N , F N are respectively, the number of true positives, false positives, true negatives and false negatives. The table I synthesizes the results obtained for each class. If we increase the size of the TS to 30%, we obtain an average TPR of 98.7% with 0.5% FPR. For a TS of 50%, we obtain an average TPR of 99.3% with 0.3% FPR. We deduct that the size of the TS is important for the performance of the classification, this is why we update the TS of our system with the new ECG data captured.

Finally, we use the Receiver Operating Characteristic (ROC) curve to represent the fraction of True Positive Rate (TPR) vs. False Positive Rate (FPR) with different threshold. Figure 9 shows the ROC curve of each class. The results obtained show that the average rates of T P R and F AR for all the 4 classes are T P R=96.1% and F AR=1.3% respectively.

V. CONCLUSION

In this paper, we present a new WBAN system for ECG monitoring and prediction of cardiac anomalies. We designed our system in five layers and detailed in this paper the application layer that performs three main tasks : pre-process captured ECG signal using a digital high-pass filter and Undecimated Wavelet Transform (UWT) to remove various noises, extract 9 main ECG features and classify ECGs using a Bayesian Network Classifier as normal or abnormal among MI, PVC or PAC anomaly. We applied our proposed approach using 3 medical databases of real ECG data and experimental results show that our system achieves an average True Positive Rate of 96.1% with an average False Alarm Rate of 1.3%. We are currently working on an end-to-end prototype using market sensors to validate our experimental results in a user environment.

Fig. 2 .

 2 Fig. 2. A typical one-cycle ECG (Heartbeat)

  The proposed system is composed of of 5 layers : -A SensorLayer composed of ECG sensors able to capture, sample and transmit ECG wirelessly. Unlike the 12-lead ECG used at hospitals, these sensors are designed for ambulatory applications and generally captures a subset of 2 or 3 leads. -A CommunicationLayer where the digitized ECG is received wirelessly from sensors via bluetooth and alerts are sent to healthcare professionals along with the captured ECG and the position of the patient. These tasks are achieved by a Smartphone playing the role of a gateway. -An ApplicationLayer designed as a mobile application running on a Smartphone that pre-processes the ECG, extracts ECG features (P , Q, R, S, T ) and analyzes the extracted features using a Bayesian Network model, in order to predict the presence of PAC, PVC or MI anomaly. The model is built based on a training set of real ECG data and applied on each new ECG features extracted to classify it and decide whether to raise an alarm or not. This model is updated with the new data whenever a new ECG is received. -A DataLayer containing databases used to save the patient data (general information, medical background, ECG captured, features extracted, alerts raised, etc.

Fig. 3 .

 3 Fig. 3. Design steps of the application layerThe pre-processing stage aims to remove different types of noise present in the digitized ECG received from sensors. These noises overlap with the cardiac components in the frequency

  : -P amp : Amplitude of the P wave of ECG. -P dur : Duration of the P wave of ECG. -QRSamp : Amplitude of the QRS complex of ECG. -QRSdur : Duration of the QRS complex of ECG. -T amp: Amplitude of the T wave of ECG. -T dur: Duration of the T wave of ECG. -P Rdur: Duration of the P R Interval -ST amp: Amplitude of ST Segment. -QT dur: Duration of the QT Interval. These 9 parameters are used as the random variables of our Bayesian Network (BN) Classifier discussed bellow.

Fig. 6 .Fig. 7 .

 67 Fig. 6. ECGs after Features Extraction


Samuel Lellouch

Laurent Sanchez-Palencia

Localization transition in weakly interacting Bose superfluids in one-dimensional quasiperdiodic lattices

Keywords: numbers: 03.75.-b, 05.30.Jp, 05.70.Ln

We study the localization of collective pair excitations in weakly-interacting Bose superfluids in one-dimensional quasiperiodic lattices. The localization diagram is first determined numerically. For intermediate interaction and quasiperiodic amplitude we find a sharp localization transition, with extended low-energy states and localized high-energy states. We then develop an analytical treatment, which allows us to quantitatively map the localization transition into that of an effective multiharmonic quasiperiodic system.

Quasiperdiodic systems, which are formed of a small number of incommensurate sinusoidal components, constitute an appealing intermediate between disordered and periodic systems. Such structures are basic models for a wide variety of physical systems. They appear naturally in the growth of certain crystals [1] or as a result of charge-density waves [2]. They also describe twodimensional lattice electrons in perpendicular magnetic fields [3]. Moreover, they can be created on purpose in solid crystals [4], photonic crystals [5], and ultracoldatom optical lattices [6][7][8]. In quasiperiodic systems, the lack of translation invariance can induce localization of linear waves, similarly as the phenomenon of Anderson localization in disordered systems [9]. In quasiperiodic systems, however, the quasi-repetition of finite patterns radically changes the localization picture. For instance, in a one-dimensional disordered system, any quantum particle is localized with an energy-dependent localization length [10,[START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF]. In contrast, for a quasiperiodic system made of a single incommensurate sinusoidal modulation of a main periodic lattice, there is a localization transition for some critical strength of the quasiperiodic component beyond which the states are localized with a localization length that is independent of the energy [START_REF] Aubry | [END_REF].

The extension of the concept of localization to interacting quantum systems is attracting a considerable attention as regards phase diagrams [13], many-body localization transitions [14], and localization of collective excitations [15][16][17]. These issues have been first investigated for purely disordered systems and extensions to quasiperiodic systems are just starting. So far, most studies focused on the phase diagram of one-dimensional bosons in quasiperiodic lattices at zero temperature [18][19][20], finite temperature [21], and infinite temperature [22]. Conversely, the localization of collective excitations remains largely open. This issue is particularly important because the transport of collective excitations governs many dynamical effects in correlated quantum systems [23], for instance the propagation of correlations in recently-developed quench experiments [24].

Here we study the localization of collective pair excitations in weakly interacting Bose superfluids subjected to a one-dimensional quasiperiodic lattice. We first determine the localization diagram numerically and show that, for intermediate interaction and quasiperiodic amplitude, there is a sharp localization transition. This nontrivial transition separates bands of states that are extended at low energy and localized at high energy. We then develop an analytical treatment, which allows us to reproduce the numerical results accurately and to quantitatively map the localization transition into that of an effective multiharmonic quasiperiodic system. Finally, we discuss experimental observability and possible extension of our results.

The starting point of our study is the Aubry-André-Hubbard Hamiltonian,

Ĥ = - j,l T j,l â † j âl + j V j â † j âj + U 2 j â † j â † j âj âj , (1)
which governs the low-energy physics of interacting bosons in one-dimensional (1D) quasiperiodic lattices. In Eq. ( 1), âj and â † j are the bosonic annihilation and creation operators at the lattice site j. The first term represents quantum tunneling with the hopping matrix T , which includes nearest-neighbor tunneling, T j,j±1 = t and T j,l = 0 for |j -l| > 1, as well as the homogeneous on-site term T j,j = -2t, for convenience. Within this convention, the free-particle spectrum, ε k = 4t sin 2 (k/2), is centered on ε = 2t with the band edges ε = 0 and ε = 4t. The second term represents the on-site quasiperiodic potential modulation, V j = ∆ cos(2πrj + ϕ), where ϕ is a phase, ∆ is the quasiperiodic amplitude, and r is an irrational number. The third term represents on-site repulsive interactions with the interaction energy U > 0.

In the weakly interacting regime with high occupation number per lattice site (n ≫ U/t, with n the mean density), we can rely on mean field theory [25]. A similar approach has been presented elsewhere for disordered systems in continuous [17,26,27] or lattice [28] spaces, and we just outline it here. The density background n j is first determined by minimizing the classical energy functional, obtained by replacing the operator âj by the real-valued field φ j ≡ √ n j in Eq. (1). It yields the Gross-Pitaevskii equation (GPE),

µφ j = -T φ j + V j φ j + U φ 3 j , (2)
where the term T φ j is a shortcut for the hopping matrix contribution t(φ j+1 -2φ j + φ j-1 ) and µ is the chemical potential. The collective pair excitations of the Bose superfluid, which are represented by two fields u j and v j , are then found by expanding Hamiltonian (1) up to second order in the Bogoliubov operator δn j /2 √ n j + i √ n j δ θj , where δn j and δ θj are the density and phase fluctuation operators, and diagonalizing the resulting quadratic Hamiltonian. The excitation energy E and wave functions u j and v j are the solutions of the Bogoliubov-de Gennes equations (BdGEs)

-T +V j -µ+2U n j U n j -U n j T -V j +µ-2U n j u j v j = E u j v j . (3)
Equations ( 2) and ( 3) form the complete set to determine the elementary excitations of the Bose fluid in the quasiperiodic lattice. We first solve Eqs. ( 2) and ( 3) numerically in the 1D quasiperiodic lattice. The number of lattice sites is chosen to be a Fibonacci number F p and r is taken as the ratio F p-1 /F p . It allows us to use periodic boundary conditions and a good approximation of an incommensurate ratio ( √ 5 -1)/2 [START_REF] Aubry | [END_REF]. In practice, we use F p = 610, which yields r = ( √ 5 -1)/2 ± 0.000002 [START_REF]We have checked that the results do not depend on the system size by using other numbers of[END_REF]. The density background is computed by solving the GPE using imaginary-time propagation with a Crank-Nicolson scheme [START_REF] Press | Numerical Recipes, 3rd Edition[END_REF]. The good numerical convergence of the imaginary-time propagation of the GPE is a delicate point for the subsequent determination of the collective excitations using Eq. ( 3). The convergence criterion applies to the effective, imaginary-time-dependent chemical potential µ(τ ) ≡ -2 d dτ ln j n j . We have checked that the density profile is unaffected when the precision threshold varies from 10 -8 to 10 -15 . The same holds when the imaginary time step ∆τ used in the propagation varies from 0.01/t to 0.5/t. Moreover, the density profile precisely agrees with the perturbative expansion of the GPE solution implemented up to order 50 [31]. All together, the precision on the density profile n j is of the order of 10 -8 for all results presented here. The excitations are then computed by exact diagonalization of the matrix in Eq. (3) using the Lanczos algorithm for sparse non-Hermitian eigenproblems [START_REF] Press | Numerical Recipes, 3rd Edition[END_REF].

The numerical results are summarized on the diagram in Fig. 1(a). It displays three different regimes. For weak quasiperiodic amplitude ∆ and strong interaction U , the density background is fully connected and all excitations are extended ('extended regime'). For a given interaction strength U and tunneling t, the density modulations increase with the quasiperiodic amplitude ∆. Above a critical value of ∆ c , the density profile gets fragmented ('fragmented regime'), which yields the upper boundary on the diagram. The fragmentation condition is chosen .

to be the minimal value of ∆ such that at least one lattice site has a density lower than 0.01 atom per site. We have checked that varying this arbitrary threshold down to 0.001 yields insignificant changes of the fragmentation boundary. Moreover the latter is in good agreement with the experimental observation of Ref. [7]. In the fragmented regime, the density profile is cut in disconnected pieces. It corresponds to trivial localization, a case that we disregard in the following. Notice that in the limit U → 0, we recover the critical value ∆ c = 2t, which is the localization transition of the noninteracting Aubry-André model. The most interesting regime appears for intermediate quasiperiodic amplitude ('ext-loc regime').

In this regime, although the density background is fully connected, we find a localization transition of the collective excitations. Remarkably enough, they are the highenergy excitations that are exponentially localized over a few lattice sites [see Fig. 1(b)] while the low-energy excitations are extended over the whole system [see Fig. 1(c)]. This transition is sharp as exemplified in Figs. 1(b) and (c), which correspond to two excitations of consecutive energies for U n = 1.75t and ∆ = 3.3t.

In order to characterize the localization transition, we compute two Lyapunov exponents for the excitations, which correspond to the two Bogoliubov wave functions, γ u (E) ≡ -lim j→∞ ln |u j |/j and γ v (E) ≡ -lim j→∞ ln |v j |/j. They are extracted from fits in the tails of the logarithm of the wave functions u and v. tinguishable and hereafter we omit the wave function index u or v. In the 'ext-loc' regime the Lyapunov exponent curves clearly show the transition, separating extended (γ = 0) and localized (γ > 0) states [32]. The excitation spectrum splits in several bands separated by minigaps, a general feature in quasiperiodic systems [START_REF] Aubry | [END_REF][START_REF] Sokoloff | [END_REF][34][35][36]. The transition generally lies in one of the minigaps. To determine the mobility edge E c , we thus rely on fits of the γ(E) curves with several fitting functionals [linear, γ(E) ∼ E -E c ; power-law, γ(E) ∼ E α -E α c ; and logarithmic, γ(E) ∼ ln(E/E c )]. The result is found to be almost independent of the fitting functional and thus provides a reliable estimate of the mobility edge. Figure 2(b) shows the mobility edge versus the quasiperiodic amplitude for various interaction strengths. The errorlike bars represent the edges of the minigap containing the mobility edge. The uncertainty on the fitted mobility edge is smaller than these bars.

In order to interpret those results, we now turn to an analytical treatment of the localization problem. The main difficulty relies on the fact that localization in quasiperiodic systems occurs for strong quasiperiodic amplitude ∆ [START_REF] Aubry | [END_REF]. For this reason, the lowest-order perturbation theory, which proved successful for 1D disordered systems [16,17], fails here [START_REF]We have found that lowest-order perturbation theory as used for 1D disordered, interacting Bose gases in Refs[END_REF]. To overcome this issue, we develop an approach based on a generic expansion in harmonics of the quasiperiodic potential. The structure of the GPE (2) shows that the field φ j takes the form of a series of harmonics of the quasiperiodic potential. The density field n j thus reads n j = (µ -Ṽj )/U , where

Ṽj = p≥1 A p cos[p(2πrj + ϕ)] (4)
is a multiharmonic quasiperiodic field, the coefficients of which can be computed iteratively [31]. Using the energy-dependent linear transform [17] g

± j = ±ρ ±1/2 E (u j -v j ) + ρ ∓1/2 E (u j + v j ) ( 5
)
where ρ E = 1 + (µ/E) 2 + µ/E, the BdGEs (3) exactly rewrite

-(ρ -1 E E + T )g + j + V j - 3 + ρ 2 E 1 + ρ 2 E Ṽj g + j = 2ρ E Ṽj 1 + ρ 2 E g - j (6) (ρ E E -T )g - j + V j - 1 + 3ρ 2 E 1 + ρ 2 E Ṽj g - j = 2ρ E Ṽj 1 + ρ 2 E g + j . ( 7
)
The solution of these equations is significantly simplified by noticing that the lattice-space Green function of the operator -T + ρ E E is of width t/ρ E E and amplitude 1/(ρ

E E + 2t) 1 -[2t/(ρ E E + 2t)] 2 .
Hence, for ρ E E ≫ t, this operator can be replaced by the local operator ρ E E+2t in Eq. (7). It is then straightforward to write the expression of g - j as a function of g + j and of the potentials V j and Ṽj . Inserting this expression into Eq. ( 6) we find a closed equation for g + j ,

-T g + j + V E j g + j = Eρ -1 E g + j , (8)
with the effective potential

V E j ≃ V j - 3 + ρ 2 E 1 + ρ 2 E Ṽj - 2ρ E 1 + ρ 2 E 2 Ṽ 2 j ρ E E + 2t + V j - 1 + 3ρ 2 E 1 + ρ 2 E Ṽj . (9)
Using exact diagonalization of Eq. ( 8) with the potential (9) around energy E, we have checked that the Lyapunov exponents and the localization transition given by our effective model coincide with those found using direct diagonalization of the BdGEs (3). It validates the effective model ( 8)-( 9) and the approximation

-T + ρ E E ≃ 2t + ρ E E used above.
In this model, the quantity V E j is a multiharmonic periodic potential of spacing 1/r incommensurate with that of the main lattice, which is unity. Such systems are known to exhibit in general an energy-dependent mobility edge with low-energy extended states and high-energy localized states [START_REF] Sokoloff | [END_REF][34][35][36]. This holds except in the particular case of self-dual models, among which the Aubry-André model is a celebrated example [START_REF] Aubry | [END_REF]. Self-duality requires a specific relation between the amplitudes of the pth harmonics and of the tunneling rate to the pth neighbors. The latter does not apply in our case since tunneling is strictly restricted to the first neighbors. It qualitatively explains the localization transition of the collective excitations reported here.

Localization properties in quasiperiodic systems can be further inferred from locator perturbation theory [START_REF] Sokoloff | [END_REF].

Here the localization criterion roughly corresponds to the convergence of the self-energy in the thermodynamic limit, which reads D(E) > 1 where D(E) is the so-called localization function. In the case of Eq. ( 8), it reads

D(E) = exp r 1/r 0 dx ln Eρ -1 E -2t -V E (x) t . (10)
Equation ( 10) can, in principle, be applied to the full effective potential V E (x). To obtain analytical results, it is, however, worth truncating the infinite series of harmonics in V E . Keeping only one harmonic is not sufficient to capture the physics even qualitatively, since it would unphysically restore duality and change the universality class of the localization transition. On the other hand, beyond two, the number of harmonics does not change the universality class. We may thus restrict ourselves to the two lowest-order harmonics, which are generated in first instance in second-order perturbation theory. As we shall see, this order of expansion turns out to be sufficient for a significant part of the localization diagram of Fig. 1. It yields the effective two-harmonic potential

V E j ≃ ∆ (0) E + ∆ (1)
E cos(2πrj + ϕ) + ∆

E cos[2(2πrj + ϕ)] (11) with the amplitudes

∆ (0) E = 3+ρ 2 E 1+ρ 2 E ∆ 2 4U n (f r -f 2 r ) - 2ρ 2 E (1+ρ 2 E ) 2 ∆ 2 f 2 r ρ E E +2t (12)
∆

E = ∆ 1 - 3 + ρ 2 E 1 + ρ 2 E f r (1)

∆

(2) It predicts the correct localization transition with collective excitations that are extended at low energy and localized at high energy, and a mobility edge that is in very good agreement with the full numerical result. Perturbation theory beyond second order generates high-order harmonics in the effective potential and renormalizes the amplitudes ∆ (p) . For U n 3t, we find that they induce negligible effects and do not significantly affect the prediction for the mobility edge. For higher values of U , however, second-order perturbation theory is not sufficient to accurately estimate the density background, and higher-order terms should be included.

E = 3 + ρ 2 E 1 + ρ 2 E ∆ 2 4U n f 2 r 2 + f r - 3 2 f 2 r f 2r - 2ρ 2 E (1 + ρ 2 E ) 2 ∆ 2 f 2 r ρ E E + 2t . ( 14
In summary, we have shown that the collective excitations of lattice Bose superfluids subjected to a singleharmonic quasiperiodic potential undergo a nontrivial localization transition with extended low-energy states and localized high-energy states. Therefore the interactions change the universality class of the localization transition, in striking contrast with the purely disordered case [16,17]. In the quasiperiodic case the transition can be understood as the result of the scattering of the excitations from the potential and the density background, which contains an infinite series of harmonics of the potential. It could be observed in ultracold-atom experiments, using for instance spectroscopy techniques, which give direct access to the excitations [START_REF] Stenger | [END_REF], or in quench experiments, which generate collective excitations that govern the propagation of experimentally-observable correlations [24,39]. It could also be observed in photonic crystals, which can combine quasiperdiodic structures [5] and photon nonlinearities [40]. Finally, it would be interesting to study the counterpart of the localization transition discussed here in Fermi superconductors, which may directly apply to electronic quasicrystals.

Expansion of the effective potential V E j

The effective potential

V E j ≃ V j - 3 + ρ 2 E 1 + ρ 2 E Ṽj - 2ρ E 1 + ρ 2 E 2 Ṽ 2 j ρ E E + 2t + V j - 1 + 3ρ 2 E 1 + ρ 2 E Ṽj , (S16)
which appears in the Eq. ( 8) of the paper, can as well be expanded in powers of the quasiperiodic amplitude by expanding the denominator and using the previously obtained expansions for Ṽj . Up to second order, we get

V E j ≃ V j - 3 + ρ 2 E 1 + ρ 2 E Ṽ (1) j - 3 + ρ 2 E 1 + ρ 2 E Ṽ (2) j - 2ρ E 1 + ρ 2 E 2 ρ E E + 2t Ṽ (1) j 2 , ( S17
)
where Ṽ (1) and Ṽ (2) are given in Eqs. (S11) and (S14). It yields the two-harmonic effective potential

V E j ≃ ∆ (0) E + ∆ (1)
E cos(2πrj + ϕ) + ∆

E cos[2(2πrj + ϕ)] (S18)

with the amplitudes ∆

(1)

E = ∆ 1 - 3 + ρ 2 E 1 + ρ 2 E f r (S19) ∆ (2)
E = 3 + ρ 2 E 1 + ρ 2 E ∆ 2 4U n f 2 r 2 + f r - 3 2 f 2 r f 2r - 2ρ 2 E (1 + ρ 2 E ) 2 ∆ 2 f 2 r ρ E E + 2t (S20) ∆ (0) E = 3 + ρ 2 E 1 + ρ 2 E ∆ 2 4U n (f r -f 2 r ) - 2ρ 2 E (1 + ρ 2 E ) 2 ∆ 2 f 2 r ρ E E + 2t . (S21)

Figure 1 .

 1 Figure 1. (Color online) Numerical results. (a) Localization diagram as a function of the interaction strength and the quasiperiodic amplitude. It displays three regimes: (i) 'extended regime' where the density background is connected and all excitations are extended; (ii) 'fragmented regime' where the density background is fragmented; and (iii) 'extended-localized regime' where the density background is connected and the excitation spectrum shows a delocalization-localization transition with exponentially localized high-energy states and extended low-energy states. (b)-(c) Typical excitation wave function u in the localized (b) and extended (c) regimes, plotted in semilogarithmic scale and for the 150 first lattice sites (similar plots are found for the v wave functions). The two panels correspond to two excitations with consecutive energies above (b) and below (c) the mobility edge for U n/t = 1.75 and ∆/t = 3.3..

  Figure 2(a) displays those Lyapunov exponents versus the excitation energy E, for fixed interaction and disorder strengths. The Lyapunov exponents γ u and γ v are indis-

Figure 2 .

 2 Figure 2. (Color online) (a) Lyapunov exponents of the Bogoliubov wave functions u and v, for U n/t = 1.75 and ∆/t = 3.3, 3.5, 3.7. The excitation spectrum is banded and displays a sharp localization transition separating extended (γ = 0) and localized (γ > 0) states. (b) Mobility edge as a function of the quasiperiodic amplitude ∆/t as extracted from power-law fits to the numerical γ(E) curves [solid lines on panel (a)]. Error-like bars correspond to the edges of the minigap containing the mobility edge. The dotted, red line shows the analytical prediction of the locator theory applied to the effective model (8) with the potential (11).

) and f r = 1 1 +

 1 2t sin 2 (πr)/U n [31]. The results of the locator perturbation theory applied to the two-harmonic potential (11) is shown in Fig 2(b).

We are grateful to Giovanni Modugno and Massimo Inguscio for insightful discussions. This research was supported by the European Research Council (FP7/2007-2013 Grant Agreement No. 256294), the Ministère de l'Enseignement Supérieur et de la Recherche, and the Institut Francilien de Recherche sur les Atomes Froids (IFRAF). We acknowledge the use of the computing facility cluster GMPCS of the LUMAT federation (FR LUMAT 2764) and HPC resources from GENCI-IDRIS (Grants 2013057143 and 2014057143).

Supplemental Material

This supplemental material aims at providing some details about the harmonic structure of the density profile and the effective potential V E j , as well as the series expansion used to determine them.

Series expansion of the density background

To determine the density background n j ≡ φ 2 j , we solve the Gross-Pitaevskii equation (GPE) together with the normalization condition n = 1 L Σ j n j where L is the number of sites and n is the averaged density. To do so, we perform a series expansion in powers of the quasiperiodic potential V j . In the absence of an external potential, we have φ j =

√ n and the chemical potential µ = U n. In the presence of an external potential, we then write

)

where the superscripts denote increasing orders in the quasiperiodic amplitude ∆, and φ

Notice that the chemical potential has to be expanded also to fulfill the normalization condition. Inserting the expansions (S1) and (S2) in the GPE [Eq. (2) of the paper] and the expansion (S3) in the normalization condition, we get

Then, collecting all the terms of same order p in the quasiperiodic amplitude yields

Equations ( S6) and (S7) can then be used to compute all φ (p) j and µ (p) at any order p iteratively. The iteration process works as follows. Given all φ (k) j and µ (k) at orders k < p, we calculate φ (p) j as a function of µ (p) from Eq. (S6) by inverting the operator 1 -T /2U n. The quantity µ (p) is then found by inserting this expression for φ (p) j into Eq. (S7).

Having determined φ (p) j , we then find the density field using Eq. (S3), whose expansion in powers of the quasiperiodic amplitude writes

This procedure is completely general and can be applied to any external potential V j . In the case of the quasiperiodic potential V j = ∆ cos(2πrj+ϕ), the above iterative process is fully algebraic because the operator 1-T /2U n in Eq. (S6) can be analytically inverted at any order (see below). We have implemented this expansion up to order 50 and found excellent agreement with the direct numerical solution of the GPE (2). It provides a cross-check of the precision of the numerical solution and of the convergence of the present analytical expansion.

Analytical expansion in the case of a quasiperiodic potential

We now give some explicit formulas for the lowest order terms and discuss the harmonic structure of the density background n j . As in the paper, we generically write the density field n j n j = (µ -Ṽj )/U, (S9)

where the field Ṽ includes terms of all orders.

At first order, Eq. (S6) reduces to -T φ

(1)

= -V j + U nµ (1) . It is straightforward to solve it in Fourier space where the operator T is diagonal. It yields φ

where ε 0 k = 4t sin 2 (k/2). Inserting this expression into Eq. (S7), we find µ (1) = V k=0 /U n = 0 and φ

. Remarkably, since the quasiperiodic potential contains only one spatial frequency, V k = ∆(e iϕ δ k,+2πr + e -iϕ δ k,-2πr )/2, one can immediately get back to real space and write

where

. Hence, to lowest order, the density profile is quasiperiodic field. It follows the modulations of the quasiperiodic potential with a reduced amplitude since f r < 1. The factor f r is a remainder of the nonlocal operator 1 -T /2U n in the l.h.s. of Eq. (S6), which reduces to an algebraic operation in the case of a quasiperiodic potential. Then, Eqs. (S8) and (S9) yield the first order term of the field Ṽj . It reads Ṽ

where φ

(1) j is given by Eq. (S10), i.e.

The next orders are found following the same process, which remains algebraic to any order in the case of the quasiperiodic potential. To second order, it yields the term

and a negative shift on the chemical potential,

The field Ṽj is then given at second order by Ṽ (2)

)] where φ (1) , φ (2) and µ (2) are given by Eqs. (S10), (S12), and (S13), i.e.

Hence, the second-order terms φ More generally, it is straightforward to show recursively that the terms of order p, φ (p) j and Ṽ (p) j , contain the pth harmonics of the quasiperiodic potential, cos[p(2πrj + ϕ)], as well as all lower harmonics of same parity. In particular, a constant term in φ j and a correction to the chemical potential µ appear only at even orders. Hence, the field Ṽj takes the multiharmonic quasiperiodic form Ṽj = p A p cos[p(2πrj + ϕ)],

(S15)

where the amplitude A p of the pth harmonics is a power series of order p, A p ∼ α p (∆/2U n) p + α p+2 (∆/2U n) p+2 + ....