Twisto-Electrochemical Activity Volcanoes in Trilayer Graphene

In this work, we develop a twist-dependent electrochemical activity map, combining a low-energy continuum electronic structure model with modified Marcus–Hush–Chidsey kinetics in trilayer graphene. We identify a counterintuitive rate enhancement region spanning the magic angle curve and incommensurate twists in the system geometry. We find a broad activity peak with a ruthenium hexamine redox couple in regions corresponding to both magic angles and incommensurate angles, a result qualitatively distinct from the twisted bilayer case. Flat bands and incommensurability offer new avenues for reaction rate enhancements in electrochemical transformations.

Enhancing electrochemical reaction rates is critical for chemical transformations [1], energy conversion, and storage [2].Electrochemical reactions at electrodeelectrolyte interfaces are often controlled by modifying the substrate and thereby tuning the adsorption energy to reach peaks of activity "volcanoes" [3].Another lessstudied approach to enhancing electrode-electrolyte reactions is through modifying the electronic density of states (DoS) of the electrode such that larger numbers of states coincide with the transitions of the redox couple [4,5].For instance, it has been shown that in single-molecule reactions with metallic electrodes (such as gold or copper), rates can increase with increasing surface DoS near the Fermi level [6,7].
Interfacial chemical reactivity of few-layer graphene is significant for nanoscale electrochemical devices and electrocatalysis [8][9][10].It offers the possibility of controllable enhancement of the DoS by modifying the twist angle [11,12].In a recent work, we showed that twisted bilayer graphene (tBLG) can be used to control the reaction rates with a canonical redox couple, producing an activity volcano of an order of magnitude enhanced rates relative to the untwisted Bernal stacked graphene [13].However, the reaction rate maximum was confined close to the magic angle (MA) 1.1 • , which would present challenges in fabrication of practical functional devices.
In this work, we analyze the electrochemical activity contour in twisted trilayer graphene (tTLG), which has unique electronic and mechanical relaxation properties associated with incommensurability (i.e. the induced moiré cells of the individual pairs of layers cannot be made to coincide as integer multiples) [14,15] compared to tBLG.We use a general momentum-space model that incorporates these properties to compute the DoS as a function of twist [14].We then map the electron transfer rates given by the overlap of DoS with the redox couple states and electron/hole occupations under an appropriate theory framework, which is the so-called Marcus-Hush-Chidsey model [5] (explained in more detail below), modified for the case of low-dimensional electrodes.Finally, we modulate temperature and redox couple parameters to adjust enhancement and magnitude shifts in the activity volcano.
We show that the rate enhancement is not necessarily localized at the MA in twisted trilayer systems, which can allow better flexibility over a range of electrochemical and redox couple parameters.We find that the rates can increase at angles away from the MA where nonoverlapping bands from incommensurability increase the number of interacting states at higher energies.We can "trade off" the rate enhancement at these two regions by varying the temperature, where both features become equally significant near the standard condition, forming an activity volcano.We find retention of activity at optimal incommensurate angles over MA at higher redox couple potentials.Our analysis shows that low reorganization energy and high temperature augments the participating solvent states and improves the reaction rate magnitude.
A general scale for electrochemical reaction rates is given by the heterogeneous rate constant (k 0 ), typically computed via the Butler-Volmer (BV) equation [16,17].At higher rates, the linearity assumption in BV breaks down and Marcus theory is needed, which utilizes a second-order free energy expansion [17][18][19][20][21]. Chidsey introduced a further modification to account for electron and hole occupation functions, known today as the Marcus-Hush-Chidsey (MHC) model [22,23].The MHC model assumes a constant, energy-independent density of states.This assumption is relaxed in its extended versions [5,24] for better predictions at high overpotentials [25] and applications like reactivity of graphene edge states and defects [26,27] or lithium stripping and electrodeposition [5].
For materials with sharply varying DoS such as flatband twisted systems, it is crucial to incorporate DoS effects on the interfacial kinetics with an electrolyte.Due to the large number of states at the magic angle, twisted graphene has the potential for exceptionally high rate enhancement comparable to bulk graphite [13].We propose that a significant enhancement of the reaction rate takes place when the redox couple energy levels line up with the energy of electrode states whose density has been increased by the tTLG.
We study the electron transfer dynamics between twisted trilayer graphene in contact with a redox couple (R), approaching from the z direction, as illustrated in the schematic (Fig. 1a).We use ruthenium hexamine (Ru +3 (NH 3 ) 6 , RuHex hereafter) in KCl solvent as our reference redox couple whose formal potential is closest to the charge neutrality point (CNP) of layered graphene [13].
The extended MHC theory [5,26] predicts the reaction rate given the reorganization energy (λ), applied overpotential (η), temperature (T ) and DoS (D(ϵ)) of the electrode, where f is the Fermi-Dirac distribution.The integrand has three terms: redox couple states ("Marcus-like" terms), electron or hole occupation probabilities, and the DoS (illustrated in Fig. 1(b)).Assuming area-normalized rates, the proportionality constant is the squared energyindependent electronic coupling term between the electrode and the redox couple [7,27].Calculation of electronic coupling would require an electronic structure of the ground-state atomic configuration, which is not feasible for the large moiré supercells characteristic of smallangle twisted multilayer graphene.
In low-dimensional electrodes with energy dependent DoS [28], we account for the relatively small and potential-dependent quantum capacitance C q , which creates a dynamic doping effect [26,29], by shifting the argument of the DoS by an energy of eV q as indicated in Eq. 2. In bilayer and trilayer graphene, the dielectric capacitance is high enough that this effect can be neglected [30], as the voltage drops mainly across quantum and the electrical double layer (EDL) capacitance.At equilibrium, eV q subtracts from the formal energetic difference between the redox couple and electrode potentials (E 0 ), while the rest of the voltage drops across the EDL, V dl = E 0 − V q .Assuming the rigid-band approximation, these quantities are obtained from the following coupled system of equations [13,31,32], where ∆E is the difference in the Fermi levels, F t (Eq.4) is the thermal broadening function, and C q and C dl are aggregated quantum and EDL capacitances, respectively.Previous calculations on tBLG [13] suggest that a greater fraction of E 0 goes into into V q near the CNP from flat electronic bands.This contributes to a small offset in the rate enhancement away from the MA, also observed in tBLG [13].
The product of the first two terms (redox couple states and Fermi function) in the integrand of Eq.2 almost approximates a Gaussian filter being convolved with the DoS.As illustrated in Fig. 1c, when this filter is closer to the flat-band energy, the reaction rate is higher due to a greater overlap with the DoS.A number of control parameters can be tuned that change the shape and position of the overlapping terms in the integrand.Fig. 1(b) highlights the tunable parameters of the solvent (E 0 , λ) and the electrode (θ 12 , θ 23 , η).The temperature (T ), broadens both the redox couple states and the Fermi function.
Using this model, we have shown a qualitative agreement with experiments in twisted bilayer graphene (tBLG) and RuHex, where the rate enhancement is observed near the MA (1.1 • ) [13].Notably, the rate enhancement is an order of magnitude (∼ 10×) higher than the theoretical prediction, possibly due to the twist dependence of the electron coupling constant not included in the model.Modifications to the coupling constant like those studied on graphene edges [33,34] and metallic electrodes [35,36] can help narrow the gap between theory and experiment.
In the trilayer system, we employ a low-energy momentum-space continuum model for electronic structure and DoS calculations [14].The model includes outof-plane relaxation, which is known to affect the DoS more than the in-plane components.Additionally, when the twist angles are both ≥ 1 • , the relaxation magnitude is small and the quantitative difference between [37,38].The model utilizes one twist angle between each adjacent pair of layers (θ 12 and θ 23 , measured in opposite senses but both defined to be > 0, see Fig. 1a).
The degree of twist modulates the van Hove singularity (VHS) separation and associated DoS peaks (Fig. 1b).Fig. 2a shows the magnitude of VHS DoS peaks at a range of twist angles (reproduced from Zhu et al. [14]).The high-magnitude peaks asymptotically converge to tBLG MA (1.1 • ) when the other layer is decoupled (twisted at a large angle >4 • ).Plugging the tight-binding DoS into the MHC model, the reaction rate contour (k 0 ) for RuHex (Fig. 2b) contrasts starkly with the DoS peak contour.The color magnitudes are shown relative to the reference ABA system with dispersive bands for direct visualization of the rate enhancement.The comparison assumes no change in electronic coupling with the redox couple between twisted and untwisted systems.Unlike tBLG, the rate enhancement is not confined to the MA curve, but spans a sizeable triangular area (volcano) as indicated in Fig. 2b.
In this region, the value of k 0 is ∼ 2.1× higher than untwisted ABA, which is similar to the theoretical enhancement reported for tBLG (2.2×) [13] with respect to the Bernal stacked AB graphene.The magnitude of k 0 in the tTLG volcano is 3.2× higher than tBLG at 1.1 • , as expected from the flatter DoS of the former (comparing Fig. 3 and S4).Experimentally, this could translate to more than an order of magnitude enhancement in tTLG compared to tBLG.We also show the volcano peak across a 45 • arc passing over the enhancement region (Fig. 2b inset).
The activity volcano starts at the MA curve (line AD, Fig. 2b) where the DoS peaks are the sharpest.Going horizontally from A to C, the VHS separate towards higher energies (Fig. 3), similar to tBLG.In tTLG however, there is an additional change in the band hybridization from commensurate to incommensurate angles, as shown in Fig. 1b in [14].These changes are frequent from A to B, where B is at the dominant (2,1) moiré harmonic line (Fig. 1b [14]).After B, there is a significant region of incommensurate angles until the (1,1) harmonic, in which many competing harmonics coexist, and non-overlapping bands occur over different energy ranges.These bands show an increased DoS in the energy range ±0.25 eV (Fig. 3), where the DoS has an optimal overlap with the MHC filter (Fig. S1).Consequently, the DoS area within ±0.25 eV increases from B to C (Fig. S6), producing a high value for k 0 in this region.After C, the incommensurate states fall outside the filter range, whereby k 0 disappears quickly.These non-overlapping bands are not found in polytypes of tTLG with unbroken commensurability (Fig. S4) like "monolayer-twistbilayer" (MtB) and alternating twist (AtA), and are expected to enhance rates only near the MA.Compared to tBLG, these polytypes have flatter bands from correlated electronic phases (Fig. S4) and are predicted to show an order of magnitude higher value of k 0 in the experiments.
Notably, the vertex of the volcano (D) starts at the point where the (2,1) harmonic line meets the MA curve [14].Therefore ∆ABD is enhanced by the flat bands, and ∆BCD by non-overlapping bands from incommensurate angles.Hence we observe a shift in maximum rate enhancement region from DoS peak (A at low T ) to DoS area (C at high T ).Near the standard conditions (250-300 K), the whole area is kinetically enhanced because both factors are equally significant.The dominance of either factor can be switched as the temperature changes.
Assuming weak T dependence of λ, DoS and electronic coupling, the absolute magnitude of rates increases rapidly with temperature as the Marcus factor increases.However, there is a shift in the relative rate enhancement from the MA curve, ∆ABD (Fig. 4a) to ∆BCD at higher temperatures (Fig. 4b).Also, the maximum rate relative to the untwisted system decreases from 3× at 100 K to 2× at 600 K.This is due to the small filter width at lower temperatures [27] (Fig. S1), enabling an efficient overlap with the flat bands at the MA.At higher temperatures the filter overlap with DoS is greater in a higher energy range (±0.25eV), which is accomplished by non-overlapping bands in ∆BCD highlighted at 600 K.Note that at 100 K the intensified region is at a small offset from the MA due to the tight overlap with the DoS offset from non-zero quantum capacitance.Hence with a given redox couple, temperature modulations can highlight different features of the DoS in reaction rate contours.
The redox couple parameters (E 0 and λ) can be adjusted to study the sensitivity and possible rate improvement.Due to the low energy range of the DoS, we do not show the effect of shifting the filter peak away from E 0 = ±0.3eV.Higher values of E 0 would reduce the overlap with the flat bands as reflected by the gradual shift in the rate enhancement region away from the MA curve at E 0 = 0.2 eV and 0.3 eV (Fig. S5).Incommensurate angles, however, retain high activity near C due to the non-overlapping states filling the band gaps.We observe a shift in the rate enhancement to higher incommensurate angles corresponding to the purple region of the dominant moiré harmonic plot (Fig. 1b of [14]).Hence, previously inactive redox couples in tBLG like CoPhen (Cobalt Phenanthroline) and FcMeOH (Ferrocene Methanol) [13] (E 0 ∼0.3 eV) may achieve comparable activity with RuHex at the incommensurate angles.
Tuning the reorganization energy, however, does not produce the same shift in enhancement between MA and the incommensurate region.As shown in Fig. S2, the filter peak reduces with λ, but unlike with temperature, there is no effect on the width.Hence, the rate magnitude decreases, but relative to the untwisted system, there is almost no change to the volcano with increasing λ (Fig. S3).These observations suggest the possibility of engineering pairs of redox couples and solvents with low λ and high T for high values of k 0 .To boost the effect from flat bands, employ low E 0 (like RuHex) and T which selectively enhance the twist angles near the MA curve.In comparison, the rate enhancement at 100K is higher (∼3) and are closer to the magic angle curve than at 600K.Small filter width at 100K magnifies the effect of flat band singularities, while larger width at 600K promotes higher area under the DoS from non-overlapping bands.
In conclusion, we have explored the role of incommensurability in identifying regions of unexpected rate enhancement in tTLG, which was otherwise locked at the magic angles in commensurate systems like tBLG, MtB and AtA systems.Incommensurability introduces nonoverlapping bands that increase the number of interacting states at higher energies, thus magnifying equilibrium rates away from the MA curve.Temperature variations can adjust rate enhancements between the MA and the incommensurate angles; low temperature contracts the MHC filter thereby favouring flat bands at the MA, and high temperature diffuses the filter favouring greater DoS area at incommensurate angles.With the RuHex redox couple, we observed a high activity volcano spanning both regions at room temperature.Incommensurate angles in twisted trilayer can retain activity with other redox couple like CoPhen and FcMeOH, unlike commensurate systems.The role of twist dependence on the electron coupling constant will be explored in future studies.This analysis shows the immense flexibility offered by tri-layer systems in maintaining high activity for a wide range of electrochemical processes.The product of the Marcus and Fermi functions inside the rate integrand is shown at different temperatures (fig. 1) and reorganization energies (fig.2).This product overlaps with the density of states to give the overall rate integrand, and is referred to as the MHC filter in the main text.The filter peak and width increases with temperature.The optimal energy range for high overlap of the density of states with the MHC filter is ±0.25eV.
Increasing the reorganization energy reduces the peak but the energy spread of the product is not affected.The effect of this is exemplified by the almost identical equilibrium rate maps for two reorganization energies (λ = 0.2eV, 1.0eV) at standard values for RuHex (T =295K, E 0 =-0.07V,η=0) as shown in fig. 3.
The DoS peaks are flatter in tTLG than in tBLG (fig.3 main text vs fig.4(a)), the equilibrium rate therefore is higher (×2) in the former system at its MA.Like tBLG, polytypes of tTLG, MtB and AtA are commensurate at all twist angles, and do not show filled band gaps from non-overlapping states at higher energies (∼ ±0.25 eV) away from their MA (fig. 3 main text vs fig.4).As a result, maximum activity in commensurate systems is localized at the MA.Comparison of equilibrium rate maps at higher values of E 0 is shown in fig. 5.As the magnitude of E 0 increases, the overlap with the flat bands decreases, which shifts the rate enhancement region to higher incommensurate twist angles, while retaining activity at the incommensurate angles (fig.5).The rate enhancement region covers the incommensurate angles marked by the purple region of the dominant moiré harmonic plot shown in fig.1b

FIG. 1 .
FIG. 1. Schematic of the reaction system (a) with the trilayer graphene twisted by θ12, θ23, exchanging electron from a redox couple (RuHex) approaching from z.The extended Marcus-Hush-Chidsey model (b) has three terms in the integrand; the Marcus term, Fermi function, and the DoS of the electrode as shown from left to right.All tunable parameters are labelled in (b).Marcus term and the Fermi function illustrated with the same color in (b) multiply to give either oxidation (red) or reduction (blue) filters, identical at equilibrium as indicated in yellow (c).The filter overlaps with the DoS (black) shifted by quantum capacitance voltage Vq.The equilibrium rate is proportional to the area under the product of the three terms.

FIG. 2 .
FIG. 2. Maximum DoS peak (reproduced from Zhu et al. [14]) (a).(b) k0 rate map for RuHex (E0 = −0.07eV, λ = 0.82 eV, η = 0 eV) as a function of the two independent twist angles in tTLG system.The "magic angle curve" is shown by the black line.A number of weak hotspots in k0 (b) occur near the diagonal due to numerical artifacts in the DoS.The solid white triangle encloses the activity volcano, also marked by A, B, C and D. Inset shows variation of k0 across a 45 • arc as indicated.All rates are relative to the untwisted ABA system.Maximum rate enhancement in the volcano is 2.1× over ABA.

FIG. 4 .
FIG. 4. Equilibrium rate constant as a function of temperature; formal potential E0 = −0.07eV and reorganization energy λ = 0.82 eV set for RuHex solvent at (a) T = 100K and (b) T = 600K.The solid black line marks the magic angle curve for tTLG.In comparison, the rate enhancement at 100K is higher (∼3) and are closer to the magic angle curve than at 600K.Small filter width at 100K magnifies the effect of flat band singularities, while larger width at 600K promotes higher area under the DoS from non-overlapping bands.

FIG. 1 .
FIG. 1.The product of Marcus and the Fermi functions at different temperatures, which acts as an energy filter of the density of states.

FIG. 3 .
FIG. 3. Equilibrium rate maps for (a) λ = 0.2 eV and (b) 1.2 eV for RuHex solvent.The rate enhancement region is almost unaffected since the integrand has a constant energy spread with λ.