Rui Li
email: rui.li@orange.com

Bertrand Decocq
email: bertrand.decocq@orange.com

Anne Barros
email: anne.barros@centralesupelec.fr

Yiping Fang
email: yiping.fang@centralesupelec.fr

Zhiguo Zeng
email: zhiguo.zeng@centralesupelec.fr

Petri Net-Based Model for 5G and Beyond Networks Resilience Evaluation

Keywords: Petri Net, 5G networks, B5G, resilience, availability, modeling, simulation

The promise of telecommunication networks to deliver more demanding and complex applications requires them to become more flexible and efficient. To achieve better performance, telecommunication networks adopt technologies such as NFV (Network Function Virtualization). However, this evolution also brings more potential risks to the telecommunication network. Reliability and resilience are becoming critical for service delivery in the networks. To answer to service requirements of high level availability and reliability, a model with a global view of infrastructure, virtual network elements, and network layer structure is required. Toward this end, this paper presents a Petri Net method to model 5G and beyond telecommunication networks. We introduce an extended Petri Net to model physical infrastructure, virtual infrastructure, network services, their behaviors, and dependencies. We present a simulation result on network availability estimation. This result shows the potential of the Petri Net-based model to be applied to a complex telecommunication system resilience assessment.

I. INTRODUCTION

Telecommunication networks are becoming indispensable for modern production and living. Facing the diverse and high requirements from a broad vertical industry, 5G and beyond networks are expected to be both efficient and reliable for service delivery. Keeping such systems at good performance during their whole life cycle is essential for service providers and operators. Although 5G has been under development for the last ten years, its resilience has not been studied enough.

5G is more service-oriented than 4G by offering a transition from a "horizontal" service delivery model toward a "vertical" one [START_REF] Banchs | A 5G mobile network architecture to support vertical industries[END_REF]. To better meet the different requirements and high demand, Network Function Virtualization (NFV) should be introduced for both RAN (Radio Access Network) and CN (Core Network) to better adjust the network configuration to the requirements. Thus, the resilience of the telecommunication networks is no longer an issue only for infrastructure but also for virtual elements and service delivery [START_REF] Li | Complexity in 5G Network Applications and use cases[END_REF].

This study builds a Petri Net-based model to describe 5G and beyond networks. This model could be applied to communication service availability [START_REF] Sayrac | 5G Common Terminology[END_REF] (the ability to allow correct operation of the application) analysis, communication service reliability (the measure of continuous correct service delivery) estimation, and the resilience [START_REF][END_REF] (the ability to provide and maintain an acceptable level of service in the face of various faults and challenges) evaluation of the system.

While still at an early stage, this paper introduces a case study regarding network virtualization characteristics for testing Network Function Virtualization (NFV) self-healing in the dysfunctional mode and analyzing the network availability.

The paper is structured as follows. First, we present related work concerning telecommunication resilience analysis and modeling in section II. Then in section III, we focus on the telecommunication network model and, in particular explaining how an extended Petri Net models the virtualization characteristics. A case study on self-healing and the results are given in section IV. Section V concludes the work with some remarks and outlines the future works.

II. RELATED WORK

Recently, some research regarding the 5G and beyond telecommunication network resilience has been carried out. These studies mainly focus on network resilience optimization, only a few on resilience assessment. In optimization, while solving linear programming problems is still the mainstream of the research as found in [START_REF] Hößler | Applying reliability theory for future wireless communication networks[END_REF]- [START_REF] Harutyunyan | Latency-Aware Service Function Chain Placement in 5G Mobile Networks[END_REF], other methods such as Shortest Paths [START_REF] Qu | A reliability-aware network service chain provisioning with delay guarantees in NFVenabled enterprise datacenter networks[END_REF], Divide and Conquer [START_REF] Carlinet | Minimum-cost virtual network function resilience[END_REF] are also applied. The complexity of solving such a problem grows with the number of constraints. However, numerous new constraints on resource allocation and network management will be required if the virtualization layer is considered. In assessment, methods such as Reliability block diagram [START_REF] Mauro | Comparative performability assessment of SFCs: the case of containerized IP multimedia subsystem[END_REF], Markov chain [START_REF] Mosahebfard | Modelling the admission ratio in NFV-based converged optical-wireless 5G networks[END_REF], [START_REF] Farooq | Continuous time Markov chain based reliability analysis for future cellular networks[END_REF] are addressed. These methods fail to model a complex telecommunication network by considering infrastructure, virtualization, network layers, and their dependencies. A model that can capture all network elements and their relationship for resilience assessment and optimization is still missing.

As a critical technique for the 5G and beyond, Network Virtualization has been particularly studied in some works. An availability model and analysis of a virtualized system based on Virtual Machines(VMs) are introduced in [START_REF] Kim | Availability modeling and analysis of a virtualized system[END_REF]. The authors in [START_REF] Khan | A performance modelling approach for SLA-aware resource recommendation in cloud native network functions[END_REF] present a performance modeling approach that goes into the microservice level to estimate the effect of resource configuration on the Quality of Service (QoS). Although VMs are already widely used in virtualization, containerization, as a novel and lightweight virtualization method, is believed to be a promising solution for 5G and beyond. However, rarely is container-based virtualization modeled. Because it makes the telecommunication network modeling even more complex [START_REF] Li | Complexity in 5G Network Applications and use cases[END_REF].

Some recent works draw attention to Petri Net-based model to study telecommunication networks thanks to its convenience in modeling discrete event systems. In [START_REF] Rui | Petri Net-based reliability assessment and migration optimization strategy of SFC[END_REF], a Petri Net model is proposed for Service Function Chain (SFC) reliability assessment. However, the model does not take into consideration the risks of infrastructure failures. A Queuing Petri Net model is applied in [START_REF] Schneider | Specifying and analyzing virtual network services using queuing Petri Nets[END_REF] to evaluate the QoS of a video streaming service. The dysfunctional state of the system is not yet considered in their study. In [START_REF] Mauro | Availability modeling and evaluation of a network service deployed via NFV[END_REF], the authors apply a Petri Net model to describe the probabilistic behaviors of network service. In [START_REF] Tola | On the resilience of the NFV-MANO: an availability model of a cloud-native architecture[END_REF], an extension of Petri Net is applied only to model the NFV MANO framework to analyze its availability. Finally, the authors in [START_REF] Medel | Modelling performance and resource management in kubernetes[END_REF], introduce a Petri Net-based performance model for containerized applications deployed by Kubernetes. However, these works are limited and do not propose a comprehensive perspective by considering NFV characteristics, infrastructure behaviors, and QoS in functional and dysfunctional mode.

Inspired by the related work, this paper proposes an extended Petri Net approach to fill the gaps in related work. This extended Petri Net is given to model SFC, containerizationbased NFV elements, and infrastructure layer of the network. It is also capable to estimate the performance and the resilience.

III. PETRI NET-BASED NETWORK MODEL

A. Timed Stochastic Colored Petri Net 1) Petri Net mathematical representation: Petri Net is also known as Place / Transition net. It is a widely used technique tracking systems' states, dynamics, and constraints. As well defined in [START_REF] Murata | Petri Nets: Properties, analysis and applications[END_REF], the Petri Net is a 5-tuple N = P, T, F, W, M 0 , where P is a finite set of places often representing the different states of a system. Places are graphically presented in circles. T is a finite set of transitions representing the state-changing process. Transitions are graphically presented in rectangles or squares. F is a finite set of arcs with F ⊆ (P × T ) ∪ (T × P ). W is a multi-set of arcs (P × T ) ∪ (T × P ) → N assigning the weight to inputs and outputs of a transition. M is the marking of the Petri Net graph and M 0 = P → {m 1 , m 2 , ..., m |P | }, therefore, assigning the initial marking of the graph. Tokens of the graph describe the dynamic and concurrent activities of systems. The marking in Petri Net records the token number of each place.

2) Extensions of Petri Net: The classical Petri Net is not directly applicable to telecommunication systems. Some statechanging processes in such a system could be stochastic, timedependent, and require additional information. Some extensions of Petri Net can help modeling a complex network better.

One of the most important extensions of Petri Net is Stochastic Petri Net [START_REF] Marsan | A class of generalized stochastic Petri Nets for the performance evaluation of multiprocessor systems[END_REF]. It includes a new set R = {r 1 , r 2 , ..., r |T | }, representing the firing rate of each transition. This extension could be applied to describe a failure process in the telecommunication network.

In order to describe a time-dependent process, for instance, the packet transmission, Timed Petri Nets [START_REF] Popova-Zeugmann | Time petri nets[END_REF] are introduced. A new set D : T → Q + 0 associates each transition with a specific non-negative number to represent the time factor.

Colored Petri Net attaches a value to a token. It indeed distinguishes different kinds of tokens that a place holds. This extension adds the following items [START_REF] Jensen | Basic Concepts[END_REF], [START_REF] Liu | Modeling workflow processes with colored Petri Nets[END_REF]:

1) Σ is a finite set of non-empty types, called color sets.

2) C is a function P → Σ defining the type of tokens allowed in a place. 3) G: T → B associate the transition with a precondition g (Boolean expression). The transition will be fired only when g returns true value. 4) E is an arc expression function defined from F into expressions such that ∀a ∈ F : C(E(a)) = C(p). 5) I is an initialization function mapping place p ∈ P with an expression such that I(p) is associated to C(p). Combining the extensions as aforementioned, we use a Timed Stochastic Colored Petri Net (TSCPN) to describe the 5G system. Such a TSCPN is a multi-tuple: T SCP N = Σ, P, T, F, W, m 0 , C, G, E, I, R, D TSCPN is then applied to describe the different parts of the telecommunication system.

B. Composition of the network system

We divide the 5G and beyond networks into three layers. The first layer is the service layer. In this layer, the network service is delivered by steering packets between a set of functions called the service function chain. We consider a network service where all its network functions engaged are virtualized. The service delivery is presented by a series of Virtual Network Functions(VNFs) connected via virtual links.

The second layer is the NFV elements layer. A virtual link in this layer is based on a physical transport network, and a VNF is a virtual functional building block hosted on a physical server. Unlike physical elements, virtual elements may have an unfixed size and an unfixed number of replicas.

The third layer is the infrastructure layer. Physical machines and physical links belong to this layer.

C. Basic hypothesis 1) Service function chain:

We consider an End-to-End service with an ordered SFC. In this system, user equipment sends service request packets to the SFC in the network. We assume that every packet conveys a same size of data, and its SFC always follows the same order of VNFs.

An SFC is a series of VNFs connected by links. The transport network is considered as a perfectly reliable system. We only consider a fixed time delay spent on the transmission link between a user equipment and VNF, and between different VNFs.

2) VNF and Virtualization: A VNF is, in fact, an application that consists of several microservices. An example of 3 VNFs and their microservices are shown in Fig. 1. Each microservice is considered as a sub-function of VNF. We assume that if a request packet needs multiple microservices, they should be pursued in a given order. A packet cannot consult two different microservices at the same time. There are many ways of virtualization. In this paper, we choose to model the deployment of these microservices in containers. Kubernetes is used as the system for automating deployment and managing containerized applications.

Pods are the smallest deployable units in Kubernetes. A pod is one or a cluster of containers with shared storage and network resources. We assume that only one container is deployed on a pod. For each container, it corresponds to a microservice the VNF supplier predefines. Pods are running on Kubernetes nodes. All these nodes are physical machines.

3) Infrastructure and resources: The infrastructures used to deliver an end-to-end function are physical links and physical servers. Each physical machine has a certain amount of CPU, storage, and network resources. Pods can only be hosted on the server with enough resources.

4) Orchestration and management: Kubernetes is an enabler for the orchestration and management of containerized applications. To evaluate the system resilience under failure, we consider self-healing operation. Kubernetes can regularly detect the healthiness of the pods or nodes. In case of failure, they will be terminated, and new ones will be created. Other operations such as auto-scaling, in which Kubernetes detects particular indicators and change the deployment manners accordingly, will be studied in future work. 2) Virtual Network Functions: We expand the VNF i process transition t VNF i process as depicted in Fig. 3. This transition takes a token from the place p Start of VNF i , and after processing, it returns the token to the place p End of VNF i . Inside a VNF i, there are n microservices embedded in containers. The token color set is Σ = {'Packet', 'Packet list'}. The latter refers to a list of packets to be treated. When a packet gets into this i-th VNF, an immediate transition t service selection routes them to the microservice k it looks for with the help of guard function G. Then this packet is inserted to the 'Packet list' token at the place p MS queue k . Microservice k transition t MS k-transition is enabled as long as there is at least one microservice pod with enough capacity left to treat this packet. After the treatment, the packet arrives at p End of MS in VNF i . Two assertion transitions associated with guard functions will check if other microservices should be consulted before leaving current VNF. 3) Microservices: Microservices level Petri Net explains the transition t MS k-transition in detail. The microservices applications are in the form of containers and they are embedded into pods. A new token color type, 'Pod', is added. These tokens are stored at p available pods and p failed pods .

The microservice is modeled with two transitions, as shown in Fig. 4. The transition t MS bounding couples the first packet of the packet list token with an available pod. A Boolean guard function B associated with this transition checks if the pod is eligible to provide service for the packet. To complete the task, the packet will borrow a certain computation resources from the pod. These resources are seen as homogeneous by assumption. A function D assigns the processing time of the microservice to this transition. An output packet token will be sent to p bounded packets after this duration. The transition t MS process verifies the completion of the service. If a pod fails during the packet treatment, the packet needs to redo the same microservice. Otherwise, the borrowed resource will be released, and the packet will try the next microservice. 4) Failure and Self-healing: Numerous failures could happen in a telecommunication system. In this study, we mainly consider pod failure, one of the most common failures that affect the a vitrulized system's performance.

Due to space, only pod software failure is explained, as shown in Fig. 5. A stochastic transition t Stochastic pod failing process connects the state change of a pod token. We assume that all pods in our system are identical, and thus, they have the same mean time to failure (MTTF). We also assume that a pod failure is in accordance with the exponential distribution X ∼ Exp(λ), and with a constant rate of λ = MTTF -1 .

Kubernetes launches a liveness probe once in a while to detect the healthiness of pods. This time interval is called periodsecond of the probe. If a pod is unhealthy, Kubernetes starts the self-healing by terminating the pod and creating new one. The transition t Pods termination consumes the failed token after a graceful termination time. Transition t Pods creation will create a pod containing the same microservice on an available node with enough resources. We introduce place p Available node containing a new token color, 'Node', to represent the nodes.

IV. NETWORKS AVAILABILITY ESTIMATION

The first step of our work is to apply the model to estimate the system resilience by looking at the virtualization and infrastructure layer without mapping them to telecommunication services. By doing so, the availability of the network to provide services to the packets is estimated. Two major failures, the physical failure on nodes and the software failure on pods, are identified as the main risks to the system. When a failure occurs, the Kubernetes Master will do self-healing to ensure availability. We consider a system with one VNF that consists of two microservices. We assume that these microservices container pods are deployed on the same Data Center. For load balancing reasons, each microservice initially has three identical pod replicas. Other parameters are given in Table . I. After 20000 simulation iterations, the results converge well. It took from half to two hours (depending on detection interval) to run these 20000 simulations in CPN tools on a personal computer equipped with Windows 10, 2.10 GHz CPU, and 8GB memory. Indeed, the computation time is proportional to the number of pods and inversely proportional to the periodsecond.

V. CONCLUSION

This paper presents a Petri Net-based model to analyze the performance and resilience of 5G and beyond networks. This model divides a telecommunication system into multiple layers and proves its ability to describe new features of 5G and beyond. The results of the Monte-Carlo simulation on VNF self-healing show the prospects of this model on telecommunication network availability analysis.

The results remain optimistic since other risks such as network failure or maintenance are not fully considered. In addition, more precise parameters need to be collected from our experts and suppliers. For the next step, the auto-scaling case study will be carried out to complete the service-level reliability and resilience analysis, and to see how networks adapt to different packet traffic. We also intend to expand the case study from one single VNF to an SFC and apply the model to simulate a real use case from the verticals.

Fig. 1 .

 1 Fig. 1. VNFs and their microservices.

D. Telecommunication network modeling 1 )

 1 5G Service Function Chain: This is the top layer Petri Net which represents the process of an SFC containing m VNFs as presented in Fig. 2. This pipeline style Petri Net consists of a set P of 2(m + 1) places and a set T of 2m + 1 transitions. It signifies the progress of packet processing.

Fig. 2 .

 2 Fig. 2. Petri Net of SFC example.At this level, a token refers to a packet that conveys a request message that needs to traverse the SFC. The color function

Fig. 3 .

 3 Fig. 3. Petri Net of VNF decomposition.

Fig. 4 .

 4 Fig. 4. Petri Net of Mirco-service treatment.

Fig. 5 .

 5 Fig. 5. Pod Self-healing process

Fig. 6 .

 6 Fig. 6. Microservice availability with immediate detection.

Fig. 7 .

 7 Fig. 7. Overall availability under immediate, 2, 5, and 10 seconds detection.

95% confidence interval [99.9993520304%, 99.9993526071%]

95% confidence interval [99.9987192945%, 99.9987202301%]

A microservice is considered available at time t, if the token quantity of such microservice at p available pods is greater than the desired replica quantity. The uptime of a microservice is the duration of time that a microservice is available. Then the average availability of a microservice i can be calculated as:

microservice i uptime total simulation time In the first situation, we assume the self-healing detection is immediate, i.e., a failure on pod or node can be detected with no delay. We simulate the microservice behavior over 50 years. The average value of microservices' uptime over 20000 simulations is taken as the final result. We assume that the two microservices are from the same VNF supplier and are managed by the same Kubernetes Master. The results in Fig. 6 show that if the desired replica quantity is three, then the availability of a single microservice is 99.9996712% (5 nines). If the desired replica quantity is one (one pod is enough, but three initial pods bring high redundancy), then the availability of this microservice can achieve up to 9 nines. The overall availability for the VNF (at least three available replicas for both microservice 1 and 2) is 99.9993523% 1 .

In the second situation, the effect of self-healing probe frequency t p on system availability is studied. The result is shown in Fig. 7. We compare the overall VNF availability for t p varying from 0 to 10 seconds. The longer the probe periodsecond, the lower the overall availability. The availability drops from 99.9993523% (5 nines) to 99.9987198% 2 (4 nines) by changing immediate detection to 10 seconds. Thus, the telecommunication network can consume less energy while satisfying the availability requirement by wisely optimizing the periodsecond if allowed, according to this result.


Q Bizot

O Politano

V Turlo

F Baras

Molecular dynamics simulations of nanoscale solidification in the context of Ni additive manufacturing

Keywords: additive manufacturing, molecular dynamics simulations, nickel, directional solidification, grain growth

Solidification is a key step in additive manufacturing technology because it determines the microstructure and performance of the final product. Our work provides the nanoscale description of relevant solidification processes for polycrystalline Ni by means of molecular dynamics simulations. In particular, we focus on the thermal effects associated with the characteristic non-stationary conditions of additive manufacturing. Directional solidification and homogeneous nucleation are investigated as a function of operating parameters that control the temperature gradient and the cooling rate. We show that a planar solid/liquid interface propagates at a constant speed, in the presence of a temperature gradient between the melt pool and the solidified region. By adding cooling to the melt pool, a columnar-to-equiaxed transition can be captured at the nanoscale. If the undercooling is not sufficient to promote nucleation, a strong instability of the interface develops, forming protusions. The behaviors observed at the nanoscale are interpreted in terms of the classical theories of solidification and nucleation.

Introduction

Metal additive manufacturing (AM) is a promising technology that enables complexshaped parts with tailored properties to be produced using layer-by-layer deposition. This technology uses a local high-power heat source, such as a laser or electron beam, to melt metallic material in the form of powder or wire. The melted material solidifies as the heat source moves away from the molten region [START_REF] Kim | Literature review of metal additive manufacturing defects[END_REF][START_REF] Lewandowski | Metal Additive Manufacturing: A Review of Mechanical Properties[END_REF][START_REF] Tepylo | Laser-Based Additive Manufacturing Technologies for Aerospace Applications[END_REF][START_REF] Anant Pidge | Additive manufacturing: A review on 3 D printing of metals and study of residual stress, buckling load capacity of strut members[END_REF]. The main advantage of AM compared to conventional manufacturing is that it is possible to design complex shapes in one step, while reducing waste, and eliminating assembly time and cost. Innovative metal AM has already found numerous applications in the automotive, aerospace, electronics, and biomedical industries [START_REF] Vafadar | Advances in Metal Additive Manufacturing: A Review of Common Processes, Industrial Applications, and Current Challenges[END_REF][START_REF] Bhavar | A review on powder bed fusion technology of metal additive manufacturing, Additive Manufacturing Handbook[END_REF][START_REF] Ahn | Direct metal additive manufacturing processes and their sustainable applications for green technology: A review[END_REF][START_REF] Gisario | Metal additive manufacturing in the commercial aviation industry: A review[END_REF].

The columnar-to-equiaxed transition (CET) is central to the control of the solidified microstructure. In alloy solidification, the main parameters governing CET are the temperature gradient, the velocity of the solidification front, the cooling rate, and constitutional undercooling. The critical gradient condition for fully steady-state equiaxed growth with a given density of inoculants has been estimated by Hunt [START_REF] Hunt | Steady state columnar and equiaxed growth of dendrites and eutectic[END_REF] and is commonly used as a rough estimate of CET. In order to deal with the specific characteristics of AM, extensive experimental work (see, for example, [START_REF] Bermingham | Promoting the columnar to equiaxed transition and grain refinement of titanium alloys during additive manufacturing[END_REF]) and theoretical studies have been carried out to predict favorable CET conditions. Modeling AM processes requires a multiscale approach: the macroscopic scale to record thermal history, the mesoscopic scale to predict melt pool dynamics, and the microscale to describe microstructure formation. Process modeling based on finite element methods or computational fluid dynamics often integrates experimental measurements to predict process defects such as cracks and pores as well as to estimate surface roughness of the final parts [START_REF] Tan | Microstructure modelling for metallic additive manufacturing: a review[END_REF]. At the microscale, stochastic models such as cellular automata and kinetic Monte Carlo can provide micrographs of grain structures close to those observed experimentally. Such models nevertheless leave aside many physical phenomena, relying on effective parameters that are not transferable, but have to be fitted for each new material of interest [START_REF] Spittle | Columnar to equiaxed grain transition in as solidified alloys[END_REF]. In contrast, phase field modeling provides a much more comprehensive description of solidification at the microscale that involves the formation of dendrites, microsegregation, and precipitation [START_REF] Kurz | Progress in modelling solidification microstructures in metals and alloys. Part II: dendrites from 2001 to 2018[END_REF]. Phase-field can also be combined with cellular automaton to study competitive growth of columnar dendritic grains under temperature gradient [START_REF] Dorari | Growth competition between columnar dendritic grains -The role of microstructural length scales[END_REF]. The multiscale phase field method can also be used to describe heterogeneous nucleation, grain selection, and epitaxial growth to assess the role of AM parameters in CET [START_REF] Liu | Insight into the mechanisms of columnar to equiaxed grain transition during T metallic additive manufacturing[END_REF]. However, at nanoscale, few models report CET, except for a recent work by Kurian et al. [START_REF] Kurian | Selective laser melting of aluminum nano-powder particles, a molecular dynamics study[END_REF], which studied the micro-selective laser melting process using molecular dynamics (MD) simulations in order to investigate melting and solidification of a randomly-distributed aluminum nano-powder bed.

Many solidification processes have nevertheless been investigated at the nanoscale by means of MD simulations. The pioneering work of Lu and Spzunar [START_REF] Lu | Molecular-dynamics simulation of rapid solidification of aluminum[END_REF] described the structural changes induced by rapid solidification. Fast quenching leads to a nonequilibrium glassy state, while a slow cooling rate results in a crystalline state [START_REF] Tian | Molecular dynamics simulation for cooling rate dependence of solidification microstructures of silver[END_REF][START_REF] Pei | The rapid solidification of Ti3Al: a molecular dynamics study[END_REF][START_REF] Hou | Cooling rate dependence of solidification for liquid aluminium: a large-scale molecular dynamics simulation study[END_REF]. The characteristics of a solidification front have been measured by MD in two-phase metallic systems under isothermal conditions for different orientations of the solidliquid interface in order to obtain the kinetic coefficients [START_REF] Celestini | Measuring kinetic coefficients by molecular dynamics simulation of zone melting[END_REF][START_REF] Sun | Kinetic coefficient of Ni solid-liquid interfaces from molecular-dynamics simulations[END_REF]. The case of directional solidification has been studied using non-equilibrium MD simulations for Lennard-Jones systems [START_REF] Celestini | Nonequilibrium molecular dynamics simulation of rapid directional solidification[END_REF], Al-Cu alloys [START_REF] Mahata | Effects of solidification defects on nanoscale mechanical properties of rapid T directionally solidified Al-Cu Alloy: A large scale molecular dynamics study[END_REF], and CrNi-alloyed steels [START_REF] Bahramyan | Nano-scale simulation of directional solidification in TWIP stainless steels: A focus on plastic deformation mechanisms[END_REF]. Large-scale MD simulations with several million atoms have been used to study the very first steps of solidification in various metallic systems. Homogeneous nucleation has been studied at constant undercooling temperature or by applying a constant cooling rate on a melted metal [START_REF] Shibuta | Heterogeneity in homogeneous nucleation from billion-atom molecular dynamics simulation of solidification of pure metal[END_REF][START_REF] Mahata | Understanding homogeneous nucleation in solidification of aluminum by molecular dynamics simulations[END_REF]. Athermal heterogeneous nucleation through grain refiners in under- cooled metal has also been investigated by MD simulations [START_REF] Fujinaga | Molecular dynamics simulation of athermal heterogeneous nucleation of T solidification[END_REF].

Solidified metal

In AM, thermal effects are predominant (see Fig. 1). Indeed, the melt pool is subjected to radiation and Marangoni convection at the same time, inducing an effective cooling rate on the melt side. In addition, the temperature gradient between the melt pool and the solidified region produces thermal transfer by conduction on the solid side. At least two different solidification processes coexist: directional solidification due to the temperature gradient, and homogeneous nucleation in an undercooled liquid. Heterogeneous nucleation can also be observed through the presence of leftover grains and grain refiners. In the present work, we developed MD simulations of pure nickel submitted to typical non-stationary thermal conditions specific to AM in order to study the columnar-to-equiaxed transition. Among metals and alloys, nickel-based alloys are commonly used in AM [START_REF] Bandyopadhyay | Recent developments in metal additive manufacturing[END_REF]. As this study is restricted to pure metal, the microstructure will result exclusively from thermal effects occurring during solidification processes. The first representative model (model A) was set up to describe the rapid directional solidification of a textured Ni substrate. The second model (model B) was used to investigate the effect of the cooling rate on the evolution of the microstructure. This approach allowed us to assess the effects of process parameters on the resulting microstructure at the nanoscale. In addition, MD simulations were compared to predictions based on the classical theory of solidification of pure metals and on classical nucleation theory.

Methods

The simulations were performed using the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software [30], and the Embedded Atom Method (EAM) potential, developed by Purja Pun and Mishin for Ni-Al systems [START_REF] Purja Pun | Development of an interatomic potential for the Ni-Al system[END_REF]. This potential has been extensively used to study reactivity in Ni-Al nanofoils [START_REF] Politano | Reaction front propagation in nanocrystalline Ni/Al composites: A molecular dynamics study[END_REF][START_REF] Politano | Molecular Dynamics Studies in Nanojoining: Self-Propagating Reaction in Ni/Al Nanocomposites[END_REF]. In addition, this potential is known to reproduce well the physical and thermodynamical properties of nickel [START_REF] Turlo | Comparative study of embedded-atom methods applied to the reactivity in the Ni-Al system[END_REF] (See We considered a pseudo 2D system in a simulation box of 100 nm × 1.4 nm × 200 nm. As shown in Fig. 2, the system was divided into three regions along the z direction: a lower region located between 0 and 20 nm kept at a fixed temperature and named the substrate, a middle region from 20 to 50 nm, and an upper region that was effectively affected by the laser beam. The y direction was slightly larger than twice the cut-off distance (r cut ) of the potential to avoid multiple interactions. Periodic boundary conditions were applied in the x and y directions, while the boundary condition was nonperiodic and shrink-wrapped in the z direction. The system of 25 grains was built with the Atomsk software [START_REF] Hirel | Atomsk: A tool for manipulating and converting atomic data files[END_REF]. The procedure used to create this polycrystalline sample is analogous to that used to obtain Voronoi cells, which has already been applied in previous works [START_REF] Politano | Reaction front propagation in nanocrystalline Ni/Al composites: A molecular dynamics study[END_REF][START_REF] Perron | Oxidation of nanocrystalline aluminum by variable charge molecular dynamics[END_REF]. Each grain had its [010] axis aligned with the y direction and could rotate along the y axis. The equations of motion were integrated using the velocity Verlet algorithm with a timestep of 1 fs. The atomic positions were visualized with OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO-the Open Visualization Tool[END_REF]. The local atomic environment (i.e. fcc, hcp, or unknown) was determined using Polyhedral Template Matching (PTM) [START_REF] Larsen | Robust structural identification via polyhedral template matching[END_REF].

Two representative models were prepared: with and without cooling. Different temperature values of the substrate were considered in order to explore a broad range of microstructures. Typically, the first model (A) allowed us to understand rapid direc-tional solidification in a temperature gradient, without cooling. In the second model (B), a cooling rate was applied to force the system to cool down at a controlled rate, as observed in AM solidification processes (see [START_REF] Bahramyan | Nano-scale simulation of directional solidification in TWIP stainless steels: A focus on plastic deformation mechanisms[END_REF]).

In model A, the system was first equilibrated at T sub for 200 ps in the NPT (isothermalisobaric) ensemble, followed by another equilibration in the NVT ensemble for 100 ps. After thermalization, the two extremities of the sample were thermostated in an NVT ensemble for 2 ns. The laser region extended from 50 to 200 nm (see Fig. 2) and its temperature was set at 1750 K. The temperature of the substrate was set at T sub . Note that, in the central part of the sample, atoms evolved freely as their equations of motion were integrated into the NVE ensemble. After 2 ns, the hot thermostat was removed but the cold thermostat was kept until the end of the simulation. We considered 2 different temperatures for the substrate: T sub = 300 K and 1000 K, in order to consider two very different gradients.

In model B, the system was first equilibrated at T sub for 200 ps in the NPT (isothermalisobaric) ensemble, followed by another equilibration in the NVT ensemble for 100 ps. After thermalization, the two extremities of the sample were thermostated as in model A. The laser region was set at 1750 K whereas the substrate was set at T sub for 2 ns. Then the laser thermostat was removed and replaced by a temperature ramp corresponding to cooling rates ranging from 75 K/ns to 750 K/ns (see Table 2 for more details). The substrate thermostat was kept until the end of the cooling process. At the end of the cooling ramp, which is when the whole system reaches the substrate temperature, T sub , the two thermostats were removed, and the simulation of the entire system was carried out in the NVE ensemble. Six different temperatures of T sub were considered: 800 K, 900 K, 1000 K, 1100 K, 1200 K, and 1300 K.

Molecular dynamics simulations enable characterization of solidification processes by means of thermal effects and nucleation. In order to evaluate quantities such as critical nuclei, nucleation rate, or thermal fluxes, it is necessary to estimate structural parameters (latttice parameter), thermodynamic properties (melting temperature, solid/liquid interfacial energy, latent heat of fusion, heat capacity, etc.) and thermal transport (heat conductivity). For this purpose, we developed separate MD simulations (see Turlo et al. [START_REF] Turlo | Comparative study of embedded-atom methods applied to the reactivity in the Ni-Al system[END_REF] for more details). The parameters for Ni are reported in Table 1. The structural order of the Ni crystalline phase is fcc. The temperature evolution of the lattice parameter a 0 was evaluated, and its value at melting point, a(T m ), is reported. The melting temperature was determined with the isothermal two-phase method. The heating-cooling cycle procedure was used to obtain the other thermodynamic properties: the bulk heat of fusion L V and the heat capacity C p (see Appendix A). The heat capacity, which is proportional to the slope of the relative dependence of the enthalpy on temperature, was in good agreement with experimental data. Reasonable agreement between experimental value and MD calculation was also obtained for the latent heat of fusion. ). Experimental values are given for information only. Most of these values are discussed in Thurnay [START_REF] Thurnay | Thermal properties of transition metals[END_REF]. Temperature-dependent parameters are evaluated at the bulk melting temperature T m . † Cohesive energies are evaluated at 0 K.

be understood if one considers that the conduction of heat observed in this MD study with the classical potential accounts only for phonon transport and does not account for electronic transport. In the present approach, phonon thermal conductivity was computed with non-equilibrium molecular dynamics (NEMD), see Appendix B for details.

In the case of our Ni system, the MD value led to an effective characteristic length of heat conduction that was shorter than the real one by a factor D exp th /D MD th = 4.2. This can be considered as an asset in simulations because we can investigate thermal effects on MD scales, which allows the refinement and validation of existing theories and the development of new theories for nanoscale systems [START_REF] Turlo | Modeling self-sustaining waves of exothermic dissolution in nanometric Ni-Al multilayers[END_REF]. After the theoretical model has been refined, experimental values for material properties can be inserted to generate reliable predictions, which can be further validated experimentally. In addition, such theoretical models solved analytically or numerically are much more computationally efficient, which allows us to cover much greater length and time scales, well beyond the MD capabilities.

The solid-liquid interfacial energy γ sl was determined using the critical nucleus method [START_REF] Xia | Molecular dynamics studies on the correlation of undercoolability and thermophysical properties of liquid Ni-Al alloys[END_REF]. Details are given in Appendix C.

Results

In AM, laser heating imposes dynamic temperature conditions that affect the solidification process. In order to understand their influence on microstructure formation, we first considered the role of a temperature gradient between the solidified material and the melt pool (model A). Next, we introduced the effect of global cooling (model B).

Model A

After the complete melting of the upper region through laser heating, the system was divided into 3 regions: a polycrystalline substrate from 0 to 20 nm maintained at T sub , a melt pool at 1750 K (> T m (Ni)), and a solid region between the two (middle region). The resulting temperature profile at t = 0 ns, after the passage of the heat source, is depicted in Fig. 3, showing that the steep temperature gradient in the solid part was established between 20 and 80 nm. The abrupt decrease in the fraction of fcc atoms in Fig. 3 corresponds to a solid/liquid interface. In the specific initial microstructure of this simulation, the system was composed of three columnar grains at the interface, as shown in the snapshot at t = 0 ns.

After the release of the thermostat from the melt pool, solidification proceeded with propagation of a planar solid/liquid interface in the z direction. We observed that the melted region shrank while the solidified region increased. This progressively induced a smoother temperature gradient in the solidified region. As one can see in Fig. 3, the fraction of fcc atoms in the solid region is not strictly equal to 1 because of grain boundaries and point defects. Thus, the thickness of the front was measured as the zone with a fraction of fcc atoms between 0.02 and 0.95, which never exceeded 10 nm. The position of the solidification front, associated with the decrease in the fraction of fcc atoms, was defined as the point where the fraction of fcc atoms is equal to 0.8. This intersection corresponds to 1650 K, which is slightly lower than the melting point of Ni (1701 K). The resulting position of the solid/liquid interface is plotted in Fig. 4 as a function of time for two distinct substrate temperatures of 300 K and 1000 K. The crystallization velocity was found to be strongly dependent on the temperature of the substrate but with common features, including a change in slope after some time. By analyzing the regions with linear time dependence on the position of the solid/liquid interface, we extracted the corresponding solidification velocities, listed in Fig. 4. During a transient period (before 3 ns for T sub = 300 K and 5 ns for T sub = 1000 K), the instantaneous velocities were larger than the stationary ones by 2.5-3 m/s. In the stationary regime, the solidification front velocities were reduced to 11.5 m/s for T sub = 300 K and 5.5 m/s for T sub = 1000 K, suggesting a linear relationship between substrate temperature and steady-state solidification front velocity. Extrapolating the linear trend set by these two points to zero velocity would result in a temperature of 1642 K, close to the melting temperature of Ni (1701 K). Solidification in pure metals is entirely controlled by thermal processes [START_REF] Porter | Phase Transformations in Metals and Alloys[END_REF]. In the absence of convection, the latent heat of solidification released during solidification is transported by conduction away from the solid/liquid interface. Conduction can take place through either the solid or the liquid, depending on the temperature gradients at the interface. For solid growth at a velocity v with a planar interface, the balance between heat flow through the solid, heat flow through the liquid, and latent heat released at the interface reads where G s = dT dz sol is the local temperature gradient in the solid, G ℓ = dT dz liq is the local temperature gradient in the liquid, and L V is the latent heat of fusion. In Model A, the liquid is superheated, and temperature gradients at the interface in the liquid and the solid are both positive. If we assume that thermal conductivity is very similar in the liquid and in the solid κ sol = κ liq = κ, eq. ( 1) gives an estimate of the velocity v:

κ sol G s = κ liq G ℓ + vL V (1)
v = κ L V [G s -G ℓ ] (2)
Since the local gradient in the melt pool is weak, the velocity mainly depends on the temperature gradient in the solidified region. The relationship (2) explains why the velocity is greater at the beginning of the solidification process when the gradient is steeper (see Fig. 4). We also understand why the crystallization velocity increases for decreasing substrate temperatures, since the temperature gradient is larger in this case.

In MD studies, the crystallization velocity is usually estimated at constant temperature using the two-phase method [START_REF] Purja Pun | Development of an interatomic potential for the Ni-Al system[END_REF]. The crystallization velocity as a function of temperature and interface orientation was evaluated using this method (see Appendix D). The velocity depends on the orientation of the Ni interface: v(001) > v(101) > v(111) and the degree of undercooling. The difference in velocity as a function of orientation becomes smaller as it approaches the melting temperature. In the case of directional solidification, investigated in the present work (model A), we estimated the solidification velocity as a function of orientation in the case of an initial system composed of a single grain. The velocity was almost the same for the three orientations: v(001) ∼ v(101) ≳ v(111). This is not surprising because the solid/liquid interface corresponds to a temperature value close to melting point. It is also supported by the fact that the velocities measured in the polycrystalline samples (see Fig. 4) are in the same range as the crystallization velocity measured by the two-phase method between 1600 K and 1700 K. It is also important to note that the interface (111) was the least mobile and, thus, was not considered in our polycrystalline sample, as crystal rotations around the direction [010] did not allow us to establish such an interface. The crystallization velocity in the case of directional solidification (model A) in polycrystalline samples is fixed by the temperature gradient rather than by the temperature value at the solidification front or by the grain orientation. This observation corroborates the fact that the front remained flat in the case of a polycrystalline system.

T sub (K) α (K/ns)

Model B

In Model B, a progressive global cooling was imposed in order to mimic the thermal conditions experienced during AM. As in the previous model, A, the system was composed of a substrate maintained at a given temperature, T sub , a melt pool, and a solidified region between the substrate and the melt pool (middle region). Both the solidified and melted regions were affected by the progressive cooling. Two parameters define the thermal conditions: the temperature of the substrate, T sub , and the cooling rate α mea-sured in K/ns. All the conditions considered in the present work are reported in Table 2. The temperature profile and the associated fraction of fcc atoms are depicted in Fig. 5 at different times for T sub = 1000 K and α = 375 K/ns. Time t = 0 ns corresponds to the time at which the heat source was released. We observed a solidified region ranging from 15 to 85 nm between the substrate and the solidification front. The solidification front separating the solid and liquid regions was located at 1650 K, close to the melting point. A temperature gradient of about 8.75 K/nm was established in this solid region, whereas most of the liquid region was at the melting temperature of Ni. Because of the imposed cooling, the liquid went from a superheated to a supercooled state. The temperature profile at 0.8 ns showed a peak reflecting the exothermicity of the solidification process, whereas most of the liquid region was at a temperature lower than the melting point of Ni. The solidification front was located precisely at the peak of the temperature front. From 0.8 ns to 1.6 ns, the solidification front propagated forward, resulting in columnar grain growth. At 1.6 ns, a small peak of fcc phase was observed in the liquid region, reflecting the appearance of a nucleus of crystalline phase. After 1.6 ns, grains appeared uniformly and randomly in the liquid region. At 2 ns, the exothermic peak in the temperature profile was less pronounced because of the homogeneous nucleation and growth of the crystalline grains, increasing the overall temperature in this previously liquid region. After 2.4 ns, the upper region was completely filled with randomly oriented and randomly distributed Ni grains. By comparison with model A, the cooling imposed in model B induced a drastic change in the microstructure, with the coexistence of columnar and equiaxed grains (see Fig. 5b).

In the case of a substrate at T sub = 1000 K, the influence of the cooling rate α was evaluated to be in the range of one order of magnitude from 75 K/ns to 750 K/ns. We observed that the cooling rate modified the time at which the first nucleus appeared (t 2 in Table 2). The higher the cooling rate, the faster nucleation occurred. Consequently, the zone of equiaxed grains ∆L z was thicker for greater cooling rates. The position of the solidification front as a function of time is reported in Fig. 6. The typical curve is convex, in comparison with Model A (see Fig. 4). This means that the instantaneous velocity progressively increased with the degree of undercooling of the melt pool. Furthermore, greater cooling rates correspond to faster solidification. After first nucleus appearance (black diamond in Fig. 6), columnar grain solidification ceased in favor of equiaxed grains.

The only difference between model A and model B was the imposed cooling rate. In model A, the solidification developped in a superheated fluid and the front kept its planar shape. In contrast, in model B, the solid grew in a supercooled liquid and the heat was conducted away from the front in the liquid due to the negative gradient G ℓ < 0, as shown in Fig. 6, in the case of a cooling equal to 375 K/ns. Equation (1) still holds and the velocity reads

v = κ L V [G s + |G ℓ |] (3)
The negative temperature gradient in the liquid G ℓ significantly decreased until it reached a minimum when the first nucleus appeared. This observation explains why the instantaneous velocity increased during the associated solidification stage. After that stage, the heat released by the formation of solid grains in the liquid thwarted the cooling, which became less efficient. In Fig. 5 at t = 2 ns, we noted that the peak corresponding to solidification was less pronounced, when compared to previous times. Each grain for- mation released heat in the liquid region and, consequently, the melt pool temperature increased. Figure 7 shows the instantaneous solidification front velocity of columnar grains in representative cases of Table 2. For a given cooling rate (α = 375 K/ns), the substrate temperature, T sub , has a limited impact on front velocity. This is due to the fact that G s in eq. ( 3) is the local temperature gradient evaluated on the left side of exothermic peaks (see Fig. 5a) rather than the global temperature gradient established in the solidified region. In the case of columnar and equiaxed grain microstructure, velocity increased as a function of time. This is related to the increase in magnitude of the temperature gradient in the liquid G ℓ during cooling. The increase in v is thus more pronounced for greater cooling rates. Due to the formation of equiaxed grains in the melt, the front associated with columnar grains either stopped or slowed down before stopping. For a very low cooling rate (75 K/ns), the velocity was almost constant over 4 ns before a slight increase. The front reached the system edge before the degree of undercooling (∆T = T m -T = 450 K) was sufficient for homogeneous nucleation. The crystallization velocity is also a function of the degree of undercooling (see Appendix D for more details):

v = K T m -T T exp - Q k B T ( 4
)
where Q is the activation enthalpy, k B is the Boltzmann constant, and K is a constant. According to eq. ( 4), the velocity increases as a function of the degree of undercooling T m -T . In contrast, the formation of equiaxed grains, which occurred at around the same degree of undercooling (500-600 K) for all systems (see Table 2), led to a drastic change in the local temperature gradient due to the release of heat associated with nucleation. The nucleation provoked a slowing down of the velocity at a characteristic time corresponding to the target degree of undercooling. As a consequence, the local temperature gradient in the liquid reached about the same minimum value (-20 K/nm) whatever the cooling rate or substrate temperature (see Fig. 6). The velocity at columnar-to-equiaxed transition was in the range of 60-65m/s, corresponding to this specific degree of undercooling.

During the propagation of solidification from 0 to 1.6 ns, the front progressively lost its planar shape, developing a characteristic rounded grain shape at the interface. This is also reflected by the greater thickness of the front. We observed a typical microstructure of columnar and equiaxed grains for most of the parameters considered, except for a high substrate temperature value (see Table 2). Homogeneous nucleation began when the degree of undercooling of the melt pool reached a value in the range 500 K -600 K, whatever the substrate temperature T sub and the cooling rate α. Is it possible to interpret this observation in terms of Classical Nucleation Theory (CNT)? In CNT, the rate of homogeneous nucleation I reads [START_REF] Mahata | Understanding homogeneous nucleation in solidification of aluminum by molecular dynamics simulations[END_REF]

I = I 0 exp - A T 3 (∆T ) 2 (5)
with

A = 16πγ 3 SL T 4 m 3k B L 2 V (6) (a) (b)
Figure 8: Normalized nucleation rate I/I 0 (eq. ( 5)) (a) and critical radius (b) as a function the degree of undercooling.

where I 0 is a coefficient that depends on the interface temperature and free energy. The factor A is a constant that depends on the solid-liquid surface energy and latent heat of solidification. The normalized rate I/I 0 is plotted in Fig. 8a as a function of the degree of undercooling for the parameter values corresponding to the MD simulations (see Table 1). The nucleation rate strongly depends on the degree of undercooling. No nucleation was expected below ∆T = 500 K. The critical temperature T cr corresponds to the maximum of the nucleation rate: T cr = 3T m /5 = 1026 K, i.e. ∆T = 684 K. This estimate is in good agreement with the range of undercooling where nucleation was observed in MD simulations. The experimental estimate is that the critical nucleation temperature is between 0.5 and 0.6 times the melting temperature.

According to CNT, the critical radius r * of a spherical nucleus depends on the degree of undercooling ∆T = T m -T :

r * = 2γ sl L v T m ∆T (7)
The corresponding number of atoms in a critial nucleus is

N * = 4 × 4π 3 r * 3 /a 3 0 ( 8
)
where the factor 4 stands for the number of fcc atoms/unit cell. In simulations, the degree of undercooling ∆T , at which nucleation appears, is in the range [500-600 K]. The critial radius r * is thus between 0.66 nm and 0.75 nm, corresponding to N * between 100 and 160 atoms. The critical radius predicted by CNT with MD parameters is plotted in Fig. 8b for the degree of undercooling corresponding to the appearance of the first nucleus. We note that the critical nucleus is almost the same in all simulations and smaller than the size of the simulation box along the y axis, thus avoiding any simulation artifacts. Nucleus formation dynamics is of primary importance in the final microstructure after solidification. In order to describe the dynamics, we analyzed the instantaneous system state at different times after the first nucleation during the growth process. Nuclei were identified by deleting all non-fcc/liquid atoms and then using the cluster analysis available in Ovito. A nucleus is defined here as a set of connected atoms, each of which is within the reach of the other atoms in the same cluster. The cutoff distance was set at 2.5 Å, which is the distance between nearest neighbors in the fcc structure. The results are presented in Fig. 9. From 1.6 ns, the number of clusters increased, reaching its maximum at 1.8 ns. A maximum of fifteen clusters was observed. Note that the number of clusters does not correspond to the rate of nucleation (eq. ( 5)): it corresponds to the cumulative nucleation rate minus the clusters that disappeared due to coalescence. The mean cluster size remained lower than 5000 atoms. It corresponds to small clusters that are larger than the critical size. Note that the critical nucleus size decreases during system cooling, which ends at 2.5 ns (see eq. ( 7)). We also followed the size of the largest cluster. At the very beginning, there is competition between clusters of similar size and the largest one may differ from time to time. At 1.75 ns, the largest cluster was clearly identified. It grew significantly because of adhesion of Ni atoms that diffused in the liquid and came into contact with existing clusters. After 1.8 ns, the number of clusters decreased following coalescence between isolated clusters or with columnar grains. Coalescence was associated with an increase in mean cluster size. Fluctuations in cluster numbers correspond to a detach/attach process between 2 or more adjacent grains. After 2 ns, the number of clusters decreased and their mean size stabilized. During the growth process of equiaxed grains, columnar grains extended in the direction of solidification (z > 0), but their growth progressively slowed due to the growth of equiaxed grains. Fluctuations in the size of the largest cluster reflect the attachment/detachment of the grain with columnar grains. At 2.3 ns, all grains were connected. At 2.5 ns, most of the system was solidified, with the coexistence of columnar and equiaxed grains. Grain boundaries were occupied by atoms with no well-defined structure (termed unknown atoms) and can be considered disordered grain boundaries.

Figure 5b shows that nuclei appeared randomly in the melt. In the front view provided by the snapshots, we first see the formation of round nuclei. These initial critical spherical nuclei (r < 0.7 nm) grew in the form of disks due to their extent being limited in the y direction. Each grain is made of fcc (green) and hcp (red) atoms. As the difference in cohesive energy per atom for fcc and hcp is very small (less than 0.03 eV; see Table 1), any thermal fluctuation may cause the formation of hcp Ni atoms, leading to the formation of nanotwinned grains. Let us now consider the case of a substrate at 1300 K (see Fig. 10). Here we considered a system twice as long in the z direction (L z = 426 nm). At t = 0 ns, laser heating was removed, and the region ranging from 50 nm to 426 nm was submitted to a cooling rate of 375 K/ns. The temperature profile at t = 0.4 ns depicted in Fig. 10a corresponds to a liquid region at 1600 K and a solid region submitted to a temperature gradient. The limit between the solidified region and the liquid pool is delineated by the abrupt decrease in the fraction of fcc atoms. In about 1 ns, the melt pool cooled down to the substrate temperature T sub = 1300K. This corresponds to an undercooled melt pool. The degree of undercooling ∆T = 410K did not allow spontaneous nucleation (see Fig. 8a). The nucleation rate is almost zero, contrasting with the formation of nuclei observed for a substrate temperature at 1000 K. When cooling ended, the system evolved into adiabatic conditions (NVE ensemble). Figure 10a depicts the temperature profiles and fraction of fcc atoms at different times under these specific conditions. At 1.4 ns, the temperature peak due to the release of latent heat was evacuated into the solid and liquid parts. At 2.2 ns, this peak became higher, while the temperature gradient in the solidified region became smoother. For later times, this tendency became more pronounced, with a smooth gradient in the solid part and a sharp one in the liquid around the solidification front.

In solidification theory [START_REF] Porter | Phase Transformations in Metals and Alloys[END_REF], it is well known that when a solid grows in an undercooled liquid, a planar solid/liquid interface is unstable. At 2.2 ns, we observed that protrusions located in the central columnar grain developed at the interface to form a cellular structure with fingers (snapshot in Fig. 10). At later times (for instance, 3 ns and 5 ns), these fingers extended, with a marked bottom and tip. In Fig. 10b, we plotted the temperature profiles and the fraction of fcc calculated in a slice containing the tip or the bottom of the finger. A striking fact is that the thickness of the solidification front between the tip and the bottom increased significantly between 3 and 5 ns. In Fig. 10b, we see that the heat released at the tip can only be conducted in the liquid (Fig. 10b circle in red) when the local gradient is very steep. Local temperature T tip ∼ 1450 K is lower than the melting point. The heat released at the finger bottom is conducted in both solid and liquid parts. The local temperature T bottom ∼ 1698 K is close to the melting point, T m (Fig. 10b -circle in red). The positions of the tip and the bottom are plotted as a function of time in Fig. 11a. As expected, the velocity of the tip v tip is greater than that of the bottom: v bottom . The difference in position between the tip and the bottom increases as a function of time, indicating the progressive extent of the fingers observed in the snapshots. Heat is conducted away from the tip in the liquid, due to the local steep temperature gradient in the liquid part. Due to the Gibbs-Thomson effect [START_REF] Porter | Phase Transformations in Metals and Alloys[END_REF], the solidification velocity at the tip reads

v tip = κ liq L V 1 r ∆T 0 1 - r * r (9)
where ∆T 0 is the degree of undercooling, r * is the corresponding critical radius, and r is the characteristic size of the tip. Given the parameters (see Table 1), the theoretical value of the tip velocity is given in Fig. 11b. The measured value of the tip is close to 50 m/s, corresponding to a characteristic size of 6 nm, in good agreement with our observations. As shown in Table 2, the microstructure observed after solidification depends on both substrate temperature, T sub , and cooling rate, α. The value of T sub modifies the steepness of the temperature gradient in the solid region between the substrate and the solidification front. This temperature gradient becomes smoother as the solid/liquid interface moves forward. In addition, since T sub is the temperature of the melt pool after completion of cooling, its value gives the degree of undercooling. In the case of low values, T sub ≤ 1200 K, the final structure corresponds to the coexistence of columnar and equiaxed grains. For values T sub ≥ 1300 K, we observed columnar grains and the development of protrusions associated with solidification front instability. In this case, the degree of undercooling was never high enough for homogeneous nucleation, since the rate of homogeneous nucleation is almost zero.

Conclusions

In the present work, molecular dynamics (MD) simulations were used to provide in situ observations of solidification during additive manufacturing (AM). The thermal conditions to which the system is submitted during AM are dynamic. The temperature gradient between the solidified region and the melt pool is constantly evolving. Thus, conduction of heat away from the solid/liquid interface takes place in non-stationary conditions. In addition, thermal radiation and convective cooling of the melt pool induce a progressive cooling of the liquid region. Molecular dynamics simulations allow us to reproduce these conditions in a nano-box system corresponding to the close vicinity of the solidification front.

We first analyzed the case of directional solidification in a non-stationary temperature gradient. The steepness of the temperature gradient is directly related to the substrate temperature T sub . A low T sub value increases the temperature gradient in the solid region (G s ), leading to rapid propagation of the solid/liquid interface in the direction of the gradient. The heat release by conduction in the solidified region determines the front characteristics (velocity and stability). The microstructure is a typical columnar structure. The solid/liquid interface remains flat, expect at grain boundaries, where grain boundary grooving was observed.

When the cooling of the melt pool is taken into account, the situation is completely different. Heat is conducted away from the interface in both solidified and melted regions, leading to an intrinsically unstable solid/liquid interface, with the development of rounded columnar grains. In addition, when the undercooling induced by the cooling was sufficient, we observed spontaneous nucleation of seeds in the melt pool and growth of equiaxed grains. These randomly oriented grains impede further propagation of columnar grains.

The faster the cooling, the faster the equiaxed grains appeared, leading to a large region with equiaxed grains. The propagation velocity of columnar grains is directly related to the temperature gradients from either side of the interface, in the solid (G s ) and liquid (G l ) regions. A greater cooling rate promotes the propagation of columnar grains by increasing their velocity. So the same operating parameter (cooling rate) influences the two competing solidification modes (directional solidification and homogeneous nucleation). The maximum degree of undercooling ∆T max depends on the substrate temperature: ∆T max = T sub -T m . If ∆T max is not in the appropriate range, the nucleation rate vanishes. In this case, the columnar grains continued to grow, but thermal instability led to grain fingering, very similar to dendrite formation.

Although the length scale is nanometric, molecular dynamics simulations reproduce the complex behavior associated with solidification during additive manufacturing. With complementary simulations, we were able to evaluate the structural and thermodynamic parameters of Ni associated with the EAM potential [START_REF] Purja Pun | Development of an interatomic potential for the Ni-Al system[END_REF]: latent heat of fusion, solid/liquid interfacial energy, heat conductivity, ... (See Appendices A-C). This study allowed us to compare the behaviors at nanoscale to those expected in the framework of the classical solidification theory. The main conclusion is that predictions of the solidification theory of pure metals still hold at the nanoscale. Another interesting conclusion of the present work is the observation of CET in a pure metallic system. Here, the CET is driven exclusively by thermal conditions. The next step is to consider systems with inoculants in the melt pool, in order to further control the CET, as is usual in experiments. In addition, through the detailed characterization of relevant parameters, we could consider supplementing this study with simulations at mesoscopic scales [START_REF] Kavousi | Quantitative prediction of rapid solidification by integrated atomistic and phase-field modeling[END_REF][START_REF] Kavousi | A temperature-dependent atomisticinformed phase-field model to study dendritic growth[END_REF]. slightly differs from the experimental estimation (2.53 10 9 J/m 3 ). Nevertheless, this indicates good adequacy of the potential to determine the latent heat of fusion. Phonon thermal conductivity was computed with non-equilibrium molecular dynamics (NEMD) [START_REF] Juan | The impact of the thermostats on the non-equilibrium computer simulations of the interfacial thermal conductance[END_REF]. We considered a periodic simulation box elongated in the z direction, with a cross-section of 20 × 20 lattice units (i.e. 7.04 nm × 7.04 nm). This system was equilibrated at T in an isothermal-isobaric ensemble (NPT) over 0.2 ns. Then it was divided into 20 bins. The first bin defines the hot region maintained at T + 50 K, while the middle bin (i.e. 11) is the cold region thermalized at T -50 K. Two distinct Langevin thermostats were used. The Langevin damping parameter was set equal to 0.1 ps. The equations of motions of the atoms outside the thermostated regions were integrated into the NVE ensemble. A run of 0.28 ns was long enough to obtain a steady state temperature profile between heat sink and heat source. The energies added/removed by the two thermostats and the temperature profile between the two regions were monitored and averaged for 0.1 nanoseconds. Thermal conductivity was computed using Fourier's law :

Appendix B

κ = - J dT /dz ( 10
)
where J is the heat flux and dT /dz is the thermal gradient between the two thermostated regions. In the case of solids, finite size effects have to be taken into account. The distance between the two thermostats limits the phonons' mean free path that lowers thermal conductivity compared to bulk values. This was overcome by computing thermal conductivity with samples of different lengths. Bulk thermal conductivity was then estimated by extrapolating reciprocal thermal resistivity 1/κ versus reciprocal system length (1/L z ). Note that, in metallic melts (or liquids) the phonon's mean free path for heat transport is much shorter than the distance between the two thermostats since it is meanly due to collisions between close neighbors. That lead us to simulate between 40 and 100 lattice units (i.e. 14.08 nm to 35.2 nm) along the z-direction and nine temperature values were considered: T = 300 K, 500 K, 800 K, 1100 K, 1300 K, 1500 K, 1900 K, 2100 K, and 2300 K (see Fig. 13). ii One half of the sample was frozen and the other half (i.e between 0 and L x /2) was melted in the NPT ensemble for 1 ns at a temperature T sup = 3000 K, much higher than the melting temperature estimated with the EAM potential, in order to ensure complete melting and produce a liquid phase.

iii Finally, the simulation was carried out in the NPT ensemble for 2.5 ns at the target temperature T for the whole system. If the target temperature is below the melting temperature (T m = 1701 K), the interface moves toward the right and the system crystallizes (Table 15a). If the target temperature T is above the melting temperature, a melting front propagates to the left and the system melts completely.

We considered 7 target temperatures below the melting temperature: 1100 K, 1200 K, 1300 K, 1400 K, 1500 K, 1600 K, and 1700 K. In order to estimate the solid/liquid interface velocity, the front position was followed and compared to the position of a reference atom in the solid phase (Table 3). We observed that the solid/liquid interface velocity depends on the degree of undercooling. The Wilson-Frenkel model was used to approximate temperature dependence (Eq. 4). In addition, the crystallization velocity is consistent with the sequence v (100) > v (110) > v (111), reported in the literature. Table 3: Left: Crystallization velocity (in m/s) calculated with the two-phase method at different temperatures and for different interface orientations. Right: Wilson-Frenkel parameters obtained from the fit by eq. ( 4).

Figure 1 :

 1 Figure 1: Schematic representation of the melt pool and thermal processes. The inset corresponds to the simulated system.

Figure 2 :

 2 Figure 2: Schematic representation of the model. The laser region is represented in red between 50 nm and 200 nm. The substrate region is represented in blue between 0 and 20 nm.

Figure 3 :

 3 Figure 3: System (model A) with a substrate maintained at T sub = 1000 K. (a) Temperature profiles (solid lines) and the fraction of fcc atoms (dashed lines) at different moments in time. (b) Snapshots of the sytem at 0, 4, and 8 ns. Color coding: fcc atoms are in green, and unknown/liquid atoms in light gray.

Figure 4 :

 4 Figure 4: Position of the solid/liquid interface as a function of time in the polycrystalline system (Model A) with a substrate at 300 K and 1000 K.

Figure 5 :

 5 Figure 5: System (Model B) submitted to a cooling rate α = 375 K/ns with a substrate maintained at T sub = 1000 K. (a) Temperature profiles (solid lines, left axis) and fraction of fcc atoms (dashed lines, right axis) at different times during the global cooling. (b) Snapshots of the sytem at 0.3, 1.6, 2, and 2.4 ns. Color coding: fcc atoms are in green, hcp in red, and unknown/liquid atoms in light gray.

Figure 6 :

 6 Figure 6: Position of the solid/liquid interface as a function of time for a substrate at 1000 K; different cooling rates are shown as solid lines. Local temperature gradients in the solid G s and the liquid G ℓ around the solidfication front at a cooling rate of 375 K/ns are shown as dash-dotted and dashed lines, respectively.

Figure 7 :

 7 Figure 7: Instantaneous velocity as a function of time for different T sub and different cooling rates α.

Figure 9 :

 9 Figure 9: Cluster analysis: number of clusters, mean cluster size, and size of the largest cluster as a function of time. Size is measured in terms of the number of atoms in the cluster (T sub = 1000 K and α = 375 K/ns)

Figure 10 :

 10 Figure 10: System (Model B) submitted to a cooling rate α = 375 K/ns with a substrate maintained at T sub = 1300 K. (a) Temperature profiles (solid lines, left axis) and fraction of fcc atoms (dashed lines, right axis) at different times during global cooling. (b) Enlarged view of temperature profiles (solid lines, left axis) and fraction of fcc atoms (dashed lines, right axis) restricted to the bottom and tip of the protrusion, at 5 ns. The slice over which spatial binning was calculated is shown in the snapshot (c) Snapshots of the system at 0.4, 2.2, 3, and 5 ns. Color coding: fcc atoms are in green, and unknown/liquid atoms in light gray.

Figure 11 :

 11 Figure 11: (a) Position of the finger tip and bottom as a function of time. (b) Theoretical estimate of tip velocity as a function of the characteristic size of the finger (see eq. 9).

Figure 13 :

 13 Figure 13: (a) κ as a function of the simulation box length for the considered temperatures. Symbols are the numerical results. Lines are the corresponding fits. (b) Bulk thermal conductivity results as a function of temperature compared with the values reported by Levchenko et al. [48] and Turlo et al.[START_REF] Turlo | Comparative study of embedded-atom methods applied to the reactivity in the Ni-Al system[END_REF] . These two authors computed κ with the Green-Kubo approach, using Mishin 2004 (M04)[START_REF] Mishin | Atomistic modeling of the γ and γ ′ -phases of the Ni-Al system[END_REF] and Purja Pun and Mishin 2009 (M09)[START_REF] Purja Pun | Development of an interatomic potential for the Ni-Al system[END_REF] potentials for Ni.

Figure 14 :

 14 Figure 14: (a) Snapshot of the solid particle, in green, surrounded by the liquid, in gray. (b) Inverse of the critical radius as a function of the degree of undercooling. Red dots are the simulation results. Black curve is the fit used for the interfacial energy estimate

Figure 15 :

 15 Figure 15: (a) Schematic representation of the system. (b) Solid/liquid interface velocity of the three orientations (100), (110), and (111) as a function of the degree of undercooling.

Table

Table 1 :

 1 Parameters of pure Ni computed with the EAM-09 interatomic potential and compared to ex-

	The huge

perimental values or thermodynamic calculations ([a] estimated in the review of Jiang and Lu

[START_REF] Jiang | Size dependent interface energy and its applications[END_REF]

; [b] Thermo-Calc; [c] Assael et al.

[START_REF] Assael | Reference correlations for the thermal conductivity of liquid copper, gallium, indium, iron, lead, nickel and tin[END_REF]

Table 2 :

 2 morphology t 1 (ns) t 2 (ns) ∆T (K) ∆L z (nm) Summary of MD results for Model B as a function of the temperature of the substrate T sub and the cooling rate α expressed in K/ns. The morphology is named 'C' for columnar grains, 'C+E' for columnar and equiaxed grains, 'C+D' for columnar and dentrites; t 1 is the time at which cooling ended, t 2 is the time at which the first nucleus appeared, ∆T is the degree of undercooling at which first nucleation appeared, ∆L z is the width of the zone with equiaxed grains.

	800	750	C+E	1.27	1.00	551	61.5
		375	C+E	2.53	1.73	527	47.8
		250	C+E	3.8	2.55	525	27.8
	900	375	C+E	2.27	1.60	527	55
	1000	750	C+E	1.00	0.84	545	79.4
		375	C+E	2.00	1.57	528	60.9
		250	C+E	3.00	2.34	530	39.4
		187.5	C+E	4.00	2.99	522	34.4
		150	C+E	5.00	3.72	520	15.4
		75	C	10.00	-	-	-
	1100	375	C+E	1.73	1.50	533	66.1
	1200	375	C+E	1.47	1.46	528	70.9
	1300	750	C+D	0.6	-	-	-
		375	C+D	1.2	-	-	-

  Q (J) 1.63 10 -20 2.41 10 -20 3.65 10 -20

	T (K) v (100) v (110) v (111)
	1100	90.95	59.98	39.53
	1200 1300 1400	75.68 63.07 45.41	52.56 42.55 32.34	35.05 30.20 24.38	K	100 488.1	110 531.2	111 775.1
	1500	26.80	22.51	18.85
	1600	15.53	13.09	12.08
	1700	3.41	2.71	4.89

Preprint submitted to MaterialiaDecember 9, 2022

Acknowledgments

The use of computational facilities at the Computing Center of the University of Bourgogne, DNUM-CCUB, is gratefully acknowledged. Part of this work was also carried out using HCP resources from GENCI-CINES and GENCI-TGCC (Grant 2020 -A0090912032 and Grant 2021 -A0110912032). The authors thank Carmela Chateau-Smith for the careful reading of the manuscript.

Appendix A

In order to estimate the latent heat of fusion, L V , we considered a bulk system of 5.28 nm × 5.28 nm × 5.28 nm at zero pressure, thermalized at T = 100 K. Periodic boundary conditions were imposed in all directions. The sample was then heated with a temperature ramp of 50 K over 50 ps in the NPT ensemble, and then relaxed at constant temperature in the NPT ensemble over 50 ps, and in the NVT ensemble over 10 ps. Quantities were measured in the NVE ensemble over 10 ps. The corresponding average heating was ∼ 450 K/ns. This procedure was repeated until complete melting of Ni at 3000 K. Note that the melting temperature of the bulk is higher than the thermodynamic melting temperature (T m ), because the crystal is perfect, without any free surfaces or imperfections. The system was then cooled with a temperature ramp similar to that imposed for the heating process. A rough estimation of the melting temperature,

is deduced from the temperature for the maximum degree of superheating (T + ) and supercooling (T -) [START_REF] Luo | Solid-Liquid Interfacial Energy and Melting Properties of Nickel under Pressure from Molecular Dynamics[END_REF].

The latent heat of fusion is determined at the melting point of the bulk system (1701 K) as the difference between the enthalpy of the solid phase and the enthalpy of the liquid phase. The latent heat of fusion is estimated at 0.22 eV, equivalent to 21.23 10 6 J/kmol (Fig. 12a). The experimental reported value in the range between 17.03 10 6 J/kmol and 18.90 10 6 J/kmol [START_REF] Thurnay | Thermal properties of transition metals[END_REF] is slightly lower than the value measured in MD. The volume per atom (V at ) is estimated at 11.55 Å 3 (Fig. 12b) at the melting temperature of Ni (1701 K, value of the potential [START_REF] Purja Pun | Development of an interatomic potential for the Ni-Al system[END_REF]). The molar volume is V m = V at × N A = 0.006955 m 3 /kmol. Note that the molar volume estimated with the potential at ambient temperature 0.006612 m 3 /kmol is close to the experimental value 0.006590 m 3 /kmol. Finally, the latent heat of fusion, L V , expressed in appropriate units, is equal to 21.227 10 6 /V m = 3.052 10 9 J/m 3 . This value

Appendix C

In order to estimate solid-liquid interfacial energy γ SL , we built a pseudo 2D system with a solid Ni cylinder surrounded by liquid Ni. Periodic boundary conditions were imposed in all directions. The entire sample was thermalized at a target temperature in the NPT ensemble for 200 ps, followed by the NVT ensemble for 100 ps. The liquid phase around the cylinder was introduced by heating the liquid region in the NPT ensemble for 200 ps and the NVT ensemble for 100 ps to 2500 K (Fig. 14a). The entire sample was then maintained in the NPT ensemble for 5 ns at the target temperature. Several radii were considered: 1.79 nm, 3.58 nm, 7.17 nm, 14.34 nm, and 17.93 nm. Given a particle radius, we chose a large set of target temperatures in order to find its melting point.

According to the classical nucleation theory, since there is no interaction between solid and liquid for the flat face of the cylinder, the Gibbs free energy for a cylinder can be expressed as ∆G = -πr 2 h∆g V + 2πrhγ SL , where h is the cylinder height, r the cylinder radius, and ∆g V is the Gibbs free energy per volume unit. The critical radius r c is obtained when the first derivative of the Gibbs free energy is zero. Moreover, ∆g V can be approximated by ∆g V ∼ L V ∆T /T m . The critial radius r * is then expressed as r * = (γ SL × T m )/(L V × ∆T ), where ∆T = T m -T the degree of undercooling at which the cylinder melts. By considering the values reported in Table 1, the simulation results were fitted by this expression in order to obtain the solid-liquid interfacial energy γ SL = 284 mJ.m -2 (Fig. 14b).

Appendix D

Two-phase systems were considered, with a solid/liquid interface perpendicular to the x direction. The crystallization velocity of Ni was determined for three orientations (100), (110), and (111) of the solid/liquid interface. The following procedure was considered: