Retrieve Search Browse Feedback HelpHome Retrieve Search Browse Feedback Help
Previous article Next article Up to abstract
Calculations of the condensate fraction and momentum distribution of liquid 4He are carried out using variational ground-state wave functions containing two-body (Jastrow) and three-body correlations. These wave functions give a satisfactory description of the equation of state of the liquid at T=0. The calculations are performed within the scheme of hypernetted-chain equations, using the scaling approximation for evaluating the contributions of the elementary diagrams. Results are reported at densities rho =0.365 sigma -3, 0.401 sigma -3, and 0.438 sigma -3, and compared with momentum distributions obtained from experimental data and Green's-function Monte Carlo calculations.
©1985 The American Physical Society.
PACS: 67.40.-w
The theoretical and experimental studies of condensate fraction and momentum distribution of atoms in liquid ^{4He}are a longstanding problem of fundamental interest. These studies started four decades ago after London^{1}proposed a connection between the lambda transition and the Bose Einstein condensation. In recent years, a variety of experimental techniques^{2-4}have been employed to measure the momentum distribution in helium liquids. At the equilibrium density, analysis of the experimental data^{2,3}gives condensate fraction n_{0}app 9 ndash 13^_% at T=0. The value of n_{0_}obtained from the Green's-function Monte Carlo (GFMC) calculations^{5}using the HFDHE2 interatomic potential of Aziz et al.^{6}is found to be 9%, which is at the lower limit of the experimental value. The value of the average kinetic energy per atom has also been inferred from the experimental data.^{3}The experimental value of 14.0 +- 0.5 K is in reasonable agreement with that obtained in GFMC (14.5 K) and variational^{7}(14.8 K) calculations.
The momentum distribution n(k) is generally obtained by calculating the one body density matrix rho ( r_{{1}1 prime} )_, and taking its Fourier transform.^{8}It is known^{9}that to obtain an accurate evaluation of the rho ( r_{{1}1 prime} )_, with chain summation methods, it is necessary to sum the various series of elementary diagrams. Fabrocini and Rosati^{10}used the method of interpolating between the hypernetted chain (HNC) and Percus Yevick (PY) equations to calculate n(k). In this paper we generalize the scaling method (HNC/S) developed earlier^{11,7}for the calculation of two- and three particle distribution functions, to sum the series of elementary diagrams in the calculation of rho ( r_{{1}1 prime} )_.
During the course of this work Puoskari and Kallio^{12}(PK) have also calculated the rho ( r_{{1}1 prime} )_ with a very similar HNC/S method. They have used Jastrow variational wave functions
PSI_{J}= ^sum from {i < j} ^f ( r_{ij}) , (1.1)
obtained with the Lennard Jones potential, whereas we use the more realistic wave functions
PSI_{v}= ^prod from {i < j} ^f ( r_{ij}) ^prod from {i < j < k} ^f_{3}( r_{ij}, r_{jk}, r_{ki}) (1.2)
containing three particle correlations. Our wave functions are obtained from the Aziz potential^{6}which gives a better equation of state than the Lennard Jones potential.^{5}
The HNC/S equations of the rho ( r_{{1}1 prime} )_ of a Jastrow wave function are given in Sec. II. The small differences between PK and our formalism are discussed at the end of this section. The results of HNC/S calculations for the McMillan Jastrow wave function are reported; these are in excellent agreement with the exact Monte Carlo results.^{13}
The HNC/S equations for wave functions containing two- and three body correlations are given in Sec. III. The n(k) calculated from the best available variational wave functions (those of Ref. 7) is reported in Sec. IV. Our values of n_{0_}are similar to those obtained with the GFMC method, and a little below those obtained from the analysis of neutron scattering data. The density dependence of n_{0_}and n(k) is discussed.
The one particle density matrix is defined as
rho ( r vec_{{1}1 prime} ) =N {int^PSI^{star}( r vec_{{1}prime} , r vec_{2},..., r vec_{N}) PSI ( r vec_{1}, r vec_{2},..., r vec_{N}) d^{3}r_{2}... d^{3}r_{N}}over {int^| PSI ( r vec_{1},..., r vec_{N}) |^{2}^d^{3}r_{1}... d^{3}r_{N}} , (2.1)
where N is the total number of particles in the liquid and r vec_{{1}1 prime} = r vec_{1}- r vec_{{1}prime}_. In homogeneous liquids (N -> inf, at constant density) the rho ( r_{{1}1 prime} )_ is a function of | r vec_{{1}1 prime} |_ only. Its Fourier transform gives the occupation probability for single particle states with momentum k:
n(k) = ^int^rho (r) e^{{i}k vec cdot r vec} d^{3}r . (2.2)
In Bose systems n(k) is written as
n(k) = N n_{0}delta_{k,0}+ ^int^[ rho (r) - rho (r -> inf ) ] e^{{i}k vec cdot r vec} d^{3}r , (2.3)
where n_{0_}is the condensate fraction, and the second term represents the momentum distribution for k > 0.
By expanding the many body integrals in Eq. (2.1), in powers of the functions
h( r_{ij}) = f^{2}( r_{ij}) -1 ,
zeta ( r_{ij}) = f( r_{ij}) - 1 , N (2.4)
the following irreducible structure for the one particle density matrix has been obtained:^{8}
rho (r) = rho n_{0}exp [ N_{ww}(r) + E_{ww}(r)] . (2.5)
Here rho is the number density of the liquid and N_{ww}(r)_ is the sum of nodal diagrams given by the convolution
N_{ww}( r_{{1}1 prime} ) = rho ^int mark ^ d^{3}r_{2}[ g_{wd}( r_{12}) -1]
lineup times [ g_{wd}( r_{{1}prime 2} ) - N_{wd}( r_{{1}prime 2} ) -1] . (2.6)
The functions g_{wd}(r)_ and N_{wd}(r)_ are obtained by solving the HNC equations
g_{wd}( r_{12}) = f( r_{12}) exp [N_{wd}( r_{12}) + E_{wd}( r_{12})] , (2.7)
N_{wd}( r_{12}) = rho ^int mark ^ d^{3}r_{3}[ g_{wd}( r_{13})-1]
lineup times [ g_{dd}( r_{32}) - N_{dd}( r_{32})-1] . (2.8)
The functions g_{dd}(r)_ and N_{dd}(r)_ are the familiar pair distribution and nodal functions denoted by g(r) and N(r) in Ref. 7.
The condensate fraction n_{0_}is given by^{8}
n_{0}= exp ( 2R_{w}- R_{d}) , (2.9)
where
R_{w}= mark ^ rho ^int^d^{3}r [ g_{wd}(r) - 1- N_{wd}(r)- E_{wd}(r)] L
lineup - cf12 rho ^int ^ d^{3}r [ g_{wd}(r) -1] [ N_{wd}(r) + 2E_{wd}(r) ] + E_{w} .
The R_{d_}is obtained by replacing all the w subscripts in the above equation by d.
The E_{ww}(r)_, E_{wd}(r)_, E_{dd}(r)_, E_{w_,}and E_{d_}represent the contribution of elementary diagrams. There is no analytic method available to evaluate these functions. The lowest order (HNC/0) approximation is characterized by setting all these functions equal to zero. In the HNC/4 approximation the E apos sare calculated from the four point elementary diagrams shown in Fig. 1. The dashed straight and wavy lines in these diagrams represent, respectively, g_{dd}-1_ and g_{wd}-1_. The solid and open circles represent internal and external points. The contributions of these diagrams are given by E_{ww,4}( r_{{1}1 prime} ) = cf12 rho^{2}^int mark ^ (g_{wd}-1)_{12}( g_{wd}-1)_{13 }(2.10)
lineup times ( g_{wd}-1)_{{1}prime 2} ( g_{wd}-1)_{{1}prime 3}
lineup times ( g_{dd}-1)_{23}d^{3}r_{2}d^{3}r_{3} , (2.11a)
E_{wd,4}( r_{12}) = cf12 rho^{2}^int mark ^ (g_{wd}-1)_{13}( g_{wd}-1)_{14
}
lineup times ( g_{dd}-1)_{23}( g_{dd}-1)_{24
}
lineup times ( g_{dd}-1)_{34}d^{3}r_{3}d^{3}r_{4} , (2.11b)
E_{dd,4}( r_{23}) = cf12 rho^{2}^int mark ^ (g_{dd}-1)_{24}( g_{dd}-1)_{25
}
lineup times ( g_{dd}-1)_{34}( g_{dd}-1)_{35
}
lineup times ( g_{dd}-1)_{45}d^{3}r_{4}d^{3}r_{5} , (2.11c)
E_{w,4}= rho over 3 ^int^( g_{wd}-1)_{12}E_{wd,4}( r_{12}) d^{3}r_{12} , (2.11d)
E_{d,4}= rho over 3 ^int^ ( g_{dd}-1)_{12}E_{dd,4}( r_{12}) d^{3}_{12} . (2.11e)
In the HNC/S method^{11}we assumed that
E_{dd}(r) approx (1+ s_{dd}) E_{dd,4}(r) . (2.12)
This approximation is based on the observation that the contributions of higher order elementary diagrams, containing five or more points, have approximately the same spatial behavior as that of E_{dd,4}(r)_. The scaling constant s_{dd_}is determined by fulfilling identities that are valid when the elementary diagrams are summed to all orders.
We note in Fig. 2 that the structure of g_{wd_}is similar to that of g_{dd_,}and the topological structure of five- and higher body E_{ww_}and E_{wd_}diagrams is identical to that of E_{dd_}diagrams. In fact E_{wd_}( E_{ww})_ diagrams are obtained from the E_{dd_}diagrams by replacing all the g_{dd}-1_ bonds from one (both) external point by g_{wd}-1_ bonds. We thus expect the higher body wd and ww elementary diagrams to have the same spatial behavior as that of E_{wd,4_}and E_{ww,4_,}respectively, and approximate the total E_{wd_}and E_{ww_}as follows:
E_{wd}(r) approx (1+ s_{wd}) E_{wd,4}(r) ,
E_{ww}(r) approx (1+ s_{ww}) E_{ww,4}(r) . N (2.13)
The one point elementary diagrams give a factor exp (2E_{w}- E_{d})_ in the expression for n_{0_}[Eq. (2.9)]. This factor is neglected in all the earlier work. In HNC/4 approximation (2E_{w,4}- E_{d,4})_ is of order 0.01, and thus this factor can be safely neglected. We can, of course, scale the E_{d,4_}and E_{w,4_}to approximate the total E_{d_}and E_{w_.}The E_{d_}and E_{w_}diagrams are obtained by dressing the E_{dd_}and E_{wd_}diagrams with g_{dd}-1_ and g_{wd}-1_, respectively, and integrating over one particle. However, the symmetry factors of the one- and two point diagrams are generally different, and hence the scaling factors of the one- and two point diagrams are different, but related. This is very similar to the scaling of three point Abe diagrams A, and that of the E_{dd_.}In Ref. 11 we argued that the ratio of s_{dd_}to s_{a_,}the scaling constant for Abe diagrams, is app 2. Similar arguments suggest that the ratios s_{d}/ s_{dd_}and s_{w}/ s_{wd_}should be app case 3 over case 2. We, hence, approximate the one point elementary diagrams as follows:
E_{d}approx (1+ case 3 over case 2 s_{dd}) E_{d,4} ,
E_{w}approx ( 1+ case 3 over case 2 s_{wd}) E_{w,4} . N (2.14)
We find that 2E_{w}- E_{d_}is app 0.03, and thus not very important.
The kinetic energy and normalization identities are used to determine the scaling constants s_{dd_,}s_{wd_,}and s_{ww_.}These are
T_{roman}MD = T_{roman}JF = T_{roman}PB , (2.15)
T_{roman}MD = ^int^ {d^{3}k} over {( 2 pi )^{3}rho} {hbar^{2}}over 2m k^{2}n(k) , (2.16)
T_{roman}JF = - {hbar^{2}}over 4m rho ^int^g_{dd}(r) ( {del^{2}f} over f - {( del f )^{2}}over f^{2}right ) d^{3}r , (2.17)
T_{roman}PB = mark ^ - {hbar^{2}}over 2m rho ^int^ g_{dd}(r) {del^{2}f} over f d^{3}r
lineup - hbar^{2}over 2m rho^{2}^int mark ^ g_{3}( r_{12}, r_{13}, r_{23}) {del vec_{1}f( r_{12}) cdot del vec_{1}f( r_{13})} over {f( r_{12}) f ( r_{13})}
lineup times d^{3}r_{12}d^{3}r_{13} , (2.18)
for the kinetic energy per atom. The subscripts MD, JF, and PB denote momentum distribution, Jackson Feenberg, and Pandharipande Bethe expressions. The normalization condition
int^ {d^{3}k} over {(2 pi )^{3}rho} n(k) =1 , (2.19)
is equivalent to
rho ( r_{{1}1 prime} =0) = rho . (2.20)
In Refs. 11 and 7 the identity T_{roman}PB = T_{roman}JF_ is used to determine the s_{dd_.}Here we use T_{roman}MD = T_{roman}JF_, and normalization Eq. (2.20) to determine the s_{wd_}and s_{ww_.}The numerical work is simplified by assuming that n(k) is exponential for k > 3.5 A ang^sup -1,
n(k > 3.5 Aang^sup -1 ) = exp [ alpha (k-3.5)]n(k=3.5) . (2.21)
The constant alpha is determined by fitting the ln [n(k)] in k=3 to 3.5 A ang^sup -1 interval to a straight line. This approximation has been verified to be good up to app 4 A ang^sup -1\p (see Fig. 7 in Sec. IV). The k > 3.5 A ang^sup -1 region gives less than 5% (1%) contribution to T_{roman}MD_ [normalization integral (2.19)].
The results obtained with the McMillan correlation:
f(r) = exp [- cf12 (b/r)^{5}] (2.22)
with b=1 sigma, at the equilibrium density, are shown in Table I and Fig. 3. The importance of the elementary diagram contributions is obvious. The close agreement between the HNC/S and Monte Carlo results of Refs. 13 and 14 suggests that the scaling approximations for the elementary diagrams are quite accurate.
The values of the three scaling factors for the McMillan Jastrow wave function at equilibrium density are found to be
s_{dd}= 2 , s_{wd}= 1 , s_{ww}= 1 , (2.23)
Puoskari and Kallio^{12}use both the two component mixture and Fantoni's formalism used here to calculate the rho ( r_{{1}1 prime} )_. At any level of approximation the mixture formalism and Fantoni's rho ( r_{{1}1 prime} )_ are proportional to each other. The only difference is that in mixture formalism the n_{0_}is calculated by using the normalization condition (2.20) in Eq. (2.5), whereas Fantoni calculates it independently by Eq. (2.9). PK also use scaling constants s_{dd_,}s_{dw_,}and s_{ww_}(their kappa_{{alpha}beta}_ equal 1+ s_{{alpha}beta}_ in our notation), and determine them from T_{roman}JF = T_{roman}PB_, T_{roman}MD = T_{roman}JF_, and T_{roman}MD_(mixture) =T_{roman}MD_ where T_{roman}MD_(mixture) is the kinetic energy obtained with n(k) from mixture formalism. This procedure is identical to ours because n(k)(mixture) is proportional to n(k), and so T_{roman}MD_(mixture)=T_{roman}MD_ is identical to the normalization condition. Thus we do not find that the mixture formalism offers any simplification. PK neglect the contribution of one body elementary diagrams E_{d_}and E_{w_;}we include them, but find that they are small.
A significant improvement in the variational energy of liquid helium is obtained by including three body correlations in the wave function.^{7,15}The wave function (the J+T denotes Jastrow plus triplet) is taken as
PSI_{{roman}J+T} = ^prod from {i < j} ^f( r_{ij}) ^prod from {i < j < k} ^ f_{3}( r vec_{ij}, r vec_{ik}) , L (3.1)
f_{3}( r vec_{ij}, r vec_{ik}) = exp ( - cf12 ^sum from cyc sum from l=0,2 ^ xi_{l}( r_{ij}) xi_{l}( r_{ik}) P_{l}( r vec_{ij}cdot r vec_{ik}) right ) . L
The l=1 term of f_{3_}gives the dominant contribution, the l=0 term gives a small contribution, and the l=2 term has negligible effect. ^{7}
HNC equations for the distribution functions of the J+T wave function have been discussed in Ref. 7. The HNC equations for the density matrix are obtained in an analogous way by replacing the E_{xy}, xy =dd_, wd, and ww as follows: (3.3) E_{xy}= C_{xy}+ E_{xy^{g}+}E_{xy^{t} . }(3.2)
Here C_{xy_}are three body elements given by
C_{dd}( r_{ij}) = rho ^int mark ^ [ f_{3^{2}(}r vec_{ia}, r vec_{ja}) -1]
lineup times g_{dd}( r_{ia}) g_{dd}( r_{ja}) d^{3}r_{a} , (3.4)
C_{wd}( r_{1j}) = rho ^int mark ^ [f_{3}( r vec_{1a}, r vec_{ja}) -1]
lineup times g_{wd}( r_{1a}) g_{dd}( r_{ja}) d^{3}r_{a} , (3.5)
C_{ww}( r_{{1}1 prime} ) = 0 . (3.6)
E_{xy^{g_}is}the sum of elementary diagrams having only g_{xy}-1_ bonds, and E_{xy^{t_}is}the sum of elementary diagrams having one or more three body correlations. The E_{dd,4^{t_}diagrams}are given in Fig. 1 of Ref. 7, and E_{wd,4^{t_,}E}sub ww,4^{t,}and E_{w,4^{t_}diagrams}are given in Figs. 4, 5, and 6, respectively. In these diagrams a wiggly line triangle 1jk with 1 as an external point represents
[ f_{3}( r vec_{1j}, r vec_{1k})-1] g_{wd}( r_{1j}) g_{wd}( r_{1k}) g_{dd}( r_{jk}) ,
whereas a plain triangle ijk represents
[ f_{3^{2}(}r vec_{ij}, r vec_{ik})-1] g_{dd}( r_{ij}) g_{dd}( r_{ik}) g_{dd}( r_{jk}) .
As in Ref. 7, a cross on a side ij of the triangle indicates that the factor g_{wd_}or g_{dd_}is to be omitted.
The equation (2.10) for R_{x},x=w,d_ becomes
R_{x}= mark ^ rho ^int ^d^{3}r ( g_{xd}-1 - N_{xd}- C_{xd}- E_{xd^{g}-}E_{xd^{t})
}
lineup - cf12 rho ^int mark ^ d^{3}r ( g_{xd}- 1) ( N_{xd}+ 2 E_{xd^{g}+}2E_{xd^{t1})
}
lineup - rho ^int^ d^{3}rg_{xd}( E_{xd^{t2}+}cf12 C_{xd}) + E_{x^{g}+}E_{x^{t} ,
}
where E_{xy^{t2}(r)_}[ E_{xy^{t1}(r)]_}are sums of E_{xy}(r)_ diagrams that contain [do not contain] a three body correlation connecting the two external points. For example, diagrams (4.1), (4.3), and (4.4) of Fig. 4 contribute to E_{wd,4^{t1_}while}the rest contribute to E_{wd,4^{t2_.}
}The bonds g_{dd}(r)-1_ and g_{wd}(r)-1_ do not change much when three body correlations are added to the wave function. Thus following Ref. 7 we assume that the scaling factors s_{xy_}for the E_{xy^{g_}do}not change by switching on the three body correlations. The f_{3_}is set to one and the s_{xy_}are calculated as discussed in the preceding section.
The contribution of E_{xy^{t_}diagrams}is generally smaller than that of the E_{xy^{g_}diagrams.}The E_{dd^{t}(r)_}changes sign^{7}at r app 0.5 sigma, and it is unlikely that a scaling approximation for the total E_{dd^{t_}is}valid. Hence in Ref. 7, and in the present work all five or more point E_{xy^{t_}diagrams}are neglected. We also neglect E_{ww,4^{t_}for}reasons discussed later. Thus in practice the f_{3_}is switched on after completing the Jastrow calculation, and the HNC equations for J+T are solved with E_{xy^{g_}and}E_{x^{g_}given}by the scaling approximations and E_{dd^{t}=}E_{dd,4^{t}, E}sub wd^{t}= E_{wd,4^{t} , }(3.7)
E_{w^{t}=}E_{w,4^{t}, E}sub d^{t}= E_{d,4^{t};} E_{ww^{t}=}0 . N (3.8)
At equilibrium density the calculation described above, carried out with the wave function of Ref. 7, gives T_{roman}MD = 14.8_ K, in close agreement with T_{roman}JF =14.72_ K, and normalization rho ( r_{{1}1 prime} =0) = 0 rho_. It thus satisfies the known identities to a reasonable accuracy. However, if we take
E_{ww^{t}=}E_{ww,4^{t} , }(3.9)
the T_{roman}MD = 18.8_ K and rho ( r_{{1}1 prime} =0)=1.065 rho_. The kinetic energy identity is particularly violated, and so it appears that it is better to neglect E_{ww^{t_}than}to approximate it with E_{ww,4_.}Crude numerical studies of E_{ww,5^{t_}(r=0)}show that at least at r=0, E_{ww,5^{t_}is}of the order of magnitude of E_{ww,4_}but opposite in sign.
In this section we present the results obtained with realistic J+T wave functions of Ref. 7 at three densities. In addition to the triplet correlations these wave functions contain an optimized pair correlation having the asymptotic behavior
f( r -> inf ) = 1 - mc over {2 pi^{2}hbar rho} 1 over r^{2} , (4.1)
a consequence of the long wavelength phonons.^{16}Here c is the velocity of zero sound. It may be easily verified that this asymptotic behavior of f implies the following asymptotic conditions for the nodal functions:
N_{dd}( r -> inf ) = mc over {pi^{2}hbar rho} 1 over r^{2} , (4.2)
N_{wd}( r -> inf ) = N_{dd}/2 , (4.3)
N_{ww}( r -> inf ) = N_{dd}/4 , (4.4)
and the n(k),
n(k -> 0) = n_{0}mc over {2 hbar} 1 over k . (4.5)
The results of our calculations are given in Tables II and III, and Figs. 7 ndash9. The results of the full calculation are labeled J+T, while those labeled J are obtained on switching off the triplet correlation. The GFMC results are taken from a calculation by Whitlock and Panoff.^{17}In all these results n(k) is normalized such that
int^n(k) d^{3}k=1 . (4.6)
In general we find that the triplet correlation by itself has little effect on the n(k). The n(k) is seen to decrease exponentially for k > 3 A ang^sup -1 in Fig. 7.
The kn(k) obtained from the neutron scattering data^{3,18}is compared with theoretical results in Fig. 8. Both the experimental and the GFMC n(k) do not have the correct k -> 0 asymptotic behavior. The difference between GFMC and J+T results has to be attributed to (i) the approximations in the use of J+T wave function, and those in the HNC/S calculation; and (ii) the finite box size in the GFMC simulation. The latter effect is particularly manifested at small k. The difference between theory and experiment may be mostly due to the inadequacy of the Aziz potential, or the impulse approximation used in relating the n(k) to neutron scattering cross sections at large momentum transfer. The use of impulse approximation for analysis of scattering from hard core liquids has been recently criticized. ^{19}There certainly is more than qualitative agreement between theory and experiment, marred by significant differences at k=0.5 and 2.3 A ang^sup -1. The density dependence of the J+T n(k) is given in Table II. The n(k) becomes broader as the density is increased.
The condensate fraction and the kinetic energies are given along with the scaling constants, in Table III. At rho = 0.438sigma^{-3}the T_{roman}MD_ is app 10% larger than the T_{roman}JF_ indicating increased importance of the neglected E_{xy^{t_}diagrams.}
The theoretical and experimental condensate fractions are compared in Fig. 9. At equilibrium density and GFMC and J+T values of n_{0_}are identical, but they are app 20% below the values deduced from neutron scattering experiments. The density dependence of the J+T n_{0_}is in crude agreement with that of Ref. 2. On the other hand the experiments of Wirth et al.^{20}have shown no density dependence of n_{0_,}while Mook^{21}finds a much stronger decrease in n_{0_}with rho.
The authors wish to thank Dr. R. O. Hilleke, Dr. R. M. Panoff, Dr. D. L. Price, Dr. V. F. Sears, Dr. R. O. Simmons, Dr. P. E. Sokol, and Dr. P. A. Whitlock for communicating their results. This work was supported by the U.S. Department of Energy, Division of Materials Sciences, under Contract No. DE AC02 76ER01198.
Previous article Next article Up to abstract
Home Retrieve Search Browse Feedback Help