Control of proton transport and hydrogenation in double-gated graphene

Device fabrication

Apertures 10 µm in diameter were etched in silicon nitride substrates (500 nm SiN_x) using photolithography, wet etching and reactive ion etching, as previously reported¹. Source and drain electrodes (Au/Cr) were patterned using photolithography and electron-beam evaporation. Mechanically exfoliated monolayer and bilayer graphene crystals were transferred over the apertures and on the electrodes (Extended Data Fig. 1a). We selected crystal flakes with a rectangular shape, with their long side being several tens of micrometres, to form a conducting channel between the source and drain electrodes. The width of the flake was chosen to be only a couple of micrometres wider than the aperture (Extended Data Fig. 1b), which ensured that a whole cross-section of the flake could be gated with the two gates effectively. The cross-section of the flake became the active area of the device, with the non-gated areas acting as electrical contacts to the gated section. An SU-8 photo-curable epoxy washer with a hole 15 μm in diameter was transferred over the flake and on the source and drain electrodes³ with the hole in the washer aligned with the aperture in the silicon nitride substrate (Extended Data Fig. 1b). The polymer seal ensured that the electrodes were electrically insulated from the electrolyte. The electrolyte used was 0.18 M HTFSI dissolved in polyethylene glycol (number-averaged molecular mass, M_n, of 600)¹⁴, never exposed to ambient conditions. This electrolyte was drop-cast on both sides of the device in a glovebox containing an inert gas atmosphere. For reference, we measured devices prepared in the same way, except HTFSI was substituted for LiTFSI in the electrolyte. Palladium hydride foils (around 0.5 cm²) were used as gate electrodes. The device was placed inside a gas-tight Linkam chamber (HFS600E-PB4) filled with argon for electrical measurements.

Transport measurements

For electrical measurements, a dual-channel Keithley 2636B sourcemeter was used to bias both gates. The applied voltages yielded two proton current signals, top (I_t) and bottom (I_b) channel current, also recorded with a Keithley sourcemeter. These two signals quantified the two halves of the proton-transport circuit: proton transport from one gate electrode towards graphene, and then from graphene towards the other gate electrode. Because of the proton permeability of graphene, these two currents were effectively identical (Extended Data Fig. 3), differing by only about 1%. For this reason, it is sufficient to use only one of them to unambiguously characterize proton transport in the device, I. A second Keithley 2614 A sourcemeter was used both to apply drain-source bias (V_ds) and to measure the electronic conductance (σ_e).

To confirm that the two gates are independent, we connected the electrolyte in the top channel with a reference electrode (PdHx foil, the same size as the gate electrodes), and monitored its potential, V_t^ref, with a Keithley 2182 A nanovoltmeter (Extended Data Fig. 2a). In the first experiment, we swept V_b for various fixed V_t. Extended Data Fig. 2b shows that V_t^ref = V_t for all values of V_b within the experimental scatter of less than 4 mV. This demonstrates that sweeping V_b does not affect V_t. In the second experiment, we swept V_t for various fixed V_b. This measurement also showed that V_t^ref = V_t (Extended Data Fig. 2c), which demonstrates that a fixed V_b does not affect V_t either. These experiments therefore demonstrate that the gates are independent from each other.

To obtain maps of the proton and electronic systems, we measured I and σ_e simultaneously as a function of V_t and V_b. The maps were obtained using software that allowed us to control V_t − V_b and V_t + V_b as independent variables. We swept V_t − V_b (for a fixed V_t + V_b) at a rate of 10 mV s⁻¹ for each gate and stepped V_t + V_b with intervals of 10 mV. The maximum V_t or V_b applied was ±1.4 V, which resulted in a maximum V_t + V_b and V_t − V_b of ±2.8 V. We normally did not apply V bias beyond these voltage ranges to avoid damaging the devices.

Raman spectroscopy

For Raman measurements, the graphene devices were left in the same gas-tight Linkam chamber (HFS600E-PB4) used for electrical measurements. It has an optical window. The Raman spectra of devices were measured as a function of applied V bias using a 514 nm laser. The background signal from the electrolyte was removed for clarity, resulting in relatively weak Raman spectra for pristine graphene (Extended Data Fig. 5). After hydrogenation, a strong D peak appeared and the intensity of the G peak increased while the 2D peak became broader, in agreement with previous work¹⁴. The density of adsorbed hydrogen atoms in hydrogenated graphene was estimated from the ratio of peak-height intensities of the D and G bands^40,41, I_D/I_G. In our devices, I_D/I_G ≈ 1, which corresponds to a distance between hydrogen atoms in graphene of L_D ≈ 1 nm. An equivalent analysis using the integrated-area ratio in these peaks⁴², A_D/A_G ≈ 2, yields L_D < 1.2 nm. Both estimates yield a density of hydrogen atoms of around 1 × 10¹⁴ cm⁻², in agreement with previous reports on hydrogenated graphene^13,14.

Electrolyte characterization

To characterize the limiting conductivity of our devices, we measured devices similar to those described above but in which the aperture in the silicon nitride substrate was not covered with graphene (‘open hole device’). Extended Data Fig. 9 shows that the I–V characteristics of these open-hole devices were linear in all the V-bias range used in this work. This demonstrates that the field effect we observed in graphene devices did not arise from changes in the electrolyte conductivity at high V, consistent with the known large electrochemical window of this electrolyte (4–5 V; ref. ⁴³).

To characterize the capacitance of the electrolyte, we patterned two gold electrodes on a silicon nitride substrate using photolithography and electron-beam evaporation. The electrodes were connected in an electrical circuit and a polymer mask was used to cover all the electrodes, except for an active area that was exposed to the electrolyte. The area of the electrodes differed by a factor of around 50, which ensured that the total capacitance was dominated by the smaller one and allowed us to observe differences, where present, in the response of the devices under positive and negative potentials⁴⁴. Cyclic voltammetry (CV) measurements with scan speeds in the range 1–40 mV s⁻¹ were performed over the voltage range −0.1 V to 0.1 V. Extended Data Fig. 6a shows that the CV curves displayed no redox peaks or asymmetry between the positive and negative voltage branches. The area-normalized capacitance of the electrolyte, C, could then be obtained from the CV curves from the expression^44,45 C = (A × ΔV × ν )⁻¹ ∫ I dV, where A is the active area of the electrode, ΔV is the voltage range in the CV, I is the current and v is the scan speed. Extended Data Fig. 6b shows the extracted C as a function of v. For the smallest scan rate (1 mV s⁻¹), we found C ≈ 30 µF cm⁻². This value decreases with v increases, as expected. Because our measurements use v = 10 mV s⁻¹, we used the value obtained at such v, C ≈ 20 µF cm⁻², in our estimates involving C.

Estimation of E and n

The Debye length in our electrolyte (0.18 M salt and solvent dielectric constant, ε_r ≈ 10)^46,47 can be estimated as λ_D ≈ 0.3 nm. Given this and the relatively large gate bias used in this work, the electrical potential across the graphene–electrolyte interface in our devices dropped almost entirely across the Stern layer. Hence, each of the gate potentials can be described using a parallel plate capacitor model, which we used to derive the relations between the gate potentials and E and n, as shown in refs. ^{25,26,27,48,49}.

To derive the relation between V_t + V_b and n, we note that if only one gate (top) operates, the charge induced is neC⁻¹ = (V_t − V_t^NP), where the superscript NP marks the neutrality point and e is the elementary charge constant. An equivalent relation holds for the top gate. Hence, the total charge from both gates is given by the addition of their contributions: neC⁻¹ = (V_t + V_b) − Δ^NP, where Δ^NP ≡ V_t^NP + V_b^NP. To consider the quantum capacitance of graphene, we note that \({\mu }_{e}={\hbar v}_{{\rm{F}}}\,\sqrt{{\rm{\pi }}n}\), where v_F ≈ 1 × 10⁶ m s⁻¹ is the Fermi velocity in graphene and ħ is the reduced Planck constant. This changes the relation to⁵⁰:

From the estimate of C above, we get (V_t + V_b) − Δ^NP ≈ (0.8 × 10⁻¹⁴ V cm²) n + (1.16 × 10⁻⁷ V cm) n^1/2. Note that this description is accurate only if the Fermi energy of the system is outside a bandgap³¹. Hence, we use it only when graphene is conductive. After the hydrogenation transition, we cannot assess n, as indicated by the breaks in the top axes in Fig. 2.

To derive the relation between V_t − V_b and E, we note that if only one gate (top) operates, the electric field induced by the gate is \({E}_{{\rm{t}}}={en}_{{\rm{t}}}{(2\varepsilon )}^{-1}={C(2\varepsilon )}^{-1}({V}_{{\rm{t}}}-{V}_{{\rm{t}}}^{{\rm{NP}}})\), where ε is the dielectric constant of the solvent. Note that the electric field points in the direction between graphene and its corresponding electrical double layer, which we define as +x for the top gate. An equivalent relation holds for the bottom gate, except that E_b points in the −x direction (towards its corresponding electrical double layer). The total electric field in graphene is then:

This yields E ≈ 1.13 × 10⁹ m⁻¹ (V_t − V_b).

Analytical model of proton transport and hydrogenation in double-gated graphene

We used an analytical model to illustrate how the gate voltages affect proton transport and hydrogenation in double-gated graphene. The energy barrier for proton transport through the centre of the hexagonal ring in graphene is modelled using a Gaussian function: V_p = G₀ × exp(−x/W)², where G₀ = 0.8 eV is the barrier height determined experimentally in the low-electric-field limit¹ and W = 0.5 Å is the barrier width (Extended Data Fig. 7b). For the hydrogenation process, the proton is directed towards the top of a carbon atom in graphene.

The potential energy profile for the hydrogenation process consists of two parts: the energy barrier (V_Hb) and the adsorption well (V_Ha). The energy barrier is modelled with a Lorentzian-type function: V_Hb = V₀ [((x − |x₀|)/d₀)³ + 1]⁻¹, where V₀ = 0.2 eV is the barrier height, x₀ = 1.7 Å is the distance between the barrier and graphene, and d₀ = 0.4 Å is the barrier width. The third power in the denominator models long-range van der Waals interactions. The adsorption well is modelled with a Lorentzian: V_Ha = V₁ [((x − |x₁|)/d₁)² + 1]⁻¹, where V₁ = −0.8 eV is the well depth, x₁ = 1.1 Å is the distance between the well and graphene, and d₁ = 0.25 Å is the width of the well. Note that the well is modelled to be strongly repulsive at x = 0 to capture the repulsion between the carbon and hydrogen atoms at very short distances. The parameters for these functions are taken from DFT calculations^14,51 and the total potential for the hydrogenation process is then V_Hb + V_Ha (Extended Data Fig. 7a).

The gate potential profiles (V_t and V_b) are modelled with a Guoy–Chapman–Stern model²⁴, using dielectric constant ε_r = 10 for the solvent, an electrolyte concentration of 0.18 M and a Stern-layer thickness of 0.4 nm. The resulting gate potentials drop almost exclusively over the Stern layer (Extended Data Fig. 7) and, as a result, the graphene–electrolyte interface behaves as a capacitor, as discussed above. The qualitative findings of our model are relatively insensitive to the specific parameters of the Guoy–Chapman–Stern model if the Stern layer exceeds 0.3 nm. The superposition of the potentials for each of the processes with the gate potentials model the behaviour of the devices.

To illustrate the role of the gates in the hydrogenation process, we set them to yield n = 1.2 × 10¹⁴ cm⁻² but E = 0. Extended Data Fig. 7a shows that this distorts the potential energy profile for hydrogenation, such that the hydrogenation barrier is now easily surmounted by incoming protons, which become trapped in the adsorption well and hydrogenate graphene. To illustrate the role of E in the proton-transport process, we set the gates to produce large E = 1.7 V nm⁻¹ but n = 0 cm⁻². Extended Data Fig. 7b shows that this distorts the potential energy profile for proton transport, such that the barrier is now easily surmounted by a proton moving in the direction of the electric field (from the −x to the +x). To illustrate the role of electron doping in proton transport, we set the gates to give large n = 1 × 10¹⁴ cm⁻² but E = 0.67 V nm⁻¹. Extended Data Fig. 7c shows that this also distorts the potential energy profile for the incoming protons, resulting in facilitated transmission over the barrier. The model illustrates that the distortion of the energy profile for incoming protons due to E and n in these devices is comparable to the barrier height. For this reason, these variables dominate the response of the devices, and previously identified effects, such as strain and curvature, should have a secondary role.

Electrochemical description of the hydrogenation process

The transport data are described using the variables E and n. However, it is equivalent to describe the system using the electrochemical potential of electrons in graphene with respect to the NP, μ_e, instead of n. Indeed, one important property of graphene is that n and µ_e are related by the formula \({\mu }_{e}={\hbar v}_{{\rm{F}}}\,\sqrt{{\rm{\pi }}n}\), where v_F ≈ 1 × 10⁶ m s⁻¹ is the Fermi velocity in graphene. This relation is fundamental, arising from the density of states in the material, and holds exactly in experimental systems^52,53. Moreover, this relation is valid independently of whether the material is gated or not. Hence, the top x axis in the hydrogenation map in Fig. 2a can be re-expressed in terms of µ_e, which illustrates that the hydrogenation process is driven by μ_e. This is consistent with the well-established notion that electrochemical charge-transfer processes are driven by this variable. Note that although the relation between n and μ_e is fixed, applying a gate voltage to graphene shifts both variables^31,54. However, these variables are not independent, as discussed above. To determine their dependence on the gate voltage, we need to establish the electrostatic gate capacitance, C (Extended Data Fig. 6). When C is determined, the dependence of n (or μ_e) on the gate voltage is described by equation (1) above.

Hydrogenation transition

It is instructive to compare our results with previous work on plasma-hydrogenated graphene^13,55. In those earlier studies, plasma-hydrogenated graphene typically displayed around 100-times higher electronic resistivity than in the non-hydrogenated state. By contrast, in our work and in ref. ¹⁴, this factor is about 10⁴, yielding an insulating state that was mostly insensitive to the gate voltage. There are at least two possibilities for this difference. The first is that the hydrogen-atom densities obtained by the different methods are different. Indeed, although the Raman spectra of the current work and ref. ¹³ yield I_D/I_G ≈ 1, this could arise from a hydrogen-atom density of less than 10¹² cm⁻² or around 10¹⁴ cm⁻², because of the bell-shape^40,41 of the graph of I_D/I_G against defect density. To decide which one applies, it is therefore necessary to look for further evidence of disorder in the spectra. The Raman spectra in ref. ¹³ displayed a sharp 2D band, which is typical of ordered samples and indicates that the hydrogen density was likely to be less than 10¹² cm⁻². This contrasts with the 2D band in our spectra, which is smeared, consistent with a figure of around 10¹⁴ cm⁻². The second possibility is that both systems had the same hydrogen density. In this case, the higher resistivity could arise from a more disordered hydrogen-atom distribution. Indeed, hydrogen atoms in plasma-hydrogenated graphene are known to cluster^33,56, which reduces the number of effective scattering centres proportionally to the number of atoms in the cluster. The reduction could be considerable because the scattering radius around each hydrogen atom extends to second neighbours (nine carbon atoms)^33,56,57. The electrochemical system could be less prone to clusters, perhaps because the electrolyte stabilizes the proton as it adsorbs on graphene, making the reaction more likely to happen than in a vacuum, thus yielding a more random distribution.

Another difference between the two hydrogenation methods is their reversibility. According to ref. ¹³, the plasma-hydrogenation process could be almost completely reversed by annealing the material in an argon atmosphere. However, a D band was still notable after annealing and some of the electronic properties of graphene were not fully recovered¹³. This imperfect reversibility was attributed¹³ to the presence of vacancy defects introduced during the plasma exposure. In both our work and in ref. ¹⁴, the hydrogenated transition is fully reversible, with no D peak apparent in the Raman spectra of dehydrogenated samples.

Another difference with plasma-hydrogenated samples is that electrochemical hydrogenation allows the dependence of the transition on n to be studied. This has revealed that the transition is sharp. We attribute this sharpness to a percolation-type transition⁵⁸ triggered both by the high density of adsorbed hydrogen atoms in the samples and to the carrier scattering associated with them^59,60. We propose that the insulating state in the samples is therefore a consequence of their high disorder, as suggested previously^59,60, rather than a bandgap. This is consistent with experimental studies reporting that a bandgap in plasma-hydrogenated samples typically requires either the patterned distribution of hydrogen atoms⁶¹ or a much higher hydrogen-atom density⁶² than in the samples in this work.

DFT calculations of graphene hydrogenation

Graphene hydrogenation was simulated using the Vienna Ab initio Simulation Package (VASP)^63,64,65,66. Electron–ion interactions were modelled using the projector augmented wave method, and the exchange correlations of electrons were modelled with the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation functional⁶⁷. Spin polarization was considered and the van der Waals interactions were incorporated by using the Grimme’s DFT-D3 method⁶⁸. Initial crystal-structure relaxation was performed with a force criterion of 0.005 eV Å⁻¹ and an electronic convergence of 10⁻⁶ eV, accelerated with a Gaussian smearing of 0.05 eV. The energy cut-off was set at 500 eV, and Monkhorst–Pack k-point mesh with a reciprocal spacing of 2π × 0.025 Å⁻¹ was implemented, which ensured energy convergence to 1 meV. We constructed a cubic simulation, consisting of a 4 × 7 orthogonal supercell with 112 carbon atoms placed at the centre in the z direction (perpendicular to the 2D plane) and with a vacuum slab to prevent interactions between adjacent periodic images. After relaxation, the energy barriers for a proton to be adsorbed on top of a carbon atom under vacuum conditions were calculated by ab initio molecular-dynamics simulations using the microcanonical ensemble and the same convergence criteria as mentioned above. We used a time step of 0.1 fs and a minimal initial kinetic energy for the proton in the direction perpendicular to the 2D layer, as previously reported^1,69. A dipole correction was implemented to study the influence of an external electric field perpendicular to the 2D layer (in the z direction)⁷⁰. Owing to the periodic boundary conditions, this dipole is repeatedly inserted in all the simulation boxes in the z direction, yielding a constant electric field in the direction perpendicular to graphene⁷⁰.

Extended Data Fig. 8 shows the calculated potential energy curves for the proton–graphene system. The curves were calculated as a function of distance between the proton and the top of a carbon-atom site with a fully relaxed lattice. The potential energy curves display a minimum (adsorption well) at around 1.14 Å (C–H bond) and a small adsorption barrier around 2 Å, in agreement with previous studies¹⁴. We find the electric field distorts the potential energy profile for hydrogenation, favouring the process in agreement with the analytical model. For reference, we also performed calculations using the non-local optB88-vdW and the hybrid functional HSE06. These resulted only in minor differences (<0.1 eV) in the hydrogenation-barrier height compared with PBE.

DFT calculations of proton transport through graphene

The DFT calculations of proton transport through graphene were performed using VASP^63,64,65,66 and the plane-wave self-consistent field (PWscf) package with Quantum Espresso (QE). We used the optB88-vdW⁷¹ functional, with a 3 × 3 × 3 Γ-centred k-point grid, a 1,000 eV energy cut-off with hard pseudopotentials^72,73, and a force-convergence criterion of 0.03 eV Å⁻¹. We used a 4 × 4 unit cell with a vacuum separating periodically repeating graphene sheets of 12 Å for pristine graphene and around 23 Å for hydrogenated graphene. The zero electric field energy profiles were computed using the climbing-image nudged elastic band method⁷⁴ with VASP. Charged cells were used to describe the protons in the simulations with a uniform compensating background. In the model, proton transfer was simulated from a water molecule on one side of graphene to another one on the opposite side. Using these two water molecules minimizes spurious charge transfer from the graphene sheet to the proton, as confirmed with a Bader⁷⁵ charge analysis. To incorporate the electric field, we modelled the system using QE. Here, we used the optB88-vdW functional^{71,76,77,78,79}, a 3 × 3 × 3 Γ-centred k-point grid and a 600 Ry energy cut-off. We confirmed that the VASP zero electric field energy barriers were reproduced within around 15 meV in QE. The electric field in QE was simulated as a saw-like potential added to the ionic potential, together with a dipole correction implemented according to ref. ⁸⁰. The saw-like potential increased in the region from 0.1 a₃ to 0.9 a₃, where a₃ is the lattice vector perpendicular to the graphene sheet, which was placed at the centre of the cell (0.5 a₃), then decreased to 0 at a₃ and 0. The discontinuity of the sawtooth potential was placed in the vacuum region. The electric field was applied in the perpendicular direction to the graphene basal plane (the z direction). For reference, we also performed calculations using the PBE-D3 functional, which gave comparable results.

We first calculated the energy profile for proton transport through graphene in the absence of an electric field and for two different levels of hydrogen-atom coverage of the lattice (0% and 20%). The choice of 20% hydrogenation was to take into account the fact that adsorbed hydrogen atoms typically form dimer structures consisting of two hydrogen atoms per eight-carbon-atom sublattice^33,56,81, which correspond to a local lattice coverage of about 25%. In agreement with ref. ¹⁸, we observed that the energy barrier for pristine graphene reduced by around 30% for 20% hydrogen-atom coverage. The barrier at zero field we found, Γ₀ ≈ 3.1–3.4 eV for the different functionals, is larger than the typically found values⁷ of Γ₀ ≈ 1–2 eV because in our approach the computed proton trajectory involved a chemisorption state, as described previously¹⁸. However, we note that the absolute values of the barriers in these simplified models are not especially informative, as discussed in ref. ⁶⁹. These models aim to provide only qualitative insights into the influence of E and hydrogenation in proton transport through graphene. Next, we computed the energy profiles along the same pathway used in the zero-E calculations, but now including a perpendicular electric field, E, along the direction of motion of the proton. Extended Data Fig. 10 shows the energy profiles along the reaction path for the two different levels of hydrogenation of the lattice for various electric fields. Regardless of the extent of hydrogenation, we observed a roughly linear barrier reduction when the electric field was switched on, achieving an approximately 20% reduction with E at around 1 V nm⁻¹.

Logic and memory measurements

For logic and memory measurements, we defined V_t and V_b as the IN1 and IN2 signals, respectively, and, guided by the maps of the devices, we systematically explored their proton and electronic responses to different input signals. To test the stability of the memory states as a function of time, the electronic system was pre-programmed into a conducting (dehydrogenated) or insulating (hydrogenated) state applying V_t + V_b = −2.8 V and +2.8 V, respectively. The retention of the insulating state was measured for more than a day with a constant IN1 = IN2 = 0 V, and a reading in-plane V_ds of 0.5 mV was applied for 20 s every 1,000 s. During logic-and-memory measurements, the electronic system was pre-programmed into a conducting or insulating state as described above. We then applied the input signals. The optimal parameters were found to be 0 V and +1.0 V for both IN1 and IN2 signals, because this yields high E but low n and thus enables strong modulation of the proton channel with minimum disruption of the electronic memory state. We found that in these measurements, the potentials at which graphene became hydrogenated were larger than in our transport maps. We attribute this to the fact that the fast sweeping of the gates may be altering the composition of the electrochemical double layer, probably resulting in lower concentrations of protons in the graphene–electrolyte interface and thus requiring higher potentials to hydrogenate graphene given the short timescales of this measurement. To implement the logic-and-memory application, the input signals were applied as a function of time in squared waveform patterns. Low and high gate voltages were defined as the logic inputs 0 and 1, respectively, yielding continuous cycles of different input combinations (00, 01, 11, 10).