Alexis Bonnin

Clément Diboune

Nassim Zahzam

Yannick Bidel

Malo Cadoret

Alexandre Bresson

New concepts of inertial measurements with multi-species atom interferometry

come

Introduction

Cold atom interferometry has allowed in the last decades the development of extremely sensitive and accurate inertial sensors for measuring the gravity acceleration [1], Earth's gravity gradient [2] or rotations [3]. They appear very promising for a wide range of applications such as inertial navigation [4], geodesy [5], natural resource exploration [6], or fundamental physics [7][8][9][10][11][12][13]. To address most of these applications, an important research effort is made to make these cold atom instruments more adapted to withstand operational constraints [14,15,19]. In addition to improve these instruments in terms of compactness and robustness, it is of major concern to extend their dynamic measurement range. When shotto-shot acceleration variations are small compared to the dynamic measurement range, defined here as one half of the interferometer fringe spacing, cold atom gravimeters allow to retrieve the gravity acceleration g with a resolution of ≈ 10 -9 g at a repetition frequency of few Hz. This dynamic measurement range lies in the range 0.2 -3 µg for state-of-the-art cold atom gravimeters [16][17][18],

Correspondence to: nassim.zahzam@onera.fr much less than the acceleration level of roughly 50 µg that could be encountered in a typical laboratory environment, 50 mg that could be encountered in a plane [15] and 0.3 g that could be encountered in a ship [19]. To allow the operation of these atomic sensors in a laboratory environment, they are typically used on a passive [20,21] or active vibration isolation platform [22][23][24] or set directly on the ground by combining the atomic signal with that of a conventional accelerometer [15,19,25,26]. To our knowledge, only few onboard demonstrations of inertial measurements have been reported whether in an elevator [21], on a plane [15], or on a boat [19]. For the elevator demonstration, the interrogation time of the atom interferometer is reduced significantly to increase the fringe spacing, and therefore the dynamic measurement range, but at the cost of a drastic loss of the sensitivity measurement. For the plane and boat demonstrations, the atom interferometer is coupled to a conventional accelerometer which helps to determine, for each measurement cycle, the fringe index on which the atom interferometer is operating. In this last configuration, the use of a conventional accelerometer often limits the performance of the instrument either because of the lower performance intrinsic noise of the conventional accelerometer, its non-linearity, the lack of knowledge of its transfer function or because of a mismatch between the localization of the conventional accelerometer measurement point and that of the atom interferometer.

As an alternative or in complement to vibration isolation platforms and conventional accelerometers, we propose and study in this article new concepts of an atomic accelerometer dedicated to onboard applications. These concepts could be interesting in the future for increasing the dynamic measurement range of the atomic instrument and allowing the atomic interferometer to operate always in its linear regime [25], leading to a maximum sensitivity whatever the measured acceleration. The instrument configurations that will be presented do not rely on external devices making our instrument selfindependent and only relying on atom interferometry. This work is in straight line with the one reported earlier by Sorrentino et al. [27], where two separated atomic clouds of 87 Rb are simultaneously operated in a gravity gradiometer configuration, that opened the way to this kind of self-independent atomic sensors. Here, we study innovative concepts based on atom interferometry with different atomic species which allows simultaneous operation of several atomic clouds at the same spatial position. Up to now, the few experiments dealing with multi-species atom interferometry were dedicated to test the universality of free fall [28][29][30][31][32][33] and their interests to improve single species atomic inertial measurements were never discussed to our knowledge. The principle of the experiment is the following. The atomic species are trapped and cooled simultaneously in the same vacuum chamber, sharing the same laser beams and constituting a whole unique embedded atomic cloud. The atomic species are then free falling simultaneously, feeling therefore the same inertial forces, and during which they are interrogated through atom interferometry by the same laser beam. At the end of the free fall, each atomic species is detected and provide complementary acceleration measurements that allow to improve performances reached by standard single species atom interferometers. Using different atomic species in the same vacuum chamber is of great interest since each atom-light interaction for manipulating one atomic species should not disturb significantly the other atomic species. Indeed each laser line dedicated to manipulate one atomic species is seen far of resonance regarding the other atomic species.

2 Presentation of the multi-species concepts

Introduction

In a general case, the phase shift ∆Φ accumulated between both paths of a Mach-Zehnder type interferometer is proportional to the acceleration of the atoms a along the interferometer's laser direction of propagation [1]:

∆Φ = k eff • a • T 2 (1)
The interferometer is characterized by three laser pulses of effective wave-number k eff ≈ 4π λ , where λ is the wavelength of the laser, and separated by the interrogation time T . Experimentally, the measured quantity at the output of the interferometer is the transition probability P to be in a particular atomic state which is a sinusoidal function of the interferometric phase:

P = P 0 - A 2 cos (∆Φ + ∆Φ op ) (2)
where P 0 and A are respectively the offset and the amplitude of the fringes. ∆Φ op is a very well known and controlled phase shift added by the operator. For instance, in an atom gravimeter, a frequency chirp α is added to the interferometer's laser frequency to compensate the Doppler shift induced by gravity g such as

∆Φ op = 2παT 2 = -k eff • g • T 2 .
The sinusoidal response of the atom interferometer is responsible of two important drawbacks when the phase shift ∆Φ is extracted.

First, when shot-to-shot variation of the acceleration signal induces phase shifts larger than π, then the measured transition probability P could correspond to a large number of acceleration values (see Fig. 1):

â = s k eff T 2 arccos 2 P 0 -P A - ∆Φ op k eff T 2 + n 2π k eff T 2 (3)
where n ∈ Z and s = ±1 correspond to the fringe indexes of the atom interferometer at each cycle. For onboard applications, a conventional accelerometer is typically used to determine these fringe indexes at each cycle [15,19] or could be used in a real time phase compensation loop [26].

Secondly, the sensitivity of the interferometer depends on the value of the total phase shift and consequently on the measured acceleration. When ∆Φ + ∆Φ op = nπ (n being an integer), i.e when the interferometer works at the top or at the bottom of the fringes, the sensor becomes highly insensitive because the derivative of the measured probability P with respect to ∆Φ is zero (see Fig. 1). In a low vibration environment, cold atom gravimeters are usually forced to operate at mid-fringe by adjusting ∆Φ op at each cycle to ensure a maximum sensitivity of the instrument [1]. For vibration levels exceeding the measurement range of the atom interferometer, a solution consists in compensating the vibrations in real time thanks to a conventional accelerometer allowing the atom interferometer to stay at mid-fringe [26]. In this article, we propose two configurations to mitigate these two issues.

First concept

The first configuration relies on the implementation of two simultaneous atom interferometers with two different interrogation times T 1 and T 2 [see Fig. 2(a)]. The atom interferometer with the highest interrogation time T 2 presents a greater sensitivity but conversely a smaller measurement range. On the other hand, the atom interferometer with a lower interrogation time T 1 benefits from an increased measurement range, by a factor (T 2 /T 1 )

2 , but at the expense of a degraded sensitivity, following the same scaling factor. Associating both interferometer signals allows to benefit simultaneously from a high sensitivity and a high measurement range. Note that, for a given measurement range fixed by the atom interferometer with T 1 , this configuration allows also to benefit from a higher accuracy thanks to the interferometer with the highest interrogation time T 2 . Indeed, the impact of a significant number of systematic effects scales inversely with the interrogation time.

Second concept

The second proposed configuration relies on the implementation of two simultaneous atom interferometers operating in phase quadrature [see Fig. 2(b)]. The linear ranges from both interferometers can be coupled and a fully linear response gravimeter is obtained allowing to get rid of insensitive acceleration regions where |dP/d∆Φ| ≈ 0.

Experimental demonstrations

Experimental setup

We report here preliminary experimental demonstration of these two configurations. The experimental setup is similar to the one described in Ref. [34]. With this setup, approximatively 6 × 10 8 , respectively 8 × 10 8 , atoms of 87 Rb, resp. 85 Rb, are simultaneously loaded from a background vapor into a 3D Magneto-Optical Trap (MOT) in 250 ms. The atoms are then further cooled down in an optical molasses phase allowing the atoms to reach a temperature of ∼ 2 µK. The atoms are then selected in the magnetic insensitive state |F = 1, m F = 0 for 87 Rb and |F = 2, m F = 0 for 85 Rb thanks to a microwave pulse and then dropped over a distance of ∼ 6 cm. For both atomic species, a Mach-Zehnder type interferometer is implemented, consisting in three equally spaced Raman laser pulses π/2 -π -π/2. Note that the interrogation times can be varied for each species to a maximum value of T = 47 ms. The Raman laser beam is retro-reflected on a mirror which acts as the reference for the inertial measurement. Since both atom interferometers share the same Raman laser beam, the two measurement axis are therefore rigorously the same. Finally, the atomic probabilities are successively measured for each isotope by fluorescence detection. The whole sequence is performed at a repetition rate of 4 Hz. The instrument is installed on a vibration isolation platform to clearly resolve the atomic fringes and minimize the interferometer phase noise which will allow us to analyze more simply the potential of the presented concepts. Despite an apparent complexity, the dual species experiment stays rather simple to implement and not significantly more complex than a single species instrument. All the required laser lines necessary for both isotopes manipulation are synthesized thanks to phase modulation of a single laser source [34]. During the cooling stage, the carrier frequency is tuned on the 87 Rb cooling transition. The 85 Rb cooling line and both repumper lines are generated by three modulation frequencies injected into the phase modulator as following: f 85 cooling = f carrier + 1.126 GHz, f 87 repumper = f carrier + 6.568 GHz, f 85 repumper = f carrier + 1.126 + 2.915 GHz. Phase modulation is also used for generating the laser lines during the interferometric sequence. Both Raman pairs are generated by directly injecting the Raman difference frequencies associated to both isotopes (i.e. 6.834 GHz for 87 Rb and 3.035 GHz for 85 Rb), making the carrier frequency common to both Raman pairs. The Raman pair corresponding to 87 Rb is red detuned by 0.59 GHz with respect to the excited hyperfine state F = 2 and therefore the one corresponding to 85 Rb is red detuned by 1.86 GHz with respect to F = 3. Compared to the one species instrument, the dual species instrument requires only one additional microwave source fitting in a standard 3U 19 inch rack.

First concept

The first dual-species operating mode that we study in this article allows to increase the measurement range of a single species atom accelerometer. For the experimental realization of this operating mode, we set the interrogation time for 87 Rb to T 87 = 20 ms and the one for 85 Rb to T 85 = 47 ms. The central π-pulse is kept common for both isotopes (see Fig. 3). This feature allows to keep the global symmetry of the entire interferometric sequence in order to cancel out roughly impact of light shifts on one atomic species due to the laser lines used to address the other isotope.

The dual atom interferometer output should be impacted by light shifts, first in the same way as the single species interferometers, due to the interaction of the atoms with their associated Raman laser (Type I), and secondly, by the interaction of the Raman laser of one atomic species with the other atomic species (Type II). In this article, we focus on proof-of-principle experimental demonstrations of the presented concepts without studying light shift or other potential effects that could impacts the dual interferometer accuracy or stability. Nevertheless, the one photon light shift of type I and II should be efficiently rejected by the well-known reversing technique of the Raman laser wave vector [1]. For the concept of interferometers with two different interrogation times, light shift of type II should be further eliminated due to the symmetry of the embedded Mach-Zehnder schemes, sharing the same central pi Raman pulse. For this concept, the interferometer with lower T is not impacted since the π/2 laser pulses of the other atomic species are temporally outside the interferometer and the central pi pulse should impact both paths of the interferometer in the same way. The interferometer with higher T should be impacted only weakly since the effect of the first π/2 pulse on both interferometers paths should be canceled by the effect of the last pi/2 pulse. Concerning the two-photon light shift, the two-photon detuning of is too large to expect significant contributions at least regarding terrestrial applications.

Moreover, in such configuration, the effective measurement points of both interferometers are temporarily In order to highlight the interest of such configuration to allow larger acceleration variation measurements, we have extracted the value of α from the two experimental interferometric signals thanks to Eq.( 3):

α = - k eff a 2π + s 2πT 2 arccos 2 P 0 -P A + n T 2 (4)
The extracted value of α corresponds to α and a can be considered constant and equal to the gravity acceleration g. Note that applying this frequency chirp α is equiva- from the frequency noise of the laser beam, which could lead to a bad determination of the 85 Rb interferometer fringe indexes, as it's shown on Fig. 5. For the 85 Rb interferometer noise, the main contribution comes from the microwave phase noise allowing generation of the second Raman laser line. In Fig. 5, n 85 [Fig. 5(a)] and s 85 [Fig. 5(b)] correspond to the fringe indexes of the 85 Rb signal used in Eq. 3 as n and s. The green circles give for each measurement point the values of the fringe indexes of the 85 Rb signal to determine unambiguously α. These indexes are determined by using the 87 Rb signal. To evaluate the efficiency of this determination and the impact of the 87 Rb interferometer noise on it, we have plotted on Fig. 5 in orange line the fringe indexes values that would be determined with a noiseless 87 Rb signal, using a sine function fitting the 87 Rb output.

We can see that the determination of the fringe indexes is not perfect which limits the gain in sensitivity. In the ideal case of a noiseless 87 Rb signal, the sensitivity gain reaches a factor 3.8, which is still lower than what could be expected from the ratio (T 85 /T 87 ) 2 ∼ 5.5 coming from the T 2 dependance of the atom interferometer scale factor. This reduced gain results mainly from measurement points near the fringe extrema where the atom interferometer sensitivity is significantly reduced. The noise on the experimental data is here a strong limitation of this new concept and it prevents us from further reducing the interrogation time of the 87 Rb interfer-ometer to increase the measurement range. Indeed, the shorter the interrogation time, the worse the determination of the fringe index: when T is decreased, the impact of the atomic signal noise is increased with an acceleration noise scaling as T -2 . We can see for instance on Fig. 6 the poor efficiency of the fringe indexes determination with T 85 = 47 ms and T 87 = 10 ms due to the noise of the 87 Rb signal. In this case, there is no sensitivity gain with the dual atom interferometer compared to a single species interferometer while, in the ideal situation of a noiseless 87 Rb signal for the determination of the 85 Rb fringe indexes, the expected gain sensitivity thanks to a dual species interferometer reaches now a factor 17. Even removing measurement points near the extrema does not lead to an improvement of the sensitivity gain because the noise of the 87 Rb signal leads now to determination errors exceeding one fringe of the 85 Rb signal.

We performed numerical simulations to deduce the influence of the Signal-to-Noise Ratio (SNR) on the sensitivity gain provided by the dual atom interferometer, with T 85 = 47 ms and T 87 = 20 ms, compared to a single species atom interferometer with T 87 = 20 ms. For the simulation, we have considered noise affecting only the amplitude of the fringes (e.g. detection noise), equal SNR for the two different atom interferometers and a random gaussian noise σ P centered on 0 directly added to the transition probability. The SNR is defined as A/σ P where A is the fringe amplitude defined in Eq.2.

In our experiment, the SNR is typically of ∼ 50. The range of α variation is taken similar to the one chosen for the experimental demonstration with T 85 = 47 ms and T 87 = 20 ms. As for the experimental results, in the dual atom interferometer, the 87 Rb signal is only used to determine the fringe indexes of the 85 Rb interferometer. The simulation results are reported on Fig. 7. We find no gain for low SNR, meaning that the atom interferometer with low interrogation time provides equivalent performance compared to the dual atom interferometer. However, as the SNR is further increased, the dual AI provides better performance reaching an optimal gain of roughly a factor 3.8 for a SNR ∼ 350. For higher SNR, the gain decreases. This behavior is attributed to the presence of fringe extrema in the 85 Rb atom interferometer signal with higher interrogation time. Extraction of α near those fringe extrema lowers the global sensitivity of the dual atom interferometer whereas the single atom interferometer, with a lower interrogation time, has no fringe extrema and keeps then a higher global sensitivity. Note that for the experimental SNR of ∼ 50, the simulations foresee a sensitivity gain of ∼ 2.5, in agreement with that obtained experimentally.

Second concept

In the more specific issue of sensitivity loss due to measurement points near fringe extrema, we have also studied a second operating mode which relies on two perfectly simultaneous interferometers working in phase quadrature. Both effective Raman wave-vectors need to be the same in order to guarantee equal scale factors and thus equal fringe spacings. For the implementation of this configuration, a differential phase φ d = π/2 is experimentally created by adding during the π-pulse a controlled phase jump on the microwave source used to generate the Raman laser pair of the 85 Rb. Note that for this concept, the time sequence of Raman laser pulses is the same for both species. The effect of light shift on one atomic species due to laser dedicated to the other atomic species should be efficiently rejected by reversing the sign the Raman laser wavevector.

Fig. 8 shows the two normalized interferometric signals in phase quadrature obtained with 87 Rb and 85 Rb for T = 47 ms. Whatever the acceleration undergone by the atoms, corresponding to a given value of α, there is always one of the two interferometers which is working in its linear range (blue-grey circles). Considering a fringe amplitude A = 1, this linear range corresponds to transition probabilities P such as 2-

√ 2 4 < P < 2+ √ 2
4 . Selecting from both atom interferometer outputs and at each cycle the measurement corresponding to this linear range gives birth to a fully linear sensor. We can also notice that the correlation between both signals leads to extend the measurement range by a factor 2 compared to a single-species sensor. This is clearly visible on Fig. 9 showing the extraction of α for a one atomic species interferometer (red circles) and for the phase shifted dual species interferometer (blue-grey circles). In order to estimate the gain in sensitivity provided by this configuration compared to a single-species sensor, we performed a numerical simulation. The interferometric phase is uniformly distributed on [0, π] and a random gaussian noise centered on 0 is directly added to the transition probability. The interferometric phases are then extracted thanks to Eq.(3).

The gain in sensitivity is given in TABLE 1. We notice that the higher the Signal-to-Noise Ratio, the higher the gain in sensitivity is. Indeed, when the SNR increases, i.e when the noise decreases for a given fringe amplitude, the uncertainty on the acceleration estimation decreases much more for points at mid-fringe than for points at the top or at the bottom of the fringes. Typically in our experiment the SNR is ∼ 50, and we expect then a sensitivity gain of ∼ 2.7. Putting this analysis into perspective, we can assume higher performance atom interferometer exhibiting SNR approaching 10 4 [35] for which a significant gain in sensitivity, higher than a factor 8, could be reached using this phase shifted dual species interferometer mode.

Conclusion

In conclusion, we have presented and made first analysis of new concepts of cold atom inertial sensor using multi-species atom interferometry. These original configurations allow an improvement of the measurement range or sensitivity compared to standard configurations of a single species atomic accelerometer. For the first configuration, involving interferometers with different interrogation times, the low interrogation time atom interferometer noise limits clearly the achievement of ambitious performances. The second configuration remains very promising to get a fully linear atom interferometer and should be analyzed experimentally more deeply. Note that even if the use of a classical accelerometer is still necessary for the operation of a cold atom accelerometer, this last configuration seems still very interesting to implement. Indeed, this would improve the sensitivity of the global instrument and mitigate the required performance noise of the classical accelerometer. Manipulating simultaneously different atomic species in the same instrument offers new perspectives to be further studied. Using three atomic species such as 87 Rb, 85 Rb and 133 Cs [36] would allow for instance to implement simultaneously the two concepts presented in this article, to develop an inertial sensor without deadtimes [37] or a simultaneous multiaxis inertial sensor. A strong point of using different atomic species in the same instrument compared to a single species experiment is the ability to increase the number of complementary measurements without disturbance that could arise form spontaneous emission, the lasers manipulating one species being far of resonance relative to the other species.

Fig. 1

 1 Fig.1Representation of a typical output signal from an atom accelerometer. The measured signal P corresponds to a large number of possible acceleration values a. The interferometer operates at its maximum sensitivity in a restricted linear range, away from the fringe extrema where |dP/d∆Φ| = 0. Considering an intrinsic noise of the interferometer σP , the associated acceleration noise when operating on the fringe extrema σ a is larger than the associated acceleration noise when operating in the linear regime σa.

Fig. 2

 2 Fig. 2 (a) Representation of the output signals from two atom interferometers operating with two different interrogation times T1 and T2 = 3T1. The atom interferometer operating with T1 benefits from a 9 times higher measurement range. (b) Representation of the output signals from two atom interferometers in phase quadrature. Whatever the acceleration, at least one atom interferometer is operating in the linear regime.

Fig. 3 3 .

 33 Fig. 3 At the top, the global interferometric time sequence for interrogating both isotopes. At the bottom, experimental interferometric signals for two simultaneous interferometers with two different interrogation times: T87 = 20 ms (black squares) and T85 = 47 ms (red circles). The interferometric fringes are scanned thanks to a frequency chirp α applied on the Raman laser.

Fig. 4 Fig. 5 For

 45 Fig. 4 (a) Non-ambiguous extraction of α from the experimental signals of Fig.3 where T87 = 20 ms and T85 = 47 ms. The frequency chirp, corresponding to a constant acceleration applied on the instrument, can be extracted from interferometric signals of 87 Rb alone (black squares), of 85 Rb alone (red circles) or from the hybridization of both signals (bluegrey circles). The true value of the applied frequency chirp is given in green line. The 85 Rb signal alone is not sufficient to extract α unambiguously. (b) Residue coming from the comparison of the extracted value of α with its true value. Standard deviation of the residue coming from the dual species signal (blue-grey circles) is 2.6 times smaller than the one of residue coming from the single species 87 Rb interferometer (black squares).

Fig. 6 For

 6 Fig. 6 For T85 = 47 ms and T85 = 10 ms, determination of the 85 Rb signal (red circles) fringe indexes n85 [green circles in (a)] and s85 [green circles in (b)] with the help of the 87 Rb signal (black squares). In orange line, the same fringe indexes have been determined using a noiseless 87 Rb signal corresponding to the fitted curve of the 87 Rb signal (black line). (c) The residue coming from the comparison of the extracted value of α with its true value. Standard deviation of the residue coming from the dual species signal (blue-grey circles) is roughly equal to the one of residue coming from the single species 87 Rb interferometer (black squares).

Fig. 7

 7 Fig. 7 Estimation of the sensitivity gain provided by a dual atom interferometer, with T85 = 47 ms and T87 = 20 ms, compared to a single atom interferometer with T87 = 20 ms. The sensitivity gain is estimated relative to the SNR of the atom interferometers where we assume SNR85=SNR87=SNR.

Fig. 8

 8 Fig. 8 Experimental realization of a fully linear atom interferometer. Coupling the 85 Rb (red circles) and the 87 Rb (black squares) interferometer outputs phase shifted by π/2, with T85 = T87 = 47 ms and selecting the adequate signal in each phase range gives rise to a linear response sensor (blue-grey circles) with an constant optimal sensitivity. Both atomic signals have been normalized to better highlight the linear response of the dual atom accelerometer.

Fig. 9

 9 Fig. 9 Extraction of α from the experimental signals of Fig.8. The frequency chirp, corresponding to a constant acceleration applied on the instrument, can be extracted from interferometric signals of 85 Rb alone (red circles), of 87 Rb alone (not represented on the figure but similar to the 85 Rb signal with an offset due to the phase quadrature) or from the hybridization of both signals (blue-grey circles). The true value of the applied frequency chirp is given in green line. The dashed lines show the measurement range associated to the single species 85 Rb interferometer (in red) and the measurement range associated to the dual species 85 Rb/ 87 Rb interferometer in phase quadrature (in blue). For the dual species signal, only data in the linear range of each single atom interferometer have been kept, excluding therefore low sensitive data near fringe extrema.

Table 1

 1 Estimation of the sensitivity gains obtained by coupling the linear ranges of two interferometers in phase quadrature compared to a single species interferometer, as a function of the Signal-to-Noise Ratio (SNR).

	SNR	Sensitivity Gain
	10 2	2.7
	10 3	4.8
	10 4	8.5


Jingjing Pan

Cédric Le Bastard

Member, IEEE Yide Wang

Meng Sun
email: msun@shmtu.edu.cn

Time-Delay Estimation

Keywords: Ground Penetrating Radar (GPR), Support Vector Regression (SVR), Forward-Backward Linear Prediction (FBLP), Time-Delay Estimation (TDE), low snapshots

  L'archive ouverte pluridisciplinaire

G ROUND penetrating radar (GPR) is a common tool for subsurface sensing in the field of civil engineering, defense, agriculture and environment [START_REF] Daniels | Ground penetrating radar[END_REF]- [START_REF] Benedetto | Civil engineering applications of ground penetrating radar[END_REF]. It allows nondestructive probing and therefore gains much interest both in media parameters estimation and in buried targets localization [START_REF] Li | Enhanced GPR signal for layered media time-delay estimation in low-SNR scenario[END_REF]- [START_REF] Giannakis | A realistic FDTD numerical modeling framework of ground penetrating radar for landmine detection[END_REF].

In civil engineering, GPR is used to survey horizontally stratified media, for example, roadways. The information of the vertical structure of the stratified media can be extracted from radar profiles by means of echo detection and amplitude estimation. Echo detection provides the time-delay estimation (TDE) associated with each interface, whereas amplitude estimation is used to retrieve the wave speed within each layer. However, TDE has several signal processing challenges in practical GPR applications. Problem arises when the backscattered echoes are too close to each other, which requires the use of signal processing algorithms with high temporal resolution. Subspace methods and linear prediction (LP) are more appropriate than the conventional FFT-based methods for estimating time-delays of overlapping echoes [START_REF] Bastard | Thin-pavement thickness GPR with high-resolution and superresolution methods[END_REF]. Furthermore, the performance of signal processing algorithms degrades when the number of available snapshots used for covariance matrix estimation is limited [START_REF] Gaudes | Robust array beamforming with sidelobe control using support vector machines[END_REF]- [START_REF] Xin | Linear prediction approach to direction estimation of cyclostationary signals in multipath environment[END_REF]. In addition, the backscattered echoes obtained from GPR are usually highly correlated (even coherent) with each other, which significantly adds up to the difficulties in TDE. Decorrelation techniques [START_REF] Sun | Time-delay estimation using ESPRIT with extended improved spatial smoothing techniques for radar signals[END_REF], [START_REF] Shan | On spatial smoothing for directionof-arrival estimation of coherent signals[END_REF], [START_REF] Qu | Time-delay estimation for ground penetrating radar using ESPRIT with improved spatial smoothing technique[END_REF] are attractive, especially with subspace methods, like multiple signal classification (MUSIC) [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] and estimation of signal parameters via rational invariance technique (ESPRIT) [START_REF] Roy | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF].

It is worth noting that LP methods have a high temporal resolution. They are able to perform spatial smoothing implicitly [START_REF] Shan | On spatial smoothing for directionof-arrival estimation of coherent signals[END_REF], which is of great importance with coherent signals.

In GPR applications, there has been a considerable number of research making use of LP in the detection of buried objects [START_REF] Van Genderen | Imaging of stepped frequency continuous wave GPR data using the Yule-Walker parametric method[END_REF], [START_REF] Ho | A linear prediction land mine detection algorithm for hand held ground penetrating radar[END_REF]- [START_REF] Ratto | Analysis of linear prediction for soil characterization in GPR data for countermine applications[END_REF]. LP is used to predict the next GPR signal from the previous observations. An object appears when the measured signal is different from the prediction. There are also applications about autoregressive model (AR, one-side LP) in the estimation of parameters of subsurface materials, for example, the estimation of time-delay [START_REF] Bastard | Thin-pavement thickness GPR with high-resolution and superresolution methods[END_REF] and soil permittivity [START_REF] Ratto | Estimation of soil permittivity through autoregressive modeling of time-domain ground-penetrating radar data[END_REF]. Nevertheless, the performance of LP is limited when the observation records are short [START_REF] Van Genderen | Imaging of stepped frequency continuous wave GPR data using the Yule-Walker parametric method[END_REF].

In essence, LP methods are about the estimation of a weight vector. A promising solution to the estimation problem due to the short observation records is the support vector regression (SVR) [START_REF] Vapnik | Statistical learning theory. Adaptive and learning systems for signal processing, communications and control[END_REF]. SVR is the regression form of support vector machine (SVM). Based on the principle of structural risk minimization, SVR has fantastic generalization ability. It is originally used only in the real domain and later extended into complex form [START_REF] Gaudes | Robust array beamforming with sidelobe control using support vector machines[END_REF], [START_REF] Gonnouni | A support vector machine MUSIC algorithm[END_REF], [START_REF] Bouboulis | Complex support vector machines for regression and quaternary classification[END_REF]. In the literature, SVR has been integrated with several linear signal processing algorithms. A support vector autoregressive method (AR-SVR) is proposed for frequency estimation problems [START_REF] Xu | Applications of support vector machines in electromagnetic problems[END_REF], [START_REF] Rojo-Álvarez | Support vector machines framework for linear signal processing[END_REF]. The authors in [START_REF] Rojo-Álvarez | Support vector method for robust ARMA system identification[END_REF] combines autoregressive moving average (ARMA) with SVR for system identification applications. The theory of AR is closely related to the forward linear prediction (FLP) [START_REF] Haykin | Adaptive filter theory[END_REF]. However, forward-backward linear prediction (FBLP) performs better than one-side prediction methods (FLP and backward linear prediction, BLP) [START_REF] Chen | Kalman-based spatial domain forwardbackward linear predictor for DOA estimation[END_REF]. Therefore, we propose to combine the theory of SVR with FBLP together in TDE. The objective function can be transformed from complex to real domain like in [START_REF] Gaudes | Robust array beamforming with sidelobe control using support vector machines[END_REF]. In this paper, the function is directly formulated in complex domain with Wirtinger's calculus [START_REF] Bouboulis | Complex support vector machines for regression and quaternary classification[END_REF], [START_REF] Bouboulis | Extension of Wirtinger's calculus to reproducing kernel Hilbert spaces and the complex kernel LMS[END_REF].

A c c e p t e d m a n u s c r i t p t

The remainder of the paper is organized as follows. Section II presents the signal model. Section III introduces the theory of FBLP and the proposed SVR based FBLP method. Section IV presents some simulation results. Section V shows the performance of the proposed method with experimental data. The conclusion is drawn in Section VI.

In the following, (.) T , (.) * , and (.) H denote transpose, conjugate, and complex conjugate transpose operations, respectively. R and C denote real and complex sets, respectively. Re(z) and Im(z) denote the real and imaginary parts of z, respectively. Vectors and matrices appear in boldface lowercase letters and boldface capital letters, respectively.

II. SIGNAL MODEL

The signal model used in [START_REF] Bastard | Thin-pavement thickness GPR with high-resolution and superresolution methods[END_REF], [START_REF] Qu | Time-delay estimation for ground penetrating radar using ESPRIT with improved spatial smoothing technique[END_REF], [START_REF] Li | An efficient algorithm for time delay estimation[END_REF] with stratified media is adopted in this paper. Thus, the backscattered echoes are the time-shifted and attenuated replicas of the transmitted signal. The observation is performed with M discrete frequencies. The frequency f m is defined as f m = f 1 + (m -1)∆f , with f 1 the beginning of the bandwidth and ∆f the frequency shift, m = 1, . . . , M . The received signal at frequency f m can be expressed in frequency domain as

g(f m ) = K k=1 e(f m )s k e -j2πfmτ k + n(f m ) (1)
where K is the number of backscattered echoes, which can be primary and multiple reflection echoes. s k denotes the amplitude of the kth backscattered echo. τ k is the kth timedelay corresponding to the kth echo. At frequency f m , e(f m ) represents the radar pulse in frequency domain while n(f m ) is an additive white Gaussian noise with zero mean and variance σ 2 . In order to use the principle of LP methods, a whitening procedure [START_REF] Bastard | Thin-pavement thickness GPR with high-resolution and superresolution methods[END_REF], [START_REF] Sun | Time-delay estimation using ESPRIT with extended improved spatial smoothing techniques for radar signals[END_REF] is applied. In the following, each echo is divided by the radar pulse. Thus, (1) becomes:

r m = g(f m ) e(f m ) = K k=1 s k e -j2πfmτ k + n(f m ) e(f m ) . (2)
With L independent snapshots, the data sample at frequency f m can be written as

r m = [r m (1), r m (2), . . . , r m (L)] T . (3)
III. METHODOLOGY

A. Forward-backward linear prediction

In theory, LP estimates the unknown samples with a linear combination of the known observations by minimizing the mean square prediction error. FBLP uses the observations from both forward and backward observation sequences. If the order of the prediction filter equals P , the prediction equation can be modeled as [START_REF] Wang | Spatial spectrum estimation theory and algorithm[END_REF], [START_REF] Tufts | Estimation of frequencies of multiple sinusoids: Making linear prediction perform like maximum likelihood[END_REF]: . . . . . . . . .

            r P r P -1 . . . r 1 . . . . . . . . . r M -1 r M -2 . . .
r * M -P +1 r * M -P +2 . . . r * M                ω 1 . . . ω P    =             r P +1 . . . r M . . . . . . r * 1 . . . r * M -P             . ( 4
)
The other LP methods, FLP or BLP, can be expressed with the first or second half (above or below dotted line) of ( 4), respectively. In a more compact way, (4) can be written as

Xω = y (5)
where

y ∈ C N T ×1 , X ∈ C N T ×P and ω ∈ C P ×1 , N T = 2(M -P )L.
The weight coefficient vector ω can be estimated by [START_REF] Wang | Spatial spectrum estimation theory and algorithm[END_REF]:

ω = R -1 r (6)
where R = X H X/L and r = X H y/L. The inversion of R requires the inequality constraint: N T P . Therefore, we have the constraint on the number of snapshots

L, L ≥ P/[2(M -P )].
The order of the prediction filter P should satisfy:

K P M -K/2 [29].
The time-delay of each echo can be estimated by searching the peak positions in the power spectrum density (PSD) of FBLP. With the weight vector ω, the PSD is defined as

P (t) = 1 a H (t) 1 -ω 2 (7)
with a(t) = [e -2jπf1t , e -2jπf2t , . . . , e -2jπf P +1 t ] T , the mode vector.

B. Proposed method: FBLP-SVR

The key issue in LP methods is to estimate the weight vector ω. However, in scenarios where the number of available snapshots is insufficient, the performance of LP methods suffers great deterioration. The principle of SVR works not only with limited number of samples but also with complex data. Therefore, we propose to combine SVR with FBLP in this paper, which is called FBLP-SVR in the following.

We can see that ( 5) is a typical form of SVR in complex domain. In this paper, SVR is directly formulated in complex domain to deal with complex data in a natural way [START_REF] Gonnouni | A support vector machine MUSIC algorithm[END_REF], [START_REF] Bouboulis | Complex support vector machines for regression and quaternary classification[END_REF]. The research in [START_REF] Bouboulis | Complex support vector machines for regression and quaternary classification[END_REF] deals with complex-valued nonlinear regression problems by exploiting the widely linear estimation method in complex domain, which is very instructive in finding the solution of [START_REF] Sun | Time-delay estimation using ESPRIT with extended improved spatial smoothing techniques for radar signals[END_REF]. In regression problems, the main objective of SVR is to find a hyper plane to fit the data within a deviation less than a given value . The -insensitive loss function is used here [START_REF] Martínez-Ramón | Support vector machines for antenna array processing and electromagnetics[END_REF]. The primal objective function is to minimize the prediction error both structurally and empirically. Since the A c c e p t e d m a n u s c r i t p t empirical errors are complex-valued, the residuals should be minimized both in their real and imaginary parts. Therefore, the optimization task can be expressed as [START_REF] Gonnouni | A support vector machine MUSIC algorithm[END_REF], [START_REF] Bouboulis | Complex support vector machines for regression and quaternary classification[END_REF], [START_REF] Martínez-Ramón | Support vector machines for antenna array processing and electromagnetics[END_REF]:

min ω,ξ r , ξr ,ξ i , ξi ( 1 2 ω 2 + C N T n=1 (ξ r n + ξr n + ξ i n + ξi n )) (8) s. t.                Re(y n -x H n ω) + ξ r n , n = 1, . . . , N T Re(-y n + x H n ω) + ξr n , n = 1, . . . , N T Im(y n -x H n ω) + ξ i n , n = 1, . . . , N T Im(-y n + x H n ω) + ξi n , n = 1, . . . , N T ξ r n , ξr n , ξ i n , ξi n 0, n = 1, . . . , N T
where x n is the nth column of X H and y n = y(n), n = 1, . . . , N T . ξ r n and ξr n stand for positive and negative errors in the real part of y n while ξ i n and ξi n are for the corresponding imaginary part. The value C controls the trade-off between the structural and empirical errors.

The method of Lagrangian multipliers is employed to find the solutions of ( 8) by introducing a dual set variables. Therefore, we have the primal-dual objective function:

L pd (ω) = 1 2 ω 2 + C N T n=1 (ξ r n + ξr n + ξ i n + ξi n ) + N T n=1 a n (Re(y n -x H n ω) --ξ r n ) + N T n=1 ân (Re(-y n + x H n ω) --ξr n ) + N T n=1 b n (Im(y n -x H n ω) --ξ i n ) + N T n=1 bn (Im(-y n + x H n ω) --ξi n ) - N T n=1 η n ξ r n - N T n=1 ηn ξr n - N T n=1 λ n ξ i n - N T n=1 λn ξi n ( 9
)
where a n , ân , b n , bn , η n , ηn , λ n , and λn are Lagrangian multipliers.

(8) and ( 9) are real-valued functions defined on complex variables. The constructed Lagrange function has a saddle point which minimizes over the primal variables and maximizes over the dual variables. In order to exploit the Karush-Kuhn-Tucker (KKT) theorem, the Wirtinger's calculus [START_REF] Bouboulis | Complex support vector machines for regression and quaternary classification[END_REF], [START_REF] Bouboulis | Extension of Wirtinger's calculus to reproducing kernel Hilbert spaces and the complex kernel LMS[END_REF] is carried out for the complex variable ω. Besides, the gradients of real variables are computed in the traditional way. Then, we deduce:

                 ∂L pd ∂ω * = 1 2 ω -1 2 N T n=1 ((a n -ân ) + i(b n -bn ))x n = 0 ∂L pd ∂ξ r n = C -a n -η n = 0, n = 1, . . . , N T ∂L pd ∂ ξr n = C -ân -ηn = 0, n = 1, . . . , N T ∂L pd ∂ξ i n = C -b n -λ n = 0, n = 1, . . . , N T ∂L pd ∂ ξi n = C -bn -λn = 0, n = 1, . . . , N T . (10
) Substituting ( 10) into (9), we have the following maximization task:

max a,â,b, b - 1 2 a - â b -b T Re(XX H ) -Im(XX H ) Im(XX H ) Re(XX H ) a - â b -b -1 T (a + â + b + b) + Re(y T )(a -â) + Im(y T )(b -b) s. t. 0 a n , αn , b n , bn C, n = 1, . . . , N T ( 11
)
where 1 is a all-one column vector with

N T elements, a = [a 1 , . . . , a N T ] T , â = [â 1 , . . . , âN T ] T , b = [b 1 , . . . , b N T ] T , and b = [ b1 , . . . , bN T ] T .
In [START_REF] Qu | Time-delay estimation for ground penetrating radar using ESPRIT with improved spatial smoothing technique[END_REF], the quadratic optimization function will reach its maximum with respect to the Lagrange coefficient vectors a, â, b, b, which can be computed using quadratic programming solvers. Then, the weight ω can be obtained according to [START_REF] Shan | On spatial smoothing for directionof-arrival estimation of coherent signals[END_REF], ω =

N T n=1 ((a n -â n )+i(b n -bn ))x n .
In the calculation, a small identity term γI is added in the cost function in case of ill-conditional inaccuracies [START_REF] Martínez-Ramón | Support vector machines for antenna array processing and electromagnetics[END_REF]. Compared with the standard FBLP, the combined FBLP-SVR needs extra computation. However, the number of snapshots in the considered situations being small, the increase of computational burden is insignificant. Qualitative evaluations will be shown in the simulation part.

IV. SIMULATION RESULTS

A. Simulation Settings

In order to evaluate the performance of the proposed FBLP-SVR method, three simulations are carried out in this section. The frequency range of the step frequency radar is 1.0 -4.0 GHz, with M = 21 frequency samples. The underground structure is assumed to have three layers, as shown in Fig. 1. The medium is composed of three layers, Layer 1, Layer 2, and Layer 3. The relative permittivities and thicknesses of the layers are listed in Table I. In the simulation, four echoes (S 0 , S 1 , S M 1 , and S 2 ) are considered. S 0 , S 1 , and S 2 are the primary echoes, with the corresponding time-delays [τ 0 , τ 1 , τ 2 ] = [6.67, 6.95, 7.89] ns. S M 1 is the multiple echo within the first layer with time-delay τ M 1 = 7.24 ns. The time resolution is determined by B∆τ , where B is the bandwidth and ∆τ the time shift between two echoes [START_REF] Bastard | Thin-pavement thickness GPR with high-resolution and superresolution methods[END_REF]. If the product B∆τ is greater than 1, the echoes are distinguishable by conventional FFT based methods. In this simulation, the first three echoes (the first and second primary echo, the multiple echo) are overlapped. The third primary echo is not overlapped with the others. The four echoes are coherent. SNR is defined as the ratio between the power of the last primary echo and the noise variance. The work in [START_REF] Wang | Spatial spectrum estimation theory and algorithm[END_REF] indicates that the smaller the order of prediction filter P is, the worse the estimation performance of LP will be. Therefore, the order P in this section is set to 19. In this situation, there should be at least 5 snapshots when using the standard FBLP so that the covariance matrix is reversible. In contrast, SVR based method does not require such constraint.

In all simulations, the SVR-related parameters are set as = 0, C = 1, and γ = 10 -6 . In our extensive simulation experiences, the value C and γ are insensitive parameters while should be small values. There may be other analytic parameter selection methods or other pairs of parameter settings, but the proposed FBLP-SVR shows its robustness with = 0, C = 1 and γ = 10 -6 at different scenarios, as in [START_REF] Gonnouni | A support vector machine MUSIC algorithm[END_REF].

{ r0 , r1 , r2 , r3 } {1, 3, 8, 9} {H 1 , H 2 } mm {25, 50}

B. Power Spectrum Density

In the first simulation, the normalized PSD of the proposed FBLP-SVR is compared with the standard FBLP and MUSIC. We assume that the top four highest peak locations in the spectrum allow estimating the four time-delays. The simulation is conducted with 5 and 100 snapshots to show the influence of the number of snapshots on the estimation performance. SNR = 20 dB. The results are depicted in Fig. 2. The vertical dashed lines are located at the true values of timedelays. In the scenario with 5 snapshots, the proposed FBLP-SVR can detect the three primary echoes, but the amplitude of PSD of the multiple echo is very weak. The PSD of the standard FBLP has false peaks and it can not correctly detect the multiple echo. When there are 100 snapshots, FBLP and FBLP-SVR perform similarly. Both of them are capable of detecting the three primary echoes and the multiple echo. The performance of FBLP and FBLP-SVR is enhanced with more snapshots. Unfortunately, MUSIC fails to detect the echoes not only with 5 snapshots but also with 100 snapshots, due to the fact that the echoes are totally correlated.

In Fig. 2, the multiple echo has little impact on the three primary echoes, even when the first two primary echoes are overlapped with the multiple echo. Therefore, the estimation of the multiple echo is excluded in the following.

C. Performance versus the number of snapshots

In the second simulation, we evaluate the accuracy of the proposed FBLP-SVR as a function of the snapshots number L, L ∈ [1,. . . ,50]. Since the standard FBLP requires L 5, its results are shown only from 5 snapshots. SVR based FLP (FLP-SVR) is also considered in this part to make a comparison with FBLP-SVR. Other settings are the same as the first simulation. For each number of snapshots, the methods are evaluated with 500 Monte Carlo trials. The performance is assessed with relative-root-mean-square-error (RRMSE), which is defined as [START_REF] Sun | Time-delay estimation using ESPRIT with extended improved spatial smoothing techniques for radar signals[END_REF] RRM SE =

1 J J j=1 (ẑ j -z) 2 z ( 12
)
where ẑj denotes the estimated time-delay for the jth run; z is the true value; J is the number of Monte Carlo trials. The RRMSEs of the three methods versus the number of snapshots are illustrated in Figs. 345. With the increasing of the number of snapshots, the RRMSEs of FBLP, FBLP-SVR FLP-SVR decrease. FBLP-SVR has lower RRMSEs than FLP-SVR at different numbers of snapshots. When the number of snapshots is large (more than 15 in this case), FBLP and FBLP-SVR achieve similar results. But when the number of snapshots is limited, the proposed FBLP-SVR has the best accuracy and outperforms the one-directional FLP-SVR and traditional FBLP.

The execution time, by a computer equipped with a processor unit (CPU) of 2.7 GHz and 16 GB of RAM, is used to get a rough idea about the computational burden of the proposed method. In the comparison, three methods (the traditional FBLP, FLP-SVR, and FBLP-SVR) are evaluated with L = 5 during 500 Monte Carlo trials. The average execution time of one trial using the traditional FBLP, FLP-SVR, and FBLP-SVR are 0.0674 s, 0.0816 s, and 0.0899 s, respectively. In view of the average execution time, these three methods are not time consuming. The combination of FBLP (or FLP) with SVR slightly increases the computational complexity. However, the proposed FBLP-SVR greatly improves the estimation

D. Performance versus SNR

In the third simulation, the methods are applied to GPR data with different SNRs ranging from 0 dB to 30 dB. Only one snapshot is considered in this simulation. The standard FBLP can not work with single snapshot hence it is not presented in the comparison. The performance of FBLP-SVR and FLP-SVR is tested with a Monte Carlo process, which consists of 500 independent runs of the methods. Figs. 678show the RRMSEs of the three primary echoes using FBLP-SVR and FLP-SVR as a function of SNR. It can be seen that the RRMSEs of the three primary echoes continuously decrease as SNR increases. The estimation results of FBLP-SVR have lower RRMSEs than those of FLP-SVR. The RRMSE of the methods depends on the echo amplitude, that is, the larger the echo amplitude is, the smaller RRMSE will be. The RRMSE difference between FBLP-SVR and FLP-SVR is the smallest for the third primary echo (in Figs. 5 and8), which might lie in the fact that the third echo is not overlapped with the other echoes.

V. EXPERIMENT

In this section, the proposed method is tested with two experimental databases.

A. Laboratory experiment

In the first experiment, a PVC slab is probed in laboratory by a monostatic step frequency radar within frequency bandwidth f ∈ [1.6, 3] GHz in far-field, with M = 71 frequency samples. The height of the antenna is 70 cm. The PVC slab has a thickness of 4 cm and relative permittivity r = 2.97 + 0.015j. In this case, the product B∆τ is about 0.64. Thus, the two backscattered echoes overlap with each other. The radar pulse is measured with a metal plane [START_REF] Bastard | Thin-pavement thickness GPR with high-resolution and superresolution methods[END_REF]. The data are acquired with single snapshot. The measured GPR data are preprocessed with temporal filtering and data whitening, like in [START_REF] Bastard | Thin-pavement thickness GPR with high-resolution and superresolution methods[END_REF], [START_REF] Sun | Time-delay estimation using ESPRIT with extended improved spatial smoothing techniques for radar signals[END_REF]. After that, the proposed FBLP-SVR and FBLP are applied to TDE. The order of prediction filter P is set to 70. Fig. 9 illustrates the obtained results. The backscattered echoes are overlapped and correlated with each other and there is only one snapshot. The proposed FBLP-SVR can still accurately detect the two echoes while FBLP fails because of the limitation on the number of snapshots (L P/[2(M -P )] = 35). The relative error of FBLP-SVR in Fig. 9(c) is about 2.7%.

B. Field experiment

Secondly, an experiment is conducted to probe a pavement of IFSTTAR fatigue carousel [START_REF] Simonin | Progress in monitoring the debonding within pavement structures during accelerated pavement testing on the fatigue carousel[END_REF]. The pavement consists of two layers of asphalt. The relative permittivities of the two layers are very close. The thickness of the first layer is about 5 cm. The measurement is implemented by a quasi-monostatic step frequency radar with transmitter (Tx) and receiver (Rx) close to each other. The distance between Tx and Rx antennas is constant during the B-scan. The B-scan is composed of 21 traces (A-scans). During the measurement, the far-field condition is verified. Preprocessing techniques (filtering the air wave, data whitening) are performed before applying the proposed algorithm [START_REF] Bastard | Thin-pavement thickness GPR with high-resolution and superresolution methods[END_REF], [START_REF] Sun | Time-delay estimation using ESPRIT with extended improved spatial smoothing techniques for radar signals[END_REF]. By the inverse Fourier transform, we have the B-scan obtained from the experimental data without air wave with a large frequency band (f ∈ [0.8, 10.8] GHz), shown in Fig. 10. There are two echoes and the time shift ∆τ between them is about 1.07 ns. Fig. 11 shows the B-scan obtained by the proposed FBLP-SVR using only one snapshot. As expected, the two echoes are well resolved by FBLP-SVR over 21 A-scans. The frequency band used in the estimation is f ∈ [3.77, 4.42] GHz with 27 frequency elements. The prediction order P is set to 25. Thus, the product B∆τ is about 0.7, which means that the two echoes are overlapped. The standard FBLP can not work since the limitation on the number of snapshots is L 7 and there is only one snapshot in the experiment. Fig. 12 illustrates the results obtained by the proposed method and FBLP at the 5th trace of the B-scan (5th A-scan of Fig. 10). The proposed method is able to detect the two echoes as in shown Figs.

2 HFig. 1 .

 21 Fig.1. Stratified layers from a simulated setup. H i and ri are the thickness and relative permittivity of Layer i (i ∈ [0, 1, 2, 3]), respectively. S i and τ i are the reflected echo and time-delay from interface i, respectively. S M 1 and τ M 1 are for the multiple echo within the first layer.

Fig. 2 .

 2 Fig. 2. PSD of FBLP-SVR, FBLP and MUSIC in the estimation of timedelays with (a) 5 snapshots and (b) 100 snapshots.

Fig. 3 .Fig. 4 .

 34 Fig. 3. RRMSEs of TDE versus the number of snapshots, first primary echo.

Fig. 5 .

 5 Fig.5. RRMSEs of TDE versus the number of snapshots, third primary echo.

Fig. 6 .Fig. 7 .

 67 Fig. 6. RRMSEs of TDE versus SNR, first primary echo.

Fig. 8 .

 8 Fig. 8. RRMSEs of TDE versus SNR, third primary echo.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Processing of GPR measurements. (a) Raw data. (b) Time-filtered data. (c) PSD of FBLP-SVR and FBLP (B∆τ ≈ 0.64).

  11 and 12. But FBLP shows no peak at the position of true values. In addition, the time-delays between the two echoes for the 21 A-scan are calculated. The mean and median of the estimated time shifts are 1.058 ns and 1.110 ns, respectively. The estimated time-delays are close to the real values.

A c c e p t e d m a n u s c r i t p t A c c e p t e d m a n u s c r i t p t

  r M -P . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	r * 2	r * 3	. . . r * P +1

ACKNOWLEDGMENT

The authors would like to thank the China Scholarship Council (No. 201506040050) for funding part of this work. The authors would like to thank Dr. David Guilbert for the data from the first experiment and Dr. Jean-Michel Simonin, Dr. Vincent Baltazart and Dr. Xavier Dérobert for the data from the second experiment. This work may contribute to COST Action TU1208 Civil Engineering Applications of Ground Penetrating Radar. The authors also wish to thank the referees for their helpful comments.

VI. CONCLUSION

In this paper, a SVR based FBLP is proposed. The proposed method takes advantage of the good properties of FBLP and SVR. FBLP has high resolution and can decorrelate coherent signals directly. SVR is robust with small samples. The performance of the performed method is validated with numerical and experimental data, in coherent scenarios with both overlapping and non-overlapping signals and limited snapshots. The proposed FBLP-SVR outperforms the traditional FBLP and FLP-SVR methods, especially when the number of snapshots is low. Furthermore, FBLP-SVR is applicable with only one snapshot.


Piotr Gwiazda
email: pg-wiazda@mimuw.edu.pl

Benoît Perthame
email: benoit.perthame@sorbonne-universite.fr

Agnieszka Świerczewska-Gwiazda

A two species hyperbolic-parabolic model of tissue growth

Keywords: 2010 Mathematics Subject Classification. 35B45, 35K57, 35K55, 35K65, 35Q92, 76N10, 76S99, Porous medium system, Darcy's law, Cross-diffusion systems, Reactiondiffusion, Hyperbolic-parabolic systems

Models of tissue growth are now well established, in particular in relation to their applications to cancer. They describe the dynamics of cells subject to motion resulting from a pressure gradient generated by the death and birth of cells, itself controlled primarily by pressure through contact inhibition. In the compressible regime we consider, when pressure results from the cell densities and when two different populations of cells are considered, a specific difficulty arises from the hyperbolic character of the equation for each cell density, and to the parabolic aspect of the equation for the total cell density. For that reason, few a priori estimates are available and discontinuities may occur. Therefore the existence of solutions is a difficult problem.

, our strategy is to ignore compactness on the cell densities and to prove strong compactness on the pressure gradient. We improve known results in two directions; we obtain new estimates, we treat higher dimension than 1 and we deal with singularities resulting from vacuum.

Introduction

The topic of modeling tissue growth has recently progressed with various inputs from physics and mechanics [START_REF] Preziosi | Multiphase modelling of tumour growth and extracellular matrix interaction: mathematical tools and applications[END_REF][START_REF] Byrne | A two-phase model of solid tumour growth[END_REF][START_REF] Goriely | Differential growth and instability in elastic shells[END_REF][START_REF] Ranft | Fluidization of tissues by cell division and apoptosis[END_REF]. Models are now used for image-based prediction of cancer growth [START_REF] Bresch | Computational modeling of solid tumor growth: the avascular stage[END_REF][START_REF] Ribba | A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents[END_REF].

They describe the dynamics of cell number density subject to motion resulting from a pressure gradient generated by the death and birth of cells, itself controlled primarily by pressure through contact inhibition. In the compressible regime, pressure results from a combination of the cell densities and controls both the motion through Darcy's law and birth and death of cells according to a finding in [START_REF] Byrne | Individual-based and continuum models of growing cell populations: a comparison[END_REF] and commonly used since then. Models with a single type of cells have been studied recently by many authors, as well as their incompressible limit [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF][START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF][START_REF] Kim | Uniform convergence for the incompressible limit of a tumor growth model[END_REF][START_REF] Kim | Free boundary problems for tumor growth: a viscosity solutions approach[END_REF][START_REF] Maury | Congestion-driven dendritic growth[END_REF]. More general formalisms using incompressibility conditions also occur in two phase flows, and they appear, e.g., in oil recovery [START_REF] Alt | Nonsteady flow of water and oil through inhomogeneous porous media[END_REF][START_REF] Alt | On nonstationary flow through porous media[END_REF] where each phase has its own pressure. Models may also contain several "phases", and have also been widely established and studied [START_REF] Byrne | Modelling solid tumour growth using the theory of mixtures[END_REF][START_REF] Sherratt | A new mathematical model for avascular tumour growth[END_REF][START_REF] Cancès | Incompressible immiscible multiphase flows in porous media: a variational approach[END_REF][START_REF] Kim | Porous medium equation to Hele-Shaw flow with general initial density[END_REF][START_REF] Bertsch | A nonlinear parabolic-hyperbolic system for contact inhibition of cell-growth[END_REF]. For instance, a specific question is to understand when segregation occurs [START_REF] Bertsch | On interacting populations that disperse to avoid crowding: the case of equal dispersal velocities[END_REF][START_REF] Carrillo | Splitting schemes & segregation in reaction-(cross-)diffusion systems[END_REF].

Here, we consider the following compressible two cell population model, that we state in the full space for the sake of simplicity,

   ∂ t n 1 -div[n 1 ∇p] = n 1 F 1 (p) + n 2 G 1 (p), x ∈ R d , t ≥ 0, ∂ t n 2 -div[n 2 ∇p] = n 1 F 2 (p) + n 2 G 2 (p), (1)
with

n := n 1 + n 2 , p = n γ , γ > 1. (2)
We assume that there is a value P H > 0 (the name homeostatic pressure was coined in [START_REF] Ranft | Fluidization of tissues by cell division and apoptosis[END_REF]) such that the smooth functions F i , G i , describing the division/death rates of cells, satisfy the properties

F (p) := F 1 (p) + F 2 (p) ≤ 0, G(p) := G 1 (p) + G 2 (p) ≤ 0, ∀p ≥ P H . (3)
We also assume that the initial data n 0 1 , n 0 2 , n 0 = n 0 1 + n 0 2 satisfy

n 0 1 ≥ 0, n 0 2 ≥ 0, p 0 := (n 0 1 + n 0 2 ) γ ≤ P H , (4)
n 0 (1 + |x| 2 + | ln(n 0 )|) ∈ L 1 (R d ), (5)
∇p 0 ∈ L 2 (R d ), ∆p 0 ∈ M loc (R d ), ∆p 0 -∈ L 2 loc (R d ), (6)
where M loc (R d ) refers to the vector space of locally bounded measures. At some point, we will also need the restrictions that γ is large enough when d ≥ 5 and that near p = 0 some cancelation occurs, namely

γ > 2 - 4 d , sup 0≤p≤P H |F (p) -G(p)| 2 p 1/γ ≤ C H . (7)
In words, the total proliferation rates of cells n 1 and n 2 are the same when p ≈ 0.

A specific difficulty arises from the hyperbolic character of the equation for each cell density n i , and to the parabolic aspect of the total cell density n. For example, it is known that solutions n 1 , n 2 may have discontinuities. For that reason, the existence of solutions is a difficult problem by lack of strong a priori estimates. Also we cannot hope for strong solutions in general. Here, we establish the existence of weak solutions. In opposition to the method used in the recent paper [START_REF] Carrillo | Splitting schemes & segregation in reaction-(cross-)diffusion systems[END_REF], our strategy is to ignore compactness on the cell densities and to prove strong compactness on the pressure gradient. Therefore, we improve known results in two directions; we treat higher dimension than 1 as in [START_REF] Carrillo | Splitting schemes & segregation in reaction-(cross-)diffusion systems[END_REF] and we deal with vacuum while [START_REF] Bertsch | A nonlinear parabolic-hyperbolic system for contact inhibition of cell-growth[END_REF] only considers uniformly positive and smooth solutions.

Theorem 1 (A priori estimates) With the assumptions (3)-( 6), the following estimates hold true for all T > 0 with constants C(T ) which only depend on the bounds in the above assumptions

n(x, t) ≥ 0, R d n(x, t) dx ≤ Ce Ct , p(x, t) ≤ P H , (8)
T 0 R d |∇p| 2 p 1-1/γ dxdt ≤ C(T ). ( 9
)
Assuming also [START_REF] Bouchut | Uniqueness, renormalization, and smooth approximations for linear transport equations[END_REF], we have for t ∈ (0, T ), and with φ(

•) ∈ C 2 Comp (R d ) a localizing function, R d |∆p(t)| 2 -φ(x) dx ≤ C(T ), T 0 R d |∆p(t)| 3 -φ(x) dxdt ≤ C(T ), (10)
R d |∆p(x, t)| φ(x) dx ≤ C(T ). ( 11
)
Since our framework includes the Barenblatt solutions, see [START_REF] Vázquez | The porous medium equation. Mathematical theory. Oxford Mathematical Monographs[END_REF], we know that these estimates are sharp in the sense that ∆p may be a singular measure supported by the free boundary. Note that p being bounded, the estimate (9) also gives an L 2 t,x bound on ∇p. Another a priori estimate is also available, which we do not use in the subsequent results, and that we postpone to the Appendix.

As a consequence of the estimates in Theorem 1, we establish the following stability result Theorem 2 (Stability of weak solutions) Assume (3) and [START_REF] Bouchut | Uniqueness, renormalization, and smooth approximations for linear transport equations[END_REF] and that the family of initial data satisfies, with uniform bounds, the assumptions (4)- [START_REF] Bertsch | A nonlinear parabolic-hyperbolic system for contact inhibition of cell-growth[END_REF]. Then, the corresponding weak solutions n ε i , with the above bounds true, satisfy after extraction of subsequences,

n ε i n i , in L ∞ (0, T ) × R d -w * , i = 1, 2, n ε → n, p ε → p, in L q (0, T ) × R d , 1 ≤ q < ∞, ∇p ε → ∇p, in L 2 (0, T ) × R d
and n 1 , n 2 , p satisfy, in the weak sense, the system (1)-( 2) with initial data n 0 1 , n 0 2 .

Finally, these two results lead us to the existence theorem, which is the main result of the current paper.

Theorem 3 (Existence of weak solutions) With the assumptions of Theorem 2, there exists a weak solution n 1 , n 2 , p ∈ L ∞ ((0, T ) × R d ) to the system (1)-(2), i.e., for i = 1, 2

T 0 R d -n i ∂ t ψ + n i ∇p.∇ψ -n 1 F i (p) + n 2 G i (p) ψ dxdt = R d n 0 i ψ(0) dx (12)
holds for all

ψ ∈ C 1 Comp (R × R d ) and relations (2) hold a.e. in (0, T ) × R d .
The main observation is that, while our problem is of hyperbolic nature, we can take advantage of informations coming from the parabolic equation on n

∂ t n -div[n∇p] = n 1 F (p) + n 2 G(p) =: nR(c 1 , c 2 , p), x ∈ R d , t ≥ 0, (13)
where, following [START_REF] Carrillo | Splitting schemes & segregation in reaction-(cross-)diffusion systems[END_REF], we define

c i = n i n ≤ 1 and c i (x, t) = 0 when n(x, t) = 0, R = c 1 F (p) + c 2 G(p) ∈ L ∞ . (14)
Next, multiplying equation [START_REF] Byrne | A two-phase model of solid tumour growth[END_REF] with p (n), we compute that p satisfies

∂ t p -|∇p| 2 -γp∆p = γpR. (15)
It is also useful for later purpose to state the equation for the c i 's

∂ t c i -∇p.∇c i = c 1 F i (p) + c 2 G i (p) -c i R. ( 16
)
To obtain the equation for c i we multiply the equation for n i with 1/n and add it to the equation for n multiplied by -n i /n 2 . Indeed, observe that

- 1 n div[n i ∇p] + n i n 2 div[n∇p] = - 1 n ∇n i .∇p - n i n ∆p + n i n 2 ∇n.∇p + n i n 2 n∆p = -∇( 1 n n i ).∇p
and the remaining terms are immediate.

The rest of the paper is devoted to the proofs of these three theorems which we perform in the three next sections. Some remarks and open problems are commented in the conclusion.

Proof of Theorem 1

The first estimates come from the balance law expressed by the equation (1) and from the maximum principle for equation [START_REF] Cancès | Incompressible immiscible multiphase flows in porous media: a variational approach[END_REF]. One easily gets, by integrating (13) over R d and using the Gronwall inequality, that

R d n(t, x) dx ≤ R d n 0 (x) dx exp(t R ∞ ). ( 17
)
To show the uniform bound on p we multiply [START_REF] Cancès | Incompressible immiscible multiphase flows in porous media: a variational approach[END_REF]

) 2 + + [γη (p) -(p -P H ) + ]|∇p| 2 -γ∆η(p) = γp(p -P H ) + R(p). (18)
We integrate over R d and observe that as γ > 1, thus γη (p) -(p -P H ) + ≥ 0 and

1 2 d dt R d (p -P H ) 2 + dx ≤ γ R d p(p -P H ) + R(p) ≤ 0, (19)
where the last inequality follows from ( 3) and ( 4).

Before passing to the next estimate, on |∇p| 2 , let us observe that the second moment of n is bounded. Indeed, multiplying (13) with x 2 Φ L , where Φ L (|x|) is a radially symmetric smooth function, which vanishes outside the ball of radius L + 1, equals to 1 on the ball of radius L, with ∇Φ L and ∆Φ L bounded uniformly in L; and integrating by parts over R d gives

d dt R d x 2 nΦ L dx + γ γ + 1 R d n γ+1 (x 2 ∆Φ L + 4x∇Φ L + 2Φ L ) dx = R d x 2 nRΦ L dx. (20)
Since n γ ≤ P H , we can furthermore obtain

d dt R d x 2 nΦ L dx ≤ C R d x 2 nΦ L dx+ R d nΦ L dx+ {L≤|x|≤L+1} nx|∇Φ L | dx+ {L≤|x|≤L+1} nx 2 |∆Φ L | dx, (21)
where the constant C depends on R ∞ , γ and P H . As n ∈ L 1 , we can claim that the second term on the right-hand side is bounded. For the moment let us assume that n vanishes sufficiently fast at infinity, what we will prove later. Then, by the Lebesgue dominated convergence theorem we can pass to the limit in terms containing Φ L and using that n vanishes for large x and both ∇Φ L and ∆Φ L are uniformly bounded, we show that the last two terms on the right-hand side vanish. We complete the estimate by applying the Gronwall inequality.

The estimate (9) comes from the entropy relation. We multiply [START_REF] Byrne | A two-phase model of solid tumour growth[END_REF] with Φ L ln n, where Φ L is the same truncation function as above, integrate over R d , and find

d dt R d n(ln(n) -1)Φ L dx + 1 γ R d p -1+1/γ |∇p| 2 Φ L dx - 1 γ -1 {L≤|x|≤L+1} p 1-1 γ (ln(p) -1)∆Φ L dx = R d n ln(n) RΦ L dx. (22)
It is easy to observe that if the function n, and thus also p vanishes sufficiently fast, then again the integral over the annulus vanishes as L → ∞ and we conclude [START_REF] Bresch | Computational modeling of solid tumor growth: the avascular stage[END_REF]. For that purpose we recall also that a control of the second moment in x is used here to control the negative values of n ln(n). Indeed, observe that

R d n| ln(n)| dx = R d n ln(n) dx - {x∈R d :n<1} n ln(n) dx ≤ R d n ln(n) dx + 2 R d n|x| 2 dx + c. ( 23
)
And the above inequality allows us to get

R d n| ln(n)| dx ≤ R ∞ T 0 R d n| ln(n)| + 2 R d n|x| 2 + R d n dx + R d n 0 (ln(n 0 ) -n 0 ) dx + c. ( 24
)
We complete the estimate (9) using the Gronwall lemma.

Finally, the fundamental estimates (10) come from Aronson and Benilan's method [START_REF] Aronson | Régularité des solutions de l'équation des milieux poreux dans R N[END_REF][START_REF] Vázquez | The porous medium equation. Mathematical theory. Oxford Mathematical Monographs[END_REF] for the porous media equation with several adaptations. Firstly, and this s a new feature here, we weaken their estimate to L 2 rather than L ∞ . Secondly, we need to localize the estimate in space. Thirdly, we adapt the functional under consideration using also the idea from [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF], and we do not work directly with ∆p but with

w = ∆p + R, ∂ t p = |∇p| 2 + γpw.
We compute

∂ t ∆p = 2(∂ ij p) 2 + 2∇p∇∆p + γ∆(pw), ∂ t R = R c 1 ∂ t c 1 + R c 2 ∂ t c 2 + R p ∂ t p = F (p)[∇c 1 .∇p + ξ 1 ] + G(p)[∇c 2 .∇p + ξ 2 ] + R p [|∇p| 2 + γpw] (25)
where

ξ i := c 1 F i (p) + c 2 G i (p) + c i R, i = 1, 2,
are the right-hand sides from equations [START_REF] Carrillo | Splitting schemes & segregation in reaction-(cross-)diffusion systems[END_REF]. Therefore, since c 1 + c 2 = 1, we find

∂ t w = 2(∂ ij p) 2 + 2∇p∇∆p + γ∆(pw) + [F (p) -G(p)]∇c 1 .∇p + R p [|∇p| 2 + γpw] + Bdd 1
with "Bdd" terms which are bounded in L ∞ ,

Bdd 1 := F (p)ξ 1 + G(p)ξ 2 ,
but this may change from line to line. Since

∇p.∇∆p = ∇p.∇(w -R) = ∇p.∇w -div(R∇p) + R(w -R)
this is also

∂ t w ≥ 2 d (∆p) 2 + 2∇p∇w -2div(R∇p) + 2R(w -R) + γ∆(pw) + [F (p) -G(p)]∇c 1 .∇p + R p [|∇p| 2 + γpw] + Bdd 1 . (26)
The negative part, that we denote by |w| -, therefore satisfies, with sgn -:= 1 I {w<0} ,

∂ t |w| -≤ - 2 d |w| 2 -+ 2∇p∇|w| -+ 2sgn -div(R∇p) + 2(1 - 2 d )R|w| -+ γ∆(p|w| -) -sgn -[F (p) -G(p)]∇c 1 .∇p -sgn -R p |∇p| 2 + γR p p|w| -+ Bdd (27) using that -2 d (w -R) 2 = -2 d w 2 + 4 d Rw -2 d R 2
, where the last term we include within bounded terms that we still gather in Bdd.

We reorganize this inequality as (here the parameter α > 0 can be chosen as small as we wish)

∂ t |w| -≤ -( 2 d -α)|w| 2 -+ 2∇p∇|w| -+ γ∆(p|w| -) + 2sgn -div(R∇p) -sgn -[F (p) -G(p)]∇c 1 .∇p -sgn -R p |∇p| 2 + Bdd .
Notice that, above, we have applied the Young inequality to the terms 2(1

-2 d )R|w| -and γR p p|w| -, that is 2R|w| -(1 - 2 d ) ≤ α 2 |w| 2 -+ c(α)R 2 , R p γp|w| -≤ α 2 |w| 2 -+ c(α)|R p γp| 2 ,
and thus the term Bdd is given by

Bdd := F (p)ξ 1 + G(p)ξ 2 - 2 d R 2 + c(α)(R 2 + |γR p p| 2 ).
We need to localize and use a nonnegative, compactly supported, smooth test function Φ to compute

d dt R d |w| 2 - 2 Φ dx ≤ I + II. ( 28
)
The terms II are those with ∇Φ which are certainly better because they contain one less derivative in the unknowns.

For the difficult term we have, after several integrations by parts, in particular to eliminate derivatives in c 1 which are the worse,

I ≤ -( 2 d -α) R d Φ|w| 3 -dx - R d Φ|w| 2 -∆p dx -γ R d Φp|∇|w| -| 2 dx + γ 2 R d Φ|w| 2 -∆p dx -2 R d ΦR∇|w| -.∇p dx + R d Φc 1 [F (p) -G(p)]|w| -∆p dx + R d Φc 1 [F (p) -G(p)]∇|w| -.∇p dx + C R d Φ|w| -|∇p| 2 dx + R d Φ|w| -Bdd dx . ( 29
)
where C is constant which here takes into account c 1 (F -G ) as well as R p and which is changing from line to line below. The linear and quadratic terms in |w| -are not a problem because the dominant term contains a cubic power of |w| -. We observe about the third-power terms that

R d Φ|w| 2 -∆p dx = R d Φ|w| 2 -(w -R) dx = - R d Φ|w| 3 -dx - R d Φ|w| 2 -R dx. Because R = G(p) + c 1 [F (p) -G(p)],
we also have

I ≤ -( γ 2 + 2 d -1 -α) R d Φ|w| 3 -dx -γ R d Φp|∇|w| -| 2 dx -2 R d ΦG∇|w| -∇p dx + R d Φc 1 [F (p) -G(p)]|w| -∆p dx - R d Φc 1 [F (p) -G(p)]∇|w| -.∇p dx + C R d Φ|w| -|∇p| 2 dx + (1 - γ 2 ) R d Φ|w| 2 -R dx + R d Φ|w| -Bdd dx , (30)
it holds

R d Φc 1 [F (p) -G(p)]|w| -∆p dx = - R d Φc 1 [F (p) -G(p)]|w| 2 -dx - R d Φc 1 [F (p) -G(p)]|w| -R dx
and estimating further (the second term of the above we include already in

R d Φ|w| -Bdd dx), I ≤ -( γ 2 + 2 d -1 -α) R d Φ|w| 3 -dx -γ R d Φp|∇|w| -| 2 dx + 2 R d ΦG|w| -∆p dx + 2 R d ΦG (p)|w| -|∇p| 2 dx + 1 2 R d Φpc 2 1 |∇|w| -| 2 dx + 1 2 R d Φ [F (p) -G(p)] 2 p |∇p| 2 dx + C R d Φ|w| -|∇p| 2 dx + (1 - γ 2 ) R d Φ|w| 2 -R dx - R d Φc 1 [F (p) -G(p)]|w| 2 -dx + R d Φ|w| -Bdd dx . (31)
We arrive at the final form, using the constant C H in [START_REF] Bouchut | Uniqueness, renormalization, and smooth approximations for linear transport equations[END_REF],

I ≤ -( γ 2 + 2 d -1 -α) R d Φ|w| 3 -dx -(γ - 1 2
)

R d Φp|∇|w| -| 2 dx + C R d Φ|w| 2 -dx + C H 2 R d Φ |∇p| 2 p 1-1/γ dx + C R d Φ|w| -|∇p| 2 dx + R d Φ|w| -Bdd dx . ( 32
)
where a constant C standing next to the integral R d Φ|w| -|∇p| 2 dx takes into account also 2G . The first two terms on the right-hand side are the "good terms". Indeed, they have a good sign, unless γ > 2(1 -2 d ), which is automatically satisfied as we assume γ > 1 in d ≤ 4 and in higher dimensions implies higher requirement of the exponent γ as stated in [START_REF] Bouchut | Uniqueness, renormalization, and smooth approximations for linear transport equations[END_REF]. The difficult term is at the highest order

R d Φ|w| -|∇p| 2 dx = - R d Φ[∇|w| -p∇p + |w| -p∆p] dx. ( 33
)
which is under control. Indeed,

R d Φ∇|w| -p∇p dx ≤ (γ - 1 2 )
R d Φp|∇|w| -| 2 dx + c(γ) R d Φp|∇p| 2 dx.
Here, the first term on the right hand side just cancels the second "good term"and the second is bounded.

For the second term of the right-hand side of (33) we have

R d Φ|w| -p∆p dx = R d Φ|w| -p(w -R) dx = - R d Φ|w| 2 -p dx + R d Φ|w| -pR dx
where both these terms are under control to give the final estimate

I ≤ -( γ 2 + 2 d -1 -α) R d Φ|w| 3 -dx + C. (34)
The terms containing gradient of Φ are collected in

II II = R d ∇p|w| 2 -∇Φ dx -γ R d ∇p|w| 2 -∇Φ dx -γ R d p∇(|w| 2 -)∇Φ dx -2 R d R∇p|w| -∇Φ dx + R d [F (p) -G(p)]c 1 ∇p|w| -∇Φ dx + 2 R d ∇ΦG|w| -∇p dx - R d ∇Φ|w| -p∇p dx (35)
and they all do not bring additional difficulties. Therefore, using the inequality (28) and the negative sign in the right hand side of (34), we obtain the a priori estimates announced in [START_REF] Bubba | Incompressible limit for a two species model of tissue growth in one space dimension[END_REF]. The L 1 bound for ∆p in [START_REF] Byrne | Growth of necrotic tumors in the presence and absence of inhibitors[END_REF] is a simple consequence because Φ∆pdx is bounded, therefore Φ|∆p| + dx is controlled by Φ|∆p| -dx which itself is controlled thanks to [START_REF] Bubba | Incompressible limit for a two species model of tissue growth in one space dimension[END_REF].

Proof of Theorem 2

The goal here is to explain the main compactness argument which is used to pass to the limit in an approximate sequence. As in [START_REF] Carrillo | Splitting schemes & segregation in reaction-(cross-)diffusion systems[END_REF], the compactness in time is a major issue.

Weak convergence of the quantities n ε i follows from the bound in L ∞ . The strong convergence of p ε follows from compactness by Sobolev injections. Indeed, on the one hand, we control 9) because the pressure is bounded by P H . On the other hand, we may win time compactness by the Lions-Aubin Lemma using equation [START_REF] Cancès | Incompressible immiscible multiphase flows in porous media: a variational approach[END_REF], which also reads

T 0 R d |∇p ε | 2 dxdt from (
∂ t p ε = (γ -1)|∇p ε | 2 + γ 2 ∆p ε2 + γp ε R ε , ( 36
)
and the space compactness on p ε together with the known bounds provide time compactness. Therefore the expression n ε = (p ε ) 1/γ shows that we may also extract a sub-sequence of n ε which converges.

The strong compactness for ∇p ε is more involved. It mainly relies on the second estimate ( 10) which provides the space compactness, still by the Sobolev embedding theorems. Indeed, the control in L 1 loc of ∆p ε is enough for compactness of ∇p ε , a fact which can be inferred from the representation formula for the solution of the Laplace equation.

For time compactness of ∇p ε , we write, using (36),

∂ t ∇p ε = ∇ (γ -1)|∇p ε | 2 + γp ε R ε + γ 2 ∇∆p ε2 .
Again, we know the local space compactness of ∇p ε from the previous paragraph, the right hand side is a sum of space derivatives of bounded functions, therefore we may apply the Lions-Aubin compactness argument and find that ∇p ε is compact in space and time.

To pass to the limit in the equations is now easy. All the nonlinear terms, that are

n ε i ∇p ε , n ε i F j (p ε ), n ε i G j (p ε ),
have limits as products of weak limits of n ε i by strong limits of p ε and ∇p ε . This completes the proof of Theorem 2.

At this stage, let us point out that our strategy differs deeply from that in [START_REF] Carrillo | Splitting schemes & segregation in reaction-(cross-)diffusion systems[END_REF] based on BV estimates for the quantities c ε i in one dimension. This estimate is somehow sharp since examples with discontinuities on the n ε i are known. Also the method for time compactness is very different since [START_REF] Carrillo | Splitting schemes & segregation in reaction-(cross-)diffusion systems[END_REF] use a control of the Wasserstein distance.

Proof of Theorem 3

We already have a priori estimates and a weak sequential stability result, thus to complete the existence proof we need to construct an approximate system compatible with these estimates. We do that in two steps. Firstly, we make positive the initial data and prove a control from below by a (small) Gaussian. Secondly, we introduce a uniform parabolic regularization.

First step. A regularized problem with a positive control from below. We show that the function

n(t, x) = c exp - |x| 2 2 -ct (37)
is a subsolution to equation ( 13) if we choose c sufficiently large. Since ∂ t n = -cn, ∇n = -xn and ∆n = dn + |x| 2 n, we may insert (37) into the equation for n and as we search for a subsolution, we change equality to inequality. We obtain

-cn -γ(γ + 1)n γ+1 |x| 2 + γn γ+1 -Rn ≤ 0 (38)
which holds true choosing c large enough so that the inequality is satisfied

-c -γ(γ + 1) exp - |x| 2 2 -ct γ |x| 2 + γ exp - |x| 2 2 -ct γ + R ∞ ≤ 0. ( 39
)
It is now a matter of standard estimates, [START_REF] Vázquez | The porous medium equation. Mathematical theory. Oxford Mathematical Monographs[END_REF], to obtain that if we start with the specific initial condition larger than n 0 , we will call it n 0 δ := n 0 + δ exp( -|x| 2 2 ), with δ > 0, then the solution to the problem will be larger than the subsolution n given by (37) with c large enough.

Thus our first approximation step is to replace an initial data n 0 by n 0 δ as it is introduced above. In a consequence the corresponding solution n δ , as well as p δ are locally bounded away from zero, we call these bounds n δ and p δ . The solution has a regularity L q (0, T ; W 2,q loc (R d )), see [START_REF] Vázquez | The porous medium equation. Mathematical theory. Oxford Mathematical Monographs[END_REF] for details and the method in [START_REF] Bertsch | A nonlinear parabolic-hyperbolic system for contact inhibition of cell-growth[END_REF].

Note that an analogue maximum estimate can be proven to provide a bound from above and justify that n vanishes at infinity, what we announced earlier.

Second step. A uniformly parabolic approximation. We consider the system of equations, which consists of a parabolic equation for n and hyperbolic equations for c 1 and c 2 . Thus we construct a parabolic approximation of the equation for c i , i = 1, 2. Let ε > 0,

∂ t c ε i -∇p ε .∇c ε i -ε div[p ε ∇c ε i ] = c ε 1 F i (p ε ) + c ε 2 G i (p ε ) -c ε i R(p ε ). ( 40
)
Note that all the quantities are for simplicity labelled only with ε, but they depend both on ε and δ, i.e. p ε := p ε,δ as well as the other quantities. We proceed now as follows: We solve a parabolic system consisting of ( 13) and (40) with initial data p 0 δ and

n 0 i,δ n 0 δ , completed with the relation p ε = (n ε ) γ .
The equations (40) allow to observe the crucial property, which c i possessed and which was the only information on these quantities used in a priori estimates. Indeed, adding the equations on c ε i , we keep the fundamental relationship c ε 1 + c ε 2 ≡ 1 thanks to the definition of R in ( 14), since initially

c ε 1 + c ε 2 = 1.
Finally, with the fully parabolic framework at hand, it is in the folklore of the domain to obtain the existence of the coupled problem between n ε and c ε i . Next, we notice that all the a priori bounds used to pass to the limit are true. Multiplying with

c ε i |w ε | -Φ and integrating over (0, T ) × R d gives R d |c ε i | 2 |w ε | -Φ dx + T 0 R d ∇p ε .∇c ε i c ε i |w ε | -Φ dxdt + ε T 0 R d p ε |∇c ε i | 2 |w ε | -Φ dxdt + T 0 R d p ε ∇c ε i ∇|w ε | -Φ dxdt = T 0 R d [c ε 1 F i (p ε ) + c ε 2 G i (p ε ) + c ε i R(p ε )]|w ε | -c ε i Φ dxdt + R d |c ε i (0)| 2 |w ε | -Φ dx. ( 41
)
We observe that

T 0 R d ∇p ε .∇c ε i c ε i |w ε | -Φ dxdt = 1 2 T 0 R d ∇p ε .∇|c ε i | 2 |w ε | -Φ dxdt = - 1 2 T 0 R d ∆p ε |c ε i | 2 Φ dxdt - 1 2 T 0 R d ∇p ε |c ε i | 2 ∇Φ dxdt - 1 2 T 0 R d ∇p ε |c ε i | 2 ∇|w ε | -Φ dxdt. (42)
Consequently, we find that

ε T 0 R d p ε |∇c ε i | 2 |w ε | -Φ dxdt + T 0 R d p ε ∇c ε i ∇|w ε | -Φ dxdt ≤ C T 0 R d |∆p ε |Φ dxdt + C T 0 R d |∇p ε |∇Φ dxdt + C T 0 R d ∇p ε |c ε i | 2 ∇|w ε | -Φ dxdt + C T 0 R d |w ε | -c ε i Φ dxdt + R d |c ε i (0)| 2 |w ε | -Φ dx.
(43) Therefore the first integral on the right-hand side can be again estimated by

T 0 R d |∆p ε |Φ dxdt ≤ T 0 R d (|w ε | -+ R)Φ dxdt.
This approximation step will affect our a priori estimates on the level of using computations [START_REF] Maury | Congestion-driven dendritic growth[END_REF], as additional terms related with parabolic approximation, which are estimated above, will appear. Taking these into account, we may pass to the limit as in Section 2, first with ε → 0 and no major difficulty arises. Thus we obtain a limit system for n and c i 's, but still we lack the information whether the equations for n i are satisfied in distributional sense. To recover this we multiply the equations for c ε i with n ε and add the equation for n ε multiplied with c ε i . Let us then define n ε i := c ε i n ε and observe that this operation will lead us to equations for

n ε i ∂ t n ε i -div[n ε i ∇p e ] -εdiv[p ε ∇c ε i ]n ε = n ε 1 F i (p ε ) + n ε 2 G i (p ε ) (44)
The only term that needs to be discussed is εdiv[p ε ∇c ε i ]n ε . To show that this term vanishes in a limit observe that

R d div[p ε ∇c ε i ]n ε φdx = - R d p ε ∇c ε i .∇[(p ε ) 1 γ ]φdx - R d p ε ∇c ε i (p ε ) 1 γ ∇φdx = - 1 γ R d (p ε ) 1 γ ∇p ε .∇c ε i φdx - R d (p ε ) 1 γ +1 ∇c ε i .∇φdx = 1 γ R d (p ε ) 1 γ ∆p ε c ε i φdx + 1 γ R d ∇((p ε ) 1 γ ).∇p ε c ε i φdx + 1 γ R d (p ε ) 1 γ ∇p ε c ε i ∇φdx + R d (p ε ) 1 γ +1 c ε i ∆φdx + R d ∇[(p ε ) 1 γ +1 ]c ε i ∇φdx. ( 45
) Since p ε and c ε i are bounded, then the first term on the right-hand side is bounded due to [START_REF] Byrne | Growth of necrotic tumors in the presence and absence of inhibitors[END_REF]. The boundedness of the second term is provided by [START_REF] Bresch | Computational modeling of solid tumor growth: the avascular stage[END_REF]. For the third and fifth term we use Young's inequality and argue with boundedness of ∇p ε in L 2 t,x . The fourth term is obvious. Thus after letting ε → 0 this error term will vanish. Finally we let δ → 0 and complete the proof.

Conclusion and perspectives

We have proposed a strategy to prove existence of weak solutions for a two species model of tumor invasion. It relies on the extension of the Aronson-Benilan regularizing effect for porous media equations which provides estimates of the Laplacian of the pressure. The most important limitation so far is a combined condition on the two bulk growth terms and it is an open question to remove it. A route in this direction could be to use the energy type estimate given in Theroem 4 in the appendix.

A question which we do not handle here is the strong compactness on the n ε i in the stability result of the approximation process. The bounds on ∆p are too weak for the L 1 theory in [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] and are boarder line to apply the compactness theorems in [START_REF] Ambrosio | Transport equation and Cauchy problem for BV vector fields[END_REF][START_REF] Bouchut | Uniqueness, renormalization, and smooth approximations for linear transport equations[END_REF] which require that D 2 p is a bounded measure.

The extension to more than two species, with the present strategy, requires combined conditions on the three growth terms which read, in the case of three species for instance, c 1 F (p)+c 2 G(p)+c 3 H(p) ≤ Cp 1/γ whenever the nonnegative c i satisfy c 1 + c 2 + c 3 = 1. Then, the analysis goes through without major changes.

There are other questions which arise in this area and that we leave open. One of them concerns the 'incompressible limit' γ → ∞ which has attracted much attention recently [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF][START_REF] Kim | Uniform convergence for the incompressible limit of a tumor growth model[END_REF][START_REF] Craig | Congested aggregation via Newtonian interaction[END_REF][START_REF] Kim | Porous medium equation to Hele-Shaw flow with general initial density[END_REF][START_REF] Kim | Free boundary problems for tumor growth: a viscosity solutions approach[END_REF] because of its relation to congested traffic [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF][START_REF] Brasco | Congested traffic dynamics, weak flows and very degenerate elliptic equations [corrected version of mr2584740[END_REF][START_REF] Maury | Congestion-driven dendritic growth[END_REF][START_REF] Maury | Handling congestion in crowd motion modeling[END_REF]. Clearly the bounds provided here are not enough to investigate this question. However the one dimensional case is under investigation [START_REF] Bubba | Incompressible limit for a two species model of tissue growth in one space dimension[END_REF] based upon arguments from [START_REF] Carrillo | Splitting schemes & segregation in reaction-(cross-)diffusion systems[END_REF]. Another question is about different mobilities, see [START_REF] Lorenzi | On interfaces between cell populations with different mobilities[END_REF][START_REF] Byrne | Individual-based and continuum models of growing cell populations: a comparison[END_REF][START_REF] Byrne | Growth of necrotic tumors in the presence and absence of inhibitors[END_REF], where the parabolic aspects of the equation for n = n 1 + n 2 do not apply.

A Additional a priori bounds

Another remarkable estimate can be obtained for solutions of the system (1)-( 2). We give it here for the sake of completeness. It can be interpreted as some kind of energy because the kinetic energy is given by E K = n |v| 2 2 = p 1/γ |∇p| 2 2 .

Theorem 4 (Energy type a priori estimates) With the assumptions (3)-( 6), the following estimates hold true with constants C(T ) which only depend on the bounds in the above assumptions. For α * = 2 γ , we control 2p . Therefore, we obtain the inequalities (46) and (47).

2 ∇p) -p α+1 2 |∇p| 2 2p ] 2 + C div(p α+1 2 ∇p) -p α+1 2 |∇p| 2

 2222222 To create a negative square, we use the special value of α given byα * = 2 γand the right-hand side is controlled as -[div(p α+1

  with (p -P H ) + . Observe that for any η ∈ C

		2 it
	holds ∆η(p) = ∆pη (p) + η (p)|∇p| 2 , which allows us to handle the highest order term with η (p) =
	p(p -P H ) + . Thus we get
	1 2	∂

t (p -P H

  p α |∇p| 2 2 = p α ∇p.∇[|∇p| 2 + γp∆p + γpR] + αp α-1 |∇p| 2 2 [|∇p| 2 + γp∆p + γpR]. [-p α ∆p -αp α-1 |∇p| 2 + αp α-1 |∇p| 2 2 ] [|∇p| 2 + γp∆p + γpR] [-p α ∆p -αp α-1 |∇p| 2 2 ] [|∇p| 2 + γp∆p + γpR]but p α ∆p is not a good quantity. So the right-hand side has to be rewritten (divide it by γ)

	Therefore, we find
		d dt R d	p α |∇p| 2 2	=		R d
				d dt R d	p α |∇p| 2 2	=		R d
				-[p	α+1 2 ∆p + αp	α+1 2	|∇p| 2 2p	] [p	α+1 2	|∇p| 2 pγ	+ p	α+1 2 ∆p + p	α+1 2 R] =
	-[div(p	α+1 2 ∇p)-	α + 1 2	p	α+1 2	|∇p| 2 p	+αp	α+1 2	|∇p| 2 2p	] [p	α+1 2	|∇p| 2 pγ	-	α + 1 2	p	α+1 2	|∇p| 2 p	+div(p	α+1 2 ∇p)+p
									0	T	R d	div(p	α * +1 2	∇p) -p	α * +1 2	|∇p| 2 2p	2	≤ C(T ),	(46)
												p α * |∇p(t)| 2 dx ≤ C(T )	∀t ∈ (0, T ).	(47)
										R d

Proof. These two estimates, (46) and (47) come together and require some elaborate computations. We write

∂ t ∇p = ∇[|∇p| 2 + γp∆p + γpR], ∂ t |∇p| 2 2 = ∇p.∇[|∇p| 2 + γp∆p + γpR],

∂ t

from the National Science Centre (Poland§ B.P. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 740623) ¶ This work has been initiated when this author was visiting the Banach center in Warsaw during Simons Semester partially supported by the Simons -Foundation grant 346300 and the Polish Government MNiSW 2015-2019 matching fund.