E D Ngounou Ntoukam

V C Kamla

Y Mualla

I H Tchappi

S Galland

J C Kamgang

S Emvudu Wono

Towards a multi-agent model to prevent damage caused by cocoa mirids to cocoa pods

Keywords: Mirid, Agent-based modeling, Cocoa, GAIA methodology Mirid, Modélisation à base d'agents, Cacao, Méthodologie GAIA

Agriculture in general and cocoa farm in particular is one of the primary sources of income for several African countries. Although quite useful for the development of Africa, cocoa production is faced with several challenges like the diseases of cocoa trees. Mirid (Sahlbergella singularis) is the worst pest affecting cocoa production across west Africa causing 30% to 40% of damage on cocoa production. It is difficult to control the mirids' population and its impact on the cacao production. Related works on mirids studied the mirids population as a macroscopic entity without considering local interactions of individual mirids. These models did not consider the behaviours and the environment of mirids which are essential due to their complexity. To overcome these drawbacks, this paper proposes an agent-based model of mirids considering the individual interactions of mirids within the environment. The proposed model consists of agents representing the mirids through their life cycle in cocoa farm environments. It is based on the biological and ecological partial knowledge found in the literature, and is built according to the GAIA methodology. The main assets of this model are the reproduction of the dynamics of the mirids on a cocoa pod and damage caused by them at the micro level for human users (entomologists, agricultural engineers, farmers and decision makers) to thoroughly understand the phenomena and the emergent properties outcomes.

Résumé

L'agriculture en général et la culture du cacao en particulier est l'une des principales sources de revenus de plusieurs pays africains. Cependant, la production du cacao est confrontée à plusieurs défis tels que les maladies et l'attaque des insectes nuisibles. Les mirides (Sahlbergella singularis) sont les insectes les plus préjudiciables à la culture du cacao en Afrique. Avec 30 à 40 % de dommages sur la production, il est difficile de contrôler la population des mirides ainsi que son impact sur la production de cacao. Les

travaux précédents sur les mirides, étudient la population de cet insecte selon le point de vue macroscopique sans tenir compte des interactions locales entre les mirides. Ces modèles ne prennent pas en compte les comportements individuels des mirides et leurs interactions avec l'environnement qui sont essentiels en raison de leurs complexités. Pour remédier à ces insuffisances, le présent article propose un modèle de mirides à base d'agents prenant en compte les interactions individuelles des mirides dans leur environnement. Le modéle proposé est constitué d'agents représentant le cycle de vie du miride sur une cabosse de cacao. Il est basé sur les connaissances biologiques et écologiques du miride, et est construit selon la méthodologie GAIA. Ces modèles ainsi que leurs implémentations sont destinés à des utilisateurs tels que les entomologistes, les ingénieurs agronomes, les agriculteurs et les décideurs. Cela leur permettra de mieux comprendre les phénomènes émergents dues aux mirides ainsi que leurs propriétés. 1 Introduction

With more than 72% of worldwide cocoa production, Africa (Ivory Coast, Ghana, Nigeria and Cameroon) is currently the main producer of cocoa [START_REF] Organization | Icco quarterly bulletin of cocoa statistics[END_REF]. Cocoa is cultivated by smallholder with lower livelihood, yet it is an important source of revenue for these farmers.

African production faces numerous difficulties like : ageing of farmers, ageing of plantations, non-respect of good agricultural practices, presence of diseases and pests [START_REF] Fountain | Baromètre du cacao[END_REF]. Cocoa Swollen Shoot Virus (CSSV) and cocoa Black Pod Rot (BPR) are the two diseases that most seriously affect African cocoa production [START_REF] Babin | Isolation and characterization of twelve polymorphic microsatellite loci for the cocoa mirid bug sahlbergella singularis[END_REF]. Additionally, Sahlbergella singularis and Distantiella theobroma are the two pests that most negatively impact production of cacao in all African producer countries [START_REF] Babin | True Bug (Heteroptera) Impact on Cocoa Fruit Mortality and Productivity[END_REF]. Between these two species, Sahlbergella singularis is the most common and most harmful [START_REF] Babin | True Bug (Heteroptera) Impact on Cocoa Fruit Mortality and Productivity[END_REF]. These two species, also called mirid bugs, feed principally on cocoa (Theobroma cacao). Mirids feed on pods and shoots and they cause the drying up of leaves and a quick destruction of the cocoa canopy. The immature (five instars) and adult stages of these sucking insects feed on the sap of young branches, by injecting a digestive saliva [START_REF] Babin | Contribution à l'amélioration de la lutte contre le miride du cacaoyer Sahlbergella singularis Hagl.(Hemiptera : Miridae). Influence des fac-teurs agro-écologiques sur la dynamique des populations du ravageur[END_REF].

Due to the considerable thickness of the pod cortex, damages caused by mirids are relatively superficial on pods and there is no direct link between the intensity of feeding lesions on mature pods and bean quality [START_REF] Babin | True Bug (Heteroptera) Impact on Cocoa Fruit Mortality and Productivity[END_REF]. Nevertheless, young fruits (less than 10 cm in length and less than 2 months old) also called cherelles, can be seriously damaged by the mirid's feeding lesions. These feeding lesions can increase the abortion of young fruits and thus impact production [START_REF] Babin | Isolation and characterization of twelve polymorphic microsatellite loci for the cocoa mirid bug sahlbergella singularis[END_REF].

The main strategy to fight against mirids is the use of insecticide. Although alternative cocoa pest control methods have been developed to fight against mirids. They include cultural management [START_REF] Mahob | Field evaluation of synthetic sex pheromone traps for the cocoa mirid Sahlbergella singularis (Hemiptera : Miridae)[END_REF], varietal management [START_REF] Anikwe | Evaluation of resistance in selected cocoa genotypes to the brown cocoa mirid, Sahlbergella singularis Haglund in Nigeria[END_REF], as well as semi-chemical management with the use of synthetic sexual pheromone traps [START_REF] Anikwe | Evaluation of resistance in selected cocoa genotypes to the brown cocoa mirid, Sahlbergella singularis Haglund in Nigeria[END_REF][START_REF] Sarfo | Design and placement of synthetic sex pheromone traps for cacao mirids in ghana[END_REF]. However, their effectiveness remains limited. Cocoa is cultivated in tropical agro-forests which are characterized by complex associations of multifunctional and uneven-aged trees and crops [START_REF] Ngo Bieng | Diversity and spatial clustering of shade trees affect cacao yield and pathogen pressure in Costa Rican agroforests[END_REF]. In such environments, it is not obvious to count mirid individuals (immatures and adults) or evaluate the impact of their actions on cocoa plant. Agent-based approaches are suitable to model complex systems [9, 27? ]. They seem to be the most relevant option to model the dynamic behavior of mirids at the micro level in order to understand mirids behavior and thus forecast its evolution and adopt an efficiency control strategy. The aim of this work is therefore to propose an agentbased model for the mirids dynamic. To this end, an initial agent-based model is proposed based on the mirid's life cycle. We make some representations of a mirid agent on a cocoa pod to present how it progresses. In order to be logical and clear, we follow the GAIA methodology [START_REF] Zambonelli | Developing multiagent systems : The gaia methodology[END_REF]. In fact, GAIA methodology is wellsuited to design and build systems in complex and open environments at both micro and macro levels [START_REF] Cossentino | Handbook on agent-oriented design processes[END_REF] such as the motion and evolution of mirids in a cocoa farm. This work is ongoing and the paper focuses only on the interactions of mirids in the environment. In the future, other interactions will be considered such as the interaction with human users (entomologists, agricultural engineers, farmers and decision makers). Therefore, we shall make sure that the target people will be able to understand our models and simulation results. To do so, the simulation should explain to the users about the behaviours of all agents and their interactions. The paper is organized as follows : Section 1.1 highlights a brief overview on theoretical models of mirids. Section 2 presents our agent-based model of the life cycle of mirids based on the GAIA methodology. Section 3 discusses the need for explainability and how to achieve it and section 4 concludes the work and gives some prospective directions.

Related Work

It is important to understand mirid's life cycle and to have information about its biology and ecology. The life cycle of a mirid is of three phases : egg, nymph and adult. Mirids develop on the cocoa pod. Females insert eggs in cocoa living tissues, and egg development requires healthy tissues to be completed [START_REF] Babin | Rearing method and life-table data for the cocoa mirid bug Sahlbergella singularis Haglund (Hemiptera : Miridae)[END_REF]. The embryonic development of a mirid is between 15 and 20 days in average with minimum 9 days and 21 days maximum [START_REF] Babin | Rearing method and life-table data for the cocoa mirid bug Sahlbergella singularis Haglund (Hemiptera : Miridae)[END_REF]. After the embryonic phase, the egg becomes a nymph. On average, the life-cycle duration (eggs to adults) of mirids is 40 days with a minimum of 36 days [START_REF] Babin | Effects of cocoa tree phenology and canopy microclimate on the performance of the mirid bug Sahlbergella singularis[END_REF] and a maximum of 50 days [START_REF] Babin | Rearing method and life-table data for the cocoa mirid bug Sahlbergella singularis Haglund (Hemiptera : Miridae)[END_REF]. The rate of hatching in the laboratory is between 30% and 82%, but it could be up to 96.53%. During the nymph stage, the mirid moves on the cocoa tree to feed of the pods and the shoots. The rate of survival of nymphs is around 68% [START_REF] Babin | Contribution à l'amélioration de la lutte contre le miride du cacaoyer Sahlbergella singularis Hagl.(Hemiptera : Miridae). Influence des fac-teurs agro-écologiques sur la dynamique des populations du ravageur[END_REF]. The sex ratio is 1 female for 0.7 male. After the last phase of nymph. a female mirid (S. singularis) starts to mate after 6 to 10 days while the male mirid starts after 2 days [START_REF] Babin | Rearing method and life-table data for the cocoa mirid bug Sahlbergella singularis Haglund (Hemiptera : Miridae)[END_REF]. After the mating, a female mirid is considered to be mature, and it can lay her eggs in many pods at different times. To lay an egg, a female mirid needs 4 to 8 days, on average, after the mating. More than 72.1% of female mirids become mature (able to lay eggs) [START_REF] Mahob | Field evaluation of synthetic sex pheromone traps for the cocoa mirid Sahlbergella singularis (Hemiptera : Miridae)[END_REF]. In average, the fecundity of female mirids is 52.5 eggs, and the fecundity period lasts, on average, for 16 days [START_REF] Tapi | Models for Miridae, a cocoa insect pest. Application in control strategies[END_REF]. The daily survival of mirids adult is around 98.14% for immature females, 92.8% for mature females, and 93% for males [START_REF] Tapi | Models for Miridae, a cocoa insect pest. Application in control strategies[END_REF]. To obtain a new generation, 50 to 60 days are needed on average [START_REF] Mahob | Field evaluation of synthetic sex pheromone traps for the cocoa mirid Sahlbergella singularis (Hemiptera : Miridae)[END_REF]. In order to understand and forecast the mirid population dynamics, several studies have been done using a variety of approaches and techniques : development of a rearing method [START_REF] Babin | Rearing method and life-table data for the cocoa mirid bug Sahlbergella singularis Haglund (Hemiptera : Miridae)[END_REF], sheeting method [START_REF] Kouamé | Variations saisonnières des populations de mirides du cacaoyer dans la région de l'indénié-djuablin en côte d'ivoire[END_REF], phe-romone trapping [START_REF] Sarfo | Design and placement of synthetic sex pheromone traps for cacao mirids in ghana[END_REF][START_REF] Mahob | Field evaluation of synthetic sex pheromone traps for the cocoa mirid Sahlbergella singularis (Hemiptera : Miridae)[END_REF] and mathematical modelling [START_REF] Tapi | Models for Miridae, a cocoa insect pest. Application in control strategies[END_REF]. Babin et al. [START_REF] Babin | Rearing method and life-table data for the cocoa mirid bug Sahlbergella singularis Haglund (Hemiptera : Miridae)[END_REF] developed a rearing method for cocoa mirids in order to assess the performance using the calculation of life-table parameters and to investigate the relationship between rearing performance, especially female fecundity, and climatic conditions. Sarfo et al., [START_REF] Sarfo | Design and placement of synthetic sex pheromone traps for cacao mirids in ghana[END_REF] studied the potential of sex pheromone trapping as an alternative environmentally-acceptable method for managing the mirids. Based on the behavioral responses of the mirids to pheromones in traps, parameters were standardized for efficient performance of the traps.

In general, field studies like those of Sarfo et al. [START_REF] Sarfo | Design and placement of synthetic sex pheromone traps for cacao mirids in ghana[END_REF] and Babin et al. [START_REF] Babin | Effects of cocoa tree phenology and canopy microclimate on the performance of the mirid bug Sahlbergella singularis[END_REF] focus on finding correlations between observed data (mirid populations) and explanatory variables (shade, temperature, rain, pheromone, etc.). Linear regression helps to determine the relationship between variables. Such statistical models can help to explain or forecast relationships between parameters. The limits of such field studies are that they are expensive, and do not allow to test many hypotheses. Additionally, the quality of the acquired data is sometimes poor [START_REF] Mahob | Field evaluation of synthetic sex pheromone traps for the cocoa mirid Sahlbergella singularis (Hemiptera : Miridae)[END_REF]. Moreover, the laboratory data and field data are not always accordance due to the fact that, in the field, it is plausible that not all parameters that influence mirids are considered. Field data are often site-specific, and extrapolation of their results to other regions is not obvious [START_REF] Duarte | The limits to models in ecology[END_REF][START_REF] Grimm | Mathematical models and understanding in ecology[END_REF]. In field studies, the combination of many parameters can also pose a severe challenge. A theoretical modelling approach or representation of such an environment could be useful to explore how the mirids live and which parameters are important to them.

Tapi et al. [START_REF] Tapi | Models for Miridae, a cocoa insect pest. Application in control strategies[END_REF] developed a generic mathematical model for mirid populations to better predict its dynamics in a plot under different management strategies. Their first model [START_REF] Tapi | Models for Miridae, a cocoa insect pest. Application in control strategies[END_REF] is a compartmental model with constant or periodic parameters. It is developed based on the mirid life cycle. They have done a theoretical and sensitivity analysis to provide numerical simulations. Following this, Tapi et al. [START_REF] Tapi | Models for Miridae, a cocoa insect pest. Application in control strategies[END_REF] developed a second model (a time-delayed model) based on the previous one. The second model was better in terms of considering the developmental time at different stages. After a brief theoretical analysis, they provided numerical simulations with/without control efforts. Mathematical models allow predicting the long-term dynamics of the studied systems. However, the variability between individuals is difficult to represent [START_REF] Huynh | Coupling equation based models and agent-based models : example of a multi-strains and switch SIR toy model[END_REF]. Moreover, mirids live in traditional cocoa agro-forests which are complex systems where interactions and components are multiple and different [START_REF] Malézieux | Mixing plant species in cropping systems : concepts, tools and models : a review[END_REF]. Although, mathematical models proposed in the literature are easy to prove, they can be difficult to explain to nonmathematicians. Therefore, they are not the best tools to communicate with certain target groups, such as farmers [START_REF] Moglia | A review of Agent-Based Modelling of technology diffusion with special reference to residential energy efficiency[END_REF]. The design of each model is founded upon approximations of the pest's biology and behavior on a population level [START_REF] Grimm | Mathematical models and understanding in ecology[END_REF]. The models presented in the literature are often dedicated to specific farms, and local interactions of individual mirids are not represented or incorporated explicitly into the population dynamics. In order to better understand the mirid dynamics, it would be useful to consider the local interactions of mirids in space and time. Thus, an agent-based approach seems to be suitable to deal with mirid motion and local interactions. Agent-based modeling is the computational study of social agents interacting in an autonomous manner as evolving systems. It allows the study of complex adaptive systems and facilitates investigations into how macro-phenomena develops from micro-level behavior among heterogeneous sets of interacting agents [START_REF] Janssen | Internet encyclopedia of ecological economics : Agent-based modelling[END_REF]. Mirids are governed by their biological preferences and limitations as documented in the literature. Therefore, and by representing mirids as individual agents, the resulting relationships between these agents can be used to predict population dynamics of the pest. Additionally, this prediction is more realistic over space and time (based on local interactions) than the ones with the population aggregated mirid models. In the context of our work where users are in the loop, we think that there is a need for a simplified way to explain the system results to users, and to interact with them and acquire feedback. In our opinion, simulation and visualization are a great way to explain the system and its results, while questionnaires and direct feedback are important to make sure the user in the loop understands the process.

Multi-agent Model of Mirids

The modelling of our system is based on the second version of GAIA proposed by Zambonelli et al. in [START_REF] Zambonelli | Developing multiagent systems : The gaia methodology[END_REF]. They have extended the version of Wooldridge et al. [START_REF] Zambonelli | Developing multiagent systems : The gaia methodology[END_REF][START_REF] Wooldridge | The gaia methodology for agent-oriented analysis and design[END_REF] towards social issues and open systems. GAIA prescribes following an ordered sequence of steps, an ordered identification of a set of organizational models, and an indication of the interrelationships guiding developers towards the development of a Multiagent System (MAS). The GAIA methodology starts first with the analysis phase to collect and organize the specifications which are the basis for the design of the computational organization [START_REF] Zambonelli | Developing multiagent systems : The gaia methodology[END_REF]. Then, the architectural design phase defines the most proper organizational structure for the MAS, the topology of interactions in the MAS, and the control regime of these interactions. Finally, the detailed design is responsible for identifying, eventually, the agent model and the services model, which act as guidelines for the actual implementation of agents and their activities. In order to emphasize the miscroscopic behavior of a mirid and for reasons of simplicity, we assume that the mirids live and grow on one cocoa pod. However an extension of this work could consider several cocoa pods in a 3D dimension. A cocoa pod constitutes the environment like showed in Figure 1. In the mentioned figure, black points represent the places where mirids have fed themselves. Mirid's position is represented in a 2D plan. This paper focuses on the modeling of mirids interactions through the environment without considering external interactions such as humans.

Analysis

Environment model : For simplicity reasons, we assume that the environment in our model is restricted to the cocoa pod where the mirids live, feed themselves, grow, mate and die. Table 1 presents the environment model for our work. We have identified some global variables that influenced the dynamics of the system. Climate parameters and shade impact the growing of the mirid and the

Preliminary role model :

The aim of this model is to describe the whole activity and procedure that can be done by a role. We identify a role in our system (MIRID). Table 3 presents the preliminary role model of our system (emphasized represented activities).

To execute the MIRID role, it is important to know at which position the mirid is situated in the plan. Its actions can impact the environment by destroying a part of the pod or by occupying a set of points. A liveness property of the MIRID role specifies that a mirid, after its birth, must grow. After certain progression, it can start moving and feed itself. When it is mature, it can mate and if it is a female mirid it can lay eggs.

At the end, it will die eventually. In our system, the role MIRID must guarantee that we have more than two mirids (at least both genders). For mating, female and male must be at the same position. To move in the system, there should be an empty position inside the whole system.

Preliminary interaction model :

The preliminary interaction that can be observed here is the mating between a male mirid and a female one. This interaction starts with the stigmergy 1 communication. The female mirid emits a pheromone on the pod. When the male mirid perceives the pheromone, it will move to look for the female. When the position of the female mirid and the male mirid is the same, the mating starts. Figure 2 presents the preliminary interaction in our model.

Figure 2 -The mating interaction protocol

Organizational rules : According to the authors of GAIA version 2 in [START_REF] Zambonelli | Developing multiagent systems : The gaia methodology[END_REF], the main goal of organizational rules is to present general relationships between roles, between protocols, and between roles & protocols. Organizational rules have two properties : the liveness property which defines how the dynamics of the organization should evolve over time and the safety property which defines the time-independent global invariants of the organization that must be respected. The safety property of our system is to ensure the continuity of the mirid species. For this, the mirid role is played indefinitely. The liveness property of our system is to ensure that the population evolves over time. For this, the mirids must mate to create a new generation of mirids. Table 4 shows the different organizational rules of our system.

MIRID

To lay an egg or to grow, there is a need for mating safety rules

M irid ω

The role MIRID can be played indefinitely tions of an agent leave signs in the environment. These signs are perceived by itself and other agents, which determine their next actions [START_REF] Parunak | How to calm hyperactive agents[END_REF]

Architectural Design

Organizational structure : From the biological description of the mirid on a cocoa pod, we define the organizational structure of our system. A peer relationship is the type of relations inside our system. It means that there is an indirect conflict on food. All agents of our system collaborate to keep the continuity of their population. Whatever the mirid (i or j) of our system, it collaborates with the others. The male and female mirids also collaborate ( --→ peer). The lack of competition is due to the abundance of food. This can be simply expressed by Equation 1.

∀(i, j) ∈ M irid 2 , i = j ⇒ i --→ peer j (1)

Role model :

After having the organizational structure, it is possible to complete the role model. This model presents the internal state of agents in the system. State-transition of the role is showed in Figure 3. It defines the specification of activities and services of the basic role and the organizational roles according to the above adopted organizational structure. Interaction model : The goal of this part is to identify new protocols or interactions deriving from the adopted organizational structure. Some interactions, notably laying eggs (Figure 4) and feeding (Figure 5), result from the organizational topology and the identification of agents.

Detailed Design

Agent model : Building an agent model consists, essentially, of establishing a link between the roles and the agent. It indicates how many instances of each agent will exist in the system, and which roles each agent will implement. It permits us to specify how many instances of each class will appear in the MAS using the cardinality standard [START_REF] Cossentino | Handbook on agent-oriented design processes[END_REF]. Our system has only one role named MIRID which can be played once or many times by the agent MIRID. Figure 6 shows our agent model.

Service model : It identifies the services associated to each agent or class, with the role to be played by each agent class. Generally, the protocols, activities and liveness expressions of the role allow us to find different services that will be implemented by an agent. Agent performs at least one service, it is necessary to document the properties of services . Table 5 specifies our service model.

Discussion

The previous section presents our agent-based model of the mirid on a cocoa pod, which is a complex system. In order to implement such a complex system, agent-based simulation frameworks have been developed to allow transferring real-world scenarios into executable models (i.e. simulating the activities of mirids in a digital environment) [START_REF] Mualla | Comparison of agent-based simulation frameworks for unmanned aerial transportation applications[END_REF]. To implement the model using agent-based simulation, some requirements should be respected :

-Specify the characteristics of our environment.

The size of the pod cocoa depends on the cocoa variety and it extremely varies due to climate fac- tors.

-Define the surface of the cocoa that should be considered according to the cocoa variety. The surface of the cocoa increase with time. -Define the position of our initial mirid. When the sun appears, mirids concentrate on the shade-less areas, they disperse themselves when the night comes. For the moving of mirids, we shall define a set of points where the mirids can go to if they are empty. Mirids move only when the sun appears and they disappear for mating or if there is an enemy around them. -For the reproduction, the female mirid emits a pheromone that will attract the male mirid. -Define the probability of death of mirids at any stage of the development. All those steps of implementation will help us implement our model using agent-based simulation. Our models, simulations and interfaces should be explainable to human users (entomologists, agricultural engineers, farmers and decision makers) in the cocoa domain. Explaining to users how intelligent systems come to their conclusions is an area of research with a history in expert systems and recommended systems [START_REF] Gilbert | Explanation and dialogue[END_REF][START_REF] Kass | The need for user models in generating expert system explanations[END_REF]. Recently, as intelligent systems are getting increasingly complex. eXplainable Artificial Intelligence (XAI) reemerged to push for interpreting the process of the machine learning mechanisms, and understanding the decisions of robots and agents [START_REF] Anjomshoae | Explainable agents and robots : Results from a systematic literature review[END_REF][START_REF] Calvaresi | Explainable multi-agent systems through blockchain technology[END_REF]. This work is on progress, and in the future we want to add the different actions of different human users (entomologists, agricultural engineers, farmers and decision makers). We need to measure how the regularity of their actions can impact the damage caused by mirids. Therefore, the main goal of our work is to provide tools (simulators and models) for the users. The tools make a recommendation, decision, or action. Then, they provide an explanation to the user that justifies its recommendation, decision, or action. Entomologists will use these tools for the preparation of their experimentation on the fields. It will help them have an idea of the type of results that they could have. The usefulness of our tools for agricultural engineers and farmers is that it will help in the pest surveillance and the forecasting of optimal timing for insecticide spraying. Farmers of cocoa are mostly poor, so it is important to them to have optimal timing for insecticide spraying. At the end of our simulations, we could estimate the damage caused by the mirids to the cocoa pods. This type of information is very important for decision makers and the government. Explainability has to be handled in two steps in our system (Figure 7) :

1. From the cocoa literature, some works have been done on the mirid, especially by entomologists.

We use these works to build our agent-based models. After the modelling, implementation outputs are interpreted by graphs. At this step, we need the interaction of entomologists to validate our models. They will be able to evaluate if our results are close to what they often observe in the fields. With the feedback provided by the entomologists, the model can be revisited and better explainable interfaces can be built. 2. All outputs of previous step concern directly the farmers and decision makers. Interfaces are necessary because they allow individual tracking of each agent in the system. This individual follow-up would better explain emerging behaviors, which is not the case for macroscopic modeling. Interfaces are useful to forecast the future according to parameters that they would have entered previously. Farmers could also use those outputs to have an optimal plan for the different interventions (insecticide spray) on the field. Also the entomologist could need this outputs to have a first idea of the results that they could have for their experimentations. The ultimate goal for this work will be to explain to users the behavior of an agent-model and simulation of the mirids in a cocoa farm. This explanation will help the users understand the behavior of the mirids, and hence, help in their decision making process to prevent the damage caused by mirids to the cocoa pods in cocoa farms.

Conclusion

In this paper, an initial agent-based model representing the population dynamics of mirids interacting on a cocoa pod is presented. The model is built based on the literature and by different steps of GAIA methodology. According to our knowledge of the literature, our work is the first agent-based model that is developed to study mirids pest. This approach will greatly improve In terms of integrating pest management strategy, our model aims to create a comprehensively designspecific, detailed, and integrated system to control specific pest species within specific crop zones. With enough refinement and insight, the behavior of mirids will possibly be evaluated at a specific zone. Next areas of our work could be : (i) the simulation and validation of our model ; (ii) the evaluation of damage caused by the mirids ; (iii) the integration of the impact of human actions (spread of insecticide) on the mirids population ; (iv) how to auto-regulate the population of mirids by the use of ants for biological control of insects ; (v) testing our model and simulation with entomological research data.

Figure 1 -

 1 Figure1-Mirid on a cocoa pod[START_REF] Babin | Rearing method and life-table data for the cocoa mirid bug Sahlbergella singularis Haglund (Hemiptera : Miridae)[END_REF]

Figure 3 -

 3 Figure 3 -State-Transition diagram of the role MIRID

Figure 4 -Figure 5 -

 45 Figure 4 -The lay eggs interaction protocol

Figure 6 -

 6 Figure 6 -Agent model

Figure 7 -

 7 Figure 7 -Workflow of our future work

Table 1 -

 1 The environment model

	Action	Environment Abstraction	Description
			The cocoa pod
	Reads	Agents position	identifies position occupied by the
			mirid
			Part of the surface
	Changes	Surface destroyed	of the pod is des-troyed by the mirid
		Position updated The mirid moved

Table 2 -

 2 Global variables of the environment

	Variables	Description
	Relative humidity	The climate gives the current relative humidity of the system
	Temperature	The climate gives the current temperature of the system
		The climate gives the current
	Shade	shade received by a part of the
		system

Table 3 -

 3 The MIRID functional role schema

	Role name : MIRID
	Description : This role is in charge of representing
	the birth, growth, mating, feeding, and death of a
	cocoa mirid
	Protocol and activities :
	move, feed, grow, be born, mate, lay egg, die
	Permissions :
	Reads : PosX, PosY
	Changes : Occupied_area,	an area of cocoa
	pod occupied by the mirid
	Destroyed_area, an area destroyed by
	the feeding of mirids
	Responsibilities :
	Liveness : M IRID = (to_be_born.grow +
	.move.feed.[mate].[ lay_egg]).die) n
	Safety : numberMirid > 2,
	miridMale.PosX=MiridFemale.PosX,
	miridMale.PosY=MiridFemale.PosY,
	N ew_Destination ∈ set_points_of _cocoa_pod,
	Content(N ew_Destination) = ∅
	cocoa pod. Table 2 presents the description of these
	global variables.

Table 4 -

 4 Organizational rules

	Organizational Rules Description
	Liveness rules

Table 5 -

 5 The service model

	Services Inputs		Pre-	Post-
				conditions	conditions
	Initiali-	Pod_area	Mirids num-
	sation	and	Mi-	ber ≥ 2 and
		rids		at least one
		number	female and
				one male
	Moving Destination	Destination	Occupied posi-
		position	position	is	tion
				free
	Death	Actual po-	Mirids num-	Mirids number
		sition		ber > 2	decreased
	Feeding Feed area	∃	eatable	Destroyed
				area	area

This term has been introduced by the French biologist Pierre-Paul Grassé in 1959 for describing the termite behavior. It is defined as : "Stimulation of the workers by the work that they performed." This term expresses the notion that the ac-

Acknowledgments

This work is supported by AUF(Agence Universitaire de la Francophonie), and partially by the Regional Council of Bourgogne Franche-Comté (RBFC, France) within the project UrbanFly 20174-06234/06242.


Giacomo Di Benedetto

Riccardo Scorretti

Arnaud Bréard

Damien Voyer

Fabien Ndagijimana

Kassem Jomaa

Melina Bouldi

Giacomo Di Benedetto

Arnaud Breard

Exposure of workers to the electromagnetic field of a RF welding machine

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

Introduction

INRS estimates around 100 000 radio-frequency (RF) industrial machines in France, which generate a high electromagnetic field [START_REF] Demaret | Intégrer le risque «rayonnements électromagnétiques» dans le document unique d'évaluation des risques professionnels[END_REF]. For this reason, the main objective of EXTI 1 (PNREST Anses, 2015/2 RF/02) is to develop an effective tool for easily obtaining relevant information on workers exposure to RF electromagnetic fields and, if necessary, to propose appropriate measures to reduce exposure. In the present work, the effect of human body immersed in an RF electromagnetic field was studied by simulating a preliminary virtual phantom built with concentric spherical layers of different materials.

Materials and methods

A shielded plastic welding machine (3 kW RF output power) was used as test device. The distance between the machine and the nearest wall was greater than 1m. The magnetic field was measured in different points in front of the electrode at a sampling frequency of 250 MHz. With these preliminary measurements we have characterized the frequency spectrum of the RF field, which has been used for numerical simulations.

A RF model implemented for finite element simulations in HFSS was developed. The model was built starting from the CAD of the plastic welding machine electrodes. For simplicity, in these preliminary simulations a multilayer sphere (Table I) has been used instead of a realistic anthropomorphic computational phantom. For living tissues the dielectric properties in [2] were used. The source is simulated as a lumped port between the electrodes with a power of 3 kW and at 2 different frequencies (27.12 MHz and 108.48 MHz). For each phantom and frequency the electric field peak value (in V/m) is shown in Table I, while for the most interesting simulations the electric field vectors are shown in Figure 1.

Results and conclusion

Some of the simulated configurations are depicted in Fig. 1. It is found that the presence of the phantom (#2,4,5) modifies to a large extent the electric field, with respect of the case where the phantom is absent (#1). Moreover, in presence of any phantom it is observed that the electric field is nearly perpendicular to the external surface of the phantom, which can be easily explained by the high permittivity and conductivity of the phantom with respect of air (Table I). In the simulated frequency range (27.12 -108.48 MHz) the peak value of the electric field exhibits a low dependence (about 10%) on the dielectric parameters of the material which compose the phantom, as long as they have the same order of magnitude of biological tissues. In particular, the internal structure of the phantom seems to play a quite minor role. The main factor which explain the observed differences between the phantoms #2,3 and #4,5,6 is the high relative permittivity of the skin, compared to the one of PVC.

On the base of these preliminary simulations, we aim for building an instrumented phantom composed of PVC filled (at least partially) with a tissue-mimic liquid. In fact, it is well known that the human body modifies significantly the electric field distribution. As our purpose is to check the compliance with basic restrictions, the electric field must be measured in presence of a phantom which mimics the human body, in order to obtain realistic data to be used for numerical dosimetry. Probes for measuring the normal component n⋅E of the electric field should be placed on the surface of the phantom. The knowledge of n⋅E, measured in realistic situations, should allow to numerically compute the electric field inside the body, so as to characterize the human exposure of workers to the RF field (more detailed simulations will be shown in the full paper).

  , ɛ = 166, σ = 0.33), blood (r = 120, ɛ = 127, σ = 1.16) 492.92 257.20 3 Skin (r = 130, ɛ = 166, σ = 0.33), blood (r = 120, ɛ = 127, σ = 1.16), brain (r = 60, ɛ = 221, σ = 0.57) r = 130, ɛ = 2.7, σ = 0), blood (r = 120, ɛ = 127, σ = 1.16) 437.68 232.54 5 PVC (r = 130, ɛ = 2.7, σ = 0), blood (r = 120, ɛ = 127, σ = 1.16), PVC (r = 90), vacuum (r = 80, ɛ = 1, σ = 0) r = 130, ɛ = 2.7, σ = 0), blood (r = 120, ɛ = 127, σ = 1.16), PVC (r = 60), vacuum (r = 50, ɛ = 1, σ = 0)

Table I -

 I Electric field peak value for different modelling configurations. For each material the radius (r) in mm, the

	relative permittivity (ɛ) and the conductivity (σ) in S/m at 27 MHz are shown in brackets
	# Model	E (V/m) at	E (V/m) at
		27.12 MHz	108.48 MHz
	1 No phantom	235.97	109.75
	2 Skin (r = 130

EXTI : EXposition des Travailleurs aux champs électromagnétiques Industriels.


J P N Watson

A Gardeisen

Confirmation of rabbit domestication in the 4th century AD at Ambrussum: evidence of morphological change in the pelvis

come

During the excavation of a well, Puits 103, at the site of Ambrussum, near Villetelle (Hérault, France) (see map in Figure 1), numerous rabbit bones were found in a context that suggested that complete rabbit carcasses had been thrown into the well within a short space of time [START_REF] Fiches | Quatre puits de l'agglomération routière galloromaine d'Ambrussum (Villetelle, Hérault)[END_REF][START_REF] Gardeisen | Des clapiers à Ambrussum au IVe s. de notre ère? Un nouvel indice d'élevage de lapins dans le Midi de la France[END_REF]. The minimum number of individuals was 35. All parts of the skeleton were well represented apart from the very small bones such as phalanges and there was no sign of butchery or burning.

A broad spread of ages was represented, suggesting that the bones accumulated from the successive throwing of animals into the well as they died, rather than the result of a sudden disease striking the entire population. This even spread of ages suggested that the rabbits were being kept in captivity (Gardeisen, 2012, 301).

To test whether there was any sign of morphological change that might be the result of domestication, whether by deliberate selection or not, six measurements were taken on the coxal bone of the pelvis (os coxae) in those specimens that had the acetabulum fused or in the process of fusion. The minimum number of individuals with the acetabulum fully fused was 13 and the minimum number of individuals with the acetabulum still in the process of fusion was 6. The measurements are described in Appendix I and Figures 2 and3. These measurements had previously been taken on modern wild rabbits from various locations in France and the Iberian Peninsula as well as on rabbits from three archaeological sites in the Peninsula, and unexpected differences had been found between rabbits of different geographical origins [START_REF] Watson | Shape differences in the pelvis of the rabbit, Oryctolagus cuniculus (L.), and their genetic associations[END_REF].

When the measurements were taken on the Ambrussum specimens it was immediately seen that the shape of the ischium is very different in the larger Ambrussum specimens and in fact the minimum diameter of the shaft, measurement "I", is clearly bimodal in the adults (Figure 4). At first it was thought that this bimodality might represent a difference between the sexes, but the difference is much greater than that seen between the sexes in the specimens from Roissy, the modern sample in which the sexual dimorphism is greatest (Figure 5).

In the study by Watson and Davis, the principal components of the measurements were found and subsequently a linear discriminant function based on the 2nd to 6th principal components was calculated. The first principal component was then plotted against the discriminant score. This is an established technique for separating the variation in the measurements into a sizecomponent, represented by the first principal component, and a shapecomponent, represented by the discriminant score (see e.g. Burke et al., 2003). Figure 6 shows the mean first principal component plotted against the mean discriminant score for the five largest samples of modern wild rabbits from the subspecies Oryctolagus cuniculus cuniculus and the rabbits from three Late Palaeolithic phases at Cova Matutano [START_REF] Olària | Cova Matutano (Vilafamés, Plana Alta, Castellón). Un modelo ocupacional del magdaleniense superiorfinal en la vertiente mediterránea peninsular[END_REF][START_REF] Watson | Análisis Preliminar de los Mamíferos del Sector 3 de Cova Matutano (Vilafamés, Castellón)[END_REF], using Roissy as one of the groups for the discriminant function and Matutano A as the other.

Figure 7 shows the individual specimens from Ambrussum compared with the modern wild rabbits as a whole. It can be seen that the specimens from Ambrussum with a thick ischial shaft (i.e. a large "I") were both larger overall and more "Iberian" in shape than the mass of modern wild rabbits, but so were two specimens with a slender ischial shaft. This group will be referred to as the "Large" group from now on, while the specimens that fall amongst the modern wild rabbits will be referred to as the "Small" group. The very large wild specimen to the left of the main part of the "Large" group is in fact an extreme outlier from Roissy in France; the measurements have been retaken and are correct.

Figure 8 shows the individual specimens from Ambrussum compared with the Magdalenian ones from Matutano. Although the "Large" Ambrussum group overlaps with the largest specimens from Matutano, it is on the edge of the mass of Matutano rabbits as a whole. Figure 9 shows the result of adding the means for the two Ambrussum groups to Figure 6. It shows that although the "Small" group fits in well with the wild rabbits the mean for the "Large" group plots out a long way above those for Cova Matutano.

In accordance with Bergmann's rule (Bergmann, 1847;[START_REF] Davis | Of metapodials, measurements and music -eight years of miscellaneous zooarchaeological discoveries at the IPA, Lisbon. O Arqueólogo Português[END_REF], which states that mammals from colder environments tend to be larger, one explanation for the large size of the "Large" Ambrussum rabbits might be an origin in a cold part of the Peninsula, which would not be impossible, given that Ambrussum lies on the main line of communication between the Peninsula and Italy, the Via Domitia.

The other explanation would be that the difference in shape is the result of domestication or the processes leading to it. This is supported by the observation that measurement "I" is disproportionately large in all but four of the "Large" Ambrussum specimens compared with measurement "H" when the Ambrussum specimens are compared with the modern wild ones (Figure 10). Although this difference is much less marked when the Ambrussum specimens are compared with the Magdalenian ones from Matutano, it is still noticeable (Figure 11).

It should be stressed that both measurements "H" and "I" are very dependent on age up to the point at which the acetabulum is fully fused. After this, growth appears to slow sharply. Figure 12 shows the Ambrussum specimens classified according to their stage of development [START_REF] Gardeisen | Des clapiers à Ambrussum au IVe s. de notre ère? Un nouvel indice d'élevage de lapins dans le Midi de la France[END_REF]. There is little difference between Gardeisen's Stades 4 and 5 but the Stade 3 specimens are markedly smaller. For this reason only adult specimens, i.e. those in Stades 4 and 5, are included in the study, but even so the possibility that differences in "H" and "I", especially "H", might correspond to differences in age needs to be borne in mind. However, it is clear that the Stade 3 specimen with a value of 37 for "I" (3.7 mm) must belong with the "Large" group, not only because of its position in Figure 12 but also because of its behaviour in the multivariate analyses.

Nevertheless, given that measurement "I" is so large in the "Large" Ambrussum specimens, it is difficult to escape the conclusion that the differences between them and all the wild specimens studied, both modern and ancient, are the result of domestication and subsequent genetic selection, whether deliberate or not.

What is surprising is the direction of the change in shape. Figure 13 shows the Ambrussum specimens together with the modern wild ones, as in Figure 10, but with the addition of specimens from two domestic breeds, the Gris de Vienne and the Papillon [START_REF] Watson | Shape differences in the pelvis of the rabbit, Oryctolagus cuniculus (L.), and their genetic associations[END_REF]. These modern breeds are very much larger and so the scale of the figure has had to be altered. The striking feature of the figure is that there is a disproportionate increase in measurement "H" in the modern domestic breeds, which is exactly the opposite to what is observed for the "Large" Ambrussum specimens.

This difference between the "Large" Ambrussum specimens and the modern domestic breeds is also seen in the plot of the the mean first principal component against the mean discriminant score. The two domestic breeds plot out on the opposite side to the Ambrussum specimens and have mean discriminant scores similar to those for the French wild specimens (Figure 14). This is partly the consequence of the differences in the relationship between measurements "I" and "H" shown in Figure 10, but it is also influenced by differences in the ilium.

That the modern domestic rabbit breeds resemble the French wild rabbits in shape is not surprising given that they are thought to be descended from French wild stock and have been shown to resemble modern French wild rabbits genetically [START_REF] Carneiro | The Genetic Structure of Domestic Rabbits[END_REF]. Studies of the mitochondrial DNA show that the B1 type, widespread in France, predominates in all the main domestic breeds and that no exclusively Iberian types occur in them [START_REF] Queney | Different Levels of Human Intervention in Domestic Rabbits: Effects on Genetic Diversity[END_REF]. But what then might be the genetic affinities of the "Large" Ambrussum rabbits?

This unexpected discovery of a presumably domestic strain with a morphology different from presentday domestic breeds raises the question of whether there could have been two separate episodes of domestication (in the sense of morphological change), the welldocumented one in the 11th12th centuries AD leading to the spread of warrens ("garennes") across France (Callou, 2003, 255256) and another one in the 4th century AD or earlier. In the intervening centuries the strain or strains developed could have died out or could simply have been absorbed into the much larger later expansion of domesticated rabbits.

The specimens in the "Small" group from Ambrussum also raise many questions. Although they are clearly differentiated from the "Large" group on the shape of the ischium, there are also important similarities, namely their rather "Iberian" shape and the large overall size of some of them. It needs to be asked whether they might not represent the population from which the "Large" group had been bred or had mutated. The smaller ones are indistinguishable from modern samples from Navarre but very different from the modern sample from Tour du Valat, which is less than 50 kilometres away from Ambrussum on the far side of the Camargue. Unfortunately there is no modern material available from the Hérault apart from a single wild rabbit from Lattes, which falls into the middle of the "Small" Ambrussum group. It could be that the modern wild rabbits are different from those in the Camargue and closer in shape to those from Navarre and that the wild rabbits of the 4th century AD were similar.

The "Large" coxal bones from Ambrussum, although they form a striking cluster, include both left and right bones. Some of these, if not all, are paired bones from the same individual, so that the actual number of individuals represented is only between 5 and 8. Since the rabbits all appear to have been deposited in the well within a short space of time, it is quite possible that all these individuals were closely related. It might be thought, therefore, that these bones were just the result of a chance mutation that went no further, were it not for the fact that a very similar bone occurs in a 4th century AD level at the nearby site of Lunel Viel (Raynaud, 1990). The value of measurement "I" is 39, which puts it into the group with a thick ischial shaft (see Figure 4). Unfortunately measurement "F" is slightly damaged, but the damage is not enough to affect the measurement by more than three tenths of a millimetre, so upper and lower limits for the value have been used and the first principal component and discriminant function calculated in each case. The result is shown in Figure 15 and it can be seen that the specimen clearly falls into the "Large" group. The site at Lunel Viel is only 7 km away from Ambrussum, so there is still a faint possibility that the animal was closely related to the "Large" Ambrussum rabbits, but it seems more likely that it represents a general phenomenon and that large rabbits were being deliberately bred in the eastern Languedoc in the 4th century AD. Figure 9. The same as Figure 6 but with the addition of the means for the "Large" and "Small" specimens from Ambrussum. It can be seen that the "Large" specimens not only are considerably larger than the wild rabbits, both modern and Magdalenian, but also have a shape that is at the opposite extreme to that of the French wild rabbits from Roissy, Versailles and Tour du Valat.

Figure 10. Measurement "I" plotted against measurement "H" for the modern wild specimens (small grey dots), the "Small" Ambrussum specimens (grey circles) and the "Large" Ambrussum specimens (black filled circles and black open circles). It can be seen that the "Small" specimens are indistinguishable from the wild ones but that most of the "Large" specimens stand apart from both, if the outlier from Roissy is excluded, with the exception of the two "Large" specimens with a slender ischial shaft. There are fewer "Small" Ambrussum specimens shown than in Figures 6 and7 because in four cases circles are superimposed.

Figure 11. Measurement "I" plotted against measurement "H" for the Magdalenian specimens from Cova Matutano (grey squares), the "Small" Ambrussum specimens (grey circles) and the "Large" Ambrussum specimens (black filled circles and black open circles). It can be seen that the "Small" specimens are indistinguishable from the Matutano ones but that in the "Large" specimens measurement "I" is greater in the majority of cases although there is considerable overlap, Figure 12. Measurement "I" plotted against measurement "H" for the three growth stages recorded in the Ambrussum specimens [START_REF] Gardeisen | Des clapiers à Ambrussum au IVe s. de notre ère? Un nouvel indice d'élevage de lapins dans le Midi de la France[END_REF]. The triangles indicate specimens in Stade 3 (juveniles), the circles those in Stade 4 (subadults) and Stade 5 (adults). The "Small" specimens are in grey and the "Large" ones in black.

Figure 13. Measurement "I" plotted against measurement "H" for the modern wild specimens (small grey dots), the "Small" Ambrussum specimens (grey circles), the "Large" Ambrussum specimens (black circles), specimens of the Gris de Vienne breed (black triangles) and specimens of the Papillon breed (grey triangles). Figures 1 and2 show the approximate positions in which the measurements are taken. All the measurements are taken with a vernier gauge with "knifeedge" blades (i.e. they should have reasonably narrow edges, although that is only really important for F1).

The formal definitions are as follows: H = minimax diameter of ischial shaft, between the acetabulum and the spina ischiadica. I = absolute minimum diameter of ischial shaft (usually in the same place as H, but sometimes much more caudal). J = minimum height of the lateral side of the acetabulum, i.e. the minimax. This is only just possible in lagomorphs, one of the blades of the vernier gauge lying almost flat against the back of the acetabulum (i.e. the medial side). It also needs a certain delicacy of touch because the edge of the acetabulum is easily damaged. In addition the acetabulum needs to be fully fused. E = absolute minimum diameter of the shaft of the ilium. This is very oblique in rabbits, at about 45 degrees to the axis of the ilium. It should not enter the area of fusion with the sacrum (which it very occasionally tends to do). F = minimax diameter of the shaft of the ilium, usually the absolute minimax but sometimes a local minimax. The measurement is near to or equivalent to the dorsoventral height in most cases, but taken as a minimax. However, in some cases, when the spina iliaca dorsalis caudalis extends a long way caudally, this is only a local minimax and the absolute minimax is further round, more or less parallel to F1. Sometimes it is difficult to take because the depression on the ventral side of the spina iliaca ventralis caudalis is very shallow and the bone then has to be placed right between the tips of the blades of the vernier gauge in order to find the minimax. F1 = minimax diameter of the shaft of the ilium with one blade of the vernier gauge placed between the edge of the acetabulum and the spina iliaca ventralis caudalis. One needs to be careful to avoid damaging the spina.

For detailed drawings and nomenclature see Barone et al. (1973, 24).

Ambrussum:

No.

Side

  ., Pavaux, C., Blin, P.C. & Cuq, P. (1973). Atlas d'anatomie du lapin. Paris: Masson. Bergmann, C. (1847). Ueber die Verhältnisse der Wärmeökonomie der Thiere zu ihrer Grösse. Göttinger Studien, Abt 1., 595708. Burke, A., Eisenmann, V. & Ambler, G.K. (2003). The systematic position of Equus hydruntinus, an extinct species of Pleistocene equid. Quaternary Research 59, 459469. Callou, C. (2003). De la garenne au clapier: étude archéozoologique du lapin en Europe occidentale. Mémoires du Muséum national d'Histoire naturelle, 189. Paris: MNHN.

Figure 1 .

 1 Figure 1. Map showing the location of the site of Ambrussum.

Figure 2 .

 2 Figure 2. Drawing showing the approximate positions in which measurements H, J and F1 are taken. (N.B. H, J and F1 are not necessarily in parallel planes; they just happen to have that property in the specimen illustrated.)

Figure 4 .

 4 Figure 4. Histogram of measurement "I" for the adult specimens from Ambrussum, showing the marked bimodality. The horizontal axis is in tenths of a millimetre. The specimens classified as "adult" are those belonging to Stades 4 and 5 as defined in Gardeisen & Valenzuela Lamas (2004, Figure 9).

Figure 5 .

 5 Figure 5. Histograms of measurement "I" for the modern specimens from Roissy. The upper histogram shows the male specimens and the lower histogram the female ones.

Figure 6 .

 6 Figure 6. Plot of the mean first principal component against the mean discriminant score for various samples of modern rabbits and for three Magdalenian faunal phases at Cova Matutano. Only adult specimens are included, i.e. those in Stades 4 and 5 as defined by Gardeisen & Valenzuela Lamas (2004).

Figure 7 .

 7 Figure 7. The individual values for the first principal component plotted against the individual discriminant scores for the adult Ambrussum specimens and for modern wild rabbits from Roissy, Versailles, Tour du Valat, Peralta and Tudela, i.e. those belonging to the subspecies Oryctolagus cuniculus cuniculus. The black filled circles are the Ambrussum specimens with a thick ischial shaft. The black open circles and the grey filled circles are the Ambrussum specimens with a slender ischial shaft. The small grey dots are the modern wild rabbits.

Figure 8 .

 8 Figure 8. The individual values for the first principal component plotted against the individual discriminant scores for the adult Ambrussum specimens and for the Magdalenian specimens from the faunal phases A, B and C at Cova Matutano. The black filled circles and the black open circles are the "Large" Ambrussum specimens. The grey filled circles are the "Small" Ambrussum specimens. The paler grey squares are the Cova Matutano specimens.

Figure 14 .

 14 Figure14. The mean first principal component plotted against the mean discriminant score, as in Figure9, with the addition of the specimens of the Gris de Vienne and Papillon breeds.

Figure 15 .

 15 Figure 15. The same as Figure 7 but with the specimen from Lunel Viel added (black triangles). The two triangles close together joined by a line are the maximum and minimum estimates.

Appendix 2: List of measurements.

Column 1 shows JPNW's reference number. Column 2 shows the side; l = left; r = right. Column 3 shows the agegroup according to AG's Stades. Column 4 shows the group the specimen has been assigned to : L = "Large"; S = "Small". Columns 5 to 10 show the measurements in millimetres.

To calculate the first principal component and the discriminant score for a specimen directly from the raw measurements in millimetres, the following formulae can be used: PC1 = 0.753F1 + 1.150F + 1.080E + 0.974J + 1.022H + 1.327I 37.6209 Score = 95.75F1 + 25. 65F + 75.64E 36.61J + 94.27H 64.97I + 155.109


C Langlois

T Douillard

H Yuan

N P Blanchard

A Descamps-Mandine

B Van De Moortèle

C Rigotti

T Epicier

Crystal Orientation Mapping via ion channeling: an alternative to EBSD

Keywords: electron backscattered diffraction, orientation mapping, ion imaging, focused ion beam, texture, channeling

A new method, which we name iCHORD (ion CHanneling ORientation Determination), is proposed to obtain orientation maps on polycrystals via ion channeling. The iChord method exploits the dependence between grain orientation and ion beam induced secondary electron image contrast. At each position of the region of interest, intensity profiles are obtained from a series of images acquired with different orientations with respect to the ion beam. The profiles are then compared to a database of theoretical profiles of known orientation. The Euler triplet associated to the most similar theoretical profile gives the orientation at that position. The proof-of-concept is obtained on a titanium nitride sample. The potentialities of iCHORD as an alternative to EBSD are then discussed.

Introduction

Electron backscatter diffraction (EBSD) is routinely employed as a characterization tool coupled with a scanning electron microscope (SEM) to obtain individual grain orientations, local texture, point-to-point orientation correlations, and phase identification and distributions on the surfaces of bulk polycrystals [START_REF] Schwartz | Electron backscatter diffraction in materials science[END_REF]. As commercial EBSD systems are now readily and widely available, the technique is no longer reserved to specialized research laboratories. Consequently EBSD has become a tool of choice for process development and quality control. With the goal of full automation in mind, efforts are being made to increase the speed of the technique [START_REF] Yang | Adaptive orientation imaging microscopy[END_REF] to obtain orientation maps.

Several avenues are currently being explored, such as increasing the sensitivity and speed of the digital camera [START_REF] Maitland | Advances in High-Speed EBSD Orientation Mapping[END_REF] as proposed by the manufacturers, as well as improving the indexation algorithms for treating and indexing patterns [START_REF] Britton | Measurement of residual elastic strain and lattice rotations with high resolution electron backscatter diffraction[END_REF]. In this context, any new technique that can ease the implementation of orientation mapping at an industrial scale would be welcome. The precision of the indexation is also the subject of intense research efforts. Indexing algorithms tend to be more and more precise, revisiting the classical Hough transform used to extract the diffraction lines from the Kikuchi patterns [START_REF] Tao | Errors, artifacts, and improvements in EBSD processing and mapping[END_REF] [START_REF] Maurice | A 3D Hough transform for indexing EBSD and Kossel patterns[END_REF]. Other approaches are currently explored, using for instance the pattern matching technique, which consist in comparing directly the Kikuchi patterns to a database of simulated ones. This idea was proposed originally by Wright et al. in 1991 [7] and was recently implemented again by Chen et al. with the benefits of modern computational power [START_REF] Chen | A Dictionary Approach to EBSD Indexing[END_REF]. The same technique has been used in Transmission Electron Microscopy to obtain orientation maps using precession diffraction [START_REF] Rauch | Automatic crystal orientation and phase mapping in TEM by precession diffraction[END_REF][START_REF] Rauch | Rapid spot diffraction patterns idendification through template matching[END_REF].

We propose in this paper a new technique for orientation mapping based on the well-known channeling contrast phenomenon [START_REF] Reimer | image formation and microanalysis[END_REF] [START_REF] Zaefferer | Theory and application of electron channelling contrast imaging under controlled diffraction conditions[END_REF]. As the channeling contrast is emphasized in ionic images [START_REF] Phaneuf | Applications of focused ion beam microscopy to materials science specimens[END_REF], the experiments were conducted in a FIB-SEM instrument to acquire ionic image series as well as comparative EBSD orientation maps. According to the channeling theory for ions and electrons, the intensity of induced secondary electrons (SE) at a given position on the surface of a polycrystalline sample will change if the orientation of the primary incident beam with respect to the sample surface is changed [START_REF] Yahiro | Crystallographic orientation contrast associated with Ga+ ion channelling for Fe and Cu in focused ion beam method[END_REF]. An intensity profile can then be drawn as a function of the sample orientation, from which information about the orientation of the crystalline lattice at this position can be retrieved. This idea has been suggested several times in the literature [START_REF] Yahiro | Crystallographic orientation contrast associated with Ga+ ion channelling for Fe and Cu in focused ion beam method[END_REF][START_REF] Silk | Determination of lattice orientation in aluminium alloy grains by low energy gallium ion-channelling[END_REF][START_REF] Veligura | Channeling in helium ion microscopy: Mapping of crystal orientation[END_REF], but, to our knowledge, no comprehensive method currently exists that is capable of yielding full orientation maps on common polycrystalline samples, such as pure metals or metallic alloys. A relevant contribution was recently brought by Veligura et al. [START_REF] Veligura | Channeling in helium ion microscopy: Mapping of crystal orientation[END_REF] but orientation maps were obtained only for particles deposited in epitaxy on a substrate, with only a single degree of freedom in rotation for the lattice orientation of the different grains.

This paper is divided into several sections. In Section 2, the approach we used to obtain ionic image series is outlined, with a comparison between a tilting image series and a rotation image series. Section 3 deals with the experimental details and Section 4 addresses the question of the different possible ways to recover the crystallographic orientation associated with an experimental intensity profile. The results and performances obtained with our method on a TiN sample are exposed in Section 5. The following discussion (Section 6) addresses several important points such as the spatial resolution, the time to obtain an orientation map and the applicability of this technique to non-cubic materials as well as multi-phased materials. We conclude this study in Section 7, with an outline of our key results.

Intensity profiles as a signature of the orientation

As mentioned in the introduction, the channeling effect has been theoretically modeled both for ions and electrons [START_REF] Giannuzzi | Comparison of channeling contrast between ion and electron images[END_REF]. These models include a significant role of lattice plane orientation with respect to the charged particle beam and channeling contrast has frequently been applied to enhance microstructural details, especially in SEM [START_REF] Levi-Setti | Ion channeling effects in scanning ion microscopy with a 60 keV Ga< sup>+</sup> probe[END_REF] [19] [13] [15, 16] [20]. The work of Yahiro et al. clearly illustrates the variation of the SE signal under a Ga + probe, with minima appearing when the beam arrives parallel to low index crystallographic planes [START_REF] Yahiro | Crystallographic orientation contrast associated with Ga+ ion channelling for Fe and Cu in focused ion beam method[END_REF]. However, if the aim of an experiment is to recover the crystallographic orientation of the sample, acquiring a tilt series of images is clearly not the best choice for symmetry reasons. Usually, the tilt axis lies in the sample surface, let say along the e x direction as represented in the drawing of Figure 1 a and b which represent respectively the reference frame, common to the sample and the stage, as well as its position in a stereographic projection. Within a tilt series, the ion or electron beam originally perpendicular to the sample surface along e z direction crosses the stereographic projection linked to the reference frame along a plane passing through the center of the projection, containing e y and e z directions. If we assume that the signal received only depends on the intersections between the beam and the crystallographic planes or directions, the same signal evolution will be obtained for two different crystals with orientations linked by mirror symmetry with respect to the geometric plane in which the beam travels. An example is given in Fig. 1c, showing two stereographic projections corresponding to such mirror orientations, with the beam path represented by a vertical dashed line. A beam tilted around the e x direction will cross the same planes at exactly the same angular positions on the two stereographic projections, giving the same intensity variations. This kind of situation does not allow an unambiguous determination of the crystallographic orientation from the obtained intensity profile and additional experiments would then be required. On the contrary, if we give a fixed tilt to the sample, for arguments sake 40°, and then rotate the sample around its tilted normal (e z direction), the beam will circulate along a cone figured by the dotted circle in the stereographic projections of a given grain (see Figure 1c). A series of image acquired during such a rotation of the sample (around a rotation axis not parallel to the beam) then results in unique intensity profiles for a given orientation because there is no other stereographic projection that would result in intersections with the same planes at the same angular positions, with the same sense of rotation. In this work, we used rotation image series to obtain the orientation of the crystal for each point of the region of interest (ROI).

Concerning the initial tilt, because of the necessity to cross enough planes to obtain a significant signature, we have chosen rather arbitrarily 40°. It is worth noting that a correction of beam deflection during the image scan has to be included if working at low magnifications [START_REF] Giannuzzi | Comparison of channeling contrast between ion and electron images[END_REF]. In the present work, the magnification is 500X, with an horizontal field of view of 230 µm. In these conditions, the beam deflection is low enough to ignore the effect.

Experimental setup

Experiments were carried out on a ZEISS NVision40 FIB-SEM machine in order to obtain both ioninduced electron images and reference orientation maps, acquired with an Oxford Instruments F +TM EBSD system (fast camera NordlysII and Channel 5 software). Concerning the ionic source, the microscope is equipped with a Seiko Zeta column and a Ga liquid metal ion source. FIB images were acquired at 30 kV and a current of 3nA. This quite high current was chosen to have a probe size of 170 nm that matched the pixel size in the images. Matching the probe size and the image pixel size allows milling the surface evenly, without localized craters during the acquisition of an image series [START_REF] Orloff | Fundamental limits to imaging resolution for focused ion beams[END_REF]. Moreover, the dwell time per pixel was 0.1 µs, which limits milling, amorphisation and Ga implantation. With these conditions, the ion channeling does not fade out, even after completing the whole series of 360 images. A titanium nitride (TiN) sample was used to acquire ion image series and EBSD maps. Titanium nitride has a face-centered cubic structure with a TiN motif localized at the origin of the crystallographic cell, and in the middle of each face. The sample does not present any strain or any precipitate in its microstructure. The surface of the sample was prepared as for classical EBSD: mechanical polishing, and vibratory polishing using colloidal silica as the final preparation step.

Computational approach

Acquisition and processing of image stacks

The raw images acquired in the experimental conditions exposed in Section 3 were numerically treated within the Fiji software [START_REF] Schindelin | Fiji: an open-source platform for biological-image analysis[END_REF] in order to get an aligned stack of images. Firstly, all the images being acquired with a tilt of 40°, they are compressed in the direction perpendicular to the tilt axis and have to be rescaled. Then, as we chose to make a rotation series, each image must be corrected for its own rotation angle. Finally, residual drifts in X-Y directions were corrected [START_REF] Thevenaz | A pyramid approach to subpixel registration based on intensity, Image Processing[END_REF] and the largest common area in all the images of the series was finally cropped. The final result is a stack of images reflecting the same area, with identical dimensions, but showing variations of contrast in each pixel due to the channeling phenomenon. At any position (x, y) in the image, we obtain an intensity profile by representing the intensity variation along the depth of the stack, the number of the slice being linked to the rotation angle. Figure 2 illustrates the processing from aligned images to intensity profiles.

Linking intensity profiles and crystallographic orientation

Each intensity profile that can be extracted on each pixel of the ROI is a signature of the crystallographic orientation at that position. Several approaches are possible to recover the orientation from the intensity profiles.

The first one would be to analyze a given profile and to extract the significant information, like the plane family associated to each peak, as well as the angular separation between the peaks. Knowing the experimental tilt of the sample (here 40°) during the acquisition, this could be enough to determine geometrically the entire stereographic projection, and hence to come back to the three Euler angles. However, we observed that the peaks sometime correspond to the contribution of several planes, when the beam comes parallel or near a zone axis. Therefore, the value of an intensity minimum cannot be exclusively assigned to a given plane crossed by the beam. The second approach is a "brute force" technique which consists in constituting a theoretical database of the profiles corresponding to all the possible orientations for the material. Then, a search is carried out in the database to find the closest theoretical profile as compared to the experimental one, according to a given metric. Since the Euler angle triplet of this theoretical profile is known by construction, the orientation will be determined. This approach has been successfully used for determining orientations in a TEM by matching experimental precession diffraction patterns to a database of theoretical patterns [START_REF] Rauch | Automatic crystal orientation and phase mapping in TEM by precession diffraction[END_REF][START_REF] Rauch | Rapid spot diffraction patterns idendification through template matching[END_REF]. In the present case, the key point is to construct the theoretical database.

This theoretical database must contain the intensity profiles associated to orientations represented by Euler triplets. To construct an intensity profile for a given orientation, one needs to know i) all the planes that will be crossed or closely approached by the beam ii) the shape of the peaks, iii) the loss amplitude attributed for each crystallographic plane, and iv) an adjustable factor to scale the calculated profiles to the intensity range of the experimental images. In our experimental setup, the planes giving a contribution in the profile are simply the ones that have a normal direction elevated by an angle less than (90° -40° -ε) with respect to the surface of the sample. Here, the parameter ε is a critical angle describing the limit out of which the beam does not "see" the plane.

Normally, each family of planes has its own critical angle [START_REF] Giannuzzi | Comparison of channeling contrast between ion and electron images[END_REF], but for simplicity's sake, the same ε was taken for all the planes. A Gaussian curve was used to model the intensity loss observed when the beam crosses a crystallographic plane perpendicularly. Therefore, for each position of the beam (i.e. each abscissa position in the intensity profile) the angles between the beam and the planes considered are calculated, and the contributions of each plane to the intensity loss are added. It is worth noting that, with this method, when the beam crosses a plane perpendicularly, the loss peak is quite narrow whereas when the beam crosses a plane at a grazing incidence, the peak is wider. Moreover, due to this "critical angle", even if the beam does not geometrically cross a plane, it can be close enough to the plane (angle to the plane between 0 and ε ) to give a significant channeling.

By trial and error using the Graphical User Interface presented in Figure 3, the value for the Gaussian curve full width at half max has been fixed to 7°. Then, the intensity loss associated with each plane (i.e the depth of the Gaussian curve associated to a plane) must be determined. In this feasibility experiment, we can measure this parameter accurately by acquiring an EBSD map on the same area to construct the stereographic projection associated to a given experimental profile obtained in a given position. Knowing both the stereographic projection and the experimental profile, one can index the peaks in the profile and determine which intensity loss is associated to which family of lowindex planes. Figure 3 presents the Graphical User Interface (GUI) that has been coded to help determine the intensity loss assigned to each crystallographic plane from the EBSD data and experimental profiles. The values of intensity losses for the three main plane families are, in arbitrary units: 1000 for {100}, 550 for {110} and 358 for {111} (arbitrary units). A global scale factor has been fixed to 10 to account for the fact that the experimental data have not been normalized. It was found that considering only the {100}, {110} and {111} planes families was enough to obtain a good fit to the experimental profiles. All theoretical profiles can be stored in a database together with their corresponding Euler triplets. The link between an experimental profile, acquired during the sample rotation, and the crystallographic orientation can be established by searching in this database to find the theoretical profile that is the most similar to the experimental one and to retrieve the associated Euler triplet. As for the acquisition, performance is an issue since the database that has to be generated can be very large, especially if the step used for each Euler angle is small, and because in addition the search through the database has to be repeated for each position in the area of interest. Fast selection and access to stored data have been intensively investigated for several decades in the database research field. For instance, data structures as B trees and B+ trees [START_REF] Elmasri | Fundamentals of Database Systems sixth Edition Pearson Education[END_REF] have been designed to offer sub-linear access costs with respect to the size of the database. This means that, when using these kinds of data organization, if the size of a database is multiplied by K, then the time needed to retrieve an object according to a simple selection criterion will also increase, but less than by a factor of K.

In the context of intensity profiles, each object is a profile described as a sequence of N intensity values (N being the number of images in the image series) and the selection process is to find among a set of theoretical profiles the one that is the most similar to a given experimental profile. Each profile is a point in a space having N dimensions, and the problem can be reformulated as the search for the theoretical profile that is the closest to the experimental profile using a distance in this N dimensional space. So, all theoretical profiles are stored in a R-tree [START_REF] Guttman | R-trees: a dynamic index structure for spatial searching[END_REF], that is a data structure designed to support the fast retrieval of multi-dimensional information from databases. Informally, an R-tree is a tree representing a decomposition of the data space in hyper-rectangles that are bounding boxes of sets of points. At the root of the tree, is a hyper-rectangle H that encompasses the whole set, S, of the points stored in the database. Each child of the root is a node storing a hyperrectangle included in H, that is a bounding box of a subset of S. This refinement is repeated along the branches, and the points are stored at the leaf level. This structure is used to find efficiently the nearest neighbor of a given point in space [START_REF] Roussopoulos | Nearest neighbor queries[END_REF], with the Euclidean distance as metric.

To generate the theoretical profile database, the Euler angle steps were initially chosen as 3°, 2° and 2° respectively for  1 ,  and  2 . Owing to the cubic symmetry of the titanium nitride TiN lattice, the Euler angle  1 varies from 0° to 360°, and  and  2 vary from 0° to 90° [START_REF] Randle | Introduction to texture analysis: macrotexture, microtexture and orientation mapping[END_REF]. With this mapping of the Euler space, there is a maximum disorientation of 2.7 degrees between adjacent orientations in the database. This results in a database of 243,000 theoretical profiles.

To execute the database searches, the method takes direct advantage of the presence of several processor cores by distributing the set of experimental profiles over the cores, each process working with its own copy of the R-tree (containing the theoretical profiles). Since no communication is required between the processes, the execution time still scales linearly with respect to the number of experimental profiles, as if a single core were used. This is confirmed by Figure 4a that reports the whole search execution time using 30 cores when changing the number of experimental profiles to process (all profiles described in a space having 90 dimensions, simulating a case where the image series contains 90 images). The benefits of using a R-tree are underlined Figure 4b by the curve depicting the search time when varying the size of the database. This figure shows that the search time scales nicely in a sub-linear way with respect to the size of the theoretical profile databases (execution times reported for a single core, a set of 80 experimental profiles to process, using 90 dimensions per profile).

Results and performances

Full database search

We used the EBSD map acquired on the same area as a reference to evaluate the correctness of the iChord indexation. In the context of a 'proof-of-concept', we evaluate that a disorientation lower than 5° between the iChord orientation and the EBSD orientation at the same position can be considered as acceptable. To compare the results of the iChord method and the EBSD results (considered as a reference), the iChord and the EBSD data have been treated to have a common color code representing the Euler angles. The treatment, for a given orientation from EBSD or iChord maps, consist in choosing among the three symmetry-equivalent Euler triplets present in the Euler space (phi1 0°-360° ; PHI 0°-90° ; phi2 0°-90°) the one with the largest phi1 value [START_REF] Wenk | The representation of orientation distributions[END_REF]. These three equivalent triplets in the Euler space result from a rotation about the three-fold axis of the cubic sample. Then, the color code used in Figure 5a and 5b is a proportional correspondence between the three Euler angles and values from 0 to 255 respectively for the three Red, Green and Blue (RGB) channels of the 8 bit color images. This allows comparing visually the iChord and EBSD orientations.

For a quantitative comparison, the disorientation between homologous pixels in the iChord and the EBSD maps must be computed. To have a pixel per pixel correspondence between the iChord and the EBSD maps, it is necessary to apply an affine transform to the EBSD to correct from spatial distorsions between the two images. After this affine transformation, a position on the sample surface is represented by the same (X,Y) coordinates both in the iChord map and the EBSD map. It is worth noting that this affine transformation does not alter the values of the Euler angles stored in the EBSD maps. The iChord map is presented on Figure 5a, the EBSD map after affine transformation is presented on Figure 5b, and the disorientation map is presented on Figure 5c. On Figure 5c, dark blue corresponds to 0° of disorientation, and white corresponds to a disorientation of 5° or higher.

of the article). An affine transform has been applied to the EBSD map to correct for the spatial distortions between the ion images and the EBSD image scan.

It is clear from Figure 5c that all the grains have been indexed with a disorientation under 5° compared to EBSD results. The mean disorientation over all the grains is 3.1°. We can conclude that, globally, our method allows the crystallographic orientation of large numbers of grains to be determined via a series of images showing channeling contrast.

Decreasing the acquisition time

A potential application of this method could be to rapidly obtain information about the texture of a material. To minimize the acquisition time, the rotation series should consist of the least possible number of images. The number of images can easily be reduced, for example from 360 to 90 in two different manners: (i) a rotation step of 2° and a rotation range of 180°, or (ii) a rotation step of 4° and a rotation range of 360°. The results of these two "virtual" experiments (images were simply extracted from the same set of images as previously) are presented in Figure 6. Clearly, restricting the rotation range reduces the accuracy of the orientation map, with more incorrect indexations (Figure 6a). On the contrary, taking a rotation step of 4° but keeping the full rotation range of 360° gives a disorientation map (Figure 6b) identical to the one obtained by using the full set of images.

Figure 6 (color online): (a) disorientation between the averaged EBSD map and an iCHORD map using a reduction of dimension of (step 2° : range 180°) -90 images used, see section 5.2. (b) disorientation between the averaged EBSD map and an iCHORD map using a reduction of dimension of (step 4° :

range 360°) -90 images used, see section 5.2. The max disorientation found among the grains in (b) is 4.8°.

One can wonder how far we can decrease the number of images while still preserving a good final accuracy of the orientation map. As a second attempt, we tried a setting of (rotation step : rotation range) equal to (12° : 360°), which means only 30 images were needed. With these conditions, the orientation map obtained, presented on Figure 7, exhibits a slightly higher proportion of incorrect indexation compared to the (1° : 360°) map. A full statistical study would be necessary to evaluate the effect of this degradation on the texture characterization. In the present proof-of concept experiment, we did not perform this analysis simply because it will remain a case-by-case study without any universal tendency.

Figure 7 (color online): disorientation between the averaged EBSD map and an iCHORD map. Dark

blue corresponds to 0° disorientation and white corresponds to a disorientation of 5° or higher. Grain boundaries have not been taken into account in this map. The dimensional reduction in this case was (step 12° : range 360°) -see section 5.2., which means that only 30 images were needed to obtain this disorientation map.

Discussion

We have proposed a new method based on channeling contrast in order to measure the crystallographic orientation of polycrystals by acquiring a series of images from a tilted sample in a FIB-SEM. Using a model sample, it has been proven that the same information as the EBSD technique can be obtained. It means that all the data treatments usually made after an EBSD experiment in order to characterize the statistical distribution of grain orientation (pole figure, texture, etc.) can be attained with this method. The main difference with respect to EBSD is the way to obtain the Euler angle triplet (phi1, PHI, phi2) on every point of a ROI. Moreover, the acquisition and processing time are potentially small, with only a few tens of images being necessary and a very small processing time (once optimized) due to its nature -a search by similarity in a database. If the time to obtain texture information is the key-point, there is also another way to speed up the processing, by considering a mean orientation for each grain. From the image series, it is in principle easy to obtain the contour of each grain by an edge detection and averaging z-projections over the whole image stack, and consequently to deduce the positions of each grain center. Interestingly, whereas this operation is hazardous when using a single 2D image (generally because of a poor signal-to-noise ratio and/or the lack of contrast between some adjacent grains), it is quite straightforward when considering the rotation stack of images as acquired in the present method. Figure 8 presents an example for the present case with edge detection on single images and a summation of such images over the whole stack. Figure 8a, b, c andd show respectively: the original image stack (first image represented), an example of edge detection on a single image of the stack, the sum of all edge-detected images of the stack (360 images) and the sum of only 10 edge-detected images from the stack. Some details become hardly visible when comparing the sum of 360 edge-detected images and the sum of only 10 edge detected images. Summing only ten images is the lowest limit to be able to detect clearly all the grain boundaries. A mean indexation for each grain can then be determined by indexing a few points around its center. In our case, the ROI contains a few hundred grains, which means that only a few hundred indexing operations are necessary, instead of indexing all pixels of the image (264 × 186 = 49104 pixels). This approach would not be possible with classical EBSD because a complete scan has to be done first before knowing the contours of all the grains. However, some interesting solutions have been proposed to determine the position of the grain boundaries using the forward-scatter detectors installed around the EBSD screen (FSD diodes) [START_REF] Day | A comparison of grain imaging and measurement using horizontal orientation and colour orientation contrast imaging, electron backscatter pattern and optical methods[END_REF], the EBSD screen itself as a detector [START_REF] Wright | Electron imaging with an EBSD detector[END_REF] or using a scanning mesh refinement technique [START_REF] Yang | Adaptive orientation imaging microscopy[END_REF][START_REF] Wu | Mapping the mesoscale interface structure in polycrystalline materials[END_REF].

Another improvement brought by this new method, as compared to EBSD, might be the lateral spatial resolution of the orientation map. The small depth of penetration of the ion beam as well as the geometry of acquisition are favorable concerning the spatial resolution. On the contrary, a drawback of using a relatively high atomic number ion beam such as Ga, is surface sputtering as well as possible ion implantation and amorphisation effects. In this case, using an He+ beam would be a benefit [START_REF] Veligura | Channeling in helium ion microscopy: Mapping of crystal orientation[END_REF]. However by correctly choosing the ion beam conditions (energy, dose, beam overlap), these effects will be negligible during the acquisition time of the azimuthal series.

Another important aspect is the possibility to implement this method as a stand-alone technique.

As mentioned in section 4.2, we used an EBSD map of the same ROI to determine the intensity loss associated to each plane. This does not mean that an EBSD pre-calibration will be systematically required when applying the method to any new material. Intuitively, one may suggest that the intensity drop off should be directly linked to the scattering power of the corresponding plane family. Indeed, it was easy to check afterwards that using coefficients, for each plane, proportional to their structure factor allows pertinent intensity profiles to be obtained. This implies that theoretical or more appropriately simulated profiles could thus be constructed without the help of an EBSD map. It further means that two materials with the same symmetry (i.e., cubic as tested so far) but with different motif and Bravais lattice could be differentiated.

To totally evaluate the general interest of this new method, it will be necessary to further investigate (i) the spatial resolution of the technique, (ii) the accuracy of the indexation, (iii) if it can be applied to non-cubic materials, (iv) if it allows distinguishing between two (or more) crystallographic phases in the sample, (v) if it can be applied even if the channeling contrast is not very intense (for instance for heavily deformed materials). These goals are not in the scope of the present investigation, but will be the subject of further work.

A final question remains concerning the applicability of this technique to the SEM. Although the availability of FIB instruments is increasing, SEM instruments remain by far more largely installed in most laboratories and industries dealing with materials, in research and development or quality / control units. The possibility to obtain orientation maps only by imaging, without supplementary hardware would then be much welcomed. Among the challenges to take up with electrons, one has to consider the smaller intensity of the channeling contrast, the larger depth of penetration in the sample, or the implementation of the theoretical background necessary to construct the profile databases.

Conclusion

This work has demonstrated the proof-of-concept that orientation maps can be obtained by other means than EBSD and X-Ray diffraction. The new method presented here relies on the acquisition of a series of ion-induced SE images in a FIB, presenting a contrast due to the channeling effect. No additional hardware is necessary for this method. The first results, obtained on a TiN sample, show that an orientation map with a mean orientation error of less than 5° can be obtained, by comparison with an EBSD map acquired on the same area, and considered as a reference.

Preliminary tests have shown that only a few tens of images of the same ROI are necessary to compute an orientation map. The key idea of our method is the ability to compute a database of theoretical intensity profiles for any orientations. An experimental intensity profile extracted from the image series at a position (X, Y) is then compared to the database, with the best match giving the most probable orientation at this position.

Experimentally, a large number of microscope parameters have to be refined to determine the best acquisition conditions, in relation with the studied materials. For instance, the acceleration voltage, the detector type or the working distance were not varied during our experiments, therefore more work needs to be done to understand their role. Theoretically, our program to generate the intensity profile for each orientation is quite simple. Further experiments will be carried out to evaluate how this new method can be used (i) for non-cubic materials, (ii) for multi-phased materials and (iii) with an electron beam instead of an ion beam.

Aside from the possibility to obtain orientation maps, this work also provides two interesting ways to improve the characterization and understanding of microstructures by means of ion or electron microscopy. At first, acquiring a rotation series of images may significantly ease the detection of grain boundaries in polycrystals by simply averaging all the images as depicted in figure 8. Secondly, it has been shown that the crystallographic structure factor of the planes under the beam seems to drive their loss in intensity, when the beam arrives parallel to the corresponding plane. This is an alternative view that could complete the "particle" theory usually invoked to explain the channeling by ions [START_REF] Levi-Setti | Ion channeling effects in scanning ion microscopy with a 60 keV Ga< sup>+</sup> probe[END_REF].

Figure 1 (

 1 Figure 1 (color online): (a) schematic of the sample together with its orthonormal reference frame (e x

Figure 2 (

 2 Figure 2 (color online): Diagram representing the successive process operations from the raw images

Figure 3 (

 3 Figure 3 (color online): Graphical User Interface displaying the crystallographic orientation (the three

Figure 4 :

 4 Figure 4: (a) graph showing the indexing time (seconds) of images of different sizes (in pixels). The

Figure 5 (

 5 Figure 5 (color online): (a) Orientation map obtained with the iCHORD method presented in this work.

Figure 8 :

 8 Figure 8: (a) first image of the rotation image series. (b) edge detection on a single image of the

Acknowledgments: the CLYM (www.clym.fr) is gratefully acknowledged for the access to the FIB instrument which was financially supported by the CNRS, the Région Rhône-Alpes, the 'GrandLyon' and the French Ministry of Research and Higher Education.