Potential-pH diagrams considering complex oxide solution phases for understanding aqueous corrosion of multi-principal element alloys

The potential-pH diagram, a graphical representation of the thermodynamically predominant reaction products in aqueous corrosion, is originally proposed for the corrosion of pure metals. The original approach only leads to stoichiometric oxides and hydroxides as the oxidation products. However, numerous experiments show that non-stoichiometric oxide scales are prevalent in the aqueous corrosion of alloys. In the present study, a room temperature potential-pH diagram considering oxide solid solutions, as a generalization of the traditional potential-pH diagram with stoichiometric oxides, is constructed for an FCC single-phase multi-principal element alloy (MPEA) based on the CALculation of PHAse Diagram method. The predominant reaction products, the ions in aqueous solution, and the cation distribution in oxides are predicted. The oxide solid solution is stabilized by the mixing free energy (or mixing entropy) and the stabilizing effect becomes more significant as the temperature increases. Consequently, solid solution oxides are stable in large regions of the potential-pH diagram and the mixing free energy mostly affects the equilibrium composition of the stable oxides, while the shape of stable regions for oxides is mostly determined by the structure of the stable oxides. Agreements are found for Ni2+, Fe2+, and Mn2+ between the atomic emission spectroelectrochemistry measurements and thermodynamic calculations, while deviations exist for Cr3+ and Co2+ possibly due to surface complexation with species such as Cl− and the oxide dissolution. By incorporating the solution models of oxides, the current work presents a general and more accurate way to analyze the reaction products during aqueous corrosion of MPEAs.


INTRODUCTION
Degradation by aqueous corrosion is ubiquitous for materials under various environments, causing significant damage to structural and functional engineered components in most industrial sectors 1,2 . The predominant electrochemical and chemical reactions governing material degradation processes can be represented by the potential-pH diagram, a pioneering approach advanced by Marcel Pourbaix (hence also named Pourbaix diagram) 3,4 . In this diagram, the regions for thermodynamically predominant reaction products (e.g., aqueous complexes, oxides, hydroxides, oxy-hydroxides, etc.) are plotted as a function of the pH of water (or aqueous solution) and electrode potential (E) based on the minimization of the Gibbs energies of formation for all the possible reactions, in order to distinguish the E-pH-T conditions for immunity, active corrosion, and passivity of a given element [2][3][4][5][6] .
Originally, the E-pH diagram was proposed for the predominant reactions between pure metals and an aqueous solution containing dissolved H + , OH − , H 2 , and O 2 3,4 . The solid reaction products were limited to the metal-hydrogen-oxygen (M-H-O) type of compounds, hence referred to as M-H-O type E-pH diagram [2][3][4][5][6] . To construct the E-pH diagrams by calculation, independent reactions involving solid and aqueous products are listed with their reaction thermodynamics, and the equilibrium of reactions correspond to the boundaries between predominant species. Such diagrams cannot be applied to predict the dissolution of an alloy composed of multiple metallic elements 6 . A common way to predict the dissolution and the predominant reaction products of an alloy using such diagrams is to superimpose several independent E-pH diagrams for pure metals. However, this method cannot yield an accurate prediction of the predominant reaction products in such situations 6,7 .
To overcome this obstacle, the E-pH diagrams for binary Fe-Cu 8,9 , Fe-Cr 8,10 , Fe-Ni 7 , and Ni-Cr alloys 7 were constructed focusing on the effect of multi-metallic oxides in comparison with E-pH diagrams for pure metals. It was found that the passivity regions for alloys were larger due to the formation of multimetallic oxides 7 . For example, the NiFe 2 O 4 spinel phase in the corrosion of Ni-Fe alloys cannot be predicted by simply superimposing the E-pH diagrams of Ni and Fe, since NiFe 2 O 4 is not present in either of them 6,7 . Using the thermodynamic data for oxides from binary alloys, the E-pH diagram for Fe-Cr-Ni alloy was constructed, which consequently resembled the combination of E-pH diagrams for binary alloys 7 . Moreover, the heat capacities of ions were estimated by the Criss-Cobble method 11 , which is only applicable to charged aqueous complexes. The E-pH diagram for Fe-Cr-Ni alloy was revised with updated thermodynamic data for spinel oxides (i.e., NiCr 2 O 4 , FeCr 2 O 4 , NiFe 2 O 4 ) 12 and the modified Helgeson-Kirkham-Flowers model was adopted to further incorporate the uncharged aqueous species 13,14 . However, in these E-pH diagrams, the metallic elements were assumed to be in their standard states with unit activities and the thermodynamics resulting from alloying was not considered [7][8][9][10]12 .
The E-pH diagram is exceedingly complex for multicomponent alloys (more than three alloy elements), mostly due to a large number of possible oxides as reaction products. Moreover, in the case of aqueous solutions containing various salts, species, and phases (e.g., sulfur-and chlorine-containing species) not encountered in the M-H-O type E-pH diagrams are expected to be stable under certain E-pH conditions. The lack of thermodynamic data for these possible reaction products has become the major obstacle in constructing the E-pH diagram for complex systems. With the advent and popularization of first-principles calculations based on the density functional theory (DFT) in materials research, the free energy of the predominant species can be calculated given their crystal structures without experimental input 15,16 . However, the accuracy of this approach in determining phase boundary is inevitably limited by the inherent uncertainty in energies from first-principles calculations, i.e., up to ±0.5 eV per atom for certain oxides 15 . Moreover, the difficulty arises when attempting to model the aqueous solution or disordered non-stoichiometric solid phases in DFT, which requires exceptionally large simulation cells and poses computational challenges 17 .
In recent decades, thermodynamic databases for multicomponent systems based on the methodology called CALculation of PHAse Diagram (CALPHAD) has been widely applied to materials design. In this method, a thermodynamic model is fitted to the thermochemical data for a phase or equilibrium for phase transformations from experiments. This enables optimization of the corresponding thermodynamic parameters so that an analytical expression of the Gibbs energy as a function of temperature and composition can be obtained 18 . This is a very active field of materials research and has been successfully applied to the design of various materials 19 . Employing the Redlich-Kister polynomials 18 , the excess Gibbs energy due to the non-ideal thermodynamic interactions between components in complex systems can be efficiently represented without significantly increasing the number of fitting parameters 18 . Moreover, thermodynamic parameters for sub-systems can be directly adopted to describe the thermodynamics of higher-order systems, making it convenient to further develop the database for multicomponent systems (e.g., commercial alloys) 18,19 . The Gibbs energy of formation can be obtained with CALPHAD-type thermodynamic databases for chemical-free energy and the Nernst potential for the energy change associated with charge transfer 6,7 for a large number of phases in multicomponent systems. This significantly increases the capability for the analysis of complex systems, such as in the case of aqueous corrosion of commercial alloys 20,21 . Such approaches have been adopted to analyze the aqueous corrosion of stainless steels 20-24 , SAM40 amorphous steel coatings 25,26 , Nibased superalloys 1,27 , and multi-principal element alloys (MPEAs, or high entropy alloys, HEAs) 28,29 .
The previous E-pH diagrams for alloys are calculated considering only stoichiometric oxides along with other compounds, mostly due to the complexity in thermodynamics of nonstoichiometric oxides. However, non-stoichiometric oxides are increasingly prevalent in the aqueous corrosion of many alloys as the chemistry of the metallic materials becomes more complicated. Non-stoichiometric corundum (Fe,Cr) 2 O 3 and spinel (Fe,Cr, Ni) 3 O 4 were characterized in the oxide scales by transmission electron microscopy (TEM) instruments equipped with energydispersive spectroscopy (EDS) for 316 and 304 stainless steels in pressurized water at high temperatures 30 30,31 . Although attempts were made to model the nonstoichiometric spinel in the E-pH diagram, the non-stoichiometric oxides were treated as a mechanical mixture of the stoichiometric endmembers, therefore, the solid solution nature of oxides was not considered 32,33 . In the corrosion of Ni-22%Cr-6%Mo (wt. %) alloy in aqueous solution at pH 4 with K 2 S 2 O 8 and Na 2 SO 4 , nonstoichiometric halite (Ni,Cr)O was observed, which is attributed to the kinetic "solute capture" effect due to the fast passivation (compared with the diffusion of elements) process 34 . A recent review showed the formation of various complex oxides in aqueous corrosion and high-temperature oxidation of refractory HEAs 35 . Recently, a non-stoichiometric oxide solid solution was also observed by 3D-atom probe tomography (3D-APT) and X-ray photoelectron spectroscopy (XPS) during passivation of 38Ni-21Cr-20Fe-13Ru-6Mo-2W (at. %) high entropy alloy in sulfate solutions 36 . A solid solution oxide was also observed on a 38Ni-20Fe-22Cr-10Mn-10Co (at.%) MPEA passivated in NaCl at pH 4 and −0.25 V vs. saturated calomel electrode (SCE) 37 . Therefore, the solution nature of oxides as a generalization of stoichiometric ones should be incorporated in the construction of E-pH diagrams to predict the predominant reactions observed during the corrosion of such complex alloys. Calculation of E-pH diagrams with an oxide solid solution is more challenging since the nonstoichiometric oxides modeled as solid solutions would require more thermodynamic interaction parameters on certain sublattices compared with the stoichiometric oxides. In addition, the numerical process of minimizing the Gibbs energy becomes even more demanding due to the extra thermodynamic parameters in thermodynamic models for solutions.
In the current work, we investigated the predominant reactions in heterogeneous interaction systems between 10 −2 mol of FCC single-phase MPEA Ni 38 Fe 20 Cr 22 Mn 10 Co 10 (at.%, denoted as Cr22 hereafter) 37 and 1 kg of water containing 0.6 mol NaCl at 25°C and 1 atm ambient pressure. The non-stoichiometric oxides were modeled by the compound energy formulism (CEF or sublattice model) in CALPHAD approach 18 . The effects of mixing of cations on stability of oxides, the E-pH diagrams, the molality (i.e., moles of solute per kg of solvent, denoted as 'm') of predominant ions in aqueous solution, and stability of solid oxide solution phases were analyzed. The calculations were validated by the experimental results of the dissolved ions in an aqueous solution.

RESULTS
Thermodynamic modeling of oxides The free energy of stoichiometric corundum, i.e., Cr 2 O 3 , Fe 2 O 3 and Mn 2 O 3 , and oxide solution (Cr, Fe, Mn) 2 O 3 are calculated to reveal the stability of non-stoichiometric oxides vs. their stoichiometric counterparts. For stoichiometric phases, the Gibbs energy is modeled by Eq. (2) in the "Methods" section and is only dependent on the temperature as shown in Fig. 1a. The nonstoichiometric oxides can be modeled by the CEF with the mixing of alloy elements on certain sublattices. In Fig. 1b, the surface for the mixing free energy of corundum oxide solution (Cr, Fe, Mn) 2 O 3 is shown. The mixing of Cr and Fe stabilizes (Cr, Fe, Mn) 2 O 3 while mixing with Mn is not favored on the cation sublattice. As temperature increases, the mixing free energy decreases significantly due to the role of mixing entropy and mixing of Cr, Fe, and Mn are expected on the sublattice for cations. If the temperature is high enough, the role of chemical interactions becomes less important and the system will approach the ideal mixing limit, indicated by the fact that the composition corresponding to the minimum of ΔG mix approaches the equimolar limit on the cation sublattice (Fig. 1c). Therefore, the mixing of alloy elements will lead to a decrease in mixing free energy, which stabilizes the corundum oxide solution. And the stabilizing effect becomes much more significant at high temperatures.

E-pH diagram with stoichiometric oxides
The E-pH diagram of Cr22 MPEA, and various compounds (e.g., oxide, hydroxide and oxy-hydroxides) at 25°C and 1 atm ambient pressure is calculated for the heterogeneous interaction between 10 −2 mol of alloy (i.e., N(Ni) = 3.8 × 10 −3 mol, N(Fe) = 2 × 10 −3 mol, N(Cr) = 2.2 × 10 −3 mol, N(Mn) = 1 × 10 −3 mol, N(Co) = 1 × 10 −3 mol) and 1 kg H 2 O containing 0.6 m of NaCl. The standard hydrogen electrode (SHE) is chosen as the reference for the electrode potential. To assess the thermodynamic stability of the alloy, the pseudo-binary phase diagram of Ni (60−x) Fe 20 Cr x Mn 10 Co 10 (at.%) is calculated, which reveals a large single-phase region for FCC phase from 905 to 1611 K for x = 22 ( Supplementary Fig. 1). Hence, only the FCC phase is considered for the Cr22 MPEA in the calculations of E-pH diagrams. To reveal the complete stability regions of oxide phases, both gas phases and gaseous species dissolved in aqueous solution (i.e., H 2 and O 2 ) are excluded from the calculation, which otherwise would be overlaid with the reaction products and covered by gas-phase regions of predominance existing outside of the 1.23 V vs. SHE, a wide region of water stability.
The calculated E-pH diagram for Cr22 MPEA with stoichiometric oxides is shown in Fig. 2. The lower and upper dashed lines denote the pH and potential where water decomposes to 1 atm gaseous H 2 and O 2 , respectively, and are referred to as hydrogen evolution reaction (HER) and oxygen evolution reaction (OER). Immunity is defined as the scenario where the metal is stable and none of the metallic elements is oxidized 7 . Hence, neither solid corrosion product is formed nor do the metallic elements dissolve into an aqueous solution in the immunity region. In the current work, an E-pH condition is considered immunity if no solid reaction product is formed and the largest concentration of all the dissolved ions for each applicable alloy elements is <10 −6 m 6-10,12 (see region 1 in Fig. 2). Other commonly used criteria for critical ion concentration included 10 −6 gram-atom per liter 3,4 , 10 −8 m 12 , and ion activity of 10 −6 5 , etc. Various critical ion concentrations of different magnitudes were also adopted to examine its effect on phase boundaries [3][4][5][6] . Note that the choice of critical ion concentration for the definition of immunity is arbitrary and modifications might be more appropriate for particular applications 3,4 . With this criterion, the critical potential for the onset of dissolution at low pH, i.e., the upper boundary of the immunity region, is predicted to be −1.31 V vs. SHE (line between region 1 and 2 in Fig. 2). This critical potential between immunity and active dissolution is determined by the chemical potential of the most active elements in the alloy, i.e., Mn in Cr22 MPEA. Similarly, the highest potential under which the dealloyed FCC phase can be stable (between region 2 and 3 in Fig. 2) is determined by the most inert alloy element, i.e., Ni in Cr22 MPEA. The right boundary of the immunity region is determined by the formation of Cr 2 O 3 corundum at high pH (between region 1 and 14 in Fig. 2). With increased potential at low pH, the alloy elements dissolve into the aqueous solution and are stable as ions. As such, no solid species is formed even as all the metallic elements are dissolved (green region in Fig. 2), which corresponds to active corrosion.
Passivity is defined as the region in E-pH diagram where solid reaction products are thermodynamically stable. Note that thermodynamics alone cannot determine whether the solid reaction products will stay in a stable and protective form for a given time or not, thus an experimental input is required to support the thermodynamic prediction 5,29 . In Fig. 2, regions 4-29 predict various thermodynamically stable solid compounds including oxides, hydroxides, and oxy-hydroxides. Two corundum phases Fe 2 O 3 and Cr 2 O 3 are predicted to be stable in Fig. 2, with Fe 2 O 3 being stable over regions 5-10 and 15 while Cr 2 O 3 would be stable over 9-11, 19-20, 22-23, 26, and 29. Note that Fe 2 O 3 corundum is mainly stable at high potentials while Cr 2 O 3 is stable at low and medium potentials. Multiple spinels are thermodynamically stable over regions [16][17][18][19][20][21][22][23][24][25][26][27] and would coexist as a mixture. In the calculation, the stable region of NiO and Ni(OH) 2 always coincide with each other, which is due to the closeness between the free energy of reactions that produce NiO and Ni(OH) 2 at room temperature 38 . Considering that NiO and Ni(OH) 2 usually coexist as complex corrosion products in experiments [38][39][40][41] , the region for NiO and Ni(OH) 2 is labeled "NiO/Ni(OH) 2 ".
Dissolved gaseous species (i.e., H 2 and O 2 ) in the aqueous solution can also play important roles in corrosion processes. For (a) (c) example, water or proton reduction in the stable region for H 2 affects the initiation of stress corrosion cracking 42 and reduction of dissolved oxygen as a cathodic half-cell reaction occurring below OER line in the stable region for water and hydroxyl enable spontaneous corrosion 43 . To analyze the effects of dissolved gaseous species on the stability of metallic species outside the stable region for water, the E-pH diagram for Cr22 MPEA is also calculated with the H 2 and O 2 gas phases excluded while dissolved gaseous species are included ( Supplementary Fig. 2). This corresponds to the situation where the reactions producing gaseous species are slow due to higher kinetic barriers than other reactions, hence the gas phases can be excluded 44 . The diagram with gaseous species in aqueous solution is the same as the one without gaseous species (Fig. 2) within the stable region for water (between the potentials for HER and OER). Above the potential for OER, water is thermodynamically unstable and the transition to O 2 species (i.e., 2H 2 O ! O 2 þ 4H þ þ 4e À ) is favored. Since the gas phases are excluded in the analysis, the aqueous solution will be oversaturated with dissolved O 2 species above OER. As a result, several horizontal boundaries of oxides at high potentials are predicted, which agree with the features in the calculated E-pH diagrams for C-22 alloy with gaseous species in aqueous solutions 27 . Below the potential for HER (lower dashed line in Supplementary Fig. 2), the aqueous solution will be saturated with H 2 species and the calculations do not converge if the concentration of H 2 species in the aqueous solution becomes unrealistically large (red region in Supplementary Fig. 2).

E-pH diagram with non-stoichiometric oxides
To incorporate the solution nature of the non-stoichiometric oxides, the CEF with multiple sublattices is adopted, where the site fractions of elements on each sublattice are optimized to find the equilibrium composition of the corresponding phase. However, the sublattice model introduces a large number of site fractions as the degrees of freedom (DOFs) in the Gibbs energy minimization process. A large number of DOFs pose challenges to the numerical convergence and make it more difficult to find the equilibrium state (defined by the global minimum in free energy  Supplementary Fig. 3). This makes it difficult to automatically track the phase boundaries of non-stoichiometric oxides on E-pH diagram using POLY-3 module.
With the free energy minimization in POLY-3 module, the equilibrium state can be determined for a given set of conditions (namely temperature, pressure, amount of water, amount of the alloy, and concentrations of dissolved species), which corresponds to a single point on the E-pH diagram. To manually track the boundaries for predominant reactions, the equilibrium information from many combinations of E-pH conditions are collected. The unreasonable results due to numerical non-convergence are excluded, and then the phase boundaries are determined accordingly. To test the validity of this approach, the variations of predominant reaction products as a function of pH at 0.5 V vs. SHE and as a function of potential at pH = 4 are analyzed. At 0.5 V vs. SHE and pH lower than 2, no oxide is predicted to be stable as shown in Fig. 3a, b. Here, the sum of concentrations of cations in the aqueous solution reaches the corresponding amount of elements in the alloy, i.e., 3.8 × 10 −3 m, 2 × 10 −3 m, 2.2 × 10 −3 m, 1 × 10 −3 m, 1 × 10 −3 m for Ni 2+ , Fe 2+ , Cr 3+ , Mn 2+ , and Co 2+ , respectively. Under such conditions, the aqueous solution is the only stable phase and all the alloy elements are dissolved,  indicating active corrosion of the alloy. As pH increases, solid solution corundum and spinel containing various elements but not all, as well as stoichiometric Ni 6 MnO 8 and NiMnO 3 , are predicted to be stable, meanwhile dissolved Cr and Ni hydroxide ions become stable in aqueous solution. If the potential is low enough at pH = 4, only the FCC solid solution, other than the aqueous solution, is predicted to be stable. As such, the concentrations of all the ions in the aqueous phase are below a threshold value of 10 −6 m, indicating thermodynamic immunity to corrosion [6][7][8][9][10] . As potential increases, the concentration of Mn 2+ surpasses 10 −6 m at the critical potential of −1.31 V vs. SHE, and dealloying by the dissolution of Mn begins. Although Mn starts to dissolve under such conditions, the other alloy elements (namely Ni, Fe, Cr, and Co) remain unoxidized in the alloy, thus the FCC phase with a different composition from the bulk Cr22 MPEA is still stable (Supplementary Fig. 3). At further increased potentials, oxide solution corundum and spinel, in addition to Fe oxyhydroxide (FeOOH), are predicted to be stable. Note that the critical potential for the onset of active corrosion (i.e., −1.31 V vs. SHE) in Fig. 3 for non-stoichiometric oxides is the same as reported in Fig. 2 since the metal is directly oxidized to aqueous cations in a low pH region and none of the oxides are stable. Hence, the solution nature of oxides does not affect certain electrochemical reactions at relatively low potentials. In case of heterogeneous distribution of alloy elements or phase separation in the parent alloy, preferential dissolution could lead to local curvatures in morphology, which depend on the scale of the compositional inhomogeneities and may change the local potential 46,47 . Such behavior cannot be analyzed using the current approach with the thermodynamics of bulk phases.
With the results from many single-point (E-pH) equilibrium calculations, the regions of stability on the E-pH diagram for each non-stoichiometric oxide are manually tracked and their borders are converted to phase boundary lines, as shown in Fig. 4. Compared with Fig. 2, the regions for immunity, dealloying of the FCC phase, and active corrosion are almost the same as Fig. 4. However, the phase boundaries in the passivity region of Fig. 4  , are found to be stable over large portions of the E-pH diagram. The findings of corundum and spinel are broadly consistent with XPS and 3D-APT results on Cr22 where Ni, Cr, Fe, Mn, Co are detected in the oxide formed during passivation in non-stoichiometric proportions with binding energies that cannot rule out corundum and spinels 37 . The oxide solution-phase modeled by CEF is a general case and reduces to the pure stoichiometric version if the alloying elements are limited. Comparing the stable regions for stoichiometric with nonstoichiometric oxides (e.g., Cr 2 O 3 + Fe 2 O 3 in Fig. 2 vs. corundum solid solution in Fig. 4), the phase boundaries are moved although the shapes of stable regions are not changed significantly for both corundum and spinel solid solution phases. Therefore, the mixing of alloy elements on the sublattices of oxide structures is important in determining the stability and composition of oxide scales during passivation.
The concentrations of the dominant ions in the aqueous solution can also be determined from the single-point equilibrium calculations. The variations in the concentrations of several ions are shown in Fig. 5. Numerical tests demonstrated that ion concentration <10 −9 m is no longer reliable with current settings for equilibrium calculations (see "Methods" section). Hence, the distribution of ions with concentration varying from 10 −2 to 10 −9 m are shown. The distribution maps for most of the ions (e.g., Ni 2+ , ) follow similar shapes as the threshold values change, while those for Mn 2+ and Co 2+ change due to the formation of other oxides or hydroxides. Note that, distribution maps for some ions (especially for Co 2+ ) shows irregular noisy data, which are due to the numerically nonconvergent single-point calculations. For practical purposes, ions with concentration <10 −6 m will be difficult to detect and may be ignored in the analysis [6][7][8][9][10] . Therefore, the threshold value of 10 −6 m is used to show the variation of the dominant ions in the aqueous solution as a function of potential and pH, shown in Fig. 6. At two regions on the diagram, there are no ions in the aqueous solution, namely the immunity region, and low potential/high pH region (i.e., labeled 32 in Fig. 6). The former is expected as there are no alloy elements dissolved in an aqueous solution within the immunity region. For the latter, the numerical calculations show that all the metallic elements either remain in dealloyed FCC phase or form solid stable species such as oxides, hydroxides, and oxy-hydroxides. Other than these two pH domains (i.e., regions 1 and 32 in Fig. 6), the predominant ions in aqueous solution oxidized from alloy are shown. At regions with low potential or pH, cations for metallic elements with various valence states are dominant. Other ions, namely chromate (VI) anions, manganate (VII) anion, and hydroxide ions with positive (e.g., Cr(OH) 2+ ), negative (e.g., Ni (OH) 3 − ), or neutral complexes (e.g., Ni(OH) 2 ), becomes prevalent with increasing potential and pH.

Cation distribution in oxides
The non-stoichiometric spinel, corundum, and cubic M 2 O 3 modeled by CEF are predicted to be stable under certain conditions for Cr22 MPEA, shown in Figs. 7, 8 and Supplementary  Fig. 4, respectively. Despite some shifts in phase boundaries, the regions for non-stoichiometric oxides overlap with combined regions of their stoichiometric counterparts. The molar fractions   (denoted as x) of alloy elements in spinel and corundum are shown in Figs. 7 and 8 Fig. 4). Note that, different from corundum and spinel oxide solutions, mixing of O with cations is predicted on the sublattice for O in cubic M 2 O 3 oxide solution at medium potentials so that the sum of the molar fractions of Ni and Mn is larger than 40%.
The current calculations employing the CALPHAD-type databases only included the thermodynamics of bulk phases. Therefore, it cannot predict the phase stability as altered or modified by interfaces/surfaces, defects, strain energy, or kinetic processes 18  Co Fig. 7 The distributions of Ni, Fe, Cr, Mn, and Co in the spinel phase at various potential and pH conditions are shown, respectively. Molar fractions of alloy elements in the spinel phase on E-pH diagram.  thermodynamics of equilibrium bulk phases are necessary to account for the effects due to surface energy (or size effects), surface termination, water adsorption, substrate, etc 48 . In addition, the metastable phase defined by local free energy minimums will not appear in the calculated diagram since the global minimum of Gibbs energy is searched for each single-point equilibrium calculations. In the current analysis for Cr22 MPEA, stoichiometric Ni 2+ based halite is predicted to be stable while solid solution halite (Ni,Cr)O was observed in the experiments on Ni-22Cr-6Mo (wt.%) alloy by Yu et al. 34 . To further understand the reason for the stoichiometry of halite, the solubility of elements in halite phase in equilibrium with other phases in Ni-Fe-Cr-Mn-Co-O system is plotted as a function of temperature ( Supplementary Fig. 5). At 300 K, the equilibrium molar fraction of Mn and Fe are 0.13 % and 0.015 %, respectively, while those of Co and Cr are zero. Therefore, non-stoichiometric halite is not stable in E-pH diagram of Cr22 MPEA at 25°C due to the fact that other cations do not tend to mix with Ni (II) on the cation sublattice of NiO halite.

Experimental validation
To evaluate the current thermodynamic analysis with an experiment, the atomic emission spectroelectrochemistry (AESEC) [49][50][51][52] was introduced to measure the total quantity of dissolved alloy components in aqueous solution, which was shown to be consistent with the 3D-APT and XPS analysis 36,37,52,53 . AESEC measurements of dissolved alloy components were performed in 0.1 M (i.e., mol per liter or molarity) NaCl, pH = 4 solution at applied potentials from 0.00 to 0.34 V vs. SHE. The dissolved elements from a series of potentiostatic experiments were quantified as integration of the dissolution rate (mol s −1 ) vs. time (s) curve, then divided by the total volume of electrolyte (cm 3 ), assuming the density of electrolyte as 1.00 g cm −3 . Note that the experimentally measured amount of each element is calibrated based on the assumption that Ni dissolved completely from the MPEA to the solution, which agrees with the current thermodynamic calculations at pH = 4 (Figs. 3d and 5).
Considering that the effect of critical ion concentration on the boundaries of predominant species is usually examined in logarithmic scale [3][4][5][6] , the comparison between AESEC measurements and thermodynamic calculations is shown with logarithmic scale in Fig. 9. To quantify the extent of agreements in between, the relative error is defined as |log 10 C AESEC − log 10 C calc |/log 10 C AESEC , with C as the concentration of ions (in molality). The calculated concentrations for Ni 2+ agree well with the AESEC measurements (<0.1% relative error), supporting the assumption of preferential complete Ni 2+ dissolution. The disagreements for Fe 2+ at 0.00 and 0.14 V vs. SHE are <3.0%, while the deviations are obvious at 0.24 and 0.34 V vs. SHE (i.e., larger than 71 %). Note that large errors are found around the transition between competing Fe species (i.e., Fe 2+ and (Fe,Cr) 2 O 3 ), and thus is likely to be caused by the insufficient accuracy of the thermodynamic data. Mn 2+ showed a good agreement in that the largest relative error is~15 %. Compared with Ni 2+ , Fe 2+ , and Mn 2+ , large deviations are observed for Cr 3+ and Co 2+ , with relative error larger than 100%. In the thermodynamic analysis, the low concentrations of Cr 3+ and Co 2+ in the aqueous solution are due to the formation of solid oxide solutions, i.e., spinel (Co,Fe)Cr 2 O 4 with~1.3 at.% Fe and 13.0 at.% Co and corundum (Fe,Cr) 2 O 3 with~5 at.% Cr at pH = 4. Recently, at pH = 4 and −0.25 V vs. SCE (−0.006 V vs. SHE), solid solution products (i.e., consistent with spinels and corundum) were reported in the passive film on Cr22 MPEA 37 . Factors not captured in the E-pH prediction apparently prevail in experiments and mostly hydroxides such as Cr(OH) 3 and possibly Cr containing corundum and spinel prevail as indicated by XPS 37 . Current thermodynamic analysis indicates that Co and Cr are predicted to form mainly CoCr 2 O 4 at pH = 4 from 0.00 to 0.34 V vs. SHE (Fig. 7). This is not evident in experiments as Co 2+ dissolution is kinetically fast in Cl − and does not appear in oxides formed. In AESEC experiments, Cr showed a non-congruent dissolution (i.e., lower than the proportional alloy composition) over most of the potential domain investigated in this work, indicating that Cr is present in the oxide film 37 . Nevertheless, dissolved Cr and Co cations in the solution were~1 or 2 orders of magnitude higher than thermodynamic prediction as shown in Fig. 9.

DISCUSSIONS
In some early E-pH diagrams for alloys, the dissolution potentials for the alloying elements were calculated separately assuming unit activities (i.e., standard states) for alloy elements, and upper boundary of the immunity region was determined by the dissolution potential of the most active elements at low pH 7-10,12 . However, this approach cannot predict the change of equilibrium composition of the dealloyed phase with increasing potentials since the thermodynamics of the actual alloy phase is not considered. In the current work, the thermodynamics of the alloy phase is treated as a solid solution modeled by CALPHAD 18  . By adopting thermodynamics of the alloy phase, the remaining metallic elements in the alloy and the dissolved species sum up to moles of the alloy with nominal composition. Therefore, the equilibrium composition of dealloyed FCC phase and the dissolved species can be predicted consistently in this way.
Due to the mixing of metallic elements on certain sublattices, the oxide solution phases can be more stable than stoichiometric ones and the composition of the oxides become nonstoichiometric. Moreover, the boundaries in E-pH diagram with non-stoichiometric oxides become less densely spaced (Fig. 4) since boundaries for stoichiometric ones are merged due to the mixing of cations. At low temperatures (25°C), the role of mixing free energy (i.e., −TS with S being the mixing entropy) is limited compared with enthalpy, which is determined by the crystalline structure of oxides. Hence, the E-pH ranges of the stable regions of oxide solution phases are determined by their crystalline structures and do not change significantly due to mixing. For an alloy with different composition, the oxide solutions will form at the same (E, pH) conditions as that for Cr22 MPEA as long as the necessary metallic elements are present in the system, although the equilibrium composition of oxide solid solutions may change inevitably due to the amount of available alloy elements. As temperature increases, the energetic contribution from mixing free energy increases (Fig. 1c) and the oxide solution phases become more stable at high temperatures. For MPEAs, mixing entropy becomes more significant due to an increased number of principal elements. Therefore, the solid solution oxides are expected to be more prevalent in the passivation layer of MPEAs at high temperatures.
Current tests show clear differences between the predicted amount of ions and that from AESEC experiments for Co 2+ and Cr 3+ . This could be due to the surface complexation resulted from the formation/dissolution of Cr-based oxide 54 . Both the depletion and enrichment of Cr from the oxide have both been observed for an equimolar CrFeCoNiMn MPEA tested in 0.1 M NaCl 55 and 0.1 M H 2 SO 4 56 , respectively. It is also reported that Co-based oxide dissolves more rapidly than other oxides in acidic media 57 . Based on the thermodynamic analysis and AESEC measurements, we postulate that CoCr 2 O 4 spinel forms first, and then dissolves rapidly in acidic media following the reaction 58 : Further, Cr 2 O 3 dissolves afterward adding to surface complexation, leading to the phenomenon that Co dissolves congruently while Cr forms oxides, but a large amount of Cr 3+ was also detected in the aqueous solution eventually. Note that the interaction parameters in thermodynamic models are evaluated from the thermal properties of equilibrium bulk phases 18 . Hence, the effects of surface adsorption, local curvatures in morphology, defects, interfaces, microstructures, and strain energy cannot be incorporated which will inevitably contribute to the deviations between thermodynamic calculations and experiments along with kinetic factors.

E-pH diagram calculations
The thermodynamic descriptions of stoichiometric and non-stoichiometric oxides are determined by the theoretical models for the free energy, e.g., Gibbs energy for isobaric processes, which is the case for aqueous corrosion under atmospheric conditions. For a stoichiometric phase θ (i.e., fixed composition), the molar Gibbs energy (G θ m ) within a certain temperature range is only a function of temperature 18 , where b i is the molar fraction of i-th element and P i b i H SER i represents the sum of the enthalpies of the elements in their standard element reference (SER), i.e., the stable state at 298.15 K and 1 bar. The reference term is needed because there is no absolute value of the enthalpy of a system. The coefficients in Eq. (2) represent the temperature dependence of free energy and are usually obtained by fitting Eq. (2) to the thermal properties of θ phase from experiments or first-principles calculations. To deal with the solution nature of non-stoichiometric phase θ, the molar Gibbs energy is represented by CEF (or sublattice model) 18,45 , where I 0 is a constituent array of zeroth-order (i.e., one constituent on each sublattice) specifying one constituent in each sublattice and P I0 Y ð Þ is the corresponding product of the constituent fractions specified by I 0 . G 0 I0 is the compound energy parameter representing the Gibbs energy of formation of the compound I 0 with SER as reference states. The factor a s is the number of sites on sublattice s and y s i denotes the site-fraction of component i on sublattice s. P Ik Y ð Þ is the k-th order product of the constituent fractions specified by I k and L Ik is the corresponding temperature-dependent interaction energies. The three terms on the right-hand-side of Eq. (3) represent the contributions from the unreacted mixture, ideal mixing entropy, and non-ideal interactions, respectively.
The Thermo-Calc software is adopted to analyze the electrochemical reactions and construct the E-pH diagram for heterogeneous interaction system between 10 −2 mol of FCC single-phase Cr22 MPEA and 1 kg of water containing 0.6 m NaCl at 25°C and 1 atm ambient pressure. The thermodynamics of Cr22 MPEA, aqueous solution (and activity of ions), stoichiometric compounds, and oxide solid solutions are described by the TCHEA3, TCAQ3, SSUB6, and TCOX8 databases 59 , respectively. The gas phase for E-pH predominance outside the region for stable water is suppressed in all calculations to reveal the stability regions for the predominant reaction products at low and high potentials. The gaseous species dissolved in aqueous solution is excluded with the exception for Supplementary Fig. 2 to study the case where the aqueous solution is supersaturated with H 2 and O 2 molecules. The Gibbs energy minimization is performed using POLY-3 module with the maximum number of iterations as 20,000, required accuracy for convergence as 10 −6 , smallest possible fraction as 10 −20 . Additional numerical tests are performed to check whether the equilibrium results are the global minimums, not local ones, in Gibbs energy surface. Despite some small differences, the diagrams with and without checking for global minimums agree with each other, hence the global minimum check is disabled to save computing time.
In the analysis with oxide solution phases modeled by Eq. (3), a large number of additional DOFs are needed to account for the site fractions of cations on sublattices, which leads to challenges to achieve numerical convergence at certain conditions. To avoid this issue and manually map the E-pH diagram, the equilibrium states under many set of conditions on the grids of pH and potential axis are computed with the intervals of pH and potential as 0.125 and 0.025 V, which corresponds to a total of 13673 data points on the E-pH diagram. The equilibrium results under these conditions are analyzed and the non-convergent data points are excluded in the analysis. The remaining data are used to map the E-pH diagram, concentrations of ions in aqueous solution, and molar fractions of alloy elements in oxide solution phases.

AESEC experiment
The Cr22 MPEA was produced and characterized indicating a single-phase FCC structure 37 (also see the phase diagram in Supplementary Fig. 1). The sample was degreased with ethanol in an ultrasonic bath for 10 min, rinsed with deionized (DI) water (18.2 MΩ cm) then dried with flowing Ar. The sample surface was ground with the Si-C paper (P4000) under DI water then again dried with flowing Ar. 0.1 M NaCl electrolyte was prepared with DI water. The final pH was adjusted to 4.0 by adding 1.0 M HCl solution. Prior to the experiment, the electrolyte was deaerated by flowing Ar for 30 min and maintained during experiments.
The total quantity of dissolved alloy components is obtained by AESEC technique [49][50][51][52] . This technique measures simultaneously all the dissolved species in the electrolyte by coupling an Ultima 2C Horiba Jobin-Yvon inductively coupled plasma atomic emission spectrometer (ICP-AES) with a specially designed flow cell. The principle and quantification method are described in detail elsewhere 49 . The specimen of interest is in contact with the electrolyte in the electrochemical flow cell where a reference electrode and a counter electrode (Pt foil) are located in a separated compartment with a porous membrane. An SCE (−241 mV vs. SHE) was used as a reference electrode and then converted to SHE for comparison with the predicted thermodynamic data. A Gamry Reference 600 TM potentiostat was utilized to control the potential during the experiment. A constant potential was applied when the flowing electrolyte was brought to the sample. Prior to each potentiostatic experiment, a cathodic potential at −1.06 V vs. SHE was applied for 600 s to reduce the pre-existing oxide. The total amount of each alloy component in the electrolyte was calculated by integrating the dissolution rate vs. time curve.

DATA AVAILABILITY
The authors declare that the data supporting the findings of this study are available within the paper and its Supplementary Information files. The simulation output files are available upon reasonable request. They are not publicly available due to the very large file sizes. Parameters of the input files are described in the computational methods.

CODE AVAILABILITY
The Thermo-Calc Software is proprietary software available for purchase at https:// www.thermocalc.com/. The input files for Thermo-Calc console mode and the Python scripts used to process output files of Thermo-Calc are available in Data set 1.