Carolina Baruffi

Alphonse Finel

Yann Le Bouar

Brigitte Bacroix

Oguz Umut

Alphonse Baruffi

Yann Finel

Brigitte Le Bouar

Oguz Bacroix

Salman Umut

Atomistic

O U Salman

Atomistic simulation of martensite microstructural evolution during temperature driven β → α transition in pure titanium

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

Introduction

Martensitic transformation (MT) is a particular subclass of solid-to-solid structural phase transformations observed in many metals and alloys [START_REF] Bhattacharya | Microstructure of Martensite: Why it Forms 735 and how it Gives Rise to the Shape-memory Effect[END_REF]. In 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. Often, a complex microstructure governed by the symmetry of the different 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 | Transmis-745 sion electron microscopy investigation of microstructures in lowhysteresis 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 755 patterning mechanisms of martensites by phase-field and lagrangian methods[END_REF][START_REF] Vattré | Polymorphism of iron at high pressure: A 3D phase-field 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 highstrength [START_REF] Saito | Multifunctional alloys obtained via a Dislocation-Free plastic deformation mechanism[END_REF][START_REF] Zhang | Complexion-mediated martensitic phase transformation in titanium[END_REF]. In particular, materials of strong interest for 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], such as titanium, zirconium and their alloys, undergo martensitic phase transformation. The present paper specifically deals with martensitic phase transition in pure titanium. However, the observations and conclusions of our work are potentially valid, at least in general terms, also for elements such as zirconium and other alloys undergoing the BCC→HCP transition. Titanium exists as Hexagonal Close Packed (HCP) phase 1 alphonse.finel@onera.fr (α-phase) at room temperature and atmospheric pressure. On raising the pressure, while keeping the temperature constant, it transforms to a hexagonal phase (ω phase) 25 around 2 GPa pressure [START_REF] Banerjee | Phase transformations: examples from titanium and zirconium alloys[END_REF]. When the temperature is raised at atmospheric pressure conditions, the HCP structure transforms to a Body-Centered Cubic (BCC) structure (β phase) at 1155 K, stabilized by lattice vibrations [START_REF] Heiming | Phonon dispersion of the bcc phase of 795 group-iv metals. ii. bcc zirconium, a model case of dynamical precursors of martensitic transitions[END_REF]. It is clear that, during any conventional transforma-30 tion route (e.g., metal forming [START_REF] Cabus | The influence of grain morphology on texture measured after phase transformation in mul-800 tiphase steels[END_REF]) or 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 | 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]), the transition between the three possible phases can generally occur multiple times, combined with other metallurgical and mechanical phenomena such as plasticity [START_REF] Conti | A variational model for reconstructive 815 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] Furuhara | Variant selection in heterogeneous nucleation on defects in diffusional phase transformation and precipitation[END_REF] 26] or recrystallization [START_REF] Jedrychowski | Investigation of SIBM driven recrystallization in alpha zirconium based on EBSD data and monte carlo modeling[END_REF]. This fact induces significant modifications of the material microstructure and, as a result, of its mechanical properties. Optimizing these properties requires a clear understanding of the correlations between the processing parameters and the mechanisms 40 involved in the microstructure formation and evolution. Despite the large number of numerical and experimental studies, the debate is still open on topics such as the relationship between initial and final phases [START_REF] Cabus | Orientation relationship 830 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 45 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][START_REF] Furuhara | Variant selection in heterogeneous nucleation on defects in diffusional phase transformation and precipitation[END_REF], the exact kinetic and sequence of transformation events [START_REF] Srivastava | Evolution of microstructure during fabrication of zr-2.5 wt pct nb alloy pressure tubes[END_REF][START_REF] Farabi | Fiveparameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF], the structure of interfaces between different variant domains [START_REF] Banerjee | Substructure in titanium alloy martensite[END_REF][START_REF] Matsuda | Crystallography and morphology of antiphase boundary-like structure induced by martensitic transformation in ti-pd shape memory alloy[END_REF][START_REF] Otsuka | Physical metallurgy of ti-ni-based shape memory alloys[END_REF]. The present work aims to partially fill this gap of knowledge in the case of the temperature-driven BCC→HCP phase transition in pure titanium and, in particular, performs a computational investigation on martensite microstructural evolution during the transition and its final morphologies growing under different stress conditions. To our knowledge, this is the first atomistic study investigating such processes in pure titanium. We employ atomistic modeling to simulate the transition under different stress conditions. Contrary to microstructural models written at the mesoscale, this modeling technique does not either exclude any local mechanism or preselect a specific kinetic pathway. Thus, it is well suited to gain insight into the details of the transition. In particular, we used a recently developed atomistic approach based on the overdamped Langevin 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], which we here apply to simulate a displacive solid-state phase transition.

As opposed to experiments, which are mostly limited to the observation of final microstructures, modeling allows to follow the entire microstructural evolution during the transition. We analyse the sequence of nucleation events that subsequently lead to long-stage microstructures emergence under different stress conditions. The analysis of martensite morphologies resulting form the transformation confirms the importance of different stress conditions in the variant selection process and, consequently, on the appearance of different defects (e.g., intervariant boundaries and antiphase defects).Two types of interfaces are mainly observed. The first ones are a consequence of the shuffle degeneracy within each orientational variant. As they do not involve any elastic accomodation, these interfaces are essentially wavy. The second ones result from the impingement between different orientational variants and, consequently, display plane morphologies due to longrange elastic relaxation.

Methods

Modeling approach

We employed 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] to model the BCC→HCP transition under different stress conditions. With this method, particle positions are treated as stochastic variables that follow a first-order in time dynamics that do not explicitly incorporate high-frequency vibrations of the crystalline grid (phonons), which limits the time scale of classical Molecular Dynamics to a few nanoseconds [START_REF] Paul | Molecular dynamics simulation, elementary methods. by jm haile[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. While this approach has been widely used in studying the dynamics of soft matter systems as well as in biomolecular simulations [START_REF] Ando | Multiple time step brownian dynamics for long time simulation of biomolecules[END_REF][START_REF] Manghi | Hydrodynamic effects in driven soft matter[END_REF][START_REF] Ma | From generalized langevin equations to brownian dynamics and embedded brownian dynamics[END_REF], it has, to our knowledge, never been employed for the simulation of crystalline materials. In this section, we report the main equations used in the model. The 105 reader may refer 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. In the proposed approach, 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 110 to a cartesian coordinate. Correspondingly, the dynamics involves only the first derivatives of x n i and reads as

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

115 Gaussian noise η n i (t) such that η n i (t) = 0, η n i (t)η m j (t ) = δ nm δ ij δ(t -t ).
δ ij 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]. We 120 hypothesize that the coefficients ν and B are independent from particle positions and related by a fluctuationdissipation relation, i.e., B = √ 2k B T ν -1 . This guarantees that, in the long-time limit t → ∞, the distribution probability P ({x n i }) generated by Eqs. (2.1) converges to a 125 steady state characterized by the Boltzmann distribution

t → ∞ : P ({x n i }) → A exp - Φ({x n i }) k B T . (2.2)
Supplemented by the constraint that the particles stay within a pre-defined simulation box, the dynamics represented by the set of equations (2.1) is valid in the (N V T ) 130 thermodynamical ensemble, i.e., the number of particle N , the volume V and the temperature T are fixed. To deal with applied stress conditions, we extended the model to the (N P T ) ensemble, where P stands for the first Piola-Kirchhoff tensor. We present now briefly the stochastic 135 dynamics required for this (N P T ) 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 the shape of the simulation box. Next, in order to couple F to the degrees 140 of freedom associated to the atomic positions, we introduce 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 = F L 0 , where L 0 is 145 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 ,

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  
η 2 = 3 2 a a 0 , η 3 = c √ 2a 0
, where a 0 and (a, c) are the lattice parameters of the BCC and HCP lattices.

as

dx n i dt = -ν -1 ∂ H ∂ xn i + 2k B T ν -1 η n i (t) i = 1, 2, 3 ; n = 1, ...N , dF ij dt = -γ -1 ∂ H ∂F ij + 2k B T γ -1 ξ ij (t) i, j = 1, ..., 3 , (2.4)
where ξ ij (t) is a white Gaussian noise with ξ ij (t) = 0 and ξ ij (t)ξ kl (t ) = δ ik δ jl δ(t -t ), γ 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 (N P T ) 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. We note that this extended Hamiltonian is equal to the usual enthalpy 2 (Φ + V 0 P F ) supplemented by an extra logarithmic term, which ensures that the Langevin dynamics converges towards the correct thermodynamic equilibrium [START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF]. In principle, this logarithmic term should include a normalization volume unit V B that would make the argument of the logarithmic term dimensionless but, as long as we are considering kinetics and the associated driving forces, this normalizing term plays no role 3 and, therefore, is not included in Eq. (2.5). Altogether, the use of the extended Hamiltonian within Eq. (2.4) guarantees that the dynamics converges to the correct equilibrium state 2 Note that our sign convention is such that P ij = δ ij P where P > 0 corresponds to a system under hydrostatic pressure. 3 In principle, the logarithmic term in Eq. (2.5) should be written as log (V 0 /V B det F ) where the quantity V B , which has the dimension of a volume, is given by V B = Λ 3 where Λ is the Broglie wavelength. This term is reminiscent of the quantum and, therefore, discrete nature of the Hamiltonian and should be incorporated if we want to ensure a proper normalization of the entropies and energies. However, as it enters only through a constant term that disappears through derivation, this normalizing term does not enter the driving forces and plays no role within the dynamics. This is the reason why we do not include it in Eq. (2.5), even though it would make the argument of the logarithmic term dimensionless. associated with the (N P T ) ensemble in which an extern-175 ally applied first Piola-Kirchhoff stress controls the system. We stress that the appearance of the first Piola-Kirchhoff stress is simply a consequence of the fact that this is the stress measure conjugated to the deformation gradient F . The deformation gradient itself emerges because, within 180 an atomic-scale approach that naturally uses atomic coordinates that refer to a fixed reference state, it is most convenient to use a Lagrangian description within which the degrees of freedom of the fluctuating simulation box are simply the entries of the deformation gradient F (see

185
for example [START_REF] Miller | Molecular dynamics at constant cauchy stress[END_REF]). Finally, we mention that, even though it is not required for the implementation and integration of the kinetic equations, we could define an instantaneous internal first Piola-Kirchhoff stress whose statistical average adopts a virial form which, at equilibrium, is automatic-190 ally equal to the externally applied Piola-Kirchhoff stress (see Appendix A for details). We conclude this paragraph with a general remark on the proposed overdamped Langevin dynamics. Such dynamics, also known as Brownian dynamics, have already been 195 used for the modeling of a number of spatio-temporal processes, in which heavy particles (such as biomolecules) interact with a bath of light particles (e.g. solvent molecules). In the limit of vanishing ratio m/M , where m and M are the masses of the light and heavy particles, it 200 can be argued that the heavy particles follow a Markovian dynamics, which permits to exclude from the simulation the light particles, which are of no direct interest (see for example [START_REF] Erban | From molecular dynamics to brownian dynamics[END_REF]). The present situation is different, as the atomic species that constitute our materials cannot 205 be separated into subclasses with different masses. However, our aim is to show that, because of the randomness generated by the chaotic character of the phase space trajectories, the initial deterministic dynamics, which is of second order in time, can be replaced by a first order in 210 time stochastic dynamics. Of course, an exact formulation of the overdamped Langevin equations should proceed through an explicit coarse-graining procedure over the initial Newtonian dynamics, which requires the identification of a characteristic time over which time averages 215 can be performed while preserving the characteristic times associated to the dynamics processes we are interested in. This coarse-graining procedure would naturally lead to a coarse grained potential Φ cg ({x n i }) that will differ from the original potential Φ ({x n i }), as phonons will be adiabatically embedded into Φ cg ({x n i }), together with explicit expressions for the noise terms and for the viscosity coefficient ν, whose knowledge is needed to access to the time scale of our first order in time dynamics. Because of the intrinsically anharmonic character of the initial potential, we anticipate that the coarse-grained potential will be significantly softer than the original one, allowing to use large time step when the dynamics is numerically discretised. In this paper, we propose a simplification that consists in replacing the coarse-grained potential by the 230 original one. We therefore cannot associate a real time scale to our simulations. However, we show below that our first-order in time dynamics does reproduce the real nature of the dynamics in terms of reaction pathway and observed microstructures. This asserts that a first-order in time out-of-lattice dynamics can be safely used to simulate spatio-temporal processes even though the underlying atomic species cannot be separated into subclasses associated to different characteristic times.

Numerical implementation

The model described in section 2.1 has been implemented in a Fortran code of our own. We referred to the work of Goedecker [START_REF] Goedecker | Optimization and parallelization of a force field for silicon using openmp[END_REF] to optimize interatomic forces computation by parallel computing. To time integrate Eqs. (2.4), we used an explicit predictor-corrector method [START_REF] Burrage | Numerical methods for secondorder stochastic differential equations[END_REF][START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF].

Dimensionless equations are obtained by introducing an energy unit E 0 and a time unit t 0 = ν/E 0 , where ν is the viscosity term that enters into the dynamics of the scaled atomic coordinates. When needed, we display our results with reference to the dimensionless time τ = t/t 0 .

For simulation in the (N V T ) ensemble, the time integration is performed using an explicit predictor-corrector method with a dimensionless time step ∆τ = 10 -5 . For the (N P T ) ensemble, the simulation box viscosity γ is set equal to 0.0855 ν and the dimensionless time step is ∆τ = 10 -7 . More details on the numerical integration of Eqs. (2.4) are given in Appendix B.

Simulation setup and interatomic potential

To clarify the influence of different external conditions on martensite microstructure, we simulate the BCC→HCP 260 transition in both thermodynamic ensembles (N V T ) and (N P T ). We use periodic boundary conditions. In the (NPT) ensemble, we consider stress-free conditions by setting to zero the first Piola-Kirchhoff stress, allowing 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 equilibrate BCC titanium at 1400 K and then quench it at 700 K. We perform the quenching by an instantaneous rescaling of the temperature parameter that fixes the noise amplitudes in Eqs. 2.4.

275

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 (a 0 = 3.417 Å). The total number of atoms is 93312.

But first, to perform Molecular Dynamics or Langevin atomic simulations, a relevant interatomic potential is needed.

280

We have considered two empirical atomic potentials for titanium from the literature that could be relevant because they were in particular developed to study the BCC→HCP transition. The first potential, of the EAM type, proposed by Mendelev et al. [START_REF] Mendelev | Development of an interatomic potential for the simulation of defects, plasticity, and phase transformations in titanium[END_REF] and referred to as Ti-1 EAM, 285 was fitted to reproduce the HCP stacking fault energy, the BCC-HCP transformation temperature (T ∼ 1150K) and the melting temperature. The second potential, of the MEAM type, proposed by Henning et al. [START_REF] Hennig | Classical potential describes martensitic phase transformations between the α, β, and ω titanium phases[END_REF] describes the structure and energetics of α, β and ω phases 290 in Ti. Optimization of the parameters is performed using a database of density-functional calculations and yields an accurate potential as verified by comparison to experimental and density-functional data for phonons, surface and stacking fault energies, and energy barriers for ho-295 mogeneous martensitic transformations. In addition, the elastic constants, phonon frequencies, surface energies, and defect formation energies closely match density-functional results even when these were not included in the fitting procedure. The authors have also verified using Molecular 300 Dynamics that the equilibrium phase diagram is in close agreement with experimental measurements.

We tested these EAM and MEAM potentials by performing preliminary simulations with classical Molecular Dynamics using the simulator LAMMPS [START_REF] Plimpton | Fast parallel algorithms for Short-Range molecular dynamics[END_REF]. We ob-305 tained the following results: when the MEAM-type potential is used, we were able to observe a stable BCC phase transforming into HCP upon quenching, coherently with previous Molecular Dynamics simulations proving its ability to reproduce the whole temperature-pressure phase dia-310 gram of titanium [START_REF] Hennig | Classical potential describes martensitic phase transformations between the α, β, and ω titanium phases[END_REF]. On the other hand, when using the EAM potential we did not observe phase transition after cooling, although we were able to get a stable BCC structure at high temperatures. We increased the simulation duration up to 1 nanosecond and we tested different sim-315 ulation box sizes. However, the transition did not occur in any of the two thermodynamic ensembles. The possible reasons for that could be: i) the simpler functional form of the EAM potential compared to the MEAM. The lack of any angular dependency in the embedding term describ-320 ing the electron density makes the EAM much cheaper than the MEAM from a computational point of view but it impacts its ability to model materials with strong bond directionality metals with partially full-d shell, ii) the presence of a high energy barrier for the nucleation of the HCP 325 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 by following a previous work on many-body force field imple-mentation [START_REF] Goedecker | Optimization and parallelization of a force field for silicon using openmp[END_REF].

Variant and phase identification

To characterize the microstructure formation and evolution, we need to identify the different crystal structures and, especially, the different variants.

To identify the different HCP variants, we use a deformation gradient map representing the local lattice distortion. Indeed, when a material undergoes a martensitic transformation, several energetically equivalent variants differing in their relative crystallographic orientation appear [START_REF] Pitteri | Continuum models for phase transitions and twinning in crystals[END_REF]. 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 to represent the actual lattice sites' displacements [START_REF] Banerjee | Phase transformations: examples from titanium and zirconium alloys[END_REF]. Before explaining the procedure used to identify F , we first recall that a simple homogeneous deformation gradient cannot fully describe the martensitic transformation from BCC to HCP: 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 sublattice of the deformed lattice. In the present situation, this shuffling consists in translating every second basal plane of the hexagonal lattice obtained after 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-bodycentered modification into the hexagonal-close-packed modification of zirconium[END_REF] states that the following crystallographic planes and directions are parallel:

(110) bcc (0001) hcp ; [ 111] bcc [ 2110] hcp , (2.6)
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 .

(2.7)

The two mechanisms differ only in that the Burgers mechanism requires a rotation of ±5.

k = 1, . . . , 6 : D (k) 2 n = m∈Ωn ∆r nm (t * )-F (k)
n ∆, r nm (0) 2

(2.8) where Ω n is the neighborhood set that is associated with a given (110) plane. The local deformation gradient F n is 390 defined as the one that, among the six tensors F lim . To monitor the evolution of the phase fraction of each 400 phase without distinguishing variants, we use the Polyhedral Template Match analysis (PTM) [START_REF] Larsen | Robust structural identification via polyhedral template matching[END_REF] implemented in the software OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simula-930 tion data with ovito-the open visualization tool[END_REF]. This method classifies crystal structures according to the topology of the local atomic environment. It provides a flexible tool for structural identi-405 fication even in the 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 410 structural template has been set equal 0.14 Å.

Results

Simulations in the (NPT) ensemble with stress-free

boundary conditions After equilibrating the system in the BCC state at 415 T=1400 K, we drop the temperature down to T=700 K. After a short relaxation, the system transforms into a HCP structure. In Fig. 1a, we report the evolution of the BCC and HCP volume fractions during the transition. Almost no BCC 420 phase is left after the transformation has been completed. A non-negligible residual fraction of atoms (≈ 0.20) exhibits crystallographic structure different than HCP, which suggests that some defects are generated. In Fig. 1b, we report the evolution of variant fractions as 425 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 (1) and U (5) . Afterwards, the structure coarsens further and forms 430 a single variant U (5) domain (Fig. 2). Snapshots taken during the transition (Fig. 3) show that, when the microstructure is coarsening, HCP domains with same c axis orientation but different shuffling directions (referred to as a couple of "anti-variants" [START_REF] Gao | Diffuse scattering as an indicator for martensitic variant selection[END_REF]) come into contact and generate an antiphase boundary (see the inset Fig. 3b). This boundary is equivalent to a stacking fault when the 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 a HCP structure without any remaining BCC domain, similarly to what is observed in the previous (N P T ) simulation. Fig. 4a shows the evolution of the BCC and HCP volume fractions while Fig. 4b shows the evolution of the six variant volume fractions. The overall transformation proceeds through different stages. We first observe a nucleation stage, that extends up to point A in Fig. 4a and4b, during which local HCP fluctuations emerge. We underline that the length scale of these fluctuations is too small to allow the identification of different variants with the procedure outlined in section 2, thus leading to an apparent incompatibility between Fig. 4a and Fig. 4b. Indeed, the procedure used to identify HCP variants relies on neighborhoods 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, HCP nuclei 1 grow rapidly (from point A to point B) and stabilizes to 460 a quasi-stationary stage during which the six HCP variants reach finite volume fractions that are roughly constant (from point B to point C). Then, the system enters a stage during which the volume fraction of three HCP variants increases at the expense of the three others (from 465 C to D). After this growing stage (from point D on),the system stabilizes in a microstructure that consists of only three variants. The selected variants share the [111] BCC direction in the parent phase, i.e., a 11 20 HCP dense direction. 470 Fig. 5 shows the final 3-variant microstructure along a plane orthogonal to the common dense HCP direction, within a color map of the c axis orientation and a schematic indication of inter-variant misorientation. The three variants organize around several triple junctions by form-475 ing three interfaces which show the structure of the type I {10 11} twin boundaries (twinning plane K 1 = {10 11}, shear direction η 1 = 1 210 ). The strains associated with the 3 variants forming the 3-plate morphology respect the twinning equations [START_REF] Bhattacharya | Microstructure of Martensite: Why it Forms 735 and how it Gives Rise to the Shape-memory Effect[END_REF][START_REF] Bowles | The crystallography of martensite transformations i[END_REF]. As later discussed in the Dis-480 cussion part, this 3-plate geometry is not fully compatible with the formation of three {10 11} twins so further strain is required for its accomodation. In Fig. 6, we show four snapshots of the microstructural evolution during the transition. First, stable nuclei of all 485 the six variants appear (a), and all the different HCP domains develop (b). At this point, two stable triple junctions (indicated by arrows) are already formed and lead to the final 3-plate morphology after the subsequent microstructure coarsening.

We repeated the simulation in the (N V T ) ensemble several times and changed the random noise term. The time evolution of variant fractions (Fig. 7) show that in each case the system behaves similarly and, after the nucleation of all the possible variants, progressively selects a triplet. In all the simulations, the selected triplets share a 111 BCC direction i.e., a 11 20 HCP direction. In terms of microstructure, the selected triplet always organize in the 3-plate morphology. Only in one case, shown in Fig. 8, the selected variants form two laminates consisting of parallel twins along {10 11} HCP plane. At the crossing point between the laminates, an FCC domain appears. This FCC domain shares coherent interfaces with the neighboring HCP variants; these sharp interfaces consist in a one-layer thick transition from a {111} FCC plane to a HCP basal plane.

Discussion

The results of simulations performed in the (N V T ) and (N P T ) ensembles give highlights on the origin of different defects experimentally observed in martensite and confirms how deeply local mechanical constraints influence martensite morphology. In both thermodynamic ensembles, at the really beginning of the transition, all variants appear. This is expected, since we consider situations in which all variants are energetically equivalent and therefore equally probable to appear 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, only part of them end up to form the final microstructure. We observed two distinct microstructural evolutions, depending on the applied boundary conditions. For simulations in the (N P T ) ensemble with P = 0, the microstructure coarsens and forms a single variant domain with antiphase defects where anti-variants domains come into contact (Fig. 3). We recall that a couple of anti-variants are two variants with the same orientation but different shuffling directions [31]. On the other hand, for simulations in the (N V T ) ensemble, a triplet of variants with a common 11 20 HCP direction is systematically selected and form stable triple junctions that drive the overall microstructural evolution 530 (Fig. 6). In this case, the final microstructure is richer in interfaces and, consequently, has a higher energy with respect to the mono-variant domain obtained for simulations in the (N P T ) ensemble.

As discussed in section 3, in almost all the simulations 535 performed in the (N V T ) ensemble, the 3 selected variants cluster in a 3-plate geometry around the common dense direction. This morphology has been experimentally observed in pure titanium [START_REF] Farabi | Fiveparameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF], Ti-Nb shape memory alloys [START_REF] Chai | Self-accommodation in ti-nb shape memory alloys[END_REF], zirconium alloys [START_REF] Srivastava | Self accomodation morphology of martensite variants in zr 2.5 wt% nb alloy[END_REF], and agrees with predic-540 tions based on the phenomenological theory of martensite [START_REF] Bowles | The crystallography of martensite transformations i[END_REF]. The 3-plate cluster minimizes the overall mesoscopic shape strain after transition [START_REF] Farabi | Fiveparameter 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] because the three selected variants are self-accommodating [START_REF] Miyazaki | The shape memory mechanism associated with the martensitic transformation in ti-ni alloys-i. self-accommodation[END_REF][START_REF] Pitteri | Continuum models for phase transitions and twinning in crystals[END_REF]. Consequently, this morphology is strongly favored when 545 strain energy minimization is dominant in driving the microstructure evolution. The morphology is experimentally observed at different length scales (micrometer [START_REF] Farabi | Fiveparameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF] and sub-micrometer scale [START_REF] Srivastava | Self accomodation morphology of martensite variants in zr 2.5 wt% nb alloy[END_REF]) and numerically reproduced by us using a simulation domain of nanometer scale. This sug-550 gests that the dominant driving force is the elastic relaxation that largely dominates the interface energy. Furthermore, experiments show how these 3-variant triangles are formed in regions delimited by big martensite laths, which are supposed to originate from first nucleation events. The domain are then progressively filled by smaller and smaller triangles [START_REF] Srivastava | Self accomodation morphology of martensite variants in zr 2.5 wt% nb alloy[END_REF][START_REF] Farabi | Fiveparameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF]. These observations together with our simulation results suggest that: i) the 3-variant cluster formation is mainly driven by elastic relaxation, which manifests itself at different length scales. We here showed that, at the nanometer scale, the elastic relaxation is still dominant. Consequently, mechanisms driving the formation of larger microstructures can be easily studied by analysing the evolution of small BCC domains undergoing transition, ii) this 3-variant cluster formation is directly related to martensite nuclei forming under a situation of local confinement. The agreement of our simulations with experiments indicates that simple fixed-volume conditions are well adapted to reproduce this state of local constraints experimented by real systems.

We mentioned previously that we occasionally observe the appearance of domains with FCC crystallography. In this case, the 3 variants form two laminates with an FCC domain at the crossing point. Low-energy coherent interfaces are formed between the HCP {0001} basal planes and the {111} FCC planes. While in the simulation showing a 3-plate morphology the three variants have 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. 4). The possible presence of FCC phase after transition has little experimental evidence [START_REF] Nishiyama | Transmission 950 electron microscope study of the martensites in a titanium-3 wt% iron alloy[END_REF] and could be an artifact due to the potential. Nevertheless, it has been also reported in other Molecular Dynamics simulations of martensitic phase transition in zirconium [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 960 simulations of the martensitic phase transition process[END_REF]. Subsequent mechanical/thermal treatments, leading to further evolution of the microstructure, could then lead to its progressive extinction in favour of the lower energy HCP phase.

We go now on with the analysis of the defects generated during the transition under different stress conditions. Our 590 simulations give insight on the origin of different kind of interfaces that are found in martensite in pure titanium and corroborates hypothesis based on experimental observations.

The single-variant microstructure rising from simulations 595 in the (N P T ) ensemble (i.e., mimicking the absence of any local constraint in the surroundings) shows antiphase boundaries separating HCP domains with the same orientation of the c axis but different shuffling directions (a couple of anti-variants [START_REF] Gao | Diffuse scattering as an indicator for martensitic variant selection[END_REF]). Similarly, single-variant 600 martensite plates containing antiphase boundary networks has been experimentally observed in titanium alloys [START_REF] Banerjee | Substructure in titanium alloy martensite[END_REF] as well as in shape-memory alloys [START_REF] Otsuka | Physical metallurgy of ti-ni-based shape memory alloys[END_REF][START_REF] Matsuda | Crystallography and morphology of antiphase boundary-like structure induced by martensitic transformation in ti-pd shape memory alloy[END_REF]. These early experimental studies led to the hypothesis that these interfaces origin from the nucleation, growth and subsequent 605 impingement of martensite domains, i.e., they are a direct consequence of the randomness of shuffling displacements during the transition [START_REF] Banerjee | Substructure in titanium alloy martensite[END_REF][START_REF] Matsuda | Crystallography and morphology of antiphase boundary-like structure induced by martensitic transformation in ti-pd shape memory alloy[END_REF]. The microstructural evolution observed in our simulations (see Fig. 3), where one variant domain disappear leading two anti-variants do- Figure 10: Histograms of off-diagonal coefficients for a single variant U (4) in the (N V T ) ensemble at three different time steps (point A,B,C in the HCP fraction curve of Fig. 4). Theoretical values are shown in dotted line.

K 1 = {10 11}, shear direction η 1 = 1 210 ) [62, 1]
), as shown in inset of Fig. 5. This triple junction is not fully consistent with the inter-variant misorientation expected from the Burgers orientation relationship and, from a purely geometrical point of view, requires some further strain to be present, because the {10 11} pyramidal plane form an angle of 61.5 • with the basal HCP plane. The stability showed by these triple junctions during the microstructural evolution (see Fig. 6) together with the observation of this specific interface arrangement in all our simulations suggests that the energy cost linked to the additional strain required to form three {10 11} symmetric boundary is negligible when compared to the energy gain in forming three low energy coherent interfaces [START_REF] Wang | Atomic structures of symmetric tilt grain boundaries in hexagonal close packed (hcp) crystals[END_REF]. This corroborates experimental observations in pure titanium 630 [START_REF] Wang | Effect of self-accommodation on α/α boundary populations in pure titanium[END_REF][START_REF] Farabi | Fiveparameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF], titanium alloys [START_REF] Beladi | Variant selection and intervariant crystallographic planes distribution in martensite in a ti-6al-4v alloy[END_REF] and zirconium alloys [START_REF] Srivastava | Self accomodation morphology of martensite variants in zr 2.5 wt% nb alloy[END_REF]. In particular, a recent study on grain boundary plane distribution in pure titanium subjected to temperature driven martensitic transformation, which completes previous experimental observations on grain boundary axis angle displicitly observed the three-variant cluster in triple junction morphology [START_REF] Farabi | Fiveparameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF]. We conclude our discussion by analyzing the numerically computed stretch tensors used to identify variants. Fig. 9 shows the histograms of the diagonal and off-diagonal components of U n in the simulations in the (N P T ) and (N V T ) ensemble, with dotted lines showing the theoretical values predicted by the Mao/Burgers mechanisms. We note that for the simulation in the (N P T ) 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. 10 compares the histograms of the off-diagonal coefficients for the simulation in the (N V T ) ensemble at three different time steps t 1 , t 2 , 655 t 3 , corresponding to the end of the nucleation stage, the quasi-stationary regime and the final stationary state (see Fig. 4). Differently from the histograms of Fig. 9, we report the value for only one of the three variants forming the final microstructure, the U (4) . The comparison between 660 histograms highlights how the deviation from the theoretical strain values (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 particularly 665 evident for the off-diagonal coefficient U 13 . Although further investigations are needed, we can conclude that these deviations are related to the additional textural evolution arriving in fixed-volume conditions when the first nucleation stage has ended and specific variants domains start 670 growing around the stable triple junctions to form the final 3-plate morphology.

Conclusions

To summarize, this work is the first numerical study at the atomic scale of the microstructural evolution of pure ti-675 tanium undergoing temperature-driven BCC→HCP transition under different stress conditions. For this purpose, we performed a set of extended atomistic simulations with suitable empirical interatomic potential, and we analyzed the microstructural evolution during the transition as well 680 as the final microstructures. Our main results concern the analysis of the atomistic mechanisms inducing the formation of different defects experimentally observed in martensite in pure titanium and the assessment of the influence of macroscopic constraints 685 on these defects and on the final martensite morphology. When no constraints are present, i.e., when the crystal is allowed to change its shape, a simple mono-variant domain decorated by wavy antiphase boundaries forms. Our simulations confirm previous experimental hypothesis that 690 trace back the origin of these interfaces in the growth and subsequent impingement of variant domains with same orientation but different shuffling directions [START_REF] Banerjee | Substructure in titanium alloy martensite[END_REF][START_REF] Matsuda | Crystallography and morphology of antiphase boundary-like structure induced by martensitic transformation in ti-pd shape memory alloy[END_REF]. In contrast, a poly-variant microstructure develops when local constraints prevent a free deformation of the environ-695 ment surrounding the growing martensite nuclei. This microstructure shows a specific 3-variant morphology which has been extensively documented in experiments on titanium, zirconium and their alloys and which allow the minimization of the strain energy [START_REF] Farabi | Fiveparameter 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] Beladi | Variant selection and intervariant crystallographic planes distribution in martensite in a ti-6al-4v alloy[END_REF]. This micro-700 structure develops around stable triple junctions between variant that are formed at the beginning of the transition after the first nucleation stage. The characterization of interfaces in this microstructure confirms a strong preference for the formation of boundaries along the {10 11} 705 HCP pyramidal plane as experimentally documented [START_REF] Farabi | Fiveparameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF]. Also, we observed the possible appearance of an FCC phase after transition, although further studies are needed to check the absence of any artifact due to the use of empirical inter-atomic potentials [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 960 simulations of the martensitic phase transition process[END_REF].

Finally, we stress that this is the first time that the overdamped Langevin dynamics, which has been mostly applied in field of soft matter and bio-molecular simulation, is successfully applied to simulate a fully 3D displacive solid phase transition.This is an important step towards the use of a first-order in time dynamics. The full application of this modelling tool would require its proper derivation through coarse-graining, which will adiabatically eliminate phonons through their incorporation within a coarse-grained potential. This potential will be much softer than the initial one, allowing to use much larger time steps than those required when using the original one. A natural extension of our work will be to investigate how final microstructures here obtained may 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 730 BCC-HCP transition, formulated using finite strains [START_REF] Vattré | Polymorphism of iron at high pressure: A 3D phase-field model for displacive transitions with finite elastoplastic deformations[END_REF][START_REF] Denoual | A phase field approach with a reaction 970 pathways-based 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].

Piola-Kirchhoff stress. We also discuss the link between this instantaneous first Piola-Kirchhoff stress and an instantaneous Cauchy stress, even though the prior computation of a Cauchy stress is not required for our (N P T ) calculations.

As explained in the text, within the (N P T ) ensemble, the internal degrees of freedom associated with the fluctuating box are the entries of the deformation gradient F whose dynamics are given by Eq. (2.4) which, for the sake of completeness, we recall here

1000 dF ij dt = -γ -1 ∂ H ∂F ij + 2k B T γ -1 ξ ij (t) i, j = 1, ..., 3 , (5.1)
where ξ ij (t) is a white Gaussian noise such that ξ ij (t) = 0, ξ ij (t)ξ lm (t ) = δ ij δ lm δ(t -t ) and η n i (t)ξ lm (t ) = 0. δ ij and δ lm are Kronecker symbols and δ(t -t ) stands 1005 for the Dirac-delta distribution. As explained in the main text, the driving forces must be computed with the extended Hamiltonian H given by Eq. (2.5)

H = Φ({F ij L 0 j xn j }) + V 0 P ij F ij -N k B T ln (V 0 det F ) , (5.
2) which must be considered as a function of the scaled co-1010 ordinates {x n i } and of the deformation gradient F . The scaled coordinates are related to the initial atomic coordinates by

xn i = H -1 ij x n j , i = 1, 2, 3, (5.3)
where the matrix H is defined by H = F L 0 , where L 0 is a diagonal matrix containing the length of the orthogonal vectors that define the initial simulation box. The driving force that enters Eq. (5.1) is given by

∂ H ∂F ij = ∂Φ ∂F ij + V 0 P ij -N k B T F -T ij = n,l ∂Φ ∂x n l ∂x n l F ij + V 0 P ij -N k B T F -T ij .
(5.4) Using Eq. ( 5.3), it is trivial to show that

1015 ∂x n l ∂F ij = δ li x n k F -T kj , (5.5)
which inserted in Eq. (5.4), leads to

∂ H ∂F ij = n ∂Φ ∂x n i x n k F -T kj + V 0 P ij -N k B T F -T ij . (5.6)
Introducing the virial tensor V, defined by

V ij = - n ∂Φ ∂x n i x n j , (5.7)
1020 the driving forces become

∂ H ∂F ij = -{N k B T δ ik + V ik }F -T kj + V 0 P ij . (5.8)
We now define the instantaneous internal first Piola-Kirchhoff stress as

P inst ij = 1 V 0 (N k B T 1 + V)F -T , (5.9
)

1025
where 1 is the identity matrix. Eq. (5.8) becomes

∂ H ∂F ij = V 0 (P ij -P inst ij ), (5.10)
The kinetic equation (5.1) then reads

dF ij dt = -γ -1 V 0 (P ij -P inst ij ) + 2k B T γ -1 ξ ij (t), (5.11)
At equilibrium, the statistical average of the l.h.s. of this 1030 equation is equal to zero (by "statistical average", we mean a time average over a sufficient long time window). As, by definition, the statistical average of the noise term is also equal to zero, we get

P inst ij = P ij , (5.12)
1035 where X stands for the statistical average of X. This equation invites us to define an internal first Piola-Kirchhoff stress P int as the statistical average of the instantaneous stress P inst

P int = P inst = 1 V 0 (N k B T 1 + V)F -T . (5.13) 1040
Once the equilibrium is reached, this internal first Piola-Kirchhoff stress, which is defined unambiguously only at equilibrium because it relates on a time average, equilibrates exactly the imposed Piola-Kirchhoff stress:

P int = P . (5.14) 1045
We note that the internal stress defined in Eq. (5.13) adopts a virial form, as it is common to any internal stress defined at the atomistic scale and, therefore, linked to interatomic forces. We also note that the numerical computation of this Piola-Kirchhoff stress is straightforward 1050 and does not require the prior computation of any other stress, such as a Cauchy stress, even tough an internal Cauchy stress could also be independently defined and related, through usual relations, to our internal first Piola-Kirchhoff stress (see below). Finally, we also note that this 1055 internal first Piola-Kirchhoff stress is not a local quantity defined for each point within the simulation box but rather a global quantity that is associated with the whole system.

Therefore, its computation does not require any recipe to define and compute numerically a local stress, such as the 1060 ones proposed by Hardy [START_REF] Hardy | Formulas for determining local properties in molecular-dynamics simulations: Shock waves[END_REF].

Up to this point, the first Piola-Kirchhoff stress P int has been directly defined in the (N P T ) ensemble within which the deformation gradient F is a fluctuating quantity.

It has naturally emerged within the kinetic equations associated with the box shape and is such that, at equilibrium, it equilibrates the externally applied first Piola-Kirchhoff stress. We show now that, as expected, it may also be associated, by conjugacy, with the deformation gradient F .

To show this, we need to introduce the canonical (N F T ) statistical ensemble, in which F is fixed, and to compute the associated free energy F whose derivative with respect to F , when properly averaged, will lead to P int . Within the (N F T ) ensemble, the partition function Z is 1075 given by

Z = 1 Λ 3N N ! Π n,i dx n i e -βΦ({x n i }) , (5.15)
where Φ({x n i }) is the interatomic potential and β = 1/(k B T ). The term Λ 3N , where Λ is the de Broglie wavelength, is reminiscent of the quantum and therefore discrete nature of 1080 the problem. It appears in the classical limit of quantum mechanics and ensures a proper normalization of entropy and free energy. Our aim here is to compute the derivation of the free energy F with respect to the deformation gradient F . Therefore, we must introduce the scaled co-1085 ordinates {x n i } which are related to the initial coordinates {x n i } through the deformation gradient F :

x n i = (F L 0 ) ij xn j .
The partition function then becomes

Z = 1 Λ 3N N ! (det F L 0 ) N Π n,i 1
0
dx n i e -βΦ({F L 0 xn }) .

1090

Taking the derivative of the free energy F = -k B T log Z with respect to the deformation gradient F , we get:

- ∂F ∂F ij = k B T    N F -T ij + Π n,i 1
0 dx n i -β ∂Φ ∂Fij e -βΦ Π n,i 1 0 dx n i e -βΦ    = k B T N F -T ij + -β ∂Φ ∂F ij N F T = k B T N F -T ij + n -β ∂Φ ∂x n i x n k F -T kj N F T .
(5.16)

where X N F T is the statistical average of X within the ensemble (N F T ). We now define the canonical first Piola-Kirchhoff stress as the stress conjugated to the deforma-tion gradient F4

P cano = - 1 V 0 ∂F ∂F ij .
(5.17)

1095

Using Eq. (5.16), we get:

P cano = 1 V 0 (N k B T 1 + V) F -T N F T , (5.18)
where, as it does not fluctuate within the (N F T ) ensemble, the constant term N k B T F -T has been included within the statistical average and where the virial tensor 1100

V has been defined in Eq. (5.7).

Comparison of Eq. (5.13) which gives P int , the internal first Piola-Kirchhoff stress within the (N P T ) ensemble, and Eq. (5.18) which gives P cano , the canonical first Piola-Kirchhoff stress defined within the (N F T ) ensemble, shows 1105 that P int and P cano differ only through different statistical averages associated with their respective statistical ensembles. Obviously, for any observable X, a statistical average within the (N P T ) ensemble may be split into a statistical average at a fixed F followed by an average over 1110 the fluctuations of F , which, with obvious notations, leads to X = X F N F T .

(5. [START_REF] Heiming | Phonon dispersion of the bcc phase of 795 group-iv metals. ii. bcc zirconium, a model case of dynamical precursors of martensitic transitions[END_REF] where, as in Eq. 5.12 and 5.13, X refers to the statistical average of X within the (N P T ) ensemble. Thus, 1115 P int , defined in the (N P T ) ensemble, is related to P cano , defined within the (N F T ) ensemble, by

P int = P cano F = - 1 V 0 ∂F ∂F F N F T .
(5.20)

In conclusion, the internal Piola-Kirchhoff stress defined in Eq. (5.13) is, as expected, related by conjugacy to the 1120 deformation gradient F .

The introduction of the first Piola-Kirchhoff stress in the model used here is simply a consequence of the fact that it is the stress measure related by conjugacy to the 1125 deformation gradient F , which, within a lagrangian setup, is the degree of freedom associated to the fluctuating box. Of course, we could also introduce other stress measures, even though this is not needed to integrate our Langevin dynamics. As example, we could define an instantaneous 1130 Cauchy stress in such a way that it is related to the instantaneous first Piola-Kirchhoff stress through the usual relation:

σ inst 1 = 1 det F P inst F T , (5.21)
which, with Eq. (5.9), leads to:

1135 σ inst 1 = 1 V (N k B T 1 + V), (5.22)
whose statistical average in the (N P T ) ensemble leads to the definition of an internal Cauchy stress:

σ int 1 = 1 V (N k B T 1 + V) .
(5.23)

The question now arises as to whether this internal Cauchy stress is equal, at equilibrium, to an external Cauchy stress. Naturally, we would like this external Cauchy stress σ to be related to the applied first Piola-Kirchhoff stress by the usual relation:

  n . Finally, polar decomposition leads to the local stretch tensor U n and to a local stretch deformation map. The non-affine displacement D (k) 2 n quantifies the degree at which an affine transformation can describe 395 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

Figure 1 :

 1 Figure 1: Evolution of BCC and HCP volume fractions (a) and variant fractions (b) in the (N P T ) ensemble. The symbols A,B,C,D indicates the times selected for displaying the microstructures in Fig. 2.

Figure 2 :

 2 Figure 2: Microstructure evolution in the (N P T ) ensemble. Only atoms within a thin slice of the simulation box and 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.The microstructures a,b,c,d correspond to the reduced times 0.01, 0.04, 0.07, 0.2, respectively. These times are indicated in Fig.1

Figure 3 :

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

Figure 4 :

 4 Figure 4: Evolution of BCC and HCP volume fractions (a) and variant fractions (b) in the (N V T ) ensemble.

Figure 5 :

 5 Figure 5: Final microstructure obtained at 700 K in the (N V T ) 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 6 :

 6 Figure 6: Evolution of the microstructure in the (N V T ) 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.4a), 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.4a) in a 3-plate morphology (final stable state). Crystallographic directions refer to the parent BCC phase.

610

  mains to come in contact, confirms these hypothesis. For simulations performed in the (N V T ) ensemble showing the 3-variants plate morphology, we verified that the grain boundaries separating the different variant domains are all {10 11} type I twin boundaries (twinning plane 615

Figure 7 :

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

Figure 8 :

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

Figure 9 :

 9 Figure 9: Histograms of Un diagonal and off-diagonal coefficients for a)-b) simulations in the (N P T ) and c)-d) (N V T ) ensemble. The dotted line shows theoretical values.

  (5.13) and(5.14), which say that, at equilibrium in the (N P T ) ensemble, the average of the instantaneous first Piola-Kirchhoff stress is equal to the applied first Piola-Kirchhoff stress, and Eq. (5.21) for the definition of the instantaneous Cauchy stress, we get:

Figure 11 :

 11 Figure 11: a) Six possible orientations of the (110) bcc planes which may transform into (0001) hcp planes during the BCC-HCP transformation. Note that the central atom of the BCC cubic cell is not shown for the sake of clarity. b) Example of a neighbor set Ωn (colored in black, containing 14 atoms) for a given atom n (colored in red) for one of the six orientations.

Table 1 :

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

  . For each of the six (110) BCC planes transforming in the final (0001) HCP plane, we define a specific set of neighboring sites, see Appendix C. Then, for each 385 atom n, we identify six local deformation gradients F . . . , 6, that minimize the following local descriptors:

	(k) n ,
	k = 1,
	26 • around the [0001] HCP
	axis in order to obtain the proposed direction correspond-
	ence [58]. Consequently, whereas the Mao relationship
	generates only 6 HCP orientation variants, the Burgers
	mechanism generates 12 HCP lattices. However, as they

differ only by rotations, the two mechanisms are associated with exactly the same six stretch tensors U (k) , k = 1, . . . , 6. These tensors are listed in Tab. 1. We define a procedure meant to identify these local stretch tensors

U (k) n

Note that the sign convention used here is implied by the enthlapy definition (Φ + V 0 P F ) used to introduce the extended Hamiltonian in Eqs. (2.5) and (5.2).

tribution [START_REF] Wang | Effect of self-accommodation on α/α boundary populations in pure titanium[END_REF], highlights a strong anisotropy in the grain boundary plane distribution with most of the grain boundaries terminating on {10 11} pyramidal planes [START_REF] Farabi | Fiveparameter intervariant boundary characterization of martensite in commercially pure titanium[END_REF]. The authors report that these boundaries are associated with symmetric tilt 60 • [ 11 20] inter-variant boundaries and ex-Appendix A

985

We show below that within the statistical ensemble (N P T ), there is a natural way to define and compute an internal instantaneous first Piola-Kirchhoff tensor that adopts a virial form and whose statistical average at equilibrium is automatically equal to the externally applied

Because of the coupling between the fluctuations of det F , σ inst and F , the r.h.s of this equation cannot be further simplified. Therefore, strictly speaking, the statistical average of the instantaneous Cauchy stress as defined in 1155

Eqs. (5.21)- (5.22) is not equal to the applied Cauchy stress defined in Eq. (5.24)

As an alternative, we could define an instantaneous Cauchy stress in the following way:

(5.27)

Using Eqs. (5.13), (5.14) and (5.24) we immediately see that the statistical average of this instantaneous stress is now equal to the applied Cauchy stress:

(5.28)

1165

However, we note that, strictly speaking, the definition of the instantaneous Cauchy stress given in Eq. (5.27) is not entirely satisfactorily because its numerical application requires the prior knowledge of the statistical averages of det F and F . As a final comment, we note that 1170 in the thermodynamic limit of an infinite system, fluctuations may be neglected (provided the system is not going through a phase transition): definitions (5.21) and (5.27) become equivalent and Eq. (5.26) becomes and equality.

Appendix B: integration scheme 1175

Our Langevin dynamics is defined by stochastic equations with white noise (see Eqs. (2.4)) that display the following generic form:

where a({X i }) is a drift term, B a noise amplitude and the 1180 differential dW (t) denotes an infinitesimal increment of the Wiener process W i (t). The integral form of Eq. (5.29) is:

where the first term in the r.h.s. is a Riemann integral and the second term is a stochastic integral. To numerically 1185 evaluate Eq. (5.30) we used an explicit predictor-corrector method which results in the following numerical scheme [START_REF] Burrage | Numerical methods for secondorder stochastic differential equations[END_REF][START_REF] Baruffi | Overdamped langevin dynamics simulations of grain boundary motion[END_REF]:

where the finite increment

1190 can be calculated as ∆W i (t) = √ ∆tξ(t) with ∆t ∈ R and ξ(t) taken from a normal distribution with unit variance. Xi is evaluated using the explicit Euler method:

Appendix C: numerical procedure for variant iden-

tification

The numerical procedure used to identify the different variants 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 1200 following procedure as described in [START_REF] Wu | The vari-980 ant 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 mainframe axis. For every atom n, we consider six possible sets of 14 neighbors by taking six different configurations defined based on the six cubic cells which 1205 can deform into the orthorhombic one, as schematically illustrated in Fig. 11. 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 1210 of the deformed lattice are needed. These displacements consist in an alternate shuffling of {110} planes along the 110 directions. The overall deformation of the lattice does not describe the shuffling so, to numerically evaluate the local strain, half of the atoms of the original BCC lat-1215 tice 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 (k) n , 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 n . Based on this assignment, we label the atom n as belonging to the corresponding variant.