Jean-Marc Leyssale
email: leyssale@lcts.u-bordeaux1.fr

Jean-Pierre Da Costa

Christian Germain

Patrick Weisbecker

Antoine Villesuzanne

Gérard Vignoles

Nanostructure of pyrocarbons : from HRTEM imaging to 3D atomistic models

Keywords: Pyrolytic Carbons, Nanostructure, atomistic modeling

Introduction:

Atomistic reconstruction methods are nowadays extremely useful tools to build realistic models of materials and correlates their properties to their atomistic structure [START_REF] Dove | Reverse Monte Carlo methods[END_REF]. However most of them, like the well-known Reverse Monte Carlo method, consists in building a model able to reproduce some orientationally averaged structural quantities like the Pair Distribution Functions (PDF) of the material, which is usually not discriminatory enough to yield a unique solution for systems having pronounced anisotropic nanotextures. For instance, two drastically different models have been reconstructed from the same PDF of a nanoporous carbon [START_REF] Petkov | Local structure of nanoporous carbons[END_REF]. The first one is a perfectly isotropic pack of small Polycyclic Aromatic Hydrocarbons (PAH) [START_REF] Acharya | Simulation of nanoporous carbons: A chemically constrained structure[END_REF] while the second one is a perfect stack of identical infinite graphene sheets containing a few percent of non-hexagonal rings [START_REF] Smith | A simple model describes the PDF of a non-graphitizing carbon[END_REF].

High Resolution Transmission Electron Microscopy is the technique of choice to describe the nanotexture of aromatic carbons [START_REF] Oberlin | Pyrocarbons[END_REF] and we have recently introduced a new method, the Image Guided Atomistic Reconstruction (IGAR) method, based on a combination of image analysis and synthesis and molecular dynamics simulations, aiming at building realistic atomistic representations of dense nanotextured carbons [START_REF] Leyssale | An Image Guided Atomistic Reconstruction of pyrolytic carbons[END_REF]. In what follows we briefly recall the methodological details of this approach and how we use it to build models of a rough laminar pyrocarbon, as prepared and after a partial graphitization through heat treatment. The models are then discussed in terms of nanotexture (HRTEM images), local average atomic structure (PDFs) as well as in terms structural defects (sp3 hybridations, non-hexagonal rings. Finally, from these structural observations we end up this paper by proposing some new, high stability, defective graphene sheets.

Image Guided Atomistic Reconstruction of pyrocarbons:

The IGAR procedure starts from some selected lattice fringe HRTEM images of the materials like those of a rough laminar pyrocarbon as prepared and after partial graphitization shown respectively on Figures 1a and1d. From these images, first (mean, variance, skewness and kurtosis) and second (autocorrelations) order multi-resolution 2D statistical descriptors are collected. Then, these descriptors are extended to 3D according to the orthotropy of the material, in our case the statistical equivalence of every plane by rotation around the main stacking direction (vertical plane on Figure 1). From these 3D statistics, the 2D image synthesis method of Portilla and Simoncelli [START_REF] Portilla | A Parametric Texture Model Based on Joint Statistics of Complex Wavelet Coefficients[END_REF] is easily extended to 3D and 3D HRTEM-like image blocks of the AP and HT pyCs can be obtained (Figures 1b and1e respectively). These blocks are then filled with carbon atoms at the suited density and a simulated annealing simulation is launched to cool the system from the liquid state (8000 K) to room temperature. In this simulation, molecalur dynamics is used to trace the atomic motions under the combine effects of an interatomic potential U REBO , the second generation REBO potential [START_REF] Brenner | A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons[END_REF], and of an external potential U HRTEM defined by the 3D image block:

U HRTEM = k Im I(r i ) i=1 N " .
In this equation, i is an index running over the atoms, I(r i ) is the grey level, ranging from 0 (white) to 1 (black) at position r i and k Im = 2 eV, is an energy conversion factor. This potential has for effect to bring the atoms to sit preferentially on the light areas of the image block while UREBO ensures that the atoms make proper chemical bonds between each other. As discussed in ref [START_REF] Leyssale | An Image Guided Atomistic Reconstruction of pyrolytic carbons[END_REF], a non-uniform cooling sequence with the slowest quench rate around the melting point of carbon is used. The models obtained for the two pyrocarbons are shown on Figures 1c and1f. They are very similar on the "chemical" point of view with around 98 % of three-fold atoms, forming mainly hexagonal rings (around 89 %) with 5 to 6 % of pentagonal and heptagonal rings. Both models also show a structure made of hexagonal rings domains linked together by high defects density areas and screw disclocations. Actually, the main difference between the two model materials is on the nanotexture; the HT pyC model showing extended and flat carbon sheets while the AP one shows smaller and more curved sheets.

Figure 1 :

 1 Figure 1: Illustration of the IGAR reconstruction of an "as prepared" (AP) and a "heat treated" (HT) pyrocarbons. (a): Filtered experimental HRTEM image of the AP pyC; (b) Synthesized 3D HRTEM-like image of the AP pyC; (c) Reconstructed atomistic model of the AP pyC; (d,e,f): equivalent to (a,b,c) for the HT pyC. Blue: Bonds bet w e en atoms belonging only to hexagonal rings; Orange: Other bonds.

Figure 2 :

 2 Figure 2: HRTEM images of the AP (a) and HT (b) simulated from the atomistic pyC models.

Figure 2

 2 Figure2shows the HRTEM images simulated; using the NCEMMS software[START_REF] Kilaas | Interactive software for simulation of high resolution TEM images[END_REF], from the atomistic models of the AP (Fig.2a) and HT (Fig2b) pyCs. Comparing these images with the filtered experimental images of Figs1a and 1dconfirms the ability of the procedure to generate a given nanotexture (small and curved fringes for the AP pyC; long and flat for the HT pyC).

Figure 3 :

 3 Figure 3: Pair Distribution Functions (PDFs ) of an AP pyC. Blue : Neutron scattering; Red: AP pyC model.

Figure 3

 3 Figure 3 shows the reduced atomic pair distribution functions (PDFs) G(r), rG(r) and r 2 G(r) measured on an AP pyC powder using Neutron scattering and the corresponding functions computed from the AP model. As can be seen, the model perfectly reproduces the experimental G(r) although the two first peaks are a bit too narrow, probably because of too rigid C-C bonds in the potential model. Plots of rG(r) and r 2 G(r) highlight the differences between the

Figure 4 :

 4 Figure 4: Snapshot of a carbon sheet taken from the HT pyC model. Blue: purely hexagonal C-C bonds; Orange: other C-C bonds; Red: C5-C7 pairs and chains.

Figure 4

 4 Figure4shows a carbon sheet extracted from the HT pyC model. The intra-sheet structure is made of misoriented nanometric graphenic domains (blue) connected by defective grain boundaries (orange) having a high density of pentagonal (C 5 ) and heptagonal (C 7 ) rings. These non hexagonal rings are found in similar proportions (between 5 and 6 %) and form C 5 -C 7 pairs and pair lines (red) explaining the flat shape of the sheets despite this high amounts of defects. Finally, many inter-sheet connections exist in the materials through screw dislocations as can be seen on the top right side of Figure4. The structure of the AP pyC is very similar and the main differences are that graphenic domains are smaller, C 5 -C 7 pairs are more scattered in the structure (shorter lines) and that the density of screw dislocations is higher.

Figure 5 :

 5 Figure 5: A ne w, stable, defective carbon sheet. a: Snapshot taken from the HT pyC model sho wing (green rectangle) a defect made of t wo C7, one inverted C6 and t wo C5 rings (2C7iC62C5). b: Carbon sheet built by periodically repeating the green rectangle.

Finally, we show on Figure 5a a particularly interesting class of defect observed in our HT pyC model. It is made of two adjacent heptagons (2C 7 ) connected with an inverted hexagon (iC 6 ) which is surrounded by two pentagons (2C 5 ). This motif allows two build the crystalline sheet shown on Figure 5B. Preliminary DFT calculations on this structure shows that it is more stable than haeckelite [10], the most stable carbon sheet found so far after graphene.


Maxime Berar

Françoise Tilotta

Joan Alexis Glaunès

Yves Rozenholc

Francoise M Tilotta

Joan A Glaun

Craniofacial reconstruction as a prediction

Keywords: Facial Reconstruction, Anatomical Landmarks, Surface Registration, Geodesics, Latent Root Regression

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In forensic medicine, craniofacial reconstruction refers to any process that aims to approximate the morphology of the face from the shape of the skull [START_REF] Wilkinson | Computerized forensic facial reconstruction: A review of current systems[END_REF]. It is usually considered when confronted with an unrecognisable corpse and when no other identification evidence is available. This reconstruction may hopefully provide a route to a positive identification. In recent years, computer-assisted techniques have been developed following the evolution in medical imaging and computer science. As presented in the surveys in [START_REF] Buzug | Special issue on computer assisted craniofacial reconstruction and modelling: Editorial[END_REF][START_REF] Clement | Computer-Graphic CranioFacial Reconstruction[END_REF][START_REF] De Greef | Three-dimensional cranio-facial reconstruction in forensic identification: latest progress and new tendencies in the 21st century[END_REF], new approaches are now available with reduced performance timeline and operator subjectivity.

Reproducing the stages of manual methods, the first machine-aided techniques fitted a skin surface mask to a set of interactively-placed virtual dowels which were on the digitised model of the remains [START_REF] Evenhouse | Computer-aided forensic facial reconstruction[END_REF][START_REF] Vanezis | Facial reconstruction using 3-D computer graphics[END_REF][START_REF] Shahrom | Techniques in facial identification: computer-aided facial reconstruction using a laser scanner and video superimposition[END_REF]. Latter techniques have moved away from the manual techniques and two kind of methods can be distinguished based on the representation of bone and soft-tissue volumes. The first techniques aim to conserve the continuous nature of the skull and soft-tissue surfaces. Estimates of the face are obtained by applying space-deformations to couples of identified bone and soft tissue surfaces, called reference surfaces. These deformations are learned between the surface of the dry skull and the surfaces of the reference skull. They can be parametric (e-g B-splines) [START_REF] Kermi | A non-linear registration method guided by b-splines freeform deformations for three-dimensional facial reconstruction[END_REF][START_REF] Vandermeulen | Computerized craniofacial reconstruction using CT-derived implicit surface representations[END_REF], implicit using variational methods [START_REF] Mang | A multi-modality computer-aided framework towards postmortem identification[END_REF][START_REF] Mang | Soft-tissue segmentation in forensic applications[END_REF] or volumetric [START_REF] Quatrehomme | A fully three-dimensional method for facial reconstruction based on deformable models[END_REF][START_REF] Nelson | The application of volume deformation to three dimensional facial reconstruction: a comparison with previous technique[END_REF]. Depending on the method, the final estimated face can be either the deformed face whose reference skull is the nearest to the dry skull [START_REF] Quatrehomme | A fully three-dimensional method for facial reconstruction based on deformable models[END_REF][START_REF] Nelson | The application of volume deformation to three dimensional facial reconstruction: a comparison with previous technique[END_REF] or a combination of all the deformed soft-tissue surfaces [START_REF] Vandermeulen | Computerized craniofacial reconstruction using CT-derived implicit surface representations[END_REF][START_REF] Tu | Automatic face recognition from skeletal remains[END_REF].

The second type of approach chooses to represent individuals using a common set of points. As the position of the corresponding points for all the individual can be summarized as variables in a table, the main idea is then to use statistics to decipher the relation between the skull and the soft-tissue. The common set of points can either be anatomical landmarks [START_REF] Vanezis | Facial reconstruction using 3-D computer graphics[END_REF][START_REF] Claes | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: methodology and validation[END_REF] or semi-landmarks located following a point correspondence procedure [START_REF] Berar | 3D semi-landmarks-based statistical face reconstruction[END_REF][START_REF] Kähler | Reanimating the dead: Reconstruction of expressive faces from skull data[END_REF]. Semi-landmarks are defined as points that do not have anatomical meanings but that match across all the samples of a data set under a reasonable model of deformation [START_REF] Bookstein | Landmark methods for forms without landmarks: morphometrics of group difference in outline shape[END_REF]. The larger the set of points, the closer this surface representation is to a real surface.

Apart from the practical constraint of the number of anatomical landmarks that an expert can define and extract, there is no justification to a chosen number of points. Indeed, the information given by the position of a skull anatomical landmarks is double. First, there is the geometric information given by the coordinates of the point. Then, the "anatomic" information is provided by the measuring of tissue thickness made at this point. This information is available for a limited number of points, due to the difficulties in localisation and in thickness computation. However, the geometric information given by the position of the point can be completed by automatic methods of semi-landmarks extraction. The second part of the data analysis framework consists in learning relationships between the soft-tissue variables and the bone variables. Facial reconstruction consists of predicting the positions of the soft-tissue points knowing the positions of the set of skull landmarks. In this regression problem, the skull points are entries of a regression model while the face points are the model output. Typically, a linear modelisation is used to fit this regression based on Principal Components Analysis (PCA) for example [START_REF] Tu | Automatic face recognition from skeletal remains[END_REF][START_REF] Claes | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: methodology and validation[END_REF][START_REF] Berar | 3D semi-landmarks-based statistical face reconstruction[END_REF], following the work made in the statistical atlas field. However, several linear regression methods have been developed, some sharing the use of PCA. For example latent variables regression methods such as Latent Root regression [START_REF] Vigneau | A new algorithm for latent root regression analysis[END_REF][START_REF] Gunst | A comparison of least squares and latent root regression estimator[END_REF] are designed to take the presence of variable colinearity into account, namely the positions of the skull landmarks and of the face semi-landmarks.

In this paper, we propose a facial reconstruction technique using Latent Root Regression, before comparing the results given to those obtained by a PCA model. In conjunction, we question the number of skull landmarks necessary. Basing our first experimentation on skull anatomical landmarks extracted by an expert, we will iteratively add supplementary mathematical skull semi-landmarks following the point correspondence technique described in [START_REF] Wang | Shape-Based 3D Surface Correspondence Using Geodesics and Local Geometry[END_REF], which relies on the geodesic paths between the landmarks to define new landmarks.

The paper is organised as follows. The material and method are presented in Section 2: Section 2.1 presents the material on which this study has been carried out. Section 2.2 and Section 2.3 focus on resolving the point correspondence problem, describing the two methods used to obtain the two subject-shared description of the bone and soft-tissue surfaces. Section 3 presents the statistical methods used: the building and use of a statistical shape model and the multivariate Latent Root Regression method. Section 4 shows the results obtained by the different models and discusses the influence of the number of skull landmarks and of the statistical method chosen.

Material and method

2.1. Material. This study was performed using whole head and skull surface meshes extracted from whole head CT scanners for a project on facial reconstruction at Paris Descartes University. The head CT-scan database of healthy people build for this project is composed of several type of data: head CT-scans, triangulated and closed surfaces covering the skull and the face, soft-tissue measurements at predefined skull landmarks [START_REF] Tilotta | Statistical Facial reconstruction by tree functional regression on surfaces[END_REF][START_REF] Tilotta | Construction and analysis of a head ct-scan database for craniofacial reconstruction[END_REF]. These anatomical skull landmarks were manually located on each CT Scan according to classical methods of physical anthropology (13 midpoints and two sets of 13 lateral points). This database, and the processes performed on it, are described in detail in [START_REF] Tilotta | Construction and analysis of a head ct-scan database for craniofacial reconstruction[END_REF] In the framework of this study, we focus on a group of 47 European female patients aged from 20 to 40 years. As properties such as age, gender and ancestry can be determined by anthropological examination [START_REF] Reichs | Forensic anthropology in the 1990s[END_REF], we choose to build statistical models on a homogeneous database following these criteria. In our case, the 85 subjects of the database were distributed according to sex and age group (20-40 years, 40-65 years) with age groups determined so as to take ageing into account. Women aged from 20 to 40 years correspond to the largest group of the database. Body weight is another factor that affects facial form [START_REF] Starbuck | The affect of tissue depth variation on craniofacial reconstructions[END_REF]. In our case, none of the 47 women has a Body Mass Index (BMI) superior to 30. 8 subjects are underweight (BMI < 19), 34 subjects correspond to a normal weight (BMI: 19-24) and 5 subjects are overweight (BMI: 14-30).

The entries of our database will consist of left or right halves of each surface. The skull and the face do not have symmetric shapes, but the relationship between these shapes do not depend on the side of the head. The plan minimising the distances to the anatomical midpoints has been chosen as an artificial boundary between the right and left part of the shapes. The next step is to extract points which correspond to the same places on the different individuals, with respect to this symmetry constraint. 2.2. Point correspondence procedure for the bone surface. The anatomical landmarks located by the expert (Figure 1,A) establish a first correspondence between the skulls. Following the scheme presented in [START_REF] Wang | Shape-Based 3D Surface Correspondence Using Geodesics and Local Geometry[END_REF], we define a set of triangular connections between these anatomical landmarks. For each pair of connected points, we can extract a geodesic between these points. Geodesics are defined to be the shortest path between points on the curved spaces of the shape surfaces (see Figure 1,B). As the shape surface between two landmarks is different from a sphere, theses geodesics are unique with high probability. At this step, a gross template of curbs on the surface between the landmarks is built. We then can define new landmarks as the midpoints of each geodesic. These landmarks are used to build new geodesics as seen in Figure 1,C and a more dense triangulation is then derived. The process is iterated to form a dense geodesic triangulation with associated semilandmarks. The obtained structures form meshes, which share the same structure for each individual, and implicitly solve the point correspondence problem. Moreover, the defined structure is symmetric: the two entries (left and right) of the database share a common substructure and set of midpoints (Figure 1,D). Due to numerical instabilities, two methods of geodesics computation on surface meshes have been used: The Surashky algorithm [START_REF] Surazhsky | Fast exact and approximate geodesics on meshes[END_REF] and the Fast Marching Algorithm algorithm [START_REF] Sethian | Level sets methods and fast marching methods[END_REF] implemented by Peyre in the Geowave library. For two iterations of this procedure, the results show three sets of skull landmarks for each individual: the 26 original landmarks (13 midpoints), a second set composed of 55 points added by the first iteration of the procedure (10 midpoints) then completed with 164 points (24 midpoints). For more iteration, the point correspondence step is confronted with degenerated triangles -i.e. segments or points -in the zones where the anatomical landmarks are close to each other.

2.3.

Point correspondence procedure for the soft-tissue surface. For the soft-tissue surface, no landmarks are located. Moreover measures of tissue thickness are not provided. Instead of facial points analogous to the anatomical skull points, we extract for each individual a set of semi-landmarks which are neither really dense nor sparse. Working on the "half" surfaces previously defined, the point correspondence procedure registers a reference mesh (see Figure 2,A) on the individual soft-tissue surface (see Figure 2,B) resulting in a deformed reference mesh (see Figure 2,C). The mesh used as the reference mesh corresponds to the face region of the head mesh modelled by F. Pighin [START_REF] Pighin | Resynthesizing facial animation through 3D model-based tracking[END_REF], where the density of vertices is greater in zones with high bending than in zones with low bending. The registration is made computing an elastic deformation between the reference mesh and soft-tissue surfaces of the database. The deformed meshes of each entry of the database have the same number of vertices (1741 for the mesh of a half face). The assumption of semi landmarks is then made: each vertex of the deformed reference meshes matches the same point for every individual. The 3D to 3D meshes matching algorithm used is a modified version of Szeliski algorithm [START_REF] Szeliski | Matching 3-D anatomical surfaces with non-rigid deformation using octree splines[END_REF].

This dissimilarity between the soft-tissue surface meshes and the reference meshes have consequences. The distances from the vertices of the deformed reference mesh to its associated soft-tissue surface mesh are null. But, the distances from the vertices of the soft-tissue mesh surface to the deformed reference mesh are not null, as seen in Figure 3. The highest distances (larger than 3 mm) correspond to parts of the soft-tissue surfaces which do not have corresponding regions in the reference surface. Other distances correspond to regions like the forehead or the cheeks where the registration suffers from the low density of the semi-landmarks reference mesh. Vertices with no direct counterparts can be as far as 2 mm from the surface defined by the deformed reference mesh. Measures of the error introduced during this point correspondence are influenced by the large distances generated by the lack of correspondence on the boundaries. The individual mean of the distances ranges from 0.54 mm to 2.66 mm. However, the median which correspond to the middle sample of the ordered samples, is by definition less influenced by the presence of a small number of high value outliers. Individual distance medians range from 0.17 mm to 0.34 mm. The mean median of distances is 0.22 mm (with standard deviation of 0.04 mm).

Statistical Methods

xi = [S x 1 S y 1 S z 1 • • • S x N S y N S y N ], (1) ỹi = [F x 1 F y 1 F z 1 • • • F x L F y L F z L ]. (2)
Two geometricaly averaged templates x and ȳ are computed and the data centred :

x i = xi -x, (3)
y i = ỹi -ȳ. (4)
The datatables X and Y, of size n × 3N and n × 3L, encompass the variables corresponding to the n centred samples x i and y i of the learning database. Consider the matrix Z = [XY] formed by merging data tables X and Y and perform Principal Component Analysis on Z. The result of this PCA is a correlation-ranked set of statistically independent modes of principal variations a j . These principal modes or components are in turn vectors of 3D coordinates (of size 3(N + L)) defined as linear combinations of points position. These vectors capture the variations observed over all subjects in the database and are the eigenvectors of the matrix t ZZ associated to the eigenvalues λ i sorted such as

λ 1 ≥ • • • ≥ λ i ≥ • • • ≥ λ N +L ≥ 0. λ i a i = t ZZa i . ( 5
)
Every entry z i in the database can now be represented as a weighted linear combination of these principal components :

z i = m j=1 b ij a j , (6)
where b ij is the weight attached to sample i and component j. Each eigenvector a i obtained with PCA for which λ i > 0 can be decomposed as the juxtaposition of two vectors a i = [v i w i ], with v i ∈ ℜ 3N and w i ∈ ℜ 3L . Each part x i and y i associated to the entry z i can be represented in the model sharing the same weights b ij and the 3D coordinates v i and w i of the vectors of the principal components :

x i = xi -x = m j=1 b ij v j , (7)
y i = ỹi -ȳ = m j=1 b ij w j . (8
0j ||x 0j-1 -b 0j v j || 2 , (9)
x 0j = x 0j-1 -b 0j v j , (10)
with x 00 = x 0 .

The solution is given by the projection of the residual centred skull x 0j landmarks upon the vector

v j : b0j = < x 0j , v j > ||v j || 2 . ( 11
)
and the facial reconstruction obtained is:

ỹi = ȳ + m j=1 b0j w j . (12)
3.3. Latent Root Regression. Latent Root Regression (LRR) is also a linear regression method. The multiresponse linear regression model for centred data is defined as :

Y = XB + E, ( 13
)
where B is a 3N × 3L matrix of regression coefficients and E is a noise matrix of size n × 3L. The elements of the matrix E are assumed to be centred with EE = 0 with covariance matrix Σ. Given a new set of centred skull landmarks x 0 , an estimate of y 0 will be:

ŷ0 = x 0 B. ( 14
)
where B is an estimate of B.

LRR is similar to Partial Least Square (PLS) regression or Principal Component Regression (PCR). It originates from chemiometrics were a small number of predictors must predict a great number of output. It is then particulary adapted to the ratio between the 245 skull landmarks and the 1741 face points. Multiresponse Latent Root Regression gives an estimate of the regression matrix coefficients using Principal Components Analysis. Single response Latent Root Regression [START_REF] Gunst | A comparison of least squares and latent root regression estimator[END_REF] uses the same vectors v i as the statistical shape model to estimate B. As these vectors are not necessarily orthogonal, an iterative procedure build upon the first latent variable is needed (see Vigneau and Qannari [START_REF] Vigneau | A new algorithm for latent root regression analysis[END_REF] for details). It results in a sequence of r orthogonal vectors ṽi which enables us to compute regression coefficients B, following the formula :

B = r i=1 ṽi t ṽi t ṽi t XXṽ i t XY. ( 15
)
The prediction of the vector ŷ0 is then given by

t ŷ0 = t x 0 r i=1 ṽi t ṽi t ṽi t XXṽ i t XY. (16)

Results

4.1.

Validation. The validation of the proposed statistical methods for craniofacial reconstruction is obtained by a leave-one-out cross-validation procedure. In turn, a sample is removed from the database and used as a test case. The remaining samples are used to create the statistical model. The reconstruction procedure is applied to the test sample and a reconstruction error computed. Each sample is composed of two couples, left and right, of soft-tissue and bone surfaces while skull points of each couple are used as separate entries. The resulting predicted surfaces for the half of the face can then be compared with the skin surface of the test case representing the home truth.

A quantitative reconstruction error is performed by calculating the distances between every point on the reconstructed skin surface and the measured skin surface of the test case. Other criteria for assessing the resemblance between the reconstructed and the original face surfaces can be chosen. On the one hand, exterior experts can judge the resemblance between the two shapes offering a subjective evaluation of the results [START_REF] Claes | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: methodology and validation[END_REF][START_REF] Quatrehomme | Assessment of the accuracy of three-dimensional manual craniofacial reconstruction: a series of 25 controlled cases[END_REF]. However, this kind of evaluation needs the reconstruction process to be more complete than the prediction of the shape of surfaces, for example by including a virtual make-up step. On the other hand a selection of specific measurements can be made between landmarks of the bone and tissue surfaces (for examples of such anthropological distances see [START_REF] Quatrehomme | Assessment of the accuracy of three-dimensional manual craniofacial reconstruction: a series of 25 controlled cases[END_REF]). As our methodology chooses to use semilandmarks instead of landmarks located by an expert, such a type of measurement will in our case be very dependent from the point correspondence step. However, the local repartition of the distance error can be computed for each semi-landmark providing distance maps which give the prediction error for each semi-landmark. Another point-to-surface error can be computed: the distance between the points of the skin surface of the test case and the reconstructed skin surface. As it has been observed during the evaluation of the mesh-matching procedure, the dissimilarity in densities of the meshes overestimates the error and the median of the distances is a good way to take these problems into account. Subsequently, reconstruction errors will denote the mean distance between every point on the reconstructed skin surface and the measured skin surface of the test case.

Results

. The leave-one-out cross validation procedure described previously was performed on the currently available database resulting in 47 test cases. As the database is composed of half parts of the bone and skin surface, we can use as much as 91 components for the prediction of the soft-tissue surface. A first degree of freedom is taken away by the data average. The other limiting factor of the number of components is the number of skull points. The numbers of coordinates of the known points are 65 (N 0 ), 220 (N 1 ) and 688 (N 2 ), which in the first case is less than the size of the learning base. The reconstruction errors for the first 64 components are presented in Figure 4 for the varying number of points and following the statistical methods presented in sections 3.2 and 3.3. For the first method, the reconstruction errors decrease greatly for the first parameters until a minimal value is achieved, then the reconstruction errors increase lightly. The minima of the mean reconstruction errors are 2.07 ± 0.78 mm (N 0 ), 1.95 ± 0.73 mm (N 1 ), 1.90 ± 0.72 mm (N 2 )) and are achieved for 10 components (N 0 ,N 1 ) or 12 components. For LRR models, reconstruction errors follow another behaviour: reconstruction errors increase greatly after the optimal number of components. The minima of the mean reconstruction errors are 1.69±0.40 mm achieved for 7 components (N 0 ), 1.44±0.36 mm achieved for 14 components (N 1 ) and 1.39 ± 0.26 mm achieved for 15 components (N 2 ). For all models, the increase in the number of skull points corresponds to a decrease in the reconstruction error. It also results in the increase in the number of selected components. The numbers of skull semi-landmarks are defined by the initial configuration of the mesh linking the anatomical landmarks. Improvements to the point correspondence procedure of the skull semi-landmarks will enable more points to be obtained as well as trying to find the optimal number of skull semi-landmarks. The minima of the reconstruction error for the LRR models are a half millimetre below the minima obtained for the PCA models. The selection of the number of components for the model has not the same consequences for the two methodologies. The reconstruction error for the maximum number of components in PCA models is near the minimal error whereas for LRR models the validation of the number of components must be carried out. 6. The resulting spatial maps of the quality of the reconstruction procedure can be seen in Figure 7. The distance card from the predicted points toward the original surface has a mean error of 1.08 mm (left). For this case, the highest reconstruction errors are located on the superior part of the forehead and on the lips. This part of the forehead corresponds to a zone with no skull landmarks. The position of the lips is influenced by the mouth overture, where variations occur during the acquisition. Similar remarks can be made concerning the more localised error in the regions of the nose and of the eyes.

For the optimal component number, we calculate mean and standard deviation for each predicted point of the mask. The mean of the point errors correspond to the reconstruction error (1.39 mm) and the standard deviation by points is then 1.02 mm. The resulting spatial maps of the quality of the reconstruction procedure for the optimal number of parameters can be seen in Figure 8 for the local mean and Figure 9 for the local standard deviation. The facial area with the highest reconstruction error (3 mm) corresponds to the boundary between the neck and the chin and is an artifact of the point correspondence step. There is no explicit correspondence to fix the limits of the mask in these zone. Others regions where the prediction error is high (2 mm to 3 mm) are the masseter region and the tip of the nose. These regions have few skull-landmarks and the bones do not support the softtissue for a large part of the cheeks. The regions with the smallest errors (0.5 mm to 1 mm) are concentrated towards the middle of the face, a part where the number of skull landmarks is high and where the inter-subject correspondence between the face meshes is more constrained. Figure 10 shows the distance map obtained for the best PCA model. The repartition of the error is similar with a greater amplitude.

Discussion

We compare our results to those of [START_REF] Claes | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: methodology and validation[END_REF] and [START_REF] Vandermeulen | Computerized craniofacial reconstruction using CT-derived implicit surface representations[END_REF]. Among reconstruction techniques, the technique described in [START_REF] Claes | Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: methodology and validation[END_REF] is close to the PCA model reconstruction techniques described here, with a supplementary deformation phase after the statistical prediction. The study is conducted on a database of 118 samples and use 57 landmarks for prediction. The number of face points obtained is likely to be high as the face meshes presented seem to tend towards surfaces. The reconstruction error corresponds roughly to our point-to-surface reconstruction errors. The mean reconstruction error observed is 1.14 mm with a standard deviation of 1.04 mm. The highest reconstruction errors (4 mm) are located in the chin and eye regions. Then, the region of the cheeks and the nose contains reconstruction errors of up to 2 mm, with a notable exception of the tip of the nose which is predicted separately. Regarding the smaller database and difference in the points correspondence step, our method achieves similar results with a generally simpler methodology. Concerning the distribution of the error on the face, we share similar zones in the cheeks and the nose whereas the chin and eye regions show smaller errors. The quality of prediction in the latter regions is more susceptible to be influenced by the point correspondence step of the soft-tissue surfaces.

The technique developed in [START_REF] Vandermeulen | Computerized craniofacial reconstruction using CT-derived implicit surface representations[END_REF] is based on the use of continuous surface and the study conducted on 20 samples. The mean reconstruction error is 1.9 mm with a standard deviation of 1.7 mm. The largest reconstruction errors (2-3 mm on average) occurs on the nostrils and masseter region. We surpass those results, however based on a smaller database. We can add that the regions with large reconstruction errors coincide once more.

A local method of facial reconstruction proposed by [START_REF] Tilotta | Statistical Facial reconstruction by tree functional regression on surfaces[END_REF] combines prediction obtained on surface patches, delimited by landmarks. The study has been performed on two regions: the nose region and the chin region. Local mean reconstruction errors of 1.40 ± 0.25 mm (nose) and 1.51 ± 0.67 mm (chin) occur for our method. The results presented in this report outperform these estimations with a mean reconstruction error of 0.99 mm, which motivates us to consider more local procedure in the reconstruction process.

Conclusion

We proposed a statistical method for 3D computerised forensic facial reconstruction, based upon Latent Root Regression. In conjunction, we have evaluated the influence of the number skull landmarks used. This method was based on an homogeneous database described in [START_REF] Tilotta | Construction and analysis of a head ct-scan database for craniofacial reconstruction[END_REF]. As with most techniques based on a data analysis framework, the presented method relies on the use of a common set of points for the description of the individuals. In this set of points, anatomical skull landmarks are completed iteratively by points located upon geodesic curbs linking the anatomical landmarks. It enables us to artificially increase the number of skull landmarks. Facial landmarks are obtained using a mesh-matching algorithm between a common reference mesh and the individual soft-tissue surface meshes.

The regression problem of predicting the positions of the semi-landmarks of the soft-tissue surface knowing the positions of the skull landmarks has been resolved using two methodologies: the building of a linear regression model, based on PCA and the use of Latent Root Regression. The accuracy of the reconstructions made by the methods was measured by leave-one-out cross-validation. First these results show the impact of an increased number of skull landmarks for the accuracy of the error computation of the reconstruction for both statistical methods. On this same set of couples of points, the LRR model shows better results whatever the number of components and the number of points. These results were discussed in regard to other facial reconstruction methods and performs comparatively to these other methods.

Some extensions can be proposed to the reconstruction method. First of all, having a larger database will increase the flexibility of the model. The more surface examples the model has, the better it learns the relationship between the two surfaces. Secondly, improvements of the correspondence steps will greatly benefit the reconstruction process. Improvement in the geodesics point correspondence method will help define an optimal number of skull landmarks. Improvements in the point correspondence step for the face will reduce the disparity observed in the computation of the distances between the reconstructed and the reference surface. Lastly an automatic extraction of the anatomical landmarks from the skull would make the complete reconstruction pipeline automatic.

Figure 1 .

 1 Figure 1. Extracted landmarks and semi-landmarks for the skull ([22])

Figure 2 .Figure 3 .

 23 Figure 2. Mesh-matching points correspondence procedure : (A) reference mesh (B) Soft-tissue surface (C) Deformed reference mesh

3. 1 .

 1 Statistical Frameworks. The variables xi are obtained from the positions of the N skull points. The variables ỹi are obtained from the positions of the L soft-tissue points :

) 3 . 2 .

 32 Best Model Fit. For facial reconstruction, we seek the best model fit of the PCA-model, i.e. the weights b 0j of the model the nearer from the centred measured skull landmarks x 0 = x0x, which are the only points known. The problem is resolved finding successively each weight b 0j and removing its contribution to the skull landmarks : b0j = arg min b

Figure 4 .

 4 Figure 4. Mean reconstruction errors

Figure 5 .

 5 Figure 5. histogram of the individual reconstruction errors for N2 points

Figure 6 .Figure 7 .

 67 Figure 6. Example of reconstructed complete face (left) and reference (right)

Figure 8 .

 8 Figure 8. Mean of the local error. The color plot shows the distribution over the surface from 0.0 mm (blue) to 3.0 mm (red).

Figure 9 .

 9 Figure 9. Standard deviation of the local error.The color plot shows the distribution over the surface from 0.0 mm (blue) to 3.0 mm (red).

Figure 10 .

 10 Figure 10. Mean of the local error for PCA model. The color plot shows the distribution over the surface from 0.0 mm (blue) to 3.0 mm (red).


Jean-Marc Leyssale

Jean-Pierre Da Costa

Christian Germain

Patrick Weisbecker

Gérard L Vignoles

An atomistic reconstruction of the nanostructure of pyrolitic carbons guided by HRTEM data

A new method for the generation of atomistic models of dense nanotextured carbons is presented. This method is based on the statistical analysis of HRTEM images and their threedimensional extension through image synthesis under constraint. The resulting 3D images then serve as an external potential bringing the atoms to settle preferentially on the white areas during a conventional simulated annealing simulation. Application of this method to the case of two laminar pyrocarbons, differing in their degree of disorder, highlights the promising nature of this approach.

INTRODUCTION

Atomistic reconstruction methods are nowadays well-established tools for linking experimental characterization data to the atomic scale structure of matter [1]. Most of them are based on the reproduction of orientation-averaged structural features like the pair distribution functions (PDFs) so that these methods are usually very efficient for isotropic systems. However, when dealing with materials displaying a neatly anisotropic nanotexture, like most dense graphene-based carbons, the success of such approaches relies strongly on the initial guess of the atomic structure [2]. In crude terms, this means the failure of these methods for that class of materials. For instance, two PDF-based computer reconstructions of the same material [3] show a drastically different nanotexture: one, quite isotropic, is an amorphous carbon with a dominating sp2 character [START_REF] Acharya | MRS Fall Meeting Symposium NN proceedings[END_REF] while the other one is a strongly anisotropic stack of faulted graphene sheets [5].

Nowadays, the nanotexture of turbostratic carbons, like for instance the pyrocarbon matrices present in many industrial C/C composites [6], is usually described from High Resolution Transmission Electron Microscopy (HRTEM) lattice fringe images [7]. Although the information provided by this technique is mainly qualitative, recent developments in image analysis methods allow a finer description of the nanotexture to be drawn [8][9][10][START_REF] Costa | An image analysis approach for the structural characterization of pyrocarbons[END_REF].

We propose here a new atomistic reconstruction strategy based on the information provided by such HRTEM images. The two steps of this method, i.e. (i) synthesis of 3D "HRTEM-like" images from 2D experimental HRTEM images under an orthotropy assumption and (ii) molecular dynamics annealing simulations using the synthetic 3D images as a complementary energy term to the interatomic potential, are described in this what follows. This new reconstruction strategy is then applied to two HRTEM images of a rough laminar pyrocarbon. The first image is typical of an as prepared material whereas the second one concerns a heat treated, partially graphitized material.

THEORY Image analysis/synthesis

To synthesize 3D HRTEM-like images, we use an original pyramidal analysis/synthesis scheme derived from a standard 2D approach [START_REF] Portilla | [END_REF]. The novel procedure breaks down into three steps:

(i)
2D analysis: decomposition of the HRTEM image sample into a set of multi-resolution subbands using the steerable pyramid decomposition [START_REF] Portilla | [END_REF]; statistical analysis of the subbands; production of a pyramidal collection of first order (mean, variance, skewness, kurtosis) and second order (autocorrelation coefficients) 2D statistics. All the analysis is performed on experimental HRTEM images after application of a radial and directional band-pass filter to remove low-frequency (background gradients) and high-frequency (noise) artifacts. This filtering also has the benefit of making the images circular (or periodic), a suitable property for image synthesis and atomistic reconstruction methods.

(ii)

2D to 3D statistical inference: production of a pyramidal collection of 3D target statistics.

(iii) 3D HRTEM-like image synthesis: a 3D data block, initially random, is alternately decomposed using a 3D pyramidal decomposition, modified to meet the 3D target statistics, reconstructed, re-decomposed, etc.

The main difficulty of this new approach is to obtain 3D target statistics from the analysis of 2D experimental images. This is trivial for first order statistics (mean, variance, skewness and kurtosis), the main difficulty is in determining target 3D autocorrelation coefficients r 3D (x,y,z) from the reference 2D autocorrelation coefficients r 2D (u,v). Such statistics describe the dependence relationship between two pixels in a texture (respectively two voxels in a solid texture). In the case of dense pre-graphitic carbons, texture anisotropy comes from the property of orthotropy. The latter can be stated as follows: all 2D views containing the direction of orthotropy, the stacking direction of graphene sheets, have the same statistical properties. For instance the direction of orthotropy on the two HRTEM images of Let (O, x, z) be the 2D Cartesian coordinate system of the image sample. The labels x and z refer respectively to the horizontal and vertical directions. Let y be the direction orthogonal to x and z so that the system (O, x, y, z) is orthonormal. 3D autocorrelation coefficients are related to 2D coefficients through polar/spherical coordinates: In other words, r 3D (x,y,z), can thus be deduced from r 2D (u,v) by interpolation.

Once inferred from 2D, the target statistics have to be imposed on a 3D random block through an ad hoc multi-resolution framework. This framework must be based on a pyramidal decomposition adapted to the anisotropic nature of the 3D textures addressed. Each 3D subband must relate to one specific 2D subband so that the imposition of the autocorrelation coefficients of the 2D subband on the 3D subband makes sense. Let B 3D (ω) be any of the frequency filters used in the 3D pyramidal decomposition to extract a specific subband. Let ω denote a 3D frequency, (ω z ,ω y ,ω z ) its expression in Cartesian coordinates and (ω ρ ,ω θ ,ω φ ) its expression in spherical coordinates. We have just seen above that the 3D autocorrelation coefficients do not depend on θ. For the same reason, the 3D expressions of the frequency filters used for the 3D multiscale decomposition must not depend on ω θ , i.e.:

B 3D (ω) = B 3D (ω x ,ω y ,ω z ) = B 3D sph (ω ρ ,ω θ ,ω ϕ ) = B 3D sph (ω ρ ,0,ω ϕ ) = B 2D pol (ω ρ ,ω ϕ ) = B 2D (ω ρ cosω ϕ ,ω ρ sinω ϕ ). ( 2
)
In other words, the 3D frequency filters can be deduced from the 2D frequency filters by interpolation.

This analysis/synthesis method is applied in what follows to the filtered experimental images of two rough laminar pyrocarbons : "as prepared (AP)" (Fig 1a ) and "heat-treated (HT)" (Fig 1b). Figs 1a (and1c) and 1b (and 1d) span (in average) on respectively 16 and 17 fringes. The steerable pyramid decomposition is performed with 3 levels and 4 orientations. Autocorrelation coefficients are taken on 7x7 neighborhoods.

Atomistic reconstruction

Starting from the 3D synthetic images, cubic simulation boxes of compatible sizes (number of fringes multiplied by the suited interlayer spacing d 002 ) are created and randomly filled by carbon atoms to generate a system at the suited density. No heteroatoms are considered in this work. Typical pyrocarbon values of d 002 and of the density, respectively 3.5 Å and 2.1 g/cm3, are adopted. This corresponds to 18507 and 22198 atoms for respectively the AP and HT pyrocarbon images of Figs 1a and1b.

After a short Hard-Sphere Metropolis Monte Carlo simulation (to ensure that no C-C distance is below 1.3 Å), random initial velocities corresponding to a high temperature (8000 K) are assigned to the system and a realistic interatomic potential, the second generation REBO potential [START_REF] Brenner | [END_REF] is switched on. The REBO energy is defined as

[ ] ∑∑ = > - = N i i j ij A ij ij R ij c REBO r V b r V r f V 1 ) ( ) ( ) ( (3)
.

where N is the number of atoms, r ij is the distance between atoms i and j, f c (r) is a switching function that goes smoothly from unity for r = 1.7 Å to zero for r = 2 Å, V R (r) and V A (r) are respectively the repulsive and attractive potentials and b ij is an empirical bond order term, depending on both the environments of atoms i and j, and which properly takes into account the different bonding possibilities of carbon. This potential, which forces the system to adopt a proper carbon bonding, is complemented by an external energy term V HRTEM to impose the nanotexture of the 3D image to the system. We chose here V HRTEM to be a simple linear function of the greyscale level:

V HRTEM = k im I(r i ) i=1 N ∑ (4).
where k im is a proportionality factor and I(r i ) is the greyscale level (ranging from 0 to 1 for respectively white and black voxels) corresponding to the position r i of atom i.

A simulated annealing procedure is then launched using molecular dynamics together with the stochastic thermostatting method of Andersen [START_REF] Allen | Computer Simulation of Liquids[END_REF]. During this simulation the system is "slowly" cooled from its initial temperature down to 0 K by adjusting the target temperature of the thermostat at each timestep. Such a procedure mimics the dynamics of a carbon system cooled from the liquid state in the presence of an external field (the 3D image potential) bringing the system to a given nanotexture. Tuning the relative importance of the two terms of the potential by changing the value of k im allows to get a good balance between finding low energy structures and obtaining the suited nanotexture.

Preliminary simulations performed on smaller images (13 fringes height, not presented here) showed that a value of k im = 2 eV, approximately one fourth of the sp2 carbon energy, is very appropriate. For the cooling rate, these simulations revealed that a non uniform cooling (fast cooling at high and low temperatures and slow cooling at intermediate temperatures) was the optimal solution. In this work we adopted the following cooling sequence: 125 K/ps (T > 7000 K), 25 K/ps (7000 K ≥ T > 6000 K), 5 K/ps (6000 K ≥ T > 5000 K), 1 K/ps (5000 K ≥ T > 4000 K), 5 K/ps (4000 K ≥ T > 3000 K), 25 K/ps (3000 K ≥ T > 2000 K) and 125 K/ps (2000 K ≥ T). This temperature ramp, with the slower cooling rate (1 K/ps) around the carbon melting point gives lower final energies and higher sp 2 fractions than for instance the much more computationally expensive uniform cooling at 2 K/ps. After completion of the annealing, and in order to avoid any bias that the image potential could introduce, this term is switched off and the REBO potential is replaced by the AIREBO potential of Stuart [START_REF] Stuart | [END_REF]. This potential is a modification of the REBO potential, incorporating in an adaptive manner the van der Waals interactions occurring between non-bonded atoms (interactions between different graphitic sheets). At that point the system is relaxed for some picoseconds at a constant temperature of 300 K before calculation of properties. 1a and1c) appear to be much less extended and more curved than those of the HT sample (Figs 1b and1d).

DISCUSSION

The atomistic configurations obtained for the two materials at the end of the above described reconstruction procedure are shown in Fig. 2 As can be seen on this figure, the two atomistic models show a well pronounced laminar nanotexture, typical of pyrocarbons. Carbon bonding in these two models is almost entirely of sp 2 nature for both models (97 and 99% for respectively AP and HT pyCs), the remaining carbon atoms being mainly sp 3 (with some traces of sp) hybridized. In order to better characterize the materials, we performed an analysis of ring statistics, limited to three-coordinated atoms, according to the "shortest path ring" (SPR) algorithm of Franzblau [17]. We found that, on average, one carbon atom belongs to respectively 2.53 and 2.81 SPR for respectively the AP and HT pyCs (in graphite each atom belongs to three SPR for comparison). Moreover, a similar repartition of rings is found in the two materials with respectively 88.6 and 88.9 % of sixmember rings (C 6 ), 5.5 and 5.4 % of 5-member rings (C 5 ) and 5.7 and 5.5 % of 7-member rings (C 7 ) for respectively the AP and HT models. In both cases a low fraction (0.2 %) of eightmember rings (C 8 ) has also been observed. Although these fractions of non-hexagonal rings (around 11 %) are relatively high, we expect that lowering the annealing rate will increase the fraction of hexagonal rings, at least in the case of the HT material as heat treatment is supposed to anneal some structural defects.

To further unravel the nanostructural features of the reconstructed models, we performed a search for "pure C 6 " rings , namely those C 6 rings whose atoms belong only to C 6 rings. In Although no quantitative results are provided in this paper, visual comparison of Figs 2a and2b shows that the HT model presents a much higher degree of order than the AP one. Indeed, pure C 6 domains appear to be more extended in HT than AP. This leads to quite straight carbon sheets in the HT model as compared to the undulated ones observed in the AP model. Even if the fractions of non hexagonal rings are similar in the two materials, we can see on Going back to the very beginning of this work, we show again in Fig. 3 the HRTEM images of the AP (Fig. 3a) and HT (Fig. 3d) pyCs together with those simulated from the AP (Figs 3b and3c) and HT (Figs 3e and3f) atomistic models using the NCEMSS HRTEM simulation software of Kilaas [START_REF] Kilaas | Interactive software for simulation of high resoluton TEM images[END_REF]. c andd) is obvious and confirm that the atomistic models contain most of the nanotextural information present in the HRTEM images. Note that a high-frequency signal, corresponding to the presence of atoms, is observed in the simulated images (Figs 3b and3d) and was removed from the experimental images (Figs 3a and3b) due to filtering.

As a last check up of our atomistic models, we plot on Fig. 4 the reduced Pair Distribution Functions (PDF) computed for the AP and HT models. The reduced PDF or G(r) is defined as : ).

G(r) = 4πρ 0 r g(r) -1 [ ] (5
where r is the interatomic distance, ρ 0 is the number density and g(r) is the atomic pair distribution function defining the normalized probability of finding two atoms distant by r :

( ) 6).

∑∑ = ≠ - = N i i j ij r r N r g 1 0 1 ) ( δ ρ (
Even though it is theoretically possible to measure the PDF of pyrocarbons, the difficulty to obtain these materials (usually deposited on fibrous performs) in bulk quantity, required to perform X-ray or Neutron scattering experiments, probably explains why this function has never been published for such materials. Nevertheless, we also show, for comparison purposes, the experimental G(r) of a nanoporous disordered carbon [3] and the computed G(r) of an ideal hexagonal graphite lattice on Fig. 4. We easily observe that the two pyC models have almost identical PDFs. This is not surprising as the PDF is an orientation-averaged measure of local order, and the two models have very similar structures in the investigated length range. When compared to the nanoporous carbon and to graphite we see that the models possess the typical structure of dense carbons, with all intra-and inter-layer peaks presents. However, when compared for instance to the nanoporous carbon PDF measurements of Petkov et al. [3], we see that some peaks in the pyrocarbon PDFs are wider and less pronounced (see for instance the third neighbor peak at 2.9 Å or the peak at 6.6 Å). This can be due to the relatively high fraction (around 11 %) of non-hexagonal rings in our models while Petkov's PDF has been modeled with atomistic configuration comprising zero [START_REF] Acharya | MRS Fall Meeting Symposium NN proceedings[END_REF] to 5% [5] of such rings.

CONCLUSIONS

In this paper we have proposed a new method, based on HRTEM data, for the atomistic reconstruction of nanotextured carbons. This method has been applied to two nanotextures typical of rough laminar pyrocarbons. The first one corresponds to an as prepared material and the second one to a heat treated pyrocarbon. This new image guided method has led to relatively satisfactory models. Indeed, the nanotextures of the two models, as well as the differences between them, are fairly compatible with what we know from the experimental samples (namely, the initial HRTEM images). Moreover, on the chemical point of view, the atomistic models are typical of sp 2 carbons and present well defined aromatic domains, well stacked along the orthotropy axis. However, at that point, the number of non-hexagonal rings (11 %) is slightly larger that what we would expect (a few percents) for such dense pre-graphitic carbons. The latter point is confirmed by the observation of the PDF, showing larger peaks than for instance the nanoporous carbon studied by Petkov et al. [3]. Nevertheless, we have reasons to believe that this "chemical quality" can be neatly improved by lowering the rate of the simulated annealing process or by including in this process, together with the 3D image, a penalty term based on an experimental PDF, if available. Such an approach has been used by Jain et al. [START_REF] Jain | [END_REF] in the simpler case of low density porous carbons. Furthermore, introduction of H atoms in the simulated structures will also help in this direction.

Fig 1 (

 1 Fig 1a and 1b) is the vertical direction.

Figure 1 :

 1 Figure 1: Filtered experimental HRTEM images of an "as prepared (AP)" (a) and a "heat-treated (HT)" (b) rough laminar pyrocarbon. Subfigures (c) and (d) show the 3D HRTEM-like images synthesized from target statistics collected from subfigures (a) and (b) respectively. As can be seen in Fig 1, the synthetic 3D images of Figs 1c and 1d exhibit anisotropy properties similar to those observed in the HRTEM samples of Figs 1a and 1b. Particularly, straight and distorted fringe lengths and disorientations observed on Fig 1c (respectively 1d), though not evaluated quantitatively, seem to be quite similar to the ones observed on the HRTEM sample of Fig 1a (respectively 1b). More precisely, the fringes observed for the AP sample (Figs.1a and 1c) appear to be much less extended and more curved than those of the HT sample (Figs1b and 1d).

Figure 2 :

 2 Figure 2: Reconstructed atomistic models of "as prepared (AP)" (a) and "heat treated (HT)" (b) rough laminar pyrocarbons. Blue sticks: bonds between carbon atoms in pure C 6 rings; Orange sticks: other C-C bonds.

Fig 2 ,

 2 pure C 6 rings are shown with blue sticks while other C-C bonds are shown in orange.

  Fig 2 that such rings are rather involved in inter-sheets connections in the AP model (Fig 2a) while they are mainly present as intra-sheet disorder in the HT model (Fig 2b). The important fraction of non 2009 MRS Fall Meeting Symposium NN proceedings. November 30 -December 3, 2009, Boston, MRS Symp. Proc. (2010) vol. 1231, pp. 50-59 hexagonal rings in rather flat sheets (Fig 2b) is explained by the segregation of such rings in some disordered but flat domains with neighboring C 5 and C 7 rings like in the well-known Stone-Wales defects.

Figure 3 :

 3 Figure 3: Filtered experimental HRTEM images of an "as prepared (AP)" (a) and a "heat treated (HT)" (b) rough laminar pyrocarbon. Simulated HRTEM images obtained from the reconstructed atomistic models of the AP (c) and HT (d) pyCs.

  2009 MRS Fall MeetingSymposium NN proceedings. November 30 -December 3, 2009, Boston, MRS Symp. Proc. (2010) vol. 1231, pp. 50-59

Figure 4 :

 4 Figure 4: Room temperature Reduced Pair Distribution Function of reconstructed atomistic models of "as prepared (AP)" (full circles) and "heat treated (HT)" (straight line) pyrocarbons. Also displayed are the experimental PDF of a nanoporous carbon [3] (dashed line and empty squares) and the computed PDF of ideal hexagonal graphite (dotted line).

ACKNOWLEDGMENTS (OPTIONAL)

We thank R. Vitti for his contributions to the early stages of this project. Dr R. J.-M. Pellenq is gratefully acknowledge by JML and GLV for useful discussions and advices about this work.


N Berbel

R Fernández-García

I Gil

B Li

A Boyer

S Bendhia

Experimental verification of the usefulness of the nth power law MOSFET model under hot carrier wearout

In this paper the usefulness of the nth power law MOSFET model under Hot Carrier Injection (HCI) wearout has been experimentally demonstrated. In order to do that, three types of nFET transistors have been analyzed under different HCI conditions and the nth power law MOSFET model has been extracted for each sample. The results show that the model can reproduce the MOSFET behavior under HCI wearout mechanism. Therefore, the impact of HCI on circuits can be analyzed by using the nth power law MOSFET model.

Introduction

The nth power law MOSFET model, also namely Sakurai-Newton model, is a wellknown analytical model, which can reproduce the MOSFET behavior with only six parameters [START_REF] Sakurai | A simple MOSFET model for circuit analysis[END_REF]. The model is expressed in the following equation:

V DSAT =K•(V GS -V TH ) m I DSAT =B•(V GS -V TH ) n (1)
where V GS and V DS are the gate-to-source and the drain-to-source voltages, respectively. V TH denotes the threshold voltage, V DSAT the drain saturation voltage and I DSAT the drain saturation current. Parameters K and m control the linear region characteristics whereas B and n determine the saturated region characteristics. λ is related to the finite drain conductance in the saturated region. The simplicity of the model allows an analytical processing. The model (or models derived from it) has been extensively used in order to evaluate the behavior of several CMOS circuits [START_REF] Sakurai | A simple MOSFET model for circuit analysis[END_REF], [START_REF] Rossello | An analytical charge-based compact delay model for submicrometer CMOS inverters[END_REF], [START_REF] Fernandez | DC broken down MOSFET model for circuit reliability simulation[END_REF] and [START_REF] Nishant | Extended-Sakurai-Newton MOSFET Model for Ultra-Deep-Submicrometer CMOS Digital Design[END_REF] due to its simplicity. It should be pointed out that the nth power law MOSFET parameters correspond to an analytical non-physical model and it is only valid for voltage gates above the voltage threshold.

Hot Carriers Injection (HCI) is the phenomenon in solid-state or semiconductor electronic devices where either an electron or a "hole" gains sufficient kinetic energy to overcome a potential barrier necessary to break an interface state [START_REF] Acovic | A review of hot-carrier degradation mechanisms in MOSFETs[END_REF]. HCI does not always lead to circuit failure [START_REF] Li | SRAM circuit-failure modeling and reliability simulation with SPICE[END_REF], but provokes undesired effects such as the reduction of the speed at which the transistor can operate. Therefore, to determine how sensitive are the electronic circuits to HCI, some accurate models are required. In this sense, MOSFET models, which take into account the effects of HCI, and that could be implemented in circuit simulators are required. Equivalent electrical circuits for FETs under HCI have been already presented, which can be included in SPICE simulators. In this approach, additional components are needed [START_REF] Li | Compact modeling of MOSFET wearout mechanisms for circuit-reliability simulation[END_REF]. An alternative approach is the description of the HCI MOSFET using models like, for example, BSIM, whose parameters are extracted after the HCI stress. This approach would be the most accurate, but a large number of parameters and a complex extraction process are required. In this work, and alternative approach has been followed, and the usefulness of nth power law MOSFET model under HCI effect has been demonstrated experimentally. Moreover, the impact of HCI wearout on each Sakurai-Newton model parameters has been experimentally obtained.

Experimental

To analyze the impact of the HCI in the nMOS, a Keithley 2602A Source-meter, a digital multimeter and an acquisition board have been used, all of them controlled with a laptop (Fig. 1). First of all, it is necessary to obtain the stress conditions. The V DS stress voltage (V STRESS ) should be equal to the drain-to-source voltage breakdown minus 0.5 V. The V GS stress voltage is obtained by measuring the substrate current (I SUB ) whereas the V GS is swept with the V DS set to V STRESS . The V GS stress voltage is defined as the gate voltage corresponding to the maximum I SUB value on the I SUB -V GS curve. Table 1 summarizes the different stress voltages applied to the different type of nMOS transistors. Fig. 1. Set-up to obtain the impact of the HCI on nMOS. Table 1. Type of transistor measured and stress voltages applied.

V DS | STRESS (V) V GS | STRESS (V)

High voltage long (W = 10 μm L = 1.52 μm) 7 2.6 5 1.9

High voltage short (W = 10 μm L = 0.38 μm) 7 2.05 5 1.6

Low voltage long (W = 10 μm L = 0.4 μm) 3 0.97 2.5 0.9

Results

Fig. 2 and Fig. 3 show the ID-VDS and the ID-VGS curve characteristic, respectively, for several stress time. As it is observed, the drain saturation current is reduced due to the effect of the HCI. Fig. 4 and Fig. 5 plot the fitting of the measurements and the Sakurai-Newton MOSFET model for the fresh component and after 700 s of stress. A good correlation between the mathematical model and the measurements has been obtained. Therefore, the impact of HCI wearout on each nth power law MOSFET model has been evaluated and it is shown below.

Conclusions

In this paper, the usefulness of the Sakurai-Newton MOSFET model under HCI wearout has been experimentally demonstrated for different nFET transistors and different stress conditions. The impact of HCI on each parameter has been also analyzed. The results show that the models can be useful in order to fit the MOSFET behavior under HCI stress and it can be included in a circuit simulator to predict the impact of HCI on digital circuit behavior, where the sub-threshold behaviors are not so important.

Fig. 2 .

 2 Fig. 2. I D -V DS for several stress time and several V GS .

Fig. 4 .

 4 Fig. 4. Comparison between the measurements and the Sakurai-Newton MOSFET model for the fresh component.

Fig. 5 .

 5 Fig. 5. Comparison between the measurements and the Sakurai-Newton MOSFET model for the t STRESS = 700 s.

Fig. 6 Fig. 6 .

 66 Fig.6depicts the increase of the finite drain conductance (λ), which produces a change in the slope of the drain current (in the saturated region), as can be observed in Fig.2. Fig.7plots the threshold voltage shift, which is more significant on the high voltage short transistor, rather than on the high voltage long and low voltage long. The VTH enhancement provokes a drain current decrease, as can be observed on Fig.3. The n parameters shift is plotted in Fig.8, which corresponds to a value less than 10% for the three types of nMOS transistors. On the other hand, B parameter shift is plotted in Fig.9, which decreases as long as the stress time increases. This parameter is related with the transconductance of the device. Fig.10and Fig.11show the K and m parameters shift, related with control of the linear region. Both parameters are modified due to the change of the linear ID-VDS curve characteristic, as seen in Fig.4and Fig.5.

Fig. 7 .

 7 Fig. 7. Variation of the threshold voltage versus the stress time.

Fig. 8 .

 8 Fig. 8. Variation of the n parameter versus the stress time.

Fig. 9 .

 9 Fig. 9. Variation of the transconductance (B) versus the stress time.

Fig. 10 .

 10 Fig. 10. Variation of the m parameter versus the stress time.

Fig. 11 .

 11 Fig. 11. Variation of the K parameter versus the stress time.

Acknowledgements

This work has been supported by the Spain-MICINN under Projects TEC2010-18550 and AGAUR 2009 SGR 1425. The authors want to thank to B. Vrignon from Freescale Semiconductors for sample provision.


Binhong Li
email: binhong.li@insa-toulouse.fr

Alexandre Boyer

Sonia Ben Dhia

Christophe Lemoine

S Bendhia

Ageing effect on electromagnetic susceptibility of a phase locked loop

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

During their lifetime, integrated circuits (ICs) can be affected by failure mechanisms such as electro migration, gate oxide break down or hot carrier injection [START_REF] White | Microelectronics Reliability: Physics-of-Failure Based Modeling and Lifetime Evaluation[END_REF], mainly activated by harsh environmental conditions such as high or low temperature and electrical overstress. Even if failure mechanisms do not compromise the circuit operation, IC intrinsic degradations can have a significant impact on performances, such as noise sensitivity, operating frequency… or Electromagnetic Compatibility (EMC). [START_REF] Weir | Ultra-Thin Gate Dielectrics: They Break Down, But Do They Fail?[END_REF] The concept of electromagnetic reliability (EMR) [START_REF] Li | Characterization of the electromagnetic robustness of a nanoscale CMOS integrated circuit[END_REF] has been developed to study and evaluate the impact of circuit ageing on EMC. Analog circuits as PLL are very sensitive to electromagnetic interferences (EMI) [START_REF] Redoute | Design of EMI resisting analog integrated circuits[END_REF]. Changes of transistors physical characteristics (such as threshold voltage and mobility) and mismatches due to intrinsic degradations could affect the susceptibility of analog circuits. In this paper, we focus on the evolution of the susceptibility of a CMOS Voltage-Controlled Oscillator (VCO) to EMI coupled on the power supply after accelerated ageing, as power supply fluctuations of VCO is the dominant noise source in PLL. In the section II the principles of EMR characterization is described. EMC measurement methods and accelerated life tests are given in section III. Section IV presents the device under test (DUT). In the fifth part, the immunity levels measurement before and after accelerated ageing test are compared and the EMI-induced failures are detailed. Finally, experimental results are summed up and discussed.

EMR Characterization principles

The characterization of the EMR of a circuit consists in measuring its emission and/or susceptibility level before and after an ageing procedure which accelerates intrinsic degradation mechanisms. As illustrate in Fig. 1, it aims at extracting the EMC level drifts induced by ageing as accurately as possible.

Once the EMC and the ageing procedure set-ups have been completely defined, they have to be validated experimentally. The validation step must prove that ageing procedure can only have an impact on the DUT. The EMC validation step must ensure that the EMC measurements are enough repeatable to extract precise information concerning the EMC drifts related to ageing effects. Measurement repeatability and uncertainties set a limit for the consistence of EMC level drifts.

At the end of the measurement campaign, the EMC measurement results obtained before and after ageing are processed to extract statistical data concerning the EMC level drifts, such as the case and the mean drifts. These data are required to predict the risk that a component becomes incompliant after ageing and, if necessary, readjust the EMC margins at the design level.

Direct Power Injection Set-up

The direct power injection (DPI) [START_REF]IEC 61967, Integrated Circuits, Measurement of Electromagnetic Emissions, 150 KHz -1 GHz[END_REF] is a standard and efficient method to couple a conducted disturbance to a pin of a circuit (Fig. 2). It consists in superimposing RF disturbances to a low frequency signal through a decoupling network. Uncertainties of the method are acceptable to ensure consistency of EMR results.

Susceptibility is characterized over the band 10 MHz to 1 GHz, the maximum forward power is set to 45 dBm. Three susceptibility criteria are considered: a static margin of the output amplitude set to 20 % of the power supply Vdd, a dynamic margin set to 20% of the output signal period (8ns) and a power current limited to 40 times the nominal current.

Burn-in Test Conditions

According to the standard AEC-Q100 [START_REF]Stress test qualification for integrated circuits[END_REF], accelerated life tests are required to qualify integrated circuits and guarantee the quality and the robustness of automotive applications. These tests consist in applying an overstress conditions during a short time in order to accelerate the damage rate for relevant degradation mechanisms. As ICs embedded in automotive applications often suffer from high temperature, High Temperature Operating Life (HTOL) which consists in applying overstresses conditions.

Ten components are tested and placed during 408 hours in a climatic chamber which regulates the temperature at 150°c +/-2.5°c. An external source provides a power supply 10% higher than the normal supply voltage. Figure 3

Device under test

The PLL is based on a delay controlled ring oscillator designed in 0.25µm SMARTMOS 8 technology from Freescale Semiconductor (Toulouse -

Experimental results analysis

Validation of the experimental set-up

First, the repeatability of EMC test bench is evaluated for validation purpose by repeating the measurement 10 times on the same sample. The test bench repeatability is 0.27dB which is highly acceptable since the repeatability errors are inferior to the uncertainty of the DPI test bench, estimated to 2 dB. Conducted EMI injections are performed on the power supply pins of the three PLL sub-blocks. Fig. 5 presents the susceptibility threshold for injection done on three separated power rails: phase comparator, VCO, and divider.

PLL Susceptibility Analysis

Injection on PLL power supply

The comparison among three immunity curves shows that the VCO is the most sensitive block in PLL to the power supply noise.

VCO Susceptibility Analysis

VCO's output frequency generally depends on power supply voltage. The variation on power supply can change the propagation delay, thus changing the oscillation frequency in VCO. The sensitivity of VCO to the power supply voltage can be extracted by measuring the output frequency of the VCO for different power supply voltage values (from 2V to 5.5V). Fig. 6 shows the VCO output frequency as a function of its static power supply

VCO output sensitivity

VCO Power supply [V] VCO output Freq [MHz]

Fig. 6. Dependence of VCO output frequency to power supply voltage.

Any RFI coupled on VCO power supply can modulate the VCO output frequency at the RFI rhythm, resulting in a jitter added to the PLL output signal.

Power supply noise amplitude[mV] VCO output frequency[MHz]

Fmax-Fmin Fmean Fig. 9. PLL output spectrum without (on the left) and with RF injection Fig. 9 plots the impact of a RFI on VCO power supply to PLL output frequency. RFI induce a considerable jitter and tends to increase the VCO output frequency.

PLL immunity Level Drift

With evaluating ageing effect on power supply immunity of three sub-blocks, the perturbation on VCO power supply exhibits a significant evolution of PLL susceptibility threshold after accelerated ageing. Fig. 9 presents the susceptibility threshold of the 6th sample for injection done on phase comparator and VCO power supply pins before and after the HTOL test. The two figures show that VCO immunity is more sensible to accelerate ageing. Its immunity level reduction can reach up to 10 dB over a large frequency range. A statistical analysis shows that among 10 samples after HTOL the maximum reduction is 15.88 dB; while the mean reduction for all 10 samples over all frequency is 2.57dB with a standard deviation 3.64dB. Moreover, the process dispersion among 10 samples after HTOL changes from 1.64dB to 4.83 dB. This result indicates that the VCO and thus the PLL becomes more sensitive to power supply fluctuations after ageing. After ageing, a shift of the power supply voltage range ensuring the PLL locking has been measured (fig. 11).

VCO output sensitivity before and after ageing

Discussion

After HTOL ageing tests, the functionality of PLL on all the tested samples is preserved. The current consumption and the performances are unchanged after ageing, proving the excellent robustness of the Freescale SMARTMOS 8 technology. However, acceleratedageing test tends to decrease significantly the immunity levels of VCO and can strongly affect the failure risk of an application which uses this PLL. Moreover, the dispersions in terms of susceptibility between components tend to increase after ageing, indicating that the intrinsic degradation mechanisms follow different evolution in each sample. As HTOL mainly accelerates TDDB and NBTI which can degrade transistor parameters, which could constitute the source of the immunity level variation after ageing. For example, the degradation of transistor transconductance (g m ) and mobility (µ) could directly affect the power supply rejection ratio and delay time of gates composing the VCO. Electrical simulations are currently leaded to confirm that the immunity reduction is induced by both a degradation of power supply noise rejection and VCO switching characteristics, provoked by changes of MOS transistor physical parameters. A model including the ageing parameters is going to be studied.

Conclusion

In this paper, an original study of the impact of accelerated ageing tests on PLL was presented. The results show that the functionality of circuit was preserved, but the electromagnetic susceptibility of VCO is affected by accelerated ageing test based on high temperature and electrical overvoltage conditions. The immunity level is reduced significantly over a large frequency range, which can increase interference risks with neighbor circuits. The correlation between the immunity level drift and the intrinsic degradation mechanisms needs to be studied in detail to understand the origins of the drift and predict in which extent the susceptibility threshold can decrease.

Fig. 1 .

 1 Fig.1. General description of the electromagnetic robustness of a component 3. Experimental set-up and device under test

Fig. 2 .

 2 Fig. 2. DPI setup, RF injected on power supply rail of VCO

  describes the accelerated life test set-up. The HTOL stress mainly accelerates wearout failure mechanisms as Time Dependent Dielectric Breakdown (TDDB) and Negative bias temperature instability (NBTI).

Fig. 3

 3 Fig. 3 Accelerated life test set up The experimental set-up uses three boards (Tab 1). Dedicated sockets are used on the EMC and ageing boards to move the circuit under test. With this set-up, only the circuits under test are aged.

  dedicated to automotive applications. The PLL is divided in three sub-blocks: VCO, fraction 1/4 frequency divider and phase comparator. The PLL operates at 24 MHz. Power supplies of these sub-blocks have been separated to identify their own sensitivities to RFI. (Fig.4) (a) (b) Fig. 4. (a) VCO schematic and (b) PLL schematic, it three separated power rails for its three sub-blocks

Fig. 5 .

 5 Fig. 5. PLL susceptibility threshold during conducted injection on power supply pins of the three sub-blocks (Divider, VCO, Phase comparator)

Fig. 7 .Fig. 8 .

 78 Fig. 7. PLL output signal without (on the left) and with RF injection

Fig. 9 .

 9 Fig. 9. PLL susceptibility threshold with an injection on (a) Phase comparator and (b) VCO power rails before and after ageing

Fig. 10 .Fig

 10 Fig.10. Dependence of VCO output frequency to power supply voltage before and after ageing

Number of samples EMC tests Aging stress Test set-up optimization (Time, cost, accuracy, repeatability) Experimental set- up validation EMC level measurements Aging stress Data processing Maximum EMC level drifts Mean EMC level drifts EMC level dispersion drift Individual and statistical analysis Outputs EMC failure rate (risk of non compliance for the full lifetime) Readjustment of EMC margin Experiments EMC set-up validation Aging procedure validation

Table 1 -

 1 Three dedicated test boards

		with a RF Socket to measure
	EMC board	susceptibility level before and after
		accelerated ageing
	Ageing board	placed in the climatic chamber with the sockets special for ageing
	Monitoring	supply and monitor components under
	board	test during accelerated ageing

Acknowledgements

The authors would like to thank Freescale Semiconductor (Toulouse -France) and ESEO (Ecole supérieure d'Electronique de l'Ouest, ANGERS -France) for their technical support and their precious advices. This study has been done with the financial support of French National Research Agency (ANR) and the regional council of Midi-Pyrénées.