Electronic and superconducting properties of hydrogenated graphene from first-principles calculations

We performed first-principles calculations on two hydrogenated graphene systems with different hydrogen coverages, C8H2 and C50H2, to analyze their electronic and superconducting properties. Our results show that their electronic properties are highly correlated to the hydrogenation positions. If the two hydrogen atoms are attached to the same sublattice, the final system will be ferromagnetic. Otherwise, it will maintain nonmagnetic rather than anti-ferromagnetic. Moreover, the distance between the doped hydrogens can trigger the movement of Dirac points, and even annihilate Dirac points when the distance is close to the maximum. We further studied their superconducting properties by applying hole doping and tensile strains. The results show that the superconducting transition temperature T c increases with more holes and reaches its maximum of about 20.2 K at the critical doping level (x c = 0.17 holes/cell). Our results show that the superconductivity mainly originates from the coupling between the out-of-plane lattice vibration modes and the electronic p z orbitals of carbon atoms. The increase of T c can be attributed to the stronger coupling between the electrons and the low-frequency phonon. However, the application of biaxial and uniaxial tensile strain will depress the superconductivity because of the modulation of the low-frequency phonon. It is worthy to note that weak anharmonicity exists in the hydrogenated graphene systems. This work provides a systematic study on tuning the superconductivity of hydrogenated graphene.


Introduction
Superconductivity in layer compounds has come into focus in recent years due to their intriguing properties, for instance, the single unit cell FeSe, few-layers NbSe 2 , and twisted graphene, etc [1][2][3]. In addition, some elemental superconductors, T-graphene and 2D boron, has been predicted under high pressure [4,5]. Graphene, as a typical two-dimensional (2D) material, has become a hot title for its promising properties [6][7][8][9][10][11] since it was mechanically exfoliated from graphite in 2004 [12]. Although it is a zero-bandgap [13] semimetal that is not positive for superconductivity [14][15][16][17], it can still be tuned into a superconductor via certain methods. For example, moderate charge doping could modify the band structures by shifting the Fermi level away from the Dirac point, making graphene a superconductor [18,19]. Further studies also predicted that the superconducting transition temperature T c of hole-doped graphene could be raised to 30 K by applying extra tensile strain [20]. The intercalation of metal atoms in graphene can lead to charge transfer, thus can also be used to shift the Fermi level and induce superconductivity [21,22].
Hydrogenation is another way to modify the electronic properties of graphene. Different results have been achieved by changing the doping amount or the adsorption position of hydrogen atoms [23][24][25][26]. For example, fully hydrogenated graphene shows a wide bandgap with annihilated Dirac points [27,28], while half-hydrogenated graphene shows anti-ferromagnetism (AFM) instead of ferromagnetism (FM) [29,30]. Moreover, it has been discussed how the band-gap opening changes under different hydrogen coverages [31]. Recently, the coexistence of FM and superconductivity has been theoretically predicted in C 6 H 1 and C 6 H 5 systems [32]. This indicates that superconductivity might exist in hydrogenated graphene. However, the coexistence of FM and superconductivity does not obey the Bardeen-Cooper-Schrieffer (BCS) theory, which applies to conventional superconductors only [33,34]. What's more, systematic studies about the influence to electronic properties from the different hydrogen coverages and the distance between the adsorption hydrogen atoms are still missing, which can only be possibly achieved in calculations with large supercells.
In this article, we perform first-principles calculations to study the electronic and superconducting properties of graphene doped with two hydrogen atoms. The reason why we choose two hydrogen atoms is to include doping on two different sites of graphene while decreasing the number of configurations for convenience. Here we only consider the situations in which the two hydrogenated atoms are attached to different sides of the graphene monolayer since they usually have lower energy [35]. To study how hydrogenation positions affect the electronic properties, we adopt 2 × 2 × 1 and 5 × 5 × 1 supercells of graphene, which include three and eleven inequivalent doping methods, respectively. The results show that superconductivity cannot appear after hydrogenation. With moderate hole doping, however, superconductivity emerges and can be further engineered by tensile strain. Due to the presence of hydrogen atoms, we further analyze the anharmonicity in this system. To determine the mechanism of electron-phonon coupling, we also study the lattice vibration modes in the strong coupling zone, together with the electronic band structures and density of states (DOS).

Calculation details
In this article, we calculate the electronic properties of C 8 H 2 and C 50 H 2 via the Vienna ab initio simulation package (VASP) [36] and the superconducting properties of C 8 H 2 via the Quantum ESPRESSO (QE) [37] program separately using the plane-wave basis and the projected augmented wave (PAW) method [38]. We adopt the pseudopotential with the generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) [39] for exchange-correlation potential throughout the calculations. The phonon spectrum is calculated using density functional perturbation theory (DFPT) [40] as implemented in QE. The criteria for the convergence of total energy and force are 10 −7 eV (10 −7 Ry) and 10 −4 eV Å −1 (10 −6 Ry Bohr -1 ) in VASP (QE). In VASP, the cutoff energy for the plane-wave basis set is 600 eV, and the k-point mesh is 12 × 12 × 1 (4 × 4 × 1) in structure optimization and 48 × 48 × 1 (48 × 48 × 1) in the DOS calculation for C 8 H 2 (C 50 H 2 ). We choose the same mesh for the DOS calculation to ensure that the selected k-points are dense enough to locate the Dirac points at the Fermi level. In QE, the energy and charge density cutoffs are 80 Ry and 480 Ry, respectively, and the k-point mesh for structure optimization is 12 × 12 × 1. The superconducting calculation requires a denser k-point mesh to contain all the k-points and q-points. We take meshes of 4 × 4 × 1 and 48 × 48 × 1 for the q-point and denser k-point respectively.
The tensile strain ε is described by ε = l−l 0 l 0 × 100%, where l 0 and l are the lattice constants before and after the strain is added. To relax hole-doped structures with tensile strain, we first add holes into the pure system and relax both lattice and atom positions by decreasing the valance electrons. To deal with the divergence of Coulomb potential in Fourier transform, a compensating uniform background is applied as implemented in VASP and QE. Doping level x describes the number of holes added into the unit cell. After fully relaxing the structure with a certain x, we add ε on both sides of the lattice and relax again. This time, we only relax the atom positions. The final structures will be used in further calculations.
QE uses the modified Allen-Dynes formula [41] to calculate T c , which can be written as: where ω log , λ and μ * represent the logarithmic average of the phonon frequency, the electron-phonon coupling parameter and retarded Coulomb pseudopotential, respectively. We choose μ * = 0.1 as regular. λ is calculated throughout the Brillouin zone as [41,42]: is the Eliashberg spectral function calculated as: where N q is the number of q-points in generating the phonon mesh and N(E f ) is the density of states at the Fermi level. ω q j and γ q j are the phonon frequency and phonon linewidth of a certain phonon mode j with momentum q, respectively. ω log is computed by this formula: Phonon anharmonicity can be computed using a newly released Python package stochastic self-consistent harmonic approximation (SSCHA) developed by Monacelli et al [43]. The supercell of doped C 8 H 2 is enlarged by 4 × 4 × 1 as the q-point mesh. Five hundred configurations are generated during each self-consistent cycle by moving the atoms according to the modes from the trial or updated dynamical matrices. The newly calculated dynamical matrices, once they achieve convergence, can be read by QE for a new electron-phonon calculation.

Electronic properties of C 8 H 2
To construct hydrogenated graphene, we attached hydrogen atoms vertically to different carbon atoms in the supercell of graphene. The initial distance between the carbon and hydrogen atoms is set to 1.11 Å, which is identical to the distance within fully hydrogenated graphene [27]. To preserve the 2D character, the simulating vacuum space between two graphene layers is set to 20 Å along the c direction. Considering the lattice periodicity, some configurations are equivalent to each other, so we can only find three different arrangements for two hydrogen atoms in a 2 × 2 × 1 supercell. There are two carbon atoms in the original supercell that will be attached to two hydrogen atoms, and we call the distance between these two carbon atoms the doping distance in this article, which is labeled d. The C-C bond length in the original graphene is labeled a, which equals 1.42 Å. Thus, we can mark these three different configurations with d = a, d = √ 3a and d = 2a, as shown in figures 1(a)-(c). The relaxed C-H bond lengths are 1.114, 1.131, and 1.116 Å. Figure 1(d) shows the first Brillouin zone used in the following calculations. It is similar to the hexagonal shape in graphene but does not preserve the perfect sixfold rotation symmetry. As a result, the original K and M points located in the corner and middle of the edge, respectively, become ordinary points. Nevertheless, we still choose one of the triangle paths, for example, K 1 → Γ → M 1 → K 1 for band structure calculations, which is accurate enough to explain some phenomena. However, there should be inevitable conflicts between the band structures and DOS calculations in some cases. Here, we emphasize that more accurate results rely on DOS instead of band structures.
The doping distance, lattice parameter, energy difference compared to lowest energy state, and space group of C 8 H 2 are listed in table 1. When d = a, the two hydrogen atoms are attached to the adjacent two carbon atoms, one of them (labeled C 1 ) bonds with one hydrogen atom and forms a sp 3 -hybridized structure. This carbon atom rises slightly from the original plane, making the nearest carbon atom C 2 naturally lower, which is positive for attracting the other hydrogen atom. As a result, this system has the lowest energy. Moreover, the C 2 atom, when attached by hydrogen atoms, does not show the same hydrogenation progress as the C 1 atom because the C 1 -C 2 bond rotates slightly, making the C-H bonds in the final C 8 H 2 tilt, which is different from other cases. Additionally, the tilt of C-H bonds depresses the expansion of lattice parameters. The original pure graphene supercell has a 1 = a 2 = 4.936 Å, and the configuration with d = a has a similar value, as shown in table 1. For the other two configurations, however, the lattice constants and volume expand. Among these three different configurations, only case d = √ 3a shows FM with a total magnetic moment of 2 μ B /cell, and the other two cases show nonmagnetism (NM). This can be identified from the band structure and DOS calculations. As we can see in figure 2(a), the calculated band structure without spin polarization exhibits a flat band near the Fermi level, which causes a giant peak in DOS near the Fermi level and indicates Stoner instability [44]. To lower the total energy, the bands prefer to split into spin-up and spin-down parts and open a bandgap, thus making the system FM. The bonding modes of p z orbitals in carbon atoms are responsible for the formation of the FM state. In the original graphene, the p z orbitals are bonded into delocalized π bonds due to sp 2 hybridization. By adding one hydrogen atom, the attached carbon atom changes into sp 3 hybridization, and the specific p z orbital is bonded by a σ bond with the hydrogen atom. This breaks the symmetry of π bonds and results in the giant peak of DOS near the Fermi level, which changes this system into a FM state. When adding the other hydrogen atom, whether the final system maintains FM depends on the position of the attached carbon atom, as analyzed later in section 4.
The Dirac points in perfect graphene appear at the K and K points. Although these two Dirac points are maintained at the Fermi level in the case with d = a, the hydrogenation changes the point group from D 6h to C 2h , thus displacing the Dirac point from the K points to drift along the K 1 → Γ path, as shown in figure 2(c). Moreover, the Dirac cones show little anisotropy, as in other studies, which was attributed to different hopping parameters in the tight-binding (TB) model [35]. In section 4, we will further discuss the movement of Dirac points, which depends on the doping positions and directly determine the space group of these configurations.
Bandgap opening also occurs in the case d = 2a. These two hydrogen atoms are placed diagonally in a C 6 ring and are the farthest in these three doping configurations. This structure has an indirect bandgap of approximately 3 eV, which is a wide bandgap insulator. The transition from direct zero-bandgap semimetal to indirect wide-bandgap insulator makes this configuration much more different from perfect graphene. In figure 2(d), the DOS has two giant peaks located within the valance and conduction band zones. Compared to the case d = a, the distance of hydrogenation doping positions increases, and the DOS exhibits more localization on both sides of the Fermi level, thus opening a wide bandgap.

Electronic properties of C 50 H 2
In the section before, we discuss the situation of attaching two hydrogen atoms to the 2 × 2 × 1 supercell of graphene and identify three possible configurations of C 8 H 2 . We analyze their electronic properties and obtain some interesting results related to the doping distance d. In C 8 H 2 , the lowest energy state d = a maintains Dirac points and drifts slightly. When d increases to 2a, the Dirac points disappear, and a wide bandgap emerges. Moreover, only the configuration with d = √ 3a shows FM. We decide to enlarge the graphene supercell to 5 × 5 × 1 for further study of how doping distance d influences the final structure. C 50 H 2 possesses eleven different d values from d = a to d = 5a, as shown in table 2.
Among these eleven configurations, four show FM, while the rest show NM. This originates from hydrogenation on different sublattices. We can label these two sublattices as A and B, respectively. We first The red (blue) lines in (b) represent bands occupied by spin-up (down) electrons. keep one carbon atom fixed at sublattice A (or B, as the same) and then label the other atom. As shown in figure 3(a), the blue and red positions indicate that the final system exhibits FM and NM, respectively. The arrangement shows that if two hydrogen atoms attach to the same sublattice, the final configuration will present FM; otherwise, NM. This can be described with Lieb's theorem [45]: graphene has two identical sublattices, and once the hydrogen atoms are doped in, the total spin moment of the ground state is given by , hydrogenation occurs on carbon atoms from the same sublattice, which results in a FM structure with a total spin moment of 2 μ B /cell, while in other cases, the attached carbon atoms are from different sublattices, thus exhibiting NM (or AFM, which is not observed through our calculations and will be discussed in following paragraph). This phenomenon is also predicted and observed by other calculations and experiments [46][47][48][49]. The spin density is also calculated as shown in figure 3(b). The attached carbon atoms have little contribution to the total spin, while the nearest carbon  atoms and the hydrogen atoms have a significant contribution to the spin density. The contribution from the p z orbital of carbon and the s orbital of hydrogen makes the final system FM.
Previous work [50,51] discussed the existence of AFM in hydrogenated graphene and they draw different conclusions on the same configuration, as referred in this article, d = √ 7a. Boukhvalov et al [50] proved that it is NM while Ferro et al [51] stated that it is AFM. However, they all showed that the formation of π bond is crucial to determine the final magnetic state. Here, we performed the spin polarized calculations to include AFM situations and found out no local spin moments survive in all the S = 0 configurations. To check the π bond formation, we calculate charge density differences of d = √ 7a, d = √ 13a and d = 5a. As shown in figure 4, the formation of π bond still occurs even at the largest distance, which is contrary to reference [51]. Although the π bond is weak at d = 5a, it is still significantly meaningful to show that electrons can hop from the attached graphene hexagon to a farther hexagon, which is contrary to reference [51] where electrons can only hop inside the attached hexagon. In addition, as shown in table 2, nonmagnetic configurations d = √ 7a and d = √ 13a are naturally more stable than FM configurations d = 3a and d = √ 12a in our calculated total energy, which is different to the results in reference [51]. Furthermore, figures 3 and S12 in reference [47] also showed that, even at a larger doping distance, the S = 0 configurations are still NM and have lower energy than FM configurations. Therefore, our results support the conclusion of reference [50] that the d = √ 7a configuration is NM. Nevertheless, our results can not exclude that with more hydrogen atoms or larger distance than reference [47] AFM might appear, which is beyond the scope of this work.
We further performed the band structure and DOS calculations, as shown in figure 5. Four FM structures always have a flat band across the Fermi level under non-spin polarized calculations and experience spin splitting under spin-polarized calculations, as shown in figures 5(a) and (b) with d = √ 3a as an example. The NM structures present more specific rules. The lowest energy state d = a in figure 5(c) is also a zero-bandgap semimetal with Dirac points. When the d value enlarges, the system still keeps the Dirac points at the Fermi level until d = √ 19a it becomes an insulator with a bandgap, as shown in figure 6(a). This reveals that the annihilation of Dirac points and the creation of bandgap are highly related to the maximum of d. The reason why Dirac points change into bandgap also resembles C 8 H 2 . When d = a, the DOS equals zero at Fermi level. As d increases near the maximum, the peaks of DOS near both sides of the Fermi level will be highly localized, thus opening the gap and triggering the topological transition. From figure 6(b), we also know that the maximum bandgap in C 50 H 2 is approximately 1 eV compared to 3 eV in C 8 H 2 . This means that as the graphene supercell enlarges, the maximum bandgap will decrease. In addition, the zero bandgap cases show an interesting movement of Dirac points. For the cases d = a, d = 2a and d = 4a in which the space groups are C2/m, the Dirac points always appear along the calculated k-point path, as shown in figures 5(c) and (d). For the other two cases, their space groups are P-1, and the Dirac points are still maintained, as shown in figures 6(c) and (d), but cannot be observed along all the triangle paths, which means they are located inside the triangle zone surrounded by the calculated path. We also choose the path K 1 → Γ → M 1 → K 1 for comparison. In figures 6(c) and (d), the DOS reaches the zero point at the Fermi level, although the band structures still have a bandgap. It can be explained detailly by TB theory. The space groups are directly decided by the hydrogenation positions, which will further affect the hopping parameter of carbon atoms judged by virtual crystal approximation (VCA) [52,53]. In VCA, the doped supercell of graphene can be transformed into a graphene-like two-atom primitive cell with three  [54,55]: where D x,± and D y,± are the x and y components of these two Dirac points and t i (i = 1, 2, 3) represents the three nearest-next hopping parameters. It should be mentioned that we can set the connection line of any pair of Dirac points as the x-axis, such as K i K i (i = 1, 2, 3). For the three C2/m cases, we have t 1 = t 2 = t 3 ; thus, the Dirac points drift along the x-axis. The difference is the movement direction. In the original hexagonal cell, there are six identical positions. When the space group transfers to C2/m, to maintain the periodicity and central symmetry, one pair will move inside the cell and the other two pairs will move outside, and vice versa. As shown in figure 7, there are two types of Dirac point movements in these three cases. For the case d = a and d = 4a, one pair moves in along K 1 K 1 and the other two pairs move out along K 2 K 3 and K 2 K 3 , which are parallel to K 1 K 1 . For the case d = 2a, K 2 K 2 becomes the x-axis, and this pair moves out, leaving other pairs moving in on the edge. The two P-1 cases with t 1 = t 2 = t 3 , however, will have two Dirac points inside the area surrounded by the calculated path, which will not be detected in the band structure calculation.

Superconducting properties of C 8 H 2 with doping or tensile strain
None of the different structures of C 8 H 2 or C 50 H 2 show any evidence for conventional superconductivity judging from the band structures and DOS because they are not metallic. The FM states do not have C 3v symmetry and cannot be identical to the unconventional superconductivity described in previous research [32]. Whether the coexistence of FM states and superconductivity could occur in this case needs to be further studied. Here, we focus on the NM configurations, which are either a zero-bandgap semimetal or a wide bandgap insulator. This indicates that superconductivity cannot appear in such pure systems. However, it is still hopeful to obtain superconductivity in doped systems. We focus on C 8 H 2 next for convenience of calculation. For the case d = a, it maintains two Dirac points as the original graphene but in different positions of the Brillouin zone compared with pure graphene. Previous research shows that graphene can be doped into a superconductor, and the transition temperature T c will rise if tensile strain is added [20]. Then, we infer that a similar phenomenon can also happen in this case. However, two reasons prevent us from calculating further. First, the point group is C 2h , which is not convenient for comparison with the original graphene. Second, the maximum N(E f ) doping is rather small and not promising. Naturally, we choose the case d = 2a, of which the point group is D 3d and the DOS near valance band maximum (VBM) is highly localized. Here, we only focus on hole doping because the system is dynamically unstable even when doping in a few electrons. By doping moderate holes, the Fermi level shifts downward, and N(E f ) will change rapidly. We ignore the subscript of high symmetry points in the following calculations because each triangle path should be the same under D 3d symmetry. First, it is necessary to check out the dynamic stability of the case d = 2a. Figure 8(a) shows the phonon dispersion curves and phonon density of states (PHDOS) of this structure. No imaginary frequency emerges, which indicates that the pure system is dynamically stable. We also see that there are two nearly identical high-frequency branches of approximately 2800 cm −1 contributed by the vibration of hydrogen atoms with weak dispersion. Then, we start to dope holes into this system. As heavy doping always makes the systems unstable, we test the critical doping level x c in this system, and the results show that x c = 0.17.
We further calculate electron-phonon properties when doping levels x = 0.05, 0.10, 0.15 and 0.17. The results show that the transition temperature could rise to 17.3, 18.9, 19.9 and 20.2 K with μ * = 0.10, respectively. The corresponding N(E f ), λ, ω log and T c are listed in table 3. The N(E f ) value does not change too much due to the smooth trend of the DOS near the first singularity. Figures 8(b) and (c) shows the phonon spectrum, PHDOS and Eliashberg spectrum function α 2 F(ω) of x = 0.05 and 0.15. As the doping level x increases, the final configuration has a larger λ and smaller ω log . However, the increase in λ takes advantage of the decrease in ω log ; thus, the final T c increases. We can also learn from the dimensionless factor γ qj ω qj that electron-phonon coupling concentrates on three parts: the low-frequency zone from 0 to 500 cm −1 , the medium-frequency zone from 1300 to 1500 cm −1 and two high-frequency branches. The contribution of the medium-frequency zone and two high-frequency branches is rather small due to the inverse relation. Therefore, the low-frequency zone contributes significantly when inducing superconductivity, especially along the path M → K. When more holes are doped, γ qj ω qj will become larger, triggering a stronger coupling, most of which is from the low-frequency zone. This can explain a larger λ and smaller ω log .
In the following, we discuss how tensile strain modulates the transition temperature. Recently, it is found that large bandgap materials such as diamond, silicon and silicon carbide can be metallized and become superconductor under biaxial shear deformation [56,57]. Besides, an amazing robust in T c was predicted when adding uniaxial strain to solid molecular hydrogen [58]. Here, we first added 5% and 10% biaxial tensile strain to the doped system with x = 0.15 and calculated their electron-phonon properties as above. Figures 8(d) and (e) and table 3 show how biaxial tensile strain affects T c . Under 5% and 10% biaxial tensile strain, λ decreases to 0.78 and 0.65, respectively, while ω log increases to 328.0 and 489.2 K. The decrease in λ takes advantage of the increase in ω log ; thus, the final T c decreases to 14.4 and 13.8 K, just opposite from the hole doping. This phenomenon can be explained from two aspects: on one hand, from table 3, we can see that N(E f ) decreases more than 60% when biaxial strain is applied; one the other hand, the depression of electron phonon coupling in low-frequency zone affects significantly to total T c . Then we add 5% uniaxial tensile strain to the doped system with x = 0.05 and x = 0.15. We chose the full Brillouin path Γ → M 3 → K 3 → M 2 → K 1 → M 1 → Γ in phonon linewidth calculation for accuracy. Compared to the pure system, N(E f ) decreases about 30% while electron phonon coupling in low-frequency zone does not change too much, as shown in figures 8(f) and (g). Under 5% uniaxial strain, when x = 0.05 and 0.15, λ decreases about 5% to 1.07 and increases about 15% to 1.56, respectively. However, ω log decreases to 174.9 and 117.7 K which lowers T c to 13.4 and 13.9 K. Comparing biaxial and uniaxial strain, for example, at Table 3. Doping level x, tensile strain ε and calculated electron-phonon properties of C 8 H 2 . The row with (bi), (uni) and (an) represents the results with the biaxial, uniaxial strain and anharmonic effect.
x (holes/cell) ε  x = 0.15, although they all tend to lower N(E f ), depress the superconductivity and trigger a similar T c , the reasons are different. From the perspective of γ qj ω qj , biaxial strain weakens electron phonon coupling because of obvious phonon softening in low-frequency zone while uniaxial strain does little effect on phonon softening. The large anisotropy induced by uniaxial strain changes the symmetry of Brillouin zone and further changes the distribution of frequencies and γ qj ω q j in low-frequency zone. For example, in figures 8 (f) and (g), we can see that near K 3 and K 1 point, both of which should be the same in pure system, γ qj ω qj is different. The phonon frequencies at M 3 point also differ from that at M 1 or M 2 point. This anisotropy induced change of distribution significantly affects PHDOS and α 2 F(ω) in low-frequency zone compared to isotropic system, not to mention ω log , and further modulates T c .
It has been found that the transition temperature of hydrogen-rich material H 3 S and LaH 10 cannot be exactly predicted unless taking into account the anharmonic effect [59][60][61]. The presence of hydrogen atoms indicates that these systems may exhibit anharmonic effects. We used SSCHA to examine the anharmonicity in doped C 8 H 2 . To guarantee that the final system is still dynamically stable, we choose doping level x = 0.05 for safety. From figure 8(h), we can see that the anharmonicity softens two high-frequency branches, while the distribution of γ qj ω qj does not show an obvious difference from the harmonic situation. The reason is that the contribution of the coupling in the low-frequency zone contributes the majority of λ, similar to the harmonic situation. The insignificant difference in γ qj ω qj between harmonic and anharmonic situations indicates that the corrected coupling caused by anharmonicity is rather weak. Under the anharmonic potential, λ and ω log change from 1.12 to 1.03 and 210.9 to 235.4 K, and final T c from 17.3 to 17.2 K. The results indicate that hydrogenated graphene shows weak anharmonicity because the number of hydrogen atoms is rather small, and unlike the caged structures, the two hydrogen atoms on both sides of the single graphene layer are not strong enough to effectively modify the harmonic potential. It should be mentioned that the appearance of anharmonicity does not necessarily rely on the hydrogen atoms and the dimensions of the systems. In other systems without hydrogen atoms [62] and even monolayer structures [63], a strong anharmonic effect still appears.
To further reveal the electron-phonon coupling mechanisms between electron orbitals and phonon vibrations, we calculate the phonon vibration modes at the K point in the low-frequency zone where the size of magenta circles is the largest in figure 8(b) which means strong coupling occurs at this point. We choose the case with doping level x = 0.05 without loss of generality. As shown in figures 9(a)-(c), there are three vibration modes at the K point that contribute significantly to coupling. The vibration mainly comes from the out-of-plane mode of the carbon atoms. The projected band structure and DOS (figure 9(d)) indicate that the electrons around the Fermi level are mainly from the p z orbitals of C atoms. The contribution from the hydrogen atoms near the Fermi level is particularly small, which further indicates weak anharmonicity. The origin of superconductivity comes from two parts: hydrogenation and hole doping. Hydrogenation helps localize the p z orbitals of carbon atoms near the Fermi level. By doping moderate holes, the VBM moves across the Fermi level. The coupling between out-of-plane vibration modes and the p z orbitals of carbon atoms triggers superconductivity.

Conclusion
In this work, we study the electronic and superconducting properties of hydrogenated graphene with two different coverages: C 8 H 2 and C 50 H 2 . In C 8 H 2 , there are three different configurations, and they exhibit different properties. For instance, the case d = a is a zero-bandgap semimetal, while the configuration with d = 2a is a wide bandgap insulator, and the case d = √ 3a shows FM with a spin moment of 2 μ B /cell. For C 50 H 2 , eleven configurations are found together with interesting regularity. First, when two hydrogen atoms are attached to different sublattices, the final structures will be NM rather than AFM due to the preservation of π bond; however, when they are attached to the same sublattice, they turn into FM with a spin moment of 2 μ B /cell. We further determined that magnetism mainly originates from the hydrogen atoms and the nearest carbon atoms of the attached atoms. Second, as d increases near the maximum in the NM cases, the DOS near the Fermi level will be maximally localized, the Dirac points will disappear and the bandgap will be opened, triggering the topological transformation. Third, the positions of Dirac points are highly related to the space group, which is directly decided by the hydrogenation position. This can be explained by different hopping parameters in TB theory using VCA.
In addition, by doping moderate holes into the case d = 2a in C 8 H 2 , we can tune the system into a superconductor. The superconducting transition temperature T c will rise from approximately 17.3 K at doping level x = 0.05 to approximately 20.2 K at the critical doping level x c = 0.17 due to the stronger coupling from the low-frequency phonon. After adding biaxial and uniaxial tensile strain ε to the doped system, T c will decrease because of the modulation of phonons in low-frequency zone: biaxial tensile strain softens phonons and uniaxial tensile strain makes phonon dispersion unsymmetrical. Due to the presence of hydrogen atoms, we study the anharmonic effect at x = 0.05. The calculated T c is approximately 17.3 and 17.2 K in harmonic and anharmonic situations, respectively. The reason for the weak anharmonicity can be attributed to the small number of hydrogen atoms and the absence of the caged structure, together with a small contribution near the Fermi level. The superconductivity is obtained by hydrogenation and hole doping, which will strengthen the coupling between out-of-plane phonon modes and p z orbitals of carbon atoms. Therefore, by applying hole doping and tensile strain on hydrogenated graphene, we can obtain superconductivity and achieve effective modulation of the superconducting properties, which is robust against the anharmonic effect. This work combines the hydrogenation and superconductivity in graphene and shows the significance of further study in doping graphene into a superconductor.