Sankalp Gilda
email: sankalp.gilda@gmail.com

Antoine De Mathelin
email: antoine.de-mathelin-de-papigny@michelin.com

Sabine Bellstedt
email: sabine.bellstedt@uwa.edu.au

Guillaume Richard
email: guillaume-gu.richard@edf.fr

Unsupervised Domain Adaptation for Constraining Star Formation Histories

The prevalent paradigm of machine learning today is to use past observations to predict future ones. What if, however, we are interested in knowing the past given the present? This situation is indeed one that astronomers must contend with often. To understand the formation of our universe, we must derive the time evolution of the visible mass content of galaxies. However, to observe a complete star life, one would need to wait for one billion years! To overcome this difficulty, astrophysicists leverage supercomputers and evolve simulated models of galaxies till the current age of the universe, thus establishing a mapping between observed radiation and star formation histories (SFHs). Such ground-truth SFHs are lacking for actual galaxy observations, where they are usually inferred -with often poor confidence -from spectral energy distributions (SEDs) using Bayesian fitting methods. In this investigation, we discuss the ability of unsupervised domain adaptation to derive accurate SFHs for galaxies with simulated data as a necessary first step in developing a technique that can ultimately be applied to observational data. The code and the data used for the experiments conducted in this work are available on Github 1 .

Introduction

In recent times, many transfer problems have arisen, and various methods exist to solve them. We can broadly classify these into three categories. First, supervised domain adaptation, where only a few labeled target data are available [START_REF] Dai | Boosting for Transfer Learning[END_REF][START_REF] De Mathelin | Adversarial weighting for domain adaptation in regression[END_REF]Motiian et al. 2017a,b). Second, semi-supervised domain adaptation, where in addition to these labeled data, a large amount of unlabeled data is available [START_REF] Kumar | Co-regularization based semi-supervised domain adaptation[END_REF][START_REF] Saito | Semi-Supervised Domain Adaptation via Minimax Entropy[END_REF][START_REF] Tzeng | Simultaneous Deep Transfer Across Domains and Tasks[END_REF]. Finally, unsupervised domain adaptation, where only unlabeled data is available in the target [START_REF] Ganin | Domain-adversarial Training of Neural Networks[END_REF][START_REF] Huang | Correcting Sample Selection Bias by Unlabeled Data[END_REF][START_REF] Richard | Unsupervised Multi-source Domain Adaptation for Regression[END_REF][START_REF] Saito | Maximum classifier discrepancy for unsupervised domain adaptation[END_REF][START_REF] Sugiyama | Direct Importance Estimation with Model Selection and Its Application to Covariate Shift Adaptation[END_REF]. Researchers have shown that adding target labels can increase the performance of the model (Motiian et al. 2017a). This setting is also the most encountered in practice: [START_REF] Cortes | Adaptation Based on Generalized Discrepancy[END_REF], as it is often possible to label at least a few target samples.

Figure 1: Radiation spectrums with corresponding star formation histories can be recorded from cosmological simulation to form a dataset (X, y). This dataset is used to learn a machine learning model. Because of the domain shift between real radiation spectrums and simulated ones, domain adaptation is used to adapt the learned machine learning model to the real observations. However, there are cases in which getting supervised data is not possible. Such is the case of the present astrophysics problem of SFH prediction. Here, we aim to learn the history of our Universe through the evolution of its galaxies' masses throughout their individual histories. We try to obtain this history from the radiation that reaches us on earth -the Spectral Energy Distribution (SED), a function of the brightness of a galaxy with the wavelength of observation.

As galaxies evolve over billions of years, it is impossible to completely capture any single galaxy's SED as a function of time. Astrophysicists then generate models of light through physically-motivated mechanisms, which allow building artificial star formation histories with the corresponding radiation. One of the objectives of these artificial data sets is to find the correspondence between the radiation and the SFH. We try to find a function that can robustly perform this mapping. Classically, Bayesian models have been developed : PROSPECT (Robotham et al. 2020), MAGPHYS (da Cunha, Charlot, and[START_REF] Da Cunha | A simple model to interpret the ultraviolet, optical and infrared emission from galaxies[END_REF], CIGALE [START_REF] Noll | Analysis of galaxy spectral energy distributions from far-UV to far-IR with CIGALE: studying a SINGS test sample[END_REF], BAG-PIPES [START_REF] Carnall | Inferring the star formation histories of massive quiescent galaxies with BAGPIPES: evidence for multiple quenching mechanisms[END_REF], PROSPECTOR [START_REF] Johnson | Bd-J/Prospector: Initial Release[END_REF], and numerous others, more recently MIRKWOOD, [START_REF] Gilda | mirkwood: Fast and Accurate SED Modeling Using Machine Learning[END_REF] showed that deep learning approaches have the potential to infer SFHs from SEDs.

A significant challenge with applying all the above models arises due to the 'domain shift' between simulated and real galaxies. Even a deep neural network (DNN) trained only on the synthetic data sets will falter when predicting the SFHs of real galaxies. Therefore, in this work, we investigate the use of unsupervised domain adaptation as an alternative method of extracting SFHs for individual galaxies without being limited by the simplifications of physical models and parametrizations of conventional SED-fitting techniques. We turn toward cosmological, hydrodynamic simulations, including EAGLE (Schaye et al. 2015a), ILLUSTRISTNG [START_REF] Nelson | First results from the IllustrisTNG simulations: the galaxy colour bimodality[END_REF][START_REF] Pillepich | First results from the IllustrisTNG simulations: the stellar mass content of groups and clusters of galaxies[END_REF] and SIMBA [START_REF] Davé | SIMBA: Cosmological simulations with black hole growth and feedback[END_REF]. These simulate a volume of the Universe from shortly after the Big Bang to the present day, recording the evolution of dark and baryonic matter over time. Galaxies formed within simulations can be studied as a proxy for real galaxies and can be compared directly to observations. They can thus serve to augment our understanding of the evolution of the Universe.

Data

Our training and test data sets consist of SEDs from three state-of-the-art cosmological galaxy formation simulations, where the true physical properties are known -including their true star formation histories. SIMBA [START_REF] Davé | SIMBA: Cosmological simulations with black hole growth and feedback[END_REF], EAGLE (Schaye et al. 2015b;[START_REF] Schaller | The EAGLE simulations of galaxy formation: the importance of the hydrodynamics scheme[END_REF][START_REF] Mcalpine | The EAGLE simulations of galaxy formation: Public release of halo and galaxy catalogues[END_REF], and ILLUSTRISTNG [START_REF] Vogelsberger | Introducing the Illustris Project: simulating the coevolution of dark and visible matter in the Universe[END_REF], with 1,688, 4,697, and 9,633 samples respectively, together constitute a diverse sample of galaxies with realistic growth histories. Specifically, we select galaxies at a redshift of 0this corresponds to simulated galaxies at the "present" epoch. Each SED is comprised of 20 measurements of the "brightness" of a galaxy at different wavelength, in the form of flux density (with units of Jansky) and are the co-variates or features that we train/learn on. The outputs/labels are star formation history (SFH) time series vectors with 29 scalar elements, in units of solar mass per year (M yr -1 ).

The SFH for a galaxy informs us about the net stellar mass generated as a function of time -sum of masses of all stars born, less the sum of stellar mass lost in stellar winds as stars age and eventually die. In other words, SFHs are plots of the star formation rates (SFRs) of galaxies against the lookback time, which for galaxies at z = 0 extends from 0 to the very age of our Universe (∼ 13.8 Gyrs). By stacking (adding) together SFHs of several thousand galaxies from a simulation, and dividing by the volume of the simulated box (i.e. the size of the simulated universe), we can derive the cosmic star formation rate density (CSFRD), which is a wellstudied global property of our Universe. Astrophysicists aim to tune various physical properties in any given simulation until the the CSFRD plot derived from it matches closely to the most-widely accepted one derived from observations (Madau and Dickinson 2014). 2We sequentially train on any two of these simulated data sets, and predict on the third, thus giving us three sets of

Methodology

We consider the problem of prediction of star formation history (SFH) where the learner has access to a data set X ∈ R n×p encoding the radiations of n galaxies with p the number of filters/wavelength-bins, and a data set Y ∈ R n×T giving the corresponding SFH for each galaxy with T is the time length of the history. Each SED (row of X) is comprised of p = 20 measurements of the "brightness" of a galaxy at different wavelength, in the form of flux density (with units of Jansky) (see Figure 2). Below are our pre-processing steps:

1. First, we create three sets of experiments. For each, we use two galaxies in the training and validation sets (with a 9:1 split) and the third galaxy in the test set. 2. Second, we normalize each SFH time series (each row of Y ) by its sum and store the resultant normalized SFH (SFH norm ) and the sum (SFH sum ) separately. This step is needed because of the large dynamic range of the various star formation histories (see Appendix D): SFH curves have a large variety of scales (some increasing to more than 100 whereas others never increase over 0.1. By scaling each SFH time series, we make the learning of the curve trend easier. The learning of the SFHs is now decoupled in the learning of SFH norm and SFH sum . 3. Third, we further ease the learning of the SFH norm by reducing the curves to their first 3 Kernel-PCA components [START_REF] Soentpiet | Advances in kernel methods: support vector learning[END_REF]. The choice of this decomposition method is motivated by extensive experimentation that showed that Kernel-PCA beat both linear PCA (Jolliffe and Cadima 2016) and discrete wavelet transform [START_REF] Shensa | The discrete wavelet transform: wedding the a trous and Mallat algorithms[END_REF] in their ability to recreate the original time series successfully. For each of the three experiments, we provide the Kernel-PCA with a wide range of hyperparameters and pick the ones that can recreate the original SFH time series back, judged according to the DILATE loss metric (Le Guen and Thome 2019) (see Appendix A and Figure 4). The DILATE similarity metric between two time series is defined as an equally weighted average of the dynamic time warping (DTW) and the temporal distortion index (TDI) similarity scores. The number of principal components (3) was chosen by selecting the smallest set that explains at least 80% variance in the validation set. We refer to these kernel-PCA components as SFH kPCA . 4. Fourth, we normalize the input features (the columns from X) via log-scaling, and follow this up by standard scaling normalization. In Figure 2 we visualize the log-scaled flux densities, in units of Janskies, for all three simulations, in 3 out of 20 filters. 5. Finally, we derive KLIEP weights for the training samples in all three experiments.

After performing these steps, we apply a domain adaptation method to correct the shift between the source and target input distributions. We choose the method KLIEP [START_REF] Sugiyama | Direct Importance Estimation with Model Selection and Its Application to Covariate Shift Adaptation[END_REF]), an instance-based method that reweights the sources in order to minimize the KL-divergence between the The feature names at the top of each plot are the names of the filters, each centered at a different wavelength, in which photometry was simulated. We notice that for the three features shown here, all three simulations share the same support, which justifies our decision of using KLIEP. two domains. Instance-based approaches have been widely used to handle regression domain adaptation issues [START_REF] Cortes | Domain adaptation and sample bias correction theory and algorithm for regression[END_REF][START_REF] Huang | Correcting Sample Selection Bias by Unlabeled Data[END_REF][START_REF] Mansour | Domain Adaptation: Learning Bounds and Algorithms[END_REF][START_REF] Sugiyama | Direct Importance Estimation with Model Selection and Its Application to Covariate Shift Adaptation[END_REF], and are particularly robust to negative transfer [START_REF] De Mathelin | Adversarial weighting for domain adaptation in regression[END_REF]. We also visually observe on the marginal distributions for the 20 input filters that all domains have the same support in the feature space, which is the framework considered by KLIEP. Finally, KLIEP has the critical advantage of proposing an unsupervised selection procedure to select a relevant bandwidth.

For each of SFH kPCA and SFH sum , we train a 4 layer feed-forward DNN with drop-out and 256 nodes in each layer. We use ReLu as the activation function, and Adam [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF] as the optimizer, with learning rate of 1e -3 . We train for 200 epochs, with early stopping to prevent over-fitting.

Results

There are two main ways of assessing the SFH outputs derived via our machine learning implementation: 1. Comparing the derived SFH to the true SFH for individual galaxies in the test set. 2. Comparing the predicted ΣSFH to the true ΣSFH for each simulation. True ΣSFH is the sum of SFH for all galaxies in a simulation and is a critical metric enabling us to verify the correctness of input physics in a simulation. We know from observations (Madau and Dickinson 2014) that star formation in the observable Universe peaked about 2 Gyrs after its formation, and all hydrodynamical simulations must produce galaxies that satisfy this observation. By ensuring that the ΣSFH curve inferred from our neural networks matches those from the underlying simulations, we ground our predictions in science while simultaneously enabling appropriate tuning of model architectures, loss functions, and hyperparameters in case of mismatches. Such comparison also enables us to assess any systematic effects in modeling when training using one simulation and comparing predictions on another.

In Figure 3, we plot the true and derived ΣSFH curves for the three distinct test datasets (SIMBA, ILLUSTRISTNG, and EAGLE), where the training data consists of samples from the other two simulations. In these cases, the resulting sample SFHs (and thus ΣSFH) are most deviant from their ground truth vectors. This is unsurprising, as the different simulations have intrinsically different sample SFHs; this is a known difference between different simulations (as shown in Figure A1 of [START_REF] Bellstedt | Galaxy And Mass Assembly (GAMA): A forensic SED reconstruction of the cosmic star formation history and metallicity evolution by galaxy type[END_REF]. Tables 1 and2 show the five metrics tested within this work (MAE, RMSE, BE, DTW, and TDI) for average predictions of SHF, and for ΣSFH, respectively. Figures 7, 8, and9 show examples of true and derived SFHs for individual galaxies within the simulations using KLIEP. Based on the five metrics, each row shows an example of the best-performing galaxy output on the left and the worst-performing galaxy on the right. The first thing to note here is that in most cases, the best-and worst-performing galaxies are different when we use different metrics to pick them. This discrepancy highlights the fact that each metric compares time-series differently, and hence it is difficult to pick one metric as the loss function to minimize.

Another observation made from these three figures is how well our technique can reproduce the stochastic nature of the simulated galaxies' SFHs. As star formation can be an incredibly stochastic process (as is clear from examples such as the top-left panel in Figure 7 and the top-left panel in Figure 8), SFHs can regularly fluctuate between high and low values. In general, we find that such stochastic SFHs are poorly recovered. For the sake of galaxy property analysis, the accurate recovery of overall SFH trends is more important than the recovery of individual star formation rate (SFR) epochs. In the bottom-right panel of Figure 7, for example (galaxy index 333) it can be seen that there is a star formation event early on at ∼ 10 -12 Gyr, and then a secondary star formation event from ∼ 2 Gyr to the present day. Recovering these two main epochs is more crucial than correctly recovering the individual SFR peaks within each epoch. One way of achieving this is to temporally smooth the SFHs of individual galaxies prior to training and testing. We can justify such smoothing of simulated features given that we would never expect the derived SFHs for observed galaxies (the ultimate aim of this work) to reproduce such short-scale stochastic features. As an example of this, see Figure C1 of [START_REF] Robotham | ProSpect: generating spectral energy distributions with complex star formation and metallicity histories[END_REF], where "good" fits to the SFHs from the semi-analytic model SHARK [START_REF] Lagos | Shark: introducing an open source, free, and flexible semi-analytic model of galaxy formation[END_REF]) from the SED-fitting code PROSPECT do not recover the stochastic SFHs. Results presented here are only one part of a continuing effort to apply sophisticated data-driven methods to SED-fitting, a crucial first step for extracting galaxy properties. In future work, we will explore advanced multi-source domain adaptation approaches [START_REF] Zhao | Adversarial Multiple Source Domain Adaptation[END_REF][START_REF] Richard | Unsupervised Multi-source Domain Adaptation for Regression[END_REF] which should more robustly account for the presence of multiple source-domain datasets as well as adaptively reweight the training samples depending on the prediction task, thus providing higher quality results. We are also actively analyzing and attempting to correct systematic biases induced by our various design choices. How do our SFH inferences for galaxies of low mass compare to those of high mass? Is our modeling sufficient to accurately infer SFH from both old and young galaxies, unlike traditional parametric approaches which suffer at earlier epochs in the Universe's history? Answers to such questions will enable us to discover and correct potential biases in our predictions.

Finally, we also plan on vastly increasing the size of our simulated data to account for the immense diversity in the physics of galaxy formation and evolution modelingdifferent initial mass functions (IMFs) (for example [START_REF] Salpeter | The Luminosity Function and Stellar Evolution[END_REF][START_REF] Kroupa | On the variation of the initial mass function[END_REF][START_REF] Chabrier | Galactic Stellar and Substellar Initial Mass Function[END_REF], nebular emission models (for example [START_REF] Bruzual | Stellar population synthesis at the resolution of 2003[END_REF], mass-metallicity relationships (such as those presented by [START_REF] Tremonti | The Origin of the Mass-Metallicity Relation: Insights from 53,000 Star-forming Galaxies in the Sloan Digital Sky Survey[END_REF]Jimmy et al. 2015;[START_REF] Lara-Lopez | Galaxy and mass assembly (GAMA): the connection between metals, specific SFR and hi gas in galaxies: the Z-SSFR relation[END_REF], among others. While this task has a relatively longer horizon-generating and saving new simulations is a compute-intensive task-it is crucial to undertake before data-driven methods such as the one proposed in this paper become trusted by astronomers.

B KLIEP Reweighting

C Extended Results

RMSE

Figure 2 :

 2 Figure2: Kernel density estimate plots of the log of flux densities (FD) for the first three features. The feature names at the top of each plot are the names of the filters, each centered at a different wavelength, in which photometry was simulated. We notice that for the three features shown here, all three simulations share the same support, which justifies our decision of using KLIEP.

Figure 3 :

 3 Figure 3: Global SFH predictions for the three experiments. The curves correspond to the sums over all SFH or predicted SFH.

Figure 4 :

 4 Figure 4: Comparison of SFH reduction approaches for SIMBA. On the left, the evolution of DILATE loss between the true sfh and the reconstructed signal after one of the three transformations: DWT, PCA, kernelPCA. On the right, reconstructed signal for one of the SIMBA SFH.

Figure 5 :

 5 Figure 5: KLIEP brings different domains closer by reweighing source samples. Top: Scatter plot of the first two PCA components in the input space. Bottom: KDE plots of log of the first two features.

Figure 6 :

 6 Figure 6: Global SFH predictions for the three experiments. The curves correspond to the sums over all SFH or predicted SFH.

Figure 7 :

 7 Figure 7: Individuals galaxies for SIMBA from target and source specific clusters.

Figure 8 :

 8 Figure 8: Individuals galaxies for ILLUSTRISTNG from target and source specific clusters.

Figure 9 :

 9 Figure 9: Individuals galaxies for EAGLE from target and source specific clusters.

Table 1 :

 1 Forecasting results for individual SFH. The first column gives the target domain : ILLUSTRISTNG, EAGLE or SIMBA, the two other domains are used as source domains. The computed metrics between each individual SFH prediction and the corresponding groundtruth are averaged in the target domain. 16 th , 50 th , and 84 th quantile values are drawn from the 50 predictions per galaxy, corresponding to the 50 neural networks used. Leading method between UDA and no-UDA (baseline) are shown in bold.

	+0.08 -0.03	5.13 +0.87 -0.89

Table 2 :

 2 Forecasting results for ΣSFH (total star formation history). The first column gives the target domain : ILLUSTRISTNG, EAGLE or SIMBA, the two other domains are used as source domains. Predictions for all samples in the target domain are summed and compared to the true ΣSFH according to the different metrics. 16 th , 50 th , and 84 th quantile values are drawn from the 50 predictions per galaxy, corresponding to the 50 neural networks used. Leading method between UDA and no-UDA

(baseline) are shown in bold. Metrics are scaled (MAE × 1000, RMSE × 1000, BE × 1000, DTW × 1000 for readability.

Table 3 :

 3 Forecasting results, with ILLUSTRISTNG and SIMBA as source domains and EAGLE as the target domain. Predictions for all samples in EAGLE have been averaged, and 16 th , 50 th , and 84 th quantile values are drawn from the 50 predictions per galaxy, corresponding to the 50 neural networks used.

	RMSE	MAE	BE	DTW	TDI
	Base KLIEP MDD DANN DeepCORAL 0.94 +0.23 1.16 +0.1 -0.07 1.15 +0.08 -0.08 1.3 +0.16 -0.11 1.18 +0.08 -0.13 -0.24 KMM 1.2 +0.12 -0.12 CORAL 0.83 +0.24 -0.1	0.84 +0.11 -0.06 0.85 +0.07 -0.09 1.03 +0.11 -0.09 0.86 +0.06 -0.08 0.86 +0.24 -0.22 0.91 +0.08 -0.13 0.63 +0.14 -0.08	-0.4 +0.05 -0.08 -0.4 +0.05 -0.04 -0.37 +0.09 -0.11 -0.4 +0.1 -0.11 0.81 +0.29 -0.35 -0.4 +0.22 -0.17 0.41 +0.1 -0.1	3.58 +0.48 -0.5 3.52 +0.36 -0.32 4.11 +0.57 -0.44 3.59 +0.47 -0.54 2.6 +0.65 -0.59 3.77 +0.87 -1.01 1.68 +0.96 -0.41	0.47 +0.09 -0.11 0.47 +0.17 -0.15 0.73 +0.2 -0.2 0.46 +0.04 -0.07 0.39 +0.13 -0.16 0.56 +0.21 -0.2 0.21 +0.14 -0.04

Table 4 :

 4 Forecasting results for ΣSFH (total star formation history) with ILLUSTRISTNG and SIMBA as source domains and EAGLE as the target domain. Predictions for all samples in EAGLE have been added, and 16 th , 50 th , and 84 th quantile values are drawn from the 50 predictions per galaxy, corresponding to the 50 neural networks used. Metrics are scaled (MAE × 1000, RMSE × 1000, BE × 1000, DTW × 1000 for readability.

	RMSE	MAE	BE	DTW	TDI
	Base KLIEP MDD DANN DeepCORAL 2.19 +0.14 0.93 +0.03 -0.02 0.92 +0.02 -0.02 0.92 +0.02 -0.02 0.94 +0.03 -0.02 -0.17 KMM 1.02 +0.04 -0.04 CORAL 0.95 +0.02 -0.01	0.63 +0.03 -0.02 0.6 +0.02 -0.02 0.61 +0.02 -0.01 0.65 +0.03 -0.02 1.74 +0.15 -0.17 0.68 +0.04 -0.03 0.58 +0.02 -0.02	0.09 +0.05 -0.03 0.02 +0.03 -0.02 0.1 +0.05 -0.05 0.11 +0.06 -0.04 1.06 +0.18 -0.18 0.09 +0.09 -0.08 -0.22 +0.03 -0.02	3.14 +0.13 -0.11 3.25 +0.1 -0.12 3.3 +0.14 -0.23 3.15 +0.13 -0.14 9.81 +0.87 -0.97 3.47 +0.31 -0.11 4.29 +0.09 -0.09	4.6 +0.38 -0.94 3.63 +0.88 -0.95 4.34 +0.49 -0.44 4.5 +0.56 -1.0 7.55 +0.22 -0.2 4.55 +0.76 -0.51 4.4 +0.39 -0.62

Table 5 :

 5 Forecasting results, with ILLUSTRISTNG and EAGLE as source domains and SIMBA as the target domain. Predictions for all samples in SIMBA have been averaged, and 16 th , 50 th , and 84 th quantile values are drawn from the 50 predictions per galaxy, corresponding to the 50 neural networks used.

	RMSE	MAE	BE	DTW	TDI
	Base KLIEP MDD DANN DeepCORAL 1.91 +0.31 0.38 +0.09 -0.07 0.34 +0.05 -0.08 0.29 +0.12 -0.07 0.42 +0.12 -0.07 -0.32 KMM 0.38 +0.08 -0.11 CORAL 0.49 +0.05 -0.06	0.3 +0.06 -0.05 0.29 +0.05 -0.08 0.22 +0.1 -0.05 0.33 +0.1 -0.05 1.8 +0.3 -0.31 0.31 +0.06 -0.09 0.41 +0.04 -0.05	0.14 +0.09 -0.04 0.04 +0.04 -0.04 0.16 +0.08 -0.08 0.19 +0.11 -0.07 1.8 +0.3 -0.31 0.15 +0.14 -0.14 -0.37 +0.06 -0.03	1.17 +0.39 -0.27 0.82 +0.22 -0.14 0.98 +0.47 -0.28 1.38 +0.55 -0.33 9.38 +1.8 -1.84 1.11 +0.38 -0.43 2.02 +0.35 -0.38	0.29 +0.1 -0.07 0.34 +0.16 -0.13 0.24 +0.13 -0.08 0.28 +0.13 -0.06 2.12 +0.4 -0.93 0.27 +0.15 -0.09 0.96 +0.31 -0.23

Table 6 :

 6 Forecasting results for ΣSFH (total star formation history) with ILLUSTRISTNG and EAGLE as source domains and SIMBA as the target domain. Predictions for all samples in SIMBA have been added, and 16 th , 50 th , and 84 th quantile values are drawn from the 50 predictions per galaxy, corresponding to the 50 neural networks used. Metrics are scaled (MAE × 1000, RMSE × 1000, BE × 1000, DTW × 1000 for readability.

	RMSE	MAE	BE	DTW	TDI
	Base KLIEP MDD DANN DeepCORAL 0.5 +0.04 0.3 +0.03 -0.01 0.27 +0.01 -0.01 0.3 +0.02 -0.01 0.3 +0.02 -0.01 -0.03 KMM 0.33 +0.05 -0.03 CORAL 0.55 +0.25 -0.09	0.23 +0.02 -0.01 0.2 +0.01 -0.01 0.23 +0.02 -0.01 0.22 +0.02 -0.01 0.41 +0.04 -0.03 0.25 +0.04 -0.03 0.45 +0.21 -0.07	0.04 +0.05 -0.04 0.02 +0.03 -0.02 0.02 +0.04 -0.03 0.02 +0.05 -0.05 0.01 +0.08 -0.05 0.05 +0.06 -0.04 0.32 +0.12 -0.12	1.0 +0.08 -0.03 0.94 +0.02 -0.02 0.98 +0.08 -0.04 1.01 +0.08 -0.07 2.08 +0.25 -0.16 1.1 +0.13 -0.08 2.08 +1.3 -0.48	5.13 +0.87 -0.89 3.19 +0.82 -0.41 4.92 +0.85 -0.58 4.99 +0.72 -0.9 6.74 +0.68 -0.34 4.83 +0.67 -0.75 6.49 +0.88 -1.24

Table 7 :

 7 Forecasting results, with SIMBA and EAGLE as source domains and ILLUSTRISTNG as the target domain. Predictions for all samples in ILLUSTRISTNG have been averaged, and 16 th , 50 th , and 84 th quantile values are drawn from the 50 predictions per galaxy, corresponding to the 50 neural networks used.

	RMSE	MAE	BE	DTW	TDI
	Base KLIEP MDD DANN DeepCORAL 0.68 +0.44 0.72 +0.34 -0.13 0.63 +0.2 -0.15 0.79 +0.36 -0.2 0.76 +0.27 -0.19 -0.17 KMM 0.85 +0.57 -0.32 CORAL 3.85 +1.92 -0.8	0.61 +0.26 -0.17 0.49 +0.15 -0.13 0.64 +0.26 -0.17 0.61 +0.23 -0.12 0.58 +0.46 -0.2 0.67 +0.43 -0.27 3.23 +1.44 -0.74	0.34 +0.49 -0.34 0.21 +0.25 -0.18 0.24 +0.42 -0.3 0.21 +0.43 -0.51 0.06 +0.8 -0.51 0.44 +0.55 -0.37 3.04 +1.11 -1.14	2.64 +1.37 -1.08 2.09 +1.17 -0.64 2.51 +1.55 -0.67 2.57 +1.01 -0.77 1.93 +1.67 -0.66 2.87 +2.89 -1.34 18.75 +12.24 -6.1	0.17 +0.37 -0.08 0.12 +0.08 -0.05 0.15 +0.24 -0.06 0.19 +0.4 -0.09 0.38 +0.46 -0.28 0.16 +0.17 -0.08 0.4 +3.29 -0.38

Table 8 :

 8 Forecasting results for ΣSFH (total star formation history) with SIMBA and EAGLE as source domains and ILLUS-TRISTNG as the target domain. Predictions for all samples in ILLUSTRISTNG have been added, and 16 th , 50 th , and 84 th quantile values are drawn from the 50 predictions per galaxy, corresponding to the 50 neural networks used. Metrics are scaled (MAE × 1000, RMSE × 1000, BE × 1000, DTW × 1000 for readability.

The observed CSFRD (especially at z > 2) is still a topic of active research currently, as the amount of star formation obscured by dust in the early Universe is still unknown. source-and target-domain data and results.


C Pardanaud
email: cedric.pardanaud@univ-amu.fr

M I Rusu

C Martin

G Giacometti

P Roubin

Y Ferro

A Allouche

M Oberkofler

M Koeppen

T Dittmar

M Köppen

Ch Linsmeier

Hydrogen retention in beryllium: concentration effect and nanocrystalline growth

Keywords: Hydrogen storage, plasma-wall interaction, Raman microscopy, beryllium hydride, crystalline growth

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The problematic of hydrogen isotope retention in beryllium addresses main issues related to plasma wall interactions in tokamaks [START_REF] Neu | Preparing the scientific basis for an all metal ITER[END_REF] and to hydrogen storage [START_REF] Sakintuna | Metal hydride materials for solid hydrogen storage: A review[END_REF].

Plasma wall interactions (erosion, retention, material modification, dust formation,...) occurring in tokamaks like the future ITER reactor (devoted to demonstrate the feasibility of obtaining energy from a magnetically confined D+T plasma) can lead to tritium retention, which is a major safety issue, due to its radioactivity [START_REF] Roth | Plasma-wall interaction: Important ion induced surface processes and strategy of the EU Task Force[END_REF][START_REF] Roth | Itpa Sol, Recent analysis of key plasma wall interactions issues for ITER[END_REF]. Hydrogen isotope retention in beryllium can also play a role in the life time estimation of the ≈ 700 m 2 ITER's beryllium inner walls because of erosion [START_REF] Borodin | Determination of Be sputtering yields from spectroscopic observations at the JET ITER-like wall based on three-dimensional ERO modelling[END_REF]. This explains why several studies based on ion irradiation of Be samples have been devoted to D behavior to ensure that tritium retention will not be a limiting issue in ITER operations. For deuterium kinetic energies in the range 15-60 eV, a D/Be ratio up to 0.7 has been reached in forming co-deposited Be layers and by varying substrate temperature and deposition rate [START_REF] De Temmerman | An empirical scaling for deuterium retention in co-deposited beryllium layers[END_REF][START_REF] De Temmerman | Revised scaling equation for the prediction of tritium retention in beryllium co-deposited layers[END_REF]. Older data also focusing on co-deposits reported values in the range D/Be=0.1-0.4 [START_REF] Mayer | Codeposition of hydrogen with beryllium, carbon and tungsten[END_REF][START_REF] Causey | Codeposition of deuterium with beryllium[END_REF][START_REF] Baldwin | Composition and hydrogen isotope retention analysis of co-deposited C/Be layers[END_REF]. Values in the range D/Be=0.2-0.4 have been measured under ion implantation of 20 keV/D with 10 19 D cm -2 [START_REF] Zhao | Dynamical saturated concentration of deuterium in a beryllium foil studied by low energy D(d,p)T reaction[END_REF]. For 9 keV/D implantation, a storage of D as D 2 molecules has been reported in bubbles tending to interconnect in the material [START_REF] Alimov | Depth distribution of deuterium atoms and molecules in beryllium implanted with D ions[END_REF]. Studies of 0.6 and 1 keV/D ions impinging on polycrystalline beryllium samples have shown that at fluencies ≈ 2 10 17 cm -2 , there is a saturation of the D content with a D/Be ratio in the range 0.2-0.3 [START_REF] Oberkofler | Deuterium release from implanted beryllium and beryllium oxide[END_REF][START_REF] Reinelt | Retention mechanisms and binding states of deuterium implanted into beryllium[END_REF], close to values reported in [START_REF] Haasz | Deuterium retention in beryllium, molybdenum and tungsten at high fluences[END_REF]. The diversity of results given above indicates impinging energies are not the only parameter driving the D/Be ratio. The thermal desorption spectroscopy of such samples displayed a single D 2 peak at 900 K for fluences lower than 0.7 10 17 cm -2 . A second peak was shown to appear at 750 K at a fluence of ≈1 10 17 cm -2 .

Two additional peaks appeared between 450 and 500 K for fluences higher than 1.2 10 17 cm -2 , a structural modification being involved to explain their existence. A comparative study of thermal desorption of D implanted in Be samples with two different crystalline orientation and in a polycrystalline Be sample has been done, revealing the role of both crystal orientation and D diffusion along grain boundaries under thermal annealing [START_REF] Piechoczek | Deuterium trapping and release in Be(0001), Be(11-20) and polycrystalline beryllium[END_REF]. It should be noted that none of all these studies have reported the formation of beryllium hydride in the form BeH 2 , except in [START_REF] Reinelt | Retention mechanisms and binding states of deuterium implanted into beryllium[END_REF] where it was suggested.

To act efficiently as an energy carrier, hydrogen should be absorbed and desorbed easily in specific materials, developed to fit with these two constraints [START_REF] Fichtner | Nanotechnological aspects in materials for hydrogen storage[END_REF][START_REF] Vajeeston | Predicting New Materials for Hydrogen Storage Application[END_REF]. MgH 2 is the prototype of the metallic hydrides having an interesting gravimetric capacity (7.6 wt.%), but its decomposition temperature is too high compared to others [START_REF] Van Den Berg | Materials for hydrogen storage: current research trends and perspectives[END_REF]. Beryllium hydride, BeH 2 being one of them, has been envisaged very recently, from a theoretical point of view, as a promising and more efficient hydrogen storage media than MgH 2 [START_REF] Zdetsis | Ab initio theoretical investigation of beryllium and beryllium hydride nanoparticles and nanocrystals with implications for the corresponding infinite systems[END_REF]. It has been synthesized in the past in a crystalline form [START_REF] Brendel | CRYSTALLINE BERYLLIUM HYDRIDE[END_REF][START_REF] Senin | STRUCTURE, AND PROPERTIES OF BERYLLIUM HYDRIDE[END_REF], in an amorphous form [23], and transitions between crystalline and amorphous forms have been studied [START_REF] Maienschein | THERMAL CRYSTALLINE-TO-AMORPHOUS TRANSITION OF BERYLLIUM HYDRIDE[END_REF].

In this study we report on the formation of BeD 2 crystalline structures (i.e. D/Be=2) by using plasma irradiation, and on their characterization by using nuclear reaction analysis (NRA), atomic force microscopy (AFM), optical microscopy, and quantum calculations compared to Raman microscopy. It was shown previously that this latter technique is sensitive to the way hydrogen isotopes are bonded and organized in materials and thus that it can be relevant for both fusion [START_REF] Pardanaud | Raman microscopy as a defect microprobe for hydrogen bonding characterization in materials used in fusion applications[END_REF] and hydrogen storage [START_REF] Reed | Recent applications of Raman spectroscopy to the study of complex hydrides for hydrogen storage[END_REF] applications. This work is organized as follows: experimental and computational details are given in Section 2. Results obtained from the computational study are presented in Section 3. Results obtained from experimental measurements are presented in Section 4. A discussion highlighting the consistency of all these results is presented in Section 5. Conclusions are given in Section 6.

Experimental and computational details

Ion implantation and sample characterization

The experimental study focuses on a 1.5 cm diameter polycrystalline beryllium sample that has been bombarded by a relatively homogeneous beam of 6 keV D 3+ (2 keV per deuterium) at a fluence of 1.9 10 17 D cm -2 . At the center of the sample, roughly 30% of the samples surface contains 1.0 10 17 deuterium atoms. On the other parts of the sample, the number of implanted deuterium atoms decreases rapidly down to ≈0 at the edges. Sample implantation and in situ ion-beam analysis are performed in the ultra high vacuum experiment ARTOSS [START_REF] Linsmeier | ARTOSS -A new surface science experiment to study the hydrogen inventory in multi-component materials[END_REF]. Detailed information on experimental setup, sample and procedure are published elsewhere [START_REF] Oberkofler | Properties of nitrogen-implanted beryllium and its interaction with energetic deuterium[END_REF]. In brief, a polycrystalline Be sample is cleaned by cyclic sputtering with Ar + ions and subsequent annealing to 1000 K. This Raman spectra were recorded on several locations of the sample, using a Horiba-Jobin-Yvon HR LabRAM apparatus (with x100 objective, with a numerical aperture of 0.9) in the backscattering geometry. λ L = 514.5 and 632.8 nm were used. The laser power was kept at ∼ 1mW µm -2 to prevent damages. Spectra were recorded with various exposure times to check that samples do not evolve under laser irradiation. We also used, when indicated in the text, the Raman imaging mode. It consists in moving the sample under the laser beam by 0.5 µm steps, and recording spectra at all the locations probed. Optical images obtained with the 100x objective used for making Raman spectra were recorded on several locations and some images were thresholded and binarized to determine the surface coverage of structures detailed below. Atomic force microscopy (AFM) was performed on some zones of the sample using the tapping mode of an AFM-NTMDT solver. The tip radius was smaller than 10 nm, the vertical and horizontal resolutions were ≈ 1 nm and ≈ 10 nm, respectively.

Computational details

The calculations were performed within the framework of the spin-polarized gradient-corrected density functional theory (DFT) using the Perdew-Burke-Ernzerhof functionals (PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF]) and including a semi-empirical pairwise London damped dispersion correction (PBE-D2 [START_REF] Grimme | Accurate description of van der Waals complexes by density functional theory including empirical corrections[END_REF][START_REF] Meijer | A density-functional study of the intermolecular interactions of benzene[END_REF]). A plane-wave basis set was used with kinetic energy cutoff for wavefunctions of 70 Rydberg (952 eV). The calculations were carried out using the Quantum-Espresso suite [START_REF] Giannozzi | QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials[END_REF]. The Vanderbilt ultrasoft pseudo potentials are taken from the GBRV database [START_REF] Garrity | Pseudopotentials for high-throughput DFT calculations[END_REF], all the electrons (core + valence) are included in the DFT calculation i.e. one e -for H and 4 e -for Be, this justifies the large energy cutoff. Within this methodology, the Be-H bond length of the BeH 2 molecule is calculated to 1.340Å, in good agreement with experimental data (1.332 to 1.342 Å [START_REF] Shayesteh | Infrared emission spectra of BeH2 and BeD2[END_REF]); the vibrational frequencies are ν 1 =707.8, ν 2 =1990.6 and ν 3 =2225.9 cm -1 , in very good agreement with results based on many body Coupled Cluster CCSD(T) / ccpV6Z calculations [START_REF] Koput | Ab initio prediction of the potential energy surface and vibrationrotation energy levels of BeH2[END_REF], ν 1 =713.5, ν 2 =1985.2 cm -1 and ν 3 =2171.5 cm -1 and experimental results ν 3 =2178.9 cm -1 [START_REF] Shayesteh | Infrared emission spectra of BeH2 and BeD2[END_REF]. These frequencies are plotted in Figure 1.The height of the peak is proportional to the number of modes and the peaks are convoluted by 20 cm -1 half-width The calculated lattice parameter of the hexagonal-compact unit-cell containing two beryllium atoms is a = 2.28 . It is in very good agreement with the experimental one 2.29 Å, and in better agreement than previous DFT-GGA results [36] [37]. This improvement is the consequence of adding London dispersion forces; it also improves the bulk modulus determined to 120 GPa (exp: 110 ~127 GPa). The phonon spectrum at the Γ point of the Brillouin zone is made of three frequencies 469.2, 469.2 and 723.9 cm -1 , with deviations lower than 7% as compared with experimental results [START_REF] Stedman | PHONON SPECTRUM OF BERYLLIUM AT 80 K[END_REF].

Studying trapping of deuterium atoms in beryllium requires supercells to avoid interactions between the repeated images from cells to cells. The 2x2x2 supercell (16 Be atoms) is the compromise we choose between the size of the system and the computational cost required.

Computational studies

We first calculated the phonon spectra of a BeD 2 molecule in a beryllium crystal. BeD 2 was built by adding two deuterium atoms to a beryllium of the crystal network while keeping the cell parameters unchanged. Two configurations were considered for BeD 2 : they are shown in inset of Figure 1 and will be referred as linear and bent.

The linear configuration is the most stable one. More precisely, the linear BeD 2 geometry is such that θ D-Be_D = 169° and the two non-equivalent Be-D distances are 1.44 Å and 1.51 Å. The phonon spectrum of the linear BeD 2 is shown in Figure 1; it is markedly different from that of the molecule in the gas phase. The bent configuration is 0.1eV less stable than the linear one. The geometry is such that θ D-Be_D = 71° and the non-equivalent Be-D distances are 1.51 Å and 1.61 Å. The phonon spectrum is shown in Figure 1. It's specific feature is to exhibit two low frequencies at 627 and 653 cm -1 , which would make the bent configuration in principle identifiable by Raman measurements.

The lattice system of crystalline beryllium hydride is orthorhombic. The beryllium atoms are in a tetrahedral-like symmetry; their environment is described in Figure 3 where angles and bond-lengths distributions are plotted. Figure 3 Finally, and in order to rationalized phonon spectra, we increased the D concentration starting from one BeD 2 linear molecule in the 2×2×2 supercell (D/Be=0.125) to the full D/Be=2 ratio maintaining the hexagonal symmetry and keeping cell parameters frozen. As the D concentration increases, the angle and bond-lengths distributions come closer and closer to that of beryllium hydride (Figure 3)

while the linear geometry is lost. This result tends to show that beryllium adopts a tetrahedral-like environment as the D/Be ratio increases. The calculated phonon spectra (Figure 2) indicate that increasing the D concentration produces the progressive enlargement of the spectrum with a shift toward the high frequencies, while the low frequencies remains more or less unchanged. However, the symmetry of crystalline beryllium hydride is orthorhombic, and the inter-beryllium atoms distance is different as compared to the hexagonal previous case. Such a difference induces a high stress in the hexagonal structure. We therefore relaxed the cell parameters of the hexagonal structure for two D/Be ratios (0.5 and 1) within the standard thermodynamic conditions (P=1 bar).

This relaxation has the effect to make appear low frequency in the phonon spectrum as can be seen in Figure 4. As a conclusion, high frequencies without low frequencies in spectra must be considered as the signature of highly constrained structures. On the contrary, low frequencies along with the full spectrum we calculated make the crystalline beryllium hydride identifiable by Raman spectroscopy.

Experimental studies

The influence of the number of implanted D atoms on the pristine material can be seen in Figure 5-a, which displays the spectral window where the beryllium stretching mode lies. The band, at 459 cm -1 for the pristine sample, downshifts and broadens when N D increases. This is usually interpreted as due to the creation of defects in the material induced by ion implantation. However, as it is known that Be-D bonds can be formed easily when D is implanted into a beryllium crystal, and that this results in a volume increase per beryllium atom [START_REF] Allouche | Quantum Modeling of Hydrogen Retention in Beryllium Bulk and Vacancies[END_REF], one can conclude, according to [START_REF] Gouadec | Raman spectroscopy of nanomaterials: How spectra relate to disorder, particle size and mechanical properties[END_REF],

that this downshift is induced by tensile stress in the material, exactly what is expected theoretically (see part III of this work). Note that a study of the effect of external compressive stress on a beryllium crystal has reported the upshift and the broadening of the band related to the beryllium stretching mode with the pressure increase [START_REF] Evans | X-ray diffraction and Raman studies of beryllium: Static and elastic properties at high pressures[END_REF]. Figure 5-b displays the 200-1800 cm -1 spectral window, with a zoom on the vertical axis. Contrary to what is expected by the Raman selection rules of a perfect Be crystal, which predict only the beryllium stretching mode to be active (found at 459 cm -1 in our study), the spectrum of the pristine sample is not flat elsewhere, especially in the 750-1100 cm -1 region where large bumps are clearly seen. We do not want to interpret in details their origin here but we can notice that it is not related to the presence of D in the material, as they are also displayed in the D free spectrum, nor to the presence of a surface oxide layer because it does not fit with the 338, 684, 1082 cm -1 frequencies reported for BeO [START_REF] Morell | Anharmonic interactions in beryllium oxide[END_REF]. When increasing N D , several bands appear at 255, 278, 617, 907, 1397, ≈1420, 1696, 1724 cm -1 . Most of these bands have their wavenumbers close to the one calculated for the BeD 2 crystal. It proves first that these bands are due to a Be-D bond, second: the D/Be ratio should be close to 2, and third: the material probed is crystalline, as the experimental and calculated spectra are very similar. Another argument suggests that it is crystalline: the full-widths-at-half-maximum (FWHM) of the Raman bands are close to 10 cm -1 , contrary to what is generally measured for amorphous beryllium hydride (generally hundreds of cm -1 ) [START_REF] Sampath | Vibrational dynamics of amorphous beryllium hydride and lithium beryllium hydrides[END_REF]. The intensities of these bands seem to be related to the amount of D trapped in the material as they increase with N D .

Having a crystalline structure subsequently to a high fluence implantation is counterintuitive and needs to have a closer insight. optical microscope images that the dendrites coexist with circular domes, the latter being the dark spots visible in Figure 6. From that we can conclude that domes appear at lower deuterium N D than dendrites. Dendrite heights are in the range 50-200 nm, whereas domes are smaller. These two kinds of structures are clearly distinguishable in Figure 7-a, where dendrites appear dense and bulky and domes are regularly distributed on the surface. These two kinds of structures are not so clearly distinguishable in Figure 7-b where dendrites appear sparsely connected, and domes are not distributed regularly, as if they were constituting units of the nascent dendrites.

Figure 8 shows a closer look on the domes (Figure 8-a) and dendrites (Figure 8-b) by means of AFM. The domes are smoother than the dendrites, the latter can even be composed of two collinear branches. This observation suggests that dendrites are more complex objects than domes.

Figure 8-c and 8-d display some dome and dendrite profiles, respectively. The larger heights can be as high as 75 and 160 nm for domes and dendrites, respectively. The dome average profiles exhibit a parabolic average shape with a 50 nm height and 600 nm bases. The average profiles of the dendrites are more rectangular with a 90 nm height and a 1 µm base. The rectangularity is just the consequence of summing up single and double dendrite branches, and do not reflect the shape of single dendrite profiles, but it is useful for mean height estimation on all the sample.

Figure 9 shows Raman imaging in the vicinity of a dendrite, in the zone of the sample where N D =3 10 17 D cm -2 . In Figures 9-c and9-d the intensities of the band close to 459 cm -1 due to the Be stretching mode, and the band at 1397 cm -1 due to a Be-D mode are displayed, respectively. These two bands are more intense on the dendrite than off the dendrite. In Figure 9-e, we compare the spectra on and off dendrite. These two spectra have been obtained by summing all the spectra over the two areas and normalizing to the band at 459 cm -1 , the threshold distinguishing between on and off dendrite being fixed at a value of 100 in the Figure 9-c map. The 5 bands identified in Figure 5-b (and pointed by a star in Figure 9) as due to Be-D modes in a BeD 2 crystal are much more intense on dendrites than off dendrites. The spatial correlation between the dendrites and the BeD 2 Raman signature allows us to conclude that the dendrites, formed under ion bombardment, are a high reservoir of trapped deuterium. From 0 to 1.5 10 17 D cm -2 , ν decreases linearly by ≈ 3 cm -1 and the FWHM increases linearly by ≈ 6 cm - 1

, for both laser wavelengths used. In this interval, the dendrite coverage is close to zero, meaning that most of D atoms are implanted in a ≈ 50 nm depth layer, according to the implantation depth profile calculated using SRIM [START_REF] Ziegler | The Stopping and Range of Ions in Matter[END_REF] (given in the insert of Figure 10), a free access Monte Carlo computer program based on the binary collision approximation. For N D = 1.5 10 17 D cm -2 , and taking into account the ≈ 50 nm of implantation depth, it gives a value D/Be ≈ 0.24, which is compatible with what was previously measured [START_REF] Reinelt | Retention mechanisms and binding states of deuterium implanted into beryllium[END_REF]. For N D > 1.5 10 17 D cm -2 , ν and FWHM both reach a plateau.

The dendrite coverage starts to increase linearly from ≈ 0 at N D ≈ 2 10 17 D cm -2 and reaches up to 10% close to 3 10 17 D cm -2 , without tendency for saturation in our experimental conditions. Note that the values of ν and FWHM recorded using the 632.8 nm laser are sligthly lower than the one recorded using the 514.5 nm laser. This may be due to a resonance effect, as suggested previously in [START_REF] Evans | X-ray diffraction and Raman studies of beryllium: Static and elastic properties at high pressures[END_REF] and references therein, to explain the slight discrepencies found in the literature.

Discussion

The comparison of vibrational spectra obtained from quantum calculations of amorphous and crystalline Be including D, and Raman spectra obtained for the dendrites formed on the D implanted Be sample indicates that the D/Be ratio in dendrites is close to 2 and the structure is probably crystalline and not amorphous. On the zone of the sample where 10% of the surface is covered by dendrites, the number of D retained in such structures is N = (3 10 17 -2 10 17 D cm -2 )= 1 10 17 D cm -2 . Knowing that dendrites are 90 nm in height and cover 10 % of the surface, the density, expressed in g cm -3 , is:

ρ=(N × (2 + 9/x)) / (5.4×10 15 )
where x=D/Be. Assuming a density of 1.23 g cm -3 , obtained from quantum calculations for the BeD 2 crystal (see Table 1), x is found equal to 2.1, confirming in very good agreement the assumption of the crystalline structure suggested by the comparison of calculated and measured spectra.

Knowing first that dendrites appear when N D is higher than 2 10 17 D cm -2 , and second that the 450 and 500 K TDS peaks of [START_REF] Oberkofler | Deuterium release from implanted beryllium and beryllium oxide[END_REF][START_REF] Reinelt | Retention mechanisms and binding states of deuterium implanted into beryllium[END_REF] (with samples implanted in similar conditions than ours) which appear for fluences higher than 1.2 10 17 D cm -2 , we propose that the two peaks result from the thermal degradation of dendrites.

The formation of circular domes and dendrites is reported here, but not fully explained. We only suggest that increasing tensile stress in the material (suggested by both Raman measurements and quantum calculations) due to volume increase subsequently to Be-D bond creation may play a role in the formation of the BeD 2 crystal. It is also evident that the impinging ion kinetic energy, tuning the implantation depth in the material, is a key parameter in creating this stress, similarly to what was found on bombarded graphite: the height of nanometer sized domes was observe to increase when a higher amount of underlying amorphous carbon (generating stress) was created, using higher implantation energies [START_REF] Pardanaud | In-plane and out-of-plane defects of graphite bombarded by H, D and He investigated by atomic force and Raman microscopies[END_REF]. In the literature, other structures have been obtained on Be samples after irradiation: cone-like shapes under plasma exposures in the PISCES-B plasma facility [START_REF] Doerner | Impact of surface morphology on sputtering during high-fluence plasma exposure[END_REF], interconnected labyrinths with hydrogen bubbles at higher energies [START_REF] Chernikov | Gas-induced swelling of beryllium implanted with deuterium ions[END_REF]. In ITER, only taking into account the charge exchange process leads to anticipate at least a homogeneous 4 × 10 17 cm -2 s -1 flux bombarding the Be walls [START_REF] Carpentier | Modelling of beryllium erosion-redeposition on ITER first wall panels[END_REF]. Extrapolating our result to this flux, we guess that the initial Be surface state in ITER will be completely transformed in a rough surface in only few seconds.

These surface modifications under bombardment can significantly influence the local sputtering yield [START_REF] Doerner | Erosion of beryllium under high-flux plasma impact[END_REF], which is of importance for extrapolating the lifetime of the future ITER reactor beryllium walls, future work should focus on the understanding of formation mechanisms leading to these structures (and for example explaining the route leading from dome formation to dendrite formation).

However, to our knowledge the observation of such dendrites was not yet reported in the JET tokamak with ITER-like wall configuration where such fluences can be found.

The eventual role of particle (sputtered Be, hydrogen,...) mobility on the surface may also be considered for explaining these dendrites formation in the framework of the diffusion limited aggregation model [START_REF] Witten | Diffusion-limited aggregation[END_REF]. Stability of the sample after being placed in air for several months (time spent between the different analyses) should also be addressed.

In parallel, but for hydrogen storage purposes, ion assisted techniques should be used, systematically varying the relevant parameters, to optimize the volume of BeD 2 formed, to increase the possibility of storing hydrogen isotope in beryllium.

Conclusion

After deuterium ion irradiation of a polycrystalline beryllium target (2 keV/D), we have shown that deuterium is trapped up to a saturation value of D/Be≈0.24 when N D =1.5 10 17 D cm -2 , with an increasing tensile stress developing inside the material, and saturating at this N D value. For higher N D values, we have observed that a higher amount of D can be trapped in dendritic structures that start to grow at N D =2.0 10 17 D cm -2 and cover for example 10% of the surface when N D =3 10 17 D cm -2 . By using a multi-technique approach (nuclear reaction analysis, atomic force microscopy, optical microscopy, and quantum calculations compared together with Raman microscopy), we conclude that these dendrites are crystalline and composed of BeD 2 . These structures can be suitable for storing hydrogen, even if a safety issue concerning the use of Be needs to be performed. As they might play a role in hydrogen isotope retention, future experiments should be devoted to search for such dendrites in the JET tokamak, or in the future ITER tokamak. )

  sample is exposed to a mass separated beam of D 3 + ions at 2 keV/D with a fluence of 1.9*10 17 D cm -2 in the central area of 5*8 mm 2 . The amount of deuterium trapped is measured on several locations of the sample surface by means of the D( 3 He,p)He nuclear reaction with a lateral resolution on the order of 1 mm. Raman spectra were also recorded in the vicinity of some of these regions. The numbers given in that figure, and referred to as N D later in the text and figures, are the absolute numbers of implanted D atoms measured by means of this nuclear reaction. For comparison with results from the literature we give N D in deuterium per cm 2 .

  Gaussian functions. Figures 2 to 4 are plotted the same way. As a consequence, these Figures cannot be directly compared to Raman spectrum , except for the bottom of figure 2 where selection rule are also taken into account.

Figure 6

 6 shows a survey of the surface morphology at different locations of the sample, obtained by means of optical microscopy. In that figure we display some images showing that in zones where N D <1.5 10 17 D cm -2 , some dark spots can be evidenced, their concentration depending on N D . In zones where N D >1.5 10 17 D cm -2 dendrites appear together with the already existing dark spots. The dendrite shape, contrast, size and concentration vary depending on N D , too. For example, low contrast dendrites in Figure6-d seems to be composed of disconnected 1-2 µm size branches, whereas high contrast dendrites in Figure6-g are composed of 5 (or more) µm size interconnected branches. Figures 7-a and 7-b show height measurements of these structures, such as those found in Figures 6-d and 6-g, by means of AFM. It appears more clearly than in the

Figure 10 displays

 10 Figure 10 displays the evolution with N D of two spectroscopic parameters related to the Be stretching mode (wavenumber, ν, and FWHM) and compared to the dendrite surface coverage (obtained by thresholding + binarizing ≈ 4500 µm 2 optical images, as shown in the Figure 10-c insert).

Figure 1 .

 1 Figure 1. Calculated vibrational frequencies obtained in dressing the numerical values by a Gaussian function of the BeD 2 molecule in beryllium and in gas phase (same convention in Figures2 and 4).

Figure 2 .Figure 3 .Figure 4 . 6 Figure 5 .

 23465 Figure 2. Vibrational frequencies associated to deuterium trapping in the beryllium crystal as a function of the deuterium concentration (D/Be atomic ratio), only the modes involving D are represented. In the lower part of the figure the calculated full beryllium hydride crystal Raman active frequencies are shown.

Figure 6 .

 6 Figure 6. Optical microscope images of various regions on a deuterium-implanted beryllium sample. (a) Pristine. (b) and (c) In zones of the sample where N D <2.0 10 17 D cm -2 . (d) to (g) In zones of the sample where N D >2.0 10 17 D cm -2 .

Figure 7 .

 7 Figure 7. Height of the deuterium implanted sample in zones containing N D >2.0 10 17 cm -2 , obtained by AFM. (a) and (b) represent typical morphologies found in that zone.

Figure 8 .

 8 Figure 8. Height profile analysis of the domes and dendrites. (a) and (b) are AFM images representing

Figure 9 .Figure 10 .

 910 Figure 9. Raman imaging in the vicinity of a dendrite. (a, b) optical microscope images. (c) Intensity of

  reveals a single bond length at 1.41Å and two specific angles at 96° and 116° not so far from the 109° of a perfect tetrahedron. The phonon spectrum of the crystalline Beryllium hydride (BeD 2 ) is shown in the lower part of Figure2. It is computed taking into account selection rules for Raman spectroscopy. As a consequence, this spectrum is directly comparable with

Raman measurements. It is characterized by low frequencies components appearing below 500cm -1 , which again would make beryllium hydride experimentally identifiable.

Table 1 .

 1 Volume and density evolution as a function of deuterium concentration at the pressure of 1 atmosphere.

	D/Be		Density
			(g.cm -3 )
		0	1.95
	Amorphous	0.5 1	1.77 1.58
		2	1.33
	Crystal	2	1.23

Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement number 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The authors were granted access to the HPC resources of [CCRT/CINES/IDRIS] under the allocation i2014080509 made by GENCI (Grand Equipement National de Calcul Intensif).