Puzzles of Surface Segregation in Binary Pt–Pd Nanoparticles: Molecular Dynamics and Thermodynamic Simulations

: Up till now, there have been extremely contradictory opinions and inadequate results concerning surface segregation in binary platinum–palladium (Pt–Pd) nanoparticles, including the problems regarding segregating components, as well as the size and temperature dependences of segregation. Taking into account such a situation, we investigated the surface segregation in Pt–Pd nanoparticles by combining atomistic (molecular dynamics) and thermodynamic simulations. For molecular dynamics experiments, the well-known program LAMMPS and the embedded atom method were employed. In the course of the atomistic simulations, two different sets of param-eterizations for the Pt–Pt, Pd–Pd, and Pt–Pd interatomic interaction potentials were used. The thermodynamic simulation was based on solving the Butler equation by employing several successive approximations. The results obtained via atomistic simulation and thermodynamic simulation on the basis of the Butler equation were compared with each other, as well as with predictions that were based on the Langmuir–McLean equation and some experimental data. Both simulation methods (atomistic and thermodynamic) predicted the surface segregation of Pd, which diminishes with the nanoparticle size and with increasing temperature. Our simulation results do not conﬁrm the predictions of some authors on surface segregation inversion, i


Introduction
Particles of nanometer size can exhibit properties that differ from those of their corresponding bulk materials. Nowadays, metallic nanoparticles (NPs) are considered an important type of nanomaterials that have, in some cases, unique physical and chemical properties; metallic NPs provide for many applications in the fields of optics, optoelectronics, chemical production, biosensors, catalysis, and other directions of nanotechnology [1,2]. In particular, plasmonic applications should be noted [2]. Metallic and metal-based NPs can also demonstrate unique tuned magnetic properties [3][4][5][6]. In addition, NPs play a role in the building blocks of fabricated nanostructured materials; for example, via powder metallurgy technologies that are based on NP sintering [7]. In the last decade, applications of metallic NPs in biotechnology and medicine have also been reported [8][9][10]. However, nanocatalysis, including that in fuel cells, seems to be the main area of application for metallic NPs [11][12][13].
Regarding monometallic NPs, gold NPs have been discussed in the literature in great detail [14,15]. Among monometallic catalysts, Pt seems to be most effective as it is employed both as separate active atomic sites in organic compounds, and as the NPs in some metallic structures, i.e., the constituents of heterogeneous catalysts [16][17][18]. According to [16], the automotive industry is the main consumer of the noble metals platinum (Pt), palladium (Pd), and rhodium. These metals are used as components of the automotive catalysts that are applied to decrease the emissions of carbon monoxide, unburnt hydrocarbons, and nitrogen oxides in exhaust gases.
(i) Bond strengths, i.e., the relative strengths of the A-A, B-B, and A-B bonds between metal atoms. The components that form stronger bonds should tend to the core positions, while stronger A-B bonds will favor mixing (we believe that the bond strengths are determined, first of all, by the binding energies of components A and B); (ii) Surface energies: the metal with the lower surface energy (and the lower surface tension) will segregate to the NP surface; (iii) Atomic sizes that relate to the release of strains in the NPs with smaller atoms that tend to the core positions; (iv) Specific electronic effects for certain metals (i.e., the directionality of bonds, according to [27]); (v) Charge transfer effects: the most electronegative atom tends to occupy a surface site.
The last two factors should be taken into account in certain specific cases only, e.g., in the case of a 38-atom Ag-Au nanocluster (see Chapter 8 in monograph [47]). Therefore, if the atomic size mismatch is small enough, the two first factors affecting the surface segregation seem to be considered basic. In addition, all the above factors relate to each other. In particular, the strengths of bonds determine the surface energies (i.e., the surface tensions) of components. According to the thermodynamic simulation method that was based on solving the Butler equation [48] considered in Section 4, the surface segregation in the A-B nanoalloys was determined, first of all, by the difference (σ On the one hand, the Pt-Pd system may be considered simple enough for investigation. In reality, Pt and Pd have similar properties, including an absence of the atomic size mismatch and a simple type of phase diagram that corresponds to the unlimited mutual solubility of components. Only the first two inter-related factors of segregation (among the five listed above) may be of importance, and these factors can be adequately taken into account in molecular dynamics (MD) simulations. On the other hand, properties of Pt, Pd, and Pt-Pd NPs are studied much less frequently than properties of NPs of metals with Metals 2023, 13, 1269 3 of 20 lower melting temperatures. In particular, as noted in ref. [49], experimental data on the size dependence of the melting temperatures of Pt and Pd NPs seem to be absent in the literature. Experimental data on the segregation of components in Pt-Pd NPs are quite rare [50][51][52], and relate only to a few NP compositions and sizes or to quite narrow size ranges. For example, experimental data [50] relate to the Pt-Pd NP size region from 1.8 to 4.2 nm, while experimental data [52] relate to Pt-Pd NPs of 3.8, 4.5, and 4.8 nm in diameter. However, all the available experimental data predict the surface segregation of Pd.
At the same time, some theoretical models and atomistic simulations predict that the surface segregating component in Pt-Pd NPs can change at a critical size. For example, in [53], a critical value α of the NP diameter d is introduced, below which the surface segregation of Pd should be switched to the surface segregation of Pt. Then, MD results [54] predict that the segregation in Pt-Pd NPs of 2 nm in size is practically absent, whereas in NPs 4 nm in size, the surface segregation of Pd switches to the surface segregation of Pt when the molar fraction of Pd is smaller than 0.6. Experimental data and accurate density functional theory (DFT) results [55] (see also the references cited in [56]) for small Pt-Pd nanoclusters do not corroborate such a transition. Thus, it would be reasonable to investigate the size dependence of surface segregation in Pt-Pd NPs using the embedded atom method (EAM), i.e., the main alternative to other interatomic potentials used for metals including the Sutton-Chen potential [57] employed in [54].
Notably, the authors of paper [53] concluded that surface segregation in Pt-Pd NPs increases with decreasing the NP size. This statement contradicts the other conclusion about the inversion of the surface segregation made in the same paper. The temperature dependence of surface segregation in binary metallic NPs, including those composed of Pt-Pd, is also a controversial issue. In the present study, the surface segregation in Pt-Pd NPs, including their size and temperature dependences, was investigated using MD and thermodynamic simulations. The thermodynamic results obtained in the frames of thermostatics relate to the equilibrium surface segregation, whereas the MD simulations reproduce kinetics of segregation at the time scales available in the MD experiments. In general, the obtained MD and thermodynamic results agree with each other, although some kinetic effects were also revealed in our MD simulations.
In ref. [49], we revealed some problems related to parameterizations of the interatomic interaction potentials that are employed for atomistic simulations in the frames of EAM. This refers to the EAM parameterizations proposed by Zhou and his co-authors [58] for 12 metals, including Pt and Pd. As noted in ref. [49], we found that parameterizations [58] inadequately predict a higher melting temperature for Pd than for Pt. For this reason, in ref. [49], we recalculated the embedding functions for Pt and Pd to adequately predict the relationship between the bulk melting points of Pt and Pd. In the present paper, our MD results obtained by employing both EAM parameterizations ( [49,58]) are compared. The abovementioned drawback of parameterizations [58] made it possible to establish whether the relationship between the melting temperatures of components A and B noticeably affects the surface segregation in A-B nanoalloys.
It is also worth noting that for the last decade, the scientific community has shown interest in trimetallic [34,59,60] and even quaternary [61] NPs containing Pt and Pd as two of the three or four components. As we believe, the above open questions concerning binary Pt-Pd NPs should be solved before proper detailed atomistic simulations of ternary and quaternary Pt-and Pd-based nanoalloys can be started.

Approaches to Molecular Dynamics Simulation
For our MD simulations, the well-known open program LAMMPS (version-23 June 2022, Sandia National Laboratories, Albuquerque, NM, USA) [62], employing EAM to describe the interatomic interactions in Pt-Pd NPs, was used. LAMMPS makes it possible to perform parallel computing on graphics processing units (GPUs), which is one of the advantages of this program. The standard velocity Verlet integrator [63] and the Nosé-Hoover thermostat [64] were used in our MD simulations with the time step t 0 = 1 fs usually being chosen in contemporary MD simulations. We have found that the value t * damp = t damp /t o = 10 of the reduced damping parameter, i.e., of the relaxation time t damp , seems to be optimal for our MD simulations. Indeed, our preliminary MD simulations showed that for molecular (Lennard-Jones) systems and for metals at constant temperatures, t * damp = 100 may be suitable. However, our MD simulations included heating and cooling NPs, and a damping parameter of 100 does not ensure proper controlling temperature during these processes. No significant difference was noted in the results of NP relaxation at the final constant temperatures obtained for t * damp = 10 and t * damp = 100. Although a more systematic study on selecting the damping parameter in isothermal MD simulations would be necessary, this problem is beyond the scope of our study.
As mentioned in the Introduction, two EAM sets of parameterizations were used: ours, [49], and that proposed in [58]. In both cases, the algorithm in [58] was employed to generate the Pt-Pd potential on the basis of Pt-Pt and Pd-Pd potentials. Both sets of parameterizations were verified on the bulk Pt and Pd phases, Pt and Pd NPs in our previous paper [49]. According to Tables 1 and 2, both sets of parameterizations predict close and adequate values of densities of bulk Pt and Pd in solid and liquid states, but parameterizations [58] inadequately predict a higher bulk melting point T    Unlike ab initio methods, MD simulations reproduce kinetics of the process in the systems under study. The kinetics of segregation should strongly depend on the NP size and temperature. Due to parallel calculations on the GPUs, LAMMPS significantly expands the range of sizes of simulated objects and the time of the reproducible MD evolution. No doubt, in MD simulations, the choice of the interatomic potential and its parameterization largely determines the reliability of the simulation results. EAM and the tight-binding potential [67] can be considered as the two main and well-proven types of many-body potentials commonly used in simulations of metallic systems. In ref. [49], we verified our parameterizations and additionally verified the parameterizations [58] which were previously tested by the authors of [58] themselves.
The performance of these potentials in describing the Pt-Pd interaction is of special importance. The validity and reliability of the Pt-Pd interaction potential was not substantiated in [49]. However, according to [58], an algorithm proposed by the authors of this paper and employed in the present study provides a reasonable approximation of the crossinteraction potentials, including Pt-Pd. In particular, this algorithm adequately predicts the heats of solutions. In ref. [49], we performed more accurate calculations of the embedding functions for Pt and Pd, which determined qualitatively correct relationships between the melting temperatures of Pt and Pd. Obviously, a higher correctness of the embedding functions should not lead to a lower correctness of the Pt-Pd parameterization; our present work can be considered as a further validation of the EAM Pd-Pt interaction potential.
There is no doubt that the decrease in the NP size affects the electronic configuration of Pt and Pd atoms, and this effect cannot be fully taken into account in the framework of EAM. However, as mentioned in Introduction, accurate DFT calculations for small Pt-Pd nanoclusters [55,56] confirm our MD results, at least indirectly.
The initial configuration of a binary A-B NP corresponded to a spherical fragment of the bulk fcc lattice with the random distributions of the A and B components. According to [68,69], Pt NPs and Pt-based nanoalloys can grow as truncated octahedral, octahedral, and even tetrahedral structures. However, when subsequent melting, cooling, and annealing the initial configuration are employed to promote segregation, the shape and structure of the initial configurations are less important. After placing the initial configuration into the center of the 60 nm × 60 nm × 60 nm simulation box, the particle was relaxed (annealed) at a sufficiently low temperature, T 0 , of 300 K, and then uniformly heated up to a temperature T 1 , knowingly higher than the melting points, T (∞) m , of both components. Such a procedure made it possible to reach the melted state of the NP. Then, the NP was gradually cooled at a cooling rate of . T = 0.12 TK/s, low for MD experiments, down to a temperature, T (at which then the segregation parameters would be determined), and finally, the NP relaxed for about 100 ns at the chosen fixed temperature. The subsequent binary NP melting, cooling, and annealing before registration of the segregation were used to facilitate the rearrangement of the initially non-equilibrium particle structure, i.e., to promote the potential trend to the spontaneous surface segregation of one of the NP components. Segregation was studied when the particle structure and its properties had not changed after relaxation (annealing), i.e., when the particle state could be treated as equilibrium or close to it. The relevance of such an approach will be confirmed by agreement between our MD simulation results and the results of the thermodynamic simulation, as well as some experimental and DFT results of other authors.

Results of Molecular Dynamics Simulation
As noted in the previous subsection, the initial configurations of binary Pt-Pd NPs corresponded to spherical fragments of the fcc lattice with random and approximately uniform distributions of atoms of both components. A snapshot and a central cross-section of one of the initial configurations containing N = 5000 atoms (40% Pd, particle size 5.3 nm) are presented in Figure 1. The initial radial distributions of the local molar fractions χ Pt (r) and χ Pd (r) of Pt and Pd, respectively, are shown in Figure 2.
The results corresponding to the final MD configurations are presented in Figures 3 and 4. These results were obtained using EAM parameterizations [58]. The figures clearly demonstrate the surface segregation of Pd in the final configurations of binary Pt-Pd NPs. The thickness of the surface outer layer, δ, can be chosen arbitrarily, for example, as the thickness of the surface monolayer or of a layer that depicts a noticeable surface segregation of one of the components. In our case, it is advisable to visually determine the spherical surface of radius r c = 2.4 nm as a geometric boundary between the particle core and its surface layer (shell). The thickness of the surface layer, δ, will be equal to r 0 − r c ≈ 0.3 nm, which corresponds to the surface monolayer enriched with Pd atoms (r 0 is the particle radius).  The results corresponding to the final MD configurations are presented in Figures 3 and 4. These results were obtained using EAM parameterizations [58]. The figures clearly demonstrate the surface segregation of Pd in the final configurations of binary Pt-Pd NPs. The thickness of the surface outer layer, , can be chosen arbitrarily, for example, as the thickness of the surface monolayer or of a layer that depicts a noticeable surface segregation of one of the components. In our case, it is advisable to visually determine the spherical surface of radius = 2.4 nm as a geometric boundary between the particle core and its surface layer (shell). The thickness of the surface layer, , will be equal to 0 − ≈ 0.3 nm, which corresponds to the surface monolayer enriched with Pd atoms ( 0 is the particle radius).  The results corresponding to the final MD configurations are presented in Figures 3 and 4. These results were obtained using EAM parameterizations [58]. The figures clearly demonstrate the surface segregation of Pd in the final configurations of binary Pt-Pd NPs. The thickness of the surface outer layer, , can be chosen arbitrarily, for example, as the thickness of the surface monolayer or of a layer that depicts a noticeable surface segregation of one of the components. In our case, it is advisable to visually determine the spherical surface of radius = 2.4 nm as a geometric boundary between the particle core and its surface layer (shell). The thickness of the surface layer, , will be equal to 0 − ≈ 0.3 nm, which corresponds to the surface monolayer enriched with Pd atoms ( 0 is the particle radius).  Figure 1) after cooling the melted NP down to = 300 K and subsequent annealing at this final temperature.   Figure 3. A snapshot (a) and a central cross-section (b) of the same Pt-Pd NP (corresponding to the initial configuration shown in Figure 1) after cooling the melted NP down to = 300 K and subsequent annealing at this final temperature.     According to Figures 5 and 6, surface segregation decreases with the particle size; this result agrees with our MD results obtained previously for other nanoalloys [69][70][71]. As one can see, both EAM parameterizations ( [58] and [49]) predict practically the same results for the ( ) (̅̅̅̅̅ ) dependences, except that our potentials predicted a slightly higher surface segregation of Pd. Thus, inadequately higher melting temperature of Pd (in comparison with Pt) predicted by parameterizations [58] for Pt and Pd potentials does not noticeably affect the surface segregation in binary Pt-Pd NPs. Our MD and thermodynamic results on the temperature dependence of the surface segregation in Pt-Pd NPs are presented in Section 4. According to Figures 5 and 6, surface segregation decreases with the particle size; this result agrees with our MD results obtained previously for other nanoalloys [69][70][71]. As one can see, both EAM parameterizations ( [49,58]) predict practically the same results for the χ (s) Pd (χ Pd ) dependences, except that our potentials predicted a slightly higher surface segregation of Pd. Thus, inadequately higher melting temperature of Pd (in comparison with Pt) predicted by parameterizations [58] for Pt and Pd potentials does not noticeably affect the surface segregation in binary Pt-Pd NPs. Our MD and thermodynamic results on the temperature dependence of the surface segregation in Pt-Pd NPs are presented in Section 4.

Approach to Thermodynamic Simulation
It is hardly possible to logically draw a definite line between the conventional thermodynamic calculations and the thermodynamic simulation. However, following our publications [70][71][72], we use the term "thermodynamic simulation" when an equation or a set of equations is solved by employing some sequential approximations or when exact values of some necessary input parameters are not known. In this study, we used a thermodynamic simulation method based on solving the Butler equation [48], derived in 1932. Much later, this equation was additionally justified and verified by Kaptay on liquid and solid alloys, including nanoalloys [73] and nanograined metals [74].
Following our papers [70][71][72], we arbitrarily divided a binary NP under consideration into a central area (core) and a surface layer (shell) surrounding the core. In particular, the outer surface monolayer can be considered as the NP shell. The goal of the thermodynamic simulation is to determine the equilibrium distributions of components in the NP core and shell. In our case, thermodynamic simulations may be performed by employing several thermodynamic models and approximations. In refs. [70][71][72], we decided to distinguish between two thermodynamic models-of the unlimited source for the segregating component and of the limited source-as well as between two sequential approximations-of the ideal solution (the partial excess Gibbs energies of mixing ∆G m,i of components are zero) and approximation, where relevant non-zero values of ∆G m,i are taken into account. The notation ∆G m,i corresponds to the partial molar excess Gibbs energy of mixing for component i. In the frames of our two-cell NP model, the particle core is treated as the source for the segregating atoms.
the specific partial surface free energies, σ i , of all the particle components (including components i and j of a multicomponent system) should be equal to each other (see additional comments in [73,74] and in our papers [70][71][72]). The specific partial surface free energy, σ i , of component i can be presented as the next function of the relative molar fraction χ i of this component: where R is the molar gas constant, ω i is the partial molar surface of component i, i.e., the partial surface area in m 2 mol −1 , ω m,i are the partial molar excess Gibbs energies of component i in the particle core and the shell, respectively. Hereafter, superscript c and s are used to indicate the particle core and the surface shell, respectively. In the present paper, we have omitted the derivation of Equation (2) since a very detailed derivation of an analogue of this equation but for a nanograin surrounded by the grain boundary was given in [74]. In (2), the specific partial surface free energy, σ i , of component i features instead of the partial free energy of the component in the grain boundary. In addition, in (2), the partial molar excess Gibbs energies, ∆G

Results of Thermodynamic Simulation
The results of our thermodynamic simulations of the surface segregation in binary Pt-Pd NPs of various sizes are shown in Figure 7. These results corresponding to T = 300 K were obtained when using the ideal solution approximation (∆G m,B = 0), i.e., when solving the system of Equations (A6) written for the case of the model of the limited source for the segregating component. In general, these results are consistent with the results of our atomistic simulations predicting the surface segregation of Pd, which decreases with the decreasing NP size. The values of the specific surface free energy, σ (0) , for Pt and Pd in the solid state (necessary when calculating the segregation coefficient (A4)) were taken from ref. [75]. Figure 8 demonstrates similar results, but were obtained when taking into account that the excess Gibbs energies were nonzero. Considering the excess Gibbs energies resulted in a significant decrease in the surface segregation (by 40% on average), but the surface segregation of Pd was not reversed into the surface segregation of Pt even at the smallest NP size. K were obtained when using the ideal solution approximation (∆ ,A = ∆ ,A = ∆ ,B = ∆ ,B ( ) = 0), i.e., when solving the system of Equations (A6) written for the case of the model of the limited source for the segregating component. In general, these results are consistent with the results of our atomistic simulations predicting the surface segregation of Pd, which decreases with the decreasing NP size. The values of the specific surface free energy, (0) , for Pt and Pd in the solid state (necessary when calculating the segregation coefficient (A4)) were taken from ref. [75].   coefficient (A4)) were taken from ref. [75].

Discussion
Our thermodynamic results obtained within the framework of thermostatics were related to equilibrium surface segregation, whereas MD simulations reproduce kinetics of segregation at the time scales available in the MD experiments. However, the results of the MD and thermodynamic simulations were consistent with each other and predicted the surface segregation of Pd in binary Pt-Pd NPs. As noted in the Introduction, the surface segregation of Pd was experimentally confirmed [50][51][52], but these experimental results related to Pt-Pd NPs of some certain compositions and sizes. Obviously, this is due to a number of difficulties associated with synthesis and following investigation of NPs in arbitrary wide ranges of their sizes and compositions. As a result, most of the available experimental data and theoretical results, as well as the results of atomistic and ab initio simulations, confirm the presence of Pd in the surface layers of Pt-Pd NPs, albeit only qualitatively. For example, the authors of [76] found that there was a high number of Pd atoms on the surface of synthesized Pt-Pd NPs. Then, the MD and DFT results [77,78] show that Pt atoms are predominantly localized in the Pt-Pd NP core, whereas Pd atoms on the NP surface. Two other interesting indirect confirmations, [79,80], may also be noted. The authors of [79] experimentally found that the sintering of Pt and Pd NPs results in the formation of Pt-Pd NPs with Pt-rich cores and Pd-rich shells. According to [80], synthesized Pd@Pt dendrimer-encapsulated NPs invert into ones with Pt-rich cores, and "primarily Pd in the shells". We believe that, in both cases, the effect of Pd enrichment in the surface layers of synthesized NPs reflects a trend to the spontaneous surface segregation of Pd.
Usually, surface segregation in binary alloys [81], including Pt-Pd alloys, is interpreted as a result of a great positive value of the heat of segregation Q. According to [81], for Pd atoms on the surface of the bulk Pt-Pd alloy Q = +27.7 kJ/mol. In our case, Q is the work necessary to exchange a surface Pd atom and a core Pt atom. Our thermodynamic simulation approach (Section 3) does not involve the heat of segregation; the surface segregation of Pd follows, first of all, from a lower (in comparison with Pt) value of the specific surface free energy σ (0) of Pd: the difference σ Pd between these quantities affects the segregation coefficient (A4). The depletion effect [19], i.e., the effect of the limited source for the segregating component and the effect of the non-zero excess Gibbs energy of mixing ∆G m , result in noticeable diminishing of χ (s) Pd . However, in our atomistic and thermodynamic simulations, the surface segregation of Pd was never inverted into the surface segregation of Pt.
As for theoretical predictions of segregation in binary Pt-Pd NPs, papers [53,82] are of special interest: the authors discussed a number of important and interesting aspects of predicting the properties of Pt-Pd bimetallic nanocatalysts and processes in binary metallic NPs, including the surface segregation of one of the components. In particular, in ref. [82], the next two segregation rules were formulated, which are partially alternative to the segregation factors [19] listed in the Introduction: (i) If the bulk melting temperature of component A is larger than that of component B, then element A will segregate to the surface; (ii) If the bulk melting temperatures of both components have more or less the same magnitude, then the segregated element will be determined by a lower surface energy.
In [53], the authors additionally noted that the first rule should be used if the melting temperatures of components differ by more than 10%. For the Pt-Pd pair, this criterion is equal to 11-12%, depending on the choice of metal in the denominator of the percentage. However, the authors of paper [53] decided that for the Pt-Pd alloy, the second rule should be employed, although such a choice contradicts their above criterion. It is also worth mentioning that the second rule [53,82] coincides with our conclusion drawn in [70,71] about the σ (0) difference as the main factor affecting surface segregation in binary A-B nanoalloys.
In [53], an intriguing but highly controversial (as we believe) conclusion was also drawn: at a certain characteristic value, α, of the NP diameter, d, an inversion should take place of the surface segregation. In other words, according to [53], at d < α, Pt atoms should segregate to the surface of Pt-Pd NPs instead of Pd NPs. In [53], parameter α was referred to as the form factor, since its value depends on the NP shape. The very conclusion about the existence of such a critical NP size was made on the basis of the next relationship: where Q is the heat (energy) of the segregation of Pd in NPs of diameter d, and Q (∞) is the heat of segregation for the surface of the corresponding bulk alloy. According to (3), Q > 0 for d > α and Q < 0 for d < α. Thus, the negative sign of parameter Q was interpreted in [53] as the reason of inversion in surface segregation. In turn, Equation (3) was written in [53] as a consequence of the Thomson formula for the size dependence of the melting temperature, T m , of NPs. Derivation of the Thomson formula is discussed in our paper [83]. A more general Thomson-type equation can be written as follows [84,85]: Equation (3) was written in [53] by assuming that: However, for d < α, Equation (4) predicts negative values for the melting point of NPs, which is physically inadequate. Thus, in [53], the conclusion about the inversion of the surface segregation in binary Pt-Pd NPs was based on the erroneous assumption (5).
In ref. [54], another scenario of surface segregation inversion in binary Pt-Pd NPs was predicted from the results of the MD simulations. The authors of [54] used not EAM, but the Sutton-Chen potential [57]. According to [54], at T = 353.15 K, there was practically no segregation in NPs with a size of 2 nm (χ At the same time, we have no reason to completely deny the possibility of the surface segregation inversion in binary metallic NPs. At least in principle, one should admit that some segregation effects cannot be predicted by only employing atomistic and thermodynamic simulations. Maybe, some specific effects can be predicted by ab initio calculations or discovered experimentally. From this methodological point of view, paper [86] is of interest where ab initio simulation results predicted that small additives of carbon can switch the surface segregation of Cu into the surface segregation of Ni in some Cu-Ni slabs. However, in [86] the third component was put into play.
The problem of the size dependence of the surface segregation in binary nanoalloys is also no less controversial than the problem of the surface segregation inversion. Presumably, there are no experimental data on this topic on the Pt-Pd nanosystem. However, for binary NPs of other fcc metals, e.g., for Au-Pt nanoalloys, the results of atomistic Monte Carlo simulation [43] predict a decrease in the surface segregation with the NP size.
The authors of paper [53] theoretically predicted an opposite effect, i.e., the effect of the surface segregation increasing with the decreasing NP size (see Figure 7 in [53]). However, this conclusion is inconsistent with the following known Langmuir-McLean-type relation [81]: from which the authors of [53] proceeded by quite reasonably replacing χ in Formula (6), proposed for the surface of the bulk alloy, by χ (c) : Here, Q * = Q/RT is the reduced heat of segregation. If, following [53], Formula (3) is also involved, then the heat of segregation Q and, consequently, the surface segregation itself will decrease with the particle size; at d → α , χ (s) will tend to χ (c) , i.e., the segregation effect will completely disappear. Our MD results and the results of our thermodynamic simulations based on the discussed solutions of the Butler equation demonstrate a gradual decrease in the surface segregation with the particle size. However, the results of atomistic simulations (Figures 5 and 6) are less unambiguous: they show only a tendency towards a decrease in segregation with the particle size. Indeed, according to these figures, the minimum value of Pd also corresponds to the minimum NP size, i.e., to NPs containing 140 atoms. However, Obviously, such a partial disagreement of our MD results with the results of the thermodynamic simulations can be explained by the fact that, despite the MD experiments combining successive stages of melting, cooling (quenching), and annealing, the final state of the largest NPs is not completely in equilibrium. This explanation is consistent with a conclusion made in paper [87] devoted to MD simulations of the sintering of small Au nanoclusters beginning with NPs containing only 147 atoms. Preliminarily, NPs (used as initial configurations in the subsequent simulation of sintering) were relaxed (annealed) for approximately 1 µs at temperatures of 400 and 500 K. The results of the relaxation of various structural isomers are interesting and, in our opinion, unexpected: at T = 500 K, the same final equilibrium state of NPs was achieved. However, at T = 400 K, even a long MD experimental annealing time was not sufficient to reach completely equilibrium states.
In accordance with both thermodynamic simulation algorithms (based on the Butler and Langmuir-McLean equations), the main reason for the decrease in the surface segregation with the NP size is the depletion effect [19]. When the size of a binary NP decreases, the stock of the segregating component in its core is depleted. The unlimited source for the segregating component model becomes more and more inadequate. Thus, the depletion effect should be taken into account in any version of the thermodynamic simulation by the second equation of system (A6), i.e., by the mass balance equation for the segregating component. Sections 2 and 3 present the results related to a sufficiently wide range of the Pt-Pd NP sizes (1.6-6.5 nm), but to the same temperature of 300 K, which is low for metallic NPs. Figure 9 demonstrates the results of our atomistic and thermodynamic simulations for Pt-Pd NPs of the same size (d = 5 nm, N = 5000), but annealed at two different temperatures: 300 and 1000 K. The results shown in Figure 9 were obtained using our EAM parameterizations [49] (the results obtained using parameterizations [58] did not noticeably differ from those presented in Figure 9). One can see that the MD results are in good agreement with our most precise results of the thermodynamic simulation corresponding to the most adequate model of the limited source for the segregating component and to approximation, taking into account the excess Gibbs energy of mixing. Figure 9 shows that both alternative simulation approaches (atomistic and thermodynamic) predict a decrease in the surface segregation of Pd with increasing temperature.
Our conclusion about decreasing the surface segregation in binary NPs with increasing temperature is consistent with the available experimental data and theoretical predictions, although just for Pt-Pd NPs there are no experimental data or atomistic simulation results from other authors. Theoretically, for a thin layer of a binary alloy, the effect of a decrease in the surface segregation with a decreasing layer thickness and increasing temperature was predicted in [88]. In addition, for the bulk binary alloys and binary metallic NPs, there are some experimental data and theoretical estimations confirming the effect of decreasing equilibrium surface segregation with increasing temperature. For example, a decrease in the segregation of Cu on the surface of the bulk Cu-Ni alloy was experimentally observed in [89], and on the surface of a thin Cu-Ni film in [90]. According to [43], the atomistic Monte Carlo simulation of binary Au-Pt NPs also revealed (at T > 800 K), diminishing the surface segregation with increasing temperature.
In the frames of our approach to the thermodynamic simulation, the conclusion regarding the decrease in surface segregation with increasing temperature is quite obvious and unambiguous. Indeed, the segregation coefficient (A4) decreases with increasing temperature. At T → ∞ , the segregation coefficient tends to unity, and χ (s) B tends to χ B . perature was predicted in [88]. In addition, for the bulk binary alloys and binary metallic NPs, there are some experimental data and theoretical estimations confirming the effect of decreasing equilibrium surface segregation with increasing temperature. For example, a decrease in the segregation of Cu on the surface of the bulk Cu-Ni alloy was experimentally observed in [89], and on the surface of a thin Cu-Ni film in [90]. According to [43], the atomistic Monte Carlo simulation of binary Au-Pt NPs also revealed (at > 800 K), diminishing the surface segregation with increasing temperature. = 300 K, model of the limited source, and taking into account the non-zero excess Gibbs energy; curve 5 corresponds to the results for = 1000 K (the same most exact model and non-zero excess Gibbs energy); the dashed straight line (6) corresponds to the limiting case of no segregation. Dots present MD results averaged over 5 independent MD experiments and obtained for temperatures = 300 K (•) and = 1000 K (■) employing our EAM parameterizations (data from [49]).
In the frames of our approach to the thermodynamic simulation, the conclusion regarding the decrease in surface segregation with increasing temperature is quite obvious and unambiguous. Indeed, the segregation coefficient (A4) decreases with increasing temperature. At → ∞, the segregation coefficient tends to unity, and ( ) tends to ̅̅̅.

Conclusions
In accordance with our results presented in this paper, two different approaches based on MD and thermodynamic simulations predict the surface segregation of Pd in Pt-Pd NPs of all selected sizes and compositions. Both methods predict a decrease in the equilibrium surface segregation with the NP size. At the same time, we did not observe any inversion of the surface segregation at a small characteristic value of the NP diameter. We also found that the surface segregation decreased with increasing temperature.

Conclusions
In accordance with our results presented in this paper, two different approaches based on MD and thermodynamic simulations predict the surface segregation of Pd in Pt-Pd NPs of all selected sizes and compositions. Both methods predict a decrease in the equilibrium surface segregation with the NP size. At the same time, we did not observe any inversion of the surface segregation at a small characteristic value of the NP diameter. We also found that the surface segregation decreased with increasing temperature. Employing two EAM parameterizations [49,58], we established that the ratio between the binding energies (which is the same for parameterizations [49,58]) is a much more significant factor than the ratio between the melting temperatures of components. In turn, the ratio between the binding energies should qualitatively be the same as between the surface energies of components.
Exemplifying on Pt-Pd NPs, we have confirmed a conclusion made previously, that the difference between the surface energies (surface tensions) of components is the main factor leading to surface segregation [71]. This conclusion coincides with the second segregation rule formulated in [53,82]. At the same time, we believe that the first segregation rule [82], according to which the component with the higher melting point should segregate to the surface, is not obeyed. Obviously, the second rule [82] relates to the equilibrium surface segregation. However, if one of the components of a binary NP or of a bulk body has a lower surface energy but a higher melting temperature, the melting temperature will play the role of the main kinetic factor damping the thermodynamic trend to the surface segregation of this component. In ref. [72], we proposed to distinguish between the thermodynamic and kinetic stability of the core-shell bimetallic nanostructures, which correspond to the limiting case of the pronounced surface segregation. One cannot be entirely sure that the final state of solid NPs in MD experiments is completely in equilibrium. However, good agreement between our MD and thermodynamic results indicate that the approach used in our MD experiments, combining the preliminary melting, solidification, and subsequent annealing of NPs, ensures their final equilibrium states, or at least states close to equilibrium.
The thermodynamic simulation results relate to equilibrium segregation which, in some cases, cannot be achieved in real systems, especially at low temperatures. The accuracy and reliability of thermodynamic results seem to be lower for small nanoclusters (d < 1-2 nm) as, in particular, the employed thermodynamic simulation approach does not take into account the size dependence of the surface energies of components.
We believe that our paper will serve as a stimulus for further experimental and theoretical studies of regularities and mechanisms of segregation not only in Pt-Pd NPs, but also in other binary metallic alloys, including nanoalloys, which are of interest from a scientific point of view, and for numerous potential practical applications.

Data Availability Statement:
The main simulation results were obtained by employing the open access program LAMMPS. The authors are not able now to freely share their own additional software for processing MD results (evaluation of the radial distributions of the local molar fractions of the NP components), as it is a proprietary of their organization (Tver State University). However, after the paper is published, the authors will be able to provide the trial software to interested readers. In addition, we propose to make the toolkit available on the GitHub website after the developed software has been registered.

Appendix A. Algorithms of the Thermodynamic Simulation of the Surface Segregation in Binary NPs Based on Solving the Butler Equation
In the frames of the model of the infinite reservoir for the segregating component, we assumed that χ can be rewritten as follows: We have also taken into account that, according to [73,74], the equality ω i ≈ ω The coefficient: is referred to as the segregation coefficient for component B [73,74]. The model of the infinite reservoir for the segregating component is relevant for sufficiently large NPs, and most adequate for the surface boundary of the bulk alloy. In turn, for NP sizes smaller than 10 nm, the particle core cannot be quite adequately interpreted as an infinite reservoir for the segregating component. Therefore, the average molar fractions: for the particle core and its shell, respectively, where ν A and ν B are the numbers of moles of components A and B, ν = ν A + ν B is the total number of moles, N A and N B are the numbers of A and B atoms, respectively, and N = N A + N B is the total number of atoms in the particle.
As already noted, our approach to the thermodynamic simulation does not put any restrictions on the NP shell thickness. However, the fraction of the surface atoms ξ = N (s) /N should be preassigned as one of the input parameters. In the surface monolayer approximation, the ξ parameter is approximately equal to the reduced (dimensionless) degree of dispersity D * : The degree of dispersity, D, is defined as the surface-to-volume ratio for all particles of a disperse system. For the monodisperse system, D = 3/r 0 . In (A5), n = N/(4/3)πr 3 0 is the average density (concentration) of atoms in the particle, n (s) = N (s) /4πr 2 0 is the average surface density, and r * 0 = r 0 /a, where the atomic diameter, a, can be interpreted as the average value of the atomic diameters a A and a B of the components. The corresponding reduced densities n * and n (s) * are defined as na 3 and n (s) a 2 , respectively. A loose but adequate approximation ξ D * corresponds to n * n (s) * 1. Of course, Equation (A5) does not exclude some further exact approximations for the ξ parameter.
Evaluation of χ The solution in question reduces to solving a quadratic algebraic equation for χ (c) B . Metallic NPs at moderate temperatures are practically non-volatile. Thus, the average molar fraction χ B of component B in the particle will not depend on time, i.e., the condition χ B = const should be fulfilled. The second equation of set (A6) has a simple physical meaning: it is the mass balance equation for the B component in the modelling particle. In the present paper, subscript B is used for Pd atoms and subscript A is used for Pt atoms.
When the ideal solution approximation is not employed, the partial molar excess Gibbs energies, ∆G m,A and ∆G m,B , can be found via the integral excess Gibbs energy ∆G m : In turn, we calculated the integral excess Gibbs energy, ∆G m , using the Redlich-Kister equation [91] where A 1 , A 2 , A 3 , etc., are the constants of expansion into the (χ A − χ B ) powers.
The Redlich-Kister equation can be rewritten in the following form [73]: where L 0 = RTA 1 and L 1 = RTA 2 are the redefined A 1 and A 2 constants. Kaptay [74] introduced the exponential multiplier exp(−T/T 1 ) to provide the tending of ∆G m to zero at high temperatures (T 1 is a very high characteristic temperature of 3000 K). Then, in [92], Kaptay drew the conclusion that the Pt-Pd alloy is adequate to the regular solution model, because the properties of Pt and Pd are very similar. Accordingly, the expansion coefficient, L 1 , can be assumed to be zero, and for the expansion coefficient, L 0 , Kaptay proposed the next equation [92]: L 0 = 24500e (−T/3000) , J/mol.