Magnetism of topological boundary states induced by boron substitution in graphene nanoribbons

Graphene nanoribbons (GNRs), low-dimensional platforms for carbon-based electronics, show the promising perspective to also incorporate spin polarization in their conjugated electron system. However, magnetism in GNRs is generally associated to localized states around zigzag edges, difficult to fabricate and with high reactivity. Here we demonstrate that magnetism can also be induced away from physical GNR zigzag edges through atomically precise engineering topological defects in its interior. A pair of substitutional boron atoms inserted in the carbon backbone breaks the conjugation of their topological bands and builds two spin-polarized boundary states around. The spin state was detected in electrical transport measurements through boron-substituted GNRs suspended between tip and sample of a scanning tunneling microscope. First-principle simulations find that boron pairs induce a spin 1, which is modified by tuning the spacing between pairs. Our results demonstrate a route to embed spin chains in GNRs, turning them basic elements of spintronic devices.

Graphene nanoribbons (GNRs), low-dimensional platforms for carbon-based electronics, show the promising perspective to also incorporate spin polarization in their conjugated electron system. However, magnetism in GNRs is generally associated to localized states around zigzag edges, difficult to fabricate and with high reactivity. Here we demonstrate that magnetism can also be induced away from physical GNR zigzag edges through atomically precise engineering topological defects in its interior. A pair of substitutional boron atoms inserted in the carbon backbone breaks the conjugation of their topological bands and builds two spin-polarized boundary states around. The spin state was detected in electrical transport measurements through boron-substituted GNRs suspended between tip and sample of a scanning tunneling microscope. First-principle simulations find that boron pairs induce a spin 1, which is modified by tuning the spacing between pairs. Our results demonstrate a route to embed spin chains in GNRs, turning them basic elements of spintronic devices.
In spite of being a diamagnetic material, graphene can develop a special class of magnetism via the polarization of its π-electron cloud. Such π-paramagnetism is less localized than the more conventional d-or fmagnetism, and can interact over longer distances. Magnetic graphene nanostructures thus offer promising perspectives for a la carte engineering of interacting spin systems with applications in quantum spintronics devices [1][2][3][4]. The vision of graphene π-paramagnetism has been recently boosted by the development of onsurface synthesis (OSS) as a versatile bottom-up route. In OSS, nanoscale graphene flakes with customized shape and composition are fabricated over a metal substrate through the steered reactions between designed organic precursors [5,6]. Solid evidence of magnetism in flakes with segments of zigzag edges has been revealed in scanning tunneling spectroscopy experiments [7][8][9][10].
Substituting one carbon atom of the graphene lattice by heteratoms has been commonly identified as a potential route to induce magnetism [12,13]. A representative case is the doping of graphene with substitutional boron atoms ( Fig. 1), because it can be idealized as the removal of one electron from the conjugated bipartite lattice plus the energy upshift of a p z state. However, boron atoms do not induce any spin imbalance around, but simply behave as a point potential [12]. A prerequisite for the emergence of π-paramagnetism is that the point defect also causes a sufficiently large rupture of the conjugated electron system, for example by completely removing lattice sites or saturating p z orbitals [14][15][16], resulting in the localization of radical states. Here we show that inserting a pair of boron atoms in the carbon lattice of graphene nanoribbons (GNRs) enables a magnetic ground state. Density Functional Theory (DFT) simulations ( Fig. 1) show that, while magnetism is completely absent around a pair of such substitutional B atoms in different sublattices of extended graphene, in a 7-carbon-wide armchair GNR (7AGNR) the boron pair builds up a net magnetic moment of 2µ B (two Bohr magnetons). The spin polarization, shown in Fig. 1(a), decays towards the pristine segments with the characteristic shape of the 7AGNR end states [17] (see Supplementary Information (SI), [11] for a comparison). In fact, the spin cloud emerges from the rupture of the conjugated system imposed by the 2B-doped ring and the two neighbouring Clar sextets (green in Fig. 1(a)).  This moiety behaves as a highly reflective barrier for valence band electrons [18,19], thus inducing localized end states associated with the termination of the topological 7AGNR valence band [20]. This striking result offers the vision of combining band topology of nanoribbons [20][21][22] and heteroatoms for shaping spin textures in graphene ribbons.
We used the STM tip to pick individual 2B-7AGNRs from one end (cross in Fig. 2(b)) and lift them off to lie free-standing between tip and sample [7,29]. In this configuration, we measured the (two-terminal) transport through the suspended 2B-7AGNR as a function of tip retraction distance z to confirm the presence of a magnetic state. At the initial stages of suspension (2B unit still on the surface), the current passing through the ribbon showed a weak exponential decrease with retraction distance z (Fig. 2(c)), indicating a small tunneling decay constant through the elevated pristine GNR [29]. However, at a certain retraction length z p the current exhibited a pronounced peak, returning afterwards to the previous exponential decay. The current peak and its vertical position z p were reproduced for several retraction/approach cycles of the same ribbon, and appeared in all ribbons with a single 2B pair studied. In every case, the value of z p correlated with the distance between the 2B site and the contacted GNR-end (see SI [11]), thus indicating that the current peaks were caused by the detachment of a 2B moiety from the surface.
To find out the origin of the anomalous current increase, we measured current-voltage (I − V ) characteristics at z positions around z p . The differential conductance (dI/dV ) spectra (Fig. 3) show the sudden appearance of a narrow zero-bias dI/dV resonance at z p = 2.7 nm, which gradually decreases its amplitude with tip retraction, and disappears for distances above z > 2.9 nm. The resonance remained pinned at zero bias in all the z range observed, and its narrow linewidth reached a maximum value of Γ HWHM ≈ 6.0 ± 0.4 mV at z p ( Fig. 3(a)). The resonance vanished abruptly when the tip was approached below z p , but it was recovered by increasing z back above the z p = 2.7 nm onset. From its The GNR zigzag termination on the surface holds a spin-polarized radical state, absent at the contacted end due the bond formed with the tip's apex [7].
narrow line shape and fixed zero-bias alignment, we conclude that the resonance is a manifestation of the Kondo effect [31][32][33]. A Kondo-derived resonance appears in dI/dV spectra when a spin polarized state weakly interacts with the conduction electrons of an underlying metal [34]. Thus, its observation here shows that the 2B moieties abruptly develop a net spin around the point of their detachment from the surface.
To correlate these observations with the 2B-induced spin polarization predicted in Fig. 1, we performed DFT simulations of a finite 2B-7AGNR suspended between a model gold tip and the surface of a Au(111) slab (see Methods). Figure 4 shows the relaxed atomic structures of the GNR-junctions before, while, and after detachment of a 2B unit, and includes the computed constant spin density isosurfaces. Before 2B-detachment from the surface (z < z p , Fig. 4(a)), the intrinsic magnetism around the 2B units is quenched: the PDOS in the regions 1 and 2 around each boron atom is broad and spin unpolarized, contrasting with the clear spin polarization of free ribbons (shown as dashed plots). This is caused by the strong hybridization of the B atoms with the gold surface [19,24], which appear 0.6Å closer to the surface than the carbon backbone.
The detachment of the 2B moieties from the metal surface causes the emergence of a net spin polariza-tion, clearly reflected in their PDOS and spin density isosurfaces ( Fig. 4(b,c)). At the intermediate snapshot of Fig. 4(b), only one of the two B heteroatoms is detached from the surface, and the ribbon hosts a net spin S = 1/2 extending towards the free-standing segment (region 4). For the fully detached 2B case (Fig. 4(c)), both regions around the two B atoms (regions 5 and 6) are spin polarized, recovering the S = 1 state of the isolated 2B-7AGNR ( Fig. 1(b)). Since the Kondo effect is msot common for a spin 1/2 near a metal substrate, the most probable scenario is the one pictured in Fig. 4(b). There, the Kondo effect is caused by the spin 1/2 of region 4 interacting weakly with the surface when the first boron atom is detached, and vanishes when the second boron is cleaved ∼2Å higher, in consistency with the B-B distance. The experimental results are thus consistent with the spin polarization around 2B moieties in the free-standing 2B-7AGNRs.
Although the Kondo signal vanishes quickly with retraction, DFT finds that the S = 1 state of the free ribbon remains, and is clearly favored over an anti-parallel alignment by ∼ 14 meV per isolated 2B pair. The presence of a triplet state is striking; the two spin clouds at each side of the 2B center extend symmetrically over opposite sub-lattices of the 7AGNR, what usually favors an antiparallel kinetic exchange [16]. A detailed analysis re- veals that the hopping matrix elements between the two localized states at the sides of one 2B unit are very small (t intra ∼ 18 meV, see SI [11]). Consequently, the moiety formed by a 2B-doped ring surrounded by two Clar sextets is a very stable element that blocks conjugated electrons from hopping across. This explains the presence of a magnetic state because the borylated element acts as a barrier for valence band electrons of the 7AGNR segments [18], and induces spin polarized boundary states due to the non-trivial topology of this band [20]. Additionally, the 2B barrier also disconnects the boundary states at each side, and hinders the (anti-parallel) kinetic exchange between them. The stabilization of the triplet configuration is then the result of the weak direct overlap between both spin-polarized boundary states through the 2B barrier, which, due to the tiny hopping between them, dominates the exchange interaction and induces the ferromagnetic alignment of the spins according to Hund's rule.
The synthetic route to produce 2B-7AGNRs also allows incorporating multiple 2B units in one GNR to explore their interactions. We studied GNRs with two consecutive 2B moieties like in Fig. 5(a,b), spaced by 1.2 nm. Transport experiments through these GNR junctions as a function of tip-sample distance (Fig. 5(c)) also reveal deviations from an exponential decay with z, but now showing two peak features at retraction distances z p1 ≈ 1.60 nm and z p2 ≈ 3.05 nm. These val-ues are related to the positions of the 2B units (nominally ∼ 2.0 and ∼ 3.2 nm from the contact point, respectively). A map of (normalized) differential conductance as a function of bias and z (Fig. 5(d)) shows that both current features are also caused by narrow zerobias dI/dV resonances appearing at ranges z < 1.9 nm and 2.9 nm < z < 3.5 nm, respectively. Their line shape (Fig. 5(e)) is similar to the resonances observed for the single 2B case, and can also be attributed to Kondo states, which reflect the emergence of spin-polarization in the ribbon as each 2B unit is detached from the surface.
Although these results apparently suggest that each 2B behaves as an independent spin center, DFT simulations of free (2B) 2 -7AGNRs ( Fig. 5(b)) find that the two singly-occupied boundary states between neighboring 2B elements interact strongly and open a large hybridization gap [19], forcing them into a closed shell configuration. As a consequence, the spin polarization vanishes between two 2B sites, but persists outside this region as two uncompensated spin 1/2 clouds (Fig. 5(b)) with barely no preferred relative spin alignment. From our electronic structure calculations [11,19] we can characterize this hybridization by a relatively large effective hopping term t inter between boundary states of neighboring 2B units, which contrasts with the weak hopping t intra across each 2B unit. In fact, the electronic structure close to the Fermi level of a sequence of bory-lated units can be mapped onto the Su-Schrieffer-Heeger (SSH) model [35], characterized by two alternating hoppings along a 1D wire. Since t inter > t intra , an alternative way to understand the spin-polarized states in Fig. 5(b) is as zero-energy topological modes of a very short SSH chain. These simulations allows us to predict that a S=1 spin chain will emerge for larger inter-2B spacing, when both hopping terms become smaller than the Coulomb charging energy U of the boundary states [11].
We also simulated a (2B) 2 -7AGNR suspended between tip and sample to explore how these interactions survive in the experimental geometry. While detaching the first 2B site from the surface is equivalent to the single 2B case of Fig. 4, the spin polarization between the two 2B moieties is absent even when the second 2B remains partially hybridized with the surface (see Fig. 5(f)). Only when the last boron atom is detached, the terminal spin shows a 1/2 spin cloud outwards, being this the state responsible for the more extended Kondo effect observed in the experiment above z p2 . These results confirm the spin polarization predicted for boundary states appearing at the interface between (2B) 2 chains and pristine 7AGNR segments [20], which in essence are zero-energy modes of the 7AGNR valence band of similar nature than those created by a single 2B unit at every side.
The peculiar spin polarization of single 2B units and dimers is a remarkable consequence of the large and longrange exchange interactions present in GNRs [7,20]. Our results indicate that the spins should survive in freestanding GNRs for sufficiently low concentrations of 2B units forming either a chain of S = 1 spins, or for sufficiently dense doping S = 1/2 spins at the boundary between pristine and borylated segments. Tuning the separation between 2B units appears to be a promising strategy to control spin polarization through a change in the correspondence between inter-2B and intra-2B interactions, offering the perspective of engineering spin chains of different kinds. Furthermore, the stronger sensitivity of substitutional boron heteroatoms to chemical bonding endows these systems with ideal properties to manipulate complex spin states in chains. References 8

Sample preparation
To prepare the Au(111) substrate we perform a cleaning cycle by sputtering the crystal with Neon ions for 10 minutes and afterwards annealing the crystal in UHV conditions (p < 10 −9 mbar) at T = 740 K for 15 minutes. We prepare our samples by simultaneously evaporating DBBA and B 2 -DBTA molecular precursors on the clean Au(111) surface [1][2][3]. Afterwards an Ullmann reaction and cyclodhydrogenation are induced by successive annealing the sample to 200 • C for few minutes and flashing to 400 • C for 30 seconds.
The samples are prepared in situ and transferred directly into a homebuild STM kept with liquid Helium at ≈ 6 K. We use a PtIr tip with an atomically sharp gold termination. The tip is prepared by controlled indention into the clean Au(111) surface.

Lifting procedure
The lifting is performed as following. After stabilization of the tip above the clean Au(111) surface with a fixed set of parameters (I = 30 pA, V b = −1000 mV), we open the feedback loop of the control electronics. Afterwards, we position the tip above the termination of the GNR and approach towards the surface with a bias V b = 25 mV. At some point during the approach a sudden increase in the current is observed. This corresponds to the formation of a tip-molecule bond and defines z = 0 during each experiment. After the tip-molecule bond is formed we are able to partially lift the GNR from the surface.
Custom cleaving trajectories are employed for the different molecules, depending on the strength of the tip-molecule bond. In the most robust cases we cleave the molecule by simple retraction along the z-axis. In less strongly bound cases, we need to move additionally in the xy-plane, following a trajectory along the direction of the GNR-backbone on the surface. However, we find that variations in the lifting trajectory do not influence at which z 1 we observe the peak in I and the Kondo-resonance in the differential conductance spectra.

Normalization of differential conductance spectra
Each individual spectrum taken at one z 0 was normalized following G norm (V, z 0 ) = , where G(V, z 0 ) is the measured dI/dV -signal at z 0 and G min,z 0 (G max,z 0 ) is the minimal (maximal) value of G(V, z 0 ) in the analysed bias interval. The normalization is only applied to the data presented in Fig. 4D.

Details on ab-initio calculations
First-principles electronic structure calculations are performed using density functional theory as implemented in the SIESTA package [4,5], where the valence-electron wave functions are expanded using a linear combination of numerical atomic-orbitals as a basis set and the core electrons are replaced by norm-conserving Troullier-Martins pseudopotentials [6]. A double-ζ plus polarization (DZP) basis set is adopted for the surface Au atoms, where an extended basis is considered for the top atomic layer optimized for the description of the Au(111) surface [7]. The tip model is defined by a gold rod with its axis along the Au(100) crystalline direction and its valence-electrons described using a DZP basis set for s-orbitals and single-ζ non-polarized (SZ) basis for d-orbitals (i.e. DZP-SZd). Such geometry and basis set has been shown to provide a reasonable description of the contact, with a smooth local density of states at the tip apex around the Fermi level [8]. For the atoms defining the borylated ribbons a double-ζ non-polarized (DZ) basis set is employed, where the orbital radii are defined using a 30 meV energy shift [5]. Dispersive interactions between the borylated ribbon and the metallic surface are considered using the van der Waals density functional by Dion et al. [9] with the modified exchange by Klimeš, Bowler and Michaelides [10]. The real-space grid for integrations is defined using a 300 Ry energy cutoff [5]. The smearing of the electronic occupations is defined by an electronic temperature of 300 K with a Fermi-Dirac distribution. The self-consistency cycle is stopped when variations on the elements of the density matrix are lower than 10 −5 as well as 10 −5 eV for the Hamiltonian matrix elements.
The freestanding 2B-GNRs, as the one shown in Fig. 2A in the main text, are computed within a simulation cell that includes 20 and 40Å of empty space in the directions along and perpendicular to the ribbon, respectively, to avoid interactions with the periodic replicas. Geometry optimizations are performed using the conjugate gradient (CG) method until all forces are lower than 20 meV/Å. For the calculations involving the suspended borylated GNRs, the Au(111) surface is represented by a 3-layer thick slab within a simulation cell of dimensions 46.16 × 24.98Å 2 along the periodic directions. A H passivation is employed at the reverse side of the slab in order to prevent spurious effects due to interaction between surface states belonging to the top and bottom surfaces of the slab [11]. A 1×2 k-point mesh is used to sample the bidimensional Brillouin zone. The slab defined in this way comprises 640 atoms (7680 orbitals), where the atoms from top Au layer are allowed to relax with the CG method until forces are lower than 20 meV/Å (the H passivation and the two bottom Au layers are kept fixed).

Lift of GNRs with 2B sites at different positions
In Fig. 1 we present two further examples of I-z retraction curves of two 2B-GNRs with the 2B site at 2.0 nm and 3.6 nm from the contact point, respectively. The peaks in the current are observed at z p = 1.2 nm and z p = 2.4 nm, in agreement with the different 2B to contact point distances. In both cases the peak in the current is caused by a zero-bias Kondo resonance, like in the case presented in the main text.

Density Functional Theory simulations of the edge states and boron-induced inner states in 7AGNR
The boundary states induced by the substitutional boron pairs exhibit the same shape, symmetry, and decay towards the pristine segments as the 7AGNRs end states, localized in the zigzag terminations. Fig. 2 below presents a comparison between two 7AGNRs end states facing each other (shown in the top) and the 2B-induced bonding and anti-bonding states (bottom). Notice that, interestingly, for 2B centers the "bonding" state is higher in energy than the "anti-bonding" linear combination.

Minimal model of an isolated 2B center in 7AGNR
A single 2B moiety in a 7-armchair graphene nanoribbon (7AGNR) is found to build up a total spin S=1 due to the predominant ferromagnetic interaction between the two singly occupied boundary states (B-states) localized at each side. The purpose of this section is to formulate the minimal model (depicted in Fig. 3) capable of explaining the low-energy electronic structure and the magnetism of isolated 2B pairs.
The effective Hamiltonian describing an isolated 2B pair also can be seen in Fig. 3 in the limit U t intra . Exact and mean field solutions of this model are readily obtained and can be compared with the relative energies and the band structures of different DFT solutions. In particular, we use information from DFT simulations (both non-spin-polarized and spinpolarized with either ferromagnetic or anti-ferromagnetic spin alignment) to estimate the following values of the effective parameters describing the system according to the model of Fig. 3: the charging energy of each boundary state U∼260 meV, the effective hopping matrix element between both boundary states t intra ∼18 meV, and the Coulomb exchange interaction J intra ∼17 meV.
While the Coulomb exchange J intra favours parallel alignment, the kinetic exchange (KE) term, given by 4( t 2 intra U ) (Fig. 3) tends to align spins antiferromagnetically. The small value of t intra relative to U translates into a small KE interaction, of the order of ∼1.2 meV, i.e., around 14 times smaller than the Coulomb exchange within each 2B pair. Therefore, the ground state of an isolated 2B pair corresponds to a total spin S=1 as expected from Hund's rules.

Electronic and magnetic interactions between neighbouring 2B centers
We now consider the case of periodic doping of 7AGNRs with 2B centers. The distance between two 2B sites is given by the number n of 7AGNR unit cells defining the periodicity of the system. Fig. 4 compares the computed hopping matrix elements between B-states at each side of a given 2B pair (t intra ) and those between two neighbour 2B pairs (t inter ). These values are obtained by fitting to the Su-Schrieffer-Heeger (SSH) model [12] the B-derived in-gap states appearing in the non-spin-polarized DFT band structures in Fig. 5 of Ref. [3] for 2B-(7AGNR) n with n = 3, . . . , 13. The fitting is rather accurate for n > 3, while for n=3 there is a somewhat larger uncertainty due to the interaction between the lower B-derived state and the valence band of the 7AGNR near the Γ-point. As mentioned above we obtain a very small effective hopping parameter between the boundary states across the diboron moiety t intra , which amounts to 18 meV for n ≥ 5.
The effective electron hopping between proximal B states t inter is much larger than t intra for n < 7, and decreases exponentially with the separation, resulting residual above n ≥ 7. The spacing between diboron sites in the results of Fig. 3 of the manuscript corresponds to n=3. The inter-diboron hopping is very large in this case (∼350 meV > U), resulting in a closed-shell state with no spin polarization in the spacing between them. The small electron hopping across a 2B site protects the outer B states, which maintain each a net spin polarization of S=1/2, as detected in Fig. 3.
The larger hopping between neighbour diboron sites than inside one, i.e., t inter > t intra , resemble a non-trivial bulk-boundary correspondence, as described by the SSH model. According to this model, for the smallest spacings between diborons (n=3, 5) boundary states are expected at the terminations since t intra < t inter . However, if we consider Coulomb correlations, only the case n=3 shows end states because it fulfills U < t inter . The n=5 case, for which already U > t inter , results in a chain with a spin texture formed by S=1 spins localized at the 2B sites and aligned antiferromagnetically along the chain due to the large kinetic exchange associated with t inter . Electronic interactions between B-derived states from DFT calculations: (a) Schematic model of a 2B-AGNR with intra-2B and inter-2B hopping terms indicated. J intra represents the direct Coulomb exchange resulting in ferromagnetic alignment of the spins within each pair. (b) Computed values of intra-2B and inter-2B hopping terms from the analysis of nonspin-polarized DFT band structures of the periodic systems [3]. The inter-2B spacing given by n, the number of 7AGNR unit cells between 2B sites. The correspondence with the Su-Schrieffer-Heeger (SSH) model is evident. The inset shows the exponential decay of t inter and magnifies the crossing point where the topology of the non-spin-polarized band structure of the system changes. However, as discussed in the text, the non-polarized solution is not relevant in the limit where t inter < U that clearly holds for n ≥ 5.