Bayesian Data Assimilation of Temperature Dependence of Solid–Liquid Interfacial Properties of Nickel

Temperature dependence of solid–liquid interfacial properties during crystal growth in nickel was investigated by ensemble Kalman filter (EnKF)-based data assimilation, in which the phase-field simulation was combined with atomic configurations of molecular dynamics (MD) simulation. Negative temperature dependence was found in the solid–liquid interfacial energy, the kinetic coefficient, and their anisotropy parameters from simultaneous estimation of four parameters. On the other hand, it is difficult to obtain a concrete value for the anisotropy parameter of solid–liquid interfacial energy since this factor is less influential for the MD simulation of crystal growth at high undercooling temperatures. The present study is significant in shedding light on the high potential of Bayesian data assimilation as a novel methodology of parameter estimation of practical materials an out of equilibrium condition.

However, interfacial parameters estimated from MD simulations are largely limited to those of static interfaces in equilibrium conditions. For example, the capillary fluctuation method is mostly applicable at melting point [18,19,28], where the solid-liquid interface does not propagate. Moreover, temperature dependence of interfacial properties is mostly deduced in the CNT-based method under the assumption of temperature dependence of thermodynamic parameters [27,28]. Therefore, it is not straightforward to discuss the change in solid-liquid interfacial energy out of the equilibrium condition.

Phase-Field Model
In the phase-field model, morphology of microstructure is defined by spatial distribution of the phase-field variable, φ. Here, φ = 1 and φ = −1 represent solids and liquids, respectively, and the solid-liquid interface is defined by the continuous change between −1 and 1. In this study, a phase-field model for the isothermal solidification of pure metal [55] is employed following our previous study [46], which is described as follows: Here, λ = a 1 W 0 /d 0 is the coupling constant with a 1 = 5 √ 2/8, W 0 is the interface thickness, and d 0 is the capillary length represented by d 0 = σ 0 (T m c p /∆H 2 ). β 0 is the kinetic coefficient. u int is the dimensionless undercooling at the interface defined as where ∆T is the undercooling temperature, ∆H is the latent heat, and c p is the specific heat. τ (n) and W (n) are relaxation time and interface width, respectively. a c (n) and a k (n) are functions based on the interface orientation, respectively, and are expressed as follows.
ε c and ε k represent the strength of anisotropy of interfacial energy and the kinetic coefficient.
(n x , n y ) is the normal vector at the interface obtained by

State Vector and System Model
Data assimilation is a method to combine a simulation model with observation data to estimate states and/or parameters of the system [47]. In the simulation model, solution of the governing equation is deterministically derived once initial and boundary conditions are given. However, it is, in general, impossible for the simulation to reproduce phenomena in nature precisely since there is a discrepancy between simulation results and true dynamics due to the incompleteness of the simulation model, fluctuations in the phenomenon, unknown boundary conditions, and so on. Therefore, uncertainty of time evolution of states is considered by introducing state variables with a probability distribution, which is called a system model. That is the basic concept of the data assimilation.
In the simulation, time evolution of a system is calculated by discretizing the governing equation in time and space in general. Physical quantities are represented by values at the representative points (mostly grid points). A column vector consisting of all representative values i,j on the grid point (i, j) at time t is called state vector x t : Each element of the state vector is called a state variable. A system model, which derives the state vector x t at time t from the state vector x t-1 at time t − 1, is defined as where f t is the simulation model at time t and v t is the system noise representing the imperfectness of simulation model. The observation model, which represents the relationship between observation data y t and the state vector x t at time t, is defined as Observation noise w t includes measurement error and imperfectness of the simulation model.

Ensemble Kalman Filter (EnKF)
In data assimilation, state variables follow a probability density function. State variables are iteratively updated by integrating observation data into existing probability density functions based on the Bayesian inference approach [47]. This procedure is called filtering. There are various types of filtering approaches, including Kalman filter [56], ensemble Kalman filter (EnKF) [57,58], and particle filter [59]. In the present study, EnKF was employed. In the EnKF, many simulations with different parameters are performed in parallel, which represents the probability distribution function. Figure 1 shows schematic image of the calculation procedure of EnKF, which consists of cycles of prediction and filtering steps. Multiple simulations (called ensembles) are independently performed in the prediction step to obtain predictive state variables and parameters, which are then corrected based on the observation data. Mathematical expression of the EnKF is summarized in Appendix A. Expected values of parameters are obtained from the ensemble average at each filtering step as sequential data with respect to time. The estimated value of each parameter is then obtained by a time average of the expectation value in the time period after the convergence.

Ensemble Kalman filter (EnKF)
In data assimilation, state variables follow a probability density function. State variables are iteratively updated by integrating observation data into existing probability density functions based on the Bayesian inference approach [47]. This procedure is called filtering. There are various types of filtering approaches, including Kalman filter [56], ensemble Kalman filter (EnKF) [57,58], and particle filter [59]. In the present study, EnKF was employed. In the EnKF, many simulations with different parameters are performed in parallel, which represents the probability distribution function. Figure 1 shows schematic image of the calculation procedure of EnKF, which consists of cycles of prediction and filtering steps. Multiple simulations (called ensembles) are independently performed in the prediction step to obtain predictive state variables and parameters, which are then corrected based on the observation data. Mathematical expression of the EnKF is summarized in Appendix A. Expected values of parameters are obtained from the ensemble average at each filtering step as sequential data with respect to time. The estimated value of each parameter is then obtained by a time average of the expectation value in the time period after the convergence.

Calculation Procedure of Data Assimilation
In this study, kinetic coefficient β0, interface energy σ0, and their anisotropy parameters εk and εc during the growth of a single crystal under isothermal conditions were estimated by the EnKF procedure. The time evolution equation of the phase-field method (Equation (1)) was used as ft in the system model (Equation (A1)), which describes the time evolution of state vectors from t − 1 to t. Equation (1) was discretized in a standard finite different scheme with second order accuracy in the space and it was solved in an explicit Euler scheme. The calculation system was divided into 100 × 100 grid points. State variable vector xt and observation vector yt are given as:

Calculation Procedure of Data Assimilation
In this study, kinetic coefficient β 0 , interface energy σ 0 , and their anisotropy parameters ε k and ε c during the growth of a single crystal under isothermal conditions were estimated by the EnKF procedure. The time evolution equation of the phase-field method (Equation (1)) was used as f t in the system model (Equation (A1)), which describes the time evolution of state vectors from t − 1 to t. Equation (1) was discretized in a standard finite different scheme with second order accuracy in the space and it was solved in an explicit Euler scheme. The calculation system was divided into 100 × 100 grid points. State variable vector x t and observation vector y t are given as: Nanomaterials 2021, 11, 2308 5 of 14 A total of 10,004 variables, including phase-field variables φ i,j at all lattice points (i, j) and four parameters to be estimated (β 0 , σ 0 , ε k , ε c ), were used as state variables in state variable vectors. The phase-field variables at each grid point in the observation data were used as the observation variables in the observation vector. Tables 1 and 2 show the parameters used for the phase-field simulation and for the EnKF data assimilation, respectively. We created 100 phase-field simulations using independent state vectors and optimized the simulations based on the observed data by alternately executing prediction by the system model and filtering by Equation (A6). Equation (A1) was used as a system model, where Equation (1) was employed as a nonlinear operator f t , and Equation (A2) was used as an observation model. The observation noise w t was set as a random number vector generated from Gaussian distribution according to covariance matrix R φ which is an identity matrix in the shape of 10,000 × 10,000. Observation matrix H t in the observation model is given as: where the dimension of the matrix is 10,000 × 10,004. In this study, atomic configurations from MD simulations of the growth of a single crystal of nickel were employed as observation data, the preparation of which is described in the next section.  Table 2. Parameters for EnKF data assimilation.

Parameter Symbol Value
Ensemble number

Molecular Dynamics Simulation for Observation Data
Prior to the data assimilation by EnKF, observation data were prepared by MD simulation. MD simulation was performed using large-scale atomic/molecular massively parallel simulator (LAMMPS) [61]. An embedded atom method (EAM) potential fitted by Purja Pun and Mishin [62,63] was employed to obtain the interatomic potential for nickel. Representative properties are listed in Supplementary Materials. In the present study, nickel was employed as a representative face-centered cubic (fcc) metal. It is known that interfacial anisotropy of fcc metals is generally larger than that of body-centered cubic (bcc) metals [64]. Velocity-Verlet method was used to integrate the classical equation of motion with a time step of 1.0 fs. Nose-Hoover thermostat and barostat [65,66] were employed to control temperature and pressure. Note that there exists a large temperature distribution in the calculation system during crystal growth due to the release of latent heat with a conventional thermostat [46,67], even though average temperature of the system is steadily controlled at the target temperature. Therefore, Langevin thermostat [68] was also employed to keep the temperature uniform over all areas of the system, which is essential for a precise comparison with the isothermal phase-field model. Figure 2a shows the initial configuration of the calculation system for MD simulation, which was prepared as follows. The liquid structure was prepared by heating an fcc crystal of nickel consisting of 260 × 260 × 10 unit cells (2,704,000 atoms) at 2000 K with the canonical (i.e., the number of atoms, volume, and temperature constant) ensemble for 10 ps. Separately, a solid nucleus was prepared as an octagonal cutout from the fcc crystal with four {100} and four {110} facets. The solid nucleus was then inserted into the liquid structure while omitting all liquid atoms located within 2.5 Å from a solid atom to avoid unexpected proximity between liquid and solid atoms at the interface. Energy minimization was performed for the combined structure. The prepared structure was then relaxed at 1455, 1480, 1505, and 1530 K for 800 ps with the isobaric-isothermal (i.e., the number of atoms, pressure, and temperature constant) ensemble and growth behavior of the solid nucleus was investigated. Periodic boundary condition was employed in all directions. Note that the melting point of Ni of this EAM potential was approximately 1680 K [60], which was estimated by the convergence temperature technique [69,70]. That is, the temperature 1455, 1480, 1505, and 1530 K corresponds to the undercooling temperature ∆T = 225, 200, 175, and 150 K. Atomic configurations from MD simulation were analyzed by polyhedral template matching (PTM) [71] implemented in the Open Visualization Tool (OVITO) [72] to identify local atomic structures (i.e., solid or liquid). Purja Pun and Mishin [62,63] was employed to obtain the interatomic potential for nickel. Representative properties are listed in Supplementary Materials. In the present study, nickel was employed as a representative face-centered cubic (fcc) metal. It is known that interfacial anisotropy of fcc metals is generally larger than that of body-centered cubic (bcc) metals [64]. Velocity-Verlet method was used to integrate the classical equation of motion with a time step of 1.0 fs. Nose-Hoover thermostat and barostat [65,66] were employed to control temperature and pressure. Note that there exists a large temperature distribution in the calculation system during crystal growth due to the release of latent heat with a conventional thermostat [46,67], even though average temperature of the system is steadily controlled at the target temperature. Therefore, Langevin thermostat [68] was also employed to keep the temperature uniform over all areas of the system, which is essential for a precise comparison with the isothermal phase-field model. Figure 2a shows the initial configuration of the calculation system for MD simulation, which was prepared as follows. The liquid structure was prepared by heating an fcc crystal of nickel consisting of 260 × 260 × 10 unit cells (2,704,000 atoms) at 2000 K with the canonical (i.e., the number of atoms, volume, and temperature constant) ensemble for 10 ps. Separately, a solid nucleus was prepared as an octagonal cutout from the fcc crystal with four {100} and four {110} facets. The solid nucleus was then inserted into the liquid structure while omitting all liquid atoms located within 2.5 Å from a solid atom to avoid unexpected proximity between liquid and solid atoms at the interface. Energy minimization was performed for the combined structure. The prepared structure was then relaxed at 1455, 1480, 1505, and 1530 K for 800 ps with the isobaric-isothermal (i.e., the number of atoms, pressure, and temperature constant) ensemble and growth behavior of the solid nucleus was investigated. Periodic boundary condition was employed in all directions. Note that the melting point of Ni of this EAM potential was approximately 1680 K [60], which was estimated by the convergence temperature technique [69,70]. That is, the temperature 1455, 1480, 1505, and 1530 K corresponds to the undercooling temperature ΔT = 225, 200, 175, and 150 K. Atomic configurations from MD simulation were analyzed by polyhedral template matching (PTM) [71] implemented in the Open Visualization Tool (OVITO) [72] to identify local atomic structures (i.e., solid or liquid).   Figure 2b shows snapshots of atomic configuration during growth of the solid nucleus from the undercooled melt of ∆T = 200 K. The nucleus grew preferentially in <100> directions and it became a rhombic-like structure, which means that fourfold symmetry appeared. The obtained atomic configuration was not suitable for the observation data for data assimilation with phase-field simulation. Therefore, atomic configuration from the MD simulation was converted into a phase-field profile in line with our previous studies [73,74]. A cross-section (90 × 90 nm 2 ) of the MD simulation cell was divided into two-dimensional difference grids of 100 × 100. After assigning all atoms in the closest grid, the majority of local atom configurations (i.e., solid or liquid) for assigned atoms were employed as the phase-field variables of each grid point (solid: 1, liquid −1). Since this voxel structure had no interfacial thickness, it was relaxed by solving the phase-field equation without the curvature effect [73] to obtain the phase-field profile with diffuse interface. This conversion procedure was carried out for the time series of atomic configuration of MD simulation with 10 ps interval. Obtained phase-field profiles were employed as the observation data for data assimilation. The observation data between 300 and 600 ps at 10 ps interval were used in the data assimilation in following the sections to avoid the initial relaxation period of solid nucleus in the MD simulation.

Results and Discussion
Now, parameter estimation of β 0 , σ 0 , ε k , and ε c from the dataset of MD simulation at 1480 K was performed. Figure 3a shows time change of the estimated values of four parameters. Estimated values of three parameters, β 0 , ε k , and σ 0 converged to certain values with decreasing the variance. On the other hand, variance of the estimated value of ε c did not decrease during the estimation although the estimated value itself came close to a certain value. It was expected that accuracy of estimation of ε c was lower than those of the other parameters. Figure 3b shows snapshots of observation data from the MD simulation and representative results of the estimated structure. The crystal shape in the observation data and that of the representative ensemble member were in good agreement. The same procedure was performed for other datasets of MD simulations at 1455, 1480, and 1505 K. These results are summarized in the Supplementary Materials. In general, results for the other temperatures agreed with that of 1480 K. That is, estimated values of three parameters, β 0 , ε k , and σ 0 converged to certain values with decreasing variance, while accuracy of estimation of ε c was again low compared to the other parameters. Figure 4 shows the temperature dependence of estimated values of four parameters from observed data of 1455, 1480, 1505, and 1530 K. Estimated values of the last filtering step are plotted in the figures. Estimated values of β 0 , ε k , and σ 0 decreased with increasing temperature. The negative temperature dependence of σ 0 agreed with our previous estimation of the solid-liquid interfacial energy of bcc-Fe by EnKF [46] and some reports in the literature [42][43][44]. Bayesian inference derived the most probable values of solid-liquid interfacial energy at various temperatures from the results of the MD simulation without any prior knowledge. It is guaranteed that the phase-field model employed in this study reproduces the Gibbs-Thomson effect properly. Therefore, the parameters derived in this study were appropriately within the range where the Gibbs-Thomson effect is valid. One possible reason of discrepancy from some studies of positive temperature dependence may be due to the effect of interface curvature [37,41]. However, it was difficult to find the physical origin of negative temperature dependence directly from our result. The degree of temperature dependence of β 0 was smaller than that of σ 0 , Incidentally, β 0 was nearly independent of the temperature within the examined temperature range in our previous study for bcc-Fe [46]. Moreover, it was difficult to find a clear trend in the temperature dependence of ε c since the accuracy of estimation of ε c was lower, as described above.
The estimated values of σ 0 ranged from 0.27 to 0.38 J/m 2 , which basically overlapped experimental and theoretical values of σ 0 for Ni at melting point (0.255 [13], 0.284 [75], 0.306 [76], and 0.325 J/m 2 [14]). Regarding the kinetic coefficient, β 0 took values between 0.0035 and 0.0037 s/m. β 0 can be converted into the interfacial mobility µ by the following relation, µ = c p /β 0 ∆H [55]. Using the values of c p and ∆H in Table 1, β 0 = 0.0035 m/s was converted into µ = 0.418 m/sK. This is within the range of reported values, 0.18-0.45 m/sK, which were derived from MD simulations with planar solid-liquid interfaces of Ni [77]. It is convincing that both the solid-liquid interfacial energy and interfacial mobility estimated in this study were consistent with previous reported values from various methodologies. On the other hand, temperature dependence of ε k took the opposite trend to our previous estimation of bcc-Fe, which was the positive temperature dependence. This difference might come from the difference in the strength of anisotropy. That is, a strong anisotropy appeared in the crystal structure of fcc-Ni in this study, whereas a weak anisotropy appeared in that of bcc-Fe structure [46]. Further study is needed to discuss the anisotropy in interfacial mobility.  other hand, temperature dependence of εk took the opposite trend to our previous estimation of bcc-Fe, which was the positive temperature dependence. This difference might come from the difference in the strength of anisotropy. That is, a strong anisotropy appeared in the crystal structure of fcc-Ni in this study, whereas a weak anisotropy appeared in that of bcc-Fe structure [46]. Further study is needed to discuss the anisotropy in interfacial mobility.  As described above, the accuracy of estimation of ε c was low compared to the other parameters in the four-parameter estimation. Therefore, ε c was separately estimated while fixing the other parameters at estimated values of β 0 = 0.00275 [s/m], ε k = 0.338, and σ 0 = 0.277 [J/m 2 ]. For this one parameter estimation, state variable vector and the observation matrix were modified as follows.
Observation vector y t was the same as Equation (12). The observation matrix is given as: where the dimension of the matrix is 10,000 × 10,001. Figure 5a shows the time change of the estimated values of ε c starting from different initial distributions at 1505 K. The other conditions were the same as those of the four-parameter estimation. Two estimations did not converge to the same value. That is, it was not successful in obtaining the converged value of the estimation for ε c , even from the procedure of one parameter estimation. Furthermore, the effect of the anisotropy parameter on growth morphology of crystal structure was examined by a phase-field simulation with fixed parameters. Three parameters were fixed to the values estimated above at 1455 K and three values of ε c = 0.008, 0.010, and 0.012 were employed. The other parameters were the same as listed in Table 1. Figure 5b shows the phase-field profile after 30,000 step simulations. The morphologies of the structures did not change significantly with respect to ε c . Therefore, it is considered that the effect of anisotropy of solid-liquid interfacial energy on the crystal structure is very small under the condition of large undercooling temperature. In that condition, growth velocity is very fast and the anisotropy in interfacial mobility is dominant. In other words, it is difficult to estimate the parameters of a less influential factor in the framework of the present study. The anisotropy parameter ε c may be estimated when a near equilibrium structure of the crystal is employed as observation data, which will be investigated in the next step.
Observation vector yt was the same as Equation (12). The observation matrix is given as: where the dimension of the matrix is 10,000 × 10,001. Figure 5a shows the time change of the estimated values of εc starting from different initial distributions at 1505 K. The other conditions were the same as those of the four-parameter estimation. Two estimations did not converge to the same value. That is, it was not successful in obtaining the converged value of the estimation for εc, even from the procedure of one parameter estimation. Furthermore, the effect of the anisotropy parameter on growth morphology of crystal structure was examined by a phase-field simulation with fixed parameters. Three parameters were fixed to the values estimated above at 1455 K and three values of εc = 0.008, 0.010, and 0.012 were employed. The other parameters were the same as listed in Table 1. Figure 5b shows the phase-field profile after 30,000 step simulations. The morphologies of the structures did not change significantly with respect to εc. Therefore, it is considered that the effect of anisotropy of solid-liquid interfacial energy on the crystal structure is very small under the condition of large undercooling temperature. In that condition, growth velocity is very fast and the anisotropy in interfacial mobility is dominant. In other words, it is difficult to estimate the parameters of a less influential factor in the framework of the present study. The anisotropy parameter εc may be estimated when a near equilibrium structure of the crystal is employed as observation data, which will be investigated in the next step.

Conclusions
In the present study, temperature dependence of solid-liquid interfacial properties during the crystal growth was investigated by data assimilation with EnKF. It is the advantage of the methodology in this study that both the interfacial energy and mobility out of equilibrium condition could be estimated simultaneously from Bayesian inference. Negative temperature dependence was found in the solid-liquid interfacial energy, kinetic coefficient, and anisotropy of kinetic coefficient from simultaneous estimations of four parameters. However, the anisotropy parameter of the solid-liquid interfacial energy

Conclusions
In the present study, temperature dependence of solid-liquid interfacial properties during the crystal growth was investigated by data assimilation with EnKF. It is the advantage of the methodology in this study that both the interfacial energy and mobility out of equilibrium condition could be estimated simultaneously from Bayesian inference. Negative temperature dependence was found in the solid-liquid interfacial energy, kinetic coefficient, and anisotropy of kinetic coefficient from simultaneous estimations of four parameters. However, the anisotropy parameter of the solid-liquid interfacial energy did not converge during the four-parameter estimation and it did not converge even in the subsequent single parameter estimation. Since the anisotropy parameter of the solid-liquid interfacial energy did not affect the morphology of the crystal in the phase-field simulation with fixed parameters, it is difficult to estimate the parameter for the less influential factors in the observation phenomena. In other words, we can find less influential factors in the observation data without any prior knowledge for target phenomena. In summary, it is significant that this study showed the high potential of data assimilation as a methodology of parameter estimation in the out of equilibrium condition. The overlap between atomistic and continuum simulations, which was achieved owing to recent progress in highperformance computing, creates new research concepts and fields. We call this cross-scale modeling, as an evolution from conventional multi-scale modeling [30].
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/nano11092308/s1, Table S1. Representative properties of the EAM potential for Ni employed in this study. Figure S1: Estimation of four parameters (kinetic coefficient β 0 , interfacial energy σ 0 , and their anisotropy parameters, ε k and ε c ) using observation data of molecular dynamics (MD) simulation at 1455 K. Figure S2: Estimation of four parameters (kinetic coefficient β 0 , interfacial energy σ 0 , and their anisotropy parameters, ε k and ε c ) using observation data of molecular dynamics (MD) simulation at 1505 K. Figure S3: Estimation of four parameters (kinetic coefficient β 0 , interfacial energy σ 0 , and their anisotropy parameters, ε k and ε c ) using observation data of molecular dynamics (MD) simulation at 1530 K.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
The authors acknowledge Takuya Fujinaga for the kind support in performing molecular dynamics simulation in Section 4.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Mathematical Expression of Ensemble Kalman Filter
In the EnKF, system and observation models are defined as w t is given as Gaussian distribution with zero mean and a covariance matrix of R t . In the prediction step, time evolution of states vector is calculated as where superscript (i) represents ith state in N ensembles. Subscript a|b represents the state vector of time a, which is filtered at time b. Probability distribution function after the prediction step is given as where δ is the Dirac delta function, and y 1:t−1 represents observational data from time t = 1 to time t − 1. In the filtering step, states vector is updated as [44,52]: K t is ensemble approximation of Kalman gain, given as V t|t−1 and R t are sample covariance matrixes of state vector and observation error, given as: Probability distribution function after the filtering step is given as