Costanza Simoncini
email: costanza.simoncini@univ-rennes1.fr

Yan Rolland

Valery Morgenthaler

Hervé Saint-Jalmes Krzysztof Jurczuk

Pierre-Antoine Eliat

Marek Kretowski

Johanne Bezy-Wendling

Hervé Saint-Jalmes

Blood Flow Simulation in Patient-Specific Segmented Hepatic Arterial Tree

Keywords: Computational fluid dynamics, hepatic artery, image processing, liver tumor, phase contrast MRI, radioembolization

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

Introduction

Hepato-cellular carcinoma (HCC) is the fourth cause of mortality in the world. Treatment options are limited, less than 15% of patients are candidates for surgery and half of them for chemotherapy or radiotherapy [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF]. A viable option to avoid side effects is the development of local treatments such as Selective Internal Radiation Therapy (SIRT) [START_REF] Grouphealth | Clinical Review Criteria. SIRT (Selective Internal Radiation Therapy) Therasphere and SIR Sphere for Unresectable Hepatocellular Carcinoma[END_REF], which brings millions of radioactive microspheres directly into arterial tumor blood supply through a catheter, inducing far less embolizing syndrome than its equivalent in chemotherapy, and comparable (even better in an end-stage HCC) benefits [START_REF] Andreana | Treatment of hepatocellular carcinoma (HCC) by intra-arterial infusion of radio-emitter compounds: trans-arterial radio-embolisation of HCC[END_REF]. The hepatic vascular supply is double: 25% comes from the hepatic artery and 75% from the portal vein. HCC is mainly vascularized by the hepatic artery. When the vessels irrigating the tumor are embolized, the arterial vascular supply of healthy liver can be also strongly reduced if the tumor targeting is not optimized. Therefore, if the portal supply does not compensate enough, post-embolisation syndrome can be severe.

Administration of SIRT requires a very precise knowledge of the hepatic arterial tree: currently, there is no available tool providing physicians the position where the catheter should be placed in the arterial tree, for the optimization of tumor targeting. In order to identify the best injection point, the physician realizes a cone-beam CT angiography to visualize the arteries geometry. When a therapy delivery point is chosen, biodegradable mock-spheres (Technetium 99m Tc albumin macro-aggregated, MAA) are injected and a scintigraphy acquired one hour later reveals their distribution. The distribution of radioactive microspheres is thus approximately predicted from only one preliminary injection of 99m Tc MAA, and for a same injection site, their behavior can actually be sometimes different [START_REF] Garin | Utility of Quantitative Tc-MAA SPECT/CT for yttrium-Labelled Microsphere Treatment Planning: Calculating Vascularized Hepatic Volume and Dosimetric Approach[END_REF]. Then, the simulation of several microspheres injection scenarii is of great interest for the clinician.

In the current work we aim at modeling arterial blood flow in a patient-specific hepatic arterial tree. This is a necessary step towards the patient-specific simulation of the distribution of the injected microspheres transported by blood. If the arteries segmentation is fine enough to detect the principal vessels irrigating the tumor, and if the simulation is initialized with patient-specific values of pressure or velocity, then the concentration of microspheres in the tumor can be more accurately estimated by existing appropriate CFD (computational fluid dynamics) methods.

It has been shown [START_REF] Basciano | Computer modeling of controlled microsphere release and targeting in a representative hepatic artery system[END_REF] that the distribution of microspheres at a bifurcation is not simply proportional to the blood flow distribution in the two descendant vessels, and that it can significantly depend on the vessels geometry, which impacts on the flow streamlines and consequently on the particles trajectories. Blood flow characteristics like pressure and velocity can be deduced from angiography and through CFD if and only if we provide appropriate BC (boundary conditions) like pressure and/or velocity at the inlet and outlets of the vascular tree. This is why, to simulate the flow, we plan to combine these anatomic information with hemodynamics information extracted from PC MRI.

Of course, the vasculature and perfusion in tumors and healthy tissue is a key point to estimate the final microspheres (and thus dose) distributions. Nevertheless, the first step is to have the best estimation of velocity and pressure in the vessels upstream from such vasculature. Clearly the fluid and microspheres propagation in small vessels is far more complex to analyze, and for this reason we cannot expect to have access to a precise vessels network shape: based on angiographic data we can only expect, at this level, to segment the tumor regions, and this is not in the scope of this paper. Let us only indicate here that coupling real data and simulated data based on image characteristics could help solving this scale limitation, that could be achieved using a computational model of the finer vascular network, like the one we previously developed [START_REF] Kretowski | Physiologically based modeling of 3-D vascular networks and CT scan angiography[END_REF][START_REF] Jurczuk | Silico Modeling of Magnetic Resonance Flow Imaging in Complex Vascular Networks[END_REF].

A complete model of SIRT will let the radiologist simulate and quantify the distribution of microspheres in tumor and in healthy tissue, starting from any potential injection points and with different injection modalities, including dose amount and injection site and velocity. Vessels enhancement from 3D angiography will be described in Section 2.1. Next, the creation of a triangular mesh for numerical simulation on the extracted geometry and the simulation setup are described in Section 2.2. Section 3 presents the principal results we obtained concerning blood flow and the preliminary validations we performed. As shown in Section 4, the procedure we developed opens the way to a patient-specific simulation of microspheres transport by blood flow, able to predict their distribution in liver depending on injection modalities and hemodynamic characteristics.

Method

Vessels extraction is performed on cone-beam CT 3D angiography. Data sets are up to 512 × 512 × 374 voxels, and voxel size is 0.463 × 0.463 × 0.463 mm. Images are acquired a few dozens of seconds after the injection of iodine-based contrast media, during the so-called arterial phase, i.e. when contrast agent is only in the arteries and is still not spread in all hepatic parenchyma. This acquisition gives a detailed portrait of the arterial tree until vessels with a diameter of around 0.5 mm.

All procedures involving human participants were in accordance with the ethical standards of the institutional research committee and with the 1964 Declaration of Helsinki. Approval by the institutional Medical Ethics Review Board and proper informed consent were obtained.

Patient-specific vasculature extraction

Arteries segmentation is performed by a first automatic enhancement of vessel structures using a 3D multiscale Frangi's algorithm [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF]. Methods based on the analysis of the Hessian matrix are indeed shown to provide good segmentations results in most vascular cases, and the multiscale character of Frangi's filter is particularly suitable for hepatic vasculature, whose range of diameters is relatively large [START_REF] Jimenez-Carretero | 3D Frangi-based lung vessel enhancement filter penalizing airways[END_REF]. The algorithm is based on the idea that a given voxel belongs to a vessel if the local second derivative of the contrast function is low in a direction (i.e. along the direction of the vessel) and high in the two orthogonal directions (i.e. towards the walls of the vessel). This property can be exploited through an analysis of the eigenvalues (λ 1 , λ 2 , and λ 3 ) of an estimation of the Hessian matrix. In every voxel such matrix is obtained from second order finite differences, computed on the grey level image previously smoothed with a 3D Gaussian kernel. The size of this kernel is a scale parameter: a larger scale parameter corresponds to a more efficient smoothing. In particular, three geometrical characteristics are calculated from the eigenvalues in order to determine the tubular-like character of a given structure.

R A = |λ 2 | |λ 3 | indicates if it is closer to a plate-like or a line structure. R B = |λ 1 | √ |λ 2 λ 3 |
reflects its divergence from a blob-like form, and S = 3 j=1 λ 2 j is used to distinguish vessels information from background noisy structures (the background noise becomes, after smoothing, a low and approximately constant grey level image, leading to small second order differences values). Frangi's filter computes the product of 3 monotonic nonlinear functions of, respectively R A , R B , and S. The sensitivities of the filter to R A , R B , and S are tuned separately by three parameters. The efficient size of the smoothing parameter, i.e. of the scale parameter, depends on the vessel diameter. Therefore two other parameters have to be defined in this method, in order to take into account the size variability of the vessels: a scale range (depending on the diameters of the smallest and the largest vessel visible in the acquisition), and the scale ratio giving the step used to go through this range. We identified a set of parameters convenient for cone-beam CT angiography, given in Table 1. Vessels segmentation is then accomplished by a thresholding of the Frangi's enhanced image, followed by a partitioning of the resulting image into connected regions with a union-find algorithm (function

Blood flow simulation

We achieved patient-specific blood flow simulation in ANSYS Fluent (ANSYS, Inc., Canonsburg, PA). The conversion from an irregular segmented binary volume to a more regular meshed geometry, from which stable CFD simulations can be performed, is not straightforward and will be described below. Next, blood flow simulation setup will be outlined.

Mesh generation

In order to implement blood flow simulation with a CFD software, a closed smooth mesh is needed, where every surface corresponding to an input or output of blood flow must be separately defined. The binary volume resulting from the Frangi-based segmentation is thus used to create a triangular mesh corresponding to the surface of the segmented vessels. Input and outputs surfaces selection is achieved by a manual cut of the mesh at every vessel extremity. The result is a mesh describing a surface which is open at the entrance of blood and at every vessel end. A number of surfaces, called cap, corresponding to the number of mesh holes, needs to be created in order to perfectly close the open mesh. Such caps can be created thanks to the properties of the mesh: indeed, mesh consists of only triangles. When the mesh is clipped by a plane, new mesh nodes are created on the intersection between the mesh triangles and the plane, transforming some of the triangles into quadrilaterals. We developed an algorithm that first looks for all the mesh nodes which belong to a triangle or a quadrilateral, and which belong to only three facets or less. Such property guarantees that they are on the cut border of the mesh. All the identified border points then need to be distinguished between every cap. A k-means clustering algorithm [START_REF] Arthur | K-means++: The Advantages of Careful Seeding[END_REF] accomplishes this step minimizing the average squared Euclidean distance into every group. Cap surfaces are created by a simple preliminary triangularization: for every set of nodes the barycenter is calculated, points are then classified in a trigonometric order, and triangles are finally defined by linking every point to its following neighbor and to the barycenter. Finally, quadrilaterals resulting from mesh clipping are split in two triangles in order to redefine a proper triangular mesh.

Afterwards, the surface corresponding to vessels wall and those corresponding to caps are imported in ANSYS Fluent meshing. A remeshing is performed using wrapping in order to obtain a smooth closed surface and mesh quality compatible with fluid simulations. Triangles minimum and maximum sizes can be adapted to the segmented vessels in order to smooth the surface or to better describe bifurcations and narrow vessels by choosing, respectively, a larger or smaller triangle size.

Different parameters for minimum and maximum sizes of the mesh elements were tested. A minimum size of 0.3 mm and a maximum one of 2.5 mm were chosen and give a good, smooth, representation of the geometry. In order to capture the behavior of the boundary layer with sufficient precision and reduce numerical instabilities, multiple layers of prisms are capping the artery walls. A volume mesh is generated on the basis of the surface mesh and passed to the ANSYS Fluent Solver, based on the finite volumes method.

Simulation setup

Blood is preliminarily considered as an incompressible Newtonian fluid. Its simulated viscosity value is defined as 0.004 Pa s, and its density as 1060 Kg /m 3 [START_REF] Levick | An Introduction to Cardiovascular Physiology[END_REF]. Inlet condition at the proper hepatic artery is given at mean arterial blood velocity, which is defined to be 30 cm/s over the entire vessel section. Different values are available in the literature [START_REF] Dyvorne | Quantification of hepatic blood flow using a high-resolution phase-contrast MRI sequence with compressed sensing acceleration[END_REF][START_REF] Tziafalia | Echo-Doppler measurements of portal vein and hepatic artery in asymptomatic patients with hepatitis B virus and healthy adults[END_REF][START_REF] Hubner | Hepatic arterial blood flow velocities: assessment by transcutaneous and intravascular Doppler sonography[END_REF][START_REF] Kito | Doppler sonography of hepatic arterial blood flow velocity after percutaneous transhepatic portal vein embolization[END_REF]. In such studies blood flow is measured in either the common, left or right hepatic artery, and the average velocity varies between 11 and 50 cm/s.

Due to the lack of an agreement in the literature concerning blood velocity values in hepatic vasculature, we measured it in the proper hepatic artery using a phase contrast (PC) MR sequence realized on a 1.5 T MR system (Aera, Siemens healthcare, Erlangen, Germany). Due to the retrospective nature of this study, it was not possible to acquire PC MRI on the same patients whose angiographies were analyzed, but the sequence was performed on a healthy volunteer in order to acquire in vivo velocity values. Phase contrast MRI has been used and improved since the eighties for blood flow measurements [START_REF] Bryant | Measurement of flow with NMR imaging using a gradient pulse and phase difference technique[END_REF][START_REF] Dyvorne | Quantification of hepatic blood flow using a high-resolution phase-contrast MRI sequence with compressed sensing acceleration[END_REF][START_REF] Bollache | Numerical modeling of arterial pulse wave propagation to characterize aortic hemodynamic: Validation using magnetic resonance data[END_REF]. It is a gradient-echo sequence presenting bipolar gradients aiming to encode blood velocity in any given direction. Thus, image contrast is proportional to blood velocity in the chosen spatial direction. In this study a 3D PC sequence was used, providing a time-averaged velocity in the three spatial directions without requiring too long acquisition times, prohibitive due to clinical constraints. The parameters of the sequence we acquired are given in Table 2. Velocity encoding value, corresponding to the maximum velocity that is to be acquired, was chosen as 80 cm/s in order to cover all the possible spectrum of mean physiological velocities in abdomen. The acquisition was performed the axial plane and velocity in each of the three directions was encoded. The acquisition on the volunteer shows a velocity of 15 cm/s, which is Pressure-Linked Equations) scheme [START_REF] Karki | Pressure based calculation procedure for viscous flows at all speedsin arbitrary configurations[END_REF]. The convergence is verified by checking the mass and momentum conservation in between inlet and outlets. Calculation are stopped when a difference of 10 -7 for continuity and 10 -4 for momentum is reached. To enhance convergence, a Least-Square Gradient Method was chosen. Such method is comparable to the one of the node-based gradient on skewed and distorted unstructured meshes and is less expensive to compute [START_REF] Rausch | Spatial Adaptation of Unstructured Meshes for Unsteady Aerodynamic Flow Computations[END_REF].

Results

The procedure was applied to the cone-beam angiographies of 4 patients suffering from an advanced stage HCC, as can be seen in Figures 1 and2.

Extraction of patient's data

Cone-beam angiographies of 2 of the 4 patients are shown in Fig. 1a and Fig. 2a Min. face area (mm 2 ) 2.5 e -3 1.1 e -3 1.9 e -3 1.1 e -3

Max. face area (mm 2 ) 0.56 0.298 0.589 0.47

Flow simulation results

Blood flow simulation was performed on a portable computer equipped with the quad-core Intel Core i7-2820QM CPU (2.3 GHz) and 16 GB of RAM. Computational time of simulation varied between 15 and 45 minutes on 3 cores. Table 4 shows area-weighted average pressure and velocity values for all the 4 patients, as well as segmented vessels radii. Data are given, for each patient, for the blood inlet surface (proper hepatic artery) and for the outlets average. and 2. Velocity results for the same patients are shown in Fig. 1c and2c, where blood pathlines are colored based on velocity magnitude.

Validation

Pressure and velocity values at the entrance of the proper hepatic artery are coherent with literature [START_REF] Ficher | The Liver[END_REF], whereas blood flow characteristics farther in the hepatic arterial tree (i.e. for smaller in current study can vary approximately between 70 and 100 mmHg [START_REF] Bura-Rivière | Physiologie et exploration de la circulation artérielle[END_REF]. Information about blood velocity at small hepatic arteries is also hard to determine exactly. According to [START_REF] Bura-Rivière | Physiologie et exploration de la circulation artérielle[END_REF], mean blood velocity for arteries whose diameter varies between 2 and 6 mm, can vary between 20 and 50 cm/s.

To address this hindrance, we performed preliminary MRI phase contrast acquisition on healthy volunteer which gave the results shown in Table 5 concerning mean blood velocity until a vessel diameter of 2.2 mm.

In order to justify the choice of the output boundary conditions of the simulations, the ratios between diameter and blood flow in vessels of similar sizes were also observed. It is well known that vessels of different sizes (aorta, arteries, arterioles,...) have different blood flow behaviors, this is why only vessels with similar sizes were compared. Even if the velocity measurements on PC MRI were possible only on vessels much larger that the outputs of the simulation, we verified that our choice of boundary conditions was reasonable. The ratios between diameter and blood flow are similar in all considered vessels, except in the left hepatic artery. Nevertheless, the diameter of this vessel is smaller and, consequently, the measurement less reliable.

Discussion and conclusion

We developed a procedure to extract a patient-specific hepatic arterial tree from cone-beam CT 3D angiography, until a vessel diameter of around 0.5 mm. We realized blood flow simulation, the results of which appear to be coherent with the literature, even if flow precise characteristics at such level are very little known and vary significantly between different individuals. Our MRI measurements are not statistically representative but are validated by an experienced radiologist and by known mean velocities at the aorta and inferior vena cava levels [START_REF] Gabe | Measurement of Instantaneous Blood Flow Velocity and Pressure in Conscious Man with a Catheter-Tip Velocity Probe[END_REF]. Their interest is double: information about blood velocity at smaller arteries is crucial for the definition of boundary conditions, by letting us tuning output conditions that are currently little known, in particular tumors blood flow. Output conditions would be adjusted in order to obey to some data fidelity relative to velocity measured on PC MRI. Moreover, it will let us validate blood flow results.

Different approaches, not considered in this paper, exist for the estimation of boundary con- We plan to undertake a clinical study on blood flow measurements, optimizing phase contrast MRI sequence. A better image resolution can be achieved by synchronizing the acquisition to breathing and by choosing a velocity encoding value tuned to the patient's peculiarities (obesity, cirrhosis, HCC can have an influence on blood flow). Phase contrast sequence is not currently included in the clinical protocol for HCC, but such clinical study will allow us to apply the described procedure to patients who have undertaken both cone-beam CT angiography and MRI phase contrast, which was not possible for the current study due to its retrospective nature. Moreover, new 4D PC MRI techniques, not yet available on the MR system used for this study, should allow time-resolved acquisitions of three dimensional PC MRI, synchronized to the cardiac cycle [START_REF] Stankovic | 4D flow imaging with MRI[END_REF]. It should be noticed though, that time constraints in the clinical workflow do not actually encourage this approach.

A more precise description of variable blood viscosity and pulsatility will also be performed.

Moreover, before SIRT, patients can undertake chemoembolization, vessels coiling or clipping whose metal can induce artifacts in acquisitions and reduce the quality of the segmentation, as for Patient 1, where the diameters of tumor vessels near the clips are overestimated. Results are inevitably dependent on segmentation, which is currently validated by placing a mask on the original image, and whose tuning is still operator-dependent. In particular, the choice of the low threshold needed before the classification into connected components influences the final segmentation especially concerning smaller vessels, whose diameter may be overestimated. In the future an adaptive thresholding method like the Otsu's method or the Kuwahara filter will be used in order to extract more realistic vessels diameters and too complex tumor vessels may not be segmented but simulated by the vascular model [START_REF] Kretowski | Physiologically based modeling of 3-D vascular networks and CT scan angiography[END_REF][START_REF] Jurczuk | Silico Modeling of Magnetic Resonance Flow Imaging in Complex Vascular Networks[END_REF] instead.

Hepatic artery diameters are, for the four patients, larger than the average, which is assessed at 3.8 ± 0.8 in [START_REF] Tziafalia | Echo-Doppler measurements of portal vein and hepatic artery in asymptomatic patients with hepatitis B virus and healthy adults[END_REF] for healthy patients. This can be explained by the presence of cirrhosis and HCC, which can induce a more important blood flow in liver arteries and, in turn, a larger diameter of the proper hepatic artery. The patient-specific blood flow simulations described in this paper represent a preliminary, compulsory step that will allow us to to perform, under the same boundary conditions and taking into account every patient's peculiarities, an accurate numerical simulation of SIRT microspheres transport in the patient's hepatic arterial tree, thanks to the Fluent particles transport tools. In particular, we envisage a computation of particles trajectories based on a Lagrangian particles tracking approach, that would take into account all the forces acting on the microspheres as fluid drag force and gravity. On the other side, the influence of microspheres on blood flow can be neglected due to the large ratio between the diameters of vessels and microspheres and the low microspheres density. Such simulations have a crucial importance for physicians in the aim of tumor targeting.

Figure 1 :Figure 2 :

 12 Figure 1: Cone-beam angiography and blood flow simulation results on segmented arterial tree for Patient 1

  ditions. For example, Aramburu et al. [2] introduce a perfusion simplified 2 parameters model in order to complete a 0-D model version of the flow in the imaged tree. The parameters identification requires the knowledge of synchronised values of pressure and flow rate at the inlet and at different phases of the cardiac cycle (systolic and diastolic phases).

Table 1 :

 1 Frangi's filter parameters for vessels enhancement Matlab, where voxels are considered to belong to the same connected components if they are linked by at list one of their 26 neighbors). In order to reduce image noise, the classification into connected components requires a very low threshold, whose value is currently manually chosen based on the enhanced image histogram. The output of the this step corresponds to a 3D image where every connected component is labeled by a different grey level. The connected component corresponding to the arterial tree is then visually identified and extracted. This allows us to delete every enhanced structure that does not belong to the hepatic artery.

	Scale range Scale ratio α	β	c
	[1, 5]	0.5	0.5 0.9 500
	"bwlabeln" in

Table 2 :

 2 Parameters of the 3D PC MRI sequence acquired for blood quantification in the hepatic artery

	Acquisition type	3D
	Slice thickness	1.3 mm
	Pixel size	1.6 mm × 1.6 mm
	Repetition time	51.75 ms
	Echo time	4.56 ms
	Flip angle	15
	Acquisition plane	axial
	Velocity encoding	80 cm/s

coherent with measurements in

[START_REF] Dyvorne | Quantification of hepatic blood flow using a high-resolution phase-contrast MRI sequence with compressed sensing acceleration[END_REF]

. Output boundary conditions are given on blood flow percentage, as to be proportional to every outlet vessel diameter. Classical no-slip conditions are enforced on vessel walls.

Navier-Stokes Equations are finally solved in ANSYS Fluent. Due to the relatively low velocities, a pressure-based solver was used together with a segregated SIMPLE (Semi-Implicit Method for

  through a Volume Rendering technique. Patient's tumor and the vessels that irrigate it are clearly visible.It can be noticed in Fig.1that the vessels supplying blood to the tumor are more numerous and more irregular, due to tumor angiogenesis. Moreover, metallic clips near the patient's tumor make the image particularly noisy and complicate the segmentation of the tumor vessels.Vessels segmentation results in arterial trees with 20 to 50 outlet vessels. Radius of the proper hepatic artery varies between 2.3 and 3.6 mm. The smallest radius of segmented vessels is 0.49 mm. Details about the obtained mesh are given in Table3for all the patients.

Table 3 :

 3 Mesh characteristics for the hepatic arteries of the four patients

		Patient 1 Patient 2 Patient 3 Patient 4
	Number of outlets	24	20	29	50
	Number of faces	2133552	968791	1353653	3995926
	Number of nodes	288175	130792	187330	540732

Table 4 :

 4 Radius, area-weighted average pressure and velocity results in segmented arterial trees of the 4 patients

	Patient 1 Patient 2 Patient 3 Patient 4

Table 5 :

 5 Blood flow velocity measured on phase contrast MRI acquisition on a 27 years old healthy volunteer

		Diameter Mean blood velocity Ratio between diameter
				and flow ratio
	Abdominal aorta	1.5 cm	21 cm/s	-
	Inferior vena cava	1.1 cm	19.4 cm/s	-
	Celiac trunk	7 mm	35 cm/s	-
	Common hepatic artery	3.1 mm	15.5 cm/s	44
	Gastroduodenal artery	3 mm	12 cm/s	59
	Right hepatic artery	2.8 mm	12.9 cm/s	59
	Left hepatic artery	2.2 mm	9.8 cm/s	98
	Splenic artery	3 mm	13.7 cm/s	51

Acknowledgment

The authors would like to thank Jean-Jacques Bellanger, Jean-Claude Nunes and Michel Rochette for their precious help in this work.

Funding: This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Disclosure: Each author has materially participated in the research and article preparation.

All authors have approved the final article. The authors have no conflict of interest. V. Morgenthaler is employed by ANSYS.


J Byggmästar
email: jesper.byggmastar@helsinki.fi

E A Hodille

Y Ferro

K Nordlund

Analytical bond order potential for simulations of BeO 1D and 2D nanostructures and plasma-surface interactions

Keywords: interatomic potential, beryllium, oxygen, molecular dynamics

An analytical interatomic bond order potential for the Be-O system is presented. The potential is fitted and compared to a large database of bulk BeO and point defect properties obtained using density functional theory. Its main applications include simulations of plasma-surface interactions involving oxygen or oxide layers on beryllium, as well as simulations of BeO nanotubes and nanosheets. We apply the potential in a study of oxygen irradiation of Be surfaces, and observe the early stages of an oxide layer forming on the Be surface. Predicted thermal and elastic properties of BeO nanotubes and nanosheets are simulated and compared with published ab initio data.

Introduction

The plasma-facing materials in future fusion reactors must be able to withstand extreme conditions due to high thermal loads, ion fluxes, and neutron irradiation. Beryllium has been selected as the first-wall material of ITER. The properties contributing to the choice of Be include a low atomic mass (and hence low plasma contamination), low tritium retention, and its oxygen gettering ability [1]. The behaviour of Be under fusion-relevant conditions has therefore been extensively studied in experimental and modelling work [2,3,4,5]. However, little is known about the effects of oxygen on the structural and erosion properties of the Be walls [6,7]. The high affinity for oxygen is likely to lead to a growing plasmafacing oxide layer if exposed to air, with subsequent erosion and sputtering of oxide molecules. Addressing the consequences of the otherwise desirable oxygen gettering ability of Be is therefore crucial.

Beryllium oxide has also recently gained considerable interest within the nanomaterials community, due to the partly covalent bonding nature of the otherwise ionic compound. Ever since the discovery of single layered carbon nanostructures, nanotubes and nanosheets of carbon and other materials have been intensively studied for their exceptional electronic and mechanical properties [8,9]. The possibility of producing carbonlike sp 2 -bonded BeO nanostructures has been recently proposed [10], and been the subject of a number of first-principles studies [10,11,12,13,14,15,16,17].

Molecular dynamics (MD) simulations are a useful tool for studying atomic processes on time and length scales inaccessible to experiments and density functional theory, such as sputtering mechanisms or radiation damage production. The accuracy of MD simulations is directly defined by the interatomic potential used to describe the interactions between individual atoms. Developing accurate many-body interatomic potentials for different materials is therefore an essential part of atomistic modelling, and no many-body potential has to our knowledge been developed for the Be-O system. A many-body potential for the Be-O system will allow simulations of both fusion-relevant plasma-surface interactions with oxide layer formation and erosion as well as simulations of mechanical and thermal stability of BeO nanostructures.

The potential formalism used here was originally developed by Tersoff [18,19], using the concept of bond order as discussed by Abell [20]. Despite being originally developed to model covalently bonded materials, the potential has been shown to closely resemble the embedded atom method potentials used for metals [21,22]. The bond order formalism therefore allows modelling of differently bonded materials, from pure metals to e.g. metal and semiconductor compounds. The potential function used in this work, referred to as the analytical bond order potential (ABOP), is similar to the original form by Tersoff, but slightly extended and rewritten. The potential is capable of modelling bond breaking and forming, and has been previously applied to a wide variety of materials, including pure metallic, metal-carbonhydrogen, and metal-oxygen systems [23,24,25,26,27,28,29]. In this work, we extend previously developed potentials for the fusion-relevant Be-W-C-H materials [25,28,30,31] to include the Be-O interaction.

The article is structured as follows. In section 2, we briefly describe the functional form of the ABOP, and the strategy used for fitting the potential to the constructed database of material properties. In section 3, we summarise the accuracy of the fitted properties of the potential, and discuss results from testing of the potential by calculating properties not included in the fitting process. Finally, we apply the ABOP in simulations of O irradiation of a Be surface, as well as investigations of the thermal stability and elasticity of BeO nanotubes and nanosheets. Concluding remarks are given in section 4.

Methods

Potential formalism

The functional form and its parameters has been extensively discussed previously (see e.g. Ref. [23]), and will only be briefly presented here. The total potential energy in the ABOP is expressed as a sum over all atomic bonds, given by (1) where the repulsive and attractive functions form a simple Morse-like potential [32] V

V = i j>i V ij = i j>i f C (r ij )[V R (r ij )-b ij V A (r ij )],
R (r ij ) = D 0 S -1 exp -β √ 2S(r ij -r 0 ) (2)
V A (r ij ) = SD 0 S -1 exp -β 2/S(r ijr 0 ) .

D 0 and r 0 are the dimer bond energy and length. The interaction range is typically restricted to the nearestneighbour shell by the cutoff function

f C (r) =        1, r ≤ R -D 1 2 - 1 2 sin π 2D (r -R) , |R -r| ≤ D 0, r ≥ R + D. (4)
Here, R is the cutoff distance and D defines the width of the cutoff region, where the potential energy smoothly approaches zero. All many-body interactions are incorporated in the bond order function, which is written in the symmetric form

b ij = b ij + b ji 2 , (5)
where

b ij = (1 + χ ij ) -1/2 . (6)
Furthermore, χ ij is given by

χ ij = k( =i,j) f C (r ik )g ik (θ ijk )ω ijk exp [α ijk (r ij -r ik )] , (7)
and the angular function by

g ik (θ ijk ) = γ ik 1 + c 2 ik d 2 ik - c 2 ik d 2 ik + (h ik + cos θ ijk ) 2 . (8)
In simulations involving high kinetic energies, such as irradiation studies, the repulsive part of the potential given above is unphysically weak and must be modified to describe repulsive short-range interactions more accurately. A common approach is to join the universal repulsive Ziegler-Biersack-Littmark potential [START_REF] Ziegler | The Stopping and Range of Ions in Matter Treatise on Heavy-Ion Science[END_REF], V ZBL (r ij ), with the original potential V ij as

V ij = F (r ij )V ij + [1 -F (r ij )] V ZBL . (9)
Here, F (r) is the Fermi function

F (r) = 1 1 + exp[-b f (r -r f )] , (10)
which ensures a smooth transition between the original potential and the repulsive ZBL potential. The parameters b f and r f , defining the transition range and distance, are chosen so that the equilibrium properties are nearly unchanged for all considered coordination numbers.

Computational details

The interatomic potential was fitted using the code tulip [START_REF] Henriksson | TULIP: Program for fitting an interatomic potential to data[END_REF]. Molecular dynamics simulations were carried out with the codes parcas [START_REF] Nordlund | [END_REF]36] and lammps [37]. The time step required for energy conservation in the MD simulations was 0.6 fs, tested at room temperature. The irradiation simulations employ an adaptive time step according to the selection criterion in [38]. Structural optimisations, defect energies, and the phonon dispersion in the ABOP were calculated using ase [39] and lammps [37]. All density functional theory (DFT) calculations were carried out using the Quantum Espresso code [40].

Fitting database

Beryllium oxide crystallises in the wurtzite structure at atmospheric pressure and room temperature [41,42]. A large amount of experimental data is available for characterising BeO wurtzite: the lattice constants [43,44,45,46], the cohesive energy [START_REF]CRC Handbook of Chemistry and Physics[END_REF], the bulk modulus [43], and the elastic constants [43,[START_REF] Bentle | [END_REF]49,50]. However, some other hypothetical Be-O structures in which the beryllium atom has different coordination numbers have to also be considered. These structures will enlarge the fitting database and ensure a better transferability of the interatomic potential to any kind of structural environment.

Consequently, DFT calculations were carried out for the Be-O dimer, a linear chain of Be-O units, a graphene-like 2D structure in which Be-O units form hexagons in a plane, and the 3D wurtzite, zinc blende, rock salt, and ceasium chloride structures. The corresponding cohesive energies, lattice constants, bulk moduli and elastic constants were determined. In addition, the formation energies of different defects in BeO wurtzite were calculated as well as the solution energies of oxygen in pure hcp Be. Surface properties were also considered by investigating the adsorption of an oxygen on the Be hcp (0001) surface. Additionally, phonon properties of BeO wurtzite were determined to assess the reliability of the ABOP.

The energy of the oxygen molecule in its ground state is required to compute the cohesive, formation, and adsorption energies discussed in this work. The ground state of the O 2 molecule is spin triplet, which consequently requires spin-polarised calculations. In addition, point defects such as the Be vacancy in BeO wurtzite displayed magnetic properties, as previously pointed out in the literature [51,52,53], which also requires spin-polarisation.

The DFT calculations have been performed with the GGA-PBE exchange and correlation functional [54] and ultra-soft pseudo-potentials [55]. Cutoff energies of 52 Ry (707 eV) and 416 Ry (5658 eV) were used for the truncation of the expansion of the wave functions and the electronic density, respectively. The k-point sampling has been converged to 24 × 24 × 24 for the wurtzite unit cell containing 2 Be atoms and 2 O atoms. The phonon calculations were performed in the framework of the Density Functional Pertubation Theory [56] with the same cutoff energies and k-point sampling as above.

The k-point samplings for the others structures were scaled as the inverse of the size of the box. The dimer, linear, and graphene-like structures exhibit nonperiodic structures in at least one direction. Thus, in these directions, the size of the box was increased to be more than 20 Å to prevent interactions across the periodic boundaries. In these directions, only one kpoint was used.

The formation energies of defects in the BeO wurtzite lattice were calculated in a 3 × 3 × 3 supercell (54 Be atoms and 54 O atoms) with a k-point sampling of 6 × 6 × 6. Relaxations of both the atomic positions and the volume of the cell were carried out.

The calculations of an O impurity in the Be hcp lattice and on the Be hcp (0001) surface were performed with a cutoff energy of 50 Ry (680 eV) for the wave functions and 400 Ry (5440 eV) for the electronic density. The formation energy of O was calculated in a 5 × 5 × 4 Be supercell (200 Be atoms). For these calculations, the k-point sampling was 4 × 4 × 5 (20 × 20 × 20 for the Be hcp unit cell with 2 Be atoms). The positions of all atoms as well as the volume of the cell were relaxed. For the surface calculations, a 3 × 3 supercell was used in the a and b directions (parallel to the (0001) surface). Eleven hexagonal Be layer were considered in the c direction (perpendicular to the surface) with a 24 Å vacuum layer to prevent interactions between the studied surface and the bottom of the repeated supercell in the c direction. The k-point sampling used for these calculations was 6 × 6 × 1. The calculations were done by relaxing all atoms except the two bottom layers.

Fitting methods

With reference data for BeO phases in a wide range of coordination numbers from the DFT calculations, the fitting of the potential parameters followed a similar strategy as previously outlined in e.g. Ref. [23]. The cutoff range of Be-O interactions was chosen to be between the first and second nearest neighbour Be-O shells for all structures included in the fitting process. When the interaction range is restricted to nearest neighbours only, the energy per bond, E b , follows the Pauling relation [START_REF] Pauling | The Nature of the Chemical Bond and the Structure of Molecules and Crystals: An Introduction to Modern Structural Chemistry[END_REF], given in the ABOP by

E b = D 0 exp -β √ 2S(r b -r 0 ) , (11)
where r b is the bond length. D 0 and r 0 are typically chosen to give the correct dimer energy and bond length, β can be calculated from the ground state frequency of the dimer, and S fitted to give the correct bond energies for the higher coordinated structures. However, the bond energy of the dimer given by DFT and experiments is weak compared to the bond energies of the higher coordinated structures, and could not be used to fix D 0 . Both D 0 and r 0 were therefore included, together with S and β, when fitting equation 11 to satisfy the bond energies of all structures studied by DFT. Care had to be taken so that the interaction range of the underlying Be-Be potential (version II from Ref. [28]) did not interfere with the above fitting strategy. Including both nearest neighbour Be-O and Be-Be interactions effectively changes the above Pauling relation for Be-O structures. The nearest neighbour Be-Be bond in the ground state of BeO (wurtzite), is 2.7 Å according to experimental data [43]. The cutoff radius of the Be-Be potential is 2.685 Å, resulting in Be-Be bonds being included after only a small compressive strain. Consequently, fitting the ABOP to the experimental lattice constants leads to a poor description of elastic properties, as the energy cost for straining the structure is highly asymmetrical around the equilibrium.

In order to avoid such unphysical behaviour, the ABOP was fitted to the (slightly overestimated) lattice constants and cohesive energies obtained in the DFT calculations.

The parameters of the angular function g ik (θ ijk ) in equation 8 were numerically fitted using the reference data obtained by DFT. The parameter h, defining the optimal angle between the atoms ijk, was allowed to vary between the equilibrium angles of the lowest energy phases, wurtzite (109.5 • ; h = 0.33), and the graphene-like sheet (120 • ; h = 0.5).

Results and discussion

Fitted properties

Table 1 shows the properties of different BeO phases given by the ABOP and compared to data obtained by the DFT calculations and available experimental data. The potential parameters are given in table 2. The structural properties of the ground state (wurtzite) phase are well reproduced by both DFT and the ABOP, including all elastic constants. However, it should be noted that the experimental values for the elastic constants vary significantly from study to study. The lattice constants obtained in the DFT calculations are slightly overestimated compared to the experimental data, which is a known feature of the GGA-PBE functionals [START_REF] Haas | [END_REF]. Nevertheless, as noted previously, the ABOP was fitted to reproduce the DFT lattice constants of the different phases, in order to get good elastic properties.

As a consequence of restricting the interactions to nearest neighbours, the ABOP formalism cannot reproduce the subtle energy difference between the wurtzite and zinc blende structures. The geometry of the nearest neighbour shell is identical in both structures, and the energy and bond distances will [43] b Calculated from the room temperature enthalpy of formation from Ref. [START_REF]CRC Handbook of Chemistry and Physics[END_REF] c Ref. [START_REF] Bentle | [END_REF] d Ref. [49] e Ref. [50] therefore be exactly the same in the ABOP. The cohesive energy and bond distance of the graphenelike BeO sheet are also in good agreement with DFT, and more importantly, the difference in cohesive energy between the wurtzite and graphene phases is in excellent agreement with DFT.

The overall good agreement between the ABOP and DFT across a wide range of coordination numbers generally ensures a good transferability to bonding geometries and coordination numbers not found in the fitted structures. In the following sections, we test this assumption by calculating properties not included in the fitting database.

Tested properties

BeO bulk properties.

When fitting an interatomic potential to properties of a select amount of crystalline phases, it is important to confirm that the Table 2: Parameters for the Be-O ABOP. The Be-Be and O-O parameters are taken from earlier studies (Be-Be version II from Ref. [28]), and the Be-O parameters are fitted in this study.

Be-Be [28] O Upon heating to temperatures above the melting point, we found that BeO melts into a very low-density liquid phase with large regions of empty space between the clusters of atoms. The melting point is commonly estimated by creating a liquid-solid interface, and determining the temperature at which the molten and crystalline phases are in equilibrium. However, due to the melting behaviour of BeO, this was not possible. We therefore carried out simple heating simulations with different heating rates to determine an approximate melting temperature. Simple heating of a single-crystalline system in MD is known to overestimate the melting point due to the the lack of nucleation points [59], and due to the extremely high heating rates required for MD time scales. The observed melting point for all tested heating rates was about 3000 K. This should be considered an upper limit predicted by the ABOP, and is therefore in good agreement with the experimental melting point of 2851 K [START_REF]CRC Handbook of Chemistry and Physics[END_REF].

We also determined the thermal expansion coefficients of wurtzite BeO and compared the obtained results to experimental data. The experimental data from Kozlovskii et al. [60] are for the temperature range 293 K-1823 K. Slack and Bartram [61] reviewed values in the temperature range 300 K-2200 K, which is the range on which we focused our study.

Both studies reported the linear expansion, i.e. the variation of the sample in a single direction and not the variation of the volume of the sample. Since the crystal we are simulating is a single crystal, unlike the experiments, we have easily access to the linear thermal expansion coefficients, α i∈{x,y,z} (in K -1 ), of the different directions of the lattice. The appropriate coefficient to compare to the experimental data on polycrystalline BeO is then the directionally averaged value [61]

α = 1 3 (α x + α y + α z ). ( 12
)
To obtain the thermal expansion, a cuboid box of wurtzite BeO, containing 15840 atoms was used. The temperature of the system was quickly increased from 0 K to 2000 K in 10 ps, kept constant for about 40 ps so that the volume of the box relaxed to zero pressure, and finally the temperature was decreased with a rate of 10 K/ps. The thermal expansion was determined by looking at the variation of the box lengths in the x, y and z directions that give the different α i∈{x,y,z} . The comparison with the experimental data is reported in figure 1. The simulation data with the ABOP is overall in good agreement with the experimental values, especially above 1000 K. At lower temperatures, the thermal expansion is slightly overestimated. In the range 1300-2200 K, the experimental values reported by Slack and Bartram seem to be constant. In the ABOP, α also has this behaviour, but for the range 1000-2000 K. Above 2000 K, the lengths of the box exhibits large changes that lead to a quick increase of the thermal expansion.

In the simulation, the values of α x and α y are the same, except at high temperatures. Indeed, wurtzite has hexagonal symmetry, and consequently the x and y directions are not symmetrically equivalent. There is also a clear anisotropy in the thermal expansion in the z direction compared to the two other directions, which is expected for a wurtzite structure.

Figure 2 shows the phonon dispersion of the wurtzite BeO phase obtained with the ABOP and compared to DFT and experimental results [50]. Our DFT results are in good agreement with the experimental measurements.

The lower acoustic branches are well reproduced in the ABOP. The subtle differences in the energies of the acoustic branches between the ABOP, DFT, and experiment are likely due to the corresponding differences in elastic constants, as seen in table 1. However, the optical branches at around 90 meV in DFT and experiments are strongly overestimated in the ABOP, with the corresponding branches at around 140-150 meV energies. The overestimation of the optical branches can be attributed to the lack of long-range

α x α y 1 3 • (α x + α y + α z ) α z

ABOP Exp Kozlovskii Exp Slack and Bartam

Figure 2:

Phonon dispersion of wurtzite BeO calculated with the ABOP and compared with DFT calculations and experimental data from Ref. [50]. coulombic interactions in the ABOP [START_REF] Ashcroft | Solid State Physics[END_REF], and was previously also observed in the ZnO ABOP [29]. 3 shows the bond lengths and cohesive energies of a few selected Be-O molecules predicted by the ABOP and compared with our DFT calculations and experimental results [START_REF] Chupka | [END_REF]. The properties of the dimer BeO are directly determined by the potential parameters D 0 and r 0 . As discussed earlier, choosing D 0 to reproduce the (BeO) 2 and (BeO) 3 rings were therefore relaxed in DFT and with the ABOP. The relaxed shape and bond angles are different in the DFT and ABOP simulations, but the bond lengths and energies are in reasonable agreement. As the number of atoms in the molecules is increased, the energy and bond length tend towards those of an infinite BeO chain, as seen in table 3. The qualitatively good agreement between the ABOP and DFT for molecules is promising for applying the ABOP in studies of surface irradiation, where sputtering of various molecules is expected.

Be-O molecules. Table

Interstitial oxygen atoms in Be.

Formation energies for oxygen in the common high-symmetry interstitial sites in hcp Be were calculated with the ABOP, and compared to DFT calculations. The formation energy is defined as

E f = E D -n Be µ Be -n O µ O , (13)
where E D is the total energy of the system containing the defect, and n Be (n O ) the number of Be (O) atoms.

The chemical potentials at 0 K are given by the cohesive energies as µ Be = -3.624 eV and µ O = E O2 /2 = -2.583 eV in the ABOP. Figure 3 shows an illustration of the different interstitial sites with the adopted naming and abbreviations. The calculated formation energies are given in table 4. The most stable interstitial position according to DFT is the the basal tetrahedral (BT) site. However, as can be seen in table 4, our DFT data are about 1.5 eV higher in energy than the DFT data by Middleburgh et al. [64] and Zhang et al. [65]. On the other hand, the relative stabilities between the BT, NBO and BO interstitial sites are in good agreement (within a margin of 0.06 eV) between the three DFT studies. Thus, the reason for the difference in formation energy is clearly that the reference energy for oxygen is not the same in all three studies. To investigate the effect of the oxygen reference, we calculated the energies of oxygen atoms and molecules in different states: the spin triplet (ground state) and a state with no spin. Using the spin triplet state, we obtained the formation energies reported in table 4. However, using the state with no spin, the formation energy of O in BT becomes -2.73 eV, in NBO -1.31 eV, and in BO -0.08 eV. These values are much closer to the values reported by Zhang et al. and Middleburgh et al., indicating that their values have been calculated using a state with no spin for the oxygen reference. Clearly, no spins exist in the ABOP formalism, and the only possible reference state of oxygen is the cohesive energy of the O 2 molecule given by the ABOP parametrisation.

In the ABOP, the most stable site is a position above the BT site, close to the NBT site (denoted r-NBT in table 4). However, an oxygen interstitial at the NBT site relaxes to the BT site, as opposed to the close-by ground state r-NBT position, in agreement with DFT observations [64,65]. The order of stability for the remaining interstitial sites predicted by the ABOP is in agreement with DFT results, although the differences in energy vary, and e.g. the stability of the NBO site is strongly overestimated.

Despite predicting the incorrect ground state oxygen interstitial site, the ABOP still qualitatively reproduces the diffusive properties obtained by DFT. The energy barrier for migration between adjacent BT sites, through an intermediate NBO site, is about 1.6 eV according to DFT calculations by Zhang et al. [65]. Although the lowest energy migration path in the ABOP is between the predicted ground state (r-NBT) and the close-by BT site along the [0001] axis (see figure 3), further migration through the bulk requires exiting the tetrahedral surrounding, following paths similar to those studied earlier in DFT. The migration energy of the BT-NBO-BT path is about 2.3 eV in the ABOP, higher than the 1.6 eV obtained by DFT [65]. In comparison, the vacancy migration in Be is about 1.3 eV in the ABOP, also higher than the DFT value of 0.72 eV reported in [64]. The self-interstitial migration energy in Be is 0.73 eV according to the ABOP [28], compared to 0.64 eV given by DFT [64]. Due to the relatively high O migration energy and the fact that the ABOP correctly predicts the O interstitial to be trapped in a tetrahedral atomic environment, even though the exact position inside the tetrahedron contradicts DFT results, can therefore be considered acceptable for most practical use in atomistic simulations. Oxygen is essentially immobile during MD time scales, which is important for oxidation and deposition simulations. Oxidation of Be surfaces is experimentally known to be limited by Be atoms migrating up through the oxide layer [66], as opposed to oxygen migrating into the Be bulk, which due to the high oxygen migration barriers is likely to be reproduced by the ABOP. We will return to this discussion in section 3.3.1, where we study oxide layer formation and growth using the ABOP.

3.2.4. Adatom energetics on Be. Additionally, we calculated adsorption energies of oxygen atoms on a Be (0001) surface with the ABOP, and compared to DFT calculations. The adsorption energy of a single oxygen atom is defined as

E ad = E D -E S -µ O , (14)
where E D is the total energy of the surface slab with the adsorbate, E S the total energy of the clean surface slab, and µ O the chemical potential of oxygen (-2.583 eV in the ABOP). Figure 4 shows the common adsorption sites on a hcp (0001) surface. The corresponding adsorption energies and the distances to the surface layer are given in table 5. There is a noticeable difference in the formation energies between DFT and the ABOP. However, due to the different possibilities of treating µ O in DFT, as discussed earlier, emphasis should be put on the order of stability when comparing the ABOP with DFT, and not the absolute values. The relaxed distances from the surface layer are in good agreement between DFT and the ABOP, except for the fcc site. The hcp and fcc sites have identical nearestneighbour symmetry, and therefore similar adsorption energies. However, the ABOP predicts the fcc site as the preferred adsorption site over the hcp, in contrast to the DFT results. This is a consequence of the interaction range of the ABOP and cannot be corrected without sacrificing the fitted bulk BeO properties. For the hcp site, the second-nearest Be neighbour (the atom directly below the hcp site in the first sub-surface layer) is barely within the cutoff radius of the ABOP after relaxation, while the same atom is outside of the cutoff sphere for the fcc site. This is also evident from the difference in the relaxed distances from the surface given in table 5. The functional form of the ABOP leads to a weakening of the nearest neighbour bonds due to the second nearest neighbour for the hcp site, resulting in a higher energy than for the fcc site.

Point defects in BeO.

The chemical potentials of the elements can vary depending on the chemical environment. The formation energy for a neutral point defect can then be expressed in the more general form used in Ref. [67], as where ∆µ is thermodynamically restricted to the range -∆H f < ∆µ < ∆H f , and ∆H f is the enthalpy of formation of BeO in its ground state phase (-4.91 eV/f.u. in the ABOP). The limits of ∆µ correspond to O-rich and Be-rich conditions. Here, we use ∆µ = 0 in all calculations. Formation energies of simple neutral point defects were calculated with DFT and compared to values predicted by the ABOP, as listed in table 6. The formation energies for single vacancies in BeO are strongly underestimated in the ABOP, particularly the Be vacancy. The stability of point defects in metal oxides, such as BeO or ZnO, are generally dependent on the electronic structure and charge state of the point defect [68,69,70], effects that the ABOP cannot capture. Good accuracy in describing point defects in BeO can therefore not be expected from the ABOP. In the case of the Be vacancy, our DFT calculation predicts a magnetisation of 2.2 µ B /cell induced by the Be vacancy, which could explain the large difference in the value obtained with the ABOP. Similarly, the Be antisite induces a 2.05 µ B /cell magnetisation, and consequently the difference in the formation energy predicted by the ABOP and DFT is large (around 5 eV for both the Be vacancy and antisite). Interestingly, the same 2-3 eV difference as for the O interstitial formation energies in pure Be (table 4), is observed between DFT and the ABOP for defects in BeO with an excess of one O atom (O antisite and O vacancy).

E f = E D - 1 2 (n Be + n O )µ BeO - 1 2 (n Be -n O )(µ Be -µ O ) - 1 2 (n Be -n O )∆µ, (15)
In all of these cases, DFT gives higher formation energies than the ABOP. This difference could be due to the different manner the O 2 molecule is considered, inducing a shift between the DFT and ABOP values, as discussed previously.

The lowest energy Be and O interstitial sites in BeO predicted by the ABOP were identified by creating 200 BeO systems with a randomly placed O or Be interstitial. All systems were relaxed at 1000 K and subsequently quenched to 0 K, allowing the interstitials to find a nearby energy minimum. The ground state O interstitial configuration in the ABOP (E f = 6.75 eV) was an oxygen atom in an octahedral site on the Be-basal plane, but with strong relaxation of the geometry of the surrounding atoms. This complex configuration is likely an artefact of the potential. However, only 4 out of the 200 random simulations ended in this configuration. The majority of cases produced an O-O dumbbell with E f = 8.08 eV or E f = 8.57 eV, depending on the orientation. The formation energy obtained by DFT for an O interstitial in a Be-basal octahedral, corresponding to the closest high-symmetry site to the ABOP ground state configuration, is 9.50 eV. The O-O dumbbell interstitial was found to be significantly more stable in DFT, with E f = 6.00 eV.

For the Be interstitial, the ABOP predicts a slightly off-plane O-basal octahedral site as the preferred interstitial position, with E f = 7.50 eV. The corresponding DFT calculation resulted in relaxation of the Be atom in the opposite direction, with a final formation energy of 9.72 eV. The second lowest energy site in ABOP was an off-axis tetrahedral site with a slightly higher formation energy of 7.73 eV.

Note that only values for the lowest energy sites are listed in table 6, and in the case of an O interstitial, the DFT and ABOP values do not correspond to the same interstitial configuration. Furthermore, while the preferred sites in the ABOP certainly correspond to the global minimum energy sites for interstitials, that is not necessarily true for the DFT values, as we only considered a few selected sites. The calculated DFT values only allow a qualitative justification of the magnitudes of interstitial formation energies in BeO predicted by the ABOP.

Applications

Oxygen irradiation of Be.

As an application of the Be-O ABOP, we studied the oxidation of beryllium by incident oxygen ions. Simulations of O ion irradiation at normal incidence on a (0001) Be surface were carried out. The simulations were performed using an initial cell containing 4704 atoms. The energies of the ions were 1, 10, and 100 eV, and the temperature of the Be target was 300 K, 1000 K,

M = 2.2 µ B /cell b M = 2.05 µ B /cell
and 1500 K. One thousand cumulative simulations at a flux of 2.6×10 28 m -2 s -1 were done for all energies and temperatures. Between every single bombardment, the cell was randomly shifted in the x and y directions in order to model a uniform distribution of the ions on the surface. Periodic boundaries were used in the x and y directions. The temperature of the simulation box was kept constant by controlling the temperature of the borders of the simulation cell. The two bottom layers were kept frozen to mimic an infinite lattice. 5 ps separated each bombardment. In ITER, the incident flux of plasma particles on the Be first wall is estimated to be about 10 19 -10 21 m -2 s -1 [71]. Assuming a high 10% O concentration in the plasma, one could expect an oxygen flux of around 10 18 -10 20 m -2 s -1 on the Be plasma-facing components. This is obviously a value much lower than the one used in the MD simulations, but given the short time scale available, this problem would only be resolved by increasing the surface by 10 orders of magnitude, which is not realistic at the MD length scale.

Snapshots of the irradiated (0001) Be surfaces at 300 K are shown in figure 5a-c. In the early stage of the irradiation, after around 10-100 deposited ions, the deposition region is still almost hexagonal Be. In the 100 eV irradiation case, the O atoms mainly stop in the r-NBT and NBO sites, which are the most stables sites according to the ABOP (table 4). In the lower energy cases, and especially in the 1 eV case, the ions do not have enough energy to overcome the energy barrier at the surface, and consequently remain on the surface mainly in the hcp and fcc hollow sites (table 5). After 1000 O ions, the region where the O ions stop is no longer a hexagonal Be structure, and an amorphous oxidised layer with a rough surface has been formed, as seen in figure 5. The thickness of the oxide layer decrease as the energy of the O ions decrease. For 100 eV ions, a large amount of O ions are implanted below the surface while for 1 eV ions, the O ions are almost exclusively stuck on top of the surface. The roughness of the oxide surface is most pronounced for the irradiation with 100 eV ions. These ions have enough energy to cause sputtering from the surface, and also induce damage deep below the surface (figure 5a). For 100 eV ions at 300 K, several O atoms are also seen much deeper than the oxidised layer. These deep atoms are results of channelling [72], where ions travel through the material in lines of low atomic density. Channelled ions are also observed in the irradiations at 1000 K and 1500 K, but closer to the surface. This could be explained by the increased thermal motion of the target atoms at the higher temperatures, effectively making the channelling paths more narrow.

To further analyse the process of oxidation in our MD simulations, atom profiles (Be and O), time evolutions of the limits of the oxide layer as well as potential energy per atom profiles were plotted as functions of depth. The limits of the oxide layer are constructed by looking at the depth of the deepest oxygen atom and the one closest to the vacuum limit. Thus, in the 100 eV case, channelled ions might be taken into account. These plots are presented in figure 6a-c for the three temperatures and the three energies. The initial position of the (0001) Be surface is indicated with a vertical line, in order to see where the oxidised layer is formed. The oxide layer and its growth can be divided into two zones: the one below the initial surface and the one above the initial surface. To explain the growth of the oxide layer toward the bulk, one could propose either implantation of O ions or diffusion of O through Be. Similarly, to explain the growth toward the vacuum, migration of Be through the oxide layer could be envisaged, which can be caused either by transfer of recoil energy or diffusion of Be through the layer. The latter process has been used to explain the growth mechanism of oxide layers in experiments by Roth et al. [66,6]. The following analyses of figure 6 is an attempt to discriminate which of these mechanisms limits the growth rate of the oxide layer.

For 100 eV ions (figure 6a), the oxide layer has a thickness of about 35 Å at all temperatures (excluding the channelled ions).

The layer is equivalently distributed below and above the initial Be surface. This evidences the fact that part of the oxidised layer is formed by implantation of the O ions. However, as can be seen in the Be profiles, some Be atoms are also extracted from below the surface to construct a sparse oxide layer above the initial Be surface. Looking at the time evolution of the oxide layer limit below the initial Be surface, one can see straight lines at the same depth for all temperatures within an approximation of a few Å. This depth corresponds to the maximum implantation depth of the O ions, meaning that the growth of the oxide layer toward the bulk is mainly due to implantation of the ions. On the other hand, the limit of the oxide layers above the initial Be surface keeps growing as the implantation time increases. For the three temperatures considered, the growth is almost identical (within a margin of a few Å), indicating that the limiting process does not seem to be temperature-dependent. Thus, diffusion of Be atoms is unlikely to be the limiting oxidation mechanism for 100 eV ions, instead the transfer of kinetic energy from the O ions to the Be atoms seems to be the main growth mechanism. Looking at the profiles of potential energy per atom, one can see that in the oxidised region, the potential energy drops from about -3.5 eV (pure Be) to a minimum of -5.0 eV. The minimum potential energy is obtained for a O/Be ratio of 1, meaning that stoichiometric BeO is preferentially formed in the oxide layer. However, the potential energy does not reach that of perfect wurtzite (horizontal line in the lower plots) even after cooling of the oxidised structure to 0 K. The time scale for formation of a crystalline wurtzite BeO layer is likely beyond reach for MD, and the oxide layer is still in an amorphous BeO phase (see figure 5a).

For 10 eV ions (figure 6b), the oxide layer has a thickness between 11.0 Å (300 K), and 20.2 Å (1500 K), i.e. an increase in thickness of about 80%. Unlike the 100 eV case, most of the oxidised layer is above the initial Be surface, and only a few atoms are implanted below the surface. As previously mentioned, the limit of the oxide layer below the initial Be surface is almost a straight line at the maximum depth of implantation of the ions, with no clear time evolution toward the bulk or the surface. On the other hand, the limit of the layer above the initial surface keeps growing with time. In addition, a clear temperature effect can be seen that explains the thickness increase between 300 K and 1500 K. In this case, the limiting growth mechanism is a thermally activated process, i.e. the diffusion mechanism of oxide layer growth seem to be significant even at MD time scales for this range of small thicknesses. The effect of temperature is also clearly seen in the potential energy profile in figure 6b. The oxide layer formed at a high temperature is more stable by 0.34 eV/atom, bringing the oxide layer closer to a wurtzite configuration (in terms of potential energy). Again, the minimum potential energy is obtained for a O/Be ratio of 1 at all temperatures.

The observations discussed above for the 10 eV ions are also true for the 1 eV case (figure 6c), but now even more pronounced. At 300 K, the resulting layer after 1000 bombardments can hardly be called an oxide layer, and is essentially a complete layer of adsorbed O atoms on top of the Be surface: the lower limit of the oxide layer almost overlaps with the position of the initial surface, indicating no implantation. The thickness of the layer is only 3.7 Å, which is less than the lattice constant of the wurtzite in the [0001] direction (table 1). In addition, the upper limit of the layer stays constant after about 2 ns (400 incident O ions) and the total number of adsorbed atoms stays approximately equal to 130 O atoms. Consequently, the decrease in the potential energy at 300 K is only 0.59 eV in the O-adsorbed layer compared to the pure Be bulk. At 1000 K and 1500 K, the situation is completely different. The thickness of the oxide layer is significantly increased, and is 14.7 Å at 1500 K, which is comparable to the thickness of the oxide layer for the 10 eV irradiation at the same temperature, and represents an increase in thickness of about 300%. Looking at the limits of the oxide layer, one can conclude that this increase occurs mainly above the initial surface with a clear temperature dependence. In addition, as in the case of the 10 eV irradiation, the potential energy of the oxide layer decreases substantially as the temperature is increased. The oxide layer at 1500 K is 0.84 eV/atom more stable than the layer at 300 K. Thus, the strong temperature dependence even at a very low energy of the incident ions suggests that a thermally activated process, such as diffusion of atoms through the oxide layer, is the origin of the growth of the oxide layer at higher temperatures, as proposed by Roth et al. [73,66,6] To summarise the observations, at 100 eV, no clear temperature effect is observed. Thus, the thickness of the oxide layer is mostly controlled by the implantation of O ions (toward the bulk of the material) and the migration of Be due to energy transfer from the incident O ions to the Be atoms. Below 10 eV, the growth of the oxide layer is thermally activated and mainly directed toward the vacuum. It suggests that the diffusion of Be atoms through the oxide layer is the limiting process, which is in line with experimental observations [73,66,6]. Here, it is worth noting that the oxidation toward the substrate might be underestimated due to the overestimation of the energy barrier for O migration in Be. To quantify this possible underestimation, we can use the Wert and Zener theory [74] to estimate the pre-exponential factor of the diffusion coefficient of O in Be (in our cases ranging around 5-6×10 -7 m 2 s -1 ). Thus, assuming the characteristic diffusion length l diff is estimated as l diff = √ D O t (where D O is the diffusion coefficient of O in Be), the diffusion length in these MD simulations can be estimated for both the 1.6 eV (DFT) and the 2.3 eV (ABOP) migration barriers. For the 1500 K case, where the highest migration depth can be expected, one obtains l DFT diff ≈ 1 Å and l ABOP diff ≈ 2 pm. In both cases, the characteristic length of migration is much smaller than the resulting thickness of the formed oxide layer. Thus, oxygen migration into the Be bulk does not affect the growth rate of the oxide layer, and cannot be the rate-limiting growth mechanism, despite the overestimation of the migration energy in the ABOP.

Finally, the irradiation simulations allow us to investigate the sputtering of beryllium by oxygen ions at different temperatures and incident energies. Figure 7 shows the evolution of the sputtering yields (number of sputtered Be atoms per incoming O ions) as a function of incident O energy for the three different temperatures used in the simulations. In The MD data are compared to the sputtering yield given by BCA using SRIM [START_REF] Ziegler | Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms[END_REF]. The error bars are the standard errors obtained after repeating the 1000 O irradiation series four times for each temperature and energy. addition to the sputtering yields obtained by the MD simulations, SRIM [START_REF] Ziegler | Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms[END_REF] simulations are completed to compare the MD results to results obtained in the binary collision approximation (BCA). For these calculations, we used SRIM version 2013.00 in the Monolayer collision steps/Surface sputtering mode, with a density of 3.01 g/cm 3 and default values for the other parameters. These were the displacement energy of 25 eV for Be and 28 eV for O in BeO, lattice binding energy of 3 eV for both elements, and surface binding energies of 3.38 eV for Be and 2 eV for O. For each energy, 100 000 incoming O ions were simulated.

The SRIM simulations only give non-zero sputtering yields for incident energies above 30 eV. In addition, the sputtering yield calculated by SRIM is underestimated by one order of magnitude compared to the MD-calculated values. This large discrepancy is explained by the strong chemical sputtering of Be by O ions. Chemical sputtering is dominant in the range of ion energies and temperatures considered here. At 1 eV and 10 eV, the sputtering yield is even purely chemical. Thus, from these simulations, it is clear that to accurately tackle the issue of beryllium sputtering by oxygen ions in tokamaks, MD is needed. The products of the chemical sputtering are mainly BeO and BeO 2 , even though a few larger molecules like Be 2 O 3 or Be 3 O 4 are eventually observed once a O/Be ratio of 1 is obtained in the oxidised layer. The chemical erosion for 100 eV irradiation is found to be dependent on the state of the oxidised layer. In the first stage of the irradiation, only single Be atoms are eroded from the surface, but once the oxidised layer has a O/Be ratio of about 1, molecules are eventually sputtered. From figure 7, it is also clear that the erosion is much more pronounced at higher energies, and that the temperature, by increasing the motion of the beryllium atoms on the surface, increases sputtering from the surface. Thus, for 1 eV irradiation, no sputtering is observed at 300 K while at 1000 K and 1500 K, significant erosion of the material is observed.

3.3.2.

Thermal stability and elasticity of BeO nanosheets and nanotubes. We applied the Be-O ABOP in an investigation of the thermal stability and elasticity of BeO nanosheets and nanotubes. Carbon graphene has been shown to be thermally very stable with an estimated melting point of around 4500 K [START_REF] Los | [END_REF].

To estimate the melting point of the corresponding BeO sheet with the ABOP, we carried out MD simulations with different heating rates. The system consisted of 800 atoms with periodic boundaries in the planar dimensions. The BeO sheet was heated up at zero pressure with three different heating rates, 10 K/ps, 5 K/ps, and 2.5 K/ps. The observed melting point was about 3300 K and did not significantly change with the different heating rates. More sophisticated methods have been used to obtain more accurate estimates of the melting point of (carbon) graphene [START_REF] Los | [END_REF], but this is beyond the scope of this work. The obtained melting point of 3300 K should therefore be considered an upper limit predicted by the ABOP.

Snapshots from the melting simulation are seen in figure 8.

During high temperatures, breaking of individual Be-O bonds was common, with neighbouring Be and O atoms moving in opposite outof-plane directions. However, individual bond breaking rarely resulted in any permanent damage, as the bond was generally reformed quickly. At temperatures close to the observed melting point, bond breaking often lead to chains or entire hexagons of BeO being temporarily detached from the nanosheet, connected only by chains of Be and O atoms as seen in figure 8a. Figure 8b shows the onset of the melting process, where the first permanently damaged regions are formed. Melting occurred by the entire nanosheet separating into chainand ring-like BeO structures with two-and three-fold coordinations, as seen in figure 8c.

The stability of the BeO nanosheet close to the observed melting point was assessed in a constant temperature simulation at 3000 K. Previously, Wu et al. [15] studied the thermal stability of a BeO nanoribbon with ab initio MD, and observed that the nanoribbon remained intact throughout a 5 ps simulation at 2000 K. Using the ABOP, we can extend the time scale to nanoseconds. The BeO sheet remained stable throughout the 1 ns simulation time, altough two BeO molecules were released, leaving vacancies in the monolayer.

The elasticity of nanotubes predicted by the ABOP was compared with published DFT results of the Young's modulus. The Young's modulus for 3D structures is typically calculated as the second derivative of the energy with respect to the strain at equilibrium, divided by the equilibrium volume. Due to the ambiguity of defining the shell thickness when calculating the volume of a nanotube, the Young's modulus for nanotubes is conventionally defined [77] using the equilibrium area of the cylindrical nanotube shell S 0 , as

Y = 1 S 0 ∂ 2 E ∂ε 2 ε=0 , (16)
where S 0 = 2πLR, L is the length, and R is the radius of the nanotube. Baumeier et al. [11] calculated the Young's moduli for both armchair, expressed by chiral indices (n, n), and zigzag, (n, 0), BeO nanotubes as a function of the diameter. The corresponding Young's moduli for both armchair and zigzag nanotubes were calculated using the ABOP. The length of the nanotubes was about 50 Å, and periodic boundaries were used in the length dimension. The potential energy as a function of strain at 0 K was extracted by statically straining the nanotube around the equilibrium length, and carrying out conjugate gradient energy minimisations at each strain value. The Young's modulus was then obtained from a second order polynomial fit to the strain-energy curve according to equation 16.

Figure 9 shows the Young's modulus as a function of diameter calculated with the ABOP and compared to DFT results by Baumeier et al. [11]. The ABOP slightly overestimates the Young's moduli compared to the DFT data, but the overall trend is similar. The zigzag type nanotubes are elastically softer than the armchair nanotubes, and the Young's modulus rapidly tends towards the value of the infinite graphene-like sheet, in agreement with DFT. Recently, Baima et al. [12] also calculated the Young's modulus of zigzag BeO nanotubes in DFT, and obtained values slightly higher than the results of Baumeier et al. given in figure 9. The calculated Young's modulus of the infinite nanosheet with the ABOP is 0.163 TPa nm, compared to the DFT values 0.13 TPa nm by Baumeier et al., and 0.137 TPa nm by Baima et al.

Finally, we used the ABOP to estimate the temperature-dependence of the Young's modulus of BeO nanotubes. Figure 10 shows the obtained data for an armchair (7, 7) and a zigzag (12, 0) nanotube. The Young's moduli were calculated from a linear fit to the stress-strain response of the stepwise strained [11]. The limiting values for the infinite graphene-like BeO sheet are given at the right. nanotube, which proved to be more accurate than using the potential energy at high temperatures. Each data point in figure 10 is an average of three simulations. The armchair and zigzag nanotubes show identical elastic softening as a function of temperature, and since the Young's modulus for both types quickly becomes independent of diameters (figure 9), a similar behaviour can also be expected for BeO nanotubes of any size. Figure 10 shows that the Young's modulus drops rapidly in the 0-100 K range, after which the it decreases linearly to about 73% of the 0 K value at 2000 K. A rapid drop in the Young's modulus has previously been observed for carbon nanotubes at around 1100 K [78], similar to the low-temperature behaviour observed here, and was hypothesised to be due to the onset of a subtle change in the atomic structure. In our simulations, however, no clear structural changes could be identified. Nevertheless, in reality, zero-point vibrations would likely dominate over thermal vibrations at low temperatures and result in a different behaviour than predicted by MD.

Conclusions

An analytical bond order potential for the Be-O system was fitted and tested against a large database of density functional theory results. The potential shows promising capabilities for simulations of the fusionrelevant Be-O surface interactions. The structural, elastic, thermal, and phonon properties of BeO are well reproduced, as well as the energetics of simple Be-O molecules. Point defect and diffusion properties in Be and BeO obtained by DFT are qualitatively reproduced by the potential. The potential was applied in simulations of the early stages of oxidation of a Be surface by irradiation at different energies and temperatures. We observed the formation of an oxide layer on the Be surface and investigated the possible growth mechanisms. Additionally, we found that chemical sputtering of different Be-O molecules is significant even at very low ion energies. Finally, the potential was shown to be suited for simulations of BeO nanotubes and nanosheets. We provided an estimate of the melting point of BeO nanosheets as well as for the temperature-dependence of the Young's modulus of BeO nanotubes.

Figure 1 :

 1 Figure1: Thermal expansion coefficients, α i∈{x,y,z} , obtained from the variation of the simulation box lengths, and the directionally averaged coefficient α = 1/3(α x +α y +α z ) compared with experimental data[60,61].

Figure 3 :

 3 Figure 3: Common interstitial sites in the hcp lattice. BT: basal tetrahedral, BO: basal octahedral, BC: basal crowdion, NBT: non-basal tetrahedral, NBO: nonbasal octahedral, and NBC: non-basal crowdion.

Figure 4 :

 4 Figure 4: Adsorption sites on a hcp (0001) surface.

Figure 6 :

 6 Figure6: Profiles of O and Be atoms (upper figures), time evolution of the limits of the oxide layer (middle figures), and potential energy per atom (lower figures) as functions of depth, after bombardment of 1000 O ions on a (0001) Be with energies 100 eV (a), 10 eV (b), and 1 eV (c). The initial position of the (0001) surface is indicated by the black vertical lines, and the black horizontal lines correspond to the potential energy of the perfect wurtzite BeO structure. In some cases, particles have been sputtered during the 1000th irradiation, and are not taken into account in the plots. In these cases, the potential energy goes to 0 eV (and marks the position of the surface of the oxide layer).

Figure 5 :

 5 Figure 5: Snapshots of the beryllium surfaces bombarded at 300 K by 1000 O ions with energies 100 eV (a), 10 eV (b), and 1 eV (c). The grey circles correspond to beryllium atoms and white circles to oxygen atoms.

Figure 7 :

 7 Figure 7: Sputtering yields of Be by O ions at different temperatures and energies. The MD data are compared to the sputtering yield given by BCA using SRIM [75]. The error bars are the standard errors obtained after repeating the 1000 O irradiation series four times for each temperature and energy.

  (a) 3200 K, 1.282 ns. (b) 3280 K, 1.316 ns. (c) 3310 K, 1.333 ns.

Figure 8 :

 8 Figure 8: Snapshots (side and top views) of the melting of a BeO nanosheet simulated with the ABOP. Be atoms are coloured grey and O atoms white.

Figure 9 :

 9 Figure 9: Young's modulus of armchair (n, n) and zigzag (n, 0) BeO nanotubes as a function of the diameter, calculated with the ABOP and compared with DFT results [11]. The limiting values for the infinite graphene-like BeO sheet are given at the right.

Figure 10 :

 10 Figure 10: Young's modulus as a function of temperature, normalised by the 0 K value, as predicted by the ABOP for both armchair (7, 7) and zigzag (12, 0) BeO nanotubes.

Table 1 :

 1 Properties of different existing and hypothetical BeO phases obtained with DFT and the ABOP, compared with available experimental results. The listed properties with units are; r 0 : bond length ( Å), E coh : cohesive energy per atom (eV/atom), a, c: lattice constants ( Å), z: internal coordinate of wurtzite, B: bulk modulus (GPa), B : pressure derivative of the bulk modulus, c ij : elastic constants (GPa).

	DFT	Experiment ABOP

c , 147.7 d , 155 e 164 c 66 145.3 152 c , 167.0 d , 185 e 183 BeO zinc blende (F 43m, no. 216)

  -O [29] Be-O

	D 0 1.03571	5.166	6.783
	r 0	2.07880	1.2075	1.41
	β	1.3	2.3090	1.6525
	S	1.88982	1.3864	1.8754
	γ c	0.3940568 0.035608 1.4 8.19587 × 10 -7 0.82595 89.3894
	d	0.27443	0.046496 0.821
	h	0.7606934	0.45056	0.488
	R	2.535	2.1	2.5
	D	0.15	0.2	0.2
	α	0.0	0.0	0.0
	ω	1.0	1.0	1.0
	b f	15.0	12.0	15.0
	r f	0.8	0.5	0.8

desired ground state phase is correctly predicted by the potential, and that no phase not considered in the fitting process is lower in energy. The ABOP was therefore subjected to extensive annealing and heating/quenching MD simulations of wurtzite BeO, and no recrystallisations into other phases were observed.

Table 3 :

 3 Bond lengths ( Å) and cohesive energies (eV/atom) of different Be-O molecules predicted by the ABOP and compared with DFT results and experimental data[START_REF] Chupka | [END_REF].

	DFT Exp. [63] ABOP

experimental or DFT dimer energy is not compatible with the cohesive energies of the bulk BeO phases, and the bond energy of the Be-O dimer is therefore overestimated in the ABOP. (BeO) n molecules were experimentally predicted to form closed ring-like structures

[START_REF] Chupka | [END_REF]

.

Table 4 :

 4 Formation energies (eV) of O interstitial atoms in Be. The different interstitial sites with abbreviations are illustrated in figure3.The values from Ref.[64] are calculated as the reported accommodation energies plus half the dissociation energy of an O 2 molecule (2.56 eV).

	Site		DFT		ABOP
		This work Ref. [64] Ref. [65]
	BT BO BC NBT	-0.85 1.82 --unstable unstable unstable -2.50 -2.37 -4.09 0.21 0.39 -1.03 ---2.85
	NBO NBC r-NBT	0.59 --	-1.13 --	-0.87 -0.41 -	-4.06 -2.90 -4.71

Table 5

 5

	: Adsorption energies and relaxed distances
	from the surface layer for O atoms on a Be (0001)
	surface. See figure 4 for the positions of the different
	sites.
		DFT	ABOP
	Site	E ad (eV) d ( Å) E ad (eV) d ( Å)
	fcc hcp bridge on top -0.74 -3.10 -3.54 unstable 0.91 0.82 1.52	-6.22 -5.91 -4.72 -2.11	0.59 0.85 0.71 1.47

Table 6 :

 6 Formation energies (eV) of point defects in BeO. The Be vacancy and antisite both lead to a magnetisation of the cell in DFT, as indicated below the table.

		DFT ABOP
	V O	4.42	1.62
	V Be	5.47 a	0.55
	O Be 9.57	11.64
	Be O 12.93 b	7.34
	O i	6.00	6.75
	Be i	9.72	7.50
	a

Acknowledgements

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The authors of this work were granted access to the HPC resources of IDRIS and CINES under the allocation A0020806612 made by GENCI (Grand Equipement National de Calcul Intensif) and to the Marconi Supercomputer at CINECA SuperComputing Application and Innovation Department, Bologna, Italy.

Grants of computer capacity from CSC -IT Center for Science, Finland are gratefully acknowledged. JB acknowledges financial support from Svenska Kulturfonden (Arvid och Greta Olins fond).