domain: math info

text:
Benoit Celse

J Becker

D Sinoquet

V Costa Ifpen

PARAMETER FITTING: WHICH ALGORITHM TO CHOOSE ?

Keywords: Optimization algorithms, Derivative Free Optimization fit, nonlinear model, Hydrocracking, kinetic model Mackie 2015

Introduction

Hydrocracking (HCK) is one of the most versatile petroleum refining processes. It usually converts a heavy, low quality feedstock (VGO: Vacuum Gas Oil) into lighter, valuable transportation fuels, contributing significantly to the overall profitability of the refinery [START_REF] Charon-Revellin | Kinetic Modeling of Vacuum Gas Oil Hydrotreatment using a Molecular Reconstruction Approach[END_REF]. A robust kinetic model allows the optimal process design and operating conditions to be chosen to maximize the desired cuts and product characteristics. Hydrocracking of VGO residue performed in a two-step process: 1) a hydrotreatment step in the first reactor, which serves mainly to remove nitrogen and sulfur from the feed; 2) a hydrocracking step in the second reactor, which performs the main hydrocracking reactions on a zeolite-type catalyst. A kinetic model is defined by its structure and its parameters, which are estimated from collected data. The main difficulty is the parameter fit on real data. The aim of this paper is to compare several optimization algorithms on a continuous lumping model of the hydrocracking step. It is structured as follows:

 Case study description,  Description of the chosen optimization algorithms,  Results.

Case study: hydrocracking

The experimental runs presented in this study were performed in a pilot plant at IFP Energies Nouvelles, Solaize, France. The hydrocracking step was performed on a commercial zeolite cracking catalyst. The plant consists of a number of fixed beds, up-flow reactors, designed to mirror the operating conditions in industrial hydrocracking units. A series of mass balances with different operating conditions is thus taken from each experimental run. Each mass balance corresponds to a single experimental point.

Analyses were performed on the feedstocks, the liquid and gaseous effluents. The most relevant measurements in this study were the feed sulfur and nitrogen contents, the partial pressures of NH 3 , H 2 S, and H 2 gasses, as well as the simulated distillation (SIMDIS). A calibration database, consisting of 29 mass balances, was compiled. This database was used for the identification of the empirical parameters in a continuous lumping model originally developed for the first reactor but extended to the second one [START_REF] Becker | Hydrotreatment modeling for a variety of VGO feedstocks: A continuous lumping approach[END_REF]. The range of the main operating conditions, temperature and liquid hourly space velocity (LHSV) is classical: T[370;400°C], LHSV[0.5;3 h -1 ]. The conversion of the 370°C fraction (X370+) is between 50 to 90%w/w.

Optimization method tested

The parameter fit problem is formulated as the following least-square minimization problem This optimization problem belongs to the class of derivative free or black-box optimization problems. Indeed, the simulator associated with the kinetic model computes the simulated data to be compared with experimental ones but does not provide the associated derivatives with respect to the parameters. Classical optimization methods, namely gradient based methods, require those derivatives: therefore, when they are not available as simulator outputs, they are usually estimated by finite difference scheme, which requires NP simulations for each gradient computation (NP being the number of parameters). The associated computational cost may thus become too high when the number of parameters increases. Moreover, the tuning of the perturbation step may be cumbersome in practice when numerical noise is present in the simulation.

An alternative is derivative free optimization methods which have become very popular with the emergence of adapted trust region methods [START_REF] Powell | The NEWUOA software for unconstrained optimization[END_REF][START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF]. SQA, developed by [START_REF] Langouët | Optimisation sans dérivées sous contraintes. Deux applications industrielles en ingénierie de réservoir et en calibration des moteurs[END_REF] at IFPEN, is a trust region method coupled with interpolating quadratic models. It has been applied successfully on several industrial applications [START_REF] Font | Derivative free optimization method and physical simulations coupled with statistical models for transient engine calibration[END_REF][START_REF] Langouët | A derivative free optimization method for reservoir characterization inverse problem[END_REF] and has shown better performances than gradient based methods with finite difference estimate of derivatives. This method is an extension of the method proposed by Powell in [START_REF] Powell | The NEWUOA software for unconstrained optimization[END_REF] to nonlinear constrained problems. The main principle of the method is the following: in order to save simulations, quadratic interpolating models are used as surrogate of the simulator responses. These quadratic models are minimized and updated thanks to additional simulations performed along the optimization process iterations. For leastsquare formulations, the quadratic models approximate the residuals (differences between experimental and simulated data) instead of the single objective function, in order to improve the accuracy of the approximation.

In the following section, SQA method is compared to the Gauss-Newton method implemented in Port library (DN2FB method) [START_REF] Dennis | An Adaptive Nonlinear Least-Squares Algorithm[END_REF]. For this latter method, the Jacobian matrix (derivatives of simulated data with respect to parameters) is estimated by finite differences.

Results

The performances of the SQA and Port algorithms for the Hydrocracking model parameter identification have been compared. The continuous lumping model requires a total of NP=28 experimental parameters to be identified. A target function, based on the yield structure (standard cuts: <150°C, 150-250°C, 250-370°C, >370°C and C1-C4 gas) is defined. The same target function and initial parameters were used for both algorithms.

The SSQ is decreased by an order of magnitude, from 1.7x10 9 to 4.1x10 8 by the Port algorithm and to 4.2x10 8 by the SQA algorithm. The evolution of the sum of square residuals with the number of function evaluations for the SQA and Port (DN2FB) algorithms are shown in figure 2. Both methods attain a minimum. The gradient base method (DN2FB) reached a minimum after around 240 function calls. The response-surface based method (SQA) reached a minimum after 190 function calls. The decrease of the SSQ is very gradual for the DN2FB algorithm. The SQA algorithm first constructs the interpolation of the response surface for a total of 2xNP+1 = 57 function calls. This is characterized by an exploration phase around the initial value of the SSQ. Once the response surface has been constructed, the algorithm rapidly descends towards the minimum (about 60 function calls). The baseline of the SSQ does not decrease significantly after this step. A number of oscillations can be observed, as the algorithm refines the search for the local minimum. Both of the local optimization algorithms, DN2FB and SQA are found to converge to a local minimum with the same square residual. The trust-region method (SQA) shows better performance on the test case than the gradient based-method with finite difference approximation (DN2FB). Once the response surface has been constructed (i.e. exploration around initial point in Figure 1), the SQA algorithm rapidly attains the local minimum. The DN2FB algorithm decreases the residuals much more gradually. Furthermore, the parameter space is explored much more thoroughly by the SQA algorithm, which is illustrated by the large spikes in Figure 1. The finite difference algorithm perturbs only slightly the parameters around each iterate for gradient estimation. This feature makes SQA less likely to be caught in a local minimum.

Conclusion

This paper compares two optimization methods for fitting kinetic model parameters. In this example, the SQA method proves to be more efficient than classical gradient method.

The fact that two different sets of parameters yield to the same sum of square residuals suggests the presence of local minima. Adding chemical/physical a-priori information can help to constrain the system and to remove this under-determination.

Figure 1 :

 1 Figure 1 : Sum of Square Residuals (SSQ) with function evaluations for SQA and Port (DN2FB) algorithms Bibliography


Frédéric Lançon
email: frederic.lancon@cea.fr

Nina Gunkelmann

Damien Caliste

Jean-Luc Rouvière

Incommensurate grain boundary in silicon and the silver-ratio sequence

A scheme is proposed to solve the structure of incommensurate interfaces, starting from highresolution images of electron microscopy, supplemented by adapted simulation techniques and complemented by theoretical tools. Direct silicon bonding is a way to produce artificial interfaces, in particular incommensurate ones. We focus on a technology-driven tilt grain boundary in silicon. While the Fibonacci sequence, linked to the golden ratio, is a prototype of the quasicrystalline structures, a silver-ratio sequence allows us to analyze this incommensurate interface. The four-fold coordination of the Si atoms is kept at the interface.

INTRODUCTION

Direct wafer bonding 1,2 is a pathway to produce grain boundaries that may not occur naturally. In particular, attractive structures for semiconductor technology can be produced by bonding together a Si{011} wafer and a Si{100} one, with the [001] direction of the first wafer aligned with the [110] direction of the second one [3][4][5][START_REF] Signamarcheix | Crystallographic orientation engineering in silicon-on-insulator substrates[END_REF][START_REF] Kato | Structural change of direct silicon bonding substrates by interfacial oxide out-diffusion annealing[END_REF] . Indeed, in this system with mixed orientations, the hole mobility and the electron mobility have opposite enhancement in the wafers, leading to important technological advancements [START_REF] Sung | High performance cmos bulk technology using direct silicon bond (dsb) mixed crystal orientation substrates[END_REF] .

However, the resulting 90° 1 1 0 tilt boundaries are incommensurate, since the two aligned directions [001] and [110] have periodicities with a ratio √ 2 between them. Because of this lack of periodicity of the grain boundaries, solving their atomic structures is not simple, though necessary to fully understand their properties. Several attempts have already been conducted 5,[START_REF] Kariyazaki | Molecular simulation on interfacial structure and gettering efficiency of direct silicon bonded (110)/(100) substrates[END_REF] , but without considering their incommensurate nature.

Although more complex than the periodic cases, the theoretical background to describe incommensurate grain boundaries was quickly provided as an extension of the quasicrystal understanding [START_REF] Rivier | A botanical quasicrystal[END_REF][START_REF] Gratias | Hidden symmetries in general grain boundaries[END_REF][START_REF] Sutton | Irrational tilt grain boundaries as onedimensional quasicrystals[END_REF][START_REF] Deymier | Experimental evidence for a structural unit model of quasiperiodic grain boundaries in aluminium[END_REF][START_REF] Sutton | Irrational interfaces[END_REF][START_REF] Sutton | Interfaces in Crystalline Materials[END_REF] . However, nowadays, high-resolution electron microscopes, as well as mature and precise simulation techniques [START_REF] Lejaeghere | Reproducibility in density functional theory calculations of solids[END_REF][START_REF] Genovese | Daubechies wavelets as a basis set for density functional pseudopotential calculations[END_REF][START_REF] Lançon | Towards simulation at picometer-scale resolution: Revisiting inversion domain boundaries in GaN[END_REF] , can be pursued to give accurate characterizations and quantitative analyses of the incommensurate grain boundaries.

In this work, we focus on the 90° 1 1 0 tilt boundary, but the approach is general to incommensurate interfaces in covalent materials. Our purpose is not to completely explain the atomic structure of the 90° 1 1 0 tilt grain boundary by energy considerations [START_REF] Cherkashin | Application of the O-lattice theory for the reconstruction of the high-angle near 90 degree tilt Si(110)/(001) boundary created by wafer bonding[END_REF] , but rather to extract the local atomic structures from high-resolution experiments. Indeed, through pattern matching, the characteristics of the interface are highlighted and twodimensional coordinates are extracted. In principle, for a non periodic interface, a large area of the material should be analyzed to determine the subtle crystallinesite distortions at the interface. In our case however, two elemental units can be selected and arranged to produce infinite sequences relative to the silver ratio 1 + √ 2, thanks to the quasicrystalline framework (see section VI). From the 2D coordinates extracted from the experimental images of these units, we introduce a two-step energy minimization to retrieve the tridimensional atomic coordinates of the interface elements (see section IV). The atomic structures are then refined by nowadays precise DFT calculations (density functional theory).

Finally a complete atomic description of this incommensurate grain boundary in silicon is given. Note that a similar but metallic incommensurate boundary has been described with the same level of details [START_REF] Gautam | Atomic structure characterization of an incommensurate grain boundary[END_REF] . However, we address here a different class of incommensurate grain boundaries where the structural unit model [START_REF] Sutton | Irrational interfaces[END_REF][START_REF] Sutton | Interfaces in Crystalline Materials[END_REF] is well adapted to describe them.

I. GRAIN-BOUNDARY DESCRIPTION

A 90° 1 1 0 tilt grain boundary in a cubic material is an interface between two crystalline grains having a common crystallographic direction 1 1 0 but where one grain is rotated by 90°with respect to the other one, around this common axis. Because the cubic cell of side a is projected into a rectangle of sides a and √ 2a onto a {1 1 0} plane, the two grains have no coincidence site lattice in this plane. If there was one, the site coordinates of this lattice expressed respectively in both 2D rectangular cells would imply relations na = m √ 2a with n and m being integers. The grain boundary is thus called incommensurate.

The orientation of the boundary itself is fixed by the grain bonding. In this study, the interface is parallel to the common axis and to one of the sides of the 2D rectangular cells (see figure 1). More precisely, [1 1 0] is the common axis of the two grains and taken as our z axis. This axis is therefore a commensurate direction of the interface. The interface contains the [110] direction of one grain, let's say grain I, and thus contains the [001] direction of grain II. This direction is the incommensurate direction of the interface and is chosen to be our x axis. The y axes is perpendicular to the interface.

II. APPROXIMANTS

The periods of grains I and II along the incommensurate direction (x) are denoted s and , respectively (s for short and for long). Grains I and II are, respectively, the lower and upper grains in figure 1. In the ideal incommensurate grain boundary, s = a/ √ 2 and = a, where a is the cubic lattice parameter of silicon. However, the atomic configurations considered for computation usually correspond to periodic approximants of the incommensurate grain boundary. In the same way that a sequence of rational numbers n s /n can tend to an irrational number, a sequence of periodic approximants can approach better and better an incommensurate structure. In the context of this article, a "n s :n approximant" is a periodic boundary where n s periodic cells of grain I are facing n cells of grain II. The best matches between the two grains are obtained when the ideal lengths n and n s s have close values. This corresponds precisely to the strong condition to be verified by n s and n for n s /n to be among the best rational approximations of the ideal √ 2 here:

∆ [n s , n ] ≡ n √ 2 -n s < q √ 2 -p (1)
for any rational number p/q so long as q < n . With this condition, the best approximations are the continuedfraction convergents of √ 2. These convergents are 1/1, 3/2, 7/5, 17/12, and more generally p k /q k = (2p k-1 + p k-2 )/(2q k-1 + q k-2 ) at step k > 2. In our calculations, the period P of the n s :n approximant is the arithmetic mean of the corresponding ideal grain periods:

P = a (n + n s / √ 2)/2 (2)
The residual strains in grains I and II in the approximant are then respectively

ε x,I = (P -n s s)/(n s s)= ( n n s - 1 √ 2 )/ √ 2 (3) ε x,II = (P -n )/(n )= ( n s n - √ 2)/(2 √ 2) = -ε x,I + O ( n s n - √ 2) 2
and thus the minimal strains are also obtained by the best rational approximations n s /n of √ 2. Note that a convergent n s /n of √ 2 corresponds to the convergent n /n s of 1/ √ 2.

III. ELECTRON MICROSCOPY RESULTS

Figure 2 is an electron-microscopy (STEM) image of the grain boundary observed along the common [1 10] direction of the two grains. The experimental details are given in reference [5]. An atomistic model of the periodic 7:5 approximant of this grain boundary had already been proposed 5 . Indeed, a repeating pattern can be recognized with a pseudo periodicity corresponding to n = 5 periods of grain II and n s = 7 periods of grain I. However, even if 7/5 is only 1 % smaller than √ 2, this would correspond to a tensile strain in grain I and a compressive strain in grain II of ±5 × 10 -3 (equation 3), impossible for the material to sustain macroscopically. In fact, a similar but longer pattern can also be identified in the microscopy images. It corresponds to a 10:7 approximant, i.e., n s = 10 and n = 7. While 10/7 does not belong to the best approximations of √ 2, this approximant has about the same strain amount than the 7:5 one, but with an opposite sign for each grain. Thus, combining them together, as in figure 2, can reduce the strains to ±9 × 10 -4 . Indeed, 17/12 = (7+10)/(5+7) is the next best rational approximation after 7/5. In the following, we will consider the 7:5 and the 10:7 approximants as the structural units of the grain boundary. Their 3D structures will be solved by combining image analyses and numerical energy minimizations. The general way to combine them in a sequence of best approximants and thus to describe the infinite grain boundary will be discussed in section VI.

IV. FROM 2D TO 3D Pattern recognition

The resolution of the experimental images is high enough to deduce with a good precision the x and y atomic coordinates of the Si atomic columns perpendicular to the image plane. To automatically extract a large number of coordinates, a pattern recognition technique [START_REF] Paciornik | A pattern recognition technique for the analysis of grain boundary structure by HREM[END_REF] based on cross-correlation functions has been used. From the experimental image, three small portions representative of three material characteristics are selected: the individual atomic column and the two possible orientations of Si dumbbells composed of two columns. Three crosscorrelation functions are first calculated between the image and each of the three templates. As described in [21], three smoother patterns are obtained by averaging image sections located at high correlation maximums (see figure 3). A second set of correlation functions is then calculated. Taking advantage of the three-dimensional nature of the color image coding, we introduce here a way to gather these three functions and get a direct visualization of the information. For all pixels, each color component red, green, or blue, is associated with the value of one of the cross-correlation functions. Figure 3 shows that the different constituents of the grain boundary are now clearly identified. Then, isolated columns, or pairs of them, are associated with Gaussian functions that are fitted on the initial image using the locations of the correlation maxima as initial guesses. To improve the recognition of the atom positions, the low intensity around each image spot was assigned to zero thanks to an apodization function. The 2D Si positions finally extracted from the experimental image are also shown in Figure 3.

Stillinger-Weber atomic interaction

The 2D atomic coordinates, say x and y, extracted from experimental images need to be completed by the z coordinates, corresponding to the hidden component perpendicular to the image plane. We introduce here a technique involving, in two steps, two energy-minimization calculations. These calculations are done with the widely used and successfully tested Stillinger-Weber potential [START_REF] Frank | Computer simulation of local order in condensed phases of silicon[END_REF] and with the parametrization of [23].

Far from the interface, the crystal structures of both grains are known and in particular the Si dumbbells correspond to atoms with two different altitudes z along the [1 1 0] direction, common to both grains. The crystal pe- riod along this direction is s = a/ √ 2. We first assign a constant altitude z = 1 /2 s to all atoms except for the Si dumbbells far enough from the interface to clearly belong to crystal sites. For these dumbbells, we set z = 1 /4 s and z = 3 /4 s to their respective two atoms. The isolated Si columns in Figure 3 correspond to atoms linked by bonds parallel to [1 1 0] in order to be tetrahedrally bonded. This is incompatible with a single-period structure along this direction and instead the fourfold coordination symmetry implies a doubling of the period. To get the initial atomic configuration, the atomic layers of two periods are therefore stacked together along z, resulting in a periodicity of 2s. Periodic boundary conditions are set along the z-direction and free surface boundary conditions are set for x and y.

Doing a mere energy minimization starting from this configuration would be impossible because of the short inter-atomic distances resulting from the unrealistic identical z value settings. The configuration would be unpredictably destroyed during the minimization. In a first step, we instead minimize the energy taking the z coordinates as the only variables. Consequently, only the forces along z are involved and they are initially zero by symmetry except for the few atoms in the crystal grains that have been given different z values. At the beginning of the minimization, the surrounding of these few particularized atoms start to move, and then all the z coordinates are progressively revealed. When the z coordinates were first assigned, it made no difference which atom in a dumbbell is up and which atom is down as long as they are coherently displaced in the crystal. However, we have here two crystals and so we have a choice for the second grain once the first one is set. Considering the two possibilities, we have actually constructed two initial configurations that both have had their energy minimized with respect to the z coordinates. Despite these two possibilities, only one correct configuration is obtained. Depending on the cases, either both configurations converge to the same final configuration, or one can be easily discarded because of a residual antiphase boundary inside one of the grains.

As a second step, a regular energy minimization with respect to all coordinates is performed. This regularizes the crystalline structure of the grains since the experimental x and y data where extracted with some random fluctuations around the equilibrium atomic positions.

With the result obtained from an entire experimental image, two slices are cut from the final atomic configuration. The first one corresponds to a 7:5 approximant, the second one to a 10:7 approximant. For each structure, the slice width is chosen to be its ideal periodicity P introduced above (Eq 2). Subsequent energy minimizations are performed to adjust the atomic positions to the new boundary conditions, but now with periodic boundary conditions for the x and z directions and free surface boundary conditions for y. As shown in Figure 4, the final atomic configurations are in very good agreement with the corresponding parts of the experimental images. The coordinate files are available in Supporting Information 24 . Our numerical method has been able to resolve the pentagon and heptagon structures as well as the connectivity of the atoms corresponding to individual spots in the images. All the atoms are now tetracoordinated.

V. FIRST-PRINCIPLE CALCULATIONS

BigDFT calculations

The approach presented in the previous section for recovering the atomic z coordinates can in principle be done with first-principal calculations in the density functional theory (DFT) framework. This is not possible in this study because of the large number of atoms to consider. However, the interface structure is far from the regular lattice structure of silicon. For instance there are pentagonal and heptagonal Si loops instead of the usual hexagonal loops. Thus, to confirm the stability of the atomic configurations found with the effective interaction potential, we perform DFT calculations with the two approximants. This will also give us reliable interface energies. In this section, the configurations stable with Stillinger-Weber interactions are now the starting points for computing the electronic and atomic structures.

We have used BigDFT code [START_REF] Genovese | Daubechies wavelets as a basis set for density functional pseudopotential calculations[END_REF] , which has the particularity of expressing the electronic wave functions with a set of wavelet basis functions. This characteristic makes it possible to control the computational precision in a systematic way, as with plane-wave basis functions, but furthermore gives a local decomposition of the wave functions. Besides the strong compression of the data and the excellent efficiency for parallel calculations both made possible by this locality, a variety of boundary conditions can be used, from isolated molecules to periodic crystals. This point is important here because we need to consider periodic conditions along x and z, but free boundary conditions along y. Indeed, the configurations are periodic approximant along x, the direction z of the interface is commensurate, and the y free boundary conditions correspond to a material slab with surfaces parallel to the interface. As usual, the perturbation by these external surfaces on the bulk electronic structure is minimized here by terminating the dangling bonds with hydrogen atoms. For the same purpose, no surface reconstruction is considered and the external Si atoms are set on their perfect lattice sites. Keeping these atoms fixed on both surfaces would be possible but would preclude any relative global shift of the grains, which is an important parameter of the grain-boundary structures. Instead, free-bloc conditions are used, i.e., all H atoms and the Si atoms in the two outermost external atomic layers at each surface are fixed all together, but their center of gravity follows the resultant of the forces applied to them. This feature has been implemented in the BigDFT code 25 for this study.

To check the independence of the results on the slab width, two different sizes along the y direction are considered. In the 5:7 approximant there are 48 H atoms and either 312 or 408 Si atoms depending on the slab width. In the 10:7 approximant there are 68 H atoms and either 440 or 576 Si atoms. Pseudopotentials [START_REF] Goedecker | Separable dualspace gaussian pseudopotentials[END_REF][START_REF] Hartwigsen | Relativistic separable dual-space Gaussian pseudopotentials from H to Rn[END_REF] are used to simulate the core electrons. The Perdew-Burke-Ernzerhof (PBE) exchange correlation functional [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] is used in the calculation. The choice of the code, of the pseudopotentials, and of the PBE functional has been assessed in a general study [START_REF] Lejaeghere | Reproducibility in density functional theory calculations of solids[END_REF] and gives high precision results. The k-point grid for calculations with the Si cubic cell is a 4 × 4 × 4 Monkhorst Pack mesh. A 6 × 4 × 6 mesh has been used for the orthogonal cell with paramerters a/ √ 2, a and a/ √ 2 (see Fig. 1). The real-space grid spacing between the wavelet-function centers, h grid , is set to 21 pm corresponding to a high length resolution. With these conditions, the calculated value of the stable Si lattice parameter is 0.5466 nm in agreement with the nowadays DFT results [START_REF] Lejaeghere | Reproducibility in density functional theory calculations of solids[END_REF] . For the grain-boundary calculations, the model periodicity along the z axis corresponds to two [1 1 0] periods and thus three k-points are used in this direction. Because of the large value of the periodicity along x and the free-boundary conditions along y, no grid of k-points is necessary for both of these axes.

The density functional calculations consist in minimizing the energy with respect to both the atomic positions and the electronic structure. The energy minimization is stopped when forces on the atoms are lower than 10 meV/ Å. In the final configurations, all the Si atoms at the interface have kept their covalent bonding. The electronic density of the relaxed 7:5 approximant is shown in figure 5 and the coordinates of both approximants are available in Supporting Information 24 .

The interface energy γ I is estimated by removing the other energy contributions from the total potential energy E pot of the configurations: bulk energy Si , the residual elastic energies due to the residual strains (Eq. 3), and hydrogenated-surface energies γ H{1 0 0} and γ H{1 1 0} . To calculate the surface energies we also consider slab geometries for each of the two grain orientations (see appendix and its figure 11). Along y, the thickness of the slabs are similar to those of the grain-boundary configurations. Free surfaces are hydrogenated too. Complete density functional calculations are done for all slabs. Now by considering the same geometry for each grain, but with additional atomic layers to thicken the slabs, we can extract the bulk energy from their total energy differences. This way, we finally get three estimates of the bulk energy, one for each of the grain orientations plus one calculated merely from the cubic cell. The maximum discrepancy between them is only 2 meV per atom. This procedure is actually most useful to obtain all together bulk energy plus elastic energies of the slightly constrained grains. Indeed, the elastic constrains on the slabs are particular, with fixed strains along x and z (corresponding to imposed periodic conditions) and null stress along y (free-bloc conditions). Taking into account these auxiliary calculations, the computed interface energies γ I are 0.652 J/m² and 0.667 J/m², respectively for the 7:5 and 10:7 approximants. The discrepancy between the thickest and thinnest slabs is less then 0.007 J/m². An estimate of γ I for the infinite incommensurate grain boundary, and not approximants, will be given in sec-tionVI and table I.

Another quantity, the adhesion energy γ A is the gain per surface unit between the energy E pot of the configuration and the energy it would have if the two grains would be separated. So here, the higher energy γ A , the more stable is the interface. In our case, this last energy is deduced from slab calculations with a clean free surface on one side and hydrogenated surface on the other side. No complex surface reconstructions are taken into account, except the inter-plane relaxations, which are direct outputs of the calculation. The adhesion energy γ A is equal to 3.26 J/m² for the 7:5 approximant. Though very close, the gain corresponds to a lower energy than the reported value γ A = 3.2 J/m² for an alternate structure of the 7:5 approximant [START_REF] Kariyazaki | Molecular simulation on interfacial structure and gettering efficiency of direct silicon bonded (110)/(100) substrates[END_REF] . For the 10:7 approximant, γ A is found equal to 3.21 J/m². The energy γ A for separating the grains can be compared to the energies 4.65 J/m² and 3.38 J/m² obtained in this work for cleaving Si crystal at 2 fixing the ratio N S /N L = √ 2. The lower limit of the strip goes through the point (1, 0) and the upper limit through (0,1). The strip is half-bounded with the upper excluded here. The square edges and the sites of the square lattice inside the strip are selected. Once projected onto a straight line, they correspond respectively to the tiles and vertices of the tiling. If projected onto a line of slope √ 2 identical to the slope of the strip, the length of the segments would be in the ratio P L /P S = 1/ √ 2. To get the ratio P L /P S = √ 2, they are projected onto a line of slope 1/ √ 2. The equivalence with the cut method is also sketched: a line segment is periodically set on every site of the lattice (for clarity only those selected by the strip are drawn). The vertices of the tiling are now the intersection points between all segments and the line of slope √ 2 going through the origin. To get the correct length ratio, the segments are perpendicular to the line with slope 1/ √ 2. The common length of the segments is the width of the strip along this direction. The white disks inside the strip correspond to the best rational approximations

N S /N L of √ 2, i.e. |N S /N L - √ 2| < |p/q - √ 2|
, ∀ q≤N L . Among them, those with a central black dot correspond to the continued-fraction convergents, which are solutions of the stronger condition Eq 1.

(1 0 0) and (1 1 0) surfaces respectively. These values are similar to those found in [9].

VI. FROM APPROXIMANTS TO A COMPLETE DESCRIPTION OF THE INTERFACE

Hyper-space description and substitution rules

The Fibonacci tiling, or Fibonacci word, is a prototype of 1D quasicrystals [START_REF] Levine | Quasicrystals. I. definition and structure[END_REF] . It is regularly presented to illustrate the strip and projection method [START_REF] Katz | Quasiperiodic structures obtained by the projection method[END_REF][START_REF] Torikai | Binary self-similar one-dimensional quasilattices: Mutual local-derivability classification and substitution rules[END_REF] , the cut method [START_REF] Bak | Icosahedral crystals from cuts in six-dimensional space[END_REF][START_REF] Janssen | Aperiodic crystals: A contradictio in terminis?[END_REF] or the inflation-rule method [START_REF] Sutton | Irrational tilt grain boundaries as onedimensional quasicrystals[END_REF][START_REF] Levine | Quasicrystals. I. definition and structure[END_REF][START_REF] Lançon | Binary tilings and quasi-quasicrystalline tilings[END_REF] . These methods generate or describe the quasicrystalline structures. The Fibonacci tiling is closely related to the golden ratio, i.e. the irrational number ϕ = 1 + √ 5 /2. Here, to combine together replicas of the two approximants described in the previous sections and get an incommensurate grain boundary, we derive similar properties for a tiling now based on √ 2 and associated with the silver ratio δ S = 1 + √ 2 instead of the golden ratio ϕ. To comply with the usual notations and broaden the scope of the discussion, we denote S = [7 : 5] and L = [10 : 7] the two approximants now regarded as structural units, with respective numbers N S and N L in a given sequence (S for small and L for large).

First of all, why a grain boundary should be associated with such a tiling related to √ 2 ? The answer is given by equations 3. They show that the sequence of structural units should maintain, locally and globally, the ratio n s /n the closest as possible to √ 2 in order to minimize the strain and thus the stress in the grains. It can be easily tested that for a sequence of N S and N L units with N S + N L → ∞, the limit of

n s /n = (7N S + 10N L )/(5N S + 7N L ) is √ 2 when the limit of N S /N L is itself √ 2.
The lengths P L and P S of the two units are identified with their periods along x (see Eqs. 2) and the ratio P L /P S is equal to √ 2. The strip and projection method is a method of choice to generate an incommensurate tiling (or sequence) of S and L units. To get a √ 2-tiling, a strip of slope √ 2 is used to select points of coordinates (N L , N S ) from a 2D simple square lattice. This is illustrated in figure 6 with the strip derived from the elementary lattice cell. This correspond to a very simple algorithm using the deviation ∆ introduced in Eq. (1). Starting from no units, we add them one by one with the following rule: from N S and N L units, add one new S unit if ∆ [N S + 1, N L ] < 1, if not add a L unit. In the last case, the numbers satisfied the relation ∆ [N S , N L + 1] ≤ √ 2. These two inequalities are linked to the strong condition for the best rational approximation seen in section II. Fig. 7 shows the relation between the strip-method and the residual strain due to the lattice misfit at the interface.

In these √ 2 tilings, the units can be grouped into larger tiles L 1 and S 1 , which have a length ratio of √ 2 too. This gives the substitution rules:

L → L 1 ≡ SLS (4a) S → S 1 ≡ LS (4b)
where, in the context of this article, L 1 and S 1 are periodic grain boundaries corresponding to juxtapositions of elementary periods of the grain boundaries L ≡ L 0 and S ≡ S 0 . Grain boundary S 1 corresponds to the convergent structure [17 : 12] and L 1 to the semi-convergent structure [24 : 17]. Starting from an initial unit and iterating the substitution rules (4) generates periodic tilings of increasing period. Equivalent by construction, an alternative method consists of using the recurrence relations

L n =S n-1 L n-1 S n-1 (5a) S n =L n-1 S n-1 (5b)
and generating two sets of sequences in parallel. Since at step zero the ratio P L /P S of the tile lengths is √ 2, these relations show that the length ratio of L n and S n is still √ 2 for any generation n, and that these lengths increase by a factor δ S = 1 + √ 2 after each substitution. The original substitution rules leading to the Fibonacci tiling are L → LS and S → L, with the transformation of the numbers N S and N L from one generation to the other strongly linked to the golden ratio ϕ. For the √ 2-tiling, the matrix M that relates the column vector (N L ; N S ) before and after one substitution is

M = 1 1 2 1 (6)
The two eigenvalues of M are the silver ratio δ S = 1 + √ 2 and its conjugate 1 -√ 2 = -1/δ S . Since they are the solutions of a quadratic equation, i.e. the characteristic equation (1x) 2 -2 = 0, with a module greater than unity for the first one and smaller than unity for the second one, the eigenvalue δ S is a Pisot number. This property is coherent with the fact that we are constructing a quasicrystalline tiling [START_REF] Lançon | Binary tilings and quasi-quasicrystalline tilings[END_REF][START_REF] Godrèche | A simple example of a non-Pisot tiling with five-fold symmetry[END_REF][START_REF] Godrèche | Indexing the diffraction spectrum of a non-pisot self-similar structure[END_REF] . The vector 1, √ 2 is an eigenvector of M corresponding to the largest eigenvalue δ S . Therefore the ratio N S /N L and thus the ratio n s /n tend to √ 2 when the number of successive substitutions tends to ∞. Thus, in both grains and away from the interface, the residual strain due to the lattice parameter mismatch actually tends to zero [see Eq. ( 3)].

Note 1. The √ 2-sequence described in this work is very closely related to the octonacci sequence [START_REF] Sire | Geometric study of a 2D tiling related to the octagonal quasiperiodic tiling[END_REF] [START_REF] Vink | Fitting the stillinger-weber potential to amorphous silicon[END_REF] . The energy of the atoms with this potential is color coded (online version) from blue (bulk energy) to red and yellow (highest energy). Note that the image contrast is correlated with the atom energies: when a high energy has been found the contrast is lower, probably indicating structure perturbations like a higher vacancy concentration or larger atomic moves.

(different orders in the tile concatenations in ref. 37, 39, and 40). Indeed, both sequences have the recurrence relation

S n = S n-2 S n-1 S n-1 (8)
deduced here from the equations (5a) and (5b). One sequence is related to the other by choosing L = LS and S = S. So, getting one or the other is a matter of the choice of the two basic tiles, L in the √ 2-sequence being smaller than L . Since the elementary blocks L and S are directly extracted from the electron micrographs, the √ 2-sequence appears as a natural choice here. Note 2. For simplicity in the Fibonacci tiling, the length ratio of the L and S tiles are usually set to ϕ, i.e. the same value as the tile proportion N L /N S . That way, in the strip-and-projection method on a square 2Dlattice, the strip selects the tiles and a projection orthogonal to it gives directly the tiles. Here, in the √ 2-tiling, the length ratio is fixed to √ 2 by the periodicity of the approximants L = [10 : 7] and S = [7 : 5] observed by microscopy, while their proportion N L /N S is the inverse 1/ √ 2 to minimize the strain in the Si crystals.. Thus, the strip and the direction of the projection are not orthogonal as shown in figure 6. Equivalently in the cut method, the segment motif of the 2D-unit cell is not orthogonal to the strip direction, as shown in figure 6 too. A similar result could be obtained from a rhombic 2D-lattice.

Note 3. While the infinite quasicrystalline √ 2-tiling corresponds to the slope √ 2, finite sequences L n or S n can be obtained from the square 2D-lattice with slopes N S /N L . That way, the strip-and-projection method leads to periodic 1D-lattices with L n or S n as the unit cells. Since in our case, the tiles L and S are periodic approximants at the atomic scale, we can build atomistic structures with periodic boundary conditions along the interface, These structures are thus well suitable for [7:5] , γ I, [10:7] , and γ I,[577:408] , respectively, for the periodic interfaces [7:5], [10:7], and [577:408]; γI for the incommensurate interface. Note: to compute these energies the bulk terms have been removed, in particular the elastic energy of the periodic approximants. While γ I,[7:5] < γI < γ I, [10:7] , the total energy is lower for the incommensurate interface, which has no long range elastic stress.

Atomic interactions

γ I, [7:5] simulations at the atomic level. Note 4. From Fig. 7 we can deduce that a random sequence of L and S units would increase the material strain and thus its elastic energy. This would also be the case, for instance, for sequences generated by the rules L → LSLSL and S → SLSS whose substitution matrix has the same eigenvectors than matrix M (Eq. 6) leading to the same tile ratio N S /N L = √ 2. However its largest eigenvalue, 3 + √ 2, is not a Pisot number; the smallest, 3 -√ 2, is larger than unity. As a consequence [START_REF] Lançon | Binary tilings and quasi-quasicrystalline tilings[END_REF][START_REF] Godrèche | A simple example of a non-Pisot tiling with five-fold symmetry[END_REF][START_REF] Godrèche | Indexing the diffraction spectrum of a non-pisot self-similar structure[END_REF] , the extension perpendicular to the strip diverges in the square-lattice of Fig 7 and would thus correspond to a material with a larger strain.

Application of the substitution rules

With the recurrence rules and the atomistic configurations of the structural units, large models of the grain boundary can now be generated. Figure 8 shows part of a 41:29-approximant overlaying the experimental image of figure 2 in a very satisfactory way (sequence LSSLS equivalent to S 2 by periodicity). These larger models can be easily relaxed with the Stillinger-Weber potential [START_REF] Vink | Fitting the stillinger-weber potential to amorphous silicon[END_REF] , which we started with in section IV. The agreement is satisfactory too and figure 8 shows a correlation between the highest atomic energies and the fuzziness of the experimental spots.

The interface energy γ I of the incommensurate grain boundary can now be estimated too (see appendix). Indeed, neglecting the correlation effects between structural units, the interface energies per 2D unit cell of the 7:5 and 10:7 approximants (section V) can be weighted by their respective concentrations, √ 2/δ S and 1/δ S . The 2D-unit cell areas √ 2 a P 7:5 and √ 2 a P 10:7 are averaged with the same respective weights (the z periodicity, √ 2 a, corresponds to two [1 1 0] periods because of the interface reconstruction; the x periodicity corresponds to P in equation 2). Dividing the weighted energy average by the weighted area average yields to γ I =0.660 J/m². Calculated in the same way, the calculated adhesion energy of the incommensurate grain boundary is γ A = 3.23 J/m².

For reference, we also have computed the interface energy γ I obtained with the Stillinger-Weber potential [START_REF] Vink | Fitting the stillinger-weber potential to amorphous silicon[END_REF] . This value, as well as the value obtained with the original parametrization by Stillinger and Weber 22 are given in table I. A recent potential parametrization [START_REF] Pizzagalli | A new parametrization of the stillinger-weber potential for an improved description of defects and plasticity of silicon[END_REF] has also been tested and the result is included in table I. Since much larger configurations can be considered using Stillinger-Weber interactions than by DFT calculations, the interface energy of the 577:408 approximant has also been calculated for these interatomic potentials. This approximant corresponds to the sequence S 5 made of 41 S and 29 L units. It has a period P equal to 221.6 nm along the incommensurate direction. This gives the opportunity to check the validity of partitioning the whole grain boundary into S and L units to compute the interface energy, while neglecting the correlation effects. Indeed, the relative difference is only 0.2% between γ I and γ I,[577:408] taking into account long range elastic stresses as discussed in [START_REF] Sutton | Interfaces in Crystalline Materials[END_REF] (see table I). With the large 577:408 approximant and the Stillinger-Weber potential [START_REF] Vink | Fitting the stillinger-weber potential to amorphous silicon[END_REF] , we can also test the energy difference resulting from alternative sequences of S and L units. We have tested a random sequence of 41 S and 29 L and the extreme case where the S and L units are segregated. In both cases, the energy minimizations lead to a deformation of the layer due to the strain. The interface is not flat anymore and the buckling is larger for the segregated sequence than for the random one. To avoid this artifact, the minimization has been redone with constrains. The y and z coordinates of the first layer at each surface have been fixed, while the x coordinates were unconstrained to let the strain field adjust to the S and L stacks along x. The incommensurate sequence has the lowest energy, while the random and the segregated ones have a surplus of 1 mJ/m² and 7 mJ/m², respectively. These results are coherent with the strain analysis illustrated in figure 7. The random and the segregated sequences have excursions far from the ideal S and L concentrations, the segregated one being the extreme case.

Hull plots

Sharing the same cubic misorientation than the present Si grain boundary, the 90° 1 1 0 tilt grain boundary in gold is therefore also an incommensurate grain boundary. For gold, it has been shown that this grain boundary has a frictionless -or superglide -property [START_REF] Gautam | Atomic structure characterization of an incommensurate grain boundary[END_REF][START_REF] Frédéric Lançon | Quasicrystalline gold interface with a frictionless property[END_REF][START_REF] Lançon | Superglide at an internal incommensurate boundary[END_REF] . Such a property was first introduced for a unidimensional theoretical model of atoms in a periodic potential field [START_REF] Aubry | The new concept of transitions by breaking of analyticity in a crystallographic model[END_REF][START_REF] Quémerais | Ground-states and growth structures for one-dimensional models with periodic to quasiperiodic transitions[END_REF] . Then, it was theoretically studied and experimentally found in several studies on friction [START_REF] Lançon | Aubry transition in a real material: Prediction for its existence in an incommensurate gold/gold interface[END_REF][START_REF] Frédéric Lançon | Quasicrystalline gold interface with a frictionless property[END_REF][START_REF] Sacco | Free sliding in lattices with two incommensurate periodicities[END_REF][START_REF] Hirano | Atomistic locking and friction[END_REF][START_REF] Shinjo | Dynamics of friction: superlubric state[END_REF][START_REF] Cieplak | Molecular origins of friction: The force on adsorbed layers[END_REF][START_REF] Weiss | Dry friction in the Frenkel-Kontorova-Tomlinson model: Static properties[END_REF][START_REF] Sørensen | Simulations of atomic-scale sliding friction[END_REF][START_REF] Hirano | Observation of superlubricity by scanning tunneling microscopy[END_REF][START_REF] He | Adsorbed layers and the origin of static friction[END_REF][START_REF] Müser | Conditions for static friction between flat crystalline surfaces[END_REF][START_REF] Ko | Friction anisotropy at Ni(100)/Ni(100) interfaces[END_REF][START_REF] Müser | Statistical mechanics of static and low-velocity kinetic friction[END_REF][START_REF] Dienwiebel | Superlubricity of graphite[END_REF][START_REF] Socoliuc | Transition from stick-slip to continuous sliding in atomic friction: Entering a new regime of ultralow friction[END_REF][START_REF]Superlubricity[END_REF][START_REF] Brazda | Experimental observation of the Aubry transition in two-dimensional colloidal monolayers[END_REF] .

The frictionless property is unlikely for this silicon interface, because of its reconstruction extension with faceting clearly visible on the calculated-structure image Only one atomic layer perpendicular to z is not shaded by a transparent mask. Far from the interface, the atoms selected by the two lowest horizontal lines are representative of a crystalline atomic row, i, of grain I. In the same way, the two highest horizontal lines correspond to an atomic row, ii, of grain II. At the interface, the horizontal lines correspond to an atomic row, iii, that has been chosen for calculating its hull plots. Any atomic row could be selected, however this one is the closest to the interface with atoms having the same first-neighbor environment than those of grain I. That way, the influence of grain II on grain I can be analyzed. The opposite would be possible by choosing an atomic row in grain II. To simplify the discussion, row i has been chosen to select crystallographically equivalent lattice sites than those of row iii. The atomic configuration has been relaxed with the Stillinger-Weber potential [START_REF] Vink | Fitting the stillinger-weber potential to amorphous silicon[END_REF] to get atomic energies. These energies are color coded with two successive color gradients. From zero to 15 meV, the first scale shows the energy range in rows iii. From 15 meV to 550 meV, the second scale highlights the high energy spots at the interface.

(figure 9) as well as on the experimental images (figure 2 and in references [5, 6, and 19]). However, Aubry hullfunction analysis [START_REF] Aubry | The new concept of transitions by breaking of analyticity in a crystallographic model[END_REF][START_REF] Aubry | Devil's staircase and order without periodicity in classical condensed matter[END_REF] is the way to get a clear signature of this dynamical property from the static structure of an incommensurate atomic configuration. The hull function characterizes the atomic relaxation from a periodic lattice when the atoms are subject to a periodic potential having a different and incommensurate periodicity. Depending on the strength of the interactions, two regimes exist separated by an Aubry transition corresponding to a "breaking of analyticity" 45 when going from the unpinned states to the pinned one. Since the original unidimensional Frenkel-Kontorova case, the hull function has been generalized to more complex systems [START_REF] Gautam | Atomic structure characterization of an incommensurate grain boundary[END_REF][START_REF] Weiss | Dry friction in the Frenkel-Kontorova-Tomlinson model: Static properties[END_REF][START_REF] Aubry | Effective discommensurations in the incommensurate ground states of the extended Frenkel-Kontorowa models[END_REF] .

At a grain boundary, each crystal row parallel to the incommensurate direction is perturbed by the other crystal and we consider the displacements of the atoms from their ideal lattices positions. Thus, we need these lattices positions, and a reference can be selected far from Hull plots of the first atomic row of grain I at the incommensurate interface (see row iii of the lower grain in figure 9). The lattice parameter a is also the periodicity of grain II. Top: atomic relaxation ∆x along the incommensurate direction. Bottom: relaxation ∆y perpendicular to the interface. The ∆y origin corresponds to the average of the y coordinates in row iii. These hull plots are calculated with 17 Si atoms of the [17:12] approximant (large circles) and 577 atoms of the [577:408] approximant (small dots). Function discontinuities cannot be revealed by inspecting the [17:12] results. However, note that the atomic relaxations of the [17:12] approximant are representative of those of the [577:408] one, the small-model data following the trends of the large-model ones. The discontinuities and the apparent fuzziness of the [577:408] curves are intrinsic to this type of incommensurate interface. The curves could be continuous for another material, but here they are discontinuous everywhere. It is a pinned-state characteristic [START_REF] Lançon | Aubry transition in a real material: Prediction for its existence in an incommensurate gold/gold interface[END_REF] .

the interface. This is illustrated in figure 9 where rows i and ii are representative of the lattice grain I and II respectively.

To represent the hull modulation function ∆x (see figure 10), the atomic displacements are plotted versus the place of their ideal sites relative to the perturbation, i.e. the other lattice periodicity. In grain I with period s along x, the ideal coordinate of an atom with a lattice index j ∈ Z, is x 0 j = j s + α. The shift α is the x coordinate of one arbitrary atom of row i, the origin being set here on one arbitrary atom of row ii. The relaxation displacement of an atom in row iii, i.e. its modulation from the ideal lattice site, is ∆x j = x jx 0 j . This dis-placement ∆x j is plotted versus χ j = (x 0 j modulo ), where is the periodicity of grain II. The system being invariant by a shift of grain I because of the grain II periodicity, the modulo function wraps the complete relaxation behavior of the atoms in a function defined on [0, [. The relaxation ∆y j perpendicular to the interface versus χ j is also plotted in figure 10. Here, ∆y j = y jȳ, where ȳ is the average of the y j in row iii. Using a lattice reference rather than ȳ for calculating ∆y j is also possible and would only result in a shift of the origin for ∆y j .

The discontinuous nature of the hull functions is clearly visible in figure 10. An infinitesimal glide of grain I relative to grain II, would correspond to an infinitesimal change of the coordinates x 0 j and consequently to the quantity χ j (∈ [0, [). Discontinuities of ∆x(χ) would lead to finite jumps of the atomic displacements, and thus the system would have to overcome energy barriers. Therefore, this Si grain boundary has a pinned behavior and not a superglide one. The "broken analyticity" is particularly visible in the ∆y hull plot.

VII. CONCLUSION

A general description of an incommensurate grain boundary in silicon has been given. While this paper is focused on a precise structure, it illustrates a general methodology for semi-conductors different from the analysis of incommensurate boundaries in metals [START_REF] Gautam | Atomic structure characterization of an incommensurate grain boundary[END_REF][START_REF] Lançon | Aubry transition in a real material: Prediction for its existence in an incommensurate gold/gold interface[END_REF][START_REF] Frédéric Lançon | Quasicrystalline gold interface with a frictionless property[END_REF] . This differences are mainly related to the behavior of the atomic bindings. Because of the covalent bonds, the boundary here is a quasicrystalline sequence of two elementary units with opposite internal strains: the 7:5 approximant S and the 10:7 approximant L. First, these units have been selected since they are present in the experimental electron-microscopy images. Second, the unit S corresponds to the ratio 7/5, one of the first best rational approximants of √ 2, while the unit L is a kind of complementary unit to S because of its almost opposite residual strain. Unit L corresponds to the semiconvergent ratio 10/7 of √ 2. So, in the material, these two units can be combined together in a suitable way to reduce the strain and stress resulting from the incommensurate periodicities of the crystals on both sides of the boundary. For instance, concatenating L with S gives LS, which corresponds to the next best rational approximant 17/12.

For both elementary units, we have conducted a detailed analysis of the 2D-images obtained by electron microscopy. A pattern recognition technique has been applied on high resolution images and has given reliable in-plane coordinates of the atomic spots. Taking advantage of the three components of color coding, the three correlation functions of the characteristic patterns -two orthogonal dumbbells and one single spot -can be combined together to give a comprehensive image of the grain-boundary. The three-dimensional atomic co-ordinates have been resolved with atomic-forces calculations using a Stilling-Weber potential and a two-step method. A first energy minimization is conducted, restricted only to the atomic coordinates perpendicular to the plane view. Once the third coordinates are recovered, the second step is an usual minimization that rectifies coherently all the coordinates of the structure. Finally, the structure has been refined with electronic structure calculations. The resulting S and L structures are fourfoldcoordinated silicon networks and are perfectly periodic in the two dimensions of the boundary.

To give a complete atomic description of the incommensurate boundary, the elementary units should form a sequence that reduces the interface stress. With tools adapted for quasicrystals, the atomic structures of S and L units can be combined together. In particular, we have given inflation rules related to the silver ratio and suitable for S and L tilings. Si atoms are still fourfold-coordinated in the resulting interfaces.

Suitable for incommensurate structures, the hull function of the grain boundary has been calculated, thanks to the large systems that can be constructed with the inflation rules. It demonstrates that no hypofriction property can be expected here.

This paper describes a structure corresponding to a perfect incommensurate interface minimizing the residual stress in the grains. However, defects of this ideal structure should be present in real materials, depending on the elaboration process. This is indeed what has been observed by electron microscopy and a rich variety of situations can be imagined. It is however out of the scope of this article to describe them here. bulk energy atom and the free surfaces ones γ surfaces to remove their contributions from E total :

γ I S I = E total -N total atomγ surfaces S I where S I is the area of the grain boundary and N total is the number of atoms in the configuration. In this paper, we determine γ I of an incommensurate grain boundary from the contributions of its periodic grain-boundary units. These units are slightly constrained due to the incommensurate periodicities of the grains, while the incommensurate grain boundary itself has no intrinsic long range strain. Thus, to carefully calculate γ I , we must generalize the previous equation. The configurations taken into account are shown schematically in figure 11.

For instance, from the total energy E total of the 7:5 con-figuration ("a" in figure 11) we calculate :

γ I,[7:5] S I,[7:5] = E total -N I h,i -N II j,k γ h,i S I,[7:5]γ j,k S I, [7:5] where h,i and γ h,i are the strained bulk and surface energies calculated from configurations "h" and "i" in figure 11. Similarly, j,k and γ j,k are derived from configurations "j" and "k". Area S I, [7:5] is the interface area of the 7:5 unit. In the same way, γ I, [10:7] can be calculated with the energies of configurations "b" and "l" to "o".

FIG. 1 .

 1 FIG. 1. Crystalline grain orientations on both side of the grain boundary. The Miller indices are relative to the cubic cell. The plan view is perpendicular to [1 1 0], which is the common direction of both grains (z-axes here). The shades of the spheres indicate their altitude z = 0 or z = a/(2 √ 2). The light-shade rectangles of sides a and a/ √ 2 represent the silicon orthorhombic cells with the lattice vectors corresponding to [001], [110], and [1 1 0]. The cubic cell is projected onto (1 1 0) as a rectangle of sides a (side of the cube) and √ 2a (diagonal of the face). The vector [001] of one of the crystal is parallel to the vector [110] of the other one, and thus the grain boundary is incommensurate. In this study, the interface is parallel to the xz-plane.

FIG. 2 .

 2 FIG. 2. Scanning transmission electron microscopy (STEM)image of the incommensurate grain boundary. Si atoms, individual or in adjacent pairs (dumbbells), are all resolved in this image. More precisely, each spot corresponds to an atomic column perpendicular to the plane view. Vertical lines indicate a possible decomposition into 7:5 and 10:7 approximants of the grain boundary (thinner and thicker strips, respectively). Note the good atomic match between atomic patterns at the different line locations making possible periodic boundary conditions for each approximant, or combinations of several 7:5 and 10:7 approximants. value of /s, i.e.√ 2 here:

  FIG. 3. (Color online) Illustration of the pattern recognition technique. (Top-left) Small part of an image, with grain II on top of grain I. (Top-right) Three patterns corresponding to grains I, II, and to isolated atomic columns. (Bottomleft) Color image representing the three correlation functions corresponding to the three patterns. Each correlation function represents one component of the RGB image coding (redgreen-blue). A gamma correction of 4 has been applied to enhance the highest values of the correlations, i.e., a power-law function with exponent 4 is applied to the components (normalized between 0 and 1). (Bottom-right) x and y atomic positions obtained by fitting, in the first image, individual or paired Gaussian functions, initially located at the maximum of the correlation functions. The individual atomic columns are distinguished by the color cyan.

FIG. 4

 4 FIG. 4. (Color online) Electron microscopy images overlaid by the atomic configurations of the 5:7 and 10:7 approximants stable with the Stillinger-Weber interactions (respectively top and bottom images). For the bottom image, two images has been combined to get a sharper contrast. The light vertical lines indicate the periodicity along the interface. The coordinates z perpendicular to the images are color coded from red (lower z) to yellow (higher z). All Si atoms are tetracoordinated. Atoms connected by bonds orthogonal to the images are surrounded by cyan circles.

FIG. 5 .

 5 FIG. 5. Cross-section of the valence-electron density and ball-and-stick model of the 5:7 approximant. White and black in the grayscale represent respectively the highest and the lowest electron density. Left: Oblique view showing the location of the cross-section, which is perpendicular to the commensurate direction of the interface (common [1 10] direction). Along this direction, the periodicity of the structure is √ 2 a, where a is the cubic lattice parameter of silicon. Right: Plane view. The horizontal axis is the [001] direction of the upper grain and the [110] direction of the lower one. It is the opposite for the vertical axis. Along the incommensurate direction, P is the periodicity of the approximant. The logarithmic scale of the density plot is indicated by the gray-scale bar. The average density is Ne/V = 196 nm -3 . The lattice parameter a = 547 pm gives the length scale. The free-surfaces are passivated by hydrogen atoms.

FIG. 6 .

 6 FIG.6. Strip and projection method to generate a √ 2-tiling from the square lattice. The gray strip has a slope equal to √ 2 fixing the ratio N S /N L = √ 2. The lower limit of the strip goes through the point (1, 0) and the upper limit through (0,1). The strip is half-bounded with the upper excluded here. The square edges and the sites of the square lattice inside the strip are selected. Once projected onto a straight line, they correspond respectively to the tiles and vertices of the tiling. If projected onto a line of slope √ 2 identical to the slope of the strip, the length of the segments would be in the ratio P L /P S = 1/ √ 2. To get the ratio P L /P S = √ 2, they are projected onto a line of slope 1/ √ 2. The equivalence with the cut method is also sketched: a line segment is periodically set on every site of the lattice (for clarity only those selected by the strip are drawn). The vertices of the tiling are now the intersection points between all segments and the line of slope √ 2 going through the origin. To get the correct length ratio, the segments are perpendicular to the line with slope 1/ √ 2. The common length of the segments is the width of the strip along this direction. The white disks inside the strip correspond to the best rational approximations N S /N L of √ 2, i.e. |N S /N L -√ 2| < |p/q -√ 2| , ∀ q≤N L . Among them, those with a central black dot correspond to the continued-fraction convergents, which are solutions of the stronger condition Eq 1.

FIG. 7 .

 7 FIG.7. Density plot of the strain εx,II in grain II for sequences made of N L tiles L and N S tiles S due to the misfit between the periodicities along the incommensurate direction of the two crystalline grains (see equations Eq. 3). Strain εx,I in grain I is the opposite. The red line (bold polygonal curve) represents S3 = LSSLSLSSLSLS (Eqs. 4), a sequence that minimizes the strain field locally and on average.

FIG. 8 .

 8 FIG. 8. Top: Overlay of the experimental image by a sequence of S and L tiles, i.e. respectively the [7 : 5] and [10 : 7] approximants, each one calculated by DFT. The vertical lines indicate the tile boundaries. Bottom: Overlay of the experimental image by a [41 : 29] approximant corresponding to the sequence S2, now relaxed with the Stillinger-Weber potential[START_REF] Vink | Fitting the stillinger-weber potential to amorphous silicon[END_REF] . The energy of the atoms with this potential is color coded (online version) from blue (bulk energy) to red and yellow (highest energy). Note that the image contrast is correlated with the atom energies: when a high energy has been found the contrast is lower, probably indicating structure perturbations like a higher vacancy concentration or larger atomic moves.

FIG. 9 .

 9 FIG.9. Atomic-configuration visualization highlighting the calculation of the hull functions ∆x and ∆y plotted in figure 10. Only one atomic layer perpendicular to z is not shaded by a transparent mask. Far from the interface, the atoms selected by the two lowest horizontal lines are representative of a crystalline atomic row, i, of grain I. In the same way, the two highest horizontal lines correspond to an atomic row, ii, of grain II. At the interface, the horizontal lines correspond to an atomic row, iii, that has been chosen for calculating its hull plots. Any atomic row could be selected, however this one is the closest to the interface with atoms having the same first-neighbor environment than those of grain I. That way, the influence of grain II on grain I can be analyzed. The opposite would be possible by choosing an atomic row in grain II. To simplify the discussion, row i has been chosen to select crystallographically equivalent lattice sites than those of row iii. The atomic configuration has been relaxed with the Stillinger-Weber potential[START_REF] Vink | Fitting the stillinger-weber potential to amorphous silicon[END_REF] to get atomic energies. These energies are color coded with two successive color gradients. From zero to 15 meV, the first scale shows the energy range in rows iii. From 15 meV to 550 meV, the second scale highlights the high energy spots at the interface.

  FIG. 10.Hull plots of the first atomic row of grain I at the incommensurate interface (see row iii of the lower grain in figure9). The lattice parameter a is also the periodicity of grain II. Top: atomic relaxation ∆x along the incommensurate direction. Bottom: relaxation ∆y perpendicular to the interface. The ∆y origin corresponds to the average of the y coordinates in row iii. These hull plots are calculated with 17 Si atoms of the[17:12] approximant (large circles) and 577 atoms of the [577:408] approximant (small dots). Function discontinuities cannot be revealed by inspecting the[17:12] results. However, note that the atomic relaxations of the[17:12] approximant are representative of those of the [577:408] one, the small-model data following the trends of the large-model ones. The discontinuities and the apparent fuzziness of the [577:408] curves are intrinsic to this type of incommensurate interface. The curves could be continuous for another material, but here they are discontinuous everywhere. It is a pinned-state characteristic[START_REF] Lançon | Aubry transition in a real material: Prediction for its existence in an incommensurate gold/gold interface[END_REF] .

FIG. 11 .

 11 FIG. 11. Schematic view of the configurations constructed with grains I and II to compute the interface energy γI. The solid lines represent the periodic boundary conditions. The dashed lines show the interfaces and the free surfaces of types (001) or (110). Configurations a and b are the 7:5 and 10:7 approximants, respectively. The tensile and compressive residual strains are shown in blue and red colors, respectively. These strains are due to the incommensurability between grain I and grain II periodicities. For taking into account their effects on the elastic energy and on the surface energies, these strains have been applied to configurations h to o. More precisely, the strains present in the 7:5 approximant (configuration a) have been used in configurations h to k, while those of the 10:7 approximant (configuration b) have been used in l to o. The different layer widths, for instance between configurations n and o, allow us to separately determine the strained-bulk energy and the strained-surface energy. Indeed, only the numbers of bulk atoms are different, while the free-surface areas are identical. Configuration c is a simple crystalline cell used to compute the reference bulk energy. Without applied strain too, configurations d to g provide us with two other ways to compute the bulk energy, as well as the values of the surface energies γ {0 0 1} and γ {1 1 0} .

  Finally, taking into account the respective frequencies c S = √ 2/δ S and c L = 1/δ S of units S = 7:5 and L = 10:7, the interface energy of the incommensurate grain boundary is γ I = c S γ I,[7:5] S I,[7:5] + c L γ I,[10:7] S I,[10:7] c S S I,[7:5] + c L S I,[10:7]

TABLE I .

 I Interfaces energies (J/m²) obtained in this work with DFT calculations and different parametrizations of the Stilling-Weber potential: γ I,

ACKNOWLEDGMENTS

The authors thank Alain Marty for useful discussions and suggestions.

SUPPLEMENTAL MATERIAL

The atomic coordinates of the S and L units are available as Supplemental Material 24 . The files are in plain ascii format and their content are self-explanatory. All coordinates are in Å unit and given in a Cartesian system. Besides the comment lines, the first two lines give the computer-box lattice parameters a, b, and c on the form a x , b x , b y , c x , c y , and c z , while a y = a z = b z = 0. Then, the following lines give for each atom the coordinates x, y, z, and the element name.

Appendix: Interface-energy calculation

To extract the interface energy γ I of a grain boundary from the total energy E total of a configuration made of two crystalline grains, one must separately compute the † Jean-Luc.Rouviere@cea.fr 1