Research ArticleChemistry

Multiple heteroatom substitution to graphene nanoribbon

See allHide authors and affiliations

Science Advances  13 Apr 2018:
Vol. 4, no. 4, eaar7181
DOI: 10.1126/sciadv.aar7181


Substituting heteroatoms into nanostructured graphene elements, such as graphene nanoribbons, offers the possibility for atomic engineering of electronic properties. To characterize these substitutions, functionalized atomic force microscopy (AFM)—a tool to directly resolve chemical structures—is one of the most promising tools, yet the chemical analysis of heteroatoms has been rarely performed. We synthesized multiple heteroatom-substituted graphene nanoribbons and showed that AFM can directly resolve elemental differences and can be correlated to the van der Waals radii, as well as the modulated local electron density caused by the substitution. This elemental-sensitive measurement takes an important step in the analysis of functionalized two-dimensional carbon materials.


Among two-dimensional (2D) materials, since the first systematic mechanical exfoliation from graphite (1), graphene has been one of the most attractive materials due to its unique electronic and mechanical properties. To realize applications such as electronics, chemical sensing, photocatalysis, and battery electrodes, ideal graphene itself is, however, insufficient because of the absence of a band gap. Structural and/or electronic functionalization is a key concept in controlling the properties, and many efforts have been directed to, for example, fabricating ribbon-like structures (2), making holes (3), and doping with other elements (4). On-surface chemical reaction with a designed molecule (5) has proven a promising technique to synthesize structurally controlled graphene, forming so-called graphene nanoribbons (GNRs) (6). In these studies, halo-substituted molecules are deposited on a metal surface and consequently heated to conjugate each other. After cyclodehydrogenation at higher temperature, GNRs can be synthesized, with the width and edge structure controlled by the molecular backbone (7). Furthermore, if a heteroatom is incorporated into the precursor, doped GNRs can be synthesized, allowing for the position and density of doped elements to be controlled. So far, substituting a single heteroatom [nitrogen (8, 9), boron (10, 11), sulfur (12), oxygen (13), and so on] into a GNR has been demonstrated.

The development of functionalized tips in atomic force microscopy (AFM) has proven to be a key stage in the evolution of the technique (1416), with the CO tip offering the ability to systematically resolve the inner structures of molecules. This high-resolution AFM is particularly beneficial for the investigation of single and assembled molecules on surfaces (1723), and even the steps of surface chemical reactions (24, 25) or tip-induced reactions (2628). Furthermore, it is sensitive enough to resolve small differences in the bond orders in polycyclic aromatic hydrocarbons (PAHs) (29) and can be also used for atomic force spectroscopy (3032). However, the use of AFM in the chemical analysis of atoms embedded in PAHs is in its early stage, despite its importance in controlling the chemical and electronic properties. For AFM on graphene, nitrogen-doping (3335) and boron-doping (36) have been investigated. However, it has been shown that AFM is also quite sensitive to environmental conditions such as moiré patterns (37). For instance, our previous work on boron-doping into GNR showed that the contrast of incorporated boron atoms mainly relates to the adsorption height difference (30 pm) rather than the elemental difference (10). Therefore, synthesizing atomically flat but multiple heteroatom-doped PAHs on surfaces is of particular importance to provide an ideal playground for this topic, alongside the intrinsic interest in graphene band gap engineering.

Here, we demonstrate the discrimination of boron, carbon, and nitrogen atoms embedded at the center of an n = 7 GNR with a CO-functionalized tip of AFM at low temperature. The multiple heteroatom-doped GNR was in situ synthesized by on-surface chemical reactions, and its structure and electronic properties were investigated by AFM and scanning tunneling microscopy (STM), as well as density functional theory (DFT) calculations. The combined study shows that each element can be resolved, and the attractive interaction force increases in the order of nitrogen, carbon, and boron. This can be correlated with the radii of each species and the relative onset of repulsive interactions rather than the difference in charge states.


To embed heteroatoms at equivalent atomic sites in a nanographene structure, we used an approach using on-surface chemical reactions with designed precursor molecules (5). Similar to previous work on periodic boron-doping into GNRs (10, 11), we used 5-(10-chloroanthracen-9-yl)-10-(10-iodoanthracen-9-yl)-5,10-dihydrodibenzo[b,e][1,4]azaborinine, which is composed of three anthracene units conjugated with single bonds (Fig. 1A). The center unit has a C4BN cyclic ring as the boron and nitrogen atoms locate closer to the adjacent anthracene units following the iodine and chlorine substitutions, respectively. The precursor was deposited onto a Au(111) surface and subsequently was annealed at 200°C to induce the conjugation and finally at 400°C to induce the cyclodehydrogenation (6). In this manner, an n = 7 BN-GNR can be synthesized. Note that because of the asymmetric structure of the precursor molecule, two different reaction paths can be expected (Fig. 1, B and C).

Fig. 1 Synthesis of multiple heteroatom-substituted GNR.

(A) Schematic drawing of 5-(10-chloroanthracen-9-yl)-10-(10-iodoanthracen-9-yl)-5,10-dihydrodibenzo[b,e][1,4]azaborinine. (B) On-surface reaction paths of long- and (C) short-periodicity boron-nitrogen GNRs (BN-GNR). Gray and purple triangles indicate boron- and nitrogen-centered membered rings, respectively. (D) STM topography of the synthesized BN-GNR on Au(111). (E and F) Corresponding AFM images. Measurement parameters: Vtip = −200 mV and I = 0.8 pA in (D). Vtip = 0 mV and A = 60 pm in (E) and (F).

Following the synthesis process, an STM topography shows the BN-GNRs synthesized on Au(111) (Fig. 1D). The length of BN-GNR was usually shorter than that of the pristine GNR (6), as well as B-GNR (10). This may relate to the fact that the large dimeric precursor, terminated with Cl atoms as the major product (Fig. 1B), requires more energy to diffuse on the surface. Consequently, the probability to cause defects, such as termination by hydrogen or cleavage of C–N or C–B bonds, presumably becomes higher (figs. S1 and S2). Nevertheless, we can synthesize a longer BN-GNR with a lower doping concentration by codepositing 10,10′-dibromo-9,9′-bianthryl (Sigma-Aldrich) and 5-(10-chloroanthracen-9-yl)-10-(10-iodoanthracen-9-yl)-5,10-dihydrodibenzo[b,e][1,4]azaborinine (fig. S3). Figure 1 (E and F) shows AFM images in which the inner structures of the BN-GNRs are resolved. A dark triangular feature is most significantly pronounced and would be in good agreement with the previous measurement on B-GNR (10) if one takes account of the fact that only one boron exists at the center of the doped anthracene instead of two (Fig. 1A). Thus, one might readily conclude that the boron atom is at the center of the dark triangular feature; however, as we show in the following, this is incorrect.

To verify this assignment of the dopants, we performed DFT calculations combined with simulated AFM. In the previous study of GNR (10), the source of the dark contrast seen on the B-B dopant was due to a displacement of the borons toward the gold surface, and the chemical and electronic properties of boron had a minor influence. In the case of the BN-GNR, the B-N dopant does not show this displacement and the ribbon appears flat (relative adsorption height to carbons, B = −2 pm and N = +4 pm; Fig. 2A). Therefore, the elemental character is dominant in the AFM contrast—it is actually the nitrogen site that is seen as dark (Fig. 2B and fig. S4), and the adjacent bright site corresponds to the boron atom (indicated by blue arrows in Fig. 1, E and F). The different orientations of B-N dopant relate to the different reaction paths (Fig. 1, B and C).

Fig. 2 AFM results of BN-GNR.

(A) Simulated atomic structure from above (left) and to the side (right), with carbon atoms shown in cyan, gold in yellow, boron in blue, nitrogen in green, and hydrogen in white. (B) Simulated AFM images at decreasing tip height. (C to J) Series of high-resolution AFM images taken at different Z distances. (K) Force spectroscopic measurements taken at different elements in the BN-GNR. (L) Close view. (M) Calculated frequency change (Δf) over labeled atomic sites as a function of tip-surface distance (note that this is referenced to the fixed part of the tip, and the effective distance is about 4.0 Å smaller once the CO molecule is included). Measurement parameters: Vtip = 0 mV and A = 60 pm.

Figure 2 (C to J) shows a series of AFM images taken at different Z heights. At the largest separation, a weak and bright contrast appears at the boron atom site, as indicated by blue arrows. By setting the tip closer to the BN-GNR, six-membered carbon rings appear (Fig. 2E). Because no significant C–N bond is observed at the distance, the dark triangle is pronounced in the image. At Z = 32 pm, three C–N bonds also appear, yet the tip-sample distance for the other sites is already closer than the optimal (Fig. 2G). Thus, the distortion of the C6 rings is observed particularly at the armchair edge, relating to the asymmetric CO tip deflection caused by the lateral force field. At a smaller separation (Z = 10 pm), the centers of the C5B ring are filled with the bright contrast, as indicated by red arrows, due to the strong repulsive interactions (Fig. 2I) (29). Similar features are also observed at the C6 and C5N rings at closest distance—note that this rich transition of the contrast occurs in a small Z distance range (~100 pm). The features seen in AFM images as a function of the Z distance are reproduced in the simulation images shown in Fig. 2B. In particular, at close approach, the appearance of C–N bonds is weakly seen, and the filling of the C5B rings and the apparent bond distortions are visible. Taking into account that CO-tip AFM is dominated by short-range repulsive interactions, the nature of the contrast in the images shows that the onset of short-range repulsion occurs first at B sites, then C, and finally N.

To quantitatively investigate the tip-sample interaction, we conducted force spectroscopic measurements at the different element sites, but all locating at the centers of the flat BN-GNR (Fig. 2K). The insets show the high-resolution AFM image and the corresponding chemical structure. Note that the long-range nonsite-dependent macroscopic interaction between tip and substrate was removed by subtracting with a curve taken on a bare Au(111), yet the van der Waals interaction between tip and GNR remains in the presented curve. Furthermore, no significant site-dependent long-range van der Waals contribution is seen because the tails of frequency curves merge perfectly. Nevertheless, because all elements were substituted at equivalent sites in the symmetric structure of n = 7 GNR, detected differences readily relate to the elemental contributions. No significant difference is seen in the tails of the attractive force curve (above 200 pm), which is rarely seen in the covalent bonding (38, 39), and atomic pairwise van der Waals interactions (31). Figure 2L shows the close view of the frequency shift curves. Before the turning point, a small difference can be seen (the boron atom has a lesser negative frequency shift), which corresponds to the contrast development in Fig. 2C. In this region where the repulsive interaction starts affecting the total interaction, clear differences are seen in the corresponding AFM images (Fig. 2E). We fitted the turning points with a polygonal curve and found that the nitrogen site has the most negative frequency shift (−2.49 Hz), followed in order by carbon (−2.14 Hz) and boron sites (−1.90 Hz). Note that the value for the carbon sites is a mean of three points. The simulated curves, calculated using the DFT structure and Hartree potential (see the Supplementary Materials for details), show a qualitatively similar behavior with the nitrogen site significantly deeper and boron and carbon sites closer to each other (Fig. 2M; the difference, of about a factor of 2, in absolute magnitude of frequency shift seen in simulations reflects that we do not know the experimental atomic tip structure exactly but do not need to for interpreting the results). Despite that the differences of the extracted attractive forces were rather small (<10 pN) (fig. S5), AFM can resolve the difference of B, C, and N in the PAHs.

The CO-terminated tip is known to have a complex charge state, arising from the sum of Smoluchowski effect of the metal tip and the electrostatic field of the CO molecule attached (40), with each contribution having their own decay length. For a neutral or weakly charged tip, as used here, for a separation larger than the gap between boron and nitrogen, no significant contrast would be developed because the interaction with negatively and positively charged atoms balances. However, a naïve interpretation would expect that atomic-scale contrast should appear if the negatively charged CO tip approaches to the sites close enough, distinguishing the substitutional atoms by their charge. At a glance, these mechanisms would be captured in the experimental and theoretical calculations (Fig. 2, B and E), yet the calculated charge analysis gives the opposite trend—B and N are positively and negatively charged, respectively (fig. S6). Therefore, the local charge state does not reflect to the AFM image contrast once the heteroatoms substitute to PAHs, and the short-range repulsive contributions dominate—increasing the charge on the tip soon results in images that deviate significantly from experiment. In general, this can simply be interpreted in terms of the van der Waals radius of the elements (B, 192 pm; C, 170 pm; and N, 155 pm) (41), such that repulsive interactions dominate at larger distances over the larger B atom and closer on the smaller C and N. Therefore, once heteroatoms are embedded in the graphitic material, the out-of-plane interaction is dominated by the van der Waals interaction. However, the coincidence in the tails of the attractive force cannot be represented only with the difference of the radii, and the effective atomic C6 coefficients are also affected (fig. S7), which is also seen in the charge analysis (fig. S6). In this manner, AFM with a CO tip can resolve the differences of the radii of heteroatoms if they are substituted into flat PAHs on surfaces, and no significant local charge transfer effect occurs.

Although the variation of the repulsive forces is rather small, significant differences among B–C, C–C, and C–N bond lengths in the AFM image are seen. Boron, nitrogen, and carbon are labeled, as shown in Fig. 3A. Here, we measured the apparent bond lengths at the center, as summarized in Table 1. We found that the C–N bond is the longest (196 pm), followed in order by C–C (142 pm) and B–C (135 pm). However, the actual bond length obtained by the DFT calculation has a different trend. This inconsistency can be rationalized by the fact that AFM with a CO tip is more sensitive to the bond status (order) (29). The N–C bond has a lower electron density than C–B, whereas the N atom has a greater electron density than the B atom (Fig. 3B). This complex of the local electron density relates to the donation and acceptance of the electron by the B and N substitution to PAH. In this manner, a significant elemental contribution appears as the resultant distortion of the six-membered ring in PAHs. This analysis is valid even if the heteroatom-substituted GNR is no longer flat as in B-GNR (10).

Fig. 3 Apparent bond analysis.

(A) Analysis of the bond lengths in the AFM image. (B) Calculated valence electron density (at an isocontour of 1.0 eÅ−3).

Table 1 Summary of apparent bond lengths.
View this table:

We finally investigated the electronic structure of BN-GNR by STM and scanning tunneling spectroscopy (STS). Figure 4A shows the STM topography taken near the Fermi level (−2 mV), in which the corrugation at the center is seen. The corresponding AFM image indicates that the bright and dark spots in the STM topography relate to the boron and nitrogen sites, respectively (Fig. 4B). dI/dV spectra were taken with a clean Au tip above six different sites, as well as above clean Au(111) site for reference (Fig. 4, C and D). Below the Fermi energy, the contribution originating from the Shockley surface states of Au(111) is seen in all curves, relating to the coupling to the substrate. This modulation obscures the clear valence band edge (VBE), yet the corresponding constant height dI/dV maps indicate that the VBE is around −0.3 eV below the Fermi energy (Fig. 4, E to L). By contrast, a clear conduction band edge is seen at 2.0 eV. These values are upshifted from the pristine GNR by 0.4 eV, although the band gap of 2.3 eV remains the same, which is related to the fact that the substituted atoms locate at the center of n = 7 BN-GNR (42). Furthermore, several peaks relating to the nitrogen donor and boron acceptor band energy levels can be seen in the band gap (1.2 and 1.4 eV)—seen also in the simulated density of states (DOS) (Fig. 4N). Our finding will aid in synthesizing and analyzing functionalized ultrathin carbon materials with AFM-based on-surface chemical reaction and will motivate the fabrication of forthcoming GNR-based devices.

Fig. 4 Electronic structure of BN-GNR.

(A) STM topography and (B) the corresponding AFM image. STS taken at the center in (C) and the edge in (D). CBE, conduction band edge. Six different positions were measured as labeled in (A). (E to L) Series of constant height dI/dV maps measured with a lock-in amplifier (root mean square amplitude = 14 mV and frequency = 512 Hz) at the different bias voltages. Arrows indicate the positions of the STS measurements. (M) Simulated STM topographies at a bias of −1.0 V (left) and +1.5 V (right). (N) Calculated DOS.


In summary, we demonstrated boron- and nitrogen-substituted GNRs by on-surface chemical reaction. Both elements were substituted at the center of n = 7 GNRs. In addition to the synthesis and investigation of the developed boron and nitrogen band structures, we present the capability of AFM with a CO-functionalized tip to directly resolve the elemental difference in the graphene. The combination with our theoretical calculation has allowed us to discriminate the elements embedded in large PAHs via the van der Waals radii, as well as the local electron density caused by the substitution.



All measurements were performed with a commercially available Omicron low-temperature STM/AFM system, operating in ultrahigh vacuum at 4.8 K. We used a tuning fork with a chemically etched tungsten tip as a force sensor (43). The resonance frequency and the mechanical quality factor are 24,803.5 Hz and 18,648, respectively. The high stiffness of 1800 N·m−1 realizes a stable operation with small amplitude of 60 pm (44). The frequency shift, caused by the tip-sample interaction, was detected with a commercially available digital phase-locked loop (OC-4, Nanonis; and HF2-LI and HF2-PLL, Zurich Instruments) (45). For the STM measurements, the bias voltage was applied to the tip while the sample was being electronically grounded. The tip apex was ex situ sharpened by milling with focus ion beam. The tip radius was less than 10 nm. A clean gold tip was in situ formed by indenting to the Au sample surface and applying a pulse bias voltage between tip and sample several times. For AFM, the tip apex was terminated with a CO molecule, which was picked up from the surface (46). Clean Au(111) surfaces were in situ prepared by repeated cycles of standard Ar+ sputtering (4 × 10−6 mbar, 1000 eV, and 15 min) and annealing at 730 K. In this experiment, we used 5-(10-chloroanthracen-9-yl)-10-(10-iodoanthracen-9-yl)-5,10-dihydrodibenzo[b,e][1,4]azaborinine that was deposited on the surfaces from crucibles of Knudsen cell by heating at 220°C. The temperature of the substrate was kept at room temperature. 2D frequency shift mapping was measured by a series of Z distance measurement of the frequency shift. Measured images were analyzed using the WSxM software (47).

Theoretical calculations

All first-principles calculations in this work were performed using the periodic plane-wave basis Vienna Ab initio Simulation Package code (48, 49) implementing the spin-polarized DFT. To accurately include van der Waals interactions in this system, we used the optB86B+vdW-DF functional (5052), selected on the basis of previous work showing that it provides a sufficiently accurate description for all subsystems involved in the measurement. We also calculated adsorbed molecular structure using vdW functionals D3 (53) and TS (54) but observed no significant differences in the structure compared to optB86B. Projected augmented wave potentials were used to describe the core electrons (55), with a kinetic energy cutoff of 550 eV (with PREC = accurate). Systematic k-point convergence was checked for all systems, with sampling chosen according to system size and a mesh of 3 × 1 × 1 being used for the final production runs. This approach converged the total energy of all the systems to the order of millielectron volt. The properties of the bulk and surface of Au, graphene, and the isolated molecular structure were carefully checked within this methodology, and excellent agreement was achieved with experiments.

For calculations of the infinite doped ribbon, a surface slab of 9 × 6 × 5 in terms of the Au unit cell was used, with a vacuum gap of at least 1.5 nm and the upper three layers of Au, and the ribbon was allowed to relax. All adsorption energies were calculated by subtracting the individual components of the system, in the same unit cell, from the total energy of the final system. Bader charge analysis was used to estimate charge transfer in the simulations (56).

For calculated AFM images, we used our implementation of the model developed by Hapala et al. (57, 58). The ribbon structure was taken from DFT simulations, and the electrostatic potential was extracted from the Hartree potential (59). We also used a quadropole dz2 tip model (60), with an effective charge of −0.2 and a tip lateral spring constant of about 0.25 N·m−1, which proved to give the best agreement with experiment. All other parameters are the same as intended by Hapala et al., and the simulated AFM scan was performed at a resolution of 5 pm (in all directions), with a force tolerance criterion of 4 × 10−6 eV·Å−1. The 3D force field was subsequently converted into a frequency shift image using the experimental parameters.

Synthesis of precursor molecule

5-(10-chloroanthracen-9-yl)-10-(10-iodoanthracen-9-yl)-5,10-dihydrodibenzo[b,e][1,4]azaborinine was synthesized from compounds N-(2-bromophenyl)-N-phenyl-10-chloroanthracen-9-amine 2 (0.459 g, 1.0 mmol) and 9,10-dibromoanthracene 3 (1.01 g, 3.0 mmol) in 8.6% yield (59.3 mg, a pale yellow solid). The synthetic details are described in the Supplementary Materials. Infrared (neat): 3051 (Ar-H), 1599, 1564, 1543, 1479, 1468, 1437, 1408, 1344, 1254, 1213, 1169, 1144, 1040, 1024, 939, 926, 899, 746, 667, and 646 cm−1; mp: >300°C; 1H NMR (nuclear magnetic resonance) (CS2/CDCl3 = 2/1, 500 MHz) δ 6.40 (d, J = 9.2 Hz, 2H, e), 6.90 (t, J = 7.2 Hz, 2H, g), 7.21 to 7.25 (m, 4H, f, j), 7.42 (dd, J = 1.2, 8.9 Hz, 2H, d), 7.45 (ddd, J = 1.2, 6.0, 8.9 Hz, 2H, c), 7.51 (ddd, J = 1.2, 6.3, 9.2 Hz, 2H, k), 7.59 (dd, J = 1.7, 7.2 Hz, 2H, h), 7.62 (d, J = 1.2, 8.6 Hz, 2H, i), 7.66 (ddd, J = 1.7, 6.0, 8.9 Hz, 2H, b), 8.55 (d, J = 9.2 Hz, 2H, l), and 8.71 (d, J =8.6 Hz, 2H, a); 13C NMR (CS2/CDCl3 = 2/1, 126 MHz) δ 117.0 (2C), 120.7 (2C), 123.6 (2C), 124.2 (2C), 125.7 (2C), 127.0 (2C), 127.2 (2C), 127.4 (2C), 127.9 (2C), 129.7 (2C), 129.9 (2C), 131.2 (2C), 131.9, 133.4 (2C), 133.7 (2C), 133.8, 133.9 (2C), 135.1, 137.7 (2C), and 146.4 (2C). The NMR signal of the carbon α to the boron was not observed. 11B NMR (CS2/CDCl3 = 2/1, 160 MHz) δ 55.0; high-resolution mass spectroscopy (direct analysis in real time) mass/charge ratio [M + H]+ calculated for C40H25B1I1Cl1N1 691.0742; observed 691.0742.


Supplementary material for this article is available at

Supplementary Materials and Methods

fig. S1. Examples of various couplings.

fig. S2. Small boron-nitrogen codoped GNR.

fig. S3. Boron-nitrogen codoped GNR with a low-doped density.

fig. S4. Simulated BN-GNR structures, spectra, and electronic parameters.

fig. S5. Force extraction.

fig. S6. Simulated charge analysis.

fig. S7. Simple model calculation analysis based on Lennard-Jones potential as E(Z) = 4ε[ (σ/Z)12 − (σ/Z)6].

fig. S8. Synthesis of BN-GNR precursor.

fig. S9. Chemical structure of 9-bromo-10-chloroanthracene.

fig. S10. Chemical structure of N-(2-bromophenyl)-N-phenyl-10-chloroanthracen-9-amine (2).

fig. S11. 1H NMR spectrum of 2 in CDCl3 at 25°C.

fig. S12. 13C NMR spectrum of 2 in CDCl3 at 25°C.

fig. S13. Chemical structure of 5-(10-chloroanthracen-9-yl)-10-(10-iodoanthracen-9-yl)-5,10-dihydrodibenzo[b,e][1,4]azaborinine (1).

fig. S14. 1H NMR spectrum of 1 in CS2/CDCl3 (=2/1) at 45°C.

fig. S15. 13C NMR spectrum of 1 in CS2/CDCl3 (=2/1) at 45°C.

fig. S16. 11B NMR spectrum of 1 in CS2/CDCl3 (=2/1) at 45°C.

References (6163)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: J.T. and A.S.F. acknowledge use of the CSC, Helsinki for computational resources. Funding: This work was supported, in part, by Japan Society for the Promotion of Science KAKENHI grant number 15K21765, the Japan Science and Technology Agency (JST) Precursory Research for Embryonic Science and Technology for a project of “Molecular technology and creation of new functions,” a Grant-in-Aid for Scientific Research on Innovative Areas “π-System Figuration” (17H05164), the Swiss National Science Foundation, and the Swiss Nanoscience Institute. A.S.F. has been supported by the Academy of Finland through its Centres of Excellence Program (project no. 915804) and EU project PAMS (contract no. 610446). Author contributions: S.K. conducted experiments and carried out data analysis. S.N. and T.H. designed and synthesized the precursor molecule. J.T. and A.S.F. conducted theoretical calculations. S.K. and A.S.F. contributed to writing the manuscript. All authors commented on the manuscript. S.K. and T.H. planned the project. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article