Cu@Phosphorene as a Promising Catalyst for CO₂ to Formic Acid Conversion: A Mechanistic DFT Approach

1. Introduction

As industrialization and population growth increase, the issue of environmental pollution from organic contaminants, coupled with the rapid rise in greenhouse gas emissions, is steadily becoming a global concern [1]. The burning of fossil fuels leads to an increased accumulation of the greenhouse gas, CO₂, which subsequently causes major environmental problems [2]. An excessive amount of CO₂ emissions can lead to harmful changes in the climate, such as global warming and rising sea levels [3]. Therefore, designing and developing effective strategies for capturing and converting CO₂ into valuable chemicals is essential. Such methods can reduce the greenhouse effect of CO₂ and are highly sought after [4,5]. The process of reducing CO₂ levels and converting it into useful substances such as formic acid, methanol, or methane holds considerable value from both economic and environmental perspectives [6,7].

In recent years, single-atom catalysis [8] has brought about a significant transformation in the fields of CO oxidation [9,10,11], CO₂ reduction [3,5,12], and Hydrogen Evolution Reaction [13,14,15]. To achieve optimal catalytic performance and durability, it is essential to fully utilize metal atoms embedded within the support structure [8]. In single-atom catalysts (SACs), each active site is entirely exposed, allowing for maximum atom utilization. As a result, the size of the atoms has a notable impact on the catalyst’s stability and reactivity [16,17,18,19]. Due to their unique characteristics, SACs serve as an innovative type of catalyst with a wide range of applications. Specifically, SACs that feature atomically dispersed metal sites on their supports have the potential to maximize metal utilization [20,21,22].

In the realm of diverse supports, two-dimensional materials emerge as a strategic option due to their larger surface area and multiple anchoring sites for metal atoms, which facilitate the practical fabrication of high-capacity single-atom catalysts [23,24]. Various 2D materials such as MoS₂ [25,26,27,28], Graphene [28,29,30,31,32,33], metal organic frameworks (MOFs) [30,31,33,34], TMDs [31], Graphdiyne [35], and MXenes [36] have received extensive attention in the field of catalysis (e.g., water splitting, CO oxidation, CO₂ reduction, etc.). Phosphorene is a two-dimensional phosphorus atomic layer with a honeycomb-like structure and a high surface-to-volume ratio. Since being exfoliated from black phosphorus in 2014, phosphorene has attracted significant research interest in synthesis, fundamental properties, and potential applications. Current interest in phosphorene is mainly driven by studies of its physical properties, such as tunable bandgap, strong in-plane anisotropy, and high carrier mobility, along with its promising applications in electronics, optoelectronics, spintronics, and energy conversion and storage devices [37]. Numerous investigations, both theoretical and experimental, have explored phosphorene's electrocatalytic properties [38,39]. Phosphorene can serve as an electrocatalyst for hydrogen and oxygen evolution reactions [40].

Phosphorene has several advantages over other 2D supports, such as MXenes and TMDs. The binding between the MXene surface and the captured CO₂ is too strong, which can slow down the hydrogenation process and hinder product desorption [41]. In the case of TMDs, they suffer from poor selectivity, often leading to unwanted by-products, and their stability is compromised under harsh reaction conditions such as high temperatures or high pressures, which are often required for CO₂ reduction. This can degrade the catalyst over time. Additionally, their energy-intensive thermal catalysis and slow kinetics in CO₂ reduction limit overall efficiency [42]. Phosphorene is the most stable allotrope of phosphorus. It is derived from white phosphorus when subjected to a temperature of 200 °C under a substantial pressure of around 1.2GPa [43,44]. Previous studies have demonstrated that copper is a well-established electrocatalyst for reducing CO₂ to C₂ products, as it facilitates intermediate binding, promotes C–C bond formation, supports multiple oxidation states, and exhibits favorable surface dynamics and electronic properties. As a result, Cu-doped phosphorene shows significant potential as an electrocatalyst for CO₂ reduction reactions (CRR) [45,46].

The conversion of CO₂ into formic acid through hydrogenation has attracted significant interest because of its many industrial uses. Formic acid is an important raw material in the organic chemical industry. It is used in the perfume industry, as a mordant in dyeing, as a neutralizer in tanning, and as a sanitizer and preservative in hygiene facilities. Its role as a fundamental chemical in the industry has enabled the production of formate esters, which are key in making various organic byproducts like aldehydes, ketones, carboxylic acids, and amides [47].

Esrafili et al. theoretically investigated the mechanistic pathway for the catalytic hydrogenation of CO₂ over platinum- and nickel-doped graphene [3]. In another study, Sirijaraensre and Limtrakul employed DFT computations and discovered that graphene doped with copper holds potential as a catalyst for converting CO₂ to formic acid using H₂. To date, numerous experimental studies have investigated the chemical and electrochemical conversion of CO₂ into formic acid. However, several aspects, particularly the reaction mechanism and the role of the catalyst, are still unknown. Meanwhile, theoretical studies can help in understanding the transformation process, the methods of adsorption, the levels of activation energies, the properties of transition states, and the role of the catalyst [48,49].

Phosphorene's air sensitivity has historically limited its use in catalytic and electrical devices. To address this issue, many passivation techniques have been developed. Encapsulating phosphorene with materials such as graphene, hexagonal BN (hBN), AlOx, TiOx, or SiO2 creates stable van der Waals heterostructures that maintain phosphorene's electronic structure while preventing oxidation [50]. Covalent P–C bonds formed with aryl diazonium salts, AIBN, or iodonium salts passivate phosphorus’s reactive sites, improving its resistance to oxidation. This stability is confirmed by XPS, FTIR, and SEM analysis [51]. Elemental doping has proven effective in enhancing both the stability and performance of black phosphorus (BP)–based systems. Theoretical density functional theory (DFT) studies predict that dopants like Cu, Te, and S tend to chemisorb on BP's surface or intercalate into its lattice due to surface tension, Coulombic forces, and π–π interactions. These dopants introduce extra electrons that help stabilize reactive phosphorus atoms without changing their electronic properties.

Elemental doping, such as with Cu, Te, or S, improves BP's stability and performance by stabilizing reactive P atoms without altering the electronic structure. DFT predictions have been confirmed through experiments; for example, Cu-doped BP exhibited high mobility and stability [52]. Te-doped BP FETs remained air-stable [53], and Ni nanoparticles prevented degradation by passivating P sites. Both Ni nanoparticles (Ni NPs) and Te-doping significantly enhance BP's durability against environmental damage through surface passivation [54,55]. Recent experimental research has shown promising progress in Cu-based catalysts for CO₂ electroreduction. Li et al. developed facet-engineered Cu₂O@Cu₂S nanostructures that achieved an FE of 43.9% for ethanol at −0.65 V vs. RHE in flow cells [56]. Although these methods involve surface modification and heterostructuring, they still face challenges related to maintaining active site stability during long-term electrochemical operation. In contrast, our Cu@phosphorene single-atom catalyst (SAC) design offers atomic-level Cu dispersion and robust Cu-phosphorene coordination, leading to enhanced catalytic durability.

Recent DFT studies, such as Zhang et al. [57], have shown that Cu-doped phosphorene favors CO₂ reduction over hydrogen evolution by stabilizing key intermediates (*COOH and *OCHO) more than *H. Their work also indicated improved selectivity toward HCOOH due to easier desorption compared to CO. These findings support our choice of Cu@Phosphorene as a promising catalyst. Therefore, in the present work, we have used density functional theory (DFT) for the mechanistic study of CO₂ reduction into formic acid using Cu@Phosphorene as a SAC and studied its stability, interaction energies, and catalytic performance for the CO₂ reduction reaction.

2. Computational Methods

All calculations for this study were performed using the Gaussian 16 (Revision B.01) suite [58]. The optimized molecules were visualized with CYLView [59]. The calculations employed Adamo’s hybrid version of the Perdew, Burke, and Ernzerhof functional (PBE0) [60,61] and Ahlrichs' triple-ζ basis set def2-TZVP [62]. Grimme’s empirical D3 correction with Becke-Johnston damping (D3BJ) [63,64,65] was applied to account for dispersion effects accurately. Frequencies were computed to confirm that the optimized structures are true minimum energy structures on the potential energy surface, as indicated by the absence of imaginary frequencies. The Berny algorithm, implemented in Gaussian 16, was used to identify the geometries of transition states.

A phosphorene sheet was modeled as a puckered honeycomb-like structure consisting of 46 phosphorus atoms, which has been shown to be sufficiently large to accurately represent the local surface environment and electronic properties of 2D phosphorene while remaining computationally feasible. Similar cluster sizes have been effectively employed in our previous DFT studies to investigate adsorption and catalytic processes on phosphorene-based systems [10,66]. Formation energy (E-form) for each species was calculated using the formula:

$${\mathit{E}}_{\mathit{a}\mathit{d}\mathit{s}}={\mathit{E}}_{\mathit{C}\mathit{u}-\mathit{P}\mathit{h}\mathit{o}\mathit{s}\mathit{p}\mathit{h}\mathit{o}\mathit{r}\mathit{e}\mathit{n}\mathit{e}+\mathit{g}\mathit{a}\mathit{s}\mathit{m}\mathit{o}\mathit{l}\mathit{e}\mathit{c}\mathit{u}\mathit{l}\mathit{e}}-{\mathit{E}}_{\mathit{C}\mathit{u}-\mathit{P}\mathit{h}\mathit{o}\mathit{s}\mathit{p}\mathit{h}\mathit{o}\mathit{r}\mathit{e}\mathit{n}\mathit{e}}-{\mathit{E}}_{\mathit{g}\mathit{a}\mathit{s}\mathit{m}\mathit{o}\mathit{l}\mathit{e}\mathit{c}\mathit{u}\mathit{l}\mathit{e}}$$

Where E_{Cu-Phosphorene} represents the electronic energy of metal-embedded phosphorene substrate, E_ads is the metal atom energy acting as SAC, E_Phosphorene represents the energy of the pristine phosphorene sheet. Moreover, the activation energy (E_a) of a reaction has been calculated as:

$${\mathit{E}}_{\mathit{a}}={\mathit{E}}_{\mathit{T}\mathit{S}}-{\mathit{E}}_{\mathit{I}\mathit{S}}$$

The difference between the TS and IS (initial state) of a reaction is termed the activation energy of that reaction. E_ads indicates the adsorption energy of different species, i.e., CO₂, H₂, and HCOOH, after they adsorb over the catalyst’s surface. A negative value of E_ads manifests that the reaction has been completed exothermically, while an endothermic reaction shows a positive value of E_ads.

The charge analysis was conducted employing the same theoretical approach. i.e., PBE0-D3BJ/def2-TZVP with Hirshfeld charge analysis [67] (NBO analysis) and the values of charges are given in Figure 9 and Figure 10. The Multiwfn software (version 3.80) [68] was used to perform and analyze non-covalent interaction (NCI) and Quantum Theory of Atoms in Molecules (QTAIM). The resulting graphs were then visualized and interpreted using the VMD software.[69]

A kinetic study was performed at temperatures up to 225°C (~500K) to investigate the relationship between the rate constant and temperature, using the equation below[70].

$$K^{TST}={\displaystyle \frac{K_{b}T}{h}}\times {\displaystyle \frac{{-∆G}_a}{K_{b}T}}$$

Where kb, T, and h are the Boltzmann constant (8.6173303 ×10^-5 eV K ^-1), the temperature of the reaction, and Planck's constant (6.582119514 × 10^-16 eV s), respectively, while ΔG is the activation energy barrier in Gibbs free energy. AIMD simulations were performed for 10 ps at 500K with a 0.5 fs step size using ORCA 5.0.3 [71] at the same level of theory as for optimizations, as studied by Zhengy and Gao et al. [72,73] to assess the catalyst's stability under experimental conditions. All calculations in this study were carried out in the gas phase to provide fundamental insights into the intrinsic catalytic behavior of Cu@Phosphorene. This approach is widely adopted for establishing reaction mechanisms and active site characteristics. However, we acknowledge that real catalytic environments involve solvent, pH, and pressure effects, which may influence the reaction energetics. Future studies incorporating solvation models and pH-dependent pathways will help extend these findings to more realistic conditions.

3. Results and Discussion

3.1. Geometry Optimization of Cu Embedded Phosphorene

This study aims to reveal the catalytic mechanism of CO₂ reduction on a Cu-anchored phosphorene surface (referred to as Cu@Phosphorene). To build a detailed understanding, we first examine the adsorption properties and interaction energies of CO₂ and H₂ molecules on the Cu@Phosphorene system. The optimized structure of Cu@Phosphorene is shown in Figure 1.

To evaluate the structural stability, the Cu atom was placed at various potential adsorption locations on the phosphorene surface—specifically, the cavity, middle, and terminal sites—as previously described in our work [10]. Among these, the cavity site proved to be the most favorable, showing the most negative interaction energy, which indicates high thermodynamic stability.

Spin state analysis showed that the Cu atom favors a doublet ground state, with an interaction energy of −1.89 eV, compared to the less favorable quartet state, which has a positive interaction energy of 0.64 eV. This confirms that a single Cu atom on phosphorene is most stable in the doublet spin state.

To further clarify the bonding between the Cu atom and the phosphorene substrate, we used Bader's quantum theory of atoms in molecules (QTAIM). This analysis provided essential electronic descriptors, including electron density (ρ), Laplacian of electron density (∇²ρ), potential energy V(r), Lagrangian kinetic energy G(r), total energy density H(r), and bond interaction energies between Cu and P atoms. These parameters give valuable insight into the electronic structure and bonding nature of the catalytic interface.

There are three bond critical points (BCPs) between Cu metal and phosphorene support, as mentioned below.

i.

CP 88 (Cu--P27)

ii.

CP 95 (Cu--P10)

iii.

CP 101 (Cu--P23).

Table 1. QTAIM analysis parameters.

QTAIM analysis BCPs (Cu@Phosphorene)	Density (ρ)	Laplacian (∇2ρ)	Lagrangian kinetic energy G(r)	Potential energy V(r)	Total energy density H(r)	Bond interaction energies
Cu-P27	0.06124	0.15367	0.05039	−0.06237	−0.01197	−19.56
Cu-P10	0.07704	0.2149	0.07078	−0.08783	−0.01705	−27.56
Cu-P23	0.05922	0.1468	0.048119	−0.05952	−0.01140	−18.67

¹ Tables may have a footer.

Table 1. QTAIM analysis parameters.

QTAIM analysis BCPs (Cu@Phosphorene)	Density (ρ)	Laplacian (∇2ρ)	Lagrangian kinetic energy G(r)	Potential energy V(r)	Total energy density H(r)	Bond interaction energies
Cu-P27	0.06124	0.15367	0.05039	−0.06237	−0.01197	−19.56
Cu-P10	0.07704	0.2149	0.07078	−0.08783	−0.01705	−27.56
Cu-P23	0.05922	0.1468	0.048119	−0.05952	−0.01140	−18.67

¹ Tables may have a footer.

The values of these parameters show that there is a partial covalent bond present between Cu and the phosphorene sheet. Hence, this type of bonding shows that the Cu@Phosphorene sheet is stable. From all of these critical points, Cu-P10 shows the strongest interaction as shown in Figure 2.

3.2. Adsorption of Different Species on Cu-Embedded Phosphorene

The adsorption behavior of different species, including CO₂, H₂, CO₂+H₂, 2H_2, and 2CO₂ on Cu@Phosphorene has been studied to confirm the competitive adsorption of different species and to confirm that the adsorption of the species involved in the mechanism will be favorable energetically. Table 2 presents the bond lengths and adsorption energies of the species, which facilitate the estimation of their adsorption behavior over a metal catalyst for CO₂ reduction.

In Figure 3 (complex A), the adsorption of the H₂ molecule over Cu@Phosphorene is evaluated to identify the most stable configuration. The calculated (E_ads) value of H₂ over Cu@Phosphorene is −0.06 eV Figure 3. The bond length of the adsorbed H₂ molecule has significantly increased (from 0.60 to 0.79 Å). In complex B, the adsorption of two H₂ molecules on the Cu@Phosphorene substrate is investigated. The adsorption energy of this complex is −0.15 eV. One of the H₂ bond lengths increases from 0.60 to 0.79Å while the other increases up to 0.76Å. In complex C, the adsorption of individual CO₂ and H₂ molecules on a Cu@Phosphorene substrate is examined as part of the study on CO₂ reduction. The adsorption energy of H₂-CO₂ when they were adsorbed simultaneously on Cu@Phosphorene was −0.26 eV, and the bond distances between Cu@Phosphorene and H₂ and CO₂ were 2.26 Å and 2.85 Å, respectively. The adsorption energy of CO₂ on Cu@Phosphorene was −0.19 eV, with the bond distance between adsorbed CO₂ and Cu@Phosphorene as 2.28 Å (Figure 3, complex D).

On the other hand, when two CO₂ molecules were simultaneously adsorbed on the surface of the Cu@Phosphorene, the calculated adsorption energy was −0.59 eV, and their distances from Cu@Phosphorene were 2.26 and 2.27 Å, respectively. These adsorption energy values depict that the adsorption of two carbon dioxide molecules simultaneously adsorbed on the surface of Cu@Phosphorene appear more favorable than that of other reacting species. These results also indicated that CO₂ gets adsorbed on the catalytic surface in the initial state, and it also helped to avoid the potential hindrance of the active sites (Figure 3, complex E).

3.3. Reaction Mechanism for CO₂ Reduction Using Cu@Phosphorene

In previous studies, two mechanisms have been employed for the CO₂ reduction reaction: the Biomolecular mechanism (BM) and the Ter-molecular mechanism. There are a number of reports of CO₂ reduction, such as Pt/Ni atoms embedded in graphene support [74] and catalytic reduction of CO₂ to HCOOH using Ti-graphene nano-flakes [3]. In this study, we focused specifically on the bimolecular and trimolecular pathways for CO₂ reduction to formic acid because these routes represent essential, thermodynamically feasible molecular reactions on the Cu@Phosphorene surface under non-electrochemical conditions. Although other mechanistic pathways—such as electrochemical CO₂ reduction and side reactions like CO formation—are also relevant to real-world catalysis, they add complexities such as applied potential, solvation effects, and competing proton-coupled electron transfer steps. These aspects are outside the current scope but are important for a broader understanding of catalytic selectivity and will be explored in future work. Our decision to concentrate on BM and TER mechanisms allows for a clearer understanding of the inherent reactivity and selectivity of the Cu@Phosphorene catalyst toward HCOOH formation.

3.3.1. CO2 Reduction via Bi-Molecular (BM) Mechanism

The following is a proposed bimolecular reaction mechanism for the catalytic reduction of CO₂ and the production of formic acid:

$$\mathit{C}{\mathit{O}}_2+\mathit{C}\mathit{u}@\mathit{P}\to {\mathit{C}\mathit{O}}_2····\mathbf{C}\mathbf{u}@P$$

$${\mathit{H}}_2+{\mathit{C}\mathit{O}}_2\dots .\mathit{C}\mathit{u}@\mathit{P}\to {\mathit{C}\mathit{O}}_2/{\mathit{H}}_2····\mathbf{C}\mathbf{u}@\mathbf{P} \mathbf{(IS1/TS1)}$$

$${\mathit{C}\mathit{O}}_2/{\mathit{H}}_2····\mathbf{C}\mathbf{u}@\mathbf{P}\to \mathit{H}-\mathit{C}\mathit{u}@\mathit{P}-\mathit{C}\mathit{O}\mathit{O}\mathit{H} \mathbf{(IM1/TS2)}$$

$$\mathit{H}-\mathit{C}\mathit{u}@\mathit{P}-\mathit{C}\mathit{O}\mathit{O}\mathit{H}\to \mathit{H}\mathit{C}\mathit{O}\mathit{O}\mathit{H}····\mathbf{C}\mathbf{u}@\mathbf{P} \mathbf{(FS2)}$$

The optimized geometries and energy profiles of the initial state (IS), transition state (TS), and final state over Cu@P are depicted in Figure 4.

According to the proposed BM mechanism, the reduction of CO₂ to formic acid occurs in two steps: first, an intermediate called OCOH is formed, and then this intermediate is transformed into formic acid. The configuration of bimolecular adsorption of carbon dioxide and hydrogen gases on the Cu-adsorbed phosphorene substrate is the initial step in the reduction of CO₂. As previously mentioned, CO₂ first adsorbs on the Cu@phosphorene surface, followed by H₂ co-adsorption, depending on the resulting adsorption energies, as shown in Table 2. This configuration is chosen as the starting point for a reaction, i.e., (IS-1). In this state, the bond length of an adsorbed H₂ molecule is 0.80 Å, which increases to 1.10 Å in TS1. The H₂ molecule approaches the O of the CO₂ molecule on the Cu@Phosphorene in TS1, facing a barrier height of 1.55 eV; vibrational frequency analysis of TS1 shows one negative frequency of −1113.85 cm⁻¹, indicating the cleavage of H1-H2 and formation of the O1-H1 bond. Additionally, the carbon atom in CO₂ shifts from sp hybridization to sp², and the O−C−O angle decreases from 177.6° in IS1 to 138.8° in TS-1. This step results in the hydrogenation of the O atom in the CO₂ molecule and the formation of the OCOH intermediate, releasing 0.44 eV of energy from TS1 to the intermediate (IM). The formation of IM is challenging because it involves crossing a moderately large energy barrier of 1.55 eV (Figure 4).

In the IM1 stage, the bond between H-H breaks, and the first Hydrogen atom, H1, attaches to the oxygen of CO₂ to form an OCOH intermediate. The distance between the oxygen of CO₂ and the attached H1 in TS1 is 1.16 Å, which decreases to 0.974 Å in IM1. Additionally, the distance between the second Hydrogen, H₂, and the Cu metal adsorbed on phosphorene changes from 1.58 Å to 1.59 Å from TS1 to IM1. In the second step, H₂ shifts from the Cu@phosphorene to the carbon of OCOH via TS2. During this process, the H₂–C distance decreases from 2.65 Å to 1.81 Å, and the Cu-H distance shortens from 1.59 Å to 1.55 Å. Frequency analysis reveals a single imaginary frequency of −599.02 cm⁻¹, confirming the formation of TS2. This vibration corresponds to both the dissociation of the Cu-H bond and the formation of the C-H bond. The activation energy for this step is 0.44 eV. Ultimately, HCOOH forms on the Cu atom of Cu@phosphorene, releasing 1.74 eV of energy. The energy of the FS2 is −0.31 eV relative to IS-1. Furthermore, the transition from IM1 to FS2 is an exothermic process.

Figure 5. Optimized geometries of the reaction configurations via a Bi-molecular mechanism for the reduction of CO₂ over Cu@Phosphorene. All distances are in Å. The yellow, grey, red, white, and brown colors represent the Cu, C, O, H, and P atoms, respectively.

Figure 5. Optimized geometries of the reaction configurations via a Bi-molecular mechanism for the reduction of CO₂ over Cu@Phosphorene. All distances are in Å. The yellow, grey, red, white, and brown colors represent the Cu, C, O, H, and P atoms, respectively.

3.3.2. CO₂ Reduction via Ter-Molecular (TER) Mechanism

Catalytic hydrogenation of carbon dioxide and the formation of formic acid are thought to occur via the Ter-molecular reaction process as explained below:

$$2{\mathit{C}\mathit{O}}_2+\mathit{C}\mathit{u}@\mathit{P}\to 2{\mathit{C}\mathit{O}}_2······\mathit{C}\mathit{u}@\mathit{P}$$

$${\mathit{H}}_2+2{\mathit{C}\mathit{O}}_2······\mathit{C}\mathit{u}@\mathit{P}\to 2{\mathit{C}\mathit{O}}_2/{\mathit{H}}_2······\mathit{C}\mathit{u}@\mathit{P} \mathbf{(IS3/TS3)}$$

$$2{\mathit{C}\mathit{O}}_2/{\mathit{H}}_2······\mathit{C}\mathit{u}@\mathit{P}\to \mathit{C}\mathit{O}\mathit{O}\mathit{H}···\mathit{C}\mathit{u}@\mathit{P}···\mathit{C}\mathit{O}\mathit{O}\mathit{H} \mathbf{(FS3)}$$

$$\mathit{C}\mathit{O}\mathit{O}\mathit{H}···\mathit{C}\mathit{u}@\mathit{P}···\mathit{C}\mathit{O}\mathit{O}\mathit{H}+{\mathit{H}}_2\to {\mathit{H}}_2-\mathit{C}\mathit{O}\mathit{O}\mathit{H}···\mathit{C}\mathit{u}@\mathit{P}···\mathit{C}\mathit{O}\mathit{O}\mathit{H} \mathbf{(IM4)}$$

$${\mathit{H}}_2-\mathit{C}\mathit{O}\mathit{O}\mathit{H}···\mathit{C}\mathit{u}@\mathit{P}···\mathit{C}\mathit{O}\mathit{O}\mathit{H}\to \mathit{H}\mathit{C}\mathit{O}\mathit{O}\mathit{H}···\mathit{C}\mathit{u}@\mathit{P}···\mathit{C}\mathit{O}\mathit{O}\mathit{H}/\mathit{H} \mathbf{(IM5/TS5)}$$

$$\mathit{H}\mathit{C}\mathit{O}\mathit{O}\mathit{H}···\mathit{C}\mathit{u}@\mathit{P}···\mathit{C}\mathit{O}\mathit{O}\mathit{H}/\mathit{H}\to \mathit{H}\mathit{C}\mathit{O}\mathit{O}\mathit{H}···\mathit{C}\mathit{u}@\mathit{P}···\mathit{H}\mathit{C}\mathit{O}\mathit{O}\mathit{H} \mathbf{(FS5)}$$

The ter-molecular mechanism (Figure 6) proposed here involves the co-adsorption of two carbon dioxide molecules on the Cu@Phosphorene followed by the adsorption of H₂. Then, the pre-adsorbed CO₂ molecules activate the H₂ molecule. The energy profile associated with the ter-molecular process on the Cu@Phosphorene sheet is displayed in Figure 6.

The initial state of this process (IS3) involves H₂ and two CO₂ molecules co-adsorbed on a Cu atom of the Cu@Phosphorene sheet. The Cu-C distances on both sides are 3.05 and 3.03 Å, respectively. The H-H bond is stretched by 0.82 Å, allowing hydrogen atoms to move closer to the co-adsorbed CO₂ molecules. The predicted activation energy to form the transition state TS3 is 1.39 eV.

Figure 6. The energy profile for a Ter-molecular mechanism of CO₂ reduction over Cu@Phosphorene. The blue colored values indicate the energy relative to the initial state while pink colored values show the energy difference of each step.

Figure 6. The energy profile for a Ter-molecular mechanism of CO₂ reduction over Cu@Phosphorene. The blue colored values indicate the energy relative to the initial state while pink colored values show the energy difference of each step.

Figure 7. Optimized geometries of the reaction configurations via a Ter-molecular mechanism for the reduction of CO₂ over Cu@Phosphorene. All distances are in Å. The yellow, grey, red, white, and brown colours represent the Cu, C, O, H, and P atoms, respectively.

Figure 7. Optimized geometries of the reaction configurations via a Ter-molecular mechanism for the reduction of CO₂ over Cu@Phosphorene. All distances are in Å. The yellow, grey, red, white, and brown colours represent the Cu, C, O, H, and P atoms, respectively.

At FS3, two COOH fragments are generated over a Cu@Phosphorene sheet, with Cu-C lengths of 2.46 Å and 2.48 Å. Additionally, the formation of these intermediates is thermodynamically feasible and exothermic at ambient temperature. The second hydrogen molecule then adsorbs onto the Cu atom, forming the first intermediate structure, i.e., IM4 in this mechanism. The hydrogen (H₂) molecule dissociates, and the H3 atom moves toward the OCOH intermediate's carbon atom, creating the first formic acid molecule. It should be noted that the predicted energy of activation for the hydrogenation of the OCOH group is less than that of the bimolecular route. Finally, the H4 atom reacts with the second OCOH molecule, resulting in the formation of a second formic acid molecule on the Cu@Phosphorene sheet. The energy of the FS5 state is −3.95 eV compared to the IS3. Moreover, the step from IM5 to FS5 is exothermic and can proceed under normal temperature and pressure conditions.

Based on the analysis of activation and desorption energies, the termolecular mechanism is more suitable for our research. The low activation energy of 1.39 eV in the termolecular process facilitates the initiation of the reaction, resulting in faster overall reaction rates. While the bimolecular mechanism has a lower desorption energy of 0.68 eV, allowing for quicker product release, its higher activation energy of 1.54 eV slows down the process of reaching the transition state, reducing overall efficiency. Therefore, the termolecular mechanism offers better catalytic efficiency and faster reaction kinetics, making it the preferred pathway for CO₂ reduction in this research. The desorption energies of HCOOH suggest that the product, HCOOH (formic acid), can be easily released from the catalyst surface. This reduces the risk of catalyst poisoning since CO₂ initially adsorbs and then regenerates the active site when converted to HCOOH.

We analyze the stability of key intermediates and the potential for catalyst deactivation by calculating the interaction (adsorption) energies of the intermediates involved in both the bimolecular and termolecular CO₂ reduction pathways on Cu@Phosphorene. In the bimolecular mechanism, the intermediate IM1, with an interaction energy of −0.85 eV, adsorbs to the Cu@Phosphorene surface. In the termolecular mechanism, intermediates IM4 and IM5, with interaction energies of −0.95 eV and −1.04 eV, respectively, also adsorb to the Cu@Phosphorene surface. These values indicate stable adsorption and effective interaction with the catalyst surface, which are within an optimal range and are strong enough to support reaction steps, but not so strong as to promote premature product desorption or cause catalyst poisoning.

Our findings are further supported by the study of Zhang et al. [57], who examined the selectivity between the CO₂ reduction reaction (CRR) and the hydrogen evolution reaction (HER) by calculating the free energy changes for forming key intermediates (*COOH, *OCHO, and *H). Their results showed that forming *COOH or *OCHO is more thermodynamically favorable than *H, indicating increased selectivity for CRR over HER with Cu doping. Additionally, the same study analyzed product desorption behavior and found that HCOOH desorbs more easily from Cu-doped phosphorene than CO, while pristine phosphorene favors CO desorption. These results support the conclusion that Cu doping not only promotes CRR over HER but also enhances selectivity toward HCOOH. Although these electrochemical pathways were not explored in the present work, we reference them here to provide a more comprehensive picture of the catalytic potential of Cu-doped phosphorene.

To evaluate potential side reactions and selectivity, we investigated the kinetic feasibility of the CO formation pathway. Specifically, we computed the transition state for the rate-determining step, following the mechanism outlined in Ref. [12]. The calculated activation energy was found to be 7.01 eV, which is considerably high. Vibrational analysis confirmed the nature of the transition state, showing a single imaginary frequency of –1147.91 cm⁻¹ corresponding to the hydrogen-transfer motion. This result suggests that the CO formation pathway is kinetically unfavourable and unlikely to proceed under the studied conditions, even prior to reaching the *CO intermediate.

Additionally, we evaluated the adsorption energies of the relevant products. CO binds to the Cu@Phosphorene surface with a moderate interaction energy of –0.76 eV, whereas HCOOH shows a stronger adsorption energy of –0.94 eV. This thermodynamic preference further supports the selectivity toward formic acid, indicating that HCOOH is both kinetically and thermodynamically favoured over CO on the Cu@Phosphorene catalyst.

So, by taking in consideration the high barrier to CO formation, the weak interaction of CO, the thermodynamic instability of *H as reported by Zhang et al. [57], and the strong stabilization of HCOOH confirm that Cu@Phosphorene exhibits high selectivity for formic acid while minimizing undesirable side products.

3.4. Ab Initio Molecular Dynamics (AIMD)

Figure 8 shows snapshots from AIMD simulations conducted for 10 ps with a 0.5 fs time step at 500 K for the Cu@Phosphorene system to evaluate its thermal stability. The simulations reveal that the Cu@Phosphorene complex remained structurally stable throughout. The Cu atom initially settled on the phosphorene sheet, with minimal migration observed during the simulation. AIMD snapshots were taken at 0, 0.25, 5, and 10 ps. In the Cu@Phosphorene complex, the Cu atom moved slightly during the early steps but quickly stabilized, maintaining its position after 0.25 ps. This small fluctuation suggests that the Cu atom is quite stable on the phosphorene surface at such high temperatures. This finding confirms that the Cu@Phosphorene complex can maintain its structure and strong Cu-P bonds even at elevated temperatures, demonstrating high thermal stability. No noticeable deformation was observed as the temperature increased.

3.5. NBO Analysis

Literature shows that charge transfer occurs when two moieties interact due to differences in electronegativity [10]. Charge transfer analysis from NBO calculations revealed that Cu@Phosphorene acts as an electron donor with a charge of 0.417 e− (here Cu also gains electrons from the surface atoms of Phosphorene).

In the bimolecular mechanism, NBO analysis is performed to examine charge differences during the absorption or desorption of different moieties from the catalyst (phosphorene) surface, which causes charge redistribution in the catalyst. As shown in Figure 9, electron density changes at each step. In IS1, the positive charge on Cu metal indicates its electropositive (electron-donating) nature, suggesting that electron density transfers from the metal to the phosphorene sheet. When the hydrogen molecule approaches Cu in TS1, the electron density of Cu decreases from 0.165 e− to 0.021 e−. This reduction results from charge transfer from the metal to hydrogen. Meanwhile, another hydrogen interacts with the oxygen of CO₂, indicated by an increase in the electronegative character of O from −0.514 e− to −0.657 e−. After hydrogen adsorbs onto the metal in IS2, charge redistribution occurs, and the negative charge on the metal shows that electron density transfers from phosphorene to the metal, then further to the hydrogen atom. Additionally, the increase in the electronegative character of Cu from −0.124 to −0.210 e− in TS2 is due to charge redistribution as electropositive carbon approaches the catalyst. In the final state (FS2), the positive charge on Cu indicates hydrogen desorption from its surface, and the product is released as formic acid.

Figure 9. NBO Analysis of CO₂ Reduction via Bi-molecular Mechanism. Here, the yellow, grey, red, white, and brown colours represent the Cu, C, O, H, and P atoms, respectively.

Figure 9. NBO Analysis of CO₂ Reduction via Bi-molecular Mechanism. Here, the yellow, grey, red, white, and brown colours represent the Cu, C, O, H, and P atoms, respectively.

In the IS3 state, as two molecules of carbon dioxide approach, the charge on the Cu atom is 0.095 e−, indicating charge transfer from CO₂ to the Cu atom (Figure 10). Meanwhile, in TS3, this charge decreases to +0.037 e−, showing that Cu is gaining electrons. In IS3, the charges on the oxygen atoms directly attached to Cu are −0.529 e− and −0.620 e−, but in TS3, these negative charges increase to −0.596 and −0.655, respectively, suggesting that bond formation may occur between CO₂ and the Cu atom. In FS3, the charge on Cu is +0.471 e− because it is directly bonded to an electronegative atom. In the next step, IM4, another H₂ molecule is introduced, and the charge on Cu is +0.357e−, indicating charge transfer toward Cu, meaning it is gaining electrons. In IM5, this charge increases to −0.153 e−, with Cu directly attached to H and OCOH, as H transfers electron density to Cu due to its electropositive nature. The bond between Cu and H breaks, and COOH interacts with H and HCOOH. In the final step, the charge on Cu is +0.385 e−, which shows that HCOOH desorbs from the catalyst surface.

Figure 10. NBO Analysis of CO₂ Reduction via Ter-molecular Mechanism. Here, the yellow, grey, red, white, and brown colours represent the Cu, C, O, H, and P atoms, respectively.

Figure 10. NBO Analysis of CO₂ Reduction via Ter-molecular Mechanism. Here, the yellow, grey, red, white, and brown colours represent the Cu, C, O, H, and P atoms, respectively.

3.6. NCI Analysis

NCI (non-covalent interactions) analysis was used to examine the nature and importance of weak interactions. The NCI plot is based on two variables: electron density (ρ(λ2)), which ranges from 0.05 to –0.05, and the reduced density gradient (RDG).

Hydrogen bonding is represented by the negative ρ(λ2) value which is blue, while steric hindrance is defined by the positive ρ(λ2) value (red colour). The green region between hydrogen bonding and repulsive forces represents the weak van der Waals interactions and ranges from 0.01 to -0.01. The NCI approach of the several states that are involved in the bi-molecular (i.e., IS1, TS1, IM1, TS2, and FS2) and the ter-molecular mechanism (i.e., IS3, TS3, FS3, IM4, TS4, 1M5, TS5, and FS5) is given in Figure S1.

Cu@Phosphorene acts as a heterogeneous catalyst, with the reaction occurring through weak surface interactions. In the bi-molecular mechanism, based on the 3D NCI plots (Figure S1), non-covalent interactions such as hydrogen bonding and van der Waals forces facilitate the approach of CO₂ and H₂ molecules to Cu@Phosphorene at the IS1 state. At TS1, a hydrogen molecule is activated, and one hydrogen atom exhibits a strong interaction with the CO₂ molecule. This is illustrated in the NCI scatter plot, where the electron density increases compared to IS1, indicating that TS1 is less stable. An OCOH moiety develops on the catalytic surface in 1M1. In TS2, the 3D NCI plot displays a blue ring-like iso-surface between the OCOH group and the hydrogen atom, indicating bonding. The NCI scatter plot also highlights van der Waals interactions between OCOH and Cu@Phosphorene, as well as hydrogen bonding between OCOH and the hydrogen atom. In FS2, HCOOH forms, shown by a prominent disc-shaped green iso-surface between HCOOH and Cu@Phosphorene in the 3D NCI plot. The NCI scatter plot further indicates that FS2 is more stable than TS2.

According to the 3D NCI plots (Figures S2 and S3), in the ter-molecular mechanism, non-covalent interactions such as hydrogen bonding and van der Waals forces cause both H₂ and two CO₂ molecules to approach the Cu@Phosphorene at the IS3 state. At TS3, the hydrogen molecule is activated, and hydrogen atoms form bonds with the Cu atom and CO₂ molecules. NCI scatter plots reveal van der Waals interactions with surface atoms and increased electron density, making this state less stable than IS3. At FS3, two OCOH groups are formed, with green disc-like iso-surfaces appearing on both sides, indicating bonding between Cu of the Cu@Phosphorene and OCOH. In 2D NCI, the electron densities are in the 0.01 to -0.01 (ρ(λ₂)) range in the RDG plot of TS3 and FS3.

At IM4, another H₂ molecule is added, leading to interactions between H₂ and OCOH moieties. At TS4, hydrogen atoms move apart, now displaying a prominent blue iso-surface between the H atom and OCOH groups. At IM5, one formic acid molecule is formed, and 2D NCI of FS5 indicates sufficient van der Waals forces and hydrogen bonding between HCOOH and the sheet, stabilizing the product. Conversely, there is minimal rotation of OCOH, and the H4 atom interacts with OCOH, shown by a strong blue disc indicating hydrogen bonding. At TS5, the 2D NCI-RDG plot reveals high electron density in areas of hydrogen bonding and van der Waals interactions, reducing the stability of the transition state. In FS5, steric hindrance appears along with non-bonding interactions seen in both 2D NCI-RDG and 3D plots, promoting the desorption of the second HCOOH.

3.7. Kinetic Analysis

The kinetic analysis of the reaction has been done by using the following equation:

$${\mathbf{k}}^{\mathit{T}\mathit{S}\mathit{T}}={\displaystyle \frac{{\mathit{k}}_{\mathit{B}}\mathit{T}}{\mathit{h}}}\times \mathbf{exp}\left({\displaystyle \frac{-\mathit{\Delta }{\mathit{G}}_{\mathit{a}}}{{\mathit{k}}_{\mathit{B}}\mathit{T}}}\right)$$

The reaction’s rate is calculated over the temperature range from 298K to 500K. By applying traditional transition-state theory (TST), k^TST (rate constant of reaction) is determined. Then the natural logarithm of k^TST is taken, and temperature (1000/K) is calculated.

Temperature has been shown to have a major impact on the k^TST, as the temperature rises, the value of the rate increases constantly. The rate constant has a smaller value at room temperature, but when there is an increase in temperature, its value shows an increase, indicating that increasing the temperature boosts the reduction of CO_2, as shown in Figure 11.

Our Cu@Phosphorene single-atom catalyst (SAC) exhibits superior structural stability compared to existing studies. Many reported catalysts, such as metal-oxide or doped-carbon systems like Cu₂O@CuBr or Cu-based SACs [75], undergo structural degradation under similar or even less demanding operating conditions. Lv, Z., et al. used Cu₂O_Br as a catalyst for CO₂ reduction and observed higher Faradaic efficiencies (FEs) for CO and C₂H₆ production with Cu₂O_Br at 0.4 V and 0.6 V vs. RHE. [76] Similarly, Liu, G., et al. used a CuBi2O4-Bi2O3 photocathode for CO₂ reduction, achieving a Faradaic efficiency (FE) for CO (FEC) that reaches up to 86.03% [77]. However, our study focuses on the selective reduction of CO₂ to formic acid. Ab initio molecular dynamics (AIMD) simulations at 500K show that Cu remains firmly anchored to the phosphorene even at this high temperature. This indicates a very high level of thermal stability and strong metal-support interaction in the Cu@Phosphorene catalytic system, making it a highly promising catalyst.

4. Conclusions

The present study highlights the promising potential of Cu@Phosphorene as an efficient and thermally stable single-atom catalyst for CO₂ reduction. Density functional theory calculations and mechanistic insights show that the termolecular pathway is energetically more favorable for forming formic acid. Complementary analyses, including NBO, QTAIM, and NCI, confirm strong interactions and effective charge transfer during catalysis. Additionally, AIMD simulations verify the thermal stability of the Cu@Phosphorene complex, further supporting its practical use. These findings not only expand the understanding of CO₂ conversion mechanisms at the atomic level but also aid in developing cost-effective catalysts for sustainable carbon management, aligning with global climate action goals.