Modeling of surface phenomena of liquid Al–Ni alloys using molecular dynamics

This work presents a study on the surface tension of liquid Aluminum–Nickel (Al–Ni) alloys. Obtaining adequate values of surface tension for this system is not a simple task as these alloys present the formation of atomic clusters with short-range order at certain compositions, which dramatically influences surface tension. The Compound Forming Model predicts the influence of these clusters on surface tension, but experimental limitations have obstructed its validation due to deficient thermodynamic data. This work attempts to overcome some of these limitations by using Molecular Dynamics (MD). By comparing the obtained results from MD simulations with those of an equivalent system without clusters, it was possible to infer the role of the atomic clusters on Al–Ni surface tension. It was found that these clusters increase surface tension by decreasing the Al content at the surface. They achieve this reduction in Al content at the surface by trapping Al atoms and hindering their travel to the surface.

www.nature.com/scientificreports/ Another approach to obtain the surface tension of an ideal alloy was proposed by Guggenheim 14 : with the same meaning of γ , γ i , α , c b i , T , and k B as explained above. For more details about the model suggested by Guggenheim, see supplementary material.
Bernard and Lupis 15 used the Sessile Drop method to measure the surface tension of the Gold-Silver (Au-Ag) binary system. Their results showed that the Au-Ag system could be described as an ideal solution. The main difference between ideal and regular solutions is that the interactions between A-A, A-B, and B-B are no longer considered to be equal. Nonetheless, the interactions between dissimilar (i.e., A-B) atoms are considered weak, and thus a random configuration is energetically most favorable 13,16 . The Quasi-Chemical Approximation (QCA) for a regular solution can be used to model the relationship between surface tension and composition 17 . Based on this model, the relationship between surface tension and composition of regular alloys is described as: where p and q are surface coordination fractions, see the Quasi-Chemical Approximation section in the supplementary material for more details about this model. Plevachuk et al. 18 measured surface tension of Bismuth-Tin (Bi-Sn) alloys at 550 K, using the Sessile Drop method. There was a good agreement between their results and the QCA model. For ideal and regular solutions, the surface tension is expected to decrease monotonously by increasing the component's content with the lowest surface tension. Therefore, the composition versus surface tension curve presents a characteristic concave shape.
Real solutions are different from regular solutions as the interaction between dissimilar atoms becomes stronger. Due to this, a random configuration is no longer the most energetically favorable one, and the atoms will tend to form groups 13 . The Compound Forming Model for strong interactions (CFM) detailed in 19 and 20 is then relevant when the tendency is to form A-B groups. Its main characteristic is that it includes the presence of A η B υ complexes (η and υ are the number of A and B atoms), or clusters, in the melt. Evidence of these short medium range order structures has been found using X-Ray diffraction 7,8,21 . For example, Brillo et al. 21 studied the local short-to-intermediate range structure in liquid Al-Ni and Aluminum-Copper (Al-Cu) alloys. Apart from the regular peaks in the structure factor curves, a distinct pre-peak was observed, which indicated the presence of Al-Ni and Al-Cu clusters at certain compositions. Donatella et al. 22 used the large drop method to measure surface tension and density of Al-Ni alloys as a function of composition and temperature. They observed that the behavior of Ni-rich alloys could be well-described by the CFM model. Das et al. 23 , who also performed neutron scattering tests and Molecular Dynamics (MD) studies on Al-Ni melts at 1795 K, confirmed a presence of pre-peaks in the structure factor, indicating the formation of Al-Ni clusters. It can be concluded that the Al-Ni system shows behavior that deviates largely from the ideal solution. The presence of these clusters could completely modify the surface tension of Al-Ni alloys. In the same work, Das et al. 23 studied the influence of clusters on Al-Ni surface tension by comparing their results to the QCA and CFM. Their findings suggested that the presence of clusters increases surface tension by decreasing the surface Al content. This study, however, assumes a priori that (1) the CFM is indeed an appropriate model for Al-Ni, (2) the data from thermodynamic databases are adequate, and (3) the effect of the high reactivity of Al-Ni system is negligible. These assumptions may be questionable, particularly because CFM does not provide surface tension predictions with good accuracy. Additionally, the accuracy and validity of thermodynamic databases for these uses has been questioned before 24 . Finally, even though experimental data was obtained using the Oscillating Drop technique, the high reactivity of Al-Ni can lead to surface contamination such as by oxygen and vaporization 25 . The limitations mentioned are accompanied by natural restrictions in the experimental setup. Due to these restrictions, to the best of the author's knowledge, there are no published studies that provide precise information about, for example, the formation, size, or lifetime of clusters, as well as their influence on surface tension. To overcome these limitations, a simulation method, like Molecular Dynamics (MD), has great potential for predicting physical-chemical properties and understanding complex phenomena. MD has been used extensively to predict the surface tension of pure [26][27][28] and binary metals 29,30 . For example, Kunwar et al. 29 predicted surface tension of pure Al and interfacial energies for Al, Ni, and Al 3 Ni 2 using MD. Calvo 30 analyzed the surface tension of Ag-Au alloys by MD with a tight binding force field. The behavior was close to that of an ideal solution; the surface tension decreases as the content of the element with the lowest surface tension increases. Despite having the same trend as found experimentally, there was a difference of up to 30% between experimental and MD values. Shin et al. 31 studied Aluminum-Iron-Nickel (Al-Fe-Ni) binary alloys using a ReaxFF force field. This study predicted the surface composition of a single alloy. Their findings showed the segregation of Al atoms at the surface, i.e., higher Al content at the surface than the bulk, which is in agreement with the CFM model. Their study does not include any surface tension measurement. As in every other MD study, enough care should be taken to select a reliable force field. As Webb and Grey proposed 32 , a charge gradient correction coupled with the appropriate EAM potential can provide surface tension values that are more accurate. However, this does not exclude the existence of EAM force fields capable of accurately predicting surface events in the absence of charge gradient correction.
In this current study, MD was used, together with the Butler, Guggenheim, QCA, and CFM, to understand how the formation of Al-Ni atomic clusters with short-range order modifies the surface tension of liquid Al-Ni www.nature.com/scientificreports/ alloys. Different systems of Al-xNi (x is a ratio of Ni content in the system) at 2000 K were investigated in terms of surface composition determination and cluster behavior, such as their formation and dissociation. The Embedded Atom Method (EAM) potential developed by Zhou et al. 33,34 is used to describe the multi-particle interactions due to the advantages and improvements with respect to other EAM potentials.

Materials and methods
Force field. Interactions between Al and Ni can be modeled by a few force fields, namely, EAM developed by Baskes et al. 35 and Mishin et al. [36][37][38] , EAM + charge transfer ionic (CTI) developed by Zhou et al. 33,34 , ReaxFF 31,39 , and COMB 40 . The choice of a force field was first based on its reliability in predicting the surface tension of pure Al. It was then calculated using these 7 force fields and compared to experimental results (see Force field section in the supplementary material). The most accurate one was EAM + CTI introduced by Zhou, which was then selected to study the interactions between Al and Ni. Down to a concentration of 30% in Al, this force field also predicted surface tensions of the alloys, which are in agreement with experimental observations as shown hereafter. In our previous publication 26 , we showed that this force field is also the most accurate one in predicting the self-diffusion coefficient.
Simulation setup. Firstly, it is desirable that the obtained results are independent of the model characteristics. In this work, the effect of initial configuration (see Figs. A.1 and A.2 in the supplementary material) and system size (see A.3 in the supplementary material) were tested to evaluate their influence on the surface tension of 50 at% Al-50 at% Ni system. It was concluded that the findings of this work were mostly independent of the prescribed simulation parameters (see the simulation setup in the supplementary material). As a result, the following procedure was used to simulate Al-Ni systems. Firstly, an FCC slab made of 50 400 Al atoms with dimensions 80.5 × 80.5 × 124.0 Å 3 and lattice spacing of 4.05 Å was constructed. Secondly, Al atoms were randomly substituted with Ni atoms until the desired composition of an alloy was reached. The studied compositions were Al (pure Al), Al-0.1Ni (10 at% Al-90 at% Ni), Al-0.2Ni (80 at% Al-20 at% Ni), Al-0.3Ni (70 at% Al-30 at% Ni), Al-0.4Ni (60 at% Al-40 at% Ni), Al-0.5Ni (50 at% Al-50 at% Ni), Al-0.6Ni (40 at% Al-60 at% Ni), and Al-0.7Ni (30 at% Al-70 at% Ni). Notice that the total number of atoms (sum of Al and Ni atoms) was the same for all compositions. Thirdly, an FCC crystal was heated up using NVT ensemble (constant number of elements, volume, and temperature) and Nose/Hoover thermostat from 300 to 2000 K at a heating rate of 0.01 K fs −141 . The dimension of the simulation box in the z-direction was adjusted so that the Al-Ni structures could expand freely when melting. Periodic boundary conditions were applied to the simulation box in all directions. Finally, after 0.6 ns when the systems were at an equilibrium state (Fig. A.4 in the supplementary material), data were collected to measure the target properties. Figure 1 shows the preparation procedure for the Al-0.5Ni system. The Velocity-Verlet algorithm 42 was applied to solve Newton's equation of motion at each time step of 1 fs. The simulations were performed using the LAMMPS software developed by Sandia National Laboratories 43 . The Open Visualization Tool (OVITO) 44  www.nature.com/scientificreports/ Properties measurement. Composition profile. The system was divided into slices of 1 Å in the z direction, and the number of Al and Ni atoms was used to obtain the Al content at each slice. The properties were calculated by averaging them over the last 0.4 ns.

Surface identification.
In an alloy-vapor system, the density ( ρ ) smoothly transitions from the bulk value of the alloy, ρ l , to the bulk value of the vapor, ρ v . Based on 45 , the position and thickness of the surface were determined by fitting the density profile to Eq. (5).
where z 0 is the position of the surface and ω is the surface thickness. For this study, ρ v = 0 as the phase in contact with the alloy is the vacuum. Therefore, the alloy-vapor interface was found at the position where ρ = ρ l /2.
Surface tension. According to Irving and Kirkwood, surface tension (γ ) can be calculated by slicing the system and adding the difference between normal and tangential components of the pressure tensor at each slice 46 . This method is commonly known as "the mechanical approach" and employs Eqs.
where P N (z) and P τ (z) are the normal and tangential components of the pressure tensor, respectively. Pressure includes the kinetic term and virial term. The kinetic term originates from the kinetic energy, whereas the virial term originates from the pairwise forces between atoms. The virial term was computed as described by Thompson et al. 47 . dz is the slice thickness, which in this work corresponds to 0.2 Å. We collected data when stable surface tensions were obtained (see Ideal and non-ideal solutions. Obtaining chemical potentials or similar properties from MD is quite a laborious process. If the system presents strong interactions, as in Al-Ni, it may not be possible to calculate them [56]. As a consequence, the Butler model was used to obtain an activity ratio, defined as: a R Al could be measured using Eq. (1): Besides, the compositions at the surface and the bulk are used to define a concentration ratio according to Eq. (11).
The Al activity and concentration ratios,a R Al and c R Al , were used to infer whether a specific alloy behaved as an ideal or non-ideal solution. An alloy was said to behave similarly to an ideal solution when its a R Al was equal to the corresponding c R Al , and non-ideal otherwise 13 .
Cluster identification. The methodology used to determine the presence of Al-Ni clusters in each alloy was divided into several steps. Firstly, the radial distribution function (RDF) was computed. Then, the Al-Ni bond length was determined from the first peak position of the partial Al-Ni RDF (Fig. 2). In this study, the Al-Ni bond length was determined to be 2.68 Å, consistent with experimental and numerical studies 33, 48-50 . Next, the pairwise distance between all the atoms was calculated, and the pairs with a distance between them corresponding to the Al-Ni bond length were selected. Finally, groups or clusters were identified by picking which of these selected pairs had atoms in common.

Results and discussion
Surface tension. Figure 3 shows the surface tension of Al-Ni systems obtained by MD (blue curve) and equivalent ideal one (red curve) using the Guggenheim equation (Eq. 3). Surface tension obtained by MD was measured by averaging over 400,000 steps (0.4 ns) when the system reached almost equilibrium. Therefore, the fluctuations of pressure in this condition were minor. For instance, the surface tension value for the Al-0.5Ni system was 1.30640 ± 0004 N m −1 when the error bar was derived using standard deviation. Error bars were less www.nature.com/scientificreports/ than 0.05% of the absolute values. Two crucial conditions must be satisfied. Firstly, we should obtain surface tension values comparable with experimental results. Second, the behavior of a real solution should be mimicked since Al-Ni alloys are not ideal or regular solutions. Figure 3 compares MD results (blue curve) with other experimental data (gray curves). The surface tension values are within the experimental range 1, 7-10 , illustrating the high accuracy of the MD results. Not only surface tension of Al-Ni alloys but also surface tension and surface energy of pure Al and Ni using this force field are within experimental range. In our previous publications 26, 29 , we showed that this force field is accurate in prediction of surface tension and surface energy of pure Al and Ni at different temperatures. The predictions obtained from MD are also compared to the ideal and regular behavior (red curve). The latter shows a concave shape, while the MD results show a characteristic non-concave shape, typical of a real solution. These two observations indicate that we can study the surface tension behavior of Al-Ni systems using MD. Data from the MD simulations should provide novel insight into the surface phenomena of these alloys. This is addressed in more detail in the following sections.
Surface identification and composition. By fitting the density profile with Eq. (5), z 0 and ω were obtained for each alloy. The Al content at the surface was then determined using these two quantities with the corresponding composition profile.  www.nature.com/scientificreports/ As seen from Table 1, all alloys featured Al surface segregation, i.e., higher Al content on the surface than in bulk. As there is no available experimental data for the surface composition of liquid Al-Ni alloys, it is reasonable that the comparison between MD and experiments can be replaced by the comparison between MD and the CFM or QCA model. Care should be taken, though, because the surface determined by these thermodynamic models may not accurately reflect the actual surface. While, in contrast, in MD study, the surface composition of Al-Ni alloys was literally directly measured. For these reasons, it was considered best to compare the general trends between c s Al obtained from MD and c s Al from CFM/QCA rather than the numerical values. Surface composition predicted by CFM and QCA were extracted from 25 . Figure 4 shows that the general trend of surface segregation at all compositions obtained from MD matches the one predicted by CFM and QCA. However, higher Al segregation at the surface of the alloys is predicted by both thermodynamic models (CFM/QCA), especially at higher Ni content. As mentioned before, the accuracy and validity of these models are under question due to the inadequate thermodynamic databases. Additionally, CFM does not provide spectacular surface tension predictions, which are inferred from the results of 25 . On the other hand, MD results depend on the force field, and care should be taken when interpreting simulation results. These deficiencies can explain the difference between MD results and the thermodynamic models in Fig. 4. Cluster identification. Cluster identification provided interesting results to unravel the relationship between composition and surface tension of Al-Ni alloys. It was observed that clusters constantly formed and dissociated, with life from 2 to 10 ps, and sizes for all alloys ranged from 6 to 12 atoms, as shown in Fig. 5. An example of a cluster forming and dissociating in Al-0.5Ni is shown in Fig. 6.
Influence of composition on surface tension. As mentioned before, one of the characteristics of real solutions is the possible formation of atomic clusters with short-range order. As observed from Fig. 7, there is an obvious general trend between composition and the number of clusters: the number of clusters decreases linearly as the Al content increases. Quite interestingly, the number of Al-Ni clusters is not the same for all Al-Ni alloy compositions, and, consequently, one could suspect that not all the alloys have the same degree of deviation from the ideal solution. This can be verified by comparing a R Al with c R Al , since their difference is proportional to the magnitude of deviation from an ideal solution (Fig. 8). It should be mentioned that a R Al was calculated using the Butler's equation (Eq. 10). Figure 8 shows that Al-0.2Ni and Al-0.1Ni behavior are close to the like ideal solutions. As such, the activity ratio is approximately equal to its concentration ratio. They were expected to have no or few Al-Ni clusters, consistent with the result shown in Fig. 7. Contrastingly, Al-0.7Ni, Al-0.6Ni, Al-0.5Ni, Al-0.4Ni, and Al-0.3Ni Table 1. Surface composition (Al content at the surface) of Al-Ni alloys. www.nature.com/scientificreports/ showed varying amounts of deviation from the ideal solution. Moreover, the amount of deviation followed the same trend as the number of clusters, i.e., a larger deviation from the ideal state corresponded with a higher number of clusters. The effect of clusters on Al-Ni surface tension was studied by comparing the surface tension obtained with MD with an equivalent ideal one using the Guggenheim equation (Eq. 3), as seen in Fig. 9. It was observed that by increasing the number of clusters, the magnitude of deviation from an ideal solution increases ( Fig. 7 and 9). The presence of clusters decreased Al surface content, i.e., decreased Al surface segregation. This corroborated the observation by Novakovic et al. 25 that clusters increase surface tension. However, in contrast with these authors, the MD clusters were actually characterized and not merely obtained from a thermodynamic model assumed a priori to be correct.   www.nature.com/scientificreports/ As an attempt to understand how clusters influence the surface tension, the surface obtained with MD and an equivalent ideal solution using Eq. (2) were compared in Fig. 10, where the presence of clusters decreased the Al surface content, i.e., decreased Al segregation. This corroborated once again the observation by Novakovic et al. 25 .
The influence of surface Al content on surface tension can be understood as follows: In any A-B binary alloy, one element will always have lower surface tension and will attempt to segregate onto the surface to lower the energy of the system. If, somehow, the segregation of the element with lower surface tension is hindered, then, naturally, the surface tension of the system increases. In an Al-Ni system, Al is the element with the lowest surface tension and thus will attempt to segregate to the surface 24 . The ideal solution model shows the unhindered Al segregation (black curve in Fig. 10) and its corresponding surface tension (black curve in Fig. 9). However, in the Al-Ni system, at compositions with < 70 at.% Al content, Al segregation will be hindered by clusters, thus reducing the surface Al content (blue curve in Fig. 10) and consequently increasing the surface tension (blue curve in Fig. 9). Therefore, it was suggested that clusters increase surface tension by decreasing Al surface segregation.   www.nature.com/scientificreports/ Once it was proposed that clusters alter surface tension by hindering Al surface segregation, the remaining question was, why do clusters affect surface segregation? The suggested answer to this question was based on two observations: the Al atoms do not remain at the surface, and the Al atoms trapped in clusters travel less than their free counterparts (the term trapped atoms makes reference to atoms entrapped by a cluster while the term free atoms denotes atoms not trapped in clusters). These two observations, and their role in explaining how clusters modify the surface composition, are addressed using Al-0.5Ni as an example. To study whether the atoms remained at the surface, a histogram with the number of times that a specific atom appeared on the surface during the last 10,000 time steps is shown in Fig. 11. As can be seen from this figure, the specific atoms on the surface change with time. This is caused by the continuous motion of atoms across the slabs.
The characteristics of the movement, however, were not the same for all atoms. Figure 12 compares the trajectory during the last 10,000 time steps of one of the trapped Al atoms with the trajectory of three random free Al atoms in the same time interval. It can be seen that trapped Al atom traveled less distance than their free counterparts. This observation was generalized to the whole system by comparing the mean square travel distance (MSTD) of the free and trapped Al atoms during the same last 400,000 time steps. This value of both groups was given by Eqs. (12) and (13), respectively, where (t) is the position of the atom at time step t, and (t 0 ) is the reference position. It should be reminded that the clusters were forming and dissociating, thus the Al atoms' status as free or trapped atoms varied during the studied time length. Therefore, it was decided to use the position at the previous time step as the reference position to avoid confusing free and trapped behavior. Figure 11. Histogram of Al atoms on the Al-0.5Ni surface. www.nature.com/scientificreports/ As can be seen from Fig. 13 trapped Al atoms traveled less than their free counterparts. Notice the MSTD shown in this figure does not increase with time because the reference position was the previous atom position.

Conclusion
Binary Al-Ni alloys have been, and continue to be, the focus of extensive research. In particular, the surface tension of liquid Al-Ni alloys is of quite an interest since this property influences the final microstructure as well as the presence of defects. Despite this interest, the relationship between composition and surface tension of liquid Al-Ni alloys is not yet fully understood. This is due to the fact that, at some compositions, Al and Ni strongly interact, forming Al-Ni atomic clusters with short-range order. Then, the thermodynamic models normally applied to many alloys to predict surface tension do not work so well. This difficulty is further compounded by experimental complications and limitations, which cause that some key properties useful to understand the Al-Ni surface tension behavior can only come from thermodynamic models.
This work used Molecular Dynamics in an effort to understand the complex relationship between composition and surface tension of liquid Al-Ni alloys. This simulating technique has the great advantage that it allows for direct prediction of atomistic behavior.
By comparing the results from the simulations with those of an equivalent ideal solution, it was possible to infer the role of Al-Ni clusters playing in surface tension. It was found clusters appear in alloys that show a large deviation from the ideal solution, a larger deviation from the ideal state corresponded with a higher number of clusters.
It was proposed that these Al-Ni clusters increase surface tension by decreasing the surface Al content. They achieve this by entrapping Al atoms and hindering their travel, which in turn decreases their probability of reaching the surface. Indeed, there is an observable relationship between the relative amount of Al atoms trapped in clusters and the difference in Al content between the surface alloy and its corresponding ideal alloy. This entrapment of Al atoms leads to an increase in surface tension because this property is related to the amount of Al at the surface. Since Al has a lower surface tension than Ni, it tends to segregate onto the surface to lower the overall energy of the system. Therefore, fewer Al atoms at the surface conduces to an increase in surface tension compared to the ideal solution. As not all Al-Ni alloys present the same amount of deviation from the ideal state, the influence of clusters on surface tension will not be the same. This leads to the non-monotonous, non-concave curve that is characteristic of this system.
In summary, although surface tension only depends on what happens at the surface, it may seem that in Al-Ni alloys, the bulk indirectly plays a role in surface tension. This is because the bulk composition will influence cluster formation, which will modify the surface composition, which in turn, will modify surface tension.
In the future, this method could be applied to other binary systems in order to investigate their behavior and compare it with ideal, regular, and real solutions.

Data availability
The datasets generated and/or analysed during the current study are not publicly available as the data also forms part of an ongoing study but are available from the corresponding author on reasonable request.