Computational study of H₂ generation from BH₃ and BH₃⁻ with HO˙ radical

Abstract

Hydrogen formation via reactions involving boron hydrides has attracted increasing attention due to its relevance to clean energy and radical chemistry. In this study, density functional theory at the M06-2X/6-311++G(d,p) level, combined with transition state theory (TST) including Wigner tunneling corrections, is employed to investigate the mechanisms and kinetics of hydrogen generation from reactions of neutral BH₃ and anionic BH₃⁻ with hydroxyl radicals (HO˙) in both the gas phase and aqueous solution. The results reveal reaction pathways involving stabilized pre-reactive intermediates followed by hydrogen transfer steps, with stronger donor acceptor interactions in the anionic system facilitating H–H bond formation. Kinetic analysis shows pronounced temperature dependence and tunneling effects. In the gas phase, a negative temperature dependence is observed, arising from a dominant pre-equilibrium associated with a strongly stabilized intermediate. In aqueous solution, the overall kinetics are governed by the combined effects of pH-dependent speciation, diffusion, and temperature. These findings provide key insights into the mechanisms of hydrogen formation from boron hydrides and identify the critical factors controlling reactivity in both gas and liquid environments.

---

1. Introduction

The growing demand for sustainable and clean energy has intensified interest in alternative energy carriers, among which hydrogen (H₂) is widely recognized as a promising candidate due to its high gravimetric energy density and environmentally benign combustion products.^1–3 In this context, the development of efficient and controllable hydrogen generation strategies from chemical hydrides has attracted considerable attention, as these materials offer high hydrogen storage capacities and tunable reactivity under relatively mild conditions.^4–7 Among them, boron-based hydrides have emerged as particularly attractive systems owing to their structural versatility, electron-deficient nature, and favorable thermodynamic profiles for hydrogen release.^8,9

Borane (BH₃) and its derivatives represent prototypical electron-deficient compounds, in which the vacant p orbital on the central boron atom enables strong interactions with nucleophiles and radical species.^10,11 This intrinsic feature renders boron hydrides highly reactive toward oxidative and radical environments, making them potential candidates for radical-mediated hydrogen evolution.^12,13 In particular, hydroxyl radicals (HO˙), as one of the most reactive transient species in both aqueous and atmospheric chemistry, can efficiently interact with boron hydrides through hydrogen abstraction or transfer pathways.^14–16 Despite this potential, the detailed molecular mechanisms governing hydrogen formation in the BH₃ + HO˙ system, including the nature of hydrogen transfer steps and the role of transient intermediates, remain incompletely understood.

From a mechanistic perspective, hydrogen evolution from boron hydrides may proceed via multiple pathways, including stepwise hydrogen abstraction and coupled proton–electron transfer processes.^17,18 Among these, double hydrogen atom transfer (DHAT) has recently been identified as an efficient route for H₂ formation in radical-mediated systems, involving the transfer of two hydrogen atoms either concertedly or in a closely coupled sequence. However, the existence, feasibility, and relative importance of DHAT pathways in the BH₃ + HO˙ system have not been systematically explored. Furthermore, the role of charge in modulating reactivity particularly the comparison between neutral BH₃ and its anionic form (BH₃⁻) remains largely unclear, despite its potential significance in tuning electronic structure and reaction barriers.

In addition to intrinsic mechanistic factors, environmental conditions such as solvation, pH, and temperature play critical roles in determining the kinetics of hydrogen evolution. In aqueous media, solvation can significantly alter the stability of intermediates and transition states,¹⁹ while the relative populations of BH₃ and BH₃⁻ are governed by acid–base equilibrium, leading to strong pH-dependent speciation. Moreover, for fast radical reactions in solution, diffusion can impose an upper limit on the observed rate constants, such that the overall kinetics reflect a complex interplay between intrinsic reactivity, chemical speciation, and mass transport. Despite its importance, a quantitative understanding of how these factors collectively influence hydrogen formation in boron hydride systems remains elusive.

Computational chemistry, particularly density functional theory (DFT), provides a robust and versatile framework for investigating reaction mechanisms and kinetics at the molecular level.²⁰ Modern hybrid meta-GGA functionals, such as M06-2X, have demonstrated reliable performance in describing thermochemistry, noncovalent interactions, and reaction barriers, making them well-suited for radical-mediated systems.^21,22 When combined with continuum solvation models and transition state theory (TST), including tunneling corrections, DFT enables the prediction of rate constants and their dependence on temperature and environment.¹⁹ Furthermore, coupling DFT with transition state theory (TST) and Eyring based kinetic formalisms allows for the prediction of rate constants and their dependence on temperature and environmental conditions.^23,24

Previous theoretical studies have primarily focused on hydrogen release from more complex boron containing systems, such as ammonia borane and related compounds, highlighting the importance of solvent effects, cooperative interactions, and multi-hydrogen transfer pathways.^4,5 In contrast, the simpler BH₃ system has received comparatively little attention, particularly in the context of radical induced hydrogen evolution. Moreover, a systematic comparison between neutral and anionic pathways, as well as a quantitative evaluation of pH-dependent kinetics, is still lacking in the current literature.

In this context, the present work aims to provide a comprehensive theoretical investigation of the reaction mechanisms and kinetics of hydrogen formation from the interactions of BH₃ and BH₃⁻ with hydroxyl radicals. Using density functional theory at the M06-2X/6-311++G(d,p) level, in combination with PCM solvation and kinetic modeling, we (i) characterize the potential energy surfaces of the BH₃ + HO˙ and BH₃⁻ + HO˙ reactions in both gas phase and aqueous solution, (ii) identify key intermediates and transition states associated with hydrogen evolution, with particular emphasis on double hydrogen atom transfer mechanisms, and (iii) evaluate the rate constants as a function of temperature and pH, including diffusion and tunneling effects, with direct comparison to available experimental data. By integrating detailed mechanistic analysis with quantitative kinetic modeling, this study provides new insights into the fundamental factors governing hydrogen evolution in boron hydride systems. The results not only clarify the role of charge, solvation, and environmental conditions in modulating reactivity but also offer valuable guidance for optimizing hydrogen generation processes, contributing to the development of efficient and sustainable energy technologies.

2. Computational methods

All electronic structure calculations were carried out using the Gaussian 09 (ref. 25) program package. Geometries of all stationary points, including reactants, intermediates, transition states, and products, were fully optimized, employing the M06-2X hybrid meta-GGA functional in combination with the 6-311++G(d,p) basis set.²⁶ This level of theory is well established for accurately describing thermochemical properties and noncovalent interactions in main-group systems, and is therefore appropriate for the BH₃ and BH₃⁻ + HO˙ reactions considered herein. The optimized Cartesian coordinates are provided in the SI (Table S1).

Harmonic vibrational frequency calculations were performed at the same level to verify the nature of each stationary point, with minima exhibiting no imaginary frequencies and transition states characterized by a single imaginary frequency. Zero-point energy (ZPE)²⁷ corrections were included in all reported energies. Open-shell species were treated using the unrestricted M06-2X formalism to ensure a proper description of spin states.²⁸ Transition states were located using the QST2/QST3 procedures within the STQN framework,²⁹ and their connectivity was confirmed by intrinsic reaction coordinate (IRC) calculations linking each transition state to the corresponding reactant and product minima (Fig. S1).^30,31 Solvent effects were incorporated using the conductor-like polarizable continuum model (PCM)³² to account for aqueous-phase conditions during both geometry optimizations and frequency analyses. To gain further insight into the electronic structure, topological analyses based on the Electron Localization Function (ELF)³³ and the Atoms-in-Molecules (AIM)³⁴ theory were performed using the Multiwfn 3.7 (ref. 35) program.

Rate constants were evaluated using conventional transition state theory (TST) at 298.15 K and 1 atm (ref. 36) according to the standard Eyring–Polanyi equation, where the activation Gibbs free energy (ΔG^≠) includes the necessary zero-point energy corrections:

where σ is the reaction symmetry number,^37,38 κ represents the tunneling correction (evaluated using the Wigner barrier),³⁹ k_B is the Boltzmann constant, ℏ is Planck constant, ΔG^≠ Gibbs free energy of activation. All rate constant calculations were performed using the Eyringpy package.⁴⁰

To accurately reproduce experimental macroscopic kinetic data, this basic gas phase formulation was subsequently corrected for both diffusion controlled limitations and solution pH effects.^40,41 In the liquid phase, bimolecular encounters are physically constrained by the rate of reactant diffusion through the solvent cage. Accordingly, the final apparent rate constant (k_app) for each specific reaction channel was determined using the Collins–Kimball equation:⁴²

where k_D is the steady-state diffusion-limited rate constant derived from the Smoluchowski expression k_D = 4πR_ABD_ABN_A, in which the mutual diffusion coefficient D_AB is modulated by solvent viscosity (η) via the Stokes–Einstein relation.

Furthermore, for reactants involving acid–base equilibria, the overall observed rate constant (k_over) was evaluated as a function of environmental pH to account for speciation shifts. This macroscopic rate coefficient was computed as the weighted sum of the kinetic contributions from all co-existing forms and their corresponding pH-dependent molar fractions (f_i):

3. Results and discussions

3.1 Validation of computational method and benchmark comparison

To validate the reliability of the M06-2X/6-311++G(d,p) computational level, the calculated energetic, structural, and vibrational properties of BH₃ and BH₃⁻ were directly compared with available experimental data (Table 1) as well as high-level theoretical results, including MP4/6-31+G(d,p) and CCSD(T)/CBS. In the gas phase, the optimized B–H bond length of BH₃ is 1.188 Å (see Fig. 1), which is in excellent agreement with experimental data (1.187 Å (ref. 43)), MP4 results (1.187–1.187 Å (ref. 44)), and high level CCSD(T)/CBS⁴⁵ calculations (1.187 Å). The deviation from the most accurate CCSD(T)/CBS value is less than 0.001 Å, clearly demonstrating the high structural accuracy of the chosen DFT method. For the BH₃⁻ anion, the calculated B–H bond length (1.205 Å) also agrees well with the reference value (1.207 Å), and correctly reproduces the expected bond elongation due to occupation of antibonding orbitals upon electron attachment.

Table 1 Computed energies, B–H bond lengths, and vibrational frequencies of BH₃ and BH₃⁻ at the M06-2X/6-311++G(d,p) level in the gas phase and PCM solvation. Literature data are provided in parentheses

Species	Energies (hartree)	B–H length (Å)	Mode (cm^−1)
BH3-G	−26.592182 (−26.508684 (ref. 44))	1.18755 (1.1869 (ref. 44)) (1.18739 (ref. 43)) (1.18674 (ref. 45))	ω1 = 2752.7 (2575.2;^45 2668 (ref. 44))
			ω2 = 1174.2 (1170.8;^45 1223 (ref. 44))
			ω3 = 2627.6 (2525.1;^45 2783 (ref. 44))
			ω4 = 1217.6 (1226.9;^45 1297 (ref. 44))
BH3-W	−26.592936	1.1875	ω1 = 2744.3
			ω2 = 1156.7
			ω3 = 2626.8
			ω4 = 1210.8
BH3^−-G	−26.588449 (−26.494303 (ref. 44))	1.20478 (1.2071 (ref. 44))	ω1 = 2485.7 (2525 (ref. 44))
			ω2 = 577.3 (563 (ref. 44))
			ω3 = 2561.9 (2606 (ref. 44))
			ω4 = 1140.1 (1244 (ref. 44))
BH3^−-W	−26.685188	1.20478	ω1 = 2553.2
			ω2 = 325.8
			ω3 = 2476.1
			ω4 = 1127.4

---

Fig. 1 Comparison of optimized geometries of BH₃ and BH₃⁻ in the gas phase and aqueous solution (PCM) at the M06-2X/6-311++G(d,p) level. Selected structural parameters are given in Å and degrees.

The vibrational frequencies further confirm the robustness of the computational approach. For BH₃, the symmetric stretching mode is calculated at 2752.7 cm⁻¹, compared to experimental and high level theoretical values of 2575.2 cm⁻¹ (CCSD(T)/CBS)⁴⁵ and 2668 cm⁻¹ (MP4).⁴⁴ Other modes, including the bending and asymmetric stretching vibrations (1174.2, 2627.6, and 1217.6 cm⁻¹), are also in good agreement with CCSD(T)/CBS results (1170.8, 2525.1, and 1226.9 cm⁻¹) and experimental data (1223 and 1297 cm⁻¹).⁴³ Importantly, the agreement with CCSD(T)/CBS data considered the benchmark for small molecules provides strong evidence for the reliability of the predicted vibrational properties.

For the BH₃⁻, the calculated vibrational modes (2485.7, 577.3, 2561.9, and 1140.1 cm⁻¹) are also consistent with the MP4 reference values (2525, 563, 2606, and 1244 cm⁻¹). The slight underestimation of high frequency modes and the more pronounced deviation in the low-frequency bending mode are expected due to anharmonic effects and the increased flexibility of the anionic structure. Nevertheless, the overall agreement confirms that the method captures both bonding strength and vibrational characteristics accurately. In terms of energetics, the calculated total energies are consistent with previously reported MP4 values, and although MP4/6-31+G(d,p) is a higher order post Hartree–Fock method in terms of formal electron correlation treatment, modern hybrid meta-GGA functionals such as M06-2X have been shown to provide comparable accuracy for thermochemical properties at significantly lower computational cost. The consistency between the M06-2X results and both the MP4 and CCSD(T)/CBS benchmarks observed here supports this conclusion.

The inclusion of solvent effects via the PCM model results in only minor changes in bond lengths and vibrational frequencies, indicating that the intrinsic structural features of BH₃ and BH₃⁻ are preserved in aqueous solution. However, a noticeable stabilization of the total energy is observed, particularly for the anionic species, due to electrostatic solvation effects. Overall, the excellent agreement between the present calculations and experimental, MP4, and especially CCSD(T)/CBS reference data confirms that the M06-2X/6-311++G(d,p) level of theory provides a reliable and balanced description of both neutral and anionic boron hydrides. This validation establishes a solid foundation for subsequent investigations of reaction mechanisms and kinetics in the BH₃ + HO˙ and BH₃⁻ + HO˙ systems.

3.2 Reaction mechanism of BH₃ + HO˙ in gas phase and aqueous solution

The potential energy surface (PES) for the reaction between neutral BH₃ and the hydroxyl radical HO˙ (Fig. 2) reveals a double hydrogen atom transfer (DHAT) mechanism leading to H₂ formation. Two distinct reaction channels, denoted as eqn (1) and (2), are identified as operative in both the gas phase and aqueous solution.

BH₃ + HO˙ → LM1 → TS1 → OHBH˙ + H₂↑ (1)

BH₃ + HO˙ → LM1 → TS2 → BH₂H˙ + H₂↑ (2)

Fig. 2 Potential energy surface and optimized geometries for the BH₃ + HO˙ → H₂ reaction at the M06-2X/6-311++G(d,p) level. Energies (kcal mol⁻¹) inside and outside parentheses refer to the aqueous phase (PCM) and gas phase, respectively. Selected bond lengths (Å) and bond angles (°) are shown for gas-phase (“G”) and aqueous (“W”) structures.

In the gas phase, the reaction is initiated by the formation of a pre-reactive complex (LM1-G), which is significantly stabilized at −48.9 kcal mol⁻¹ relative to the separated reactants (BH₃ + HO˙). This pronounced stabilization indicates a strong interaction between the electron deficient boron center and the hydroxyl radical, primarily arising from donor acceptor interactions between the lone pair of the oxygen atom and the vacant p orbital of boron.

As a consequence, two B–H bonds are elongated (1.281 Å) relative to the third one (1.189 Å), suggesting that hydrogen transfer preferentially occurs at these weakened B–H bonds. From LM1-G, the channel (1) proceeds via the first transition state (TS1-G), located at −36.7 kcal mol⁻¹, corresponding to the initial hydrogen transfer step. TS1-G is characterized by a single imaginary frequency, confirming its nature as a first order saddle point. Structurally, this transition state shows significant elongation of the B–H bonds (up to 1.512 and 1.711 Å), indicating the departure of two hydrogen atoms from the boron center, while the forming H–H bond is shortened to approximately 0.805 Å. The small H–B–H angle (28.1°) further supports the formation of an incipient H₂ molecule. The relatively low activation barrier from LM1-G to TS1-G (12 kcal mol⁻¹) suggests that this pathway is kinetically accessible, leading to the formation of the OHBH˙-G radical and H₂-G with an energy of −46.6 kcal mol⁻¹. An alternative pathway (2) proceeds via a second transition state (TS2-G, see Fig. S2), located at −0.8 kcal mol⁻¹ relative to the reactants. The imaginary frequency associated with TS2-G indicates that a DHAT mechanism is still operative; however, it differs from TS1-G in that the hydroxyl radical directly participates in the hydrogen transfer by contributing one hydrogen atom. In this structure, the O–H bond is elongated to 1.274 Å, while one B–H bond is significantly weakened (1.455 Å) and another is moderately elongated (1.347 Å). These two hydrogen atoms subsequently combine to form H₂. This distinct hydrogen transfer pattern leads to a different product channel, yielding BH₂O˙ + H₂. Overall, the reaction is highly exothermic, with the final products lying at −50.0 kcal mol⁻¹ relative to the reactants. When solvent effects are included via the PCM model, the overall reaction mechanism remains qualitatively similar to that in the gas phase, with both reaction channels (1) and (2) still operative. However, notable quantitative differences are observed. The pre-reactive complex (LM1-W) is slightly less stabilized in solution (−47.8 kcal mol⁻¹), reflecting the attenuation of specific intermolecular interactions due to solvation. The first transition state (TS1-W) is also marginally destabilized (−36.1 kcal mol⁻¹), although the relative barrier with respect to LM1-W remains nearly unchanged, indicating that the initial hydrogen transfer step is only weakly affected by the solvent environment. In contrast, the second transition state (TS2-W) exhibits a significant increase in energy, shifting from −0.8 kcal mol⁻¹ in the gas phase to +0.9 kcal mol⁻¹ in solution. This increase suggests that solvent effects impose an additional energetic penalty on the H₂ formation step, likely due to solvent reorganization and the reduced entropy associated with the formation of a more constrained transition state. The final products are also slightly less stabilized in solution (−48.0 kcal mol⁻¹) compared to the gas phase (−50.0 kcal mol⁻¹), indicating a decrease in the overall exothermicity of the reaction.

Structurally, the key features of intermediates and transition states are largely preserved in the PCM calculations, with only minor variations in bond lengths and angles, confirming that the intrinsic reaction pathway is not altered by solvation. Overall, the reaction proceeds via a DHAT mechanism, and solvent effects do not change the mechanistic pathway; instead, they slightly increase the activation barrier and reduce the thermodynamic driving force. These findings indicate that hydrogen formation from the BH₃ + HO˙ system is more favorable in the gas phase than in aqueous solution, primarily due to weakened transition state stabilization and increased entropy penalties in the solvent environment.

3.3 Reaction mechanism of BH₃⁻ + HO˙ in gas phase and aqueous solution

Following the mechanistic analysis of the neutral BH₃ + HO˙ system, it is essential to examine the corresponding reaction involving the anionic species BH₃⁻ in order to elucidate the role of charge on reactivity and hydrogen evolution. The presence of an extra electron significantly alters the electronic structure of the boron center, enhancing its nucleophilicity and modifying its interaction with the hydroxyl radical, thereby leading to distinct mechanistic features. In this context, the reaction of BH₃⁻ with HO˙ is investigated in both the gas phase and aqueous solution.

BH₃⁻ + HO˙ → LM2 → TS3 → BH₂O˙⁻ + H₂↑ (3)

BH₃⁻ + HO˙ → LM2 → TS3 → LM3 → TS3a → BH₂O˙⁻ + H₂↑ (4)

BH₃⁻ + HO˙ → LM2 → TS4 → OHBH˙⁻ + H₂↑ (5)

In the gas phase, two reaction pathways, eqn (3) and (5), are identified for the BH₃⁻ + HO˙ system, whereas in aqueous solution the reaction proceeds via eqn (4) and (5). As illustrated in Fig. 3, the corresponding potential energy surface (PES) exhibits a mechanistic profile markedly different from that of the neutral BH₃ system, primarily due to the increased electron density at the boron center, which strengthens its interaction with the hydroxyl radical.

Fig. 3 Potential energy surface and optimized geometries for the BH₃⁻ + HO˙ → H₂ reaction at the M06-2X/6-311++G(d,p) level. Energies (kcal mol⁻¹) in parentheses/without parentheses refer to the aqueous phase (PCM) and gas phase, respectively. Selected bond lengths (Å) and bond angles (°) are shown for gas-phase (“G”) and aqueous (“W”) structures.

The reaction is initiated by the formation of a strongly stabilized pre-reactive complex, LM2, with energies of −109.7 kcal mol⁻¹ in the gas phase (LM2-G) and −109.1 kcal mol⁻¹ in solution (LM2-W). This exceptionally large stabilization reflects a strong ion radical interaction, arising from electrostatic attraction and charge assisted hydrogen bonding. Structurally, LM2 (Fig. 4) is characterized by the formation of a B–O interaction (1.511 Å), which is longer than that in the neutral complex LM1, accompanied by perturbation of the B–H bonds (1.232 Å), indicating significant structural reorganization upon electron enrichment of BH₃⁻. The corresponding structure in solution remains largely unchanged, with only minor geometric variations. From LM2-G, the reaction proceeds via two competing pathways. The dominant pathway involves TS3-G, located at −72.3 kcal mol⁻¹ in the gas phase, corresponding to an activation barrier of approximately 37.4 kcal mol⁻¹ relative to LM2-G. The TS3-G is characterized by a double hydrogen atom transfer (DHAT) mechanism, in which each reactant contributes one hydrogen atom to form H₂-G. At this transition state, the two hydrogen atoms are separated by 2.033 Å, with O–H and H–B distances of 0.972 and 2.499 Å, respectively, while the B–O bond shortens to 1.372 Å. These structural features clearly indicate a concerted hydrogen transfer involving simultaneous hydride donation from boron and proton transfer from the hydroxyl radical.

Fig. 4 ELF maps and AIM analysis of LM1-G and LM2-G in the gas phase: (a) ELF map of LM1-G, (b) ELF map of LM2-G, (c) AIM analysis of LM1-G, and (d) AIM analysis of LM2-G.

In the gas phase, pathway (3) proceeds directly from TS3-G to form the products BH₂O˙⁻ + H₂ at −104.2 kcal mol⁻¹, indicating a highly exothermic process. In contrast, under aqueous PCM conditions, this pathway (4) becomes stepwise rather than concerted. At TS3-W (−78.8 kcal mol⁻¹), the first hydrogen transfer from BH₃⁻ occurs, as evidenced by the significant elongation of the B–H bond (3.774 Å), while the hydroxyl group remains coordinated to boron. This leads to the formation of an intermediate (LM3-W) at −88.8 kcal mol⁻¹. Subsequently, a second transition state (TS3a-W) is located at −90.3 kcal mol⁻¹, slightly lower in energy than LM3-W, indicating an almost barrierless rearrangement. In this step, the transferred hydrogen atom facilitates abstraction of a second hydrogen from the hydroxyl group, resulting in H₂ formation and yielding the final products BH₂O˙⁻ + H₂ at −103.4 kcal mol⁻¹. The overall process remains highly exothermic in solution. An alternative pathway (5) proceeds via the fourth transition state, the TS4-G (see Fig. S2) located at −39.0 kcal mol⁻¹ in the gas phase (TS4-W has −36.6 kcal mol⁻¹ in solution), corresponding to a significantly higher activation barrier (70 kcal mol⁻¹ relative to LM2). In this transition state, a DHAT mechanism is also observed; however, the hydroxyl group plays a largely spectator role, maintaining a B–OH interaction (1.41–1.43 Å), while two hydrogen atoms originating from BH₃⁻ combine to form H₂. This pathway leads to the formation of (OHBH˙⁻ + H₂)-G, with final energies of −42.4 kcal mol⁻¹ and −46.4 kcal mol⁻¹ (OHBH˙⁻ + H₂)-W. Compared to the TS3 pathway, this route is both kinetically less favorable and thermodynamically less exothermic, indicating that it represents a minor reaction channel.

The inclusion of solvent effects via the PCM model does not alter the overall mechanistic framework but introduces notable quantitative changes. The stability of LM2 is only slightly affected, indicating that solvation has a limited influence on the strong ion radical interaction. Importantly, TS3-W is further stabilized relative to the gas phase, effectively lowering the activation barrier and facilitating hydrogen transfer. This suggests that polar solvation stabilizes the charge-separated character of the transition state. In contrast, TS4-W is slightly destabilized, further increasing its already high barrier and reinforcing the dominance of the TS3 pathway. The relative stability of the final products is only marginally affected, and the ordering of reaction channels remains unchanged. Structurally, the optimized geometries of all stationary points show only minor differences between gas phase and solution, confirming that the intrinsic reaction mechanism is preserved upon solvation. Key features such as B–O bonding and the extent of B–H and O–H bond elongation at transition states remain consistent across environments.

Overall, the reaction of BH₃⁻ with HO˙ proceeds predominantly via a DHAT mechanism along the TS3 pathway, which is both kinetically and thermodynamically favored. The rate determining step is associated with TS3, and the overall reaction is strongly exothermic. Solvent effects enhance the feasibility of this dominant pathway by stabilizing the key transition state, in contrast to the neutral BH₃ system where solvation tends to increase the activation barrier. These results highlight the critical role of charge in modulating reactivity and provide important insights into the efficiency of hydrogen evolution in boron hydride systems.

3.4 Electronic structure analysis (ELF and AIM)

To further elucidate the electronic factors governing the reactivity differences between the neutral and anionic systems, a detailed topological analysis of the electron density was performed using the Electron Localization Function (ELF) and Atoms-in-Molecules (AIM) approaches. These complementary methods provide deeper insight into bonding characteristics, electron delocalization, and intermolecular interactions that cannot be fully captured by geometric and energetic analyses alone. In particular, they are highly effective in identifying electronic features associated with hydrogen transfer processes, such as bond polarization, electron redistribution, and the emergence of incipient H–H interactions.

The analysis was carried out in the gas phase, as ELF and AIM are fundamentally based on the topology of the electron density (Fig. 4), which reflects the intrinsic electronic structure of the system. As demonstrated in the previous sections, the inclusion of solvent effects via the PCM model results in only minor variations in structural parameters and does not alter the fundamental reaction pathways. Therefore, the key electronic features governing reactivity are expected to be preserved upon solvation. Moreover, the continuum nature of the PCM model does not explicitly account for solvent solute electronic interactions at the level required for meaningful topological analysis. Consequently, gas-phase calculations provide a reliable and physically meaningful representation of the bonding nature without introducing unnecessary complexity.

Based on this approach, ELF and AIM analyses were performed for the key pre-reactive complexes, LM1-G in the neutral system and LM2-G for anionic system. The ELF results reveal a clear distinction in electron distribution between the two systems. In LM1-G, the electron density is primarily localized around the boron center and along the B–H bonds, with limited delocalization toward the hydroxyl fragment, indicating predominantly polarized covalent bonding. In contrast, LM2-G exhibits a markedly higher degree of electron localization, particularly around the hydride ligands, where the ELF basins become more diffuse and extended. This behavior reflects the increased electron density associated with the extra electron, leading to partial delocalization toward hydrogen atoms involved in subsequent H₂ formation. Notably, electron density in LM2-G is redistributed toward regions between hydrogen atoms, suggesting a preorganization that facilitates early H–H bond formation. In addition, the ELF domain between boron and oxygen becomes more pronounced, consistent with a stronger B–O interaction.

The AIM analysis further supports these observations. In LM1-G, typical bond critical points (BCPs) are observed for the B–H and B–O interactions, with moderate electron density values characteristic of conventional bonding. In contrast, LM2-G shows enhanced BCP characteristics, particularly along the B–O interaction pathway, indicating a stronger and more stabilized interaction. At the same time, the electron density at the B–H, BCPs in LM2-G becomes reduced and more asymmetric, reflecting bond weakening and polarization, which are essential prerequisites for hydride transfer. The redistribution of electron density toward hydrogen atoms further supports the formation of a transient H–H interaction even at the pre-reactive stage.

Overall, the combined ELF and AIM analyses consistently demonstrate that the anionic system possesses a more reactive electronic structure. Enhanced electron donation from boron strengthens the B–O interaction while simultaneously weakening the B–H bonds, thereby promoting efficient hydrogen transfer and facilitating H₂ formation at an earlier stage of the reaction. These findings provide a clear electronic structure basis for the lower activation barriers and higher reactivity observed for the BH₃⁻ system compared to its neutral counterpart.

3.5 Kinetic analysis

Before discussing the detailed kinetic results, it is important to distinguish between gas-phase and solution-phase behavior, as the governing factors differ significantly in these environments. In the gas phase, the reaction kinetics are primarily controlled by intrinsic potential energy surfaces and pre-equilibrium effects, whereas in solution, additional factors such as solvation, pH-dependent speciation, and diffusion play a decisive role. Therefore, the kinetic analysis is presented sequentially: first for the gas phase, where intrinsic reactivity dominates, and then for the aqueous phase, where environmental effects modulate the overall reaction rate.

3.5.1 Kinetic analysis for hydrogen evolution reactions in the gas phase. Building upon the mechanistic and electronic structure insights obtained from the PES and ELF/AIM analyses, a quantitative kinetic investigation was carried out to evaluate the feasibility of hydrogen evolution under different conditions. In particular, rate constants were computed to assess the temperature dependence of the reactions in the gas phase, and to determine the extent to which quantum tunneling contributes to hydrogen transfer processes, given the light mass of hydrogen atoms involved in the DHAT mechanism. The temperature dependence of the rate constants for the reactions BH₃⁻ + HO˙ and BH₃ + HO˙ in the gas phase, including Wigner tunneling corrections, is presented in Fig. 5 and Table S2. Two sets of data are reported for each system: the conventional transition state theory rate constant, k(T), and the tunneling corrected rate constant, k_wig(T).

Fig. 5 Temperature dependence of rate constants for BH₃⁻ + HO˙ and BH₃ + HO˙ reactions in the gas phase (280–350 K), showing both conventional TST k(T) and Wigner-corrected k_wig(T) values.

For the anionic reaction (BH₃⁻ + HO˙), the rate constants are significantly higher than those of the neutral system over the entire temperature range (280–350 K). The values of k(T) decrease monotonically with increasing temperature, from approximately 4.7 × 10⁴³ (cm³ per molecule per s) at 280 (K) to 1.7 × 10³² (cm³ per molecule per s) at 350 K, indicating a negative temperature dependence. This behavior suggests that the reaction proceeds through a strongly stabilized pre-reactive complex (LM2-G), where the effective barrier is submerged relative to the reactants, leading to a rate that decreases as thermal energy increases. Upon inclusion of Wigner tunneling corrections, the corrected rate constants k_wig(T) are consistently higher than the corresponding k(T) values about 1.19 to 1.39 times, confirming the contribution of quantum tunneling to hydrogen transfer. However, the magnitude of the correction is relatively modest, indicating that although tunneling is operative, the reaction is not dominated by tunneling effects. The fitted expression for k_wig(T) shows a weak inverse temperature dependence (T^−1.7) combined with an exponential term, further supporting a barrierless or near-barrierless reaction profile with quantum mechanical contributions.

In contrast, the neutral system (BH₃ + HO˙) exhibits substantially lower rate constants, with k(T) values decreasing from approximately 2.9 × 10¹⁶ (cm³ per molecule per s) at 280 (K) to 4.2 × 10¹⁰ (cm³ per molecule per s) at 350 (K). Similar to the anionic case, a negative temperature dependence is observed, although less pronounced. This trend is consistent with a reaction pathway involving a pre-reactive complex (LM1-G), but with a higher effective barrier compared to the anionic system. The inclusion of Wigner tunneling corrections leads to a light increase in the rate constants, particularly at lower temperatures, reflecting the effect of tunneling in this system. The fitted expression for k_wig(T) characterized by a positive temperature exponent (T^0.11) and an exponential term, indicates a thermally activated process in which tunneling partially compensates for the activation barrier. A comparison between the two systems clearly demonstrates that the anionic reaction is kinetically more favorable by several orders of magnitude across the entire temperature range. This enhancement can be attributed to the increased electron density in BH₃⁻, which facilitates stronger interaction with the hydroxyl radical, weakens the B–H bonds, and lowers the effective barrier for hydrogen transfer. Consequently, the role of tunneling in the anionic system is less pronounced relative to its already high reactivity, whereas in the neutral system, tunneling plays a more significant role in enabling hydrogen transfer, especially at lower temperatures.

Overall, the results indicate that both reactions exhibit non-Arrhenius behavior with a negative temperature dependence, characteristic of complex mediated mechanisms. This behavior originates from the formation of a pre-reactive complex (LM1-G), which is enthalpically stabilized but entropically disfavored due to the association of reactants into a bound structure. As temperature increases, the growing entropic penalty destabilizes LM1-G relative to the separated reactants, leading to a decrease in the equilibrium constant K_eq. Although the intrinsic TST rate constant k_TST increases with temperature, the strong stabilization of LM1-G renders the pre-equilibrium contribution dominant, resulting in an overall decrease of the effective rate constant. The results also demonstrate that the presence of an excess electron in BH₃⁻ enhances intrinsic reactivity while reducing the relative contribution of tunneling compared to the neutral system.

3.5.2 Kinetic analysis in aqueous solution. In addition, the kinetics of the BH₃ (neutral and anionic) + HO˙ reactions were investigated in aqueous solution with explicit consideration of pH, temperature, and diffusion effects. Because the relative populations of BH₃ and BH₃⁻ are governed by acid–base equilibrium, variations in pH strongly influence the overall reaction rate. Moreover, in the liquid phase, diffusion imposes an additional kinetic constraint, particularly for fast radical reactions. As a result, the effective rate constants are not solely determined by intrinsic activation barriers but arise from the interplay between chemical speciation, thermal activation, and mass transport.

The kinetic parameters in the solution phase (Fig. 6), along with the detailed computational results, have been added to the SI as Table S3.

Fig. 6 Dependence of the rate constant (log k) on pH and temperature for the BH₃/BH₃⁻ + HO˙ reaction in aqueous solution, showing comparison between calculated results and experimental data for (a) pathway via TS1, and (b) pathway via TS3/TS3a.

The calculated rate constants for the BH₃ + HO˙ reaction in aqueous solution exhibit a pronounced dependence on pH, temperature, and diffusion, in good agreement with experimental observations.⁴⁶ As shown in Fig. 6, the log k increases steadily with pH over the range 1–9, followed by a plateau at higher pH values (pH ≥ 10). Specifically, log k rises from approximately 1.5–2.0 at low pH to about 9.5–10.0 under strongly alkaline conditions, indicating significant rate enhancement in basic media and a saturation regime at high pH. Notably, the calculated results closely reproduce experimental data, particularly at pH = 7 and T = 350 (K), where the predicted rate constant (2.1 × 10⁸ M⁻¹ s⁻¹) is in excellent agreement with the reported value (1.0 × 10⁸ M⁻¹ s⁻¹). This agreement validates the level of theory and the kinetic model, including the treatment of solvation, protonation equilibria, and diffusion.

Importantly, despite the intrinsic reactivity differences observed in the gas phase, the rate constants of the BH₃ + HO˙ and BH₃⁻ + HO˙ reactions become nearly comparable in aqueous solution. This convergence arises because the overall kinetics are governed by three coupled factors diffusion, pH-dependent speciation, and temperature. Under these conditions, the reactions approach the diffusion controlled limit, where mass transport and encounter frequency dominate over intrinsic activation barriers. Consequently, the advantage of the anionic pathway is partially diminished, leading to similar effective rate constants for both systems. Temperature further enhances the reaction rates across the entire pH range (280–350 K), with nearly parallel trends suggesting that the fundamental reaction mechanism remains unchanged. Overall, the combined influence of pH, temperature, and diffusion indicates that the reaction is favored under alkaline and high temperature conditions, while diffusion plays a decisive role in limiting and homogenizing the observed kinetics in solution.

4. Conclusions

In this study, the reaction mechanisms and kinetics of hydrogen formation from the interactions of BH₃ + HO˙ and BH₃⁻ + HO˙ were systematically investigated using density functional theory at the M06-2X/6-311++G(d,p) level. The calculated structural and vibrational properties show good agreement with available experimental and high level theoretical data, supporting the reliability of the computational approach.

PES analysis indicates that the BH₃ + HO˙ reaction proceeds via a DHAT mechanism, involving distinct transition states associated with sequential hydrogen transfer steps. The overall process is highly exothermic, yielding H₂ and oxygen-containing boron products. Inclusion of solvent effects via the PCM model does not alter the fundamental mechanism but slightly increases the activation barriers and reduces the reaction exothermicity.

Kinetic analysis reveals non-Arrhenius behavior in the gas phase, where the rate constants decrease with increasing temperature due to a dominant pre-equilibrium involving a strongly stabilized intermediate. In aqueous solution, the overall kinetics are governed by the combined effects of pH, temperature, and diffusion. While the anionic pathway is intrinsically more reactive, pH-dependent speciation and diffusion limitations lead to comparable effective rate constants for both systems.

Overall, this work provides key mechanistic and kinetic insights into boron hydride radical chemistry and highlights the roles of pre-equilibrium, solvation, and transport effects in controlling hydrogen generation efficiency.

Author contributions

Trinh Le Huyen: conceptualization, methodology, validation, data curation, formal analysis, visualization, writing – original draft, writing – review & editing. Pham Cam Nam: supervision, software, writing – review & editing.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6ra03503e.

Acknowledgements

This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 104.06-2025.12 (TLH).

References

1. G. Crabtree and M. Dresselhaus, Phys. Today, 2004, 57, 39–44.
2. L. Schlapbach and A. Züttel, Nature, 2001, 414, 353–358.
3. J. A. Turner, Science, 2004, 305, 972–974.
4. F. H. Stephens, V. Pons and R. T. Baker, Dalton Trans., 2007, 2613–2626.
5. A. Staubitz, A. P. Robertson and I. Manners, Chem. Rev., 2010, 110, 4079–4124.
6. Q.-L. Zhu and Q. Xu, Energy Environ. Sci., 2015, 8, 478–512.
7. T. Le Huyen, T. T. T. Huyen, C. C. Nam and P. C. Nam, RSC Adv., 2025, 15, 23760–23771.
8. A. Gutowska, L. Li, Y. Shin, C. M. Wang, X. S. Li, J. C. Linehan, R. S. Smith, B. D. Kay, B. Schmid and W. Shaw, Angew. Chem., 2005, 117, 3644–3648.
9. T. B. Marder, Angew. Chem., Int. Ed., 2007, 46, 8116–8118.
10. W. N. Lipscomb, Science, 1977, 196, 1047–1055.
11. H. C. Brown, Boranes in Organic Chemistry, Cornell University Press, 2019.
12. P. J. Dolan, J. H. Kindsvater and D. G. Peters, Inorg. Chem., 1976, 15, 2170–2173.
13. J. H. Morris, H. J. Gysling and D. Reed, Chem. Rev., 1985, 85, 51–76.
14. C. Sonntag, Free-radical-induced DNA damage and its repair: a chemical perspective, Springer Science & Business Media, 2006.
15. L. Wojnárovits, T. Tóth and E. Takács, Crit. Rev. Environ. Sci. Technol., 2018, 48, 575–613.
16. B. Halliwell, Toxicol. Ind. Health, 1993, 9, 1–21.
17. J. M. Mayer, Annu. Rev. Phys. Chem., 2004, 55, 363–390.
18. S. Hammes-Schiffer and A. A. Stuchebrukhov, Chem. Rev., 2010, 110, 6939–6960.
19. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3094.
20. C. J. Cramer and F. Bickelhaupt, Angew. Chem., Int. Ed. Engl., 2003, 42, 381.
21. Y. Zhao and D. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
22. Y. Zhao and D. G. Truhlar, Accounts Chem. Res., 2008, 41, 157–167.
23. H. Eyring, J. Chem. Phys., 1935, 3, 107–115.
24. D. G. Truhlar, B. C. Garrett and S. J. Klippenstein, J. Phys. Chem., 1996, 100, 12771–12800.
25. M. Frisch, G. Trucks, H. Schlegel, G. Scuseria, M. Robb, J. Cheeseman, G. Scalmani, V. Barone, G. Petersson and H. Nakatsuji, Gaussian 09, revision A.02, Gaussian Inc., Wallingford, 2009.
26. Y. Zhao, N. E. Schultz and D. G. Truhlar, J. Chem. Theor. Comput., 2006, 2, 364–382.
27. R. S. Grev, C. L. Janssen and H. F. Schaefer III, J. Chem. Phys., 1991, 95, 5128–5132.
28. A. Kumar, F. Neese and E. F. Valeev, J. Chem. Phys., 2020, 153, 094105.
29. C. Peng and H. Bernhard Schlegel, Isr. J. Chem., 1993, 33, 449–454.
30. S. Maeda, Y. Harabuchi, Y. Ono, T. Taketsugu and K. Morokuma, Int. J. Quantum Chem., 2015, 115, 258–269.
31. L. Deng and T. Ziegler, Int. J. Quantum Chem., 1994, 52, 731–765.
32. S. Miertuš, E. Scrocco and J. Tomasi, Chem. Phys., 1981, 55, 117–129.
33. A. Savin, R. Nesper, S. Wengert and T. F. Fässler, Angew Chem. Int. Ed. Engl., 1997, 36, 1808–1832.
34. R. F. Bader, Accounts Chem. Res., 1985, 18, 9–15.
35. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
36. B. C. Garrett and D. G. Truhlar, Encyclopedia of Computational Chemistry, 1998, 5, 3094.
37. A. Fernández-Ramos, B. A. Ellingson, R. Meana-Pañeda, J. M. Marques and D. G. Truhlar, Theor. Chem. Acc., 2007, 118, 813–826.
38. E. L. I. Pollak and P. Pechukas, J. Am. Chem. Soc., 1978, 100(10), 2984–2991.
39. E. Wigner, Z. Phys. Chem. B, 1932, 19, 203–216.
40. E. Dzib, J. L. Cabellos, F. Ortíz-Chi, S. Pan, A. Galano and G. Merino, Int. J. Quantum Chem., 2019, 119, e25686.
41. A. Galano and J. R. Alvarez-Idaboy, J. Comput. Chem., 2013, 34, 2430–2445.
42. F. C. Collins and G. E. Kimball, J. Colloid Sci., 1949, 4, 425–437.
43. K. Kawaguchi, J. Chem. Phys., 1992, 96, 3411–3415.
44. C. T. Wickham-Jones, S. Moran and G. B. Ellison, J. Chem. Phys., 1989, 90, 795–806.
45. M. S. Schuurman, W. D. Allen and H. F. Schaefer III, J. Comput. Chem., 2005, 26, 1106–1112.
46. D. Borchardt, J. G. Choi, K. Suzuki and S. Bauer, J. Chem. Phys., 1988, 88, 6282–6289.

---