Stefano Fortunati

Lorenzo Ortega

On the efficiency of misspecified Gaussian inference in nonlinear regression: application to time-delay and Doppler estimation

Keywords: Nonlinear regression, Misspecified Cramér-Rao bound, Mismatched Maximum Likelihood estimator, time-delay and Doppler estimation, band-limited signals

Nonlinear regression play a crucial role in various engineering applications. For the sake of mathematical tractability and ease of implementation, most of the existing inference procedures are derived under the assumption of independent and identically distributed (i.i.d.) Gaussian-distributed data. However, real-world situations often deviate from this assumption, with the true data generating process being a correlated, heavytailed and non-Gaussian one. The paper aims at providing the Misspecified Cramér-Rao Bound (MCRB) on the Mean Squared Error (MSE) of any unbiased (or at least consistent) estimator of the parameters of a nonlinear regression model derived under the i.i.d. Gaussian assumption in the place of the actual correlated, non-Gaussian data generating process. As a special case, the MCRB for an uncorrelated, i.i.d. Complex Elliptically Symmetric (CES) data generating process under Gaussian assumption is also provided. Consistency and asymptotic normality of the related Mismatched Maximum Likelihood Estimator (MMLE) will be discussed along with its connection with the Nonlinear Least Square Estimator (NLLSE) inherent to the nonlinear regression model. Finally, the derived theoretical findings will be applied in the well-known problem of time-delay and Doppler estimation for GNSS.

I. INTRODUCTION

N ONINEAR regression are one of the most-used statistical models in Signal Processing (SP) and related engineering applications. In a regression model, an observation vector

C 𝑁 ∋ x = f ( θ) + n, (1)
is characterized by i) a vector of unknown deterministic parameters θ, ii) a linear or nonlinear (continuous and differentiable) known function f, parameterized by θ and iii) an additive noise vector n. The function f generally specifies the measurement process while θ collects the quantities that need to be estimated. Regression models can be found in array processing, image processing, biomedical data analysis and even in climatic studies, just to name a few. While the definition of f, θ and of the measurement noise n depends on the particular application at hands, the inference procedures used to estimate Stefano Fortunati is with Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des signaux et systèmes, 91190, Gif-sur-Yvette & DR2I-IPSA, 94200, Ivry-sur-Seine, France.

Lorenzo Ortega is with TéSA, 7, Boulevard de la Gare, 31500, Toulouse, France & DR2I-IPSA, 40, Boulevard de la Marquette, 31000, Toulouse, France.

the parameter of interest usually share a common (although unrealistic) assumption: the entries of the noise vector n are sampled from an i.i.d. white Gaussian random process. This assumption is made to make the estimation algorithm mathematically tractable and easy to implement. In fact, it is well known that, under the i.i.d. Gaussian assumption, the optimal estimator is the (non) linear least square estimator (NLLSE). However, everyday engineering practice shows that this assumption is too simplistic since the noise process can be correlated and even non-Gaussian. The central question that we aim at answering in this paper is: how accurate can an i.i.d., Gaussian-based inference procedure be when the regression model is characterized by a correlated, generally non-Gaussian noise? In order to answer to this question, we will rely on the misspecification theory developed in [START_REF] Huber | The behavior of maximum likelihood estimates under nonstandard conditions[END_REF]- [START_REF] Vuong | Cramér-Rao bounds for misspecified models[END_REF] and recently rediscovered in [START_REF] Fortunati | Performance bounds for parameter estimation under misspecified models: Fundamental findings and applications[END_REF]- [START_REF] Mennad | Slepian-bangs-type formulas and the related misspecified cramér-rao bounds for complex elliptically symmetric distributions[END_REF] and the reference therein. Specifically, we will show that it is possible to derive the so-called Misspecified Cramér-Rao Bound (MCRB) on the estimation of θ when the assumed model is the "classical" i.i.d. Gaussian model while the true data model is a dependent and non-Gaussian one.

This paper is organized into seven distinct sections. In Section II, we present both the true and assumed nonlinear regression models. Section III introduce the calculation of the pseudo-true parameter vector for the misspecified Gaussian nonlinear regression model. Section IV derives the MCRB under quite general condition of the correlation structure of the true data generating process, while, in Sec. V, we specialize this general results to a case of practical interest in which the true data model is an i.d.d. Complex Elliptical Symmetric (CES) model with unspecified density generator. Sec. VI is dedicated to the investigation of the asymptotic properties of the NLLE under the above mentioned misspecified scenario and to its relation with the MMLE. Sec. VII provides an example of possible application of the theoretical results to the time-delay and Doppler estimation under the abovementioned misspecified scenario for GNSS applications. Our conclusion are collected in Sec. IX. Notation: Throughout this paper, italics indicates scalar quantities (𝑎), lower case and upper case boldface indicate column vectors (a) and matrices (A), respectively. Each entry of a matrix A is indicated as 𝑎 𝑖 𝑗 ≜ [A] 𝑖, 𝑗 . I 𝑁 defines the 𝑁 × 𝑁 identity matrix. The superscripts * , ⊤ and 𝐻 indicate the complex conjugation, the transpose and the Hermitian operators respectively, then A H = (A * ) ⊤ . The Euclidean norm of a vector a is indicated as ||a||.

II. NONLINEAR REGRESSION WHIT DEPENDENT

OBSERVATIONS

As discussed in the Introduction, the nonlinear regression is one of the most used statistical model in signal processing (SP) and statistics. The aim of this section is then to introduce firstly the model in its generality (i.e. the true signal model) and secondly to present its Gaussian-based, i.i.d. simplified version as it is generally assumed by SP practitioners for inference purposes.

A. True signal model

Let {𝑥 𝑘 ∈ C} +∞ 𝑘=-∞ be a sequence of scalar, complex-valued, observations characterized by the following data generating process:

𝑥 𝑘 = 𝑓 𝑘 ( θ) + 𝑛 𝑘 , -∞ < 𝑘 < +∞, (2)
where θ ∈ Θ ⊂ R 𝑝 indicates the real-valued 1 , true parameter vector and Θ is a compact subset of R 𝑝 . The functions 𝑓 𝑘 : Θ → C, -∞ < 𝑘 < +∞ are known continuous and differentiable functions on Θ. In practical applications, the sequence (2) will be observed from a finite integer

𝑁 1 ∈ Z to a finite integer 𝑁 2 ∈ Z, such that -∞ < 𝑁 1 < 𝑁 2 < +∞.
Consequently, by defining 𝑁 = |𝑁 2 -𝑁 1 + 1|, the sequence (2) can be written in a vectorial form as in [START_REF] Huber | The behavior of maximum likelihood estimates under nonstandard conditions[END_REF] where n ∈ C 𝑁 is a zero-mean complex-valued noise vector whose 𝑁 entries are assumed to be sampled from a Wide Sense Stationary (WSS) discrete random process {𝑛 𝑘 : ∀𝑘 } characterizing the measurement noise 𝑛 𝑘 in (2).

Let us now have a closer look at the statistical characterization of {𝑛 𝑘 : ∀𝑘 }. As a zero-mean WSS process, {𝑛 𝑘 : ∀𝑘 } is fully characterized by (see e.g. [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]Sec. 15.5], [START_REF] Stoica | Spectral Analysis of Signals[END_REF]Sec. 1.3]): i) its autocorrelation function 𝑟 𝑛 [𝑘 + 𝑗, 𝑘] = 𝑟 𝑛 [𝑘 + 𝑗 -𝑘] = 𝑟 𝑛 [ 𝑗] and ii) the joint probability density function (pdf) of the 𝑁 samples C 𝑁 ∋ n ∼ 𝑝 n , for any values of 𝑁. For further reference, we indicate the marginal pdf of each sample as 𝑛 𝑘 ∼ 𝑝 𝑛 𝑘 . We make the following extremely general (non-Gaussian, non-i.i.d.) assumption:

Assumption 1: Let {𝑛 𝑘 : ∀𝑘 } be a zero-mean, WSS discrete and circular complex-valued process [START_REF] Picinbono | On circularity[END_REF] such that the joint pdf of 𝑁 samples follows an unspecified pdf n ∼ 𝑝 n , ∀𝑁 admitting finite first and second order moments. Then, we assume that its autocorrelation function exists and satisfies

|𝑟 𝑛 [ 𝑗] | ≜ |𝐸 𝑝 n [𝑛 * 𝑘+ 𝑗 𝑛 𝑘 ] | = 𝑂 (| 𝑗 | -𝛾 ), 𝑚 ∈ Z, 𝛾 > 𝜚/( 𝜚 -1), 𝜚 > 1. 2 Note that the circularity of {𝑛 𝑘 : ∀𝑘 } implies that 𝐸 𝑝 n [𝑛 𝑘+ 𝑗 𝑛 𝑘 ] = 0, ∀𝑘, 𝑗.
It is worth noticing here that, as a direct consequence of this assumption: 1 We decided to work with real-valued parameters for two reasons. Firstly, in practical applications, the parameters of interest are real-valued physical quantities as e.g. time-delay and Doppler. Secondly, this choice allows us to avoid the technicalities related to the Wirtinger calculus [START_REF] Kreutz-Delgado | The complex gradient operator and the CRcalculus[END_REF] that may obscure the more important statistical concepts. It is worth stressing that this choice will not limit the generality of the derived results since any complex-valued vector can be recast in term of a real-valued vector by means of the standard isomorphism between C 𝑝 and R 2 𝑝 , i.e.

C 𝑝 ∋ 𝜃 ⇋ θ ≜ (Re( 𝜃 ) ⊤ , Im( 𝜃 ) ⊤ ) ⊤
2 Given a real-valued function 𝑓 ( 𝑥 ) and a positive real-valued function

𝑔 ( 𝑥 ), 𝑓 ( 𝑥 ) = 𝑂 (𝑔 ( 𝑥 ) ) if
= 𝑟 𝑛 [(𝑘 + 𝑗) -𝑘] = 𝑟 * 𝑛 [𝑘 -(𝑘 + 𝑗)] = 𝑟 𝑛 [ 𝑗] = 𝑟 * 𝑛 [-𝑗], [𝚺] 𝑘,𝑘 = 𝑟 𝑛 [0] = σ2 𝑛 , ∀𝑘, (3)
where σ2 𝑛 is the true and generally unknown noise power. We would like to stress that Assumption 1 is extremely general and allows for a wide range of realistic noise models [START_REF] Fortunati | Massive mimo radar for target detection[END_REF]. To be convinced of this, we can note that any (Gaussian and non-Gaussian) stable second-order stationary (SOS) ARMA, of any finite orders, satisfies Assumption 1, since the autocorrelation function of any stable SOS ARMA decays exponentially. It is well know that, by appropriately choosing the orders of the Autoregressive and of the Moving Average parts, an ARMA process can approximate the (continuous) power spectral density (PSD) of any complex discrete random processes [START_REF] Stoica | Spectral Analysis of Signals[END_REF]Ch. 3]. Moreover, a non-Gaussian ARMA can characterize the heavy-tailed behavior of realistic noise models. Another popular noise model of practical interest satisfying Assumption 1 is the Compound-Gaussian (CG) (or spherically invariant random vector (SIRV)) model [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]. In fact, any SIRV n ∈ C 𝑁 can be represented as [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]Def. 3] n = 𝑑 √ 𝜏m for some real-valued positive random variable 𝜏, such that 𝐸 [𝜏] = 1, called texture, independent of the zeromean, 𝑁-dimensional, circular, complex Gaussian random vector, called speckle, m ∼ CN (0, 𝚺), where 𝚺 is the covariance matrix given in (3).

To conclude this section, we note that the pdf of the data vector in eq. ( 1) can be expressed as function of the unspecified noise pdf 𝑝 n as:

x ∼ 𝑝 ε ≜ 𝑝 ε (x; σ2 𝑛 , θ) = 𝑝 n (𝑥 -f ( θ); σ2 𝑛 ), (5)
where

ε = ( σ2 𝑛 , θ⊤ ) ⊤ ∈ Γ ⊂ R + × R 𝑝 (6)
is the complete vector of the true parameters, where θ is the vector of the parameter of interest and R + ∋ σ2 𝑛 > 0 a nuisance parameter, i.e. a term whose estimation is not strictly required but the lack of its knowledge may have an impact on the estimation performance of θ.

B. Misspecified Gaussian, i.i.d. signal model

To do inference on the parameter vector θ, and specifically to estimate it, a common procedures among SP practitioners is to assume a simplified model describing the statistical behaviour of the observations in the place of the true data generating process in [START_REF] White | Maximum likelihood estimation of misspecified models[END_REF]. This model misspecification is dictated by two main reasons [START_REF] Fortunati | Performance bounds for parameter estimation under misspecified models: Fundamental findings and applications[END_REF]. The first one is that the autocorrelation structure, as well as the pdf 𝑝 n is generally not a-priori known and not easy to obtain from physical considerations on the random experiment at hand. Secondly, one could prefer a simplified model in order to derive estimation algorithms that are easy to implement and fast to compute.

One of the most popular simplifying assumption is to consider the noise process {𝑛 𝑘 : ∀𝑘 } as a zero-mean, White Gaussian WSS random process. This implies that its autocorrelation function can be expressed as 𝑟 𝑛 [ 𝑗] = σ2

𝑛 𝛿[ 𝑗], where 𝛿[ 𝑗] is the Kronecker delta sequence. As a consequence, the noise vector n ∈ C 𝑁 is distributed as a centered complex normal random vector with diagonal covariance matrix, i.e. n ∼ CN (0, 𝜎 2 𝑛 𝑰 𝑁 ). This simplifying assumption leads to the following misspecified statistical model for the data vector x ∈ C 𝑁 in (1):

F 𝝐 ≜ 𝑓 𝝐 | 𝑓 𝝐 (x; 𝝐) = CN (f (𝜽), 𝜎 2 𝑛 𝑰 𝑁 ), 𝝐 ∈ Γ , (7)
that is, each pdf belonging to F 𝝐 can be expressed as:

𝑓 𝝐 (x; 𝜎 2 𝑛 , 𝜽) = (𝜋𝜎 2 𝑛 ) -𝑁 𝑒 -||x-f (𝜽) || 2 𝜎 2 𝑛 . ( 8
)
The crucial question that we are going to answer in the next section is: is it possible to derive a lower bound to the Mean Squared Error (MSE) of any unbiased (or at least consistent) estimation procedure of θ, derived under the Gaussian, i.i.d., misspecified model F 𝝐 in [START_REF] Fortunati | The constrained Misspecified Cramér-Rao bound[END_REF] in the presence of dependent observations satisfying Assumption 1?

To answer to this question, we evaluate the MCRB [START_REF] Fortunati | Performance bounds for parameter estimation under misspecified models: Fundamental findings and applications[END_REF]- [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecification[END_REF] on the estimation of θ when the assumed model is F 𝝐 while the true data generating process is the (dependent, non-Gaussian) one introduced in (2). To this end, we start by evaluating the pseudo-true parameter vector 𝝐 0 ∈ Γ, i.e. the vector in Γ that minimizes the Kullback-Leibler Divergence (KLD) [4, A 1] [5, Sec. 4.4.1] between the true (and unknown) pdf x ∼ 𝑝 ε and any element 𝑓 𝝐 ∈ F 𝝐 of the assumed misspecified model in [START_REF] Fortunati | The constrained Misspecified Cramér-Rao bound[END_REF]. The vector 𝝐 0 ∈ Γ can be seen as a sort of "minimum divergence projector" of the true pdf onto the misspecified model F ε and then it characterizes the pdf 𝑓 𝝐 0 ∈ F 𝝐 closest, in the KLD sense, to the true pdf 𝑝 ε .

III. THE PSEUDO-TRUE PARAMETER VECTOR

As anticipated in the previous section, the pseudo-true parameter vector 𝝐 0 is the element in the parameter space Γ that minimizes the KLD between the true data pdf x ∼ 𝑝 ε and any (possibly) misspecified pdf 𝑓 𝝐 0 ∈ F 𝝐 [START_REF] Vuong | Cramér-Rao bounds for misspecified models[END_REF], [4, A1] and [START_REF] Fortunati | Chapter 4 -Parameter bounds under misspecified models for adaptive radar detection[END_REF]Sec. 4.4.1]:

𝐷 ( 𝑝 ε || 𝑓 𝝐 ) = 𝐸 𝑝 ε ln 𝑝 ε (x; ε) 𝑓 𝝐 (x; 𝝐) x ∼ 𝑝 ε , 𝑓 𝝐 ∈ F ε (9) 𝐸 𝑝 ε [•]
is the expectation with respect to (w.r.t.) the true model's pdf. Consequently:

𝝐 0 = arg min 𝝐 ∈Γ {𝐷 ( 𝑝 ε || 𝑓 𝝐 )} = arg min 𝝐 ∈Γ 𝐸 𝑝 ε [-ln 𝑓 𝝐 (x; 𝝐)] . (10
) From ( 8), it follows directly that:

𝐸 𝑝 ε [-ln 𝑓 𝝐 ] = 𝑁 ln(𝜋) + 𝑁 ln(𝜎 2 𝑛 ) + 𝐸 𝑝 ε ||x -f (𝜽)|| 2 𝜎 2 𝑛 . (11)
By substituting [START_REF] Stoica | Spectral Analysis of Signals[END_REF] in [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF], we have:

𝝐 0 = arg min 𝝐 ∈Γ 𝐸 𝑝 ε -ln 𝑓 𝝐 (x; 𝜎 2 𝑛 , 𝜽) = arg min 𝝐 ∈Γ 𝐸 𝑝 ε 1 𝜎 2 𝑛 ∥x -f ( θ) ∥ 2 + 𝑁 ln(𝜎 2 𝑛 ) (12)
Let us start by minimizing w.r.t. to 𝜽:

𝜽 0 = arg min 𝜽 𝐸 𝑝 ε [-ln 𝑓 𝝐 (x; 𝝐)] (13a) = arg min 𝜽 𝐸 𝑝 ε ∥x -f ( θ) ∥ 2 (13b) = arg min 𝜽 𝐸 𝑝 ε tr (x -f ( θ)) (x -f ( θ)) 𝐻 (13c) = arg min 𝜽 tr 𝚺 + f ( θ)f ( θ) 𝐻 -× (13d) × -f ( θ)f (𝜽) 𝐻 -f (𝜽)f ( θ) 𝐻 + f (𝜽)f (𝜽) 𝐻 (13e) = arg min 𝜽 ||f (𝜽) -f ( θ)|| 2 ⇒ 𝜽 0 = θ. (13f)
Remarkably, this result tells us that the pseudo-true parameter vector of interest 𝜽 0 is equal to the one θ.

Let us now minimize w.r.t. to the variance 𝜎 2 𝑛 . By using the result obtained in (13), we have:

𝜎 2 0 = arg min 𝜎 2 𝑛 𝐸 𝑝 ε -ln 𝑓 𝝐 (x; 𝜎 2 𝑛 , θ) (14)
⇒ 𝐸 𝑝 ε 𝜕 𝜕𝜎 2 𝑛 ln 𝑓 𝝐 (x; 𝜎 2 𝑛 , θ) 𝜎 2 𝑛 =𝜎 2 0 = 0 (15)
From direct calculation, we have:

𝐸 𝑝 ε 𝜕 𝜕𝜎 2 𝑛 ln 𝑓 𝝐 (x; 𝜎 2 𝑛 , θ) 𝜎 2 𝑛 =𝜎 2 0 (16a) = 𝐸 𝑝 ε - 𝑁 𝜎 2 𝑛 + 1 𝜎 4 𝑛 ∥x -f ( θ) ∥ 2 𝜎 2 𝑛 =𝜎 2 0 (16b) = 𝐸 𝑝 ε - 𝑁 𝜎 2 𝑛 + tr(nn 𝐻 ) 𝜎 4 𝑛 = - 𝑁 𝜎 2 0 + tr (𝚺) 𝜎 4 0 (16c) = - 𝑁 𝜎 2 0 + 𝑁𝑟 𝑛 [0] 𝜎 4 0 = - 𝑁 𝜎 2 0 + 𝑁 σ2 𝑛 𝜎 4 0 = 0 ⇒ 𝜎 2 0 = σ2 𝑛 (16d)
Again, eq. ( 16) tells us that the pseudo-true nuisance parameter 𝜎 2 0 equates the true one σ2 𝑛 . By collecting the results from eqs. ( 13) and ( 16), we have that the pseudo-true parameter vector equates the true one

𝝐 0 = ε ≜ ( σ2 𝑛 , θ⊤ ) ⊤ , (17)
under extremely mild assumptions, i.e. for any noise vector C 𝑁 ∋ n ∼ 𝑝 n sampled form a discrete random process {𝑛 𝑘 : ∀𝑘 } whose unspecified joint pdf has finite first and second order moments, that is it admits a zero-mean 𝐸 𝑝 n [n] = 0 and a covariance matrix 𝚺 ≜ 𝐸 𝑝 n [nn 𝐻 ] statisfying (3) and ( 4).

It can be noted that the equality in [START_REF] Fortunati | Semiparametric inference and lower bounds for real elliptically symmetric distributions[END_REF] does not requires the polynomial decrease of the autocorrelation function introduced in Assumption 1. However, we will see that this requirement will be crucial to derive asymptotic results about the efficiency of misspecifed Gaussian procedures.

IV. CLOSED FORM EXPRESSION FOR THE MCRB

The aim of this section is to provide the closed form expression of the Misspecified Cramér-Rao Bound (MCRB) for the estimation of ε under the misspecified scenario discussed in Sec. II-B. Following [START_REF] Vuong | Cramér-Rao bounds for misspecified models[END_REF], [4, Theo. 1] and [5, Theo. 4.1] and by exploiting the equality between the true and the pseudo-true paramater vectors, the MCRB is given by:

MCRB(𝝐 0 ) = MCRB( ε) = A( ε) -1 B( ε)A( ε) -1 , (18)
where:

[A( ε)] 𝑖, 𝑗 ≜ 𝐸 𝑝 ε ∇ 𝝐 ∇ ⊤ 𝝐 ln 𝑓 𝝐 (x; ε) 𝑖, 𝑗 = 𝐸 𝑝 ε 𝜕 2 𝜕 𝑖 𝜕 𝑗 ln 𝑓 𝝐 (x; 𝝐) 𝝐 =ε , (19)
[B( ε)] 𝑖, 𝑗 ≜ 𝐸 𝑝 ε ∇ 𝝐 ln 𝑓 𝝐 (x; ε)∇ ⊤ 𝝐 ln 𝑓 𝝐 (x; ε) 𝑖, 𝑗 = 𝐸 𝑝 ε 𝜕 𝜕 𝑖 ln 𝑓 𝝐 (x; 𝝐) 𝝐 =ε 𝜕 𝜕 𝑗 ln 𝑓 𝝐 (x; 𝝐) 𝝐 =ε , (20)
where 𝑓 𝝐 (x; 𝝐) ∈ F 𝝐 in [START_REF] Fortunati | The constrained Misspecified Cramér-Rao bound[END_REF].

The calculation of the matrices A( ε) and B( ε) will be performed in four steps:

1) Evaluation of the terms related to σ2 𝑛 . Through direct calculation, we have:

∇ 𝜎 2 𝑛 ln 𝑓 𝝐 (x; ε) = 𝜕 ln 𝑓 𝝐 (x; σ2 𝑛 , θ) 𝜕𝜎 2 𝑛 (21) = - 𝑁 𝜎 2 𝑛 + 1 𝜎 4 𝑛 ∥x -f ( θ) ∥ 2 𝜎 2 𝑛 = σ2 𝑛 = - 𝑁 σ2 𝑛 + tr(nn 𝐻 ) σ4 𝑛 ,
and then:

∇ 𝜎 2 𝑛 ∇ ⊤ 𝜎 2 𝑛 ln 𝑓 𝝐 (x; ε) = 𝑁 𝜎 4 𝑛 - 2 𝜎 6 𝑛 ∥x -f ( θ) ∥ 2 𝜎 2 𝑛 = σ2 𝑛 = 𝑁 σ4 𝑛 - 2tr(nn 𝐻 ) σ6 𝑛 . (22)
By taking the expectation w.r.t. the true data distribution 𝑝 ε and following the same calculation done in [START_REF] Van Der | Asymptotic Statistics[END_REF], we get:

𝐸 𝑝 ε ∇ 𝜎 2 𝑛 ∇ ⊤ 𝜎 2 𝑛 ln 𝑓 𝝐 (x; ε) = 𝑁 σ4 𝑛 - 2tr (𝚺) σ6 𝑛 = - 𝑁 σ4 𝑛 , (23)
where we used their linearity to invert the order of the expectation and trace operators. Similarly, we have that;

𝐸 𝑝 ε ∇ 𝜎 2 𝑛 ln 𝑓 𝝐 (x; ε) 2 = 𝐸 𝑝 n - 𝑁 σ2 𝑛 + n 𝐻 n σ4 𝑛 2 = 𝑁 2 σ4 𝑛 - 2𝑁tr (𝚺) σ6 𝑛 + 𝐸 𝑝 n (n 𝐻 n) 2 σ8 𝑛 = 𝐸 𝑝 n (n 𝐻 n) 2 -σ4 𝑛 𝑁 2 σ8 𝑛 . ( 24
)
Note that the term 𝐸 𝑝 n (n 𝐻 n) 2 cannot be further developed without specifying the true pdf of the noise 𝑝 n . We will further discuss this point in the next section. 2) Evaluation of the terms related to θ

From the assumed Gaussian pdf in eq. ( 8), we have:

∇ 𝜽 ln 𝑓 𝝐 (x; σ2 𝑛 , θ) = - 1 σ2 𝑛 ∇ 𝜽 ∥x -f ( θ) ∥ 2 = 1 σ2 𝑛 𝑁 2 ∑ 𝑘=𝑁 1 𝑥 𝑘 -𝑓 𝑘 ( θ) * ∇ 𝜽 f𝑘 + 𝑥 𝑘 -𝑓 𝑘 ( θ) ∇ * 𝜽 f𝑘 = 2 σ2 𝑛 𝑁 2 ∑ 𝑘=𝑁 1 Re 𝑛 * 𝑘 ∇ 𝜽 f𝑘 , (25)
where, for ease of notation, we posed ∇ 𝜽 𝑓 𝑘 ( θ) = ∇ 𝜽 f𝑘 . According to the eq. ( 20), we can evaluate the matrix B( θ) as showed in eq. ( 26) reported at the bottom of this page. It is worth noticing that, in the step (26c), we used the circularity assumption on {𝑛 𝑘 : ∀𝑘 }, i.e. 𝐸 𝑝 𝑛 [𝑛 𝑘 𝑛 𝑗 ] = 0, ∀𝑘, 𝑗 (see Assumption 1). The matrix P( θ) in (26f) has been introduced for further reference. Moreover, again thorough direct calculation, we have :

∇ 𝜽 ∇ ⊤ 𝜽 ln 𝑓 𝝐 (x; σ2 𝑛 , θ) = 1 σ2 𝑛 𝑁 2 ∑ 𝑘=𝑁 1 𝑥 𝑘 -𝑓 𝑘 ( θ) * ∇ 𝜽 ∇ ⊤ 𝜽 f𝑘 -∇ 𝜽 f𝑘 ∇ 𝐻 𝜽 f𝑘 × + 1 σ2 𝑛 𝑁 2 ∑ 𝑘=𝑁 1 𝑥 𝑘 -𝑓 𝑘 ( θ) ∇ 𝜽 ∇ ⊤ 𝜽 f𝑘 * -∇ * 𝜽 f𝑘 ∇ ⊤ 𝜽 f𝑘 = 2 σ2 𝑛 𝑁 2 ∑ 𝑘=𝑁 1 Re 𝑥 𝑘 -𝑓 𝑘 ( θ) ∇ 𝜽 ∇ ⊤ 𝜽 f𝑘 * (27a) × - 2 σ2 𝑛 𝑁 2 ∑ 𝑘=𝑁 1 Re ∇ 𝜽 f𝑘 ∇ 𝐻 𝜽 f𝑘 . ( 27b
)
Then we can introduce the matrix A( θ) as

A( θ) ≜ 𝐸 𝑝 ε ∇ 𝜽 ∇ ⊤ 𝜽 ln 𝑓 𝝐 (x; σ2 𝑛 , θ) (28) = - 𝑁 σ2 𝑛 2 𝑁 𝑁 2 ∑ 𝑘=𝑁 1 Re ∇ 𝜽 f𝑘 ∇ 𝐻 𝜽 f𝑘 ≜ - 𝑁 σ2 𝑛 K( θ),
where, again, we have introduced the matrix K( θ) for further reference. Note that the expectation of the term in (27a) is nil since

𝐸 𝑝 ε 𝑥 𝑘 -𝑓 𝑘 ( θ) = 𝐸 𝑝 ε [𝑛 𝑘 ] = 0, ∀𝑘.

3) Evaluation of the cross-terms

From the circularity of the noise process (see Assumption 1), it is immediate to verify that:

𝐸 𝑝 ε ∇ 𝜎 2 𝑛 ln 𝑓 𝝐 (x; ε)∇ ⊤ 𝜽 ln 𝑓 𝝐 (x; σ2 𝑛 , θ) = 𝐸 𝑝 ε ∇ 𝜎 2 𝑛 ln 𝑓 𝝐 (x; ε)∇ 𝜽 ln 𝑓 𝝐 (x; σ2 𝑛 , θ) ⊤ = 0 1×4 . ( 29
)
Moreover, we have that:

∇ 𝜎 2 𝑛 ∇ 𝜽 ln 𝑓 𝝐 (x; σ2 𝑛 , θ) = ∇ 𝜽 ∇ 𝜎 2 𝑛 ln 𝑓 𝝐 (x; σ2 𝑛 , θ) = 2 σ4 𝑛 𝑁 2 ∑ 𝑘=𝑁 1 Re 𝑛 * 𝑘 ∇ 𝜽 f𝑘 . (30)
Consequently, since the noise process is zero-mean, we trivially have that:

𝐸 𝑝 ε ∇ 𝜎 2 𝑛 ∇ ⊤ 𝜽 ln 𝑓 𝝐 (x; σ2 𝑛 , θ) = 𝐸 𝑝 ε ∇ 𝜽 ∇ 𝜎 2 𝑛 ln 𝑓 𝝐 (x; σ2 𝑛 , θ) ⊤ = 0 1×4 . (31)
4) Definition of the matrices A( ε) and B( ε) By collecting the previous results, we have that the matrix A( ε) in eq. ( 19) can be expressed as:

A( ε) = 𝑁 -1/ σ4 𝑛 0 1×4 0 4×1 -1 σ2 𝑛 K( θ) . (32)
Similarly, for the matrix B( ε) in eq. ( 20), we have:

B( ε) = 𝑁 (𝐸 𝑝n [(n 𝐻 n) 2 ]-σ4 𝑛 𝑁 2 ) 𝑁 σ8 𝑛 0 1×4 0 4×1 1 σ4 𝑛 P( θ) . (33)
As we can see from eq. ( 33), the matrix B( ε) is function of the matrix P( θ) in (26f) and the term

𝐸 𝑝 n [(n 𝐻 n) 2
] that depends on the autocorrelation function 𝑟 𝑛 [ 𝑗] involving summations over the index 𝑗 that can go until |𝑁 2 -𝑁 1 + 1| = 𝑁. Now, in order to provide asymptotic results on the number of observations, i.e. as 𝑁 → ∞, the norm of the matrix matrix B( ε) has to remain bounded as 𝑁 → ∞. As discussed in [START_REF] White | Nonlinear regression with dependent observations[END_REF], the polynomial decrease of 𝑟 𝑛 [ 𝑗] is needed to guarantees that B( ε) will not explode as 𝑁 → ∞. 3 Finally, the MCRB in ( 18) can be expressed as:

MCRB( ε) = A( ε) -1 B( ε)A( ε) -1 = 1 𝑁 (𝐸 𝑝 n (n 𝐻 n) 2 -σ4 𝑛 𝑁 2 )/𝑁 0 1×4 0 4×1 C( θ) , (34)
where

C( θ) ≜ K( θ) -1 P( θ)K( θ) -1 . ( 35
)
It is important to note that, due to the block-diagonal structure of MCRB( ε), the MCRB of the parameter of 3 Roughly speaking, the Assumption 1 guarantees the existence of a matrix B 0 , such that det(B 0 ) > 0 and a (B 0 -B( ε ) )a → 0 as 𝑁 → ∞, for any non-zero real vector a ∈ R 𝑝+1 .

interest vector θ can be simply obtained as:

MCRB( θ) = 𝑁 -1 C( θ). (36)
Remarkably, this results tells us that the estimation of θ is asymptotically decorrelatd from the nuisance parameter σ2 𝑛 .

A. Consistent estimation of the matrix P( θ)

Let us take a closer look to the matrix P( θ) in eq. (26e). It can be immediately noted that it depends on the apriori knowledge of the autocorrelation function of the noise

𝑟 𝑛 [ 𝑗] ≜ 𝐸 𝑝 n [𝑛 * 𝑘+ 𝑗 𝑛 𝑘 ].
However, to evaluate it, we need to know the true pdf 𝑝 n of the noise. This is in contrast with the Assumption 1 where 𝑝 n is left fully unspecified. We should then rely on a consistent estimator P 𝑁 of P( θ). Thanks to Assumption 1, deriving such consistent estimator is possible, even in presence of dependent observations. Following [START_REF] White | Nonlinear regression with dependent observations[END_REF], let us define the estimator P 𝑁 of P( θ) as:

P 𝑁 = 2 𝑁 ∑ 𝑁 2 𝑘=𝑁 1 | n𝑘 | 2 Re ∇ 𝜽 f𝑘 ∇ 𝐻 𝜽 f𝑘 + 4 𝑁 ∑ 𝑙 𝑗=1 ∑ 𝑁 2 -𝑗 𝑘=𝑁 1 Re n * 𝑘+ 𝑗 n𝑘 ∇ 𝜽 f𝑘+ 𝑗 ∇ 𝐻 𝜽 f𝑘 , (37)
where n𝑘 ≜ 𝑥 𝑘 -𝑓 𝑘 ( θ 𝑁 ) and ∇ 𝜽 f𝑘 ≜ ∇ 𝜽 𝑓 𝑘 ( θ 𝑁 ) and θ 𝑁 is a √ 𝑁-consistent estimator of the true parameter vector θ. Among all the possible consistent estimators, the best choice is the asymptotic efficient one that we are going to introduce in the subsequent section VI. The consistency of the estimator P 𝑁 in [START_REF] Lubeigt | Joint delay-doppler estimation performance in a dual source context[END_REF] is established in [START_REF] White | Nonlinear regression with dependent observations[END_REF]Theo. 3.5]:

Theorem 1: Under Assumption 1 and other technical regularity conditions (see A1a, A3, A4 and A7 in [START_REF] White | Nonlinear regression with dependent observations[END_REF]), if the correlation lag 𝑙 grows at the rate 𝑙 = 𝑜(𝑁 1/3 ) as 𝑁 → ∞, 4we have that P𝑁 is a consistent estimator of P( θ):

P 𝑁 𝑝 → P 𝑁 ( θ), (38)
where 𝑝 → indicates the convergence (element by element) in probability.

B( θ) ≜ 𝐸 𝑝 ε ∇ 𝜽 ln 𝑓 𝜽 (x; ε)∇ ⊤ 𝝐 ln 𝑓 𝝐 (x; ε) = 𝐸 𝑝 n 2 σ2 𝑛 ∑ 𝑁 2 𝑘=𝑁 1 Re 𝑛 * 𝑘 ∇ 𝜽 f𝑘 2 σ2 𝑛 ∑ 𝑁 2 𝑗=𝑁 1 Re 𝑛 * 𝑗 ∇ 𝜽 f 𝑗 𝑇 (26a) = 4 σ4 𝑛 𝑁 2 ∑ 𝑘=𝑁 1 𝑁 2 ∑ 𝑗=𝑁 1 𝐸 𝑝 n Re 𝑛 * 𝑘 ∇ 𝜽 f𝑘 Re 𝑛 * 𝑗 ∇ 𝑇 𝜽 f 𝑗 = 4 σ4 𝑛 𝑁 2 ∑ 𝑘=𝑁 1 𝑁 2 ∑ 𝑗=𝑁 1 𝐸 𝑝 n 𝑛 * 𝑘 ∇ 𝜽 f𝑘 + 𝑛 𝑘 ∇ * 𝜽 f𝑘 2 𝑛 * 𝑗 ∇ 𝑇 𝜽 f 𝑗 + 𝑛 𝑗 ∇ 𝐻 𝜽 f 𝑗 2 (26b) = 2 σ4 𝑛 ∑ 𝑁 2 𝑘=𝑁 1 ∑ 𝑁 2 𝑗=𝑁 1 Re 𝐸 𝑝 n 𝑛 𝑘 𝑛 𝑗 * ∇ 𝜽 f𝑘 ∇ 𝑇 𝜽 f 𝑗 + Re 𝐸 𝑝 n 𝑛 * 𝑘 𝑛 𝑗 ∇ 𝜽 f𝑘 ∇ 𝐻 𝜽 f 𝑗 (26c) = 2 σ4 𝑛 ∑ 𝑁 2 𝑘=𝑁 1 𝐸 𝑝 n |𝑛 𝑘 | 2 Re ∇ 𝜽 f𝑘 ∇ 𝐻 𝜽 f𝑘 + 2 σ4 𝑛 ∑ 𝑁 2 𝑘=𝑁 1 ∑ 𝑁 2 𝑗=𝑁 1 𝑗≠𝑘 Re 𝐸 𝑝 n 𝑛 * 𝑘 𝑛 𝑗 ∇ 𝜽 f𝑘 ∇ 𝐻 𝜽 f 𝑗 (26d) = 2 σ4 𝑛 ∑ 𝑁 2 𝑘=𝑁 1 𝐸 𝑝 n |𝑛 𝑘 | 2 Re ∇ 𝜽 f𝑘 ∇ 𝐻 𝜽 f𝑘 + 4 σ4 𝑛 ∑ 𝑁 2 -𝑁 1 𝑗=1 ∑ 𝑁 2 -𝑗 𝑘=𝑁 1 Re 𝐸 𝑝 n 𝑛 * 𝑘+ 𝑗 𝑛 𝑘 ∇ 𝜽 f𝑘+ 𝑗 ∇ 𝐻 𝜽 f𝑘 (26e) = 𝑁 σ4 𝑛 2 𝑁 ∑ 𝑁 2 𝑘=𝑁 1 𝑟 𝑛 [0]Re ∇ 𝜽 f𝑘 ∇ 𝐻 𝜽 f𝑘 + 4 𝑁 ∑ 𝑁 2 -𝑁 1 𝑗=1 ∑ 𝑁 2 -𝑗 𝑘=𝑁 1 Re 𝑟 𝑛 [ 𝑗]∇ 𝜽 f𝑘+ 𝑗 ∇ 𝐻 𝜽 f𝑘 ≜ 𝑁 σ4 𝑛 P( θ). (26f)
By a direct application of the Continuous Mapping Theorem and of the Slutsky's Lemma [16, Theo. 2.3 and Lemma 2.8],

we have that the matrix C( θ) in ( 35) can be consistently estimated as:

C 𝑁 ≜ K( θ 𝑁 ) -1 P 𝑁 K( θ 𝑁 ) -1 𝑝 → C( θ), (39)
that can be exploited to get a consistent estimation of the MCRB on the vector of the parameter of interest θ.

V. MCRB EXPRESSION FOR CES UNCORRELATED DATA

In order to highlight the importance and the generality of the results obtained in the previous section, let us consider the particular case, yet important in applications, where the noise process {𝑛 𝑘 : ∀𝑘 } is assumed to be white with independent and identically Complex Elliptically Symmetric (CES)-distributed samples. More formally, we assume that:

Assumption 2: Let {𝑛 𝑘 : ∀𝑘 } be a zero-mean, white WSS discrete and circular complex-valued process [START_REF] Picinbono | On circularity[END_REF] such that:

1) each sample 𝑛 𝑘 follows a CES distribution 𝑛 𝑘 ∼ 𝑝 𝑛 = 𝐶𝐸 𝑆(0, σ2 𝑛 , 𝑔) with unspecified density generator 𝑔, 2) its autocorrelation satisfies

𝑟 𝑛 [ 𝑗] ≜ 𝐸 𝑝 n [𝑛 * 𝑘+ 𝑗 𝑛 𝑘 ] = σ2 𝑛 𝛿[ 𝑗],
where 𝛿[ 𝑗] is the Kronecker delta sequence. As a direct consequence of Assumption 2, we have that:

• the joint pdf 𝑝 n of the noise vector n ∼ 𝑝 n is the product of the marginal densities, i.e. 𝑝 n (n; σ2

𝑛 , 𝑔) = 𝑁 2 𝑘=𝑁 1 𝑝 𝑛 𝑘 (𝑛 𝑘 ; σ2 𝑛 , 𝑔) • the covariance matrix of n ∈ C 𝑁 is a diagonal matrix, i.e.
𝚺 ≜ 𝐸 𝑝 𝑛 [nn 𝐻 ] = 𝜎 2 𝑛 I. It is worth stressing here the generality of an unspecified CES distribution for the noise samples. The CES ones is a wide class of non-Gaussian and heavy-tailed distributions encompassing the Gaussian, the Generalized Gaussian, the 𝑡-, the 𝐾-and the Weibull distributions as special cases [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]. Since its nominal density generator 𝑔 is left unspecified, we let the noise 𝑛 𝑘 ∼ 𝑝 𝑛 𝑘 have any possible distribution in the CES class.

From the Stochastic Representation Theorem [14, Theo. 3], each entry 𝑛 𝑘 can be represented as [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF]Theo. 3]:

𝑛 𝑘 = 𝑑 √ 𝑄 𝑘 σ𝑛 𝑢 𝑘 , (40)
where 𝑢 𝑘 is a complex univariate random variable uniformly distributed on C𝑆 ≜ {𝑢 ∈ C||𝑢| = 1}, i.e. 𝑢 𝑘 ∼ 𝑈 (C𝑆). The second order modular variate 𝑄 𝑘 ∼ Q is a positive random variable, independent from 𝑢 𝑘 with pdf 𝑝 Q (𝑞) = 𝛿 -1 𝑔 𝑔(𝑞), where 𝛿 𝑔 ≜ ∫ ∞ 0 𝑔(𝑞)𝑑𝑞 is a normalizing constant (see [14, Eq. ( 19)]). Since the density generator 𝑔 is left unspecified, it is immediate to verify that there is a scale ambiguity between σ2

𝑛 and 𝑔 itself. To avoid this problem, we impose that 𝐸 {Q} = 1. Note that, this constraint allows us to consider σ2 𝑛 as the statistical power 𝑃 of the noise 𝑛 𝑘 , (see the discussion in [14, Sec. III.C]), since from (40), we have that:

𝑃 ≜ 𝐸 {|𝑥 𝑘 | 2 } = 𝐸 {Q}𝐸 {|𝑢 𝑘 | 2 } σ2 𝑛 = σ2 𝑛 , (41)
where

𝐸 {|𝑢 𝑘 | 2 } = 1 [14, Lemma 1].
Let us now apply the general expression of the MCRB obtained in [START_REF] Niu | Wavelet Based Approach for Joint Time Delay and Doppler Stretch Measurements[END_REF] to the special case of an i.i.d. CES-distributed noise process formally characterized in Assumption 2. To this end, it is immediate to verify that the matrix A( ε) in (32) remain unchanged. Let us now focus on the matrix B( ε) in [START_REF] Jin | The Estimation of Time Delay and Doppler Stretch of Wideband Signals[END_REF]. We will proceed as follows.

1) Evaluation of the term [B( ε)] 1,1 in (24). Let us start by evaluating the term 𝐸 𝑝 n (n 𝐻 n) 2 . Under Assumption 2 and by exploiting the stochastic representation in [START_REF] Mcphee | On the accuracy limits of misspecified delay-doppler estimation[END_REF], the term 𝐸 𝑝 n [(n 𝐻 n) 2 ] can be evaluated as shown in [START_REF] Ortega | On gnss synchronization performance degradation under interference scenarios: Bias and misspecified cramér-rao bounds[END_REF], reported at the bottom of the page,where we used:

• the mutual independence between 𝑄 𝑖 ∼ Q and 𝑄 𝑗 ∼ Q and between 𝑢 𝑖 and 𝑢 𝑗 (see Assumption 2), • from the i.i.d. assumption, we have that Lemma 1]. By using this result, the term [B( ε)] 1,1 can be readily expressed as:

𝐸 𝑄 2 𝑖 = 𝐸 Q 2 , ∀𝑖, • the constraint 𝐸 [𝑄 𝑖 ] = 𝐸 𝑄 𝑗 = 𝐸 [Q] = 1, • the relations 𝐸 |𝑢 𝑗 | 2 = 1 and 𝐸 |𝑢 𝑗 | 4 = 1 from [14,
[B 𝑖𝑖𝑑 ( ε)] 1,1 = 𝑁 (𝐸 Q 2 -1)/ σ4 𝑛 (43)
2) Evaluation of the matrix B( ε) By putting in the general expression of B( θ) in ( 26) the autocorrelation function

𝑟 𝑛 [ 𝑗] = σ2 𝑛 𝛿[ 𝑗] (see Assumption 2)
, it is immediate to verify that:

B 𝑖𝑖𝑑 ( θ) = 2 σ2 𝑛 𝑁 2 ∑ 𝑘=𝑁 1 Re ∇ 𝜽 f𝑘 ∇ 𝐻 𝜽 f𝑘 ≜ 𝑁 σ2 𝑛 K( θ), (44)
where K( θ) is the matrix already defined eq. ( 28). Consequently, the matrix B( ε), under Assuption 2, can be expressed as:

B 𝑖𝑖𝑑 ( ε) = 𝑁 (𝐸 {Q 2 } -1)/𝜎 4 𝑛 0 1×4 0 4×1 1 σ2 𝑛 K( θ) . (45)
Finally, from the general expression in [START_REF] Niu | Wavelet Based Approach for Joint Time Delay and Doppler Stretch Measurements[END_REF], the MCRB for the estimation of ε under Assumption 2 can be expressed as:

MCRB 𝑖𝑖𝑑 ( ε) = A( ε) -1 B 𝑖𝑖𝑑 ( ε)A( ε) -1 = 1 𝑁 𝜎 4 𝑛 (𝐸 {Q 2 } -1) 0 1×4 0 4×1 σ2 𝑛 K( θ) -1 , (46)
and consequently, due to the block-diagonal structure of MCRB 𝑖𝑖𝑑 ( ε), the MCRB on the vector of the parameters of interest θ is given by:

MCRB 𝑖𝑖𝑑 ( θ) = σ2 𝑛 𝑁 K( θ) -1 . ( 47
)
It is worth highlighting here an interesting result: under the misspecificed scenario discussed in this section, i.e. when the data follow a CES, i.i.d. (true) model while the assumed one is a Gaussian, i.i.d., model, we have that:

A( θ) + B 𝑖𝑖𝑑 ( θ) = 0. ( 48
)
As explained in [START_REF] Fortunati | Chapter 4 -Parameter bounds under misspecified models for adaptive radar detection[END_REF]Lemma 4.1], the result in (48), along with the block-diagonal structure of MCRB 𝑖𝑖𝑑 ( ε) in ( 46), implies that the simplified Gaussian assumption does not lead to any degradation of the asymptotic estimation performance of the parameter vector of interest θ. In fact, MCRB 𝑖𝑖𝑑 ( θ) coincides with the lower bound that we can get if the true data model was an i.i.d. Gaussian one. This intriguing outcome can be explained through the semiparametric theory (refer to [17, Sec. IV;B] and [18, Sec. III.B]) that allows us to prove that the lack of knowledge of the density generator 𝑔 does not have any asymptotic impact on the estimation of θ.

Finally, if the true distribution is a Gaussian one, the term 𝐸 {Q 2 } is equal to 2 as proved in [START_REF] Mennad | Slepian-bangs-type formulas and the related misspecified cramér-rao bounds for complex elliptically symmetric distributions[END_REF]Eq. (41)] and this lead us to the classical result about the CRB on the estimation of the variance in complex Gaussian data.

VI. AN ASYMPTOTIC EFFICIENT ESTIMATOR UNDER DEPENDENT OBSERVATIONS

Let us go back now to the general misspecified nonlinear regression problem presented in Sec. II. After having derived the MCRB for the vector of the parameters of interest θ in [START_REF] Das | A New Compact Delay, Doppler Stretch and Phase Estimation CRB with a Band-Limited Signal for GenE[END_REF], the crucial question that arises is as follows: is it possible to derive, under the misspecified Gaussian model F 𝝐 in (7), a consistent estimator θ 𝑁 of θ able to achieve the MCRB, at least asymptotically?

It is well known that, under the i.i.d. case, the answer to this question is positive and θ 𝑁 is given by the Missmatched Maximum Likelihood estimator (MMLE) [START_REF] Huber | The behavior of maximum likelihood estimates under nonstandard conditions[END_REF], [START_REF] White | Maximum likelihood estimation of misspecified models[END_REF], [START_REF] Fortunati | Performance bounds for parameter estimation under misspecified models: Fundamental findings and applications[END_REF]- [START_REF] Richmond | Parameter bounds on estimation accuracy under model misspecification[END_REF]. The extension to the dependent case has been provided in [START_REF] White | Nonlinear regression with dependent observations[END_REF] where the asymptotic behaviour of the nonlinear least square estimator (NLLSE) for θ under the dependent data generating process in [START_REF] White | Maximum likelihood estimation of misspecified models[END_REF], that is:

θ 𝑁 = argmin 𝜽 ∈Θ 1 𝑁 ∑ 𝑁 2 𝑘=𝑁 1 |𝑥 𝑘 -𝑓 𝑘 (𝜽)| 2 , ( 49
)
has been investigated. Note that, when the misspecified Gaussian model F 𝝐 in ( 7) is assumed, it is immediate to verify that, for any finite 𝑁, the NLLSE coincides with the MMLE. In fact, from (8), the misspecified log-likelihood function is

𝑙 (𝜽) = -𝑁 ln(𝜋𝜎 2 𝑛 ) -||x -f (𝜽)|| 2 /𝜎 2 𝑛 , then the MMLE for θ in F 𝝐 is given by: argmax 𝜽 ∈Θ 𝑙 (𝜽) = argmin 𝜽 ∈Θ ||x -f (𝜽)|| 2 = θ 𝑁 , (50)
that is the same estimator of the one in (49). Remarkably, in [START_REF] White | Nonlinear regression with dependent observations[END_REF] it has been proved that: Theorem 2: Under Assumption 1 and other technical regularity conditions (see A1-A9 in [START_REF] White | Nonlinear regression with dependent observations[END_REF]), θ 𝑁 in (49) satisfies the following properties:

1) Consistency wrt the true parameter vector:

θ 𝑁 𝑎.𝑠. → θ, (51)
𝑎.𝑠.

→ indicates the almost sure convergence.

2) Asymptotic normality: Let us indicate as ∼

𝑁 →∞ the convergence in distribution, we have:

√ 𝑁 P( θ) -1/2 K( θ) θ 𝑁 -θ ∼ 𝑁 →∞ N (0, I). (52)
The interested reader can find the proof of Theorem 2 in the supporting material associated to this paper. Our proof extends the one provided in [15, Theo. 3.1 and 3.2] to the case of complex observations. It follow directly from (52) that the asymptotic error covariance matrix of θ 𝑁 equates the MCRB( θ) in [START_REF] Das | A New Compact Delay, Doppler Stretch and Phase Estimation CRB with a Band-Limited Signal for GenE[END_REF], i.e.:

lim 𝑁 →∞ w 𝑇 𝑁 𝐸 𝑝 ε ( θ 𝑁 -θ) ( θ 𝑁 -θ) ⊤ -× × + K( θ) -1 P( θ)K( θ) -1 w = 0, ∀w ∈ R 𝑝 /{0}. (53)
Consequently, the NLLSE in (49) (that coincides with the MMLE under misspecified Gaussian assumption) is exactly the consistent and asymptotically efficient estimator that we were looking for.

To conclude, it can be noted that, in the i.i.d. case discussed in Sec. V, eq. ( 53) simplifies to:

lim 𝑁 →∞ w 𝑇 𝑁 𝐸 𝑝 ε ( θ 𝑁 -θ) ( θ 𝑁 -θ) ⊤ -× × + σ2 𝑛 K( θ) -1 w = 0, ∀w ∈ R 𝑝 /{0}. (54)

VII. APPLICATION TO TIME-DELAY AND DOOPLER

ESTIMATION Time-delay and Doppler estimation is fundamental in a plethora of engineering domains, including communications, radar, and navigation [START_REF] Swick | A Review of Wideband Ambiguity Functions[END_REF]- [START_REF]Handbook of Global Navigation Satellite Systems[END_REF], as it serves as the initial step at the receiver [START_REF] Ricker | Echo Signal Processing[END_REF], [START_REF] Munoz | Position Location Techniques and Applications[END_REF], [START_REF] Yan | Review of range-based positioning algorithms[END_REF]. Due to its importance, understanding the achievable estimation performance in terms of MSE is of paramount practical interest. This crucial insight is typically provided by the CRB. Over the past few decades, numerous CRB expressions have been developed for time-delay and Doppler estimation problems, encompassing both finite narrow-band and wideband signals [START_REF] Van Trees | Detection, Estimation, and Modulation Theory, Part III: Radar -Sonar Signal Processing and Gaussian Signals in Noise[END_REF], [START_REF] Dogandzic | Cramér-Rao bounds for estimating range, velocity, and direction with an active array[END_REF]- [START_REF] Mcphee | Accounting for acceleration-signal parameters estimation performance limits in high dynamics applications[END_REF]. Furthermore, recent studies have explored scenarios in which the actual signal model at the receiver differs from the assumed one [START_REF] Lubeigt | Untangling first and second order statistics contributions in multipath scenarios[END_REF]- [START_REF] Ortega | On gnss synchronization performance degradation under interference scenarios: Bias and misspecified cramér-rao bounds[END_REF]. In these investigations, expressions for estimation boundaries, as determined by the MCRB, have been established. However all these prior studies share a common assumption: both the noise in the true signal model and the noise in the signal model assumed by the receiver follow a centered complex normal distribution with uncorrelated covariance matrix, i.e. a diagonal matrix. Surprisingly, despite the extensive

𝐸 𝑝 n (n 𝐻 n) 2 = σ4 𝑛 𝐸 ∑ 𝑁 2 𝑖=𝑁 1 𝑄 𝑖 |𝑢 𝑖 | 2 2 = σ4 𝑛 𝐸 ∑ 𝑁 2 𝑖=𝑁 1 ∑ 𝑁 2 𝑗=𝑁 1 𝑄 𝑖 𝑄 𝑗 |𝑢 𝑖 | 2 |𝑢 𝑗 | 2 = σ4 𝑛 𝐸 ∑ 𝑁 2 𝑖=𝑁 1 𝑄 2 𝑖 |𝑢 𝑖 | 4 + × (42a) × + σ4 𝑛 𝐸 𝑁 2 ∑ 𝑖=𝑁 1 𝑁 2 ∑ 𝑗=𝑁 1 , 𝑗≠𝑖 𝑄 𝑖 𝑄 𝑗 |𝑢 𝑖 | 2 |𝑢 𝑗 | 2 = σ4 𝑛 𝑁 2 ∑ 𝑖=𝑁 1 𝐸 𝑄 2 𝑖 𝐸 |𝑢 𝑖 | 4 + 𝑁 2 ∑ 𝑖=𝑁 1 𝑁 2 ∑ 𝑗=𝑁 1 , 𝑗≠𝑖 𝐸 [𝑄 𝑖 ] 𝐸 𝑄 𝑗 𝐸 |𝑢 𝑖 | 2 𝐸 |𝑢 𝑗 | 2 (42b) = σ4 𝑛 𝑁 𝐸 Q 2 + 𝑁 (𝑁 -1) = σ4 𝑛 𝑁 (𝐸 Q 2 -1) + 𝑁 2 , (42c)
research in this area, there is a notable absence in the literature regarding the ultimate attainable estimation performance for time-delay and Doppler (in terms of MSE) when the true signal model features a correlated non-Gaussian distributed noise. The aim if this section os then to fill this gap by relying on the theoretical results derived in the previous sections.

A. Signal model

We consider the transmitter 𝑇 to receiver 𝑅 direct transmission of a band-limited signal 𝑎 (𝑡) with bandwidth 𝐵

𝑎 (𝑡) = 𝑁 2 ∑ 𝑛=𝑁 1 𝑎 (𝑛𝑇) sinc (𝜋𝐵 (𝑡 -𝑛𝑇)) , 𝑇 = 1/𝐵, (55)
over a carrier with frequency 𝑓 𝑐 (𝜆 𝑐 = 𝑐/ 𝑓 𝑐 , 𝜔 𝑐 = 2𝜋 𝑓 𝑐 ). The transmitter is located at position 𝑷 𝑇 (𝑡) and the receiver is located at position 𝑷 𝑅 (𝑡). The distance travelled by the transmitted signal is 𝑷 𝑇 𝑅 = ∥𝑷 𝑇 (𝑡 -𝜏 0 (𝑡)) -𝑷 𝑅 (𝑡) ∥≈ (𝑷 𝑇 -𝑷 𝑅 ) 𝑐 + 𝑣 𝑐 𝑡, that is, a first order approximation where τ = (𝑷 𝑇 -𝑷 𝑅 ) 𝑐 and b = 𝑣 𝑐 with 𝑣 the relative velocity between the transmitter and the receiver. Once the baseband demodulation process has been completed, the received signal at the output of the Hilbert filter can be expressed as [START_REF] Dogandzic | Cramér-Rao bounds for estimating range, velocity, and direction with an active array[END_REF], [START_REF] Medina | Compact CRB for delay, doppler and phase estimation -application to GNSS SPP & RTK performance characterization[END_REF], [START_REF] Skolnik | Radar Handbook[END_REF]

𝑥 (𝑡; η) = ᾱ𝑎 (𝑡 -τ) (1 -b) 𝑒 -𝑗2 𝜋 𝑓 𝑐 ( b (𝑡 -τ ) ) + 𝑛 (𝑡) , (56) yielding to 𝑥 (𝑡; η) = ᾱ𝑎 (𝑡 -τ) 𝑒 -𝑗2 𝜋 𝑓 𝑐 ( b (𝑡 -τ ) ) + 𝑛 (𝑡) , (57)
under the narrowband assumption, i.e. the influence of the Doppler parameter on the baseband signal samples is omitted. The term ᾱ = ρ𝑒 𝑗 Φ represents a complex gain, while 𝑛(𝑡) is a zero-mean, generally non Gaussian, wide sense stationary (WSS) continuous random process. The discrete signal model is built from

𝑁 = |𝑁 1 -𝑁 2 + 1| samples at 𝑇 𝑠 = 1/𝐹 𝑠 = 1/𝐵, x = ᾱ𝝁( η) + n = ρ𝑒 𝑗 Φ 𝝁( η) + n, (58)
with x = (. . . , 𝑥 (𝑘𝑇 𝑠 ) , . . .) ⊤ , 𝑁 1 ≤ 𝑘 ≤ 𝑁 2 signal samples. Moreover, by posing ( η) = [ τ, b] ⊤ , we have:

𝝁( η) = (. . . , 𝑎(𝑘𝑇 𝑠 -τ)𝑒 -𝑗2 𝜋 𝑓 𝑐 ( b (𝑘𝑇 𝑠 -τ ) , . . .) ⊤ . (59)
Consequently, by defining the true vector of the parameters of interest as θ⊤ = ρ, Φ, η⊤ the signal model in eq. ( 58) follows the form in (1):

x = ᾱ𝝁( η) + n = f ( θ) + n. (60)
Finally, standard receivers assumes that the noise vector n ∈ C 𝑁 is distributed as a centered complex normal random vector with diagonal covariance matrix, i.e. n ∼ CN (0, 𝜎 2 𝑛 𝑰 𝑁 ). Note that this represents the same misspecified scenario introduced in section II-B. Specifically, we have that the pdf of the observation vector x in (60) belongs to the misspecified model in [START_REF] Fortunati | The constrained Misspecified Cramér-Rao bound[END_REF], i.e. x ∼ 𝑓 𝝐 ∈ F 𝝐 .

B. Time-delay and Doppler Closed-Form MCRB Expression for a Band-Limited Signal

It is interesting to note the likelihood between the expression obtained in previous sections and those already derived in the state of art. In particular, we may note that the matrix -A( θ) derived in [START_REF]Handbook of Global Navigation Satellite Systems[END_REF] represents the FIM of a single source conditional signal model (CSM) [START_REF] Stoica | Performances study of conditional and unconditional direction of arrival estimation[END_REF]. A compact expression of this FIM, that depends only on the baseband signal samples, was recently derived in [START_REF] Medina | Compact CRB for delay, doppler and phase estimation -application to GNSS SPP & RTK performance characterization[END_REF] as:

-A( θ) = 2𝐹 𝑠 𝜎 2 𝑛 Re QWQ 𝐻 (61)
with

W =       𝑤 1 𝑤 * 2 𝑤 * 3 𝑤 2 𝑊 2,2 𝑤 * 4 𝑤 3 𝑤 4 𝑊 3,3       , (62a)
Q =         𝑒 𝑗 Φ 0 0 𝑗 ᾱ 0 0 𝑗 ᾱ2𝜋 𝑓 𝑐 b 0 - ᾱ 0 -𝑗 ᾱ2𝜋 𝑓 𝑐 0         , (62b)
where the elements of W can be expressed w.r.t. the baseband signal samples as,

𝑤 1 = 1 𝐹 𝑠 a 𝐻 a, 𝑤 2 = 1 𝐹 2 𝑠 a 𝐻 Da, 𝑤 3 = a 𝐻 𝚲a, (63)
𝑤 4 = 1 𝐹 𝑠 a 𝐻 D𝚲a, 𝑊 2,2 = 1 𝐹 3 𝑠 a 𝐻 D 2 a, 𝑊 3,3 = 𝐹 𝑠 a 𝐻 Va.
with a, the baseband samples vector, D, 𝚲 and V defined as, a = (. . . , 𝑎(𝑛𝑇 𝑠 ), . . .

) ⊤ 𝑁 1 ≤𝑛≤ 𝑁 2 , (64a)
D = diag (. . . , 𝑛, . . .) 𝑁 1 ≤𝑛≤ 𝑁 2 , (64b)
(𝚲) 𝑛,𝑛 ′ = 𝑛 ′ ≠ 𝑛 : (-1) |𝑛-𝑛 ′ | 𝑛-𝑛 ′ 𝑛 ′ = 𝑛 : 0 (64c) (V) 𝑛,𝑛 ′ = 𝑛 ′ ≠ 𝑛 : (-1) |𝑛-𝑛 ′ | 2 (𝑛-𝑛 ′ ) 2 𝑛 ′ = 𝑛 : 𝜋 2 3 (64d)
Moreover, under the uncorrelated noise assumption, we note from (46) that, since B 𝑖𝑖𝑑 ( θ) = -A( θ), the MCRB on the estimation of θ is given by:

MCRB 𝑖𝑖𝑑 ( θ) = 2𝐹 𝑠 σ2 𝑛 Re QWQ 𝐻 (65)
In the case of correlated noise, the expression of the matrix B( θ) = 𝑁 -1 σ4 𝑛 P( θ) in ( 26) is more challenging since P( θ) involves the autocorrelation function 𝑟 𝑛 [ 𝑗]. We can distinguish between the following two cases. If 𝑟 𝑛 [ 𝑗] is a-priori known, for the application at hand P( θ) can be expressed as:

P( θ) = 2𝐹 𝑠 σ2 𝑛 𝑁 Re QWQ 𝐻 + 4𝐹 𝑠 𝑁 𝑙 ∑ 𝑗=1 Re 𝑟 𝑛 [ 𝑗]QW 𝑗 Q 𝐻 (66)
where:

W 𝑗 =        𝑤 ( 𝑗 ) 1 𝑤 ( 𝑗 ) * 2 𝑤 ( 𝑗 ) * 3 𝑤 ( 𝑗 ) 2 𝑊 ( 𝑗 ) 2,2 𝑤 ( 𝑗 ) * 4 𝑤 ( 𝑗 ) 3 𝑤 ( 𝑗 ) 4 𝑊 ( 𝑗 ) 3,3        , (67)
and the elements of W 𝑗 can be expressed w.r.t. the baseband signal samples as,

𝑤 ( 𝑗 ) 1 = 1 𝐹 𝑠 a 𝐻 𝑗 + a 𝑗 -, 𝑤 ( 𝑗 ) 2 = 1 𝐹 2 𝑠 a 𝐻 𝑗 + D 𝑗 a 𝑗 -, (68)
𝑤 ( 𝑗 ) 3 = a 𝐻 𝑗 + 𝚲 𝑗 a 𝑗 -, 𝑤 ( 𝑗 ) 4 = 1 𝐹 𝑠 a 𝐻 𝑗 + D 𝑗 𝚲 𝑗 a 𝑗 -, 𝑊 ( 𝑗 ) 2,2 = 1 𝐹 3 𝑠 a 𝐻 𝑗 + D 2 𝑗 a 𝑗 -, 𝑊 ( 𝑗 ) 3,3 = 𝐹 𝑠 a 𝐻 𝑗 + V 𝑗 a 𝑗 -.
with a 𝑗 -, a 𝑗 + , D 𝑗 , 𝚲 𝑗 and V 𝑗 defined as,

a 𝑗 + = (. . . , 𝑎(𝑛𝑇 𝑠 ), . . .) ⊤ 𝑁 1 + 𝑗 ≤𝑛≤ 𝑁 2 , (69a)
a 𝑗 -= (. . . , 𝑎(𝑛𝑇 𝑠 ), . . .) ⊤ 𝑁 1 ≤𝑛≤ 𝑁 2 -𝑗 , (69b)
D 𝑗 = diag (. . . , 𝑛, . . .) 𝑁 1 ≤𝑛≤ 𝑁 2 -𝑗 , (69c)
𝚲 𝑗 𝑛,𝑛 ′ = 𝑛 ′ ≠ 𝑛 : (-1) |𝑛-𝑛 ′ | 𝑛-𝑛 ′ 𝑛 ′ = 𝑛 : 0 (69d) V 𝑗 𝑛,𝑛 ′ = 𝑛 ′ ≠ 𝑛 : (-1) |𝑛-𝑛 ′ | 2 (𝑛-𝑛 ′ ) 2 𝑛 ′ = 𝑛 : 𝜋 2 3 (69e) with 𝑁 1 ≤ 𝑛, 𝑛 ′ ≤ 𝑁 2 -𝑗. If 𝑟 𝑛 [ 𝑗]
is not a-priori known, from (37), a consistent estimator of P( θ) can be implemented as:

P 𝑁 = 𝑁 -1 2𝐹 𝑠 Re QWQ 𝐻 ∑ 𝑁 2 𝑘=𝑁 1 | n𝑘 | 2 + + 𝑁 -1 4𝐹 𝑠 ∑ 𝑙 𝑗=1 Re n * 𝑘+ 𝑗 n𝑘 QW 𝑗 Q 𝐻 , (70)
where n𝑘 ≜ 𝑥 𝑘 -𝑓 𝑘 ( θ 𝑁 ) and θ 𝑁 is the consistent estimator of θ defined as [START_REF] Ortega | On gnss synchronization performance degradation under interference scenarios: Bias and misspecified cramér-rao bounds[END_REF] 5 :

η = arg max 𝜼 𝚷 𝝁 (𝜼) x 2 (71) 𝜌 = 𝝁 𝐻 ( η) 𝝁 ( η) -1 𝝁 𝐻 ( η) x (72) Φ = arg 𝝁 𝐻 ( η) 𝝁 ( η) -1 𝝁 𝐻 ( η) x (73)

VIII. SIMULATION AND DISCUSSION

To support our theoretical analysis, we examine the transmission and reception of a GPS L1 C/A signal [START_REF]Handbook of Global Navigation Satellite Systems[END_REF]. This signal employs a baseband signal represented by a periodic binary phase-shift keying (BPSK) Gold code with a length of 1023 chips of period 1ms. At the receiver, we set a sampling frequency 𝐹 𝑠 = 4 MHz, which is the standard rate for most commercial receivers. The GNSS receiver assumes that the noise follows a standard centered normal distribution.

Scenario 1 In a first scenario, we set a true signal model where the noise is sampled from a complex centered 𝑡distribution [5, Sec. 4.6.1.1] with 𝜐 > 1 degrees of freedom (or shape parameter) that control the level of non-Gaussianity and a scale parameter 𝜇. The second-order modular variate Q of a 𝑡-distribution is an 𝐹-distributed random variable with parameters 2 and 𝜐 i.e. Q ∼ 𝐹 (2, 𝜐) [14, Sec. IV.A]. Then, in order to meet the constraint 𝐸 {Q} = 1, the scale as to be set as 𝜇 = 𝜐 σ2 𝑛 ( 𝜐-1) where σ2 𝑛 depends on the signal to noise 5 Let 𝑆 = 𝑠 𝑝𝑎𝑛 (A), with A a matrix, be the linear span of the set of its column vectors. The orthogonal projector over 𝑆 is 𝚷 A = A A 𝐻 A A 𝐻 . ratio at the output of the match filter 𝑆𝑁 𝑅 𝑜𝑢𝑡 . The 𝑆𝑁 𝑅 𝑜𝑢𝑡 is defined as:

𝑆𝑁 𝑅 𝑜𝑢𝑡 = |𝛼| 2 a 𝐻 a σ2 𝑛 . (74)
Furthermore, in this scenario, we employ two autoregressive processes (AR) of order 1 and 6, respectively, to model the noise correlation. The poles of the process are set to 𝑝 = 0.9•𝑒 𝑗2 𝜋•0. the 𝑆𝑁 𝑅 𝑜𝑢𝑡 . The number of Monte Carlo is set to 1000 iterations and 𝜐 = 2.5. In the results one can observe that the RMSE ( √ 𝑀𝑆𝐸) of the pseudotrue parameter converges to the asymptotic estimation performance derived in Section IV. These results confirm the theoretical derivation. Moreover, the Gaussian i.d.d

√

𝐶 𝑅𝐵 has been included to quantify the performance with respect to the correlated case. It is worth to underline that the previous theoretical results are valid for any joint pdf 𝑝 n admitting finite first and second order moments and not only of the one obtained from a 𝑡-distribution.

Scenario 2 In the second scenario, we would like to illustrate the estimation performance in the case where the noise is i.i.d. To do so, we propose to set the true signal model as in the previous scenario, i.e the noise process is complex centered 𝑡-distributed (we remind that in this particular case, the noise is assumed to be i.i.d and there is not need to define any process to characterized the noise correlation). Moreover, we also set a true signal model where the noise is distributed according to a complex centered Generalized Gaussian (GG) distribution, [START_REF] Fortunati | Chapter 4 -Parameter bounds under misspecified models for adaptive radar detection[END_REF]Sec. 4.6.1.2] with exponent 𝑠 > 0 and scale 𝑏 > 0, where 𝑠 is a parameter controlling the level of non-Gaussianity. The second-order modular variate Q of a GG distribution is given by Q = 𝑑 𝐺 1/𝑠 where 𝐺 is a Gamma distributed random variable with parameter 1/𝑠 and 𝑏, i.e. 3 and 4 w.r.t. the 𝑆𝑁 𝑅 𝑜𝑢𝑡 . The number of Monte Carlo is set to 1000 iterations. In the simulation, complex centered Generalized Gaussian distributions with 𝑠 = {0.5, 1.5, 2.5} and complex centered 𝑡-distributions with 𝜐 = {1.1, 2, 3} have been used as a true model. In the results one can that the RMSE of the pseudotrue parameter converges to asymptotic estimation performance derived in Section IV. These results confirm the theoretical Note also that the

√

𝑀𝐶 𝑅𝐵 is equal to the √ 𝐶 𝑅𝐵. It is important to underline that the preceding theoretical findings are applicable to all true noise models characterized by Complex Elliptically Symmetric (CES) distributions, not limited to the Gaussian (GG) and 𝑡-distribution cases. As mentioned earlier, a thorough explanation of this phenomenon is rooted in semiparametric theory (refer to [17, Sec. IV.B] and [18, Sec. III.B]), and a comprehensive explanation will be provided in future research. For now, we restrict our discussion to this observation: the equivalence between the MCRB and the CRB holds true only when the parameters of interest parameterize the mean of the observation vectors. Conversely, if some parameter of interest is involved in the covariance matrix of the observations, this equivalence may no longer be valid.

IX. CONCLUSION

This paper focused on the performance evaluation of estimation procedures in nonlinear regression models. In particular, we were interested in analyzing the asymptotic performance of inference algorithms based on the simplistic i.i.d. Gaussian assumption in the presence of correlated and non-Gaussian noise. To this end, the related MCRB has been evaluated and the consistency and efficiency properties of the MMLE/NLLSE investigated. Under a weak condition on the rate of decay of the autocorrelation function (Assumption 1), our results show that: • The MCRB for the parameters of interest in the nonlinear regression model depends of the autocorrelation function of the noise but not on the joint pdf on the noise samples that can then be left unspecified. Moreover, the MMLE/NLLSE is consistent and efficient with respect to the relevant MCRB. • If the noise samples are modeled as zero-mean, i.i.d. CESdistributed (with unspecified density generator) random variables, the MCRB on the parameter of interest equates the CRB derived under i.i.d. Gaussian assumption. This means that the asymptotic performance of Gaussian-based MMLE, i.e. the NLLSE, are not affected by the lack of knowledge of the true non-Gaussian and heavy-tailed noise distribution. Since the i.i.d. Gaussian assumption is widely used in applications, these theoretical results are of great practical interest. Specifically, they implies that, a practitioner can continue to use the i.i.d. Gaussian-based inference procedures without any loss in asymptotic estimation performance even when the noise samples are heavy-tailed, non-Gaussian but still i.i.d. random variables. On the other hands, if the noise samples are correlated heavy-tailed, non-Gaussian random variables, the asymptotic performance of i.i.d. Gaussian-based procedures depends only on the autocorrelation function of the noise process and not on the specific joint pdf of its samples. Our theoretical findings have been then used to investigate the asymptotic performance of Gaussian procedure in time-delay and Doppler estimation for GNSS.

Future works will focus on the possibility to derive lower bounds on the performance of estimation procedures in the presence of an non-perfect knowledge of the nonlinear function characterizing the regression model.

Fig. 1 .

 1 Fig.1. RMSE of the MMLE of the time-delay considering complex centered t-dist. with 𝜐 = 2.5 and two AR processes of order 1 and 6, respectively, to model the noise correlation.

3

 3 for the order 1 process and 𝒑 =[0.5•𝑒 -𝑗2 𝜋•0.4 , 0.6• 𝑒 -𝑗2 𝜋•0.2 , 0.7, 0.4 • 𝑒 𝑗2 𝜋•0.1 , 0.5 • 𝑒 𝑗2 𝜋•0.3 , 0.6 • 𝑒 𝑗2 𝜋•0.35 ] for the order 6 process. The MMLE for the joint estimation of the time-delay and Doppler is defined in (71). The root mean square error (RMSE) results of the MMLE for the parameters of interest 𝜼 𝑇 = [𝜏, 𝑏] are shown in Figs. 1 and 2 w.r.t.

Fig. 2 .

 2 Fig.2. RMSE of the MMLE of the Doppler considering complex centered t-dist. with 𝜐 = 2.5 and two AR processes of order 1 and 6, respectively, to model the noise correlation.

Fig. 3 .

 3 Fig.3. RMSE of the MMLE of the time-delay considering complex centered GG dist. with 𝑠 = {0.5, 1.5, 2.5} and t-dist. with 𝜐 = {1.1, 2, 3}.

Fig. 4 .

 4 Fig. 4. RMSE of the MMLE of the Doppler considering complex centered GG dist. with 𝑠 = {0.5, 1.5, 2.5} and 𝑡-dist. with 𝜐 = {1.1, 2, 3}.

  the marginal pdf 𝑝 𝑛 𝑘 of each sample 𝑛 𝑘 ∼ 𝑝 𝑛 𝑘 is left fully unspecified, • the covariance matrix of n , i.e. 𝚺 ≜ 𝐸 𝑝 n [nn 𝐻 ], has the following Hermitian Toeplitz structure [11, Sec. 1.3]: [𝚺] 𝑘,𝑘+ 𝑗 = [𝚺] 𝑘+1,𝑘+ 𝑗+1 = [𝚺] *

	𝑘+ 𝑗+1,𝑘+1

and only if there exists a positive real number 𝑎 and a real number 𝑥 0 such that | 𝑓 ( 𝑥 ) | ≤ 𝑎𝑔 ( 𝑥 ), ∀ 𝑥 ≥ 𝑥 0 .

•

Given a real-valued function 𝑓 ( 𝑥 ) and a strictly positive real-valued function 𝑔 ( 𝑥 ), 𝑓 ( 𝑥 ) = 𝑜 (𝑔 ( 𝑥 ) ) if for every positive real number 𝑎, there exists a real number 𝑥 0 such that | 𝑓 ( 𝑥 ) | ≤ 𝑎𝑔 ( 𝑥 ), ∀𝑥 ≥ 𝑥 0 .

ACKNOWLEDGMENTS

This work was supported by TéSA.


Gildas Avoine

Xavier Carpent

Diane Leblanc-Albarel
email: diane.leblanc-albarel@irisa.fr

Rainbow Tables: How Far Can CPU Go?

Keywords: Cryptography, Time-Memory Trade-Off (TMTO), Rainbow Table

Rainbow tables are techniques commonly used in computer security to invert oneway functions, for instance to crack passwords, when the domain of definition is reasonably-sized. This article explores the limit on the problem size that can be treated by rainbow tables when the precomputation and the attack phases are both CPU-driven. We conclude that the bottleneck is no longer the memory as it may have been and the precomputation phase seems to have been underestimated so far. In contrast to the usual articles on rainbow tables, we offer a comparison of what can be done on different environments depending on the needs and available computing power of the users.

INTRODUCTION

Context A rainbow table is a data structure used to invert one-way functions -for example hash or encryption functions -when the domain of definition is reasonably-sized.

It is the most used variant of cryptanalytic Time-Memory Trade-Off (TMTO). A TMTO has characteristics between that of bruteforce search (which has high computation cost) and dictionary attack (which has high storage cost). Rainbow tables are particularly useful [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF] when the inversion (or attack ) is repeated many times, (2) when the precomputation is delegated (e.g., a powerful entity precomputes, but a simple laptop carries out the attack), or (3) when there is a small pre-determined window of opportunity for the attack but a possibly long time to prepare for it (lunchtime attack ).

The use of rainbow tables is divided in two phases: precomputation phase (performed once) and attack phase (performed for each inversion). Given a problem of size N (N is the number of possible solutions) and a memory of M T , the complexity of the precomputation phase is linear in N (with a factor varying, depending on parameters and optimizations -see discussion in [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF]) and the complexity of the attack phase is O(N 2 /M 2 T ). This means that a memory anywhere from O( √ N ) (below which the attack becomes slower than brute force) to O(N ) (above which the memory becomes larger than in a dictionary attack) can be used to accelerate the search with respect to brute force. Various papers discuss algorithmic and implementation improvements to both phases: variants [START_REF] Kim | Analysis of the nonperfect table fuzzy rainbow tradeoff[END_REF][START_REF] Hong | Variants of the distinguished point method for cryptanalytic time memory trade-offs[END_REF], attack phase optimizations [START_REF] Avoine | Characterization and improvement of time-memory trade-off based on perfect tables[END_REF][START_REF] Avoine | Analysis of rainbow tables with fingerprints[END_REF][START_REF] Avoine | Heterogeneous rainbow table widths provide faster cryptanalyses[END_REF], storage optimizations [START_REF] Avoine | Optimal storage for rainbow tables[END_REF], and precomputation improvements [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF].

Motivation This article aims to characterize the efficiency of rainbow tables from a practical point of view, and determine which largest N can be realistically addressed by rainbow tables, with these optimizations in mind. This is explored with a CPU-driven computation model, with reasonable assumptions about the computing environment and parameters.

We further aim to gain some insight from this work into the current practical impact of rainbow tables in computer security, and the future research for rainbow tables and associated concepts.

Contribution We define three potential types of environments for the precomputation phase. Each environment corresponds to typical needs and budget with today's technology. For each environment, we evaluate the largest N that can be undertaken and we identify the technological bottleneck that prevents going further. Experiments are provided to illustrate certain theoretical evaluations.

Two environments are super-computers or rather large computers, with at least 100 cores. The one that we call supercomputer is a computer typically belonging to the top-100 worldwide list 3 . The other one, which we call computer corresponds to a computer available for medium sized companies or academic research teams. The third environment, called the cloud environment corresponds to what could be expected on rented machines available in the G. Avoine, X. Carpent, and D. Leblanc-Albarel cloud. The latter environment could be economically interesting if precomputations are performed only once or occasionally. Using these three environments, we conclude what the three main entity types can achieve in TMTOs today. Section 6 provides approaches to improve TMTOs in the future. Knowing the space size that can be easily attacked by TMTOs today also allows the whole community to take measures to protect themselves.

Organization Section 2 summarizes basics on rainbow tables, especially the way the two phases are carried out. Section 3 briefly presents other TMTOs variants and improvements made.

Section 4 introduces the environments and scenarios considered in this article. The performed evaluations and experiments are described in Section 5. Discussions about the technological limits of TMTOs then conclude the article in Section 6.

BACKGROUND

Overview

The purpose of TMTOs is to retrieve a preimage of a given value obtained through a (one-way) function. Usually, this one-way function is denoted h with h : A → B. The searched preimage is denoted by x with x ∈ A, and the given value is denoted by y with y ∈ B. The aim is therefore to retrieve x from y with y = h(x), using precomputed tables.

Various TMTOs algorithms exist. One of the most efficient and commonly used variant is the rainbow tables variant, which is consequently the variant addressed in this paper.

Rainbow tables without any improvement, called classic rainbow tables, are considered in this section. Other variants and improvements are briefly described in Section 3, and more extensive descriptions can be found in, e.g., [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF][START_REF] Oechslin | Making a faster cryptanalytic time-memory trade-off[END_REF][START_REF] Hong | A comparison of cryptanalytic tradeoff algorithms[END_REF][START_REF] Lee | Comparison of perfect table cryptanalytic tradeoff algorithms[END_REF].

Precomputation Phase

Matrix Construction

To construct a rainbow table, a rainbow matrix of m rows and t+1 columns is firstly generated. Each matrix element is denoted x i,j , with x i,j ∈ A, 0 ⩽ i ⩽ t and 0 < j ⩽ m. Each element x i+1,j of a row j, is obtained from the previous element x i,j using a function denoted f i and named hash-reduction function, where:

x i+1,j = f i (x i,j ).
Given the hash function h : A → B, and a reduction function family4 r i : B → A, the hash-reduction function f i is defined by:

f i : A → A x i,j → r i (h(x i,j )) = x i+1,j
The successive application of f i forms a chain. To form a rainbow matrix of m rows and t + 1 columns, m chains are computed by iterating t times the function f i on m arbitrary chosen elements of A. . All the other columns are discarded.

f 0 f 1 x 0,1 -→ x 1,1 -→ x 2,1 . . . . x t,1 f 0 f 1 x 0,1 -→ x 1,2 -→ x 2,2 . . . . x t,2 . . . . . . . . . . . . f 0 f 1 x 0,m -→ x 1,m -→ x 2,m . . . . x t,m

Clean Rainbow Table

During the matrix generation, chains may merge. A merge between 2 chains occurs when all their elements are equal after a given column, i.e., when x i,j ̸ = x k,j and x i,j+1 = x k,j+1 , for 0 ⩽ j < t. Chain merging is possible because |B| > |A|.

To make the attack phase as efficient as possible, the table is trimmed to only keep rows that have not merged, i.e., chains with unique end points -this is called a clean or perfect table .  As the attack phase should be as fast as possible, clean rainbow tables are usually used. In what follows, tables will be considered clean, m 0 will denote the number of start points used to generate the matrix, and m t will denote the number of points remaining in column t after the cleaning process.

Non-Maximal Table

The number of start points, m 0 , is chosen to be as high as practically possible for a given t. When m 0 = N the table is said to be maximal. In that case m max t points remain in the last column after cleaning the table. m max t is defined in Theorem 2.1 from [START_REF] Avoine | Characterization and improvement of time-memory trade-off based on perfect tables[END_REF].

Theorem 2.1. Given t and a sufficiently large N, the expected maximum number of chains per clean rainbow table is:

m max t ≈ 2N t + 2 .
The higher the desired number of points remaining in column t, the higher the proportion of chains that need to be thrown out during the cleaning due to duplicate end points. Therefore, besides being very costly, taking m 0 = N is not worth the cost compared to the number of remaining chains at the end of the generation. Thus, in practice, non-maximal tables are generated with a number of start points much lower than N .

Given a number of start points m 0 < N , the maximality factor α with 0 < α < 1, is defined such that:

m t = α m max t
It is then possible to define a ratio r of chains to be generated in order to obtain a targeted

m t =
r ≈ α 1 -α A table created
with a typical ratio of r = 20 chains computed per chain kept for instance, results in m t being at about 95% of its maximal value.

Success Probability

A single clean table covers much of A, but not all of it (very close to 86.47% for maximal tables). Therefore, several independently-computed tables are used in the attack phase, reaching a success probability arbitrarily close to 1. For a single table with m t end points, the success probability P (t) is given by Equation (1), adapted from [START_REF] Oechslin | Making a faster cryptanalytic time-memory trade-off[END_REF].

P (t) = 1 -(1 - m t N ) t . (1)
Using ℓ tables the success probability P (t, ℓ) is given by Equation ( 2) from [START_REF] Oechslin | Making a faster cryptanalytic time-memory trade-off[END_REF].

P (t, ℓ) = 1 -(1 -P c (t)) ℓ . (2)

Attack phase

Once the precomputation phase is completed, the attack phase can be performed. The latter phase uses the precomputed rainbow tables to find a preimage of a given y ∈ B. The following process is iterated until an answer is found or after t iterations (in which case the process fails): at iteration i, a chain of length i is computed: f t (. . . f t-i+1 (r t-i (y)) . . . ). Its end point is then compared to the end points in the table. If a match is found, the corresponding chain of the matrix is rebuilt from its start point, up to the column t -i -1:

x t-i-1,j = f t-i-1 (. . . f 1 (x 1,j ) . . . ). If h(x t-i-1 ) = y
then the attack is completed and successful. If the two hashes are not equal, or if no match was found, the process proceeds to the next iteration. If the attack reaches t iterations, it fails. Proposition 2.2 adapted from [START_REF] Avoine | Characterization and improvement of time-memory trade-off based on perfect tables[END_REF], defines the number of operations needed to perform a search in the column c, i.e., after t -c iterations. The number of unique points remaining in column i is denoted m i and can be easily computed as demonstrated in [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF].

Proposition 2.2. After t -c iterations, the average number of hash operations C c needed to perform a search is:

C c = t -c t i=c 1 - m i N Proof.
A proof for maximal tables can be found in [START_REF] Avoine | Characterization and improvement of time-memory trade-off based on perfect tables[END_REF].

The result can be extrapolated to non-maximal rainbow tables.

The average attack time using rainbow tables can be deduced from Proposition 2.2 and is given by Theorem 2.2 Theorem 2.2. The average number of hash operations T required to perform an attack using rainbow tables, given a search space of size N and ℓ tables is:

T = ℓ t i=1   m t N 1 - m t N ℓ(i-1) c j=1 C t-j+1   +e -2ℓ ℓ t i=1 C i .
Proof. As Proposition 2.2, this theorem is derived from formulas given in [START_REF] Avoine | Characterization and improvement of time-memory trade-off based on perfect tables[END_REF].

During the attack phase, tables can be loaded in RAM or on disks (preferably SSDs). In the latter case, several approaches can be taken (see [START_REF] Avoine | How to handle rainbow tables with external memory[END_REF][START_REF] Kim | Analysis of the rainbow tradeoff algorithm used in practice[END_REF]).

In [START_REF] Avoine | How to handle rainbow tables with external memory[END_REF], the authors show that the attack time is not significantly impacted by the use of secondary memory (SSD, HDD etc.) rather than RAM. Furthermore, as discussed in Section 5 and 6, performing the attack phase on secondary memory rather than on RAM would not change the conclusions of this paper.

However, an attack on secondary memory is more complex. Thus, for the sake of clarity, in this paper, we consider attack phase on RAM.

The attack phase benefited from various algorithmic improvements that significantly improved its performances. All told, a fully-fledged and well-parameterized implementation can be from 5 to 10 times faster5 than a naive implementation. While important in practice, they do not change the order of magnitude of the attack phase cost. A brief summary is nonetheless presented in Section 3, but for clarity, rainbow tables without improvement on the attack will be used in this paper.

RELATED WORK

This section provides an overview of TMTO variants and improvements published so far, along with references that detail each of them.

TMTOs Variants

Hellman Table . The earliest TMTO was introduced by Martin Hellman in 1980 [START_REF] Hellman | A cryptanalytic time-memory trade-off[END_REF]. This work differs from the rainbow variant on the choice of the reduction function: only one reduction function is used in Hellman's case instead of a reduction function per column in the rainbow variant. At each iteration performed during the attack phase, only one additional application of the reduction function is needed, which speeds up considerably this step.

However, the coverage of Hellman's tables is very poor due to the large number of collisions it generates, which implies using a very large number of tables (typically t tables) to obtain a suitable coverage. The greater the number of tables used, the longer the attack time. In addition, cleaning Hellman's table is very difficult and not necessarily worth it. In the end, Hellman's tables are less efficient than rainbow tables and, therefore, are non longer used.

Distinguished Point (DP) Tables. DP variant is based on Hellman's tables, but instead of having chains of length t, chains are computed until reaching a socalled distinguished point (typically, points that have at least d bits at 0). DP tables have thus chains with variable lengths. End points are all distinguished points which make the table easy to clean. As chains of DP tables do not have predictable sizes, it alters both the precomputation phase and the attack phase. In the end, it has been established in [START_REF] Hong | Variants of the distinguished point method for cryptanalytic time memory trade-offs[END_REF][START_REF] Lee | Comparison of perfect table cryptanalytic tradeoff algorithms[END_REF], that DP tables are less efficient than rainbow tables.

Fuzzy Rainbow Tables. This variant is a trade-off between rainbow tables and DP tables. Fuzzy rainbow tables are extensively studied in [START_REF] Kim | Analysis of the nonperfect table fuzzy rainbow tradeoff[END_REF]. In this variant, each chain is built by concatenating DP chains. As shown in [START_REF] Hong | A comparison of cryptanalytic tradeoff algorithms[END_REF], except for very low success probabilities and, therefore, uninteresting cases in practice, fuzzy rainbow tables are less efficient in precomputation and attack than rainbow tables.

Precomputation Phase Improvement

Until recently, most of the improvements on TMTOs focused on the memory or time needed for the attack. In 2021, [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF] shows that the precomputation time can also be improved to cover larger spaces. The filtration method have thus been introduced in [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF], the method consists in cleaning the matrix as the generation progresses rather than at the end and allows dividing the precomputation time by 6 without increasing the attack time.

Attack Phase Improvements

Checkpoints. Checkpoints have been introduced in [START_REF] Avoine | Characterization and improvement of time-memory trade-off based on perfect tables[END_REF]. The concept consists in storing additional information -so call checkpoints -on each chain, such that this information allows reducing the attack phase time. The checkpoints, are stored along with the start points and end points. During the attack phase, when the attack chain matches an end point of the table, the checkpoints of the built attack chain are compared to the checkpoints of the matching chain. It at least a checkpoint differs, then the match is a false alarm and it is consequently useless to rebuild the chain of the table from its starting point. With optimal parameters, rainbow tables with checkpoints are faster than the classical case but require slightly more memory.

Fingerprints. This improvement has been introduced in [START_REF] Avoine | Analysis of rainbow tables with fingerprints[END_REF]. Its consists of applying in an efficient way the checkpoints improvement with a memory improvement called truncated end points. the principle is to store end points that have been truncated to reduce the memory needs for the attacks phase. Truncated end points alone increase the attack time significantly but used combined with checkpoints in the fingerprints method, they allow a speedup of about 2 compared to the classical rainbow table.

Heterogeneous Tables. This variant has been introduced in [START_REF] Avoine | Heterogeneous rainbow table widths provide faster cryptanalyses[END_REF]. It consists in using rainbow tables with different lengths. This implies that searching in the shorter tables takes less time than searching in classic tables but searching in the longer ones takes more time. This variant does not require more memory than the classical case. By wisely performing the attack phase (favoring searches in the shortest tables), this variant reduces the average attack time by up to 40% without increasing the memory. The counterpart is that the worst case is slightly longer than with classic tables.

ENVIRONMENTS AND SCENARIOS CONSIDERED

Context

In practice, the entities that perform TMTOs have different needs, purposes, and resources (e.g., available memory, price, time available, etc.). Resources are, in addition, not the same for the attack and precomputation phases.

It is therefore necessary to define the context in which each phase is performed. A company, for instance, does not have the same resources than a nation state. The memory available for the precomputation must be defined as well as the one for the attack. Other variables such as the available time should also be specified in advance.

To define the context in which a TMTO is performed, we will use the notion of environments and scenarios. An environment corresponds to the material resources available (RAM, number of cores, CPU performances etc.). A scenario corresponds to the non-material resources available (Time or money).

In this paper, we consider different environments, each of them described in 4.2 (precomputation) and in 4.3 (attack phase). These environments aim to represent different entities. For each environment, we have considered different scenarios depending on available time or money. 2) computer, which corresponds to a 128-core computer; and (3) cloud, which consists of rented computing units from on the main cloud platforms, e.g., AWS, Azure, or GCP. For each environment, the number of cores, the number of hashes/second/core, and the available memory (RAM) are provided for supercomputer and computer in Table 1 and in Table 2 for cloud.

Precomputation phase

Typically, supercomputer represents the computing power of a governmental agency, computer might be a computer owned by a university or a small to mediumsized business, and cloud illustrates an entity that rents commercial computing units in order to precompute rainbow tables.

For each environment, three scenarios are considered. As shown in Table 1, for supercomputer and computer, the scenarios depends on the available time for the precomputation phase: 1 year, 1 month, or 1 week. For cloud however, the precomputation phase is bounded to 1 month, and as shown in Table 2, the environments are defined by the budget assigned to the precomputation phase according to the different scenarios: 1 000 000 USD, 100 000 USD, or 10 000 USD

Attack phase

For the attack phase, we consider environments with a single core 7 . This makes the analysis and description of environments a little simpler, and as discussed in Section 6, does not change the overal situation.

The number of hashes per second n t , on the attack core is fixed for each environment n t = 11 000 000 8 . The idea is that the attacker usually has limited resources for the attack phase and cannot benefit from a supercomputer for this phase.

The only attribute that has a bearing on the efficiency of the attack is the size of the memory available. The different values for the attack memory M T are chosen to correspond to realistic, practical cases and influence the parameter t (given our assumption of ℓ = 4 tables). Memory available (but not necessarily used) are presented in Table 3.

The memory available for this phase must remain relatively small. To correspond to cases encountered in practice, it must be small compared to the memory available for the precomputation phase. The memory available for the attack is therefore adapted to the scenarios considered and to the environments used for the precomputation phase.

For each environment, scenarios depend on the time or money invested in the precomputation phase. The aim is to perform the attack as quickly as possible given the tables generated during the corresponding precomputation scenarios: 1 year, 1 month, 1 week of precomputation or 1M, 100K, 10K, USD invested.

EVALUATION OF THE MAXIMUM PROBLEM SIZE

Methodology

The maximum problem size N that can be addressed by a CPU-based TMTO in a given time-frame or budget can be accurately evaluated from the analytical formulas provided in [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF] (precomputation time) and [START_REF] Avoine | How to handle rainbow tables with external memory[END_REF][START_REF] Kim | Analysis of the rainbow tradeoff algorithm used in practice[END_REF] (attack time).

Precomputation Phase

The minimum precomputation time P min in seconds is given by Equation (3) derived from Theorem 2 in [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF].

G. Avoine, X. Carpent, and D. Leblanc-Albarel

The total precomputation time P in seconds is given by Equation ( 4) derived from Equations ( 5) to ( 9) in [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF], with the same speed values (called v o and v c ) used in the first environment of [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF].

P min = 2N n h v h ln(1 + r) (3)
P = M ax 1 n h v h a+1 i=1 m ci-1 (c i -c i-1 ); 1 v f n f a+1 i=1 m ci-1 + v o a+1 i=1 m ci-1 + v c n h a+1 i=1 m ci-1 (4)
The following parameters are used in these equations: a is the number of filters used during the precomputation phase; each c i is the columns of the i-th filter, with c 0 = 0 and c a = t; finally, the values n h , v h , n f and v f denote the number of computation nodes, the number of hashes per second per core, the number of filtration nodes and the number of filtration per second per core respectively. As in [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF], a single filtration core is sufficient, thus n f = 1. Values of n h and v h according to the environment and scenario are given in Table 1 and2. Values for v f are given in Table 4.

For computer and cloud environments the number of hashes per second corresponds to typical hashes value of an environment similar to computer. The filtration speed corresponds to filtration speed measured on an environment similar to computer. For supercomputer environment hashes and filtrations per second have been computed from typical FLOPS numbers on this kind of environment.

Environment

Hashes Several fixed parameters are defined as follows: the number of tables generated is ℓ = 4, the factor r = 20 is used (so m 0 = 20m max t is considered at the beginning of precomputation). Since environments considered in this paper are very close to the first environment of [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF], identical values of v c and v o are used for the estimations, i.e., v c = 0 and v o = 1.37.10 -10 .

Memory Used

The memory M T used by the ℓ = 4 tables is given by Equation ( 5), the maximum RAM needed for the precomputation phase M P is given by Equation ( 6) with m c1 the number of unique points in the column of the first filter. The factor 3 is explained by the fact that the hash table used for filtration has a load factor λ = 1.5 and that both start points and end points are stored.

As presented in [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF], filters are placed to minimize P . M T = 2ℓm t log 2 (N )

(5)

M P = 3m c1 log 2 (N ) (6)
The parameter t is chosen for each time frame (year, month, or week), determined from M T , and according to the budget (1M, 100K, or 10K USD) for the precomputation phase. As presented, in Proposition 1 from [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF], t impacts M P . Therefore, t influences the memory of the two phases. If t is too high, the attack phase will be too slow. On the other hand, if t is too small, the memory needed for precomputation and especially for storing tables may be too large.

Attack Phase

As only one core is used for the attack phase regardless of the environment used for the precomputation, the attack time can be well estimated by dividing result of Theorem 2.2 by the number of hashes per second n t performed by the CPU used for the attack. Therefore, the attack time T RAM , is given by Equation ( 7)

T RAM = T n t , (7)
with T the number of hashes needed to perform the attack given in Theorem 2.2. As the attack should be performed by an average computer, we have chosen to use n t = 11 000 000 as the number of hash per second used for the attack phase.

Results

The results of the evaluation are provided in Table 5, where 9 configurations are presented according to the environments and scenarios considered.

6. DISCUSSION 6.1. Noteworthy observations [START_REF] Avoine | Precomputation for rainbow tables has never been so fast[END_REF] The cost P and the theoretical lower bound P min are very close in most cases. In the cases where the gap between P min and P is more significant (e.g., supercomputer environment), converting some computation nodes to filtering nodes would reduce the gap significantly, as it is mainly explained by a lack of filtration speed.

Therefore, unless a significant algorithmic change is introduced, improving the efficiency of the precomputation phase would have no noticeable impact on the domain size that can be realistically attacked.

(2) For all environments considered, the largest domain sizes on which tables can be precomputed in each scenario lead to tables that are practical for the attack phase (i.e., reasonably small in memory, and reasonably fast to execute a search). In other words, within the context established in this paper on the environments, scenarios, technology, and algorithms, The abbreviations "d", "h", "m", and "s" respectively denote "days", "hours", "minutes", and "seconds".

TABLE 5: Evaluation of the maximum problem size the precomputation phase is the bottleneck to using rainbow tables on large domains. The slight overhead due to performing the attack phase on secondary memory (HDD or SSD) would thus not change this conclusion.

(3) The memory M P needed to perform the precomputation is not a bottleneck in the environments and scenarios analyzed in this paper.

For instance, in the case of computer environment, the RAM available for the precomputation phase is limited compared to the ones available in the other environments.

This limited memory available for precomputation does not significantly impact the attack times of the different scenarios as they remain fairly fast.

Depending on the variants and improvements considered, it might be expected that the memory for precomputation would be limiting, although in theory, this limitation would not necessarily exist.

Therefore, the cases considered in this paper demonstrate that the memory available for precomputation is not a limiting factor. Indeed, the attack times are fast when using realistic environments and scenarios with limited memory for precomputation.

Conclusion

Our main observation from Table 5 is that the precomputation cost P is, today, the bottleneck of the TMTO. Indeed, the attack time remains reasonably low, even for the bigger space sizes considered. The memory needed for the precomputation and the attack is affordable (for the corresponding entity that performs it). Therefore, it is the precomputation time P , and more precisely, the computing power of the adversary (number of computation nodes and/or number of hashes per second) dedicated to the precomputation that limits an increase in N .

Using the state-of-the-art precomputation algorithm, there is little room for improvements on the precomputation cost P . Indeed, the P is very close to the theoretical lower bound P min . Going further may require forgoing the CPU technology for a more efficient one, which could currently be GPU or FPGA. While these technologies have constraints in how they are put to use, they do operate several orders of magnitude faster than typical CPUs. Several articles, e.g., [START_REF] Lu | Time-memory trade-off attack on the GSM A5/1 stream cipher using commodity GPGPU[END_REF] treat this problem and websites propose implementations or programs 9to purchase or generate rainbow tables on GPU. Other contributions, e.g., [START_REF] Standaert | A time-memory tradeoff using distinguished points: New analysis & FPGA results[END_REF][START_REF] Mentens | Time-memory trade-off attack on FPGA platforms: UNIX password cracking[END_REF] focus on FPGAbased TMTOs. However, these papers address the problem with relatively small domains N (in the 2 40 to 2 50 range), deprecated one-way functions, or older GPU/FPGA models. Furthermore, they do not use recent improvements in TMTOs, e.g., filtration. Therefore, it is not possible to compare these results obtained on GPU/FGPA with the results we obtained in this paper. A deeper look into efficient GPUor FPGA-based precomputation represents interesting future work that could evaluate these technologies capabilities to deal with larger TMTOs.

In summary, our work helps quantify the vulnerable spaces to TMTOs performed on CPUs. In addition, although the precomputation phase has not been as widely studied as the attack phase, it nevertheless seems to be the main bottleneck preventing larger TMTObased attacks today.

Our conclusions are of course dependent on the assumptions made in the premise of this paper, i.e. specific to CPU-based attacks, under some reasonable computing environments and precomputation scenarios, with specific attack phase context in mind, etc. These assumptions were established with care, and a large G. Avoine, X. Carpent, and D. Leblanc-Albarel spectrum of conditions were considered. Nevertheless, a significant change in hardware technology, or a profound algorithmic modification (especially targeting the precomputation), or a notable deviation of the assumptions, could yet increase the domain size vulnerable to attacks.

The precomputation phase being a bottleneck in our observations may likewise be challenged by such changes. Indeed, the precomputation cost is by nature linear in N , whereas the attack phase is by nature quadratic in N . A situation where the attack phase is the bottleneck is therefore not far-fetched, depending on the development of technologies and research in this field.

FIGURE 1 :

 1 FIGURE 1: Rainbow matrix

  αm max t number of chains in the final clean table. Proposition 2.1 introduced in [1] defines this ratio r. Proposition 2.1. With a target of m t = α m max t unique end points, m 0 = rm max t chains need to be generated to construct a rainbow table, with:

TABLE 2 :

 2 6 . Precomputing cloud environments according to the scenarios For each precomputing environment and scenario, ℓ = 4 tables are used, and parameters α = 0.95 and thus r = 20 are chosen to reach a success probability larger than 99.95%. The aim is to maximize N while keeping the attack phase affordable.

	Scenarios	cloud Environments #Cores Hashes/Sec/Core RAM (TB)
	1M USD	13 055	11 000 000	256
	100K USD	1 279	11 000 000	20
	10K USD	127	11 000 000	2

TABLE 3 :

 3 Memory available for the attack according to the scenarios and environments

TABLE 4 :

 4 Hashes and filtrations per second per core

		/Sec/Core Filtration/Sec/Core
	supercomputer	25 000 000	51 270 585
	computer	11 000 000	15 949 709
	cloud	11 000 000	15 949 709

https://top500.org/

The purpose of a reduction function is to map a hash to an arbitrary element of the input space A. A notable example is r i (y) → (y + i) mod N for A = {0, 1, . . . , N -1}.

The exact speed-up depends on many factors and parameters. This window is a conservative estimate, based on typical scenarios and configurations.

The computing characteristics shown in Table2have been obtained by simulating on these costs on AWS (cluster EC2).

Note that a single core is considered to provide a reference value, but the attack phase can be easily parallelized to operate on several cores.

This corresponds to the number of SHA256 hashes per second measured on a AMD EPYC 7742 3.2 GHz processor.

For instance, https://www.cryptohaze.com/ offers a GPUbased rainbow cracker, and http://project-rainbowcrack.com/ has implementations of rainbow tables on GPU.

ACKNOWLEDGMENTS

A significant portion of this work was done while Xavier Carpent worked at COSIC/KULeuven (BELGIUM).