Superconductivity in hydrogenated carbon nanostructures

Abstract
We present an application of density functional theory for superconductors to superconductivity in hydrogenated carbon nanotubes and fullerane (hydrogenated fullerene). We show that these systems are chemically similar to graphane (hydrogenated graphene) and like graphane, upon hole doping, develop a strong electron phonon coupling. This could lead to superconducting states with critical temperatures approaching 100 K, however this possibility depends crucially on if and how metallization is achieved.



Introduction
The development of accurate ab-initio many body methods for superconductivity, such asÉliashberg theory [1,2] and density functional theory for superconductors [3], as well as efficient numerical methods to compute electronic interactions [4][5][6] allows prediction of superconducting properties and the critical temperature (T c ) of materials without need of empirical parameters or experimental input. This, combined with the fact that computational approaches are nowadays also able to scan the configuration space of stable materials [7][8][9][10][11] in a faster and efficient fashion (as compared to cumbersome trial-and-error process of synthesis and characterization), has revived hope of finding a, long dreamed, room temperature superconductor.
Another very interesting prediction was that of superconductivity in graphane [27]. Graphane [28][29][30][31] is a fully hydrogenated graphite layer where hydrogen atoms are bonded on the two sides of the layer attached to alternate carbons [32]. Unlike SH 3 , graphane does not superconduct in its stoichiometric structure, it has to be hole doped and this could be in principle achieved by electrostatic gating or by a chemical substitution. So far the doped configuration of graphane has not been synthesized and the prediction is not yet confirmed. Still this material is of great importance for several reasons: first graphane is stable at ambient conditions of pressure, making it of technological relevance; second it has the highest critical temperature ( 100 K) predicted and known at ambient conditions of pressure in phononic superconductors; third graphane, in spite of being an hydrogen compound, does not owe its high-T c to the high energy phonon modes of hydrogen (i.e. Ashcroft mechanism [33]) but rather to the carbon phonon modes and carbon electronic states. This opens the possibility of other high temperature phononic superconductors may exist in a broader set of materials, since at room pressure hydrogen states are usually fully compensated while carbon related states are more likely to be pinned at the Fermi level where they can contribute to superconductivity. In this work we will show that the excellent superconducting properties of graphane also occur in related structures of lower dimensionality such as fullerane, the hydrogenated buckyball (C 60 H 60 ), and hydrogenated carbon nanotubes (CH-NT). While in graphane hydrogenation occurs on both sides of the layer, for these systems we consider only exohedral hydrogenation as shown in Figure 1.
The superconducting state is described by means of density functional theory for superconductors [3] (SCDFT) and an extension used to compute the order parameter in real space [34]. This provides a natural description for large systems with low periodicity as nanotubes and buckyballs and is a strength of SCDFT over the more con-ventionalÉliashberg approach in reciprocal space [35][36][37][38].

Basics of SCDFT
SCDFT is an extension of DFT to account for the very peculiar symmetry breaking that occurs in a superconductor [55,56]. Proposed in 1988 [3] by Oliveira, Gross and Kohn was later revisited and extended in Hardy Gross' group [41,42,57] to merge with the multi-component DFT of Kreibich and Gross [58] to include nuclear motion.
The starting point of SCDFT is the non relativistic Hamiltonian for interacting electrons and nuclei: where e stands for electrons, n for nuclei and ext for external fields.
where ψ are electronic field operators and µ the chemical potential. Nuclei need to be considered explicitly (not just as source of an external potential like in conventional DFT [59]) because in most known superconductors the ion dynamics provides an essential part of the superconducting coupling: where Φ are ionic creation operators, M the mass and Z the atomic number. The Hamiltonian needs to includes an external symmetry breaking field [56] that for singlet superconductivity can be chosen as: In addition to this one should also add an external field coupling with the electronic density: and an external field that couples with the nuclei: In its modern form [41,42] SCDFT is based on the three densities: where 0 is the grand canonical density matrix. The SCDFT generalization of the Hohenberg-Kohn theorem [59] (at finite temperature [60]) states: 1. There is a one-to-one mapping between the set of densities ρ (r), χ (r, r ), There is a variational principle so that it exists a functional Ω that: where ρ 0 , χ 0 , Γ 0 are the ground state densities and Ω 0 the grand canonical potential. The fact that all observables are functionals of the densities and that H is the sum of internal interactions (Eq. (1)) and couplings with external fields (Eq. (5) + Eq. (6) + Eq. (7)) allows Ω [ρ, χ, Γ ] to be written as: that defines the universal functional F [ρ, χ, Γ ].

The Kohn-Sham system
As for conventional DFT it is useful to introduce the Kohn-Sham system [61], a non-interacting system which is minimized by the same densities of the physical (interacting) one, under the three external potentials: ∆s r, r = ∆ ext r, r + ∆xc r, r The subscript H stands for Hartree terms and xc are the exchange-correlation potentials obtained by a functional derivative of the xc functional of the theory F xc [ρ, χ, Γ ], that can be obtained from perturbation theory [41,42].
The electronic part of the Kohn-Sham equations are then derived by diagonalizing the Kohn-Sham Hamiltonian with a Bogoljubov-Valating transformation [56,62]: (14) leading to the diagonalization conditions: that are the electronic Kohn-Sham equation for SCDFT. Their mathematical form is well known in superconductivity literature as Bogoljubov-deGennes (BdG) equations [56] which are mostly used, within the BCS model, to describe superconducting structures in real space. In SCDFT these equations become exact for the calculation of the total energy and the densities: In absence of superconductivity both χ and ∆ are zero and the Kohn-Sham equations (15) and (16) become the usual Kohn-Sham equations of conventional DFT: This form is slightly more general because it would still include the full effect of temperature and ionic motion since it is still coupled with the ion dynamics via the potentials in equation (13).

Transformation to momentum space
Equation (18) can be solved in the superconducting state (i.e. keeping the non-zero χ in the functional v s [ρ, χ, Γ ]) and the corresponding eigenfunctions ϕ nk (r) can be used as a basis set to express the BdG equations in k space. Introducing the expansion: that when inserted into equation (15) and using the orthogonality of the basis set gives: which is a form of the BdG equations particularly useful for introducing approximations. At this stage the problem to solve is still very complicated and can not be tackled without introducing approximations. The key approximation is to decouple as much as possible the many degrees of freedom (and densities) of the problem: 1. Decouple electrons from ions separating static and dynamic part of the interaction, including the latter in a perturbative fashion.
2. Decouple the high energy chemical scale (responsible for bonding) from low energy pairing interactions (responsible for superconductivity).

Phonons and electron-phonon interaction
The problem of correlated electron-nuclear dynamics is enormously complex, however for systems close to equilibrium theoretical methods are able to describe accurately and efficiently the nuclear dynamics and electron-phonon coupling [4,63,64]).
A key approximation is to ignore the effect of superconductivity on the lattice dynamics and on the electronphonon interaction. As the superconducting transition is usually of type II, this approximation is exact near the critical temperature, where the superconducting density becomes infinitesimally small. This allows us to use the lattice dynamics of the normal state.
To compute phonons and the electron-phonon interaction one usually relies on conventional Kohn-Sham density functional theory, defining electron-phonon scattering matrix elements as: where k and q are the electron and phonon momenta, m and n Kohn-Sham band indices, ϕ nk the Kohn-Sham states, ν is the phonon branch, ω qν the phonon frequency and ∆V qν scf the variation in the Kohn-Sham potential due to the ionic displacement corresponding to the phonon mode. By means of density functional perturbation theory [4] these matrix elements can be computed accurately and at a reasonable computational cost for bulk superconductors. The electron-phonon interaction of the Kohn-Sham system reads: where ψ † σnk and ψ σnk are creation and destruction operators for Kohn-Sham states and b νq is a phonon field operator.
The step of approximating the dynamic part of H en withH en can be certainly justified empirically by its success in applications [4,65,66] but is theoretically not very rigorous. The most compelling justification is that if the Kohn-Sham band structure is close to the interacting one so will likely be its response to a lattice motion. Clearly if Kohn-Sham bands are far off from the interacting ones (like in strongly correlated systems) then the use of Kohn-Sham electron-phonon coupling is expected to be a poor approximation.

Band decoupling approximation
The electronic BdG Kohn-Sham equation (22) can be further simplified by assuming that the superconducting condensation is a small perturbation on the nonsuperconducting system. As already pointed out in the previous section, since the superconducting transition is usually of type II the assumption becomes exact close to T c therefore it would not affect the estimation of T c itself.
This assumption implies that the superconducting transition will not induce a structural one, therefore ∆ s (r, r) should keep the original lattice periodicity and the quantum number k in equation (18) must be maintained [57,67]. In other words the summations in equations (19) and (20) should only run over the band index n and not over k.
The summation over n means that the superconducting transition can still induce hybridization between different bands corresponding to the same k-point. However unless bands are degenerate (or close to degeneracy with respect to the energy scale set by ∆ s that is of the order 10 meV) this hybridization must be extremely small. Therefore, apart from anomalous cases, one can introduce a second and stronger approximation by ignoring this superconductivity induced band hybridization effect. Equations (19) and (20) reduce to: that implies ∆ s,nn kk → δ nk,n k ∆ s,nk . Inserting equation (25) into equation (22) one can formally solve these equations obtaining: with e φ nk = ∆ s (nk) / |∆ s (nk)| and E nk = ± ξ 2 nk + |∆ s (nk)| 2 . While the densities in equations (16) and (17) take on the simple form: The entire superconducting problem is now reduced to the construction of the matrix elements of the Kohn-Sham potential ∆ s (nk) that are obtained by the solution of the equation: Several approximations fro F xc have been proposed and tested [42,68], and all lead to a BCS-like form of the equation above: where the two kernels K and Z depend on the chosen F xc functional [42,[68][69][70][71] and contain all the key information about electronic and electron-phonon coupling.
The main point of strength of the theory is its low computational cost, as compared to Green's function methods, mainly because the gap equation (31), while being fully dynamical and including strong coupling effects, does not involve cumbersome Matsubara integration (all the complexity is absorbed in the process of functional construction). This means that the equation can be easily solved in its full k resolution and in a full energy window. Then all anisotropy effects [38,47,52,54,72] as well as high energy Coulomb interactions [70,71,73,74] can be included in a fully ab-initio manner, even when pairing is induced by exchange of spin-fluctuations [75,76].

Superconductivity of a doped carbon-hydrogen nanotube and doped fullerane
The lattice properties of the two chosen hydrogenated carbon nano-structures are collected in Figure 1. These are an hydrogenated (6,0) nanotube (CH-NT) and an hydrogenated fullerene (C 60 H 60 ). Both are considered with a full outer hydrogenation of the carbons. We have established that these structures are dynamically stable but we will not consider their thermodynamic stability [77].
Both CH-NT and C 60 H 60 ( Fig. 1) are insulating structures. The KS-GGA band gap is 2 eV and 3 eV respectively.
In both systems the top of the valence band is characterized by C-C an C-H bonding states, while the conduction band is delocalized and dispersive (Figs. 2 and 3).
The simplest way to simulate a doping is by rigidly shifting the Fermi level. Superconductivity is expected to be stronger upon hole doping as holes should form a metallic state in the stiff C-C covalent bonds [27,78]. Owing to the nature of the conduction band, electron doping would lead to far weaker superconducting coupling.
The lattice dynamics is characterized by high energy H stretching modes at about 380 meV and a continuous region below 200 meV that contains at its high end the H rocking modes and in its mid to low region carbon and acoustic modes (Figs. 2c and 3e).
The electron phonon coupling important for superconductivity (discussed in Sect. 2.1.2) can be reduced to thé Eliahberg α 2 F (ω) function by averaging over the Fermi surface: with k ≡ k + q, and N F the DOS at the Fermi level. This function is shown for the two systems and at a single doping level in Figures 2d and 3e.

Superconductivity of a doped CH-NT
We will first focus on superconductivity in the CH-NT as it results the most interesting and likely to occur. Figure 4 gives the full doping dependence of T c , coupling strength and α 2 F for this CH-NT.
T c is computed within SCDFT as described in Sections 2 and 4. In spite of the strong Coulomb interaction T c appears to be quite high ranging from 50 to 100 K depending on the doping regime. There are clearly two doping regimes that can be identified in Figure 4 a low doping corresponding to about 1 e − per unit cell that leads to a first peak in the T c curve, with a maximum T c of about 50 K and a heavy doping regime with T c up to 100 K that corresponds to about 3e − /cell. In the low doping regime the α 2 F (ω) function upper panel in Figure 4 shows three main structures at about 65 meV, 150 meV and a high energy structure above 350 meV. In the high doping regime there is also a prominent structure around 150 meV and an high energy structure at 380 meV, while the low energy structure is reduced in intensity and energy. Overall this high doping region shows both a larger coupling and a larger average effective phonon frequency.
While T c is of the same order of magnitude in these two cases there is a change in the properties of the condensate in these two regimes. The local structure of the superconducting order parameter χ and xc-potential is plotted in Figure 5. As discussed in reference [34], the xc-potential mostly gives information on the Coulomb renormalization in its full energy range (that weakly depends on the doping level), while χ (R, s = 0) shows where the attractive part of the coupling is located. At low doping this is concentrated in the C-H bonding region (Figs. 5b and 5e) while at high doping the attraction is more concentrate in the C-C bonding region (Figs. 5c and 5f), similar to what occurs in graphane [34]. In graphane the contribution from the C-H bond on the attractive channel is negligible, here it is still relevant even at this high doping regime. This means that the CH-NT is, to some extent, a hydrogenic superconductor [25,33,[79][80][81].
In graphane this rigid shift approach to doping mimics well the electronic effect of a B substitutional defect [27]. For the CH-NT we observe this approximation to be less accurate than in graphane as, upon B substitution, the top of the valence band (now at the Fermi level) acquires a strong B character that affects both position and shape of the 1D Van Hove singularity in the density of states. A more accurate approach would be to correctly account for electron-phonon coupling and T c for the C → B substituted system. However, this leads to a phonon softening  at finite q along the tube and could be an evidence of a charge-density wave instability. As the simulation of the charge density wave is computationally too demanding (a super cell along the instability would require over a hundred of atoms) so we opt for the rigid shift of the Fermi level.
A rigid shift approach also neglects the effect of phonon renormalization induced by the metallic screening. Phonon softening always leads to an increased electron phonon coupling, however it also reduces the energy window around the Fermi level where phonons can pair cooper pairs. The overall effect strongly depends on the system, for example in MgB 2 the Kohn anomaly is crucial [78,82,83] to lead to a correct T c .
In the present case the error has been estimated to be small. This estimation is done by comparing T c computed from the coupling at q = Γ (where phonons are stable) for the case of a doping by rigid shift and by B substitution. The former approach gives. 1 85 K the latter 95 K. An extremely small difference indicating that the rigid shift approach is sufficiently reliable

Superconductivity of a doped fullerane
We adopt the same approximations discussed for the CH-NT. Like in the CH-NT, doping in C 60 H 60 can be simulated by a rigid shift of the Fermi level, leading to T c predictions of the same order of magnitude as in the CH-NT. However the molecular C 60 H 60 units in the crystal are weakly interacting, as seen by the fact that the band dispersion is quite low at the top of the valence band. This raises doubts about the validity of several approximations in our superconductivity approach such as the decoupling approximation and Migdal's theorem [84]. A related problem is that rigid shift predictions are not confirmed by adopting a better description for the doping. This can be clearly seen comparing Figures 3b and 3c. Jellium doping causes a major transformation on the top of the valence band while C → B substitutions introduce localized states in the band gap that pin the Fermi level. Such states are not likely to lead to any superconductivity as their localization would correspond to an extremely strong pair-breaking Coulomb repulsion.
In short this means that the C 60 H 60 -crystal could possibly superconduct but this would depend crucially on how the doping is achieved and especially if doping is able to induce a stronger interaction between the molecules. This is certainly not obvious for hole doping. On the other hand our rigid shift calculation indicates that superconducting coupling would be negligible for electron doping.
Eventhough this rigid shift T c 50 K may not be very reliable, it is interesting to observe that also in this system the doped CH bonds hide a significant pairing strength. If not in the present geometry, that has a high formation energy [85], this may lead to superconductivity in related nanostructures [77]: larger fulleranes, hydrogenated buckyonions, systems with lower hydrogen content, different hydrogenation geometry [85] or under applied pressure [86].

Computational methods and details
All structural relaxation, normal state electronic structure, phonons and electron phonon calculations have been done with density functional theory [59,61] using the PBE [87] approximation for the exchange correlation functional as implemented in the ESPRESSO package [88]. Core states are described in the norm conserving [89,90] pseudopotential approximation. Lattice dynamics are evaluated with density functional perturbation theory [4,66,91]. A 100 Ry cutoff is used for the planewave expansion of the charge density. The brillouin zone integration for the electronic charge is done in a 2×2×20 (2×2×2) k-grid for CH-NT (C 60 H 60 ), phonons are computed in a 1×1×10 (1×1×1) q-grid. The k integration for theÉliashberg function and the solution of the SCDFT gap-equation (31) requires a very accurate sampling especially of those states close to the Fermi level, because the kernels K and Z are sharply peaked around ξ = 0. Therefore an accurate interpolation scheme is necessary. As described in references [42,47,54] a convenient approach is to use a large set of random k-points accumulated around the Fermi level and opportunely weighted. The properties of the corresponding states are then obtained by interpolation from calculations on the regular grids.
Superconducting T c is computed within SCDFT by adopting an isotropic approximation. Meaning that the coupling is first averaged (in the form of Eq. (32)). This neglects multiband effects that give a negligible contribution to the estimation of T c . A more relevant approximation, that will not be validated in this work, is to assume a conventional strength for the Coulomb interaction. In SCDFT Coulomb repulsion in the Cooper pair is included in a completely parameter free way by computing numerically the matrix elements of the screened interaction (usually in RPA) that enter the xc functional [41,42] both in a static [41,42] and a dynamical way [70,71]. However this is computationally quite expensive and we will leave this analysis to a further investigation. Here we will assume Coulomb matrix elements to be constant for all states. Their value V chosen such that µ = V*N F is equal to 0.6. This is similar to what is used in theÉliashberg-McMillan [1] approach to superconductivity. The only difference is that in the present case Coulomb interaction will still depend on energy via the density of states (while in the Morel-Anderson approach the density of states is also assumed to be constant). Therefore Coulomb renormalization occurs in the solution of the gap equation (31) and not analytically [92,93]. The chosen value of µ is relatively high, as expected for a system with a poor screening as doped insulator. Still a more realistic approximation may give a stronger Coulomb repulsion especially in the molecular C 60 H 60 structure and could reduce the estimated T c .

Summary and conclusion
We have constructed two hypothetical systems: an hydrogenated carbon nanotube and fullerane, the hydrogenated fullerene. These prove to be dynamically stable and are electronically similar to graphane. To characterize the superconducting state of such complex structures we have used density functional theory for superconductors and constructed the order parameter and the Kohn-Sham gap function in real space. Upon doping the valence bands of these systems show an exceptionally large electron phonon coupling, indicating that the prediction of a high critical temperature in graphane may be a general property of a large class of doped hydrogen carbon systems at room pressure. In particular this class certainly includes larger (and more stable) hydrogenated nanotubes [94]. Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.