GW space-time method: Energy band-gap of solid hydrogen
GGW space-time method: Energy band-gap of solid hydrogen
Sam Azadi, Arkadiy Davydov, and Evgeny Kozik Department of Physics, King’s College London, Strand, London WC2R 2LS, UK ∗ Department of Physics, King’s College London, Strand, London WC2R 2LS,UK; Physik-Institut, University of Zurich, 8057 Zurich, Switzerland † (Dated: June 12, 2020)We implement the GW space-time method, in which the Greens function G and the screenedCoulomb interaction W are efficiently represented in the real-space and imaginary-time domain atfinite temperatures, paying particular attention to controlling systematic errors of the representa-tion. Having validated the technique by the canonical application to silicon and germanium, weapply it to calculation of fundamental band gaps in hexagonal solid hydrogen with the bare Greensfunction obtained from the local density approximation (LDA) and the interaction screened withinthe random phase approximation (RPA). The results, obtained from the asymptotic decay of the fullGreens function without resorting to analytic continuation, predict that the solid hydrogen above200 GPa can not adopt the hexagonal-closed-pack (hcp) structure. The ability of the method tostore the full G and W functions in memory with sufficient accuracy is crucial for its subsequentextensions to include higher orders of the diagrammatic series by means of diagrammatic MonteCarlo algorithms. Keywords: High-pressure solid hydrogen, GW method, energy band gap
I. INTRODUCTION
The metallization of low-temperature high-pressuresolid hydrogen is a challenging and controversial prob-lem in the field of condensed matter physics, materialsscience, and high-pressure physics. It has been the sub-ject of wide range of studies since 1935 by experiment,theory and more recently computational techniques. Themain interesting exotic features of metallic hydrogenare related to the room-temperature superconductivity ,which was predicted based on the phonon mediated the-ory of superconductivity, a possibility of new phase ofsuperfluid liquid metallic ground state , and the im-portance of metallic solid hydrogen in astrophysics .Determining the metallization mechanism and metalliza-tion pressure of solid hydrogen is still an open problem.Experiments including Infrared (IR) and Ramanmeasurements indicate the existence of several low-temperature phases. Phase I, which is stable up to110 ± . Changes in the low-frequency regionsof the Raman and IR spectra imply the existence of phaseII, also known as the broken-symmetry phase, above110 ± ∼
150 GPa isverified by a large discontinuity in the Raman spectrumand a strong rise in the spectral weight of the molecu-lar vibrons . Phase IV, characterized by two vibronsin its Raman spectrum, was discovered at 300 K andpressures above 230 GPa . A new phase has beenobserved at pressures above 200 GPa and temperatureof 480 K . The experimental results indicate that H and hydrogen deuteride at 300 K and pressures greaterthan 325 GPa transform to a new phase, characterizedby substantial weakening of the vibrational Raman ac-tivity. Recent single-crystal X-ray diffraction study of solid hydrogen at pressures of up to 254 GPa divulges thecrystallographic nature of the transition from phase I tophases III, and indicates that hydrogen molecules remainin the hexagonal close-packed (hcp) crystal structure .Going beyond the hydrostatic pressure of 400 GPa,at which hydrogen is predicted to become metal ,is a very challenging task. Experimentally, it is notclear whether or not the low-temperature phase III ismetallic, although the existence of semimetallic molec-ular hydrogen at pressure above 350 GPa has beenobserved . Metallization can take place either via astructural transformation or via band gap closure. Thepossible structural transformations, which is expectedto cause the metallization, are insulator-molecular tometallic-molecular and insulator-molecular to metallic-atomic phase transitions . The recent Infra-red (IR)measurements suggest that metallization of hydrogenproceeds within the molecular solid .For ab-initio electronic structure calculations of high-pressure solid hydrogen the geometry of the system,which is experimentally not well identified for pressuresabove 120 GPa, substantially affects the results. Fromthe computational point of view, many properties ofsolid hydrogen including free energy phase diagram, crys-tal structure, lattice dynamics, and electronic band gapwere studied using density functional theory (DFT) and quantum Monte Carlo (QMC) methods .The DFT random structure searching method usingPerdew-Burke-Ernzerhof (PBE) exchange-correlation(XC) functional predicts several structures for thepressures above 100 GPa, which all adopt a form of dis-torted hexagonal symmetry. Almost all the predictedproperties by DFT based methods for high-pressure hy-drogen strongly rely on the XC approximation. Thesingle-particle DFT based methods are not accurate forelectronic band gap calculation , and going beyond DFTfor obtaining the fundamental and optical band gaps is a r X i v : . [ c ond - m a t . m t r l - s c i ] J un inevitable.Several structures with different symmetries, such as C /c , P bcn , P c , P /
22, and
Cmca , and different num-ber of atoms per primitive cell, were predicted as poten-tial candidates for the phase III . The energy differ-ence between the predicted structures, within the stud-ied pressure range, is ∼
30 meV/atom, which is be-low the chemical accuracy and beyond the accuracy ofwidely used XC-DFT approximations. All of the pre-dicted structures are obtained using only PBE-XC ap-proximation and there is no evidence if structure searchtechnique with other XC approximation produces thesame results. Also, it has been demonstrated that thePBE-XC functional is poor at describing the relative en-ergies of configurations and the properties of the molec-ular bond . Computational x-ray diffraction patternsof molecular structures at room-temperature are alsovery different from the crystal structures which are pre-dicted by structure searching based on the static enthalpyminimization . There is a lack of experimental sup-port for the predicted structures. Hence, In this workwe focuse on the hexagonal structure which is reportedby experiment .The diffusion quantum Monte Carlo (DMC) and vari-ational quantum Monte Carlo (VMC) methods can beused to estimate the excited states and energy band gap.By their definition, DMC and VMC techniques are pri-marily zero temperature ground state methods, and alsoare affected by single-particle and many-body finite-size(FS) errors . Especially, DMC calculations of excita-tions in solids is challenging because of a 1 /N effect: thefractional change in the total ground state energy dueto the presence of a one- or two-particle excitation is in-versely proportional to the number of electrons N in thesimulation cell. Additionally, because of large numberof atoms per primitive cell in almost all of the predictedstructures for solid hydrogen, controlling the FS errorsin DMC calculations for high-pressure solid hydrogen isa challenging task.A more efficient approach to excited state properties isprovided by the many-body perturbation theory withinthe GW approximation. In principle, the QP energies canbe obtained using the full self-consistent GW method inwhich G is determined self-consistently from Σ = GW,which in turn produces G . It was shown that thevalence bandwidth of the homogeneous electron gas ob-tained by full self-consistent GW method is ∼
15% widerthan the non-interacting case , and also it overestimatesthe valence bandwidth of Ge and Si crystals . The dif-ficulties of the full self-consistent GW method in gener-ating QP levels, specially in challenging materials suchas oxides MnO and NiO, can be overcome by an ansatzof quasiparticle self-consistent GW approximation (QSGW). The first principle QSGW method is a pertur-bation theory, which uses a single particle Hamiltonianwith exchange-correlation potential derived from the self-energy at each iteration step of the self-consistent loop,which makes it independent of DFT exchange-correlation potential. It was shown that the QSGW approach accu-rately predicts excited-state properties of large numberof weakly and moderately correlated systems . To in-vestigate the accuracy and efficiency of pseudopotential(PP) approximation in GW calculations for sp semicon-ductors, the results of all-electron GW method with thefull-potential linear muffin-tin orbitals and full-potentiallinear augmented plane-wave were compared with GWPP approach .Applying the finite temperature Green’s function for-malism to real systems within the ”chemical accuracy”is a formidable task, because the eigenvalue spectrumof realistic Hamiltonians is wide. In this work we im-plement the space-time GW method in which themany-body quantities including the Green’s function Gand screened Coulomb interaction W, are defined in realspace and imaginary time. In our implementation, de-pending on the computational step, we choose the mostefficient representation, either imaginary time or imagi-nary energy on the one hand and reciprocal space or realspace on the other hand. The fast Fourier transforms(FFTs) are frequently applied for switching between rep-resentations. We work on imaginary axis which helpsus to operate with smooth decaying quantities with fastconvergence. Functions on the imaginary time axis arerepresented using two different techniques. In the firstone we use sampling on a logarithmic grid, which is closeto the so-called power mesh concept , for which the FTto the imaginary frequency domain is performed by firstinterpolating a function to homogeneous grid and thenby an FFT call. The second method utilizes the Cheby-shev polynomial representation, with a framework for FTbuild in .Traditionally, the quasi-particle excitation energies areextracted by performing analytic continuation of the self-energy to the real frequency domain, which is an ill-defined numerical problem. However, the representationin terms of Chebyshev polynomials allows us to capturethe Greens function asymptotic behavior at long imag-inary times in a controlled way and use it to robustlyextract quasi-particle energies without resorting to ana-lytic continuation. To the best of our knowledge, this isthe first time that the sampling of imaginary time G inperiodic systems is completed via Chebyshev-grid.We perform ’ one-shot-GW ’, meaning that G and W areconstructed from the local-density approximation (LDA)potential, that produce the self energy Σ = G LDA W LDA .The quasiparticle energies are also approximated as aperturbation correction to the LDA from the matrix ele-ments of the diagonal parts of Σ − V LDAxc , as explained inthe next section. Our GW calculations use the pseudopo-tential (PP) approximation. The main goal of the currentstudy is to determine the direct and indirect band gapsof hexagonal solid hydrogen which has no core electron,and therefore inaccuracies and errors due to use of PPapproximation are negligible. Also, solid hydrogen is thesimplest crystal without complexities of transition metaloxides. Therefore, the ’ one-shot-GW ’ method could bethe right tool for this work. In our method the G andW functions are represented and stored in memory infull, which makes it amenable to immediate extension tohigher orders of the diagrammatic expansion by means ofstochastically sampling the series by the well-establishedDiagMC technique . Hence, our work sets the foun-dation for future methodological developments address-ing the accuracy of the GW approximation itself, forwhich the algorithmic framework (DiagMC) is alreadyestablished.We first calculate the electronic band structure and en-ergy band gaps of two well-studied systems of Si and Gecrystals, and we compare the outcome of our GW sim-ulation with experiment and previously published data.Within the GW approximation, we investigate the nu-merical accuracy of the method by considering differentimplementations of representing the objects and arriveat an optimal representation where numerical errors arewell under control by means of appropriate extrapola-tions, which we carefully perform. We then study theenergy band gap of hexagonal solid hydrogen as func-tion of pressure. We calculate the direct and indirectband gaps of hcp and P /m phases, and compare ourresults with available experiments and other many-bodytechniques. We speculate the fundamental role of latticedynamics and vibrons in reducing the band gap, and weconclude that the solid hydrogen above 200 GPa can notadopt the hcp crystal structure.The paper is organized as follows; the section II consistof two parts. The first part II A describes the conceptand detailed aspects of our implementation of the GWspace-time method. We discuss the analytic continuationfor the self-energy, mesh-discretization, determination ofthe quasiparticle (QP) energies, and our approach forenergy band gap calculation. In the second part II B, wediscuss numerical details of DFT and GW calculations.In the first part of section III we test and validate ourapproach by obtaining results for the standard cases ofSilicon and Germanium, for which there are wealth ofboth theoretical and experimental results. In the secondpart of section III, we present and explain our results forthe hexagonal solid hydrogen, and we compare our datawith experiment. The section IV concludes the paper. II. COMPUTATIONAL DETAILSA. The GW space-time method
In the GW space-time (GWST) method, the Green’sfunction (G), polarizability (P), dielectric function ( ε ),dynamically screened Coulomb interaction (W), and self-energy (Σ) are defined and stored in real space and imag-inary time, typically on appropriate grids. The non-interacting Green’s function G DF T ( rr (cid:48) τ ) is constructedas follows: G DF T ( rr (cid:48) τ ) = (cid:88) k G k ( τ ) ψ k ( r ) ψ ∗ k ( r (cid:48) ) (1) with the Green’s function in the Kohn-Sham (KS) basis G k ( τ ): G k ( τ ) = θ ( − ξ k ) f β ( ξ k ) e ξ k τ + θ ( ξ k )(1 − f β ( ξ k )) e − ξ k ( β − τ ) (2)where the text k stands for combined band/k-point in-dex k = { l, k } , and f is the Fermi distribution func-tion. Plane-wave (PW) density functional theory (DFT)is employed for obtaining the KS band energies ξ k , andwavefunctions ψ k ( r ).Next, the random phase approximation (RPA) irre-ducible polarizability is calculated in real space and Mat-subara time: P ( rr (cid:48) τ ) = G DF T ( rr (cid:48) τ ) G ∗ DF T ( rr (cid:48) , − τ ) (3)and then Fourier transformed to the reciprocal space P ( rr (cid:48) τ ) → P ( qGG (cid:48) τ ), where q and G are Brillouinzone quasimomentum and reciprocal lattice vertices cor-respondingly. The P ( qGG (cid:48) τ ) is further transformed tothe imaginary frequency domain either through the se-ries of interpolations when using logarithmic grids, or viatransforming to an intermediate Chebyshev representa-tion, which acts as a proxy between two domains (seeSec. II B).The frequency-dependent polarisability is used to con-struct the inverse dielectric function: ε − ( qGG (cid:48) , iω ) = (cid:20) δ GG (cid:48) − P ( qGG (cid:48) , iω ) | q + G || q + G (cid:48) | (cid:21) − (4)The inverse FT ε − ( qGG (cid:48) , iω ) → ε − ( qGG (cid:48) , τ ) is donein the similar way as the direct one for polarisability. Thescreened Coulomb interaction is obtained by: W ( qGG (cid:48) τ ) = (cid:88) qGG (cid:48) π Ω ε − ( qGG (cid:48) , τ ) | q + G || q + G (cid:48) | (5)where Ω is the volume of unit cell, and thenFourier transformed to the direct space W ( qGG (cid:48) τ ) → W ( rr (cid:48) τ ). From that the self-energy operator Σ( rr (cid:48) τ ) = − G ( rr (cid:48) τ ) W ( rr (cid:48) τ ) and its expectation values (matrix el-ements) Σ lm k ( τ ) = (cid:104) ψ l k | Σ( rr (cid:48) τ ) | ψ m k (cid:105) are computed.Finally the last Fourier transform for self-energy is per-formed: Σ lm k ( τ ) → Σ lm k ( iω ).We use two methods to extract the quasi-particleexciation energies. The first, traditional one relieson analytic continuation of the self-energy diagonalΣ ll k ( iω ) → Σ ll k ( ω ), which is done by fitting to a mul-tipole model . The quasi-particle (QP) energies arefound as solutions of self-consistent quasi-particle equa-tion: ξ qpl k = ξ l k + Re[Σ ll k ( ω = ξ qpl k ) − v xcl k ] (6)where v xcl k stands for exchange-correlation potential of theKS scheme.Analytic continuation is not always reliable since its ac-curacy is generally uncontrollable. However, the problem −6 −5 −4 −3 −2 −1
0 200 400 600 800 1000 1200 l og G ( τ ) τ (a.u −1 )ncheb=100 GG fit −6 −5 −4 −3 −2 −1
0 200 400 600 800 1000 1200 l og G ( τ ) τ (a.u −1 )ncheb=150 GG fit −6 −5 −4 −3 −2 −1
0 200 400 600 800 1000 1200 l og G ( τ ) τ (a.u −1 )ncheb=200 GG fit −6 −5 −4 −3 −2 −1
0 200 400 600 800 1000 1200 l og G ( τ ) τ (a.u −1 )ncheb=250 GG fit FIG. 1. (Color online) Actual and fitted G k of eq. 7 at allirreducible k-points as a function of τ for Si. Different plotscorrespond to different number of Chebyshev polynomials. of determining the band gap does not require additionalapproximations, as the energy excitations nearest to theFermi energy can be found from the asymptotic decay ofthe full Greens function at long imaginary times: G k ( τ ) → a e k e ξ e k τ + a h k e − ξ h k ( β − τ ) , (7) τ (cid:29) [ ξ e ] − , [ β − τ ] (cid:29) [ ξ h ] − where a e ( h ) k are some constants and ξ e ( h ) are the energiesof electron(hole)-like excitations. The quantity G k cor-responds to the following sum over band indices of theGreen’s function matrix: G k ( τ ) = (cid:88) lm G lm k ( τ ) (8)Finally, the matrix G lm k is obtained first by solving thematrix Dyson equation:ˆ G k ( iω ) = (cid:104) iω ˆ1 − ( ˆΣ k ( iω ) − ˆ v xc k ) (cid:105) − . (9)with the subsequent Fourier transform ˆ G k ( iω ) → ˆ G k ( τ ).’Hats’ specify objects, which are matrices in the bandindex.The nonlinear dependence on ξ parameters in eq. 7 canbe extracted as a linear one out of the slope of log[ G ( τ )].Although the GWST scheme has an advantage that thefull non-diagonal (in band index) self-energy is easily ac-cessible, in this work however, only the diagonal partis assumed to be non-zero, in what is called the band-diagonal approximation. This is known to work well forSi and Ge. An example of the fitting procedure above isshown in Fig. 1.The advantage of this technique is that the QP ener-gies are obtained directly without the uncontrolled an-alytic continuation procedure and the subsequent needfor solving the quasiparticle equation 6. We compare theenergy band gaps of bulk Si and Ge which are obtained using both techniques. In the rest of this work, the re-sults from analytical continuation of Σ and equation (6)are named G W AC , and the QP-energies that are foundfrom eq. 7 are labelled as G W GF . B. Numerical details
The single-particle Kohn-Sham orbitals for Si crystalwith lattice-parameters of a = 10 . a.u are obtainedusing the local-density-approximation (LDA) XC DFTfunctional with energy cutoff 60 Ry, and 200 unoccu-pied energy bands. For Ge, with lattice parameter of a = 10 . a.u , we used LDA Norm-Conserving pseudopo-tential with four valence electrons, which was generatedwith a Scalar-Relativistic Calculation .For the solid hydrogen, we considered two molecularstructures with hexagonal crystal structure and symme-try of P /m −
16 and hcp ( P /mmc − P /m −
16 with c/a <
1, inwhich the c and a are the lattice parameters in the z and xy − plane directions, is the best candidate for the phaseI . The recent X-ray experiment indicates the stabilityof the hcp phase with c/a > . OurDFT calculations were carried out using the latest ver-sion of the Quantum-Espresso (QE) suite of programs and the Becke-Lee-Yang-Parr (BLYP) XC-DFT func-tional. We used a basis set of plane waves with an energycutoff of 100 (Ry). For the geometry and cell optimisa-tions a 16 × × k -point mesh is employed. For theGW calculations, we used a real space grid with numberof grid points N r = 16 , which gives an average absoluteerror in wave function norm around 10 − , and recipro-cal space mesh with number of grid points of N k = 4 .In all the calculations the QE-tabulated standard norm-conserving pseudopotential is used.The main quantities in GWST are nonlocal operatorsin spatial arguments, decay by increasing the distancebetween r and r (cid:48) . Apart of W( r , r (cid:48) ), the nonlocality isshort range and it decays faster than 1 / | r − r (cid:48) | . Followingthe original nomenclature , we define the first argu-ment of a two-point function within the unit cell (UC),while the second one within the so-called interaction cell(IC). The size of interaction cell is defined by the k-pointgrid as its inverse Fourier transform domain. The crystalsymmetries are used to reduce the computational costsand data storage requirements. Logarithmic Matsubara time mesh . We use the full[0 , β ] interval for sampling imaginary time, while imag-inary frequency range is in principle unbound and thecorresponding grid needs to be truncated, with an appro-priate account of high-frequency asymptotic. The log-arithmic sampling used for time axis is denser aroundends of the time interval, aimed at resolving effects ofhigh-energy empty states. Therefore these regions aresampled well at any other numerical parameters. How-ever, any FT needs a set of interpolations to linear grids.This can introduce irregular numerical errors, especiallyin the poorly sampled Matsubara time region away fromthe [0 : β ] boundaries. This problem can be avoided byusing the so-called ”uniform power mesh” , which isa mixture of both uniform and logarithmic grids. We,however appeal to an alternative and more systematicscheme, the Chebushev polynomial representation. Chebyshev polynomial representation provides an eas-ily controllable way of storing and manipulating corre-lation functions in both time and frequency domains, aswell as the Fourier transform algebra for switching be-tween them. The concept is essentially that of represent-ing a continuous function as an expansion in an orthonor-malized basis of polynomials, with the error controlled bythe number of basis functions N ch used in the represen-tation: F ( τ ) = N ch − (cid:88) l =0 F l T l ( τ ) , (10) F ( iω ) = N ch − (cid:88) l =0 F l T l ( iω ) , (11) T l ( iω ) = (cid:90) β dτ T l ( τ ) e iωτ (12)Once the coefficients F l are obtained and the polynomialbasis functions T l ( τ ) or T l ( iω ) are tabulated, computingthe Fourier transform is trivial.This framework has been discussed and developed forGreens functions in refs. . It was shown that, for abasis truncated at the first N ch polynomials, the optimaltime mesh for sampling F ( τ ) in the calculation of F l ,must include exactly N ch points:¯ τ αn = τ (cid:20) cos (cid:18) π n + 12 N (cid:19)(cid:21) , n = 1 ..N ch − T N ch ( τ )with an advantage that T l ( τ ) functions for l < N ch arelinearly independent, thus serving as a basis in the Mat-subara time domain. The τ on the right-hand side is amap from [ − β ], i.e., τ [ x ] = β ( x + 1) /
2. In asimilar fashion, imaginary frequency sampling is chosenas roots of the Fourier transform counterpart T N ch ( iω ),with minor difference for Bosonic and Fermionic cases .The disadvantage of this framework is in the fact thatthe more Chebyshev points/polynomials are needed forresolving the effects of states far away from the Fermilevel, which typically results in steep decay around edgesof the imaginary time interval. This behaviour can makeproblematic the treatment of core states or extremelyhigh-energy empty states, if included. Having 200 bandsfor Si, we have to use 250 polynomials and only about 100points in the logarithmic mesh to achieve negligible dis-cretisation errors. However, the advantage of Chebyshevpolynomial framework is that once the result is convergedwith N ch , the error of the FT is automatically negligi-ble. Second, all time/frequency dependent functions are W E g a p ( e V ) r k−2x2x2k−4x4x4k−6x6x6k−8x8x8 k i FIG. 2. Convergence of the Γ v → X c band gap with thenumber reciprocal and real space grid points with its linear fit.The inset shows the convergence with the number of k-pointsat N r → ∞ . well represented in all regions of time/frequency axis uni-formly, which is essential for obtaining energies from theasymptotic behavior eq. 7. Finally, the Fourier transform G ( iω ) → G ( τ ) generally requires a careful treatment ofthe high-frequency tail due to the slow convergence ofthe truncation error, but in the Chebyshev polynomialrepresentation the correct asymptotic is build in auto-matically. III. RESULTS AND DISCUSSIONA. Silicon and Germanium
For accurately controlling the systematic errors, wehave investigated the convergence of calculated bandgaps as function of size of interaction cell ( N k ) andnumber of points in sampling of the unit cell N r usinglogarithmic-grid (Table I) and Chebyshev-polynomials(Table II) representation of the imaginary-time depen-dence. Table I lists the Γ v → X c energy band gapsof Si and Ge which are calculated using G W AC and G W GF approaches. We used four reciprocal spacemeshes (k-mesh) with number of mesh points of N k =2 , , , , and three real space mesh with number ofgrid points of N r = 8 , , . This convergence pat-tern is visualised for Si in Fig. 2.The Γ v → X c energy band gap of Si and Ge whichis obtained by logarithmic-grid and G W GF is slightlylarger than G W AC with the same grid. The differ-ence between ”AC” and ”GF” energy band gaps whichare obtained by Chebyshev polynomials is smaller thenthe same difference that is calculated by logarithmic-grid. For instance, comparing the largest simulationof N k = 8 k-mesh and N r = 12 , the difference be-tween ”AC” and ”GF” band gaps is 67 meV and 45meV by logarithmic-grid and Chebyshev polynomials, re- Si- G W AC k-mesh N r = 8 N r = 10 N r = 12 N r → ∞ N k = 2 N k = 4 N k = 6 N k = 8 N k → ∞ ... ... ... 1.410(5)Si- G W GF N k = 2 N k = 4 N k = 6 N k = 8 N k → ∞ ... ... ... 1.51(1)Ge- G W AC k-mesh N r = 8 N r = 10 N r = 12 N r → ∞ N k = 2 N k = 4 N k = 6 N k = 8 N k → ∞ ... ... ... 1.272(6)Ge- G W GF N k = 2 N k = 4 N k = 6 N k = 8 N k → ∞ ... ... ... 1.36(1)TABLE I. Γ v → X c energy band gap of Si and Ge obtainedusing G W AC and G W GF at k-mesh of N k = 2 , , , and real-space-mesh with number of grid points ( N r ) of 8,10, and 12 at each direction. The logarithmic grid with alinear fitting to infinite real and reciprocal system size, whereasymptotic standard errors are included, is used.Si- G W AC k-mesh N r = 8 N r = 10 N r = 12 N r → ∞ N k = 4 N k = 6 N k = 8 N k → ∞ ... ... ... 1.346(2)Si- G W GF N k = 4 N k = 6 N k = 8 N k → ∞ ... ... ... 1.402(1)TABLE II. Similar calculations to table I for Si, but Cheby-shev polynomials is used. The number of Chebyshev polyno-mials is 250. spectively. The band gaps which are obtained by bothschemes are in good agreement with previously reporteddata .Figure 3 illustrates the electronic band structures of Siand Ge which are calculated by LDA and G W . The G W AC fundamental band gap of Si and Ge is 1.39, and0.86 eV, respectively. The fundamental band gap of Siand Ge obtained by G W GF is 1.46, and 0.94 eV, respec- −10−5 0 5 10 L Γ X W E n e r g y ( e V ) Si −10−5 0 5 10 L Γ X W E n e r g y ( e V ) Ge FIG. 3. (Color online) Band structure of Si and Ge at roomtemperature obtained by LDA (dashed line) and G W (solidline). The 8 × × k-mesh and N r = 12 real-space grid isused. The LDA fundamental gap for Si and Ge is 0.54, and0.24 eV, respectively. The G W AC and G W GF fundamentalgap for Si is 1.39, and 1.46 eV, respectively. The G W AC and G W GF fundamental gap (direct gap at Γ) for Ge is 0.86, and,0.94 eV respectively. tively. The experimental energy band gap of Si and Geis 1.17 and 0.78 , respectively. The calculated Si energyband gap is in excellent agreement with the one reportedin original representation of GWST method . B. Solid Hydrogen
All the following G W QP-energies results are ob-tained using our G W GF technique. Two hexagonalstructures of hcp and P /m with c/a > c/a <
1, respectively, are studied. Recent X-ray diffractionexperiment , observed three Bragg peaks which are con-sistent with an P /mmc structure with c/a ratio closeto c/a = (cid:112) /
3, and therefore the structure of stud-ied system was reported as an isostructural hcp. TheX-ray diffraction measurement can not provide enoughinformation for determining the molecular orientation.Whereas, in DFT simulation we can observe a muchmore complicated structure which has orientational or-der with a larger primitive unit cell. The stability of P /m structure, which has the largest energy band gapamong other DFT-predicted solid hydrogen molecularstructures, was suggested by hybrid-functional DFT andDMC simulations . Comparing the GW direct andindirect band gaps of these hexagonal structures with ex-periment can yield essential information about the met-allization of them within the pressure range of phase III.The electronic band structure of hcp and P /m struc-tures, which we calculated using DFT-BLYP, are illus-trated in figure 4. The indirect DFT-BLYP energy bandgap of hcp and P /m at P ∼
100 GPa are, 0.56 and4.46 (eV), respectively. There are three high-symmetryk-points of Γ, A , and M which play the crucial rolein the process of band gap closure. We calculated theQP-energies at high-symmetry k-points of A (center of
10 0 10 20 En e r g y ( e V ) (cid:75) M K (cid:75)
A L H A L M H K (cid:75) M K (cid:75)
A L H A L M H K En e r g y ( e V ) FIG. 4. (Color online) The DFT-BLYP band structure of hcp (top panel) and P /m (bottom panel) solid hydrogenwith two and eight H molecules per primitive unit cell, re-spectively.k-point hcp P /m Gap hcp P /m Γ v -3.28 -3.17 Γ 5.28 8.75Γ c M v -3.11 -0.79 M M c A v -13.27 -4.03 A 32.95 10.43 A c L v -4.76 -4.64 L 12.34 15.56 L c G W QP-energies and op-tical band gaps of hcp and P /m structures at ∼
100 GPa.The fundamental electronic band gaps of hcp and P − /m structures are 2.75, and 6.37 eV, respectively. Energies are ineV and with respect to the E F . a hexagonal face), M (center of a rectangular face), andΓ (center of the Brillouin zone) for the top of valenceband and the bottom of conduction band. Table III com-pares the room-temperature G W QP-energies of hcpand P /m structures at pressure ∼
100 GPa. For bothstructures the optical band gap, which is larger than in-direct (thermal) gap, occurs at Γ point.By increasing the pressures the thermal gap closes be-fore the optical gap. Experimentally, it is difficult tooptically measure the thermal band gap due to small-ness of phonon-driven transition which is responsible forabsorption. Most of the previous works focused on themetallization of solid hydrogen by indirect band gap clo- −4−2 0 2 4 6 8 10 Γ M Γ M E n e r g y ( e V )
100 GPa150 GPa200 GPa 300 GPa350 GPa400 GPa −4−2 0 2 4 6 8M A E n e r g y ( e V )
100 GPa150 GPa200 GPa 300 GPa350 GPa400 GPa
FIG. 5. (Color online) G W band structure of P /m structure at 100-400 GPa. The bands are calculated for thehigh-symmetry k-points which are responsible for the bandgap closure. sure. We calculated both the fundamental and opticalband gaps at high-symmetry k-points. The fundamentalband gap of the hcp and P /m structures at pressure ∼
100 GPa is 2.75, and 6.37 eV, respectively. The in-direct energy band gap reduction per pressure is ∼ . ∼ ∼ <
250 GPa.We just estimated the static band gap closure, meaningthat the zero point (ZP) motion contributions, which canreduce the energy band gap ∼ ∼
350 GPa. For instance, a recent observa-tion by Eremets etal . predicts that at pressure of 350-360GPa and T < K , the hydrogen starts to conduct witha semimetallic behaviour . Hence, the structure of high-pressure solid hydrogen has to have a finite energy bandgap at pressures below 350 GPa, which is not the case forthe hcp. Our G W band gap results for the hcp phasepredict that this structure can not be the right candidatefor the structure of phase III of solid hydrogen. The en-ergy band gap and electronic band structure of hcp phasewere studied in reference 16, where an indirect band gapof 3.8 eV was reported at ∼
100 GPa, which is inconsis-tent with our band gap data.Figure 5 illustrates the G W band structure of P /m structure at pressure range of 100 to 400 GPa. The QP-bands are calculated between high-symmetry k-points ofΓ, M , and A which are involved in the process of metal-lization through fundamental band gap closure. Increas-ing the pressures depresses the conduction band with re-spect to the Fermi energy. The valence band maximumoccurs away from the Brillouin zone (BZ) centre Γ but atthe edge of BZ M point. Hence, transitions at the bandedge have to involve a big change in the electron wavevector. Since, optical frequency photons have a smallmomentum k vector, the interband transition requires aphonon to conserve momentum and energy. The first va-lence band maximum below the M -point takes place atthe Γ for P <
250 GPa, whereas by increasing the pres-sure it moves to A -point. For the pressures below 300GPa, the conduction band minimum takes place at thecentre of BZ Γ (Γ c ), but at the larger pressures the con-duction band minimum adopts the A k-point symmetry( A c ). Hence the optical band gap at A is decreased byincreasing the pressure, as it is illustrated in Figure 6.Static enthalpy-pressure phase diagram, which was ob-tained by many-body wavefunction-based technique dif-fusion Monte Carlo (DMC), shows that the P /m struc-ture is stable up to ∼
300 GPa where the energy bandextrema switches the momentum symmetry from Γ to A . This k-point symmetry exchange affects the natureof fundamental band gap. At lower pressures the indi-rect band transition happens through M → Γ, while athigher pressures it takes place via M → A , as it is shownin figure 6.Figure 6 shows the QP-energies, optical and funda-mental band gaps of P /m structure at the pressurewindow of 100 to 400 GPa. The QP-energies are plottedfor the first band below the Fermi energy, valence band(1 v ), and the first band above the Fermi level, conduc-tion band (1 c ), at high-symmetry k-points of Γ, M , and A . The QP-energy is almost a linear function of pressurewith negative slope. The optical and fundamental energyband gaps are also reduced, almost linearly, by increasingthe pressure, except the optical band gap at the Γ whichis almost independent of pressure. The smallest opticalband transition occurs at the BZ-edge point of M . Our G W indirect energy band gaps at 250, 300, and 350GPa agree well with the diffusion quantum Monte Carlo(DMC) energy band gap simulation , which are plottedin the figure 6. The optical band gaps, which were mea-sured by Infra-Red experiment and reported in reference21, are also shown in Figure 6. By a linear extrapolationof band gap as function of pressure, one can see that the M → Γ and M → A indirect band gaps close at 502, and484 GPa, respectively. Whereas, the linear extrapolationof optical band gap of M → M as function of pressuresuggests a zero band gap at 871 GPa.The M → M optical band gap at 400 GPa is ∼ . While,the fundamental band gap at 400 GPa, which is 1.9 eV,differs from the measured band gap by ∼ and Coupled Electron Ion MonteCarlo , is 2-2.5 eV depends on the crystal structure ofsolid hydrogen. It is not clear if the optical band gapreduction due to compressing the system is the same.Taking into account the lattice dynamic effects, we esti- FIG. 6. (Color online) Evolution of QP-energies (top panel),optical and thermal band gaps (bottom panel) as function ofpressure for hexagonal P /m structure. 1 v and 1 c indicesindicate the first valence band below the E F and the firstconduction band above the E F , respectively. The diffusionquantum Monte Carlo (DMC) data and experimental (EXP)optical band gap results are from references 36 and 21, re-spectively. mate that the fundamental energy band gap of P /m structure at 400 GPa is closed, while there is a wideopen optical band gap. The difference between GW pre-dicted indirect and direct band gaps of P /m structureat 400 GPa is 4 eV. The comparison between GW bandgap of P /m structure and experiment indicates thatthis structure may not be the right candidate for thephase III.In principle, for being consistent with experiment, itcan be speculated that the difference between indirectand direct band gaps of probable candidate for the phaseIII should goes to zero. Ideally, this may occur by havingeither large number of hydrogen atoms per primitive unitcell or orientationally distorted molecular configurationto delocalise the electronic density. Therefore, one wouldexpect an electronic bandstructure with nearly flat bandsvicinity the Fermi energy. IV. CONCLUSION
We presented and discussed our implementation of theGW space-time method which allows us to calculate theproperties of extended systems at room-temperature. Wehave calculated QP-energies using two approaches. Inthe first one, named G W AC , the analytic continuationof the self energy is employed. In our second approach G W GF , we obtained the valence- and conductance-bandenergies from the asymptotic decay of the Greens func-tion at long imaginary times. Representation of functionsof imaginary time was performed using two methods: (i)the traditional logarithmic grid and (ii) in the basis ofChebyshev polynomials. We discussed the numerical in-stabilities associated with logarithmic grid and the ad-vantages of using the Chebyshev polynomials. By com-paring and analysing the results yielded by two methods,we found that the Chebyshev representation of the imag-inary time axis improve the accuracy for given memoryusage, and also enables us to control the systematic errorsintroduced by the numerical representation of G and W.To validate our implementation, we calculated the energyband gap of two well-studied systems of crystalline Si andGe, and compared our data with experiment and estab-lished benchmark. Our robust results with controlledsystematic errors within the GW approximation, agreewell with experiment and other published works.The necessity to store the full functions G and W inmemory is an important byproduct of our GW method.It enables constructing on the fly and sampling stochas-tically the higher-order diagrams by means of the es-tablished DiagMC methodology . DiagMC allows tostochastically sum the series in terms of G and W in anumerically exact way, with or without self-consistencywhich can make calculations of materials properties fullycontrolled and reliable.We applied our G W GF technique on the problem ofmetallization of solid hydrogen. We calculated the bandstructure, direct and indirect band gaps of hexagonalsolid hydrogen within pressure range of 100-400 GPa. Weconsidered two hexagonal structures with the symmetryof hcp and P /m with four and sixteen H-atoms perprimitive unit-cell, respectively. The hcp static (not in-cluding ZP effects) fundamental band gap closes at pres-sures below 250 GPa. Including the ZPM reduces themetallization pressure of hcp structure tp pressures be-low 200 GPa. Majority of theoretical and experimentaldata indicate that the solid hydrogen energy band gapat T <
250 K remains open till ∼ GP a . Thus, theright structure of solid hydrogen at 250 GPa has to havea non-zero energy band gap, which is not true for the hcpstructure. Our QP results predict that the hcp structurecan not be the correct candidate for the phase III of solidhydrogen.Our GW fundamental band gap of P /m structureis consistent with the DMC result. However, we found ∼ P /m and the recent experiment . The main sources of this disparity are neglecting the band gap renormal-ization due to lattice dynamics, and the selected hexag-onal structure. It is possible that the realised structureof solid hydrogen in this regime is different. The bandgap reduction due to phonon contribution and ZP vibra-tion can be as large as ∼ ∼ P /m structureand the measured optical band gap of solid hydrogen at ∼
400 GPa. The comparison between the calculated di-rect band gap and experiment suggests that the phaseIII of solid hydrogen can not adopt the P /m structureeither.Thus, provided the GW approximation is adequate,the comparison of our results with experimental dataprovides evidence that rules out the possibility of re-alisation of a hexagonal structure in the phase III ofhigh-pressure solid hydrogen. This prediction is incon-sistent with the recent interpretation of X-ray difrractionmeasurements . According to energy band gap analysisand the suggestion by experiment that the metallizationof hydrogen proceeds within the molecular solid , thelikely candidate for the structure of phase III of solid hy-drogen should have a finite static (ignoring the zero pointmotion) direct energy band gap at ∼
400 GPa. Providedthe size of this gap is comparable with the reduction duezero point vibrations, realisation of a metallic-molecularphase at ∼
400 GPa would be consistent with the exper-imental data.However, in most of the DFT-predicted molecularstructures for phase III, the fundamental energy bandgap is indirect. By increasing the pressure, indirect bandgap closure takes place before direct band gap and themetalliztion occurs through indirect band gap termina-tion. In the very recent study , it has been speculatedthat within the pressure window, where the indirect bandgap is zero but the direct band gap is still open, the den-sity of states near the Fermi energy enlarges by increasingthe pressure and the system becomes a bad metal withproperties similar to a semimetal which has also beenclaimed by experiment . Coorespondingly, the othercandidates for the phase III with indirect fundamentalgap at pressures above 300 GPa have symmetries of C /c , Cmca , P c , and
P bcn with 24 (or 12), 24 (or 12), 48, and48 number of hydrogen atoms per primitive unit cell. Thedensity of energy bands and number of band structure ex-trema near the Fermi energy in these structures are large,and therefore a dense k-mesh is required in GW calcula-tion for accurately obtaining the nature of band gap andtransition between bands. These simulations are consid-erably memory demanding and expensive which call fora separate study.0
V. ACKNOWLEDGEMENT
This work was supported by the Simons Foundation aspart of the Simons Collaboration on the Many-ElectronProblem. ∗ [email protected] † [email protected] E. Wigner and H. B. Huntington, J. Chem. Phys. , 764(1935). N. W. Ashcroft, Phys. Rev. Lett. , 1748 (1968). E. Babaev, A. Sudbo, N. W. Ashcroft, Nature , 666(2004) S. A. Bonev, E. Schwegler, T. Ogitsu, and G. Galli, Nature , 669 (2004). H. K. Mao and R. J. Hemley, Rev. Mod. Phys. , 671(1994). V. L. Ginzburg, Phys. Usp. , 353 (1999). J. M. McMahon, M. A. Morales, C. Pierleoni, and D. M.Ceperley, Rev. Mod. Phys. , 1607 (2012). R. M. Hazen, H. K. Mao, L. W. Finger, and R. J. Hemley,Phys. Rev. B , 3944 (1987) Y. Akahama, M. Nishimura, H. Kawamura, N. Hirao, Y.Ohishi, and K. Takemura, Phys. Rev. B , 060101(R)(2010) M. I. Eremets and I. A. Troyan, Nat. Mater. , 927(2011). R. T. Howie, C. L. Guillaume, T. Scheler, A. F. Goncharov,and E. Gregoryanz, Phys. Rev. Lett. , 125501 (2012). R. T. Howie, I. B. Magdˇau, A. F. Goncharov, G. J. Ack-land, and E. Gregoryanz, Phys. Rev. Lett. , 175501(2014). A. F. Goncharov, I. Chuvashova, C. Ji, and H. Mao, PNAS(2019) DOI: https://doi.org/10.1073/pnas.1916385116 R. T. Howie, P. Dalladay-Simpson, and E. Gregoryanz,Nat. Mater. , 495 (2015). P. Dalladay-Simpson, R. T. Howie, and E. Gregoryanz,Nature , 63 (2016). C. Ji, and et al , Nature , 558 (2019) R. P. Dias, and I. F. Silvera, Science , 715 (2017). X. Liu, P. Dalladay-Simpson, R. T. Howie, B. Li, and E.Gregoryanz, Science , eaan2286 (2017) M. I. Eremets, A. P. Drozdov, P. P. Kong, and H. WangNature Phys. , 1246 (2019). S. Azadi, B. Monserrat, W. M. C. Foulkes, and R. J. Needs,Phys. Rev. Lett. , 165501 (2014) P. Loubeyre, F. Occelli, and P. Dumas, Nature , 631(2020). C. J. Pickard and R. J. Needs, Nat. Phys. , 473 (2007). C. J. Pickard, M. Martinez-Canales, and R. J. Needs, Phys.Rev. B , 214114 (2012). A. F. Goncharov, J. S. Tse, H. Wang, J. Yang, V. V.Struzhkin, R. T. Howie, and E. Gregoryanz, Phys. Rev.B , 024101 (2013). I. B. Magd˘au and G. J. Ackland, Phys. Rev. B , 174110(2013). I. I. Naumov, R. E. Cohen, and R. J. Hemley, Phys. Rev.B , 045125 (2013). M. A. Morales, J. M. McMahon, C. Pierleoni, and D.M.Ceperley, Phys. Rev. B , 184107 (2013). R. C. Clay, J. Mcminis, J. M.McMahon, C. Pierleoni, D.M. Ceperley, and M. A. Morales, Phys. Rev. B , 184106(2014). R. C. Clay, M. Holzmann, D. M. Ceperley, and M. A.Morales, Phys. Rev. B , 035121 (2016). S. Azadi and T. D. K¨uhne, JETP letters , 449 (2012). S. Azadi and W. M. C. Foulkes, Phys. Rev. B , 014115(2013). S. Azadi, and G. J. Ackland, Phys. Chem. Chem. Phys. , 21829 (2017). B. Monserrat, N. D. Drummond, P. D. Simpson, R. T.Howie, P. L. R´ıos, E. Gregoryanz, C. J. Pickard, and R. J.Needs, Phys. Rev. Lett. , 255701 (2018) S. Azadi, W. M. C. Foulkes, and T. D. K¨uhne, New J.Phys. , 113005 (2013). N. D. Drummond, B. Monserrat, J. H Lloyd-Williams, P.L´opez R´ıos, C. J. Pickard, and R. J. Needs, Nat. Comm. , 7794 (2015). S. Azadi, N. D. Drummond, and W. M. C. Foulkes, Phys.Rev. B , 035142 (2017). S. Azadi, and T. D. K¨uhne, Phys. Rev. B , 155103(2019). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). J. P. Perdew, et al , PNAS (2017) DOI:https://doi.org/10.1073/pnas.1621352114 G. J. Ackland, and J. S. Loveday, Phys. Rev. B ,094104 (2020). S. Azadi, and W. M. C. Foulkes, Phys. Rev. B , 245142(2019) Y. Yang, V. Gorelov, C. Pierleoni, D. M. Ceperley, and M.Holzmann, Phys. Rev. B , 085115 (2020) M. P. Surh, S. G. Louie, and M. L. Cohen, Phys. Rev. B , 9126 (1991) W. Luo, S. Ismail-Beigi, M. L. Cohen, and S. G. Louie,Phys. Rev. B , 195215 (2002) B. Holm, and U. von Barth, Phys. Rev. B , 2108 (1998) W. Ku, and A. G. Eguiluz, Phys. Rev. Lett. , 126401(2002) S. V. Faleev, M. van Schilfgaarde, and T. Kotani, Phys.Rev. Lett. , 126406 (2004) S. V. Faleev, M. van Schilfgaarde, T. Kotani, F. L´eonard,and M. P. Desjarlais Phys. Rev. B , 033101 (2006) M. van Schilfgaarde, T. Kotani, and S. V.Faleev, Phys.Rev. Lett. , 226402 (2006) T. Kotani, M. van Schilfgaarde, and S. V.Faleev, Phys.Rev. B , 165106 (2007) M. van Schilfgaarde, T. Kotani, and S. V.Faleev, Phys.Rev. B , 245125 (2006) G. Onida, L. Reining, R. W. Godby, R. Del Sole, and W.Andreoni, Phys. Rev. Lett. , 818 (1995) H. N. Rojas, R. W. Godby, R. J. Needs, Phys. Rev. Lett. , 1827 (1995) I. D. White, R. W. Godby, M. M. Rieger, and R. J. Needs,Phys. Rev. Lett. M. M. Rieger, L. Steinbeck, I. D. White, H. N. Rojas, andR. W. Godby Com. Phys. Comm. , 211 (1999) L. Steinbeck, A. Rubiob, L. Reiningc, M. Torrent, I. D.Whited, and R. W. Godby, Com. Phys. Comm. , 105(2000) E. Gull, S. Iskakov, I. Krivenko, A. A. Rusakov, and D.Zgid, Phys. Rev. B , 075127 (2018) Jia Li, Markus Wallerberger, Naoya Chikano, Chia-NanYeh, Emanuel Gull, and Hiroshi Shinaoka Phys. Rev. B101, 035144 N. V. Prokof’ev, and B. V. Svistunov, Phys. Rev. Lett. ,2514 (1998) E. Kozik, K. Van Houcke, E. Gull, L. Pollet, N. V.Prokof’ev, B. V. Svistunov, and M. Troyer, EPL (Euro-physics Letters) , 10004 (2010) K. Van Houcke, E. Kozik, N. V. Prokof’ev, B. V. Svis-tunov, Physics Procedia , 95 (2010) A. S. Mishchenko, N. Nagaosa, and N. Prokof’ev, Phys.Rev. Lett. , 166402 (2014) K. Van Houcke, F. Werner, E. Kozik, N. Prokof’ev, B. Svis-tunov, M.J.H. Ku, A.T. Sommer, L.W. Cheuk, A. Schi-rotzek, and M.W. Zwierlein, Nature Physics , 366 (2012) Y. Deng, E. Kozik, N. V. Prokof’ev, and B. V. Svistunov,EPL (Europhysics Letters) , 57001 (2015) K. Chen, and K. Haule, Nature Comm. , 3725 (2019) X. Gonze, R. Stumpf, and M. Scheffler, Phys. Rev. B ,8503 (1991) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni,I. Dabo, et al. , J. Phys.: Condens. Matt. , 395502 (2009). C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B , 785(1988). S. Azadi, R. Singh, T.D. K¨uhne, J. Comp. Chem. , 262(2018) V. Gorelov, M. Holzmann, D. M. Ceperley, and C. Pier-leoni, Phys. Rev. Lett. , 116401 (2020)71