R Zoubkoff

F Triozon

Y.-M Niquet

S Latil

Boron and nitrogen codoping effect on transport properties of carbon nanotubes

Keywords: number(s): 61.48.De, 73.22.-f, 73.23.-b, 31.15.aq

This paper reports a theoretical study of the effect of boron and nitrogen codoping on the transport properties of carbon nanotubes (CNTs) at the mesoscopic scale. A new tight-binding parametrization has been set up, based on density functional theory calculations, that enables a reliable description of the electronic structure of realistic BN-doped CNTs. With this model, we have carried out a deep analysis of the electronic mean free path (MFP) exhibited by these nanostructures. The MFP is highly sensitive to the geometry of the scattering centers. We report that the relative distance between B and N atoms in the network influences drastically the electronic conduction. Moreover, we point out that the scattering induced by small hexagonal BN domains in the carbon network is less important than the BN-pair case.

I. INTRODUCTION

Single-wall carbon nanotubes (CNTs) have attracted an intense research effort from both fundamental and applied aspects. They exhibit very promising electronic properties that are directly related to their specific quasi-one-dimensional geometry. In order to tune their electronic structure and properties, B and N chemical doping can be used. The substitutional doping of B and/or N onto the hexagonal network has been intensively studied [START_REF] Terrones | Carbon Nanotubes[END_REF].

Beyond doping, the direct synthesis of C-BN heterostructures is highly desirable, because nanotubes made from hexagonal boron nitride (h-BN) are wide gap insulators [2]. We expect that C-BN heterojunctions can act as basic nanodevices, since the electronic properties of selected regions can be tailored to act as n or p semiconductors. Many applications can be addressed in nanoelectronics, nonlinear optics, or field emission sources. Doped carbon CNTs can be produced by adding moieties containing B or N elements during the growth process [3,4]. They can also be produced by a high-temperature chemical posttreatment that includes heteroatoms onto the carbon sheet [5].

The chemical nature and the specific geometry induced by the substitutional sites have been probed by various techniques: electronic energy loss spectroscopy, x-ray photoemission spectroscopy, and scanning tunnelling microscopy. More recently, nanotube C-BN heterostructures have been synthesized [6]. In these nanostructures, boron and nitrogen species are known to form hexagonal BN domains embedded into the carbon network.

In parallel to these experimental investigations on growth and characterization, theoretical studies have been carried out on the BN domain geometrical structures and on the electronic transport properties of B-or N-doped nanotubes [7,8]. For a single substitution in the carbon honeycomb lattice, it was established that the transmission function presents minimum * sylvain.latil@cea.fr values corresponding to the quasibound state (QBS) induced by the dopants [7]. Conductance is damped for injection energy around 0.5 eV above (below) the Fermi level E F in the case of N (B) substitution. The QBSs affect the electronic conduction in broader ranges when distributions of dopants are taken into account. The conduction regimes have been analyzed and the evolution of the electronic mean free path (MFP) has been estimated, as functions of dopant density, nanotube radius, applied magnetic field, etc. [9][10][11]. Besides, the electronic structure of boron and nitrogen codoped nanotubes has been investigated, with special attention given to the variation of the electrostatic potential along the tube. It has been pointed out that electronic gaps could arise in specific situations [12].

Focusing on the electronic transport in boron-and nitrogensubstituted nanotube heterostructures, Khalfoun and coworkers reported first that the effect of the B and N elements on the transport properties is highly dependent on the relative position of the two doping species [13]. When the substitutional atoms are first nearest neighbors, the QBSs repel each other and the effect of the two dopants is annihilated, which preserves the transmission at E F , whereas the carrier diffusion by each quasibound state is restored when the two dopants are spaced, resulting in a strong reduction of the electronic transmission properties. Second, in the case of nanodomains, they showed that the mutual interaction of the N-and B-induced QBSs results in a discrete spectrum of retrodiffusion states. The energy of the resulting QBSs induced by the shape and the nature of the defects or domains is then the key parameter driving the transmission function at E F .

In this paper we present our theoretical investigations on the electronic transport properties of BN-doped single-wall CNTs. All the structures have the B x C 1-2x N x stoichiometry, with x the doping ratio. The model for the electronic structure is presented in Sec. II A. Section II B describes the techniques used to calculate the electronic transport and to extract the mean free path. We have investigated the evolution of the MFP as a function of radius and doping ratio for metallic (n,n) and semiconducting (n,0) single-wall CNTs. Our results are presented in Sec. III.

II. THEORETICAL MODELS

To describe the electronic properties of carbon nanotubes containing a selected density of combined BN defects, we need first a model for the electronic structure enabling a reliable description for a wide range of structures, with the same accuracy and electronic reference energy. On the other hand, we want to characterize electronic properties at the mesoscopic scale (μm size), hence dealing with nanotubes containing up to 10 6 atoms. Furthermore, to estimate how a given density of BN defects affects the transport, statistics on different configurations are needed. We then adopt a tight-binding (TB) model that provides a good balance between reliability and efficiency.

A. Basis set, parametrization

Up to now, most calculations at the mesoscopic scale have been carried out with the simplest tight-binding approaches, usually in an orthogonal basis set. That model gives a good description of the electronic structure for pristine nanotubes. To account for boron and nitrogen substitutions in the network, an electrostatic correction was added on the diagonal part of the H matrix only [8,9,11]. That method is unfortunately restricted to the simplest situations for which the combined influence of boron and nitrogen can be neglected: low density of defects, pure N or B doping, etc.

With a very different approach, Lopez-Bezanilla et al. [14] have studied the transport in functionalized nanotubes, with a covalent link between the functional group and the nanotube wall. In this situation, the simplest TB approach is not reliable, mainly because sp 2 and sp 3 hybridizations coexist in the simulation cell. Model Hamiltonian and overlap matrices (expanded on an orbital atomic basis set) were directly constructed from the density functional theory (DFT) Hamiltonian and overlap matrices (H and S) obtained after the self consistent field treatment. This highly reliable approach preserves the accuracy of the ab initio electronic description. This method has also been successfully applied by Khalfoun and coworkers to study the effect of combined B and N substitutions in nanotubes [13]. It is unfortunately not tractable in the present problem, since the heavy basis set employed in DFT (eight orbitals per atom with long-range hopping integral terms) leads to time-consuming calculations, which limit the statistical averages on μm length CNTs.

In our model, a nonorthogonal TB formalism, we adopt a basis set formed by a single p ⊥ atomic orbital per atom. The curvature of the nanotube is consequently neglected, because it is not an important feature for the electronic states in the vicinity of the Fermi level. 1The determination of the matrix elements of the Hamiltonian and overlap matrices (H and S) is made by fitting matrix elements coming from DFT calculations. Reference DFT calculations were performed using the SIESTA code [16], with the LDA approximation for exchange and correlation [17] prior to the TB parametrization. Core electrons were treated by norm-conserving pseudopotentials [18] and a double-ζ s, p basis set was used to solve the valence electron eigenproblem. The basis set has been optimized using the simplex downhill method [19,20]. The unit cell is a graphene-like system with periodic boundary conditions made from 600 atoms, allowing us to reasonably reduce the Brillouin zone sampling to the point.

First the analysis of overlap matrix elements in DFT calculations indicate that the overlap matrix is significant up to the third neighbors for each element considered. It indicates that the hopping integrals also could be neglected after the third neighbors. For reference, we consider a pristine supercell and then address different doping situations perturbing the initial Hamiltonian H 0 obtained for a pure C system. Isolated B or N substitution first modifies the on-site energies of the neighboring carbon atoms [8,13]. This correction is fitted by a sum of three Gaussian functions. It can then be considered as a Coulombian potential perturbing the Hamiltonian H 0 previously set. Moreover, the presence of such a heteroatom in the network also affects the hopping integral elements. Assuming that the perturbative potential smoothly evolves on the range of one atomic orbital, we used a Hückel formulation to describe the influence of the dopant on nondiagonal Hamiltonian elements:

H ij = H 0 ij + δ i + δ j 2 S ij , (1)
where δ correspond to the on-site corrections, that is, the electrostatic potential induced by the dopant. The evolution of this perturbative potential is plotted in Fig. 1 as a function of the distance from the heteroatom. On the same figure, we report the correction induced on the hopping integral between C first neighbors as a function of the mean distance of these C to the dopant. This method is very tractable since all of the TB parameters only depend on the on-site energy correction. When both B and N substitutions are present in the supercell, the perturbative potential induced on C atoms is well described by the sum of isolated corrections [13]. In this codoping situation, we also paid attention to the evolution of on-site energy of dopants perturbing each other. Figure 2 reports the evolution of the on-site correction for various BN substitution geometries. Our fits based on isolated B or N are in good agreement with the DFT calculations for all the pair geometries.

Small domains formed by a cycle of three B and three N atoms, denoted B 3 N 3 are also considered. Unfortunately, these smallest hexagonal BN domains cannot be described with the previous additive model, since the exchange part starts to be non-negligible. For this specific case, our analysis indicates that the perturbation induced onto the matrix elements on the surrounding C atoms affects the first neighbors only, as depicted in the Fig. 2(d). For that reason, the Hamiltonian elements are then directly taken from the DFT output in our parametrization. These two methods are reliable and can be compared since they refer to similar DFT calculations, performed at the same accuracy level.

B. Electronic transport properties

To calculate the conductance of our systems, we employed the recursive Green's function formalism in real space [21][22][23] generalized to nonorthogonal basis sets [24]. This order N approach allows us to study systems containing up to 10 6 atoms. The system is divided into three parts: two semiinfinite pristine electrodes (labeled L and R) connected to the channel region (the central part) with the condition that there is no direct tunneling Hamiltonian between the two electrodes. The projection of Green's function on the central region is

G C = (ES C -H C -L -R ), (2)
where the terms are the self-energies corresponding to the coupling between the central part and the L and R electrodes. They write

X = (ES CX -H CX )g X (ES XC -H XC ), (3)
where X could be R or L and g X is Green's function for the isolated semi-infinite electrode. They are computed using the TB_Sim package [25], in which the recursive Green's function technique is implemented. The Landauer-Buttiker conductance [26] can then be evaluated from the transmission function

T (E) = Tr[ L G C R G † C ], (4)
where the terms are defined by

X = i( X - † X ). (5)
In this study the semi-infinite electrodes are treated as perfect nanotubes and the central part corresponds to the doped region.

To investigate the influence of the BN substitutions on the transport properties in a large scale, we have computed the transmission function for a set of random distribution of defects. The mean conductance is determined by averaging over 200 configurations, assuming each one is equiprobable. Figure 3(a) presents the transmission function of a pristine (10,10) nanotube, a nanotube containing a single BN pair, and the average transmission for a 180-nm-long (10,10) nanotube with a doping ratio of 0.1%. From the transmission function, the MFP can be extracted from an interpolation between the quasiballistic and diffusive regimes [27]:

T (L,E) = N (E) 1 + L e , ( 6
) -1 0 1 E -E F (eV)

FIG. 3. (Color online) (a)

The transmission function of a pristine (10,10) nanotube (dashed, black), a (10,10) nanotube containing a single BN-pair substitution (dotted, red), and one with a doping ratio of 0.1% and a length of 180 nm (solid, blue). (b) The transmission calculated at E F as a function of the length of the doped region (black squares). The solid line is obtained by fitting the calculated values with Eq. ( 6), which gives an estimation of e = 335 nm in this example. with N the number of accessible transmission channels in the similar pristine nanotube for the injection energy E. L is the length of the doped region, and l e is the electronic MFP. This interpolation scheme is illustrated in Fig. 3(b). The calculated average transmission is plotted (squares) as a function of L at the Fermi energy E F . The MFP is extracted by fitting T (L,E F ) with Eq. ( 6).

III. RESULTS AND DISCUSSION

We have considered five different BN geometries in the hexagonal network. Heteroatoms can be located as first or third neighbors, or separated by two or four carbons, as depicted on Fig. 4. We have investigated, for each geometry of BN substitutions the evolution of the transport properties as a function of the doping ratio and the radius for armchair nanotubes.

Figure 5(a) presents the evolution of the MFP for a (10,10) nanotube doped with 0.1% of B and N, as a function of the energy. The MFP is nearly zero in the vicinity of the first Van Hove singularity, close to the energy of the QBS induced by the B and N substitutions. Moreover, the MFP becomes very small for the energy ranges beyond the Van Hove singularity where the group velocity decreases drastically, as previously reported [28].

The MFP presents an inverse dependance with the density of defects e = A/ρ and a linear evolution with the radius e = αR, following the theoretical predictions based on Anderson-like disorder modeling [29]. The corresponding curves are respectively presented on Figs. 5(b) and 5(c). A relevant property is that the prefactors A and α decrease when the BN distance increases, which is reminiscent of the energy scale of the QBS for the single diffusion center situation [13]. The scattering potential induced by the defect is softened when B and N atoms get closer, thus the quasibound states created by the BN defects are shifted away in the vicinity of the van Hove singularities. All these results indicate that the BN codoping behaves as a classical single diffusion center in the low doping limit. 2 Electronic transport in nanotubes containing BN domains also follows the previous laws concerning radius and concentration of defects. For that reason, the small domains B 3 N 3 can also be considered as point defects in the nanotube. On one hand, the induced scattering potential is spatially more extended than in the BN first-neighbors case and thus may affect the electronic transmission in a stronger manner. On the other hand, for a similar BN concentration the number of defects is three times lower than for the pair case. As a result the MFP for BN cycles is larger than for the BN pairs with the same BN doping ratio.

It is worth noticing that the energy formation for a given BN defect increases with the distance between the species. Besides, structural studies of BN nanodomains in C-BN heterostructures indicate that the BN and CC bonds in the network are energetically favorable [30]. The system tends to maximize these kinds of bonds, leading to the formation of BN domains [31]. 2 When B and N substitutions get farther from each other, their combined effect is lost and the obtained mean free paths can be compared with those of a fully random distribution of B and/or N atoms, as in Refs. [9,11]. The full lines are obtained by fitting the evolution with e = A/ρ law. The obtained value of A are respectively 0.9, 0.6, 0.53 and 0.35 nm for the pairs (with an increasing distance between B and N, as shown on Fig. 4) and 0.99 nm for the B 3 N 3 domain. (c) Estimation of the MFP at E F , with a doping ratio of 0.1%, as a function of the radius R. The full lines are obtained by fitting with a linear law E F = αR. For the pairs the α term is respectively estimated to 1320, 900, 750, and 510, whereas we found α =1540 for the domain. It is interesting to point out that the transmission function presents an asymmetric behavior with respect to the electronic injection energy. Indeed, the transport properties are more affected for energies above the Fermi level than below. The maximum MFP value is reached for an energy around 0.25 eV below the Fermi level. This feature is not intuitive when considering previous studies addressing single defects only. A single BN pair, or a B 3 N 3 cycle, seems to be transparent for incident electrons, and the transmission function is just affected around the Van Hove singularities [13]. Actually, the potential well induced by nitrogen is deeper than the potential hill induced by boron. As a consequence, the QBS induced by a single substitution of boron is closer to the first Van Hove singularity than the one coming from the nitrogen case, explaining this asymmetry in the MFP curve.

We have also paid attention to a semiconducting nanotube, with the zigzag symmetry. The results we obtained for a (19,0) nanotube with a 0.1% doping ratio are presented on Fig. 6(a). The evolution of the MFP as a function of the energy indicates that this asymmetric behavior is still present for the semiconducting CNTs. Energy ranges corresponding roughly to the valence and conduction bands have been selected respectively to -0.4 and +0.4 eV. The evolution of the MFP as a function of the doping ratio is depicted in Fig. 6(b) for an energy range corresponding to the valence band and Fig. 6(c) for the conducting band.

Again, the electronic mean free path is inversely proportional to the dopant density, and the prefactor decreases with the distance between the B and N atoms. For the energy range corresponding to the valence band state, a deviation to the A/ρ law is observed for high dopant concentrations. This feature is related to a deformation of the band structure induced by the dopant that shifts the valence band. For the states in the conduction band, no deviation was observed, which suggests that the conduction band is less sensitive to the dopant than the valence band. The asymmetric behavior is also more pronounced for distant BN couples, two times larger for the valence band than for the conducting band. This asymmetry tends to disappear when the distance between B and N decreases as well as for small domains.

IV. SUMMARY

In conclusion, we have investigated transport properties of carbon nanotubes with random distributions of BN pairs and small domains. The electronic mean free path decreases when doping species are dispersed, which could have important applications since the doping process in carbon nanotubes has a tendency to segregate BN into domains. We report that all considered BN geometries can be treated as classical point defects, and that the mean free path is inversely proportional to the doping ratio and proportional to the radius.

We also report an asymmetric behavior of the mean free path with injection energy. As a perspective it could be interesting to study the behavior of larger BN domains and then estimate the critical domain size over which one of the previous laws cannot be applied.

FIG. 1 .

 1 FIG. 1. (Color online) Evolution of corrections added to the C on-site energy as a function of the distance to the single atom dopant. The evolution of the first neighbors hopping integral between two C atoms are also plotted. The black points are the SIESTA matrix elements, and the solid lines correspond to the parametrization obtained by fitting the on-site energy by a sum of three Gaussians and with Eq. (1) for the hopping integrals.

FIG. 2 .

 2 FIG. 2. (Color online) The evolution of the correction added to the on-site energy for different types of combined B and N substitutions along an armchair line (sites surrounded by the rectangle): (a) nearest neighbors BN pair, (b) B and N atoms separated by 2.48 Å, (c) B and N atoms separated by 7.1 Å, and (d) a B 3 N 3 nanodomain. The black points are the directly read SIESTA matrix elements and the solid line corresponds the previously setup TB model. For the B 3 N 3 geometry, the TB model elements are not obtained by the fit, but directly taken from the SIESTA output.

FIG. 4 .

 4 FIG. 4. (Color online) Different BN substitution geometries considered in this study. The distances between B and N for the pair cases are respectively 1.42, 2.82, 3.75 and 6.15 Å.

FIG. 5 .

 5 FIG. 5. (Color online) Transport properties in a BN-doped (10,10) nanotube. (a) The evolution of the mean free path as a function of the injection energy, for the different BN defects sketched in Fig. 4. (b) Estimation of e at the Fermi energy, as a function of the BN doping ratio.The full lines are obtained by fitting the evolution with e = A/ρ law. The obtained value of A are respectively 0.9, 0.6, 0.53 and 0.35 nm for the pairs (with an increasing distance between B and N, as shown on Fig.4) and 0.99 nm for the B 3 N 3 domain. (c) Estimation of the MFP at E F , with a doping ratio of 0.1%, as a function of the radius R. The full lines are obtained by fitting with a linear law E F = αR. For the pairs the α term is respectively estimated to 1320, 900, 750, and 510, whereas we found α =1540 for the domain.

FIG. 6 .

 6 FIG. 6. (Color online) Transport properties in a BN-doped (19,0) nanotube. (a) Evolution of the mean free path as a function of the injection energy. (b) Mean-free-path estimate at -0.4 eV below the Fermi level as a function of the doping ratio for a (19,0) nanotube. The full lines correspond to a fit by a A/ρ law. The A values are 0.65, 0.41, 0.31, 0.17 nm for the BN pair and 0.97 nm for the cycle. (c) Same data and fits at +0.4 eV over the Fermi level. The A values are 0.59, 0.34, 0.21, 0.10 nm for the BN pair and 0.93 nm for the cycle.

In the present work, we pay attention to armchair and zigzag nanotubes only. However, an extension of our model to study chiral (n,m) nanotubes is possible by adding an angular dependance of the hopping elements (with respect to the nanotube axis) as shown in Ref.[15]. Such a model method allows us to predict the small gap opening in (quasi)metallic nanotubes.

ACKNOWLEDGMENTS

Prof. Luc Henrard and Dr. Cyrille Barreteau are acknowledged for fruitful discussions about this work. This work was financed by the Nano-Simulation contract coming from the CEA Transversal Nano-Science Programme, and by the French National Research Agency (ANR) through the NANOSIMGRAPHENE Project No. ANR-09-NANO-016-01.


D G Porter

M Roger

M J Gutmann

S Uthayakumar

D Prabhakaran

A T Boothroyd

M S Pandiyan

J P Goff

Divacancy superstructures in thermoelectric calcium-doped sodium cobaltate

We have grown single crystals of NaxCayCoO2 and determined their superstructures as a function of composition using neutron and x-ray diffraction. Inclusion of Ca 2+ stabilises a single superstructure across a wide range of temperatures and concentrations. The superstructure in the Na + layers is based on arrays of divacancy clusters with Ca 2+ ions occupying the central site, and it has an ideal concentration Na 4/7 Ca 1/7 CoO2. Previous measurements of the thermoelectric properties on this system are discussed in light of this superstructure. Na 4/7 Ca 1/7 CoO2 corresponds to the maximum in thermoelectric performance of this system.

I. INTRODUCTION

Sodium cobaltate (Na x CoO 2 ) has emerged as a material of exceptional scientific and technological interest due to its potential for thermoelectric applications that convert waste heat into electricity or allow solid state refrigeration [START_REF] Terasaki | [END_REF][2][3][4] . In Na x CoO 2 , the removal of sodium ions creates holes in the cobalt layers. The Seebeck coefficient and power factor both increase with x, as the concentration of holes decreases 4 . Above a concentration of x ≈ 0.85, however, phase coexistence occurs with the insulating x = 1 phase, causing a reduction in the power factor. The substitution of a divalent ion for sodium decreases by one the number of holes in the cobalt layer. This offers the possibility to lower the hole concentration and enhance further the thermoelectric performance. There are already promising results on polycrystalline samples doped with Ca and Sr, with an increase in the power factor over the comparable pure compound [5][6][7] .

Neutron diffraction studies have demonstrated the central role played by superstructure on the thermoelectric properties of Na x CoO 2 8,9 . The structure of Na x CoO 2 comprises layers of conducting cobalt-oxygen octahedra separated by layers of sodium ions. The sodium ions can sit on one of two inequivalent sites; Na 2b sites directly between successive cobalt ions, or displaced from these position on Na 2d sites. The 2d sites are preferentially occupied due to a small energy cost for sitting directly above cobalt arising from short-range repulsion. Vacancies in the sodium layer have been found to form multi-vacancy clusters at certain fractional concentrations. These clusters form by promoting sodium ions onto the 2b sites, creating rattling cages that can reduce thermal conduction through the material 10 . Coulomb interactions cause the vacancy clusters to arrange themselves into long range patterns, causing rings of superlattice peaks in diffraction images.

The electrostatic potential generated by these super-lattices templates the Coulomb landscape in the cobaltoxygen layer, generating confined electronic pathways in the layers 8 . The combination of electron confinement and rattling cage-ions creates an ideal environment for the thermoelectric effect.

Here we present the results of an x-ray and neutron diffraction study of superstructures in Ca-doped Na x CoO 2 . We show that the Ca stabilises a single superstructure over a wide range of temperature and composition, and we use the results to gain an understanding of how Ca-doping affects the thermoelectric properties.

II. METHODS

Single crystal boules of various compositions of Na

x Ca y CoO 2 were grown using the floating zone technique 11,12 . Small, high quality crystallites were cleaved from the boules and studied using x-ray diffraction with a molybdenum source, Xcalibur diffractometer (Agilent instruments) at Royal Holloway. A Cryojet5 (Oxford Instruments) was used for temperature control and precise temperatures were calibrated by replacing the sample with a thermistor of similar size. A large, single grain of Na 0.7 Ca 0.1 CoO 2 was cleaved from its boule and studied using neutron diffraction on SXD 13 at ISIS.

The hexagonal reflections from the diffraction experiments were integrated using standard instrument methods, and properties of the average structure were refined using the program Jana2006 14 . A reverse Monte Carlo (RMC) program was written to determine the vacancy ordering within the sodium layer and to establish the position of calcium ions. RMC uses simulated annealing to minimise the chi-squared difference between the calculated intensities arising from a given set of atomic positions in the unit cell, and the experimental data. The positions of Ca and Na ions are changed with respect to the Co lattice, using the Metropolis algorithm in the Canonical ensemble. Small displacements of all atoms, with respect to their original positions, are also allowed. When the fictitious temperature is slowly decreased, the configuration corresponding to the absolute minimum of this chi-square difference can be reached. Further details on this method can be found in Refs. [ 15,[START_REF] Porter | Control of electrons through patterning of superstructures in sodium cobaltate[END_REF] ].

III. DIFFRACTION EXPERIMENTS

Three compositions of the calcium doped system were grown using the floating zone technique: Na 0.65 Ca 0.1 CoO 2 , Na 0.57 Ca 0.14 CoO 2 and Na 0.55 Ca 0.2 CoO 2 . The samples were grown as large 4-6 cm boules which were cleaved to find small, high quality crystallites. These were screened with the x-ray diffractometer and high quality samples were measured with longer exposures and greater coverage at several temperatures. Figure 1 shows reciprocal space cuts in the (h, k, 0) plane from a single-crystal sample of each composition. All three samples exhibit rings of 12 satellite reflections around the principal hexagonal Bragg peaks. These additional peaks lay on the vertices of a hexagonal grid with spacing a * /7 and can be indexed using a supercell with lattice vectors:

a = 2a -b b = a + 3b c = c
These supercell vectors index one of the two possible superlattice domains, where the second domain can be generated by a reflection in either the {1, 0, 0} or {1, 1, 0}. The intensities of the superlattice reflections indicate equal populations of these domains. No broadening of the peaks was detected within instrumental resolution and, therefore, the superstructures were found to be ordered long range. The positions of these superlattice peaks are not consistent with any of the diffraction patterns observed previously in pure Na x CoO 2 , signifying that a new superstructure has been found. All three samples have additional reflections halfway between the largest Bragg reflections due to contamination by λ/2, in addition to this the Na 0.55 Ca 0.2 CoO 2 sample displays further reflections between the primary satellites. These additional reflections lay at the vertices of an a * /4 hexagonal grid and can be indexed by a supercell with lattice vectors: These propagation vectors were first observed for Na 0.5 CoO 2 [START_REF] Huang | [END_REF] , implying a coexistence between the Na 0.5 CoO 2 superstructure and the new a * /7 superstructure.

a = 2a b = a + 2b c = c
The intensities of the satellite peaks vary significantly for different l planes as shown in Fig. 2. Satellites occur primarily in 12-fold rings around the {1, 0, l} reflections, with the modulation of peak intensity within the rings varying between even and odd planes. At high l, the peak intensity is diminished due to the x-ray form factor, however peaks can still be observed and show that satellites around the {1, 1, l} reflections gain intensity, especially in the (h, k, 11) plane. The increased satellite intensity and reduced intensity in the principal reflections in the (h, k, 11) plane is particularly striking and it is related to the fact that the c-lattice parameter is approximately 11 times the distance between the cobalt and oxygen planes.

The diffraction patterns were measured as a function of temperature between 90 and 460 K. It was found that the satellite peaks do not change below room temperature, though they do eventually disappear at high temperature, as illustrated in Fig. 3. The disappearance of the satellite peaks is attributed to the sodium layers becoming disordered. We did not observe any variation of the peak widths in this temperature range within our experimental resolution. The ordering temperature increases as the calcium concentration increases, and the order parameter indicates a sharper transition for the y = 0.2 sample. The a * /4 peaks did not vary at all within this temperature range.

The peak intensities in Fig. 3 are normalised against the intensity of a nearby principal hexagonal reflection in each case. The Na 0.65 Ca 0.1 CoO 2 intensity is much weaker than the other compositions and it has a lower disorder temperature, with a fitted order parameter giving T d = 379 ± 6K. There is little difference between the initial intensity or the disorder temperature of Na 0.57 Ca 0.14 CoO 2 and Na 0.55 Ca 0.2 CoO 2 , with T d = 401 ± 2K and T d = 410 ± 10K respectively. This suggests that the a * /7 satellites in both samples originate from a phase with a common calcium composition. The additional calcium in the Na 0.55 Ca 0.2 CoO 2 sample is, therefore, likely to be part of the coexisting phase.

A boule of Na 0.7 Ca 0.1 CoO 2 was screened on SXD. A large single grain with intense, sharp principal Bragg peaks was cleaved from the boule and longer exposures of several hours were taken at 350K and 40K. It was not possible to grow large enough single grains for neutron diffraction studies of the superstructures with higher calcium concentrations.

Superlattice peaks can be distinguished in the data by their intensities having hexagonal symmetry. Satellite peaks were observed with the same propagation vectors used to index the a * /7 superlattice observed with x-rays, indicating that this phase forms in the bulk material, as shown in Fig. 4. Several shorter exposures were performed at increasing temperatures and the disorder transition observed was consistent with the Na 0.65 Ca 0.1 CoO 2 measurement performed with x-rays, see Fig. 3. The superlattice pattern did not change as a function of temperature, indicating a single phase below the ordering transition temperature.

In addition to the satellites, further peaks were observed close to the {110} principal hexagonal reflections. The d-spacing of these peaks was ≈ 3.70 Å, which matches the {220} reflection of calcium oxide (CaO). By comparing the integrated intensities for several reflections with calculated structure factors, and averaging over domains, the relative proportion of CaO in Na 0.7 Ca 0.1 CoO 2 was determined to be 0.8% with a statistical error of 0.2%. For this system to produce a hexagonal pattern of points as found in these results, the cubic [111] direction of the CaO crystal would have to be parallel to the hexagonal [001] of Na 0.7 Ca 0.1 CoO 2 , and the cubic [110] direction of CaO would be parallel to the hexagonal [100] of Na 0.7 Ca 0.1 CoO 2 . In this direction the CaO structure can be described as successive hexagonal layers of calcium and oxygen ions with an ABC stacking sequence. The hexagonal environment of the main phase allows the CaO impurity to grow epitaxially within the Na 0.7 Ca 0.1 CoO 2 boule.

IV. SINGLE CRYSTAL REFINEMENTS

The principal hexagonal reflections in the multicomposition x-ray datasets were integrated using the three-dimensional (3D) profile method. The results of the refinements, as recorded in Table I, show that the sodium occupancies are consistently lower than the nominal values. The calcium occupancies are much closer to the nominal values however, which is consistent with loss of the more volatile sodium ions during growth. For Na 0.57 Ca 0.14 CoO 2 the fitted Ca occupancy agrees with the nominal composition for the proposed Na 4/7 Ca 1/7 CoO 2 superstructure. This is consistent with the lack of any observed coexisting phases in this sample.

The thermal parameters were refined as isotropic spheres. The data in Table I are consistent with increased static disorder away from the optimum composition. The ADPs for the sodium ions on the 2d sites are comparable to those obtained for the pure system 10 , whereas those of the divalent calcium ions on the 2b sites are significantly lower. A full discussion of ADPs in terms of the lattice vibrations is beyond the scope of the present paper.

The R-factors for these refinements show that there is excellent agreement between the refined models and the experimental data. Other refinements were attempted that changed the assumption that Ca ions occupy 2b sites and Na ions sit on 2d sites. These resulted in unrealistic concentrations of calcium or negative thermal parameters, and these models were, therefore, rejected.

V. REVERSE MONTE CARLO SIMULATIONS

Intensities from the superlattice reflections from Na 0.57 Ca 0.14 CoO 2 were obtained from the room temperature x-ray data by indexing a single superlattice domain and excluding the principal hexagonal reflections. The superlattice reflections arise purely from the superstructure, whereas the principal reflections are only slightly modified, and they are also affected by other factors such as multiple scattering, or coexisting domains and phases. The propagation vectors used to index the rings of superlattice reflections converts to a real space supercell of seven hexagonal sodium cobaltate unit cells, and the concentrations from the refinements meant that four sodium ions and one calcium were added to each layer. RMC calculations were performed using this set of intensities and supercell. Each iteration of the RMC program involves discrete hopping movements of the sodium and calcium ions to vacancy sites within the plane, as well as continuous movements of any ion by small increments in any direction. Simulated annealing of the system helped to avoid a false minimum solution, initially allowing ions to hop freely then slowly decreasing the fictitious temperature and reducing the probability of allowing movements that would produce a higher χ 2 .

The RMC calculations led to the robust solution illustrated in Fig. 5, with the calcium ion on a 2b position surrounded by sodium ions sitting on 2d sites. In the second sodium/calcium layer the same pattern emerges, however, the calcium ion is translated away from the location of calcium on the first layer. This translation relates to a significant reduction in χ 2 . The fit for this model had a χ 2 per degree of freedom of 5.82 and an R w value of 12.9%. Figure 2(b) shows the calculated x-ray patterns for this solution indicating good agreement with the experimental data in Fig. 2(a) across all l planes, including the variation of the satellite peak intensity between even and odd planes. The resulting structure has the space group P 2 1 /m with unique axis c, atomic positions are given in Table II.

The RMC solution also exhibits ordering in the distortions of the cobalt plane, with cobalt ions moving away from the closest calcium (which is sitting directly above the cobalt). This form of distortion was predicted using a Coulomb model for the pure system 8 . The magnitude of the distortion is smaller for the system doped with calcium since the decreased hole concentration reduces the electrostatic potential for the multi-vacancy clusters. Figure 5(a) depicts the calcium doped divacancy structure with cobalt distortions emphasized for clarity. It was not possible to resolve the variation in the Co-O bond lengths using our x-ray data and, therefore, we attempted to use the neutron data.

RMC simulations using the neutron diffraction data were less sensitive to the superlattice, mainly due to a lower signal-to-background. The same divacancy superstructure was obtained, but the smaller scattering contrast between the calcium and sodium ions made it difficult to locate the positions of the calcium ions.

For Na 0.55 Ca 0.2 CoO 2 this divacancy phase coexists with the phase observed for Na 0.5 CoO 2 as shown in Fig. 6. In fact, the Na 0.5 CoO 2 superstructure can also be viewed as stripes of divacancy clusters 8 . In this case it was not possible to perform RMC simulations due to the presence of contamination by higher order diffraction peaks.

VI. DISCUSSION

Analysis of the superlattice observed in all samples of Na x Ca y CoO 2 indicates the presence of the superstructure with an ideal concentration of Na 4/7 Ca 1/7 CoO 2 . For pure sodium cobaltate the corresponding superstruc- ture has a sodium concentration x = 5/7. According to calculations of the ground state energies of multivacancy structures using a long-range Coulombic model with short-range repulsion, this is the only stable divacancy cluster in this composition range 8 . It is likely that the substitution of the divalent calcium at the centre of the divacancy cluster has a stabilising effect since it lowers the charge of the divacancy cluster. The fact that this superstructure does not change with temperature and forms over a range of sample compositions suggests that it is a particularly stable superstructure. Previous studies involving the doping of sodium cobaltate with calcium have concentrated on the variation transport properties relevant to the thermoelectric performance of the material and had not considered the role of the superstructure. This has led to a number of different doping regimes being attempted though none have used exactly the concentration of the superstructure determined here. Kawata et al. found that both the electrical resistivity, ρ, and the Seebeck coefficient, S, increase with increasing calcium concentration 5 . Ono et al. found that the increase in both ρ and S was linear, so that overall the power factor, P = S 2 /ρ, also increased with calcium concentration 6 . However, the power factor reaches a maximum value for Na 0.70 Ca 0.0175 CoO 2 and above this calcium concentration, P decreases. In another com- position series, Li et al. found that the power factor first increased, and then decreased for Na 0.64 Ca 0.16 CoO 2 and higher calcium concentrations 7 . In these reports the power factor is peaked very close to the optimum concentration for the ideal Na 4/7 Ca 1/7 CoO 2 superstructure. The Seebeck coefficient increases up to this composition. Above this composition, where we find coexistence with an additional insulating phase, the electrical conductivity decreases.

For a full understanding of the thermoelectric performance of this material, the effect of the divacancy superstructure on the thermal transport must also be understood, as this is a key component for calculating the thermoelectric figure-of-merit, ZT = T S 2 /ρκ, where κ is the thermal conductivity. This contribution is beyond the scope of the current work but will be addressed in a future publication.

VII. CONCLUSIONS

In summary, both x-ray and neutron diffraction of calcium-doped sodium cobaltate have revealed a previously unobserved superlattice that can be explained by a superstructure with seven unit cells. An RMC program was used to find the ordering in both sodium layers as well as distortions in the cobalt plane in a model-independent manner, reliably generating a divacancy structure with calcium ions at the centre on a 2b site, with an ideal concentration Na 4/7 Ca 1/7 CoO 2 . The same superstructure is observed across a broad temperature range and in all compositions studied, though moving away from the optimal concentration for this structure results in weaker superlattice intensities and the emergence of additional phases. Comparison with thermoelectric measurements in the literature indicates that this ideal concentration coincides with the maximum in the thermoelectric power factor for the system.

FIG. 1 .

 1 FIG. 1. (Colour online) Reciprocal space cuts through the (h, k, 0) layer for different compositions of NaxCayCoO2, measured by single-crystal x-ray diffraction. (a) Na0.65Ca0.1CoO2 shows a 12-fold ring of satellite peaks around the principal Bragg reflections which can be indexed on a hexagonal grid with spacing a * /7. (b) Na0.57Ca0.14CoO2 displays the same superlattice structure, but with stronger intensity. (c) Na0.55Ca0.2CoO2 also displays the a * /7 superlattice peaks, plus additional peaks that lay on a a * /4 lattice.

FIG. 2 .

 2 FIG. 2. (Colour online) (a) Cuts through the reciprocal space volume for Na0.57Ca0.14CoO2 at integer l values, illustrating the variation of intensity in the superlattice peaks with l. Gaps in the data are a result of missing coverage. (b) Calculated intensities from the model generated by reverse Monte-Carlo, for comparison with (a).

FIG. 3 .

 3 FIG. 3. (Colour online) Temperature dependence of the superlattice peaks for three compositions of NaxCayCoO2. Each peak intensity is normalised by a close principal hexagonal reflection, therefore indicating the relative intensity variation between compositions. Fitting an order parameter to each dataset gives disorder temperatures of T d = 379 ± 6K for Na0.65Ca0.1CoO2, T d = 401 ± 2K for Na0.57Ca0.14CoO2 and T d = 410 ± 10K for Na0.55Ca0.2CoO2. The neutron data points indicate measurements during heating and cooling, showing no sign of hysteresis behaviour and are consistent with the x-ray measurement.

FIG. 4 .

 4 FIG. 4. (Colour online)A cut in the (h, k, 0) plane for Na0.7Ca0.1CoO2 SXD data at 40K. Satellite peaks appear in the same positions as found in the x-ray data, indicating that this is a bulk phase. Black circles and squares distinguish the two domains for the superlattice, and the red diamonds mark the {220} reflections of the epitaxial CaO impurity. Other non-indexed peaks are a result detector anomalies and noise resulting at low detector angles.

FIG. 5 .

 5 FIG. 5. (Colour online) Results from the RMC calculation of Na0.57Ca0.14CoO2 x-ray data. (a) The resultant structure, where blue balls are sodium ions, purple are calcium ions, the other balls are cobalt ions, where dark green is a large negative distortion and light green is a large positive distortion. Oxygen ions have been omitted for clarity. (b) Annealing of χ 2 during the calculation. Inset: Final comparison of calculated (Ic) and experimental (Io) intensities, where the black line displays Ic = Io.

FIG. 6 .

 6 FIG. 6. (Colour online) Na0.55Ca0.2CoO2 can be described by a coexistence of the Na0.57Ca0.14CoO2 supercell in (b) and the x = 0.5 supercell in (c). The calculation of the (h, k, 0) layer in (a) gives good agreement with the xray data from Fig. 1(c). Superlattice peaks generated by the Na0.57Ca0.14CoO2 supercell are indicated with circles, whereas the peaks produced by the x = 0.5 supercell are marked by triangles.

TABLE I

 I

	0, 0, 0), O ( 1 3 , 2 3 , 0.0930(8)), Na ( 2 3 , 1 3 , 0.25), Ca (0, 0, 0.25). Co and
	O occupancies were fixed at full occupancy.

. Refinement of NaxCayCoO2 principal Bragg reflections from x-ray diffraction data at 300K. The same model was used in each case: Space group: P 63/mmc, Positions: Co (

ACKNOWLEDGMENTS

We thank D. Voneshen for his help. We are grateful for the financial support and hospitality of ISIS. This work was supported by EPSRC grants EP/J011150/1 and EP/J012912/1.


M G Silly

G Li

Y J Dappe

J Dappe

Keywords: graphene, electronic properties, SiC

Electronic properties of zero-layer graphene on 6H-SiC(0001) substrate decoupled by silicon intercalation

Introduction

A graphene sheet constitutes a two-dimensional layer, formed by sp 2 -hybridized carbon atoms arranged in a honeycomb lattice. The band structure of graphene is governed by massless Dirac particles and exhibits characteristic linear band dispersion near the K-point linked to the half-filled π bonds. [1] This singular electronic structure sparks off novel electronic properties such as room temperature quantum Hall effect, [2] Klein tunneling [3] and high carrier mobility greater than in any known semiconductor. [4] Compatible with existing electronic device technology, graphene synthesis on SiC substrate is one of the most promising approaches for uniform coverage and structural coherence at wafer scale, a key factor to maximize the mobility and uniformity essential for applications. [5,6] Moreover, the use of wide-bandgap semiconductor SiC as semi-insulating substrate avoids hazardous transfer of the graphene layer to another insulating substrate at the origin of potential decrease in graphene electronic performances. [7] In addition, electronic properties of monolayer (ML) graphene, such as doping and electronic transport, are extremely dependent on the chemical nature of the interface. In fact, graphene grown on SiC(0001) Si-face exhibits mobilities an order of magnitude lower than on SiC(000-1) C-face, although structural graphene film quality appears comparable. [8] Hence, the control of the electronic properties of graphene layer appears to be of fundamental importance for the development of graphene electronics. [6] Therefore, in order to modify the graphene electronic properties, successful species intercalation as H, [9] O, [10] F, [11] Au, [12] Li, [13] Ge [14] and Si [15,16] has shown the high versatility of graphene doping from n-doped or p-doped materials to neutral quasi-free-standing graphene layer. As one ML graphene decoupling has been intensively investigated, zero-layer graphene or buffer layer decoupling process is still controversial, because of the complex covalent bonding between the zero-layer graphene and the substrate. [17,18] Recently, a significant study has been carried out on one ML graphene decoupling by silicon intercalation. The authors have also evidenced the decoupling of ZLG into graphene by Low Energy Electron Diffraction (LEED). [15] Nevertheless, the electronic properties of the resulting decoupled graphene layer remains unknown. The determination of the electronic properties is of crucial interest for potential applications as silicon intercalation is compatible with current electronic technology.

Here, we present a UHV study of the buffer layer decoupling from the SiC substrate by Si intercalation. We have investigated the evolution of the electronic properties of the ZLG by means of high-resolution photoemission spectroscopy (HRPES) and angle-resolved photoemission spectroscopy (ARPES) induced by silicon intercalation. The ZLG exhibits the typical linear dispersion along the Γ -K direction characteristic of a decoupled graphene layer. The decoupled layer presents a slight n-doping. Ab initio calculations confirm the transition from semiconducting to metallic behavior of the terminated carbon layer.

Methods

The experiments were performed using an on-axis-oriented n-type 6H-SiC(0001) research grade substrate with polished silicon face. The SiC substrate, initially outgassed at 600 °C for 12 h in a base pressure of 5 × 10 À10 mbars, is flashed in a pressure better than 3 × 10 À10 mbars at a temperature below 1200 °C for 30 s to terminate the substrate by the 6√3 × 6√3R30°reconstructed surface. All Si depositions and annealing were made in situ and characterized by LEED, ARPES and HRPES.

Both HRPES and ARPES were performed at the soft X-ray TEMPO beamline at Synchrotron SOLEIL. [START_REF] Polack | AIP Conf. Proc[END_REF] The beamline covers the energy range between 50 and 1500 eV with a resolving power better than 10 000. The spectra were measured using a high energy resolution Scienta SES2002 photoelectron analyzer. The experiments were performed in UHV with a base pressure of 3 × 10 À10 mbars. All the presented core level spectra were normalized with respect to the surface contributions. The core level spectra were deconvoluted according to standard curve-fitting procedure.

Ab initio calculations were performed using a very efficient Density Functional Theory (DFT)-localized orbital molecular dynamic technique (FIREBALL). [START_REF] Lewis | [END_REF][21][22][23] Basis sets of sp 3 for C and Si and s for H were used with cutoff radii (in atomic units) s = 4.5, p = 4.5 (C), s = 4.8, p = 5.4 (Si) and s = 4.1 (H). [24] In this study, we have considered supercells of five layers SiC(0001) with and without Si coverages, with a lateral size roughly corresponding to a 4 × 4 unit cell of graphene. The bottom layer is saturated with hydrogen atoms. On top of each supercell, we have set one 4 × 4 ML of graphene, and we have relaxed the system using a sample of 32 k-points in the Brillouin zone, maintaining the last three layers in bulk positions. A set of 300 special k-points along the path Γ -K-M has been used for the band structure calculations on the relaxed positions

Results and discussion

After a first flash in UHV of the 6H-SiC(0001) (Si-face) substrate below 1200 °C for 30 s, the surface exhibits the typical LEED pattern of the carbon-terminated 6√3 × 6√3R30°reconstructed surface (Fig. 1a). [25,26] This reconstruction composed of a carbon ML is also called ZLG. This epitaxial graphene layer is rotated by 30°with respect to the SiC-(1 × 1) surface cell. After several Si depositions at room temperature, corresponding to a total Si amount of about 1.75 ML and annealing at 750 °C cycles, this diffraction pattern disappears, exhibiting the bulk SiC-(1 × 1) pattern (Fig. 1c). Epitaxial single-layer graphene diffraction dots should appear close to the position of the 6√3 × 6√3 symmetry. [15] However, working at high energy, the contribution of the graphene mono layer is too weak to be resolved. Notice that the decoupling temperature of the ZLG (750 °C) appears to be lower than the temperature of one ML graphene decoupling (800 °C). [15] Finally, we can also notice an increase of the background intensity that may be attributed to inhomogeneities, unreconstructed surface or randomly distributed Si lying at the surface due to an excess of silicon at the surface.

Core level photoemission spectroscopy measurements are presented in Figs 1b and1d. The ZLG C1s core level is composed of three components (Fig. 1b): the two components at higher binding energy correspond to the carbon atom of the ZLG partially covalently bonded to the SiC substrate and the component located at lower binding energy corresponds to the SiC bulk. [25] After decoupling, the C1s core level shows two main components (Fig. 1d). The one at higher binding energy exhibits the characteristic asymmetric shape of metallic species and is attributed to the decoupled graphene layer. The component at lower binding energy corresponds to the SiC bulk. This contribution is shifted to the higher kinetic energy in good agreement with literature. [15] The ARPES measurements presented in Fig. 2 show the valence band structure collected along the Γ -K direction in the graphene The modulation periodicity appears to be 3/13 of the graphene Γ -K period, corresponding to the secondshortest reciprocal lattice vector of the 6√3 × 6√3R30°superstructure. [26] Another small shoulder rising at 1.6 eV below the Fermi level betrays the presence of g2 state partially screened by a continuum at the same energy. Between 2 and 5 eV appears an intense dispersive band presenting 5/13 of graphene periodicity, also observed in cone replica phenomena. [27] Up to 13 eV, the band structure is entangled by several broken structures. We also observe a mini gap opening, characteristic of band folding at reconstruction periodicities wave vector [28] (g1 modulation). The whole structure can be seen as the graphene π-band, modulated by the n-periodic superstructure. This corresponds to the n-folding bandstructure into the 1BZ coming from the incommensurate graphene layer unit cell and 6√3 × 6√3R30°unit cells [(13 × 13) graphene unit cell]. several cycles of silicon deposition and annealing, we have obtained the typical bandstructure of a decoupled graphene layer (Figs 2c and2d), i.e. the linear dispersion and the Dirac cone. The Dirac point is located at 250 meV below E F , in good agreement with the values for graphene decoupling by Si intercalation on iridium [29] and graphene on SiC(0001) exposed to air. [30] Consequently, we have calculated an electron density for the decoupled zero graphene layer of about 4.05 × 10 12 cm À2 close to the value obtained for epitaxial graphene layer on SiC(0001) exposed to air. [6,30] To model the ZLG decoupling process and confirm the structural arrangement of the terminated carbon layer during the silicon intercalation, we have calculated the bandstructure for two Si coverages. To calculate the ZLG, we started from truncated 6H-SiC(0001) substrate terminated by the pristine Si with a graphene layer on top. Then, the system has been relaxedusing DFTuntil the equilibrium. The final system (Fig. 3a) presents a covalently bonded graphene-like structure. This configuration is electronically stable because of conservation of mesomeric effects leading to a partial delocalization of the charges. However, the 2D symmetry of the graphene layer is broken, and the metallic character of the film is lost as shown in the band structure calculation (Fig. 3b), explaining the loss in metallicity measured in ARPES.

Moreover, a constrain induced by the covalently sp 3 -bonded carbon atoms of the graphene to the silicon atoms of the substrate leads to a distortion of the final carbon layer as shown in theoretical models. [17] Consequently, the calculated C-C distance in the graphene sheet is around 1.55 Å, bearing in mind that periodical calculations also induce an artificial distortion of the graphene network. Our model explains the measured STM images [6,17] reflecting the non-flat behavior of the carbon interface. The Si-C bond distance between the buffer layer and the substrate is equal to 2.14 Å compared with 1.89 Å for the Si-C bond in bulk SiC. Then, extra silicon atoms have been added to the previous interface, and the system has been relaxed again. We have considered here 1.75 ML intercalated Si atoms, following the experimental conditions. Si atoms break the previous Si-C bonds at the interface and link to the substrate Si atom, leading to the release of the carbon atom above, as illustrated in Fig. 3c. Hence, the graphene layer became less highly bonded to the silicon layer below. Calculations including Van der Waals interactions [31,32] show that the graphene layer now lies 2.86 Å above the substrate. This distance is smaller than the interlayer spacing for graphite (3.33-3.35 Å), [33] which explains the measured doping coming from the substrate. We also obtain a C-C distance of 1.55 Å, to be compared with the 8% stretch of the graphene lattice needed to ensure commensurability between graphene and SiC unit cells. [34] The calculated bands along the Γ -K-M direction exhibit the typical linear dispersion of the Dirac cone of a free-standing graphene layer. For the sake of clarity, we only represent here the contribution to the bandstructure coming from the graphene carbon atoms (Fig. 3d). We also observe disperse bands coming from the small hybridization of graphene states with the SiC surface states. The small gap that seems to appear in the bandstructure results from the bandfolding induced by the finite size of the unit cell. This gap does not appear experimentally. The Dirac point appears to be around 250 meV below the Fermi level. Consequently, calculations show, in good agreement with experiments, that silicon intercalation leads to decouple zero-layer graphene lying on top of the silicon-terminated SiC substrate. The calculations confirm that the carbon-terminated surface on top of one ML Si initially semiconducting exhibiting no density of states at Fermi level becomes metallic for higher Si coverage, dispersive electronic states crossing the Fermi Level.

Conclusion

By measuring the electronic properties of the ZLG before and after Si intercalation, we have evidenced the transition from semiconductor to metal of the surface. The LEED pattern revealed that a total decoupling of the ZLG can be achieved at 750 °C, a lower temperature compared with the one ML graphene. This result evidences that the presence of an extra graphene layer on top of the ZLG is an impediment to Si intercalation. The resulting decoupled graphene layer exhibits a slight n-doping revealing an electron density of 4.05 × 10 cm À2 . [12] We have demonstrated that Si intercalation leads to the formation of quasifree-standing one ML graphene. Ab initio calculations corroborate the experimental data.

Figure 1 .

 1 Figure 1. LEED patterns of the 6H-SiC(0001) surface reconstruction showing buffer layer decoupling. (a) LEED pattern of the 6H-SiC(0001) 6√3 × 6√3R30°buffer layer measured at 170 eV. (c) LEED image taken at 200 eV after several cycles of Si deposition followed by annealing at 750 °C. The image presents a 1 × 1 reconstruction showing the loss of coupling between the buffer layer and the substrate. The short red and long green arrows represent respectively the reciprocal lattice vectors of the Bulk SiC B 1 → ; B 2

Figure 2 .

 2 Figure 2. Angle-resolved photoemission spectroscopy spectra presenting the evolution of the electronic properties of the carbon-terminated 6H-SiC (0001) surface during the decoupling process. (a and b) The initial 6H-SiC(0001) 6√3 × 6√3R30°buffer layer measured around K-point along the Γ -K-M direction; the white lines point out extrema of the g1 state modulation in k. (c and d) After several cycles of Si deposition and annealing, we obtain a totally decoupled one monolayer graphene.


Ryan Requist

Silvio Modesti

Pier Paolo Baruselli

Alexander Smogunov

Michele Fabrizio

Erio Tosatti
email: tosatti@sissa.it.

Kondo conductance across the smallest spin 1/2 radical molecule

Molecular contacts are generally poorly conducting because their energy levels tend to lie far from the Fermi energy of the metal contact, necessitating undesirably large gate and bias voltages in molecular electronics applications. Molecular radicals are an exception because their partly filled orbitals undergo Kondo screening, opening the way to electron passage even at zero bias. Whereas that phenomenon has been experimentally demonstrated for several complex organic radicals, quantitative theoretical predictions have not been attempted so far. It is therefore an open question whether and to what extent an ab initio-based theory is able to make accurate predictions for Kondo temperatures and conductance lineshapes. Choosing nitric oxide (NO) as a simple and exemplary spin 1/2 molecular radical, we present calculations based on a combination of density functional theory and numerical renormalization group (DFT+NRG), predicting a zero bias spectral anomaly with a Kondo temperature of 15 K for NO/ Au (111). A scanning tunneling spectroscopy study is subsequently carried out to verify the prediction, and a striking zero bias Kondo anomaly is confirmed, still quite visible at liquid nitrogen temperatures. Comparison shows that the experimental Kondo temperature of about 43 K is larger than the theoretical one, whereas the inverted Fano lineshape implies a strong source of interference not included in the model. These discrepancies are not a surprise, providing in fact an instructive measure of the approximations used in the modeling, which supports and qualifies the viability of the density functional theory and numerical renormalization group approach to the prediction of conductance anomalies in larger molecular radicals.

nanocontacts | Anderson impurity model | ballistic conductance | phase shift E lectron transport through molecules adsorbed on metallic surfaces or suspended between metal leads is the basic ingredient of molecular electronics [START_REF] Aviram | Molecular rectifiers[END_REF][START_REF] Reed | Conductance of a molecular junction[END_REF][START_REF] Joachim | Electronics using hybrid-molecular and mono-molecular devices[END_REF]. Because the highest occupied and lowest unoccupied molecular orbitals rarely lie close to the Fermi energy, electrons must generally tunnel through the molecule, making the zero bias conductance much smaller than G 0 = 2e 2 =h, the conductance quantum, whenever gating is not easily achieved, as is the case in mechanical break junctions and scanning tunneling microscopy (STM). That problem does not persist for molecular radicals, where one or more molecular orbitals are singly occupied, generally resulting in a net spin. When brought into contact with a metal, the radical's spin is Kondo screened (4), leading to a zero bias conductance that may be of order G 0 with a Fano-like anomaly below the Kondo temperature T K [START_REF] Madhavan | Tunneling into a single magnetic atom: Spectroscopic evidence of the kondo resonance[END_REF][START_REF] Li | Kondo scattering observed at a single magnetic impurity[END_REF] and no need for gating. One reason for practical interest in such anomalies is the sensitivity of the conductance to external control parameters such as magnetic fields and mechanical strain [START_REF] Parks | Mechanical control of spin states in spin-1 molecules and the underscreened Kondo effect[END_REF]. Several molecular contacts have been studied [START_REF] Scott | Kondo resonances in molecular devices[END_REF], involving both adsorbed [START_REF] Wahl | Kondo effect of molecular complexes at surfaces: Ligand control of the local spin coupling[END_REF][START_REF] Zhao | Controlling the Kondo effect of an adsorbed magnetic ion through its chemical bonding[END_REF][START_REF] Iancu | Manipulating Kondo temperature via single molecule switching[END_REF][START_REF] Jiang | Real-space imaging of Kondo screening in a two-dimensional O 2 lattice[END_REF][START_REF] Mugarza | Spin coupling and relaxation inside molecule-metal contacts[END_REF][START_REF] Müllegger | Interactions and self-assembly of stable hydrocarbon radicals on a metal support[END_REF][START_REF] Müllegger | Surface-supported hydrocarbon π radicals show Kondo behavior[END_REF][START_REF] Zhang | Temperature and magnetic field dependence of a Kondo system in the weak coupling regime[END_REF] and contact-bridging [START_REF] Parks | Mechanical control of spin states in spin-1 molecules and the underscreened Kondo effect[END_REF][START_REF] Park | Coulomb blockade and the Kondo effect in single-atom transistors[END_REF][START_REF] Parks | Tuning the Kondo effect with a mechanically controllable break junction[END_REF][START_REF] Osorio | Electronic excitations of a single molecule contacted in a three-terminal configuration[END_REF] molecules. None so far involved a radical molecule that is both simple and spin 1/2, and the Kondo anomalies were not predicted from first principles, ahead of experiment. Both the intrinsic complexity of the contact between a large molecular radical and a metal and the unavailability of quantitatively tested ab initio electronic structure-based approaches to Kondo conductance have so far restricted the theoretical work to the role of a posteriori support of STM and break junction zero bias anomaly data. It is therefore important to achieve a first-principles predictive capability of Kondo conductance anomalies across molecular radicals and ascertain its reliability. To that end, we put to work a density functional theory and numerical renormalization group (DFT+NRG) method devised and implemented earlier in our group [START_REF] Lucignano | Kondo conductance in an atomic nanocontact from first principles[END_REF][START_REF] Baruselli | Kondo effect of magnetic impurities in nanotubes[END_REF][START_REF] Baruselli | Kondo effect of magnetic impurities on nanotubes[END_REF] but not yet verified by experiment. Although other strategies with similar goals have also been proposed in refs. 23-26, we can now demonstrate the predictive power as well as the limitations of our procedure in a specific test case where direct comparison with experiment is possible.

Seeking a molecular radical Kondo system of the utmost simplicity and stimulated by the experimentally observed Kondo screening of the S = 1 radical O 2 [START_REF] Jiang | Real-space imaging of Kondo screening in a two-dimensional O 2 lattice[END_REF], we singled out nitric oxide (NO) as the smallest and simplest S = 1=2 radical molecule that might display a Kondo conductance anomaly in an adsorbed state. We have therefore applied our theoretical program to STM conductance of the NO radical adsorbed on the Au(111) surface. NO is an interesting test case, in the past only qualitatively considered as a prototype of radical-surface Kondo interactions [START_REF] Yoshimori | Does a localized spin of adsorbed NO survive on Cu(111)?[END_REF][START_REF] Jigato | The chemisorption of spin polarized NO on Ag{111}[END_REF], whereas recently it has been used instead to quench the spin state of larger molecular radicals such as porphyrins [START_REF] Wäckerlin | Controlling spins in adsorbed molecules by a chemical switch[END_REF].

Starting with DFT calculations of an idealized Au(111)/NO/ STM tip geometry as shown in Fig. 1, we determined the adsorption geometry and the persistence of a nonzero spin state carried by the 2π p orbitals of NO. A two-orbital Anderson model describing their hybridization with the gold surface Fermi sea Significance Molecular electronics, with molecules functioning as the basic building blocks of circuits, is considered by some to be one of the contenders for next-generation electronics. Radical molecules are interesting because a quantum phenomenon, the Kondo effect, involving the radical's unpaired electron and its spin, has the potential to enhance and control conductance. Examples of molecular Kondo are known experimentally in geometries including mechanical break junctions and scanning tunneling spectroscopy (STS), but so far it has not been possible to make accurate theoretical predictions. Here we predict the Kondo effect expected for the simplest molecular radical, nitric oxide, on a gold surface. Subsequently we verify that experimentally by STS. Both agreements and discrepancies offer a quantitative lesson of considerable future use.

was constructed and its parameters were adjusted to quantitatively reproduce the DFT results. The NRG solution of the Anderson model (NRG Results) determined which of the two 2π p orbitals, namely the one with odd symmetry, is eventually screened at the Kondo fixed point. A spectral function peak and a zero bias conductance anomaly were predicted with a Kondo temperature of about 15 K [half width at half maximum (HWHM) of 5 meV]. [Adopting Wilson's definition, the Kondo temperature T K;W is inferred from Δ HWHM ≈ 4:6k B T K;W [START_REF] Zitko | Energy resolution and discretization artifacts in the numeric renormalization group[END_REF], where the HWHM is extracted from the linewidth Γ K obtained from a Frota lineshape fit [START_REF] Frota | Shape of the Kondo resonance[END_REF].]

To check these predictions we carried out STM/scanning tunneling spectroscopy (STS) experiments (Experimental STM/ STS Results) for NO/Au(111). A clear Kondo anomaly was found with an experimental Kondo temperature of about 43 K (HWHM = 16 meV). The observed conductance lineshape near zero bias, a sharp dip, is also quite close but exactly complementary to the calculated Kondo spectral peak. The broad agreement with experiment and especially the discrepancies represent an important critical validation, pointing out physical variables to be included in future calculations.

First-Principles Calculations and Modeling

The initial step is a DFT calculation (Materials and Methods and SI Text) of NO adsorbed on an unreconstructed Au(111) surface, adopted here as an approximation to the face-centered cubic (fcc) domains of the "herringbone" reconstructed Au(111) surface [START_REF] Barth | Scanning tunneling microscopy observations on the reconstructed Au(111) surface: Atomic structure, long-range superstructure, rotational domains, and surface defects[END_REF]. Our calculations (Table S1) show that the only stable NO adsorption geometry is on top, with N pointing down toward a Au surface atom and, remarkably, a tilt angle of nearly 60°as illustrated in Fig. 1. NO also binds (with lower binding energy) at Au (111) bridge and hollow sites. However, in the low coverage limit these adsorption states are stable only when NO is constrained to be upright, for otherwise there is no barrier against migration to the optimal tilted on-top adsorption state. The calculated on-top adsorption energy 320 meV is weak and compares well with an estimated 400 meV from temperature-programmed desorption spectroscopy [START_REF] Mcclure | Adsorption and reaction of nitric oxide with atomic oxygen covered Au(111)[END_REF].

Consistent with that weak adsorption state, spin-polarized DFT shows the persistence of a magnetic moment close to one Bohr magneton, meaning that adsorbed NO retains its radical character as in vacuum. Both 2π p molecular orbital resonances, their degeneracy lifted in the tilted on-top adsorption state, cross the gold Fermi level (Fig. 2). This is an important point-the DFT calculation partitions the single unpaired electron approximately equally between the two 2π p orbitals, despite a significant splitting of the orbital energy levels due to the large tilt angle. Although the evidence of NO tilting in our STM images is not firmly established, its occurrence is convincingly stabilized by the resulting hybridization with the Au surface. Of the total impurity magnetization ∼55% originates from the more hybridized even 2π p orbital with lobes in the Au top -N-O plane (Fig. 2, Inset) and the remaining 45% from the less hybridized odd 2π p orbital normal to that plane.

To calculate the hybridization of the molecular orbitals with the gold surface, we use the information contained in the calculated phase shifts of gold conduction electrons scattering off the adsorbed molecule [START_REF] Lucignano | Kondo conductance in an atomic nanocontact from first principles[END_REF][START_REF] Baruselli | Kondo effect of magnetic impurities in nanotubes[END_REF]. The quantum mechanical scattering problem is solved numerically for the geometry in Fig. 1, using pwcond (34) (Materials and Methods), which provides spinresolved phase shifts for the relevant symmetry channels. Here we focus on two symmetries, namely even (e) and odd (o) under reflection across the Au top -N-O plane, corresponding to the strongly and weakly hybridized NO 2π p orbitals, respectively. Examples of resonances in the phase shifts are shown in Fig. S1. It should be noted for later discussion that the scattering problem involved only the Bloch states of bulk Au and not those of the surface.

The next step of our theoretical protocol is to define an Anderson model Hamiltonian including channels of both symmetries,

H = e e n e + e o n o + X α=e;o X k V kα c † k c α + V * kα c † α c k + H int ; [1]
where the subscript α = e=o stands for the 2π p e=o state, V kα describes hybridization with gold conduction states, and H int contains all of the interactions in the 2π p manifold. When the molecule is brought down to the surface in an upright configuration at the on-top, fcc, or hexagonal close-packed (hcp) sites, the crystal field lowers the cylindrical symmetry of the isolated molecule to C 3v . Further tilting of the molecule at the on-top site breaks C 3v symmetry and lifts the 2π p degeneracy. The splitting of the 2π p e and 2π p o levels (Fig. S2) and their unequal hybridization with gold conduction states, at first sight an irrelevant detail, is actually crucial to the Kondo physics, as we discuss below.

The molecular orbitals hybridize with both surface and bulk states. Surface state hybridization affects the electronic structure of adsorption but only weakly. This is borne out in our slab calculations, where the adsorption energy converges already for 5 layers, whereas the energies of the surface states of the top and bottom of the slab-strongly split by their mutual interactiondo not converge until ∼23 layers (35). If surface state hybridization had a significant effect on the electronic structure, one would expect the adsorption energy to converge more slowly than it does. Surface states decay exponentially inside the bulk, and therefore their contribution to the hybridization and phase shifts, and hence to the conductance anomalies, is invisible in our scattering calculations, which involve only bulk scattering channels. We checked that the ab initio estimates of the molecular orbital hybridization linewidths 2Γ α = 2π P k jV kα j 2 δðe k -e F Þ are well converged with respect to the number of layers of the slab.

The interaction Hamiltonian is now written in terms of the 2π p states as

H int = U e n e↑ n e↓ + U o n o↑ n o↓ + U eo n e n o + J H S e • S o + W ; [2]
where the U terms are intrachannel and interchannel Hubbard interactions, the J H term is their Coulomb ferromagnetic (Hund's rule) exchange interaction, and W = W eo c † e↓ c † e↑ c o↑ c o↓ + W * eo c † o↓ c † o↑ c e↑ c e↓ is a small but nonzero pair hopping term. To fix the model parameters from the first-principles DFT input, we require the scattering phase shifts of the model Hamiltonian at the Hartree-Fock (mean-field) level to reproduce the phase shifts of the ab initio scattering calculation for the full system. The physical picture behind this type of matching is the local moment regime of the Anderson impurity model; because spinpolarized DFT and Hartree-Fock provide comparable mean-field descriptions of the local moment regime, we can reliably infer the model parameters by requiring them to give the same phase shifts. With four total phase shifts from the two channels and two spin polarizations, these conditions determine the four model parameters e e , e o , U e , and U o . To fix the remaining parameters U eo , J H , and W eo , we make use of relationships between the interaction parameters that are exact for the isolated molecule (SI Text) and should still apply to the adsorbed molecule owing to its gentle physisorption. The resulting parameters (Table S2) provide all of the ingredients we need for a many body calculation of the spectral function and Kondo conductance anomaly.

NRG Results

The two-orbital Anderson model, Eq. 1, is solved by NRG [START_REF] Zitko | Energy resolution and discretization artifacts in the numeric renormalization group[END_REF][START_REF] Zitko | NRG Ljubljana[END_REF][START_REF] Bulla | Numerical renormalization group method for quantum impurity systems[END_REF]; details are given in Materials and Methods. For tilted NO adsorbed at the on-top site, the solution indicates a competition between even and odd 2π * orbitals. The 2π p e is much more hybridized with the Au surface, whereas the 2π p o has a lower bare energy, and it is thus unclear at the outset which of the two will host the Kondo resonance. We find, somewhat surprisingly, that at the fixed point the least hybridized, odd orbital is the winner, its spectral function displaying the Kondo resonance shown in Fig. 3E. By occupying the 2π p o more substantially than the 2π p e , NRG inverts the mean-field orbital occupations initially obtained in DFT. At T = 0 K, the NRG orbital occupations are n o = 0:83 and n e = 0:34, whereas the DFT occupations were n o = 0:44 and n e = 0:65. This exemplifies a situation of more general interest, showing that quantum fluctuations may cause one channel, even one predominant in mean field, to give way to another.

As a general note of caution we stress that our Anderson model is in a regime where the parameter dependence of the Kondo temperature, lineshape, and orbital population balance is rather critical. Kondo temperatures are particularly difficult to predict, given their exponential dependence on parameters. At the same time, the accuracy of the DFT electronic structure and phase shifts near e F , which determine those parameters, is limited by inaccuracies and approximations, including in particular self-interaction errors. The Kondo temperature increases by a factor of 3 if the hybridization is increased by 20% or the orbital energy splitting e e -e o is decreased by only 50 meV. Further decreasing e e -e o leads to a rather abrupt transition-via charge transfer from 2π * o to 2π p e -to a state with a broad resonant level at e F rather than a Kondo resonance. We stick here to the results predicted by our straightforward protocol, without adjustments or data fitting.

Experimental STM/STS Results

To verify the theoretical predictions just presented, we carried out STM measurements of NO adsorbed on Au(111). Submonolayer NO was dosed on a reconstructed Au(111) surface at 30 K, and the sample was then cooled to 5 K. STM images acquired with tunneling currents below 10 -11 A show that the molecules adsorb preferentially at the fcc elbows of the gold herringbone surface reconstruction as isolated units or small disordered clusters at low coverage-below about 0.05 monolayers-and form large disordered clusters at higher coverage. The images of small clusters in Fig. 3 B and C indicate that the largest proportion of molecules adsorb on top, in agreement with the theoretical prediction. The average nearest-neighbor NO-NO distance in the clusters is about 0.5 nm. After annealing at 70 K, the molecules form ordered islands with a variety of metastable lattices with unit vectors between 0.35 nm and 0.60 nm. The angle between the unit vectors of the molecular lattices and those of the gold surface is between 0°and 30°. The structures either are incommensurate with the substrate or have a large unit cell with several inequivalent molecules, suggesting that the NO-NO interaction prevails over the NO-Au interaction. The NO lattice shown in Fig. 3D is formed by units that have a dumbbell appearance in the filled-state STM images. This shape is consistent with the charge distribution of the 2π * o orbital of the tilted molecule (Fig. 2) and supports the results of the NRG calculations that predict this state is preferentially filled.

STM spectra measured on isolated NO molecules, small disordered clusters, large disordered clusters, and ordered NO lattices all display similar dip-shaped conductance anomalies centered at the Fermi level, as shown in Fig. 3F. For isolated molecules, whose mobility at 5 K was too high to acquire images with enough resolution to detect the predicted on-top adsorption state and tilt angle, the HWHM of the dip is about 12 ± 4 meV at 5 K with an amplitude that is ∼25% of the background conductance. For NO clusters and lattices, the dip HWHM is 16 ± 4 meV at 5 K and depends on the position of the molecule. In addition to the zero bias dip, shoulders at ±30 meV are visible in the spectra taken over the isolated molecule but not over clusters (Fig. 3F). We have not been able to determine whether these inelastic features are of vibrational or electronic origin, but the symmetric displacement of the shoulders with respect to zero bias, as well as their disappearance in clusters, is suggestive of a vibrational origin. By DFT calculations we find vibrational eigenmodes with frequencies of 25 meV and 42 meV in the correct range. Spin transitions and anisotropy are not expected for a spin 1/2. Hypothetically, electronic states might offer an alternative explanation, the lower Hubbard band of the 2π p o state, visible just below Fermi in the NRG spectra (Fig. 3E), corresponding to the left shoulder, and the empty 2π p e state forming a resonance above Fermi, corresponding to the right shoulder. Altogether, an electronic explanation seems less likely.

No other sharp structures were observed between -0.8 eV and 0.8 eV from the Fermi level as shown in Fig. S3 and described in SI Text. In particular, the large step at -0.45 eV corresponding to the bottom of the clean Au(111) surface state band is absent in the spectra measured over molecules, indicating that free surface state electrons are pushed out by the NO molecule, to become detectable again in spectra 1-2 nm away. Although surface states are repelled, nevertheless their tails are expected to extend up to and hybridize with the NO radical.

To further confirm the Kondo nature of the observed zero bias STM anomaly, the temperature dependence of the half width at half maximum of the dip was measured and is shown in Fig. 4. For each temperature, the lineshape was fitted with an asymmetric version [START_REF] Prüser | Mapping itinerant electrons around Kondo impurities[END_REF] of the Frota function [START_REF] Frota | Shape of the Kondo resonance[END_REF] with asymmetry parameter ϕ = π. By 70 K, the HWHM has grown to 40-50 meV and the amplitude has decreased from ∼25% to ∼10% of the background. Due to the high mobility of the isolate molecules at high temperature and the ensuing difficulty of acquiring spectra for them, all data presented for T > 5 K were taken on NO islands. The effect of the temperature-induced broadening and attenuation of the dip is intrinsic because it is about a factor of 5 larger than the effect of the change of the experimental resolution of STS at high temperature related to the broadening of the Fermi distribution of the tip, corresponding to a HWHM of about 10 meV at 68 K. The large intrinsic broadening rules out inelastic tunneling (vibrational or electronic excitations) as the main origin of the dip. The Kondo temperature ∼43 K, extracted from the low-temperature limit of the Frota linewidth parameter Γ K , is about a factor of 3 larger than our unbiased theoretical prediction. This level of discrepancy is not surprising in view of the exponential sensitivity of the Kondo temperature to the impurity parameters and of the self-interaction errors generally affecting the first-principles DFT calculations.

Actually, a theoretical underestimate of the Kondo temperature is not unexpected, because the hybridization with gold surface states, as was mentioned above, is not completely accounted for in our calculations. Despite these deficiencies, the overall, unbiased theoretical results are gratifyingly reproduced, and so is the evolution of the zero bias anomaly with temperature, as shown in Fig. 4. The experimental offset of the center of the dip from the Fermi energy, less than 2 meV, is much less than the width, indicating that the parent level of the Kondo resonance is indeed nondegenerate. That provides an indirect and yet strong confirmation that the NO molecule is tilted, for if the 2π p levels were degenerate, the Kondo peak would be shifted above the Fermi energy to satisfy the Friedel sum rule [START_REF] Choi | SU(4) Kondo effect in carbon nanotubes[END_REF]. The data do not clearly establish whether the symmetry of the Kondo orbital is odd as predicted, and that remains to be verified, e.g., by photoemission spectroscopy.

Finally, the observed Fano conductance lineshape (Fig. 3E, Inset), a q ' 0 dip rather than q ' ∞ as predicted by the impurity spectral function, indicates that direct tunneling into the molecular 2π p states (t 2 in Fig. 3E, Inset) is dominated by another tunneling channel (t 1 ). As STS indicates, the NO molecule is fully immersed in the Au(111) surface states, whose wave functions extend far out into the vacuum, reaching out to the STM tip. Because surface states do not represent a bulk screening channel, they are not included in our scattering calculations. Because the surface states are nonetheless involved and provide a free electron reservoir for Kondo screening, their local density of states near the molecule will show a dip at the Fermi energy, complementary to the calculated bare NO spectral function peak. As observed in the Kondo STM conductance anomalies of many other adsorbed magnetic atoms, and as generally discussed by Ujsaghy et al. [START_REF] Újsághy | Theory of the fano resonance in the STM tunneling density of states due to a single kondo impurity[END_REF], dI=dV spectra are strongly modified, in our case very plausibly dominated, by direct tunneling into nearby metal surface states, specularly reversing the intrinsic Kondo anomaly.

Conclusions and Outlook

We have described how the Kondo parameters and resulting zero bias STM conductance anomaly of an adsorbed molecular radical can be predicted from first principles, with features that compare reasonably well with those of a subsequent experiment. To do so, we used NO as the simplest spin 1/2 molecular radical capable of forming a Kondo screening cloud when adsorbed on a gold surface. The discrepancies between the calculated Kondo temperature and lineshape, which we did not try to amend in any way, and the measured ones are quite instructive, highlighting in particular the need to fully incorporate the metal surface states in future calculations. The delicate renormalization group flow toward the least hybridized, odd symmetry 2π p o state also appears as an instructive many body effect and remains a prediction to be verified. The magnetic field splitting of the Kondo anomaly, presently academic in view of the large value of T K , could easily be calculated if need be.

The protocol implemented here for NO/Au(111) in a specific STM geometry is of more general applicability and can be applied to different molecular radicals and different STM and break junction geometries, where the influence of structural or mechanical deformations could be explored. As an overall result, the DFT+NRG approach demonstrates enough predictive power to be useful in the a priori evaluation of the conductance characteristics of molecular radical nanocontacts, thus providing a theoretical asset of considerable technological relevance.

Materials and Methods

Density functional theory calculations in the generalized gradient approximation of Perdew et al. [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] were performed with Quantum Espresso (42), a plane wave pseudopotential electronic structure package. The NO/Au(111) adsorption geometry was determined using the slab method in a 3 × 3 surface supercell with a vacuum layer of 16.9 Å. Scattering calculations were carried out with pwcond (34) in a 3 × 3 supercell with a vacuum layer of 8.4 Å, in correspondence with the experimental tip height of ∼10 Å. Brillouin zone integrations were performed on a 6 × 6 × 2 k-point mesh with a smearing width of 0.002-0.010 Rydberg. Plane wave cutoffs were 30 Rydberg for the wave function and 360 Rydberg for the charge density. As detailed in SI Text, Hubbard interactions were applied to the d states of Au (U = 1.5 eV) and the p states of N and O (U = 1.0 eV).

Numerical renormalization group calculations were carried out with NRG Ljubljana (41), using the z-averaging technique [START_REF] Zitko | Energy resolution and discretization artifacts in the numeric renormalization group[END_REF], the full density matrix approach [START_REF] Weichselbaum | Sum-rule conserving spectral functions from the numerical renormalization group[END_REF], and the self-energy trick [START_REF] Bulla | Numerical renormalization group calculations for the self-energy of the impurity Anderson model[END_REF]. The logarithmic discretization parameter was chosen to be Λ = 4, and a maximum of 2,000 states (not counting multiplicities), corresponding to roughly 5,000 total states, were kept at each iteration. Spectral functions were obtained by log-Gaussian broadening [START_REF] Bulla | Finite-temperature numerical renormalization group study of the Mott transition[END_REF] and at finite temperature with a kernel that interpolates between a log-Gaussian at high energy and regular Gaussian at low energy [START_REF] Weichselbaum | Sum-rule conserving spectral functions from the numerical renormalization group[END_REF].

The dI/dV spectra were measured with a clean Au tip by the lock-in technique, applying a 3-meV modulation to the bias voltage and using a maximum current of 3 × 10 -11 A to avoid tip-induced modification of the adsorption geometry. NO molecules were easily displaced at higher tunneling currents at 5 K and were mobile at temperatures above 20-30 K. The gold lattice and molecules could be resolved simultaneously only once (Fig. 3 , where V pp is the peak-to-peak value of the modulation voltage and V is the bias. This function represents the response of the STS spectroscopy with a sinusoidal modulation of the bias and a lock-in amplifier to a delta function-like density of states. The effect of the thermal broadening of the Fermi distribution in the tip at high temperature was taken into account by computing the convolution of the Frota function with the derivative of the Fermi distribution. S2). Experimental error bars indicate the SD calculated from a collection of spectra on different NO molecules. Spectra were measured on at least 10 different molecules at each temperature. The error bars on the theoretical curve represent the uncertainty corresponding to a ±20% uncertainty in the hybridization linewidth; uncertainty from other parameters is not depicted. Model parameter SDs are reported in Table S2.

Fig. 1 .

 1 Fig. 1. Model STM geometry for NO/Au(111). The scattering region shown contains one NO molecule and 28 Au atoms; it satisfies periodic boundary conditions in the transverse direction and is bounded by two planes in the c direction, where the potential is smoothly matched to the potential of bulk gold. NO is shown in the optimal on-top adsorption geometry, where the molecule is bent by about 60°. Figure produced using VESTA (46).

Fig. 2 .

 2 Fig. 2. Projected density of states for NO/Au(111). Spin-up states are plotted as positive values and spin-down states as negative. The densities of states of the NO 2π * molecular orbitals (Insets) are estimated from weighted sums of the projections on N and O p orbitals. The Fermi energy is indicated by a vertical line. The 2π * states of NO and the d z 2 states of the top-site Au atom, Au top , are shifted and scaled for visibility.

Fig. 3 .

 3 Fig. 3. STM/STS images and spectra, all obtained with clean Au tips except for B and C. (A) Isolated NO molecules (arrows) and small NO clusters on Au(111) at 5 K. (B and C) Two small NO clusters formed from a few molecules (in the center) and substrate Au atoms (at the edges) imaged with a molecule on the STM tip. The Au lattice is marked by blue dots. Most of the molecules are on top. (Scale bar, 0.5 nm.) (D) Filled-state STM image ðV bias = -60 mVÞ of an ordered single NO layer obtained after annealing at 70 K. The layer is made up of rows of units with a dumbbell shape that is consistent with the shape of the filled 2π o * orbital. The arrows point to two of these units. (E) NRG spectral functions for the 2π * orbitals of NO at the on-top adsorption site calculated at 5 K with ab initio parameters from Table S2. Insets: spectral function (Left), Fano lineshapes (Upper Right), and interfering tunneling paths t 1 and t 2 (Lower Right). (F) Tunneling spectra of isolated molecules and small clusters at 5 K and clusters at 40 K and 68 K showing the Kondo dip. The spectra are shifted vertically for clarity. The solid lines are fits with the Frota function (31).

  B and C), when an unknown molecule adsorbed on the tip. The horizontal stripes in Fig. 3 B and C are caused by temporary detachments or shifts of this molecule. The broadening of the tunneling spectra caused by the finite modulation voltage is taken into account in the fit of the measured spectra by calculating the convolution of the Frota function with ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 -ðV=V pp Þ 2 q

Fig. 4 .

 4 Fig. 4. Temperature dependence of the HWHM of the Kondo resonance from experiment (red) and NRG calculations (black) based on the ab initio parameters for the on-top adsorption geometry (TableS2). Experimental error bars indicate the SD calculated from a collection of spectra on different NO molecules. Spectra were measured on at least 10 different molecules at each temperature. The error bars on the theoretical curve represent the uncertainty corresponding to a ±20% uncertainty in the hybridization linewidth; uncertainty from other parameters is not depicted. Model parameter SDs are reported in TableS2.

| www.pnas.org/cgi/doi/10.1073/pnas.1322239111 Requist et al.

K. (B and C) Two small NO clusters formed from a few molecules (in the center) and substrate Au atoms (at the edges) imaged with a molecule on the STM tip. The Au lattice is marked by blue dots. Most of the molecules are on top. (Scale bar, 0.5 nm.) (D) Filled-state STM image ðV bias = -60 mVÞ of an ordered single NO layer obtained after annealing at 70 K. The layer is made up of rows of units with a dumbbell shape that is consistent with the shape of the filled 2π o * orbital. The arrows point to two of these units. (E) NRG spectral functions for the 2π * orbitals of NO at the on-top adsorption site calculated at 5 K with ab initio parameters from Table S2. Insets: spectral function (Left), Fano lineshapes (Upper Right), and interfering tunneling paths t 1 and t 2 (Lower Right). (F) Tunneling spectra of isolated molecules and small clusters at 5 K and clusters at 40 K and

K showing the Kondo dip. The spectra are shifted vertically for clarity. The solid lines are fits with the Frota function (31). 72 | www.pnas.org/cgi/doi/10.1073/pnas.1322239111 Requist et al.

| www.pnas.org/cgi/doi/10.1073/pnas.1322239111 Requist et al.

ACKNOWLEDGMENTS. We are grateful for the use of the NRG Ljubljana code and the high-performance computing resources of CINECA. Work was partly sponsored by Contracts Progetti di Ricerca di Interesse Nazionale/ Cofinanziamento 2010LLKJBX 004 and 2010LLKJBX 007, Sinergia CRSII2136287/1, and advances of European Research Council Advanced Grant 320796-MODPHYSFRICT.

Author contributions: R.R., S.M., P.P.B., A.S., M.F., and E.T. designed research; R.R., S.M., and P.P.B. performed research; R.R., S.M., A.S., and E.T. analyzed data; R.R., S.M., and E.T. wrote the paper; and S.M. carried out experiments.

The authors declare no conflict of interest.