Thermodynamic Properties of the Superconducting State in Metallic Hydrogen: Electronic Correlations, Non-conventional Electron-Phonon Couplings and the Anharmonic Effects

Thermodynamical properties of the superconducting state in metallic hydrogen were determined on the basis of the model of two compressed hydrogen planes. We took into account both the on-site and the inter-site electronic correlations (U and K), as well as the relevant non-conventional electron-phonon coupling functions (gU and gK). We proved, within the Eliashberg formalism, that the maximum value of the critical temperature of transition into the superconducting state is about 200 K for the harmonic approximation, and about 84 K for the Morse anharmonic approximation. Omission of the electronic correlations results in a considerable overstatement of the TC value. On the other hand, the TC value is remarkably understated if the non-conventional interactions are disregarded. Other thermodynamic quantities, such as the order parameter, the jump in the specific heat value, the thermodynamic critical field, and the upper critical field, take the values for which the non-dimensional ratios RΔ, RC, RH and RH2 do not differ substantially from the predictions of the BCS theory.


Introduction
The metal-insulator transition (MIT) in the molecular solid hydrogen at high pressure (p ∼ 25 GPa) was predicted by Wigner and Huntington [1]. Experimental studies did not confirm MIT in solid hydrogen up to 342 GPa [2]. However, Weir  and the temperature equal to T = 3000 K [3]. The newest theoretical predictions based on the exact-exchange calculations fixed the metallization pressure at 400 GPa [4].
In 2004, Ashcroft suggested that the hydrogen-rich compounds could be the high-T C superconductors at lower values of pressure than hydrogen due to the chemical precompression [20]. Ten years later Drozdov, Eremets, and Troyan announced the discovery of the superconductivity in the compressed H 2 S and H 3 S compounds at 150-203 K [21,22]. It should be noticed that the existence of the superconducting state in H 2 S and H 3 S compounds had been theoretically predicted by Li et al. [23] and Duan et al. [24,25] before the above-mentioned experiments were carried out.
The superconducting state in the system of compressed hydrogen and sulphur is induced by the electron-phonon interactions. Thermodynamic properties of the considered phase can be successfully described using the Eliashberg formalism and the density functional theory (DFT) [26][27][28][29][30][31].
In 2018, the superconducting state of the highest, so far, value of the critical temperature, T C = 215 K for p = 150 GPa or T C = 260 K for p ∈ (180-200) GPa, was experimentally confirmed in the LaH 10 compound [32,33]. The properties of the superconducting state in the LaH 10 were discussed on the theoretical level in [34]. The experimental results achieved by Drozdov et al. [32] were there related to the induction of the superconducting phase in the R3m structure (T C = 206-223 K), and the experimental results obtained by Somayazulu et al. [33] associated with the superconducting state induced in the F m3m structure, for which the critical temperature can reach the value of 280 K. There exist also strong theoretical indications suggesting that the T C value could be risen even higher by doping the LaH 10 compound with scandium (T C ∼ 294 K) or yttrium (T C ∼ 274 K) [35].
In the present paper we will discuss the properties of the superconducting state in metallic hydrogen in the framework of the model of two compressed hydrogen planes. We took into account both the relevant electronphonon coupling functions, as well as the on-site and the inter-site electronic correlations. We will show that one can expect induction of the superconducting state with the maximum critical temperature of the order of 80 K around 500 GPa.

Formalism
We assumed the effective linear electron-phonon coupling for the purpose of our considerations. For the half-filled electron band ( n = 1 and μ = 0), the statistic operator takes the form [36]: where c ( †) kσ denotes the annihilation (creation) operator for the electron with momentum k and spin σ ∈ {↑, ↓}. The symbol b ( †) q is the phonon annihilation (creation) operator, and φ q = b † −q + b q . Despite the fact that the operator (1) has the form of the conventional Fröhlich Hamiltonian [37], it is actually more general because we included both the on-site (U ) and the inter-site (K) Hubbard correlations in the effective electronic dispersion relation: ε k = ε + t cos (kR) + U 2 + K , where ε represents the electronic orbital energy and t is the electronic hopping integral. Additionally, the effective electron-phonon coupling function takes the form: g kq = g ε sin(qR) 2 . The quantities g ε and g t are conventional electron-phonon coupling functions, while g U and g K denote the non-conventional electron-phonon coupling constants corresponding to the U and K energies. The functions ε k and g kq were obtained after applying the mean-field approximation to the Hamiltonian modeling the system of two compressed hydrogen planes. The initial Hamiltonian included the four-fermion terms representing the pure electronic correlations. The symbol E q denotes the phonon dispersion relation: where ω x 0 is to be determined for either the harmonic (x = H ) or the anharmonic Morse (x = M) approximation. The equilibrium parameters of the operator (1) were calculated in the variational way, which is presented in detail in the next chapter.
In order to achieve the equations describing the thermodynamics of the system, the necessary calculations were done within the Nambu spinor representation [38,39]: under the assumption that the spinors meet strictly the anticommutation rule: represents one of the Pauli basis matrices τ j . The Hamiltonian (1) in spinor representation takes the form: The matrix Green function G k (iω n ) = k | † k iω n is now expressed by the following matrix: The diagonal elements of G k (iω n ) describe the properties of the normal state for temperature values exceeding the critical temperature of transition between the superconducting and the normal state. On the other hand, the antidiagonal elements of Green function allow to determine the properties of the superconducting phase. The symbol ω n represents the fermionic Matsubara frequency: ω n = π/β (2n − 1), where β is the inverse temperature (β = 1/k B T ). The electronic Green function satisfies the Dyson equation [40]: where the bare propagator can be written as follows: The self-energy takes the form: The procedure of achieving self-consistency was applied during the derivation of the Eliashberg equations. The following self-energy distribution in Pauli basis matrix was taken into account: where Z k (iω n ) denotes the wave function renormalisation factor. From the physical point of view, it determines renormalisation of the electron band mass by the effective electron-phonon interaction. The symbol ϕ k (iω n ) represents the real part of the order parameter function. Notice that the order parameter is expressed as k (iω n ) = ϕ k (iω n ) /Z k (iω n ).
Analytical calculations lead to the Eliashberg-type system of equations [41]: .
The form of the derived equations is identical to the form of the classic Eliashberg equations [41], with (8)-(9) containing effective functions (ε k , E q , and g kq ).
Nevertheless, to analyze the system of (8)- (9), one can use all methods normally used in classic Eliashberg formalism.
In the present work we will discuss the superconducting properties of the system at the isotropic level. Let us consider the Eliashberg function being one of input parameters: where α 2 F (k, ω) = ρ (0) q g 2 kq δ ω − E q and w k = δ (ε k ) /ρ (0). The symbol δ denotes the Dirac distribution, ρ (0) is the electronic density of states at the Fermi level.
The second input parameter is the so-called Coulomb pseudopotential (μ ) occurring in the function: μ (ω m ) = μ θ (ω C − |ω m |). It models explicitly the influence of the depairing electronic correlations on the thermodynamic properties of the superconducting state. The Coulomb pseudopotential can be determined with accuracy of the first order with respect to the Coulomb potential (U ) by the formula [42]: μ = μ/[1 + μ ln (W/ω D )], where μ = ρ (0) U . Symbols W and ω D occurring in the formula denote the half-width of the electronic band and the Debye frequency, respectively. It can be easily seen that μ differs from μ in taking into account the logarithmic reduction of the Coulomb interactions. Physically, it results from the fact that the Coulomb interactions propagate faster than the phonon-mediated interactions. The ratio W/ω D is therefore a good measure of the velocity of the interaction spreading within the considered system. Notice that θ is the Heaviside function and ω C determines the cut-off frequency: ω C = 5ω D [43].
The parameters of the superconducting state for the isotropic case can be calculated from analytical formulae if both the coupling constant and the Coulomb pseudopotential take low values (λ < 1.5 and μ < 0.2). The obtained results should be in agreement with the results obtained from the isotropic Eliashberg equations. In particular, the value of the critical temperature can be estimated from the Allen-Dynes formula [44]: The symbols used in the formula (11) are explained in Table 1. Additionally, we calculated the non-dimensional ratios, which are equivalents of the universal constants defined within the BCS theory [45,46]: The values of universal ratios according to the BCS theory are equal to 3.53, 1.43, 0.168, and 0.727, respectively. Symbols occurring in formulae (12)-(15) are of the following meaning: (0) is the value of the order parameter for the temperature of zero kelvin, C (T C ) represents the jump of the specific heat at the critical temperature, C N is the value of the specific heat for the normal state, H C Table 1 Parameters λ, ω ln , and ω 2 represent the electron-phonon coupling constant, the logarithmic phonon frequency, and the second moment of the normalized weight function, respectively. Functions f 1 and f 2 are the strong-coupling correction function and the shape correction function, respectively denotes the thermodynamic critical field, and H C2 is the upper critical field. If one knows the explicit form of the Eliashberg function, the R -R H 2 ratios should be calculated from [43,47,48]:

Estimation of Hamiltonian H Parameters
The ground state parameters of Hamiltonian H were calculated on the basis of the following formula: H = E p + E e + E f , where E p = 2/R represents the energy of the proton-proton repulsion (R is the distance between the protons: R = |R|), E e denotes the lowest energy of the electron eigenstate, and E F is the term related to the pressure: The energy of electronic states was estimated using the extended Hubbard Hamiltonian (H e ), which models all twosite interactions between the nearest atoms of two planes subjected to compression [49]: where c † jσ , ( c jσ ) is the electron creation (annihilation) operator in the position representation. The spin operator is given by S j , while the scalar product of spin operators , the operator of directing the spin up and down being introduced here, namely, The parameters of our interest which determine the values of the ε k , E q , and g kq functions in the Hamiltonian H are expressed by the following integrals: The symbol j (r) denotes the Wannier function: j (r) = A φ j (r) − Bφ l (r) , where φ j denotes 1s Slater-type orbital: φ j (r) = α 3 /π exp −α|r − R j | , and α is the variational parameter. The values of the A and B parameters were selected so that the Wannier function was normalized: The atomic overlap (S) was calculated from the formula: The physical meaning of the above listed integrals and some of their properties are given hereafter. The quantity ε is the electronic orbital energy. Its value tends to 1 Ry when the interplanar distance tends to infinity. The subsequent quantity t represents the hopping integral, which tends to zero as R → +∞. U denotes the on-site Hubbard integral. It should be noticed that it has no equivalent in electrodynamics, because two conventional charges cannot occupy the same spatial point simultaneously. The symbol K stands for the inter-site Hubbard integral. This quantity corresponds to the potential energy of the Coulomb interaction (2/R). The value of K is the closer to the potential energy value, the greater is the distance between hydrogen atoms belonging to two different planes. J is the Heisenberg exchange integral. The term J ĉ † 1↑ĉ † 1↓ĉ 2↓ĉ2↑ + h.c. models the hopping of electron singlet pairs. The full interaction of the Heisenberg-Dirac exchange is given by the formula: −J 2Ŝ 1Ŝ2 − 1 2n 1n2 . Considerations rarely take into account the so-called energy of the correlated hopping V . The process to which it is applicable consists in the hopping of an electron to the neighbouring atom, while the other one hops the same way and then returns to the previously occupied site.
The performed numerical calculations let us estimate the electronic parameters of the considered system with an accuracy of six decimal places. The results are gathered in Table 2. It can be seen that the exchange integral and the energy of the correlated hopping take so small values that they can be omitted in further considerations. The onsite and the inter-site Hubbard integrals assume high values, which grow even higher as the force F increases.
Please notice that we use the concept of the compressive force in this work, but, of course, its value can be related to the value of pressure acting on the considered system. The pressure was calculated on the basis of the estimated distance between hydrogen atoms. For this purpose, the first-principles calculations were performed according to the DFT method [50] as implemented in the Quantum-Espresso package [51,52]. The calculations were conducted at the gamma point using the Rappe-Rabe-Kaxiras-Joannopoulos ultrasoft pseudopotential, the energy cut-off value of 80 Ry, and the charge density cutoff of 320 Ry. The threshold value for the total energy change, equal to 10 −12 Ry, was established in order to ensure the energy convergence. The obtained results are gathered in Table 3, which juxtaposes the values of both the compressive force and the pressure. We added also values of both the equilibrium distance between hydrogen atoms and the variational parameter (equilibrium inverse size of the orbital), as well as values of the ground state energy.
In the numerical calculations we assumed the smallest value of F equal to 3 Ry/a 0 , which corresponds to the pressure of 494 GPa. This value is well above the value of metallization pressure for hydrogen, which is about 400 GPa based on the DFT calculations [4]. In addition, the pressure of 494 GPa correlates well with the experimentally achieved pressure of 495 GPa, for which it has been reported the hydrogen metallization [53,54]. However, it should be noted that many experienced researchers [55] questions the accuracy of the results contained in the paper [53]. Below the pressure of 494 GPa our model predicts the sharp drop in the critical temperature due to persistently high values of the on-site and inter-site Hubbard integral (see Table 2). The superconducting state disappears for the pressure of 483 GPa.
The vibrational energy levels were calculated on the basis of the enthalpy H . The potential energy in the harmonic approximation is given by The quantum harmonic oscillator energy has the well-known form: E p = ω h 0 (p + 1/2), wherein ω h 0 = k h /m . The symbol m represents the reduced mass: m = m p /2 = 918.076336 (m p is the proton mass). The vibrational energy levels can also be calculated using the Morse potential (V M ), which approximates the general form of the intermolecular energy curve. The Morse potential for F = 0 can be written as where E D is the dissociation energy measured from the minimum value of V M (R), and α M represents the measure of the curvature of the potential at its minimum. In the case of F > 0, the Morse potential should be generalized to the following form:  Table 4.
The parameters, g x = δx/δR dx/dR, where x ∈ {ε, t, U, K}, were calculated in order to determine the electron-phonon interaction. The quantization was performed as follows: δR = 1/2m E d † + d , where the symbol d † , (d) denotes the phonon creation (annihilation) operator. The equilibrium values of the g x function are presented in Table 5. The g J 0 and g V 0 coupling constants were neglected in considerations because they take very low values.
To end this chapter, let us recall that the numerical packages employed for calculations based on the firstprinciples approach were tested on the hydrogen molecule and its ions. In the case of H + 2 , the very advanced numerical analysis performed by Schaad and Hicks allowed to obtain R 0 = 1.9972 a 0 and E 0 = −1.205269238 Ry [56]. We obtained the results of R 0 = 2.003296 a 0 and E 0 = −1.173013 Ry. The equally accurate estimation of R 0 and E 0 was conducted for the H 2 molecule. The results  [57,58]. Unfortunately, inconsistent data are found in the professional literature for the H − 2 ion. The early theoretical paper by Eyring, Hirschfelder, and Taylor-who used the valence bond technique with two variation parameters-reported the stable ground state with the minimum at R 0 = 3.40151 a 0 [59]. This result correlates fairly well with our one, which is R 0 = 3.476828 a 0 (E 0 = −1.947958 Ry). However, the paper [60] suggests that the H − 2 ion is not stable and undergoes the autoionization into H 2 and an electron at infinity.

Characteristics of the Thermodynamic Properties of the Superconducting State
The performed ab initio calculations made possible to determine all significant thermodynamic parameters of the superconducting state induced in the considered system, such as the critical temperature, the order parameter, the jump in the specific heat, the thermodynamic critical field, and the upper critical field.
First, we determined values of the isotropic Eliashberg function defined by the formula (10). It was the basis for determination of the dependence of the electron-phonon coupling constant λ, the logarithmic vibration frequency ω ln and the second moment of the normalized weight function ω 2 on the value of the compressive force F . We also obtained a similar dependence for the Coulomb pseudopotential μ (in both the harmonic approximation (μ H ) and the anharmonic Morse approach (μ M )), directly by virtue of its definition. The results are presented in Fig. 1a-d. As long as the full model is considered (i.e., both U and K, and both g U and g K are taken into account), the values of the parameter λ are not especially high. The considered system is essentially within the region of indirect electron-phonon coupling. Moreover, the electron-phonon coupling constant distinctly decreases with an increase in the compressive force F , approaching the values accepted by the mean-field BCS theory (λ ∼ 0.3) [45,46]. Please notice that the obtained result is found on the account of Table 3 The values of the compressive force along with the respective values of pressure. Additionally, the values of the equilibrium distance, the equilibrium inverse size of the orbital, and the ground state energy are specified  considering all the significant electron correlations modeled with both the on-site and the inter-site Hubbard integral. Neglecting the depairing electronic correlations leads to distinct overestimation of the electron-phonon coupling constant, as can be seen in Fig. 1a.
The logarithmic phonon frequency ω ln was introduced by Allen and Dynes for better estimation of the critical temperature value characteristic for the metal-superconductor phase transition [44]. It replaced the Debye frequency ω D in the McMillan formula for the critical temperature [61]. The dependence of ω ln on the force F is plotted in Fig. 1b. One should notice that the values of the logarithmic frequency (as well as of the Debye frequency) are very high. This effect results from the very small mass of the hydrogen atomic core and was first observed by Ashcroft [5]. Please notice that the relationship ω ln ∼ 1/ √ m p is valid as a first approximation. The data in Fig. 1b indicate also that the logarithmic frequency increases with an increase in system compression, what results from the stiffening of the system (both k H and k + M elastic constants increase). Omission of the electronic correlations during the system analysis leads to the significant underestimation of the ω ln valuethe misestimation opposite to that in the case of the value of λ.
The frequency ω 2 does not affect the value of critical temperature as significantly as ω ln (formula (11)). Generally, √ ω 2 is applied for determination of values of the correction functions f 1 and f 2 . Its values are shown in Fig. 1c. One can see that it grows as the force F increases. Figure 1d presents the dependence of the Coulomb pseudopotential μ on the compressive force. The Coulomb pseudopotential is a significant quantity in the theory of the phonon-induced superconducting state, because it models directly the influence of the depairing electronic correlations on the value of the critical temperature. One can easily see that the change in the compressive force affects μ only slightly (both in the harmonic and the anharmonic Morse approximation). The average value of the Coulomb pseudopotential is equal to about 0.17. This means that the depairing electron correlations is not very large, what results from the screening of the Coulomb-type interactions in the considered system. The result we came to, however, Table 5 The values of the electron-ion coupling constants at the equilibrium point for selected values of F F (Ry/a 0 ) g ε 0 (Ry/a 0 ) g t 0 (Ry/a 0 ) g U 0 (Ry/a 0 ) g K 0 (Ry/a 0 ) g J 0 (Ry/a 0 )  [62,63]. If one would like to confirm the result in such a case, he could try to calculate the value of μ starting from the formula which would include the contributions of the U 2 order [64]. However, very high values of μ may as well conceal the induction of the phase being competitive with the superconducting phase [65]. Figure 2 illustrates the dependence of the critical temperature on the applied force F . The critical temperature reaches the maximum value of about either 200 K (for the harmonic approximation and F ranging from 3 Ry/a 0 to 20 Ry/a 0 ) or 84 K (for the anharmonic approach and Fig. 2 The influence of the force F on the critical temperature value the values of F close to 3 Ry/a 0 ), but the results obtained within the anharmonic approach are more exact. For the purpose of comparison, we referred to the results reported in the literature for F ∼ 3 Ry/a 0 (i.e., within the pressure range from ∼ 400 GPa to ∼ 500 GPa). The following estimated values of critical temperature were obtained for the DFT approach within the harmonic approximation: 84 K for 414 GPa [10][11][12][13], 179 K for 428 GPa [17,19], 242 K for 450 GPa [10][11][12], and 360 K for 539 GPa [6,14]. These results correspond well with our T C value calculated within the harmonic approach, which is equal to about 200 K. However, one should remember that the anharmonic effects significantly reduce the value of T C .
It is worth noticing that the value of T C drops down while compression increases. A similar effect we found also in the H 3 S hydrogenated compound [29]. Interestingly, our results with respect to the value of T C for this compound seem to harmonize relatively well with the maximum values of the critical temperature experimentally measured for H 3 S or also LaH 10 , reported to be [T C ] H 3 S = 203 K [22] and [T C ] LaH 10 = 215 K [32] (or else [T C ] LaH 10 = 260 K [33]). The values of the critical temperature for the hydrogenated compounds, higher than those for the metallic hydrogen, result from the presence of heavier elements, which make a noticeable contribution to the Eliashberg function within the range of low frequencies [26,35].
Please take notice that either neglecting one of the channels of the depairing electronic correlations or the inaccurate determination of the value of the Coulomb pseudopotential lead to a significant overestimation of the value of T C . This effect can be often observed in analyses based on the DFT method, when the harmonic approximation is used and the μ value is underestimated. On the other hand, elimination of the non-conventional electron-phonon interaction channels (g U and g K ) leads to the underestimation of the critical temperature.
The results obtained for other thermodynamic quantities characterizing the superconducting state are presented in Fig. 3a-d. The subfigures show the values of the subsequent non-dimensional thermodynamic ratios defined by formulas (12)- (15). Deviations from the BCS theory can be observed only for the low values of the compressive force acting on the system. This range of low values of F corresponds to the region of the indirect electron-phonon coupling, as expected. So, despite the high value of critical temperature found for the metallic hydrogen, the values of its thermodynamic parameters generally good fit with the predictions of the BCS theory. Note that the omission of U and K electron interactions leads to incorrect values of the ratios R -R H 2 , significantly different from the predictions of BCS theory. From the physical point of view, the failure to take into account the electron correlation causes the significant overestimation of electron-phonon coupling constant λ and underestimation of logarithmic phonon frequency ω ln (see Fig. 1). Also underestimated or even neglected (U = K = 0) are the electron deparing effects (μ ). As a result, we obtain too high values of the critical temperature in the area of the relatively strong electron-phonon coupling, which results in strong deviations from the results of the BCS theory.

Conclusion and Discussion of the Results
The work presents an analysis of properties of the superconducting state which can be potentially induced in metallic hydrogen. The ab initio calculations were carried out for the model of two compressed hydrogen planes in order to overcome the great mathematical complexity of the problem. We took into account all physically important channels of electronic correlations (U and K). Additionally, there were considered the non-conventional electronphonon coupling functions, g U and g K , related directly to the on-site and the inter-site Hubbard interactions, respectively.
The thermodynamic properties of the superconducting state were determined within the Eliashberg formalism. It is the most advanced approach to the problem, which makes possible the quantitative description of the superconducting phase [41,43]. We determined the form of the isotropic Eliashberg function, starting from the ab initio calculations, then we found the input parameters to the Allen-Dynes formula [44]. We proved that the superconducting state can be induced in the considered system, its maximum value of critical temperature being equal to about 200 K (in the harmonic approximation) or else 84 K (in the more exact anharmonic approach). The high T C value results physically Fig. 3 The values of the nondimensional thermodynamic ratios: a R , b R C , c R H and d R H 2 versus the force F from the very high value of the logarithmic phonon frequency, which is in turn related to the very small atomic core mass of hydrogen (a single proton). The electronphonon coupling constant λ takes rather small values within the range of low values of the force F . In the case of high compression, its value falls within the range of weak electron-phonon coupling. The obtained results indicate that the effective value of the depairing electronic correlations in the metallic hydrogen is relatively low (μ H ∼ μ M ∼ 0.17) and is only weakly dependent on the compressive force F . Please notice that the calculated values of the critical temperature correlate well both with the theoretically found values of T C reported in other publications and with the maximum values of critical temperature experimentally measured for the H 3 S and LaH 10 compounds [22,32,33].
Despite the high value of the critical temperature in the metallic hydrogen, the thermodynamic parameters of the superconducting state generally take values for which the non-dimensional thermodynamic ratios R , R C , R H , and R H 2 only slightly deviate from the values predicted by the mean-field BCS theory. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.