Pascal Brault
email: pascal.brault@univ-orleans.fr

Plasma deposition of catalytic thin films: Experiments, Applications, Molecular modeling

Keywords: Plasma deposition, catalytic films, nanocluster, fuel cells, photocatalytic TiO 2, gas sensor, catalyst treatment, catalyst functionalization, molecular dynamics simulations, thin film growth

Plasma deposition of catalytic thin films is reviewed in view of highlighting some interesting features. Plasma sputtering, plasma enhanced chemical vapor deposition and plasma enhanced metalorganic chemical vapor deposition and their preferential use in various kind of catalytic films are described. Fuel cell electrodes, gas sensors and photocatalytic films are emphasized as significant applications. As example, magnetron sputtering deposition is successfuly used for growing fuel cell electrodes with high performances. Doping doped TiO2 photocatalysts are deposited using various kind of plasma depending on the expected film morphology. Gas sensors are well designed when using plasma deposition. Plasma treatment of catalysts offers a suitable altertive to thermal treatments. Finally, associated simulations, especially recent progress in molecular dynamics simulations of catalytic films growth are surveyed. This is a suitable way to understand basic mechanisms of catalytic film growth.

Introduction

Heterogeneous catalysis is of primary importance in many chemical and electrochemical processes. To be efficient a solid catalyst is usually a supported assembly of nano-sized particles, eventually dispersed on a support with a large surface area [1,[START_REF]Fuel Cell Catalysis[END_REF]. Due to environmental concerns, both the catalyst usage conditions and catalyst preparation methods will suffer from severe international rules. The development of new and environmentally friendly catalyst preparation methods is a very challenging task in the industry. This should also satisfy the fabrication cost and preparation time reductions and if possible in increasing the catalysis efficiency and selectivity. Plasma technologies have been involved in catalysis with some success [START_REF] Kizling | [END_REF][4][5] and they are ready to be popular. Most of the studies are related to ultra-fine particles production, direct catalyst deposition on a support [6,7] and also assistance of catalytic reactions [START_REF] Khacef | Proceedings of the VIIIth Hakone[END_REF][START_REF] Gorce | [END_REF]. Some applications will directly benefits from advances in plasma technologies, among them fuel cells electrodes, photocatalytic films (mainly TiO 2 and variants) and gas sensors due to the interest of the thin film geometries [10]. Catalyst location control and loading reduction, avoiding expensive noble metal catalysts (e.g. platinum), miniaturization for mobile applications or lab-on-chip, morphology of the catalytic film [11][START_REF] Tiggelaar | TRANSDUCERS -12th International Conference on Solid-State Sensors, Actuators and Microsystems[END_REF][START_REF] Randolph | [END_REF][14][15][16][17] are requirements that can be reached using plasma technologies. This article is thus reporting recent advances in the successful use of low or high pressure non thermal plasmas for catalytic thin film growth.

Plasma experiments

Catalyst chemical synthesis is currently used for producing catalyst on an industrial scale and/or being subject to research studies. The chemical routes include impregnation-reduction method, colloidal route, carbonyl route, microemulsion and electrochemical methods [18].

Due to expected environmental limitations, the use of solvents and other products which are toxic at some level, will be prohibited or limited on an undefined timescale. So searching for alternative "clean" technologies, which can also minimize the thermal budget, is an up-to-date challenging effort.

Plasma technologies are one of the candidates for building clean technological process ready to use on an industrial scale, because they are already available in many companies. The advantages of plasma lie in the ability to generate active species which are chemically (ion or neutral radicals, excited species) and/or physically (ions, electrons) active. Moreover this allows new reaction pathways to be created as well in gas phase as at surfaces (walls, substrates). Basically plasma chemistry will lower reaction thresholds, improve catalyst activity and selectivity. Surface plasma treatments are intended to catalytic thin film deposition and catalytic surface functionalization Plasma sputtering deposition (PSD) and plasma enhanced chemical vapor deposition (PECVD), including plasma enhanced metalorganic chemical vapor deposition (PEMOCVD) and combinations of these techniques are the most preferred techniques as being very attractive. Atmospheric non-thermal plasmas such as atmospheric dielectric barrier discharges are progressively used for catalyst treatment [19]. To a less extent pulsed laser deposition is also considered for catalytic film growth [20] Plasma sputtering deposition is a fairly simple technique. A plasma is created by an electrical breakdown in a chamber filled by a sputtering gas, most often argon, possibly in a combination with a reactive gas. Argon ions generated in the plasma are attracted on one or more catalyst targets (with or without a specially designed magnetic field) and are sputtering catalyst atoms to produce on substrate an elemental or compound (alloy) catalytic film; the substrate is either a flat surface or a porous substrate. Plasma sources are of various kinds: RF, DC, pulsed DC magnetron sputter source, Inductive Coupled Plasma (ICP) source (including helicon and Transformer Coupled Plasma (TCP) sources), RF, DC parallel plate source [START_REF] Grill | Cold plasma in Materials fabrication[END_REF][START_REF] Lieberman | Principles of Plasma Discharges and Material Processing[END_REF]. Main advantage is the control of the sputtered particle flux, energy distributions at the target and at the substrate [START_REF] Thomann | [END_REF][24][25][26] by simply tuning experimental conditions. An interesting feature is the sputtered particle distribution that is non-maxwellian (contrary to vacuum evaporation) with mean kinetic energies in the range 10-30 eV. The transport of sputtered particles through the buffer gas to the substrate is the way to modify these distributions that can become thermalized if the transport path is long enough and/or the pressure is high enough. Moreover high deposition rates (up to 10µm.h -1 ) are possible. Moreover excited species can exists, due collision between energetic sputtered atoms and the buffer gas atoms (e.g. argon), excited or metastable Ar created by electron collisions, which can interact with the growing film [27][28][29][30][START_REF]Reactive Sputter Deposition[END_REF][START_REF] Westwood | Sputter Deposition[END_REF].

PECVD is also a very interesting technique, which allows both deposition and treatments like functionalization and grafting. Ability of PECVD is to deposit films from molecular precursor in gaseous or liquid form at various rates. The basic mechanism is the polymerization in the gaseous phase or at the surface. This occurs due to radical creation by electron collision on precursors which overcome limitations occurring in Chemical Vapor Deposition (CVD).

Moreover, plasma generates radiation, positive ions, neutral and metastable species, electrons that may interact with immersed surfaces [33]. So new reaction pathways are possible and thus « new » materials can be deposited as thin film with easy control of the porosity. A variant of the PECVD for which precursor are metal containing carbon bonds. Precursors are liquid or solid, so heating is often required before injection in the plasma chamber [34,35].

Typical catalytic thin films

The most "simple" catalytic film is a collection of noble metal supported nanoclusters, nanocluster being the most preferred catalyst morphology [1]. Fig. 1 shows Transmission electron Microscopy (TEM) images of of Pd clusters grown onto µc-SiC obtained by ICP sputtering of a Pd wire [6,[START_REF] Thomann | [END_REF]25]. The clusters have a well monodispersed size distribution. These nanoclusters were successfully tested against hydrogenation of 1,3 butadiene [6]. More generally, a parametric study of Pd sputter deposition shows a transition from nanocluster growth [36] to continuous film of very low thickness (< 10nm), which is of potential interest for gas separation membrane [7]. A major requirement in heterogeneous catalysis is a high surface area of the active phase. This can be achieved by decorating spherical particles (e.g. carbon) with catalyst nanoclusters, which are further randomly stacked on suitable support, providing the high area active surface. Another way is to sputter catalyst target onto porous support. As example, for building a proton exchange membrane fuel cell (PEMFC) electrode, it is proposed either to deposit Pt or Pt alloys onto uncatalyzed gas diffusion layer made with Vulcan XC72 carbon particles brushed onto a carbon cloth or paper [37]. Another way is to directly grow an active layer by Ptcarbon co-sputtering [38][39][40].

Another type of catalysts is oxide particles and films. They also are well suited for plasma synthesis. A very popular oxide catalyst is TiO 2 which is especially interesting for photocatalytic applications [41]. Variants are possible with many doping species: metals (Fe, Cr, W, Ni, …), oxides (Fe 2 O 3 , WO 3 ) of light elements (N, C, S, ..) [5]. Main interest is to increase the visible light absorption for benefiting from the solar light and/or to enhance the photocatalytic activity, for example for improving antibacterial effects, pollutant/toxic molecule decomposition, water recycling, … .The anatase phase is preferred in the form of nanoclusters, nanofibers, nanotubes, …, possibly decorated with metal clusters. It can be prepared via numerous routes, especially all plasma techniques can be used and lead to various thin films morphology, those with high specific area being preferentially aimed [42][43][44][45][46][47]. When comparing on Fig. 2, the different morphologies issued from (a) PEMOCVD, (b) Magnetron sputtering and (c) PECVD, we observe that surface disorder is greater in PEMOCVD. Growth is columnar and dense with magnetron deposition and long nanocolumns (nanofibers) are obtained using PECVD. Noble metal oxides are catalysts of hydrocarbon conversion and can also be prepared by plasma sputtering [48], which leads to performance improvement.

The field of carbon nanotubes (single wall and multiwall), nanofibers, nanowalls, and other nanostructures is subject to a huge amount of studies and are partly of interest as catalysts, eg electrochemical devices [49][50][51][52][53]. Deposition of catalysts on nanofibers (see Fig. 3), nanotubes arrays is a growing fields especially as fuel cell electrodes [54][55][56][57][58][59].

A kind of substrate call us to attention: anodic aluminum oxide (AAO) templates [60] (see Fig. 4a). Briefly, when oxidizing aluminum foils (thickness around 50 -100 µm) in an electrochemical bath, a hexagonal array of cylindrical nanopores (nanochannels) is created.

The size and density of nanochannels is depending on bias and current of the electrochemical cell. Such 2D porous array can thus be used as template for various nanowires growth eg platinum, silicon, CdS. The use of plasma sputtering deposition allows for creating a well organized cluster assembly with gradient concentration inside the pores [61] as displayed in Fig. 4b. This could be useful when searching for localized catalytic effects near a boundary.

Applications

Fuel cells

A fuel cell is a power supply that converts electrochemical energy into electricity and heat by the following reactions (for a proton exchange membrane fuel cell, PEMFC): Fuel (H 2 , CH 3 OH, …) carrying hydrogen atoms that are catalytically transformed into protons at the anode side, releasing electrons and supplying an external charge (H 2  2H + + 2 e -), while protons are transported through a ion conducting solid polymer electrolyte. This electrical circuit is closed at the cathode by the catalyst (Pt) activated oxygen reduction with protons to give water and heat. The PEMFC is thus a three layer system. The electrodes (anode and cathode) have an active layer made of dispersed catalyst nanoparticle (sizes in the range 2 to 10 nm) supported on carbon VulcanXC72 particles (size in the range 20 -100nm). This active layer of usually 20 to 50 µm is coated (painted) on gas diffusion layer, eg carbon cloth or carbon paper. The catalyst load (Pt for example) can vary from 0.3 to 2 mg Pt cm -2 per electrode for a delivered power of 0.5W.cm -2 . So a 80 kW motor requires 80g of Pt, which costs 3200 €. It is thus necessary to reduce Pt quantity and increases Pt use. Many attempts using plasma have been undergone. Recently plasma sputter deposition has reached good performances and high use of platinum and Pt quantity reduction down to 0.01 mg Pt cm -2 per electrode, leading to power density of 20 kWg Pt -1 [40]. In that case, Pt nanoclusters of size 0.7-2 nm (Fig. 5) where grown during Pt-Carbon magnetron co-sputtering on carbon Vulcan XC72 electrode and are present over 400 nm inside the active layer. Similar approaches involved Pt sputtering simultaneously with C 2 H 4 PECVD or more generally hydrocarbon PECVD leading to carbon nanostructure decorated with some catalyst [54][55][56][57][58][59]62] While Pt load reduction is a challenge, the goal is also to avoid Pt and more generally noble metal catalysts. Some special molecules called metal (mainly Fe, Cu, Co, …) porphyrins and phtalocyanins lead to promising properties [63,64]. This comes from the metal bonding to four nitrogen atoms in a planar geometry (Fig. 6). Plasma treatments of catalytic films containing such molecules were shown to improve fuel cell performances using low pressure ICP plasma and atmospheric dielectric barrier discharge [65,66].

Gas sensors

Beside fuel cells, an application demanding performances of plasma technologies is the development of sensors, especially gas sensors, as needed for process controls due to new environmental management rules. As fuel cells, gas sensors are multilayered systems eventually based upon catalysis for species detection through catalytically activated reactions [11,14,[67][68][69][70][71] (see Fig. 7)

Specifically

Catalysts treatments

Catalyst efficiencies can be improved by suitable exposure to inert or reactive plasmas.

Especially, catalyst treatments using nonthermal plasmas is becoming a growing field of interest. It is due to the obtained increase in the catalyst acidity (both Brönsted and Lewis acidities), in the dispersion and in the low-temperature activity (which is largely due to the increased dispersion of the catalysts during treatment [5]). The presence of catalyst powder is modifying the high pressure plasma indicating the effectiveness of the process (Fig. 8). Major results concern the increase in the dispersion and the stability of catalyst. Reduction of catalysts is also very active. For example, the use of microwave plasma [72], atmospheric DBD [19] have resulted in higher activity and selectivity, due to change of nanocluster sizes

Molecular modeling

For better understanding of the growth of metal or oxides films, as those needed for catalytic films, it is interesting to carry out molecular simulations at the atomic level. Building efficient catalytic layer needs to understand the growth of catalytic phase, in the form of metallic and oxide clusters or films. Molecular simulations as Molecular Dynamics (MD) simulations are well suited for understanding growth phenomena at the nanoscale. Molecular Dynamics simulations are based upon the resolution of Newton equations of motion. So for an ensemble of N atoms or molecules, the positions, velocities and acceleration are always known [73,[START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF].

If dissipation occurs, through friction terms for example, Langevin-like equations have to be solved. For simplicity, we consider the cases for which Newton's equations of motion are valid. They can be written:

)) ( , ), ( ), ( ( force with the , 1 ) ( 2 1 2 2 t r t r t r V r f f m t t r N i i i i   
The only information necessary to solve this set of N equations of motion is the potential energy. Resolution of these equations is carried out using stable algorithms as Verlet velocity scheme [73]. Statistical information and materials properties can be deduced by averaging over all trajectories and compared directly to experiments. Thus MD is a very powerful method for understanding molecular interactions. For illustrating this, examples of metal cluster growth and oxide films are described.

Pt deposition on model porous carbon substrate.

Pt-Pt interaction potentials (also valuable for transition elements) is given by the second moment approximation of the tight binding theory (TB-SMA) [START_REF] Cleri | [END_REF] and for Pt-C by an efficient Lenard-Jones potential [76]. They take the forms, respectively:

2 1 0 2 0 1 2 exp 1 exp j i ij j i ij i r r q r r p A V and 6 12 4 ) ( j i j i j i ij r r r V
Calculations are performed on various model porous carbon substrates as depicted in Fig. 9.

Main information that can be gathered is the radial density function which inform on the cluster/film structure, depth concentration profile, cluster size distribution if any. As a matter of fact, experiments display a Pt concentration profile [77] that can be fitted with the stretched exponential function ρ(z,t), solution of a generalized diffusion equation:

) ( exp ) ( ) , ( 2 t B z t A t z
. As shown in Fig. 10, the calculated θ using MD simulations is in the range -1.0 to -1.7 as for the experiments [77], indicating enhanced Pt penetration in the porous electrode.

Oxide clusters and films.

MD simulations of oxide clusters or films are a bit more difficult due to long range potential describing ionic interactions. A fairly simple potential (but not necessarily always the most efficient) is given by:

ij ij ij B ij j i ij r q r C r r A r V 6 6 exp ) (
. Numerous strategies exists for solving Newton equations with long range interactions as Ewald summation, … [START_REF] Vashishta | Amorphous Insulators and Semiconductors[END_REF]. Fig. 11 displays the evolution of a TiO 2 cluster under thermal relaxation at 300 K [START_REF] Naicker | [END_REF]. Atomic positions are largely modified compared to the frozen initial condition corresponding to the bulk anatase structure. This very simple potential is also able to reproduce experimental results of plasma sputtering deposition of magnesium aluminium oxide [80][81][82]. Evolution of the film morphology is reproduced in Fig. 12 as function of Mg content [82].

Conclusion.

As a summary of this short review, plasma technologies are very efficient for designing catalyst thin films, especially plasma deposition is very well suited for controlling catalyst content and localisation in a layer, which is out of range of standard chemical means, which always provide homogeneously dispersed catalysts on a thickness not less than 5 to 10 µm.

Numerous successes have been recorded for applications in fuel cells for which plasma deposited electrodes have a very low catalyst content and a very good efficiency.

Photocatalysis (TiO 2 films and its derivatives) benefits from versatility of plasma deposition methods for allowing a wide range of photocatalysts geometry (films, nanotube, complex nanostructure, porous films). The catalytic films in gas sensor technology use plasma deposition in a same way as for microelectronics. Catalyst treatment and functionalization make use of new plasma technology of atmospheric DBD for avoiding thermal treatments. As a consequence, catalyst treatments are softer and more efficient in the way of catalyst efficiency and durability.

From the point of view of green chemistry, the use of plasma at least represents a reduction in the use of chemicals, especially the use of noble metals and precursors. It gives also the possibility also for giving up the use of toxic chemicals, especially by using plasma sputtering deposition. Moreover, the electrical energy consumption for nonthermal plasma catalyst treatments is not high.

Molecular simulations, as Molecular Dynamics, are able to provide relevant insight into basic mechanisms of cluster and film growth leading to the possibility of designing more efficient catalytic cluster assemblies or thin films.

Figure 1 :

 1 Figure 1: TEM pictures of Pd nanocluster growth on SiC by ICP sputter deposition taken different Pd loading (a) 0.1 µg Pd .cm -2 , (b) 1.0 µg Pd .cm -2 and (c) 3.0 µg Pd .cm -2 . Scale bars are 10 nm (from Ref. 6).

Figure 2 :Figure 3 :

 23 Figure 2: TiO 2 photocatalytic thin films deposited by (a) PE-MOCVD (using bis(2-N-(2hydroxy-2-methylethylimino)-4-pentanoate)titanium(IV) (from Ref. 42), (b) reactive magnetron sputtering (from Ref. 43), (c) PECVD (using Ti(OC 3 H 7 ) 4 TIPP) (from Ref. 44),

Figure 4 :

 4 Figure 4: Hexagonal array of cylindrical pores decorated by platinum clusters grown by plasma sputtering (a) top view (b) cross-section. Pt clusters appear as white zones, from Ref. 60.

Figure 5 :

 5 Figure 5: High Resolution TEM micrographs of plasma sputter deposited Pt nanoclusters on carbon particles. Dark areas are Pt clusters. Courtesy: C. Coutanceau, LACCO, Poitiers, France.

Figure 6 :Figure 7 :

 67 Figure 6: Drawing of the metal-porphyrin: cobalt tetramethoxophenylporphyrin (Co-TMPP) atomic structure

Figure 8 : 5 )Figure 9 :Figure 10 :

 85910 Figure 8: Photos of high pressure glow discharge without (a) and with catalyst powder (b). I: negative glow space; II: positive column space (catalyst powder located on the bottom, from Ref. 5)

Figure 11 : 79 Figure 12 :

 117912 Figure 11: 3 nm anatase cluster. (i) initial, (f) and final structures at 300 K. Oxygen and titanium ions are dark and light colors, respectively, from Ref. 79

Fig. 1 22

 1 Fig.1

Fig. 4

 4 Fig. 4

  Fig. 5

Fig. 6 27

 6 Fig. 6

Figure 8 29Fig

 8 Figure 8

Fig 11 32

 32 Fig 11

Fig. 12

 12 Fig. 12

  , catalytic gas sensors for methane (CH 4 ), butane (C 4 H 10 ) and hydrogen (H 2 ) are used mainly for preventing explosive accidents. The operation principle is catalytic oxidation of combustible gases. Thick film catalysts, which are normally noble metals such as platinum and palladium can be heated by a wire coil in an alumina bead. The catalytic combustion of the gases starts when the electrically heating catalyst material reaches the required temperature. The platinum resistance variation relative to a reference resistance called compensator is measured as the temperature due to oxidation of inflammable gases is increased[11,14,[67][68][69]. Sensors for other gases such H 2 S, O 2 are also studied[70,71]

Acknowledgements

I am indebted to people who contributed to the past and present catalytic thin film projects at GREMI or through collaborations or during deep discussions, which finally leads me to write such a review, especially: A.-L. Thomman, A. Caillard, A. Ennajdaoui, H. Rabat, M. Cavarroc, M. Mougenot, C. Coutanceau, S. Baranton, C. Lamy, J. Durand, S. Roualdes, G. Pourcelly, P. Buvat, C. Josserand, F. James, J.-L. Rouet, R. Boswell, C. Charles, D.B. Graves, S. Wu.


C Ducourtioux

M F Ouedraogo

UNIQUENESS OF TRACES ON LOG-POLYHOMOGENEOUS PDOS

Keywords:

We show how to derive the uniqueness of graded or ordinary traces on some algebras of log-polyhomogeneous pseudodifferential operators (PDOs) from the uniqueness of their restriction to classical pseudodifferential ones.

Introduction

We consider a closed connected riemannian manifold M of finite dimension n and a finite rank hermitian vector bundle E over M . A pseudodifferential operator acting on smooth sections of E is called classical (or polyhomogeneous) ( [START_REF] Shubin | Pseudodifferential operators and spectral theory[END_REF]) if locally its symbol is classical, i.e. it admits an asymptotic expansion in positively homogeneous components.

A pseudodifferential operator L acting on smooth sections of E is called log-polyhomogeneous if locally its symbol has the form c XXXX Australian Mathematical Society 0263-6115/XX $A2.00 + 0.00 0 AMS 2000 Mathematics subject classification: 47G30, 58J40.

1 a k (x, ξ) log k |ξ| + a k-1 (x, ξ) log k-1 |ξ| + • • • + a 0 (x, ξ), (1)
where k ∈ N and a 0 , • • • , a k are classical symbols. We call the integer k the log-degree of The assumption on A C implies the following fundamental decomposition of A (see Lemma

?? Paragraph 2) A = +∞ k=0 A C Log k Q. (2)
We assume that A C does not only contain smoothing operators. Otherwise, A C and A are algebras of smoothing operators and V. Guillemin ([4]) has shown that the L 2 trace is the unique trace on such algebras.

On an algebra, we say that a linear form τ is a trace if for any operators A and B in the algebra, we have τ (AB) = τ (BA) i.e. τ vanishes on commutators. On a graded

algebra B = k≥0 B k , a graded trace is a sequence (τ k ) k∈N of linear forms τ k on 0≤l≤k B l which vanishes on 0≤l≤k-1 B l and such that if A is in 0≤l≤k B l and if B is in 0≤l≤m B l then τ k+m (AB) = τ k+m (BA).
Under the assumption of the uniqueness of a trace τ 0 on A C we show that there exists a unique graded trace (τ 0 k ) k∈N on the whole algebra A extending τ 0 (see Theorem ??

Paragraph 3). We also prove that if there exists a trace on A extending τ 0 then this extension is unique (see Theorem ?? Paragraph 3).

Our first result applies to the Wodzicki-Guillemin residue (also called the non commutative residue) res on the algebra C provided M is of dimension n ≥ 2. It is well known that the Wodzicki-Guillemin residue res is the unique trace on C . Setting res 0 = res, the extension res 0 k coincides with the higher non commutative residue res k introduced by M.

Lesch up to a multiplicative factor.

The k-th residue of an operator L in L of log-degree k with local symbol

σ(L) = a k (x, ξ) log k |ξ| + a k-1 (x, ξ) log k-1 |ξ| + • • • + a 0 (x, ξ) is defined as res k (L) = (k + 1)! M S * M tr((a k ) -n (x, ξ))dξdx.
On L seen as a graded algebra (the grading is given by log-degrees), M. Lesch ([7]) has

shown that the sequence (res k ) k∈N is the unique graded trace. We recover the same result by an alternative approach.

Our second result applies to the canonical trace on the algebra of odd-class log-polyhomogeneous PDOs when the manifold M is odd dimensional.

According to M. Kontsevitch and S. Vishik [START_REF] Kontsevich | Geometry of determinants of elliptic operators[END_REF], a classical operator A of order m ∈ Z is of odd-class if locally the positively homogeneous components of its symbol {a m-j : j ∈ Z} are simply homogeneous, i.e. they have the property a m-j (x, -ξ) = (-1) m-j a m-j (x, ξ).

The odd-class classical PDOs form an algebra. Following [START_REF] Paycha | A Laurent expansion for regularised integrals of holomorphic symbols[END_REF] we say that a log-polyhomogeneous PDO L is of odd-class if locally, all the classical symbols a 0 , • • • , a k arising as coefficients of powers of log |ξ| in its symbol (see ( 1)) have the above property (3). Similarly to odd-class classical PDOs, one can easely check that odd-class log-polyhomogeneous PDOs form an algebra.

When M is odd dimensional, the canonical trace TR of [START_REF] Kontsevich | Geometry of determinants of elliptic operators[END_REF] is well defined on the algebra of odd-class classical PDOs. Recently, L. Maniccia, E. Schrohe and J. Seiler ( [START_REF] Maniccea | Uniqueness of the Kontsevich-Vishik trace[END_REF]) have proved that TR is the unique trace on this algebra. The canonical trace has been first extended

L.

  We denote by • L the algebra of log-polyhomogeneous PDOs acting on smooth sections of E • C the subalgebra of classical PDOs of L • Q an admissible classical PDO of positive order such that LogQ exists • A C a subalgebra of C such that [A C , LogQ] ⊂ A C ( in general we only have that the commutator of a classical PDO and a logarithm is a classical PDO, hence we only have [A C , LogQ] ⊂ C ). • A the subalgebra of L generated by LogQ and A C .

by M. Lesch ([7]) to log-polyhomogeneous PDOs of non integer orders (this means that the orders of a 0 , • • • , a k are not integers) and then by S. Paycha and S. Scott to odd-class log-polyhomogeneous PDOs when M is odd dimensional ( [10]). We refer the reader to Paragraph 4 for more details. To our knowledge, the proof of the uniqueness of TR on odd-class log-polyhomogeneous PDOs is new.


Marcel Derrien

Herve Le Gleau

Marie-Paule Raoul

Hervé Le Gléau

The use of the high resolution visible in SAFNWC/MSG cloud mask

NWCSAF (www.nwcsaf.org) consortium continuously develops and maintains a software package to extract various products from Meteosat Second Generation (MSG) imagery. Meteo-France has developed the part of this software package retrieving cloud mask and types, cloud top temperature and height. The improvement of the detection of small scale low clouds was a challenging task asked by the users of the SAFNWC/MSG cloud mask and type, for an earlier detection of convection and for a better screening of cloud-contaminated pixels prior the use of clear-sky radiances in the study of marine or land surface properties or their assimilation in Numerical Weather Prediction models. For this aim Meteo-France has developed an algorithm analysing the High Visible Resolution (HRV) of Spinning Enhanced Visible and Infrared Images (SEVIRI) and implemented it as an option in v2010 of the SAFNWC/MSG cloud mask available for users in June 2010. The paper describes the algorithm based on temporal analysis of HRV features and its application at the coarser spatial resolution of the SEVIRI cloud mask. Real examples illustrate situations where its behaviour can be visually estimated (improvements, misses and false alarms). Finally, the impact of its implementation is evaluated through a 11 month period by comparison with SYNOP and SHIP ground-based cloud covers. The increase of Cramer's coefficient (a normalized χ 2 indicator) when using HRV denotes that it is better associated with the SYNOP and SHIP groundbased observations of the cloudiness categorized in three classes, clear(N≤2), broken (3≤N≤5) or cloudy (N≥6).

INTRODUCTION

A better detection of small scale low clouds at daytime was requested by the users of SAFNWC/MSG package during the SAFNWC users workshop held at Madrid in 2005. The advantage of the better spatial resolution of the HRV channel was not exploited by the previous versions of the SAFNWC/MSG cloud detection scheme. Therefore its use was one of the improvement tasks planned during the second phase of the project. The SEVIRI cloud mask algorithm analysed in this paper is an output of this task.

HRV IMAGERY CHARACTERISTICS

SEVIRI has one high resolution (HRes) broadband solar channel (0.3-1.1 μm) and 11 lower resolution (LRes) channels; 3 narrowband solar channels (0.6, 0.8 and 1.6 μm) and 8 thermal infrared channels (3.9, 6.2, 7.3, 8.7, 9.7, 10.8, 12.0 and 13.4 μm). The spectral response functions of the narrowband channels at 0.6, 0.8 μm and HRV channel are shown in figure 1. The HRV sampling increment angle is 251.53/3 μrad in both N/S and E/W direction, in fact 1/3 of LRes bands, leading to a 1x1km2 at subsatellite point. However, one should note that the true optical resolution is lower by a factor 1.6 for both HRes and LRes (respectively 1.67km and 4.8km at the sub-satellite point, [START_REF] Schmetz | An introduction to Meteosat Second Generation (MSG)[END_REF]. The HRV measurements are performed by 9 detectors, their calibration adjustment, registration between channels, and rectification are performed before the level 1.5 images reach the users through EUMETCast. A full HRV earth disc would be an image of 11136x11136 points. Due to MSG downlink limitations, only two distinct sub-areas are transferred to the users, the upper zone covering Europe has a fixed position in the disc and the lower one is shifted to follow the daily illumination. Between 00h UTC and 14h UTC, its position is suited to monitor the Indian Ocean, between 14hUTC and 17h UTC it is shifted westwards every hour and between 17h UTC and 00h UTC it remains in a position to monitor the Atlantic Ocean. Figure 2 illustrates the transmission scheme of HRV data. Of course this alternative sectorized coverage is a drawback for the users processing regions not permanently covered by HRV, it is increased by the need of a pair of images for the HRV-based detection algorithm. The coregistration of images is known by the position of the pixel whose centre corresponds to 0 E resp. 0 N ; in HRV image it is at position (5566,5566) and (1856,1856) in LRes SEVIRI, in the native coordinates (starting at 1 for the lower right pixel of each image). This adjustment is performed by the rectification process. It is an essential prerequisite for a good fusion of HRV and other LRes SEVIRI bands in an instantaneous image and between consecutive images. Deneke and Roebling, 2010 found a minor misregistration of the HRV relative to the narrowband channels of 0.36km±0.11km in East and 0.06±0.10km in South direction for the European area that they have studied. We have neglected this effect in our scheme.

HRV-BASED DETECTION IMPLEMENTATION DETAILS

We designed it as an add-on option of the cloud mask, as some users may be restricted by hardware considerations or would like to avoid changes in cloud detection efficiency for their operational applications. The optional HRV-based algorithm is applied only to SEVIRI pixels remained as clear after the two previous steps of the cloud mask algorithm (SEVIRI spectral and textural thresholding (Derrien and Legleau, 2005) and temporal-differencing and region-growing, (Derrien and Legleau, 2010). It requires not only the HRV data of the current scene, but also the scene observed and analysed 15 minutes sooner.

HRV reflectance test over land

The assumption is that any low cloud exhibits higher HRV reflectance than underlying clear-sky surface. This test is designed to catch static small scale low clouds that are not detected by change detection technique. A LRes SEVIRI pixel is classified as cloud contaminated if not already classified as cloud by previous tests and if:

• solar elevation > 5° and Max(R hrv 3x3) > R hrv threshold Where Max(R hrv 3x3) is the maximum HRV normalized reflectance in the 3x3 HRV array covered by the corresponding SEVIRI LRes pixel. The threshold is a crude estimate of clear-sky broadband reflectances derived from maps of monthly values and an angular empirical correction. The monthly maps are obtained from a combination of black sky albedo for 3 MODIS narrow bands (0.55μm, 0. 67μm and 0.86μm) available from a NASA website dedicated to MODIS atmosphere (http://modis-atmos.gsfc.nasa.gov/ALBEDO/index.html). Black-sky albedo is the directional hemispherical reflectance computed at local solar noon estimated from MODIS Terra and Aqua measurements. The parameterization of this test has been designed using a one-year period of collocations of SYNOP and satellite data including HRV. When applied to full disk data during summer 2009 false alarms that were not observed on the training data set and too difficult to remove (not clearly depending on Rmap or scattering angle) appeared in some places. To avoid them, the threshold has been increased lately before the software delivery to the integrator. A consequence is that some rather static topographically induced clouds such as for instance valley fog not directly detected by other SEVIRI visible channels may still remain not detected by this test.

Local spatial HRV texture test over sea

A SEVIRI pixel is classified as cloud contaminated if:

• solar elevation > 10° and [SD(R hrv 3x3)/Mean(R hrv 3x3) > .08 or SD(R hrv 3x3) > 0.8%] • solar elevation < 10° and solar elevation > 5° and [SD(R hrv 3x3)/Mean(R hrv 3x3) > .16 or SD(R hrv 3x3) > 0.4%] -Mean(R hrv 3x3) and SD(R hrv 3x3) stand respectively for the mean and standard deviation computed using the 9 reflectances of HRV pixels covered by the SEVIRI low resolution pixel This test based only on HRV local texture inside SEVIRI low resolution pixel detects heterogeneities over sea, one feature detects small clouds when background is rather dark whereas the other one is more efficient for heterogeneities inside clouds or sunglint, that were not detected by previous test based on other SEVIRI visible channels.

Local spatial texture and temporal test over land

A SEVIRI pixel is classified as cloud contaminated if:

• clear (and not snowy) and -SD MIN MAX MEAN stand respectively for the spatial standard deviation, minimum value, maximum value and mean of the feature.

-R hrv is the HRV reflectance -RN hrv is the HRV reflectance normalized by the analytical formulation of solar pahthlength valid for a standard atmosphere proposed by [START_REF] Li | On the effective solar pathlength[END_REF] (in place of classical inverse cosine function).

-indices cur and prev stands for pixels from the current image or the one 15 minutes sooner.

The assumption behind this test is that a mobile or evolving bright target inside a 3x3 HRV array is a cloud. False alarms that have been problematic are clear sides of cloud shadow limits moving over bright grounds. This is why the pixels detected by this test are further processed by a neighbourhood analysis described below.

Neighbourhood filtering

Clear restoral test over land

This test is designed to remove some false alarms passing through the previous HRV-based tests, and is applied only to pixels detected by this test.

A LRes SEVIRI pixel is restored as clear if:

• Its MEAN(R hrv 3x3) is minimum among other land pixels in its SEV 3x3 low resolution neighbourhood

Cloud restoral test over land

Obviously during visual inspection of the results some clouds remained missed after the change detection test. Brighter pixels in the neighbourhood of pixels detected as clouds by the HRV change detection algorithm described in the previous sections should also be detected as cloudy. This test is designed to this purpose.

EXAMPLES OF HRV-BASED ALGORITHM IMPACT

We have analysed the spatial and temporal distribution of HRV-based detections by counting their occurrences on a monthly base. Figure 3 illustrates its impact on derived full disc cloud observation frequency.

Continental cloud patterns

Two cases are used to illustrate a typical fair-weather convection over a desertic area in Africa (figure 6) and over Europe (figure 7).

FULL EARTH DISC AND SYNOP AND SHIP CLOUD COVERS COMPARISON

We have specifically studied the HRV-based cloud detections using SYNOP over land, provided by M.Lockhoff from CM-SAF, and SHIP over sea, for the available eleven-month period (from October 2009 to August 2010). As those detections are sensitive to fractional cloudiness, the method as described in Derrien and Legleau,2005 used to assess the cloud detection accuracy of cloudy and clear events, is no longer appropriate to estimate the improvement using HRV data, as it ignores the partial cloud covers.

We have categorized cloud cover in three classes to account for every cloudiness: clear when N≤2, broken when N=3,4,5, cloudy when N≥6. Satellite-based cloud cover is determined by counting any CMa cloud contaminated pixel in the 5x5 SEVIRI LRes box centred on the ground measurement position as 100% covered. Obviously this assumption is not always true particularly for situations when the HRV is likely to detect a cloud while the coarser SEVIRI missed it. Therefore we can note that the satellite-based cloud cover is by definition biased towards the cloudy class. Three-by-three confusion matrices are built, and the resulting classifications scheme are studied. The population is restricted to cases where at least one LRes pixel in 5x5 box is detected by HRV test and where satellite zenith angle is lower than 78 degrees. For this population the two satellite-based cloud covers (noted v2010,with HRV and v2009, without HRV) are determined, their confusion matrices are given in tables 1 and 2. We note that broken ground-based cloudiness is 50% of the population when HRV test triggers at least once inside the 5x5 box (16% when all SYNOP of the period are considered). We note also our satellite-based classification trend to overestimate ground-based cloudiness, but our method counts every cloud detection as 100% covered by clouds. The Cramer's coefficient, that is in fact a normalization of the chi-square test to the data set size and the confusion matrix size, presents the advantage that its range lies between 0 (no correlation) and 1 (perfect correlation). It is defined as :

Cloud Mask Cloudy

Cloud Mask Broken

Cloud

) , ( 2 c r NxMIN V χ =
Where N is the total number of events of the confusion matrix, r and c are the numbers of rows and columns of the confusion matrix and χ 2 is the chi-square test. The table 3 shows that it increases when the satellite-based cloudiness includes the HRV detections.

CONCLUSION

The SAFNWC/MSG v2010 cloud mask at low SEVIRI resolution makes use of HRV information to better account for partial cloudiness inside SEVIRI pixels. The complementary HRV detection algorithm, applied only on pixels not flagged as cloudy by the primary cloud mask, exploits spectral, textural and temporal HRV features. By comparison with SYNOP and SHIP ground-based cloud covers over the full earth disc we have assessed an improvement when considering a three-class categorization of cloudiness (clear, broken, cloudy), of course this improvement depends on the efficiency of the initial cloud mask used in the comparison. We have also visually inspected resulting cloud masks, few observed false alarms are cloud shadows over bright backgrounds and misses may still occur over static orographic loud clouds. The quality of SAFNWC/MSG cloud mask over the lower part of the earth disc depends on the presence the HRV data in the pair of images required by the method. Therefore SAFNWC/MSG users processing such regions may decide to use it or not when configuring their applications, moreover HRV effective use is reported at the pixel level in the cloud mask quality information. The fusion of high and low resolution imagery is a challenging task of future geostationary meteorological satellites, such as Meteosat Third Generation. In the context of cloud detection, as a feedback of this work, we think that cloud shadows and terrain orientation effects in visible imagery will require a particular attention with the advances in spatial resolution.

Figure 1 :

 1 Figure 1: Spectral response functions of the MSG (Met-9) SEVIRI bands for the 0.6μm band (red) , 0.8 μm (green) and HRV (black).

Figure 2 :

 2 Figure 2: Illustration of the HRV window positions over MSG earth disc on 14 January 2010.

  A SEVIRI pixel on land at low resolution scale in the 11x11 neighbourhood of a HRV change detection is restored as cloud if at least 5 pixels inside 11x11 neighbourhood are detected by HRV change detection tests and if: • MIN(RN hrv 3x3) cur > 10% and • MAX(R hrv 3x3) > MEAN 11x11 (Max(R hrv 3x3) of the pixels detected by HRV if nb > 5 ) and • SD(R hrv 3x3) > 1.5% or MAX(R hrv 3x3)-MIN(R hrv 3x3) > MEAN 11x11 (MAX(R hrv 3x3)-MIN(R hrv 3x3)) of the pixels detected by HRV if nb > 5 )

Figure 3 :

 3 Figure 3: HRV-based detection frequency for January 2010 (left) and July 2010 (right)

Figure 4 :Figure 5 :

 45 Figure 4: Shallow convection clouds over Northern Atlantic, 8 September 2009. 17h00UTC, (120x120 LRes size). Top left: 0.8 μm reflectance; top right: v2009 cloud type; bottom left: HRV reflectance; bottom right: v2010 cloud type

Figure 6 :

 6 Figure 6: Fair-weather cumulus over Egypt near Nile delta, 9 September 2009, 10h00UTC, (120x120 LRes size). Top left: 0.6 μm visible reflectance; top right: v2009 cloud type; bottom left: HRV reflectance; bottom right: v2010 cloud type.

Figure 7 :

 7 Figure 7: Fair-weather cumulus in subsident air over western Europe, 23 June 2009, 14h00UTC, (1270x770 LRes size). Top left: enhanced 0.6 μm visible reflectance; top right: EUMETSAT RGB (12-12-9); bottom left: v2009 cloud type; bottom right: v2010 cloud type.

  Cramer's coefficients for v2009 and v2010 cloud cover classifications in three classes (cloudy,broken,clear) stratified by background for the subset where HRV detection triggers at least once in the 5x5 satellite target and with satellite zenith angle <78°

Table 1 : Confusion matrices for cloud mask v2010 and (v2.009) cloud covers at daytime for every SYNOP in the MSG earth disc from 1 October 2009 to 31 August 2010 when any HRV detection occurred in the 5x5 satellite target and with satellite zenith angle <78°

 1

	Mask	Total
	Clear

Table 2 : Confusion matrices for cloud mask v2010 and (v2009) cloud covers at daytime for every SHIP in the MSG earth disc from 1 October 2009 to 31 August 2010 when any HRV detection occurred in the 5x5 satellite target and with satellite zenith angle <78°

 2


N Mabrouk

G Deffuant

Comparison of discrete and continuous time implementations of a simple spatially explicit individual-based model

We propose to compare discrete and continuous time implementations of a simple birth-death individual-based model (IBM) to analytically derived expression of the mean and variance of the population size. Our results show that the discrete time implementation underestimate the variance, especially when increasing the size of the timestep.

Introduction

Over the last two decades, the individual-based approach has been established as an interesting technique for modelling ecological systems [START_REF] Judson | The rise of the individual-based model in ecology[END_REF][6][10] [START_REF] Grimm | Individual-based Modeling and Ecology[END_REF][3] [START_REF] Hellweger | A bunch of tiny individuals -individual-based modeling for microbes[END_REF]. With the increase of desktop computing power, it is quiet feasable today to simulate a population by simulating its constitutive individuals as they move across a landscape, interact, reproduce and die [START_REF] Dieckmann | The geometry of ecological interactions: simplifying spatial complexity[END_REF]. Such simulations have shown to be a valuable aid in studying many complex ecological systems and in reproducing realistic spatial and temporal patterns observed in real systems [START_REF] Picioreanu | Multidimensional modelling of anaerobic granules[END_REF] [START_REF] Railsback | Analysis of habitat selection rules using an individual-based model[END_REF][15] [START_REF] Pitt | An individual-based model of canid populations: modelling terrotoriality and social structure[END_REF]. However and before inferring any conclusion from these simlulations, the modeler should ensure that the computer code producing these results faithfully translates the formal rules as defined in the individual-based model and that design decision relative to temporal and spatial discretization, events execution order and state updating mode are not introducing simulation artifacts. This step is of critical importance and yet it is one of the most difficult as most IBMs have no reference solution (i.e. analytical solution) against which the simulation results can be checked.

While the use of modern compilers and the adoption of good programming practices (i.e. unit testing) can be very effective to prevent both syntically and non-syntically programming faults, the effect of the design choices made by the modeler to put the individual-based rules into a computer program is rarely assessed during IBM developpment. In this work we are concerned with the artifacts that may result from the temporal representation in individual-based models. Time is represented in IBMs used either a discrete or a continuous models. Most IBMs adopt a discrete time model where the time is subdivided into constant timesteps. The processes are then executed at each timestep based on the state of the system at the previous timestep and the state of the system is updated synchronuously at the end of the timestep [?]. This mode of time representation is very common [START_REF] Kreft | Individual-based modelling of biofilms[END_REF][?] due to the simplicity of its implementation. In the continuous representation of time, the model time no longer proceeds in time steps but from event to event. Each event (i.e. a birth of a new individual or death of an existing one), has an exact execution time and the computer simply executes the next action in a queue that is ordered by time [START_REF] Grimm | Individual-based Modeling and Ecology[END_REF]. This way of proceeding is also called dynamic scheduling [START_REF] Grimm | Individual-based Modeling and Ecology[END_REF] or discrete event model. Contrary to the discrete time model where the states of all the individuals are updated sychronuously at the end of the time step, in the continuous time model the state of an individual is usually updated immidiatly after the event relative to that individual is executed. Hence, for the continuous time model the updating scheme is said asynchronuous.

Artifacts due to the choice of the time model has been addressed by several authors mainly from the perspective of comparing synchronous and asynchronous updating modes as synchronous updating mode is associated with the discrete time model and the asynchronuous updating scheme with the continous time model. The comparison of the synchronous and asynchronous modes has been performed for cellular automaton [START_REF] Schoenfisch | Sychronous and asychronous updating in cellular automata[END_REF][?] and individual-based models [START_REF] Caron-Lormier | Asynchronous and synchronous updating in individual-based models[END_REF]. These comparisons showed the time model mode can influence the simulation results. The synchronous update of the state of the individuals has been criticized for its non-realism as many natural systems typically exhibit asynchronous dynamics [START_REF] Guo | Multi-timescale event-scheduling in multiagent immune simulation models[END_REF][?][?]. Additionnaly the synchronous mode requires from the modeler to decide about the order of execution of the event during a timse step. Axtell [?] showed that this may introduce additional artifacts and suggested that the order of execution should be randomized.

In order to assess the effect of time representation on the simulation result of an IBM, we use a simple IBM in which the individuals are subject to birth and death events. The IBM is sufficiently simple to allow analytical developpement. We propose to compare the IBM simulation results to the analytically developped solution and analyze how time representation and individuals state updating scheme can affect the IBM simulation results. Birth-death IBMs are common in ecology as these two processes are encountered in many ecological systems and are involved in the dynamic of the population size and the spatial distribution of the individuals. In a birth-death IBM, an initial population of individuals evolve in time due to individuals births and deaths. Each individual has proper birth and death probabilities per unit of time. these probabilities can be constant or may vary depending the local environment of the individual. For simplicity we consider only the case of constant density independent birth and death rates.

The paper is organized in three parts: in the first part we provide a detailed description of the IBM and the underlying algorithms representing time as discrete and continuous. In the second part we present the results of the comparison to the analytical solution and the comparison the spatial patterns. The third part is dedicated to the discussion of these results and concluding remarks.

Methods

We propose to derive a simple IBM for a (virtual) population of individuals living in a two-dimentional domain and subject to stochastic constant-rate birth and death processes. Such IBMs are valuable in showing the degree of variability that we might expect to observe in practice [?]. In this work, however, we use this IBM to illustrate how the choice of the time model (e.g continuous time versus discrete time) may affect the individual-based simulation results. The simplicity of this IBM allows the derivation of analytical expressions for the mean and variance of the number of individuals N (t) at the instant t. We propose to compare the individual-based simulations results to the analytical results.

Individual-based model description

We consider a two-dimensional spatial domain containing a population of individuals. The individuals are considered as point particles characterized by their location x = (x 1 , x 2 ) in this plane discretized with a constant space step dx. We suppose that the individuals are subject to stochastic birth and death processes. An individual in location x produces a new individual in location x ′ with a probability B(x, x ′ ) and can be removed from the system with a probability D(x). We suppose that the birth and death rates of the individuals are constant and independant from their location (though we keep the general notation where B and D are functions of the location x). The expression of B(x, x ′ ) is given by:

B(x, x ′ ) = bK ||x -x ′ || w b ( 1
)
where b is the constant density-independent probability per unit of time that an individual gives birth to a new individual.

K ||x-x ′ || w b
is a birth kernel describing the dispersion of the newly formed individuals arround the location x of the parent individual. The death rate of the individuals is given by:

D(x) = d (2)
where d is a density-independent death probability per unit of time. Because of the stochasticity of the births and deaths of the individuals each individualbased simulation represent a realisation of a stochastic process. The number of individuals at an instant t is a stochastic variable which mean and variance can be derived analytically and are given by [?]:

mean :M [N (t)] = N 0 e (b-d)t (3) variance : V [N (t)] = 2N 0 bt (4)
Note that for the particular case where b = d, the average population size remains constant M [N (t)] = N 0 while the population size grows exponontially for b > d and goes to extinction for b < d. In this work we assume that b = d.

Continuous time implementation

In the continuous time model, the time evolution of the system is defined by calculating the waiting time to the next event to occur. The time scale is then progressed with calculating the waiting time which is a real value. Gillespie [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF] proposed a Monte Carlo procedure for implementing such time model. the procedure iterates over the following steps: One of the advantages of the continuous time formulation is that designing actions concerning the order of execution of the individuals and the updating mode of the state of the system are implicitly solved in this algorithm. The events are scheduled in a randomized order and the update of the states of the individuals is asychronuous.

Discrete time implementation

The discrete time model can be implemented in different ways all having in common that time is subdivided into discrete time steps (not necessary with a constant size, but for simplicity we assume that ∆t = constant). The different implementations differ in the order of execution of the different processes (birth, death and update of the individual states). Figure 1 shows 4 examples with different orders of execution of the processes. For each process we iterate over the individuals according to the following generic algorithm : One of the difficulties that raises the discrete time model is that the order of execution of the processes can have an impact on the simulated results. This point will be illustrated by comparing the 4 configurations presented in figure 1.

Simulation results

We compare the IBM simulated results using the continuous and the discrete time representation to the analytical solution provided by equations 3 and 4. We proceed by replicating each individual-based simultion a number of times and than calculate the average and variance of N (t). The number of replicates is determined experimentally by testing 50, 150 and 250 replicates. It is expected that, in the absence of artifacts the increase of the number of replicates should yield to better matching between the analytical and the simulated results.

Comparison of the continuous time model and the analytical solution

Comparison of the discrete time model and the analytical solution

Before comparing the discrete time IBM results to the analytical equations we assess wether the order of execution of the processes has an impact on the simulated results. Figure 3 shows the average number of individuals calculated over 250 runs of the discrete time IBM using the configuration (a) and (c) (see figure ??). The results clearly show that the order of execution of the processes can habe an important impact on the simulated results. When death process is executed first the IBM predicts the rapid extinction of the population. This We proceed to the comparison of the discrete time IBM simulation (configuration (c)) to the analytical expression of the mean and variance of the population size. The simulations are replicated 50, 150 and 250 times and the results are reported in figure 4. While the mean number of individuals convergence rapidely to the analytical solution, the variance calulated by the individual-based simulation seems to be slightly lower than the analytical solution. This result is confirmed by taking a larger time steps (∆t = 5) as can be viewed from figure 5. Hence it seems that a large timestep induces an underestimation of the variance but has no effect (for this simple IBM) on the mean population size.

Comparison of the spatial patterns simulated with continuous and discrete time implementations

We assessed the effect of the time model on the spatial patterns simulated using the individual-based model. We use the radial pair correlation to caracterise the spatial pattern. The spatial pattern in this simple birth death depend to a large extent on their properties of the birth kernel. If the size of the birth kernel is small in comparison to the size of the domain, than the newborn individuals As for the variance, the discrete time model tends to underestimate the aggregation of the individuals measured with the pair correlation function (figure 6). The increase of the time step reinforces this tendency confirming that the pattern is sensitive to the choice of the timestep in the discrete time IBM.

Discussion

While a great effort has been devoted during the last years towards the development of improved software libraries for the implementation of IBMs [11] [17][?] many researcher agree that there is an increasing need to develop a good practice that guide the IBM developer during all the phases of the model construction and implementation. Often IBM developer are faced with several design decisions relative to space and time representations (discrete versus continuous), events scheduling and individuals state updating processes (synchronuous versus asynchronuous). These decisons can in many situations induce artifacts that, due to the IBM complexity and the absence of a reference solution, are difficult to assess. resentation requires also from the modeler to specify the order in which the processes are executed. This step can introduce additional artifacts. We show that by executing death events systemitically before birth events the simulated population size vanishes to zero. Similar artifacts have been reported in [?].

While the discrete time model is easier to implement and in many situations can reduce simulation time, much care should be considered in the choice of the time step in particular for IBMs where different time scales are considered. In many situations the computationnal advantage of the discrete time model is lost if the time step is set considering the rapid processes.

1 . 6 . 7 . 8 . 9 .

 16789 Set time to t = 0 2. Calculate the sums r b = N (t) i=1 b i = bN (t) and r d = N (t)i=1 d i = dN (t).The total rate at which events occur (births or deaths) is given by r(t) = r b (t) + r d (t)3. Choose a waiting time τ for the next event to occur according to τ = -1 r lnλ, where 0 < λ ≤ 1 is a uniformly distributed random number 4. Choose a birth or death event with a probability r b /r and r d /r respectively 5. Choose an individual k with a probability b k /r b (if the event is a birth) or d k /r d if the event is a death, where b k = b and d k = d are respectively the birth and death rates of the individual k Perform the selected event on the individual k Update time according to t = t + τ Update the number of the individuals N (t) Continue from step 2 until t < t end

1 .

 1 Set time to t = 0 2. iterate over the individuals and execute process 1 3. iterate over the individuals and execute process 2 4. ..

Figure 1 :

 1 Figure 1: Discrete time model: different scheduling possibilities

Figure 2

 2 Figure 2 compares the continuous time IBM simulation results to the analytical solution for 50, 150 and 250 replicates of the IBM simulation. The results show that with the increase of the number of replicates the continuous time model coverges successufully to the analytical solution.

  (a) runs = 50 (b) runs = 50 (c) runs = 150 (d) runs = 150 (e) runs = 250 (f) runs = 250

Figure 2 :

 2 Figure 2: Continuous time simulation. Number of individuals and variance over 50, 150 and 250 of the individual-based model

Figure 3 :

 3 Figure 3: Number of individuals versus time simulated using two algorithms of the discrete time implementation. (a) death events are executed first, (b) birth events are executed first

  (a) runs = 50 (b) runs = 50 (c) runs = 150 (d) runs = 150 (e) runs = 250 (f) runs = 250

Figure 4 : 1 (b) dt = 5 Figure 5 :

 4155 Figure 4: Continuous time simulation. Number of individuals and variance over 50, 150 and 250 replicates

  Time representation is a reccurent design decision in individual-based models. Our work shows that the choice of a discrete time model can introduce artifacts, mainly by underestimating the variance. The choice of the discrete time rep-

Figure 6 :

 6 Figure 6: Pair correlation function

Acknowledgement

We thank Claude Lobry, Alain Rapaport and Fabien Campillo for their helpfull suggestions. This work is supported by the SYCOM program, a grant from the French National Agency of Research (DISCO, program ANR)