V Laurain

R Tóth
email: r.toth@tudelft.nl.

M Gilson

H Garnier

Direct Identification of Continuous-Time LPV Input/Output Models

Keywords: Continuous-time models, LPV models, system identification, refined instrumental variable, Box-Jenkins models, input/output

Controllers in the linear parameter-varying (LPV) framework are commonly designed in continuoustime (CT) requiring accurate and low-order CT models of the system. However, identification of continuous-time LPV models is largely unsolved, representing a gap between the available LPV identification methods and the needs of control synthesis. In order to bridge this gap, direct identification of CT LPV systems in an input-output setting is investigated, focusing on the case when the noise part of the data generating system is an additive discrete-time colored noise process. To provide consistent model parameter estimates in this setting, a refined instrumental variable (IV) approach is proposed and its properties are analyzed based on the prediction-error framework. The benefits of the introduced direct CT-IV approach over identification in the discrete-time case are demonstrated through a representative simulation example inspired by the Rao-Garnier benchmark.

I. INTRODUCTION

The framework of linear parameter-varying (LPV) systems was introduced in the 1990s with the purpose to handle in a simple but efficient way the often nonlinear or time-varying nature of systems encountered in practice. The LPV system class forms an intermediate step between linear time-invariant (LTI) systems and nonlinear/time-varying plants as the signal relations in LPV systems are considered to be linear just as in the LTI case, but the parameters are assumed to be functions of a measurable time-varying signal, the so-called scheduling variable p : Z → P.

Here the compact set P ⊂ R n P denotes the scheduling space. The scheduling variable p represents

• Direct approaches: The methods formulate the identification of the CT model directly based on samples of the measured CT signals.

Unfortunately, transformation of DT-LPV models to CT-LPV models is more complicated than in the LTI case and despite recent advances in LPV discretization theory (see [START_REF] Tóth | Discretization of linear fractional representations of LPV systems[END_REF], [START_REF] Tóth | On the discretization of LPV state-space representations[END_REF]) the theory of CT realization of DT models is still in an immature state. Thus, it is often difficult in practice to obtain an adequate CT realization from an identified DT model. In order to illustrate the underlaying problems, consider the following simple CT-LPV model

d dt y(t) + p(t)y(t) = bu(t), (1)
where p, y and u are the scheduling, output and input signals of the system respectively and b ∈ R is a constant parameter. When approximating the derivative in DT by, for example, using the backward Euler approximation: d dt y(t k ) ≈ y(t k )-y(t k-1 )

Ts

with T s > 0 the sampling period, [START_REF] Apkarian | A convex characterization of gain-scheduled H∞ controllers[END_REF] transforms into:

y(t k ) = 1 1 + T s p(t k ) y(t k-1 ) + bT s 1 + T s p(t k ) u(t k ). (2)
This discretized model has now two p-dependent coefficients to be estimated instead of the one single constant parameter in [START_REF] Apkarian | A convex characterization of gain-scheduled H∞ controllers[END_REF]. Moreover, the dependence of the coefficients on p is not linear anymore but rational with a singularity whenever p(t k ) = -1

Ts . An alternative way to approximate derivatives in DT is to apply a forward Euler approximation: d dt y(t k ) ≈ y(t k+1 )-y(t k )

Ts

, which gives y(t k ) = (1 -T s p(t k-1 ))y(t k-1 ) + bT s u(t k-1 ).

This discretized model has only one p-dependent coefficient and the linearity of the dependence is preserved, however now the model equation is dependent on p(t k-1 ) instead of p(t k ). This socalled dynamic dependence (dependence of the model coefficients on time-shifted versions of p)

is a common result of model transformations in the LPV case and rises problems in LPV system identification and control alike (see [START_REF]Flexible model structures for LPV identification with static scheduling dependency[END_REF]). Furthermore, it is well known in numerical analysis that the forward Euler approximation is more sensitive for the choice of T s in terms of numerical stability than the backward Euler approximation [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF]. This means that (3) requires much faster sampling rate than [START_REF] Scherer | Mixed H2/H∞ control for time-varying and linear parametrically-varying systems[END_REF] to give a stable approximation of the system and it is more sensitive to parameter uncertainties which rises problems if (3) is used for estimation. Consequently it can be concluded that even for a very simple CT-LPV model, estimation of a DT model with February 2, 2011 DRAFT the purpose of obtaining afterwards a CT realization is a tedious task with many underlaying problems for which there are no general theoretical solutions available.

Unlike an indirect approach, a direct solution offers a way to efficiently overcome the previous problems. Due to the recent technological developments of sampling instruments in terms of achievable sampling rate, use of direct CT approaches in identification has recently regained interest showing a better performance than indirect approaches for both linear and nonlinear models, see e.g. [START_REF] Rao | Identification of continuous-time systems[END_REF], [START_REF] Rao | Numerical illustrations of the relevance of direct continuous-time model identification[END_REF]- [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF]. An exhaustive review of direct estimation methods can be found in [START_REF] Rao | Identification of continuous-time systems[END_REF], [START_REF] Garnier | Continuous-time model identification from sampled data: implementation issues and performance evaluation[END_REF], [START_REF] Garnier | Identification of Continuous-time Models from Sampled Data[END_REF]. Among the available identification approaches for CT-LTI inputoutput (IO) models, the interest for instrumental variable (IV) methods has been growing in the last years [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF], [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF]. The main reason of this increasing interest is that IV methods offer similar performance as extended least square (LS) methods or other prediction-errorminimization (PEM) methods (see [START_REF]Identification of continuous-time systems: direct or indirect ?[END_REF], [START_REF] Ljung | Experiments with identification of continuous-time models[END_REF]) and provide consistent results even for an imperfect noise structure which is the case in most practical applications. These approaches have been used in many different frameworks such as direct CT [START_REF]Identification of continuous-time systems: direct or indirect ?[END_REF], [START_REF] Garnier | Identification of Continuous-time Models from Sampled Data[END_REF], direct nonlinear CT [START_REF] Laurain | Refined instrumental variable methods for identification of Hammerstein continuous-time Box-Jenkins models[END_REF] or closedloop CT identification [START_REF] Gilson | Optimal instrumental variable method for closed-loop identification[END_REF], [START_REF] Gilson | Instrumental variable methods for closed-loop system identification[END_REF] and lead to optimal estimates in the LTI case if the system belongs to the model set defined.

In this paper we aim to provide the very first step towards bridging the existing gap between LPV control and identification via the introduction of a direct CT identification approach that benefits from the properties of IV methods. It was shown recently in [START_REF] Laurain | Refined instrumental variable methods for identification of LPV outputerror and Box-Jenkins models[END_REF] that in order to minimize the classical prediction error for DT-LPV models, a Multiple Input Single Output (MISO)-LTI reformulation of the data-generating system is needed. Based on a similar reformulation, the prediction error minimization problem can be clearly stated in the present CT case. Furthermore, the proposed approach extends the recent results from the CT-LTI identification framework using an IV method to face the direct CT-LPV identification problem stated here. The resulting approach not only provides the very first global LPV identification method that is able to provide consistent estimates of LPV-IO models in continuous-time, but it is also applicable in case of colored output noise and has a low computational load. Furthermore, it opens the possibility for closed-loop CT-LPV identification.

The paper is organized as follows: in Section II, the general class of CT-LPV systems in an IO representation form is introduced pointing out the main difficulties of this model class.

Additionally, a reformulation of the dynamical description of LPV data generating plants in the February 2, 2011 DRAFT considered setting is developed which makes possible the extension of CT LTI-IV methods to the LPV framework. In Section III, CT LPV-IV methods are described and analyzed, while their performance are illustrated in Section IV through a representative simulation example inspired by the Rao-Garnier benchmark. Finally in Section V, the main conclusions of the paper are drawn and directions of future research are indicated.

II. PROBLEM DESCRIPTION A. System description

Consider the data generating CT LPV system described by the following equations

S o      A o (p t , d )χ o (t) = B o (p t , d )u(t), y(t) = χ o (t) + v o (t), (4)
where d denotes the differentiation operator w.r.

a o i (p t )d na-i and B o (p t , d ) = n b j=0 b o j (p t )d n b -j . (5)
Note that a o i and b o j are functions of p at time t, which is called static dependence. In LPV system theory, a more general p-dependence of coefficients than static is required to establish equivalence of representations. In particular, the coefficients a o i and b o j need to depend also on the time derivatives of p, which is called dynamic dependence [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]. In order to simplify the upcoming discussion, we restrict our attention to static dependence. Nevertheless, the established results hold also in the case of dynamic dependence of (4) and of the proposed model structure.

In terms of identification we can assume that sampled measurements of (y, p, u) are available with a sampling period T s > 0. Hence, we will denote the discrete-time samples of these signals as u(t k ) = u(kT s ), where k ∈ Z. The basic idea to solve the noisy CT modeling problem of (4) is to assume that the CT noise process v o (t) can be considered at the sampling instances as a DT noise process filtered by a DT transfer function. In this paper, a practically general case is considered where the colored noise associated with the sampled output measurement y(t k ) is assumed to have a rational spectral density which might has no relation to the actual process dynamics of S o . As a preliminary step towards the case of a p-dependent noise, it is also assumed that this spectral density is not dependent on p, like in case of a measurement noise.

Therefore, v o is represented by a discrete-time autoregressive-moving-average (ARMA) model

v o (t k ) = H o (q)e o (t k ) = C o (q -1 ) D o (q -1 ) e o (t k ), (6)
where e o (t k ) is a DT zero-mean white noise process, q -1 is the backward time shift operator,

i.e. q -i u(t k ) = u(t k-i ). C o and D o are monic polynomials having constant coefficients. This formulation of the noise model ( 6) avoids the rather difficult mathematical problem of treating sampled CT random process in terms of a filtered piece-wise constant CT noise source (see [29], [START_REF] Johansson | Identification of continuous-time models[END_REF]). In terms of (6), y(t k ) can be written as

y(t k ) = χ o (t k ) + v o (t k ), (7)
which corresponds to a so called hybrid Box-Jenkins system concept already used in CT identification of LTI systems (see [29], [START_REF] Johansson | Identification of continuous-time models[END_REF], [START_REF] Laurain | Refined instrumental variable methods for identification of Hammerstein continuous-time Box-Jenkins models[END_REF]). Furthermore, in terms of ( 6), exactly the same noise assumption is made as in the classical DT Box-Jenkins models (see [START_REF] Ljung | System Identification : Theory for the User -Second edition[END_REF]).

B. Model structure considered 1) Process model:

The process model is denoted by G ρ and defined in a form of an LPV-IO representation with a static scheduling dependence:

G ρ : (A(p t , d , ρ), B(p t , d , ρ)) , (8)
where the p-dependent polynomials A and B given as

A(p t , d , ρ) = d na + na i=1 a i (p t )d na-i and B(p t , d , ρ) = n b j=0 b j (p t )d n b -j ,
are parameterized as

a i (p t ) = a i,0 + nα l=1 a i,l f l (p t ), i = 1, . . . , n a , ( 9a
) b j (p t ) = b j,0 + n β l=1 b j,l g l (p t ), j = 0, . . . , n b , (9b)
In this parametrization, {f l } nα l=1 and {g l } n β l=1 are meromorphic functions of p, with static dependence, allowing the identifiability of the model (they can be chosen for example as linearly independent functions on P). The associated model parameters are stacked column wise:

February 2, 2011 DRAFT ρ = [ a 1 . . . a na b 0 . . . b n b ] ⊤ ∈ R nρ , (10)
where

a i = [ a i,0 a i,1 . . . a i,nα ] ∈ R nα+1 , b j = [ b j,0 b j,1 . . . b j,n β ] ∈ R n β +1 and n ρ = n a (n α + 1) + (n b + 1)(n β + 1). Introduce also G = {G ρ | ρ ∈ R nρ },
as the collection of all process models in the form of [START_REF] Zhou | Essentials of Robust Control[END_REF].

2) Noise model: The noise model is denoted by H and defined as a DT-LTI transfer function:

H η : (H(q, η)) , (11)
where H is a monic rational function given in the form of

H(q, η) = C(q -1 , η) D(q -1 , η) = 1 + c 1 q -1 + . . . + c nc q -nc 1 + d 1 q -1 + . . . + d n d q -n d . ( 12
)
The associated model parameters η are stacked columnwise in the parameter vector,

η = [ c 1 . . . c nc d 1 . . . d n d ] ⊤ ∈ R nη , (13)
where

n η = n c + n d . Additionally, denote H = {H η | η ∈ R nη },
the collection of all noise models in the form of [START_REF] Tóth | Discrete time LPV I/O and state space representations, differences of behavior and pitfalls of interpolation[END_REF].

3) Whole model: With respect to a given process and noise part (G ρ , H η ), the parameters can be collected as θ = [ ρ ⊤ η ⊤ ] and the signal relations of the LPV-BJ model, denoted in the sequel as M θ , are defined as:

M θ              A(p k , d , ρ)χ(t) = B(p k , d , ρ)u(t) v(t k ) = C(q -1 , η) D(q -1 , η) e(t k ) y(t k ) = χ(t k ) + v(t k ) (14)
Based on this model structure, the model set, denoted as M, with process (G ρ ) and noise (H η ) models parameterized independently, takes the form

M = (G ρ , H η ) | col(ρ, η) = θ ∈ R nρ+nη . ( 15
)
This set corresponds to the set of candidate models in which we seek the model that explains data gathered from S o the best, under a given identification criterion (cost function).

C. Predictors and prediction error

Similar to the LTI case, in the LPV prediction error framework, one is concerned about finding a model in a given LPV model structure M, which minimizes the statistical mean of the squared prediction error based on past samples of (y, u, p). However in the LPV case, no transfer function February 2, 2011 DRAFT representation of systems is available. Furthermore, multiplication with d is not commutative over the p-dependent coefficients [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF], meaning that

d (B(p, d )u(t)) = n b j=0 ∂b j ∂p (p t )d p t d n b -j u(t) + B(p, d )d u(t) (16)
which is not equal to B(p, d )d u(t).

In the DT case, in order to define predictors with respect to models M θ ∈ M, a convolution type representation of the system dynamics, i.e. an LPV Impulse Response Representation (IRR), is used where the coefficients have dynamic dependence on p [9], [START_REF] Laurain | Refined instrumental variable methods for identification of LPV outputerror and Box-Jenkins models[END_REF]. Considering the CT case, no IRR has been developed yet and thus the same concept cannot be used to define the predictors.

1) System reformulation and prediction error:

Under the assumed noise conditions and for DT-LPV-IO models, it was shown in [START_REF] Laurain | Refined instrumental variable methods for identification of LPV outputerror and Box-Jenkins models[END_REF] that an efficient way to deal with the LPV identification problem in the PEM framework is to express the LPV system as a MISO LTI model. Therefore, based on the same idea, if the system belongs to the model set defined and with a deterministic p signal, it is possible to express the CT LPV system as a CT MISO LTI system by rewriting the signal relations of (4) as

χ (na) o (t) + na i=1 a o i,0 χ (na-i) o (t) Fo(d )χo(t) + na i=1 nα l=1 a o i,l f l (p(t))χ (na-i) o (t) χ o i,l (t) = n b j=0 n β l=0 b o j,l g l (p(t))u (n b -j) (t) u j,l (t) (17)
where g 0 (t) = 1 and the superscript (n) for a signal, like u (n) , denotes the n th time-derivative of the signal, e.g.

u (n) (t) = d n u(t). Furthermore, F (d ) = d na + na i=1 a i,0 d na-i while u (n) (t k ) represents the value of the signal u (n) (t) sampled at time instance t k .
Note that in this way, the time variation of the coefficients is transposed onto the signals

χ o i,l (t) = f l (p(t))χ (na-i) o (t), {i, l} ∈ {1 . . . n a , 1 . . . n α }, (18a)
u j,l (t) = g l (p(t))u (n b -j) (t), {j, l} ∈ {1 . . . n b , 1 . . . n α }. (18b)
Therefore, the process part of the LPV-BJ model is rewritten as a Multiple-Input Single-Output

(MISO) system with (n b + 1)(n β + 1) + n a n α inputs {χ o i,l } na,nα i=1,l=1 and {u j,l } n b ,n β j=0,l=0
. By using ( 17), ( 14) can be rewritten in terms of the sampled output signal y(t k ) as

y(t k ) = - na i=1 nα l=1 a o i,l F o (d ) χ o i,l (t k ) + n b j=0 n β l=0 b o j,l F o (d ) u k,j (t k ) Go(χo,u,t k ) + v o (t k ), (19)
February 2, 2011 DRAFT which is a sampled LTI representation of the system defined in (4).

Given the assumption that v o (t k ) = H o (q)e o (t k ) and C(q -1 ) is a monic polynomial, ( 6) can be rewritten in the form

v o (t k ) = e o (t k ) + ∞ i=1 h i e o (t k-i ). ( 20
)
This shows that the knowledge of {v o (τ )} τ ≤t k-1 implies the knowledge of {e o (τ )} τ ≤t k-1 . Therefore, by using the traditional approach [START_REF] Ljung | System Identification : Theory for the User -Second edition[END_REF], the prediction of v o (t k ) is considered as the condi-

tional expectation of v o (t k ) based on {e o (τ )} τ ≤t k-1
which is according to [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF]:

v(t k ) = v(t k | t k-1 ) = E{v o (t k ) | {e o (τ )} τ ≤k-1 } = ∞ i=1 h i e o (t k-i ), ( 21
)
where E is the expectation operator. Assuming that H o has a stable inverse such that e o (t k ) =

H -1 o (q)v o (t k )
, then the classical one-step-ahead predictor can be given as [START_REF] Ljung | System Identification : Theory for the User -Second edition[END_REF]

v(t k ) = v(t k | t k-1 ) = v o (t k ) -e o (t k ) = (1 -H -1 o (q))v o (t k ). (22)
Consequently, for the considered LPV system formulated as in [START_REF] Laurain | Refined instrumental variable methods for identification of Hammerstein continuous-time Box-Jenkins models[END_REF], the one-step-ahead predictor of y(t k ) (defined as the conditional expectation

ŷ(t k | t k-1 ) on {y(t i )} i≤k-1 , {u(t i ), χ o (t i )} i≤k )
is given by

ŷ(t k ) =H -1 o (q)G o (χ o , u, t k ) + (1 -H -1 o (q))(y(t k )), ŷ(t k ) =H -1 o (q) - na i=1 nα l=1 a o i,l F o (d ) χ o i,l (t k ) + n b j=0 n β l=0 b o j,l F o (d ) u k,j (t k ) + (1 -H -1 o (q))(y(t k )). (23)
2) Prediction Error Model: Using the same idea as in Subsection II-C1, the LPV model from ( 14) can also be expressed in a MISO LTI form [START_REF] Laurain | Refined instrumental variable methods for identification of LPV outputerror and Box-Jenkins models[END_REF]:

y θ (t k ) = - na i=1 nα l=1 a i,l F (d , ρ) χ i,l (t k ) + n b j=0 n β l=0 b j,l F (d , ρ) u k,j (t k ) + H(q, η)e(t k ). (24)
Therefore, similarly to the LTI case, the one-step-ahead prediction error can be expressed and defined as [START_REF] Ljung | System Identification : Theory for the User -Second edition[END_REF]:

ε θ (t k ) = y(t k ) -ŷθ (t k ), (25)
February 2, 2011 DRAFT where ŷθ (t k ) is the one-step-ahead predictor based on the model ( 14) written as in [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF] and is defined as (see ( 23)):

ŷθ (t k ) = H -1 (q, η) - na i=1 nα l=1 a i,l F (d , ρ) χ i,l (t k ) + n b j=0 n β l=0 b j,l F (d , ρ) u k,j (t k ) + 1 -H -1 (q, η) y(t k ). ( 26
)
3) Prediction error minimization:

Denote D N = {y(t k ), u(t k ), p(t k )} N k=1 a data sequence of S o .
Then to provide an estimate of θ based on the minimization of ε θ , an identification criterion W (D N , θ) can be introduced, like the least square criterion

W (D N , θ) = 1 N N k=1 ε 2 θ (t k ), (27)
such that the parameter estimate is

θN = arg min θ∈R nρ+nη W (D N , θ). ( 28
)

4) CT filtering and sampled data:

The hybrid representation of the model ( 14) is a combined operation of CT filtering and DT filtering which implicitly appears in the formulation of [START_REF] Gilson | Optimal instrumental variable method for closed-loop identification[END_REF].

In order to clearly define the coexistence of DT and CT filtering in [START_REF] Gilson | Optimal instrumental variable method for closed-loop identification[END_REF], a detailed investigation and discussion about the assumptions and the structure of the model are needed. In this paper, we considered the practically feasible situation such that only sampled measurements of the CT signals (y, p, u) are available. In order to apply a CT filter on sampled data one can either interpolate the samples to obtain a continuous-time signal and apply the CT filter on this reconstructed signal or use a numerical approximation, i.e. DT approximation of the considered system. This is a common problem for simulation of continuous-time systems. For simulation purposes, DT approximation of the system can efficiently be dealt with by using powerful numerical algorithms available [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF]. Note that to derive an accurate DT approximation of the system, it is often sufficient in terms of the classical discretization theory to assume that the sampled free CT signals of the system are restricted to be constant in the sampling period [START_REF] Goodwin | Control System Design[END_REF],

which has also been shown in case of LPV systems with static dependence [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]. This provides the hypothesis, also used in [29], [START_REF] Johansson | Identification of continuous-time models[END_REF], that if CT (p, u) are piecewise constant between two samples, then the trajectory of y is completely determined by its observations at the sample period T s k.

Therefore, under these inter-sampling conditions, the following operation is well-defined [START_REF] Garnier | Identification of Continuous-time Models from Sampled Data[END_REF]:

February 2, 2011 DRAFT (F (d )y) (t k ) = F (d )y(t k ), (29)
Under this assumption, and considering that a CT filter can only be applied to sampled data through numerical approximation, the usual filter properties such as commutativity holds between a DT filter and the numerical approximation of a CT filter. Nevertheless, it is important to notice that the numerical approximation method used for the evaluation of a CT filter does not have any impact on the coefficients to be estimated which remain, in terms of ( 26), the coefficients of the parsimonious CT model.

D. Identification problem statement

Based on the previous considerations, the identification problem addressed in the sequel can now be defined.

Problem 1: Given a CT-LPV data generating system S o defined as in (4) and a data set D N collected from S o . Based on the hybrid LPV-BJ model structure M θ defined by ( 14), estimate the parameter vector θ using D N under the following assumptions:

A1 S o ∈ M, i.e. there exits a G o ∈ G and a H o ∈ H such that (G o , H o ) is equal to S o . A2 In (9a-b) {f l } nα l=1 and {g l } n β l=1 are chosen such that (G o , H o ) is identifiable. A3
u(t k ) and p(t k ) are not correlated to e o (t k ).

A4

D N is informative with respect to M.

A5

S o is globally BIBO stable, i.e. for any trajectory of p : R → P and any bounded input signal u, the output of S o is bounded [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF].

III. REFINED INSTRUMENTAL VARIABLE FOR LPV SYSTEMS

Based on the MISO-LTI formulation [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF], it becomes theoretically possible to achieve optimal PEM using linear regression [START_REF] Laurain | Refined instrumental variable methods for identification of LPV outputerror and Box-Jenkins models[END_REF]. This allows to extend the Refined Instrumental Variable (RIV) approach of the LTI framework to provide an efficient way of identifying hybrid LPV-BJ models.

A. Linear Regression for CT LPV-BJ models

Using the LTI model ( 14) reformulated as in [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF], y(t k ) can be written in the regression form:

y (na) (t k ) = ϕ ⊤ (t k )ρ + ṽ(t k ), (30) where,
ϕ(t k ) = [ -y (na-1) (t k ) . . . -y(t k ) -χ 1,1 (t k ) . . . -χ na,nα (t k ) u 0,0 (t k ) . . . u n b ,n β (t k ) ] ⊤ ρ = [ a 1,0 . . . a na,0 a 1,1 . . . a na,nα b 0,0 . . . b n b ,n β ] ⊤ February 2, 2011
DRAFT and ṽ(t k ) = F (d , ρ)v(t k ). The extended regressor in [START_REF] Johansson | Identification of continuous-time models[END_REF] contains the noise-free output terms {χ i,k }. Therefore, by momentary assuming that {χ i,l (t k )} na,nα i=1,l=0 are known a priori, the prediction error ε θ (t k ) for ( 30) is given in terms of (25) as:

ε θ (t k ) = D(q -1 , η) C(q -1 , η) 1 F (d , ρ) F (d , ρ)y(t k ) -- na i=1 nα l=1 a i,l χ i,l (t k ) + n b j=0 n b l=0 b j,l u j,l (t k ) . (31)
In the given context, the filters D(q -1 ,η) C(q -1 ,η) and F (d , ρ) in ( 31) commute (see Section II-C4). The later allows us to rewrite the one-step-ahead prediction error [START_REF] Ljung | System Identification : Theory for the User -Second edition[END_REF] associated with [START_REF] Johansson | Identification of continuous-time models[END_REF] as

ε θ (t k ) = (F (d , ρ)y f ) (t k ) - na i=1 nα l=1 a i,l χ f i,l (t k ) + n b j=0 n β l=0 b j,l u f k,j (t k ), (32)
where y f (t k ), u f j,l (t k ) and χ f i,l (t k ) represent the outputs of a hybrid prefiltering operation, involving the CT and DT filters (see [START_REF] Young | Refined instrumental variable identification of continous-time hybrid Box-Jenkins models[END_REF]):

Q c (d , ρ) = 1 F (d , ρ)
and

Q d (q -1 , η) = D(q -1 , η) C(q -1 , η) . (33)
In other words:

y f (t k ) = D(q -1 , η) C(q -1 , η) 1 F (d , ρ) y (t k ) . (34)
Based on [START_REF] Goodwin | Control System Design[END_REF], the associated linear-in-the-parameters model takes the form [START_REF] Young | Refined instrumental variable identification of continous-time hybrid Box-Jenkins models[END_REF]:

y (na) f (t k ) = ϕ ⊤ f (t k )ρ + ṽf (t k ), (35)
where

ϕ f (t k ) = [ -y (na-1) f (t k ) . . . -y f (t k ) -χ f 1,1 (t k ) . . . -χ f na,nα (t k ) u f 0,0 (t k ) . . . u f n b ,n β (t k ) ] ⊤ ṽf (t k ) = Q d (q -1 , η)Q c (d , ρ)ṽ(t k ) = e(t k ).

B. The refined instrumental variable approach

Under the assumption that both the inverse noise model Q d (q -1 , η) and the CT filter

Q c (d , ρ)
and consequently {χ i,l (t k )} na,nα i=1,l=0 are known a priori, traditional parametric estimation methods from the LTI framework could provide efficient estimates of ρ and η. However, in a practical situation, Q d (q -1 , η) and Q c (d , ρ) are unknown and need to be estimated as well.

Furthermore, it is important to notice here that the regressors in [START_REF] Campi | An iterative identification method for linear continuous-time systems[END_REF] and (30) contain some time-derivatives of y and u which, in the assumed framework considering sampled data, can only be approximated. It is well-known that the approximation of derivatives requires a low pass filtering of y and u. The most commonly used filters for this purpose are Poisson's filters, February 2, 2011 DRAFT or state-variable filters [START_REF] Garnier | Continuous-time model identification from sampled data: implementation issues and performance evaluation[END_REF]. The drawback of these filters is that they require the choice of a design variable. However, in the proposed approach F (d , ρ) achieves this stable low-pass filtering directly. Therefore, it is a particular strength of the presented reformulation [START_REF] Campi | An iterative identification method for linear continuous-time systems[END_REF] is that the estimated filter F (d , ρ) is not only used for the minimization of the prediction error but it also provides the filtering for the approximation of the time derivatives. In order to estimate the parameter vector in [START_REF] Campi | An iterative identification method for linear continuous-time systems[END_REF] without the prior knowledge of Q d (q -1 , η) and Q c (d , ρ), the refined instrumental variable (RIV) method is proposed due to the following reasons:

• RIV methods lead to optimal estimates in the LTI case if S o ∈ M (see [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF], [START_REF] Young | Refined instrumental variable identification of continous-time hybrid Box-Jenkins models[END_REF]).

• In a practical identification scenario, G o ∈ G might be fulfilled due to first principle or expert's knowledge, however, it is commonly fair to assume that H o / ∈ H. In such case, RIV methods have the advantage of providing consistent estimates whereas methods such the extended LS are biased and more advanced PEM methods need robust initialization [START_REF] Ljung | Experiments with identification of continuous-time models[END_REF].

• The RIV algorithm has been successfully used for models with similar hybrid structure, like in the case of linear models [START_REF] Rao | Numerical illustrations of the relevance of direct continuous-time model identification[END_REF], [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF] and nonlinear ones [START_REF] Laurain | Refined instrumental variable methods for identification of Hammerstein continuous-time Box-Jenkins models[END_REF].

Aiming at the extension of the RIV approach for the estimation of hybrid LPV-BJ models, consider the relationship between the process input and output signals as in [START_REF] Johansson | Identification of continuous-time models[END_REF]. Based on this form, the extended-IV estimate is given as [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF]:

ρXIV (N ) = arg min ρ∈R nρ 1 N N k=1 ζ f (t k )ϕ ⊤ f (t k ) ρ- 1 N N t=1 ζ f (t k )y (na) f (t k ) 2 W ,
where ζ(t k ) is the instrument, x 2 W = x T W x, with W a positive definite weighting matrix and the filtered variables ζ f , ϕ f and y f are constructed using a stable prefilter. If G o ∈ G, the extended-IV estimate is consistent under the following two conditions 2 :

C1 Ē{ζ f (t k )ϕ ⊤ f (t k )} is full column rank. C2 Ē{ζ f (t k )ṽ f (t k )} = 0.
Moreover it has been shown in [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF] and [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF] that the minimum variance estimator can be achieved if: C3 W = I. 2 The notation Ē{.} = limN→∞ 1 N P N t=1 E{.} is adopted from the prediction error framework of [START_REF] Ljung | System Identification : Theory for the User -Second edition[END_REF].

February 2, 2011 DRAFT

C4

ζ is chosen as the noise-free version of the extended regressor in [START_REF] Johansson | Identification of continuous-time models[END_REF] and is therefore defined in the present LPV case as:

ζ(t k )= -χ (na-1) (t k ) . . . -χ(t k ) -χ 1,1 (t k ) . . .-χ na,nα (t k ) u 0,0 (t k ) . . . u n b ,n β (t k ) ⊤ C5
G o ∈ G and n ρ is equal to the minimal number of parameters required to represent G o with the considered model structure.

C6

The used hybrid filter is chosen as the filter chain [START_REF] Young | Refined instrumental variable identification of continous-time hybrid Box-Jenkins models[END_REF].

While conditions C1, C2, C3 and C5 are quite straight-forward to fulfill (see [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF]), both the construction of a suitable instrument that fulfills C4 and of an optimal filter fulfilling C6 are not trivial in practice. The RIV algorithm involves an iterative (or relaxation) algorithm in which, at each iteration, an 'auxiliary model' is used to generate the instrumental variables (which guarantees C2), as well as the associated prefilters. This auxiliary model is based on the parameter estimates obtained at the previous iteration. Consequently, if convergence occurs, C4

and C6 are fulfilled. Thus, the RIV is a suitable method to i) efficiently estimate ρ in [START_REF] Campi | An iterative identification method for linear continuous-time systems[END_REF] when

S o ∈ M and ii) consistently estimate ρ in a practical situation when H o / ∈ H.
Nonetheless, it has to be added that even in the CT LTI case, the convergence of the iterative CT RIV algorithm has not been proven so far and is only empirically assumed [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF]. An additional concern is that even if conditions C1-C6 are fulfilled, assuming that the properties of the approach established in the LTI case apply to the estimation of the reformulated LPV model would mean that the noise-free output terms are a priori known. Therefore, even if the presented method considerably lowers the variance in the estimated parameters, the optimality of the estimates cannot be guaranteed.

C. The LPV-RIVC and LPV-SRIVC (Simplified RIVC) Algorithms

Based on the previous considerations, the iterative scheme of the LPV-Refined Instrumental Variable for Continuous-time models (LPV-RIVC) as well as the simplified version (LPV-SRIVC) can be given in the considered hybrid LPV framework.

1) The LPV-RIVC Algorithm: The following algorithm is designed for hybrid LPV-BJ models.

Algorithm 1 (LPV-RIVC):

Step 1 The usual initialization for CT-RIV algorithm is a DT model estimate issued from an LS method or a DT-RIV algorithm. In the LPV case however, the transformation

(ρ (0) ) ⊤ (η (0) ) ⊤ ] ⊤ is given. Set τ = 0.
Step 2 Compute an estimate of χ(t k ) via numerical approximation of

A(p t , d , ρ(τ) ) χ(t) = B(p t , d , ρ(τ) )u(t),
where ρ(τ) is estimated in the previous iteration. Based on M θ(τ) , deduce χ(t k ) which is bounded according to Assumption A5. Moreover in terms of Assumption A3, χ(t)

is not correlated to the noise.

Step 3 Compute the estimated continuous-time filter Qc (d , ρ(τ) ):

Qc (d , ρ(τ) ) = 1 F (d , ρ(τ) ) , ( 36
)
where F (d , ρ(τ) ) is as given in [START_REF] Rao | Numerical illustrations of the relevance of direct continuous-time model identification[END_REF].

Step 4 Use the CT filter Qc (d , ρ(τ) ) as well as χ(t k ) in order to generate the estimates of the derivatives which are needed later to construct the regressor:

{ Qc (d , ρ(τ) ) χ(i) } na-1 i=0 , { Qc (d , ρ(τ) )ŷ (i) } na-1 i=0 , { Qc (d , ρ(τ) )u j,l (t k )} n b ,n β j=0,l=0 , { Qc (d , ρ(τ) ) χi,l (t k )} na,nα i=1,l=0 .
Step 5 Compute the estimated discrete-time filter:

Qd (q -1 , η(τ) ) = D(q -1 , η(τ) ) C(q -1 , η(τ) ) ,
Step 6 The needed filtered signals {u f j,l (t k )} n b ,n β j=0,l=0 , y f (t k ) and {χ f i,l (t k )} na,nα i=1,l=0 are computed by applying the DT filter Qd on the estimated derivatives obtained in Step 4.

Step 7 Build the filtered estimated regressor φf (t k ) and, in terms of C4, the filtered instrument ζf (t k ) as:

φf (t k ) = -y (na-1) f (t k ) . . . -y f (t k ) -χf 1,1 (t k ) . . . -χf na,nα (t k ) u f 0,0 (t k ) . . . u f n b ,n β (t k ) ⊤ , ζf (t k ) = - χ(na-1) f (t k-1 ) . . . -χf (t k ) -χf 1,1 (t k ) . . . -χf na,nα (t k ) u f 0,0 (t k ) . . . u f n b ,n β (t k ) ⊤ .
February 2, 2011 DRAFT Note that the generation of χ(t) using Step 2 guarantees that the instrument is not correlated to the noise and therefore implies C2.

Step 8 The IV optimization problem can now be stated in the form

ρ(τ+1) (N ) = arg min ρ∈R nρ 1 N N k=1 ζf (t k ) φ⊤ f (t k ) ρ - 1 N N k=1 ζf (t k )y (na) f (t k ) 2 . ( 37
)
where the solution is obtained as

ρ(τ+1) (N ) = N k=1 ζf (t k ) φ⊤ f (t k ) -1 N k=1 ζf (t k )y (na) f (t k ).
The resulting ρ(τ+1) (N ) is the IV estimate of the process model associated parameter vector at iteration τ + 1 based on the prefiltered input/output data.

Step 9 An estimate of the noise signal v is obtained as

v(t k ) = y(t k ) -χ(t k , ρ(τ) ). (38)
Based on v, the estimation of the noise model parameter vector η(τ+1) follows, using in this case the ARMA estimation algorithm of the MATLAB identification toolbox (an IV approach can also be used for this purpose, see [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF]).

Step 10 If θ (τ +1) has converged or the maximum number of iterations is reached, then stop, else increase τ by 1 and go to Step 2.

2) The LPV-SRIVC Algorithm: Based on a similar concept for the estimation of CT LPV-OE models, the so-called simplified LPV-RIVC (LPV-SRIVC) method can also be developed. This method is based on a model structure [START_REF] Tóth | On the discretization of LPV state-space representations[END_REF] with C(q -1 , η) = D(q -1 , η) = 1 and consequently in this case Q d (q -1 , η) = 1. Therefore the LPV-SRIVC algorithm remains the same as the LPV-RIVC algorithm except Step 5, 6 and 10 of Algorithm 1 are skipped. Naturally, the LPV-SRIVC does not lead to the statistically optimal PEM for hybrid LPV-BJ models, however it still leads to consistent estimates. Moreover, a CT LPV-OE models does not involve any DT filtering and consequently, their structure is fully CT unlike for hybrid LPV-BJ models.

IV. SIMULATION EXAMPLE

As a next step, the performance of the proposed algorithms are presented on a representative simulation example. It is important to note that to the best of the authors' knowledge, the presented method is the first approach able to handle the case of colored output measurement noise for CT LPV models. Therefore, the results obtained for the presented algorithm cannot be compared to other algorithms.

February 2, 2011 DRAFT

A. Data generating system

The system taken into consideration is inspired by a benchmark example proposed by Rao and Garnier in [START_REF] Rao | Numerical illustrations of the relevance of direct continuous-time model identification[END_REF]. It has been widely used since then to demonstrate the performance of direct continuous-time identification methods [START_REF] Rao | Numerical illustrations of the relevance of direct continuous-time model identification[END_REF]- [START_REF] Laurain | Refined instrumental variable methods for identification of Hammerstein continuous-time Box-Jenkins models[END_REF], [START_REF] Ljung | Initialisation aspects for subspace and output-error identification methods[END_REF], [START_REF] Campi | An iterative identification method for linear continuous-time systems[END_REF]. In order to create a CT LPV system on which the strength of direct CT identification can be demonstrated, a "moving pole" is considered. A particular feature of LPV systems is that they have an LTI representation for every constant trajectory of p. Such an LTI representation describes the so-called frozen behavior of the system and can be expressed in a transfer function form. In terms of the frozen concept, the "moving pole" means that a particular pole of these frozen transfer functions of S o is a function of p. This phenomenon often occurs in mechatronic applications such as for instance, wafer scanners [START_REF] Wassink | LPV control for a wafer stage: Beyond the theoretical solution[END_REF]. In our case, the Rao-Garnier benchmark inspired "moving pole" LPV system is a fourth order system with non-minimum phase frozen dynamics and a p-dependent complex pole pair. It is defined as follows:

S o                    A o (d , p) = d 4 + (2ζ 2 ω 2 (p) + 2ζ 1 ω 1 ) d 3 + (ω 2 1 + ω 2 2 (p) + 4ζ 2 ζ 1 ω 2 (p)ω 1 ) d 2 + (2ζ 2 ω 2 (p)ω 2 1 + 2ζ 1 ω 2 2 (p)ω 1 )d + ω 2 2 (p)ω 2 1 B o (d , p) = -T ω 2 2 (p)ω 2 1 d + ω 2 2 (p)ω 2 1 H o (q) = 1 1 -q -1 + 0.2q -2 (39)
where T = 4 [s], ω 1 = 20 [rad/s], ζ 1 = 0.1, ζ 2 = 0.5. The slow frozen mode ω 2 is p-dependent and chosen as: ω 2 = 2 + 0.5p. Notice that the frozen behavior (p is fixed to a constant trajectory) of S o for p = 0 corresponds exactly to the Rao-Garnier benchmark defined as

G RG (d ) = -T d + 1 ( d 2 ω 2 1 + 2ζ 1 d ω 1 + 1)( d 2 ω 2 2 (0) + 2ζ 2 d ω 2 (0) + 1) . ( 40
)
Using the given numerical values, S o takes the following form

S o              A o (d , p) = d 4 + a o 1 (p)d 3 + a o 2 (p)d 2 + a o 3 (p)d + a o 4 (p) B o (d , p) = b o 0 (p)d + b o 1 (p) H o (q) = 1 1 -q -1 + 0.2q -2 (
b o 0 (p) = -6400 -3200p -400p 2 , b o 1 (p) = 1600 + 800p + 100p 2 . ( 42c
)
The Bode plot of 20 frozen behaviors of S o is depicted in Figure 1 for 20 fixed values of the scheduling variable p equally distributed from -1 to 1 where the consequence of the moving low frequency mode can be clearly observed. To obtain data records for identification purposes, the input signal u is chosen as a uniformly distributed sequence U(-1, 1), while the scheduling variable is chosen as p(t) = sin(πt).

Furthermore, the sampling period is chosen as T s = 1ms, and the simulation time is T max = 10s which, considering the currently available acquisition possibilities, is a fair assumption.

B. Model structures

In the sequel both the LPV-RIVC and the LPV-SRIVC algorithms are studied for the identification of the data generating system S o . The proposed LPV-RIVC method is applicable to the hybrid LPV-BJ model and assumes the following model structure:

M LPV-RIVC              A(d , p) = d 4 + a 1 (p)d 3 + a 2 (p)d 2 + a 3 (p)d + a 4 (p) B(d , p) = b 0 (p)d + b 1 (p) H(q) = 1 1 + d 1 q -1 + d 2 q -2
(p) = b 0,0 + b 0,1 p + b 0,2 p 2 , b 1 (p) = b 1,0 + b 1,1 p + b 1,2 p 2 . ( 43c
)
while the LPV-SRIVC method is applicable to CT LPV-OE models and assumes the following model structure:

M LPV-SRIVC            A(d , p) = d 4 + a 1 (p)d 3 + a 2 (p)d 2 + a 3 (p)d + a 4 (p) B(d , p) = b 0 (p)d + b 1 (p) H(q) = 1 (
and P x is the power of signal x. The MC results obtained using both algorithms are presented in Table I. It can be clearly seen that the estimates are unbiased which conforms to the theory.

The standard deviation for the nominal part of the coefficients (a * ,0 (p) and b * ,0 (p)) remains low whereas it raises considerably for coefficients a * ,2 (p) and b * ,2 (p).

In Figure 2, the 200 simulated model outputs are plotted together with validation data set (generated under the same excitation conditions as the one used for estimation). It can be seen that despite the large variance in the estimated parameters a * ,2 (p) and b * ,2 (p), the simulated outputs remain close to the true noise-free output signal (considering the level of noise corresponding to a SNR = 20dB). It appears therefore, that these parameter values have a low contribution to the observed output signal under these excitation conditions. Thus, in terms of minimization of the squared prediction error, their role is less significant which results in a relatively large variance of their estimates under noisy conditions. This fact underlines that experiment design is needed to minimize the variance of the parameters estimates. However, there is a lack of input design methods suited for the identification of LPV systems (see [START_REF] Laurain | Refined instrumental variable methods for identification of LPV outputerror and Box-Jenkins models[END_REF], [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF] and references therein) as the concept of persistency of excitation is not well understood yet for LPV models. A procedure for input design for CT-LTI models such as in [START_REF] Campi | An iterative identification method for linear continuous-time systems[END_REF] might be suitable to solve the estimation problem of the a * ,2 (p) and b * ,2 (p) parameters but optimal input design is not investigated for the considered system here.

C. Comments

Note that trying to identify a discrete-time model for the considered CT system is quite a tedious task and requires the following critical issues to be addressed:

• In the LTI case, it has been shown that direct CT identification methods are better suited for the identification of the Rao-Garnier benchmark [START_REF]Identification of continuous-time systems: direct or indirect ?[END_REF]. In the presented LPV context, this means that any frozen behavior of S o (for any fixed trajectory of p) is better identified using a direct CT identification method.

• Using the trapezoidal integration method [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF], the system described in (1) is discretized as:

y(t k ) = 2 -
By analyzing this simple discretization scheme, it is clear that applying it to (39) results in i) the augmentation of the number of parameters to be estimated, ii) non-linearin-the-parameters dependence on p some dynamic dependence on p (appearance of p(t k-1 ), p(t k-2 ) . . . terms).

Estimating such a DT model, where the dependences on p are non-linear-in-the-parameters, is hardly feasible using the existing parametric methods as it requires a customized nonlinear optimization approach. Alternatively, non-parametric methods can be applied or the model can be approximated by a simplified model structure. The latter relaxes the assumption G o ∈ G: in other words, the time-varying property of the LPV system (in this case, a "moving pole") would not anymore be directly linked to the coefficients of the considered DT model and consequently, could not be clearly identified.

V. CONCLUSION

Due to the lack of methods dedicated to the case of direct continuous-time identification of LPV models, there exists a clear gap between the available identification approaches and the practical needs of control synthesis. In order to bridge this gap, a novel method has been proposed in this paper for the identification of hybrid LPV-BJ models and CT LPV-OE models with a p-independent noise process.

Fig. 1 .

 1 Fig. 1. Bode plot of the frozen behaviors of the true LPV system for 20 values of the scheduling variable p (between -1 to 1)

  44) with a 1 (p), a 2 (p), a 3 (p), a 4 (p), b 0 (p), b 1 (p) as given in (43a-f). Note that to demonstrate the achievable performance with the proposed methods we assume that information about the plant in terms of model order and structural dependency is known priori. In terms of identification, the model M LPV-RIVC corresponds to the case S o ∈ M while the model M LPV-SRIVC corresponds to the more realistic practical assumption G o ∈ G and H o / ∈ H. Therefore, 17 parameters are to be estimated by the LPV-SRIVC algorithm and 19 by the LPV-RIVC algorithm. To provide representative results, a Monte Carlo (MCs) simulation of N MC = 200 random realizations is used with a Signal-to-Noise Ratio (SNR) of 20dB where SNR = 10 log P χo P vo ,

Consequently, the qualityFig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. 200 simulated model outputs together with a noise-free and a single MCs output

  The presented method is based on a particular MISO-LTI reformulation of the data equations which enables the use of Refined IV-based methods for LPV IO models in the error-predictionminimization framework. The proposed algorithm has been tested on a representative numerical simulation example inspired by the Rao-Garnier benchmark. The presented example has shown February 2, 2011 DRAFT

  t. time, i.e. d = d dt , p : R → P is the scheduling variable with p t = p(t), χ o is the noise-free output, v o is a quasi stationary noise process with bounded spectral density and it is uncorrelated to p. A o and B

o are polynomials in d with coefficients a o i and b o i that are meromorphic functions 1 of p with no singularity on P A o (p t , d ) = d na + na i=1

  DT model into a CT model is not trivial. Consequently, the initial estimate proposed for the LPV-RIVC algorithm is an LTI-RIVC estimate of M θ , i.e. θ(0) = [

	of a
	February 2, 2011	DRAFT

  (p) = a 1,0 + a 1,1 p, . a 2 (p) = a 2,0 + a 2,1 p + a 2,2 p 2 ,(43a)a 3 (p) = a 3,0 + a 3,1 p + a 3,2 p 2 , a 4 (p) = a 4,0 + a 4,1 p + a 4,2 p 2 ,

	a 1 (43b)
	b 0
	where
	February 2, 2011	DRAFT

TABLE I MONTE

 I CARLO SIMULATION RESULTS WITH ADDITIVE COLORED MEASUREMENT NOISE FOR SN R = 20 DB

			LPV-RIVC	LPV-SRIVC
	Name True Value	mean	st. dev.	mean	st. dev.
	a 1,0	5	4.99	0.058	4.99	0.059
	a 1,1	0.25	0.249	0.081	0.249	0.082
	a 2,0	408	407.94	1.27	407.95	1.29
	a 2,1	3	2.93	1.47	2.92	1.47
	a 2,2	0.25	0.268	2.35	0.251	2.37
	a 3,0	416	415.64	14.22	415.61	14.19
	a 3,1	108	107.72	11.20	107.76	11.13
	a 3,2	1	1.14	30.62	1.17	30.40
	a 4,0	1600	1598.9	44.58	1599.1	44.67
	a 4,1	800	799.31	27.22	799.38	27.13
	a 4,2	100	101.31	79.85	101.08	80.01
	b 0,0	-6400	-6396.4	47.95	-6396.4	48.03
	b 0,1	-3200	-3195.2	65.72	-3195.2	65.76
	b 0,2	-400	-398.92	70.40	-399.14	70.30
	b 1,0	1600	1593.8	366.9	1593.8	365.7
	b 1,1	800	805.02 219.77	802.43	221.7
	b 1,2	100	104.86 749.61	102.95	743.9
	d 1	-1	-0.999	0.0096	X	X
	d 2	0.2	0.202	0.0094	X	X

  T s p(t k-1 ) 2 + T s p(t k ) y(t k-1 ) + T s b 2 + T s p(t k ) u(t k ) + T s b 2 + T s p(t k ) u(t k-1).

A function f is called meromorphic if f = g h where g, h are holomorphic (analytic) functions and h is not the zero function.February

2, 2011 DRAFT

that the proposed procedure is robust to noise and can reasonably well estimate the system in case of an imperfect noise model. Furthermore, it was motivated that in the given LPV framework and for relatively complicated systems, a direct CT estimation method is an attractive approach for capturing the true time-varying nature of the studied system. In this paper, only the case of p-independent noise models has been investigated. Even if refined IV-based methods are theoretically unbiased for p-dependent noise models, the investigation of scheduling dependent noise models remain as a topic for future research.


Stanislas Rohart

Florent Tournus

Véronique Dupuis

Chemical order and size effects on the magnetic anisotropy of FePt and CoPt nanoparticles

Keywords: numbers: 75.75.+a, 75.30.Gw, 75.50.Bb, 75.50.Cc

We investigate the consequence of the dimension reduction on the magnetic anisotropy of FePt and CoPt nanoparticles. Using an extension of the magnetic anisotropy model of Néel, we show that, due to a statistical finite size effect, chemically disordered clusters can display a magnetic anisotropy energy (MAE) as high as 0.5 × 10 6 J/m 3 , more than one order of magnitude higher than the bulk MAE. Concerning L10 ordered clusters, we show that the surface induces a reduction of the MAE as compared to the bulk, due to the symmetry breaking at the cluster surface, which modifies the chemical order.

I. INTRODUCTION

The understanding and the control of the nanostructure magnetic properties is a challenging field of research, for both fundamental and applied physics. At the nanoscale, finite size effects and surface contributions become dominant and lead to a modification of the properties as compared to the bulk, such as the enhancement of the orbital magnetic moment and of the magnetic anisotropy energy (MAE) at the less coordinated atoms. [1][2][3][4] However, technological applications of such structures for ultrahigh density magnetic storage media are limited by thermal fluctuations of the magnetization direction. To overcome this so-called superparamagnetic limit, intense efforts are devoted to the investigation of materials with a high MAE. In such a view FePt [5][6][7][8][9][10] and CoPt [11][12][13][14][15][16][17] have become some of the most studied materials as they display a very large magnetocrystalline anisotropy (over 5×10 6 J/m 3 ) when ordered in the L1 0 phase.

The effect of the L1 0 ordering in FePt and CoPt on the magnetic anisotropy in bulk crystals or thin films has been widely studied on both experimental and theoretical points of view. [18][19][20][21][22] In the case of nanostructures, as the surface modifies the local chemical environment, it is expected to have potentially a large impact on the MAE. Although many experimental studies had been performed [5][6][7][8][9][10][11][12][13][14][15][16] , it is only recently that a clear impact of the dimension reduction on the MAE has been experimentally evidenced for CoPt nanoparticles. 17 By using model samples, and combining various techniques, the authors have indeed been able to reliably determine the MAE of CoPt clusters with a diameter around 3 nm: while the anisotropy of chemically disordered particles was found to be higher than that of the bulk phase, it was the opposite for chemically ordered particles.

In this paper, we use an extension of the phenomenological anisotropy model of Néel to investigate the effect of chemical ordering and size reduction in both FePt and CoPt clusters. We first show that this simple model is in good agreement with bulk and thin film results. Then, taking advantage of its simplicity, which allows us to perform an extensive study of FePt and CoPt nanoparticles, we use this model to investigate the MAE of clusters. Taking into account finite size effects by mean of a statistical analysis on a large number of nanostructures, we demonstrate the impact of L1 0 chemical ordering on the MAE and show how the surface and the finite size modify the results as compared to the bulk. Contrary to the results on pure Fe and Co clusters 4,23 , here the surface induces a lowering of the MAE for L1 0 ordered clusters.

Using a simple phenomenological pair-interaction model to calculate an effect as subtle as the magnetic anisotropy may seem a quite hazardous approach as compared to more elaborate models, like ab-initio calculations. It is indeed well known that a theoretical determination of the MAE is a hard task: because it corresponds to very slight variations in the total energy as a function of the magnetization direction, highly converged results are required. Moreover, there is nowadays no available derivation of the Néel model from first principles. Nevertheless, this model, which is based on symmetries, has been found to be in good agreement with experiments on magnetic nanostructures 4,23,24 . In particular it can remarkably well account for the prevailing effect of lowcoordinated atoms on the MAE. This let us think that this empirical model must be valid to a certain point, and that it can give reliable tendencies regarding the evolution of the MAE in nanostructures with various parameters (size, shape, configuration...). Besides, let us note that a Néel expression for the local anisotropy is also used in atomic scale magnetic simulations based on a Heisenberg-type model. 25,26 . Unfortunately, even it can be quite successfully compared to some ab-initio calculations and experiments (see further in the paper), we cannot demonstrate the validity of the Néel model in the precise case under study: there is indeed no available abinitio calculations of the MAE for comparable CoPt or FePt nanoparticles. However, in order to study the effect of size reduction and of chemical ordering in particles of a few nanometers (which is precisely the size of interest for ultra-high density magnetic storage applications), there is no other method that can reasonably be used, and that would be more precise. We are then convinced that investigating such complex systems, which are out of reach for first-principle calculations for the moment, will give valuable physical insights. Knowing the results predicted by this simple and crude model will then enable further comparison with more advanced calculations in the future. It is also important to check if the predictions of the Néel model are correct in the case of alloy nanoparticles, and we hope that our work will stimulate experimental studies.

II. MAGNETIC ANISOTROPY OF A BIMETALLIC ALLOY

The magnetocrystalline anisotropy (MCA) arises from the coupling of the magnetization with the crystal lattice through the spin-orbit interaction 27 . Then, for a given magnetic atom, it is directly linked to its local environment. For a whole structure, it is the signature of the system symmetries. It can be described using a phenomenological pair model, first introduced by L. Néel 28 . In its original form, for a single element crystal, the MCA energy is written as a pair interaction over the first neighbors. The magnetic anisotropy energy of a dimer is expressed, using the lowest order term as 28

E pair = -L (m • e) 2 ( 1
)
where m is the magnetization unit vector and e the unit vector along the bond. The prefactor L is the anisotropy parameter, which has to be determined and is generally fitted on bulk experiments.

To consider a given structure, we need two more hypothesis. First we suppose that the magnetization direction is homogeneous over the whole structure (macrospin approximation), which is a quite suitable approximation for nanostructures with dimensions close to the exchange length. [START_REF] Hubert | Magnetic domains[END_REF] Second, we also suppose that the L parameter is the same in the bulk and at the cluster surface, and can be extracted from bulk measurements. Even if one could envisage a variation of the anisotropy parameter with the different environment in a structure, in fact the crude model of Néel with a uniform L extracted from bulk measurements has already shown a good agreement (qualitative at least but also quantitative) for such systems. 4,23,24 As the purpose of our study is the prediction of the main tendencies, we will make this assumption in order to restrict the number of parameters in the model. On the basis of these two hypothesis, the total MCA energy, E a , is evaluated by performing a sum over all the nearest neighbor pairs. This leads to a quadratic form, that can be diagonalized into the orthonormal basis set (e 1 , e 2 , e 3 ). The normalization condition on m further simplifies the expression to a biaxial form

E a /V = -K 1 (m • e 1 ) 2 + K 2 (m • e 2 ) 2 + K (2)
where V is the volume of the structure, K 1 and K 2 are positive constants and K a constant that can be omitted. Note that, in general cases, the new basis is different from the one defining the atomic lattice. The direction e 1 (respectively e 2 ) corresponds to the easy (respectively hard) magnetization direction, with the anisotropy constant K 1 (respectively K 2 ). This shows that this model, which is limited to second degree terms in m, always leads to a biaxial second order magnetic anisotropy displaying two degenerated energy minima. At the equilibrium (at zero temperature), the magnetization lies along e 1 and K 1 V corresponds to the minimum energy needed to switch the magnetization from one minimum to the other one. Note that the value of K 2 does not influence this switching energy. For a bulk cubic crystal, this model would give a zero anisotropy, as the second order expression is not compatible with the cubic symmetry (for a more accurate description of the MCA in this case, one should use higher order terms in eq. 1, see ref. 24). However, this crude model is suitable to describe surface anisotropy where the cubic symmetry is locally broken.

To extend this model to an alloy containing two types of atom, named A and B in the following, three types of neighbor pair have to be considered: 22 two between identical atoms and one between A and B atoms. We introduce three different anisotropy parameters: L AA , L BB and L AB respectively for A-A, B-B and A-B pairs. The MCA energy is then evaluated by summing eq. 1 over all the nearest neighbors pairs, taking into account the chemical nature of the two atoms.

This model allows us to describe the effect of L1 0 chemical ordering on the MCA. Neglecting any tetragonal distortion, the atomic site lattice considered is face centered cubic (FCC), and we focus on equiatomic alloys. In this case, the L1 0 phase corresponds to an alternation of pure A and pure B planes in the [001] direction, leading to a chemical anisotropy along the c axis. Oppositely, in the chemically disordered phase named A1, the atoms are randomly distributed on the lattice sites. The magnitude of the chemical order is usually quantified by the long-range order (LRO) parameter S, which varies between 0 (totally disordered alloy) and 1 (perfectly ordered alloy) [START_REF] Nix | [END_REF] . Defining two sublattices labelled α and β, corresponding to successive planes in the [001] direction, S is given by

S = 2(n α -1/2) ( 3
)
where n α is the fraction of α sites occupied by atoms A.

Note that changing the sign of S is equivalent to inverting the convention for the choice of the α and β planes. As this convention is arbitrary, the sign is meaningless and only the absolute value of S has to be considered. Conversely, the probability p for each lattice site to be occupied by an atom of the correct type (A in a α plane and B in a β plane) is directly related to S as p = n α .

To calculate the MCA of a bulk crystal of given order parameter, we consider the probability for a nearest neighbors pair to be a A-A, B-B or A-B bond and we calculate its average magnetic anisotropy parameter L . Three different nearest neighbors pairs with equivalent L are obtained: those lying in the α planes, those lying in the β planes and those linking atoms from two different planes. The corresponding L parameters are respectively L α , L β , and L αβ :

L α = p 2 L AA + (1 -p) 2 L BB + 2p(1 -p)L AB L β = (1 -p) 2 L AA + p 2 L BB + 2p(1 -p)L AB L αβ = p(1 -p)L AA + p(1 -p)L BB + [p 2 + (1 -p) 2 ]L AB .
For a chemically disordered crystal (p = 1/2), the probability to have a A or B atom on a given site is identical, and L α = L β = L αβ : the cubic symmetry is recovered and the MCA is zero. As soon as a chemical order arises in the crystal (S > 0), the cubic symmetry is broken and the second order MCA term has a finite value. Summing the different bond contributions in a FCC unit cell of volume V cell , the MAE is obtained as

E a,cell = 2 (L α + L β -2L αβ ) (m • e c ) 2 = 2 (L AA + L BB -2L AB ) S 2 (m • e c ) 2 , (4)
with e c the unit vector aligned along the c axis of the crystal. It corresponds to a uniaxial anisotropy, along the [001] direction (c axis), and coincides with the chemical anisotropy direction. The anisotropy constant for the bulk alloy is then given by

K(S) = -2(L AA + L BB -2L AB )S 2 /V cell (5)
and has a simple quadratic form. This evolution can be understood looking back at our model: as a pair interaction is considered, the model involves products of the occupation probabilities for a given site, which are linear in S. If 2L AB > L AA + L BB , the anisotropy constant is positive, meaning that the c axis is the easy magnetization axis with the anisotropy constant K 1 = K(S).

We now apply this model to the case of XPt crystals, where X = Fe or Co. In these alloys, a strong second order MCA accompanies the chemical ordering into the L1 0 phase, and originates both from a local chemical environment change and from a tetragonal distortion of the FCC lattice. However, this slight deformation has a much weaker influence on the MCA 21 and we neglect it in the following, as it was supposed in the previous development. In order to apply our model, we need to determine the parameters L XX , L PtPt and L XPt . As the magnetic moment on the Pt atoms is small in these crystals 8,18 , we neglect the magnetic anisotropy of the Pt-Pt bonds by taking L PtPt = 0. A straightforward relation is then obtained between L XX and L XPt by fitting equation 5 on the experimental MAE value of the L1 0 ordered crystal. Consequently, there only remains one adjustable parameter and we finally estimate L XX from measurements on pure samples using the original Néel model 23 : taking into account the magnetostriction parameters, the model gives L XX (r) = L XX (r 0 ) + (r -r 0 ) dLXX dr (r 0 ), where r and r 0 are respectively the bond lengths in the alloy and in the pure crystal. We have determined the parameters for both FePt and CoPt (see table I). The L XX parameter for X = Fe and X = Co has been calculated using the values given in ref. 23

E a Σ = 1 s at 3L CoPt - 3 2 L CoCo = 2K S ( 6
)
where Σ is the film surface and s -1 at is the surface atomic density. K S is the usual surface anisotropy at the Co/Pt interface 32 (the factor 2 on the right hand side of eq. 6 accounts for the two Co/Pt interfaces in the film). Using the values of table I, we find with our model K S = 1.17 mJ/m 2 , close to the value of 1.15 mJ/m 2 reported in ref 32 .

III. CHEMICAL ORDER IN NANOPARTICLES

In order to investigate the relationship between chemical order and magnetic anisotropy, we first need to precisely describe chemical order in nanostructures. Indeed, due to the finite size (and consequently the surface), the description is not as simple as in the bulk. First, the numbers of α and β sites are not necessarily equal like in the bulk. The order parameter has to be adapted, so that a particle with a perfect alternation of A and B planes displays S = 1. Using f α the proportion of α sites, we use a more general definition of S: 30

S = n α -x 1 -f α (7)
where n α is, as before, the fraction of α sites occupied by a A atom, and x is the atomic concentration of A. Whereas for the bulk f α = 1/2 and eq. 3 is recovered, it varies depending on the shape and the size of a nanoparticle. Moreover, note that inverting the convention for α and β planes changes f α (for example, the truncated octahedron cluster, as the one represented in fig. 2, would display f α = 0.507 or 0.493 depending on the convention), but of course do not change the absolute value of S. Besides, a physical inversion of the atomic type on each lattice site will change S into -S, meaning that the physical chemical order is unchanged. However, in this case, the inversion leads to a cluster with a different stacking (in particular, the surface atoms are not the same) and possibly with different magnetic properties, as it will be discussed later.

To obtain a complementary description of the chemical order we can also use a short-range order (SRO) parameter, since for small size systems short-range ordering may dominate and control the nanoparticle properties. Although the SRO and LRO descriptions are equivalent for an infinite bulk (where there exists a direct relation between p, the probability for each site to be occupied by an atom of the correct type, SRO and LRO parameters 33 ) it is not true in nanoparticles. As our anisotropy model is limited to the first neighbors, we also limit the SRO parameter to the nearest neighbors and use the description of Cowley. 33 With p l,m,n the probability for an atom i on a given site to have, in the (l, m, n) direction, a nearest neighbor of the same kind, a local order parameter can be written as

a l,m,n = 1 - 1 -p l,m,n x
where x is, as before, the atomic concentration in A. In an infinite AB FCC crystal and for a given atom, 12 of these parameters can be calculated. However, regarding the L1 0 ordering, they are all equivalent (note that this This parameter is directly linked to the LRO parameter as σ = S 2 . Note that this quadratic relation is understood by the fact that the LRO description is based on the lattice sites occupation whereas the SRO description is based on nearest-neighbor pair correlations.

In a nanoparticle, the SRO parameter has to be adapted due to the surface, as it was the case for S. For each atom i, we calculate a local order parameter σ i as:

σ i = n i -x i 1 -f i ( 8
)
where n i is the proportion of nearest-neighbors (NN) of the same chemical type as atom i, in the (a, b) plane, x i it the proportion of NN of the same chemical type as atom i, among all its NN, and f i is the proportion of NN in the (a, b) plane among all the NN (if atom i is surrounded by 12 NN like in the bulk, f i = 1/3; if it is at a surface site, f i may vary). The SRO parameter σ is then defined as the average of σ i over the whole atoms. Of course, this definition of σ is identical to the Cowley-type one for the infinite bulk.

An important feature, in the case of nanoparticles, is that the finite size allows us to explicitly take into account the statistical distribution of configurations. Indeed, the atomic arrangement for a given ordering probability p is not unique, resulting in a statistical distribution of the chemical order parameters S and σ. Then, the direct relation between S and p and the one between σ and S, established for an infinite system, are no more valid.

In order to investigate the relation between p, σ and S in nanoparticles, we have randomly generated a great number of nanoparticles with a given value of p. For each cluster, the true value of S and σ is calculated and we then obtain a statistical distribution of chemical order parameters (and of cluster composition), related to the small atom number in the cluster. Note that the choice of the direction defining the α and β sublattices is completely arbitrary, so that the [100] as well as the [010] or the [001] direction can be chosen to compute S and σ. In order to have a physical description of the chemical anisotropy, we always choose the directions that gives, on one hand the highest value of |S|, and on the other hand, the highest value of σ. For reasons explained further, we restrict our study to nanoparticles with a truncated octahedron shape.

We first focus on the chemical order in clusters where each site has an equal probability to be occupied by a A or B atom (i.e. p = 1/2). For each cluster, we have calculated both LRO and SRO parameters. Contrary to the bulk where a unique S and σ correspond to a given p, clusters display distributions of LRO and SRO, which are shown in fig. 3. As explained before, this is related to the small size of the particle and illustrates the distribution of configurations. This effect decreases when the diameter of the cluster increases, as the standard deviation of these parameters scales as the square root of the total number of atoms. Consequently, finite average values are found which means that the clusters cannot be considered as perfectly chemically disordered. The 201 atoms clusters display σ = 0.037 and |S| = 0.094 but note also that some clusters have a non-negligible probability to display chemical order parameter as high as σ = 0.08 or |S| = 0.15. Surprisingly, there is no detectable correlation between LRO and SRO for these clusters. This means that, for these disordered clusters, these parameters probe the chemical order in different and complementary ways and none of them is found to be better than the other one.

In partially ordered clusters (p = 1/2), the situation is quite different. Although a statistical distribution still exists, its width decreases with increasing p, which indicates that the number of configurations decreases. A clear correlation appear between LRO and SRO and their means follow the bulk relations: S ≈ 2p -1 and σ ≈ S 2 .

In conclusion, we insist on the importance of the statistical distribution of configuration for nanoparticles. Clusters generated with a given p present a distribution of chemical order parameters, and conversely, a variety of configurations correspond to clusters with the same order parameter. This means that, for small particles, chemical order parameters are not sufficient to fully characterize the chemical arrangement of atoms in the cluster: it is impossible to sum up an entire chemical configuration with only one or two numerical parameters. A statistical analysis is then essential to correctly investigate the link between chemical order and the physical properties in nanoparticles.

IV. MAGNETIC ANISOTROPY OF NANOPARTICLES

We now use this anisotropy model to investigate the consequence of the size reduction in FePt and CoPt nanoparticles, with respect to their chemical order. We focus on nearly spherical clusters, as it corresponds to the usual shape experimentally observed. Moreover, as the host lattice for the FePt and CoPt crystal is FCC, we have chosen to restrict our study to clusters with a perfect truncated octahedron (TO) shape (see figure 2), which is the equilibrium shape in this case 6,34,35 . The magnetic anisotropy of pure Co TO clusters has already been studied using the Néel model 23 . It has been shown that, due to the symmetries of the TO shape, both surface and volume MCA are zero (the MCA of the facets compensates) and that a surface anisotropy only arises when additional facets are present. Therefore, the MAE of perfect TO clusters, studied with our model will explicitly reveal the effects of size reduction and chemical ordering. Moreover, for the sake of simplicity, and in order to avoid any effect coming from the variation of the magnetic element concentration, we have only considered clusters with a composition around 50 % Pt. The dipolar interaction is also neglected in our study, since it is expected to play a significant role only for non-spherical clusters (shape anisotropy). 36 Concerning the magnetic anisotropy, whereas the bulk FePt or CoPt MCA is only uniaxial, the cluster finite size implies that this is no more true and the more general biaxial form (see eq. 2) is always found. However, since the main point of interest is the minimum energy needed to switch the magnetization from an easy magnetization direction to the opposite one, we restrict our discussion to K 1 .

A. Chemically disordered nanoparticles

We first focus on the case of fully disordered clusters (i.e. when p = 1/2, each site having an equal probability to be occupied by a X or a Pt atom). For each cluster, we have calculated K 1 by summing eq. 1 over all the nearest neighbors bonds. We obtain a statistical distribution which is displayed on figure 4 for different types of clusters. Just like the S and σ distributions, the K 1 distribution is related to the small size of the cluster and its standard deviation scales as the square root of the total number of atoms in the cluster. Note also, that as for |S| and σ, the mean value of K 1 decreases when the cluster size increases, indicating that for a large number of atoms the effect of statistical inhomogeneities is smoothed out. Besides, we have observed that the |S| and K 1 distributions are statistically independent which means that no detectable correlation exists between them. At first sight, this may seem surprising. In fact these two quantities are very different: S probes the long-range chemical order whereas K 1 is sensitive to the local chemical order around the magnetic atoms. A perfectly L1 0 ordered bulk structure with an antiphase boundary (i.e.

an inversion in the X and Pt planes stacking sequence) exemplifies this difference: whereas S will be zero, the magnetic anisotropy will almost be unaffected. A small S is then not incompatible with a significant MCA. Note however that a strong correlation exists between |S| and K 1 when clusters display a higher chemical order, as it will be discussed later.

The local character of the magnetic anisotropy let us expect a stronger link between K 1 and the SRO parameter, because they both depend on NN correlations. Indeed, we obtain a noticeable correlation between K 1 and σ. However, although the correlation coefficient between K 1 and σ is significantly higher than the one between K 1 and S, and moreover increases with the particle size, it remains small. In particular, it corresponds to an insignificant variation of K 1 with σ, as compared to the K 1 statistical distribution (this can be seen in figure 7 for σ ∼ 0). In small disordered particles, the configuration distribution is responsible for a partial "decorrelation" between the SRO and the magnetic anisotropy. Although the SRO parameter is more suited than the LRO one, it still fails to establish a definite link between chemical order and magnetic anisotropy.

We have also investigated the spatial distribution of the cluster easy magnetization axis direction (e 1 vector). Whereas the L1 0 ordering induces an anisotropy axis along the c axis, we do not find such a result for disordered clusters. The e 1 direction distribution satisfies the cluster cubic symmetry and the three a, b and c axes are found to be equivalent. However, the distribution is inhomogeneous and a probability maximum occurs in the [110] direction (and the other equivalent directions), as shown in fig. 5. This means that the easy magnetization axis is rarely oriented towards the {001} or {111} facets of the TO cluster but generally points towards the edges connecting two {111} faces. Moreover, this direction also corresponds to a maximum of the MAE: in the case of the 201 atom FePt cluster presented in fig. 5, the mean MAE is about 0.35 × 10 6 J/m 3 for an easy axis aligned along the [110] direction while it is about 0.20 and 0.25 × 10 6 J/m 3 respectively in the [001] and [111] directions.

B. Chemical ordering

To study the chemical ordering effect in nanoparticles, clusters with different p values were generated, leading in each case to a finite width distribution for |S| and σ (except for p = 1). In this way, a large number of configurations is obtained allowing us to investigate the relation between |S|, σ and K 1 . The results are represented in fig. 6 and7 revealing the statistical nature of the link between S, σ and K 1 . These Figures are analogous to the bulk K(S) curve but also shows the finite width distributions.

Although in the bulk, simple relations have been found between K 1 and S (quadratic variation) or σ (linear vari- ation), the relation is more complex for clusters. For small |S| or σ, the mean K 1 is higher than the one of the bulk and we find the result previously described: versus S, the K 1 distribution is constant (mean value and width) and versus σ, a very small increase of the mean K 1 occurs. In this case, the magnetic anisotropy is dominated by the small size effect, in particular by the statis- tical distribution of the atoms on the lattice sites. Above a threshold (S th or σ th respectively for LRO and SRO parameters), K 1 starts increasing and the distribution narrows (note that the threshold is the same for the two order parameters as σ th = S 2 th ). K 1 is proportional to |S| 2 or to σ, similarly to the bulk, which means that the chemical order is sufficient to induce a significant effect on the MAE. However, the increase is smaller than for the bulk and the perfectly ordered L1 0 cluster magnetic anisotropy appears to be weaker than for the bulk. As it will be discussed later, this is due to a surface effect. The case of CoPt clusters is even more complex, as the magnetic anisotropy displays two distinct limits for perfectly ordered clusters. The difference is 0.32 × 10 6 J/m for the 201 atoms CoPt clusters. This corresponds to two different L1 0 ordered clusters where the Co and Pt atoms have been exchanged. As previously discussed in section III this does not changes the order parameter but, while in one case the top and bottom facets are made of Co atoms, in the other case they are made of Pt atoms, the magnetic properties are different (the MAE is higher in this case). This effect also exists for FePt clusters but is weaker (the energy difference is only 0.2 × 10 4 J/m 3 ) and it is not visible on the figures. Finally, as the cluster size increases, the threshold decreases and the K 1 value for L1 0 ordered clusters increases, meaning that the behavior is getting closer to the bulk one.

Our calculations have also revealed that, as the chemical order increases, the easy magnetization direction gets closer to the c axis. When p is high enough (which corresponds to |S| above the threshold mentioned earlier), the easy axis is contained in a cone around the c axis whose summit angle decreases as p increases, as shown in figure 8.

As a matter of fact, our results do not show any particular difference between the description using a SRO or a LRO parameter. In fact, as long as the chemical order is significant, a coherence is developed over the whole particle (statistically, a significant LRO goes with a sig- The anisotropy cone is defined as the one, around the c axis, containing 90% of the easy magnetization axes.

nificant SRO) and both parameters are highly correlated with the magnetic anisotropy. The SRO description is found to be slightly more convenient only for disordered clusters. On the other hand, the theoretical variation of the magnetic anisotropy with the LRO parameter may be more useful for comparisons with experiments, or with more advanced calculations, 5,9,10,19,21 since the LRO parameter is usually used to quantify the degree of order in CoPt and FePt alloys. Moreover, while there is to our knowledge no reported experimental measurement of a SRO parameter, for CoPt or FePt nanoparticles as those considered in our study, it is possible to determine their individual (using for instance transmission electron microscopy) or their mean LRO parameter. [37][38][39][40] V. DISCUSSION

To apprehend the joint effects of surface and of small size in XPt clusters (X = Co,Fe), we focus on two extreme cases: fully disordered and perfectly L1 0 ordered clusters. Theses systems show a completely different variation of the MCA versus the cluster size, as it can be seen in figure 9. In the limit of large diameters, bulk results are recovered: disordered clusters have a zero MCA and L1 0 ordered clusters display MCA of 7 and 5 × 10 6 J/m 3 , respectively for FePt and CoPt. For a finite diameter, chemically disordered (resp. L1 0 ordered) clusters display a higher (resp. lower) MCA than the bulk one.

For perfectly ordered clusters (S = 1), the variation of K 1 versus the number N of atoms in a cluster can be adjusted by this expression

K 1 (N ) = K V + K a N 1/3 (9)
where K V is the bulk MAE, at the corresponding S, and the second term is a small size correction. This small size correction can be read in terms of a surface magnetic anisotropy because it scales with the surface to volume ratio Σ/V . It could be alternatively written as K S Σ/V , with K S the surface magnetic anisotropy. Accordingly, the total anisotropy energy for a particle can then be written, in an usual way, as the sum of a volume and a surface contribution: E a = K V V + K S Σ. However, in the present case K S is negative, which is very unusual. The negative sign of K S indicates a lowering of the MCA at the smallest sizes, as the chemical order is broken at the surface: due to missing neighbors, the contribution of the XPt bonds is reduced and the total MCA decreases. We find K S = -0.42 and -0.33 mJ/m 2 for the FePt clusters, and -0.32 and -0.23 mJ/m 2 for the CoPt clusters. For each type of clusters, two values for K a (and thus for K S ) are found due to the non-equivalence between the two L1 0 stacking possibilities, which indicates that the surface effects depend on the type of facets: the clusters which have the highest anisotropy have pure Pt upper and lower [001] facets. This small size effect cannot be neglected for the smallest clusters: for less than 2000 atoms (about 3.5 nm in diameter), the decrease is more than 10 % as compared to the bulk. It is worth to mention that the absolute value of this surface magnetic anisotropy is comparable to the highest one reported for pure clusters embedded in various matrices [41][42][43] . For chemically disordered clusters (S = 0) equation 9 cannot be used to account for the variation of K 1 with N . Interestingly, in this case the evolution corresponds to the expression

K 1 (N ) = K V + K b N 1/2 (10)
indicating that in this case, the small size contribution does not originate from a surface effect. Here, the effect scales with the inverse square root of the number of atoms, which is reminiscent of a statistical effect. Indeed, as discussed before, for small disordered clusters the statistical distribution of the various chemical configurations is responsible for a concomitant (but uncorrelated) enhancement of the mean long range order parameter and of the mean anisotropy constant.

In the case of intermediate S values, the evolution of K 1 can be adjusted by a compination of the two preceeding expressions:

K 1 (N ) = K V + K a N 1/2 + K b N 1/3 , (11)
where the two constants K a and K b respectively account for the statistical/configuration effect and the surface effect. The results of the fits are shown in table II for different values of |S|, both for FePt and CoPt clusters. The relative weight of these two contributions gives us information on the origin of the anisotropy variation with size. We find that the magnitude of the surface contribution (i.e. the absolute value of K b ) increases with the

The predictions of the empirical Néel model can be of great help to analyze experimental magnetic measurements. Indeed, according to the present results, we can expect the chemically disordered particles to have a quite high anisotropy (much higher than that of the bulk phase), whereas it should be in fact impossible to obtain small chemically ordered nanoparticles having an intrinsic MCA as high as that of the bulk L1 0 phase. This is consistent with our very recent experimental findings on CoPt nanoparticles. 17 In addition, our calculations show that the anisotropy constant is falling off rapidly when a particle is not fully ordered, which could also be responsible of the experimental observation of a lower MAE than expected. This point may be particularly relevant, since it is predicted that the most favorable configurations for small particles correspond to S < 1. 44 We emphasize that the results reported in this paper have been obtained by completely letting aside any energetic consideration for the various chemical arrangements. Nevertheless, our model could easily be combined with a structural optimization method, either to determine the MAE of a few particular structures predicted to be the most stable ones, or to determine the MAE distribution of a statistical ensemble of structures, taking into account their probability of existence.

Another interesting feature predicted by our model is the significant MAE distribution, even for a given particle size and a fixed chemical order parameter, which is due to the existence of a statistical set of chemical configurations. As discussed elsewhere, 45 this can be one of the major source of MAE distribution in realistic assemblies of clusters (with a distribution of composition, size, shape...). This feature, which is a particularity of magnetic alloys, should be directly visible in experimental measurements on well-defined samples.

VI. CONCLUSION

In this study, we have extended the anisotropy model of Néel 28 to describe the magnetocrystalline anisotropy in mixed FePt and CoPt crystals. This phenomenological model allows us to describe the effect of L1 0 ordering on the MCA and appears to be in good agreement with previously reported results on bulk and thin films. We have used this model to investigate the specific modifications occurring in nanosized clusters of FePt and CoPt with a truncated octahedron shape. In chemically disordered clusters, we have shown that there always exists a partial chemical order, due to the small number of atoms, and that, oppositely to the bulk, the MCA is non zero. On the contrary, for L1 0 ordered clusters, we have shown that surface induces a strong diminution of the magnetic anisotropy, due to the broken L1 0 periodicity at the cluster surface.

Our finding points a lack in this type of ordered alloys as materials for high-density magnetic storage applications, where a high switching energy K eff V is needed to stabilize the magnetization direction. In pure or chemically disordered clusters, the anisotropy constant increases when the size decreases, thanks to a positive surface MAE (see e.g. ref. 3,4,43), which attenuates the size variation of the switching energy. On the contrary, in the L1 0 ordered CoPt and FePt clusters, K eff decreases when the size decreases due to the negative surface MAE, which worthen the size variation of the switching energy. Here, surface is a drawback that may be compensated by playing on appropriate surface effects. Changing the surface anisotropy sign could be achieved using the high interface magnetic anisotropy that adds to surface magnetic anisotropy when the cluster surface is in contact with an appropriate non ferromagnetic material, like Pt 42 or CoO 46 . Surface engineering, by embedding the clusters in the matrix or by synthesizing core-shell clusters with the non ferromagnetic shell around the magnetic core, appears to be complementary approaches to these L1 0 ordered materials, in order to overcome the superparamagnetic limit in nanoscale clusters.

FIG. 1 :

 1 FIG.1: Variation of the anisotropy constant K for a bulk FePt crystal versus the long-range chemical order parameter S. The results of our model (continuous line) are compared to experiments performed by S. Okamoto et al.31 and to firstprinciple calculations performed by J.B. Staunton et al.21

FIG. 2 :

 2 FIG. 2: Sketch of a L10 ordered truncated octahedron cluster containing 2406 atoms. The two colors represent the two types of atoms. Note that the topmost and bottommost planes contain the same kind of atom. Consequently, two types of cluster are possible: either the dark atoms correspond to A-type atoms, or they correspond to B-type atoms. These two clusters display different fα factors, but have the same physical LRO parameter.

FIG. 3 :

 3 FIG. 3: Distribution of the SRO (σ) and LRO (|S|) parameters for disordered clusters containing 201 atoms (circles) and 2406 atoms (triangles).

FIG. 5 :

 5 FIG. 5: (color online) Statistical distribution of the easy magnetization axis for a 201 atom FePt cluster, represented in a 3D view along the [111] direction. The color scale represents the probability density for the easy magnetization axis to point in the corresponding direction, red (resp. blue) areas indicating the most (less) probable directions.

FIG. 7 :

 7 FIG. 7: (color online) Variation of the anisotropy constant K1 versus σ for a 201 atoms FePt (a) and CoPt (b) cluster. The color scale represents in the (σ,K1) plane, the density of probability for a cluster with chemical short range order parameter |σ| to have a MCA K1 (the corresponding color bar is shown in inset). The dashed line corresponds to the bulk behavior.

FIG. 8 :

 8 FIG.8: Variation of the anisotropy cone angle θani (summit half angle) versus p for 201 and 2406 FePt atoms clusters. The anisotropy cone is defined as the one, around the c axis, containing 90% of the easy magnetization axes.

TABLE I :

 I and a bond length of respectively 2.72 Å for Fe and 2.71 Å for Co. In order to Parameters LXX and LXPt of the Néel model for FePt and CoPt, determined as explained in the text from the L10 ordered bulk MAE and the bulk Fe and Co Néel model parameters 23 .

	XPt	KL1 0	LXX	LXPt
		(J/m 3 )	(meV)	(meV)
	FePt	7.10 6	0.23	0.73
	CoPt	5.10 6	0.91	0.89

check the accuracy of our model, we compare some predictions with bulk and thin film results. The variation of bulk FePt MAE versus S is compared to experiments

31

and first principle calculations

21

, and, as it can be seen from fig.

1

, shows a good agreement. The model is also compared to thin film experimental results. The MAE of a Pt/Co/Pt(111) multilayer with a Co thickness of one monolayer is

Acknowledgments

We gratefully acknowledge V. Repain and Y. Nahas (MPQ, Paris, France) for stimulating discussions. We acknowledge support from the European Community (STREP SFINx no. NMP2-CT-2003-505587).

chemical order parameter of the particles, while the behavior of the statistical contribution is less intuitive: K a is first increasing with S and then finally vanishes when S reaches 1. Note also that, since the surface and the configuration contributions to the MAE have opposites signs, there exists a range of S values where there is almost no evolution of K 1 with the particle size.