Effects of heteroatom doping on hydrogen uptake in tungsten oxide

Abstract

Redox-active transition metal oxides (TMOs) that can undergo proton-insertion coupled electron transfer (PICET) are promising candidates for catalyzing molecular conversion reactions which require the transfer of hydrogen atoms (or the thermochemical equivalent, H⁺, e⁻). Herein, we studied the effects of isovalent (Mo⁶⁺) and aliovalent (V⁵⁺ and Nb⁵⁺) heteroatom doping on the electrochemical PICET behavior of monoclinic tungsten oxide (WO₃). Cyclic voltammetry in aqueous acidic electrolytes shows that the addition of redox-active heteroatoms (Mo⁶⁺ and V⁵⁺) leads to systematic shifts in redox couple half-wave potentials (E_1/2), broadening, and an overall decrease in the current response. Conversely, the non-redox active heteroatom (Nb⁵⁺) only reduces the current response with no observable peak-current broadening. This broadening is attributed to changes in the proton binding affinities of oxygen in different chemical environments, i.e., bridging different pairs of redox-active transition metal cations. We determined the hydrogen bond dissociation free energy (H BDFE) values to elucidate the thermodynamic effect of heteroatom substitution. Density functional theory calculations reveal a differentiation in the hydrogen binding and oxygen vacancy formation energies between heteroatom doped structures. The PICET-induced structural phase transitions of the pristine and doped samples were further probed with operando electrochemical X-ray diffraction (EC-XRD) and with ex situ chemical reduction. The broadening of the potential-dependent current response with increased heteroatom doping manifests in the operando EC-XRD results as prolonged structural regions where multiple hydrogen bronze phases exist and the appearance of cubic bronze phases at lower degrees of reduction compared to pristine WO₃.

---

1 Introduction

Early transition metal oxides (TMOs), such as those based on titanium^1,2 and tungsten,^3,4 can take up hydrogen atom equivalents (i.e., proton/electron pairs) via proton-insertion coupled electron transfer (PICET) reactions from acidic environments.⁵ This property can be utilized for a variety of electrochemical and catalytic applications, including in aqueous energy storage,^5–10 water electrolysis,^4,11 and heterogeneous hydrogen transfer catalysis.^12,13 All of these require the transfer of protons and electrons across a phase boundary (i.e., electrochemical interface, TMO/liquid, or TMO/gas). In TMOs, the thermochemical properties associated with PICET can be tuned across a variety of length scales, from the chemical composition and crystal structure to defect concentration, particle morphology, and particle size.^2,14,15 This high degree of tunability enables various strategies for controlling proton-coupled electron transfer (PCET) reaction kinetics and selectivity.

At the atomic level, reversible PICET involves the formation and homolytic cleavage of covalently bonded hydrogen atoms within the host structure. In TMOs, the inserted protons are localized to the oxide ligands, forming hydroxyl (M–OH) groups while the electrons populate conduction band states associated with the transition metal cations.⁵ The fundamental steps of PICET involve (1) mass transport of protons in the liquid electrolyte; (2) transfer of an electron to the solid host; (3) proton-coupled electron transfer (PCET) at the interface; and (4) mass transport of the proton in the solid state (Fig. 1a). The interfacial PCET reaction may proceed sequentially or in a concerted fashion and is represented by a reaction square scheme, as shown in Fig. 1b, using the example of PICET into WO₃. The energy change associated with forming the O–H bond, which requires both PT and ET steps, is represented as the diagonal reaction within the square scheme. Also termed the bond dissociation free energy (BDFE), this parameter can be quantified from the pH-dependence of the PICET half-wave potential (E_1/2).¹⁶ By understanding the factors which govern such PICET reactions on TMOs, we can better design materials for energy and catalytic applications.

Fig. 1 Fundamental steps of PICET (a) and reaction square scheme for PICET into WO₃ (b). Visual representations of heteroatom doping into WO₃ with a qualitative comparison of the ionic radii for the heteroatom dopants explored here (c) and projected electronic structures (d) of the WO₃ host upon heteroatom doping.

Matson et al. showed that heteroatom doping of polyoxotungstate clusters, which are molecular-scale analogs of extended tungsten oxides, allows for control over the PCET thermochemistry.^17–20 For example, the Lindqvist-type cluster, (TBA)₂[W₆O₁₉], undergoes two ET-reactions, with no measurable dependence on the electrolyte pK_a, suggesting that the reduction is not coupled to proton transfer (PT). Vanadium doping of this cluster to form (TBA)₃[VW₅O₁₉] led to a change in the redox behavior, wherein the two ET-reactions of (TBA)₂[W₆O₁₉] gave way to a single redox-couple in (TBA)₃[VW₅O₁₉] associated with PCET at the V-site.²⁰ The Keggin-type cluster, (TBA)₃[PW₁₂O₄₀], undergoes four PCET reactions and the BDFEs could be increased by >20 kcal mol⁻¹ by substituting a single W-site for V.¹⁹ The complete substitution of W⁶⁺ with Mo⁶⁺, forming the isostructural Keggin-type polyoxomolybdate cluster, (TBA)₃[PMo₁₂O₄₀], was also found to increase the BDFE by 20 kcal mol⁻¹.¹⁸ These molecular studies emphasize the impact of heteroatom doping and the importance of transition metal identity, dopant location, and orbital hybridization on PCET, serving as an inspiration to investigate the role of heteroatom dopants on PICET in WO₃.

In extended TMOs, the incorporation of heteroatom dopants can tune the electronic structure by changing the distribution of the electron density near the Fermi level. This effect can be achieved through the introduction of either isovalent or aliovalent transition metal dopants. The ability (or lack thereof) for the dopant ion orbital states to hybridize with the oxide p-orbitals perturbs the relative energies of conduction band states, changing the energetic requirements for electron transfer and delocalization.²¹ This compositional tuning not only modifies the electronic energy landscape of the TMO but can also impact PICET. Both isovalent and aliovalent doping can alter metal–oxygen bond lengths and angles, which in turn can modify the stable proton insertion sites. Miu et al. reported that PICET into WO₃ could be accurately predicted by analyzing the acid–base properties of the oxide host, further demonstrating a pair-wise correlation between hydrogen binding, Bader charges, and electronic density of states band centers, which were used as descriptors for the oxide Lewis acidity/basicity.²² Unique to aliovalent doping is the possibility for secondary compositional changes, such as the formation of vacancies or interstitial ion inclusions needed to compensate for the change in oxidation state (Fig. 1c), which could further alter the proton insertion behavior. Considering this, we propose two extreme hypothetical scenarios to describe the influence of heteroatom doping on PICET in TMOs: (1) limited hybridization between the dopant and host results in discreet PICET reactions such as what is seen in polyoxometalate clusters or (2) extensive hybridization induces systemic changes in the PICET reactions associated with the host transition metal (Fig. 1d). Coupled with these atomic-scale changes is the possibility of heteroatom doping to alter particle size or morphology, which could dramatically affect the quantity or availability of reactive surface sites to participate in PCET.

Here, we investigate how isovalent (Mo⁶⁺) and aliovalent (V⁵⁺ and Nb⁵⁺) heteroatom doping affects the thermodynamics and kinetics of proton-insertion coupled electron transfer (PICET) in monoclinic (γ) WO₃. Structural changes upon doping were characterized by XRD and Raman spectroscopy, while SEM and BET analyses assessed the microstructure and surface area. The distribution of the heteroatoms within the sample particles was assessed with depth-dependent XPS measurements. Electrochemical PICET behavior was studied using cyclic voltammetry at varying scan rates and pH in acidic aqueous electrolytes. Operando XRD and ex situ chemical reduction were used to probe the PICET-induced structural dynamics. Redox-active dopants (Mo⁶⁺ and V⁵⁺) shifted the onset of structural transitions, whereas the non-redox-active dopant (Nb⁵⁺) suppressed PICET and increased the hydrogen evolution onset potential. DFT calculations were further used to calculate the hydrogen binding energies (HBEs), oxygen vacancy formation energies (OVFEs), and net charge of the metal atoms within the doped systems. The HBEs were found to be dependent on the identity of the transition metals which coordinate to the proton accepting oxygen. These findings demonstrate that heteroatom doping provides a tunable handle for modulating hydrogen bonding and PCET thermochemistry in transition metal oxides.

2 Experimental methods

2.1 Chemicals

All chemicals were used as received. Na₂WO₄·2H₂O (99+%), concentrated HCl (12.1 M, Certified ACS Plus), Na₂MoO₄·2H₂O (99+%), concentrated HNO₃ (15.6 M, Certified ACS), NaVO₃ (96%), NbCl₅ (99.9%), TEMPO (98+%), toluene (Certified ACS), acetylene black (99.9+%), acetone (ACS Certified), ethanol (anhydrous, histological grade), concentrated H₂SO₄ (18 M, Certified ACS Plus), and NaOH (Certified ACS) were purchased from ThermoFisher. H₂PtCl₆ (8% by weight in H₂O), H₃PO₂ (50% by weight in H₂O), lithium (99.9% trace metals basis), 1 M LiPF₆ in ethylene carbonate (EC) : dimethyl carbonate (DMC) (50 : 50 battery grade), and n-methyl pyrrolidone (NMP, anhydrous, 99.5%) were purchased from Millipore Sigma. Poly(vinylidene fluoride) (PVDF; Kynar HSV 900) was purchased from Arkema.

2.2 Materials synthesis

Monoclinic tungsten oxide dihydrate (WO₃·2H₂O) was synthesized using the acid precipitation method developed by Freedman.²³ 50 mL of a 1 M Na₂WO₄·2H₂O aqueous solution was added dropwise to 450 mL of a 3 M HCl solution that was stirred at 300 rpm at room temperature. This mixture was allowed to stir for three hours, and the precipitate was collected via vacuum filtration and washed with DI water until the pH of the rinse solution reached ∼6. Monoclinic tungsten trioxide (γ-WO₃) was formed by thermally dehydrating the as-prepared WO₃·2H₂O at 300 °C for 5 hours in air in a box furnace.

The heteroatom doped samples (Mo⁶⁺, V⁵⁺, and Nb⁵⁺:WO₃) were prepared by modified acid-precipitation methods. For the Mo-containing compounds, first a fixed amount of Na₂WO₄·2H₂O was added to a 20 mL scintillation vial with a magnetic stir bar. Next, a stoichiometric amount of Na₂MoO₄·2H₂O was added to the vial along with 2.5 mL of DI H₂O. The precursor masses and solvent volumes used for the Mo-doped sample syntheses are given in Table S1. The sample was allowed to stir until all the solids were completely dissolved. This solution was then added dropwise to a separate vial containing 10 mL of a 6 M HNO₃ solution in DI H₂O. Upon addition, a solid precipitate formed. The sample was allowed to continue stirring for three hours to promote crystallinity and ensure reaction completeness. The powder product was then separated via vacuum filtration, washed multiple times with water (pH of rinse solution ∼6) and once with acetone, and then dried at 60 °C in air. The powder was then ground for 10 minutes with an agate mortar and pestle, then loaded into an alumina crucible and annealed in air at 300 °C for 5 hours in a box furnace.

For the V-containing compounds, it was found that excess precursor reagent, NaVO₃, was required to promote V-incorporation into the WO₃ host phase. The NaVO₃ and Na₂WO₄·2H₂O precursors were dissolved in 5 mL of DI H₂O in a 20 mL scintillation vial under constant heating and stirring at 80 °C in an oil bath. Once fully dissolved, the solution was allowed to cool to room temperature before being pipetted into a separate vial containing 10 mL of a 6 M HNO₃ solution in DI water. Upon addition, a solid precipitate formed. The sample was allowed to continue stirring for three hours to promote crystallinity and ensure reaction completeness. The powders were separated and annealed following the same procedure as that for the Mo-containing samples. Precursor masses and solvent volumes used for the V-doped sample syntheses are also given in Table S1.

For the Nb-containing compounds, NbCl₅ was weighed separately and added into a 20 mL scintillation vial inside a nitrogen-filled glove box (MBraun, O₂ < 0.5 ppm, H₂O < 0.5 ppm). The vial was capped and removed from the glovebox and the powder was subsequently dissolved in 1 mL of ethanol. Separately, Na₂WO₄·2H₂O was dissolved in 2 mL of DI H₂O. Once both solids were fully dissolved in their respective solvents, the solutions were pipetted simultaneously into a separate vial containing 10 mL of a 6 M HCl solution in DI H₂O. HCl was used instead of HNO₃ for safety considerations associated with the reaction of HNO₃ with ethanol. Upon addition, a solid precipitate formed. The sample was allowed to continue stirring for three hours to promote crystallinity and ensure reaction completeness. The powders were separated and annealed following the same procedure as that for the Mo-containing samples. A table of precursor masses and solvent volumes used for the Nb-doped sample syntheses are also given in Table S1.

2.3 Physical characterization

Powder X-ray diffraction (PXRD) data were collected on a PANalytical X'Pert Pro powder diffractometer (45 kV, 40 mA, sealed Cu X-ray tube, Kα₁, Kα₂ λ = 1.5406 Å, 1.5444 Å) equipped with an X'Celerator position sensitive detector in Bragg–Brentano geometry. Pawley refinements of the XRD data were performed using Topas Academic (V7).²⁴ Crystal structure images were generated using Vesta 3.²⁵ For direct comparison with synchrotron diffraction measurements all angular ranges were converted to the scattering vector, Q, in reciprocal space (Å⁻¹) using the equation Q = (4π/λ) × sin (2θ/2). X-ray fluorescence (XRF) spectroscopy was performed on a Rigaku SuperMini200 sequential benchtop WDXRF spectrometer. For XRF measurements, approximately 1 g of powdered sample was enclosed within a Mylar film pouch before loading into the sample holder. The XRF results were corrected for elemental results associated with the Mylar film. Raman spectroscopy was performed using a WiTEC alpha 300R confocal Raman spectrometer with a laser wavelength of 532 nm and a 100× objective lens. Scanning electron microscopy (SEM) was performed on a Hitachi SU8700 equipped with an Everhart–Thornley detector. Samples were measured with an accelerating voltage of 2 keV with no applied stage bias. Brunauer–Emmett–Teller (BET) N₂ sorption isotherms were obtained on a TriStar Micrometrics setup at 77 K. All samples were degassed at 373 K under a constant flow of N₂ overnight.

X-ray photoelectron spectroscopy (XPS) measurements were performed at the Ambient Pressure X-ray Photoelectron Spectroscopy (APXPS) beamlines 9.3.1 and 9.3.2, at the Advanced Light Source (ALS). The Mo-doped sample was measured at Beamline 9.3.1, while the V-doped sample was measured at Beamline 9.3.2. Beamline 9.3.1 is a tender X-ray facility (photon energy range: 2.0–6.0 keV) equipped with a bending magnet and a Si (111) double-crystal monochromator.²⁶ Measurements were conducted at photon energies of 3400, 4000, and 5000 eV in normal emission geometry, with a pass energy of 100 eV using an R4000 HiPP-2 Scienta analyzer. The incident beam spot size was ∼350 µm in diameter, and the sample manipulator provided four-axis control (x, y, z, and θ). Beamline 9.3.2 is a soft X-ray facility covering photon energies from 200 to 900 eV.²⁷ Measurements were performed at 570, 610, 730, 810, and 890 eV, with a pass energy of 100 eV for all energies except for 570 eV, where 50 eV was used. The X-ray spot size was 1.0 mm in diameter, incident at 15° relative to the sample surface.

All XPS measurements at Beamline 9.3.2 were carried out at room temperature under ultra-high vacuum (UHV, ∼10⁻⁹ Torr) conditions, whereas measurements at Beamline 9.3.1 were conducted under high-vacuum conditions (HV, ∼10⁻² Torr). The binding energy (BE) scale was calibrated against the Au 4f_7/2 core-level peak (BE = 84.0 eV). Samples were prepared by the drop-cast method. Briefly, 30 mg of the as-prepared powder was dispersed in 6 mL of a DI water/ethanol mixture (1 : 2 volume ratio; 2 mL DI water and 4 mL ethanol) and ultrasonicated for 30 min to achieve uniform dispersion. The resulting suspension was then drop-cast onto an Au-coated Si wafer and dried at room temperature.

2.4 Chemical reduction and hydrogen transfer reactions

Hydrogen tungsten bronze materials (H_zMo⁶⁺:WO₃ and H_zV⁵⁺:WO₃) were prepared via a solution-based reduction process previously developed by our group.¹³ Briefly, 0.2375 g of the oxide materials were added into a 20 mL scintillation vial along with 312 µL of H₂PtCl₆ (5% mass loading with respect to TMO), 5 mL of DI water, a magnetic stir bar, and 300 µL of H₃PO₂. The reaction mixtures were capped and heated between 95 °C and 100 °C in an oil bath with constant magnetic stirring at 500 rpm for 2 hours, after which they were allowed to cool to room temperature naturally with continued stirring. During the heating process, the yellow powders gradually turned blue indicating reduction of the transition metal ions. After the reaction, the samples were transferred to 50 mL centrifuge tubes and centrifuged at 4500 rpm for 5 minutes. After this, the supernatant solution was discarded and the solid was resuspended in 30 mL of H₂O before being centrifuged again. This process was repeated two more times with the final wash using acetone. The samples were then dried at room temperature under vacuum overnight. After drying, samples were stored in a N₂-filled glovebox (MBraun, <0.5 ppm O₂ and <0.5 ppm H₂O) to prevent oxidation.

For hydrogenation reactions involving TEMPO, approximately 50 mg of the hydrogen bronze powders were added into a 20 mL scintillation vial along with a magnetic stir bar and 7 mL of a 0.05 M solution of TEMPO in toluene. The contents were capped and allowed to magnetically stir overnight at room temperature. After stirring, the solutions were transferred to 10 mL centrifuge tubes and centrifuged at 4500 rpm for 5 minutes to separate the solid from the TEMPO solution.

UV-Vis-NIR spectroscopy data of the supernatant TEMPO solutions were collected on an Ocean-Insights HDX-XR spectrometer equipped with an HL-2000-HP-FHSA light source, a SQUARE-ONE cuvette holder, and QP600-1-VIS-NIR fiber optic cables. For measurements, 2.5 mL of the TEMPO supernatant solutions were added to a quartz cuvette. Scans were collected with 100 ms integration times and averaged over 15 scans with a boxcar width of 5. A calibration curve was constructed using solutions with known concentrations of TEMPO (0.10, 0.07, 0.05, 0.03, and 0.01 M) in toluene. The change in absorbance intensity at 470 nm was fit to a linear regression trend following Beer's law, absorbance = εbc, where ε is the molar absorption coefficient of TEMPO (slope of the trend line), b is the pathlength of the cuvette used (1 cm) and c is the concentration of the solution. The absorbance spectra collected for the TEMPO reactions with the bronze phases could then be used to determine the extent of conversion of optically active TEMPO to inactive TEMPOH.

2.5 Electrode preparation

For electrochemical measurements, slurry-cast composite electrodes on carbon paper (Fuel Cell Earth) substrates were prepared from the as-synthesized powder samples. The carbon paper substrates were cleaned by O₂ plasma treatment for five minutes before the composite slurries were deposited on the substrate. First, 10 mg of PVDF was dissolved in 500 µL of NMP. Next, 80 mg of the powder sample was ground together with 10 mg of acetylene black for five minutes with acetone used as a grinding medium. The solvent was evaporated, and the homogeneous powder mixture was added to the PVDF/NMP solution. The resulting mixture was then stirred for six minutes using a Thinky Mixer (Thinky U.S.A.) at 2000 rpm to create a viscous ink. 7 µL of the slurry ink were then pipetted onto half of a 2 cm² carbon paper electrode giving an approximate mass loading of 1 mg of active material per 1 cm² of electrode area. The electrodes were then heated overnight at 120 °C under vacuum to remove the NMP solvent. For electrochemical analysis in non-aqueous electrolytes, the slurry-cast composite electrodes were prepared in the same manner but were deposited onto stainless steel mesh substrates (SS304 woven, Fisher Scientific).

2.6 Electrochemical characterization

Cyclic voltammetry (CV) experiments were performed on a BioLogic MPG2 potentiostat. The three electrode cells used for electrochemical analysis consisted of the as-prepared carbon paper working electrodes suspended vertically in a three-neck flask with a Pt wire counter electrode (Sigma-Aldrich), a Ag/AgCl reference electrode (4 M KCl, Pine Research Instrumentation), and 20 mL of a 0.5 M H₂SO₄ (pH = 0.30) solution in DI water serving as the electrolyte. The pH of the electrolyte solution was adjusted by the dropwise addition of a 1 M NaOH solution in DI H₂O. A pH probe (Mettler Toledo, FiveEasy) was used to measure the pH of the electrolyte solution. The electrodes were allowed to soak in the electrolyte within the cell for 30 minutes prior to cycling to ensure proper wetting of the carbon paper substrate. CV experiments were conducted from −0.4 V to +0.5 V vs. the Ag/AgCl reference (unless otherwise noted) with scan rates of 100, 50, 10, 5, and 1 mV s⁻¹. The specific capacity and subsequent proton content were determined by integrating the current response (corrected for active material mass loading) with respect to time. The potential scale was converted to the reversible hydrogen electrode (RHE) scale using the equation

E(vs. RHE) = E(Ag/AgCl) + 0.059 × (pH) + 0.197

Non-aqueous CV experiments were performed in an argon filled glove box (MBraun, O₂ < 0.5 ppm, H₂O < 0.5 ppm) using a BioLogic VMP3 potentiostat. Measurements were conducted in three electrode cells consisting of the as-prepared stainless steel working electrodes suspended vertically in a three-neck flask with Li metal serving as both the counter and reference electrodes and 20 mL of a 1 M LiPF₆ in EC/DMC serving as the electrolyte. CV experiments were conducted from 4.0 V to 1.5 V vs. the Li-metal reference with a 1 mV s⁻¹ scan rate. The specific capacity and subsequent Li-content were determined by integrating the current response (corrected for active materials mass loading) with respect to time. The potential range was converted to the RHE scale using the following equation:

E(vs. RHE) = E(vs. Li/Li⁺) − 3.040

2.7 Operando structural characterization

Operando electrochemical X-ray powder diffraction (EC-XRD) measurements were conducted at the Stanford Synchrotron Radiation Lightsource (SSRL) using BL 2-1 in Bragg–Brentano geometry with an incident energy of 17 keV and a Pilatus 100 K area detector was used to collect the area diffraction patterns. The 2D diffraction patterns were compiled and integrated using a Python script developed at SSRL beamline 2-1. Diffraction patterns were collected from 10° to 20° 2θ with a fixed incident angle (θ) of 4°. For direct comparison with laboratory diffraction measurements all angular ranges were converted to the scattering vector, Q, in reciprocal space (Å⁻¹) using the equation Q = (4π/λ) × sin (2θ/2). All operando measurements were performed in an in situ cell following previously described protocols.^4,28 The cell was composed of a 5.0 mm thick PEEK electrode holder with a 1/32″ thick fluorosilicone rubber gasket (McMaster-Carr) for preventing electrolyte leaks and a polyimide (Kapton) film to provide a low background window suitable for X-ray transmission. The in situ cell was assembled and then placed in a helium (He) chamber. The full assembly was mounted onto the diffractometer. For operando EC-XRD measurements, slurry-cast composite electrodes were deposited onto titanium mesh (gauze, 100 mesh, ThermoFisher Scientific) substrates with dimensions of 40 mm × 5 mm. For these electrodes, the mass loading was increased to approximately 5 mg cm⁻² to improve the diffraction statistics. For measurements, cyclic voltammetry was performed in the in situ cell at a scan rate of 0.5 mV s⁻¹ for two cycles from −0.38 V to +0.4 V (vs. Ag/AgCl) using a BioLogic SP150 potentiostat in a three-electrode configuration with carbon paper serving as the counter electrode, a leakless Ag/AgCl reference electrode (Pine), and a 0.5 M H₂SO₄ solution in DI H₂O serving as the electrolyte.

2.8 Density functional theory (DFT) calculations

DFT calculations were performed starting with the pristine, monoclinic WO₃ unit cell (64 atoms) with a 2 × 2 × 2 supercell (256 atoms), which served as the basis for both the undoped and doped systems. Doping was achieved by substituting two tungsten (W) atoms with two identical dopant atoms (e.g. Mo, V, or Nb), resulting in doped WO₃ supercells (3% dopant substitution). For each doped configuration, two spatial arrangements of the dopants were considered: a “near” configuration, where the dopants were separated by a single oxygen atom, and a “far” configuration, where the dopants were further spatially separated. Hydrogen intercalation was then modeled by inserting a single hydrogen atom into the supercell. In the “near” doped structures, the hydrogen atom was bound to the bridging oxygen atom between one dopant and one tungsten atom (this is the site with the lowest energy), while in the “far” doped structures, the hydrogen atom was bound to the oxygen positioned between two tungsten atoms neighboring one dopant (this is also the site with the lowest energy). To further explore defect effects, an oxygen vacancy was introduced into the non-hydrogenated structures by removing the aforementioned oxygen atoms. Subsequently, an atomic hydrogen was added to these structures positioned as described in S4. All structures were created using the Atomic Simulation Environment (ASE) library in Python and the Avogadro software.^29,30 All generated structures were geometrically optimized while keeping the lattice parameters fixed (at the lattice of the optimized, pristine WO₃).

The geometry optimizations of the structures were carried out using CP2K,³¹ which uses a mixed Gaussian and plane-wave approach: the wavefunctions are expanded in a Gaussian basis set, while the electronic density is represented on a plane-wave grid.³² The Perdew–Burke–Ernzerhof (PBE)³³ exchange–correlation functional was used to describe electron exchange and correlation effects.³⁴ Long-range dispersion interactions were accounted for with Grimme's D3 correction.^35,36 Triple-ζ valence polarized (TZVP) basis functions [IB5] were applied to oxygen and hydrogen atoms, while double-ζ valence polarized (DZVP) basis functions [IB6] were used for tungsten, niobium, vanadium, and molybdenum.³⁷ Core electrons were treated with Goedecker–Teter–Hütter (GTH) pseudopotentials.^38,39 Geometry relaxations were performed using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm,⁴⁰ with a maximum force criterion of 0.01 eV Å⁻¹, a self-consistent field (SCF) convergence threshold of 10⁻⁷ Hartree and a kinetic energy cutoff of 550 Ry.

Following structural optimization, Bader charge analysis was performed to calculate the net charges of the atoms of interest.⁴¹ Additionally, hydrogen binding energies (HBEs) were computed using the following expression:

HBE = E_S,H − E_S − 1/2E_H2,

where E_S,H is the total energy of the hydrogen-intercalated structure, E_S is the total energy of the corresponding non-intercalated structure, and E_H2 is the energy of an isolated H₂ molecule. Oxygen vacancy formation energies (OVFEs) were also calculated using the following equation, which accounts for the formation of water as a reference:

OVFE = E_S,VAC − E_S + E_H2O − E_H2,

where E_S,VAC is the total energy of the structure containing an oxygen vacancy and E_H2O is the energy of an isolated H₂O molecule.

3 Results and discussion

3.1 Physical characterization of heteroatom doped tungsten oxides

XRF was used to assess the incorporation of the heteroatom dopants into WO₃ from the modified acid precipitation syntheses (Table S1). The three series of heteroatom doped WO₃ compositions (Mo⁶⁺, V⁵⁺, and Nb⁵⁺:WO₃) were all found to be isostructural to monoclinic (γ)-WO₃ based on the XRD and Raman results (Fig. S1). There were more noticeable structural changes in the V⁵⁺:WO₃ series due to the smaller ionic radius of V⁵⁺ compared to W⁶⁺, Mo⁶⁺, and Nb⁵⁺. Representative XRD patterns and Raman spectra for compositions with 6% dopant substitution for all three series are shown in Fig. 2a and b. SEM (Fig. S2) and BET (Fig. S3) measurements indicate similar plate-like morphologies for all compositions with the heteroatom doped samples having smaller overall particle sizes and larger surface areas compared to the pristine WO₃. Depth-dependent XPS measurements (Fig. S4 and S5) were used to probe the distribution of the redox-active heteroatoms (Mo⁶⁺ and V⁵⁺) within the samples. The higher the XPS photon energy, the further the probing depth within the sample. These measurements show that Mo⁶⁺ tends to aggregate near the particle surface while V⁵⁺ tends to aggregate more in the sub-surface and bulk of the particle. The XPS measurements of the V⁵⁺ sample also show the appearance of oxygen vacancies, consistent with the decrease in the average transition metal oxidation state (Fig. S6). A more in-depth discussion on the physical characterization of these materials is provided in Discussion S1.

Fig. 2 Structural characterization – PXRD (a) and Raman spectroscopy (b) – and electrochemical characterization (c) for each heteroatom dopant at a 6% concentration compared to pristine WO₃. CV results for the Mo-series at all dopant concentrations (d); the inset figure shows the shift of the primary oxidation feature (O2) with the Mo⁶⁺ dopant concentration.

3.2 Electrochemical characterization of heteroatom doped tungsten oxides

We studied the effects of heteroatom doping on PICET in WO₃ with cyclic voltammetry (CV) in an acidic aqueous electrolyte (0.5 M H₂SO₄, pH = 0.30; Fig. 2c and d). We previously reported the electrochemical PICET reactions of WO₃.⁴ Briefly, the two sets of redox-couples observed in the CV are associated with PICET-induced structural transitions from monoclinic WO₃ to a tetragonal (O1/R1) and cubic (O2/R2) hydrogen bronze phase, H_zWO₃, respectively. These PICET reactions occur prior to the onset of the hydrogen evolution reaction (HER), a competing PCET reaction. For the heteroatom dopants, Mo⁶⁺ and V⁵⁺ can both undergo PICET reactions in their binary oxide phases (MoO₃ and V₂O₅) while Nb⁵⁺ (Nb₂O₅) does not undergo PICET prior to the onset of hydrogen evolution.⁵ For the two redox-active dopants, three main changes are observed as a function of dopant amount: (1) the potentials of the redox couples shift, (2) the peak currents decrease, and (3) the redox peaks become increasingly broader (Fig. 2c, d and S7). The Nb-doped series also shows decreased peak currents but there are only marginal, albeit negative, shifts in peak potential and markedly less peak broadening (Fig. 2c and S7). Unique to the Nb-doped series was an increase in HER activity at more positive potentials compared to WO₃ (Fig. S7b), which made it impossible to resolve the second redox-couple (O2/R2). Because of this, the cathodic turnover potential was limited to −0.15 V (vs. Ag/AgCl) for slow scan rate sweeps (1 mV s⁻¹ and 5 mV s⁻¹). The increased HER activity of a related system, Nb-doped W₁₈O₄₉, was achieved previously.⁴²

For the Mo-doped series, we observed a positive shift in peak potentials (Fig. 2d) suggesting that Mo-doping leads to a change in PICET thermodynamics (i.e., the BDFE for hydrogen). A positive shift indicates that more energy is required to homolytically cleave the electrochemically formed O–H bond.¹⁶ For the most doped sample (18% Mo⁶⁺:WO₃), we observed a +32 mV (+0.74 kcal mol⁻¹) shift in redox peak potential compared to WO₃. The V-doped series showed a gradual shift towards more negative potentials with 1% V⁵⁺:WO₃ having a small positive shift in the O2/R2 redox-couple and no shift for the O1/R1 couple compared to WO₃ (Fig. S7a). For the most doped sample (6% V⁵⁺:WO₃), we observed a −20 mV (−0.46 kcal mol⁻¹) shift in redox peak potential compared to WO₃. A decrease in peak potential would indicate a lower BDFE or less energy required to cleave the electrochemically formed O–H bond. The potential shift observed for the Nb-doped series was significantly smaller than that for the two redox-active series, showing only a −4 mV (−0.09 kcal mol⁻¹) shift at 16% Nb-doping (Fig. S7b). For all three dopant series, we observe a decrease in the number of electrochemically inserted protons/electrons compared to WO₃ (Fig. S7c). More discussion on the scan rate dependent measurements and kinetic analysis can be found in Discussion S2, Fig S8 and S9.

The broadening of the peak-current response is only observed for redox-active dopants (Mo⁶⁺ and V⁵⁺) and increases with their amount. For PICET with WO₃, inserted protons are localized at the bridging oxygens, forming a hydroxyl (M–OH) group. The electrochemical potential associated with PICET is thus influenced by the coordination environment of the oxygens. In doped WO₃, assuming a random distribution of dopants, we would expect three distinct oxygen environments: (1) bridging between two tungsten sites (W–O–W), (2) bridging between a tungsten site and a dopant site (D, i.e., W–O–D), and (3) bridging between two dopant sites (D–O–D). We expect that these different environments should lead to unique binding sites for inserted protons, although the extent of this difference may not be substantial enough to lead to isolated redox associated with the dopant, as observed in polyoxotungstate molecular clusters.^17,19

3.3 PICET thermochemistry of heteroatom doped tungsten oxides

We determined the hydrogen bond dissociation free-energies (H BDFEs) and pH-dependent potential shifts associated with the current peaks of the PICET reactions from CV experiments in pH adjusted electrolytes (0.5 M H₂SO₄ + NaOH; Fig. S10 and S11). For a 1 H⁺ per 1 e⁻ reaction, the pH dependence of the half-wave potential should follow a Nernstian trend with a shift of 59 mV per pH unit.^5,16 The BDFE(O–H) is determined from the intercept of the Nernstian trend using eqn (1).

BDFE(O–H) = 23.06 E⁰_1/2 + 52.8 kcal mol⁻¹ (1)

where E⁰_1/2 is the half-wave potential (vs. RHE) at pH = 0. It was not possible to determine the BDFE(O–H) for all heteroatom doped samples due to the current broadening or, in the case of Nb-doping, increased HER activity. The results for WO₃ and selected heteroatom doped samples are shown in Fig. 3 for the O1/R1 and O2/R2 redox couples, where the inset figure shows the pH dependence of the half wave potentials. In all samples, the pH dependence was >59 mV per pH unit. Such “super-Nernstian” behavior in TMOs has been noted before and attributed to > 1 H⁺ per e⁻ transfer or increasing metal–oxygen bond strength with pH.¹⁴ The doped samples showed slightly higher pH shifts than WO₃. Multi-electron redox is possible for both Mo⁶⁺ and V⁵⁺ transition metal cations. However, the decrease in capacities achieved with doping suggests that this isn't operational under the given conditions. For the V⁵⁺- and Nb⁵⁺-doped samples, the presence of oxygen vacancies could facilitate the adsorption of water molecules or hydronium ions at the particle surface. For both O1/R1 and O2/R2 couples, heteroatom doping shows a fine increase of the BDFE(O–H) relative to WO₃. Of this subset of materials studied, the largest increase occurs with 6% Mo⁶⁺ doping, with 0.37 kcal mol⁻¹ and 0.51 kcal mol⁻¹ higher than those for WO₃ for O1/R1 and O2/R2, respectively. This is similar, though milder, to the modulation of BDFE(O–H) in polyoxotungstate clusters (20 kcal mol⁻¹). The presence of the single heteroatom V dopant in [VW₅O₁₉]³⁻ could “switch on” PCET reactivity [BDFE(O–H) = 64 kcal mol⁻¹] relative to the homometallic tungsten cluster by increasing the basicity of the bridging oxo sites.¹⁹ We hypothesize that the increased BDFE(O–H) upon heteroatom doping of WO₃ are also due to the increased basicity of bridging oxygen sites albeit to a lesser extent.

Fig. 3 Bond dissociation free energies [BDFE(O–H)] determined for the two sets of redox peaks observed in the CV experiments (O1/R1 and O2/R2) with respect to the half-wave potential (E_1/2) of the redox couples. The BDFEs for the O2/R2 couple of the Nb-doped samples could not be determined due to the increased HER activity. The inset figure shows the pH dependent potential shifts for all compounds measured with respect to the expected Nernstian shift value shown with a dashed line.

3.4 Operando electrochemical characterization of PICET-induced structural transitions

The electrochemically induced structural transitions of the heteroatom doped samples were probed with operando electrochemical X-ray powder diffraction (EC-XRD) measurements. The PICET induced structural transitions for WO₃ were discussed in depth previously and our results are consistent with these previous reports (Fig. 4a). The two redox features observed are associated with a transition from monoclinic (P2₁/n) to a tetragonal bronze phase (O1/R1, P4/nmm) and a cubic bronze phase (O2/R2, Im3̄). The pseudocubic (a_p) lattice parameter was used to evaluate the structural changes during reduction. The a_p parameter was used previously to establish lattice parameter trends in perovskite materials where the transitions are associated with changes in the octahedral tilting pattern.^4,13,43 The equations used to calculate the a_p-lattice parameter for the different structural phases are given in the SI (eqn (S1)–(S3)). The a_p-lattice parameter increased upon PICET into WO₃. Towards the apex of the redox transitions (which shows as current peaks in the CVs), WO₃ has a biphasic mixture consisting of the two hydrogen bronze phases (Fig. 4b). This biphasic region exists only for a couple of scans at each redox event.

Fig. 4 Surface contour plots showing the diffraction pattern evolution over the course of the operando EC-XRD measurements for WO₃ (a), 18% Mo⁶⁺:WO₃ (c), and 1% V⁵⁺:WO₃ (e). Trend in pseudocubic, a_p, lattice parameter with respect to the proton content for WO₃ (b), 18% Mo⁶⁺:WO₃ (d), and 1% V⁵⁺:WO₃ (f). Full circles indicate the reduction sweep, and empty circles indicate the oxidation sweep for the first cycle.

For 18% Mo⁶⁺:WO₃ (Fig. 4c and d) which has one of the broadest current responses, operando EC-XRD shows more hysteresis between the structural transitions. Upon reduction, there is a prolonged biphasic region spanning 0.15 ≤ H_z ≤ 0.35, suggesting a more gradual transition than what was observed for WO₃. Upon oxidation, the tetragonal phase is only observed very briefly within the biphasic region (0.15 ≤ H_z ≤ 0.20) and quickly devolves into the monoclinic structure. This can be visualized as the smearing of the diffraction peaks upon oxidation and from the trend in the pseudocubic axis. The cubic phase appears at significantly lower degrees of reduction (H_z = 0.15) compared to the undoped sample (H_z = 0.40).

For 1% V⁵⁺:WO₃ (Fig. 4e and f) where the redox couples are still well defined, the structural transition sequence is more similar to the step-wise process observed in WO₃ as opposed to the gradual transition observed in 18% Mo⁶⁺:WO₃. However, it is again observed that the cubic phase appears at lower degrees of reduction for this sample (H_z = 0.20). We expect that the hydrogen binding affinities of the different oxygen chemical environments are enabling the structural transition to the cubic bronze phase in an effort to optimize the hydrogen bonding configuration and transport pathways within the structure. This effect would be further influenced by the amount of heteroatom dopants resulting in a greater distribution of less favorable W–O–D and D–O–D bonding environments. Similar operando measurements were performed for 6%V⁵⁺:WO₃ and the results are discussed in Discussion S3 and Fig. S12. Previously, we pointed out that the broadened current responses observed for the more doped samples made it difficult to determine the BDFE values for the distinct PICET reactions. The operando EC-XRD measurements presented here further rationalize this difficulty by showing the increased coexistence of the two bronze phases in samples containing higher dopant amounts.

3.5 Ex situ characterization of heteroatom doped hydrogen bronzes

We used chemical reduction to further assess the structures of the heteroatom doped hydrogen bronzes (H_zMo⁶⁺:WO₃ and H_zV⁵⁺:WO₃; Fig. 5). We previously reported a solution-based method utilizing H₂PtCl₆ and H₃PO₂ to prepare bulk powders of Pt-decorated (Pt@) H_zWO₃ with the equilibrium structure and proton/electron content that would exist for an electrode polarized at 0 V vs. RHE.¹³ All reduced samples display a diffraction pattern consistent with the cubic bronze phase. We previously discussed that the undoped material (H_zWO₃) should have a tetragonal structure with P4mm space group symmetry. The appearance of the cubic phase for this sample is a consequence of increased preferred orientation effects hiding the tetragonal axis. This stems from the smaller particle sizes used here compared to those used in our prior report.¹³ To verify the cubic structure of the heteroatom doped samples, we conducted identical chemical reduction reactions on Mo-doped and pristine WO₃ samples that were sintered at 500 °C for 8 hours in air to increase the average particle size (Fig. S13). For the undoped sample, the pattern displays additional peak splitting consistent with the previously reported tetragonal cell, suggesting that preferred orientation leads to the appearance of the cubic phase for this sample. For the Mo-doped samples, the diffraction patterns do not change, and the cubic structure is maintained in the larger particle size samples (Fig. S13).

Fig. 5 PXRD patterns for the chemically reduced tungsten oxides for the Mo-doped series (a) and the V-doped series (b).

We further evaluated the extent of chemical reduction and hydrogenation ability of the heteroatom doped samples using a chemical oxidant, 2,2,6,6-tetramethyl-1-piperidinyloxy (TEMPO). We previously used these reactions with TEMPO to demonstrate the ability of Pt@H_zWO₃ to perform hydrogen atom transfer (HAT) reactions and to determine the proton content of the bronze phase.¹³ The conversion of TEMPO to TEMPOH can be monitored via UV-Vis absorbance spectroscopy due to the change in optical response of the organic molecule upon HAT. The molar absorptivity of TEMPO in toluene was determined to be 10.22 M⁻¹ cm⁻¹ through linear fitting of the concentration dependent absorbance intensity curve (Fig. S14). The bronze reactions were allowed to stir overnight in a 0.05 M TEMPO in toluene solution after which the supernatants were collected fo r UV-Vis measurements. After the reaction with the bronzes, the absorbance profiles for all of the TEMPO supernatant solutions show a decrease in intensity, verifying that HAT had occurred (Fig. S15). All bronze materials were found to transfer between 0.18 and 0.28 H⁺/e⁻ per mol of TMO (Fig. S15c). The pristine sample took approximately 4 hours to return to its oxidized state, WO₃, as determined by the visual color change from blue back to yellow. Both series of doped compounds maintained some degree of blue coloration after the overnight reaction. This may suggest that only a portion of the inserted protons are able to react with TEMPO. Given that the undoped sample was fully oxidized within 4 hours, we may expect that the most accessible protons come from protonated oxygen sites which bridge two W-centers and that protonated W–O–D and D–O–D sites are less accessible. The reduction extent of the heteroatom doped bronzes determined from the reactions with TEMPO (Fig. S15c and d) is consistent with the previous operando XRD results showing that the cubic bronze phase exists at lower degrees of protonation/reduction (Fig. 4).

3.6 Density functional theory (DFT) calculations

To further understand the role of heteroatom dopants in the PICET insertion site energetics from an atomistic level perspective, we performed DFT calculations on bulk (periodic in 3 dimensions) structures. From the experimental results (see Discussion S1) we can be confident that the proton insertion and oxygen vacancy formation (for aliovalent dopants) occurred within the bulk rather than merely at the surface of the materials, which dictates the use of periodic bulk oxide structures for the DFT calculations. First, we performed Bader charge analysis to calculate the net charges of atoms of interest in the sequence D–O–M–O/X–W–O–W, where D represents a dopant (or tungsten only for the pristine structure), M represents either a dopant or a tungsten atom (in near and far conformations respectively), and X denotes an oxygen vacancy. Regarding the stability of the far vs. the near configurations, the calculations showed a negligible difference in total energy, with the near configurations being less than 0.3 eV more exothermic. We further computed the hydrogen binding energies (HBEs) for all structures, both with and without dopants and oxygen vacancies. In the vacancy-containing systems, we also calculated the HBEs for hydrogen atoms located between the two metal atoms (i.e., at the vacancy site as a metal hydride). Rather than calculating hydrogen interactions at different degrees of intercalation,²² we focused on single hydrogen atom binding energy to obtain clear trends on the binding behavior of different sites created in WO₃ upon doping and to avoid combinatorial complexity arising from different hydrogen binding configurations in the cell. Finally, we determined the oxygen vacancy formation energies (OVFEs) for each structure. As summarized in Table S2, the net charge of the oxygen (which binds hydrogen after PICET) in the unprotonated structures is practically unaffected by heteroatom doping. The only exception to this is the near-configuration of the V-doped structure, where the first oxygen's net charge increases from −1.02e (in the pristine structure) to −0.81e. Nevertheless, we cannot observe any correlation between the net charge of oxygen and HBE. In general, no obvious correlation was found between the charge distribution on the atoms of interest and HBE.

Regarding the HBEs, we notice that they are mostly thermoneutral (values clustered near zero) for the vacancy-containing structures (see Fig. S16). In contrast, vacancy-free structures (see Fig. S16) display very different behavior: the isovalent-doped (Mo-doped) and pristine systems exhibit mildly exothermic HBEs (i.e., close to thermoneutrality), whereas the aliovalent-doped (Nb- and V-doped) systems show highly exothermic values. For hydrogen bound at the vacancy site (see Fig. S16), the trend is reversed; the isovalent-doped systems exhibit very exothermic HBEs, while the aliovalent-doped ones are endothermic. These three scenarios are illustrated in Fig. 6, where HBE is plotted as a function of OVFE.

Fig. 6 Hydrogen binding energy (HBE) as a function of oxygen-vacancy formation energy (OVFE) for different local hydrogen environments. (a) HBE vs. OVFE for three adsorption configurations: hydrogen bound to the middle oxygen in the D–O–M–O–W–O–W chain (no oxygen vacancy) (blue). Hydrogen bound to the left oxygen in the D–O–M–X–W–O–W chain (with an oxygen vacancy) (red). Hydrogen occupying the oxygen-vacancy site between two metals in the D–O–M–X–W–O–W chain (green). Linear fits are shown to illustrate trends. (b) Differences in HBE for experimentally observed dopant environments relative to the pristine structure. For doped systems, in which two vacancy-containing configurations are possible, the values shown represent the average HBE relative to the pristine case. (c) Schematic representations of the atomic configurations corresponding to the three cases in (a) and (b) for both near and far configurations. Cyan: dopant metal (V, Nb, or Mo), red: O, pink: W, and white: H. The bottom right legend depicts the three cases in all figures.

For vacancy-free structures, HBEs show a strong linear correlation with OVFEs. A similar trend is observed for hydrogen positioned between the two metals in the vacancy-containing structures (i.e. as a hydride). In contrast, when hydrogen is bound to an oxygen adjacent to a vacancy, the coefficient of determination (R²) is lower, because the HBE values are practically thermoneutral (slope close to 0). Another interesting result observed is the separation of the doped structures with respect to their OVFEs. All aliovalent doped structures have significantly negative (exothermic) OVFEs whereas the isovalent and pristine ones exhibit much positive OVFEs, regardless of the position of the two dopants. The DFT results suggest that the oxidation state of the dopants affects HBE the most, since the differentiation is between aliovalent and isovalent doped structures. Net charges of all atoms in the sequence D–O–M–O/X–W–O–W do not seem to play any significant role in the HBEs.

Previously, we hypothesized that the broadening of the peak-current response in the CV measurements was attributed to the distribution of oxygen chemical environments due to the introduction of the heteroatoms. The DFT calculations validate this hypothesis by showing changes in the HBEs based on the identity of the two metals which the oxygen atom bridges. For the isovalent-doped samples, both far and near configurations show more exothermic hydrogen binding compared to the pristine structure (see Fig. 6b). Experimentally, we observe a positive shift in BDFE(O–H) which agrees with the more favorable hydrogen binding on the Mo-doped oxide. For aliovalent doped samples (V and Nb), the culmination of different available sites (i.e. on an oxygen site or at the vacancy site) leads to an overall decrease in HBEs (weaker binding) for the different configurations, consistent with the shifts to more negative potentials observed experimentally. Between the Mo- and V-doped calculations, the relative change in HBEs indicates that V-doping would have a more drastic impact on proton binding site distribution. This matches our experimental results showing similar degrees of peak-current broadening for 6%V⁵⁺:WO₃ and 18%Mo⁶⁺:WO₃. For the Nb-doped series, it was computationally found that the oxygen near the vacancy is the most favorable hydrogen insertion site. The neighboring vacancy sites stabilize hydrogen as hydrides with binding energies that are less exothermic than the corresponding V-sites. As a result, the Nb-doped structure is the only oxide that stabilizes vacancies and exhibits more balanced energetics of hydrogen interaction with both bridging oxygens and neighboring vacancies (less energetic variations). This may potentially explain the increased HER activity observed experimentally, involving vacancy-sites where protons couple with hydrides on the surface, since there is no evidence of the HER occurring in the bulk.⁴ These differences in hydrogen binding site energies may also have an impact on the solid-state proton transport, restricting diffusion to the lower energy pathways. Overall, the computational results match exceptionally well with the experimental results, highlighting the importance of dopant configuration within the crystal structures.

While PICET is sensitive to the local oxygen environment, cation insertion-coupled electron transfer is less localized. Intercalated alkali cations such as Li⁺ reside in interstitial void spaces governed by the overall crystal structure and are largely agnostic to the specific chemical coordination environments of the oxide ligands. Therefore, cycling the heteroatom doped WO₃ samples in a non-aqueous alkali cation-containing electrolyte should not result in the peak current broadening observed in the acidic aqueous electrolyte. To test this hypothesis, we performed CV of WO₃ and the 6% doped samples in a Li⁺-containing electrolyte (Fig. 7). Our results support this hypothesis by showing only marginal changes in current magnitude with no significant peak broadening during electrochemical Li⁺ insertion. The 6% V-doped sample is the only composition that shows a minor shift of the redox couples to more negative potentials, which can be attributed to the previously discussed increased structural strain (i.e., more different interstitial void-site volume). These experiments verify that the peak broadening observed in the aqueous electrolyte is attributed to the increased distribution of oxygen environments.

Fig. 7 Comparison of CVs at 1 mV s⁻¹ collected in 0.5 M H₂SO₄ (top) and 1 M LiPF₆ in EC/DMC (bottom) for WO₃ and the 6% heteroatom doped compositions.

4 Conclusions

We report the effect of isovalent (Mo⁶⁺) and aliovalent (V⁵⁺ and Nb⁵⁺) doping in monoclinic (γ) WO₃ on the thermodynamics of electrochemical PICET. Mo- and Nb-doped samples show no significant change in the crystal structure, while V-doping leads to broadening of the diffraction peaks consistent with increased structural strain due to the change in ionic radii. Electrochemical cycling of samples containing redox-active dopants (Mo⁶⁺ and V⁵⁺) in 0.5 M H₂SO₄ shows systematic shifts in the half-wave potentials for the two distinct redox couples (O1/R1 and O2/R2) and a broadening of the potential-dependent current response that increases with the dopant concentration. The incorporation of a non-redox-active dopant (Nb⁵⁺) considerably lowers the current response associated with PICET while increasing HER activity, with no substantial shift in redox-peak potential. The BDFE(O–H) values and pH-potential shifts for the doped and undoped samples were determined from pH-dependent CV. The sequence of PICET-induced structural transitions was probed with operando EC-XRD conducted in an acidic electrolyte. The broadening of the current response in the CV experiments corresponds with changes in the degree of reduction at which the hydrogen bronze phases appear and the duration of the bronze structural regimes. The presence of the heteroatom dopants changes the chemical distribution of bridging oxygen sites, finely altering the kinetics and thermodynamics of electrochemical PICET. DFT calculations support this hypothesis by showing changes in hydrogen binding energies relative to WO₃ based on the identity of metal atoms with which the oxygen/vacancy site bridges, showing more of an impact in V-doped systems than in Mo-doped systems. DFT calculations further revealed a linear correlation between the hydrogen binding energies and oxygen vacancy formation energies for both doped and undoped structures. This correlation was observed in cases without oxygen vacancies and when an inserted proton was positioned between two metal atoms in the presence of an oxygen vacancy. These findings suggest that the oxygen vacancy formation energy may serve as a reliable descriptor for hydrogen binding energy. The combination of experimental and computational results presented here highlights the importance of location and distribution of dopants/defects in the thermodynamics of hydrogen atom transfer reactions involving redox-active metal oxides.

Author contributions

N. P. H.: conceptualization, investigation, formal analysis, methodology, writing, review & editing; N. E. P.: investigation, formal analysis, methodology, writing, review & editing; J. R. P.: investigation, formal analysis, methodology, writing, review & editing; G. M.: funding acquisition, supervision, resources, writing, review & editing; E. J. C.: funding acquisition, supervision, resources, writing, review & editing; V. A.: funding acquisition, resources, conceptualization, supervision, writing, review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data from this work are openly available at https://doi.org/10.5281/zenodo.17402112.

Supplementary information (SI): amounts of reagents used for the synthesis of the heteroatom doped tungsten oxides (Table S1). Physical characterization discussion (Discussion S1). XRD and Raman spectra of all heteroatom doped samples (Fig. S1). SEM images of the heteroatom doped and undoped WO₃ samples (Fig. S2). BET sorption measurements and determined surface areas of the heteroatom doped and undoped WO₃ samples (Fig. S3). Synchrotron XPS characterization results for Mo⁶⁺ (Fig. S4) and V⁵⁺ (Fig. S5 and S6). CV results at 1 mV s⁻¹ for the aliovalent doped series in 0.5 M H₂SO₄ (Fig. S7). Discussion on scan rate dependence (Discussion S2). Scan rate dependent CV experiments for all compositions (Fig. S8). B-value analysis of the scan-rate dependent CV measurements (Fig. S9). Variable pH CV measurements and analysis of undoped and Mo/V-doped compositions (Fig. S10). Variable pH CV measurements and analysis of the Nb-doped compositions (Fig. S11). Equations used to calculate the pseudocubic lattice parameters (eqn (S1)–(S3)). Discussion on 6% V⁵⁺:WO₃ operando EC-XRD measurements (Discussion S3). Surface plots and pseudocubic lattice parameter trends for 6% V⁵⁺:WO₃ from operando EC-XRD measurements in 0.5 M H₂SO₄ (Fig. S12). PXRD patterns for the chemically reduced Mo-doped samples with larger particle sizes (Fig. S13). Concentration dependent absorption measurements and calibration curve for TEMPO in toluene (Fig. S14). TEMPO chemical oxidation and UV-Vis results (Fig. S15). DFT result discussion (Discussion S4). DFT calculations for net charges, HBEs, and OVFEs (Table S2). Visualization of the atoms of interest for the 3 cases of proton binding (Fig. S16). See DOI: https://doi.org/10.1039/d5sc08564k.

Acknowledgements

We thank Dr Molleigh Preefer, Dr Nicholas Strange, and Dr Johanna Nelson Weker at the SSRL for assistance during operando EC-XRD measurements. We also thank Vincenzo Musico and Seongbak Moon at NC State for SEM and BET measurements, respectively. We thank Prof Ellen Matson at the University of Rochester, Prof James McKone at the University of Pittsburgh, and Dr Saeed Saeed for thought-provoking discussion and guidance. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under award number DE-SC0023465. E. J. C. was partially supported by the Condensed Phase and Interfacial Molecular Science Program (CPIMS), in the Chemical Sciences, Geosciences, and Biosciences Division of the Office of Basic Energy Sciences of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Characterization was performed in part at the Analytical Instrumentation Facility (AIF) at North Carolina State University, which is supported by the State of North Carolina and the National Science Foundation (award number ECCS-2025064). The AIF is a member of the North Carolina Research Triangle Nanotechnology Network (RTNN), a site in the National Nanotechnology Coordinated Infrastructure (NNCI). Use of the Stanford Synchrotron Radiation Lightsource, SLAC National Accelerator Laboratory, was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Contract No. DE-AC02-76SF00515. This research used resources of the Advanced Light Source, which is a DOE Office of Science User Facility under contract no. DE-AC02-05CH11231. This work was supported by computational time granted from the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility – ARIS – under project ID pr017025. We acknowledge, also, the EuroHPC Joint Undertaking for awarding this project access to the EuroHPC supercomputer LUMI, hosted by CSC (Finland) and the LUMI consortium under project 465001836.

References

1. Y. S. Kim, S. Kriegel, K. D. Harris, C. Costentin, B. Limoges and V. Balland, Evidencing Fast, Massive, and Reversible H⁺ Insertion in Nanostructured TiO₂ Electrodes at Neutral pH. Where Do Protons Come From?, J. Phys. Chem. C, 2017, 121(19), 10325–10335,  DOI:10.1021/acs.jpcc.7b02395.
2. N. Makivić, J. Y. Cho, K. D. Harris, J. M. Tarascon, B. Limoges and V. Balland, Evidence of Bulk Proton Insertion in Nanostructured Anatase and Amorphous TiO₂ Electrodes, Chem. Mater., 2021, 33(9), 3436–3448,  DOI:10.1021/acs.chemmater.1c00840.
3. J. B. Mitchell, W. C. Lo, A. Genc, J. Lebeau and V. Augustyn, Transition from Battery to Pseudocapacitor Behavior via Structural Water in Tungsten Oxide, Chem. Mater., 2017, 29(9), 3928–3937,  DOI:10.1021/acs.chemmater.6b05485.
4. M. A. Spencer, N. P. Holzapfel, K. E. You, G. Mpourmpakis and V. Augustyn, Participation of Electrochemically Inserted Protons in the Hydrogen Evolution Reaction on Tungsten Oxides, Chem. Sci., 2024, 15(14), 5385–5402,  10.1039/d4sc00102h.
5. N. P. Holzapfel, V. Augustyn and V. Balland, Fundamentals of Proton-Insertion Coupled Electron Transfer (PICET) in Metal Oxides for Aqueous Batteries, ACS Energy Lett., 2025,(10), 1143–1164,  DOI:10.1021/acsenergylett.4c03076.
6. H. Jiang, W. Shin, L. Ma, J. J. Hong, Z. Wei, Y. Liu, S. Zhang, X. Wu, Y. Xu, Q. Guo, M. A. Subramanian, W. F. Stickle, T. Wu, J. Lu and X. Ji, A High-Rate Aqueous Proton Battery Delivering Power Below −78 °C via an Unfrozen Phosphoric Acid, Adv. Energy Mater., 2020, 10(28), 1–8,  DOI:10.1002/aenm.202000968.
7. W. Xu, K. Zhao, X. Liao, C. Sun, K. He, Y. Yuan, W. Ren, J. Li, T. Li, C. Yang, H. Cheng, Q. Sun, I. Manke, X. Lu and J. Lu, Proton Storage in Metallic H_1.75MoO₃ Nanobelts through the Grotthuss Mechanism, J. Am. Chem. Soc., 2022, 144(38), 17407–17415,  DOI:10.1021/jacs.2c03844.
8. Y. Xu, X. Wu and X. Ji, The Renaissance of Proton Batteries, Small Struct., 2021, 2(5), 1–16,  DOI:10.1002/sstr.202000113.
9. S. Wu, H. Guo and C. Zhao, Challenges and Opportunities for Proton Batteries: From Electrodes, Electrolytes to Full-Cell Applications, Adv. Funct. Mater., 2024, 2405401,  DOI:10.1002/adfm.202405401.
10. H. Guo and C. Zhao, An Emerging Chemistry Revives Proton Batteries, Small Methods, 2024, 8(6), 1–23,  DOI:10.1002/smtd.202300699.
11. M. A. Spencer, J. Fortunato and V. Augustyn, Electrochemical Proton Insertion Modulates the Hydrogen Evolution Reaction on Tungsten Oxides, J. Chem. Phys., 2022, 156(6), 064704,  DOI:10.1063/5.0082459.
12. R. Kosydar, M. Kołodziej, E. Lalik, J. Gurgul, G. Mordarski and A. Drelinkiewicz, The Role of Hydrogen Bronzes in the Hydrogenation of Polyfunctional Reagents: Cinnamaldehyde, Furfural and 5-Hydroxymethylfurfural over Pd/H_xWO₃ and Pd/H_xMoO₃ Catalysts, Int. J. Hydrogen Energy, 2022, 47(4), 2347–2365,  DOI:10.1016/j.ijhydene.2021.10.162.
13. N. P. Holzapfel, M. Chagnot, P. S. Abdar, J. R. Paudel, E. J. Crumlin, J. R. Mckone and V. Augustyn, Solution-Phase Synthesis of Platinum-Decorated Hydrogen Tungsten Bronzes for Hydrogen Atom Transfer from Oxides to Molecules, Chem. Mater., 2024, 36, 11684–11696,  DOI:10.1021/acs.chemmater.4c02814.
14. J. M. Mayer, Site Heterogeneity and Broad Surface-Binding Isotherms in Modern Catalysis: Building Intuition beyond the Sabatier Principle, J. Catal., 2024, 439, 115725,  DOI:10.1016/j.jcat.2024.115725.
15. N. G. Altınçekiç, C. W. Lander, J. Yu, A. Roslend, Y. Shao and H. Noh, Proton, Electron, and Hydrogen-Atom Transfer Thermodynamics of the Metal–Organic Framework, Ti-MIL-125, Are Intrinsically Correlated to the Structural Disorder, J. Am. Chem. Soc., 2025, 147(38), 34777–34790,  DOI:10.1021/jacs.5c10498.
16. J. M. Mayer, Bonds over Electrons: Proton Coupled Electron Transfer at Solid-Solution Interfaces, J. Am. Chem. Soc., 2023, 145(13), 7050–7064,  DOI:10.1021/jacs.2c10212.
17. M. Dagar, A. De, Z. Lu, E. M. Matson and A. E. Thorarinsdottir, Implications of Charge and Heteroatom Dopants on the Thermodynamics and Kinetics of Redox Reactions in Keggin-Type Polyoxometalates, ACS Mater. Au, 2025, 5(1), 200–210,  DOI:10.1021/acsmaterialsau.4c00136.
18. Z. Lu and E. M. Matson, Consequences of Substituting Framework Metals in Polyoxometalates on H-Atom Uptake: Proton-Coupled Electron Transfer at a Keggin-Type Polyoxomolybdate, Chem. Commun., 2025, 61(51), 9254–9257,  10.1039/d5cc01894c.
19. Z. Lu, M. Dagar, J. R. McKone and E. M. Matson, Location of Dopant Dictates Proton-Coupled Electron Transfer Mechanism in Vanadium-Substituted Polyoxotungstates, Chem. Sci., 2025, 16(16), 6736–6743,  10.1039/d4sc08452g.
20. D. Shiels, Z. Lu, M. Pascual-Borràs, N. Cajiao, T. V. Marinho, W. W. Brennessel, M. L. Neidig, R. J. Errington and E. M. Matson, Vanadium Substitution Dictates H Atom Uptake at Lindqvist-Type Polyoxotungstates, Inorg. Chem., 2024, 63(49), 23304–23316,  DOI:10.1021/acs.inorgchem.4c03873.
21. A. Klein, Transparent Conducting Oxides: Electronic Structure-Property Relationship from Photoelectron Spectroscopy with in Situ Sample Preparation, J. Am. Ceram. Soc., 2013, 96(2), 331–345,  DOI:10.1111/jace.12143.
22. E. V. Miu, G. Mpourmpakis and J. R. Mckone, Predicting the Energetics of Hydrogen Intercalation in Metal Oxides Using Acid-Base Properties, ACS Appl. Mater. Interfaces, 2020, 12(40), 44658–44670,  DOI:10.1021/acsami.0c11300.
23. M. L. Freedman, The Tungstic Acids, J. Am. Chem. Soc., 1959, 81(15), 3834–3839,  DOI:10.1021/ja01524a009.
24. A. Cohelo, TOPAS-Academic, Powder Diffr., 2007, 312–317.
25. K. Momma and F. Izumi, VESTA 3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data, J. Appl. Crystallogr., 2011, 44, 1272–1276.
26. S. Axnanda, E. J. Crumlin, B. Mao, S. Rani, R. Chang, P. G. Karlsson, M. O. M. Edwards, M. Lundqvist, R. Moberg, P. Ross, Z. Hussain and Z. Liu, Using “Tender” X-Ray Ambient Pressure X-Ray Photoelectron Spectroscopy as A Direct Probe of Solid-Liquid Interface, Sci. Rep., 2015, 5, 1–12,  DOI:10.1038/srep09788.
27. M. E. Grass, P. G. Karlsonn, F. Aksoy, M. Lundqvist, B. Wannberg, B. S. Mun, Z. Hussain and Z. Liu, New Ambient Pressure Photoemission Endstation at Advanced Light Source Beamline 9.3.2, Rev. Sci. Instrum., 2010, 81, 053106,  DOI:10.1063/1.3427218.
28. S. Saeed, S. Fleischmann, T. Kobayashi, Z. Jusys, E. Mamontov, N. C. Osti, N. P. Holzapfel, H. Song, T. Wang, S. Dai, D. Jiang and V. Augustyn, Oxide Acidity Modulates Structural Transformations in Hydrogen Titanates during Electrochemical Li-Ion Insertion, J. Am. Chem. Soc., 2024, 146, 28795–28808,  DOI:10.1021/jacs.4c08063.
29. A. H. Larson, J. J. Motensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dulak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsberg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiotz, O. Schutt, M. Strange, K. S. Thygesen, T. Vegge, L. Villhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, The Atomic Simulation Environment — a Python Library for Working with Atoms, J. Phys. Condens. Matter, 2017, 29, 273002.
30. M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek and G. R. Hutchison, Avogadro: An Advanced Semantic Chemical Editor, Visualization, and Analysis Platform, J. Cheminf., 2012, 4(17), 17,  DOI:10.1186/1758-2946-4-17.
31. T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, CP2K: An Electronic Structure and Molecular Dynamics Software Package -Quickstep: Efficient and Accurate Electronic Structure Calculations, J. Chem. Phys., 2020, 152(19), 194103,  DOI:10.1063/5.0007045.
32. J. Vandevondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach, Comput. Phys. Commun., 2005, 167(2), 103–128,  DOI:10.1016/j.cpc.2004.12.014.
33. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77(18), 3865–3868,  DOI:10.1103/PhysRevLett.77.3865.
34. J. P. Perdew, W. Yang, K. Burke, Z. Yang, E. K. U. Gross, M. Scheffler, G. E. Scuseria, T. M. Henderson, I. Y. Zhang, A. Ruzsinszky, H. Peng, J. Sun, E. Trushin and A. Görling, Understanding Band Gaps of Solids in Generalized Kohn-Sham Theory, Proc. Natl. Acad. Sci. U. S. A., 2017, 114(11), 2801–2806,  DOI:10.1073/pnas.1621352114.
35. Y. Lee, T. Lee, W. Jang and A. Soon, Unraveling the Intercalation Chemistry of Hexagonal Tungsten Bronze and Its Optical Responses, Chem. Mater., 2016, 28(13), 4528–4535,  DOI:10.1021/acs.chemmater.5b03980.
36. Y. Hinuma, H. Hayashi, Y. Kumagai, I. Tanaka and F. Oba, Comparison of Approximations in Density Functional Theory Calculations : Energetics and Structure of Binary Oxides, Phys. Rev. B, 2017,(96), 1–24,  DOI:10.1103/PhysRevB.96.094102.
37. J. VandeVondele and J. Hutter, Gaussian Basis Sets for Accurate Calculations on Molecular Systems in Gas and Condensed Phases, J. Chem. Phys., 2007, 127(11), 114105,  DOI:10.1063/1.2770708.
38. C. Hartwigsen, S. Goedecker and J. Hutter, Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn, Phys. Rev. B:Condens. Matter Mater. Phys., 1998, 58(7), 3641–3662.
39. S. Goedecker, M. Teter and J. Hutter, Separable Dual-Space Gaussian Pseudopotentials, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54(3), 1703–1710.
40. J. D. Head and M. C. Zerner, A Broyden-Fletcher-Goldfarb-Shanno Optimization Procedure for Molecular Geometries, Chem. Phys. Lett., 1985, 122(3), 264–270,  DOI:10.1016/0009-2614(85)80574-1.
41. E. Sanville, S. D. Kenny, R. Smith and G. Henkelman, Improved Grid-Based Algorithm for Bader Charge Allocation, J. Comput. Chem., 2007, 28(5), 899–908,  DOI:10.1002/jcc.
42. H. Guo, L. Pan, M. Gao, L. Kong, J. Zhang, A. Khan, N. A. Siddiqui and J. Lin, Electronic Structure Modulating of W₁₈O₄₉ Nanospheres by Niobium Doping for Efficient Hydrogen Evolution Reaction, Chem.–Eur. J., 2024, 30, e202403043,  DOI:10.1002/chem.202403043.
43. N. W. Thomas, The Compositional Dependence of Octahedral Tilting in Orthorhombic and Tetragonal Perovskites, Acta Crystallogr., Sect. B: Struct. Sci., 1996, 52(1), 16–31,  DOI:10.1107/S0108768195006100.

---