Observation of Kekulé vortices around hydrogen adatoms in graphene

Fractional charges are one of the wonders of the fractional quantum Hall effect. Such objects are also anticipated in two-dimensional hexagonal lattices under time reversal symmetry—emerging as bound states of a rotating bond texture called a Kekulé vortex. However, the physical mechanisms inducing such topological defects remain elusive, preventing experimental realization. Here, we report the observation of Kekulé vortices in the local density of states of graphene under time reversal symmetry. The vortices result from intervalley scattering on chemisorbed hydrogen adatoms. We uncover that their 2π winding is reminiscent of the Berry phase π of the massless Dirac electrons. We can also induce a Kekulé pattern without vortices by creating point scatterers such as divacancies, which break different point symmetries. Our local-probe study thus confirms point defects as versatile building blocks for Kekulé engineering of graphene’s electronic structure.


INTRODUCTION
Real-space topological defects in crystals exhibit exotic electronic properties 11,12 , especially when combined with the reciprocal-space topological phase hosted by the bulk 13,14 .In two-dimensional hexagonal lattices, a vortex in the Kekulé order parameter is of particular interest for charge fractionalisation without breaking timereversal symmetry 4,5 .It is a two-dimensional extension of the charge fractionalisation at domain-wall solitons in the one-dimensional spinless model of polyacetylene 15,16 .Introducing the twofold spin degeneracy, a charged topological defect carries zero spin, while the neutral defect carries spin 1 ⁄2.The fractionalisation mechanism is therefore manifested by this unusual spin-charge relation.Similar to the dimerisation in one-dimensional polyacetylene, the Kekulé order in graphene corresponds to the √ 3 × √ 3R30 • unit cell tripling, with a distinct bond order within one of the three equivalent hexagonal rings.These three degenerate states define an angular order parameter space as shown in Figure 1a.A Kekulé vortex of winding 2π corresponds to the alternation of the three Kekulé domains upon encircling a singularity, which could be implemented in optical 17,18 and acoustical 19,20 metamaterials.In graphene, recent experiments have demonstrated that a Kekulé bond order emerges from the intervalley coherent quantum Hall ferromagnet states in the zeroth Landau level [21][22][23] .These states can host skyrmionic topological excitations appearing as Kekulé vortices [23][24][25] .However, these schemes require a strong magnetic field and the Kekulé vortex without breaking time-reversal symmetry remains out of reach.At zero magnetic field, theory shows that a missing electronic site provides a fractionalisation mechanism analogous to that of the Kekulé vortex. 26,27Such point defects can scatter electrons from one valley to another, thereby connecting the two valleys and offering a way to stabilize the Kekulé bond order [6][7][8][9][10] .We now demonstrate that they also induce Kekulé vortices in graphene.
Figure 1b shows a scanning tunneling microscopy (STM) image of graphene with a single chemisorbed hydrogen adatom (see Methods for experimental details and Extended data Fig. 1 for another example.).In this image, we have filtered out the intense signal of the close-tozero-energy state induced by the adatom [29][30][31] , in order to highlight the result of elastic intervalley scattering of the massless relativistic electrons by the adatom 32 .The image reveals a hexagonal superlattice commensurate with that of graphene but with a unit cell three times larger, arXiv:2307.10024v1[cond-mat.mes-hall]13 Jul 2023 FIG. 1. Observation of the Berry-Kekulé vortex induced by hydrogen adatom on graphene.a, Schematic illustration showing three distinct Kekulé orders in graphene.b, STM image of graphene with a hydrogen adatom in the center for a tip bias V b =400 mV and tunneling current it = 45.5 pA.The imaged area of 7.4×7.4nm 2 has been filtered (see Extended Data Fig. 2 for raw images) to remove the low frequency signal associated to the H atom which position is given by the white disk.The bottom right inset shows selected harmonics.The colored tiling evidences three domains, each corresponding to one of the three Kekulé orders (bold hexagons) and separated from one another by domain walls along the armchair direction.The outer colour wheel (shown in the bottom left) labels the winding of the Kekulé order parameter around the hydrogen adatom.The coloured dots highlight the usual on-site quasi-particle interference signal 28 , whose 4π winding is highligted by the central colour wheel (bottom left).The pseudo-spin of incoming electrons scattering off the H atom is locked on the azimutal coordinate θr of the STM tip 28 (upper right inset).c, Filtered STM image with the Kekulé harmonics only (see inset) for comparison with the LDOS calculated using the Green's function approach.The imaged area is 4.8×4.8nm 2 .d, LDOS around a hydrogen adatom calculated using the Green's function approach and integrated between 0 and 400 meV in order to provide comparison to the experiment.The simulated image includes on-site and bond contributions (see the main text, Extended Data Fig. 3 and Supplementary Information for details).e, Clar's sextet configurations of graphene in presence of a hydrogen adatom (black dot) illustrating the emergence of a Kelulé vortex.
as emphasised by the three-colour tiling (See Extended data Fig. 2 for details).The pattern is dominated by the onsite signal on sublattice B (assuming the H adatom is on sublattice A), which alternates between the three nonequivalent B sites.The coloured dots and the inner colour wheel in Fig. 1b highlights that this onsite signal winds 4π when circling around the adatom, consistent with previous studies 28,33,34 .A closer inspection of the interference pattern further reveals the underlying signal on graphene bonds, which is consistent with the Kekulé ordering.In Figure 1b, bold coloured hexagons identify three distinct Kekulé domains separated from one another by domain walls along the armchair direction.Remarkably, the Kekulé bond order shows a 2π winding when circling around the adatom.This demonstrates that the chemisorbed H adatom induces a Kekulé vortex on the surrounding graphene bonds.

ESTABLISHING THE ELECTRONIC ORIGIN OF THE KEKUL É VORTEX
We will now show that the Kekulé vortex has a purely electronic origin related to the band topology of graphene rather than to defect-induced atomic displacements.The bond ordering of the Kekulé type relates to the electron density between carbon atoms, which are no longer equivalent due to the presence of the hydrogen adatom.Intervalley scattering resolved at the tip position r yields the following bond contribution to the local density of states (LDOS): where φ n = (2n+1)π/3 and ∆K n is a scattering wavevector between the valleys responsible for the hexagonal superlattice pattern (see Supplementary Information).This bond contribution can be understood as a two-path loop interference allowed by the overlap of neighbouring p z orbitals of carbon atoms.The vortex originates from the polar angle θ r indexing the tip orientation, which is known to be a real-space representation of the electron pseudospin θ q defined from a Dirac point in momentum space (see inset in Fig. 1b and Ref. 28).Thus, the bond contribution to the LDOS exhibits a 2π vortex centered on the adatom as a reminiscence of the topological pseudospin vortex, or equivalently the topological Berry phase, which characterizes the band topology of the massless Dirac electrons 35 .The usual onsite LDOS modulations ∆ρ A (r) and ∆ρ B (r) due to interference also contribute to the STM signal 28 .Adding all contributions leads to the energy-integrated LDOS modulations in Fig. 1d, which shows very good agreement with the filtered experimental image shown in Fig. 1c.We would like to point out that the experimentally observed vortex is additionally reproduced by both our tight-binding and density functional theory (DFT) calculations (see Methods, Fig. 2 and extended data Fig. 5).While the Fourier filtering of experimental data does not allow to separate onsite and bond contributions to the LDOS, this can be readily achieved in calculations (see Extended Data Figs. 3).We note that the vortices are not affected by the energy integration of the LDOS, since the scattering wave-vectors ∆K n are energy independent.Energyresolved STM images shown in Extended Data Fig. 4 confirm this conclusion.
The Kekulé vortex in the LDOS derives from an intrinsic topological property of the massless Dirac (that is, chiral in pseudospin 35 ) wavefunctions scattered by the adatom.Therefore, we name it a Berry-Kekulé vortex.At this point, the following two fundamental questions arise: i) Is the Berry-Kekulé vortex necessarily accompanied by a structural Kekulé distortion that would lead to opening a spectral gap? ii) If not, is the Berry-Kekulé vortex accompanied by zero-energy fractional excitations despite the absence of a spectral gap?
To answer the first question, we perform DFT calculations that intrinsically take into account lattice relaxation effects (see Methods).The simulated STM image (Fig. 2a) reproduces accurately the experiment, and in particular the 2π vortex.While the presence of the hydrogen adatom does induce minor lattice distortions that are consistent with the winding of the structural Kekulé distortion (Fig. 2b,c), the results are essentially the same as the ones provided by our analytical description in Eq. 1 and tight-binding calculations (Extended Data Fig. 5) that do not include any lattice distortion effects whatsoever.Furthermore, performing DFT calculation without any relaxation does not affect the presence of the vortex (Extended Data Fig. 6).This further confirm the electronic origin of the Berry-Kekulé vortex.
Interestingly, the Kekulé vortex can be interpreted in terms of Clar's sextet theory 36 , a set of empirical rules widely used in the chemistry community that explains graphically the stability of aromatic molecules.The sextet representation depicts the six resonantly delocalized π electrons by a circle.Clar's rules state that adjacent hexagons can never be aromatic sextets simultaneously, and that the most stable bond configuration maximizes the number of Clar's sextets.Graphene admits three equivalent resonant Clar configurations corresponding to three Kekulé orders (Fig. 1a) 37 .In the presence of a hydrogen adatom that effectively removes one site from the lattice, the circulation of delocalized electrons in the adjacent benzene rings is obstructed, which lifts the degeneracy of Clar's resonant configurations allowed in its surrounding.We then find two possible configurations shown in Fig. 1e.The only freedom lies in the positions of the double covalent bonds along the boundaries between two Kekulé domains.Clar's sextet representation leads to the same conclusion: a hydrogen adatom induces a Kekulé vortex in graphene.In supplementary infomation, we formalize the relation between these two pictures.
We will now answer question ii) with regard to fractional excitations accompanying Kekulé vortices, first for spinless electrons.The hydrogen adatom forms a covalent bond with a carbon atom in graphene lattice, eliminating the p z atomic orbital as a result of rehybridization into the sp 3 state of carbon atom 30 .In this sense, a hydrogen adatom is equivalent to a single-atom vacancy in graphene except for stronger structural reconstruction in the latter case.Such a strong potential promotes half a state from the valence band to zero energy, and half a state from the conduction band if the spectrum is particle-hole symmetric.The zero-energy state is fully polarized on the sublattice opposite to that of the removed p z orbital and presents an algebraic decay due to the gapless relativistic spectrum [29][30][31] .This quasilocalised state at zero energy has the fractional charge −e/2, which corresponds to the charge difference of the occupied valence band with and without the hydrogen adatom 26 .The number of quasi-bound states relates directly to the pseudospin real-space representation θ r encoded in the scattering wavefunctions 26 .Therefore, the topological winding of the Berry-Kekulé vortex we observe at higher energy and the number of fractional excitations at zero energy are intrinsically intertwined.
Introducing the spin degeneracy doubles the number of zero-energy states, one per spin.At half filling, this implies that one of the two spin-polarised states is fully occupied.Since each bound state is a hybrid superposi-tion of valence-conduction half states of −e/2 and +e/2 fractional charges, the quasi-bound state around the hydrogen adatom must be neutral with spin 1 /2.Remarkably, this very unusual spin-charge relation appears to be supported by two independent experimental observations, one bringing evidence of the neutral charge 38 , and the other of the magnetic moment 39 .It is also consistent with Lieb's theorem 40 .

KEKUL É ORDER NEAR DIVACANCIES
The existence of the Kekulé vortex and the quasibound states induced by the hydrogen adatom defect is related to the symmetries of this scattering center.While a hydrogen adatom breaks the C 2 symmetry and sublattice balance of graphene, a divacancy preserves these symmetries.Figure 3a shows an STM image of such divacancy in graphene.The divacancy also induces locally a Kekulé order.Importantly, this Kekulé order is exempt of any winding, with the Kelulé domain being defined by one of the three orientations of the chemical bond linking the two removed atoms.This is confirmed by DFT calculations (Fig. 3c) which show that the defect also creates significant lattice distortion originating from structural reconstruction due of the formation of two pentagonal rings (Figs.2c and Fig. 3d).
A local Kékulé order with or without vortex can therefore be induced in graphene by the specific distribution of atomic defects.Harnessing experimentally these building blocks in a novel type of defect engineering could lead to the long-awaited macroscopic Kekulé order from the cooperative effect of atomic defects 6,7 .
The samples were grown by thermal decomposition of the carbon-face SiC at temperatures close to 1150 • in ultrahigh vacuum. 41Silicon evaporation results in several graphene layers decoupled by rotational disorder.Hydrogen atoms were then deposited by thermal dissociation of hydrogen gas in a custom atomic hydrogen source as described previously. 39STM images were obtained in the constant current mode in a custom ultrahigh vacuum setup at 5 K.

Tight-binding calculations
The nearest-neighbor tight-binding Hamiltonian of graphene is expressed as = t 0 1 + e ikx + e iky 1 + e −ikx + e −iky 0 , where the nearest-neighbor hopping integral t = −2.7 eV.The TB calculation is carried out with periodic boundary conditions, with the hydrogen adatom modelled as a large on-site potential (V = 100|t| ≈ 270eV).The supercell size of 27 × 27 unit cells of graphene was used to reduce the spurious effects due to the mutual interference between the periodic images of adatoms.The information about the phase difference between wavefunction amplitudes on the nearest-neighbor sites, that is the bond order, is needed in order to describe the Kekulé order.We define the bond parameter as the LDOS between two neighboring atoms i, j.On the basis of TB model, we consider the wave function as the product of the TB eigenvector ψ and an envelope function f (r) φ(r) = n i φ n i f (r − r i ), taking the middle point between the atoms r ij = (r i + r j )/2 LDOS writes ρ(r ij ) = |f (a 0 /2)(ψ i + ψ j )| 2 = f (a 0 /2)(2Re ψ i |ψ j + ψ * i ψ i + ψ * j ψ j ), (3)   in which the inner product ψ i |ψ j governs the modulation of bond order.Therefore, we use the orbital overlap Re ψ i |ψ j as the bond-order operators in the TB calculation.In the presence of a scattering center, ψ i |ψ j is perturbed by the inter-sublattice Green's functions G AB and G BA (See SM 42 ), since the atoms connected by the bond are in different sublattices.The bond order is plotted by integrating the LDOS from 0 to 300 meV.

DFT calculations
First-principles calculations were performed using the SIESTA code 43 .We use the double-ζ plus polarization localized orbital basis set combined with the local density approximation exchange-correlation functional 44 .The energy shift for constructing the localized orbital basis functions was set to 275 meV, and the real-space cutoff to 250 Ry.The structural relaxation was performed using the conjugate gradient method.
The simulated STM images were produced using the plstm module of the SIESTA package, as a postprocessing step following the DFT calculations.The images were simulated from LDOS at a constant height of 2.5 Bohrs above the graphene plane.

FIG. 2 .
FIG.2.Kekulé vortex induced by the hydrogen adatom from first principles.a, STM image simulated using DFT with a tiling defined in the same way as in Fig.1b.The calculated LDOS was integrated between 0 and 400 meV as in the experiment.The image area is 4.2×4.2nm 2 .b, Relaxed atomic structure of graphene with a hydrogen adatom.The bond lengths are coded as the color and thickness of the bonds.The colored tiling is superposed on the graphene lattice to show the winding of the bond length.The area shown is 3×3 nm 2 .c, Calculated bond length distribution as a function of distance to the hydrogen adatom and divacancy defect.

FIG. 3 .
FIG.3.Kekulé order induced by the divacancy defect in graphene.a, An STM image of a divacancy in graphene (V b = 500 mV, it = 400 pA).The image area is 6.2×6.2 nm 2 .b, Fourier filtered image (same filter as in Fig.1b).c, DFT simulated STM image.d, Atomic structure of the divacancy defect in graphene from DFT calculations.The area shown is 3×3 nm 2 .