Luc Le Magoarou

Stéphane Paquelet

Performance of MIMO channel estimation with a physical model

Keywords: Channel estimation, physical model, MIMO

Channel estimation is challenging in multi-antenna communication systems, because of the large number of parameters to estimate. One way of facilitating this task is to use a physical model describing the multiple paths constituting the channel, in the hope of reducing the number of unknowns in the problem. The achievable performance of estimation using this kind of physical model is studied theoretically. It is found that adjusting the number of estimated paths leads to a bias-variance tradeoff which is characterized. Moreover, computing the Fisher information matrix of the model allows to identify orthogonal parameters, ultimately leading to fast and asymptotically optimal algorithms as a byproduct.

I. INTRODUCTION

M ULTIPLE-input multiple-output (MIMO) communica- tion systems allow for a dramatic increase in channel capacity, adding space to the classical time and frequency dimensions [START_REF] Telatar | Capacity of multi-antenna gaussian channels[END_REF], [START_REF] Tse | Fundamentals of wireless communication[END_REF]. This is done by using several antennas at the transmitter (N t ) and at the receiver (N r ). The capacity of MIMO systems is maximized if the channel state is perfectly known at both ends of the link, which requires a channel estimation step.

Channel estimation is deeply impacted by the transition from single antenna to MIMO systems. Indeed, it amounts to determine a complex gain for each transmit/receive antenna pair, for each subcarrier. If N f subcarriers are used, the number of real parameters to estimate is thus 2N r N t N f , which may be very large for massive MIMO systems, i.e. systems with up to several hundreds of antennas [START_REF] Rusek | Scaling up mimo: Opportunities and challenges with very large arrays[END_REF], [START_REF] Larsson | Massive mimo for next generation wireless systems[END_REF].

In a massive MIMO context, the classical least squares (LS) channel estimator is not adapted because the high dimensionality of the parameter space leads to an ill-posed problem. Therefore, in order to add prior information, it has been proposed classically to use Bayesian estimation and thus model the channel matrix as random with a known distribution, giving rise to estimators such as the linear minimum mean square error (LMMSE) [START_REF] Björnson | Massive mimo networks: Spectral, energy, and hardware efficiency[END_REF], [START_REF] Bazzi | Downlink training sequence design for fdd multiuser massive mimo systems[END_REF]. Another possibility is to use a parametric model based on the physics of wave propagation, in which the channel is expressed as a sum of p paths [START_REF] Sayeed | Deconstructing multiantenna fading channels[END_REF]. Whereas LS and LMMSE estimators have been studied extensively in terms of optimal training sequences and performance This work has been performed in the framework of the Horizon 2020 project ONE5G (ICT-760809) receiving funds from the European Union. The authors would like to acknowledge the contributions of their colleagues in the project, although the views expressed in this contribution are those of the authors and do not necessarily represent the project.

Luc Le Magoarou and Stéphane Paquelet are both with b<>com, Rennes, France. Contact addresses: luc.lemagoarou@b-com.com, stephane.paquelet@b-com.com. [START_REF] Biguesh | Training-based mimo channel estimation: a study of estimator tradeoffs and optimal training signals[END_REF], a similar study is still lacking for channel estimators based on a physical model. Contributions and organization. In this paper, the performance of MIMO wideband channel estimators based on a physical model is studied. The considered model is presented and the estimation error is decomposed into a bias and a variance terms in section II. Each term is then analyzed theoretically, beginning with the variance term in section III and the bias term in section IV, highlighting a bias-variance tradeoff piloted by the number of paths p considered by the model. The theoretical analysis also allows to nicely interpret the channel estimation error and leads to the design of a computationally efficient channel estimation algorithm with optimality properties. Finally, the mathematical developments are assessed empirically in section V, justifying the use of such physical models and showing the interest of the designed computationally efficient algorithm.

Note that this paper is based on some of our previous work [START_REF] Magoarou | Parametric channel estimation for massive MIMO[END_REF], [START_REF]Bias-variance tradeoff in MIMO channel estimation[END_REF], but significantly extends it in several ways. First, the wideband MIMO channel is considered here whereas only the narrowband MIMO channel was considered in [START_REF] Magoarou | Parametric channel estimation for massive MIMO[END_REF], [START_REF]Bias-variance tradeoff in MIMO channel estimation[END_REF]. Second, the assumptions required for the study of the variance are relaxed, and a much more extensive set of experiments is performed.

II. PROBLEM FORMULATION

Notations. Matrices and vectors are denoted by bold uppercase and lower-case letters: A and a (except 3D "spatial" vectors that are denoted -→ a ); the ith column of a matrix A by a i ; its entry at the ith line and jth column by a ij . A matrix transpose, conjugate and transconjugate is denoted by A T , A * and A H respectively. The trace of a linear transformation represented by A is denoted Tr(A). For matrices A and B, A B means that A -B is positive semidefinite. The linear span of a set of vectors A is denoted: span(A). The Kronecker product and vectorization operators are denoted by ⊗ and vec(•) respectively, and the Hadamard (entry-wise) product by . The identity matrix is denoted by Id. CN (µ, Σ) denotes the standard complex gaussian distribution with mean µ and covariance Σ. E(•) denotes expectation and cov(•) the covariance of its argument. S 2 denotes the 3D unit sphere.

A. Physical channel model

Let us study a wideband block fading MIMO channel whose features are:

• N t transmit antennas located at the positions -→ a t,1 , . . . , --→ a t,Nt with respect to the centroid of the transmit antenna array, whose radius is denoted R t max j -→ a t,j 2 .

• N r receive antennas located at the positions --→ a r,1 , . . . , ---→ a r,Nr with respect to the centroid of the receive antenna array, whose radius is denoted

R t max i -→ a r,i 2 . • N f subcarriers at the frequencies f c + f 1 , . . . , f c + f N f ,
where f c is the central frequency ( k f k = 0), and the bandwidth is denoted B f N f -f 1 . • P propagation paths, where β l is the complex gain of the l-th path, τ l is its delay, -→ u r,l is its direction of arrival (DoA) and -→ u t,l is its direction of departure (DoD). Making the plane wave assumption, and further assuming RrB c

1 and RtB c 1 so that the phase difference due to the antenna positions is the same for all subcarriers, the channel between the j-th transmit antenna and the i-th receive antenna at the k-th subcarrier can be expressed

h ijk = P l=1 β l e -j2π[ 1 λ --→ ar,i. --→ u r,l -1 λ --→ at,j . --→ u t,l +f k τ l ] , (1)
where λ c fc is the wavelength at the central frequency. The used system representation is illustrated in figure 1. It is possible to express the whole channel in a compact form by introducing the steering vectors,

e t ( - → u ) 1 √ Nt e -j 2π λ --→ at,1. - → u , . . . , e -j 2π λ ---→ a t,N t . - → u T ∈ C Nt , e r ( - → u ) 1 √ Nr e -j 2π λ --→ ar,1. - → u , . . . , e -j 2π λ ---→ a r,Nr . - → u T ∈ C Nr ,
and the delay vector,

e f (τ ) 1 √ N f e -j2πf1τ , . . . , e -j2πf N f τ T ∈ C N f .
Indeed, defining the characteristic vectors,

e( -→ u r , - → u t , τ ) e f (τ ) ⊗ e t ( - → u t ) * ⊗ e r ( -→ u r ) ∈ C NrNtN f ,
the channel can be expressed as a simple linear combination:

h = N r N t N f P l=1 β l e( -→ u r,l , -→ u t,l , τ l ). ( 2
)
This expression of the channel is very general and encompasses many cases encountered in practice. For example, a narrowband (single carrier) channel corresponds to take N f = 1 and e f (τ ) = 1, ∀τ , and a multiple-input singleoutput (MISO) channel corresponds to take N r = 1 and e r ( -→ u r ) = 1, ∀ -→ u r . Most MIMO channel simulators [START_REF] Sun | A novel millimeterwave channel simulator and applications for 5g wireless communications[END_REF]- [START_REF]Study on channel model for frequencies from 0.5 to 100 ghz[END_REF] are based on a similar parameterization which can be described by a vector φ ∈ R N φ containing the physical parameters:

φ (β l , -→ u r,l , -→ u t,l , τ l ) P l=1 . ( 3
)
The channel is thus described with N φ = 7P real numbers (two for the complex gain, two for each direction and one for the delay of each path), with P taking values up to several hundreds in the most widepsread models [START_REF] Sun | A novel millimeterwave channel simulator and applications for 5g wireless communications[END_REF]- [START_REF]Study on channel model for frequencies from 0.5 to 100 ghz[END_REF].

B. Channel estimation

Observations. Let us consider that in order to carry out channel estimation, N m noisy linear measurements of the channel are taken. The obtained observations are expressed

y = Mh + n, (4)
where M ∈ C Nm×NrNtN f is the observation matrix and n is the received noise, which is assumed complex gaussian: n ∼ CN (0, Σ). Once again, this way of expressing the observations is very general. For example in the case of an hybrid system [START_REF] Ayach | Spatially sparse precoding in millimeter wave mimo systems[END_REF]- [START_REF] Sayeed | Millimeter-Wave MIMO Transceivers: Theory, Design and Implementation[END_REF], if the channel is observed through N c analog combiners given by the combining matrix W ∈ C Nc×Nr , and the same training sequence of duration N s given by the matrix X ∈ C Nt×Ns is sent on N ps subcarriers whose positions are given by a row-sampled identity matrix F ∈ {0, 1} Nps×N f , then the observation matrix would take the form M = F⊗X T ⊗W H ∈ C NcNsNps×NrNtN f . On the other hand, if every subcarrier is used to send pilots and the output of every receive antenna is observed, the observation matrix would take the form

M = Id N f ⊗ X T ⊗ Id Nr ∈ C NrNsN f ×NrNtN f .
Objective. Channel estimation aims at retrieving h from the observation of y, knowing the observation matrix M and the distribution of the noise vector n. The channel estimator is denoted ĥ, and is assessed by its mean squared error MSE( ĥ) E hĥ 2 2

(5)

= h -E[ ĥ] 2 2 + E ĥ -E[ ĥ] 2 2 , (6)
where the second line is the well-known bias-variance decomposition [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. At first sight, the search space of channel estimation is of dimension 2N r N t N f (number of real numbers needed to describe the channel vector), which may be very large in massive MIMO systems (up to several thousands). For this reason, classical estimation methods such as the least squares (LS) may not be appropriate. In order to overcome this limitation, some information about the channel has to be used to regularize the problem. For example, the channel can be considered as a random vector whose distribution is known, yielding Bayesian channel estimation [START_REF] Björnson | Massive mimo networks: Spectral, energy, and hardware efficiency[END_REF], [START_REF] Bazzi | Downlink training sequence design for fdd multiuser massive mimo systems[END_REF]. Another possibility is to use a physical channel model, as the one presented in section II-A. In that case, according to (2), the dimension of the search space is equal to 7P , which can also be very large. However in practice, the number of estimated paths p is much smaller than P , so that there are only N θ = 7p parameters to estimate. In that case the estimate belongs by construction to the set M p of vectors that can be expressed by a sum of at most p virtual paths. This set is called model set hereafter and is expressed

M p = a ∈ C NrNtN f | a = p n=1 γ n e( --→ w r,n , --→ w t,n , ξ n ) ,
with α n ∈ C, --→ w r,n ∈ S 2 , --→ w t,n ∈ S 2 and ξ n ∈ R + . Such sets obey the inclusion relation M q ⊂ M q+1 . The best element of the model set is denoted

h Mp argmin g∈Mp h -g 2 . ( 7
)
It is simply the projection of the channel onto the model set M p . By construction, this projection can be written

h Mp = N r N t N f p n=1 γ n e( --→ w r,n , --→ w t,n , ξ n ), (8)
which gives optimal values for the estimated parameters, or virtual parameters

θ (γ n , --→ w r,n , --→ w t,n , ξ n ) p n=1 ∈ R N θ . (9)
The actual estimate will be denoted

ĥ = N r N t N f p m=1 α m e( --→ v r,m , --→ v t,m , ζ m ), (10)
so that the estimated parameters are

θ (α m , --→ v r,m , --→ v t,m , ζ m ) p m=1 ∈ R N θ . (11)
In summary, the channel depends on the 7P physical parameters φ given in [START_REF] Rusek | Scaling up mimo: Opportunities and challenges with very large arrays[END_REF], but the considered model depends on only 7p parameters whose optimal values make up the virtual parameters θ defined in [START_REF] Magoarou | Parametric channel estimation for massive MIMO[END_REF], and the channel estimate depends on the actually estimated parameters θ given in [START_REF] Sun | A novel millimeterwave channel simulator and applications for 5g wireless communications[END_REF]. The central question of this paper is:

How well can an estimator that yields estimates belonging to M p perform with a channel taking the form (2)?

MSE decomposition.

In order to answer this fundamental question, let us first notice that an estimate taking the form (10) cannot be better than the best element of the model set, that is hĥ

2 ≥ h -h Mp 2 .
The study of this paper focuses on estimators whose expectation E[ ĥ] is equal to the projection h Mp , so that the MSE can be decomposed as in [START_REF] Bazzi | Downlink training sequence design for fdd multiuser massive mimo systems[END_REF]:

MSE( ĥ) = h -h Mp 2 2 + E ĥ -h Mp 2 2 . ( 12
)
According to this assumption, the bias of the estimator is identified to the model bias (distance of the true channel h to the model M p ), and the variance is computed with respect to the projection h Mp . In the next sections, the two terms of this decomposition are analyzed separately.

III. VARIANCE ANALYSIS

In this section, the variance term of the mean squared error [START_REF] Jaeckel | Quadriga: A 3-d multi-cell channel model with time evolution for enabling virtual field trials[END_REF] is bounded using the Cramér-Rao Bound (CRB) [START_REF] Rao | Information and the accuracy attainable in the estimation of statistical parameters[END_REF], [START_REF] Cramér | Mathematical Methods of Statistics[END_REF], which is valid for any unbiased estimator. The case of a linear channel parameterization is first evoked in section III-A because of its remarkable similarities with the studied problem. The bound for the studied model is then derived in section III-B, and computationally efficient estimation algorithms are deduced of the particular form of the Fisher information matrix in sections III-C and III-D.

A. The linear case

The projected channel h Mp is a nonlinear function of the virtual parameters θ. Let us first look at a simplified problem in which the channel to estimate h l is linearly linked to the parameter vector θ l ∈ R N θ via the equation h l Dθ l . In that hypothetical case, the observations would read

y l = MDθ l + n.
The least-squares estimate of the parameter vector (without constraining it to be real) is then

θl = (D H M H MD) -1 D H M H y l = θ l + (D H M H MD) -1 D H M H n,
and a channel estimate can be obtained as

ĥl = D θl = Dθ l h l + D(D H M H MD) -1 D H M H n n .
This is an unbiased estimator of h l . Assuming a white noise n ∼ CN (0, σ 2 Id), the covariance of n is given by

E[n n H ] = σ 2 D(D H M H MD) -1 D H .
The variance of the estimator ĥl is thus given by

E ĥl -h l 2 2 = σ 2 Tr[D(D H M H MD) -1 D H ] = σ 2 N θ M 2 2 if D H M H M M 2 2 D = D H D.
According to the second line of the above equation, provided the observation matrix fulfills some condition (being conformal on the column space of D), the variance is proportional to the number of parameters and the noise variance, and inversely proportional to the observation matrix squared norm. As will be highlighted in the next subsection, the lower bound on the variance of ĥ is surprisingly similar to this expression, despite the non-linearities involved.

B. Cramér-Rao bound

Let us go back to the original problem and consider the virtual parameters θ defined in (9) as true parameters and θ defined in [START_REF] Sun | A novel millimeterwave channel simulator and applications for 5g wireless communications[END_REF] as their unbiased estimates. Applying the complex CRB [START_REF] Van Den Bos | A cramér-rao lower bound for complex parameters[END_REF] yields

E ĥ -h Mp 2 2 ≥ Tr ∂h Mp ∂θ I(θ) -1 ∂h Mp ∂θ H CRB, with ∂h Mp ∂θ ∂h Mp ∂θ1 , . . . , ∂h Mp ∂θ N θ ∈ C NrNtN f ×N θ
, and I(θ) ∈ R N θ ×N θ being the Fisher information matrix (FIM). Noticing that the observation y follows a complex gaussian distribution,

y ∼ CN Mh µ(θ)
, Σ , the FIM is given by the Slepian-Bangs formula [START_REF] Slepian | Estimation of signal parameters in the presence of noise[END_REF], [START_REF] Bangs | Array processing with generalized beamformers[END_REF]:

I(θ) = 2Re ∂µ(θ) ∂θ H Σ -1 ∂µ(θ) ∂θ (13) = 2Re ∂h Mp ∂θ H M H Σ -1 M ∂h Mp ∂θ , (14)
where the second line is true assuming that locally, h can be decomposed as h = h Mp + r, with r not being a function of the virtual parameters θ, so that ∂h ∂θ = ∂h Mp ∂θ . Let us now state the main result of the paper regarding the variance.

Theorem 1. Provided the model is identifiable, CRB ≥ N θ 2 M H Σ -1 M 2 , ( 15
)
with equality if the condition

C opt : ∂h H Mp ∂θ M H Σ -1 M M H Σ -1 M 2 ∂h Mp ∂θ = ∂h H Mp ∂θ ∂h Mp ∂θ is fulfilled.
The theorem is proven in appendix A, in order to keep the flow of the paper.

An important feature of this result is that the optimal CRB is proportional to the number of parameters, and consequently to the number of estimated paths p, if C opt is fulfilled. However, the condition C opt may seem a bit abstract in its current general form. In order to ease interpretation, let us consider a special case given by the following corollary.

Corollary 1. For a white gaussian noise with Σ = σ 2 Id, the bound of theorem 1 becomes

CRB ≥ N θ σ 2 2 M 2 2 , ( 16
)
with equality if C opt :

∂h H Mp ∂θ M H M M 2 2 ∂h Mp ∂θ = ∂h H Mp ∂θ ∂h Mp ∂θ .
Note that this bound is very similar to the variance obtained in the linear case in section III-A, except for the division by two (due to the fact that the parameters being real was overlooked in section III-A). This means that, at least asymptotically, the nonlinear model behaves like a linear one (provided an efficient estimator is available in the nonlinear case). Interpretations. The condition C opt is quite easily understood. It means that the observation matrix M has to preserve angles (be conformal) on the column space of ∂h Mp ∂θ . One obvious although conservative way to fulfill C opt is then to take M H M = α 2 Id, so that M is conformal on the whole space C NrNtN f (including the column space of ∂h Mp ∂θ ). This is possible only if N m ≥ N r N t N f , meaning that the number of measurements has to be greater than the dimension of the channel. Another (more clever) way of fulfilling C opt is to take M = αQU H where Q ∈ C N θ ×N θ is an unitary matrix and U ∈ C NrNtN f ×N θ has its columns forming an orthogonal basis of the column space of ∂h Mp ∂θ . That way, N m = N θ , which may be much smaller than the channel dimension N r N t N f .

Let us now look at the quantity M 2 2 when using the two aforementioned strategies with a fixed power per measurement P m . First, if each coordinate of C NrNtN f is measured K times, this yields M = √ P m (Id (1) , . . . , Id (K) ) T and N m = KN r N t N f . Second, if each basis vector of the column space of ∂h Mp ∂θ is measured K times, this yields

M = √ P m (UQ H (1) , . . . , UQ H (K) ) H and N m = KN θ . In these two cases, M 2 2 = P m K so that CRB = N θ σ 2
2PmK . The bound is thus in both cases inversely proportional to the measurements power, as well as to the number of measurements taken. However the second strategy is much more computationally efficient since it requires to take only N m = KN θ measurements instead of N m = KN r N t N f . Note that P m is strongly linked to the transmit power, for example in the case of a single subcarrier and an observation matrix of the form M = X T ⊗ Id Nr where each column of X has its squared norm equal to P t , then P m = P t .

A study of practical and resource efficient pilot designs that fulfill (exactly or approximately) C opt would be of great interest, but is beyond the scope of this paper. The condition C opt is assumed fulfilled in the remaining of this paper, yielding an optimal observation matrix.

C. Fisher information matrix

The CRB computed in the previous subsection is attained by efficient estimators [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. Maximum likelihood estimators (MLEs) are asymptotically efficient, and The Fisher information matrix (FIM) determines their asymptotic properties (it is the expected Hessian of the negative log-likelihood), and can be used to design efficient estimation algorithms. Let us compute the FIM for the considered physical model [START_REF] Magoarou | Parametric channel estimation for massive MIMO[END_REF]. It exhibits a block structure

I(θ) =   I (1,1) I (1,2) ... I (1,p) I (2,1) I (2,2) . . . . . . I (p,1) I (p,p)   ,
where, according to ( 14) assuming C opt , the off-diagonal block I (m,n) ∈ R 7×7 contains the correlations between the sensitivities of the channel to parameters of the m-th and n-th paths, and the diagonal block I (n,n) ∈ R 7×7 contains the correlations between the sensitivities to parameters of the n-th path. Let us focus here on the diagonal blocks, since as will be stated in more details in the next subsection, most estimation algorithms do not handle paths jointly but sequentially, implicitly assuming distinct estimated paths have uncorrelated effects on the channel. It remains to compute the derivatives (sensitivities)

∂h Mp ∂θ with θ defined in [START_REF] Magoarou | Parametric channel estimation for massive MIMO[END_REF]. Writing

γ n = ρ n e jφn yields ∂h Mp ∂ρ n = N r N t N f e jφn e f (ξ n ) ⊗ e t ( --→ w t,n ) * ⊗ e r ( --→ w r,n ) and ∂h Mp ∂φ n = N r N t N f jγ n e f (ξ n ) ⊗ e t ( --→ w t,n ) * ⊗ e r ( --→ w r,n ). Then, denoting ---→ b r,n,1 , ---→ b r,n,2 the two basis vectors used to describe a change in --→ w r,n and ---→ b t,n,1 , ---→ b t,n,2 those used to describe a change in --→ w t,n , ∂h Mp ∂ --→ b r,n,i = N r N t N f γ n e f (ξ n )⊗e t ( --→ w t,n ) * ⊗(a r,n,i e r ( --→ w r,n ))
with a r,n,i -j 2π λ ( --→ a r,1 .

--→ b r,n,i , . . . , ---→ a r,Nr .

--→ b r,n,i ) T , and

∂h Mp ∂ --→ b t,n,i = N r N t N f γ n e f (ξ n )⊗(a t,n,i e t ( --→ w t,n )) * ⊗e r ( --→ w r,n ) with a t,n,i -j 2π λ ( --→ a r,1 . --→ b t,n,i , . . . , ---→ a r,Nt . --→ b t,n,i ) T . Finally, ∂h Mp ∂ξ n = N r N t N f γ n (f e f (ξ n )) ⊗ e t ( --→ w t,n ) * ⊗ e r ( --→ w r,n ) with f -j2π(f 1 , . . . , f N f ) T .
Note that in section II-A the antenna locations are expressed with respect to the centroid and the frequencies with respect to the central frequency. This is not an arbitrary choice. Indeed, this carefully chosen parameterization guarantees that

1 T f = 1 T a r,n,i = 1 T a t,n,i = 0,
where 1 is a vector of the appropriate size having all entries equal to one. Taking this into account and injecting these derivatives in [START_REF] Ayach | Spatially sparse precoding in millimeter wave mimo systems[END_REF] assuming C opt yields

I (n,n) = C               1 0 0 0 0 0 ρ 2 n 0 0 0 0 0 B r,n 0 0 0 0 0 B t,n 0
0 0 0 0 ρ 2 n f 2 2               (17)
with

C = N r N t N f M H Σ -1 M 2 , 0 denotes zero vectors or matrices of appropriate size, B r,n ρ 2 n a r,n,1 2 2 a H r,n,1 a r,n,2 a H r,n,2 a r,n,1 a r,n,2 2 2
and

B t,n ρ 2 n a t,n,1 2 2 a H t,n,1 a t,n,2 a H t,n,2 a t,n,1 a t,n,2 2 2
.

The important feature of ( 17) is that the FIM I (n,n) being block diagonal, the parameters of the same path are orthogonal to each other (thanks to the chosen parameterization). Parameter orthogonality has several implications [START_REF] Cox | Parameter orthogonality and approximate conditional inference[END_REF]. It can for instance be exploited to design efficient estimation algorithms, as is the topic of the next subsection. The diagonal blocks of the FIM are not analyzed here, since they do not have any impact on the CRB as long as the model is identifiable (some interpretations regarding the diagonal terms in a simplified setting are available in [START_REF] Magoarou | Parametric channel estimation for massive MIMO[END_REF]).

D. Efficient estimation algorithms

Let us focus on algorithms aimed at obtaining an estimate taking the form of [START_REF]Bias-variance tradeoff in MIMO channel estimation[END_REF]. First, one can rewrite [START_REF]Bias-variance tradeoff in MIMO channel estimation[END_REF] as ĥ = Eα, with E (e( -→ v r,1 , -→ v t,1 , ζ 1 ), . . . , e( -→ v r,p , -→ v t,p , ζ p )) and α N r N t N f (α 1 , . . . , α p ) T . Assuming a white gaussian noise, maximum likelihood channel estimation corresponds to solve the optimization problem minimize

E,α y -MEα 2 2 , ĥ ← Eα. ( 18
)
Note that given E, the optimal vector α can be obtained as the solution of a least squares problem as α opt = (E H M H ME) -1 E H M H y, so that in the end channel estimation amounts to find an optimal E, i.e. an optimal set of p vectors {e( -→ v r,1 , -→ v t,1 , ζ 1 ), . . . , e( -→ v r,p , -→ v t,p , ζ p )}. Greedy strategy. Looking for the p vectors jointly yields a very complex optimization problem. Instead, greedy strategies have been proposed which consist in building a dictionary of characteristic vectors and applying a sparse recovery algorithm such as orthogonal matching pursuit (OMP) [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]- [START_REF] Bajwa | Compressed channel sensing: A new approach to estimating sparse multipath channels[END_REF]. This amounts to estimate the paths one by one, i.e. building the matrix E column by column. Denoting

E (k) (e( -→ v r,1 , -→ v t,1 , ζ 1 ), . . . , e( --→ v r,k , -→ v t,k , ζ k ))
the state of the matrix E at the k-th iteration, the optimal vector

α (k) ← (E (k)H M H ME (k) ) -1 E (k)H M H y is computed so that a residual r (k+1) ← y-ME (k) α (k) is used at the next iteration.
Such a strategy is summarized in algorithm 1.

Algorithm 1 Greedy channel estimation

Input: Observation y, observation matrix M, number of paths to estimate p. 1: r (1) ← y 2: for i = 1, . . . , p do

3: Estimate a characteristic vector e( -→ v r,i , -→ v t,i , ζ i ) based on r (i)
4:

E (i) ← (e( -→ v r,1 , -→ v t,1 , ζ 1 ), . . . , e( -→ v r,i , -→ v t,i , ζ i )) 5:
Update the coefficients:

α (i) ← E (i)H M H ME (i) -1 E (i)H M H y 6:
Update the residual: r (i+1) ← y -ME (i) α (i) 7: end for 8: ĥ ← E (p) α (p) Output: Channel estimate ĥ

The critical and most computationally intensive part of the above algorithm is the characteristic vector estimation (line 3), which in the maximum likelihood framework amounts to solve the optimization problem maximize

x f (x, r (i) ) |x H M H r (i) | 2 Mx 2 2
, with x being a characteristic vector. Joint estimation. One straightforward and popular way to solve this problem is given in algorithm 2 which, if testing

N ζ delays, N-→ vr DoAs and N-→ vt DoDs requires N ζ N-→ vr N-→ vt computations of the cost function f (x, r (i) ).

Algorithm 2 Joint characteristic vector estimation

Input: Current residual r (i) , observation matrix M

1: Estimate physical parameters -→ v r,i , -→ v t,i , ζ i based on r (i) -→ v r,i , -→ v t,i , ζ i ← argmax -→ vr, - → vt ,ζ f e( - → v r , - → v t , ζ), r (i) Output: Characteristic vector e( -→ v r,i , -→ v t,i , ζ i )
Sequential estimation. Let us propose a new way of estimating the characteristic vector, taking profit of the analysis of the previous subsection. The idea is to exploit the parameters orthogonality, adopting a sequential strategy presented in algorithm 3. In this algorithm, the first step (line 1) corresponds to the maximum likelihood estimation of the delay considering DoD and DoA as unknown nuisance parameters (b k being the k-th vector of the standard basis of R Nt and c l being the l-th vector of the standard basis of R Nr ), the second step (line 2) corresponds to the maximum likelihood estimation of the DoD considering a known delay and the DoA as an unknown nuisance parameter and finally the third step (line 3) corresponds to the maximum likelihood estimation of the DoA considering both delay and DoD are known. This strategy, if testing N ζ delays, N-→ vr DoAs and

N-→ vt DoDs requires N ζ N r N t + N-→ vt N r + N-→
vr computations of the cost function f (x, r (i) ). Thanks to the parameters orthogonality [START_REF] Cox | Parameter orthogonality and approximate conditional inference[END_REF], this strategy is asymptotically (at high SNR or with a great number of measurements N m ) equivalent to the joint estimation, although it may be much less complex (provided N-→ vt N t and N-→ vr N r ). Note that this is true for any estimation order (the DoD or DoA can be estimated first also, without affecting the conclusions), so the order should be chosen so as to yield the lowest complexity. This sequential strategy is empirically assessed and compared to the classical joint strategy in section V.

Algorithm 3 Sequential characteristic vector estimation

Input: Current residual r (i) , observation matrix M 1: Estimate delay ζ i based on r (i)

ζ i ← argmax ζ Nt k=1 Nr l=1 f e f (ζ) ⊗ b k ⊗ c l , r (i) 2: Estimate DoD -→ v t,i based on r (i) and ζ i -→ v t,i ← argmax - → vt Nr l=1 f e f (ζ i ) ⊗ e t ( - → v t ) * ⊗ c l , r (i) 3: Estimate DoA -→ v r,i based on r (i) , ζ i and -→ v t,i -→ v r,i ← argmax -→ vr f e( - → v r , -→ v t,i , ζ i ), r (i) Output: Characteristic vector e( -→ v r,i , -→ v t,i , ζ i )

IV. BIAS ANALYSIS

The main result of the previous section (theorem 1) indicates that the variance of the studied class of estimators is at best proportional to the number of virtual paths p. This result may lead to choose p as small as possible so as to minimize the variance. However, if p is taken too small, h Mp may become an oversimplified version of the channel h leading to a high MSE due to a high model bias hh Mp 2 2 . There is a biasvariance tradeoff. How to set p appropriately? Knowing that the number of physical paths P is in general very large (up to several hundreds), is it possible to approximate it with a few virtual paths? This section studies these questions, trying to understand the mechanisms allowing to merge a large number of physical paths into much fewer virtual paths without incurring a large bias.

Computing the bias defined in [START_REF] Jaeckel | Quadriga: A 3-d multi-cell channel model with time evolution for enabling virtual field trials[END_REF] amounts to compute the projection h Mp . Unfortunately, even considering a discretized set of candidates DoAs, DoDs and delays, this problem (which then becomes a sparse approximation problem) is NPhard [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]. The projection can be approximated numerically using sparse recovery methods, as will be done in section V. However the objective of this section is to study theoretically the bias and give an interpretable upper bound. Bound on the bias. To do so, let us consider a simple situation in which L physical paths are approximated by a single virtual path, in order to understand under which circumstances they can be merged. In that case, the physical channel is expressed

h = N r N t N f L l=1 β l e( -→ u r,l , -→ u t,l , τ l ),
and the bias is bounded as

h -h M1 2 ≤ h -h 2 for any h = N r N t N f γe( -→ w r , -→ w t , ξ),
where -→ w r , -→ w t and ξ are the DoA, DoD and delay of the approximating virtual path. This is true in particular considering the optimal coefficient

γ opt = L l=1 β l e( -→ w r , -→ w t , ξ) H e( -→ u r,l , -→ u t,l , τ l )
for which h is the orthogonal projection of the channel h onto the characteristic vector e( -→ w r , -→ w t , ξ). With this optimal coefficient, introducing e e( -→ w r , -→ w t , ξ) and e l e( -→ u r,l , -→ u t,l , τ l ) in order to lighten notations, it follows

h -h 2 = N r N t N f L l=1 β l e l -e H e l e 2 ≤ N r N t N f L l=1 β l e l -e H e l e) 2 = N r N t N f L l=1 |β l | 1 -|e H e l | 2 . ( 19
) This inequality highlights the fact that collinear characteristic vectors lead to small error. Said otherwise, physical paths with characteristic vectors collinear to the one of the approximating virtual path can be merged without accuracy loss.

Moreover, the inner product between characteristic vectors is expressed as the product of inner products between steering and delay vectors,

e H e l = [e f (ξ) H e f (τ l )][e t ( -→ u t,l ) H e t ( -→ w t )][e r ( -→ w r ) H e t ( -→ u r,l )].
Analyzing jointly the three inner products yields the following approximation bound which is the main result of this section.

Theorem 2. Provided -→ u r,l --→ w r 2 < 1 √ 2π Rr λ , -→ u t,l --→ w t 2 < 1 √ 2π R t λ and |τ l -ξ| < 1 π √ 2 B the bias is bounded by h -h M1 2 ≤ N r N t N f L l=1 |β l | 1 -(1 -x l ) 2 (1 -y l ) 2 (1 -z l ) 2 ,
where

• x l = 2π 2 (τ l -ξ) 2 1 N f N f i=1 f 2 i , • y l = 2π 2 -→ u t,l --→ w t 2 2 1 Nt Nt j=1 --→ at,j 2 2 λ 2 cos 2 ( -→ a t,j , -→ u t,l --→ w t ), • z l = 2π 2 -→ u r,l --→ w r 2 2 1 Nr Nr k=1 --→ a r,k 2 2 λ 2 cos 2 ( --→ a r,k , --→ u r,k --→ w r ).
The theorem is proven in appendix B, in order to keep the flow of the paper. Interpretations. This bound yields quite intuitive results that can be nicely interpreted:

• First of all, physical paths that are close to the approximating virtual path in the delay, DoD and DoA domains can be merged. Indeed, if the quantities

(τ l -ξ) 2 , -→ u t,l --→ w t 2 2 and -→ u r,l --→ w r 2
2 are small the l-th term of the bound is small and the l-th path can be merged with low approximation error. It is important to notice that the physical path and the virtual path have to be close in the three domains, because one domain in which they are far apart is sufficient to get almost orthogonal characteristic vectors.

• Second, the discrimination power of the system is taken into account by the bound. A larger bandwidth or larger arrays lead to more difficulty to merge physical paths. Indeed, if the quantities

N f i=1 f 2 i , Nt j=1 --→ at,j 2 2 λ 2
and Nr k=1

--→ a r,k 2 2 λ 2
increase the approximation error also increases. On the other hand, if there is for example only one antenna at the receiver, a physical path that has a close delay and DoD with the approximating virtual path can be merged irrespective of its DoA, because the system has then no discrimination power in the DoA domain. • Finally, the ability to merge paths depends also on the orientation of the transmit and receive antenna arrays. Some paths can be merged even if they are far from the virtual path in the DoD and DoA domains, if their direction differences lie in a direction in which the array is insensitive. This can be seen by the presence of the quantities cos 2 ( -→ a t,j , -→ u t,l --→ w t ) and cos 2 ( --→ a r,k , --→ u r,k --→ w r ) in the bound. for example for a ULA aligned with the z-axis, a path that has the same elevation as the approximating virtual path will be perfectly mergeable irrespective of its azimuth, because the ULA has sensitivity only in one direction (the cosine will be null in this case). The proposed bound lends itself to nice interpretations as shown above, but may be loose because of the use of the triangle inequality. A way of possibly improving the bound is evoked in appendix C.

V. EXPERIMENTS

The objective of this section is to assess empirically the mathematical developments of the previous sections. The model bias is first computed with a varying number of virtual paths. Then, the joint and sequential estimation strategies (algorithms 2 and 3) are compared, showing a bias-variance tradeoff. All experiments are done using realistic channels generated with help of the NYUSIM channel simulator [START_REF] Sun | A novel millimeterwave channel simulator and applications for 5g wireless communications[END_REF], in a millimeter wave massive MIMO context. In particular, the central frequency is set to f c = 28 GHz and the distance between transmitter and receiver to d = 30 m to obtain the DoDs, DoAs, delays, gains and phases of each path. The channel matrix is then obtained from (2) (with the total number of physical paths P being random and ranging from fifty to a hundred). All results shown in this section are averages over one hundred channel realizations.

A. Empirical evaluation of the model bias

Computing the model bias h-h Mp 2 2 requires to compute the projection onto the model h Mp = proj Mp (h). This is unfortunately a NP-hard problem even with discretized directions and delays [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]. It is nevertheless possible to compute an approximation ĥMp by using algorithm 1 called directly on h (in perfect observation conditions and without noise). The relative bias is then approximated by the relative error

h -ĥMp 2 2 h 2 2 .
For this first set of experiments, let us compute this relative error in various configurations (varying p, N t , N r , N f and the orientation of the antenna array). To do so, algorithm 1 is used in conjunction with the joint characteristic vector estimation (algorithm 2), for which the optimization problem is solved by exhaustive testing of N ζ delays, N-→ vr DoAs and N-→ vt DoDs (evenly sampling the three domains). In order to simplify the analysis, let us have a single oversampling parameter S controlling the number of tests. This amounts to take

N ζ = SN f if N f > 1, N ζ = 1
otherwise (the delay has no influence with N f = 1), and setting similarly N-→ vr and N-→ vt . The value S = 6 is taken here (it was found empirically that testing more values does not improve the result, the impact or reducing S is assessed in the next subsection). Results on average over 100 NYUSIM channel realizations are shown on figure 2. Several comments are in order:

• On all three subfigures, the error decreases quickly when the number of virtual paths p increases. This shows that even a large number of physical paths can indeed be merged into much fewer virtual paths with little accuracy loss. This justifies the use of channel estimators based on physical models taking the form of (10) for massive MIMO in the millimeter wave band. • On figure 2a, the influence of the number of antennas on the bias is assessed. To do so, let us consider an uniform linear array (ULA) aligned with the x-axis with half-wavelength separated antennas with a varying number of transmit antennas N t ∈ {16, 64, 256} at the transmitter, a single antenna receiver (N r = 1) and a single subcarrier (N f = 1). First, notice that taking more virtual paths (a greater p) leads to lower error, as expected since each physical path is then closer in average to a virtual path. This could be predicted from the bound of theorem 2, in which adding virtual paths amounts to reduce the quantity -→ u t,l --→ w t 2 2 . This conclusion holds also for figures 2b and 2c. Second, more antennas require more virtual paths to attain the same error, because the discrimination power of the antenna array is higher, as was predicted by the quantity Nt j=1

--→ at,j 2 2 λ 2
in the bound of theorem 2.

• Let us now consider a fixed number of transmit antennas (64) and add several receive antenna and several subcarriers. Results of this experiment are shown on figure 2b. It is clear that considering 12 subcarriers (with a 15 MHz spacing between adjacent subcarriers) leads to a higher discrimination power, which results in higher error for a fixed number of virtual paths. This is also the case when a ULA with 4 receive antennas is considered, but the difference is subtle, since the DoA discrimination power of an array with 4 antennas is quite low. These observations are nicely interpreted with help of the bound of theorem 2, because the quantities

N f i=1 f 2 i and Nr k=1 --→ a r,k 2 2 λ 2
are not null with several subcarriers and receive antennas, whereas they were for the previous experiment.

• Finally, changing the orientation of the antenna array also has an influence on the bias, as evidenced on figure 2c. Indeed, taking a ULA aligned with the z-axis (a vertical ULA) allows to get a lower error for a fixed number of virtual paths compared to the horizontal ULA aligned with the x-axis. This is because the physical paths generated by NYUSIM have mostly similar elevations, they differ mainly for their azimuth. This observation corresponds in the bound of theorem 2 to cos 2 ( -→ a t,j , -→ u t,l --→ w t ) that is small in most cases for the vertical ULA, since it has discrimination power only on the elevation, and not on the azimuth (this is the opposite for the horizontal ULA).

B. Comparison of joint and sequential strategies

The objective of this subsection is to compare the joint and sequential estimation strategies (algorithms 2 and 3) for several signal to noise ratios (SNR). To do so, let us consider a system with N t = 64 transmit antennas (arranged as an ULA with half-wavelength separated antennas), a single receive antenna (N r = 1) and N f = 12 subcarriers (with a 15 MHz spacing between adjacent subcarriers). When the sequential characteristic vector estimation (algorithm 3) is used, the DoD is estimated first and then the delay (the DoA being irrelevant since N r = 1). The observation matrix M is taken as the identity (the objective here being to study the channel model and not pilot sequences). The greedy channel estimation algorithm is used on noisy observations taking the form of (4) (with various noise levels) to get a channel estimate ĥ. The oversampling factor S defined in the previous subsection is taken in {2, 4, 6}. The performance is measured by the relative error hĥ 2 2 h 2

2

.

Results are shown on figure 3, where each subfigure corresponds to a different oversampling factor S. Note that the approximate bias (computed with S = 6) is shown on the three plots to serve as a reference. Several comments are in order:

• First, notice that for both characteristic vector estimation methods, the error is close to the bias for a small p, and then linearly increasing for a large p (error dominated by the variance). This is a bias-variance tradeoff. It is interesting to notice that the optimal number of virtual paths increases with the SNR (for both methods) and stays relatively small (no more than a dozen, which yields around 50 real parameters to estimate), whereas the channel vector has 768 complex entries (which yields 1536 real parameters to estimate if not using a parametric model). For a large p, the slope at which the error increases is inversely proportional to the SNR, which was predicted by corollary 1. • Second, the joint strategy (algorithm 2) performs better than the sequential strategy (algorithm 3) for all SNRs (at the optimal p). However, the difference becomes very small at high SNR. For example, the difference is approximately 3% for a SNR of -8 dB and only approximately 0.5% for a SNR of 0 dB. This is in total agreement with the parameter orthogonality between DoD and delay evidenced in section III, which implies asymptotic independence. • Finally, one can notice that increasing the oversampling factor S improves accuracy, but also increases complexity.

Indeed, the joint strategy complexity is S 2 N t N f whereas the sequential strategy complexity is SN t N f + SN f (it is roughly S times less complex). This implies that the complexity advantage of the sequential strategy increases linearly with S (it would increase quadratically with several In summary, the sequential strategy was empirically found to be almost as good as the joint strategy, whereas it is much less complex. Moreover, its complexity advantage is higher if a high accuracy is sought for. This shows that such sequential characteristic vector estimation strategies are promising for greedy channel estimation methods.

VI. CONCLUSION

In this paper, the performance of MIMO wideband channel estimators using a physical model was theoretically studied. To do so, the mean squared error (MSE) of channel estimation was decomposed in to a bias and a variance term. The variance was studied with help of the Cramér-Rao bound, shown to be proportional to the number of considered virtual paths, provided a condition on the observation matrix C opt is fulfilled. Moreover, computing the Fisher information matrix (FIM) allowed to discover that DoD, DoA and delay are orthogonal parameters, leading to the design of a computationally efficient sequential estimation algorithm (algorithm 3) with asymptotic optimality properties. Then, the bias term was bounded by a quantity lending itself to nice interpretations, depending both on the propagation properties of the channel, on the antenna arrays geometries and on the position of the subcarriers.

The aforementioned mathematical developments were then assessed experimentally. First, the bias was approximated in various configurations, verifying the predictions of the computed bound, and furthermore showing that few virtual paths are sufficient to approximate well realistic channels in the millimeter wave band. The sequential estimation algorithm was then compared to the classical joint one. It was found that from moderately high SNR it performed as well as the joint one, whereas being much more computationally efficient. This shows it is a promising alternative and should be considered for system design.

In the future, it would be of great interest to study practical designs of observation matrices that take into account practical constraints (for example hybrid systems) while fulfilling exactly or approximately the optimal observation condition C opt and requiring a reasonable amount of resources (with a number of linear measurements N m much smaller than N t N r N f ). On a more technical side, improving the bound on the bias proposed here may be possible, as mentioned in appendix C.

The normalized matrix Ā Ā 2 is symmetric positive semidefinite and has all its eigenvalues upper bounded by one so that Dividing both sides by 2 Ā 2 and noticing that the matrices A and Ā have the same eigenvalues (even though with a doubled multiplicity for Ā) so that A 2 = Ā 2 , (15) is obtained.

B. Proof of theorem 2

The inner product of characteristic vectors involved in ( 19) can be expressed as e H e l = [e f (ξ) H e f (τ l )][e t ( -→ u t,l ) H e t ( -→ w t )][e r ( -→ w r ) H e t ( -→ u r,l )],

which is a product of three inner products, that can all be expressed in the general form Noticing that C k n ≥ 0 for k even, one gets that the real part of the series is an alternating series. Its terms are decreasing in magnitude if

C k n > k! (k + 2)! C (k+2) n , ∀k, n.
This is the case for any k if it is fulfilled for k = 0, leading to the sufficient condition C 2 n < 2, ∀n. Specializing this condition to the three values C n can take allows to get the three sufficient conditions of the theorem:

• |τ l -ξ| < 1 π √ 2 B , • -→ u r,l --→ w r 2 < 1 √ 2π Rr λ , • -→ u t,l --→ w t 2 < 1 √ 2π R t λ .
If they are fulfilled, the real part of the series can be bounded,

Re 1 N N n=1 e -jCn ≥ 1 - 1 2N N n=1 C 2 n ,
which implies

1 N N n=1 e -jCn ≥ 1 - 1 2N N n=1 C 2 n .
Replacing C n and N by their respective values for the three inner products gives the theorem.

C. Improving the bound on the bias

The bound of this paper is obtained applying the triangle inequality, which may lead to a loose bound. Instead, one could write q ij = e H i e j -e H i ee H e j . Studying the properties of the matrix Q may lead to a tighter bound, but is not guaranteed to lead to as easily interpretable results.

h -h 2 = N r N t N f L l=1 β l e l -

Fig. 1 :

 1 Fig. 1: Example of system representation for the l-th path. The length of the path from the first transmit antenna to the first receive antenna is equal to the length of the red line plus the length of the blue line minus the length of the green line, which illustrates (1).

Nr = 1 ,Fig. 2 :

 12 Fig. 2: Bias assessment in several configurations with a varying number of virtual paths p.

Fig. 3 :

 3 Fig. 3: Comparison of the joint and sequential strategies for several oversampling factors S and several SNRs.

APPENDIX A. Proof of theorem 1 2

 12 Starting from[START_REF] Ayach | Spatially sparse precoding in millimeter wave mimo systems[END_REF] and denoting A M H Σ -1 M, the CRB can be re-expressed using only real matrices by introducing Tr∂h Mp ∂θ ( DT Ā D) -1 ∂h Mp ∂θ H .Moreover, since for any symmetric real matrix F and any complex matrix E,Tr EFE H = Tr Re{E} Im{E} F Re{E} T Im{E} T , Tr D( DT Ā Ā 2 D) -1 DT .

Using [ 28 ,- 1

 281 Theorem 4.3] and assuming the model is identifiable (so that the inversion is possible), it follows DT ≥ Tr D( DT D) -1 DT = N θ .

•

  taking C n = 2πf n (τ l -ξ) and N = N f yields the delay vectors inner product e f (ξ) H e f (τ l ),• taking C n = 2π λ --→ a t,n .( -→ w t --→ u t,l) and N = N t yields the transmit steering vectors inner product e t ( -→ u t,l ) H e t ( -→ w t ),• taking C n = 2π λ --→ a r,n .( -→u r,l --→ w r ) and N = N r yields the receive steering vectors inner product e r ( -→ w r ) H e r ( -→ u r,l ).This general expression allows to study these three inner products jointly.Using the series representation of the exponential e x =The antennas at both the transmitter and receiver being located with respect to the centroid of the array and the frequencies being expressed with respect to the central frequency, one gets N n=1 C n = 0, which implies that the term k = 1 of the series is null in all three cases, leading to

e H e l e 2 =

 2 β H Qβ where β = (β 1 , . . . , β L ) T and Q ∈ C L×L with


G Egoroff

P De Wever

A Cornée

6 -Just as an old tree keeps all the records of its growth and life, the Earth retains memories of its past...

The national inventory of geological heritage in France is based on the most exhaustive selection of geological sites. The following table shows the geological interests currently used in the program and the future interests that will be added to the database.

Géosite des carrières d'Orival

Intérêts

Intérêt touristique ou économique :

Site à proximité de circuits de randonnées et de la Route des Moulins.

Vulnérabilité, menaces

Vulnérabilité naturelle

Géosite des carrières d'Orival

Menaces anthropiques actuelles

Aucune.

Menaces anthropiques prévisibles

Comblement des carrières en activité.

Traçabilité Bibliographie

  https://inpn.mnhn.fr/accueil/recherche-de-donnees/inpg Concept and methods Perspective Progress Publications This diagram represents the program database and allows users to exchange program data. 4 200 geosites included in the national database 2 650 geosites validated at national level 1 088 geosites available on the national website 500 people collaborating on the project 70 institutions involved A national methodology published in French and in English An exchange standard for national use A

  Caen, prendre la D22 en direction de Creully. Après le lieu-dit Pierrepont, continuer tout droit au croisement avec la D35 et prendre le chemin à droite qui longe la Seulles. Les carrières se situent sur la droite. de nombreuses coupes remarquables des formations du Bathonien moyen, principalement le Calcaire de Creully (équivalent du Calcaire de Saint-Pierre-du-Mont), intensément exploité, surmonté par la Caillasse de Fontaine-Henry (équivalent de la Caillasse de Longues). De nature bioclastique, le Calcaire de Creully présente des litages obliques et renferme des niveaux de silex. Les structures sédimentaires de type mégarides, vagues de sable et chenaux en auge ou en Description géologique Carte(s) topographique(s) IGN à 1/25 000 BAYEUX.PORT-EN-BESSIN (1512OT) Carte(s) géologique(s) BRGM à 1Ensemble de carrières avec fronts de taille bien exposés d'une quinzaine de mètres de hauteur environ. SuperficieCet ensemble comporte une carrière exploitée et de nombreuses carrières abandonnées. L'ensemble des parties abandonnées a été classé en réserve naturelle régionale en 2008. Le calcaire exploité est une pierre de taille de renommée régionale appelée "Pierre d'Orival" ou calcaire d'Orival. L'extraction, limitée à 90000 T/an selon l'arrêté préfectoral d'exploitation du 05/02/1999 s'appliquant jusqu'au 05/02/2029, est en réalité de 20000 T/an. exceptionnellement bien exposés et témoignent des paléocourants et de différents types d'érosion. La Caillasse de Fontaine-Henry, calcaire argileux intensément bioturbé et très fossilifère, est bien exposé en sommet des carrières d'Orival. Son faciès sédimentaire particulier contraste fortement avec les faciès à litages obliques du Calcaire de Creully.Carrière en exploitation sur une partie du siteDate usage UsageLa partie classée en réserve naturelle régionale est gérée par le Conservatoire Fédératif des Espaces Naturels de Basse-Normandie (CFEN). La partie en exploitation est gérée par la société Girard et Fossez et Cie.

  taille montrent de nombreuses figures de sédimentation bien visibles : litages obliques, mégarides, chenaux, qui permettent d'étudier les paléocourants du Bathonien moyen en domaine tidal et la construction des dunes sous-marines.Intérêt(s) géologique(s) secondaire(s) : PaléontologieLa Caillasse de Fontaine-Henry renferme de nombreux fossiles stratigraphiques et de faciès. Ressources naturell Exploitation de la Pierre d'Orival (Calcaire de Creully) au moins depuis le XIVe siècle. Présence de fours à chaux de différentes générations.StratigraphieLe site abrite une des coupes de référence du Calcaire de Creully qui est l'équivalent latéral du Calcaire de Saint-Pierre-du-Mont affleurant dans les falaises du Bessin. Reconstitution de la dynamique sédimentaire de la plate-forme bathonienne. Observations de mégarides en 3 dimensions.

Example of a typical data sheet in the national inventory which allows a description and an evaluation of a geosite located in Normandy (north west of France). This geosite includes a set of quarries (old or still active), exploiting the heritage stone of Pierre de Caen. The main geological interest of this site is sedimentology.

In the database, one geological site means:

-100 fields of information to be completed (30 mandatory) -a GIS cartography -a review of existing exsitu collections of objects -a bibliographic review -an evaluation of the patrimonial interest and the need of protection

The totality of this data are in free open access

The data fulfiling the conditions of national validation are published on the INPN website (National inventory of natural heritage) which is the official tool for disseminating data concerning "Nature".

The INPN website is managed by the French National Museum of Natural History (MNHN).

For more than 10 years, dozens of books or articles on the national inventory of geological heritage have been published. The publications presented here come from a collection specifically dedicated to the geological heritage, with the support of the French Geological Society.

-English translation for the exchange standard -Create connexion with :

-collection database (https://www.recolnat.org/fr/) -natural and biodiversity database (Natura2000 sites, ... The National Geoheritage Inventory in France, from knowledge to outreach


Rim Bechaieb
email: rim.bechaieb@gmail.com

Zohra B Lakhdar

Hélène Gérard

DFT and TD-DFT Studies of Mg-Substitution in Chlorophyll by Cr (II), Fe (II) and Ni (II)

Keywords: Chlorophyll, Metal, Substitution, DFT, TD-DFT, NBO analysis

Electronic contributions in electronic transitions and energetic data associated to the Mg-substitution in chlorophyll by three transition metals: chrom (Cr 2+ ), iron (Fe 2+ ) and nickel (Ni 2+ ) have been studied theoretically using Density Functional Theory (DFT) and Time Dependent Density Functional Theory (TD-DFT) methods. The binding energies are stronger than for Mg 2+ in the case of all three cations especially in the case of Ni 2+ . The Mg-substitution process is found to be exergonic for all title elements in gas phase and in acetonitrile using both implicit and explicit models of solvation. The natural population analysis (NPA) results, which estimated by natural bond orbital (NBO) analysis, showed significant charge transfer from pheophytin ligand to the central cation. The UV-visible proprieties of the different substitution compounds have been studied using the TD-DFT method evidencing that substitution of Mg by Cr, Fe or Ni is associated to a blue shift of the Q-band for the three cations.

Introduction

The essential steps in photosynthesis, which is one of the most important energy-conversion mechanisms carried out for living systems, namely light harvesting (photon absorption) and charge separation, are mediated by green pigments named chlorophylls buried in the three dimensional structures of large specialized proteins fixed in cell or thylakoid membranes. Until recently, Chlorophyll (a) (Chla) was thought to be the principal photosynthetic pigment in reaction center of the photosynthetic units of all of cyanobacteria, algae and plants [START_REF] Blankenship | Molecular Mechanisms of Photosynthesis[END_REF][START_REF] Ke | Photosynthesis: Photobiochemistry and Photobiophysics[END_REF]. Considerable progress has been made in the past four decades in the understanding of the in vitro properties of chlorophyll and this has contributed to a better understanding of the role of Chla at the molecular level in photosynthesis [START_REF] Whitmarsh | The photosynthetic process[END_REF][START_REF] Callahan Alan | [END_REF][5][6]. The salient feature of Chla relevant for its molecular organization processes is the metal bonded in the center of the tetrapyrrole macrocycle cal led pheophytin (Pheo). It provides the binding site to nucleophilic electron donors and it functions as a switch between various types of activities. This central pocket is occupied by the magnesium cation (Mg 2+ ) in Chla (which is thus a Pheo(Mg) complex). The association of the metal with the N-atoms from the four symmetric pyrrole rings is ensured by two covalent and two coordination bonds. In parallel, mastering the potential modulation of its properties to develop new bio-inspired tools for health [7] or energy conversion [8] remains a major challenge to chemists [9]. In particular, the formation of bimetallic complexes proved to be an interesting direction. Plants easily absorb many toxic metals (M), most often referred to as "heavy metals" [10]. Once absorbed, they penetrate the plant tissues and in lower concentrations iron for example is an essential micronutrient for higher plants and algae and constituents of photosynthetic apparatus photosystem I and II. But higher external metal-concentrations in plants growing on sols polluted by heavy metals may produce a lot of damaging effects [11][12][13] among which inhibition of photosynthesis [14][15][16].

The exact processes remain unsure. The majority of the experiments are based on the most important site of chlorophyll (a) which is the Mg 2+ cation bonded in the center of tetrapyrrole macrocycle called pheophytin (Pheo). Several experimental studies propose that divalent heavy metals can replace the Mg 2+ cation of Chla to form Mg-substituted Chla complexes of general formula Pheo(M) [11,[17][18][19][20][21], the substitution resulting in the irreversible loss of photosynthetic activity [22,23]. In vivo, experiments showed a spontaneous insertion of transition metals, such as Cu and Ni [11,24,25].

From a theoretical point of view, we have recently examined the Mg-substitution by Zn and Cu, and its competition with various chelation processes at the periphery [26]. Substitution by Cu and Zn was found to be possible and the most favored process in the case of copper when implicit solvation is included. In a second study, we examined the structural, electronic and energetic data associated to the Mg-substitution in chlorophyll (a) by three major toxic pollutants: Cd 2+ , Hg 2+ and Pb 2+ . Our various approaches of free cations solvation allowed us to evidence that Mg-substitution should be possible for all title elements [27]. But the metal exchange accounting for Pheo(M) formation are still not well established for Cr, Fe and Ni. Therefore further studies are needed in order to obtain a comprehensive picture for the Mg-substitution by Cr 2+ , Fe 2+ and Ni 2+ .

We thus decided to examine theoretically the structural, electronic and energetic aspects of the substitution of the Mg center of chlorophyll (a) by Cr, Fe and Ni (Scheme 1). Quantum chemical calculations at the density functional theory (DFT) have been carried out to describe the geometries of the Pheo(M) complexes. The electronic proprieties are then probed through NPA analysis as well as modelling of the UV-visible spectra using the TD-DFT implementation at the same computational level as geometry optimization. The results are analysed and compared with those obtained from experimental studies. Scheme 1: Mg-Substitution process in Chlorophyll (a).

Computational details

Computations have been carried using the Gaussian 09 [START_REF] Gaussian | Gaussian[END_REF] program within the framework of the density functional theory at the B3LYP/6-31G* level [START_REF] Becke | [END_REF][30][31], Closed-shell system, namely Mg 2+ (3s 0 3p 0 ), but also Ni 2+ (4s 0 3d 8 4p 0 ), Cr 2+ (4s 0 3d 4 4p 0 ) and Fe 2+ (4s 0 3d 6 4p 0 ) in their singlet states were carried out using spin restricted formalism. On the other hand, open shell systems, namely Ni 2+ in its triplet state, Cr 2+ and Fe 2+ in theirs triplet and quintuplet states have been studied using a spin-unrestricted formalism. The lowest energy spin states calculated are the one retained for the following discussion. They were found to be the triple t state for Ni 2+ ion (both naked and coordinated to six molecules of acetonitrile) and the singlet for the complex formed with chlorophyll (a) (Pheo(Ni)). For Fe, the more stable spin states were calculated to be the quintuplet state for Fe 2+ ion (both naked and coordinated to six molecules of acetonitrile) and the triplet for the complex formed with chlorophyll (a). The studied Cr 2+ derivatives are always more stable in its quintuplet state. These results are in line with electronic structures obtained for porphyrin complexes [32].

They are all consistent with a simple molecular orbitals model based on the 3d orbital splitting [START_REF] Crabtree | The organometallic chemistry of the transition metals, A Wileyn New-York. 34[END_REF]. The solvated structures are in a low-field (small Δ value, see Fig. 1) high-spin octahedral like arrangement. The Pheo(M) complex corresponds to a high-field (large Δ) low-spin square planar electronic arrangement (Fig. 1). Calculations were based on a Chla model containing 73 atoms reported in Scheme 1 which was proposed in our previous work [26,27,34]. This model retains the structure of Chla except for the phytyl ester side chain which was replaced by a methyl group to reduce the number of atoms for computation.

The geometries were fully optimized in the gas phase and vibrational frequencies were computed in the harmonic approximation to confirm structural assignment to energy minima and contributions to Gibbs free energies. Empirical corrections for dispersion were evaluated through single point calculations using DFT -D3(BJ) corrections [35].

Chla and its substitution complex (Pheo(M)) can accommodate additional solvent molecules on the central cation, in addition to the tetracoordination of the Pheo ligand through binding of acetonitrile which was used as a model solvent in our previous works about the interaction of different heavy metals with chlorophyll (a) [26,27]. This solvent has been chosen as it cannot form H-bonds with the C=O groups of chlorophyll (which would be artefact toward the examined properties) and has been recently used as solvent for a detailed experimental study of transmetallation of substituted chlorophylls [36].

The Gibbs free binding energy of the first acetonitrile ligand, including all thermal corrections; BSSE, dispersion corrections and ZPE computed in the harmonic approximation, varies from slightly exothermic for Mg 2+ (-2.4 kcal mol -1 ) to slightly endothermic for Ni 2+ and Cr 2+ (+5.0 and +3.3 kcal mol -1 respectively).

As such, Chla should be considered as monosolvated, whereas Pheo(Ni) and Pheo(Cr) should be conside red as non-solvated. Nevertheless, this difference in solvation does not yield significant variations in the substitution free energies (only a 2.4 kcal mol -1 difference). As a consequence, for the sake of homogeneity, no solvation will of Pheo(M), whatever the nature of M.

We followed the same calculation strategy to model the free cations, we examined the explicit coordination of the different cations M 2+ with acetonitrile (MeCN) molecules (Eq. 1) to result in [M(MeCN)n] 2+ complexes. For each metal (Mg, Cr, Fe and Ni), the number of acetonitrile molecules ''n'' to be added was chosen as the smallest ''n'' value for which the Gibbs free energy ΔG(n) for equation 1 is positive.

Eq (1): [M(CH3CN)n] + CH3CN → [M(CH3CN)n+1] ΔG(n)
The results for the isolated cations Mg 2+ , Cr 2+ , Fe 2+ and Ni 2+ are gathered in Table 1. Gibbs free energies for coordination of an additional solvent molecule to [M(MeCN) 6 ] 2+ complex, including all thermal corrections (BSSE, dispersion and ZPE), is strongly positive (> 10 kcal mol -1 ) in the case of Mg, slightly negative in the case of Cr and Ni (-1.3 kcal mol -1 and -1.7 kcal mol -1 respectively) and optimization is even impossible in the case of iron (Fe) as decoordination of the MeCN molecule occurs during the course of the optimization. As a consequence, and for the sake of homogeneity, six-solvated cations will be used to model the structure of all free cations in acetonitrile. Modelling of the UV-Vis spectra was carried out using the TD-DFT implementation [START_REF] Gaussian | Gaussian[END_REF] at the B3LYP, CAM-B3LYP and BP86 functionals in the gas phase and at B3LYP functional with implicit solvation in acetonitrile (see Sup. Info). NPA population analyses were carried out using the NBO version [37] implemented in Gaussian.

Results and Discussion

Structural properties

In the first part, we analyzed the structural properties of the optimized geometries in the gas phase for the complexes formed by substitution of the central Mg atom by Cr, Fe and Ni. Before considering the stability of Pheo(M) complexes, some general comments can be made on their structures. The main results of the geometry optimization of the Pheo(M) monomer are reported in Table 2. At first glance, it can be seen that the metal cation, native (Mg) or substituted (Cr, Fe or Ni), is located at the center and in the plane of the chlorin ring. The value obtained for Mg is in good accord with experimental X-ray results for chlorophylls and bacteriochlorophylls in protein [38][39][40].

The substitution retains the order of M-N bond distances of Chla which is reported experimentally [38,39] and theoretically [41,42] but the average < M-N > and the out-of plane distortions vary with the metal. [a] <M-NNN> is out of the plane distortions.

[b, c, d, e] M-N bond distances between the central metal and the four nitrogen atoms (Scheme 1).

[f] The average of the M-N distances.

[g] Δ(M-N) is the difference between the longer and the shorter M-N bond.

[h] Ebind= E(Pheo(M)) -(E(M2+) + E(Pheo)).

Correlation between the geometric and energetic data for binding within Pheo(M) can be gleaned from the binding energies in Table 2.

The binding energies (Ebind) were first computed as the interaction between two charged fragments, the dianionic Pheo ligand and the dicationic metal ion according to Eq (2).

Eq (2)

:

Ebind= E(Pheo(M)) -(E(M 2+ ) + E(Pheo))
The calculated binding energies of the substitution metals are all more negative than for the native one and are quite sensitive to the nature of the metal (as a 100 kcal.mol -1 variation is observed, despite the similarity of the charges of the metal cations). They vary in the order: Mg > Fe > Cr > Ni. The trend in the binding energies parallels the <M-N> bond length in that large absolute values of Ebind are associated with shorter M-N bond distances.

Electronic properties

The population analysis of the Chla model and of the three substitution compounds Pheo(Cr), Pheo(Fe) and We next tried to evaluate the UV-visible properties of these substitution compounds (Table 3). The Q-band of chlorophyll can be quite easily reproduced computationally as it is made of a single strong transition at 579.6 nm (Qy, whereas Qx participation to the band maximum remains negligible). As this is not the case for the B-band (B1 transition at 391.6 nm and B2 at 373.1 nm, see Sup. Info.), and in line with our previous studies [26,27], we focused only on the Qy transition. A blue shift is observed when going from Mg to metals of line 4. This shift decreases when going from Cr to Ni (14.6 nm, 6.6 nm and 3 nm for Cr, Fe and Ni respectively). We concluded so, that substitution of Mg by Cr, Fe or Ni is associated to a small blue shift of the Q-transition for the three cations. Let us also mention that the oscillator strength is two time smaller in the case of Cr and Fe (open shell systems).

Better understanding of the nature of this shift can be found in a frontier molecular orbitals (FMOs) analysis [43][44][45][46]. The TD-DFT calculations reveal the main molecular orbitals contributing in the Qy electronic transition, which will be discussed in the following. The Qy transition of Chla is attributed to electronic excitations from HOMO to LUMO (85%) and from HOMO-1 to LUMO+1 (14%) (Fig. 2).

This result is in a good agreement with the observed UV-Vis absorption spectra of Chla that obtained theoretically using TD-DFT method [47][48][49].

A very similar result seems to be obtained for Pheo(Ni) as the Qy transition at 576.6 nm (f=0.246) in the UV-Vis spectrum is found to be associated to two electronic excitations, from HOMO to LUMO (89%) and from HOMO-1 to LUMO+1 (8%) (Fig. 3). Nevertheless, a close examination of different molecular orbitals presented in Fig. 2 and Fig. 3, we found a very different spatial distribution of electron. Indeed, whereas in the case of Chla, all four orbitals involved are distributed on the porphyrin ring, HOMO and LUMO+1 are metal centered in the case of Pheo(Ni), so that the transition in the latest complex is mainly (89%) metal → ligand electron transfer, and more weakly (8%) ligand → metal.

The results are more complex to decipher in the case of Cr and Fe as, since they are open shell systems and don't exhibit a strongly dominating contribution to the Qy band. As a consequence, only contribution larger than 10% will be discussed. For Cr, the Qy transition is located at 564.9 nm with oscillator strength of 0.1385 and attributed mainly to (Figure 4); HOMO-3(α)→LUMO(α) (34%), HOMO-2(α)→LUMO(α) (19%), HOMO-1(α)→LUMO(α) (26%), HOMO(α)→LUMO(α) (52%) and HOMO(β)→LUMO(β) (32%).

As evidenced, all contributions end up in the LUMO (either α or β). When looking into more details in the spatial distribution of the orbital, they are essentially ligand centred, so that all transition can be seen as intra-ligand excitations (Fig. 4). Nevertheless, in contrast with Chla, a small metal contribution is observed on most orbitals. A similar observation can be made in the case of Pheo(Fe). We found only one alpha electronic excitation from HOMO to LUMO (61%) attributed to the Qy transition at 573.7 nm with an oscillator strength 0.1266 and several β transitions: HOMO-2(β)→LUMO+1(β) (30%), HOMO-1(β)→LUMO(β) (65%), HOMO(β)→LUMO+1(β) (19%). As for Cr, most transition are between ligand centered orbitals, with a minor contribution of the metal center (Figure 5). The latter nevertheless seems most important than for Cr. As a conclusion, in most of the observed transition, the structure of the porphyrin ring is responsible for electronic transitions in Pheo(M) complexes. This result is in agreement with previous theoretical study of chlorophyll [47] and metalloporphyrin [50]. Charge transfers between metal and ligands are only observed in the case of Ni, despite the apparent similarity of the UV-visible spectrum.

Energetic properties

The thermodynamic data of the Pheo(M) coordination, cation binding and exchange (Mg-substitution) process are examined in this part.

Let us first comment the results of Mg-substitution in gas phase (Eq 3) which can be directly connected to the binding energy reported above.

Eq (3):

Pheo(Mg) + M 2+ → Pheo(M) + Mg 2+ ΔG(g)
The Gibbs free energy of the substitution process (Eq (3)) ΔG(g) is reported in Table 4. The substitution process is found to be strongly exothermic (ΔH < 0) and exergonic (ΔG < 0). The extent of the stability follows the order: Ni 2+ > Cr 2+ > Fe 2+ . However, these results cannot be used to predict the thermodynamic of the reaction in solution as it requires taking into account solvent binding, not only in chlorophyll and in its substitution complexes, but also when the cation is free in solution. [i] The Gibbs free energy of the substitution process in gas phase.

[j] Binding energy of Pheo(M) complexes with one acetonitrile molecule including BSSE and dispersion corrections.

[k] The Gibbs free energy of Pheo(M) complexes with one acetonitrile molecule including all thermal corrections (BSSE, dispersion corrections and ZPE).

[l] The Gibbs free energy of the free cation coordination with six acetonitrile molecules.

[m] The Gibbs free energy of the substitution process using explicitly solvated cations.

[n] The Gibbs free energy of the substitution process using a mixed explicit/implicit solvation model Evaluation of the substitution process is obtained by computing the Gibbs free energy using explicitly solvated cations ΔGexplicit (Eq 4, Table 4):

Eq (4): [M(MeCN)6] 2+ + Pheo(Mg) → [Mg(MeCN)6] 2+ + Pheo(M) ΔGexplicit
Mg-substitution is still favored with Cr 2+ , Fe 2+ and Ni 2+ , but it appears that the explicit solvation model in the case of acetonitrile reduces significantly the exergonicity of exchange, which can be associated with the fact that acetonitrile coordinates less strongly with Mg 2+ (-323.9 kcal mol -1 ; coordination Gibbs free energies of six acetonitrile molecules) than the other metals (-395.6 kcal mol -1 for Ni 2+ , -339.3 kcal mol -1 for Fe 2+ and -338.1 kcal mol -1 for Cr 2+ ). This is similar to the experimental hydration free enthalpies [51] of Fe 2+ (-440.2 kcal mol -1 ), Cr 2+ (-442.6 kcal mol -1 ) and Ni 2+ (-473.7 kcal mol -1 ) which are stronger than Mg 2+ (-437.8 kcal mol -1 ). Therefore, desolvation of the metal cation which is necessary for the substitution, strongly disfavors the substitution. Gibbs free energies using the explicit model of solvation (ΔG explicit) were found following the order Cr > Fe ˃ Ni.

Finally, we decided to probe another correction to the explicit model commonly proposed in the literature.

Continuum solvation effects were incorporated, leading to ΔGexplicit/implicit computed according to a mixed explicit/implicit solvation model [52].

Correction using single point PCM calculations does not modify the exergonicity of exchange, but it results an inversion in the order of the thermodynamic stabilities (ΔGexplicit/implicit): Cr > Ni ˃ Fe. This calculations lead to an important correction in all cases by increasing the Gibbs free energi es; by 2.9 kcal mol -1 , 6.1 kcal mol -1 and 9.5 kcal mol -1 for Cr, Ni and Fe respectively which indicates that implicit solvation opposes the substitution with these three metals. We examined the origin of this effect. Surprisingly, the principal corrections of Gibbs free energies using the implicit solvation were associated with the effects in the [M(MeCN) n ] 2+ complexes (correction about -135.8 kcal mol -1 in the case of Mg 2+ , -120.0 kcal mol -1 for Ni 2+ , -118.9 kcal mol -1 for Fe 2+ et -118.8 kcal mol -1 for Cr 2+ ).

Conclusions

We have studied the substitution of the central cation Mg 2+ of chlorophyll (a) by three dicationic transition metals, Cr 2+ , Fe 2+ and Ni 2+ . The substitution is found to be energetically favored and leads to compounds with very similar electronic and geometric proprieties of those of chlorophyll (a). The substituted complex does not permit the coordination of the apical ligand to the central metal and the different models of solvation, explicit or explicit/implicit, opposes the substitution in the case of acetonitrile and reduces significantly the exergonicity of exchange but it still favored in all title elements.

Fig. 1 .

 1 Fig. 1. Model 3d orbital splitting: a) octahedral-like ; b) square-planar like. Adapted from ref. 32

  Pheo(Ni) reveals that the charge on the central atom is smaller in the Pheo(Fe) and Pheo(Ni) compounds, compared to chlorophyll, but it is bigger in the case of Cr (1.437 vs 1.413 for chlorophyll (a)). Even though the central atom formally bears in all cases a +2 charge, all the cations have a computed charge significantly smaller to +2, as expected from the strong charge transfer from the Pheo macrocycle to the central cation. The charge transfer nevertheless remains significantly smaller than 1 e. A population analysis of the valence shell allows decomposing this charge transfer. In the singlet state of Pheo(Ni) the charge transfer (0.783 e) is distributed essentially on orbital 4s (0.43) with a smaller overpopulation of the d shell (8.33 instead of 8). The population analysis of the valence shell of Pheo(Cr) showed the charge transfer is equal to 0.563 e and decomposes on the 4s orbitals (0.33 e) and in a smaller extent on an overpopulation of the d shell (4.21 instead of 4 e). A similar result is obtained for Fe: 0.46e on the 4s orbitals, and overpopulation by 0.27 e on the 3d orbitals.

Fig. 2

 2 Fig. 2 Frontier molecular orbitals and associated electronic transitions of Qy band (579.6 nm) for (Chla) calculated at B3LYP/6-31G* level.

Fig. 3

 3 Fig. 3 Frontier molecular orbitals and associated electronic transitions of Qy band (576.6 nm) for (Pheo(Ni)) calculated at B3LYP/6-31G* level.

Fig. 4

 4 Fig. 4 Molecular orbitals and associated electronic transitions of Qy band (564.9 nm) for (Pheo(Cr) calculated at B3LYP/6-31G* level.

Fig. 5

 5 Fig. 5 Alpha and Beta Frontier molecular orbitals and associated electronic transitions of Qy band (573.7 nm) for (Pheo(Fe) calculated at B3LYP/6-31G* level.

Table 1

 1

		Mg 2+	Cr 2+	Fe 2+	Ni 2+
	n = 3	-39.9 (-49.1)	-48.3 (-53.1)	-37.4 (-48.8)	-44.8 (-55.9)
	n = 4	-19.2 (-27.2)	-17.3 (-28.9)	-15.9 (-23.9)	-19.8 (-31.2)
	n = 5	-11.4 (-26.4)	-7.6 (-16.7)	-6.4 (-23.6)	-15.6 (-24.3)
	n = 6	10.6 (-0.7)	-1.3 (-12.1)	Not obtained	-1.7 (-12.4)

Gibbs Free energies (ΔG(n) in kcal mol -1 ) ΔH(n) in kcal mol -1 is given in parenthesis) for the coordination of an additional acetonitrile molecule to the [M(MeCN)n] 2+ complex

Table 2

 2 Geometric (distances in Å) and Ebind (in kcal mol -1 ) parameters for Pheo(M) complexes.

		Mg 2+	Cr 2+	Fe 2+	Ni 2+
	Spin multiplicity	1	5	3	1
	<M-NNN> [a]	0.004	0.004	0.005	0.009
	M-N(I) [b]	2.031	2.027	2.015	1.979
	M-N(II) [c]	2.149	2.109	2.152	2.121
	M-N(III) [d]	2.017	2.003	2.007	1.962
	M-N(IV) [e]	2.074	2.041	2.066	2.028
	<M-N> [f]	2.068	2.045	2.060	2.022
	Δ(M-N) [g]	0.132	0.106	0.145	0.159
	E bind [h]	-721.6	-762.6	-758.8	-815.6

Table 3

 3 NPA charge, population and Qy state (transition wave length in nm, oscillator strength in parenthesis) of

	Pheo(M) complexes in their stable states
		Mg 2+	Cr 2+	Fe 2+	Ni 2+
	NPA charge	1.413	1.437	1.343	1.218
	Valence pop NPA	3s(0.26) 3p(0.32)	4s (0.33) 3d (4.21) 4p (0.02)	4s (0.37) 3d (6.27) 4p (0.02) 5s (0.01)	4s (0.43) 3d (8.33) 5s (0.01) 5p (0.02)
	Qy state	579.6	564.9	573.6	576.6
	(B3LYP (gas))	(0.24)	(0.14)	(0.13)	(0.25)

Table 4

 4 Energetic data (in kcal mol-1) for Pheo(M) complexes

		Mg 2+	Cr 2+	Fe 2+	Ni 2+
	G(g) [i]	-	-40.3	-37.0	-92.4
	EMeCN [j]	-20.4	-8.0	-3.8	-6.2
	GMeCN [k]	-7.9	3.6	7.8	3.9
	Gs (M 2+ ) [l]	-323.9 -338.1	-339.3 -395.6
	Gexplicit [m]	-	-26.2	-21.8	-20.8
	Gexplicit/implicit [n] -	-23.3	-12.3	-14.7