Terkia Braik
email: terkia.braik@im2np.fr

Nathalie Mangelinck-Noël
email: nathalie.mangelinck@im2np.fr

Guillaume Reinhart
email: guillaume.reinhart@im2np.fr

Jean-Marc Debierre
email: jean-marc.debierre@im2np.fr

Interplay between facetting and confinement in directional solidification of silicon: a direct comparison between experiments and phase-field simulations

Keywords: A1. Directional solidification, A1. Computer simulation, A1. in situ X-ray imaging, A1. Solid-liquid interface, A1. Crystal facets, B1. Photovoltaic silicon

We combine the analysis of experiments based on in situ X-ray imaging and phase-field simulations to get insight into the interplay between faceting and confinement effects in thin silicon samples during directional solidification from a monocrystalline seed. The phase-field model is parametrized by available experimental results, especially regarding the anisotropy functions for the surface energy and for the kinetic attachment coefficient. In the present work, a method is proposed to extrapolate 2D fully quantitative simulations results to 3D thin sample geometry. Two simple seed orientations along the 110 and 100 crystal directions are considered. For 110 direction, it is shown that confinement effects are weak so 2D simulations can be readily compared to the experiments. For 100 direction, much stronger confinement effects are identified and characterized by a confinement factor that is estimated from phase-field simulations of thin 3D systems. For both directions, the length of the side facets is shown to increase with the solidification velocity and an analytical expression including the confinment factor is obtained for the law relating both quantities. A very quantitative experiment-simulation agreement is obtained in each case.

Introduction

Directional solidification is one leading process to fabricate silicon for the ever increasing needs of solar energy conversion [1]. To obtain materials with properties meeting the requirements for solar cell fabrication, several types of features and defects appearing during the solidification step must be controlled. During the solidification of silicon, crystal faceting and the associated twin formation are frequent and are reported in industrial ingots where the latter is the main mechanism of the grain structure formation. The corresponding facet dynamics, twin generation mechanisms and the interaction of both with other structural defects are of current interest [START_REF] Huang | Growth and Defects in Cast-Mono Silicon for Solar Cells: A Review[END_REF][START_REF] Lan | The emergence of high-performance multi-crystalline silicon in Photovoltaics[END_REF][START_REF] Wang | Seed-Assisted Growth of Cast-Mono Silicon for Photovoltaic Application: Challenges and Strategies[END_REF][START_REF] Tsoutsouva | In situ investigation of the structural defect generation and evolution during the directional solidification of 110 seeded growth Si[END_REF][START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF]. Other main structural defects are dislocations that can be at the origin of subgrains known to be deleterious for photovoltaic properties and which interact with facets, twins and other defects in general [START_REF] Wang | Seed-Assisted Growth of Cast-Mono Silicon for Photovoltaic Application: Challenges and Strategies[END_REF][START_REF] Oriwol | Small-angle subgrain boundaries emanating from dislocation pile-ups in multicrystalline silicon studied with synchrotron whitebeam X-ray topography[END_REF]. During the directional solidification of pure silicon, for positive temperature gradients generally used in the processes, the solid-liquid interface remains stable without microstructures [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF][START_REF] Riberi-Béridot | Strain building and correlation with grain nucleation during silicon growth[END_REF]. Besides, due to triple junction conditions, the interface usually becomes faceted near the crucible walls [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF]. In stationary growth conditions, the facet lengths are observed to sensitively depend on misalignment effects that are almost unavoidable in the experimental situation. During the solidification of thin model samples, it can be anticipated that the combination of the temperature gradient misalignment angle θ m with the confinement due to a very limited sample thickness (L y = 300 µm here), results in a complex law relating the facet length to its normal growth velocity.

The ultimate goal of this paper is to shed some light on this matter by drawing a comparison between experimental investigations and phase field simulations. Experimental investigations are based on directional solidification experiments of thin model samples in which the interface dynamics and generation of structural defects is characterized by in situ X-ray imaging during solidification [START_REF] Ouaddah | X-ray Based in Situ Investigation of Silicon Growth Mechanism Dynamics-Application to Grain and Defect Formation[END_REF].

The thin interface phase-field model (TIPM) was introduced more than two decades ago by Alain Karma and his collaborators, first for the free solidification of a pure substance [START_REF] Karma | Phase-field method for computationally efficient modeling of solidification with arbitrary interface kinetics[END_REF][START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF] and later for the directional solidification of a dilute binary alloy [START_REF] Karma | Phase-field formulation for quantitative modeling of alloy solidification[END_REF][START_REF] Echebarria | Quantitative phase-field model of alloy solidification[END_REF]. Since then, this model has proved successful, notably in providing quantitative support to theoretical results [START_REF] Singer-Loginova | The phase field technique for modeling multiphase materials[END_REF]. In our previous work, a comparison of TIPM simulations with experimental results for silicon free growth [START_REF] Yang | Crystal growth and equilibrium crystal shapes of silicon in the melt[END_REF] and for directional solidification [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF] recently allowed to quantitatively parametrize anisotropy functions for both the surface energy and the kinetic attachment coefficient [START_REF] Boukellal | Equilibrium and growth facetted shapes in isothermal solidification of silicon: 3D phasefield simulations[END_REF][START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF]. In the case of silicon directional solidification, the quantitativeness of our numerical code was recently assessed by comparing two-dimensional (2D) simulations with an analytical model [START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF]. In the present study, the 3D version of the phase-field code is implemented to conform to the thin layer geometry of the experimental setup.

After a short description of the experimental context, our numerical model is briefly described. The main part of this paper concerns the comparison of the results obtained by both approaches and presents the informations provided by this combined investigation. As a first investigation step, we consider a simple situation without misalignment (θ m = 0) to determine the confinement effects due to the limited sample thickness. Based on the obtained results, the second step is to determine an accurate estimate for the kinetic coefficient β 111 , assuming a linear kinetic growth law for {111} facets. Experimental and phase-field curves relating the side facet length to the facet velocity are then compared. Finally, we establish a simple scaling relation that predicts the influence of θ m on the lengths of both side facets. A quantitative comparison is drawn between the results of experiment and simulation. Finally, conclusions are presented with a brief discussion of possible perspectives.

Experimental

The experiments are performed in the GaTSBI (Growth at high Temperature observed by Synchrotron Beam Imaging) device. It allows to study the directional solidification of silicon (Si) characterized in situ and in real-time during the solidification process. It is composed of a high temperature (up to 1800 • C) directional solidification (DS) furnace coupled with synchrotron radiation X-ray imaging techniques (Bragg diffraction imaging and radiography). A more detailed description of the device and of the X-ray imaging methods can be found in Ouaddah et al. [START_REF] Ouaddah | X-ray Based in Situ Investigation of Silicon Growth Mechanism Dynamics-Application to Grain and Defect Formation[END_REF]. The Si samples (38 mm × 5.8 mm × 0.3 mm), housed in a boron nitride crucible, are introduced inside the DS furnace. The crucible assembly can be seen for example in Figure 1b in [START_REF] Tsoutsouva | In situ investigation of the structural defect generation and evolution during the directional solidification of 110 seeded growth Si[END_REF]. The results shown in the following were obtained during the solidification of two silicon samples. They were cut from the same Si (110) wafer (floating zone process by SILITRONIX) in order to study two orientations along the growth direction: 110 and 100 . Initially, the samples were partially melted to keep a seed. Then, they were directionally solidified from the seed at a constant applied temperature gradient, G appl. = 37.5 K/cm, imposed between the two heaters of the solidification furnace and with a cooling rate of -1 K/min applied on both heaters. In the following, the samples are labelled DS 011 (green color code in the figures) and DS 001 (red color code in the figures) for seed orientation along the solidification direction 110 and 100 , respectively. The crystal axes corresponding to the two seed orientations are represented in figure 1 and superimposed on the octahedral crystal unit (see column (a)).

The facet dynamics of these experiments was studied in particular in Stamelou et al [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF]. In our experiments, both main faces of the silicon samples are in contact with the boron nitride crucible plates [START_REF] Tsoutsouva | In situ investigation of the structural defect generation and evolution during the directional solidification of 110 seeded growth Si[END_REF][START_REF] Jhang | Heterogeneous twinning during directional solidification of multi-crystalline silicon[END_REF]. As a consequence, the solid-liquid interface is as well in contact with these plates, as shown schematically in figure 2a right. A sketch of the front view of the solid-liquid interface can be seen in figure 2a left. It can be directly compared to the X-ray radiography image of the solid-liquid interface (figure 2b). Facets at the side edges of the samples are always experimentally observed at the solidliquid interface, as shown in figure 2b. We showed that they correspond to {111} facets [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF][START_REF] Ouaddah | X-ray Based in Situ Investigation of Silicon Growth Mechanism Dynamics-Application to Grain and Defect Formation[END_REF]. These {111} facets are particularly important for the establishment of the grain structure. Indeed, we showed that the grain structure in pure samples is fully driven by twinning on the {111} side facets at the solid-liquid interface [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF][START_REF] Riberi-Béridot | Strain building and correlation with grain nucleation during silicon growth[END_REF][START_REF] Ouaddah | X-ray Based in Situ Investigation of Silicon Growth Mechanism Dynamics-Application to Grain and Defect Formation[END_REF]. As such, the {111} facet dynamics was studied and their velocity as a function of the undercooling measured; details can be found in [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF]. The right side of figure 2a displays a sketch of the solid-liquid interface viewed from the side. It corresponds to the situation of sample DS 001 . Due to the projection mode of X-ray radiography, we have no information on the precise shape of the solid-liquid within the sample thickness by X-ray radiographs. However, it was shown repeatedly using X-ray topographs and EBSD measurements that the {111} facets cross the whole thickness of the samples from the front to the back boron nitride crucible plates [START_REF] Tsoutsouva | In situ investigation of the structural defect generation and evolution during the directional solidification of 110 seeded growth Si[END_REF][START_REF] Riberi-Béridot | Strain building and correlation with grain nucleation during silicon growth[END_REF][START_REF] Ouaddah | X-ray Based in Situ Investigation of Silicon Growth Mechanism Dynamics-Application to Grain and Defect Formation[END_REF]. The orientation of the {111} traversing facets is determined by the possible {111} plane family orientations corresponding to the seed crystallographic orientation [START_REF] Tsoutsouva | In situ investigation of the structural defect generation and evolution during the directional solidification of 110 seeded growth Si[END_REF][START_REF] Riberi-Béridot | Strain building and correlation with grain nucleation during silicon growth[END_REF][START_REF] Ouaddah | X-ray Based in Situ Investigation of Silicon Growth Mechanism Dynamics-Application to Grain and Defect Formation[END_REF]19]. However, the precise shape of the solid-liquid interface within the thickness can only be obtained by the phase field simulation and will be discussed in the following.

Numerical model

Phase-field model

We adopt the TIPF model defined in [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF], in which the phase field ϕ varies continuously across the interface, between ϕ = 1 far on the solid side (r 0) and ϕ = -1 far on the liquid side (r 0), r being the coordinate along the normal to the interface which is located at r = 0. This variation is only appreciable over a distance of the order of the interface width W 0 which is the governing model parameter. For a flat interface, one shows that in the stationary state,

ϕ = -tanh r W 0 √ 2 . ( 1
)
Our implementation of the TIPM makes use of the reciprocal function

ψ = √ 2 tanh -1 (ϕ) (2)
that is equal to -r/W 0 for a planar interface. This preconditioned phase field ψ was introduced by Glasner [START_REF] Glasner | Nonlinear preconditioning for diffuse interfaces[END_REF] to reduce the simulation times. For the present simulations, we further increase the code efficiency by using Graphics Processing Unit (GPU) parallel programming [START_REF] Ruetsch | CUDA Fortran for scientists and engineers: best practices for efficient CUDA Fortran programming[END_REF]. Similar implementations of the TIPM were recently used in a number of situations where the simulations could be directly compared to experimental data [START_REF] Boukellal | Equilibrium and growth facetted shapes in isothermal solidification of silicon: 3D phasefield simulations[END_REF][START_REF] Ghmadh | Directional solidification of inclined structures in thin samples[END_REF][START_REF] Bergeon | Spatiotemporal dynamics of oscillatory cellular patterns in three-dimensional directional solidification[END_REF][START_REF] Ghmadh | Physical mechanisms of the oscillatory 2 λ-O mode in directional solidification[END_REF][START_REF] Boukellal | Scaling laws governing the growth and interaction of equiaxed Al-Cu dendrites: A study combining experiments with phase-field simulations[END_REF].

We only give the main lines of our TIPM model here; more details can be found elsewhere [START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF]. In order to take into account the dependence of the surface energy on the crystal orientation, the anisotropy function a s ( n) must be defined, where n is the unit vector normal to the solid-liquid interface. Similarly, the kinetic attachment coefficient is characterize by its own anisotropy function a k ( n). Both functions are characteristic of the material under study. In the case of pure silicon, based on symmetry considerations and on experimental results, we recently proposed an analytical form for both functions [START_REF] Boukellal | Equilibrium and growth facetted shapes in isothermal solidification of silicon: 3D phasefield simulations[END_REF][START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF]. In the phase-field model, the singularities of these anisotropy functions have to be smoothed out to avoid numerical divergences. For the surface energy,

a s ( n) = 1 + δ √ cos 2 α + 2 + δ 0 sin 2 α + 2 , ( 3
)
where α is the conical angle between the normal to the interface n and the unit vector n 111 normal to a {111} facet, and where is a parameter much smaller than one that is used to avoid a sharp cusp at α = 0 [START_REF] Wang | Phase-field model of vapor-liquidsolid nanowire growth[END_REF]. For the reasonable choice = 10 -2 , the amplitudes δ 2.5 and δ 0 1.75 were obtained by fitting phase-field equilibrium shapes [START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF] to their experimental counterparts [START_REF] Yang | Crystal growth and equilibrium crystal shapes of silicon in the melt[END_REF]. For the kinetic attachment,

a k ( n) = b 0 exp - θ ∆θ 2 - φ ∆φ 2 , (4)
where θ and φ are two spherical angles separating n from n 111 and where ∆θ and ∆φ are their dispersions [START_REF] Boukellal | Equilibrium and growth facetted shapes in isothermal solidification of silicon: 3D phasefield simulations[END_REF]. The Gaussian amplitude b 0 is proportional to the physical kinetic coefficient β 111 and inversely proportional to the phasefield coefficient β 0 defined below. Moreover, it must also be adjusted to the criterion that defines the facet extension to ensure quantitative results [START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF].

Taking the interface width W 0 and the relaxation time τ 0 as the length and time units, the ψ evolution equation reads

a τ ( n) ∂ψ ∂t = √ 2 ϕ -a 1 ξ(1 -ϕ 2 )u +a 2 s ∇ 2 ψ - √ 2ϕ( ∇ψ) 2 +2a s ∇a s • ∇ψ + √ 2 (1 -ϕ 2 ) ∇ • A, (5)
where

u = c p (T -T m )/l H (6)
is the dimensionless temperature field, with c p the specific heat at constant pressure, l H the latent heat, T m the melting temperature, and a 1 0.8839 a model constant [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF]. In situ observation of the solid-liquid interface dynamics show that, in the low growth velocity regime considered here, a constant temperature gradient settles rapidly in the interface region [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF]. In addition, the latent heat released at the solid-liquid interface is rapidly dissipated by the experimental setup [START_REF] Riberi-Béridot | situ characterization by X-ray synchrotron imaging of the solidification of silicon for the photovoltaic applications: control of the grain structure and interaction with the defects and the impurities[END_REF]. Relying on these observations, we thus chose to impose a frozen linear temperature profile in the simulations,

u(z) = c p G(z -z m )/l H , (7)
using the measured value of the local temperature gradient G. In equation ( 5), the anisotropy vector A concentrates all the terms depending on the surface energy anisotropy a s ( n). Details about this formulation can be found elsewhere [START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF][START_REF] Boukellal | Scaling laws governing the growth and interaction of equiaxed Al-Cu dendrites: A study combining experiments with phase-field simulations[END_REF][START_REF] Boukellal | Tip dynamics for equiaxed Al-Cu dendrites in thin samples: Phase-field study of thermodynamic effects[END_REF]. As a function of the scale parameter ξ, the characteristic length and time read W 0 = ξd 0 , and τ 0 = a 0 (d 2 0 /D)ξ 3 , D being the Si thermal diffusivity, d 0 the capillary length, and a 0 0.5539 [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF]. The kinetic attachment then reads β( n) = β 0 a k ( n) with β 0 = a 0 W 0 /D. Following [START_REF] Boukellal | Equilibrium and growth facetted shapes in isothermal solidification of silicon: 3D phasefield simulations[END_REF] and [START_REF] Lin | Phase field modeling of facet formation during directional solidification of silicon film[END_REF], the anisotropy function a τ is taken as

a τ ( n) = a s ( n)[a s ( n) + a k ( n)]. (8)
The expression of anisotropy function a s is given in the frame of the three main cubic crystal axes, [001], [010], and [001]. In the simulations, the three components of the unit vector n normal to the solid-liquid interface are calculated in the (x, y, z) frame of the numerical domain that is also the frame of the crystal cubic directions that define the seed orientation, as shown in figure 1. Three Euler angles are then used to obtain n 100 , n 010 , n 001 from n x , n y , n z by inverse rotations.

Numerical parameters and conditions

Obtaining physically meaningful simulation results necessitates to carefully specify several numerical parameters and conditions. The first parameter is the interface width W 0 . A detailed convergence test of our phase-field code with this parameter showed that quantitative results can be obtained by decreasing it down to W 0 1 µm [START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF]. For such low W 0 values, the domain size is already very large for 2D simulations (up to 2400 × 2 × 2400 mesh points), so 3D fully quantitative simulations are out of reach. In order to compare our 2D results with experimental ones obtained in thin samples, we adopt the following strategy. First, semi-quantitative 3D simulations for higher values of W 0 are performed to determine the confinement effects in the experimental thin samples. Second, these confinement effects are taken into account to adjust the 2D simulation results to the case of thin 3D systems.

For the quantitative simulations in 2D systems, the domain size is reduced to two mesh points and mirror boundary conditions are imposed in the y (thickness) direction. For thin 3D systems, the number of mesh points along y is increased to reproduce the experimental sample thickness (L y = 300 µm) and different wetting and/or mirror boundary conditions are imposed. Concerning the contacts of the solid-liquid interface with the domain boundaries, it is modeled on the left domain side by imposing

∂ψ ∂x = σ 0 (1 -ϕ 2 ), (9)
where σ 0 is a positive constant. This wetting boundary condition was already used in [START_REF] Ghmadh | Directional solidification of inclined structures in thin samples[END_REF]. Here, we verify that, as observed experimentally [START_REF] Riberi-Béridot | Strain building and correlation with grain nucleation during silicon growth[END_REF], it does produce an almost straight facet on the side of the interface with the crystallographic angle θ f that is expected for the seed crystal orientation (see Table 1). The condition on the right domain boundary is either a mirror when we simulate only a half of the system, or the wetting condition (∂ψ/∂x) r = -σ 0 (1 -ϕ 2 ) when the whole system is simulated. The imposed initial condition is a flat interface located at the altitude z = z m of the T = T m isotherm. Accordingly, the phase field is initially set to ψ(z) = (z m -z)/W 0 .

The two major difficulties in the process of rendering the phase-field simulations fully quantitative are:

(i) to achieve an accurate parametrization of the anisotropy functions a s and a k that relies on experimental data, (ii) to ensure a satisfactory convergence of the numerical results with decreasing interface width, W 0 → 0.

Point (i) has been discussed in the previous subsection for the parameters δ and δ 0 entering a s . The method used to deal with the parameter β 111 of a k , will be discussed in the next section. For point (ii), a systematic convergence method was recently introduced for the present model. The method consists in building a scaling curve that is valid for a wide range of growth velocities and for a wide range of W 0 values. This master curve is then used to extrapolate the side facet length Λ to the W 0 → 0 limit [START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF].

Physical parameters

Simulations are performed with the physical parameters gathered in table 1. A constant temperature gradient G is considered here: it represents a reasonable average value of the local interfacial gradients estimated experimentally. The variable geometrical parameter is the space dimension (2D or confined 3D) and the variable control parameters are the vertical growth velocity V and the misalignment angle θ m .

Physical parameter

Symbol

Value melting temperature T m 1687 K (ref. [START_REF] Boukellal | Isothermal solidification of Si and Al-Cu: 3D phasefield simulations[END_REF]) capillary length d 0 1.94 × 10 -4 µm (ref. [START_REF] Boukellal | Isothermal solidification of Si and Al-Cu: 3D phasefield simulations[END_REF]) thermal diffusivity D 2.3 × 10 7 µm 2 /s (ref. [START_REF] Nishi | Thermal diffusivities and conductivities of molten germanium and silicon[END_REF]) specific heat at constant pressure c p 1.032 J/g/K (ref. [START_REF] Kagaya | Specific heats of c-Si and molten Si[END_REF]) specific latent heat l H 1650 J/g (ref. [

Results of the comparison between experiment and simulation

We first investigate the confinement effects due to the limited sample thickness.

Dimensionality and confinement effects

Most of our simulations are performed in 2D systems while experiments use thin samples (300 µm in thickness). As the average experimental facet length is 600 µm, significant confinement effects are expected. They necessarily influence the solid-liquid interface shape as a whole. As a consequence, extracting an accurate value of the kinetic coefficient from local measurements on a side facet relies on a good knowledge of the confinement effects that affect the facet shape.

Shape of the solidification front

A few simulations of thin layers are first performed. In order to keep reasonable simulation times, the interface width W 0 must be increased. As discussed in [START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF] the results of these 3D simulations are not fully quantitative in this case. The boundary conditions imposed along the y direction at the front and the back domain boundaries differ. The reason is that in experiments, one crucible wall is a flat and atomically smooth piece of boron nitride (BN), while the other wall is carved into a similar piece of material by a mechanical tool that introduces surface roughness. It is well known that the roughness of the BN surface is an important factor affecting its wetting by silicon [START_REF] Drevet | Wetting and adhesion of Si on Si3N4 and BN substrates[END_REF]. We thus expect different contact angles on front and back walls. To incorporate these features into the simulations, we impose a wetting condition like (9) at the front boundary (smooth wall) and a mirror boundary condition at the back one (rough wall).

Let us first consider the side facets that appear on both x edges of the experimental samples. The interface equation that is applicable to a facetted 2D solid-liquid interface is the Ben Amar-Pomeau (BAP) equation [START_REF] Ben Amar | Growth of faceted needle crystals: theory[END_REF]. Generalizing this equation to 3D is not straightforward but one can argue that, provided the kinetic term is much larger than the capillary one, it should approximately read ū -

β 111 V f , (10)
where ū represents the dimensionless temperature u averaged over the whole facet area and V f = V sin θ f is the normal facet velocity. In the frozen linear temperature field imposed in the simulations,

ū = - c p l H G(H + f ūΛ cos θ f ). ( 11
)
Here H is the distance from the T = T m isotherm to the facet highest point, f ū a shape factor relating the location of the facet center of mass to its height, Λ the facet length observed from some remote point along the y direction, and θ f the crystallographic facet angle (see figure 2c). It was shown previously by a 2D analytical model that H can be neglected in first approximation [START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF]. As a result, the side facet length

Λ β 111 l H c p V f f ūG cos θ f (12)
increases linearly with the facet normal velocity V f , with a proportionality factor that depends on the facet geometry through its center of mass relative height f ū.

For the DS 011 seed orientation, the phase-field simulations show that the solid-liquid interface is rough, except on both sides of the sample where wetting conditions are used to obtain two {111} facets as observed in the experiments [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF][START_REF] Riberi-Béridot | Strain building and correlation with grain nucleation during silicon growth[END_REF]. Large rectangular-shaped facets are formed (top row in figure 1), due to the crystallographic octahedral structure of silicon where the y direction is the common edge of two {111} planes. This suggests that for the DS 011 orientation, 2D domains should be rather representative of 3D ones, thus that simulations results do not depend much on the domain thickness, as illustrated in figure 1. For the DS 001 orientation, side facets still exist but the interface is also largely faceted in the sample thickness because the apex of the octahedron is now common to four {111} facets (see bottom row in figure 1). In thin samples two triangular-shaped side facets pointing upward are formed along x and one rectangular-shaped facet (thickness facet) crosses the domain in the y direction. For both orientations, the side facets are perpendicular to the direction of observation (y axis). For this reason, the facet height H f = Λ cos θ f , where Λ is the apparent facet length observed through the sample. Then the facet center of mass approximatively lies at an altitude H f /2 below the apex of the rectangular facet and at an altitude 2H f /3 for a triangular facet pointing upward. This gives the values f ū(011) 1/2 and f ū(001) 2/3 for the shape factor introduced in equation ( 11).

In the case of the DS 001 seed orientation, due to the presence of a thickness facet, confinement effects are expected to increase with the growth velocity. At low velocities, the side facet shape is somewhat intermediate between a triangle and a rectangle while it becomes fully triangular-shaped at high velocities. The facet shape factor f ū(001) increases accordingly from 1/2 to 2/3. A confinement criterion is obtained by comparing the maximum length of the thickness facet, Λ c = L y / sin θ f with the unconfined facet length Λ(V f ) deduced from equation [START_REF] Karma | Phase-field formulation for quantitative modeling of alloy solidification[END_REF]

with f ū = 1/2. Strong confinement approximately corresponds to Λ > Λ c , that is to V f > V c , with V c = Gc p β 111 l H L y 2 cot θ f . ( 13
)
Using the physical parameters listed in table 1, one obtains V c 5.80 µm/s. Strong confinement is thus predicted in the experiments because V f is always larger than 7 µm/s. In the 3D phase-field simulations, we observe at high growth velocities that strong confinement may induce the alternance of two quasi-stationary states. The more stable one corresponds to side facets that are shorter than the thickness facet, as shown in the bottom row of figure 1c. The reason why this shape is more stable is that the BAP equation ( 12) is now simultaneously verified by the side and the thickness facets because Λ and f ū adjust to give the same average undercooling -ū = β 111 V f on both facets. Alternatively, in the case of the DS 011 seed orientation, our 3D phase-field simulations show that f ū(011) 1/2 independently of V f . Confinement effects are thus negligible for this orientation (top row in figure 1) so that 2D simulations are sufficient to describe the experimental system. Our next step is to determine an accurate estimate for the kinetic coefficient β 111 of the {111} facets, assuming a linear kinetic law.

Kinetic coefficient for a linear law

The standard version of TIPM is designed to quantitatively incorporate linear kinetic effects [START_REF] Karma | Phase-field method for computationally efficient modeling of solidification with arbitrary interface kinetics[END_REF][START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF], as confirmed by recent phase-field studies of silicon solidification based on this model [START_REF] Boukellal | Equilibrium and growth facetted shapes in isothermal solidification of silicon: 3D phasefield simulations[END_REF][START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF][START_REF] Lin | Three-dimensional phase field modeling of silicon thin-film growth during directional solidification: Facet formation and grain competition[END_REF]. On the other hand, attempts made to incorporate nonlinear kinetics showed that it is much more difficult to ensure quantitativeness in this case [START_REF] Liu | Phase field modeling of crystal growth with nonlinear kinetics[END_REF]. An estimate of the kinetic coefficient β 111 was previously obtained by plotting experimental values of the side facet normal velocity V f as a function of the groove undercooling ∆T g = Gd g , where d g is the total groove depth, and by fitting them to a second order kinetic law [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF]. In the present context, we fit them to a linear law in oder to allow for direct comparison with the linear kinetics imposed in our phasefield simulations. The experimental data points and the corresponding linear adjustment curves obtained for the DS 001 and DS 011 crystal orientations are shown in figure 3a. When the appropriate shape factor f ū is applied for each experimental batch of data, the two straight lines convincingly merge (figure 3b). Indeed, this is the expected result because all the facets on which measurements are performed are portions of {111} planes, independently of the sample crystal orientation. This implies that the quantitative 2D results for the lateral facet length Λ presented in the next subsection have to be corrected by a confinement factor f Λ = (1/2)/f ū prior to compare them with the experimental 3D ones because f ū = 1/2 for a 2D facet. As discussed recently [START_REF] Braik | Quantitativeness of phasefield simulations for directional solidification of faceted silicon monograins in thin samples[END_REF], 2D simulations performed for the DS 011 orientation conform more closely to the experimental thin layer geometry because there is a natural crystallographic invariance along the thickness y direction in this case. Based on this remark, we retain the slope of the experimental data points for the DS 011 sample in order to evaluate β 111 more accurately. A linear fit of these points gives

β 111 = 24.7 ± 4.3 × 10 -6 s/µm. ( 14
)
As discussed in a recent review [START_REF] Miller | Some remarks on the undercooling of the Si (111) facet and the Monte Carlo modeling of silicon crystal growth by[END_REF], several kinetic laws (not all linear) relating the {111} facet undercooling to its growth velocity have been proposed in the literature. To our knowledge, linear kinetic laws have been observed only for global solid-liquid interfaces of multicrystalline silicon [START_REF] Riberi-Béridot | Growth undercooling in multi-crystalline pure silicon and in silicon containing light impurities (C and O)[END_REF][START_REF] Fujiwara | In situ observation of Si faceted dendrite growth from low-degree-of-undercooling melts[END_REF][START_REF] Miyahara | Effect of twin growth on unidirectional solidification control of multicrystal silicon for solar cells[END_REF]. The corresponding β estimates are 78.2, 86.4, and 104.2 (×10 -6 s/µm), respectively. These values are of the same order of magnitude as our present estimate [START_REF] Singer-Loginova | The phase field technique for modeling multiphase materials[END_REF] but they are still significantly larger. This is probably due to the much higher undercoolings, up to 7 K, reached in the case of multicrystalline solid-liquid interfaces.

We now compare experimental and phase-field data for the variations of the side facet length as a function of the facet velocity.

Facet lengths

The side facet lengths Λ obtained in our 2D phase-field simulations are corrected for 3D thickness effects and compared to our experimental results in thin samples. As just discussed, according to the facet shape obtained in the simulations of thin 3D samples, the Λ values obtained in 2D simulations have to be multiplied by a confinement factor of either f Λ (011) = (1/2)/(1/2) = 1 or f Λ (001) = (1/2)/(2/3) = 3/4. The corresponding results are displayed in figure 4. For both growth directions, we obtain a good overall agreement between the experimental data and the 2D phase-field data corrected by the confinement factor f Λ (figures 4a and 4c).

We finally extend our experiment-simulation comparison to the case of a misalignment between the gradient direction and the crystal orientation imposed along the z axis.

Misalignment effect

In the experiments, the isotherms are slightly tilted as compared to perfectly horizontal isotherms. This is mainly due to the thermal configuration of the furnace as discussed in [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF][START_REF] Riberi-Béridot | Strain building and correlation with grain nucleation during silicon growth[END_REF]. As a result, while the isotherms remain basically planar, they are often rotated around the y axis by a misalignment angle θ m of a few degrees (figure 2b). When the temperature gradient is perfectly aligned with the crystal growth axis (θ m = 0), the groove depth d g associated to a facet of length Λ 0 is obtained from equation [START_REF] Karma | Phase-field formulation for quantitative modeling of alloy solidification[END_REF],

d g = Λ 0 cos θ f β 111 l H c p V f f ūG . ( 15
)
For θ m = 0, as illustrated in figure 5, the depth of the left facet groove measured along the gradient direction is d left = Λ left cos(θ f +θ m ), while d right = Λ right cos(θ f -θ m ) for the right facet groove, the facet angle θ f being the same on both sides. According to equation [START_REF] Yang | Crystal growth and equilibrium crystal shapes of silicon in the melt[END_REF], the groove depth measured along the gradient direction is a constant for a given value of V f , and one can thus conclude that

d g = d left = d right , (16) so
Λ 0 cos θ f = Λ left cos(θ f + θ m ) = Λ right cos(θ f -θ m ). ( 17
)
To study in more details the misalignment effect, we sorted out the experimental Λ(V f ) data in two groups corresponding to left and right facets. In parallel, we performed 2D phase-field simulations with a nonzero misalignment angle θ m to numerically determine Λ left and Λ right and draw a comparison between both approaches. The corresponding data are gathered in figure 6. For both seed orientations, the curves display a very good agreement between phase-field and experiment. Comparing with figures 4b and c, one see that the experiment-simulation agreement is definitely improved by taking misalignment into account.

To go further, one can remark that equations ( 16) and ( 17) impose a simple relation between the groove depths obtained at different misalignments θ m , since they are all related to d g , the groove depth for θ m = 0. In the following, d g (V f ) is estimated from our phase-field simulations performed without misalignment, and compared to d left (V f ) and to d right (V f ) both for the phase-field and the experimental data obtained with a nonzero misalignment. According to Eq. ( 16), plotting d left (V f ) and d right (V f ) versus d g (V f ) should cause the data points to collapse onto the identity line. This should be equally true for the growth directions [011] and [001] and for different values of the misalignment angle θ m . Figure 7 fully confirms this prediction, showing that the misalignment effect is geometrically based.

Conclusions

To accurately characterize and model the physical behavior of directionallysolidified thin silicon samples, it is appropriate to restrict their study to well-defined simple situations, namely here well defined seed orientation and monocrystal growth. In the present work, these conditions have been imposed both in experiments and in numerical simulations with the concern to allow a fine and relevant comparison of their results over a large range of growth velocities. This has lead us to systematically parametrize the phase-field model by available experimental data, with a particular attention devoted to the anisotropy functions of both the surface energy and the kinetic attachment coefficient. The results of 2D simulations have shown that useful information could be obtained from a planar description but that quantitative comparison with experiments also require 3D simulations to assess the influence of the sample thickness. Identifying and quantifying the resulting confinement effect have demonstrated that they turn from very weak to high when the seed orientation is changed from DS 011 to DS 001 . In the first case, 2D simulation results have proven representative of the thin experimental layers. In the latter case, the predicted increase of the facet length with the growth velocity opposed to the limited sample thickness has shown to provoke a strong faceting-confinement interplay, and the corresponding correction factor between 2D and thin layer geometries has been quantitatively estimated. The analytical law relating the facet length to the normal facet growth velocity in 3D confined and unconfined systems has been extended to include this correction factor. The origin of the small misalignment effect observed experimentally has been identified as geometrical and the experiment-simulation agreement has improved by taking misalignment into account.

Altogether, the present results bring valuable pieces of information and strong support to experimental studies of silicon solidification through in situ X-ray imaging. Moreover they also validate the proposed 3D phase-field model that uses anisotropy functions based on experimental studies. The remarkable experiment-simulation agreement obtained in the present study opens new perspectives to perform quantitative simulations of experimental silicon solidification in systems governed by more complex physical conditions.

Figure 1 :

 1 Figure 1: Cubic crystal axes defining the two seed orientations discussed in the text, DS 011 (top row) and DS 001 (bottom row). These crystal axes are labeled x, y, z in the paper. (a) Unit octahedron of the diamond crystal structure, (b) phase-field 3D growth shape in a domain with a square base, (c) phase-field 3D growth shape in a thin domain (L y = 300 µm). The facet growth velocities imposed in the simulations are V f = 4.89 µm/s (top row) and V f = 10.39 µm/s (bottom row).

Figure 2 :

 2 Figure 2: (a) Sketch of a solid-liquid interface showing a traversing facet: front view perpendicular to the y axis (as in the experimental radiograph) and side view perpendicular to the x axis. (b) Experimental in situ radiograph: the sample width is about 6 mm. Dashed lines parallel to the side facets are added to guide the eye. (c) Notations used in the text to describe a side facet.

Figure 3 :

 3 Figure 3: (a) Experimental estimates of the {111} side facet growth velocity V f along the facet normal as a function of the total groove undercooling ∆T g = Gd g , d g being the total groove depth. (b) Same data with a confinement correction factor f ū applied to the undercooling: f ū = 1/2 for the DS 011 crystal orientation and f ū = 2/3 for DS 001 . The symbol size is of the order of the experimental error.

Figure 4 :

 4 Figure 4: Variations of the {111} side facet length Λ with the facet normal velocity V f compared between experiments (crosses) and 2D phase-field simulations (full circles): (a) DS 011 crystal orientation, (b) DS 001 crystal orientation, (c) same as (b) with a correction factor f Λ (001) = 3/4 applied to the Λ values obtained in the simulations.

Figure 5 :

 5 Figure 5: Schematic representation of the misalignment between the local temperature gradient G and the vertical crystal direction (dashed vertical lines). The angle between the {111} side facet and the vertical z axis is denoted by θ f and the angle between G and z axis by θ m .

Figure 6 :

 6 Figure 6: Variations of the {111} side facet length Λ with the facet normal velocity V f compared between experiments (squares) and phase-field simulations (lines plus circles). Right and left facet data are respectively represented by open and closed symbols. For both growth directions, the results are obtained for a misalignment angle θ m = 5 • .

Figure 7 :

 7 Figure 7: Universal curve showing the collapse of the experimental (squares) and phasefield (other symbols) values of the groove depth d measured along the gradient direction. Values of d left (V f ) and d right (V f ) obtained for nonzero misalignment are represented as functions of the phase-field groove depth d g (V f ) obtained with perfect alignment. Empty (full) symbols are used for the right (left) {111} facets. Misalignment values are θ m = 3 • (asterisks), θ m = 7 • (plusses), and θ m = 5 • (other symbols). The straight line is the identity line.

Table 1 :

 1 Physical silicon parameters used in the phase-field simulations.

	33]

Acknowledgments

We would like to thank Mathis Plapp, Kozo Fujiwara and Klaus Kassner for interesting discussions about this work.

This research was partially funded by Agence Nationale de la Recherche (grant No. 14-CE05-0046-01, CrySaLID project).