Calculations of oxide formation on low-index Cu surfaces

Chemical stability is an important property for materials exposed to air or water. It is important, therefore, to understand on a fundamental level the surface chemistry by which materials interact with their environment. Here, we focus on the low-index surfaces of copper, Cu(100), (110), and (111), which have a long history as model systems in metal oxidation studies.¹ 

Previous experimental and theoretical studies have shown that the kinetics and morphology of Cu oxidation are complex and influenced by many factors, including surface orientation, ambient temperature, surface defects, and electric fields.^2–9 Most previous studies have focused on the formation of a monolayer (ML) of surface oxides.^10–21 However, the process of oxygen transport into the Cu sub-surface is a key step for the transformation between the oxygenated Cu surfaces and bulk oxides. Based on x-ray photoelectron spectroscopy (XPS), x-ray induced Auger electron spectroscopy (XAES), and scanning tunneling microscopy (STM) characterization of Cu(100) oxidation, it was understood that subsurface oxidation leads to the growth of disordered Cu₂O.^22,23 Based upon the observed oxidation of the Cu(100) surface,^24–28 it is believed that the nucleation of Cu₂O on the Cu(110) surface is also initiated by subsurface oxygen.^29,30 Although there are fewer studies of the Cu(111) surface, STM measurements suggest that oxygen etches into the step edge to form over-layer oxide structures at high oxygen coverage.^31–34 The oxide-like over-layer can be viewed as the initial layer of a Cu₂O(111) film, which can potentially act as a template for the further growth of Cu₂O(111) layers.

Even though the oxidation of Cu low index surfaces has received considerable attention, the detailed atomistic mechanisms that govern the formation of copper oxide are not well-understood. In any viable mechanism, O atom must be supplied from O₂ dissociation, which is expected to become a bottleneck for further oxidation when the surface is saturated with oxygen. On some metals, including Al,^35–37 formation of a compact oxide layer can effectively hinder further oxidation of the bulk material. Thus penetration of oxygen through the surface oxide is critical for further oxidation. There are few studies on how molecular oxygen dissociates on O-covered Cu surfaces, or how O atoms subsequently embed into the Cu sub-surface. Studies of Cu(100) oxidation, for example, have fixed the oxygen coverage at 0.5 monolayer (ML) assuming that oxygenated surfaces are inert towards further oxidation.^10,19 On Cu(110), the oxidation process beyond c(6 × 2) reconstruction has never been reported. This lack of atomic detail is not surprising because it is difficult for experimental techniques to resolve the dynamic processes responsible for surface oxidation. It is, however, possible to understand such oxidation reaction mechanisms using a computational approach.

Here, we use density functional theory (DFT) to model oxygen dissociation and diffusion and evaluate the kinetics of subsurface oxide formation on the low-index Cu surfaces.

DFT calculations were performed using the generalized gradient approximation functional of the Perdew–Burke–Ernzerhof (PBE) form³⁸ as implemented in the Vienna ab initio simulation package.^39–42 The projector augmented wave approach^43,44 was used to model the core electrons and a plane wave basis with a cutoff energy of 350 eV was used to describe the valence electrons. Brillouin-zone integration was performed using (10 × 10 × 1), (6 × 10 × 1), and (8 × 8 × 1) Monkhorst-Pack grids⁴⁵ for the (1 × 1) surface unit cells of Cu(100), Cu(110), and Cu(111) surfaces, respectively. The Cu surfaces were modeled using periodic slabs with six layers parallel to the surfaces, where the bottom two layers were fixed in their bulk positions. For larger surface cells, correspondingly smaller grids were used to ensure an equivalent sampling in reciprocal space. Periodic images along the direction perpendicular to the surface were separated by a vacuum region of 15 Å. All of our calculations were closed-shell except for those involving free molecular or atomic oxygen, where the calculations were spin polarized. Adsorption was done on one side of the slab only. The positions of all the atoms, except those of the bottom layers, were relaxed until the force on each atom was less than 0.01 eV/Å. Reaction barriers were calculated with the climbing-image nudged elastic band (NEB) method,⁴⁶ where at least five intermediate images were included between the initial and final states. Low-energy structures of oxygen adsorption on the Cu low-index surfaces were obtained using basin-hopping global optimization.⁴⁷ The STM images are simulated by plotting the partial charge density in the energy range from −0.1 to 0.1 eV with respect to the Fermi level. The images are taken from a slice 0.2 Å above the top-most atomic layer.

To understand the relative stability of different oxidized surfaces, we calculated the formation energy per atom as a function of O concentration, to construct a convex hull. Any structure with a formation energy higher than a linear combination of two neighboring compositions is unstable with respect to phase separation into those two structures. Structures unable to decompose in this way form the convex hull. In our definition of oxygen concentration, we consider only the percentage of O atoms in the top three layers; the remaining Cu slab is regarded as a boundary condition for the surface of interest. Specifically, the formation energy of Cu_1−xO_x is calculated as

where $ECua+Cu1−xOx$ is the total energy of the adsorbate-substrate system, E_Cua is the energy of bottom three layers, E_Cu is the average energy of one Cu atom in a clean Cu slab, E_O2 is the energy of an isolated oxygen molecule, x is the percentage of O atoms in the top three layers, $EsubstrateCu1−xOx$ is the energy of Cu_1−xO_x on the three-layer Cu substrate.

The slope of any point with oxygen concentration x connecting to the end point, x = 0, is calculated as

The oxygen binding energy or chemical potential with respect to clean Cu is calculated as

Comparing Equations (2) and (3),

the slope is equivalent to the oxygen binding energy.

In the previous work, an average oxygen binding energy was used to determine phase stability. This stability criterion, however, assumes that appropriate reference states are gas phase O₂ and a clean Cu slab. This is appropriate at low oxygen coverage but it can predict incorrect stable structures at high O concentration. The convex hull construction, that we use in this present work, is appropriate for determining stable Cu–O complexes even at high O concentration; it is widely used for studying the entire phase space of alloys.^48,49

The oxygen adsorption energy E_ads is calculated as

where $ECuO2slab$ is the total energy of the Cu–O system and E_O2 is the energy of an isolated oxygen molecule. E_ref is the energy of the substrate without O₂ adsorption. Adsorption energies are defined such that a negative value indicates a thermodynamically favorable adsorption process.

The first step in copper oxidation involves O₂ dissociative adsorption. Our calculations show that dissociation of O₂ on the clean Cu(100) surface is spontaneous, which is consistent with the previous work.^11,50,51 There are three different O atom adsorption sites on the surface as shown in Fig. 1(a): the hollow (H), bridge (B), and top (T) sites. The highly coordinated hollow binding site is preferred at all oxygen coverages.^19,52,53 The structure of the unreconstructed Cu(100)-c(2 × 2) with 0.5 ML oxygen⁵⁴ and reconstructed Cu(100)-missing row (MRR) surfaces^55–57 is shown in Figs. 1(b) and 1(c), respectively. In the c(2 × 2)−0.5 ML structure, O atoms occupy the hollow sites in every other row, and the MRR reconstruction is characterized by a missing row of copper atoms in every fourth row. For clarity, the O coverage is defined as the ratio of the total number of O atoms to the number of Cu atoms per layer in the corresponding clean surface.

The basin-hopping algorithm was coupled with DFT to obtain globally stable oxide structures with large supercells with various oxygen concentrations. Structures with formation energies (Equation (1)) on the convex hull (see Fig. 2) are stable with respect to decomposition into neighboring phases. The Cu(100)−0.25 ML, Cu(100)-MRR, and Cu₂O (Cu₂O(111) surface is the most stable among low-index Cu₂O surfaces)⁵⁸ are located on the convex hull, in agreement with the previous reports.^19,24,59 The Cu(100)-c(2 × 2)−0.5 ML structure, however, which was not shown in the phase diagrams of the previous work,^19,24 is found on our hull. According to the same slope of the 0.25 ML and c(2 × 2)−0.5 ML, the c(2 × 2)−0.5 ML can be observed while 0.25 ML cannot under a certain oxygen partial pressure because of the lack of a controlled O/Cu ratio. Moreover, the existence of c(2 × 2) is well documented experimentally,^13,54 but there is no evidence for a 0.25 ML structure. Structures with high oxygen coverage above the hull are thermodynamically unstable and are unlikely to be seen experimentally. During the initial stage of oxidation, STM experiments show that the copper surface is covered by either the c(2 × 2)−0.5 ML and/or the MRR structure. Further oxidation of copper inevitably starts from either of these two surface phases, but the mechanism is not known. Thus a theoretical investigation of oxygen dissociation and diffusion on the c(2 × 2)−0.5 ML and MRR surfaces is desirable.

Oxidation of the MRR surface is a possible mechanism for Cu₂O formation. By comparing various adsorption sites and dissociation paths, we find a feasible mechanism for O₂ dissociation and O diffusion into the MRR subsurface. Fig. 3 shows (a) the energetics of our reaction path and (b) the structures of key stationary points along the path.

The O₂ adsorption energy on the MRR surface is +0.18 eV. Dissociation of the O₂ molecule goes through a saddle point (TS1) with a barrier of 0.57 eV. Since the adsorption of O₂ is endoergic, the physically relevant dissociation barrier of 0.75 eV is calculated with respect to O₂ in the gas phase. Next, the two dissociated O atoms pull one Cu to the middle of them (IM3) with a barrier of 0.15 eV, forming a stable linear O–Cu–O trimer. Finally, after some facile rearrangement of Cu and O atoms on the surface (IM4), the second adsorbed O atom diffuses to the sublayer (IM5) with a barrier of 0.35 eV. The entire process is exoergic by 0.51 eV. Thus our NEB calculation shows a reaction pathway with an overall barrier of 0.75 eV, much lower than found in the previous work,¹⁰ indicating that further oxidation of the MRR surface will be facile at 400 K (according to harmonic transition state theory, a 400 K reaction with an energy barrier of 0.75 eV occurs at 3.68 × 10³/s). O adatoms diffusing into the sublayer in IM5 indicate the onset of the bulk oxide formation.

The unreconstructed Cu(100)-c(2 × 2)−0.5 ML surface is another stable phase up to oxygen coverage of 0.5 ML. According to both experiment and theory, the MRR phase is slightly more stable than the Cu(100)-c(2 × 2)−0.5 ML phase in O rich conditions.^19,26,59 The processes of c(2 × 2) transits to MRR and further oxidation on c(2 × 2) can proceed simultaneously, so it is necessary to investigate other oxidation mechanisms of c(2 × 2) which do not involve Cu ejection and diffusion. Our calculated minimum energy path (MEP) for O₂ dissociation and O penetration on the c(2 × 2)−0.5 ML phase is shown in Fig. 4(a).

O₂ molecular adsorption (IM1) is again endoergic. To react, the O₂ molecule pushes an O atom from a surface Cu–O chain into the subsurface (IM2) overcoming an energy barrier of 0.28 eV (TS1). With O₂ in a Cu hollow site, the molecule readily dissociates into two O atoms IM3, 0.30 eV barrier. Finally, one of the dissociated O atoms diffuses into the sublayer to form a subsurface oxide unit that is stabilized by 1.43 eV (IM4) with respect to O₂ gas. The highest barrier along this path is 0.97 eV. It can be seen in Fig. 4(b) that the structure of IM4 has recovered the c(2 × 2) structure on the surface layer but that it is pushed upwards to accommodate a c(2 × 2)−0.5 ML structure in the sublayer. Following this pattern, O₂ can continue to dissociate and diffuse into the sublayer to form an oxide island.

While both the Cu(100)-MRR and c(2 × 2)−0.5 ML surfaces have significant barriers for O₂ dissociation, we have found mechanisms on both surfaces for which the subsequent oxidation of the subsurface layers is facile.

Both our calculations and previous studies show that O₂ spontaneously dissociatively adsorbs on the clean Cu(110) surface.⁶⁰ Oxygen atoms can adsorb at five distinct sites on the non-reconstructed Cu(110) surface: the hollow (H), long bridge (LB), short bridge (SB), top (T), and the shifted-hollow site (shH). The shH site is pseudo threefold coordinated and located roughly halfway between the H and SB sites. At low oxygen coverage, O preferentially binds to shH site, as shown in Fig. 5(a).¹⁹ With increasing oxygen coverage, interactions between adsorbed O atoms lead to a wider variety of adsorption sites. The complexity of such surfaces made the basin-hopping method necessary to find the ground state at each oxygen concentration and the convex hull shown in Fig. 6. Both the (2 × 1) and c(6 × 2) structures, illustrated in Figs. 5(b) and 5(c), are on the hull, in agreement with what is observed experimentally.^61–64 The (2 × 1) structure has an added row (Cu–O–Cu chains along the [001] direction) structure with every other [001] Cu atomic row absent; each unit cell consists of one Cu–O–Cu chain. For the c(6 × 2) structure, the unit cell consists of two pairs of Cu–O–Cu chains along the [001] direction connected by Cu atoms coordinating every other O atom along the chains.

Previous studies suggest that a critical oxygen coverage is required for a crossover from oxygen chemisorption to bulk oxide formation on Cu(110).²⁹ In Secs. III B 1 and III B 2, we attempt to provide some atomistic insights into the further oxidation of the c(6 × 2) and 1 ML structures. The former is the surface oxide with the highest O concentration structure on the hull and the latter one is slightly above the hull but with higher oxygen coverage.

The path from c(6 × 2) to bulk Cu₂O is likely complicated, but one process that cannot be avoided is O₂ dissociation on c(6 × 2). Fig. 7(a) shows the lowest energy reaction pathway that we calculated for this step.

O₂ spontaneously adsorbs on the c(6 × 2) surface (−0.24 eV, IM1 in Fig. 7(b)). After small surface rearrangements (IM2), O₂ dissociates (0.59 eV barrier, TS1) into two stable adsorbed O atoms (IM3). Next, we propose a transition path (0.70 eV barrier) that involves the concerted movement of three O atoms (A, B, and C, in Fig. 7(b)) laterally into a more stable site. In a final step, O atom “D” diffuses (0.72 eV barrier) to form a geometry (IM5) that is consistent with the stable structure reported in the previous DFT calculations.²⁹ The IM5 structure is 1.50 eV more stable than the initial O₂ adsorption state. Remarkably, we find that the IM5 can spontaneously adsorb O₂ with a binding energy of −0.20 eV. The O₂ dissociation barrier rises, however, to 0.92 eV. Thus the further oxidation of c(6 × 2) needs to overcome barriers for oxygen dissociation and diffusion that are prohibitive at room temperature.

Experimental STM images clearly show the sequential formation of the (2 × 1) and c(6 × 2) reconstructions under controlled oxygen pressure and temperature conditions.⁶⁵ These reports, however, focus on the morphology of the reconstructions rather than the unreconstructed portions of the surface, which were either ignored or considered to be the clean Cu substrate. Following the same logic as the Cu(100) surface, studying the oxidation of unreconstructed surface is necessary because the O₂ dissociation on the clean Cu(110) surface has no barrier while reconstruction growth from the step is not fast enough to cover the whole surface. We propose that the unreconstructed portion of the surface (i.e., not (2 × 1) or c(6 × 2)) is Cu(110) fully covered by O atoms, as shown in Fig. 8(a). This structure, herein described as 1 ML, is produced by O₂ dissociatively chemisorption on the clean Cu(110) surface that has not already reconstructed. Even though the 1 ML surface is thermodynamically unstable (above the hull in Fig. 6), kinetic limitations should make it difficult to decompose into separate c(6 × 2) and Cu₂O phases. To help make connections with the experiment, we have simulated an STM image of the 1 ML surface, along with the (2 × 1) and c(6 × 2) surfaces, as shown in Fig. 8(b). Compared with the ordered (2 × 1) and c(6 × 2) reconstructions, the disordered nature of the 1 ML surface makes it difficult to recognize STM images which might be why there is no experimental report of this structure. The work function is another way to identify the 1 ML structure. The work functions of the (2 × 1) and c(6 × 2) reconstructions are similar, as shown in Table I, and 0.3 eV higher than that of the clean Cu(110) surface. The work function of the 1 ML structure is even higher, 0.5 eV above Cu(110), making it possible to distinguish from the bare Cu surface.

In the 1 ML structure, surface Cu atoms are saturated by adsorbed O atoms, and further oxidation requires O diffusion into the sublayer. Fig. 9(a) shows our calculated MEP of the oxidation on 1 ML. Overall, the energy decreases by ∼1.1 eV upon oxidation and the highest barrier along the path is 0.31 eV. Initially, the O₂ molecule adsorbs on the surface to form IM1 with an adsorption energy of −0.20 eV. Then the O₂ molecule pushes one of the O atoms on the surface into the sublayer region (IM2) with a negligible reaction barrier. Next, O₂ dissociation occurs via TS1 with an energy barrier of 0.31 eV. In the final step, the sublayer O atom (IM3) diffuses into the third layer (IM4). The barrier for this final step is only 0.15 eV and results in a significant drop of the total energy by 0.70 eV.

In contrast to the c(6 × 2) surface, oxidation of the 1 ML subsurface is facile even at room temperature.

On the clean Cu(111) surface, the O₂ adsorption energy is −1.16 eV and the dissociation barrier is lower than 0.15 eV. Threefold hollows are found to be the sites where O atoms adsorb. Consistent with our DFT results, molecular O₂ is not detectable on Cu(111) at temperatures above 170 K.⁶⁶ O atoms can occupy the face centered cubic (FCC) hollow, hexagonally close packed (HCP) hollow, bridge (B) sites. Our calculations show that O binds most strongly to the fcc site with an adsorption energy of −1.61 eV in the low-coverage limit, which is lower than on Cu(110) or Cu(100). Many experiments have reported the formation of Cu₂O(111) on the Cu(111) surface upon oxidation.^33,34,67 Our calculated structures of the Cu(111) and Cu₂O(111) surfaces are shown in Fig. 10. The Cu₂O(111) surface consists of a three-layer repeat unit with each copper layer sandwiched between two layers of oxygen atoms. The surface unit cell is 2.35 times larger than Cu(111), giving rise to a lattice mismatch of 17.5% between the two.

Fig. 11 shows stable surface oxide structures on which we can investigate further oxidation mechanisms. Limited by the supercell size, we are unable to simulate the “44” or “29” structures that have been observed experimentally.^68,69 Common to these stable structures, however, is a Cu₂O(111)-like oxide layer that forms on the Cu(111) surface. Clean Cu(111) surfaces with low oxygen coverage (0.0625 ML, 0.125 ML, and 0.25 ML) are also on the hull. As the oxygen coverage increases to 0.5 ML, the lattice mismatch between the Cu₂O(111)-like over layer and the Cu(111) substrate destabilizes such structures. At this point, structures with less Cu in the top layer, labeled as 9Cu + 6O, 12Cu + 6O, and step + 8O in Fig. 11, become more stable, as there is more room to rearrange the top layer to better fit the substrate (Fig. 12). The step-like structures, step + 6O, step + 7O, and step + 8O, which are expected to form near step edges, lie on the hull because the top layer oxide only covers part of the surface and thus the mismatch with the substrate is minimized.

Following the same logic used to analyze the Cu(100) and Cu(110) surfaces, we start with the point on the hull with less oxygen than Cu₂O to study oxidation beyond the first layer, as well as structures above the hull with intermediate compositions. The latter are not thermodynamically stable, but the mechanism of their decomposition will correspond to Cu₂O formation.

The oxidation process of Cu(111) surface is very different from those of Cu(100) and Cu(110). First no single ordered surface oxide is seen in experiments at low temperature. Our calculated values of O₂ adsorption and dissociation on the oxidized structure mentioned above are listed in Table II. For the 0.5 ML surface, the O₂ adsorption energy is only −0.12 eV, but the dissociation barrier of 0.28 eV is low as well. The 9Cu + 6O surface adsorbs O₂ more strongly and the dissociation barrier is 0.23 eV. The dissociation barrier on the step + 8O structure is higher (0.78 eV). The 12Cu + 6O reconstruction, which has the structure most similar to one layer of Cu₂O(111), has the strongest O₂ adsorption among the oxidized surfaces, but the O₂ dissociation energy is also high (0.65 eV). Since the number of Cu atoms in the Cu(111) top layer cannot exactly match any of the above structures, a mixture of surface structures is expected. In our calculations, once the O₂ molecule dissociates, the O atoms spontaneously move to the sublayer. Thus the barriers for O diffusion on the oxidized surfaces are not listed, since O₂ dissociation is the rate-limiting step.

In this work, we have calculated oxidation mechanisms of low-index Cu surfaces at high oxygen coverages. At low oxygen coverage, the adsorption of oxygen on the surface is energetically favored and dissociation barriers are negligible. As oxygen atoms accumulate on the surfaces, O₂ dissociation becomes more difficult, especially when the surface Cu atoms are saturated by O atoms. The oxidation process beyond this surface is the key to bulk oxide formation. The first question addressed is what structure does this saturation point correspond to? The point on the convex hull with the highest oxygen concentration before Cu₂O is certainly relevant. However, this single point is always a reconstruction that requires a specific Cu/O surface stoichiometry, and thus this phase cannot cover the entire Cu surface, as the concentration of Cu adatoms will vary. Further oxidation may also occur on the residual un-reconstructed oxide surfaces, which forms spontaneously from the clean Cu surface in the presence of oxygen. The barriers for oxygen dissociation and diffusion on the three low index surfaces are summarized in Table III. The dissociation and diffusion processes must occur sequentially, whereas oxidation on the reconstructed and un-reconstructed surfaces can occur in parallel when both are present. The lowest energy dissociation process on each surface, which is also the rate-limiting step for subsurface oxidation, is highlighted in bold text. The order of oxidation rates of the low-index Cu surfaces at high oxygen coverage is (111) > (110) > (100), with (111) being the easiest to oxidize. This order is consistent with the observations by Zhou et al. Using TEM observations, they found that the nucleation rate of oxide islands at 350 °C on Cu(111) is much faster than that on the Cu(100) and Cu(110) surfaces.^7,70

To understand these results, we return to the original structures of the Cu(100), Cu(110), and Cu(111) surfaces. The Cu(100) surface has the most open structure and rapidly dissociates O₂ at hollow sites, as observed in experiments. With increasing oxygen coverage up to 0.5 ML, the surface Cu atoms become saturated by O, leading to a stable surface oxide which hinders further O₂ dissociative adsorption. The Cu(111) surface is the most close-packed and stable surface. It is then expected that the reactivity of molecular O₂ is lower on Cu(111) than on the more open Cu(100) and Cu(110), at low oxygen coverage. At high oxygen coverage, the O-induced reconstructions of the Cu(111) surface result in corrugated surface overlayers. These disordered oxidized surface structures provide a facile diffusion channel for O atoms into the sublayer.² In brief, the more stable the initial structure, the less stable the surface oxide, and the easier it is for oxygen to penetrate into the bulk.

Although there is no consensus on a single rate law describing the oxidation dynamics for copper oxidation,⁷¹ our results provide a qualitative explanation for the observed kinetics governing the initial stages of oxide island nucleation.

We have investigated the oxidation on Cu(100), Cu(110), and Cu(111) surfaces using DFT at high oxygen coverages, including O₂ molecule dissociation and O diffuse into the sub-surface. For the O/Cu(100) system, the $22×2R45°$ missing-row and the c(2 × 2) become energetically favored at 0.5 ML oxygen coverage. The saturated Cu bonds, however, make these two structures inert towards further oxidation. We find that the lowest barriers for oxygen dissociation and diffusion on the MRR and c(2 × 2)−0.5 ML are 0.75 eV and 0.97 eV, respectively, which are significantly lower than previously reported. For the O/Cu(110) system, the c(6 × 2) has the highest oxygen concentration prior to Cu₂O formation. The oxygen dissociation barrier is 0.59 eV, but O diffusion into the subsurface has a high barrier of 0.72 eV. Considering that the c(6 × 2) structure cannot fully cover the surface, another structure that we call 1 ML has been investigated. The O₂ dissociation barrier on this phase is 0.31 eV and O atoms can rapidly embed into the subsurface. No ordered structures have been observed on Cu(111) at low temperature, thus we choose the last point on the convex hull before Cu₂O, and points between them to investigate the further oxidation. The O₂ dissociation barrier is high (0.78 eV) on the step + 8O structure, but under experimental conditions, we know that the step + 8O structure cannot cover the entire substrate due to the precise stoichiometry required. We anticipate that at least some of the surface is covered by the 0.5 ML structure, where the dissociation barrier is only 0.28 eV and O can subsequently diffuse spontaneously into the sublayer region.

Comparing the various barriers of oxidation processes on three low-index Cu surfaces: the Cu(111) surfaces yield a fast oxidation at high oxygen coverage; the Cu(110) surface maintains a moderate oxidation rate throughout the process; and for Cu(100) surface, the oxidation rate is low at high oxygen coverage. These computational results are broadly consistent with experimental observations and thus provide mechanistic insight into the subsurface oxidation of copper.

We gratefully acknowledge the NSF under Grant No. DMR-1410335 and the Texas Advanced Computing Center. Xin Lian thanks the China Scholarship Council for financial support.

The authors declare no competing financial interest.