---

Practical solution to the Monte Carlo sign problem:

arXiv:nucl-th/9310026v2 22 Dec 1993Practical solution to the Monte Carlo sign problem:Realistic calculations of 54FeY. Alhassid†, D. J.

Dean, S. E. Koonin, G. Lang, and W. E. OrmandW. K. Kellogg Radiation Laboratory, 106-38, California Institute of TechnologyPasadena, California 91125 USA(August 9, 2018)AbstractWe present a practical solution to the “sign problem” in the auxiliary fieldMonte Carlo approach to the nuclear shell model.

The method is based onextrapolation from a continuous family of problem-free Hamiltonians.Todemonstrate the resultant ability to treat large shell-model problems, wepresent results for 54Fe in the full fp-shell basis using the Brown-Richterinteraction. We find the Gamow-Teller β+ strength to be quenched by 58%relative to the single-particle estimate, in better agreement with experimentthan previous estimates based on truncated bases.PACS Nos.

21.60.Cs, 21.60.Ka, 27.40.+z, 25.40.KvTypeset using REVTEX1

Recent publications [1,2] have described quantum Monte Carlo methods for exact so-lution of the nuclear shell model. The methods are based on the Hubbard-Stratonovich(HS) representation [3] of the imaginary-time many-body propagator in terms of one-bodypropagators of non-interacting nucleons moving in a fluctuating field.

Thermal averages canbe calculated, as can ground state properties; errors arise only from discretization and sta-tistical sampling, both of which can be controlled. As these computations scale much moregently with the number of single particle orbits (Ns) and/or the number of valence nucleons(Nv) than do direct diagonalization techniques, they hold great promise for treating verylarge model spaces.Unfortunately, the applicability of shell model Monte Carlo calculations has heretoforebeen limited by the “sign problem” generic to all fermionic Monte Carlo techniques [1,2,4,5].The sign of the integrand may vary from sample to sample and the net integral resultsfrom a delicate cancellation that is difficult to reproduce with a finite number of samples.The problem is well-documented (and as yet unsolved) in simulations of correlated electronsystems [4].

Except for an important, yet schematic, class of nuclear interactions [2], wehave found that all realistic nuclear shell model Hamiltonians suffer from a sign problem.In this Letter, we report a practical solution to the sign problem and present the firstrealistic calculation of a mid-fp-shell nucleus, 54Fe [6]. Our method is based on an extrap-olation of observables calculated for a “nearby” family of Hamiltonians whose integrandshave a positive sign.

Success depends crucially upon the degree of extrapolation required.We have found that, for all of the many realistic interactions tested in the sd- and fp-shells,the extrapolation required is modest, amounting to a factor-of-two variation in the isovectormonopole pairing strength.A general time-reversal invariant Hamiltonian with two-body interactions can be broughtto the formH =Xαǫ∗α ¯Oα + ǫαOα+ 12XαVαnOα, ¯Oαo,(1)where the Oα are a convenient set of one-body operators and ¯O denotes the time-reverse of2

O. For real Vα, H in Eq.

(1) is a manifestly time-reversal invariant. The auxiliary field MonteCarlo approach utilizes the HS representation of the imaginary-time many-body propagatorU = exp(−βH) as a path integral over one-body propagators in fluctuating auxiliary fields.Upon introducing Nt time slices of duration ∆β = β/Nt and complex c-number auxiliaryfields σαn (n = 1, .

. .

, Nt), we can write the canonical expectation value of an observable Oas⟨O⟩≡Tr (Oe−βH)Tr (e−βH) ≈R D[σ]W(σ)Φ(σ)⟨O⟩σR D[σ]W(σ)Φ(σ). (2)Here,theapproximationbecomesexactasNt →∞and the metric is D[σ] = Πα,n [dσαndσ∗αn∆β|Vα|/2π].

The non-negative weightis W(σ) = ζ(σ) exp(−P |Vα||σαn|2∆β), where ζ(σ) ≡Tr Uσ is the canonical partition func-tion of the one-body evolution operator Uσ ≡UNt . .

. U1, where Un = exp(−∆βhn), and theone-body Hamiltonian for the nth time slice is hn =Pα(ǫ∗α+sαVασαn) ¯Oα+(ǫα+sαVασ∗αn)Oα,with sα = ±1 for Vα < 0 and sα = ±i for Vα > 0.

The “sign” is Φ(σ) ≡ζ(σ)/|ζ(σ)| and⟨O⟩σ ≡Tr (OUσ)/ζ(σ). Both ζ(σ) and ⟨O⟩σ can be evaluated in terms of the Ns×Ns matrixUσ that represents the evolution operator Uσ in the single-particle space.The sign problem arises because the one-body partition function ζ(σ) is not necessarilypositive, so that the Monte Carlo uncertainty in the denominator of Eq.

(2) (the W-weightedaverage sign, ⟨Φ⟩) can become comparable to or larger than ⟨Φ⟩itself. In most cases ⟨Φ⟩decreases exponentially with β or with the number of time slices [5].An important class of interactions free from the sign problem (i.e., Φ(σ) ≡1) was foundin Ref.

[2]. This occurs when Vα < 0 for all α in Eq.

(1). In that case, sα = 1 for all α,so that both hn and Uσ are time-reversal invariant.

The eigenvectors of Uσ then occur astime-reversed pairs with complex conjugate eigenvalues λi, λ∗i (i = 1, . .

. , Ns/2), the grandcanonical partition function ζ(σ) = Πi|1+λi|2 is positive definite, and the canonical partitionfunction for even Nv is also positive definite.Based on the above observation, it is possible to decompose H into its “good” and “bad”parts, H = HG + HB, with3

HG =Xα(ǫ∗α ¯Oα + ǫαOα) + 12XVα<0VαnOα, ¯OαoHB = 12XVα>0VαnOα, ¯Oαo. (3)The “good” Hamiltonian HG includes, in addition to the one-body terms, all the two-bodyinteractions with Vα < 0, while the “bad” Hamiltonian HB contains all interactions withVα > 0.

By construction, calculations with HG alone have Φ(σ) ≡1 and are thus free of thesign problem.We define a family of Hamiltonians Hg that depend on a continuous real parameter g asHg = HG + gHB, so that Hg=1 = H. If the Vα that are large in magnitude are “good”, weexpect that Hg=0 = HG is a reasonable starting point for the calculation of an observable⟨O⟩. One might then hope to calculate ⟨O⟩g = Tr (Oe−βHg)/Tr (e−βHg) for small g > 0 andthen to extrapolate to g = 1, but typically ⟨Φ⟩collapses even for small positive g. However,it is evident from our construction that Hg is characterized by Φ(σ) ≡1 for any g ≤0,since all the “bad” Vα(> 0) are replaced by “good” gVα < 0.

We can therefore calculate⟨O⟩g for any g ≤0 by a Monte Carlo sampling that is free of the sign problem. If ⟨O⟩g is asmooth function of g, it should then be possible to extrapolate to g = 1 (i.e., to the originalHamiltonian) from g ≤0.

We emphasize that g = 0 is not expected to be a singular pointof ⟨O⟩g; it is special only in the Monte Carlo evaluation.In the nuclear shell model, the two-body interaction can be written in a density decom-position as [2]12XabcdXKTπEπKT(ac, bd)XM(−)MρKMT(ac)ρK−MT(bd) .Here ρKMT = ρpKM + (−)TρnKM (T = 0, 1), ρ(p,n)KM (ac) = (a†a × ˜ac)KM is the one-body densityoperator for the pair of proton or neutron orbits (a, c) coupled to angular momentum Kand its z-projection M, and π = (−)la+lc = (−)lb+ld is the parity. The matrices EπKT areconstructed from the two-body matrix elements V πJT(ab, cd) of good angular momentum J,isospin T, and parity π through a Pandya transformation.

For interactions that are time-reversal invariant and conserve parity, the EπKT(i, j) are real symmetric matrices that can be4

diagonalized by a real orthogonal transformation. The eigenvectors ρKM(α) play the role ofOα in Eq.

(1), and the eigenvalues λKπ(α) are proportional to Vα. In the Condon-Shortley[7] convention ¯ρKM = π(−)K+MρK−M so that the “good” eigenvalues satisfy sign [λKπ(α)] =π(−)K+1 [8].

To minimize the number of auxiliary fields required, we use the freedom to addan arbitrary symmetric interaction to H [2] and choose V SJT=1 = V AJT=0 so that EKT=1 ≡0.EKT=0 is then uniquely determined by the anti-symmetric part of the interaction throughthe combinationV AJT=0 + V AJT=1.To demonstrate the viability and utility of the method, we have applied it to the mid-fpshell nucleus 54Fe using the realistic Brown-Richter interaction [9]. The number of m-schemeSlater determinants describing the 6 valence protons and 8 valence neutrons moving amongthe Ns = 20 single-particle states of the 0f7/2,5/2 and 0p3/2,1/2 orbitals is 206 208≈5 · 109.For comparison, the largest model space treated by standard diagonalization techniques iscurrently 48Ti [10] where the m-scheme dimension is ≈7 · 106.Figure 1 (upper) shows the eigenvalues VKπα = π(−)KλKπ(α) of the Brown-Richterinteraction; only about half of the eigenvalues are negative.

However, those with the largestmagnitude are all “good”. It is possible to use an inverse Pandya transformation to calculatethe usual two-body matrix elements V πJT(ab, cd) for the “good” and “bad” interactions,allowing the matrix elements of HG to be compared in Fig.

1 (lower) with those of the fullinteraction. The greatest deviation is for J = 0, T = 1 (the monopole pairing interaction),where HG is about twice as attractive as the physical H. In all other channels, HG and Hare quite similar.We have performed Monte Carlo calculations for β = 2 MeV−1 using Nt = 32 (so that∆β = 0.0625 MeV−1).

For g = −1, −0.8, −0.6, −0.4, −0.2, and 0, we took approximately3300 uncorrelated samples.The computations were performed on the Intel TouchstoneDELTA 512-node parallel computer, where each node is an Intel i860 processor. Each nodeproduced and analyzed a sample in about 4 minutes, so that each value of g took about25 minutes in total.

Selected calculations for larger values of β or Nt show that we haveconverged to the true ground-state properties.5

The results for various observables are shown in Fig. 2.

The extrapolations to the physicalHamiltonian (g = 1) are done by least-squares polynomials. For each observable except⟨H⟩, the degree of the polynomial is chosen to be the lowest for which χ2 per degree offreedom is less than 1; linear or quadratic extrapolations are almost always sufficient.

For⟨H⟩, the variational principle implies the additional constraint of vanishing derivative atg = 1, in which case a quadratic or cubic polynomial is used. We have also calculatedresponse functions R(τ) = ⟨O†(τ)O(0)⟩by polynomial extrapolation of our calculationsof ln[Rg(τ)/Rg(0)] for g ≤0.

Fitting ln[R1(τ)/R1(0)] to a polynomial in τ allows us todetermine moments of the normalized strength function fO(E), such as ¯E ≡R EfO(E)dE.Our overall method was checked in detail [11] by comparison with direct diagonalizationin the sd-shell using the Brown-Wildenthal interaction [12] and in the lower fp-shell (44Ti)using the Brown-Richter interaction [9].Table 1 summarizes the extrapolated results for various observables.Note that thestatistical uncertainty in these values is proportional to the uncertainties in the MonteCarlo results for g ≤0, and so can be reduced by increasing the number of samples. Thecalculated first moment of the isoscalar quadrupole strength function, (1.25 ± 0.16) MeV,should be compared with the empirical excitation energy of the first 2+ state, 1.408 MeV.

Ourestimate for the B(E2) for the decay of this state assuming free nucleon charges (and thatthis transition has all of the strength) is (96 ± 1) e2 fm4, while effective charges (ep, en) =(1.1, 0.1)e would be required to reproduce the experimental value of 126 e2fm4.Thesecharges are significantly smaller than the (1.35, 0.35)e used in truncated calculations [13]or the (1.33, 0.64)e used in the lower fp shell [9]. The total mass quadrupole strength,⟨Q2⟩= (1482 ± 84) fm4, is significantly larger than the simple single particle estimate of380 fm4.

The total M1 strength ⟨(M1)2⟩= (14.1±0.4) µ2N is quenched relative to the singleparticle estimate of 42.55 µ2N. It is also interesting to note that the occupation numbers ofthe single particle orbits are smeared across the Fermi surface.Of particular physical interest are the Gamow-Teller operators.

Our calculations exactlysatisfy the sum rule ⟨(GT−)2⟩−⟨(GT+)2⟩= 3(N −Z) = 6. The single particle estimate for6

⟨(GT+)2⟩corresponding to the f7/2 proton →f5/2 neutron transition is 10.28 [14], so the shellmodel Monte Carlo value of 4.32±0.24 is quenched by 58%. This value is comparable to theexperimental result of 3.1±0.6 [15], but significantly smaller than previous estimates of 6.40or 6.70 based on truncated bases [13].

The additional quenching on the full space correlateswith the enhanced B(E2, 2+1 , →0+1 ), (i.e., smaller effective charge), as was surmised in [13].Direct comparison with experimental Gamow-Teller strength functions requires that weknow the energy of the daughter ground state relative to 54Fe. Since the 54Co ground stateis the isobaric analog state (IAS) of the 54Fe ground state, we find a mean (p, n) excitationenergy of ¯Ex = (6.13 ± 0.17) MeV.

This is in agreement with the systematics of Nakayamaet al. [16], which give EGT−−EIAS = 5.81 MeV, but is somewhat low relative to theexperimental value of 8.2 MeV [15].

When our calculations of the mean (n, p) excitationenergy are corrected for the Coulomb energy (including exchange) and the nucleon massdifference, we find ¯Ex = (1.24 ± 0.2) MeV, to be compared with the experimental centroidof 3 MeV [15]. A more consistent theoretical value of ¯Ex can be obtained by calculating themass differences of the A = 54 isobars within the shell model Monte Carlo [11].The method presented in this Letter is a practical solution to the sign problem for realisticshell model interactions.

A full-basis calculation of 54Fe with the Brown-Richter interactionshows the feasibility of the method, with significant quenching of the Gamow-Teller β+strength. Systematic studies of the temperature, nuclide, and interaction dependence ofthese calculations will be reported elsewhere.

Our techniques also enable the determinationof an optimal effective interaction and effective operators in a greatly enlarged model space.This work was supported in part by the NSF, Grants No. PHY90-13248 and PHY91-15574, by the DOE, Grant No.

DE-FG02-91ER40608, and by a Caltech DuBridge postdoc-toral fellowship to WEO. We thank P. Vogel and B.

A. Brown for helpful discussions, andthe Concurrent Supercomputing Consortium for a grant of DELTA time.7

REFERENCES† Permanent Address: Center for Theoretical Physics, Sloane Physics Laboratory, YaleUniversity, New Haven, Connecticut 06511. [1] C. W. Johnson, S. E. Koonin, G. H. Lang, and W. E. Ormand, Phys.

Rev. Lett.

69, 3157(1992). [2] G. H. Lang, C. W. Johnson, S. E. Koonin, and W. E. Ormand, Phys.

Rev. C 48, 1518(1993); W. E. Ormand, D. J.

Dean, C. W. Johnson, G. H. Lang, and S. E. Koonin, Phys.Rev. C (in press).

[3] J. Hubbard, Phys. Rev.

Lett. 3, 77 (1959); R. L. Stratonovich, Dokl.

Akad. Nauk.

S.S.S.R.115, 1097 (1957). [4] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis, and R. T. Scalettar,Phys.

Rev. B40, 506 (1989).

[5] For a review see W. Von der Linden, Phys. Rep. 220, 53 (1992); E. Y. Loh, Jr. and J.E.

Gubernatis, in Electronic Phase Transitions, edited by W. Hanke and Y. V. Kopaev(North Holland, Amsterdam, 1992). [6] The results presented here supersede and correct two conference reports: S. E. Koonin, inProceedings of the International Conference on Perspectives in Nuclear Structure, Copen-hagen, 1993, edited by J. J. Gaardhøje et al.

(Elsevier, Amsterdam, 1993); D. J. Dean, inProceedings of the International Conference on Computational Physics II, Beijing, 1993(International Press, Hong Kong, 1993).

[7] E. V. Condon and G. H. Shortley, The Theory of Atomic Spectra, Ch. III (CambridgeUniversity Press, 1935).

[8] If all single-particle orbits in the model space have the same parity, π = +1 and thesign rule of Ref. [2] is reproduced.

The sign rule of Ref. [2] holds in the Biedenharn-8

Rose convention, but in the calculation of the two-body matrix elements the sphericalharmonics Ylm should be multiplied by il. [9] W. A. Richter, M. G. Vandermerwe, R. E. Julies, and B.

A. Brown, Nucl.

Phys. A523,325 (1991).

[10] E. Caurier, A. Poves, and A. P. Zuker, Phys. Lett.

B256, 301 (1991). [11] Y. Alhassid, D. J.

Dean, S. E. Koonin, G. H. Lang, and W. E. Ormand, to be submittedto Phys. Rev.

C.[12] B. A.

Brown and B. H. Wildenthal, Ann. Rev.

Nucl. Part.

Sci. 38, 29 (1988).

[13] N. Auerbach, G. F. Bertsch, B. A.

Brown, and L. Zhao, Nucl. Phys.

A556, 190 (1993). [14] M. B. Aufderheide, G. E. Brown, T. T. S. Kuo, D. B. Stout, and P. Vogel, Ap.

J. 362,241 (1990).

[15] M. C. Vetterli et al., Phys. Rev.

C40, 559 (1989). [16] K. Nakayama, A. P. Galeao, and F. Krmpotic, Phys.

Lett. B114, 217 (1982).9

FIGURESFIG. 1.

Upper: The eigenvalues Vα of the Brown-Richter interaction in the fp-shell. Eigen-values for each particle-hole angular momentum K are plotted in increasing order.

Bottom: Thetwo-body matrix elements VJT=1(ab, cd) of the Brown-Richter interaction (solid circles) and its“good” part (open circles), for J ≤4. The ordering for each J is arbitary.

Plots of the VJT=0 andthe remaining T = 1 matrix elements (not shown) are similar to those shown for J ≥1.FIG. 2.

The results of the Monte Carlo calculations for 54Fe at β = 2 MeV−1 for severalobservables as a function of g ≤0. Q = Qp + Qn is the isoscalar quadrupole, Qv = Qp −Qnis the isovector quadrupole, GT+ is the Gamow-Teller operator changing a proton to a neutron,and M1 is the magnetic moment operator using the free-nucleon g-factors.

The lines are polyno-mial extrapolations; the extrapolated values and corresponding uncertainties are shown at g = 1.The extrapolation is linear for ⟨M12⟩, but quadratic for ⟨Q2⟩, ⟨Q2v⟩, and ⟨GT2+⟩. For ⟨H⟩, theextrapolation is cubic with the constraint of vanishing derivative at g = 1.10

TABLESTABLE I. Monte Carlo results for 54Fe.⟨H⟩= −55.5 ± 0.5 MeVTotal strength¯E (MeV)Isoscalar quadrupole⟨Q2⟩= 1482 ± 84 fm41.25 ± 0.16Isovector quadrupole⟨Q2v⟩= 381.3 ± 33.8 fm412.7 ± 0.2Gamow-Teller (p, n)⟨(GT−)2⟩= 10.32 ± 0.246.13 ± 0.17Gamow-Teller (n, p)⟨(GT+)2⟩= 4.32 ± 0.249.7 ± 0.2M1⟨(M1)2⟩= 14.1 ± 0.4 µ2N8.6 ± 0.7Occupation NumbersProtonsNeutrons⟨a†a⟩f7/2 = 4.92 ± 0.03⟨a†a⟩f7/2 = 6.35 ± 0.03⟨a†a⟩p3/2 = 0.56 ± 0.02⟨a†a⟩p3/2 = 0.86 ± 0.02⟨a†a⟩p1/2 = 0.11 ± 0.01⟨a†a⟩p1/2 = 0.17 ± 0.01⟨a†a⟩f5/2 = 0.41 ± 0.01⟨a†a⟩f5/2 = 0.61 ± 0.0111

---

출처: arXiv:9310.026 • 원문 보기