G Rotondo
email: rotondo@recherche.enac.fr

Dr Paul Thevenon

Dr Carl Milner

Dr Christophe Macabiau

M Felux

Rotondo Giuseppe

Dr Achim Hornbostel

Simona Mihaela

Circiu

Methodology for determining Pseudorange noise and multipath models for a multi-constellation, multi-frequency GBAS system

in 2004 and obtained in 2007 a research master at ISAE in space telecommunications. In 2010, he obtained a Ph.D. degree in the signal processing laboratory of ENAC in Toulouse, France. From 2010 to 2013, he was employed by CNES, the French space agency, to supervise GNSS research activities and measurement campaigns. Since the July 2013, he is employed by ENAC as Assistant Professor. His current

INTRODUCTION

The ILS is used for guiding aircraft on final approach during precision approach operations at almost all the major airports worldwide. The ILS whilst not expensive to install requires frequent and expensive maintenance. Furthermore, the system can only support a straight-in approach trajectory for a single runway end, such that multiple installations for one airport are required to support multiple runways [START_REF] Dautermann | Approach service type D evaluation of the DLR GBAS testbed[END_REF]. Known issues also include multipath caused by uneven ground surface or other aircrafts on the airport surface which limit the capacity by restricting separation minima. Using the Global Navigation Satellite Systems (GNSS) with the GBAS could provide the safe and reliable guidance required with greatly improved flexibility in the definition of approach tracks. Moreover, GBAS could be a more cost effective solution since only one ground subsystem installation could be used to support multiple runway approaches at a single or potentially multiple aerodromes. The GBAS enhances the core constellation by providing differential corrections and integrity monitoring. Different types of services were developed classified with the acronym GBAS Approach Service Type (GAST), with GAST C supporting the Category (CAT) I precision approach type using Single Frequency (SF) and Single Constellation (SC) position solution. GAST D is designated for the Category II/III precision approach operations utilizing the L1 C/A signal of the Global Positioning System. With the advent of GAST D, a secondary shorter smoothing filter time constant of 30s has been introduced to limit the maximum residual differential ionospheric error, whilst GAST C is based on a 100s time constant [START_REF] Dautermann | Approach service type D evaluation of the DLR GBAS testbed[END_REF]. Despite this, the primary threat to GBAS users remains the gross ionospheric differential error induced by strong gradients. GBAS like all differential navigation systems relies on strong spatial correlation of errors such as this ionospheric delay between the ground reference stations and aircraft, such that they may be mitigated through the broadcast of corrections. However, extreme ionosphere storms, causing large ionosphere differential errors must be protected against and monitoring for such threats inevitably impacts on continuity and availability, potentially resulting in a degradation of service. One means to overcome this problem is to use the signals on Multiple Frequencies (MF) to form combinations which partially or totally remove the delays caused by the ionosphere [START_REF] Circiu | Evaluation of Dual Frequency GBAS Performance Using Flight Data[END_REF] [START_REF] Konno | Design of an aircraft landing system using dual-frequency GNSS[END_REF]. Furthermore, the additional satellites available in the Multiple-Constellation (MC) environment will significantly improve performance by adding geometric redundancy. Dual frequency techniques have been investigated in previous work [START_REF] Konno | Dual-Frequency Smoothing for CAT III LAAS: Performance Assessment Considering Ionosphere Anomalies[END_REF] [5], leading to two smoothing algorithms, Divergence Free (Dfree) and Ionosphere Free (I-free) smoothing. The differences between the two algorithms relate to the level of mitigation of the ionospheric delay and the resulting noise inflation of the final observables. The I-free technique removes the ionosphere delay in its entirety but at the cost of increased noise on the observable used for positioning. This is achieved through combining both the code and phase measurements on two frequencies. The Dfree technique removes only the part of the ionospheric delay relating to the transient temporal divergence, but no increase in the standard deviation of the noise and multipath over the single frequency smoothing output occurs [START_REF] Konno | Dual-Frequency Smoothing for CAT III LAAS: Performance Assessment Considering Ionosphere Anomalies[END_REF] [START_REF] Hwang | Enhanced DifferentialGPS Carrier-Smoothed Code Processing Using Dual-Frequency Measurements[END_REF].

These techniques can thus be used to mitigate the ionosphere and provide Cat II/III services when the GAST D service would be unavailable (under ionospheric gradient conditions or under poor geometry conditions). The SESAR 15.3.7 project is developing the GAST F concept through the investigation of these processing methodologies amongst others. In order to assess accurately the performance which may be achieved, the error model for the GALILEO E1 signal, the Galileo E5a signal and the GPS L5 signals must be determined. It is important to determine this firstly at the raw pseudorange level before addressing the impact of smoothing. Furthermore, different smoothing time constants and correction update rates are being considered within the SESAR framework which will require newly characterized models than those presented within the MOPS and SARPs [START_REF] Do ; A | Minimum aviation system performance standards for the local area augmentation system (LAAS)[END_REF] [START_REF]Development Baseline SARPs Proposal[END_REF]. In addition, new GAST F constraints regarding the antenna environment for the ground installation and on the tracking configuration on the receiver may be defined which would modify the impact of noise and multipath on the measurement.

CURRENT ERROR MODEL

The error model of a GBAS is composed by different errors contributing to the total error model. The non-aircraft contribution of the residual error is given in [START_REF] Do ; A | Minimum aviation system performance standards for the local area augmentation system (LAAS)[END_REF] and [START_REF] Icao | International Standards and Recommended Practices[END_REF] by the formula:

𝑅𝑀𝑆 𝑝𝑟 𝑔𝑟𝑜𝑢𝑛𝑑 (𝜃 𝑖 ) ≤ √ 1 𝑀 (𝑎 0 + 𝑎 1 𝑒 (- 𝜃 𝑖 𝜃 0 ) ) 2 + (𝑎 2 ) 2 (1)
Where:  𝑀 is the number of the receivers in the ground subsystem.  𝜃 𝑖 is the elevation angle for the 𝑖 𝑡ℎ ranging source.  The values of 𝑎 0 , 𝑎 1 , 𝑎 2 and 𝜃 0 are given in Table 1 The main contribution to this type of error is given by noise and multipath at the ground station, other sources of error are due to residual atmospheric error due to the physical separation between the ground station and the aircraft. In GAST F new possible processing options should be used and consequentially different smoothing time constant, the evaluation of the errors before the smoothing filter in this case can provide precious information to the user. The evaluation of the multipath plus noise error on L1 has been done using following the formula of the Code-Carrier (CMC) with Divergence-Free combination of the phase measurement [START_REF] Konno | Dual-Frequency Smoothing for CAT III LAAS: Performance Assessment Considering Ionosphere Anomalies[END_REF]. This particular combination of code and dual frequency phase measurement is composed mainly of the noise and multipath error on the code measurement. The equation of the CMC is:

𝐶𝑀𝐶 𝐿1 = 𝜌 1 -𝜙 1 + 2 𝛼 (𝜙 1 -𝜙 2 ) (2)
Where:  𝜌 1 is the code measurement on L1;  𝜙 1 and 𝜙 2 are the phase measurement on L1 and L2

 𝛼 = 1 - 𝑓 1 2 𝑓 2
2 , 𝑓 1 and 𝑓 2 are the frequencies of the signals on L1 and L2. Considering the model of the code and phase measurement as:

𝜌 1 = 𝑟 + 𝑐(𝑑𝑡 𝑠𝑣 + 𝑑𝑡 𝑟𝑒𝑐 ) + 𝑇 + 𝐼 1 + 𝜀 𝜌1 (3) 𝜙 1 = 𝑟 + 𝑐(𝑑𝑡 𝑠𝑣 + 𝑑𝑡 𝑟𝑒𝑐 ) + 𝑇 -𝐼 1 + 𝑁 1 + 𝜀 𝜙1 (4) 𝜙 2 = 𝑟 + 𝑐(𝑑𝑡 𝑠𝑣 + 𝑑𝑡 𝑟𝑒𝑐 ) + 𝑇 -𝐼 2 + 𝑁 2 + 𝜀 𝜙2 (5)
Where:

 𝑟 is the true range;  dt sv/rec are the satellite and receiver clock errors;  𝑇 is the troposheric delay;  𝐼 1 and 𝐼 2 are the ionospheric delay related to the frequency;  𝑁 1 and 𝑁 2 are the phase ambiguities related to each frequency;  𝜀 is the noise and multipath error on the code or on the phase measurement, according to the frequency. Replacing the eq. ( 3), ( 4) and (5) in eq. ( 2) and removing all the terms that are common in the three models, the remaining terms are:

𝐶𝑀𝐶 𝜌1 = (𝐼 1 + 𝜀 𝜌1 ) -(-𝐼 1 + 𝑁 1 + 𝜀 𝜙1 - 2 𝛼 (-(𝐼 1 -𝐼 2 ) + (𝑁 1 -𝑁 2 ) + (𝜀 𝜙1 -𝜀 𝜙2 ))) (6)
Considering the relation between the ionospheric delay and the frequencies used:

𝐼 1 = 𝐼 𝑓 1 2 ; 𝐼 2 = 𝐼 𝑓 2 2 (7)
It is possible to compute the following relation:

𝐼 1 -𝐼 2 = (1 - 𝑓 1 2 𝑓 2 2 ) 𝐼 1 = 𝛼 𝐼 1 (8)
It is possible to replace it in eq. ( 6) and simplify the common terms in order to obtain:

𝐶𝑀𝐶 𝜌1 = (𝐼 1 + 𝜀 𝜌1 ) -(𝐼 1 + 𝑁 12 + 𝜀 𝜙12 ) = 𝜀 𝜌1 -(𝑁 12 + 𝜀 𝜙12 ) (9)
Where:

 𝑁 12 = 𝑁 1 - 2 𝛼 (𝑁 1 -𝑁 2 )  𝜀 𝜙12 = 𝜀 𝜙1 - 2 𝛼 (𝜀 𝜙1 -𝜀 𝜙2 )
In eq. ( 9) the first term 𝜀 𝜌1 represents the noise and multipath affecting the code measurement, this error is assumed to be zero mean over long period. The second term is the phase ambiguity combination, it is a constant values and it is easily removable considering that the noise and multipath are zero mean over long period. The last term is a combination of noise and multipath on the phase measurement, as for the noise and multipath on code measurement this error is assumed to be zero mean over long period; this error can be considered as negligible if compared with the same error on the code measurement.

In case of SF data is not available to compute the CMC as in eq. ( 2), one possible combination is:

𝐶𝑀𝐶 𝜌1 = 𝜌 1 -𝜙 1 = 𝜀 𝜌 + 𝜀 𝜙 + 𝑁 1 + 2𝐼 1 (10
) It is possible to see that now the CMC as computed in eq. ( 10) is affected by the ionospheric delay multiplied by 2. To remove this term from the CMC, it can be computed using the DF measurements from a nearby station; considering the models of the phase measurement and removing the common terms:

𝜙 1 -𝜙 2 = -(𝐼 1 -𝐼 2 ) + (𝜀 𝜙1 -𝜀 𝜙2 ) (11)
Knowing the relation of the ionospheric delay on different frequencies as given in eq. ( 7), it is possible to replace these relations in eq. ( 11) and it is possible to find the following relationship:

𝜙 1 -𝜙 2 = -𝐼 ( (𝑓 2 2 -𝑓 1 2 ) 𝑓 1 2 𝑓 2 2 ) + (𝜀 𝜙1 -𝜀 𝜙2 ) ( 12
)
Where 𝐼 is the ionospheric delay not related to any frequency.

Just multiplying all the terms by 𝑓 2 2 /(𝑓 1 2 -𝑓 2 2 ), it is possible to obtain:

𝑓 2 2 𝑓 1 2 -𝑓 2 2 (𝜙 1 -𝜙 2 ) = 𝐼 𝑓 1 2 + 𝑓 2 2 𝑓 1 2 -𝑓 2 2 (𝜀 𝜙1 -𝜀 𝜙2 ) (13)
In eq. ( 13) the first term,

𝐼 𝑓 1
2 is the ionospheric delay on the same frequency as the code measurement; the last one is the difference of noise and multipath on the phase measurement multiplied by a term related to the two used frequency; this second term is negligible if compared to the ionospheric delay.

It is therefore possible to evaluate the term 𝐼 1 from a dual frequency phase combination of a nearby station, and to remove the influence of the ionosphere in eq. [START_REF] Wang | Development of a Low-Cost GPSbased Attitude Determination System[END_REF].

When applying directly one of the two formula for the CMC computation, it is possible to see some large errors affecting it, which is not representing the multipath or the noise error, but is in fact due to phase ambiguity in Eq. ( 2). This effect is visible in Figure 1.

Figure 1 -Multipath and Noise Evaluation on L1

In order to have the correct evaluation of the multipath plus noise impact on the L1 and L5 signals, the following work has been done:

 Identify groups of data corresponding to a continuous tracking of the signal.  Search inside each group of data for possible cycle slip comparing the predicted phase measurement, eq. ( 14), and the real one [START_REF] Wang | Development of a Low-Cost GPSbased Attitude Determination System[END_REF]:

Φ ̂𝑘 = Φ 𝑘-1 + Φ ̇𝑘 + Φ ̇𝑘-1 2 Δ𝑡 (14)
Where  Φ ̂ is the predicted phase measurement  Φ is the real phase measurement.  Φ ̇ is the Doppler measurement  𝑘 is the epoch index  Δ𝑡 is the interval between the 𝑘 𝑡ℎ measurement and the previous one. If the absolute difference between the real phase measurement and the predicted phase is bigger than 1 cycle or 1 wavelength, a cycle slip is detected.

If no Doppler measurement are present in the data another methodology must be used. The computation of the double-differences on the phase measurement can be used to detect cycle slip [START_REF] Kannemans | Cycle slip detection for static and moving receivers[END_REF]:

𝑑 𝑖 = 𝜙 𝐼+1 -𝜙 𝑖 𝑑𝑑 𝑖 = 𝑑 𝑖 - 𝑑 𝑖+1 + 𝑑 𝑖-1 2 (15)
The cycle slip is detected by comparing the value of the real and predicted phase measurement, the threshold has been set at 1 cycle. For the case with the double-differences, the threshold is 0.5. From each CMC series, according to the group of data, the mean value of the CMC is removed in order to compensate for the possible phase ambiguity over the measurement. It is important to note that because noise and multipath are assumed to be zero mean over long period, the slices of CMC with less than 3000 samples are not taken into account. In Figure 2 it is possible to see that the cycle slips and the phase ambiguity are removed from the CMC.

Figure 2 -Multipath & Noise on L1 Corrected for the Phase Ambiguity

Once that the CMC has been computed, the modelling of the error for each satellite and for a 1° elevation bin is done.

An auto-regressive (AR) model is used because it is a more realistic representation of the noise plus multipath signal, in fact it takes into account the correlation time of the data. The of the signal just using sigma as representation of the standard deviation and considering it as white noise is a too optimistic representation for processing purposes. A generic AR model is represented by the following formula [START_REF] Jackson | Digital Filters and Signal Processing[END_REF]:

𝑥(𝑛) = -∑ 𝑎 𝑘 𝑥(𝑛 -𝑘) + 𝑝 𝑘=1 𝑏(𝑛) (16)
Where  𝑥 is the signal at different lags.  𝑎 𝑘 are the coefficients according to order model.  𝑏(𝑛) is a Gaussian noise, named the driving noise. The goal of the modelling process is to estimate the 𝑎 parameters, according to the model order, and the variance or the standard deviation of the driving noise for each elevation bin. Different Matlab functions are used in order to compute the 𝑎 parameters and the variance of the driving noise, then a comparison between the function has been done on order to check possible differences. They are all computing the same parameters, but using different methods. LPC (data, order). This function find the coefficients of Norder forward linear prediction and the variance of the driving noise. ARCOV (data,order) Estimate AR model parameters using covariance method and the variance of the driving noise. ARYULE (data,order) Estimate autoregressive (AR) allpole model using Yule-Walker [START_REF] Hayes | Statistical Digital Signal Processing and Modeling[END_REF] method and variance of the driving noise.

ARMCOV (data,order) Estimate AR model parameters using modified covariance method and compute the variance of the driving noise. ARBURG (data,order) Estimate AR model parameters using Burg [START_REF] Kay | Modern Spectral Estimation: Theory and Application[END_REF] method and compute the variance of the driving noise. These different methods all estimate the parameters of an AR model, but use different computation techniques. The following investigation aims at determining if one of these methods is more appropriate to the targeted kind of data.

RESULTS

The methodology explained in the previous chapter has been applied on two series of measurement: 1. Data collected at Braunschweig airport by DLR the 9th July 2014. The first observation available is at 00h 00' 00" and the last observation is at 23h 59' 59", the interval between the observations is 0.5 seconds. The antenna model is a Leica AR-25 choke ring antenna the stations are located at the positions given in Table 2.

Data Collected at Malaga airport GBAS Station

(courtesy of ENAIRE) the 30 th March 2014. The first observation available is at 00h 00' 14" and the last one is the day after at 00h 00' 14". The measurement are recorded at 2 Hz only for the GPS L1 C/A signal, data coming from the near monitoring station have been also used to compute the ionospheric delay. The antenna, for the GBAS data, is a Multipath Limiting Antenna (MLA).

DLR data collection

The first step is the choice of the order model, it derives from the analysis of three different criteria that provide an estimation of the prediction error power [START_REF] Ljung | System Identification: Theory for the User[END_REF] Figure 3 shows that the practical minimum order is around 5, however analyzing also the errors magnitude for the lower model order has been found that the difference of the error criteria between the second and the fifth order model is small (approximately between 6-7 % for each criteria). Therefore, a 2 nd order model seems reasonably close to represent the real signal, and has the benefit to keep the model simpler than a higher order model. In Figure 4, the comparison of the coefficients standard deviation confirms the choice of the model. The standard deviation of the a1 coefficients for the 1 st order model is larger than the other model orders, and from the second order, the standard deviation starts to be almost similar. The values of the a2 standard deviation are similar for all the order model analyzed. The standard deviation of the driving noise is similar for all the model orders analyzed.

From the previous analysis the second order model has been selected to model the noise and multipath, the equation that model the noise and multipath is:

𝑥(𝑛) = -𝑎 1 * 𝑥(𝑛 -1) -𝑎 2 * 𝑥(𝑛 -2) + 𝑏(𝑛) (20)
In the next figures the values of the 𝑎 coefficients and the standard deviation of the driving noise process will be show, the red lines represent the proposed coefficients and the proposed model for the noise standard deviation. Each Elevation Bin From the analysis of the previous figures, it is possible to assess that for the analyzed data there are no differences in the coefficients computed with the different Matlab functions and even the standard deviation is similar for all of them; for this reason in the next figures just the ARMCOV function will be considered. The 𝑎 1 coefficient model proposed is 0.41, instead the value of 𝑎 2 coefficient model is 0.25. The standard deviation of the driving process noise follows an elevation dependent exponential curves, the model proposed for it is:

𝜎 𝑑𝑟_𝑛𝑜𝑖𝑠𝑒 = 0,7 * 𝑒𝑥𝑝 - 𝐸𝑙 14,3 + 0,23 (21)
And it is represented by the red dashed line.

In order to be sure that the computed coefficients and variance of the driving noise are consistent, the real signal and the one generated by the use of the model parameters are compared in Figure 8, and an analysis if the autocorrelation function is shown in Figure 9. These plots show that the coefficients of the model reflects the properties of the real signal. The analysis of Figure 11 shows that the practical model order is at 10, but analyzing also the lower model order error magnitude, it is possible to see that the magnitude of the three criteria for the first model order is at most 3% bigger than the one provided by model order 10.

The choice of the 1 st model order seems to be appropriate in order to have a reasonably accurate model but not too complex. The analysis of Figure 12 shows that the standard deviation of the 𝑎 1 coefficients is really similar for all the model order and also the driving noise standard deviation has almost the same values for all the order models. Thanks to these two analysis is possible to select 1 as model order.

The equation of the autoregressive model, for this set of data, is:

𝑥(𝑛) = -

CONCLUSION

We have presented a proposal for the estimation of noise and multipath using a D-Free combination. In order to remove possible phase ambiguity and cycle slip on the phase a phase based cycle slip detector has been implemented, finally the CMC was modelled using an autoregressive method. The application of this methodology to two different data collections using only single ground receiver, has provided good results; however it is important to note that the results provided for the DLR data collection are not representative of a certificated GBAS ground station because they are not conform with the siting antenna and receiver constraints. One advantage of this method is that it provides more information about the noise and multipath temporal correlation than just describing the time series by its total standard deviation.

FUTURE WORKS

The next steps of this work consists in the application of the methodology to a large set of data and in particular to GPS L1 and L5 and GALILEO E1 and E5a, which will be the signal of interest of a GAST F GBAS system. The second step will be the analysis of the same signal applying the I-Free combination and, finally, the analysis of the two combinations after the smoothing filter using different smoothing time constants in order to verify their impact on the noise and multipath pseudorange error. The correct model of noise and multipath standard deviation vs. elevation will then be used for other measurement processing in the GAST F GBAS system, such as the computation of the Protection Levels.

Figure 3 - 5 Figure 4 -

 354 Figure 3 -Error Criteria for Satellite PRN 5

Figure 5 -Figure 6 -Figure 7 -

 567 Figure 5 -A1 Coefficients: Mean Values and Related Standard Deviation

Figure 8 -Figure 9 -Figure 10 -

 8910 Figure 8 -Real CMC and Generated CMC for Satellite PRN 2 and 30° Elevation Bin

Figure 11 - 3 Figure 12 -

 11312 Figure 11 -Error Criteria for Satellite PRN 3

  𝑎 1 * 𝑥(𝑛 -1) + 𝑏(𝑛) (22)The next plots will show the 𝑎 1 coefficients mean and the driving noise standard deviation computed for all the elevation angle between 10 and 90.

Figure 13 -Figure 14 -

 1314 Figure 13 -A1 Coefficients: Mean Values and Related Standard Deviation

Figure 15

 15 Figure 15 -Real CMC and Generated CMC for Satellite PRN 18 and 30° Elevation Bin

Figure 18 -Figure 19 -

 1819 Figure 18 -Real CMC and Generated CMC for Satellite PRN 2 and 30° Elevation Bin with a 20 th order model

Figure 20 -Figure 21 Figure 22 -

 202122 Figure 20 -Real CMC and Generated CMC for Satellite PRN 18 and 30° Elevation Bin with a 20 th order model

Table 1 -Non-Aircraft Elements Accuracy Requirement [6].

 1

	GAD	𝜽 𝒊 (degrees)	𝒂 𝟎	𝒂 𝟏	𝒂 𝟐	𝜽 𝟎 (degrees)
	A	>5	0.5	1.65	0.08	14.3
	B	>5	0.16	1.07	0.08	15.5
	C	>35 <35	0.15 0.24	0.84 0	0.04	15.5 //
	NOISE AND MULTIPATH EVALUATION
	The development of a new GBAS service, GAST F,
	requires the knowledge of the error models for each signal
	used over the two constellations. The model proposed for
	until now is adapted only for the GPS L1 signal and
	considers a smoothing filter time constant of 100 seconds.

Table 2 -Reference Receiver Locations

 2

	Indicator	Latitude[°]	Longitude[°]	Receiver
				Type
	BR01	52.321444	10.543339 E	Javad
		N		Delta
	BR02	52.322324	10.554618 E	Javad
		N		Delta
	BR03	52.317001	10.567265 E	Javad
		N		Delta
	BR04	52.321269	10.564586 E	Javad
		N		Delta

ACKNOWLEDGMENTS

The authors would like to thank DLR and INDRA/ENAIRE for provision of data and for supporting the work.

"©SESAR JOINT UNDERTAKING, 2013. Created by [ENAC, DLR] for the SESAR Joint Undertaking within the frame of the SESAR Programme co-financed by the EU and EUROCONTROL. The opinions expressed herein reflects the author's view only.

The SESAR Joint Undertaking is not liable for the use of any of the information included herein. Reprint with approval of publisher and with reference to source code only."

Dr. Michael Felux obtained a diploma in Technical

Mathematics from TU München in 2009since then. Currently, he is a Ph.D. candidate at TU München, Chair of Flight


Sepideh Mirrahimi
email: sepideh.mirrahimi@math.univ-toulouse.fr

Jean-Michel Roquejoffre
email: jean-michel.roquejoffre@math.univ-toulouse.fr

A class of Hamilton-Jacobi equations with constraint: uniqueness and constructive approach

Keywords: Hamilton-Jacobi equation with constraint, uniqueness, constructive existence result, selection-mutation models AMS Class. No: 35A02, 35F21, 35Q92

We discuss a class of time-dependent Hamilton-Jacobi equations, where an unknown function of time is intended to keep the maximum of the solution to the constant value 0. Our main result is that the full problem has a unique viscosity solution, which is in fact classical. The motivation is a selection-mutation model which, in the limit of small diffusion, exhibits concentration on the zero level set of the solution of the Hamilton-Jacobi equation. Uniqueness is obtained by noticing that, as a consequence of the dynamic programming principle, the solution of the Hamilton-Jacobi equation is classical. It is then possible to write an ODE for the maximum of the solution, and treat the full problem as a nonstandard Cauchy problem.

Introduction 1.Model and question

The purpose of this paper is to discuss existence and uniqueness for the following problem, with unknowns (I(t), u(t, x)):

       ∂ t u = |∇u| 2 + R(x, I) (t > 0, x ∈ R d ), max
x u(t, x) = 0

I(0) = I 0 > 0, u(0, x) = u 0 (x), (1)
where I 0 > 0 and u 0 is a concave, quadratic function. For a given continuous function I(t), u(t, x) solves a Hamilton-Jacobi equation. The unknown I may be thought of as a sort of regulator, or a sort of Lagrange multiplier, to maintain the maximum of u equal to 0. The constraint on the maximum of u(t, .) makes the problem nonstandard.

1

Existence of a solution (u, I) to [START_REF] Barles | Solutions de viscosité des quations de Hamilton-Jacobi[END_REF] is not new, the first result is due to Perthame and Barles [START_REF] Perthame | Dirac concentrations in Lotka-Volterra parabolic PDEs[END_REF] (see also for a result with weaker assumptions). An important improvement is given by Lorz, Perthame and the first author in [START_REF] Lorz | Dirac mass dynamics in a multidimensional nonlocal parabolic equation[END_REF]; they indeed notice that a concavity assumptions on R -that we also make here -entail regularity. This allows them to derive the dynamics of the maximum point of a solution u(t, x). See also [START_REF] Mirrahimi | Phénomènes de concentration dans certaines EDPs issues de la biologie[END_REF]. Both types of results rely on a special viscous approximation of (1) -see equation ( 4) below. Uniqueness, however, has remained an open problem, apart from a very particular case [START_REF] Perthame | Dirac concentrations in Lotka-Volterra parabolic PDEs[END_REF].

The main goal of the paper is to prove the missing uniqueness property; a result that we had already announced in [START_REF] Mirrahimi | Uniqueness in a class of Hamilton-Jacobi equations with constraints[END_REF]. We also provide a constructive existence proof which was not available in the previous existence results [START_REF] Perthame | Dirac concentrations in Lotka-Volterra parabolic PDEs[END_REF][START_REF] Barles | Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result[END_REF][START_REF] Lorz | Dirac mass dynamics in a multidimensional nonlocal parabolic equation[END_REF]. Two important consequences, that we will present in a forthcoming paper [START_REF] Mirrahimi | Approximation of solutions of selection-mutation models and error estimates[END_REF] (see also [START_REF] Mirrahimi | Uniqueness in a class of Hamilton-Jacobi equations with constraints[END_REF]), will be the convergence of the underlying selection-mutation model in a stronger sense than what is known, and asymptotic expansion of the viscous solution. The asymptotic expansion, which allows to approximate the phenotypical distribution of the population when the mutation steps are small but nonzero, is particularly interesting in view of biological applications. One of the main ingredients will be regularity under suitable concavity assumptions on R and u 0 , which is far from being available in general. Instead of relying on viscous approximations we will prove these results directly for the equation

u t = |∇u| 2 + R(t, x). (2)
This will allow a much easier treatment than in the usual viscosity sense.

The uniqueness result will also be helpful to develop the so-called Hamilton-Jacobi approach (see for instance Diekmann et al. [START_REF] Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF], [START_REF] Perthame | Dirac concentrations in Lotka-Volterra parabolic PDEs[END_REF][START_REF] Lorz | Dirac mass dynamics in a multidimensional nonlocal parabolic equation[END_REF] and Subsection 1.2) to study more complex models describing selection and mutations. For instance, our result would allow to generalize a result due to Perthame and the first author in [START_REF] Mirrahimi | Asymptotic analysis of a selection model with space[END_REF] on a selection model with spatial structure, where the proof relies on the uniqueness of the solution to a corresponding Hamilton-Jacobi equation with constraint.

Motivation

Model (1) arises in the limit ε → 0 of the solutions to the problem

∂ t n ε -ε∆n ε = n ε ε R x, I ε (t) (t > 0, x ∈ R d ), I ε (t) = R d ψ(x)n ε (t, x)dx, (3)
where n ε (t, x) is the density of a population characterized by a d-dimensional biological trait x. The population competes for a single resource, this is represented by I ε (t), where ψ is a given positive smooth function. The term R(x, I) is the reproduction rate; it is, as can be expected, very negative for large x and decreases as the competition increases. The Laplace term corresponds to the mutations. The small parameter ε is introduced to consider the long time dynamics of the population when the mutation steps are small. Such models can be derived from individual based stochastic processes in the limit of large populations (see Champagnat-Ferrière-Méléard [START_REF] Champagnat | Unifying evolutionary dynamics: From individual stochastic processes to macroscopic models[END_REF][START_REF] Champagnat | Individual-based probabilistic models of adaptive evolution and various scaling approximations[END_REF]). There is a large literature on the models of population dynamics under selection and mutations. We refer the interested reader, for instance to Geritz et al. [START_REF] Geritz | Evolutionary singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF] and Diekmann [START_REF] Diekmann | Beginner's guide to adaptive dynamics[END_REF] for an approach based on the study of the stability of differential systems (the so-called adaptive dynamics approach), to Champagnat [START_REF] Champagnat | Mathematical study of stochastic models of evolution belonging to the ecological theory of adaptive dynamics[END_REF] for the study of stochastic individual based models, to Raoul [START_REF] Raoul | Etude qualitative et numérique d'équations aux dérivées partielles issues des sciences de la nature[END_REF] and Mirrahimi [START_REF] Mirrahimi | Phénomènes de concentration dans certaines EDPs issues de la biologie[END_REF] for the study of integro-differnetial models.

The Hopf-Cole transformation n ε = exp (u ε /ε) yields the equation

∂ t u ε = ε∆u ε + |∇u ε | 2 + R(x, I ε ) (4)
which, in the limit ε → 0, yields the equation for u. Furthermore, I ε being uniformly positive and bounded in ε, the Hopf-Cole transformation leads to the constraint on u. One expects that n ε concentrates at the points where u is close to 0 and the function I ε appears, in the limit, as a sort of Lagrange multiplier. One has indeed

n ε (x, t) --⇀ ε→0 n(x, t) = ρ(t) δ(x -x(t))
, weakly in the sense of measures, with

u(t, x(t)) = max x u(t, x) = 0, ρ(t) = I(t) ψ(x(t))
.

This method to study (3) has been introduced in [START_REF] Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF] and then developed in different contexts. See for instance Perthame-Barles [START_REF] Perthame | Dirac concentrations in Lotka-Volterra parabolic PDEs[END_REF] (convergence to Hamilton-Jacobi dynamics for (4), Barles-Mirrahimi-Perthame [START_REF] Barles | Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result[END_REF] (the same type of result, but with nonlinear, nonlocal diffusion), Champagnat-Jabin [START_REF] Champagnat | The evolutionary limit for models of populations interacting competitively via several resources[END_REF] (nonlinear integro-differential model with several resources), [START_REF] Lorz | Dirac mass dynamics in a multidimensional nonlocal parabolic equation[END_REF] (convergence improvement by introduction of the concavity assumptions). This approach has a lot to do with the 'approximation of geometric optics' for reaction-diffusion equations of the Fisher-KPP type. See Freidlin [START_REF] Freidlin | Functional integration and partial differential equations[END_REF][START_REF] Freidlin | Limit theorems for large deviations and reaction-diffusion equations[END_REF] for the probabilistic approach, and Evans-Souganidis [START_REF] Evans | A PDE approach to geometric optics for certain semilinear parabolic equations[END_REF], Barles-Evans-Souganidis [START_REF] Barles | Wavefront propagation for reaction-diffusion systems of PDE[END_REF] for the viscosity solutions approach.

Assumptions

The assumptions we are stating below are in the same spirit (but slightly weaker) as in [START_REF] Lorz | Dirac mass dynamics in a multidimensional nonlocal parabolic equation[END_REF], where the authors noticed that this set of assumptions allowed them to work with smooth solutions, thus going quite far in the study of [START_REF] Barles | Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result[END_REF]. We believe that the results that we will prove would certainly be false if some of those assumptions were removed.

• Assumptions on R(x, I). We choose R to be smooth, and we suppose that there is I M > 0 such that (fixing the origin in x appropriately) max

x∈R d R(x, I M ) = 0 = R(0, I M ), (5)
-K 1 |x| 2 ≤ R(x, I) ≤ K 0 -K 1 |x| 2 , for 0 ≤ I ≤ I M , (6)
-2K 1 ≤ D 2 R(x, I) ≤ -2K 1 < 0 as symmetric matrices, (7)
-K 2 ≤ ∂R ∂I ≤ -K 2 , (8)
| ∂ 2 R ∂I∂x i (x, I)| + | ∂ 3 R ∂I∂x i ∂x j (x, I)| ≤ K 3 , for 0 ≤ I ≤ I M , and i, j = 1, 2, • • • , d, (9)
D 3 R(•, I) L ∞ (R d ) ≤ K 4 , for 0 ≤ I ≤ I M . (10)
• Assumptions on u 0 (.) and I 0 . We assume the existence of positive constants

L 0 , L 0 , L 1 , L 1 such that -L 0 -L 1 |x| 2 ≤ u 0 (x) ≤ L 0 -L 1 |x| 2 , (11)
-2L 1 ≤ D 2 u 0 ≤ -2L 1 . (12)
Note that this implies

|Du 0 (x)| ≤ L 2 (1 + |x|), (13)
for a large constant L 2 > 0. We also need that, for a positive constant L 3 ,

D 3 u 0 L ∞ (R d ) ≤ L 3 . (14)
Finally we assume that max

x u 0 (x) = u 0 (x 0 ) = 0, R(x 0 , I 0 ) = 0. ( 15
)
Note that the monotony assumption [START_REF] Diekmann | Beginner's guide to adaptive dynamics[END_REF] means that the growth rate decreases as the competition increases, which is natural from the modeling point of view. The concavity assumption is a technical one.

In Section 2 we will study an unconstrained Hamilton-Jacobi equation where we replace R(x, I) by R(t, x). To prove our results on this unconstrained problem we assume same type of regularity and concavity assumptions on R that we state below:

• Assumptions on R(t, x). We choose R to be smooth, and we suppose that

-K 1 |x| 2 ≤ R(t, x) ≤ K 0 -K 1 |x| 2 , for t ∈ R + , (16)
-2K 1 ≤ D 2 R(t, x) ≤ -2K 1 < 0 as symmetric matrices, (17)
D 3 R(t, •) L ∞ (R d ) ≤ K 4 , for t ∈ R + . (18)

Results and plan of the paper

Our first result concerns the unconstrained Hamilton-Jacobi equation

u t = |∇u| 2 + R(t, x) (t > 0, x ∈ R d ), u(0, x) = u 0 (x). ( 19
)
The assumptions on u 0 and R are those stated in the preceding subsection.

Theorem 1.1 Equation (19) has a unique viscosity solution u that is bounded from above. Moreover, it is a classical solution:

u ∈ L ∞ loc R + ; W 3,∞ loc (R d ) ∩ W 1,∞ loc R + ; L ∞ loc (R d ) , -max(2 L 1 , K 1 ) ≤ D 2 u ≤ -min(2L 1 , K 1 ) and D 3 u L ∞ ([0,T ]×R d ) ≤ L 4 (T ) where L 4 (T ) is a positive constant depending on L 1 , K 1 , K 4 , L 3 and T .
Let us point out that the assumption of boundedness from above is, most certainly, irrelevant. However the constraint in [START_REF] Barles | Solutions de viscosité des quations de Hamilton-Jacobi[END_REF] implies that all the solutions that we consider are bounded from above. This extra condition is thus legitimate, and will keep the length of the preliminary work to a minimum.

Theorem 1.2

The Hamilton-Jacobi equation with constraint (1) has a unique solution (u, I). Moreover we have

(u, I) ∈ L ∞ loc R + ; W 3,∞ loc (R d ) ∩ W 1,∞ loc R + ; L ∞ loc (R d ) × W 1,∞ (R).
The paper is organised as follows. In Section 2 we prove the Cauchy Problem for [START_REF] Mirrahimi | Uniqueness in a class of Hamilton-Jacobi equations with constraints[END_REF]. In Section 3 we reduce (1) to a (nonstandard) differential system. Theorem 1.2 is proved in Section 4. Section 5 is devoted to the study of a particular example.

The Cauchy problem

In this section, we prove Theorem 1.1.

We will first prove that the only solution to (2) that is bounded from above is the solution u(t, x) of the dynamic programming principle

u(t, x) = sup (γ(s),s)∈R d ×[0,t] γ(t)=x F (γ) : γ ∈ C 1 ([0, t]; R d ) , (20)
with

F (γ) := u 0 (γ(0)) + t 0 - | γ| 2 4 (s) + R(s, γ(s)) ds.
Such a solution is a viscosity solution to (1). We will prove, in addition, that it is classical and satisfies the properties claimed by Theorem 1.1.

Uniqueness. This step essentially consists in showing that a viscosity solution of (2) does not grow too wildly, which will reduce the problem to the application of classical arguments. We have already assumed boundedness from above, so let us show that u goes to -∞ at most in a quadratic fashion. Due to the assumptions ( 16) on R, we have

∂ t v ≥ 0, v(t, x) = u(t, x) + tK 1 |x| 2
in the viscosity sense, which implies that v is time-increasing -thus the needed estimate. Let usalthough this is elementary -explain why: choose T > 0 and assume the existence of 0 < s < t ≤ T such that the inequality v(t, x) ≥ v(s, x) does not hold. In other words there is x 0 such that v(s, x 0 ) > v(t, x 0 ). For ε > 0 consider the quantity

w ε (t, x) := v(t, x) -v(s, x 0 ) + ε T -t + |x -x 0 | 2 ε 2 .
For ε > 0 small enough, there is a local minimum point for w ε , called (t ε , x ε ), such that t ε is bounded away from s and T , and x ε → x 0 as ε → 0. At that minimum point the viscosity inequality implies

-∂ t 1 T -t ≥ 0, a contradiction.
So there is a large C > 0 such that

-C(1 + t)(1 + |x| 2 ) ≤ u(t, x) ≤ C.
And so, by an easy adaptation of Chap. 2 of Barles [START_REF] Barles | Solutions de viscosité des quations de Hamilton-Jacobi[END_REF], where a uniqueness result for solutions that grow at most exponentially fast is provided (see also [START_REF] Fujita | Asymptotic solutions of viscous Hamilton-Jacobi equations with Ornstein-Uhlenbeck operator[END_REF][START_REF] Fujita | Asymptotic solutions of Hamilton-Jacobi equations in Euclidean n space[END_REF]), there is at most one viscosity solution to (2) that is bounded from above.

Existence. We may thus turn to [START_REF] Mirrahimi | Approximation of solutions of selection-mutation models and error estimates[END_REF]. Let us suppose that (γ n ) 1≤n , with

γ n ∈ C 1 ([0, t]; R d ) and γ n (t) = x, is such that F (γ n ) → u(t, x)
as n → ∞. Since R and u 0 are bounded from above, we obtain that, for some constant

C t 0 | γn | 2 (s)ds < C.
Consequently, from γ n (t) = x we deduce, modifying the constant C if necessary, that

γ n W 1,2 [0,t] < C. It follows that, there exists γ ∈ W 1,2 ([0, t]; R d ), such that as n → ∞, γ n → γ strongly in C([0, t]; R d )
and weakly in W 1,2 ([0, t]; R d ). We deduce that, as n → ∞,

u 0 (γ n (0)) → u 0 (γ), t 0 R(s, γ n (s))ds → t 0 R(s, γ(s))ds, t 0 | γ| 2 (s)ds ≤ lim inf n→∞ t 0 | γn | 2 (s)ds.
We conclude that

u(t, x) = u 0 (γ(0)) + t 0 - | γ| 2 4 (s) + R(s, γ(s)) ds. (21)
We claim that such a trajectory is unique, which entails uniqueness for the Cauchy problem [START_REF] Mirrahimi | Uniqueness in a class of Hamilton-Jacobi equations with constraints[END_REF]. We note indeed that such trajectory γ satisfies the following Euler-Lagrange equation

     γ(s) = -2∇R(s, γ(s)), γ(0) = -2∇u 0 (γ(0)), γ(t) = x. (22)
However, from the concavity assumptions on R and u 0 we obtain that the above elliptic problem is coercive and hence the solution γ is unique.

Regularity. Let us denote by γ x (t) the unique solution of [START_REF] Raoul | Etude qualitative et numérique d'équations aux dérivées partielles issues des sciences de la nature[END_REF]. The function

(t, x) → (γ x (t), γx (t)) belongs to L ∞ loc (R + , W 2,∞ loc (R d )) because ∇u 0 and ∇R(t, .) are in W 2,∞ loc (R d ). Now, we have u(t, x) = u 0 (γ x (0)) + F (γ x ),
yielding (this is a classical computation):

∂ i u(t, x) = ∇u 0 (γ x (0)).∂ i γ x (0) + t 0 - γx (s).∂ i γx (s) 2 + ∇R(s, γ x (s))∂ i γ x (s) ds.
Integrating by parts and using the Euler-Lagrange equation ( 22), we get:

∇u(t, x) = - γx (t) 2 . ( 23
) Thus u ∈ L ∞ loc (R + , W 3,∞ loc (R d )).
Strict concavity. We will prove that u(t, x) is uniformly strictly concave, namely that D 2 u ≤ -2λI

in the sense of symmetric matrices, for λ = min L 1 , √

K 1 2
. To this end, we show that, for all σ ∈ [0, 1] and (x,

y) ∈ R d × R d : σu(t, x) + (1 -σ)u(t, y) + λσ(1 -σ)|x -y| 2 ≤ u(t, σx + (1 -σ)y). ( 24
)
Let γ x and γ y be optimal trajectories, solving [START_REF] Raoul | Etude qualitative et numérique d'équations aux dérivées partielles issues des sciences de la nature[END_REF], with γ x (t) = x and γ y (t) = y. Note from the choice of γ x and γ y that we have

u(t, x) = u 0 (γ x (0)) + t 0 - | γx (s)| 2 4 + R(s, γ x (s)) ds, u(t, y) = u 0 (γ y (0)) + t 0 - | γy (s)| 2 4 + R(s, γ y (s)) ds, and
u(t, σx + (1 -σ)y) ≥ u 0 (σγ x (0) + (1 -σ)γ y (0)) + t 0 - |σ γx + (1 -σ) γy (s)| 2 4 + R(s, σγ x (s) + (1 -σ)γ y (s)) ds.
Furthermore, from the concavity assumptions on R and u 0 we have

σu 0 (t, x) + (1 -σ)u 0 (t, y) + L 1 σ(1 -σ)|γ x (0) -γ y (0)| 2 ≤ u 0 (t, σx + (1 -σ)y), and
σ t 0 R(s, γ x (s))ds + (1 -σ) t 0 R(s, γ y (s))ds + K 1 σ(1 -σ) t 0 |γ x (s) -γ y (s)| 2 ds ≤ t 0 R(s, σγ x (s) + (1 -σ)γ y (s))ds.
Moreover, from the strict concavity of µ → -|µ| 2 , we obtain that

σ t 0 - | γx (s)| 2 4 ds + (1 -σ) t 0 - | γy (s)| 2 4 ds + σ(1 -σ) t 0 | γx (s) -γy (s)| 2 4 ds ≤ t 0 - |σ γx + (1 -σ) γy (s)| 2 4 ds.
We deduce that

u(t, σx + (1 -σ)y)) ≥ σu(t, x) + (1 -σ)u(t, y) + σ(1 -σ) t 0 ( 1 4 | γx (s) -γy (s)| 2 + K 1 |γ x (s) -γ y (s)| 2 )ds + L 1 |γ x (0) -γ y (0)| 2 (25)
Next we have

K 1 2 t 0 d ds |γ x (s) -γ y (s)| 2 ds ≤ K 1 t 0 |γ x (s) -γ y (s)| 2 ds + t 0 | γx -γy | 2 4 (s)ds. Writing |x -y| 2 = |γ x (t) -γ y (t)| 2 = |γ x (0) -γ y (0)| 2 + t 0 d ds |γ x (t) -γ y (t)| 2 ds we find |x -y| 2 ≤ |γ x (0) -γ y (0)| 2 + 2 K 1 t 0 |γ x (s) -γ y (s)| 2 ds + 1 2 K 1 t 0 | γx -γy | 2 (s)ds.
Combing the above line with (25), we obtain (24) for λ = min L 1 , K

Bounds on u and ∇u. The first thing to notice is a bound for (γ x , γx ). Indeed, the coercivity for [START_REF] Raoul | Etude qualitative et numérique d'équations aux dérivées partielles issues des sciences de la nature[END_REF], as well as the fact that x → ∇R(t, x) grows linearly with a constant only depending on t, implies the existence of a locally bounded constant K(t) such that

|(γ x (s), γx (s))| ≤ K(t)(1 + |x|), for all s ∈ [0, t]. (26)
This implies, modifying the constant K(t) if necessary, thanks to [START_REF] Perthame | Dirac concentrations in Lotka-Volterra parabolic PDEs[END_REF]:

|u(t, x)| ≤ K(t)(1 + |x| 2 ).
Moreover, from (23), we have

|∇u(t, x)| ≤ K(t)(1 + |x|). (27)
Semi-convexity. Because u is three times differentiable in x, then ∂ t u is locally W 2,∞ in x and equation ( 19) may be differentiated twice with respect to x. So, let e be any unit vector, we have

∂ t (∂ ee u) = 2|∇∂ e u| 2 + 2∇u.∇(∂ ee u) + ∂ ee R(t, x) ≥ 2|∂ ee u| 2 + 2∇u.∇(∂ ee u) + ∂ ee R(t, x).
Because of (27), the curve t → γ x (t) becomes a characteristic curve for the equation

∂ t v = 2v 2 + 2∇u(t, x).∇v(t, x) + ∂ ee R(t, x).
Moreover, along this characteristics, we find

d dt v(t, γ x (t)) = 2v(t, γ x (t)) 2 + ∂ ee R(t, γ x (t)).
We deduce, thanks to ( 12) and ( 17), that

∂ ee u ≥ -max(2 L 1 , K 1 ).
Bounds on D 3 u. As for the third derivative, we set v(t, x) = D 3 u(t, x), the equation for v -that we may obtain using differential quotients -is 3 Uniqueness: reduction to a differential system and properties of the function I

∂ t v -2∇v.∇u = S(
The idea is to change the constrained problem (1) by the following slightly nonstandard differential system:

     R (x(t), I(t)) = 0, for t ∈ R + , ẋ(t) = -D 2 u t, x(t) -1 ∇R x(t), I(t) , for t ∈ R + , ∂ t u = |∇u| 2 + R(x, I), in R + × R d , (28)
with initial conditions

I(0) = I 0 , u(0, •) = u 0 (•), x(0) = x 0 , such that max x u 0 (x) = u 0 (x 0 ) = 0 and R(x 0 , I 0 ) = 0. ( 29
)
Note that (28) is really a differential system because the assumptions on R imply that I(t) can implicitely be expressed in terms of x(t). And it is slightly nonstandard because x solves an ODE whose nonlinearity depends on u. The precise statement is the following Proof. Let (u, I) be a solution of (1) with initial datum (u 0 , I 0 ), the function I being continuous, and u a solution of the Hamilton-Jacobi equation in the sense of (20). Theorem 1.1 is applicable, and yields a solution u(t, x) which has at least three locally bounded spatial derivatives, locally uniformaly in time. Moreover, the D 2 u is bounded uniformly in time and in x, and finally the function u(t, .) is strictly concave. This allows a lot.

• There is, at each time, a unique x(t) maximising u(t, .) over R d . Thus the trivial identity ∇u(t, x(t)) = 0 (30) can (use differential quotients) be differentiated with respect to t, to yield that (i) x(t) is locally W 1,∞ and (ii) the (a priori less trivial) identity

∂ t (∇u)(t, x(t)) + D 2 u(t, x(t)). dx dt (t) = 0. ( 31
)
• The function u(t, x) has enough regularity so that we may take the gradient of (1) with respect to x, and evaluate the result at x = x(t). Because of (30) we have D 2 u(t, x(t)).∇u(t, x(t)) = 0 and, because of (31), we have

-D 2 u(t, x(t)). dx dt (t) = ∇R(x(t), I(t)). (32)
• The last item to take into account is the constraint u(t, x(t)) = 0, which we may (still with the use of differential quotients) differentiate with respect to time, in order to yield

∂ t u(t, x(t)) + ∇u(t, x(t)). dx dt (t) = 0, thus entailing ∂ t u(t, x(t)) = 0. This yields, from (1), R(x(t), I(t)) = 0. ( 33
)
Gathering (33), ( 1) and (32) shows that the constrained problem implies (28).

We also prove that regularity plus (28) easily implies [START_REF] Barles | Solutions de viscosité des quations de Hamilton-Jacobi[END_REF]. Let (u, I) solve (28). The first line of (1) derives immediately. To prove the second line, note that, thanks to Theorem 1.1 and the third line of (28) we deduce that u is strictly concave. It has hence, for all t ∈ R + , a unique strict maximum point, in the variable x, that we denote y(t). Following similar arguments as above, we obtain that ẏ(t) = (-D 2 u (t, y(t)) -1 ∇R(y(t), I(t)), for t ∈ R + ,

y(0) = x 0 .
Comparing this with the second line of ( 28) and ( 29) we obtain that x(t) = y(t), for all t ∈ R + . Finally, evaluating the third line of (28) at x(t) and using the first line of (28), we obtain that ∂ t u(t, x(t)) = 0. This equality together with ∇u(t, x(t)) = 0 and (29) implies that max x u(t, x) = u(t, x(t)) = 0. and the proof of Theorem 3.1 is complete.

Remark. What we have done here is nothing else than the derivation of the equation for x carried out in [START_REF] Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF]. In particular, the equilibrium property R(x, I) = 0 would hold in a more general setting than here. The new point here is that we establish, in a mathematically rigorous fashion, the equivalence between the initial problem (1) and the coupled system (28), and this equivalence holds because of all the differentiations with respect to x and t that we are allowed to make.

The proof of Theorem 1.2

We fix T > 0. To prove that (1) has a unique solution (u, I) in [0, T ] × R d , it is enough to prove that there exists a unique solution to (28)-( 29).

We prove this using the Banach fixed point Theorem in a small interval and then iterate.

The mapping Φ and its domain

A = x(•) ∈ C [0, δ]; B x 0 , r δ x(0) = x 0 ,
where B(z, r) is the ball of radius r centered at z, and δ is a positive constant such that

δ < min(µ, c(T )), with µ = min(2L 1 , K 1 ) C M d(Ω 0 , Ω), (34)
where d(A, B) is the distance between the sets A and B. The constant r δ is given by

r δ = C M δ min(2L 1 , K 1 ) , the constant C M is chosen such that |∇R(x, I)| ≤ C M , in Ω × [0, I M ], (35)
and c(T ), a constant depending only on T , will be chosen later. Note that, by the choice of δ and since x 0 ∈ Ω 0 , we obtain that B (x 0 , r δ ) ⊂ Ω.

Our theorem will be proved by the introduction of a mapping Φ : A → A, Φ(x) = y.

We will prove the following theorem.

Theorem 4.1 The mapping Φ is a strict contraction from A into itself.

To define Φ, we first need to introduce some other mappings. Let x(•) ∈ A. We define

I : A → C ([0, δ]; [0, I M ]) such that R(x(t), I[x](t)) = 0.
From ( 5), ( 8) and since x(t) ∈ Ω for all t ∈ [0, δ], it follows that I can be defined in a unique way. Next, we define the following domain

B = v ∈ L ∞ [0, δ]; W 3,∞ loc (R d ) ∩ W 1,∞ [0, δ]; L ∞ loc (R d ) | -max(2 L 1 , K 1 ) ≤ D 2 v ≤ -min(2L 1 , K 1 ), D 3 v L ∞ ([0,δ]×R d ) ≤ L 4 (T ) ,
and the following mapping

V : C ([0, δ]; [0, I M ]) → B V (I) = v,
where v solves

∂ t v = |∇v| 2 + R(x, I), in [0, δ] × R d , v(0, x) = u 0 (x), in R d . (36)
It follows from Theorem 1.1 that the above mapping is well-defined. Finally, we introduce a last mapping:

F : C ([0, δ]; [0, I M ]) × A × B → A, F (I, x, v) = y,
where y ∈ A solves

ẏ(t) = -D 2 v(t, x(t) -1 ∇R(x(t), I(t)), in [0, δ], y(0) = x 0 .
To prove that F is well-defined, we must verify that y(t) remains in B (x 0 , r δ ). We note that, since v ∈ B, we have

0 < 1 max(2 L 1 , K 1 ) ≤ -D 2 v(t, x(t)) -1 ≤ 1 min(2L 1 , K 1 ) .
We deduce, thanks to (35), that y(t) ∈ B (x 0 , r δ ) .

We are now ready to define mapping Φ:

Φ : A → A, Φ(x) = F I(x), x, V (I(x)) .
It follows from the above arguments that Φ is well-defined.

Proof of Theorem 1.2

In this section we will explain why uniqueness to (28) holds. One main technical lemma will be stated, its proof will be postponed to a special section. Let us prove that Φ is a contraction for c(T ) (and hence δ) small enough. To this end, we first prove that I is Lipschitz:

|I(x 1 ) -I(x 2 )| ≤ C x 1 -x 2 L ∞ ([0,δ]) , for all x 1 , x 2 ∈ C [0, δ]; B x 0 , r δ . (37)
We have indeed

R(x 1 , I(x 1 )) = R(x 2 , I(x 2 )) = 0. It follows that R(x 1 , I(x 1 )) -R(x 2 , I(x 1 )) = R(x 2 , I(x 2 )) -R(x 2 , I(x 1 )),
and thus

∇R (cx 1 + (1 -c)x 2 , I(x 1 )) • (x 2 -x 1 ) = ∂R ∂I (x 2 , J)(I(x 2 ) -I(x 1 )),
with c ∈ (0, 1) and J ∈ (I(x 1 ), I(x 2 )). Finally, using ( 8), ( 35) and the fact that x 1 (t), x 2 (t) ∈ Ω for all t ∈ [0, δ] we obtain (37

) with C = C M K 2 . Next, we have that V : C ([0, δ]; [0, I M ]) → B is also Lipschitz. Lemma 4.2 Let I 1 , I 2 ∈ C ([0, δ]; [0, I M ]). Then V (I 1 ) -V (I 2 ) W 2,∞ ([0,δ]×R d ) ≤ C I 1 -I 2 L ∞ ([0,δ]) δ. (38)
This is a nontrivial lemma, whose proof will be given in the next section.

From ( 7) and ( 9), and since x 1 , x 2 ∈ A and V 1 , V 2 ∈ B, we deduce that F : C ([0, δ]; [0, I M ])×A×B → A is Lipschitz with respect to all the variables with Lipschitz constant Cδ:

F (I 1 , x 1 , V 1 ) -F (I 2 , x 2 , V 2 ) L ∞ ([0,δ]) ≤ Cδ x 1 -x 2 L ∞ ([0,δ]) + I 1 -I 2 L ∞ ([0,δ]) + V 1 -V 2 W 2,∞ ([0,δ]×R d ) . (39)
Finally, we conclude from (37), ( 38) and (39) that Φ : A → A is a Lipschitz mapping with a Lipschitz constant Cδ:

Φ(x 1 ) -Φ(x 2 ) L ∞ ([0,δ]) ≤ Cδ x 1 -x 2 L ∞ ([0,δ]) .
Choosing c(T ) (and hence δ) small enough, we deduce that Φ is a contraction. We deduce from the Banach fixed point Theorem that Φ has a unique fixed point and consequently (28)-( 29) has a unique solution for t ∈ [0, δ].

To prove that (28)-( 29) has a unique solution in [0, T ] we iterate the above procedure K = ⌈ T δ ⌉ times. Let 2 ≤ i ≤ K and (x, I, u) be the unique solution of (28)-(29) for t ∈ [0, (i -1)δ]. Then, at the i-th step, we consider the same mapping Φ but as the initial data we choose

x 0 = x((i -1)δ), u 0 (•) = u((i -1)δ, •), I 0 = I((i -1)δ). ( 40
)
We claim that these initial conditions satisfy

x 0 ∈ Ω 0 , u 0 ∈ W 3,∞ loc (R d ), with -max(2 L 1 , K 1 ) ≤ D 2 u 0 ≤ -min(2L 1 , K 1 ) D 3 u 0 L ∞ (R d ) ≤ L 4 ((i -1)δ), max x u 0 (x) = x 0 , R( x 0 , I 0 ) = 0. ( 41
)
Lemma 4.3 Let (x, I, u) be the unique solution of (28)-( 29) in [0, τ ]. Then, I(t) is increasing with respect to t in [0, τ ] and

x(t) ∈ Ω 0 , u(t, •) ∈ W 3,∞ loc (R d ), with -max(2 L 1 , K 1 ) ≤ D 2 u(t, x) ≤ -min(2L 1 , K 1 ), D 3 u L ∞ ([0,τ ]×R d ) ≤ L 4 (τ ), max x u(t, x) = x(t), R(x(t), I(t)) = 0, , for all t ∈ [0, τ ] and x ∈ R d .
Proof. The regularity estimates on u are immediate from Section 2. The two last claims also follow immediately from (28)-( 29) and the arguments in Section 3. We only prove that I(t) is increasing with respect to t and x(t) ∈ Ω 0 .

Differentiating the first line of (28) with respect to t we obtain

∇R(x(t), I(t)) • ẋ(t) + ∂ ∂I R(x(t), I(t)) İ (t) = 0.
Moreover, multiplying the second line of (28) by ∇R we obtain

∇R(x(t), I(t)) • ẋ(t) = ∇R(x(t), I(t)) -D 2 u(t, x(t) -1 ∇R(x(t), I(t)) ≥ 0.
Combining the above lines we obtain

∂ ∂I R(x(t), I(t)) İ (t) ≤ 0, for t ∈ [0, τ ],
and hence, thanks to (8) we deduce

İ(t) ≥ 0, for t ∈ [0, τ ].
Therefore, I(t) is increasing with respect to t and in particular

I(t) ≥ I 0 , for t ∈ [0, τ ].
Consequently, from the first line of (28) and ( 8), we obtain that

R(x(t), I 0 ) > 0, for t ∈ [0, τ ]. It follows that x(t) ∈ Ω 0 , for t ∈ [0, τ ].
which concludes the proof.

It is then immediate that the initial data given by (40) verify (41). One can verify that the above conditions are the only properties that we have used to prove that Φ is well-defined and a contraction. Therefore, one can apply again the Banach fixed point Theorem and deduce that there exists a unique solution of (28)-(29) for t ∈ [(i -1)δ, iδ].

Proof of Lemma 4.2

(i) We first prove that

V (I 1 ) -V (I 2 ) L ∞ ([0,δ]×R d ) ≤ C I 1 -I 2 L ∞ ([0,δ]) δ. (42)
Let v 1 = V (I 1 ) and v 2 = V (I 2 ), and r = v 2 -v 1 . From (36) we obtain

∂ t r = (∇v 1 + ∇v 2 ) • ∇r + R(x, I 2 ) -R(x, I 1 ), in [0, δ] × R d r(0, x) = 0, for all x ∈ R d . (43)
Note that the above equation has a unique classical solution which can be computed by the method of characteristics. The characteristics verify

γ(t) = -∇v 1 (t, γ) -∇v 2 (t, γ). (44)
For any (t 1 , x 1 ) ∈ [0, δ] × R d , there exists a unique characteristic curve γ which verifies γ(t 1 ) = x 1 . Moreover, this characteristic curve is defined in [0, t 1 ].

The local existence and uniqueness of the characteristic curve γ is derived from the Cauchy-Lipschitz Theorem and the fact that ∇v i is Lipschitz with respect to x, for i = 1, 2. The latter property derives from the fact that 2L

1 ≤ D 2 v i ≤ -2L 1 , for i = 1, 2.
To prove that the characteristic curve γ is defined in [0, t 1 ] we must prove that γ remains bounded in this interval. However, this property follows using (44) and the fact that, since

v 1 , v 2 ∈ B, |∇v i (γ)| ≤ C 1 |γ| + C 2 , for i = 1, 2.
We are now ready to prove (38). To this end, we multiply (44) by ∇r(t, γ) to obtain

∇r(t, γ) • γ(t) = -(∇v 1 (t, γ) + ∇v 2 (t, γ)) • ∇r(t, γ).
Combining this with (43) we deduce that ∇r(t, γ(t)) • γ(t) + ∂ t r(t, γ(t)) = R(x, I 2 ) -R(x, I 1 ).

We integrate this between 0 and t 1 to find r(t 1 , γ(t 1 )) = Since D 2 v j , for j = 1, 2, is bounded, and using [START_REF] Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF] (iii) Finally we prove that

D 2 r L ∞ ([0,δ]×R d ) ≤ C I 1 -I 2 L ∞ ([0,δ]) δ. ( 46
)
Note that at every point (t, x) ∈ R + × R d , we can write mixed derivatives of the form r ξη in terms of pure derivatives:

r ξη = 1 2 (∂ 2 ξ+η,ξ+η r -∂ 2 ξξ r -∂ 2 ηη r). (47)
This implies the existence of ξ on the unit sphere of R d such that

D 2 r L ∞ ([0,δ]×R d ) ≤ 3 2 r ξξ L ∞ ([0,δ]×R d ) .
We differentiate (43) twice in the direction of ξ and obtain ∂ t r ξξ = (∇v 1,ξξ + ∇v 2,ξξ ) • ∇r + 2 (∇v 1,ξ + ∇v 2,ξ ) • ∇r ξ + (∇v 1 + ∇v 2 ) • ∇r ξξ + R ξξ (x, I 2 ) -R ξξ (x, I 1 ).

Using the above arguments, the fact that D 2 v i and D 3 v i are bounded, ( 9) and (45) we deduce that

  t, x, v) := 6v.D 2 u + D 3 R, where ∇v denotes the column of tensors (∂ 1 v, ..., ∂ d v) and v.D 2 u denotes the column of matrices (∂ 1 D 2 u.D 2 u, ..., ∂ d D 2 u.D 2 u). This is a linear equation with (thanks to the bound on D 2 u) bounded coefficients. Thus, local boundedness of v(t, .) L ∞ (R d ) holds, and this concludes the proof of Theorem 1.1.

Theorem 3 . 1

 31 Solving the constrained problem (1) is equivalent to solving the initial value ODE-PDE problem (28)-(29).

First, we

  define Ω = {x | R(x, 0) > 0} , Ω 0 = {x | R(x, I 0 ) ≥ 0} , and

t 1 0R

 1 (x, I 2 (τ )) -R(x, I 1 (τ ))dτ. It follows, using[START_REF] Diekmann | Beginner's guide to adaptive dynamics[END_REF], that|r(t 1 , γ(t 1 ))| ≤ K 2 |I 1 -I 2 |t 1 ,and hence (42).(ii) Next, we prove that∇r L ∞ ([0,δ]×R d ) ≤ C I 1 -I 2 L ∞ ([0,δ]) δ. (45)We differentiate (43) in the direction e i , for i = 1, • • • , d, to obtain∂ t r (i) = ∇v 1,(i) + ∇v 2,(i) • ∇r + (∇v 1 + ∇v 2 ) • ∇r (i) + R (i) (x, I 2 ) -R (i) (x, I 1 ),with the notation f (i) = ∇f • e i . We multiply by r (i) , sum over i and divide by |∇r| to obtain∂ t |∇r| ≤ (∇v 1 + ∇v 2 ) • ∇|∇r| + n i=1 ∇v1,(i) + ∇v 2,(i) |r (i) | + |∇R(x, I 2 ) -∇R(x, I 1 )|.

∂

  t |r ξξ | ≤ C I 1 -I 2 L ∞ ([0,δ]) δ + C r ξξ L ∞ ([0,δ]×R d ) + (∇v 1 + ∇v 2 ) • ∇|r ξξ | + K 3 |I 2 -I 1 |.Next we use the characteristic curves as previously. We multiply (44) by ∇|r ξξ (t, γ(t))| and obtain∇|r ξξ (t, γ(t))| • γ(t) = -(∇v 1 (t, γ) + ∇v 2 (t, γ)) • ∇|r ξξ (t, γ)|.We combine the above equations to obtaind dt |r ξξ (t, γ(t))| ≤ C I 1 -I 2 L ∞ ([0,δ]) + C r ξξ L ∞ ([0,δ]×R d ) .We conclude that|r ξξ (t, x)| ≤ Ct( I 2 -I 1 L ∞ ([0,δ]) + r ξξ L ∞ ([0,δ]×R d ) ).Restricting c(T ) (and hence δ) if necessary we getr ξξ L ∞ ([0,δ]×R d ) ≤ Cδ 1 -Cδ I 2 -I 1 L ∞ ([0,δ]) ,hence (46).

  we deduce that∂ t |∇r| ≤ (∇v 1 + ∇v 2 ) • ∇|∇r| + C|∇r| + K 3 |I 2 -I 1 |.The characteristic curves corresponding to the above equation verify again (44). We multiply (44) by ∇|∇r|(t, γ(t)) to obtain∇|∇r|(t, γ) • γ(t) = -(∇v 1 (t, γ) + ∇v 2 (t, γ)) • ∇|∇r|(t, γ). C|∇r|(t, γ(t)) + K 3 I 2 -I 1 L ∞ ([0,δ]) . follows that |∇r(t, x)| ≤ e Ct -1 K 3 C I 2 -I 1 L ∞ ([0,δ]) .Hence (45), modifying the constant C if necessary.

	Combining the above equations we obtain
	d dt |∇r|(t, γ(t)) ≤ It

Acknowledgements

The authors thank G. Barles and B. Perthame for interesting comments on an earlier version of this paper, which had a great influence on the present form of this work. They also thank M. Bardi for pointing out the possible tractability of the quadratic case. S. Mirrahimi was partially funded by the ANR projects KIBORD ANR-13-BS01-0004 and MODEVOL ANR-13-JS01-0009. J.-M. Roquejoffre was supported by the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 321186 -ReaDi -"Reaction-Diffusion Equations, Propagation and Modelling" held by Henri Berestycki, as well as the ANR project NONLOCAL ANR-14-CE25-0013. Both authors thank the Labex CIMI for a PDE-probability winter quarter in Toulouse, which provided a stimulating scientific environment to this project.

5 An example with quadratic u 0 and R Let us study the (instructive) example of a quadratic equation [START_REF] Barles | Solutions de viscosité des quations de Hamilton-Jacobi[END_REF]. In other words we choose

where A 0 and A 1 are positive definite matrices and b ∈ R d . The Euler-Lagrange equation for the dynamic programming principle writes

In order to translate the equivalent system (28), we need ∇R(I, x) and D 2 u(t, x). The first quantity is easily obtained:

We then solve (48), differentiate with respect to x and denote by Γ the differential of γ -the solution of (48) with respect to x, to obtain (after an elementary but tedious computation)

Notice that the result was expected because γ is linear function of x, and that B(s) is invertible for 0 ≤ s ≤ t, and the norms of B and B -1 are bounded uniformly in s and t. And so we have

and the matrix C(t) is bounded away from 0 and +∞, uniformly with respect to t. This, by the way, is not easily seen on the formula (50); the proof of Theorem 1.1 is still what one should use here. The equation for x(t) and I(t) is thus straightforward:

Thus, existence and uniqueness of (x, I) is straightforward. For large t, examination of (50) and (49) yields

This is consistent with the known behaviour of I(t) and x(t), as well as a closed form of the competition increase I(t) -I 0 .