Sintering of Alumina Nanoparticles: Comparison of Interatomic Potentials, Molecular Dynamics Simulations, and Data Analysis

Sintering of alumina nanoparticles is of interest both from the view of fundamental research as well as for industrial applications. Atomistic simulations are tailor-made for understanding and predicting the time- and temperature-dependent sintering behaviour. However, the quality and predictability of such analysis is strongly dependent on the performance of the underlying interatomic potentials. In this work, we investigate and benchmark four empirical interatomic potentials and discuss the resulting properties and drawbacks based on experimental and density functional theory data from the literature. The potentials, which have different origins and formulations, are then used in molecular dynamics simulations to perform a systematic study of the sintering process. To analyse the results, we develop a number of tailored data analysis approaches that are able to characterise and quantify the sintering process. Subsequently, the disparities in the sintering behaviour predicted by the potentials are critically discussed. Finally, we conclude by providing explanations for the differences in performance of the potentials, together with recommendations for molecular dynamics sintering simulations of alumina.


Introduction
Aluminium oxide (Al 2 O 3 ) or alumina has evoked significant research interest due to its use in technological and industrial applications, ranging from the field of aerospace [1], water treatment, filtration and cooling systems [2][3][4], medical and health care [5][6][7], to materials engineering [8][9][10]. Alumina coating on aluminium is known to improve wear resistance [11] and act as electrical and thermal insulation [12]. The material exhibits unique properties, such as e.g., high thermal conductivity, high hardness, and a good resistance to corrosion and abrasion [13].
An important technique for processing alumina is sintering, which involves consolidation and densification of powder compacts, and is driven by a reduction in total interface energy [14]. As one of the technologies that has been known and used for a long time, the process has proven to be important and beneficial for the production of a variety of ceramic structures, from powder-metallurgy parts to bulk components. Recently, sintering has evoked much interest in modern industrial production technologies such as 3D printing and additive manufacturing.
Accurate modelling of the sintering process is on the wish list of many sintering industries, mainly due to the engineering expediency it provides. However, the complex interplay of a multitude of parameters together with a lack of fundamental materials understanding at the nanoscale have proven to be a major hindrance to this end [15].
Indeed, many of the benefits and applications of alumina have their origins at the nanoscale making the material system and the sintering process ideally suited for investigation with atomistic modelling and simulations. Such simulations, which follow the trajectory of individual atoms, have now established themselves as an indispensable tool in advancing our understanding of materials behavior [16][17][18]. Sintering and concomitant processes have been investigated through molecular dynamics (MD) simulations for a large variety of materials (including alumina) and for a number of different techniques such as selective laser sintering or pressure-assisted/pressure-less sintering [19][20][21].
Atomistic modelling of materials behavior is, however, not straightforward and hinges primarily on three aspects [18]: i) accuracy of the initial structure used, ii) reflection of real-world boundary conditions, and most importantly, iii) the reliability, robustness and predictive capability of the interatomic potential used. Such interatomic potentials are particularly critical for oxide materials such as alumina due to the ionic (and partly covalent) nature of the oxygen-metal bonds, which require the inclusion of long-range effects. Although many interatomic potentials for alumina exist, each potential is developed with a different rationale in terms of applicability and transferability. The choice of the interatomic potential hence becomes crucial for the atomistic study, and is usually made by evaluating the predictive capability of the potential on certain basic material properties -e.g., lattice constants, cohesive and surface energies and elastic constants -and subsequently applied to the more complex problem at hand.
Aluminium oxide often exists in crystalline form -the corundum structure -and is a member of the hexagonal close packed (hcp) family. A single molecule consists of two aluminium and three oxygen atoms, the interactions among which are usually modelled by incorporating two-body and/or three-body interaction terms. The Coulomb-Buckingham [22] or Born-Mayer-Huggins type potentials [23] are based purely on two-body interactions, whilst the potential by Vashishta et al. [24] incorporates both two-body and three-body interactions, with the latter associated with triplets containing Al-O bonds. A common feature of such potentials is the usage of fixed charges -both nominal and effectivefor Al and O ions. To account for local heterogeneous electrostatic environment, particularly in the presence of interfaces, free surfaces and segregated ions, improvements have been proposed by incorporating variable point charges. Reactive potentials, such as reactive empirical bond order (REBO), reactive force field (ReaxFF) and charge optimised bond order (COMB) potentials, integrate many-body effects and two-body interactions through a bond order term [25]. Other approaches such as the charge-transfer ionic [26] and the second-moment tight-binding [27] models are simpler in their approach and allow for charge equilibration during the process under study.
In this work, we investigate the sintering of spherical nanoparticles using molecular dynamics simulations. Specifically, we compare four different interatomic potentialsthe Vashishta potential [24], the Coulomb-Buckingham potential of the Matsui type [28], the Born-Mayer-Huggins potential parametrised by Bouhadja et al. [23] and the charge transfer ionic potential [26]. We first present the mathematical formulations of the potentials and evaluate the predictive capability of these potentials towards basic material properties. Subsequently, the sintering process is simulated for three different temperatures by considering spherical nanoparticles of three different sizes. Sintering is modelled as isothermal, high temperature process. In other words, we neglect pressure effects that play an important role in low temperature sintering [29][30][31]. We develop a data analysis strategy that is based on a number of atomic and geometric properties to quantify and characterize the process. As a result, we find that despite exhibiting similar material properties, the investigated potentials result in very different sintering outcomes. We then discuss the characteristics used in terms of their ability to clearly quantify the progress of sintering by setting numerical thresholds.
In the following section, we start by introducing the used potentials and rewrite them in a consistent mathematical formulation. Subsequently, details of modelling, simulation, and data analysis for characterising sintering are provided. We then compare various alumina properties (such as elastic properties, surface energies, lattice parameters, etc.) that result from our atomistic simulations with data from experiments and density functional theory calculations. The actual sintering simulations are analysed in detail for all four potentials, which is concluded by a discussion of their respective suitability and appropriateness for MD simulations of nanoparticle sintering of alumina.

Interatomic potentials of alumina
In the following, four interatomic potentials are considered for the sintering study: the Vashishta potential (in the following abbreviated as Vash), the Coulomb-Buckingham potential (CB), Born-Mayer-Huggins potential (BMH), and the Charge Transfer Ionic + EAM potential (CTIE). Whilst the last three are based on two-body interactions, the first one, i.e. Vash is based on both two-body and three-body interactions. Three potentials (Vash, BMH, and CB), have the same fixed charges for the anions and cations, and one (CTIE), allows for variable charges. In the following, all the potentials used in our study are briefly introduced.
where the two-body part (U (2) ij (r ij )) of the effective potential is written as Here, H ij is the strength of the steric repulsion, Z i the effective charge (in units of the electronic charge |e|), D ij is the strength of charge-dipole attraction, W ij is the van der Waals interaction strength, η ij the exponents of the steric repulsion term, r ij = |r i − r j | the distance between the i th and j th atoms, and λ and ζ are the screening lengths for Coulomb and charge-dipole terms, respectively. The two-body interaction potential includes steric size effects of the ions, charge-transfer effects leading to Coulomb interactions, charge-dipole interactions due to the electronic polarisability of ions and induced dipole-dipole van der Waals interactions. It has a cut off at r c (r c = 6Å for alumina) beyond which the tow-body term is zero. The three-body term is as follows: where In Eq. 3, 4 and 5, B jik is the strength of the interaction, θ jik the angle formed by r ij and r ik , C jik andθ jik are constants, and Θ(r 0 − r ij ) is a step function. The three-body term becomes zero when θ jik =θ jik .
2.2. Two-body interaction potentials 2.2.1. The Coulomb-Buckingham potential The (CB) potential consists of of long range (Coulomb) and short range (Buckingham) potential terms as shown in Eq 6.
where Z i and Z j are the effective charges, A ij and ρ ij are the parameters of repulsion, r ij ‡is the distance between the i th and j th atoms, C ij is the Van der Waals constant. The exponential term in the CB potential provides a better description of strong repulsion due to the overlap of the closed shell electron clouds, which is vital in a simulation of bombardment by energetic atoms or ions, etc. When atoms lie in a crystalline environment, the distance between two atoms is small, and the short-range potential plays an important role. However, when atoms are far from the substrate (free atoms), the larger distances between atoms leads to a rapid reduction of the short-range contribution and the long-range potential has the dominant effect on the free atoms. This potential can be categorised as [22] (i) partial-charge (Al +1.4175 , O −0.945 ) Buckingham-type Matsui [28] (ii) full-charge (Al +3 , O −2 ) Buckingham-type Bacorisen potential [32] (iii) full-charge (Al +3 , O −2 ) Buckingham-type Sun potential [33] (iv) BKS (Beest-Krammer-Santen) potential; which is another CB type potential [34].
In the current study we use as CB potential the Matsui type.

The
Born-Mayer-Huggins or Tosi/Fumi potential The BMH potential consists of long range (Coulomb) and a short range interactions terms, where σ is an interaction-dependent length parameter, ρ is an ionic-pair dependent length parameter. The first term represents the long-range Coulomb interaction with charges Z i and Z j between i th and j th ions, separated by r ij ‡ The last three terms on the right hand side represent the Born repulsive, van der Waals and dipole dispersion interactions, respectively. In the current study, we use this potential parameterised by [23]. ‡ Implementation of CB and BMH potentials in LAMMPS requires cut-off distances for non-coulombic and coulombic interactions, and 10 and 15Å are used as cut-off values for those, respectively.

The Charge Transfer Ionic + EAM potential
The CTIE potential consists of contributions of the non-ionic interaction and those of the ionic interaction and charge transfer [26,35], where U CTI is the CTIE potential for ionic interaction and charge transfer, and U EAM is the non-ionic interaction. The first part becomes zero for non-electrostatic material. The electrostatic part can be expressed as: with, The pair potential is written as where φ ij (r ij ) is the pair energy between the i th and j th atom separated by r ij , F i (ρ i ) is the embedding energy to embed an atom i in a local site with electron charge density ρ i . More details are given in [26,35].

2.2.4.
A first juxtaposition of the four potentials From the above equations we observe the following: The formulations of CB and BMH are mathematically similar and become identical if D is zero in Eq. 7. They both are two-body interaction potentials. However, in the two-body interaction of the Vash potential, there are two additional terms for steric repulsion and charge-dipole interaction. The additional term for three-body interaction incorporated in Vash contains spatial and angular dependent factors, which is useful for amorphous alumina [24]. The three-body term applies to triplets Al-O-Al and O-Al-O with fixed angles. The charges are fixed and the same for anions and cations in all these potential. By contrast, CTIE is a potential of variable charges, where charges of two anions or cations are not identical. Numerical values of the parameters used in the four potentials are tabulated in the Appendix.

Methods
All simulations are performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) atomistic simulation software [36]. Molecular statics (MS) is employed to calculate and verify basic potential properties of individual potentials, and molecular dynamics is employed to determine the melting temperature and to simulate sintering. A stable time step of 1 fs is used in all MD simulations.

Material properties at 0 K
Five different properties of the interatomic potentials, viz., lattice constants, cohesive energy, vacancy formation energies, elastic constants and surface energies, are calculated using standard procedures via MS simulations. Further details are provided in thesupplementary material.

Melting temperature
The melting temperature of alumina for the different potentials is determined using standard molecular dynamics simulations via the following procedure. A 3D triclinic sample of size 33 × 29 × 39Å and an xy tilt factor of 16.7Å consisting of 1764 Al and 2746 O atoms and periodic boundary conditions in all directions is thermalised at several targeted temperatures using the isobaric isothermal (NPT) ensemble. The average atomic volume, determined via a Voronoi tessellation as implemented in the open visualisation tool Ovito, is plotted for the different target temperatures. The melting point for the system and potential is then approximated as the temperature at which a jump in the average atomic volume is observed. For more details, the reader is referred to the supplementary material. in x, y and z directions of the simulation box, respectively, is created using the atomic simulation environment tool [37]. Spherical particles of different radii are cut out of the bulk system using Atomsk [38]. The simulation box is kept periodic in all three directions to keep the de-clustered atoms in the box. The minimum distance between the sphere surface and the ends of the simulation box is kept to be at least 25Å to avoid the influence of periodic images. The energy is minimised using the conjugate gradient and FIRE [39] algorithms to maximum force-norm value of ≈ 10 −5 eV/Å. The particle is initialised with a Gaussian velocity distribution corresponding to 60% of the melting temperature, followed by a microcanonical ensemble for 30 ps to stabilise the temperature of the system. Once Figure 1: The simulation setup: A particle is cut out from bulk material (shown as semitransparent), where D x = D y = D z = 25Å is the smallest distance between particle and the surfaces. Al and O atoms are colored with blue and yellow, respectively. The particle is heated up to the target temperature T f and thermalised at T f for 10 ps each. It is then replicated and equilibrated at T f . The particles are then brought together at a distance within the potential cut off and at T f . Subsequently, the system is maintained at the T f for 1 ns during which sintering can occur under favorable conditions. the system is equilibrated at 0.3T m , it is subsequently heated up to the targeted sintering temperature in an canonical (NVT) ensemble for 10 ps. The NP is then thermalised at the target temperature for 10 ps. Subsequently, a replica of the particle is created, generating another particle at the same temperature and size, keeping a minimum gap of 50Å between them. The entire system is thermalised again at the targeted temperature by an NVT for 10 ps. The particles are subsequently brought close to each other so that the distance between them is less than the cut-off distance of those interatomic potentials that have such a cut-off. The temperature of the system is then kept at the targeted temperature in an NVT ensemble for 1 ns. Under favourable conditions, inter-particle diffusion takes place resulting in sintering of the two NPs at constant temperature.

Characterisation of sintering
To compare the performance of the four potentials in terms of sintering it is required to characterise and to quantify the time and temperature dependent sintering process. As a comparison of sintered particles is complex (e.g., in terms of geometry), we use altogether six different quantities to characterise and quantify sintering, with their evolution providing an indication on the amount of sintering.

Shrinkage ratio
The coalescence of two particles can be quantified by the shrinkage ratio S f which is defined as: where d 1 and d 2 are particle sizes (diameters), and d COM denotes the distance between the centres of mass of them. The latter is obtained as geometric mean of the atomic positions belonging to the particle under consideration. S f is initially negative when NPs are detached, and zero when they are brought into contact. If sintering takes place, the particles move closer to each other resulting in an increase in S f . When the two spheres are completely sintered, the centre-to-centre distance is ideally zero, resulting in S f being one.

Surface area
We note that during sintering, a reduction in surface area is to be expected: this reduction is driven by the movement of surface atoms which are more mobile than the atoms closer to the core of the particle. The transport of material results in reduction of surface area and thus a reduction of overall energy [40,41]. To compute the surface area, we use the construct surface-mesh modifier in Ovito [42]. The algorithm constructs a surface mesh with Delaunay tetrahedralisation of the atomic points in conjunction with an appropriate probe sphere radius, where elements are removed that do not fit in the probe sphere. The remaining triangular faces of the tetrahedra defines a 2D boundary or surface of the defined domain. For more details, the reader is referred to [43]. The virtual probe spheres fill in the empty space between atomic points without spanning the atomic points. A smaller value gives more details of mesh, and a larger one smooths out the mesh. Three values for the probe radius (2, 4 and 6Å) were tested. No difference was observed in the results for probe radii of 4 and 6Å; hence a radius of 6Å was chosen for the final analysis. The resulting surface area is then normalised by the initial surface area to facilitate comparison of different particle sizes and potentials.

Mean square displacement (MSD)
The MSD defines the average displacement of all atoms in a system with respect to their original positions. It provides a good approximation of diffusion taking place during sintering at different temperatures and is defined as where r i (t) and r i (0) are the current and initial positions of atom i, respectively.

Neck curvature
The transport of surface atoms to reduce the overall surface area results in a reduction of curvature at the neck, i.e., the region where the two particles meet. The neck curvature hence provides a good estimate of the amount of sintering. In the case of completely sintered particles, the neck is completely eliminated, resulting in a value of 0 for the curvature. We compute the curvature as follows: 3D atomic coordinates are projected onto a plane whose normal is perpendicular to the centre-to-centre axis of the two particles. The boundaries of this region are obtained from spline fitting of the projected data, resulting in two curves that define the profile of the neck. Subsequently, the local minimum/maximum of the two curves are identified, to which second order parabolas (y = ax 2 + bx + c ) are fitted and from which the curvature can be computed. This is repeated for rotations of 15 • of the sample about the centre-to-centre axis of the two NPs, resulting in a total of 24 radii of curvatures for a given simulation snapshot. The mean radius of curvatureR c is calculated as the arithmetic average of all values from which the mean neck curvature κ is obtained as κ = 1/R c . Fig. 2 provides a schematic diagram of the procedure.
3.4.5. Fraction of ions at the neck During the sintering process, the Al and O ions change their positions in the contact region, and this material transport reduces the curvature. The change in the geometry of the neck can then also be quantified by the fraction of ions forming the neck region. This is done via the following procedure: The atoms in the initial configuration are assigned an index (0 or 1) based on the NP they belong to. For each snapshot in time, the average index of the nearest neighbour atoms is computed, and those atoms with an averaged index between 0.9 and 1.1 are defined to constitute the neck region.
3.4.6. Norm of stretch tensor A disadvantage of MSD is that it may lead to erroneous results under the presence of rigid body motion, as will be evident by the results presented below. An alternative measure can be constructed by using displacements whilst accounting for the immediate neighbourhood of individual atoms: we first construct the deformation gradient tensor using the displacement vectors of individual atoms. Using the polar decomposition F i = U i R i , we can decompose F i into the right stretch tensor, U i , and the rotation tensor, R i . Exploiting the orthogonality of R i , we obtain U i as: and subsequently, the rotation tensor as The stretch tensor of individual atoms is subsequently averaged over the nearest neighbors as follows:Ū where V i is the atomic volume obtained via a Voronoi construction. For further analysis, we consider the square Euclidean norm of the averaged stretch tensor Ū i 2 as the characteristic that quantifies sintering.

Potential and material properties
Material properties such as lattice constants, cohesive energies, vacancy formation energies, elastic constants, surface energies and melting temperatures were calculated as described in Section 3 and the appendix, respectively. These properties, calculated with all four potentials along with the values from various experiments and density functional theory calculations are provided in Tab. 1. A graphical comparison of the range of properties is shown in Fig. 3.
In the following, we treat values from experiments and DFT equivalently and refer to them collectively as experiments.    For most properties, the values predicted by almost all interatomic potentials are well within or close to the range of values obtained from experiments (cf. Tab. 1). The notable exception is seemingly the BMH potential which predicts values for 0 K properties that are significantly higher than those obtained from experiments. Exceptions are also seen in the case of vacancy formation energy of Al; while Vash and CTIE predict values towards the higher values observed in experiments, values from CB and BMH are significantly beyond the range of experimental values. The latter two also predict significantly higher (numerically, lower) cohesive energy values than that observed from experiments. In case of the surface energies, we observe values that are by and large close to the range of data obtained from experiments. Only the CTIE is an exception here: energies for all the surfaces are significantly higher than that seen in experiments.
The predicted melting temperature values obtained via MD simulations with the four potentials show that only the CB potential predicts a value well within the range of experimental data. The values from all other potentials are significantly higher, with the highest value of 4200 K by the CTIE potential. To facilitate objective comparison of the sintering process with the different potentials, we use the homologous temperature of the corresponding potential.
In Fig. 3 we show the 25/75 percentile, median and total range of material parameters measured or calculated by various techniques. Outliers are detected based on 150% of the inter-quartile range. The figure shows that some of the obtained material properties exhibit significant scatter together with a strong skewness of the distribution. For example, the highest melting temperature (4200 K) is almost 210% more than the lowest melting temperature

Sintering
The results of the sintering simulations for each investigated temperature and particle size are shown as supplementary movies M1 to M3, with each movie comparing the prediction of all four potentials. Fig. 4 and 5 show the evolution of all investigated characteristics over time for all four potentials (indicated by the different line colours), for three different particle sizes (indicated by markers) and for different temperatures (sub-plots in columns).
Shrinkage ratio: Fig. 4a shows the evolution of the shrinkage ratio (S f ). For all four potentials, we generally observe an increase of S f with time, particularly at higher temperatures, implying that the two particles move closer to each other. However, the magnitude of this movement differs strongly from one potential to another, for all particle sizes and temperatures. At T= 0.6 T m , for R = 4 nm and 6 nm, S f remains negative or close to zero throughout the simulation, indicating that the two particles barely come into contact with each other. At this temperature, a non-negligible value of S f is only observed for particles with R = 2 nm and potentials CTIE, BMH and Vash; while the former two display a steady increase in S f , Vash shows an almost instantaneous increase to ≈ 5% at the start which is then maintained throughout the simulation.
The evolution of S f at T= 0.7 T m follows similar trends, however, for NPs with R = 2 nm, a non-negligible positive S f is seen for all four potentials. While a steady increase over the entire simulation time and a value greater than 40% is observed for BMH, CB and CTIE, Vash displays an increase only within the first 0.2 ns and a saturation to significantly lower value of roughly 10%. Additionally, CTIE shows a S f value of roughly 50% and 20% for R = 4 nm and 6 nm, respectively; for these particle radii, other potentials do not show any noticeable value of S f .
Significantly different trends are seen at the higher temperature of T= 0.8 T m . Both CB and BMH show very similar characteristics in the evolution of S f over time. For the three particle sizes, S f evolves towards higher values than at lower temperatures and attains values of 0.7, 0.4 and 0.2 after 1 ns, evidencing a trend of decreasing S f with increasing particle size.
In general, S f with CTIE is higher than that predicted by other potentials. Very contrasting trends are seen for the Vash potential: Whist the evolution for R = 2 nm follows that of CTIE, for R = 4 nm, S f saturates to a value of 0.1 after 0.2 ns and remains negative for R = 6 nm.
Surface area: The evolution of overall surface area (A f ) of the NPs is shown in Fig. 4b. To enable easy comparison, the values are normalised w.r.t. the corresponding initial area before the start of the sintering process. Therefore, values should start at 1 and are expected to reduce during the simulation time if sintering proceeds favourably.
By and large, evolution of A f displays trends similar to that observed for S f : when the two NPs come close to each other, a decreasing surface area over time is generally observed (compare at T= 0.6 T m : BMH, Vash and CTIE potentials and at R = 2 nm and T= 0.7 T m : BMH, CB and CTIE). For cases where no sintering it expected (evidenced by zero shrinkage ratio), A f remains close to unity.

Mean square displacement:
The evolution of MSD over time is shown in Fig. 4c. The curves clearly show an increase in MSD over time, albeit with different amounts for different sample sizes and potentials. With increasing temperature, a significantly higher MSD is observed, indicating the propensity for increased diffusion at higher temperatures. However, the trends do not complement those observed with S f and A f . For example, the MSD for  CTIE at T= 0.7 T m is greater than the MSD for all other potentials and particle sizes, which contradicts the trends observed with S f . This can be observed at T= 0.6 T m and T= 0.7 T m as well. Furthermore, the evolution of MSD with the CB potential at T= 0.6 T m shows nonmonotonic evolution with periodic crests and troughs, indicating artifacts associated with the superposition of rigid body motion.

Neck curvature:
The change in the neck curvature κ with time is shown in Fig. 5a and, as expected, decreases with time. A value of zero indicates a complete removal of the neck region. The initial value of κ for larger particles are higher than that of the smaller particles, and at the same time, the initial change in the curvature for R = 2 nm is faster than that for R = 4 nm and R = 6 nm for all temperatures. Generally, no unexpected trends are observed and the evolution is consistent with the trends observed for S f .

Fraction of ions at the neck:
The change in the accumulated fraction of ions at the neck, N i f , is shown in Fig. 5b. The observed trends here align well with those seen in the shrinkage ratio. With increasing temperature, N i f reaches values as high as 0.9, indicating that almost the entire system of two particles is now classified as the neck region. Such a scenario occurs when the two particles have completely sintered and remixing of atoms from the two particles has taken place. At lower temperatures N i f reaches substantially lower values; the only exception seemingly is for R = 2 nm with the CTIE potential at T= 0.7 T m , where a value of more than 0.9 is reached.
Norm of stretch tensor: Fig. 5c shows that the evolution of |U | 2 starts at a value of 3 since the initial stretch tensor is an identity tensor. The square of the norm is used here to facilitate comparison with MSD. Here, it is evident that there are differences between the evolution of |U | 2 and MSD: at T= 0.8 T m , the MSD after 1 ns for Vash is significantly lower than other potentials for R = 4 nm and 6 nm, whereas the |U | 2 values are very similar to those from BMH and CB for R = 6 nm. Furthermore, the non-monotonic evolution with periodic crests and troughs seen with the CB potential at T= 0.6 T m vanishes in the evolution of |U | 2 . However, for the Vash potential, the evolution of |U | 2 shows fluctuations from the mean at all temperatures. These fluctuations, also visible with the evolution of MSD.

Discussion
We study the different formulations and parameterisation of four potentials and investigate how they result in different amounts of sintering for three different temperatures and particle sizes. In the current work, six measures are used to characterise sintering and to quantify the process. It is hence important to discuss the pros and cons of these measures, since not all measures are able to quantify the sintering process to the same extent. Figure 6: Rotation distribution of atoms for two NPs with the distance from centre of masses of the corresponding particles is plotted. The insets are snapshots of atomic structure contoured with the amount of rotation at various simulation times with respect to their positions at time zero for potential CB and radius R = 4 nm at T= 0.6 T m . The distribution and the snapshots are chosen at the simulations times when MSD changes its evolution path (see Fig. 4c).
In particular, it is visible that the MSD is prone to influences from rigid body motion, resulting in a non-monotonic evolution over time as seen for instance with CB for R = 4 nm and T= 0.6 T m . We note that the simulations are performed with periodic boundary conditions and the MSD is calculated using only the atoms participating in the sintering process, i.e. atoms that detach themselves from the main particle are neglected. To understand the non-monotonic evolution in MSD, we investigate the rigid body rotation of the two particles, see Fig. 6. This is done by using the polar decomposition of the deformation gradient F = RU and expressing the rotation tensor R in terms of an angle that denotes the rotation of an individual atom from the initial configuration. It is evident that the crests and troughs in the evolution of MSD correspond to significant changes in the rotation of the two particles. This non-monotonic evolution is eliminated when the rotational part is removed from the displacement, as is the case with the stretch tensor |U | 2 Interpretation of the colour code: 0: no substantial change relative to the initial configuration 1: minimal formation of the neck region 2: significant progress of the neck region 2.5: sintering of the two particles without remixing of ions 3: sintering with complete remixing. Figure 7: Quantification of the amount of sintering by "visual inspection" for all investigated potentials, particle sizes and temperatures.
(cf. Fig. 5c). MSD, as a quantifying measure of the sintering process, must hence be used with caution. Indeed, we strongly suggest the use of |U | 2 , which essentially encapsulates MSD whilst accounting for the immediate neighbourhood of individual atoms. All six measures aim to quantify the same sintering process. Therefore, they should be able answer the questions a) which of these quantities give a definitive indication on the amount of sintering?, and b) what numerical value of the said quantity indicates complete sintering? In principle all six measures should be consistent with each other if they perform equally well. However, significant discrepancies can be observed in the evolution of these quantities. For instance, MSD (and also |U | 2 ) for R = 4 nm and R = 6 nm at T= 0.8 T m for CTIE is higher than that for all other potentials and particle sizes. This trend is not observed with other quantities, such as the fraction of ions in the neck, F ion , or the shrinkage ratio, S f . Also, it is not completely evident from the quantities themselves if we indeed have complete sintering or not.
How strongly do the six sintering measures correlate with our expectation and experience? Fig. 7 provides a quantification of the amount of sintering obtained via visual inspection. The scale of 0-3 indicates no substantial change in the initial configuration (0) to complete sintering (3).
The normalised surface area A f and the neck curvature are perhaps the most intuitive measures to set numerical thresholds that indicate complete sintering. In the ideal case of two particles being completely sintered and become one sphere, the normalised surface area should decrease to a value of approx. 0.79, and the neck curvature to zero. However, the discrete nature of the system and the approximations used in the two algorithms can result in large fluctuations due to small perturbations, making these quantities sensitive and error-prone.
The shrinkage ratio S f and the fraction of ions in the neck F ion are less intuitive, but seemingly the more robust measures to quantify sintering. In the ideal case of complete sintering, S f should evolve to a value of 1. In the way how we define the neck region, complete sintering would involve remixing of atoms from the two original NPs resulting in the entire sintered particle being classified as neck leading to a value of 1 for F ion . In practice, however, values of approximately 0.7 and 0.8 for S f and F ion , respectively, are seemingly sufficient to indicate complete sintering, while values of 0.5 and 0.35, respectively, indicate partial sintering, i.e., two particles have become a single particle, albeit without significant remixing of ions from the two particles. The latter case is merely a question of kinetics and would result in the former case with increased simulation time. We note, nonetheless, that these numerical values have no physical meaning and are mere observations from our simulations.
A focal point of the current work is the comparison of interatomic potentials for the simulation of sintering. A key component of atomistic simulations is the reliability and applicability of interatomic potentials. Empirical potentials, such as those used in the current work, are usually derived by fitting a functional to certain properties of the material under consideration. All four potentials used in the current work make use of lattice constants and the bulk modulus of alumina in their fitting procedure. Cohesive energy is used additionally for the parameterisation of CTIE and Vash potentials. In this regard, it is unsurprising that the values of these properties fall well within the range of values obtained from experiments (see Fig. 3). All other material properties listed in table 1, including the vacancy formation energies and the surface energies, are predicted values and differ significantly from one potential to another.
The decision on which interatomic potential to use is generally taken based on the accuracy of the predicted material properties, which are usually not a part of the fitting procedure. The two main properties of interest in the sintering process are the surface energies and the vacancy formation energies. The former is of lesser relevance in the current work due to the usage of spherical particles. Nevertheless, we note that all potentials but for BMH, predict surface energies within the range determined by experiments.
The values of the vacancy formation energies -both Al and O-site vacancy energies -predicted by the four potentials are, however, significantly different from one another. While the values of the O-site vacancy energies are very close to experimental values, those of Al-site vacancy energies predicted by CB and BMH potentials are significantly higher than experimental values. In general, the following trend is observed in the predicted values of the vacancy formation energies: Vash < CTIE < CB < BMH, with the lowest values predicted by the Vash potential.
Given the proximity in the values of most material properties to those obtained from experiments, and the low vacancy formation energies, we should expect the best prediction of the sintering process with the Vash potential. However, comparison of the amount of sintering (see Fig. 7) shows completely contradictory results with the least amount of sintering with Vash and furthermore, a near absence of sintering of particles of larger size (R = 4 nm and R = 6 nm) even at T= 0.7 T m and T= 0.8 T m . By contrast, all other potentials show an increased amount of sintering at T= 0.8 T m . The reduced progress in sintering for the highest particle radius of R = 6 nm reflects the slower kinetics, and is likely to result in a completely sintered state with increased time. With the Vash potential, however, an increased simulation time is unlikely to yield any changes as is evident in the characteristics presented in Fig. 4 and Fig. 5.
It hence follows that the origin of these differences lies in the functional form and treatment of the individual potentials. The Vash potential includes both two-body and three-body interactions, the latter being applied only to triplets of atoms. The former, which includes Coulomb interactions, is truncated at r cut = 6Å [24]. Such truncation is unlikely to have an effect on periodic systems as is evident in the accuracy of the material properties predicted by the potential. In the presence of free surfaces, such truncation results in significant discrepancies like pseudo size effects observed in the sintering of larger particles at higher temperatures.
As pointed out in section 2, the formulations of the BMH and CB potentials are effectively identical in the manner they are used in the current work. The BMH potential is indeed a re-parameterisation of the CB potential [23]. The key difference between the two potentials is in the effective charges (Z i ) used for Al and O ions: in the BMH potential charge values close to those in pure oxide systems [76] are used. This difference in charges has a significant effect on the material properties predicted by the two potentials, but has apparently no effect on kinetics of the sintering process. Evidently, for all particle sizes and temperatures considered in the current study, the amount of sintering predicted by the two potentials is almost the same.
Of the four potentials, the CTIE potential results in the highest kinetics of the sintering process. A plausible reason for this is the charge equilibration performed in this potential. In the presence of a heterogeneous electrostatic environment, as e.g. interfaces and free surfaces, ions are likely to have varying charges depending on their local environment. The ability to perform charge equilibration allows for different charges for different ions, depending on their local environment, which we believe is the reason for improved kinetics with the CTIE potential. A further improvement in the simulations can be introduced by activating dynamic charge transfer in the CTIE potential since this is known to enhance atomic diffusion in surface regions of nanoparticles [77].
We hence summarize that accurate treatment of Coulomb effects, together with correct effective charges of the ions, are important for effective prediction of sintering in atomistic simulations of alumina. The usage of effective charges can result in improved kinetics of the sintering process. Truncation of Coulomb effects, however, is more critical and can result in spurious size effects during sintering.

Conclusions
In this study, we perform atomistic simulations of sintering and compare four interatomic potentials for alumina. The four potentials -Vash, CB, BMH and CTIE -differ in the formulation as well as in their treatment of Coulomb interactions and individual charges. The potentials are first compared in terms of their ability to predict fundamental material properties. Three different particle sizes and temperatures are used to study the sintering process, which is subsequently characterised using six different quantities. The findings of the study can be summarised as follows: • All four potentials predict lattice constants that are well within the range of experimental values. Only Vash and CTIE predict cohesive energies that are close to experiments. Both CB and BMH predict significantly higher values of cohesive energy • The values of most elastic constants are significantly higher with the BMH potential.
With the other three potentials, most of the values are very close to those obtained from experiments.
• The surface energy values predicted by the CTIE potential are significantly higher than the experimental values. All other potentials predict values well within the range determined by experiments.
• The vacancy formation energy values of both Al and O sites follows the following trend: Vash < CTIE < CB < BMH. • All potentials show an increasing amount of sintering with increasing temperature for the nanoparticle with R = 2 nm, with complete sintering observed at T= 0.8 T m .
• Very little or no sintering is observed for the largest particle size of R = 6 nm at temperatures T= 0.7 T m and below.
• Except for the Vash potential, all other potentials predict significant amount of sintering at T= 0.8 T m for all particle sizes.
• For the quantification of the sintering process, the use of |U | 2 is recommended over MSD, since the latter is particularly prone to influences of rigid body motion of the particles.
• The shrinkage ratio (S f ) and fraction of ions in the neck (F ion ) provide the best quantification of the sintering process. Numerical values of 0.7 and 0.8 for S f and F ion , respectively, indicate complete sintering.
• The abrupt cutoff of Coulomb interactions is deemed as the reason for the spurious behaviour of the Vash potential, which shows almost no sintering for particle sizes above R = 2 nm, even at the high temperature of T= 0.8 T m .
• The CTIE potential predicts the highest amount of sintering amongst all potentials compared in the current work. The increased kinetics is attributed to the ability to perform charge equilibration within the potential.  Table 2: The potential parameters of four potentials. The parameters for CTIE are from the original source [26] and the parameters of EAM potential are from Ref. [35].

Supplementary Material
0.0.1. Method of calculating melting temperature An energetically minimized periodic alumina system is heated up to several temperatures for allowing volume expansion and equilibrated by NPT at the target temperatures for minimum of 100 ps. The steady state volume per atom is plotted against temperatures for Vash in Fig. S1 as an example.
Supplementary Figure S1: The steady state volume per atom, obtained at various temperatures for Vash potential, is plotted with temperature. The window of temperature within which the volume jump occurs is considered to be the melting temperature of the potential. The melting temperature is 2760 K for Vash potential associated with the biggest jump in volume per atom as shown in the figure.
1. Information concerning the computation of material properties at 0 K 1.1. Lattice constants, cohesive energy and vacancy formation energy A 3D alumina system of 8640 atoms (3456 Al + 5184 O) atoms, periodic in all directions are energetically minimized to calculate the energy per atom (cohesive energy) and lattice constants. The vacancy formation energy is calculated by taking a difference in energy between a minimized infinite pristine system and with a vacancy in it.

Surface energies
For surface energy calculations, infinite slabs with the desired crystallographic orientations in the free surface direction are minimised in energy. The value of the surface energy is then calculated as: where E slab is the energy of the slab, E coh is the cohesive energy, A is the free surface area, and N is the total number of atoms in the slab.
1.2.1. Identifying the appropriate slab thickness In order to minimise the influence of the thickness on the calculated surface energy, the same calculation is performed for several slab thicknesses. The calculated surface energy increases with the slab thickness, and it stabilises with respect to thickness beyond a certain value (see Fig. S2). Based on this study, around 60Å is chosen for the minimum slab thickness to calculate surface energies.
Supplementary Figure S2: Surface energy of [0001] for slabs of different thicknesses calculated with Vash potential.

Elastic constants
The system is deformed by applying a small strain in positive and negative Voigt strain directions. The elastic constants are calculated by taking the derivatives of the measured change in stress tensor with respect to strain.

Neck region
In Fig. S3 the neck region of the sintering nanoparticles is highlighted by blue to green contour for R = 6 nm particle and Vash. The details on defining a neck region is described in the main text.
Supplementary Figure S3: The neck region between two particles during sintering (highlighted with a color gradient). All other atoms are colored according to their chemical species: Al -gray, O -red.

Influence of rigid body motion of particles on MSD
Fig . S4 shows the centre-to-centre distance and particle coordinates as a function of the sintering time.
Supplementary Figure S4: The rigid body motion pf particles contribute to the MSD calculation. The centre-to-centre distance between particles remains almost unchanged. The centre of masses of the particles do not change substantially in x and z direction. The particles oscillate with respect to their centre point in the y direction resulting a local fluctuation in MSD.