Electrical conductivity of high-pressure liquid hydrogen by quantum Monte Carlo methods

0.2 MB·English

by  Fei Lin

Electrical conductivity of high-pressure liquid hydrogen by quantum Monte Carlo
methods
Fei Lin,1 Miguel A. Morales,1 Kris T. Delaney,2 Carlo Pierleoni,3 Richard M. Martin,1,4 and D. M. Ceperley1
1Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
2Materials Research Laboratory, University of California, Santa Barbara, CA 93106, USA
3CNISM and Department of Physics, University of L’Aquila, L’Aquila, Italy
4Department of Applied Physics, Stanford University, Stanford, CA 94305, USA
(Dated: January 22, 2010)
WecomputetheelectricalconductivityforliquidhydrogenathighpressureusingquantumMonte
0
Carlo. ThemethodusesCoupledElectron-Ion MonteCarlo togenerateconfigurationsofliquidhy-
1
drogen. Foreachconfiguration,correlatedsamplingofelectronsisperformedinordertocalculatea
0
set of lowest many-body eigenstates and current-currentcorrelation functions of thesystem, which
2
aresummed overin themany-bodyKuboformula togiveACelectrical conductivity. Theextrapo-

n lated DC conductivityat 3000 K for several densities shows a liquid semiconductor to liquid-metal
a transition at high pressure. Ourresults are in good agreement with shock-wave data.
J

2 PACSnumbers: 71.22.+i,72.15.Cz,02.70.Ss
2

] Liquidhydrogenathighpressurehasbeen the subject ductivity is
l
e ofintenseexperimentalandtheoreticalresearch,because
r- of its special place in the periodic table and its cosmic σ (ω)= 2πe2(1−e−β~ω) |hk| pi|ni|2×
st abundance. Understandingits behaviorunderhighpres- αα m2ωΩZ Xk,n Xi α
. sureandhightemperatureis importantforrevealingthe
at properties of giant planets such as Jupiter and Saturn. e−βEnδ(Ek−En−~ω), (1)
m
Metallization of liquid hydrogen at high pressure is of
where α = x,y,z; e (m) is electron charge (mass), Ω is
- particular interest. Experiments using shock waves have
d the volume, β is the inverse temperature, ω is the fre-
found a metallization transition[1]; at a pressure of 140
n quency, and pi is the α component of the momentum
o GPa and temperature of 3000 K, liquid hydrogen has α
operatorfor the ithelectron. |niand E aremany-body
c been reported to turn from an insulating fluid into a n
eigen-states and eigen-energies of the Hamiltonian, and
[
2 mcmet)a−l1li.cHoonweewveirth, thaeDorCeticcoanlldyu,cstuicvhitya mofetaabloliuztat2io0n00pr(oΩ- Z = ne−βEn the partition function.
ThPekeychallengeinevaluatingtheKuboformula(Eq.
v cess has not been very well understood. Mean-field den-
1)istodeterminethesumoverallmany-bodyeigenstates
8 sity functional theory (DFT) calculations [2] and Quan-
4 of the system. We compute a number of the lowest-
tumMolecularDynamics(QMD)calculations[3,4]have
2 energyeigenstates,andtherelevantmatrixelements,us-
been used to calculate the electrical conductivity in liq-
2 ingtheCFQMCmethod[6],whichwenowexplain. Con-
. uidhydrogen,butthesemethodsneglectcorrelationsbe-
9 sider the subspace spanned by a set of M many-body
0 tweenelectrons. Suchcorrelationsarelikelytobeimpor- basis states f (R), i.e., Φ (R) = M d f (R), where
9 tant for the accurate determination of transport proper- R representsjelectronic coiordinatePs.j=T1heijgjround basis
0 ties. In this letter, we propose and apply an alternative state, f (R),has SlaterdeterminantsofKohn-Shamsin-
: ab-initio method which combines the Coupled Electron- 0
v gleelectronorbitalstimes aJastrowpaircorrelationand
i Ion Monte Carlo (CEIMC) method [5] with the Correla-
X backflow corrections[8, 13]. For the excited basis states
tion Function Quantum Monte Carlo (CFQMC) method
f (R),j > 0, we use low-lying particle-hole excitations
ar [6, 7]. wjith respect to the ground state f (R).
0
Withinthissubspace,upperboundstotheexacteigen-
Within the CEIMC approach [5, 8], the electrons and
valuescanbefound,byminimizingtheRayleighquotient
protons are simulated with separate but coupled Monte
with respect to d ,
Carlo simulations, taking advantage of the separation ij
of time scales in the Born-Oppenheimer approximation.
dRΦ∗(R)HˆΦ (R)
The Born-Oppenheimer energy surface for the protons Λ = i i , (2)
is determined using reptation quantum Monte Carlo i R dRΦ∗i(R)Φi(R)
(RQMC) [9]. After the proton system reaches equilib- R
rium, we record uncorrelated samples of proton config- where Hˆ is the electronic Hamiltonian. Furthermore, an
urations, which are subsequently used to determine the improved basis set {f˜i(R)} can be obtained by applying
electrical conductivity. an imaginary-time projection to the basis set:
The many-body Kubo formula [10] for electrical con- f˜(R)=e−tHˆ/2f (R), (3)
i i
2
15000 0.6 Note that we have used Eq. (3) and inserted complete
(a) (b)
t=0.06 Ha-1 0.4 sets of basis states {|Ri}. These matrix elements are
t=0.22 Ha-1 ω) calculated all at once with RQMC.
t=0.62 Ha-1 S(
-1m)10000 tt==01..8022 HHaa--11 0.2 eigSeonlfvuinncgtiEoqn.s(t4h)aytiealrdesrtehqeumireadnyf-obrotdhyeeKiguebnvoafluorems aunlad.
c 0
σ(ω) (Ω 5000 0 (c) 10ω (eV)t2=01.02 Ha-130-16.8 Tathehaclufrroefntth-ceurprreonjtecctoirornelatitmioen tfu/2ncitnionRsQaMreCco[m6].putFeid-
E (Ha)n EEE012 -17.2 nbyallay,GtahuessδiafnunfucnticotnionindtuhreinKguthbeonfourmmeurilcaailscablrcoualdateinoend,
32 atoms, T=3000K EE34 whose width is of the same order as the observed many-
00 2 4 6 8 10 0 0.3 0.6 0.9 1.2 1.5-17.6 body energy spacing.
ω (eV) t (Ha-1) We use the ground state as guiding wavefunction dur-
ingRQMCsimulations[17]. Themethodisabletodeter-
FIG. 1: (Color online). Electrical conductivity (a), the cor-
minethelowest-lyingstatesofthesystem,typicallyfewer
responding sum rule (b), and five lowest energy levels (c) as
afunction of imaginary-time projection t for rs=1.40. Note than 50 in our simulations, after which the calculation
thedecreaseofσ(ω)atlowfrequencyisanartifactduetothe becomeslessefficientduetotheincreasedfluctuationsin
gap in a finite system. The simulation has M = 24 excited the matrix elements in Eqs. (5-6), a problem very much
states. Theresultsareaveragedover108twistangles(seeFig. relatedtothe well-knownfermionsignproblemofQMC.
3) and10 decorrelated proton configurations drawn from the An analysis shows[6] that the Monte Carlo (MC) effi-
thermal Boltzmann distribution sampled byourCEIMC cal- ciency must decrease with projection time as exp(−αt),
culations. Lines are just guides to theeye.
making convergence in t problematical: we call this the
“efficiency problem”. However, it is important to note
that as the basis is increased in size, the projected low
(a) (b) 0.6
energy states become more accurate [6]. This suggests
15000 MMM===123246 ω)S( 00..24 tthhaeteffiweciesnhcoyulids ninoctlusedveeraeslymraednyuceexdc.itations as long as
-1m)
Ω c10000 0 10 ω (2e0V) 30 400 15000
σ(ω) (500000 322 atom4s, T=36000K8 10E (Ha)n E-E110212 2M4 2436 0000....00114826 (3c6)EEE012 ---111776...418 -1σ(ω) (Ω cm)150000000 Γ10 p8o tiwntist average
ω (eV) M
T=3000K, r=1.50
s
FIG. 2: (Color online). Electrical conductivity (a), the cor- 32 atoms, t=1.02 Ha-1
0
responding sum rule (b), and first three lowest energy levels 0 2 4 6 8 10
ω (eV)
with energy difference E1−E0 in the inset as a function of
M, the number of trial wavefunctions in the basis set. Here
FIG. 3: (Color online). Electrical conductivity calculated at
t=1.02 Ha−1 for rs=1.40. theΓ point compared with thetwist-averaged one.
where t is the projection time. Replacing {f (R)} with Allsimulationsreportedherewereat3000K,thetem-
j
perature of the shock experiments measuring conductiv-
{f˜(R)} and minimizing Λ with respect to d , one ob-
j i ij ity[1]. WefirstinvestigatetheroleofRQMCimaginary-
tains the many-body generalized eigenvalue equation:
time projection t. Fig. 1 (a) shows the electrical con-
M M ductivity of liquid hydrogen as a function of t. Defining
Xj=1Hij(t)dkj(t)=Λk(t)Xj=1Sij(t)dkj(t), (4) Slim(ωω)→=∞Sπ2e(mnωe)R=0ωσ1(ω[1′0)]d.ωH′,etrheenceonisduelcetcitvriotyn sduemnsirtuy.leIins
the present method, we cannot include all many-body
where Λ (t) is the kth eigenvalue, and the Hamiltonian
k
HandtheoverlapmatricesSaredefinedwithrespectto states, so we do not expect the sum rule to be satisfied.
the projected basis set {f˜(R)} as However, it provides an indication of calculation qual-
i
ity. The main measure of appropriate projection time is
H = dR dR f (R )hR |Hˆe−tHˆ|R if (R ), (5) the convergenceof low-frequency σ(ω). We can see from
ij Z 1 2 j 2 2 1 i 1 Fig. 1 (a) that low-frequency σ(ω) curves converge be-
tween the interval t ∈ [0.82,1.02]Ha−1. Here we choose
Sij = Z dR1dR2fj(R2)hR2|e−tHˆ|R1ifi(R1). (6) the upper bound oft=1.02Ha−1, since we want to have
3
12000 The corresponding σ(ω) curve is also suppressed in the
r=1.40
32 atoms, T=3000K s
10000 rs=1.50 low frequency region. See Fig. 2 (a). We find that M =
r=1.60
s
r=1.65 24isagoodcompromisebetweenaccuracyandefficiency.
-1m)8000 esxp. WealsonoticefromFig. 2(b)thatS(ω)→0.6asω →∞
Ω c6000 forM =24;while for M =12 and36,the corresponding
σ(ω) (4000 sum approaches 0.5.
We note that our simulations are subject to finite-size
2000 errors in evaluating thermodynamic properties. Such ef-
fects can be reduced by using twist-averaged boundary
0
0 2 4 6 8 10
ω (eV) conditions [11]. In Fig. 3 we compare σ(ω) calculated
0.01 at the Γ point of the Brillouin zone (BZ) to the one av-
r=1.40 32 atoms, T=3000K
s eraged over a 6 × 6 × 6 grid in the twist angle space
3Bohr/spin)00..000068 rrrsss===111...566005 (tthheecBonZd)u.cWtiveitfiyndvatluhea.t HΓopwoeivnetr,safimniptlei-nsgizeoveeffreescttismhaatvees
OS (states/eV/0.004 nostboutsdebirevedee.dn evTnahtniiirssehlwyinieglllimonfointσa(htωea)dppaaetnndsmfaorarellarωesspufffioonrcsiaieblnlltelsyyfosrltaetrmhgees
D0.002
system in the metallic phase (e.g. at r ≤ 1.50). A
s
common procedure[3, 12, 14] to estimate the DC con-
0
0 2 4 6 8 10
E-E (eV) ductivity is to discarda few data points and extrapolate
0
the higher frequency data to ω = 0 . For systems in-
FIG. 4: (Color online). Electrical conductivity (scattered side the metallic phase one can use the Drude formula
points) along with the Drude fits (solid lines) (upper panel) σ(ω)=σ /(1+ω2τ2) to estimate the DC conductivity
DC
andDOS(lower panel) asafunction ofdensity forliquid hy-
σ and the electronic relaxation time τ.
drogen. Thered dot at ω=0 indicates theexperimentalDC DC
conductivity value after the liquid hydrogen is metallized by
5000
highpressureat3000K[1]. Thedropofelectricalconductiv-
ity to zero at zero frequency for rs = 1.40, 1.50 and 1.60 is 4000 54 atoms, T=3000K rrss==11..5605
dueto finite-size effects. r=1.75
-1m)3000 esxp.
Ω c
ω) (2000
σ(
as large t as possible before the efficiency is reduced to
1000
an intolerable level. Note that the decrease in σ(ω) at
ω <2eVisanartifactduetothe finitenumberofatoms
0
0 2 4 6 8 10
in the simulation, as we discuss below. Fig. 1 (b) also ω (eV)
shows that the best sum rule satisfaction is achieved at
t = 1.02Ha−1, where S(ω) ∼ 0.6 as ω → ∞. As a fur- FIG. 5: (Color online). Electrical conductivity (scattered
points) and Drude fits (solid lines) as a function of density
ther test of the t value, we show in Fig. 1 (c) five lowest
for liquid hydrogen. The red dot at ω = 0 denotes the ex-
energy states as a function of projection time. We see
perimental DC conductivity value after the liquid hydrogen
that indeed after t = 1.02Ha−1, the statistical variances
is metallized by high pressure [1].
increase substantially, and the estimates of the lowest
energies are pushed down due to the noise [7].
Fig. 2 shows (a) the electrical conductivity, (b) the We now proceed to study the electrical conductivity
corresponding S(ω), and (c) the first 3 energy levels at at various densities at T= 3000 K, where the metalliza-
fixed projection time as a function of M, the size of the tion of liquid hydrogen has been observed [1]. We run
basis set. As before, we choose M as large as possible CEIMC simulations at r = 1.40, 1.50, 1.60, and 1.65
s
whilestillobtainingareasonableMCefficiency. Tocheck for 32 protons and electrons. Results (scattered points)
the efficiency, we collect the same amounts of stochastic along with the corresponding Drude fits (solid lines) are
data for the Hamiltonian and overlap matrices in Eq. 5 shownonthe upper panelofFig. 4 andmany-bodyden-
and 6, solve the generalized eigenvalue equation Eq. 4, sity of states (DOS) on the lower panel, where all the
andlookatthe energygapbetweengroundstateE and energy values are calculated with respect to the ground-
0
thefirstexcitedstateE . Asuddenincreaseoftheenergy state energy of the configuration. We can see that DC
1
gap indicates the reduction of the efficiency, since noise conductivity decreases from about 8600 (Ω cm)−1 (see
can push down the most the lowestenergy level [7]. The below)atr =1.40to about3500(Ω cm)−1 atr =1.60
s s
inset of Fig. 2 (c) shows an increased slope for E −E and to nearly zero at r = 1.65. The dot at ω = 0
1 0 s
curveatM =36,indicatingsuchareductionofefficiency. in Fig. 4 shows the experimental DC conductivity after
4
metallization [1]. The conductivity points for r = 1.65 a complete understanding of liquid hydrogen electronic
s
near zero frequency has typical features for liquid semi- transport properties. It is also very promising to extend
conductors, while a 2000 (Ω cm)−1 is a typical liquid our method to other systems, such as lattice models.
metalvalue[15]. Thus,basedon32-atomsimulationswe
see a liquid metal to a liquid semiconductor transition F.L. thanks Michele Casula for helpful discussions.
in a density region of 1.60 < rs < 1.65 at T=3000 K C.P. thanks the Institute of Condensed Matter The-
forliquidhydrogenunderhighpressure,whichiscloseto ory at the University of Illinois at Urbana-Champaign
the experimental transition density of rs =1.62. Such a for a short term visit. This research was sponsored in
metal-semiconductortransitionisalsohintedintheDOS part by the National Nuclear Security Administration
figure. Going from rs = 1.40 to rs = 1.65, one can see under the Stewardship Science Academic Alliances pro-
the gradual opening of an energy gap. gramthroughDOE GrantDE-FG52-06NA26170;Minis-
We further look at the metallic behavior of the liquid terodell’Istruzione,dell’Universitae dellaRicercaGrant
hydrogenatrs =1.40. TheDrudefitgivesσDC =8620± PRIN2007 (to C.P.). Computer time was provided by
1000(Ωcm)−1andτ =2.0±0.2a.u. or3.1±0.3×10−16s. NCSAattheUniversityofIllinoisatUrbana-Champaign
Aroughestimateofelectronmeanfreepath(l)thengives and CINECA, Italy.
l ∼ 1.2 ˚A, which is about 1.5 times the average proton-
proton distance (1.4 Bohr) [8], necessary for the system
being metallic. For r = 1.50 and 1.60 we also fit σ(ω)
s
points to the Drude formula. The fitted curves become
flatterasthedensitydecreases,whichsignalsthatatr = [1] S. T. Weir, A. C. Mitchell, and W. J. Nellis, Phys. Rev.
s
1.60, liquid hydrogen has a small electronic relaxation Lett. 76, 1860 (1996); W. J. Nellis, S. T. Weir, and A.
C. Mitchell, Phys.Rev. B 59, 3434 (1999).
time, making it a bad metal, similar to liquid carbon
[2] O.Pfaffenzeller and D.Hohl,J.Phys.: Condens.Matter
[12]. 9, 11023 (1997).
Finally, we explicitly test the importance of finite-size [3] L. A. Collins et al.,Phys.Rev.B 63, 184110 (2001).
effects by performing simulations for 54-atom cells as [4] B.Holst,R.Redmer,andM.P.Desjarlais, Phys.Rev.B
shown in Fig. 5. The energy gap at r = 1.65 that is 77, 184201 (2008).
s
present in the 32 atom simulation now disappears. The [5] C. Pierleoni, D. M. Ceperley, and M. Holzmann, Phys.
liquid semiconductor [at r = 1.75 the extrapolated DC Rev. Lett. 93, 146402 (2004); C. Pierleoni and D. M.
conductivity is around 130s0 (Ω cm)−1 less than a typi- Ceperley, Chem. Phys. Chem. 6, 1872 (2005); ibid, Lec-
ture Notes in Physics 703, 641 (2006).
cal liquid metal value [15, 16], hence the name] to liquid [6] D. M. Ceperley and B. Bernu, J. Chem. Phys. 89, 6316
metal transition density is estimated to be in the range
(1988); B. Bernu,D.M. Ceperley, and W. A.Lester Jr.,
1.65≤ rs ≤1.75, which is again close to the experimen- J. Chem. Phys. 93, 552 (1990).
tal density of r = 1.62 [1]. Further increase of system [7] B. L. Hammond, W. A. Lester Jr., and P. J. Reynolds,
s
size is not possible at present, but previous CEIMC cal- Monte Carlo methods in ab initio quantum chemistry,
culations have found that the ground-state energy with Chap. 6, World Scientific, Singapore (1994).
[8] K. T. Delaney, C. Pierleoni, and D. M. Ceperley, Phys.
54and108atomsareverycloseiftwist-averagedbound-
Rev. Lett.97, 235702 (2006).
ary conditions are used. Therefore, we expect that the [9] S. Baroni and S. Moroni, Phys. Rev. Lett. 82, 4745
liquid metal transitiondensity determined above is close
(1999).
to the thermodynamic limit. A more definitive answer [10] R. Kubo,J. Phys.Soc. Japan 12, 570 (1957).
will require significant additional calculations. [11] C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E
Wehaveproposedamethodtocalculatethefrequency 64, 016702 (2001).
dependentelectricalconductivityusingCEIMCwithcor- [12] G.Galli,R.M.Martin,R.Car,andM.Parrinello,Phys.
Rev. B 42, 7470 (1990).
relatedMonte Carlo sampling and the many-body Kubo
[13] C.Pierleoni,K.T.Delaney,M.A.Morales,D.M.Ceper-
formula. The method is able to estimate some of the ley, and M. Holzmann, Comp. Phys. Comm. 179, 89
lowest-lyingenergiesandtheircorrespondingoverlapma-
(2008).
trices, and hence it is suitable for the extrapolation of [14] G. Galli, R. Q. Hood, A. U. Hazi, and F. Gygi, Phys.
DC conductivity from the AC conductivity. With this Rev. B 61, 909 (2000).
methodwestudytheDCconductivityofliquidhydrogen [15] M.Shimoji,LiquidMetals,p.381,AcademicPress,Lon-
as a function of density at 3000 K and at high pressure, don (1977).
[16] N. F. Mott, Philos. Mag. 24, 2 (1971).
andshowthe metallization. OurDC conductivity values
[17] It has been shown that an improved guiding wavefunc-
at the metallization point is in good agreementwith the
tion exists in the Bosonic excited state simulations [6].
shock-waveexperiments [1]. Furthermore,our metalliza-
We tested this approach by including excited states,
tiondensity(rs ∼1.65)isclosetotheexperimentalvalue whichareassignedwithsomeappropriateweights,inthe
ofr =1.62. Inthe future, itwillbe interestingto apply RQMC guiding wavefunction, but the method becomes
s
the method to more densities and temperatures to have unstable, possibly dueto theFermion sign problem.