Carolina Baruffi

Alphonse Finel

Yann Le Bouar

Brigitte Bacroix

O U Salman

Atomistic simulations of temperature-driven microstructure formation in pure Titanium

Keywords: titanium, martensitic phase transition, variant selection, atomistic simulations, overdamped Langevin dynamics

whether they are published or not. The documents may come

Introduction

Martensitic transformation (MT) is a particular sub-class of solid-to-solid structural phase transformations observed in many metals and alloys [START_REF] Bhattacharya | Microstructure of Martensite: Why it Forms and how it Gives Rise to the Shape-memory Effect[END_REF]. In 1 alphonse.finel@onera.fr Preprint submitted to Elsevier 6th August 2020 most general terms, the MT is a diffusionless displacive first-order phase transition. It involves a shear dominated change of shape in the underlying crystal lattice on alteration of the external conditions, i.e., temperature and/or pressure or stress. Very often, a very complex microstructure governed by the crystallographic symmetry of phases develops [START_REF] Khachaturyan | Some questions concerning the theory of phase transformations in solids[END_REF][START_REF] Schryvers | Microstructures and interfaces in Ni-Al martensite: comparing HRTEM observations with continuum theories[END_REF][START_REF] Delville | Transmission electron microscopy investigation of microstructures in low-hysteresis alloys with special lattice parameters[END_REF][START_REF] Finel | Phase field methods: Microstructures, mechanical properties and complexity[END_REF][START_REF] Salman | The role of phase compatibility in martensite[END_REF][START_REF] Salman | Modeling of spatio-temporal dynamics and patterning mechanisms of martensites by phase-field and lagrangian methods[END_REF][START_REF] Vattré | Polymorphism of iron at high pressure: A 3D phasefield model for displacive transitions with finite elastoplastic deformations[END_REF][START_REF] Salman | Origin of stabilization of macrotwin boundaries in martensites[END_REF], giving rise to exceptional mechanical properties such as shape-memory effect [START_REF] Hildebrand | An atomistic investigation of the kinetics of detwinning[END_REF], superelasticity [START_REF] Tanaka | Ferrous polycrystalline shape-memory alloy showing huge superelasticity[END_REF] and high-strength [START_REF] Saito | Multifunctional alloys obtained via a Dislocation-Free plastic deformation mechanism[END_REF][START_REF] Zhang | Complexionmediated martensitic phase transformation in titanium[END_REF].

More specifically, materials such as titanium, zirconium, and their alloys undergo martensitic phase transformation, and they are of strong interest for several applications in the nuclear [START_REF] Lemaignan | Zirconium alloys in nuclear applications[END_REF], aeronautic [START_REF] Benmhenni | Micromechanical modelling of monotonic loading of cp α-ti: Correlation between macroscopic and microscopic behaviour[END_REF] and bio-medical industries [START_REF] Leyens | Titanium and titanium alloys: fundamentals and applications[END_REF][START_REF] Brunette | Titanium in medicine: material science, surface science, engineering, biological responses and medical applications[END_REF]. These materials are characterized by a Body-Centered Cubic (BCC) structure, stabilized by vibrations at high temperature, which transforms into a Hexagonal Close Packed (HCP) structure when the temperature is lowered [START_REF] Heiming | Phonon dispersion of the bcc phase of group-iv metals. ii. bcc zirconium, a model case of dynamical precursors of martensitic transitions[END_REF].

Thus, during any conventional transformation route, e.g., metal forming [START_REF] Cabus | The influence of grain morphology on texture measured after phase transformation in multiphase steels[END_REF] or more advanced elaboration processes, e.g., additive manufacturing [START_REF] Chastand | Comparative study of fatigue properties of ti-6al-4v specimens built by electron beam melting (ebm) and selective laser melting (slm)[END_REF][START_REF] Thijs | Van Humbeeck, Strong morphological and crystallographic texture and resulting yield strength anisotropy in selective laser melted tantalum[END_REF][START_REF] Thijs | A study of the microstructural evolution during selective laser melting of ti-6al-4v[END_REF],

HCP→BCC and BCC→HCP transformations generally occur. Most importantly, both BCC and HCP are high-symmetry crystallographic phases without a group-subgroup relationship [START_REF] Conti | A variational model for reconstructive phase transformations in crystals, and their relation to dislocations and plasticity[END_REF][START_REF] Bhattacharya | Crystal symmetry and the reversibility of martensitic transformations[END_REF]. Therefore, the transformation involves necessarily large "strains" and thus, they can possibly be combined with other metallurgical and mechanical phenomena such as plasticity and induce significant modifications of the material microstructure and, as a result, of its mechanical properties. Consequently, optimizing these properties requires a clear understanding of the links between the processing parameters and the mechanisms involved in the microstructure formation and evolution. Although various microstructures are well documented experimentally, the debate is still open on topics such as the relationship between initial and final phases [START_REF] Cabus | Orientation relationship between austenite and bainite in a multiphased steel[END_REF][START_REF] Humbert | Evaluation of the orientation relations from misorientation between inherited variants: Application to ausformed martensite[END_REF][START_REF] Chen | Determination of the stretch tensor for structural transformations[END_REF], the variant selection criterion in bulk material and at grain boundaries [START_REF] Gao | Diffuse scattering as an indicator for martensitic variant selection[END_REF][START_REF] Wang | Effect of self-accommodation on α/α boundary populations in pure titanium[END_REF][START_REF] Lischewski | Nucleation and variant selection during the α-γ-α phase transformation in microalloyed steel[END_REF][START_REF] Shi | Variant selection during α precipitation in Ti-6Al-4V under the influence of local stress -a simulation study[END_REF], the exact kinetic and sequence of transformation events. . . Besides, phase transformation typically directly interacts with plasticity [START_REF] Conti | A variational model for reconstructive phase transformations in crystals, and their relation to dislocations and plasticity[END_REF][START_REF] Bhattacharya | Crystal symmetry and the reversibility of martensitic transformations[END_REF][START_REF] Otsuka | Fft-based modelling of transformation plasticity in polycrystalline materials during diffusive phase transformation[END_REF] or recrystallization [START_REF] Jedrychowski | Investigation of SIBM driven recrystallization in alpha zirconium based on EBSD data and monte carlo modeling[END_REF] and makes the understanding and prediction of the resulting microstructure even more complicated. Therefore, an atomic-scale simulation, that does not exclude any local mechanism and does not preselect a specific kinetic pathway, is crucial to gain insight into the details of the MT. In the present work, we investigate variant selection mechanisms in pure titanium through atomistic simulations with empirical potentials. In particular, we focus on the microstructure formation in pure titanium and on the influence of local mechanical constraints inducing a confinement effect on the material undergoing the transition. While many atomistic potentials based on the Embedded Atom Method (EAM) [START_REF] Mendelev | Development of an interatomic potential for the simulation of defects, plasticity, and phase transformations in titanium[END_REF], Modified Embedded Atom Method (MEAM) [START_REF] Hennig | Classical potential describes martensitic phase transformations between the α, β, and ω titanium phases[END_REF] and machine learning algorithms [START_REF] Takahashi | Conceptual and practical bases for the high accuracy of machine learning interatomic potentials: Application to elemental titanium[END_REF] have been suggested in the literature to describe plasticity and phase transition in titanium, the temperature-induced microstructure formation and evolution remain mostly unexplored in atomistic studies.

For our purpose, we test two different interatomic potentials to simulate the temperature-induced phase transformation under different stress conditions. Afterward, we analyze the morphology of the resulting microstructures through a deformation-gradient map representing the lattice distortion to elucidate the role of kinematic compatibility [START_REF] Falk | Dynamics of viscoplastic deformation in amorphous solids[END_REF][START_REF] Feng | Phase transformations and compatibility in helical structures[END_REF]. Moreover, in the simulations, we use a novel atomistic modeling approach based on the overdamped Langevin equations. This method has been recently proposed for the study of crystalline materials and successfully applied to grain boundary migration and compared with classical Molecular Dynamics [START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF][START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF]. Although mainly focused on the martensitic phase transition in titanium, the present work represents a further investigation of the potentialities of the method in describing structural evolution in crystalline materials.

STRETCH TENSORS U (k)

U (1) = 1 2      2η 1 0 0 0 η 2 + η 3 η 3 -η 2 0 η 3 -η 2 η 2 + η 3      U (2) = 1 2      2η 1 0 0 0 η 2 + η 3 η 2 -η 3 0 η 2 -η 3 η 2 + η 3      U (3) = 1 2      η 2 + η 3 0 η 3 -η 2 0 2η 1 0 η 3 -η 2 0 η 2 + η 3      U (4) = 1 2      η 2 + η 3 0 η 2 -η 3 0 2η 1 0 η 2 -η 3 0 η 2 + η 3      U (5) = 1 2      η 2 + η 3 -η 2 + η 3 0 -η 2 + η 3 η 2 + η 3 0 0 0 2η 1      U (6) = 1 2      η 2 + η 3 η 2 -η 3 0 η 2 -η 3 η 2 + η 3 0 0 0 2η 1     
= a a 0 , η 2 = 3 2 a a 0 , η 3 = c √ 2a 0
, where a 0 and (a, c) are the lattice parameters of the BCC ans HCP lattices, see Fig. 1.

Methods

Modelling approach

As mentioned in section 1, in our simulations, we use a recently introduced modeling approach describing the evolution of particle positions with an overdamped stochastic dynamics [START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF][START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF]. In contrast to Newtonian dynamics, particle positions are treated as stochastic variables which follow a first-order in time dynamics that do not explicitly incorporate high-frequency vibrations of the crystalline grid (phonons), which typically limit to few nanoseconds the time scale of classical Molecular Dynamics [START_REF] Paul | Molecular dynamics simulation, elementary methods[END_REF]. The chaotic nature of the Newtonian dynamics, which in the long time drives the system to a stochastic equilibrium state, is recovered in the first-order in time dynamics through the use of an additive noise term, carefully chosen to guarantee that the system converges to the correct thermodynamical state in the long-time limit. In this section, we report the equations used in the model, and we refer the reader to Refs. [START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF][START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF] for more details on its analytical derivation.

The main objective of the proposed approach is to avoid the time scale associated with phonons. Therefore, the configurational space is restricted to the coordinates x n i , where the upper index n = 1, ..., N refers to a particle and the lower index i = 1, 2, 3 to a cartesian coordinate. Correspondingly, the dynamics involves only the first derivatives of x n i .

dx n i dt = -ν -1 ∂Φ ∂x n i + Bη n i (t), (2.1)
where Φ({x n i }) is the potential energy between particles, ν a viscosity coefficient and B the amplitude of a white gaussian noise η n i (t) such that η n i (t) = 0, η n i (t)η m j (t ) = δ nm δ ij δ(t -t ). δ nm and δ ij are Kronecker symbols and δ(t -t ) stands for the Dirac delta distribution. Equations (2.1) represent a first-order in time stochastic dynamics, also known as overdamped Langevin Dynamics [START_REF] Nelson | Dynamical Theories of Brownian Motion[END_REF].

The coefficients ν and B are independent from particle positions and related through fluctuation-dissipation theorem B = √ 2k B T ν. This guarantees that, in the long-time limit t → ∞, the distribution probability P ({x n i }) generated by Eqs. (2.1) converges to a steady state characterized by the Boltzmann distribution:

t → ∞ : P ({x n i }) → A exp - Φ({x n i }) k B T (2.
2)

The dynamics represented by the set of equations (2.1) is valid in the (NVT) thermodynamical ensemble, i.e., the number of particle N, the volume V and the temperature T are fixed variables. To deal with applied stress conditions, we extended the model to the (NPT) ensemble, where P stands for the first Piola-Kirchhoff tensor. We present now briefly the stochastic dynamics required for the (NPT) ensemble. First we incorporate nine additional degrees of freedom into the model, which are the components of the deformation gradient F describing the change in shape of the simulation box. The particle coordinates within the simulation box are represented by scaled coordinates {x n i } related to the actual coordinates by:

xn i = H -1 ij x n j i = 1, 2, 3 (2.3)
where the matrix H is defined by H = FL 0 where L 0 is a diagonal matrix containing the length of the orthogonal vectors L 0 1 , L 0 2 and L 0 3 that define the initial simulation box. The extended overdamped Langevin dynamics reads as:

dx n i dt = -ν -1 ∂ H ∂ xn i + B η n i (t) i = 1, 2, 3 ; n = 1, ...N , dF ij dt = -γ -1 ∂ H ∂F ij + A ξ ij (t) i, j = 1, ..., 3 , (2.4)
where ξ ij (t) is a white gaussian noise, γ a viscosity associated with the new degrees of freedom F ij and H({x n i }, F) the Hamiltonian for the extended set of DOF. The extended Hamiltonian should of course be such that Eqs. (2.4) converge in the long time limit towards the thermodynamical equilibrium of the (NPT) ensemble. As shown in [START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF], this leads to:

H({x n i }, F) = Φ(F ij L 0 j xn j ) + V 0 P ij F ij -N k B T ln (V 0 det F) (2.5)
where V 0 is the volume of the initial simulation box. Finally, to numerically evaluate Eqs. (2.4), we used an explicit predictor-corrector method (see [START_REF] Burrage | Numerical methods for second-order stochastic differential equations[END_REF][START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF]).

Simulation setup and interatomic potential

We simulate the BCC→HCP transition in both thermodynamic ensembles (NVT) and (NPT) to clarify the influence of different external conditions on martensite microstructure. In the latter, we consider stress-free boundaries in each direction by setting each component of the first Piola-Kirchhoff tensor P to zero in all directions, which allows the material to change its macroscopic shape. Although real conditions experienced by a region in bulk material would be an intermediate case between these two conditions, the two extreme scenarios are useful for a global understanding of the influence of local mechanical constraints preventing a free change in shape and/or volume of the matrix around a martensite nucleus. In the simulations, we first allowed a BCC structure to reach its equilibrium state at 1400 K followed by a subsequent fast quenching at 700 K. We perform the quenching by instantaneous rescaling of the temperature parameter that fixes the noise amplitudes in Eqs. 2.4. The simulation box size is set equal to 36 × 36 × 36 a 3 0 , where a 0 is the BCC equilibrium lattice constant at 1400 K, see Fig. 1. The total number of atoms is 93312. We apply periodic boundary conditions in all directions.

Results

We take special care in the choice of the potential to model interatomic forces. Developing empirical potentials suitable to reproduce the properties of a given material is still a difficult issue in the field of atomistic modeling.

Although many-body potentials such as the EAM [START_REF] Daw | Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals[END_REF][START_REF] Daw | Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals[END_REF] have paved the way to a faithful description of several metals, modeling of materials with strong bond directionality (such as covalent materials or metals with partially full-d shell)

is still a challenge and typically requires more complex functional forms. In the present work, we performed preliminary simulations at the Molecular Dynamics scale, using LAMMPS [START_REF] Plimpton | Fast parallel algorithms for Short-Range molecular dynamics[END_REF] with two different empirical potentials, an EAMtype potential [START_REF] Mendelev | Development of an interatomic potential for the simulation of defects, plasticity, and phase transformations in titanium[END_REF] and a MEAM [START_REF] Hennig | Classical potential describes martensitic phase transformations between the α, β, and ω titanium phases[END_REF]. We obtained the following results: when the MEAM-type potential is used, we were able to observe a stable BCC phase transforming into HCP upon quenching. On the other hand, we did not observe a phase transition after cooling when the EAM potential is used, although we were able to get a stable BCC structure at high temperatures. Therefore, we increased the simulation duration up to 1 nanosecond during annealing of the material, and we tested different simulation box sizes. However, the transition did not occur in any of the two thermodynamic ensembles. Possible reasons for that could be the simpler functional form of the EAM potential compared to the MEAM (i.e., the lack of any angular dependency in the embedding term describing the electron density) or the presence of a high energy barrier for the nucleation of the HCP phase. Based on these results, we finally decided to use the MEAM potential to perform the overdamped Langevin simulations and implemented it in a parallel code that we developed by following the previous work on many-body force field implementation [START_REF] Goedecker | Optimization and parallelization of a force field for silicon using openmp[END_REF].

Variant and phase identification

To characterize microstructure formation and evolution occurring in the simulations, we use a deformation gradient map to represent the local lattice distortion. Indeed, when a material undergoes a martensitic transformation, several energetically equivalent variants differing in their relative crystallographic orientation appear. Each of these variants is associated with a stretch tensor U that can be easily identified once the local deformation gradient F is known:

we just need to use the polar decomposition F = QU (where Q is a rotation and U = U T is positive-definite), which is unique. However, due to the infinite degeneracy of the lattice groups of the parent and product phases, the identification of the deformation gradient F is not unique. Therefore, it is common to use a lattice correspondence between the parent and product phases which is meant to represent the actual lattice sites displacements. Before to explain the procedure used below to identify F, we first recall that the martensitic transformation from BCC to HCP cannot be fully described by a simple homogeneous deformation gradient: the BCC lattice is a Bravais lattice, the HCP is not. Therefore, the transformation strain that we want to identify must be supplemented by atomic displacements applied on a sub-lattice of the deformed lattice. In the present situation, this shuffling consists in translating every second basal plane of the hexagonal lattice obtained after the action of the homogeneous deformation gradient.

We now turn to the procedure used to identify the local deformation gradient. This identification requires a lattice correspondence between the parent and product phases. Two lattice correspondences, given in terms of orientation relationships, have been proposed. The mechanism given by Burgers [START_REF] Burgers | On the process of transition of the cubic-body-centered modification into the hexagonal-close-packed modification of zirconium[END_REF] states the following correspondence between crystallographic planes and directions:

(110) bcc || (0001) hcp ; [ 111] bcc || [ 2110] hcp , (3.1)
whereas the mechanism given by Mao [START_REF] Mao | Effect of pressure on crystal structure and lattice parameters of iron up to 300 kbar[END_REF] states the following correspondence:

(110) bcc || (0001) hcp ; [00 1] bcc || [11 20] hcp . (3.2)
The two mechanisms differ only in that the Burgers mechanism requires a rotation of ±5.

n . For this purpose, we note that this stretch tensor is crystallographically linked to the (110) BCC plane that will transform into the subsequent (0001) HCP plane.

Therefore, for each of the six (110) BCC planes, we define a set of neighboring sites that is specific to this (110) plane (see Appendix for more details). Then, for each atom n, we identify six local deformation gradients F (k) n , k = 1 to 6, that minimize the following local descriptors:

k = 1 to 6 : D (k) 2 n = m∈Ωn ∆r nm (t * ) -F (k) n ∆r nm (0) 2 (3.3)
where Ω n is the neighborhood set that is associated with a given (110) plane.

The local deformation gradient F n is defined as the one that, among the six tensors F

Simulations in the (NPT) ensemble with stress-free boundary conditions

After equilibrating the system in the BCC state at T=1400 K, we abruptly drop the temperature down to T=700 K. After a short relaxation, the system transforms into a HCP structure. We report in Fig. 2a the evolution of the BCC and HCP phase fractions during the transition. Almost no BCC phase is left after the transformation has completed. A non-negligible residual fraction of atoms (≈ 0.20) exhibits crystallographic symmetry different than HCP, which suggests that some defects are generated. We comment on that point below. In Fig. 2b, we report the evolution of variant fractions as a function of time. At the very beginning of the transition, all the six variants nucleate almost instantaneously. However, very quickly, most of them disappear, giving rise to a microstructure composed of the variants U (2) and U (4) . Afterwards, the structure coarsens further and a single variant U (2) domain forms (see Fig. 3). However, the snapshots shown in Fig. 4 put in evidence that when the microstructure is coarsening, HCP domains with the same c axis orientation but different shuffling directions (referred to as a couple of "anti-variants", see for instance [START_REF] Gao | Diffuse scattering as an indicator for martensitic variant selection[END_REF]) come into contact and generate an anti-phase boundary (see the inset Fig. 4b). This boundary is equivalent to a stacking fault in case the boundary plane between the two domains is parallel to the {0001} HCP basal plane.

Simulations in the (NVT) ensemble

After quenching the system down to T=700 K, the crystal structure transforms into an HCP structure without any remaining BCC domain, similarly to what is observed in the previous (NPT) simulation.

We report in Fig. 5a the evolution of the BCC and HCP volume fractions and, in Fig. 5b, the evolution of the volume fractions of the six variants of the HCP phase. As observed, the overall transformation proceeds through welldefined stages. We first observe a nucleation stage, that extends up to point A in Fig. 5a and5b, during which local HCP fluctuations emerge. However, because of the deep quench, the length scale of these fluctuations is too small to identify the nuclei of the local HCP variants. Indeed, the procedure used to identify HCP variants relies on neighborhood that extend beyond the second neighbor shell, whereas the PTM algorithm used to identify the local lattice relies on a neighboring set limited to the first two neighbor shells [START_REF] Larsen | Robust structural identification via polyhedral template matching[END_REF]. Next, as seen in Fig. 5, the nucleation stage switches rapidly between points A and B and, between points B and C, stabilizes to a quasi-stationary stage during which the six HCP variants reach finite volume fractions that are roughly constant. Then, the system enters a stage, that extends between points C and D, during which the volume fraction of three HCP variants increases at the expense of the three others. This growing stage ends at point D beyond which the system stabilizes in a microstructure that consists of only three orientational variants.

In the simulation shown, the selected variants share the [111] BCC direction in the parent phase, i.e., a 11 20 HCP direction. The resulting microstructure

Discussion

The simulation results in the (NVT) and (NPT) ensembles highlight how local mechanical constraints deeply influence the martensite morphology. In both the ensembles, at the really beginning of the transition all the variants appear, which is expected since all the variants are energitically equivalent and therefore accessible by the system due to the randomness of thermal fluctuations [START_REF] Gao | Diffuse scattering as an indicator for martensitic variant selection[END_REF]. However, when the structure further evolves, some of them are selected and constitues the final microstructure. We observed two distinct microstructural evolutions as a function of the applied boundary conditions. For simulations in the (NPT) ensemble, the microstructure coarsens and forms a single variant domain with anti-phase defects. On the other hand, for simulations in the (NVT) ensemble, a triplet of variants with a common 11 20

HCP direction is systematically selected. As discussed in section 3, in almost all the simulations performed, the 3 selected variants cluster in a 3-plate geometry around the common dense direction. The same morphology has been recently observed in pure titanium [START_REF] Farabi | Five-parameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF] and previously documented in Ti-Nb shape memory alloys [START_REF] Chai | Self-accommodation in ti-nb shape memory alloys[END_REF] and zirconium alloys [START_REF] Srivastava | Self accomodation morphology of martensite variants in zr 2.5 wt% nb alloy[END_REF]. Moreover, according to the phenomenological theory of martensite [START_REF] Bowles | The crystallography of martensite transformations i[END_REF], the 3-variant cluster minimizes the overall mesoscopic shape strain [START_REF] Farabi | Five-parameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF][START_REF] Srivastava | Self accomodation morphology of martensite variants in zr 2.5 wt% nb alloy[END_REF][START_REF] Wang | Effect of self-accommodation on α/α boundary populations in pure titanium[END_REF] so that this morphology should appear when strain energy minimization is dominant in driving the evolution of the microstructure and variant selection is governed by self-accommodation [START_REF] Miyazaki | The shape memory mechanism associated with the martensitic transformation in ti-ni alloys-i. selfaccommodation[END_REF][START_REF] Pitteri | Continuum models for phase transitions and twinning in crystals[END_REF].

Besides, the agreement with the microstructures observed in experiments shows that the fixed-volume conditions adopted in our simulations are well adapted to reproduce local constraints experimented by real systems. Furthermore, we also occasionally observe that the 3 selected variants form two laminates with a "buffer" FCC domain at the crossing point. While in the simulation showing a 3-plate morphology the variants are selected with a similar volume fraction (≈ 0.3), in this case one variant (the one participating in both the laminates) is dominant with respect to the others (see Fig. 5). The possible presence of FCC phase after transition has little experimental evidence [START_REF] Nishiyama | Transmission electron microscope study of the martensites in a titanium-3 wt% iron alloy[END_REF]. However, it has been already documented for zirconium by means of Molecular Dynamics simulations [START_REF] Morris | Molecular dynamic simulation of a homogeneous bcc → hcp transition[END_REF][START_REF] Pinsook | Simulation of martensitic microstructural evolution in zirconium[END_REF][START_REF] Ackland | Molecular dynamics simulations of the martensitic phase transition process[END_REF].

In terms of interfaces, simulations in the (NPT) ensemble (i.e., one variant morphology is not exactly compatible with the formation of three {10 11} boundaries. Indeed, the basal planes of the 3 selected variant are mutually rotated by exactly 60 • around the common 11 20 HCP direction (see Fig. 6 and9) but the {10 11} plane form an angle of 61.5 • with the basal plane. Consequently, some further accommodation is required to form triple junctions like those shown in the detail of Fig. 6. Nevertheless, analysis of frequency distributions for the misorientation angles in pure titanium [START_REF] Wang | Effect of self-accommodation on α/α boundary populations in pure titanium[END_REF][START_REF] Farabi | Five-parameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF] and in titanium alloys [START_REF] Beladi | Variant selection and intervariant crystallographic planes distribution in martensite in a ti-6al-4v alloy[END_REF] also highlight a strong preference for the formation of high angle boundaries with 60 • of misorientation around the 11 20 axis and grain boundary plane close to the {10 11}. A possible conclusion is that the energy cost linked to the additional strain required to form the triple junction in the 3-plate clusters is much lower than the energy gain in forming three coherent interfaces that, as shown in Molecular Dynamics simulations [START_REF] Wang | Atomic structures of symmetric tilt grain boundaries in hexagonal close packed (hcp) crystals[END_REF], display the smallest energy among between the other possible boundaries. Finally, we discuss the numerically computed stretch tensors used to identify variants. Fig. 10 shows the histograms of the diagonal and off-diagonal components of U n in the simulations in the (NPT) and (NVT) ensemble where dotted lines correspond to the values of the theoretical strain predicted by both the Mao and Burgers mechanisms. We note that for the simulation in the (NPT) ensemble, there is an overall good agreement between average values and theoretical predictions. On the other side, in fixed-volume conditions, the deviation is higher and more pronounced in the off-diagonal coefficients. Fig. 11 compares the histograms of the off-diagonal coefficients for the simulation in the (NVT) ensemble at three different instants t 1 , t 2 , t 3 , corresponding to the end of the nucleation stage, the quasi-stationary regime and the final stationary state (see Fig. 5). Unlike the histograms shown in Fig. 10, we report the value for only one of the three variants forming the final microstructure: U (4) . The comparison between histograms highlights how the divergence from the theoretical transformation matrix (shown in dotted line) starts developing after the nucleation stage, when almost all the BCC phase has disappeared and the HCP structure begins coarsening towards the final 3-plate morphology. This is par-ticularly evident for the off-diagonal coefficient U 13 . Consequently, although further investigations are needed, the deviation seems related to additional deformation mechanisms arriving in fixed-volume conditions after the first phase of nucleation and growth of HCP domains and linked to the textural evolution of martensite.

Conclusions

To summarize, we performed a set of extended overdamped Langevin dynamics simulations to study the microstructural evolution in pure titanium after identifying an appropriate interatomic potential through preliminary MD calculations. We then implemented an algorithm to construct a local deformation gradient map to characterize the lattice distortion that allowed us to clarify how crystal symmetry and kinematics influence defect and microstructure formation. Our main results concern the influence of macroscopic constraints on the transition such that a simple mono-variant domain composed of anti-phase domain boundaries forms when the crystal is allowed to change its shape. In contrast, a poly-variant microstructure develops when the shape of the crystal is constrained. This occurs because of the fundamental role of kinematic compatibility in structural phase transitions that imposes, at mesoscale, the formation of strain-energy minimizing structures that we were able to identify through our variant identification algorithm. A natural extension of this work will be to investigate how final microstructures obtained here influence the mechanical response of the material under external mechanical loading. This can easily be performed in our formulation through controlling the components of the Piola-Kirchhoff tensor (stress controlled) or the deformation gradient (strain controlled). Our findings can also be useful to develop appropriate mesoscale phase-field theories of BCC-HCP transition formulated using finite strains [START_REF] Vattré | Polymorphism of iron at high pressure: A 3D phasefield model for displacive transitions with finite elastoplastic deformations[END_REF][START_REF] Denoual | A phase field approach with a reaction pathwaysbased potential to model reconstructive martensitic transformations with a large number of variants[END_REF] and Landau-type theories with strain components used as the order parameter [START_REF] Shchyglo | Martensitic phase transformations in Ni-Ti-based shape memory alloys: The landau theory[END_REF].

Finally, the current work, which follows a previous numerical application lim-ited to a two-dimensional case [START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF], is the first application of the overdamped Langevin formalism to study a microstructural evolution problem in a fully three-dimensional setup. As shown above, even though it does not, on purpose, incorporate lattice vibrations, this dynamics leads to reliable results, without being limited by the time scale associated with phonons. However, as argued previously in Ref. [START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF] an explicit coarse-graining procedure over the initial Newtonian dynamics is required to obtain an exact formulation of the overdamped Langevin equations to reach much larger time-scales than MD. This will be subject of another work.

As already mentioned, the transformation from BCC to HCP cannot be fully described by a simple homogeneous deformation gradient but supplementary atomic displacements applied on a sublattice of the deformed lattice are needed. These displacements consist in an alternate shuffling of {110} planes along the 110 directions. The shuffling is not described by the overall deformation of the lattice so, to numerically evaluate the local strain, half of the atoms of the original BCC lattice must be considered in the six possible configurations.

After defining the six possible sets of neighbors, we calculate for each atom n six deformation gradients F

Figure 1 :

 1 Figure 1: BCC and HCP lattices and their lattice constants.

Figure 2 :

 2 Figure 2: Evolution of phase fractions (a) and variant fractions (b) in the (NPT) ensemble.

n , k = 1 to 6 .

 6 Finally, a polar decomposition leads to the local stretch tensor U n and, therefore, to a local stretch deformation map. The non-affine displacement D (k) 2 n quantifies the degree at which an affine transformation can describe the local change in the lattice. In the following analysis, we set a threshold D 2 lim = 6.5 Å 2 above which the calculated F n is considered not meaningful and exclude from post-processing atoms with D (k) 2 n> D 2 lim . Also, to monitor the evolution of the phase fraction of each phase without distinguishing variants, we use the Polyhedral Template Match analysis (PTM)[START_REF] Larsen | Robust structural identification via polyhedral template matching[END_REF]. This method classifies crystal structures according to the topology of the local atomic environment. It provides a flexible tool for structural identification even in presence of strong thermal fluctuations when other methods relying on interatomic distances (e.g., Common Neighbor Analysis[START_REF] Honeycutt | Molecular dynamics study of melting and freezing of small lennard-jones clusters[END_REF]) are less robust. In our analysis, the cut-off for the Root-Mean-Square-Deviation (RMSD) between the local atomic structure and the ideal structural template has been set equal to 0.14 Å.

Figure 3 :

 3 Figure 3: Simulation in the (NPT) ensemble, only atoms classified as HCP are shown and colored on the basis of the corresponding variant: a) initial nucleation stage, (b)-(c) two variants prevail, (d) final single variant domain. Crystallographic directions refer to the parent BCC phase.

Figure 4 :

 4 Figure 4: Final state in the (NPT) ensemble: atoms with crystallography different from HCP are colored in black: a) anti-variant domains (pink) separated by another variant (green), b) anti-phase defect (see the inset) formed at the boundary between anti-variant domains after microstructure coarsened. Crystallographic directions refer to the parent BCC phase.

Figure 5 :

 5 Figure 5: Evolution of the volume fraction of BCC and HCP phases (a), and of variants (b) in the (NVT) ensemble.

Figure 6 :

 6 Figure 6: Final microstructure obtained at 700 K in the (NVT) ensemble. The inset shows a triple junction detail. Atoms classified as HCP are colored according to the corresponding variant and atoms with crystallographic symmetry different than HCP in black. Crystallographic directions refer to the parent BCC phase.

Figure 7 :

 7 Figure 7: Evolution of the microstructure in the (NVT) ensemble, only atoms in a small slab normal to the [111] BCC direction and classified as HCP are shown and colored according to the corresponding variant: a) nucleation stage (from A to B in Fig. 5a), all the variants appear, b) quasi-stationary regime, two stable triple junction, highlighted by arrows, are identifiable,(c-d) the microstructure coarsens (from C to D in Fig.5a) in a 3-plate morphology (final stable state). Crystallographic directions refer to the parent BCC phase.

Figure 8 :

 8 Figure 8: Evolution of the variant fraction in four different simulations a)-d) in the (NVT) ensemble where the random noise term is different in each case.

Figure 9 :

 9 Figure 9: Final microstructure at 700 K in the (NVT) ensemble, highlighting the crossing between laminates. Crystallographic directions refer to the parent BCC phase.

Figure 10 :

 10 Figure 10: Histograms of Un diagonal and off-diagonal coefficients for a)-b) simulations in the (NPT) and c)-d) (NVT) ensemble. The dotted line shows theoretical values.

Figure 11 :

 11 Figure 11: Histograms of off-diagonal for a single variant U (4) n in the (NVT) ensemble at three different instants (point A,B,C in the HCP fraction curve of Fig. 5). Theoretical values are shown in dotted line.

n

  by following the approach proposed by Falk[START_REF] Falk | Dynamics of viscoplastic deformation in amorphous solids[END_REF] and, by polar decomposition, six strains U (k) n each one associated to a {110} BCC plane in the undeformed configuration. The local strain for atom n is then defined as the U (k) n with minimum D (k) 2n . Based on this assignment, we label the atom n as belonging to the corresponding variant.

Figure 12 :

 12 Figure 12: a) the six cubic cells which can deform in the final orthorhombic cell are highlighted with different colors (note that central atoms are not shown for sake of simplicity); b) an example of neighbor set (colored in black) for a given atom n (colored in red) for one of the six configurations considered.

Table 1 :

 1

The six transformation stretch tensors associated with the BCC→HCP displacive transformation, reported in the direct orthonormal basis aligned with the cubic directions of the BCC lattice. Coefficients η 1 , η 2 and η 3 are related to the BCC and HCP lattice parameters by η 1

Appendix

The numerical procedure used to identify the different variant is based on the definition of a local strain describing the transition from the BCC to the HCP structure. To calculate this local atomic strain we implemented the following procedure as described in [START_REF] Wu | The variant identification in molecular dynamics simulations for shape memory alloys[END_REF].

We start from a BCC structure with crystal axis 100 BCC parallel to the main frame axis. For every atom n, we consider six possible sets of 14 neighbors by taking six different configurations defined on the basis of the six cubic cells which can deform into the orthorhombic one, as schematically illustrated in Fig. 12.