Thermodynamics of CuPt nanoalloys

The control of structural and chemical transitions in bimetallic nanoalloys at finite temperatures is one of the challenges for their use in advanced applications. Comparing Nested Sampling and Molecular Dynamics simulations, we investigate the phase changes of CuPt nanoalloys with the aim to elucidate the role of kinetic effects during their solidification and melting processes. We find that the quasi-thermodynamic limit for the nucleation of (CuPt)309 is 965 ± 10 K, but its prediction is increasingly underestimated when the system is cooled faster than 109 K/s. The solidified nanoparticles, classified following a novel tool based on Steinhardt parameters and the relative orientation of characteristic atomic environments, are then heated back to their liquid phase. We demonstrate the kinetic origin of the hysteresis in the caloric curve as (i) it closes for rates slower than 108 K/s, with a phase change temperature of 970 K ± 25 K, in very good agreement with its quasi-thermodynamic limit; (ii) the process happens simultaneously in the inner and outer layers; (iii) an onion-shell chemical order - Cu-rich surface, Pt-rich sub-surface, and mixed core - is always preserved.

The validation for the parametrization for Cu-Pt interaction used in the manuscript is hereby discussed. Parameters for Pt-Pt and Cu-Cu interaction have been found in the literature and their average was initially used to model Pt-Cu interaction. They are numbered according to their position in Table 1  The sets are then validated against the energy ranking found for density functional theory level of accuracy, of a variety of relaxed structures made of 55 atoms. The architectures considered are clusters of 55 atoms with Cuboctahedron (Co), Icosahedron (Ih), Decahedron (Dh) morphology and Core-shell(cs), Core-Vertex(v), Core (C)-Edge(E)-Vertex(V) impurity chemical compositions. DFT calculations are performed using the Quantum Espresso package 6 using Ultrasoft potential Pardew-Burke-Ernzerhof ? exchange model within the local spin density approximation. An Energy cutoff = 45.0 eV and a Plane wave cutoff = 360.0A have been used. Calculation sample a single Γ Point in the Brillouin zone.
The parameters proposed in Ref. 7 , corresponding to the parametrization number 6, results one of the best performing set as shown by figure 1. In the case of mono-metallic and single impurity systems the energetic ranking is perfectly reproduced. The unique difficulty in reproducing the DFT data is observed for the case of the Pt vertex decorated clusters, hinting to a significant charge transfer in these architectures as also found in the mono-metallic Pt case 8 . Further tuning of the Cu-Pt interaction towards more Cu-or Pt-like interaction in the spirit of Reference 9 didn't result in a significant improvement of the matching between DFT and MD rankings for any of the parameters sets here investigated.

Figure 1.
Record of the ∆Energy of each configuration as calculated through DFT self-consistent level and TBSMA potentials. Energy is rescaled with respect to the bulk cohesion energy contributions. CS and V label refer to a Cu-rich nanoalloy, ICS and IV to a Pt-rich one. Cu 55 energy ranking is shown on the left, Pt 55 on the right. TBSMA potentials parameters sets are color coded. Lines are meant as guides to the eyes only. 3/17

Nested Sampling detailed explanation
Nested sampling estimate the density of states of a system by constructing a sequence of decreasing potential energy levels E i , starting from the vapor phase, each of which bounds from above a volume of configuration space χ i , with the property that χ i is approximately a constant factor smaller than the volume χ i−1 , corresponding to the level above. Each volume is sampled uniformly, and therefore each distribution will have an approximately constant fractional overlap with the ones immediately before and after, ensuring fast convergence of the sampling and allowing an accurate evaluation of phase-space integrals. In particular, the energy levels near the phase transition, where phase volumes change rapidly, will automatically be very narrowly spaced. The sequence of energy levels comprise a discretization of the cumulative density of states χ(E), which allows the evaluation of the partition function at arbitrary temperatures, where N is the number of particles of mass m, V is the volume, β is inverse thermodynamic temperature, h is Planck's constant, the density of states χ is the derivative of χ, and we labeled the factor resulting from the momentum integral as Z m . In each iteration, we start with a set of samples uniformly distributed and bounded above by energy E i , and a new energy level is obtained by recording the energy of the sample with highest energy, discarding that sample, and reconstructing the uniform distribution under the new energy limit by cloning a randomly chosen sample from the remaining ones, and de-correlating it from its source sample using a random walk, always obeying the new energy limit. The partition function can be evaluated a posteriori at any temperature and other thermodynamic variables, such as the internal energy, can be evaluated as a function of temperature, 4/17

Iterative molecular dynamics
Molecular Dynamics is a powerful method to integrate the Newton equations of motions and it encodes information on representative trajectories at finite temperatures. The iterative temperature Molecular Dynamics (itMD) consists of concatenated MD runs each of them long ∆τ, where the temperature in two successive runs is increased by the quantity ∆T . This procedure allows to analysis and report the kinetics of the phase change in terms of a unique parameter, named as temperature change rate, λ which is given by the ratio ∆T ∆τ . itMD allows to estimate caloric curves taking into account potential memory effects at low and high temperatures and it is a more realistic way to mimic the solidification and liquefaction of a system. Furthermore, it allows to detect finite temperature shape transitions whenever their characteristic time is lower than the simulation time at each temperature ∆τ. As standard output, we report the total energy, potential energy at each instantaneous temperature, as calculated using the equipartition theorem. From them we reconstruct the caloric curve per each simulation. Results are averaged over the set of same time interval (corresponding to same nominal temperature too) runs. Statistical analysis is obtained averaging over Figure 2. Caloric curve data as obtained from melting and freezing itMD procedures. Each point is averaged over at least 25 independent simulations and the error is the statistical standard deviation. The error on the temperature is too small for being reported.
the various -at least 25-independent itMD simulations performed per each λ value. Figure 2 reports the caloric curves for both melting and freezing procedures. Along the energy axis, each data point is plotted together with its standard deviation. As expected, the statistical error on the energy increases closer to the phase change point, while it is almost negligible in the liquid phase region.

5/17
2 Additional examples of cluster taxonomy within our new 2D gIh projection scheme The mapping of peculiar low-symmetry motif sampled by itMD and Nested sampling according to the 2D giant Icosahedron projection scheme may result not immediate to novel users of this framework. We recap in Figure 3 seven paradigmatic and most common cases of cluster mapping: a) Cp cluster with all its atoms in an either fcc or hcp arrangement. b) Dh cluster with its single five-fold symmetry axis, located off-centre. c) Dh cluster with its single five-fold symmetry axis, located in-centre. d) Ih cluster with three of its five-fold symmetry axes present and these intersect outside the cluster. e) Ih cluster where its six five-fold symmetry axes intersect off-centre. f) Ih cluster with a single five fold symmetry axis near the edge of the cluster, but with an additional sheet of hcp grain boundary at 60 • angle indicating the potential positioning of the additional five fold axes and forming an Ih structure. The cluster is also shown from a different angle. g) Ih cluster with several five-fold symmetry axes, some forming triangles, and one icosahedral centre. An additional snapshot where the white atoms are removed to highlight the five-fold axis triangular arrangements. Dashed lines and circles are guides to eyes.

Caloric curves of CuPt nanoalloys
Caloric curves for the itMD runs with fixed ∆τ or fixed ∆T rates and for the NS runs with both 4.8 × 10 4 and 9.6 × 10 4 MC steps are shown in Figure ??. The hysteresis loop of the itMD caloric curves gets narrower with slower rates both when fixing ∆T or ∆τ in the simulation. The shorter NS runs slightly underestimate the phase transition temperature. The itMD data shown are found by averaging temperature and energies of the system over the ensemble of independent runs with the same temperature change rate.  itMD-caloric curves (∆E is the total energy rescaled with respect to the potential energy of a perfect core-shell icosahedron) with λ equal to 5K/ns, 10K/ns and 25K/ns, shown by light, mid and dark symbols, respectively (hot colours for the melting and cold for freezing). The time interval per each temperature is fixed to 1 ns. For comparison, the NS solidification limits, arising from the seven different set-ups listed in the text, are shown by dashed lines (4.8 × 10 4 MC steps) and solid black lines (9.6 × 10 4 MC walk). Bottom panel: itMD-caloric curves at different λ fixing the temperature change at 25 K and varying ∆τ from 1 ns (darkest blue/red) to 250 ns (cyan/yellow). In the inset we report the temperature difference between melting and freezing temperatures as a function of λ . The closure of the hysteresis is evident.

Discussing the nanoalloy's taxonomy during phase changes
Let us first to discuss the taxonomy of CuPt nanoalloys during the exploration with NS. The structural characterization for all seven NS runs are shown in Figure 5. It is clear that NS explores multiple basins in a single run, however in different extent. All NS calculations explored configurations with at least one five-fold symmetry axis, and the icosahedral central atom was always found to be at the very edge of the cluster, except in run III, suggesting that the relative probability of a regular icosahedron is very low. Due to the large size of the cluster and limits in computational resources we had to compromise the sample size used in the calculations. As a result, only two of the runs, II. and VII, explored all four major structural groups and only the results from these two are presented in the manuscript.

itMD clusters taxonomy and morphological transitions analysis
In a similar manner as for NS, we analyze the taxonomy of CuPt during itMD-freezing simulations at different λ rates. Here we illustrate both Ih and p-Ih motifs where the p-Ih label highlights the so called poly-icosahedral clusters, that is, motifs with more than one icosahedral centre, as defined before. In the main text we sum over structures showing a single and more icosahedral centre. Fig. 6 reports the structural classification of the cluster onto the 2D giant icosahedral map when ∆T = 25 K and varying ∆τ over two order of magnitude, namely between 1 to 250 ns, at 700 K. We also report in Fig. 7 data for the simulations where the temperature rate ∆T is changed (25K, 10K, and 5K) but we keep fixed the time period between any temperature change, ∆τ .
itMD was found to explore multiple basins after nucleation only in the case of a freezing rate of 25K/250 ns. In other words, regardless fixing ∆τ or ∆T , if the rate λ is larger than 1K/ns, after nucleation, the clusters never undergo any major structural transformations. Structural changes are shown for nine of these runs as a function of temperature in Figure 8. This constitutes an indirect proof that the characteristic timescale for structural transitions below the melting temperature of our system, Cu 162 Pt 147 , is below 250 ns but much slower than 25 ns.  Each square represents an independent simulation. The green circle shows the nanoparticle shape at 700 K onto the 2D giant icosahedral map represented by red lines, as described in the main text. In the case of ∆τ = 250 ns, the green circle refers to a temperature close to 700 K and its cluster evolution is reported in Fig. 8.  Figure 7. Structural classification of itMD freezing simulations with λ equal to 5 K/ns (left), 10 K/ns (center), and 25K/ns (right), fixing ∆τ = 1 ns and varying ∆T accordingly. Each square represents an independent simulation. The green circle stands where a structure can be located onto the 2D giant icosahedral map, represented by red lines, at 700 K. The shape name is given in each panel. Paradigmatic examples of the temperature evolution of local atomic order parameters such as common neighbor analysis (CNA) signatures and Steinhard parameter Q6 and W6, are depicted in Figure 9, and 10. For the case of CNA signature, paradigmatic bulk signatures are analyzed: the (5,5,5) which detects pair of atoms lying on five-fold symmetry axis, the (4,2,2) which corresponds to pair of atoms along grain boundaries or twinning planes, and the (4,2,1) which characterizes FCC-bulk pairs. For (Q6,W6) pairs, we report specific coordinates corresponding to the ideal value for atoms whose neighbors present an Icosahedral, FCC, HCP, and five-fold symmetric symmetry. The solid and the liquid regions are clearly distinguished by a sharp transition of the order parameters under investigation. We note that while the CNA and (Q6,W6) values for the melted phase are representative of the ones found for all the Nested Sampling and itMD runs, their magnitude in the case of a solid structure depend on the morphological basin sampled.

12/17
Our NS simulations show an ensemble of liquid, frozen, and half-liquid/half-frozen configurations. Atomic Q6 and W6 parameters are also reported to distinguish solid and liquid atoms according to the mapping shown in Figure 10  Data are gathered from three independent NS runs. Also in this case note the correspondence between the nucleation of solid structure and the clustering of the Q6 and W6 values around characteristic values peculiar of local atomic environments in nanoclusters. The ICO, 5-fold(*), HCP, and FCC label refer to characteristic (Q6,W6) values for icosahedral centres, atoms in a (in)complete five-fold axis, or in FCC and HCP sheets respectively.

15/17 4 Chemical reordering in CuPt nanoalloys
To check if chemical re-ordering might happen on a time scale similar to shape transitions, we run additional itMD melting simulations starting from unfavorable chemical ordering, for example, a segregated one as in the case of a core-shell. Fig 11 shows the time evolution of the radial chemical distribution function of a Cu 162 Pt 147 cluster with an initial Pt core -Cu shell at 300K, 600K, and 900K. Within the time scale probed in this simulation, 100 ns , the chemical ordering is preserved for the two lower temperatures, and strongly altered, leading to the formation of the favorable onion-shell pattern, only for the 900 K run. Figure 11. Time evolution of the ratio of the Cu atoms within different 0.01A thick layers of a Cu 162 Pt 147 cluster at 300K, 600K, and 900K when starting from a perfect icosahedral morphology and a Pt core -Cu shell chemical ordering.