Marie-Claire Villeval Mars

Marie-Claire Villeval
email: villeval@gate.cnrs.fr

Experimental Economics: Contributions, Recent Developments, and New Challenges

Keywords: A120, C900, D010, D020 experimental economics, neuroeconomics, quantitative methods, field experiments

Although economics has long been considered as a non-experimental science, the development of experimental economics and behavioral economics is amazingly rapid and affects most fields of research. This paper first attempts at defining the main contributions of experiments to economics. It also identifies four main trends in the development of experimental research in economics. The third contribution of this paper is to identify the major theoretical and methodological challenges faced by behavioral and experimental economics.

Economics was traditionally considered as a non-experimental science that must rely on direct observation of the world because of a supposed lack of control for important aspects of economic decision-making [START_REF] Samuelson | Economics[END_REF] cited by Friedman and Sunde, 1994). This view has been challenged by the development of behavioral economics [START_REF] Camerer | Behavioral Game Theory: Experiments in Strategic Interaction[END_REF] that largely relies on experimental evidence. The laboratory offers a controlled environment that allows researchers to measure the impact of exogenous institutional changes while guaranteeing ceteris paribus conditions. The attribution of the Prize in Economic Sciences 2002 to V. Smith for his research works on market mechanisms and to D. Kahneman for having integrated insights from psychology into economics has consecrated this change in perspective.

Experiments are now carried out in most economic fields, including macroeconomics, and the number of publications incorporating experimental work is blooming. This enthusiasm for experiments in economics is partly justified by a willingness to test the predictions delivered by game-theoretical models and by decision theory. The repeated observation of anomalies in the lab has given rise to the conviction that behavior is not only and not always driven by money maximization but that utility is also affected by emotions and affect and that economic reasoning depends on cognitive conditions. This has provoked a need for further explorations. Another reason for the success of experimental economics can be found both in the search for the appropriate counterfactual and in the difficulties faced by researchers to get access to survey data. This paper does not aim at presenting a state-of-the-art of this method in economics but attempts at identifying how this method is contributing to the renewal of research and knowledge in economics. In a first section, it briefly recalls its major findings and emphasizes the diversity and the richness of the protocols tested. In a second section, it tries to identify new tendencies in the use of the experimental method in economics. In a last section, it presents a reflection about the theoretical and methodological challenges faced by this method.

A major contribution of experiments to economics

V. Smith has emphasized three functions of experiments in economics: testing theories, generating data and contributing to institutional design. Regarding the first two functions, major contributions have been made in the analysis of markets under perfect or imperfect information, in the understanding of individual decision-making (loss aversion, endowment effect, ambiguity aversion, law of small numbers, preference reversal, …) and in the study of social interactions (altruistic punishment, trust, fairness, inequality aversion, …). Some myths have then been challenged by the experimental research [START_REF] Bazerman | Economics Wins, Psychology Loses, and Society Pays[END_REF]: in particular, preferences are not always stable and consistent; the agents do not have perfect knowledge of their preferences; they do not always use all the available information to make decisions. This is especially striking in the analysis of incentives. Many studies have shown the impact of incentives on the individual behavior in terms of effort decisions, as predicted by the agency theory; we have however accumulated many proofs of evidence that in some circumstances, incentives may have a negative impact on the effort level chosen by the agents (incentives may crowd-out intrinsic motivation; they may be ineffective in the presence of income targeting; etc.). Indeed, testing economic theory, experimental research has identified several limitations to strategic thinking: these limitations are related to the existence of social motivations, of psychological and judgmental biases, of limited cognitive abilities and emotions. On the other side, other experiments have simultaneously confirmed the strong influence of market forces and competition on human decisions. The coexistence of these results influences the agenda of further research.

The contribution of lab experiments to institutional design is also noticeable. Experiments are able to address several public policy questions. They can be used by companies that plan to reform their organization, their strategy on the market, or their human resources policy.

Experiments provide useful tests of alternative market designs (for example, the design of auctions of broadcast bandwidths or the organization of pollution permit markets). They offer a cheap alternative to small-scale policy implementation. Progress has also been made in the analysis of various funding modes of collective goods or political campaigns (for example, [START_REF] Houser | Selling Favors in the Lab: Experiments on Campaign Finance Reform[END_REF]. Another contribution is related to the fact that experiments can help in analyzing the agents' preferences when markets do not exist (for example, approaching the relationship between valuation and choices, or analyzing tax compliance under various auditing regimes). The extent to which policy makers incorporate experimental evidence into policy design is however less clear.

One should now emphasize a fourth contribution of experimental economics: it offers academics an innovative teaching tool. Classroom experiments have expanded dramatically in the 90s, thanks to prominent contributors such as C. Holt (see his Veconlab Software: http://veconlab.econ.virginia.edu/admin.htm, and his recent book on Markets, Games, and Strategic Behavior, 2006). Leading journals, such as the Journal of Economic Literature, include a section devoted to the practical organization of classroom experiments. This evolution has enabled a more active learning process in economics, based on the self-discovery by students of fundamental economic notions, reconciling abstract mathematical modeling with empirical observations (see [START_REF] Holt | Economic Science: An Experimental Approach for Teaching and Research. Presidential Address to the Southern Economic Association[END_REF]. It certainly contributes to counteract the distaste of students for majoring in economics. Several examples of classroom experiments on voting, political institutions, and agenda effects are provided in these journals and website.

New tendencies

The development of its interactions with other social sciences, in particular psychology, characterizes the recent evolution of experimental economics. In particular, behavioral game theory can be seen as a willingness to blend « experimental evidence and psychology in a mathematical theory of normal strategic behavior » [START_REF] Camerer | Behavioral Game Theory: Experiments in Strategic Interaction[END_REF]. More and more experimental economists pay attention to the work of social psychologists, and vice-versa (de Cremer, Zeelenberg, and Murnighan, 2006). This leads to the import of insights from behavioral decision research into economic analysis, for a better understanding of motivation, social interactions and individual choices.

To a smaller extent, this dialogue has also developed with anthropology, in particular to investigate whether economic behavior is transversal to cultures and to various stages of development (see the replication on the field of standard economic games in small-scale societies by Henrich et al., 2004). Cooperation between experimental economists and political scientists should develop quite easily, regarding the analysis of voting behavior (see the numerous papers in the economics of voting by T.R. Palfrey, such as for example Levine andPalfrey, 2006, or Aragones and[START_REF] Aragones | Spatial Competition Between Two Candidates of Different Quality: The Effects of Candidate Ideology and Private Information[END_REF] or decision-making in committees (in the tradition of [START_REF] Fiorina | Committee Decisions under Majority Rule: An Experimental Study[END_REF]. Despite these fruitful interactions, it should however be noted that important methodological differences remain, especially regarding the role of monetary incentives, the status of deception and the use of context (see [START_REF] Croson | Contrasting Methods and Comparative Findings in Psychology and Economics[END_REF][START_REF] Hertwig | Economists' and psychologists' experimental practices: How they differ, why they differ and how they could converge[END_REF].

A second major evolution of experimental economics is related to the analysis of rationality and its limitations, thanks particularly to the development of neuroeconomics (brain imaging, EEG, PET and fMRI) and physionomics (skin conductance measures, heart rate, measures of hormones, study of pupil dilation). Opening the "black box of brain activity" [START_REF] Camerer | Neuroeconomics: How Neuroscience can Inform Economics[END_REF] aims at establishing the role of controlled processes (deliberation and learning) vs. automatic processes that are below consciousness in economic decision-making.

It also aims at analyzing the brain during decision-making to inform about the relationships between the ('hot') affective and ('cold') cognitive processes involved in performing an economic task. These works lead to questions basic standard assumptions in economics regarding choice and social interactions and contribute to a better understanding of the importance of limited foresight and emotions in economic decisions.

A third direction of evolution is the reinforcement of the realism of economic experiments and the diversification of protocols. Indeed, the external validity of experiments is still an issue. The sterility of the lab environment has many advantages, regarding the conditions of control and the measure of exogenous changes, but it raises some concerns regarding the possibility to extrapolate evidence collected in the lab to real settings. One can also criticize the lack of variability of socio-demographic characteristics of the standard student-subject pools or their lack of expertise. A partial answer can be found in the current development of field experiments, probably initiated by [START_REF] Levine | A Model of Agenda Influence on Committee Decisions[END_REF] in their study of the impact of agenda manipulation on voters' decisions. Whatever artefactual, framed or natural, field experiments constitute an important complement of lab experiments to allow general inferences [START_REF] Harrison | Field Experiments[END_REF].

The reinforcement of realism can also be achieved by combining experimental tests with survey data (see for example the call for contribution to the ANES Panel Study on voting behavior). We also observe a diversification of subject-pools, an increase in the sample size (facilitated by notably the development of internet-based experiments), or the use of representative samples [START_REF] Harrison | Estimating Individual Discount Rates in Denmark: A Field Experiment[END_REF]. This is accompanied by the design of less abstract tasks, the elicitation of valuation instead of induced valuation through virtual goods, the import of parameters values from the field. This evolution is however not without raising other potential concerns (loss of control for the cost of effort, risk of uncontrolled selection or recruitment biases, …).

Major theoretical and methodological challenges

Experimental research has revealed many anomalies in behavior that refute the predictions delivered by economic theory. However, as pointed out by [START_REF] Camerer | Behavioral Economics[END_REF], « the goal of behavioral economics is not just to create a list of anomalies. The anomalies are used to inspire and constrain formal alternatives to rational-choice theories ». These anomalies do not condemn economic theory per se but they certainly call for a transformation of the theory and a search for consistency. This raises of course a more fundamental question regarding the nature of the revision process. One possibility is to work on updating standard theory. For example, experiments have widely documented the importance of fairness in economic transactions. A possibility would be to incorporate additional constraints (for example an equity constraint) in the design of contracts, but the general framework of the agency theory would be preserved. An alternative is to develop a theory that would constitute an alternative to the rational choice theory and that would explain conflicting behaviors in various contexts. This option is explored by the so-called dual theories, opposing controlled, "cold", or deliberative processes to automatic, "hot", or affective processes (see [START_REF] Camerer | Neuroeconomics: How Neuroscience can Inform Economics[END_REF][START_REF] Loewenstein | We Can Do This the Easy Way or the Hard Way': Negative Emotions, Self-regulation and the Law[END_REF], or by the cognitive hierarchy models in which each player assumes that his strategy is the most sophisticated [START_REF] Camerer | A Cognitive Hierarchy Theory of One-shot Games[END_REF]).

Among the theoretical challenges that need to be addressed by further research, a preferential treatment should be given to the following questions. Should one incorporate social preferences, and non-monetary sources of motivation in the utility function and where should one put the limit in the manipulation of this function? When do market forces dominate and when do social preferences dominate? When do selfish people contaminate the socially-oriented subjects? Some suggest that in games where strategies are complements, a small number of socially-oriented subjects may change the behavior of selfish people [START_REF] Camerer | When Does 'Economic Man' Dominate Social Behavior? Science[END_REF]. Answering this fundamental question requires to understand why « non-standard » behavior (i.e. non predicted by rational-choice theory) survives, declines or develop. This calls for a development of research on the market implications of emotions and on how the market "disciplines" emotions. More largely, we need to study the relationships between emotions or affective decision-making and efficiency. Another important theoretical challenge is to determine what should be the status of agents' heterogeneity in preferences or abilities. Most experiments reveal a diversity of preferences and orientations that is not usually predicted by the theory, and that may have an impact on optimal matching in groups. Either one considers this diversity as pure empirical noise or one tries to give a status to heterogeneity in economic theory. Allowing for heterogeneity also means that sorting and self-selection matter. We should therefore also investigate the interactions between self-selection, emotions and rationality in decision-making.

From a methodological point of view, many questions regarding the collection of data remain unanswered. Do we test mean behavior or do we search for limit cases? A challenge is that a lab experiment must reproduce a theoretical model but it must be at the same time able to mimic the field context reliably. These two orientations may not be necessarily compatible: one can consider that a pragmatic attitude must prevail (the design of an experiment must depend on the research issue), but one can also consider that it is an epistemological issue (experiments aim at reproducing models and parallelism is a secondary issue). The discussion is still open. In the same vein, how can one interpret the differences revealed by the comparison of field and lab experiments? It is important to be able to determine whether they are related to differences in values, in optimization variables, or in strategies. Another point is related to the grasp of time in lab experiments, especially since we are aware of time-inconsistency in human decisions [START_REF] O'donoghue | Procrastination on Long-Term Projects[END_REF]. But how can one reproduce medium-term or long-term in the lab [START_REF] Gneezy | Putting Behavioral Economics to Work: Testing for Gift Exchange in Labor Markets using Field Experiments[END_REF], and how can one reproduce infinite-horizon games? The next point is related to the methodological implications of collected evidence on the impact of smartness and mood on behavior. It probably suggests that we should increase the control for subjects' smartness, attention, and mentalizing abilities. Similarly, a control for mood may be justified in some experiments since it may influence affective reasoning [START_REF] Kirchsteiger | Your Morals are Your Moods[END_REF]. In the same perspective, it matters to collect more data on the individual background of the experimental subjects to control for possible sample biases (for example personal wealth and social origin are likely to influence behavior in bargaining games; the number of siblings may affect decisions in public good or coordination games; etc.). Another concern may however arise: to which extent can we collect additional and peripheral information about our subjects' abilities and background without polluting behavior in the game? A fourth point is related to the measure of emotions. When is it sufficient to use self-reported measures and when should one use physiological measures?

The statistical analysis of experimental data requires also further improvements. Although some scholars still believe that non-parametric statistics constitute the best method to analyze experimental data, more and more numerous are those who consider that the analysis of experimental data benefits from the application of suitable econometric tools in addition to nonparametric statistics, in particular the econometrics of panel data. The development of « experimetrics » ( [START_REF] Bardsley | The Experimetrics of Public Goods: Inferring Motivations from Contributions[END_REF] allowed by a deeper cooperation between experimentalists and econometricians would enable scholars to improve the techniques adapted to the specific nature of the experimental data (in particular, the existence of small samples with a few independent observations but repeated observations). In particular, there is a need for the development of techniques enabling a treatment of the potential selection biases in repeated sequential games, or appropriate to the analysis of the diversity of strategies (for example, the socalled finite mixture models).

Conclusion

In sum, experimental economics constitutes a powerful method to develop knowledge in economics. It helps in identifying empirical regularities and so doing it inspires new research questions. It contributes to the renewal of economic theory by incorporating emotions, affect and cognition in the analysis of behavior. This evolution has been favored by a more open-minded attitude towards other disciplines, in particular psychology and neurosciences. This buoyant field of research fosters excitement for this method and the competition between experimental economists contributes to the rapid development of knowledge in the field of behavioral economics. On the downside of this evolution is the observation that experiments revealing anomalies raise more curiosity and attention than those that do not refute theory, but do not usually attempt to revise the theory. Another drawback is that it is more exciting for scholars to test new protocols than to replicate previous ones, since it pays less to replicate a protocol than to design a new experiment. Therefore, if one considers that replication is crucial for certifying the robustness of scientific results, it would be suitable to provide more incentives for the replication of experiments, in particular in terms of publications in high-rank economic journals.


E A Hodille

X Bonnin

Régis Bisson

T Angot

Charlotte S Becquart

J.-M M Layet

C Grisolia

Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials

come

Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials

Introduction

Tungsten (W) is a serious candidate material for plasma facing components (PFC) in fusion devices such as ITER and DEMO thanks to its thermal properties. Due to plasma/wall interaction, W-based PFC will be subject to high particle fluxes (10 20 -10 24 m²/s) of hydrogen isotopes (HIs). Fuel particles implanted in the subsurface region of PFC could diffuse into the tungsten matrix and be trapped deeper in the bulk. It causes safety issues because of the regulation-limited amount of tritium in the vessel walls as well as operation concerns due to possible uncontrolled HIs recycling fluxes that can affect global plasma stability. So, migration and trapping of deuterium have been extensively studied; see for instance reviews by Causey [1],

Skinner et al. [2] and Tanabe [3].

The development of a global tokamak wall model interacting with hydrogen isotopes by taking into account all the physical processes (particles implantation, migration, trapping, outgassing…) is necessary in order to extrapolate the fuel retention in W-based plasma facing materials for ITER. In this approach, macroscopic rate equations (MRE) models are an efficient way to investigate migration and trapping of HIs in metallic materials from nanometers to centimeters scales. The rate equation model for hydrogen diffusion including hydrogen trapping in materials was originally discussed by McNabb et al. [4] and used for a number of metals such as steel [5,6] and tungsten [7,[START_REF] Ogorodnikova | Hydrogen and Helium Recycling at Plasma Facing Materials[END_REF]9,10,11].

In the present work, a code based on a rate equations model has been developed to deal with the trapping of HIs in W. It is named MHIMS (Migration of Hydrogen Isotopes in MaterialS) [12] and it can be seen as a light version of the HIIPC models developed by Sang et al. [13] for tokamaks inventory simulations. In addition to the reasonable computational resources necessary for running MHIMS, we will detail here how MHIMS is a good tool to extract, from laboratory experiments, fundamental parameters of the HIstungsten interaction and how it can be used for experimentally relevant predictions. MHIMS is able to model particles implantation with traps creation during plasma wall interaction, particles depth profiles evolution in the material and Thermal Desorption Spectrometry (TDS) measurements. In the first part of this paper, the equations of the model are described (section II). Then, simulation results are benchmarked against well controlled laboratory experiments from the literature performed with polycrystalline tungsten samples that we consider as a reference case (section III). Using this benchmarked set of parameters, we compare MHIMS simulations of the evolution of the retention as a function of fluence and implantation temperature with a larger set of laboratory experiments. The objective of this part of our work is to test the robustness of our model and its predictive ability (section IV). Finally, predictions on deuterium outgassing at room temperature are made and their influences on the interpretation of laboratory experiments are discussed (section V).

II. Description of MHIMS Model

An earlier introduction to our model can be found in [12], but we provide here a more comprehensive and updated account. Fig. In our MRE model, HIs are split into two populations: mobile (or solute) and trapped species.

𝐶 𝑚 stands for the concentration of mobile particles and 𝐶 𝑡,𝑖 for the concentration of trapped particles in the i th trap type. In the following, the concentrations will be normalized to the metal density i.e. they are expressed in atomic fraction (at.fr.).

The temporal variation of each population is described by Equation ( 1) and (2).

𝜕𝐶 𝑡,𝑖 𝜕𝑡 = -𝑆 𝑡𝑟𝑎𝑝,𝑖→𝑚𝑜𝑏𝑖𝑙𝑒 + 𝑆 𝑚𝑜𝑏𝑖𝑙𝑒→𝑡𝑟𝑎𝑝,𝑖 (1)
𝜕𝐶 𝑚 𝜕𝑡 = 𝐷(𝑇) ⋅ 𝜕²𝐶 𝑚 𝜕𝑥² -∑ 𝜕𝐶 𝑡,𝑖 𝜕𝑡 + 𝑆 𝑒𝑥𝑡 (2)
where 𝑆 𝑚𝑜𝑏𝑖𝑙𝑒→𝑡𝑟𝑎𝑝,𝑖 is the source of mobile particles being trapped, 𝑆 𝑡𝑟𝑎𝑝,𝑖→𝑚𝑜𝑏𝑖𝑙𝑒 corresponds to trapped particle being detrapped and Sext is the exterior source of particles entering the volume. 𝐷(𝑇) = 𝐷 0 ⋅ 𝑒 -𝐸 𝐷 𝑘⋅𝑇 is the diffusion coefficient of HIs in the metal matrix (in m².s -1 ) with the energy barrier 𝐸 𝐷 represented in Fig. 1 (k is the Boltzmann constant). Dealing with HIs, the diffusion coefficient has to be mass dependent, thus as in [10,11] we used 𝐷 𝐻𝐼 = 𝐷 𝐻

√𝐻𝐼 𝑎𝑡𝑜𝑚𝑖𝑐 𝑚𝑎𝑠𝑠

for the diffusion of any Hydrogen Isotope. As a remark, if an equilibrium is considered between trapped and solute particle, the model correspond to Oriani's one [5] which envisaged the HIs trapping and migration as an apparent diffusion from trap to trap: in such model, the retention variation with fluence would vary as fluence 0.5 (i.e. limited by diffusion).

For each trap type, there are a finite number of available traps. We also assume that each trap captures only one HI atom. The trap density is noted 𝑛 𝑖 and it can evolve with space (inhomogeneous spatial distribution) and with time, since we include trap creation in this model (see further in the text) in contrast to some previous MRE models [7,10,11]. The number of solute sites is fixed. We call 𝑛 𝑠𝑜𝑙𝑢𝑡𝑒 the number of solute sites per tungsten atoms and it is considered that 𝑛 𝑠𝑜𝑙𝑢𝑡𝑒 ≫ ∑ 𝑛 𝑖

𝑖

. Therefore we consider that for each trap, its first neighbor is a solute site. It is also considered that the solute concentration is low: 𝐶 𝑚 ≪ 𝑛 𝑠𝑜𝑙𝑢𝑡𝑒 which means that for each trapped particle, there is a free neighbor solute site which permits detrapping.

Following these assumptions, 𝑆 𝑚𝑜𝑏𝑖𝑙𝑒→𝑡𝑟𝑎𝑝,𝑖 and 𝑆 𝑡𝑟𝑎𝑝,𝑖→𝑚𝑜𝑏𝑖𝑙𝑒 can be expressed as follows [START_REF] Ogorodnikova | Hydrogen and Helium Recycling at Plasma Facing Materials[END_REF]:

𝑆 𝑡𝑟𝑎𝑝,𝑖→𝑚𝑜𝑏𝑖𝑙𝑒 = 1 𝜏 𝑎 ⋅ 𝐶 𝑡,𝑖 (4)
𝑆 𝑚𝑜𝑏𝑖𝑙𝑒→𝑡𝑟𝑎𝑝,𝑖 = 1 𝜏 𝑏 ⋅ 𝐶 𝑚 𝑛 𝑠𝑜𝑙𝑢𝑡𝑒 ⋅ (𝑛 𝑖 -𝐶 𝑡,𝑖 ) (5)
Here, 𝜏 𝑎 and 𝜏 𝑏 are time constants for respectively the detrapping and trapping processes. As in [7,10,11], the detrapping time constant 𝜏 𝑎 can be expressed as a frequency term. Using the energy barrier 𝐸 𝑇,𝑖 = 𝐸 𝐷 + 𝐸 𝐵,𝑖 of Fig. 1,

1 𝜏 𝑎 = 𝜈 0 ⋅ 𝑒 - 𝐸 𝑇,𝑖
𝑘⋅𝑇 and Equation (4) becomes:

𝑆 𝑡𝑟𝑎𝑝,𝑖→𝑚𝑜𝑏𝑖𝑙𝑒 = 𝜈 0 ⋅ 𝐶 𝑡,𝑖 ⋅ 𝑒 - 𝐸 𝑇,𝑖 𝑘⋅𝑇 (6)
where 𝜈 0 is the pre-exponential factor (or attempt frequency) in s -1 . It is important to specify the value of 𝜈 0 when this set of equations is used to fit TDS experiments. Indeed, for a same detrapping energy, the simulated peak can be shifted by about 35 K if 𝜈 0 changes only by an order of magnitude. In the following, 𝜈 0 is taken equal to 10 13 s -1 accordingly with [7,10]. In [START_REF] Ogorodnikova | Hydrogen and Helium Recycling at Plasma Facing Materials[END_REF]9], 𝜈 0 is expressed as a function of the lattice constant and of the diffusion coefficient and 𝜈 0 ~ 3×10 13 s -1 for hydrogen and ν0 ~ 2×10 13 s -1 for deuterium which is of the same order of magnitude. The trapping source can also be understood as a detrapping effect, where the HI atom detraps from a solute site before falling in a trap site of type i. The trapping attempt frequency 1 𝜏 𝑏 ⁄ is usually expressed as a function of the diffusion coefficient by

𝜏 𝑏 = 𝜆 2 𝐷(𝑇)
, where λ is the distance between 2 solute sites or between a solute and a trap site. So the trapping source term (5) becomes:

𝑆 𝑚𝑜𝑏𝑖𝑙𝑒→𝑡𝑟𝑎𝑝,𝑖 = 𝐷(𝑇)⋅𝑛 𝑖 𝜆 2 ⋅ 𝑛 𝑠𝑜𝑙𝑢𝑡𝑒 ⋅ 𝐶 𝑚 ⋅ (1 - 𝐶 𝑡,𝑖 𝑛 𝑖 ) (7)
and 𝑛 𝑠𝑜𝑙𝑢𝑡𝑒 𝑛 𝑖 ⁄ is the average number of solute sites between two trap sites.

The exterior source of particles Sext is present only if a flux of ions is interacting with the materials (plasma or ion beam exposure). The flux 𝜑 of energetic ions is implanted at a certain depth following a distribution law 𝑓(𝑥) that is modeled using SRIM [15] and which depends on the implantation energy of ions. This depth implantation distribution 𝑓(𝑥) is approximated with a Gaussian function in our model. Non-implanted ions are those which are reflected: the reflection coefficient 𝑟 is also calculated with SRIM and we use it to define the exterior source of HIs:

𝑆 𝑒𝑥𝑡 = (1 -𝑟) ⋅ 𝜑 ⋅ 𝑓(𝑥) (8)
When all the particle sources expressions of Equation ( 6), ( 7) and ( 8) are inserted in Equations ( 1) and ( 2) for the temporal evolution of HIs populations, we obtain the following set of equations which are solved numerically:

𝜕𝐶 𝑡,𝑖 𝜕𝑡 = -𝜈 0 ⋅ 𝐶 𝑡,𝑖 ⋅ 𝑒 - 𝐸 𝑇,𝑖 𝑘⋅𝑇 + 𝐷(𝑇)⋅𝑛 𝑖 𝜆 2 ⋅ 𝑛 𝑠𝑜𝑙𝑢𝑡𝑒 ⋅ 𝐶 𝑚 ⋅ (1 - 𝐶 𝑡,𝑖 𝑛 𝑖 ) (9)
𝜕𝐶 𝑚 𝜕𝑡 = 𝐷(𝑇) ⋅ 𝜕 2𝐶 𝑚 𝜕𝑥 2 -∑ 𝜕𝐶 𝑡,𝑖 𝜕𝑡 + (1 -𝑟) ⋅ 𝜑 ⋅ 𝑓(𝑥) (10)
The diffusion coefficient 𝐷(𝑇) of HIs in tungsten has been measured by Frauenfelder [14] for hydrogen: 4.1 × 10 -7 ⋅ 𝑒 -0.39𝑒𝑉 𝑘⋅𝑇 (𝑚 2 . 𝑠 -1 ). This value is in agreement with other experimental data as shown in several review papers [1,2,3] and has been rationalized with theoretical calculations [16,17,18]. We therefore take the experimental value to parameterize our model.

Tungsten crystallizes with a bcc lattice. It has been shown by ab initio calculations [16] that hydrogen diffuses between tetrahedral interstitial sites. Thus, 𝜆 = 𝑎 𝑊 2⋅√2 ~110 pm with 𝑎 𝑊 the lattice constant of tungsten equal to 316 pm and 𝑛 𝑠𝑜𝑙𝑢𝑡𝑒 = 6 since there are 12 tetrahedral sites and 2 tungsten atoms in the bcc lattice.

As suggested on the left side of Fig. 1, there can be potentially a recombination barrier for atoms at the surface of the metal since they escape in the vacuum mostly in the form of HI molecules. Nevertheless, experimental measurements of HI retention in tungsten strongly suggest that desorption of HIs from the surface is far from being the rate limiting step [1,19]. So, we consider that the boundary condition is a Dirichlet boundary condition [10]: 𝐶 𝑚 (𝑥 = 0, 𝑡) = 0.

We will now turn back to the time evolution of the total trap density ∑ 𝑛 𝑖 𝑖 . In the work of Haasz and coworkers [20], the deuterium retention dependency on the deuterium ion incident flux has been studied. It was observed that the retained fraction of impinging deuterium ions was low and flux dependent until a flux threshold of 10 18 D.m -².s -1 above which it became high and flux independent (for flux up to 5×10 19 D/m 2 /s). Subsequently, the analysis of the evolution of the deuterium depth profile as a function of fluence [21,22] showed that deuterium concentration is inhomogeneous throughout the bulk and can be divided in three zones. The highest HIs concentration is found in the ion stopping range, then a decreasing tail of high HIs concentration is measured in the first µm depth and finally a HIs low concentration plateau is found. While the latter zone is related to natural (intrinsic) defects in the material, the two former high concentration zones were shown to increase much faster than expected from the increase of fluence. These observations were explained by ion-induced trap creation through two different processes [22]. The first trap creation process is operative in the ion stopping zone. It could be due to trap creation through local supersaturation of HIs in the ion stopping range (i.e.

solute/mobile concentrations exceeding the solubility limit), subsequently inducing stresses in the material and creating extended defects such as vacancy clusters and bubbles [9,22,24]. The second trap creation process occurring in the first µm results from the diffusion of the HIs located in the implantation zone and being "pushed" toward the bulk by the local stress field. The resulting increase of HIs concentration could also lead to the plastic deformation of the tungsten matrix from deuterium super-saturation [20,21,23].

In our model these two processes are taken into account through a time variation of the density of created traps 𝑛 3 expressed as:

𝑑𝑛 3 𝑑𝑡 = (1 -𝑟) ⋅ 𝜑 ⋅ [(1 - 𝑛 3 𝑛 3amax ) ⋅ 𝜂 3a ⋅ 𝑓(𝑥) + (1 - 𝑛 3 𝑛 3bmax ) ⋅ 𝜂 3𝑏 ⋅ 𝜃 𝑥 𝑝 (𝑥)] (11)
where both traps creation rates are proportional to the implanted flux 𝜑. Thus, for both type of traps, the amount of created traps is proportional to the HIs fluence. In the brackets of Equation (11), the first term corresponds to the traps creation in the stopping zone f(x) and follows the expression proposed by Ogorodnikova et al. [22] based on Duesing et al.'s work [25]. In this zone, the traps density is limited to 𝑛 3𝑎𝑚𝑎𝑥 and the creation rate 𝜂 3𝑎 is the number of traps created per implanted ion. The second term in the brackets of Equation ( 11) corresponds to traps creation in the first µm of the material. This expression was not included in previous models. The form of the expression of the second trap creation process resembles the one of the first trap creation process but with a different depth distribution. We chose an empirical depth distribution for these traps with the simple form 𝜃 𝑥 𝑝 (𝑥) = 1 𝑥 𝑝 between 0 and 𝑥 𝑝 , 0 elsewhere and xp ~ 1-3 µm [21]. This second type of created traps is density limited to 𝑛 3𝑏𝑚𝑎𝑥 and the creation rate is 𝜂 3𝑏 . In addition, to agree with the observed flux dependence by Haasz and coworkers [20], a threshold of trapping creation process is implemented: only if 𝜑 > 10 18 𝐷. 𝑚 -2 . 𝑠 -1 , does traps creation take place.

Firstly, we note that the first term in the brackets of Equation ( 11) could also account for trap creations through inelastic collisions of incident impurities with W atoms. However, in the set of experiments we chose to simulate the amount of impurities in the beam was carefully controlled thus we neglect this contribution. Secondly, the choice of the empirical depth distribution for the second term in Equation ( 11) is based on its resemblance with the erfc() found a change in surface morphology related to a change of flux, so 𝑛 𝑚𝑎𝑥 and 𝜂 𝑚𝑎𝑥 may also depend on the flux. However, this was observed for higher fluxes (10 21 -10 22 D/(m²s)) and at higher fluence (10 26 D/m 2 ) than for the works of Ogorodnikova et al. [9,22] and Haasz and coworkers [20,27,28,29] we want to simulate in the present article.

functional
Finally, we emphasize that MHIMS is designed to simulate a typical implantation -TDS retention measurement experiment by considering 3 different phases:

 The implantation period, lasting 𝑡 𝑖𝑚𝑝𝑙𝑎𝑛𝑡𝑎𝑡𝑖𝑜𝑛 (𝑠) at a sample temperature 𝑇 𝑖𝑚𝑝 (𝐾),

 The resting period between the end of the implantation and the beginning of the retention measurement. Here the samples are maintained at constant temperature 𝑇 𝑟 (𝐾) for a period lasting 𝑡 𝑟𝑒𝑠𝑡 (𝑠), the resting time,  The TDS phase during which the temperature is increased with a given heating ramp 𝛽 (expressed in 𝐾. 𝑠 -1 ).

MHIMS will be used to simulate laboratory experiments and to this end we will use published experimental parameters representative for these 3 phases. It will only deal with results obtained with polycrystalline samples.

III. MHIMS benchmarking on laboratory experiments

MHIMS is used in [12] to catch critical parameters from TDS experiments and used it in HIIPC to simulated trapping/migration/desorption of HIs during real tokamak material thermal cycle. In this section, we come back to this simulation in order to explain in more detail the results shown in [12]. We used laboratory experiments with hot-rolled W from Ogorodnikova et al. [9] since these data were obtained with an all in-situ apparatus, i.e. there was no air exposure between deuterium implantation and retention measurements with TDS. In their article, o 𝑛 3𝑏𝑚𝑎𝑥 = 1 × 10 -2 , 𝜂 3𝑏 = 2 × 10 -4 , 𝑥 𝑝 = 1 µ𝑚.

Ogorodnikova
The values of these traps densities fulfill the condition 𝑛 𝑠𝑜𝑙𝑢𝑡𝑒 >> 𝑛 𝑖 . In addition, the concentration of mobile particles is calculated to be ~1 × 10 -7 , the condition 𝐶 𝑚 << 𝑛 𝑠𝑜𝑙𝑢𝑡𝑒 is also fulfilled.

The simulated TDS measurement given by MHIMS reproduces well both the position and the shape of the experimental TDS spectra, the relative error of the integral of each spectra being less than 10%, and they do not miss the evident broad aspect of desorption peaks of deuterium in contrary to Ogorodnikova et al.'s original fit. Note that the detrapping energies of trap 1 and trap 3 in MHIMS agree well with Ogorodnikova et al. detrapping energies [9]. At first glance, the addition of a second intrinsic trap to better reproduce the experimental data of Ogorodnikova et al. may look artificial. However, intrinsic traps are usually associated with dislocations, grain boundaries and/or vacancies [9,22] which are all present in tungsten polycrystals. To these 3 types of defects could correspond at least two different binding energies, as it has been shown recently by Xiao and Geng [30]. Furthermore, ab-initio calculations have recently demonstrated that multiple hydrogen atoms can be trapped around a single defect, leading to a distribution of binding energies [18,31,32]. These theoretical predictions justify our use of more than a single intrinsic trap.

In a MHIMS simulation (fluence = 10 22 D/m²), we have access to the evolution of the deuterium density profile in the bulk of the material which allows for a better understanding of experimental observables such as a TDS peak. For example in Figure 5 of the article of Ogorodnikova et al. [9] one can see that TDS peaks shift to higher temperature as fluence is increased from 10 22 to 10 23 D.m -2 . This TDS peak behavior is well reproduced in our simulation (not shown) and can be rationalized as being due to particles trapping deeper in the tungsten bulk (Fig. 2(b)), since to a higher fluence corresponds a longer implantation at constant flux and thus a longer time for particles diffusion into the bulk. This difference in deuterium density profile shows up in TDS experiments as a delay of peak appearance i.e. the TDS peak appears shifted to higher temperature. Besides, looking at the bulk part of the depth profile (Fig. 2(b)), it can be seen that the deuterium concentration approaches the traps density. The concentration of trapped particles depends on a balance between Equations ( 6) and (7). Therefore, at a given temperature, it mainly depends on the couple (𝜈 0 , 𝐸 𝑇,𝑖 ) and the presence of mobile particles near a trap. With a large flux like the one used here, the source of mobile particles is sufficient to saturate both traps 1 and 2. The traps creation in the stopping zone can be seen in the zoom of the sub-surface part.

The density of created trap in the area reaches 2×10 -2 trap/W for a fluence of 10 22 D.m -2 and 10 -1 trap/W (= n3amax = maximum density) for a fluence of 10 23 D.m -2 which is very high. It is probable that in this zone W is deeply damaged. The other population of created traps (until 1 µm) is not visible for fluence of 10 22 D.m -2 and appears only at higher fluences (10 23 D.m -2 ).

However, the influence of these created traps is seen in the simulated TDS spectrum: they are at the origin of the tail at high temperature observed between 650 and 700 K. In this study, the traps created up to 1 µm are associated with a high detrapping energy but in [22] they are associated to low detrapping energy. The association is ambiguous due to difficulties in understanding correctly what happens during trap creation. Further experimental studies are needed to understand properly these processes. However, the retention in these created traps does not influence the observations, the discussion and the conclusion in the following.

IV. Comparison of the MHIMS model with various laboratory experiments

Variation of retention vs fluence

To test the reliability of the trapping parameters used above, the evolution of retention with fluence at two different implantation temperatures obtained with MHIMS is shown in Fig. 3 together with the experimental data from Ogorodnikova et al. [9] and Tian et al. [27]. Several fluences from 10 21 to 10 24 D/m² are simulated for implantation temperatures of 300 K and 500 K.

The MHIMS model parameterized on data from Ogorodnikova et al. at 300 K implantation temperature predicts deuterium retention as a function of fluence in excellent agreement with the data of Ogordnikova et al. measured at 300 K and 500 K. Absolute differences are at most of a factor of two. However, when comparing the MHIMS results with the 500 K implantation data of Tian et al. [30] a mismatch by a factor of 5 can be observed at low fluence. This difference is explained by the pre-implantation treatment of the tungsten samples used. The parameters of the model used here have been determined on experiments using a sample pre-annealed at 1273 K during an hour prior to deuterium ion implantation while the data of Tian et al. [27] have been obtained with the sample pre-annealed at 950 K during 30 min. As it is discussed in [9,12], this pre-annealing has an impact on the intrinsic traps density and so on the overall retention. This is especially visible at low fluences because trap creation is not the dominant retention process. At higher fluences, the difference diminishes because the amount of created traps does not depend on the pre-implantation treatment.

Finally, in Fig. 3, it is observed that the absolute variation of retention with fluence is different for 300K and 500K implantations. At 300 K, the retention varies according to the power law fluence 0.55 while at 500 K, the retention increases as fluence 0.7 . This change in the power law exponent can be interpreted as a transition from diffusion-limited retention (300 K) to trapcreation limited retention (500 K). Indeed for implantation at 500 K, detrapping from the lowest energy intrinsic trap is easy and thus the retention is mainly due to retention in the highest energy intrinsic trap and the extrinsic created trap, i.e. the fraction of retention due to extrinsic trap increases. Since the extrinsic trap is linearly dependent on the ion flux (within a certain domain, see Eq. ( 11)), the power exponent for 500 K implantation is closer to unity as compared to 300 K implantation.

Variation of retention with implantation temperature

The effect of implantation temperatures was simulated in more details by comparison with an extensive set of experimental data from the group of Haasz [27,28,29] between experiment and simulation with a retention decrease of one order of magnitude for 𝑇 𝑖𝑚𝑝 increasing from 300 K to 600 K and then a supplementary retention decrease of another order of magnitude between 600 K and 700 K. The simulated TDS measurements and the differential deuterium retention in each trap type are shown in Fig. 4(b) and Fig. 5, respectively, for all these implantation temperatures. Because trap 1 desorption occurs at ~400 K, the amount of trapped deuterium in trap 1 will decrease rapidly for implantation temperatures above 300 K. Since deuterium detraps efficiently from trap 1 and diffuses into the bulk, it will end up preferably trapped in the highest energy intrinsic trap 2. This explains the relative growth of trap 2 populations as compared to trap 1 between 300 K and 400 K as seen in Fig. 5. From 300 K to 600 K, the population in the extrinsic traps 3 stays nearly constant. It does not significantly increase because trap 3 is created near the implanted zone and when such a trap is created, it is almost immediately filled. Trap 3 population begins to decrease after 600 K as the sample temperature is high enough for detrapping (and thus desorption) to occur from this trap.

Even if the decreasing trends are observed both in experimental data and in simulation results, the normalized retention from experiments seems to be underestimated by the MHIMS simulation. To explain this discrepancy, two effects may be considered.

First, the pre-implantation treatment: authors from ref. [28,29] use similar annealing as in ref. [9],

which was used here to benchmark the set of parameters of the MHIMS model. So, we can exclude pre-implantation treatment as the reason of the observed difference.

Second, the resting time: Haasz et al. specified that the resting time of their sample is in the 16 h -72 h range [29]. As we will show in the next section, this experimental parameter is the reason for the discrepancy in Figure 4 between MHIMS simulations and experiments.

V. Effect of the resting time

Studies on deuterium retention in tungsten rarely have considered the influence of the resting time between ion implantation and the retention measurement. One of the reasons is the common assumption that any eventual desorption between implantation and TDS experiment is only due to solute particles. However, with the presence of a low energy trap like the MHIMS lowest energy intrinsic trap (trap 1), there could be significant desorption during a long (hours or more) resting time. To investigate this potential effect, several simulations were realized with different resting times, from 50 s (corresponding to the conditions of Ogordonikova et al. Reaction Analysis (NRA) measurements when it is used to extract trap density from deuterium depth profiles: if the resting time between deuterium implantation and NRA measurements is of the order of several tens of hours, then the trapped HI density in low energy intrinsic traps like trap 1 will be strongly underestimated or even not seen. Several tens of hours is actually quite a usual time delay before ex-situ NRA measurements. A final remark: the above results are important when deuterium implantation is performed at 300 K. For example, if we model with MHIMS an implantation at 500 K followed by a resting time at 300 K, no desorption during resting time is observed: if trap 1 is not filled during implantation due to high temperature implantation it will not be filled from trap 2 population at 300 K.

As an application of our observation with MHIMS of an effect of the resting time for 300 K implantation, we integrated this resting time parameter in our attempt to reproduce the published data from the group of Haasz [27,28,29], where the resting time is around 70 h -80 h.

The results are shown in Fig. 7. As compared to Fig. 4(a), it is clear that a better agreement between the model and the experimental data is achieved. Thus, it appears that the resting time between the end of ion implantation and the beginning of retention measurements may be an important parameter that any rate equation models should include, as it is the case in MHIMS.

Conclusion

MHIMS is a code based on a macroscopic rate equations model including ion-induced trap creation and resting time between the end of ion implantation and the characterization time.

It is applied to deal with retention of HIs in tungsten. We use a well-controlled (in-situ) TDS experiment from literature to benchmark this model. It is found that HIs can be trapped in three types of trap: two intrinsic traps (detrapping energy 0.87 eV and 1.00 eV) and an extrinsic trap created by ion irradiation (detrapping energy 1.50 eV).

The benchmarked MHIMS code is used to simulate the evolution of retention at difference fluences at 300 K and 500 K. These results are compared to a set of experimental data from different groups. Simulation results are in good agreement with experimental data and quantitative discrepancies are linked to difference in sample preparation. It is observed that for deuterium ions implantation at 300 K, the retention is limited by diffusion in the bulk (trapping in intrinsic trap). On the other hand, for ion implantation at 500 K, retention evolution as a function of fluence is limited by traps creation (extrinsic trap). The evolution of retention with implantation temperature in the 300 -700 K range is also well reproduced by the MHIMS model.

Similarly to experimental data, it is shown that retention drops by two orders of magnitude with a specific threshold around 600 K. This behavior is rationalized by analysis of the relative deuterium populations in the three types of traps.

Finally, the effect of outgassing at 300 K after ions implantation at 300 K is simulated.

After 300 000 s, a loss of 45 % of the initial retention is predicted together with an enhanced diffusion in the bulk. Taking into account this outgassing during resting time, the MHIMS model achieves a better agreement with experimental measurements of the implantation temperature dependency of deuterium retention. It is deduced that, if not taken into account, this room temperature outgassing would lead to misinterpretation of NRA data for trap density determination, particularly when a too long resting time is used between implantation at 300 K and NRA measurements.

  form obtained from non-steady-state diffusion. It however behaves with a steeper decay that could render the non-linear effects related to the super-saturation process. A more exact functional form should account for the intertwined influence from the incident flux and the material temperature that define the diffusion rate toward the bulk as well as the concentration in the implantation zone. In such, this is an extended theoretical study that would necessitate further experimental data for validation and thus it is out of the scope of this article. Nevertheless, Ogorodnikova et al. have shown[26] that the two HI concentration distribution components remain for deuterium implantation at 320 K and 500 K. Therefore Equation (11) should render at least qualitatively the evolution of HIs concentration in the bulk of the material even at temperature higher than 300 K. Thirdly, our choice of a fluence dependent model is guided by our goal to simulate experiments which were performed with fluxes comprised between 10 19 to 10 20 D/m²/s and fluence between 10 21 -10 24 D/m 2 . Indeed, Haasz and coworkers observed a HIs retention which is independent of the flux in the range 10 18 -10 19 D/(m 2 s) for the fluence range 10 21 -10 23 D/m 2 . It suggests that the integrated amount of trap created is in this range independent to the flux and only dependent on the fluence. One may note that Lindig et al.[24]

3  4 

 34 et al. simulate the TDS measurements for a fluence of 10 22 D.m -2 with a rate equations model consisting of 2 detrapping energies related to: an intrinsic trap (0.85 eV) and an extrinsic trap (1.45 eV) i.e. a trap created during deuterium ions implantation. We also fit this TDS measurement, simulating implantation at a flux of 2.5×10 19 D/m²/s with 200 eV/D ions (r = 0.56), a resting time of 50 s and a heating ramp up of 8 K/s [9]. Their experimental data and their model adjustment are shown in Fig. 2 (a) together with our MHIMS simulation. We present also the deuterium depth profile just after the resting time (i.e. before the TDS heating ramp begins) determined by MHIMS in Fig. 2(b). For the MHIMS simulation, three traps have been used:  trap 1: a low energy intrinsic trap: 𝐸 1 = 0.87 𝑒𝑉, 𝑛 1 = 1 × 10 -trap 2: a medium energy intrinsic trap: 𝐸 2 = 1.00 𝑒𝑉, 𝑛 2 = 4 × 10 -trap 3: a high energy extrinsic trap: 𝐸 3 = 1.50 𝑒𝑉 with the following density of induced trap: o n3=n3a+n3b o n3a(t=0s)= n3b(t=0s)=0 o 𝑛 3𝑎𝑚𝑎𝑥 = 1 × 10 -1 , 𝜂 3𝑎 = 6 × 10 -4

  while using the set of parameters benchmarked on Ogorodnikova et al. experiments [9]. In the MHIMS modeling, we used a constant fluence of 10 22 D/m² (limited by calculation time constraints which are temperature dependent) and an ion incident energy of 200 eV/D. Several simulations were made for implantation temperatures from 300 K to 700 K. The experimental data relates to different fluences and implantation energies. To compare the simulation results with experimental data, the retention is normalized over the maximum retention in each data set. The results of our simulations are presented in Figure 4. Fig. 4(a) shows the same decreasing trends

  [9]) to 300 000 s (~83 h) i.e. resting times typical of the Haasz group. The implantation was simulated at room temperature (300 K) with incident energy of 200 eV/D and a fluence of 10 22 D/m². After 300 000 s of resting time, also at 300 K, it is seen that 45 % of the initial retained particle were lost (Fig6 (a)).In Fig.6(a), the overall deuterium retention and the relative deuterium population in each of the three traps is plotted as a function of resting time, while in Fig.6(b) is shown the corresponding evolution of the deuterium depth profile. The concomitant decrease, respectively increase, of the deuterium population in trap 1, respectively trap 2 with resting time, shows that deuterium in trap 1 was detrapped and partly retrapped in trap 2. However, the decrease of the deuterium overall retention indicates that part of the deuterium detrapped from trap 1 diffused to the surface and was desorbed in the vacuum. The non-desorbed deuterium population tends to diffuse deeper in the bulk while being trapped in trap 2. After a resting time of 50 s, HIs concentration is saturated in the 1 st µm and is almost equal to the initial intrinsic trap concentration i.e. 𝐶 𝑡,1 (𝑡 = 50𝑠) + 𝐶 𝑡,2 (𝑡 = 50𝑠) ≅ 𝑛 1 + 𝑛 2 = 1.4 × 10 -3 . After a resting time of 300 000 s, almost all trap 1 sites are empty and the HIs are present up to 1.8 µm: the concentration of HIs is almost equal to 𝐶 𝑡,1 (𝑡 = 300000𝑠) + 𝐶 𝑡,2 (𝑡 = 300000𝑠) ≅ 𝑛 2 = 0.4 × 10 -3 . This result has potentially a profound consequence for the interpretation of Nuclear

Figure 1 .

 1 Figure 1. Potential energy diagram for a hydrogen atom in tungsten.

Figure 2 .

 2 Figure 2. (a) MHIMS simulated TDS spectrum compared with experimental and model results of Ogorodnikova et al. [9]. (b) MHIMS simulated depth profile of deuterium just before the beginning of the TDS measurement for 2 fluences (10 22 and 10 23 D.m -2 ). Implantation simulated at a flux of 2.5×10 19 D/m²/s with 200 eV/D ions, resting time is 50 s, the heating ramp rate is 8 K/s.

Figure 3 .

 3 Figure 3. Evolution of deuterium retention vs 200 eV/D ions fluence. Experimental data for 300 K -470 K [9] and 500 K [30]. Simulated results for 300 K and 500 K. "Fitted point" stands for the MHIMS simulation of Fig 2(a) used to calibrate the code. Implantation simulated at a flux of 2.5×10 19 D/m²/s with 200 eV/D ions, resting time is 50 s.

Figure 4 .

 4 Figure 4. (a) Evolution of (normalized) deuterium retention as a function of implantation temperature. Crosses: MHIMS simulations. Other symbols: experimental data. (b) MHIMS simulation of TDS measurements for different implantation temperature betweenn 300 K and 700 K. Implantation simulated at a flux of 2.5×10 19 D/m²/s with 200 eV/D ions, resting time is 50 s, the heating ramp up is 8 K/s.

Figure 5 .

 5 Figure 5. Evolution of retention in each trap in function of implantation temperature. Implantation simulated at a flux of 2.5×10 19 D/m²/s with 200 eV/D ions, resting time is 50 s.

Figure 6 .

 6 Figure 6. (a) Evolution of traps population as the resting time is increased (50 s -300 000 s). (b) Deuterium depth profile as a function of resting time. Implantation simulated at a flux of 2.5×10 19 D/m²/s with 200 eV/D ions.

Figure 7 .

 7 Figure 7. Same as Figure 4(a) but including in the MHIMS simulation with a 300 000 s resting time. Implantation simulated at a flux of 2.5×10 19 D/m²/s with 200 eV/D ions.

Figure 1 H 2 Fluence

 12 Figure 1

  Figure 6

  1 presents the general energy diagram of a hydrogen atom inside a metal with two types of trapping sites present in the bulk of the material. 𝐸 𝑅 + 𝐸 𝑆 is

the energy barrier of the hydrogen atom penetration into the metal matrix. 𝐸 𝑅 is the energy barrier for the H atom to overcome in order to reach the surface (a preliminary step for surface recombination). 𝐸 𝐷 is the barrier of diffusion of H in the metal matrix through solute sites. 𝐸 𝐵,1 and 𝐸 𝐵,2 are binding energies of HIs located in two different trap types present in the metal. Thus, energy barriers to come out of these traps, called detrapping energies, are respectively 𝐸 𝑇,1 = 𝐸 𝐵,1 + 𝐸 𝐷 and 𝐸 𝑇,2 = 𝐸 𝐵,2 + 𝐸 𝐷 .


Azzouz Dermoune
email: azzouz.dermoune@univ-lille1.fr

Cristian Preda
email: cristian.preda@polytech-lille.fr

Parametrizations, fixed and random effects

Keywords: General linear model, fixed effect, random effect, cubic spline, smoothing parameter, likelihood, climate change detection

HAL is

Motivation

We consider the problem of climate change detection. The years taken into account and the annual mean temperature are denoted by t 1 < . . . < t n+1 and y 1 , . . . , y n+1 respectively. In our work we consider the additive model

y(i) = s(t i ) + w i , i = 1, . . . , n + 1, (1)
where (w 1 , . . . , w n+1 ) is a Gaussian white noise with the variance σ 2 w , and s(t i ) is the true temperature at the year t i . We model the behaviour of the true temperature s by a random element of the set S 3 (t 1 , . . . , t n+1 ) of cubic splines having the knots t 1 < . . . < t n+1 .

The fixed effect E(s) is a straight line. It belongs to the null space Q -1 (0) of the quadratic form

Q(s) = t n+1 t 1 |s (t)| 2 dt. (2)
Here s , s are respectively the first and the second derivative of the map s. The random effect s -E(s) belongs to an unknown supplementary space E of Q -1 (0) relative to S 3 (t 1 , . . . , t n+1 ) and has the probability distribution exp{-Q(s)/2σ 2 s }1 E (s)ds/Z σs .

Here σ 2 s is a positive parameter which measures the dispersion of the random effect around the space Q -1 (0), and may depend on the space E, and Z σs is the partition function. We assume that the noise (w i ) is independent of s, but its variance σ 2 w may depend on the space E. A popular estimator of the temperature s is given by the following penalized estimation technique [START_REF] Wahba | Spline models for observational data[END_REF] ŝ = arg min{

n+1 i=1 |y(i) -s(t i )| 2 + λQ(s) : s ∈ H 2 }, (4)
where H 2 is the infinite dimensional space of all functions with square integrable second derivative, and λ > 0 denotes the smoothing parameter. Generalized cross-validation techniques are among the automatic methods used to estimate the smoothing parameter see, e.g., [START_REF] Craven | Smoothing noisy data with spline functions[END_REF], [START_REF] Green | Nonparametric Regression and Generalized linear models[END_REF], [START_REF] Wahba | Spline models for observational data[END_REF]. The estimator ŝ (4), for λ fixed, belongs to the set S 3 (t 1 , . . . , t n+1 ) of cubic splines and does not depend on the parametrization of S 3 (t 1 , . . . , t n+1 ). See [START_REF] Dermoune | A consistent estimator of the smoothing parameter in the Hodrick-Prescott filter[END_REF], [START_REF] Dermoune | Multivariate Extension of the Hodrick-Prescott Filter-Optimality and Characterization[END_REF] for a similar study. The concept of fixed and random effects has been applied to the analysis of longitudinal data. See, e.g., [START_REF] Diggle | Analysis of longitudinal data[END_REF], [START_REF] Laird | Random-effects models for longitudinal data[END_REF], [START_REF] Liang | Longitudinal data analysis using generalized linear models[END_REF], [START_REF] Zeger | Longitudinal data analysis for discrete and continuous outcomes[END_REF]. In [START_REF] Lai | Fixed and random effects selection in nonparametric additive mixed models[END_REF] the authors have modelled the fixed effects nonparametrically using truncated series expansions with B-spline basis. They have selected the fixed effects using lasso methodology, while the random effects are estimated using the Newton Raphson algorithm.

In our work the sum of the fixed and the random effects is a random cubic spline. The fixed effects are straight lines. The random effects are the supplementary spaces of the space of straight lines relative to S 3 (t 1 , . . . , t n+1 ). We introduce an original notion of basis separating the space of straight lines (fixed effects) from its supplementary spaces (random effects) relative to S 3 (t 1 , . . . , t n+1 ). We interpret the smoothing parameter as the signal-tonoise ratio σ 2 w /σ 2 s . We show that the estimator of the smoothing parameter is a function of such a basis. We also show that there exists an infinite number of separating bases and we propose comparison criteria between two separating bases using the partition functions and the maximum-likelihood method.

The plan of our work is the following. In Section 2, we introduce the notion of separating bases in a general setting , and show that they determine the shape of the fixed effect and the parametrization of the random effect, but are not enough to determine the shape of the random effect and the parametrization of the fixed effect. In Section 3 we return to the climate change detection. We show that the random effect is parametrized by s (t 1 ), . . . , s (t n+1 ), and the fixed effect is parametrized by two independent linear forms of the vectors (s(t 1 ), . . . , s(t n+1 )), (s (t 1 ), . . . , s (t n+1 )).

We construct four bases separating the fixed and the random effects, and calculate in each basis the maximum likelihood estimators of the fixed and the random effects and the dispersion parameters σ 2 s and σ 2 w . Finally we consider the separating basis as a parameter and we estimate it by maximizing the likelihood.

Separating Basis in a general setting

In this section the set S is any finite dimensional vector space having the dimension p, and s a random element of S. Its probability distribution is defined by the fixed space F having the dimension k < p and the measurable map q : S → R such that for all ν > 0 and for all supplementary space E of F relative to S, the partition function

Z ν := E exp{-q(s)/2ν 2 }ds < +∞.
Here ds denotes the Lebesgue measure on E. The random vector s has the mean (fixed effect) E(s) ∈ F . We say that F is the shape of the fixed effect. The random effect s -E(s) belongs to some supplementary space E of F relative to S. Knowing the space E (the shape of the random effect), s-E(s) has the probability distribution exp{-q(s)/2σ 2 s }1 E (s)ds/Z σs Here σ 2 s is a positive parameter which measures the dispersion of the random effect. Now we need to parametrize the set S in order to define properly an element s ∈ S. A parametrization of S is a one-to-one linear map

Θ : s ∈ S → θ ∈ R p .
The notation Z σs (b k+1 , . . . , b p ) means that the partition function depends on the family b k+1 , . . . , b p . However the component k j=1 θ j b j is only a candidate for the fixed effect , i.e., k j=1 θ j b j is a candidate for the mean of s.

The following proposition shows that there is an infinite number of separating bases, and the choice of a fixed effect depends on the practitioner's aim.

Proposition 2.2.

There is an infinite number of separating bases.

Proof. Starting from any basis (b 1 , . . . , b k ) of F and using the incomplete basis theorem we can construct an infinite number of bases B = (b 1 , . . . , b p ) of S, which achieves the proof.

Two parametrizations Θ 1 and Θ 2 are related by the passage matrix P:

θ (2) i = p j=1 P ij θ (1) j , i = 1, . . . , p, (6)
or equivalently their respective bases B (1) and B (2) are related by b

(2)

i = p j=1 P -1 ji b (1)
j .

Here P -1 denotes the inverse of the passage matrix P. Now we ask the following question. What is the link between two separating bases?

Proposition 2.3. Let B (1) = (b (1)
1 , . . . , b (1) p ) and B (2) = (b

(2) 1 , . . . , b (2)
p ) be two separating bases. Let Θ (1) and Θ (2) be the parametrizations of the element s in B (1) and B (2) respectively.

1) The passage matrix P, given by Θ (2) = PΘ (1) , has the following form

P(1 : k, 1 : k) P(1 : k, k + 1 : p) 0(k + 1 : p, 1 : k) P(k + 1 : p, k + 1 : p) ,
where for 1 ≤ n 1 < n 2 ≤ p, P(n 1 : n 2 , n 3 : n 4 ) is the sub-matrix (P ij : n 1 ≤ i ≤ n 2 , n 3 ≤ j ≤ n 4 ) of P and 0(k + 1 : p, 1 : k) denotes the sub-matrix of the null matrix 0(1 : p, 1 : p).

2) We have (θ

(2) k+1 , . . . , θ (2)
p ) = P(k + 1 : p, k + 1 : p)(θ (1)
k+1 , . . . , θ

p ) .

3) The components θ

(2) 1 , . . . , θ (2)
k may depend on the all parameters (θ

(1) 1 , . . . , θ (1)
p ). 4) The elements (b

(2) 1 , . . . , b (2)
k ) depend only on (b

(1) 1 , . . . , b (1)
k ). But the elements (b (2) k+1 , . . . , b (2) p ) may depend on the whole basis (b

(1) 1 , . . . , b (1) p ). Proof. 1) If θ (2) i = 0, i = k + 1, . . . , p, then from p i=1 θ (2) i b (2) i = p i=1 θ (1) i b (1)
i ,

we have

k i=1 {θ (2) i b (2)
i -θ

(1)

i } = p i=k+1 θ (1) i b (1)
i .

The left-hand side term belongs to F , and the right-hand side term belongs to the supplementary space span(b

(1) k+1 , . . . , b (1) p ) of F . It follows that θ (1) i
= 0 for all i = k + 1, . . . , p, which achieves the proof of 1). The proof of 2), 3) and 4) is a consequence of 1).

The vector s has in two separating bases B (1) , B (2) the following decompositions

s = k i=1 θ (1) i b (1) i + p i=k+1 θ (1) i b (1) i = k i=1 θ (2) i b (2) i + p i=k+1 θ (2) i b (2) i . 1) If the vector θ (2) k+1 , . . . , θ (2)
p has the PDF exp{-q( p j=k+1 θ (2) j b (2) j )/2σ 2 s }dθ (2)
k+1 . . . dθ (2) p /Z σs (b

(2) k+1 , . . . , b (2) p ), (7)
then assertion 2 of Proposition 2.3 and the formula of change of variables tell us that the vector θ

(1) k+1 , . . . , θ (1)
p has the PDF exp{-q( p j=k+1 θ

i b

(1)

j )/2σ 2 s }dθ (1)
k+1 . . . dθ (1) p /Z σs (b

(1) k+1 , . . . , b (1) p )
where

Z σs (b (1) k+1 , . . . , b (1) p ) = Z σs (b (2) k+1 , . . . , b (2) p )/|det{P(k + 1 : p, k + 1 : p)}|.
Hence, a random parameter in the basis B (2) remains random in the basis B (1) . 2) If θ

(2) 1 , . . . , θ (2)
k are fixed and θ

(2) k+1 , . . . , θ (2)
p has the PDF [START_REF] Dermoune | Estimation of noisy cubic spline using a natural basis[END_REF], then for i = 1, . . . , k, the parameter

θ (1) i = p j=1 P -1 ij θ
(2) j may depend on the random vector θ

(2) k+1 , . . . , θ (2)
p . The parameters θ

1 , . . . , θ

k are corrupted by θ (2) k+1 , . . . , θ (1)
p . Hence, a fixed effect in a basis B (2) is in general no longer fixed in another basis B (1) . Definition 2.4. We say that two separating bases B (1) , B (2) are equivalent if the components k i=1 θ

(2) i are the fixed effect of s respectively in the basis B (1) and B (2) and if the passage matrix (see Proposition2.3 assertion 2) satisfies |det{P(k + 1 : p, k + 1 : p)}| = 1.

Remark 2.5. Two separating bases B (1) , B (2) are equivalent if and only if the passage matrix (6) has the following form

P(1 : k, 1 : k) 0(1 : k, k + 1 : p) 0(k + 1 : p, 1 : k) P(k + 1 : p, k + 1 : p) , with |det{P(k + 1 : p, k + 1 : p)}| = 1.
Note that the PDF (5) will concentrate around the minimizers of the map q as Z σs (b k+1 , . . . , b p ) → 0. Now we can compare two separating bases as follows.

Definition 2.6. We say that the parametrization Θ (1) = (θ

(1) 1 , . . . , θ (1) p ) is more concentrated than Θ (2) = (θ (2) 1 , . . . , θ (2) p ) if for all ν > 0, Z ν (b (1) k+1 , . . . , b (1) p ) < Z ν (b (2) k+1 , . . . , b (2) p ).

Four separating bases of cubic splines

We consider the set S := S 3 (t 1 , . . . , t n+1 ) of C 2 cubic splines having the knots t 1 < . . . < t n+1 . We recall that an element s ∈ S is a C 2 map on [t 1 , t n+1 ] and is a polynomial of degree three on each interval [t i , t i+1 [ for i = 1,. . . , n.

More precisely let

p 1 := s(t 1 ), . . . , p n+1 := s(t n+1 ), q 1 := s (t 1 ), . . . , q n+1 := s (t n+1 ), u 1 := s (t 1 ), . . . , u n+1 := s (t n+1 ), v 1 = s (t 1 +), . . . , v n = s (t n +)
be respectively the values of s and its derivatives up to order three on the knots. We have for i = 1, . . . , n,

s(t) = p i + q i (t -t i ) + (t -t i ) 2 u i /2 + (t -t i ) 3 v i /6, t ∈ [t i , t i+1 [.
The following constraint for h i = t i+1 -t i , i = 1, . . . , n guarantees the hypothesis that s is C 2 :

p i + q i h i + u i h 2 i /2 + v i h 3 i /6 = p i+1 , (8)
q i + u i h i + v i h 2 i /2 = q i+1 , (9)
v i = s (3) (t i ) = (u i+1 -u i )/h i . ( 10
)
It is well known [START_REF] Deboor | A Practical Guide to Spline[END_REF] that S has the dimension n + 3. Hence an element s ∈ S is completely defined by n + 3 independent parameters. Note that the quadratic form ( 2) is equal to

n i=1 t i+1 t i |u i + t(u i+1 -u i )/h i | 2 dt = n i=1 (u 2 i + u i u i+1 + u 2 i+1 )h i /3 = u Q 2 u
where the column vector u = (u 1 , . . . , u n+1 ) and the invertible (n + 1) × (n + 1) matrix

Q 2 := [q(i, j)] (11)
is defined by q(1, 1) = h 1 /3, q(n + 1, n + 1) = h n /3, for i = 2, . . . , n, q(i, i) = (h i-1 + h i )/3, for i = 1, . . . , n, q(i, i + 1) = q(i + 1, i) = h i /6, and q(i, j) = 0 elsewhere. The PDF (3) means that the vector u is Gaussian, centred with the covariance matrix σ 2 s Q -1 2 . It follows for any separating basis (b 1 , . . . , b n+3 ) that

Z σs (b 3 , . . . , b n+3 ) = σ n+1 s /det(Q 2 ). ( 12
)
Any other parametrization of the random effect is just a linear transformation of (u 1 , . . . , u n+1 ), e.g., (u 1 , v 1 , . . . , v n ) is also a parametrization of the random effect. See assertion 2 of Proposition 2.3. From Remark 2.5, we derive that for all t ∈ [t 1 , t n+1 ] the fixed effect E(s) is given by E(s(t)) = α + βt, where

α = n+1 i=1 {α (0) i s(t i ) + α (1) i s (t i )} (13)
β = n+1 i=1 {β (0) i s(t i ) + β (1) i s (t i )}. (14)
The coefficients α

(0) i , α (1)
i , β

i , β

(1) i define the parametrization of the fixed effect.

If a priori information tells us that the initial values p 1 , p 2 are not random, then the parametrization Θ 002 = (p 1 , p 2 , u 1 , . . . , u n+1 ) and the corresponding basis B 002 = (b 002 1 , . . . , b 002 n+3 ) are the right choice for separating the fixed effect (p 1 , p 2 ) from the random effect (u 1 , . . . , u n+1 ). The subscript notation 002 is justified by the fact that

p 1 := s(t 1 ) := s (0) (t 1 ), p 2 := s(t 2 ) := s (0) (t 2 ), u 1 := s (t 1 ) := s (2) (t 1 ), . . . , u n+1 := s (t n+1 ) := s (2) (t n+1 ).
If a priori information tells us that (p 1 , q 1 ) are not random then the parametrizations are the right choices for separating the fixed effect (p 1 , q 1 ) from the random effect (u 1 , . . . , u n+1 ). The subscript notation 012 is justified by the fact that p 1 := s (0) (t 1 ), q 1 := s (t 1 ) := s (1) (t 1 ) and u 1 := s (2) (t 1 ), . . . , u n+1 := s (2) (t n+1 ). Similarly, the notation 0123 is justified by the fact that p 1 := s (0) (t 1 ), q 1 := s (1) (t 1 ), u 1 := s (2) (t 1 ) and

Θ 012 := (p 1 , q 1 , u 1 , . . . , u n+1 ), Θ 0123 := (p 1 , q 1 , u 1 , v 1 , . . . , v n ),
v 1 := s (t 1 +) := s (3) (t 1 +), . . . , v n = s (t n +) := s (3) (t n +). It follows for s ∈ S that s = p 1 b 002 1 + p 2 b 002 2 + n+1 i=1 u i b 002 2+i = p 1 b 012 1 + q 1 b 012 2 + n+1 i=1 u i b 012 2+i = p 1 b 0123 1 + q 1 b 0123 2 + u 1 b 0123 3 + n i=1 v i b 0123 3+i .
The passage matrix from the parametrization Θ 002 to the parametrization Θ 012 is given by

p 1 = p 1 , q 1 = (p 2 -p 1 )/h 1 -h 1 u 1 /3 -h 1 u 2 /6, u 1 = u 1 , . . . = . . . u n+1 = u n+1 .
Note that if we decide that (p 1 , p 2 ) is the fixed effect, then (p 1 , q 1 ) is not a fixed effect, because from (8) the parameter q 1 is random with the mean (p 2 -p 1 )/h 1 and the variance h 2 1 var(u 1 /3 + u 2 /6). Inversely, if (p 1 , q 1 ) is the fixed effect, then (p 1 , p 2 ) is not a fixed effect, because again from (8) the parameter p 2 is random with the mean p 1 + h 1 q 1 and the variance h 4 1 var(u 1 /3 + u 2 /6). We conclude that the parametrizations Θ 002 and Θ 012 have the same concentration σ n+1 s /det(Q 2 ) in the sense of Definition 2.6, but are not equivalent in the sense of Definition 2.4.

+ n i=1 s (t i +)b 0123 3+i , (15)
for all cubic spline s. By writing

s(t) = p i + q i (t -t i ) + (t -t i ) 2 u i /2 + (t -t i ) 3 v i /6
on each sub-interval [t i , t i+1 ), and using its smoothness, we obtain the following constraints

  p i+1 q i+1 u i+1   = M i   p i q i u i   + v i   h 3 i /6 h 2 i /2 h i   , (16)
where

M i =   1 h i h 2 i /2 0 1 h i 0 0 1   .
We denote by M i the product i l=1 M l . The equation ( 16) implies for i = 1, . . . , n, that

  p i+1 q i+1 u i+1   = M i   p 1 q 1 u 1   + i l=1 v l M i-l   h 3 l /6 h 2 l /2 h l   . (17)
We will use the notation M -1 i = 0, and M 0 i = I 3 . Let us define, for each i = 1, . . . , n, the piecewise functions,

χ i (t) = 1 [t i ,t i+1 ) (t), χ 1 i (t) = (t -t i )χ i (t), χ 2 i (t) = (t -t i ) 2 χ i (t)/2, χ 3
i (t) = (t -t i ) 3 χ i (t)/6. Now, using (17), we can announce the basis ) is a basis of Q -1 (0) [START_REF] Craven | Smoothing noisy data with spline functions[END_REF]. It is plotted in Figure 1. To compare with ( 13), ( 14

B 0123 = [b 0123 j : j = 1, . . . , n + 3]. Proposition 3.1. We have for t ∈ [t 1 , t 2 ), that b 0123 1 (t) b 0123 2 (t) b 0123 3 (t) b 0123 4 (t) = χ 1 (t) χ 1 1 (t) χ 2 1 (t) χ 3 1 (t) , b 0123 l = 0, l = 5, . . . , n + 3.

We have for

i = 1, . . . , n -1, t ∈ [t i+1 , t i+2 ), and l = 1, . . . , i, that b 0123 1 (t) b 0123 2 (t) b 0123 3 (t) = χ i+1 (t) χ 1 i+1 (t) χ 2 i+1 (t) M i , b 0123 3+l (t) = χ i+1 (t) χ 1 i+1 (t) χ 2 i+1 (t) M i-l   h 3 l /6 h 2 l /2 h l   , b 0123 i+4 (t) = χ 3 i+1 , b 0123 j (t) = 0, j = i + 5, . . . , n + 3.

Construction of the bases B 002 and B 012

Now, we are ready to construct the basis

B 012 = [b 012 1 , . . . , b 012 3+n ]. From the constraints (u i+1 -u i )/h i = v i , i = 1, . . . , n + 1, we have s = p 1 b 0123 1 + q 1 b 0123 2 + u 1 b 0123 3 + n i=1 v i b 0123 3+i = p 1 b 0123 1 + q 1 b 0123 2 + u 1 (b 0123 3 -b 0123 4 /h 1 ) + n i=2 u i (b 0123 2+i /h i-1 -b 0123 3+i /h i ) + u n+1 b 0123 3+n /h n .
It ) is the inverse of the matrix having the rows l 1 = (1, -1/h 1 , 0, . . . , 0), l 2 = (0, 1/h 1 , -1/h 2 , 0, . . . , 0), . . . , l n = (0, . . . , 0, 1/h n-1 , 1/h n ), l n+1 = (0, . . . , 0, 1/h n ). It follows that the partition function

Z σs (b 0123 3 , . . . , b 0123 n+3 ) = ( n i=1 h 2 i )σ n+1 s /det(Q 2 ). ( 18
)
We derive the basis B 002 as follows. The relation To compare with ( 13), ( 14), we observe that b 002

q 1 = (p 2 -p 1 )/h 1 -h 1 u 1 /3 -h 1 u 2 /6,
1 (t) = 1 -t/h 1 , b 002 2 (t) = t/h 1 and then s(t 1 )b 002 1 (t) + s(t 2 )b 002 2 (t) := α + βt, with α = s(t 1 ), β = {s(t 2 ) -s(t 1 )}/h 1 .

Basis B s202

Here we are interested in the parametrization Θ 202 = (s (t 1 ), s(t 1 ), . . . , s(t n+1 ), s (t n+1 )).

The corresponding basis B 202 = (ϕ 1 , . . . , ϕ n+3 ) was constructed in [START_REF] Dermoune | Estimation of noisy cubic spline using a natural basis[END_REF] and is plotted in Figure 1. The subscript notation 202 is justified by the fact that s (t 1 ) := s (2) (t 1 ), s(t 1 ) := s (0) (t 1 ), . . . , s(t n+1 ) := s (0) (t n+1 ), s (t n+1 ) := s (2) (t n+1 ). A cubic spline s is written in the basis B 202 as follows

s = s (t 1 )ϕ 1 + n+1 i=1 s(t i )ϕ 1+i + s (t n+1 )ϕ n+3 .
The basis B 202 does not contain any straight line (see Figure 1). It follows that B 202 is not a separating basis. However we are going to derive from B 202 an interesting separating basis denoted B s202 . We can show [START_REF] Dermoune | Estimation of noisy cubic spline using a natural basis[END_REF] that the quadratic form (2) equals (s (t 1 ), s(t 1 ), . . . , s(t n+1 ), s (t n+1 ))Q 202 (s (t 1 ), s(t 1 ), . . . , s(t n+1 ), s (t n+1 )) ,

where

Q 202 := [ t n+1 t 1 ϕ i (t)ϕ j (t)dt : i, j = 1, . . . , n + 3].
The singular valued decomposition (SVD) of Q 202 = OD(0, 0, c 1 , . . . , c n+1 )O tells us that the diagonal matrix D(0, 0, c 1 , . . . , c n+1 ) is defined by the null The new parametrization θ = O (s (t 1 ), s(t 1 ), . . . , s(t n+1 ), s (t n+1 ) , defines the new basis B s202 = OB 202 := (ψ 1 , . . . , ψ n+3 ) plotted in Figure 2. By construction ψ , ψ 2 is a basis of the null space Q -1 (0) of the quadratic form (2), and then B s202 is a separating basis. More precisely, the new parametrization (θ 1 , . . . , θ n+3 ) is defined by

θ 1 = n+1 j=1 o 1j+1 s(t j ), θ 2 = n+1 j=1 o 2j+1 s(t j ),
and for i = 3, . . . , n + 3,

θ i = o i1 s (t 1 ) + n+1 j=1 o ij+1 s(t j ) + o in+3 s (t n+1
). To compare with (13), ( 14), we observe that the fixed effect θ 1 ψ 1 +θ 2 ψ 2 := α+βt is such that α and β mix all the components s(t 1 ), . . . , s(t n+1 ).

Making the assumption that (s (t 1 ), . . . , s (t n+1 )) is centred, Gaussian with the covariance matrix σ 2 s Q -1 2 (11) is equivalent to saying that (θ 3 , . . . , θ n+3 ) is centred, Gaussian with the covariance matrix σ 2 s Q -1 s202 , where the diagonal matrix Q s202 = D(c 1 , . . . , c n+1 ). It follows that the partition function of (θ 3 , . . . , θ n+3 ) is equal to

Z σs (ψ 3 , . . . , ψ n+3 ) = σ n+1 s n+1 i=1 1/c i . (19)
4. Maximum likelihood estimators of the fixed effect and the dispersion parameters

Let s = 2 i=1 θ i b i + n+3 i=3 θ i b i be a random element of S 3 (t 1 , . . . , t n+1 ) in the basis B = (b 1 , . . . , b n+3
) separating the fixed effect β := (θ 1 , θ 2 ) from the random effect r := (θ 3 , . . . , θ n+3 ) . In the basis B the PDF of the random effect is centred, Gaussian with the covariance σ 2 s Q -1 B , where

Q B = [ t n+1 t 1 b i (t)b j (t)dt : i, j = 3, . . . , n + 3].
The matrix representation of (1) is the following additive mixed model y := (y 1 , . . . , y n+1 ) = Fβ + Rr + w.

Here the column vectors of the matrices

F = (b 1 (t), b 2 (t)), R = (b 3 (t), . . . , b n+3 (t)) (20) are b i (t) := (b i (t 1 ), . . . , b i (t n+1 )) .
If the noise w is white and Gaussian with the variance σ 2 w , then the data y is Gaussian with the mean Fβ and the covariance matrix

Σ(σ 2 s , σ w ) = σ 2 s G 0 + σ 2 w G 1 . Here G 0 = RQ -1 B R , G 1 = I n+1 . The -2 ln-likelihood of the data y in the basis B is equal to (y, B) = (n + 1) ln(2π) + ln[det{Σ(σ 2 s , σ 2 w )}] +trace{(y -Fβ)(y -Fβ) Σ -1 (σ 2 s , σ 2 w )}. (21)
The likelihood equations

∂ β i = 0, i = 1, 2, ∂ σ 2 w = 0, ∂ σ 2 s = 0,
are equivalent [START_REF] Anderson | Asymptotically efficient estimation of covariance matrices with linear structure[END_REF] to

β = {F Σ -1 (σ 2 s , σ 2 w )F} -1 F Σ -1 (σ 2 s , σ 2 w )y, trace{Σ -2 (σ 2 s , σ 2 w )G 2 0 } trace{Σ -2 (σ 2 s , σ 2 w )G 0 } trace{Σ -2 (σ 2 s , σ 2 w )G 0 } trace{Σ -2 (σ 2 s , σ 2 w )} σ 2 s σ 2 w = trace{(y -Fβ)(y -Fβ) Σ -2 (σ 2 s , σ 2 w )G 0 } trace{(y -Fβ)(y -Fβ) Σ -2 (σ 2 s , σ 2 w )} .
To solve the likelihood equations we use the following iterative scheme [START_REF] Anderson | Asymptotically efficient estimation of covariance matrices with linear structure[END_REF].

Starting from β(0) = (F F) -1 F y, we obtain initial σ 2 s (0), σ 2 w (0) by solving the system trace(G 2 0 ) trace(G 0 ) trace(G 0 ) trace(G 1 ) σ 2 s (0) σ 2 w (0) = trace(C(0)G 0 ) trace(C(0)) ,
where C(0) = {y -Fβ(0)}{y -Fβ(0)} . Having σ 2 s (0), σ 2 w (0), we construct

β(1) = {F Σ -1 (σ 2 s (0), σ 2 w (0))F} -1 F Σ -1 (σ 2 s (0), σ 2 w (0))y.
Having β(1), we construct σ 2 s (1), σ 2 w (1) by solving the system trace{Σ -2 (σ 2 s (0),

σ 2 w (0))G 2 0 } trace{Σ -2 (σ 2 s (0), σ 2 w (0))G 0 } trace{Σ -2 (σ 2 s (0), σ 2 w (0))G 0 } trace{Σ -2 (σ 2 s (0), σ 2 w (0))} σ 2 s (1) σ 2 w (1) = trace{C(1)Σ -2 (σ 2 s (0), σ 2 w (0))G 0 } trace{C(1)Σ -2 (σ 2 s (0), σ 2 w (0))}
, where C(1) = {y -Fβ(1)}{y -Fβ(1)} . If the number n+1 of observations y is large then this iterative scheme converges.

If the number n+1 is small, then we need N large i.i.d. copies y(1), . . . , y(N ) of y. In this case we only need the following one-step algorithm. Calculate OLSE

β N (0) = (F F) -1 F ȳ,
and obtain initial estimates of σ 2 s,N (0), σ 2 w,N (0) by solving the system trace(G 2 0 ) trace(G 0 ) trace(G 0 ) trace(G 1 )

σ 2 s (0) σ 2 w (0) = trace(C N (0)G 0 ) trace(C N (0)) ,
Here ȳ = N k=1 y(k)/N and C N (0) = N k=1 (y(k)-Fβ N (0))(y(k)-Fβ N (0)) /N . Having σ 2 s,N (0), σ 2 w,N (0), we construct

β N (1) = {F Σ -1 (σ 2 s,N (0), σ 2 w,N (0))F} -1 F Σ -1 (σ 2 s,N (0), σ 2 w,N (0))ȳ.
Having β N (1), we construct σ 2 s,N (1) and σ 2 w,N (1) by solving the system

trace{Σ -2 (σ 2 s,N (0), σ 2 w,N (0))G 2 0 } trace{Σ -2 (σ 2 s,N (0), σ 2 w,N (0))G 0 } trace{Σ -2 (σ 2 s,N (0), σ 2 w,N (0))G 0 } trace{Σ -2 (σ 2 s,N (0), σ 2 w,N (0))} σ 2 s (1) σ 2 w (1) = trace{C N (1)Σ -2 (σ 2 s,N (0), σ 2 w,N (0))G 0 } trace{C N (1)Σ -2 (σ 2 s,N (0), σ 2 w,N (0))} , where C N (1) = k=1 (y(k) -Fβ N (1))(y(k) -Fβ N (1)) /N . It is known [1] that √ N (β N (1) -β), √ N (σ 2 s,N (1) -σ 2 s , σ 2 w,N (1)
-σ 2 w ) T have a limiting normal distribution with means 0 and covariance matrices {F

Σ -1 (σ 2 s , σ 2 w )F} -1 , [ 1 2 trace{Σ -1 (σ 2 s , σ 2 w )G i Σ -1 (σ 2 s , σ 2
w )G j } : i, j = 0, 1] -1 respectively. Finally having the estimates β, σ2 w , σ2 s of β, σ 2 w , σ 2 s , the random effect estimate Remark 4.1. The fixed effect in the bases B 002 and B 012 does not coincide with the linear regression. It only reflects the linear tendency of the beginning of the period. In the basis B s202 the fixed effect concerns all the period. In the basis B s202 (and unlike the bases B 002 and B 012 ) the columns of the matrix F and the columns of the matrix R (20) are orthogonal. We can show theoretically that this fact implies that the maximum-likelihood estimator β = β(i) = (F F) -1 F y for each iteration i. It follows that the fixed effect F β coincides with the linear regression. See [START_REF] Dermoune | A consistent estimator of the smoothing parameter in the Hodrick-Prescott filter[END_REF], [START_REF] Dermoune | General Linear mixed model and signal extraction problem with constraint[END_REF] for a similar study.

Numerical applications

Real data application: climat change

In the climate change problem we are interested in the annual mean temperature observed in France form 1900 to 2015. For each year, Meteo France (Division des Études et Climatologie, Nord) provided us with the minimum, mean and maximum temperatures in France. In our application we are interested only in the mean temperature. Data is presented in Figure 3.

We illustrate the importance of the parametrization of the fixed effects by considering the bases B 002 and B s202 . The estimation of the parameters of the two models and the -2 ln(Likelihood) values are presented in Table 2 For each model, the fixed and the random effects are plotted in Figure 4. The difference between the two models is illustrated in Figure 5.

The numerical results (Figure 5) confirm the theoretical ones, namely the Remark 4.1. Note that the performance of the basis B s202 is superior to that of B 002 in terms of likelihood. As expected, from Figure 5 we see that the fixed effects in the two models are very different, whereas the estimates of the true signal are very close. In the basis B 002 the random effect compensates and corrects the fixed effect defined only by the first two observations (years 1900 and 1901). The numerical results also show that the basis B s202 is more concentrated than B 002 (

116 i=1 1/c i ≈ 1/27 < < 1/det(Q 2 )).
Appendix: Inference from fixed and random effects in a general setting Let s = k i=1 θ i b i + p i=k+1 θ i b i be a random element of S in the basis B = (b 1 , . . . , b p ) separating the fixed effect θ 1 , . . . , θ k from the random effect θ k+1 , . . . , θ p . The PDF of θ k+1 , . . . , θ p is equal to The mixed model (22) consists of three types of objects: observable random vector y (data), unobservable random vectors θ k+1 , . . . , θ p , w and unknown parameters τ := (θ 1 , . . . , θ k , σ 2 s , σ 2 w ). In general, the probability density function (PDF) of a random vector z is denoted by f τ (y). Here τ is a fixed parameter and y varies. The likelihood of τ given the data y is denoted by L(τ | y). The connection between the (PDF) and the likelihood is given by L(τ | y) = f τ (y). The joint (PDF) of the random effect (θ k+1 , . . . , θ p ) and y is equal to f τ (θ k+1 , . . . , θ p , y) = f σ 2 s (θ k+1 , . . . , θ p )f (θ 1 ,...,θ k ,σ 2 w ) (y) = L(τ | θ k+1 , . . . , θ p , y).

f σ 2 s (θ k+1 , . . . , θ p ) := 1 Z s exp{-q( p i=k+1 θ i b i )/2σ 2 s }.
These were called a joint likelihood by [START_REF] Henderson | Best linear unbiased estimation and prediction under a selection model[END_REF] in the context of Gaussian case and hierarchical likelihoods by [START_REF] Lee | Likelihood for random-effect models[END_REF].

Random effect estimation. Given the data y and an estimation τ of τ , we estimate the random effect by maximizing the joint likelihood of the random effect (θ k+1 , . . . , θ p ) and y given y and τ = τ .

L(τ | θ k+1 , . . . , θ p , y).

Acknowledgements

The research project is supported by Laboratoire Paul Painlevé USTL-UMR-CNRS 8524, the Labex CEMPI ANR-11-LABEX-0007-01 and the project ClinMine -ANR-13-TECS-0009. The authors would like to thank the Editorin-Chief and the anonymous referees for carefully reading the manuscript and giving many helpful and constructive suggestions resulting in the present version. The authors would also like to thank Meteo France Division des Études et Climatologie Nord for producing the French temperature dataset.

  Defining a parametrization Θ is equivalent to the existence of the basis B := (b 1 , . . . , b p ) of S such that for all s ∈ S s = p i=1 θ i b i := Bθ. Definition 2.1. We say that a basis B = (b 1 , . . . , b p ) separates (θ 1 , . . . , θ k ) from (θ k+1 , . . . , θ p ) (or simply B is a separating basis) if (b 1 , . . . , b k ) is a basis of the set of the fixed effect F . We will show in assertion 2 of Proposition 2.3 below, for all separating bases B = (b 1 , . . . , b p ) that the component p j=k+1 θ j b j has the PDF exp{-q( p j=k+1 θ j b j )/2σ 2 s }dθ k+1 . . . dθ p /Z σs (b k+1 , . . . , b p ).

  and the corresponding bases B 012 := (b 012 1 , . . . , b 012 n+3 ), B 0123 := (b 0123 1 , . . . , b 0123 n+3 )

3. 1 . 1 , 1 + s (t 1 )b 0123 2 +

 1112 Construction of the bases B 0123 , B 012 , B 002 First, we construct the basis B 0123 = (b 0123 . . . , b 0123 n+3 ) in which s = s(t 1 )b 0123 s (t 1 )b 0123 3

1 , b 0123 2

 12 The basis B 0123 is a separating basis because (b 0123

  + s (t 1 )b 0123 2 (t) := α + βt, with α = s(t 1 ), β = s (t 1 ).

follows that b 012 1 = b 0123 1 , b 012 2 = b 0123 2 , b 012 3 = b 0123 3 -b 0123 4 /h 1 ,

 11223341 for i = 4, . . . , n + 2, b 012 i = b 0123 i /h i-3 -b 0123 i+1 /h i+1-3 , b 012 n+3 = b 0123 3+n /h n . The passage matrix (b 0123 3 , . . . , b 0123 n+3 ) = P(b 012 3 , . . . , b 012 n+3

and s = p 1 b 002 1 + p 2 b 012 2 /h 1 , b 002 2 = b 2 /h 1 , b 002 3 = b 012 3 -h 1 b 012 2 4 are corrupted by b 012 2 .

 1222233242 The bases B 012 , B 0123 , B 002 are separating bases and are plotted in Figure1. Moreover, according to (18) if h i = 1 for all i, then the bases B 002 , B 012 , B 0123 have the same concentration.

-Figure 1 :Figure 2 :

 12 Figure 1: Representation of bases B 012 , B 0123 , B 002 and B 202 for n=7 equidistant intervals in [0, 1].

r = arg min r {

 r λr G 0 r + y -(F β + Rr) 2 } = ( λG 0 + R R) -1 R (y -F β),where λ := σ2 s σ2 w , and • denotes the Euclidean norm.

  11.0 11.5 12.0 12.5 13.0 13.5 Annual mean temperature from 1900 to 2015 in France years Mean annual temperature (Y)

Figure 3 :

 3 Figure 3: Annual mean temperatures in France from 1900 to 2015.

Figure 4 :Figure 5 :

 45 Figure 4: Mean annual temperatures in France from 1900 to 2015.

  It follows thatf τ (y) = R p-k f (θ 1 ,...,θp,σ 2 w ) (y)f σ 2 s (θ k+1 , . . . , θ p )dθ k+1 . . . dθ p ,wheref (θ 1 ,...,θp,σ 2 w ) (y) = 1 Z w exp{-G(y -Xs)/2σ 2 w }is the (PDF) of y known s and σ 2 w . Parameter estimation: Given the data y, we estimate τ by maximizing the likelihood τ → L(τ | y).

Table 1 :

 1 5.1. Simulation: inference from the bases B 002 , B 012In order to illustrate the importance of the parametrization for the estimation of the true fixed effect, we consider the parametrization B 002 and Estimation using B 012 and B 002 models.

	Parameter θ	p 1	p 2	q 1	σ 2 s	σ 2 w	λ
	True values	0	0.056	1	5	0.1	1/50
	θ012	0.0021 0.1110 0.9775 5.0213 0.0954 1/52.9100
	θ002	0.0022	0.056 0.4419 5.1343 0.0955 1/53.7634
	True values	0	0.172	1	0.5	0.1	1/5
	θ012	0.0107 0.1566 0.9742 0.4721 0.1025 1/4.6058
	θ002	0.0101 0.1867 1.1731 0.4578 0.1024 1/4.4707
	True values	0	0.1521	1	0.1	0.1	1
	θ012	-0.0012 0.1373 1.0109 0.1042 0.0998 1/1.0440
	θ002	-0.0014 0.1517 1.0823 0.1025 0.0912	1/1.12

Table 2 :

 2 . Estimation using B 002 and B s202 models for annual mean temperature.

	Basis	θ1	θ2	σ 2 s	σ 2 w	-2 ln(likelihood)
	B 002	11.486	11.491 1.577 × 10 -5 0.232	172.552
	B s202 -114.968 -58.545 2.755 × 10 -5 0.2283	165.495

by fixing p 1 (true), q 1 (true) (fixed effect) and u 1 (true), . . . , u n+1 (true) is a realization of the random effect u 1 , . . . , u n+1 . We obtain the true value

We generate a sample s 1 , . . . , s N of s 0 , and a sample w 1 , . . . , w N of the noise w. We obtain a sample y 1 , . . . , y N of the data y.

Under the B 002 parametrization and assuming that p 1 , p 2 are the fixed effect, we execute the scoring algorithm and obtain the estimator p 002 1 , p 002 2 , u 002 1 , . . . , u 002 n+1 of p 1 (true), p 2 (true), u 1 (true), . . . , u n+1 (true). We derive

as an estimator of q 1 (true). Now we estimate p 1 (true), q 1 (true), u 1 (true), . . . , u n+1 (true) in the basis B 012 and assume that p 1 , q 1 are the fixed effect. From the estimator p 012 1 , q 012 1 , u 012 1 , . . . , u 012 n+1 , we obtain

as a new estimator of p 2 (true). Table 1 presents the differences between the estimations of the parameters p 1 , p 2 , q 1 , σ 2 s and σ 2 w when the true model is the basis B 012 and B 002 while varying the signal-to-noise ratio λ = σ 2 w /σ 2 s . The results presented in Table 1 show that, for small values of λ, a bad choice of model produces large estimation errors. In our example, for λ = 1/50 the estimation of q 1 = 1 is q 012 1 = 0.9775 in the true model B 012 , whereas q 002 1 = 0.4419 in the wrong model B 002 . According to our simulations, these errors decrease with λ (see the results in Table 1 for λ = 1/5, and λ = 1).