Structural and atomic displacement evaluations of Aluminium nanoparticle in thermal annealing treatment: an insight through molecular dynamic simulations

In the present study, Molecular Dynamic (MD) simulations were applied to uncover the phase transition of Aluminium (Al) nanoparticles in annealing treatment. Parameters and properties, like potential energy, Mean Square Displacement (MSD), diffusion coefficient, crystal style and Radial Distribution Function (RDF), were discussed in the following article. In the linear heating stage (stage one), it was confirmed that the melting behaviour started at 960 K and ended at 1038 K. With the process of melting, the diffusion of Al atoms was accelerated from shell to core. In the thermostat stage (stage two), lots of bulges appeared due to the high kinetic energy and limited atomic interactions. Those bulges spread just like waves. In the slowly cooling stage (stage three), we observed the supercooling phenomenon. The solidification of melted Al nanoparticle started at 540 K, which was much lower than its melting temperature. Comparing with the original one, the cross-section of annealed nanoparticle showed some penetrating defects.


Introduction
As a typical energetic carrier, Al nanoparticles have played a critical role in applications of combustion and flame, such as nanofluids, solid propellants, thermites and gelled propellants [1]. The outstanding performance of Al powders is because of its high oxidations enthalpy, relatively low cost, high combustion temperature and environmentally benign products [2]. Therefore, researches of Al nanoparticles have attracted lots of concern. Dreizin et al investigated the mechanism of oxidation behaviours for Al nanoparticles by experimental approaches [3], while van Duin et al did similar work through reactive MD simulations [4]. Meanwhile, many other investigations were focused on functional coatings for Al nanoparticles, including carbon [5,6]; ether [7][8][9]; Nickel [10][11][12][13] and some other organic materials [14].
Thermal annealing treatment is a quite common technique for solid materials. In this field, plenty of researches have been published in recent years as well. Tianlin Huang et al characterized the annealing twins of high purity Al by electron backscatter diffusion experiments [15]. Zhengyi Jiang et al reported how rolling parameters and annealing conditions affected the microstructure and texture of 1060 Al alloy [16]. Hui Zhang et al tested the addition of Silicon in annealed AlSi20 alloy and 8009 Al alloy [17]. Verheyden et al annealed monocrystalline or cast Al microwires under protective environment, and also obtained their mechanical properties [18]. Xuming Pang et al were focused on the Al coating on P355NL1 steel and their annealing behaviours. Then they reported the microstructure of Al coating at different temperatures [19]. Nesma T. Aboulkhair et al investigated the laser-assisted melting of Al for 3D printing with an unprecedented extent of freedom [20].
Despite experimental studies, MD simulations were widely applied to uncover solid-liquid phase transitions as well. By monitoring spatiotemporal arrangements, Tran Phuoc Duy simulated the melting of bcc Fe liquidlike atoms, and discussed the presence of crystalline model [21]. Ojwang et al developed a reactive force field to predict the melting, solidification and crystal lattice of Al nanoparticles [22]. Yasushi Shibuta et al discussed the effect of particle size and temperature on the melting process of iron nanoparticle in the melted iron environment [23].
Before this study, Junpeng Liu et al have already published the melting behaviour of Al nanoparticle [24]. However, the entire annealing treatment not only includes the melting stage, but also the following cooling stage. Besides, the comparison of annealed and unannealed Al nanoparticle needs to be discussed as well. Therefore, we designed this study in order to present the entire annealing treatment. Due to the continuity of the following projects, MD simulations were applied.

MD simulations and EAM potential
With little time cost and powerful data export, MD simulations have been widely acknowledged as an alternative method to investigate the static and dynamic properties of microscopic structures and phases [25]. This type of numerical simulation is generally conducted on the solution of Newton's classical motion equations. Therefore, information like atomic forces, acceleration, motion, trajectory could be dealt as a function of time. Meanwhile, relative chemical and physical properties will be presented to predict performances of the investigated system. In this investigation, equilibrated MD simulations were performed under Large Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [26].
Because our MD simulations were going to be performed with plenty of Al atoms, interactions between those atoms ought to be predefined. The physical reliability of MD simulations' result is largely determined by those potentials. Although the interatomic potential is defined to describe the potential energy of the simulated system, not all chemical or physical phenomenon will be presented in the simulation. Therefore, the selection of interatomic potential is better to fit the demand of this study. As a widely discussed atomic potential for metal and alloy models, an Embedded Atom Method (EAM) potential was applied in this investigation. EAM potential is a typical empirical forcefield with efficient computing capability [27]. Which means that our model could be built closer to real scale and contain more atoms, then those obtained results would be more convincing. As å å In which, E: total energy of atom; F: embedding energy; ρ: atomic electron density; Φ: pair potential interaction; α and β: element types; i and j: atomic ID. The EAM force field applied in this study was developed by Mendelev et al [28]. According to their report, the potential could describe the atomic diffraction behaviours of both Al and Fe elements. As a well-performed potential file, it would provide solid-liquid phases description not only for this study, but also for our following Fe-Al bio-element investigations. Although only an Al nanoparticle is annealed this time, it's better to apply the same EAM potential between two projects for comparison.

Simulating model and parameters
In order to study the effect of annealing treatment on Al nanoparticles accurately, we constructed an Al nanoparticle with a radius of 8 nm. There were about 129 299 Al atoms in the nanoparticle. The nanoparticle was modelled by creating a spherical region and filling it up. Therefore, the original particle could be considered as ideal spheres without any defect. Then, the particle was placed at the centre of a simulation box with periodical boundaries. The box was about 24×24×24 nm 3 , while the vertical and horizontal minimum vacuum was about 4 nm long. The creation of LAMMPS datafile was accomplished by Visual Molecular Dynamics (VMD) [29]. The simulation box and the Al nanoparticle itself are presented as figure 1 below.
Once the model had been prepared, a pre-equilibrating process was performed for about 0.5 ns with a temperature of 300 K. This pre-treatment was designed due to the consideration of system stability and optimization. Then, temperature adjustment treatments were designed as three stages. Stage one is a linear heating process for targeted Al nanoparticle. The linear heating simulation started at 300 K, and ended at 1500 K with the heating rate about 1.2 K ps −1 . The peak temperature had already been confirmed higher than melting points of both studied Al nanoparticles and bulk Al material [24]. In stage two, the melted nanoparticle was fixed at the peak temperature for about 0.5 ns to further discuss its liquid states. Finally, stage three is a linear cooling process in order to make liquid Al nanoparticle solidify and return to room temperature. Similarly, the cooling treatment was controlled linearly as well, and the rate was about 0.8 K ps −1 .
Temperature control is one of the most critical points in this study. Comparing with other relevant theories, Nose [30]/Hoover [31] thermostat is believed to describe the canonical ensemble (NVT) with true dynamics [4]. Therefore, all simulations in this study applied Nose/Hoover thermostat. Additionally, the timestep of 1 fs has been assigned to all simulations. The preparation and export of snapshots were both accomplished by OVITO program [32].

Stage one: heating
The entire annealing simulation started with a linear heating process. In order to study continuous changes of Al nanoparticle in stage one, the deviation of temperature is roughly divided into four equal sub-stages. Figure 2 is the plot of temperature for stage one. Despite the initial and final temperature, two more temperature levels are labelled in figure 2. The deviation between each label is 400 K. Figure 3 is the snapshot of Al nanoparticle at these four sub-stages. It could be clearly observed that the most explicit change for Al nanoparticle happened in the range from 700 to 1100 K. Before 700 K, the crystal configuration of Al nanoparticle was not changed too much, and the initial shape of Al nanoparticle was still kept. On contrast, configurations obtained at 1100 and 1500 K were highly ununiform and irregular. Therefore, the melting point of this simulated Al nanoparticle could be roughly deduced as higher than 700 K and lower than 1100 K.
In order to further uncover the change of Al nanoparticle from 700 to 1100 K, a plot of MSD result is shown in figure 4. The MSD result generally demonstrates the extent of atomic displacement at each temperature.  There are two critical points in figure 4, one is point A (960 K) and the other is point B (1038 K). Point A states that the rate of average atomic displacement was accelerating. It means that the crystal restriction for Al atoms was collapsing, and the movement of each atom would be larger and larger. While at point B, the melting process ended and the whole system would be considered as a liquid phase. Furthermore, the melting point of Al nanoparticle is validated by a plot of potential energies. As figure 5 shows, a transition of total potential energy started when the temperature was in the range from 960 to 1038 K. Because of the failure of crystal limitation, atomic distances were enlarging in this period. Consequently, the potential energy of the whole system jumped into a higher level in such short term. Therefore, it is reasonable to confirm that the melting behaviour of this simulated Al nanoparticle started at 960 K, and then ended at 1038 K.  The plot of figure 6 is the statistic of crystal style for Al nanoparticle. Those results were obtained through the modification of 'Common Neighbor Analysis' in OVITO [32]. Because the edge Al crystal was incomplete itself, those surface Al atoms were classified as other style crystals. That's why the initial content of Face Centre Cubic (FCC) style was about 90.3%. As figure 6 shows, the percentage of FCC style declined with climbing temperature. Before 960 K, the trend of decline was presented like the parabolic curve. When the melting behaviour began, the decline curve was accelerated and reached the zero point at 1038 K. Obviously, the zero point presented the fully melting temperature as well. Additionally, the snapshot in figure 6 shows the distribution of FCC crystal in the cross-section of Al nanoparticle at 1000 K. Combining with the decline curve, the configuration was under rapid melting stage. As the snapshot shows, FCC style crystal was mainly distributed around the core of Al nanoparticle. Meanwhile, the layer of other crystal style was thicker than before. Therefore, it could be initially deduced that our simulated nanoparticle was melted from shell to core, in other words, this was an inward melting phenomenon.
Due to the limitation of crystal style classification, the migration of liquid phase boundary ought to be validated by atomic potential energy configurations. Figure 7 is the cross-sections of Al nanoparticle at 960 K; 992 K; 1022 K and 1038 K respectively. According to snapshot (a) of figure 7, the nanoparticle was actually covered by a layer of Al atoms with low potential energy, and this was also where the melting behaviour started. Then, the layer of high potential energy was thicker and thicker associated with climbing temperatures. It could be observed that the crystal in high potential energy zone was dislocated, so the inward melting migration could  also be considered as inward migration of crystal dislocation. Finally, as presented in the snapshot (d), the distribution of potential energy was uniform again at 1038 K.
Assisted by 'Coordination Analysis' function of OVITO [32], the lattice of Al nanoparticle was studied by RDF. Those RDF results are plotted as figure 8 below. At 300 K, the peak of RDF plot is narrow, high and regular, which demonstrates that the lattice of Al crystal was well and periodically distributed in the nanoparticle. Next at 960 K, this was the temperature that melting behaviour just started. So, the lattice of Al crystal was still kept, and periodical peaks of RDF plot could be observed as well. However, those peaks were broader and lower than those at room temperature. It implies that the probability of detecting an Al atom at a certain distance was much lower than before, although the lattice hadn't been collapsed. From 300 K to 960 K, those Al atoms were migrated out of their original position, but the extent of migration was not enough to destroy the periodical distribution of crystal lattice. Finally, at 1038 K and 1500 K, peaks of RDF plot were not detectable anymore, because the lattice of Al crystal had been destroyed.
Furthermore, the Al crystal could also be represented by two-dimensional distributions of Al atom. Figure 9 is the quantity of Al atoms while scanning alongside the Y-axis. Due to the presence of regular and perfect lattice, there were intervals through the nanoparticle. In figure 9, those intervals are clearly presented on the plot of 300 K. Meanwhile, the content of Al atom would be extremely high when the position Y is on the array of lattice. Then the plot of 1500 K does not show intervals and high peaks. Despite the absence of lattice, those results also imply that there was no vacuum in the melted nanoparticle. Besides, the mass centre of melted Al nanoparticle didn't change explicitly while comparing with solid nanoparticle. According to figure 9 and export of LAMMPS computation, the mass centre of Al nanoparticle at 1500 K was still around (120, 120, 120) Å. Additionally, another difference was on the edge of Al nanoparticle. It is obtained that the edge of melted Al nanoparticle was enlarged for about 6 Å. Figure 10 is the statistic of distances from edge atoms to mass centre. In which, the vertical dotted line represents the judgment of whether the nanoparticle was fully melted or not. It seems that the increase of boundary before 1038 K was regular and stable, while the boundary of fully melted Al nanoparticle was highly unstable. This instability is presented by the fluctuation of boundaries. Just like figure 3(d) shows, there were several bulges on nanoparticle surface which were highly unstable.

Stage two: thermostat
Stage two was performed under the maximum temperature of entire simulations. This section mainly describes the state of Al nanoparticle in liquid phase. Figure 11 has recorded both MSD and cohesive energy results from 0 to 200 ps. According to the red plot of figure 11, the linear climbing slop demonstrates the steady movement of liquid Al atoms as expected. Based on those MSD data, the diffusion coefficient of system could be further calculated [7]. It is obtained that the diffusion coefficient of liquid Al nanoparticle at 1500 K was about 1.364 Å 2 ps −1 . The cohesive atomic energy is another essential property to describe condensed materials, as it represents the extent of interatomic forces. According to our results, the cohesive energy of melted Al nanoparticle at 1500 K was about −2.993 Kcal mol −1 . The distribution of atomic forces on each atom is shown in figure 12. As this melted Al nanoparticle was physically isotropic in this study, we only select forces on X-axis to represent the entire distributions. It seems that the positive and negative distribution of atomic forces was symmetric and uniform, and the vast majority of Al atoms were within ±1 eV Å −1 . In this isolated computing system, there was no force-rich zone from core to shell, so the stress concentration effect would not appear for melted Al nanoparticles. It could also be deduced that the irregular surface melted Al nanoparticle was not caused by stress concentration on bulges. In order to further illustrate the deformation of melted Al nanoparticle, the cross-section of its ZY plane has been obtained by scattering points of atoms. As figure 13 shows, the radius of Al nanoparticle had been relatively enlarged, which is in agreement with figure 10. It could be measured that the largest bulges on surface were as high as 11.859 Å. However, it is observed that those bulges were not fixed at several certain sites but spread from one place to another. This is because Al atoms were still     under gravity and interatomic interactions, although they had been free from the crystal. On this occasion, energy of one bulge would be delivered to another bulge, just like waves in nature.
Finally, atomic X-velocity components while scanning alongside the Y-axis have been plotted in figure 14. The figure is comprised of scattering points and each point characterises a single atom on Y-axis. It is shown that the majority of atoms were moved at the speed of ±18 Å ps −1 . In this velocity interval, the quantity of Al atoms was so high that the interval is plotted like a rectangle. Meanwhile, there were still several points that were higher or lower than the main interval. The maximum velocity was about −24.46 Å ps −1 . Generally, the distribution of atomic velocity was not various with radial distance.

Stage three: cooling
Once the Al nanoparticle had been fully melted, the next stage was to cool the nanoparticle until room temperature, so that the microstructure could be equilibrated, and lattice would be reformed. In industry, the process of cooling is various. While in this study, we only applied the linear cooling treatment. The cooling rate has been mentioned above, which is a little lower than the heating rate. Figure 15 is the plot of potential energy as a function of declining temperature. Comparing with figure 5, figure 15 has shown an opposite developing trend. Due to the weakness of atomic thermal motion, atomic distances for the cooling system was decreasing. Therefore, the potential energy of the overall system dropped off as a function of time. The most critical point was about the solidification point of Al nanoparticle. In figure 15, the gradient of plot changes in the range from  540 K to 469 K. Which means that the solidification point of this simulated nanoparticle was in this temperature interval. Obviously, this was much lower than the melting points for bulk Al material and simulated nanoparticle. This phenomenon is so-called as supercooling. The supercooling always appears on pure liquid substances. Our simulation was designed and performed in an ideal environment. There was no impurity and containers in this system. Consequently, the lack of crystallization nucleus caused the supercooling phenomenon.
In terms of the cooling process, the most critical part ought to be the growth of crystal. Therefore, the interval of 540 and 469 K could be further discussed. Figure 16 shows four cross-sections of Al nanoparticle which have been coloured by potential energy. Those snapshots were obtained at 540; 517; 492 and 469 K respectively. In contrast with melting process, the solidification process was not developed from shell to core but spread from single point to the whole plane. In figure 16(a), it seems that the initial crystal appears at the top right corner. To some extent, the growth of crystal was started at the surface as well. Then, figures 16(b), (c) has illustrated the direction of entire development and final configuration. Due to the defect of obtained crystal, there was an area of high potential energy throughout the nanoparticle.
When the annealed Al nanoparticle was obtained, it would be worthwhile to uncover what has been changed due to the annealing treatment. Figure 17 is prepared in order to compare crystal style distributions between original and annealed Al nanoparticles. Figure 17 (a) is the cross-section of original model. Its crystal was wellorganized, as it was whittled from bulk material. Before annealing treatment, the percentage of FCC crystal was about 90.3%. While in figure 17(b), the annealed one, this value was declined until 74.4%. Meanwhile, the percentage of other style was raised from 9.7% to 20.0%, and Hexagonal Close Packed (HCP) structure appeared after annealing. Figure 18 is the distribution of Al atoms, while scanning alongside the Y-axis. It could be confirmed that the mass centre of annealed nanoparticle was not changed. Due to the presence of lattice, the figure shows gaps between peaks. However, those gaps were not as deep as in figure 9. This is because of the irregular placement of lattice and the presence of crystal defects. While combining figures 17 and 18, it seems that the surface of annealed Al nanoparticle was not as uniform as before. But the size of bulges had been much smaller than melted one.

Conclusion
This paper reports the microstructure and atomic arrangement of Al nanoparticle in thermal annealing treatment. Those obtained parameters were exported by LAMMPS based MD simulations. In the first heating stage, the nanoparticle was heated from 300 to 1500 K. We observed an inward melting behaviour started at 960 K. Associated with the happening of melting, the diffusion of Al atoms was accelerated, and the total potential energy jumped to a higher level. While at 1038 K, there was no crystal detected from then on. In the second thermostat stage, the edge of melted Al nanoparticle was highly unstable. Due to the high kinetic energy and limited atomic interactions, many bulges were observed on nanoparticle surface. There was no stress concentration found in the melted nanoparticle. In the third cooling stage, the supercooling phenomenon appeared due to the absence of crystallization nucleus in ideal simulation conditions. Therefore, the solidification temperature was confirmed at 540 K. When the temperature was reduced until room temperature, several defects of crystal were found throughout the Al nanoparticle.

Acknowledgment
This study is financially supported by Fundamental Research Funds for Central Universities (HEUCFG201815)