High-throughput computational-experimental screening protocol for the discovery of bimetallic catalysts

One of the important roles of computer simulations in materials science is to predict and/or design materials prior to experiments. In particular, first-principles calculations using density functional theory (DFT) have played a vital role in the field of catalysis^1,2,3. The computational framework enables the prediction or understanding of reaction pathways on catalyst surfaces, which is undoubtedly useful for experimentalists. However, the accurate estimation of some crucial properties such as the reaction barriers on catalytic surfaces is extremely time-consuming in first-principles calculation schemes. In this regard, obtaining a full understanding of the catalytic reaction mechanism via first-principles calculations is often considered rather inefficient. In our experience, although dependent on computing power, a time cost of several months is required for computing a full energetics picture of even one metallic catalyst system. The efficiency of the search of catalysts via first-principles calculations has therefore been questioned because carrying out experimental tests without the aid of computer simulations can often be even faster.

To maximize the efficiency of catalyst discovery using combined computational and experimental studies, the descriptors or features that bridge these two areas should be wisely chosen. Specifically, a simple but physically reasonable descriptor or feature enabling the representation of catalytic properties is critical. A representative example of such a descriptor is the d-band center theory that correlates between the d-band center⁴, i.e., the average energy of electronic d-states projected onto a surface atom of the catalyst, and the gas adsorption energy, which is widely accepted in the field of metallic catalysis. Coupled with the volcano relationship between catalytic activity and adsorption energy⁵, it is theoretically possible to predict catalytic activity from the d-band center without exploration of the full reaction mechanism. Later, along with the d-band center theory, consideration of the d-band shapes describing higher moments of the d-band as well as the sp-band properties describing local Pauli electronegativity were also suggested to capture the surface reactivity of transition metal alloys^6,7,8. These studies indicate that the electronic density of states (DOS) patterns themselves projected onto the surface atoms of catalysts can serve as an improved descriptor for the development of metal catalysts because they include more comprehensive information on not only d-states but also sp-states (both their values and shapes). Nonetheless, the full DOS patterns themselves have never been used as a descriptor in combined computational-experimental screening processes.

As electronic structures are key to determining the physical/chemical properties of materials, materials with similar electronic structures tend to exhibit similar properties. Indeed, a DFT study revealed that a solid-solution Ir₅₀Au₅₀ alloy can exhibit catalytic properties for the H₂ dissociation reaction with an activity comparable to that of Pt, and this result is attributed to the similar electronic structures (DOS patterns) of the Ir₅₀Au₅₀ alloy and Pt⁹. In another example, the similar electronic structures of Rh₅₀Ag₅₀ and Pd were experimentally verified for hydrogen storage¹⁰. The similar electronic structural characteristics enable the Rh₅₀Ag₅₀ alloy to exhibit hydrogen storage properties superior to those of Pd¹¹, although neither Rh nor Ag alone shows hydrogen storage properties⁹. Thus, bimetallic alloys with a DOS pattern similar to that of Pd, which is a well-known representative metal catalyst, should exhibit catalytic performance similar to that of Pd, which is a good hypothesis for the high-throughput development of bimetallic catalysts.

Here, we used the full DOS pattern as a key descriptor in high-throughput computational-experimental screening protocols. As a demonstration in the high-throughput search of bimetallic catalysts, we selected palladium (Pd), the prototypical catalyst for hydrogen peroxide (H₂O₂) synthesis from hydrogen (H₂) and oxygen (O₂) gases^12,13,14, as our reference material. Using DFT calculations, we first screened 4350 crystal structures of bimetallic alloys to investigate their thermodynamic stabilities (only closest-packed surfaces were considered). Next, the similarities in the DOS patterns between the alloys and Pd were quantified, and subsequently, their synthetic feasibility was evaluated. Through these processes, we finally proposed eight candidates with high DOS similarities with Pd, which are expected to show catalytic properties comparable to those of Pd. These eight screened alloys were experimentally synthesized and tested for H₂O₂ direct synthesis, and four of them (Ni₆₁Pt₃₉, Au₅₁Pd₄₉, Pt₅₂Pd₄₈, and Pd₅₂Ni₄₈) indeed exhibited catalytic properties comparable to those of Pd. In particular, Pd-free Ni₆₁Pt₃₉ is greatly outperformed prototypic Pd with a 9.5-fold enhancement in cost-normalized productivity (CNP) owing to the high content of inexpensive Ni.

High-throughput computational screening for bimetallic catalysts

Figure 1 shows a schematic diagram of our high-throughput screening protocol for the discovery of bimetallic catalysts. Based on 30 transition metals in periods IV, V, and VI, we considered a total of 435 (₃₀C₂) binary systems with a 1:1 (50:50) composition. For each alloy combination, 10 ordered phases that are available for 1:1 composition were investigated, namely, B1, B2, B3, B4, B11, B19, B27, B33, L1₀, and L1₁, leading to a screening of 4350 (435 × 10) crystal structures. Using DFT calculations, we calculated the formation energy (∆E_f) of each phase and determined which crystal structure is most stable in the bulk phase. Here, negative formation energy implies that the phase is thermodynamically favorable (or that the element combination is miscible), and vice versa. Because nonequilibrium alloyed phases can be stabilized by the nanosize effect, although the elemental combinations are immiscible in the bulk phases^11,15,16, a margin of ∆E_f < 0.1 eV was considered when screening the thermodynamic stabilities of the bimetallic systems. The alloyed structures with ∆E_f > 0.1 eV could be readily synthesized; however, such structures would be easily transformed into phase-separated structures during chemical reactions (two-phase structures of the pure elements), which is not suitable for their practical uses as catalysts. Eventually, a total of 435 binary systems and 249 alloys were filtered at the thermodynamic screening step.

There are three screening steps to identify catalyst candidates from 4350 candidates. The top (cyan color) section indicates the first step showing 30 elements and 10 crystal structures to collect binary alloys with a 1:1 composition. The middle (magenta color) and bottom (violet color) sections indicate the second and the last steps relevant to the thermodynamic stability and the DOS similarity, respectively.

Full size image

For the 249 thermodynamically screened alloys, we calculated the DOS pattern projected on the close-packed surface for each structure and compared it with that of the Pd(111) surface. To quantitatively compare the similarity of two DOS patterns, we defined the following measure:

Here, \({{{\mathrm{{\Delta}}} DOS}}_{2 - 1}\) in Eq. (1) indicates the similarity between the DOS of the alloy (DOS₂) and the DOS of the reference Pd(111) (DOS₁) in which a Gaussian distribution function, g(E; σ) in Eq. (2) is considered to compare the two DOS patterns near Fermi energy (E_F) with the high weight. We set the standard deviation, σ = 7 eV because most of the d-band centers for the bimetallic alloys are distributed from −3.5 eV to 0 eV relative to the Fermi energy.

In comparing two DOS patterns, both d-states and sp-states were considered. To clarify the importance of the inclusion of sp-states, we discuss this effect with an example of O₂ adsorption on Ni₅₀Pt₅₀(111), which is one of the elementary steps for H₂O₂ direct synthesis. In Fig. 2, the DOS patterns for d-states and sp-states are shown; the patterns are obtained by the sum of the partial DOSs of the surface Ni and Pt atoms on Ni₅₀Pt₅₀(111). Comparing the DOS patterns before/after O₂ adsorption, the change in the d-band DOS is negligible after O₂ adsorption. On the other hand, the sp-band DOS patterns change more smoothly after O₂ adsorption. This result implies that an O₂ molecule interacts more with the sp-bands of the surface Ni and Pt atoms than with the d-bands on Ni₅₀Pt₅₀(111) in O₂ adsorption. The partial DOSs of O₂ before/after O₂ adsorption on Ni₅₀Pt₅₀(111) are also compared in Fig. 2. An O₂ molecule in the gas phase has a pair of half-occupied orbitals (π_p*) in a triplet state, which leads to magnetic properties, resulting from the spin-splitting into up and down states. After O₂ adsorption, the orbitals of the O₂ molecule broaden and the spin-splitting character disappears. Thus, the magnetic moment of the adsorbed O₂ becomes zero (peroxo state), which is supportive evidence that the π_p* states are fully occupied through the interaction between the metal surface atoms and the O₂ molecule. Then, orbital hybridization between the sp-states of the Ni₅₀Pt₅₀ surface and the p-states of the O₂ molecule mainly occurs, therefore, consideration of the DOS patterns including the sp-state is required, at a minimum, for H₂O₂ direct synthesis.

The s-states (light), p-states (dark), and d-states (patterned) DOS plots are shown. Occupied states are filled in color; unoccupied states are depicted by lines. a–c The shape of states without adsorbed O₂ on the Ni₅₀Pt₅₀ (111) surface and O₂ gas phase in a vacuum. d–e The shape of states of the Ni₅₀Pt₅₀ (111) surface and O₂ of the adsorbed system.

Full size image

The DOS similarity values of the 249 alloy surfaces are shown in Fig. 3. According to the definition of the DOS similarity, as the value approaches zero, the electronic structure of the alloy becomes more similar to that of Pd(111), and the material is expected to show catalytic properties similar to those of Pd(111). In Fig. 3, we chose 17 candidates with the low DOS similarity values (\({{{\mathrm{{\Delta}}} DOS}}_{2 - 1}\)< 2.0): CrRh (\({{{\mathrm{{\Delta}}} DOS}}_{2 - 1}\)= 1.97, B2 structure), FeCo (1.63, B2), CoNi (1.71, L1₀), CoHf (1.89, L1₀), CoTa (1.96, L1₀), CoPt (1.78, L1₀), NiCu (1.11, L1₀), NiNb (1.96, L1₀), NiMo (1.87, L1₀), NiRh (0.98, L1₀), Ni-Pd (0.84, L1₀), NiIr (1.28, L1₀), Ni-Pt (1.16, L1₀), CuPd (1.51, B2), RhPt (1.67, L1₀), Pd-Pt (1.16, L1₀), and PdAu (1.43, L1₀). Here, owing to the rarity of Hf and Ta, the CoHf and CoTa alloys were excluded from the candidate list. In addition, we excluded NiCu and CrRh from the list because the large reduction potential difference between the elements in the binary system leads to the unfavorable formation of alloyed structures in which the two elements are homogeneously mixed via a wet-chemical method for nanoparticle (NP) synthesis¹⁷.

The diagonal components denote the most thermodynamically stable crystal structures (e.g., FCC, BCC, HCP, and SC) of pure metals. The occupied components under the diagonal denote the most thermodynamically stable crystal structures (e.g., B2, B11, B19, B27, B33, L1₀, and L1₁) of bimetallic alloys. The background colors above the diagonal describe the formation energies, and the occupied components above the diagonal denote the DOS similarity of each alloy compared with Pd(111).

Full size image

Experimental verification of the DFT-guided bimetallic catalysts

For the remaining 13 candidates, we tried to synthesize their alloyed NPs and confirm whether their crystal structures were identical to those (Fig. 3) predicted by DFT calculations. For NP synthesis, we used a butyllithium reduction method¹⁸. Details of the NP synthesis are found in the Methods section. Because it was assumed in Fig. 3 that the proposed alloys had a composition of 50:50, we intended to synthesize alloys with nearly the same composition. Inductively coupled plasma (ICP) analyses revealed that the compositions of the synthesized alloys (loading ratio of metal precursors = 50:50) were Au₅₁Pd₄₉, Ni₆₁Pt₃₉, Pd₅₀Cu₅₀, Pd₅₂Ni₄₈, Pt₅₈Co₄₂, Pt₅₂Pd₄₈, Rh₅₆Ni₄₄, and Rh₅₆Pt₄₄; however, CoFe, NiIr, NiMo, NiNb, and NiCo were not alloyed because the reducing power of butyllithium used as a reducing agent was likely insufficient. To confirm whether the alloyed NPs had the same crystal structures as predicted by DFT, their X-ray powder diffraction (XRD) patterns were investigated, as shown in Fig. 4. The samples showed alloyed diffraction patterns without signals from phases of the pure elements that compose each alloy. We also simulated XRD patterns of the alloys on the basis of the DFT crystal structures and found that the simulated patterns matched well with the experimental patterns. This result clearly reveals that we successfully synthesized alloyed NPs with the crystal structures predicted by DFT calculations for eight bimetallic systems. The lattice parameters for the alloyed samples are presented in Supplementary Table 1. In addition, we performed Rietveld refinements of the XRD patterns (Supplementary Figure 1). The synthesized NPs clearly have cubic-based structures as predicted by our DFT calculations, although in the Ni-Pt samples a small amount of non-alloyed Pt (pure Pt) coexists with alloyed Ni₇₂Pt₂₈. As we used the synthesized Ni-Pt samples themselves when measuring their catalytic performance, we used the ICP composition (Ni₆₁Pt₃₉).

XRDs for a Au₅₁Pd₄₉, b Ni₆₁Pt₃₉, c Pd₅₀Cu₅₀, d Pd₅₂Ni₄₈, e Pt₅₈Co₄₂, f Pt₅₂Pd₄₈, g Rh₅₆Ni₄₄, h Rh₅₆Pt₄₄ NP samples. For comparison, the red and green lines in the XRD patterns correspond to the peaks of pure elements in each bimetallic NP. The blue peaks were simulated with the crystal structures of the bimetallic system predicted by DFT calculations. Above each XRD pattern, the element mapping images acquired by STEM-EDX are included, where the first image (the most left) is the overlay image of the maps for elements in each bimetallic NP. The scale bars correspond to 20 nm. The enlarged images of the STEM-EDX are also provided in Supplementary Figures 2~9.

Full size image

It was also important to determine whether the synthesized NPs have ordered or disordered structures because the two types of structures can show different catalytic properties¹⁹. The main role of our descriptor (ΔDOS₂₋₁) in the high-throughput protocol for catalyst development is to discover potential catalysts and then transfer the information to experimentalists at high speed. For this purpose, a reference structure for comparison with alloys is required because various structures are possible for random solid-solution structures. In this work, ordered structures were used as a reference and were also validated by comparing simulated and experimental XRD patterns. From this process, we intended to predict the possibilities and tendencies of alloys for use as H₂O₂ direct synthesis catalysts. Nevertheless, because the synthesized sample could include random solid-solution structures, there is the potential for discrepancies between theory and experiment.

We further characterized the NP samples by high-angle annular dark-field scanning transmission electron microscopy (Fig. 4 and Supplementary Figure 2~9). The elemental color mapping and overlay clearly indicate the homogeneous mixing of elements in each NP. The mean diameters of the alloyed NPs were 3.4 ± 0.4 nm (Pd₅₀Cu₅₀), 8.8 ± 2.0 nm (Ni₆₁Pt₃₉), 5.0 ± 1.5 nm (Rh₅₆Ni₄₄), 11.6 ± 2.0 nm (Au₅₁Pd₄₉), 5.2 ± 1.5 nm (Pd₅₂Ni₄₈), 8.1 ± 2.2 nm (Pt₅₂Pd₄₈), 5.8 ± 1.0 nm (Pt₅₈Co₄₂), and 1.7 ± 0.2 nm (Rh₅₆Pt₄₄). In particular, for the four alloys (Ni₆₁Pt₃₉, Au₅₁Pd₄₉, Pd₅₂Ni₄₈, and Pt₅₂Pd₄₈) showing catalytic properties comparable to those of Pd (Table 1), their NP diameters were in the range of 5.2~11.6 nm, which is similar to the diameter (7.2 ± 1.6 nm) of Pd NPs (Supplementary Figure 10).

To investigate whether the catalytic performance of the eight proposed alloys is comparable to that of Pd, the alloyed catalysts were tested for direct H₂O₂ synthesis under mild conditions (20°C, 1 atm) in comparison to pure Pd (Table 1). Along with the industrial importance of H₂O₂ as an eco-friendly oxidizing agent^20,21, the H₂O₂ direct synthesis reaction can be performed under ambient temperature and pressure conditions, and Pd has been regarded as the archetypical catalyst for the chemical reaction although there is still room for performance enhancement²². In terms of H₂O₂ productivity representing the catalytic activity for H₂O₂ direct synthesis, three alloy systems, Au₅₁Pd₄₉ (505.2 mmol_H2O2 gmetal⁻¹ h⁻¹), Ni₆₁Pt₃₉ (880.1), and Pt₅₂Pd₄₈ (1225.6) exhibited performance comparable to that of Pd (552.4). In addition, Pd₅₂Ni₄₈ showed noticeable H₂O₂ productivity (218.1), although this productivity is lower than that of Pd. Here, it is notable that the four alloy systems have lower \({{{\mathrm{{\Delta}}} DOS}}_{2 - 1}\) values than the alloy systems showing poor catalytic activities. Although the Rh-Ni system has a low \({{{\mathrm{{\Delta}}} DOS}}_{2 - 1}\) value, the DOS patterns near E_F are quite different from those of Pd (Supplementary Figure 11), which can explain why it shows poor catalytic activity.

On the other hand, the orders of the \({{{\mathrm{{\Delta}}} DOS}}_{2 - 1}\) values and H₂ conversion do not show perfect consistency. However, it is not enough to represent the catalytic performance for H₂O₂ direct synthesis only with H₂ conversion. The H₂O₂ synthesis reaction involves several elementary steps on the catalyst, such as H₂ dissociation and the first and second hydrogenation of O₂ generating OOH and H₂O₂, respectively. Moreover, when assessing the catalytic performance, side reactions such as H₂O formation should be considered. In other words, the catalytic performance of H₂O₂ direct synthesis is not dominated by only one catalytic behavior but coupled with multiple behaviors that involve not only H₂ conversion but also H₂O₂ productivity and selectivity. In addition, we compared our descriptor of DOS similarity (ΔDOS₂₋₁) with the d-band center that has been widely used in metallic catalysts. Here, we considered the difference between the d-band centers of the alloys and that of Pd (|ε_d,Pd−ε_d|); the d-band centers of the bimetallic alloys are shown in Supplementary Table 4. For comparison, Pt₅₂Pd₄₈ and Pd₅₀Cu₅₀ were chosen because they show relatively high and low H₂ conversions, respectively (Table 1). The ΔDOS₂₋₁ values appropriately represent the order of H₂ conversion (Pt₅₂Pd₄₈: ΔDOS₂₋₁ = 1.16 and Pd₅₀Cu₅₀: ΔDOS₂₋₁ = 1.51), in contrast to the difference in d-band centers (Pt₅₂Pd₄₈: 0.31 and Pd₅₀Cu₅₀: 0.08). In this regard, the similarity in DOS patterns can be a good descriptor for the discovery of bimetallic catalysts and needs no calculations relevant to the reaction pathways on the catalyst surfaces, which would require a high computational cost.

Among the four alloy combinations showing high H₂O₂ productivity, the Au-Pd^{13,14,23,24,25}, Ni-Pd²⁶, and Pd-Pt²⁷ systems were previously reported as promising candidates for H₂O₂ direct synthesis, validating our approach once again. In addition, our approach led to the discovery of a catalyst, a Ni-Pt alloy, which has not yet been reported for H₂O₂ direct synthesis, although the system is well known as an oxygen reduction reaction (ORR) catalyst in fuel cells^28,29,30. As the price of Pd rapidly increases³¹, our Pd-free Ni-Pt alloy can serve as a more cost-effective option. To identify the cost-effectiveness of the Ni-Pt catalyst, we list the CNP (unit, mmol_H2O2 $metal⁻¹ h⁻¹)³² of the assessed materials in Table 1, in which the metal price [unit, U.S. dollar ($)] instead of catalyst mass is used. In terms of this metric, N_i61Pt₃₉ shows a CNP of 69.3, which is not only ~9.5 times higher than that of Pd (7.3) but also much higher than those of the other candidates (Au₅₁Pd₄₉: 7.1, Pd₅₂Ni₄₈: 5.5, and Pt₅₂Pd₄₈: 23.0). As another metric for catalytic activity, we also considered the turnover frequency of the catalysts. In terms of this metric, Ni₆₁Pt₃₉ (98.5) shows a performance superior to that of Pd (58.8). Moreover, for efficient H₂O₂ direct synthesis, the H₂ conversion and H₂O₂ selectivity properties of the catalysts should be considered. It is also interesting that Ni₆₁Pt₃₉ shows an H₂O₂ selectivity similar to that of Pd, although Ni₆₁Pt₃₉ shows higher H₂ conversion. These findings clearly demonstrate that the Ni-Pt bimetallic catalyst is very attractive for H₂O₂ direct synthesis.

When screening bimetallic catalysts by electronic structure calculations, only one composition (50:50) for the bimetallic systems was considered because DFT prediction to determine the optimized composition to provide the best catalytic performance would be challenging and time-consuming. Instead, composition optimization proceeded with the help of experiments, likely providing a faster result than a DFT prediction. Therefore, to optimize the composition of the Ni-Pt bimetallic catalysts, we also explored the effects of the Ni/Pt ratio in Ni_xPt_100−x NPs on their catalytic properties for H₂O₂ direct synthesis (Supplementary Table 2), and no significant structural changes in the Ni_xPt_100−x NPs were observed (Supplementary Figure 12 and 13). We found that the best catalytic performance was exhibited by Ni₆₁Pt₃₉ in terms of H₂O₂ productivity, H₂ conversion, and H₂O₂ selectivity.

Catalytic reaction mechanism

The Ni-Pt bimetallic catalysts were proposed based only on electronic structure calculations without any calculations relevant to the reaction pathways. Thus, to justify the experimental catalytic activity shown in Table 1, we investigated the reaction mechanisms on the catalyst surface using DFT calculations. We considered the close-packed (111) surface of the Ni₅₀Pt₅₀ slab in L1₀, where the alloy composition was considered for simplicity of the calculations although the best catalytic performances were experimentally observed for the Ni₆₁Pt₃₉ composition. The Langmuir-Hinshelwood mechanism^33,34,35 involving adsorption of both H₂ and O₂ over a catalyst surface was first explored. Our DFT calculations revealed that dissociative adsorption of H₂ occurs very favorably over Ni₅₀Pt₅₀(111) surfaces (Supplementary Figure 14), which is similar to Pd(111)^33,34. This finding implies that the hydrogenation of O₂ over Pd(111) and Ni₅₀Pt₅₀(111) determines the catalytic activity for H₂O₂ synthesis. In Fig. 5a, the overall H₂O₂ synthesis involving two sequential hydrogenations of O₂ (I→II and II→III) is exothermic by 0.42 eV on Ni₅₀Pt₅₀(111), which is more favorable than the value of 0.07 eV for Pd(111). Moreover, while the rate-determining step on Pd(111) is the first hydrogenation (I→II) to form OOH* (*denotes a surface site) with an energy barrier of 0.78 eV, the rate-determining step on Ni₅₀Pt₅₀(111) is observed at the desorption step of the generated H₂O₂*, with an energy barrier of 0.69 eV (III→V). These results indicate that higher H₂O₂ productivity is possible on Ni₆₁Pt₃₉ NPs than on Pd NPs, as observed in Table 1. On the other hand, Wilson and Flaherty³⁶ proposed another mechanism for H₂O₂ formation on Pd NPs called a heterolytic reaction mechanism, similar to the two-electron ORR, in which an electron is generated during the oxidation process of H₂ on the NP surfaces. We also investigated the mechanism on Ni₅₀Pt₅₀(111) and Pd(111) by DFT calculations where the computational hydrogen electrode method was used along with consideration of implicit water solvents (Fig. 5b). On both Pd(111) and Ni₅₀Pt₅₀(111), the ORR processes for H₂O₂ formation are thermodynamically favorable. In particular, although the first protonation is very similar (−1.03 eV for Pd and −1.00 eV for Ni₅₀Pt₅₀), the second protonation is slightly more favorable on Ni₅₀Pt₅₀(111) (−1.40 eV vs −1.58 eV). This mechanism also readily reveals more preferential H₂O₂ production on Ni-Pt surfaces than on Pd surfaces.

a The Langmuir-Hinshelwood mechanism. b The electron–proton transfer mechanism. The cyan lines are the energy diagram for Ni₅₀Pt₅₀(111), and the magenta lines are for Pd(111). The color code for the atoms is as follows: teal = Ni, coral = Pt, red = O, and green = H.

Full size image

In Fig. 5a, it is interesting that Ni₅₀Pt₅₀(111) stabilizes H₂O₂* (III state) more than Pd(111) does, whereas the energetics for OOH* (II state) are comparable on the two surfaces. To clarify this result, we investigated the binding characteristics between the adsorbates (OOH* and H₂O₂*) and catalyst surfaces (e.g., Pd(111) and Ni₅₀Pt₅₀(111)) via projected crystal orbital Hamilton population (pCOHP) analysis (Supplementary Figure 15a), in which we display the pCOHPs according to the typical method, namely drawing negative (i.e., bonding) contributions to the right and positive (i.e., anti-bonding) contributions to the left. For OOH* binding to Pd(111) and Ni₅₀Pt₅₀(111), the bonding contributions are comparable. In other words, there are three energy windows for the bonding contributions: −3 to −4 eV (1st window), −5 to −7 eV (2nd), and −8 to −9 eV (3rd). Here, the bonding energy levels for the first and third contributions are lower (stronger binding) for Ni₅₀Pt₅₀(111), whereas the bonding energy levels for the second are lower for Pd(111). On the other hand, in the case of H₂O₂* binding, all of the bonding contributions for Ni₅₀Pt₅₀(111) are lower than those for Pd(111), which explains why the H₂O₂* adsorbate is more stable on Ni₅₀Pt₅₀(111) than Pd(111). This behavior can be explained by the back-bonding character between H₂O₂* and Ni₅₀Pt₅₀(111). Alloying of Pt and Ni creates electron-rich Ni atoms in the alloy structure, which leads to back-bonding with the adsorbate (the Ni atom donates electrons to an empty orbital in the adsorbate). This effect is significantly enhanced for H₂O₂* compared to OOH* as the O atom in the adsorbate bonded to metal surfaces is more electron-deficient for H₂O₂* than OOH* due to the additional O-H bond in H₂O₂*. From the charge density difference plot shown in Supplementary Figure 15b, we find that Ni atoms on Ni₅₀Pt₅₀ and O atoms in H₂O₂* donate electrons (cyan isosurface) while the O atoms in H₂O₂* also accept electrons (yellow), indicating the back-bonding mode between the metal and H₂O₂.

In addition, we discuss the selectivity of H₂O₂ production on Pd(111) and Ni₅₀Pt₅₀(111) with DFT calculations. It is well known that H₂O production, the side reaction for H₂O₂ production, occurs by the O-O bond dissociation^27,32,33,34. Therefore, we investigated the possible O-O bond scission reactions on the Pd(111) and Ni₅₀Pt₅₀(111) surfaces by DFT calculations, where we considered the ratio of the O-O bond dissociation barrier at several steps to the corresponding hydrogenation step barrier or desorption barrier (O₂* dissociation vs O₂* hydrogenation, OOH* dissociation vs OOH* hydrogenation, and H₂O₂ dissociation vs H₂O₂* desorption) for H₂O₂ production as an estimation for the H₂O₂ selectivity (Supplementary Table 5). We assume that the H₂O production reaction occurs most frequently at the O-O bond scission step with the lowest reaction barrier among the three possible O-O bond scission reactions. Comparing the ratios for the three O-O bond scission reactions, Pd(111) shows the lowest ratio of 0.14 (unitless) at the H₂O₂* scission step, which is similar to that of Ni₅₀Pt₅₀(111) (0.12 at the OOH* scission step). This is consistent with the experimental finding that the selectivity for H₂O₂ production is similar for both Pd (63.0%) and Ni₅₀Pt₅₀ (61.2%).

This work demonstrates an efficient computational-experimental screening protocol for the development of bimetallic catalysts to replace conventional Pd catalysts for H₂O₂ direct synthesis. Using this protocol, a Pd-free bimetallic catalyst, Ni₆₁Pt₃₉, was successfully developed. Here, we note that Ni and Pt have not usually been considered in the catalysis field for H₂O₂ direct synthesis because they prefer to dissociate O₂ molecules, which prevents H₂O₂ formation, although Ni-Pd²⁵, Pt-Pd²⁶, and Pt-Ag³⁷ bimetallic catalysts have been reported. The success story of our protocol primarily results from the validity of our descriptor, the similarity of electronic DOS patterns, for catalyst screening. Moreover, one can readily extend our screening protocol by updating a reference material. For example, since the Ni-Pt system shows superior catalytic performance in this work, we can discover other catalysts by using the DOS pattern of the Ni-Pt alloy as a reference, which remains as future work. In this regard, this protocol will immediately provide opportunities for catalyst discovery for the replacement or reduction in the use of the platinum-group metals and then enable the further application of these catalysts to other chemical reactions.

Computational details for screening bimetallic catalysts

We generated 10 ordered phases for each bimetallic alloy with a 1:1 composition. Specifically, the crystal structure prototypes were obtained by NaCl (B1), CsCl (B2), zinc blende (B3), wurtzite (B4), TiCu (B11), CdAu (B19), FeB (B27), CrB (B33), AuCu (L1₀), and CuPt (L1₁) compounds from the Material Project³⁸. Here, the initial lattice parameters were replaced by values commensurate with covalent radii of two elements in the target compounds. By means of the DFT optimization process, we collected the most stable bulk phases of 4350 binary alloys. Furthermore, the formation energies (ΔE_f) of all the bulk phases were calculated as follows: ΔE_f = E[AB] – E[A] – E[B], where E[AB], E[A], and E[B] are the total energies of the bimetallic AB compound, A crystal, and B crystal in the unit cell. For thermodynamically favorable alloys, a close-packed surface of each crystal structure was used for the calculation of the DOS pattern of each alloy. In principle, the close-packed surface of L1₀ is the (111) facet, that of B2 is the (110) facet, and those of other crystal structures are the (001) facet³⁹. The DOS patterns were obtained by the sum of the partial DOSs on the surface atoms.

To investigate thermodynamic stabilities in bulk phases and the DOS patterns of surface atoms, spin-polarized DFT calculations were performed using the Vienna Ab Initio Simulation Package (VASP)⁴⁰ with projector-augmented-wave (PAW) pseudopotentials⁴¹ and the revised Perdew-Burke-Ernzerhof (RPBE) exchange-correlation functional⁴². Then, we used a plane-wave kinetic energy cutoff of 400 eV and Monkhorst-Pack k-point meshes of 8 × 8 × 8 for the bulk systems and 4 × 4 × 1 for the slab systems. The bulk of the 1 × 1 × 1 unit cell and the slab of the 2 × 2 unit cell with four layers were used. In the slab systems, a vacuum spacing of 15 Å was used to avoid interactions between slabs. All relaxations were performed until the energy change was >1 × 10⁻⁶ eV/cell and the forces on the individual atoms were <0.01 eV Å⁻¹.

H₂O₂ formation mechanism

The energy diagrams for H₂O₂ production on the Ni₅₀Pt₅₀(111) surface and Pd(111) surfaces were also calculated by the DFT method. The simulations were performed using the VASP⁴⁰ using the PAW⁴¹ pseudopotentials to describe the potential from the ionic core. For the exchange and correlation terms, we employed the RPBE functional⁴² with Grimme’s D3 dispersion correction⁴³. An energy cutoff of 400 eV and Monkhorst-Pack k-point meshes of 3 × 3 × 1 for the slab calculations were used. A large vacuum spacing of 15 Å was used to prevent interslab interactions. The top two layers and the reactant molecules were optimized until the energy change was <1 × 10⁻⁴ eV/cell and the force on each atom was less than 0.02 eV Å⁻¹. In the case of the Langmuir-Hinshelwood mechanism, the climbing image nudged elastic band⁴⁴ method was employed to investigate the transition states for the main reaction. For the electron–proton transfer mechanism, the computational hydrogen electrode method⁴⁵ with U = 0 V in addition to the water solvation energy calculated with the implicit solvent model implemented within the VASP_sol package⁴⁶ was used to construct the energy diagram on the following reactions:

We also performed pCOHP calculations⁴⁷ using the LOBSTER package⁴⁸ to understand the binding character between the adsorbates (e.g., OOH* and H₂O₂*) and catalyst surfaces (e.g., Pd(111) and Ni₅₀Pt₅₀(111)). The differential charge density (Δρ) of the H₂O₂* adsorbates on the Ni₅₀Pt₅₀ (111) surface was evaluated by the following:

where, ρ_tot is the charge density of the H₂O₂*-Ni₅₀Pt₅₀ system and ρ_ads, and ρ_slab are the charge densities of the H₂O₂* adsorbate and metal slab, respectively. The atomic coordinates for ρ_ads, and ρ_slab are equal to those in the H₂O₂*-Ni₅₀Pt₅₀ system.

Chemicals

Palladium(II) acetate [Pd(CH₃COO)₂, ≥99.9%], copper(II) acetate [Cu(CH₃COO)₂, 99.99%], nickel(II) acetylacetonate [Ni(acac)₂, 95%], rhodium(III) nitrate hydrate [Rh(NO₃)₃∙XH₂O, ~36%], gold(III) chloride trihydrate [HAuCl₄∙3H₂O, 99.9%], chloroplatinic acid hydrate [H₂PtCl₆∙XH₂O, ~38%], cobalt acetate tetrahydrate [Co(CH₃COO)∙4H₂O, 99.999%], dioctyl ether, oleylamine, and n-butyllithium solution [2.0 M in cyclohexane] were purchased from Sigma-Aldrich. Ethanol was analytical grade and was used without further purification.

Catalyst synthesis

All alloy metal catalysts and Pd catalysts were synthesized using the butyllithium reduction method. To synthesize the bimetallic XY catalysts, 0.017 mmol of X precursor and Y precursor (details are described in Table 3 of Supporting Information) were dissolved in 10 mL of dioctyl ether with 2 mL of oleylamine at 50°C. The solution was then injected into a butyllithium solution containing 15 mL of dioctyl ether and 1.2 mL of 2.0 M butyllithium in cyclohexane at 50°C. The colloids were stirred (500 rpm) for 20 min and then heated to 120°C for 1.5 h in an Ar atmosphere. The reaction mixture was further heated to 260°C for 1 h. The mixture was cooled to room temperature, and 1.25 mL of trioctylphosphine was injected to protect the colloids. The resulting NPs were washed with ethanol three times.

Catalyst characterization

Transmission electron microscopy (TEM) images and energy-dispersive X-ray spectra (EDS) were acquired using a transmission electron microscope (FEI Talos F200X) equipped with scanning transmission electron microscopy–energy dispersive X-ray spectroscopy (Bruker Super–X EDS system). The crystal structure was examined by X-ray diffraction (Rigaku Dmax 2500) with Cu Kα radiation (λ = 1.5406 Å). The samples for the TEM images and X-ray diffraction patterns were purified by centrifugation three times to remove the surfactants and/or excess reactants. Then, the NPs were dispersed in ethanol. The resulting solution was dropped on a copper grid coated with an amorphous carbon film for the TEM images and a soda-lime-silica glass for the X-ray diffraction pattern. The actual compositions of each alloy catalyst were measured via ICP atomic emission spectrometry using iCAP6500 Duo (Thermo).

Catalyst performance measurements

The direct synthesis of H₂O₂ was performed in a double-jacket glass reactor under stirring at 1200 rpm. The reaction medium (150 ml) was composed of ethanol and water (ethanol:water volume ratio = 1:4) with 0.9 mM sodium bromide (NaBr) and 0.02 M phosphoric acid (H₃PO₄). The catalyst weight used was 1 mg. The volumetric flow rate of the reactant gas stream was 22 ml/min (H₂:O₂ volume ratio = 1:10). The reaction was performed at 293 K and 1 atm for 1 h. After the reaction, the concentration of H₂O₂ was calculated using iodometric titration. The concentration of H₂ was calculated using gas chromatography (Younglin, 6500GC) equipped with a Molecular Sieve 5 A (Supelco, 60/80 mesh, 6 ft × 1/8 in, 2.1 mm) packed column and a thermal conductive detector. The catalytic properties (H₂O₂ productivity, H₂ conversion, and H₂O₂ selectivity) shown in Table 1 were calculated with the following equations: