{'fraction_non_alphanumeric': 0.040204032452405815, 'fraction_numerical': 0.011312822984443196, 'mean_word_length': 4.388860832551764, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Monte Carlo evaluation of path integrals for the nuclear shell model

We present in detail a formulation of the shell model as a path integral and Monte Carlo techniques for its evaluation. The formulation, which linearizes the two-body interaction by an auxiliary field, is quite general, both in the form of the effective `one-body' Hamiltonian and in the choice of ensemble. In particular, we derive formulas for the use of general (beyond monopole) pairing operators, as well as a novel extraction of the canonical (fixed-particle number) ensemble via an activity expansion. We discuss the advantages and disadvantages of the various formulations and ensembles and give several illustrative examples. We also discuss and illustrate calculation of the imaginary-time response function and the extraction, by maximum entropy methods, of the corresponding strength function. Finally, we discuss the"sign-problem"generic to fermion Monte Carlo calculations, and prove that a wide class of interactions are free of this limitation.


I. MOTIVATION AND ORGANIZATION
Exact diagonalizations of shell model Hamiltonians in the 0s-1d shell demonstrated [1] that the shell model can yield an accurate and consistent description for a wide range of nuclear properties in different nuclei, if the many-body basis is sufficiently large. However, the combinatorial scaling of the many-body space with the size of the single particle basis or the number of valence nucleons restricts such exact diagonalizations to small nuclei or nuclei with few valence particles [2].
In facing the general challenge of developing non-perturbative methods to describe strongly interacting manybody systems, various quantum Monte Carlo schemes have been proposed as an alternative to direct diagonalization [3]. Among them, the auxiliary-field Monte Carlo method is suitable for addressing interacting fermions [4]. This method is a Monte Carlo evaluation of the path-integral obtained by the Hubbard-Stratonovich transformation [5] of the imaginary-time evolution operator. The many-body wavefunction is represented by a set of determinantal wavefunctions evolving in fluctuating auxiliary fields. The method thus enforces the Pauli principle exactly, and the storage and computation time scale gently with the single-particle basis or the number of particles. Auxiliary field methods have been applied to condensed matter systems such as the Hubbard Model [3,6], yielding important information about electron correlations and magnetic properties.
In this paper, we discuss the application of auxiliary-field Monte Carlo techniques to the nuclear shell model. This involves a two-body interaction more general than the simple on-site repulsion of a Hubbard model. Our goal is to develop new methods for extending the applicability of the shell model, as well as to investigate the powers and limitations of auxiliary-field Monte Carlo methods for general fermion systems.
We have previously published a letter [7] with selected results for static observables in sd-and fp-shell nuclei. This paper serves to give the details of the implementation, introduce a method for calculating dynamical correlations and strength functions, demonstrate the method with simple (sd-shell) nuclei, and discuss several important issues that arise in the implementation. We also explore the limitations imposed by the negative contributions in the path-integral, referred to as the "sign problem" in the literature.
Our presentation is organized as follows. We begin in Section II by using the Hubbard-Stratonovich (HS) transformation to write the imaginary time evolution operator exp(−βĤ) as a path integral. This requires that the HamiltonianĤ be cast as a quadratic form using an appropriate set of operators. We discuss in Section III two ways in which this can be accomplished, using either particle density or pairing operators. The imaginary time evolution operator can be used to extract information about the system at finite temperature or in its ground state (β → ∞).
The formulae for obtaining static observables are presented in Section IV where methods for handling the canonical ensemble are also introduced. We then discuss in Section V the extraction of the strength functions for operators from the imaginary time response function. In Section VI, we briefly describe the computational algorithms for implementing our methods and present selected results of calculations for sd-shell nuclei. In the final section, we address the sign problem and also discuss a class of non-trivial interactions that give rise to a positive definite path-integral for some nuclei.

II. THE IMAGINARY TIME EVOLUTION OPERATOR
Given some many-body HamiltonianĤ, we seek a tractable expression for the imaginary time evolution operator: U = exp(−βĤ) . (2.1) Here, β has units of inverse energy and β −1 can be interpreted as an imaginary time. (Here and throughout, we takē h = 1 and measure all energies in units of MeV.) It is also clear thatÛ can be interpreted as the partition operator for temperature β −1 . We will refer toÛ as the evolution operator hereafter. The operatorĤ is usually a generalized Hamiltonian and might contain terms beyond the true Hamiltonian, such as −µN in the grand-canonical ensemble or −ωĴ z if we are 'cranking' the system.
There are two formalisms for extracting information from the evolution operator: the "thermal" formalism (on which we will concentrate) and the "zero-temperature" formalism (to which the thermal formalism reduces in the limit β → ∞). In the thermal formalism, we begin with the partition function Z =Tr exp(−βĤ) , (2.2) and then construct the thermal observable of an operatorÔ: Here, the traceTr is over many-body states of fixed (canonical) or all (grand-canonical) particle number. In the zero-temperature formalism we begin with a trial wavefunction ψ 0 and use the evolution operator to project out the ground state, assuming that ψ 0 is not orthogonal to the ground state. The expectation value ofÔ is then given by In this section, we describe how to writeÛ in a form that allows Eqs. (2.3) or (2.4) to be evaluated.

A. Path integral formulation of the evolution operator
We restrict ourselves to generalized Hamiltonians that contain at most two-body terms. The HamiltonianĤ can then be written as a quadratic form in some set of 'convenient' operatorsÔ α : where we've assumed that the quadratic term is diagonal in theÔ α . The meaning of 'convenient' will become clear shortly, but typically it refers to one-'body' operators, either one-particle ('density') or one-quasiparticle ('pairing').
The strength of the two-body interaction is characterized by the real numbers V α .
ForĤ in the quadratic form (2.5), one can write the evolution operatorÛ as a path integral. The exponential is first split into N t 'time' slices, β = N t ∆β, so that U = exp(−∆βĤ) Nt . (2.6) Then we perform the Hubbard-Stratonovich (HS) transformation on the two-body term for the n'th time slice to give [5], where the phase factor s α is ±1 if V α < 0 and is ±i if V α > 0. Each real variable σ αn is the auxiliary field associated withÔ α at time slice n.
(2. 13) In the limit of an infinite number of time slices Eq. (2.8) is exact. In practice one has a finite number of time slices and the approximation is valid only to order ∆β. The case of only one time slice is known as the Static Path Approximation (SPA); previous work on the SPA and its extensions can be found in Refs. [8] and [9].
Rewriting the evolution operator as a path integral can make the model space tractable. Consider the case where theÔ α are density operators. Then Eq. (2.1) is an exponential of two-body operators; it acts on a Slater-determinant to produce a sum of many Slater-determinants. In contrast, the path-integral formulation (2.8) contains only exponentials of one-body operators which, by Thouless' theorem [10], takes a Slater-determinant to another single Slater-determinant. Therefore, instead of having to keep track of a very large number of determinants (often many thousands for modern matrix-diagonalization shell model codes such as OXBASH [11]), we need deal only with one Slater-determinant at a time. Of course, the price to be paid is the evaluation of a high-dimensional integral. However, the number of auxiliary fields grows only quadratically with the size of the single particle basis while the corresponding number of Slater-determinants grows exponentially. Furthermore, the integral can be evaluated stochastically, making the problem ideal for parallel computation.

B. Monte Carlo evaluation of the path integral
Formulating the evolution operator as a path integral over auxiliary fields reduces the problem to quadrature. For a limited number of auxiliary fields, such as in the SPA with only a quadrupole-quadrupole interaction, the integral can be evaluated by direct numerical quadrature. However, for more general cases (typically hundreds of fields), the integral must be evaluated stochastically using Monte Carlo techniques.
Using the one-body evolution operator defined bŷ For the zero-temperature formalism while for the thermal formalism (canonical and grand-canonical), and To evaluate the path integral via Monte Carlo techniques, we must choose a normalizable positive-definite weight function W σ , and generate an ensemble of statistically independent fields {σ i } such that the probability density to find a field with values σ i is W σi . Defining the 'action' by the required observable is then simply (2.21) where N is the number of samples, and S i ≡ S σi , etc. Ideally W should approximate exp(−S) closely. However, exp(−S) is generally not positive and can even be complex. In some cases, Φ i may oscillate violently, giving rise to a very small denominator in Eq. (2.21) to be cancelled by a very small numerator. While this cancellation is exact analytically, it is only approximate in the Monte Carlo evaluation so that this 'sign problem' leads to large variances in the evaluation of the observable.
There are several possible schemes for both the choice of W and the sampling of the fields. We typically choose W = |exp(−S)| and generate the samples via random walk (Metropolis) methods.

III. DECOMPOSITIONS OF THE HAMILTONIAN
To realize the HS transformation, the two-body parts ofĤ must be cast as a quadratic form in one-body operatorŝ O α . As these latter can be either density operators or pair creation and annihilation operators (or both), there is considerable freedom in doing so. In the simplest example, let us consider an individual interaction term, where a † i , a i are anti-commuting fermion creation and annihilation operators. In the pairing decomposition, we write (using the upper and lower bracket to indicate the grouping) The commutator is a one-body operator that can be put directly in the one-body Hamiltonianĥ σ . The remaining two quadratic forms in pair-creation and -annihilation operators can be coupled to auxiliary fields in the HS transformation.
In the density decomposition, there are two ways to proceed: we can group (1,3) and (2,4) to get or group (1,4) and (2,3) to get Again the commutator terms are one-body operators, but now the quadratic forms are squares of density operators that conserve particle number. We refer to Eq. (3.3) as the 'direct' decomposition and Eq. (3.4) as the 'exchange' decomposition.
For any general two-body Hamiltonian, we can choose the pairing or density decompositions for different parts of the two-body interaction. Moreover, even within a pure density break-up decomposition, there is still freedom to choose between the direct and exchange formulations. Although the exact path integral result is independent of the scheme used, different schemes will lead to different results under certain approximations (e.g., mean field or SPA). The choice of decomposition will also affect the rate of convergence of our numerical result as N t → ∞, as well as the statistical precision of the Monte Carlo evaluation. Most significantly, it affects the fluctuation of Φ in Eq. (2.21) and thus determines the stability of the Monte Carlo calculation. (See Section VI below.) In the application of these methods to the nuclear shell model, it is particularly convenient to use quadratic forms of operators that respect rotational invariance, isospin symmetry, and the shell structure of the system. We introduce these in the following subsections for both the density and pairing decompositions.

A. Density decomposition
We begin by ignoring explicit isospin labels and by writing the rotationally invariant two-body Hamiltonian aŝ where the sum is taken over all proton and neutron single-particle orbits (denoted by a, b, c, d) and the pair creation and annihilation operators are given by The V J (ab, cd) are the angular-momentum coupled two-body matrix elements of a scalar potential V ( r 1 , r 2 ) defined as (independent of M) while the anti-symmetrized two-body matrix elements V A J (ab, cd) are given by To obtain the density decomposition ofĤ 2 , we perform a Pandya transformation to recouple (a, c) and (b, d) into density operators with definite multipolaritŷ whereã jama = (−1) ja+ma a jama . ThenĤ 2 can be rewritten aŝ 14) where the particle-hole matrix elements of the interaction are andĤ ′ 1 is a one-body operator given byĤ Note that adding symmetric matrix elements is equivalent to using the exchange decomposition for some parts of the interaction. The freedom in choosing the combinations of direct and exchange decomposition is then embodied in the arbitrary symmetric part of the matrix elements V N J . Introducing the shorthand notation i = (ac), j = (bd), we can write Eq. (3.14) aŝ Upon diagonalizing the matrix E K (i, j) to obtain eigenvalues λ Kα and associated eigenvectors v Kα , we can represent whereρ Finally, if we defineQ This completes the representation of the two-body interaction as a diagonal quadratic form in density operators. We then couple auxiliary fields σ KM (α) toQ KM and τ KM (α) toP KM in the HS transformation. (The latter are not to be confused with the "imaginary time" τ .) In the treatment thus far, protons and neutrons were not distinguished from each other. Although the original HamiltonianĤ 2 conserves proton and neutron numbers, we ultimately might deal with one-body operators ρ KM (a p , b n ) and ρ KM (a n , b p ) (n,p subscripts denoting neutron and proton) that individually do not do so. The one-body hamiltonianĥ σ appearing in the HS transformation then mixes neutrons and protons. The single-particle wavefunctions in a Slater determinant then contain both neutron and proton components and neutron and proton numbers are not conserved separately in each Monte Carlo sample; rather the conservation is enforced only statistically.
It is, of course, possible to recouple so that only density operators separately conserving neutron and proton numbers (ρ KM (a p , b p ) andρ KM (a n , b n )) are present. To do so, we write the two-body Hamiltonian in a manifestly isospin-invariant form, where, similar to the previous definition (3.6), the pair operator is Here ( 1 2 , t a ), etc. are the isospin indices with t a = − 1 2 for proton states and t a = 1 2 for neutron states, and (T T z ) are the coupled isospin quantum numbers. The two-body Hamiltonian can now be written solely in terms of density operators that conserve the proton and neutron numbers. Namely, Here, we defineρ KM,T asρ K M,T =ρ K M,p + (−1) Tρ K M,n , (3.28) and the E K,T are given by Note that these expressions allow less freedom in manipulating the decomposition since we have to couple proton with proton and neutron with neutron in forming the density operators. Also note that E K,T =0 (ac, bd)−E K,T =1 (ac, bd) is an invariant related only to the physical part of the interactions, (V A J,T =1 + V A J,T =0 ). We can choose all E K,T =1 to be zero in the above (by setting V S J,T =1 = V A J,T =0 ) leaving E K,T =0 completely determined by the physical matrix elements. In that case, we can halve the number of fields to be integrated. However, while introducing the isovector densities requires more fields, it also gives more freedom in choosing the unphysical matrix elements to optimize the calculation.
If we now diagonalize the E KT (i, j) as before and form the operatorŝ the two-body part of the Hamiltonian can finally be written aŝ In this decomposition, the one-body hamiltonianĥ σ of the HS transformation does not mix protons and neutrons. We can then represent the proton and neutron wavefunctions by separate determinants, and the number of neutrons and protons will be conserved rigorously during each Monte Carlo sample. For general interactions, even if we choose nonzero E K,T =1 matrix elements, the number of fields involved is half that for the neutron-proton mixing decomposition, and the matrix dimension is also halved. These two factors combine to speed up the computation significantly. In this sense, an isospin formalism is more favorable, although at the cost of limiting the degrees of freedom embodied in the symmetric matrix elements V S J .

B. Pairing decomposition
In nuclei where the pairing interaction is important, it is natural to cast at least part of the two-body interaction as a quadratic form in pair creation and annihilation operators. We demonstrate this for the case where the Hamiltonian is written in the isospin formalism. Upon diagonalizing V A JT (ab, cd) in Eq. (3.22), we can writê Separating A † A into commutator and anticommutator terms, we havê Clearly,Ĥ ′ 1 is a one-body operator that can be combined withĤ 1 . The remaining two-body term can be written as a sum of squares by defining As in the density decomposition, we can then couple the σ and τ fields to Q and P , respectively.
Note that in the pairing decomposition, the one-body Hamiltonian h(τ ) used in the path integral is a generalized one-body operator that includes density, pair-creation and pair-annihilation operators. The wavefunction is then propagated as a Hartree-Fock-Bogoliubov state, rather than as a simple Slater determinant.
In this decomposition, neutrons and protons are inevitably mixed together in the one-body Hamiltonianĥ σ (consider the Q,P terms for T = 0). In fact,ĥ σ also does not conserve the total number of nucleons; rather the conservation is only statistical after a large number of Monte Carlo samples.
For simplicity, we have described how to decompose the Hamiltonian solely in density operators or solely in pair operators. However, it is straightforward to mix the two decompositions with the choice depending on the type of interactions involved. Consider the 'Pairing plus Quadrupole' model, namelŷ where P † , P are the monopole pair creation and annihilation operator and Q is the quadrupole-moment operator, Naively, it would be most convenient to use a pairing decomposition for the pairing interaction and a density decomposition for the quadrupole interaction. This would require only 8 fields for each time slice in the HS transformation. The pure pairing or pure density decompositions would be much more complicated, as the former would require rewriting the quadrupole interactions in terms of pair operators and the latter would require rewriting the pairing interaction in terms of density operators. The numerical examples given in Section VI below will illustrate the behavior of such Hamiltonians under different decompositions.

IV. CALCULATION OF STATIC OBSERVABLES
In the previous sections, we expressed the evolution operator for a two-body Hamiltonian as a path integral over auxiliary fields in which the action involves only density and pair operators. In this and the following section, we show how to extract observables from this path integral. The merit of this method will be clear from the compact formulae involved, which require handling only relatively small matrices of dimension N s or 2N s , depending on the type of decomposition used. We derive formulae for three different approaches that use the evolution operator to obtain information about the system: the zero-temperature formalism, the grand-canonical ensemble, and the canonical ensemble. The zero-temperature and grand-canonical ensemble methods have been applied, with the density decomposition, to other physical systems such as the Hubbard Model [6]. However, we believe that the canonical ensemble treatment presented here is novel.
The zero-temperature approach can be used only to extract ground-state information. On the other hand, the grand-canonical ensemble allows finite-temperature calculations, but fluctuations in particle number can be very significant in a system with a small number of nucleons. Thus, the canonical ensemble is particularly important in nuclear systems, where the particle number is small and shell structure is prominent. The grand-canonical ensemble yields information averaged over neighboring nuclei, which can have very different properties.
The canonical ensemble is more difficult than the other two approaches in two respects. First, the canonical (fixednumber) trace of U σ is more difficult to compute than the wavefunction overlap of the zero-temperature formalism or the grand-canonical trace. Second, observables are more difficult to extract since Wick's theorem does not apply. We suggest two different methods to handle the canonical ensemble: the activity expansion and the integration over real coherent states. The coherent-state integration method can be applied to calculate canonical trace of any particle number involving only a negligible increase in computation time relative to the zero-temperature and grand-canonical approaches. Unfortunately, its utility is limited by the sign problem. On the other hand, the activity expansion, which is well-suited (and numerically stable) for calculating nuclei with a small number of valence particles or holes, is less susceptible to the sign problem. (Yet a third method for extracting the canonical trace, based on Fourier analysis, shows good "sign" statistics and is numerically stable for mid-shell calculations; it is therefore suitable when the activity expansion fails. Details of this method will be published elsewhere [12].) For each approach, we also derive the general formulae when pair operators (as well as density operators) are present in the single-particle Hamiltonian; i.e., when a pairing decomposition is used for some or all of the interaction. As far as we know, there is no known general formulae for the pairing decomposition formalism. Using the fermion coherent state formalism [13], we derive in Appendix A a set of formulae for the calculation of observables that are similar to the well-known formulae for a pure density decomposition. Thus, our methods can be extended to calculations using general one-body operators in the HS transformation by only doubling the dimension of the matrices involved.

A. Zero-temperature formalism
We begin with 'zero-temperature' observables. The trial wavefunction ψ 0 in Eq. (2.4) is, in principle, any state not orthogonal to the ground state. In practice, it is most conveniently a Slater-determinant. If we have N s m-scheme orbitals available and N v indistinguishable particles, ψ 0 is constructed from N v single-particle wavefunctions, each of which is a vector with N s components, and we write ψ 0 as a Slater-determinant of a N v × N s matrix, Ψ 0 .
In the pure density decomposition, consider the one-body Hamiltonianĥ = M ij a † i a j (sums over the indices i, j from 1 to N s are implicit). The evolved wavefunction, is then a Thouless transformation of the original determinant, the new state being represented by the matrix exp(−∆βM )Ψ 0 . We can therefore represent the product of evolution operators by N s × N s matrices Let us now consider the overlap function in Eq. (2.16). Let U σ (τ 2 , τ 1 ) be the matrix representing the evolution operatorÛ σ (τ 2 , τ 1 ). Choosing some value τ of the imaginary time at which to insert the operatorÔ, we introduce the right and left wavefunctions ψ R,L (τ ) defined by Then the required overlap in Eq. (2.16) is are the matrices representing ψ R and ψ L . Note that if there are two distinguishable species of nucleons -protons and neutrons -and we use a decomposition that conserves the numbers of neutrons and protons, then there is a separate determinant for each set of single-particle wavefunctions and the total overlap is the product of the proton and neutron overlaps.
With the basic overlap in hand, we now turn to the expectation value of an operatorÔ for a given field configuration (i.e., Eqs. (2.15) and (2.19)). By Wick's theorem, the expectation value of any N -body operator can be expressed as the sum of products of expectation values of one-body operators. Hence, the basic quantity required is (A straightforward derivation of this expansion can be found in [14].) Again, when the decomposition separately conserves proton and neutron numbers, the expectation values for proton and neutron operators are given by separate matrices.
These formulae can be extended to the case where pair operators are present in the one-body Hamiltonian; i.e., whenĥ is of the formĥ where Θ, ∆, and Λ are general N s × N s complex matrices. Our first task is to find a simple expression for ψ 0 |Û σ (β, 0)|ψ 0 whereÛ σ (β, 0) = Nt n=1 exp(−∆βĥ σn ). Using the fermion coherent state representation and Grassman algebra, we derive the expressions in Appendix A. Here we simply state the results.
If the trial wavefunction ψ 0 is a quasi-particle vacuum, such that is the matrix representing the evolution operatorÛ σ (β, 0), and M n is the 2N s × 2N s matrix representingĥ n : (4.10) Here, the evolution operatorÛ σ (τ 2 , τ 1 ) is represented by a 2N s × 2N s matrix and the many-body wavefunction is represented by a 2N s × N s matrix, independent of the number of particles present. In analogy to Eqs. (4.3-4.6), we can write To calculate the expectation value of a generalized one-body operator, we proceed as in the pure density case. Let Then the one-quasiparticle density matrix is (4.14) The final step follows from the properties of U in Eq. (A20). Note that both the overlap and the Green's function are similar to those of the density decomposition. However, the enlarged dimension of the representation matrices causes the overlap to be the square-root of a determinant rather than just a determinant. We know of no simple way to determine the sign of the square root. Computationally, it can be traced by watching the evolution of ψ 0 |Û σ (τ, 0)|ψ 0 asÛ σ is built up from 0 to β (Appendix B), although at the expense of increased computation. Also, the linear dimension of the matrix used in the calculation is twice that for the pure density decomposition case so that the pairing decomposition is more computationally demanding. However, when the interaction has a strong pairing character, it has the potential for a more effective Monte Carlo sampling, and offers greater freedom in the decomposition, which might be used to mitigate the sign problem.

B. Grand-Canonical Ensemble
For the grand-canonical ensemble, the trace in Eq. (2.18) is a sum over all possible many-body states with all possible nucleon numbers, and a chemical potential inĤ is required. For the pure density decomposition, the many-body trace is given by which can be proved by expanding the determinant [14].
For the expectation value of a one-body operator, one can recapitulate exactly the argument of the zero-temperature development and obtain We have extended these formulae to decompositions that involve pairing operators (Appendix A). The results are given in terms of the 2N s × 2N s matrices M n , U (β, 0) representing the Hamiltoniansĥ n and the evolution operator U (β, 0), (Eqs. (4.9,10)) namelŷ To motivate this formula, consider the simplest case of one time slice whereÛ = exp(−ĥ) andĥ is hermitian and in the form (4.7). Then,ĥ can be diagonalized by a HFB transformation and In the diagonal form,Tr[exp(−ĥ)] can be identified easily as due to the invariance of the determinant with respect to similarity transformations. Here, as in the overlap formulae for zero-temperature approach, the grand-canonical trace is given as the square root of a determinant, so that the evolution of the sign is important (Appendix B). The formula for the observables can be calculated from C. Activity expansion for the canonical ensemble As mentioned in the beginning of this section, the grand-canonical ensemble may lead to large fluctuations in the particle number for systems with few particles, and so is particularly ill-suited for small nuclear systems, and although the particle number does not fluctuate in the zero-temperature approach, that formalism can only give ground-state results. The canonical ensemble is therefore important for studying thermal behavior, as well as the ground state of large systems.
In the canonical ensemble, we have to find the trace ofÛ σ (β, 0) over all states with a fixed particle number N v (actually fixed proton and neutron numbers). We discuss two methods of achieving this: the activity expansion presented here and the integration over real coherent states presented in the following subsection.
Consider first the case when only density operators are present in the one-body Hamiltonianĥ. From the grandcanonical ensemble formulae (4.15), one can see that the canonical trace is just the sum of all the N v × N v subdeterminants. More explicitly, we consider the activity expansion: for some parameter λ, let (4.23) In our matrix representation,T r λNÛ σ (β, 0) = det(1 + λU ) . If Eq. (4.24) is expanded in powers of λ, the coefficient of λ Nv is just the canonical trace ofÛ σ (β, 0) over N v particles, so that det(1 + λU ) is the generating function for the canonical trace. Thus, The trick now is to find simpler expressions for ζ Nv , instead of doing the explicit sum over the determinants. To do this, write This expansion converges to the generating function because Z(β, λ) is a polynomial in λ of finite order (i.e., N v can be at most N s ). The coefficient of λ n in the exponential is readily found. For a given particle number N v , we need only calculate matrix traces up to Tr[U Nv ] and the coefficient of λ Nv can be extracted accordingly.
For one-body expectation values, using again the grand-canonical trace as the generating function for calculating observables and collecting all terms with coefficient λ Nv , we arrive at The expectation value of two-body operators, ρ ab ρ cd (σ) Nv , is nontrivial as Wick's theorem must be modified, but again one can simply collect the terms with coefficient λ Nv and obtain This form of the activity expansion works best for N v ≤ N s /2 (mostly empty model spaces). When N v > N s /2 (mostly filled spaces), it is more efficient to expand in the activity of the hole states. In this case, we define and as in Eq. (4.24),T The coefficient of λ N is the partition function for N holes, the expectation value of a one-body operator is and the expectation of a two-body operator is When pairing operators are involved, we use Eq. (4.17) for the grand-canonical trace, which becomes the generating function for the corresponding canonical trace. As an illustration, we display the formula for the expansion in particle activity,T where the definitions of S 11 , S 12 , S 21 , S 22 are obvious, and The expectation values of one-and two-body operators can then be derived as in the pure density decomposition.

D. Canonical ensemble via coherent states
We make use of the operatorĴ which can be shown [15] to be a resolution of unity in the Hilbert space of N v fermions occupying N s levels. Here, X ph are the real integration variables, |X are the real coherent states |0 is the N v -fermion state with the levels 1, . . . , N v occupied and C is a normalization constant, where dX = ph dX ph . The thermal canonical expectation is then We can compute the overlap and the observable in (4.41) as in the zero-temperature formalism, but we now have to do the Monte Carlo integration over the fields X ph as well. However, the number of these fields scales only as N 2 s , which is about the same as the number of σ fields in one time slice, so that introducing the coherent-state integration is not much more computationally expensive than the zero-temperature formalism. The advantage of coherent states is that we do not have to find Tr[U (β, 0) n ]. However, the integration of the coherent states may aggravate the sign problem, as will be discussed in Sections VI and VII.

V. DYNAMICAL CORRELATIONS
In the previous section, we discussed how to extract the expectation value of a one-body operator Ô and, by use of Wick's theorem and its extension to the canonical case, equal-time two-body operators ÂB , etc. A great deal of the interesting physics, however, is contained in the dynamical response function Ô † (t)Ô(0) whereÔ(t) = e −iĤtÔ e iĤt . In our calculations, we evaluate the imaginary-time response function Ô † (iτ )Ô(0) and from it deduce the associated strength function, fÔ(E).
In the zero temperature formalism, the strength function is where i is the ground state, while in the canonical or grand-canonical ensemble, Thus, the imaginary-time response function is related to the strength function by a Laplace transform, Recovering the strength function from R by inversion of the Laplace transform is an ill-posed numerical problem. Different methods have been proposed to surmount this difficulty [16,17]. We resort to making the best "guess" for the strength function via the Classic Maximum Entropy (MaxEnt) method [18,19], which was first introduced to recover strength functions in Monte Carlo calculations by Silver et al . [17] and has since been widely used in similar condensed matter simulations. It is in essence a least-squares fit biased by a measure of the phase space of the strength function.
In MaxEnt methods, the function to be minimized is 1 2 χ 2 − αS, where χ 2 is the usual square of the residuals, S is the entropy of the phase space (not to be confused with the action in the auxiliary field path integral), and α is a Lagrange multiplier. In Classic MaxEnt α is determined self-consistently. The method, described briefly in Appendix C, also can yield error estimates for the extracted strength function. Upon definingÔ we have an expression suitable for Monte Carlo evaluation, We now proceed to findÔ σ (τ ). For the purpose of illustration, we show the formulae for the pure density decomposition (formulae for the general case can be derived similarly). For the simplest case whenÔ = a † i orÔ = a i , it can be shown that [14] In this way, the creation and annihilation operators can be 'propagated' back to τ = 0. Any operatorÔ σ (τ ) can be first expressed in a σ (τ ) and a † σ (τ ) and therefore can be propagated back and expressed in a and a † . For example, supposeÔ = C ij a † i a j is a one-body operator. Then Thus, the response function can be measured in the same way as the static observables.

A. Monte Carlo Methods
Monte Carlo evaluation of the path integral requires a weight function. We have tried two different choices for the weight function, each with advantages and disadvantages.
One choice for the weight function is a Gaussian, with the the static mean-field solution as the centroid, and the widths given by the RPA frequencies. Thus, much of the known physics is embodied in the weight and the Monte Carlo evaluates corrections to the mean-field + RPA approximation. Further, it is simple to efficiently generate random field configurations with a Gaussian distribution.
Unfortunately, Gaussian sampling has several disadvantages. First, finding the RPA frequencies can be extremely time consuming since we have to calculate and diagonalize the curvature matrix ∂ 2 ∂σα,i∂σγ,j ∂ 2 S ∂σ∂σ at the mean-field solution. (Here, α, γ run through the number of quadratic terms in the Hamiltonian and i,j run through the number of time slices.) For any general interaction, the curvature matrix has a large dimension, N 2 s N t . Second, the Gaussian has to be modified when there is spontaneous symmetry breaking in the mean fields (such as quadrupole deformation). Otherwise, the Goldstone modes in the the RPA frequencies (e.g., zero eigenvalues corresponding to shape rotations) will destroy the normalizability of the weight function. Finally, multiple mean-field solutions well separated from each other can also pose a problem, so that multiple Gaussians may be needed.
A more satisfactory route is to choose | exp(−S)| = G(σ)|ζ(Û σ (β, 0))| as the weight function and to use a stochastic random walk such as the Metropolis algorithm to generate the fields. The expectation of an observable is then given by Eq. (2.21) where The Metropolis algorithm is free of the disadvantages for the Gaussian weight function, in that it will eventually sample the entire space where | exp(−S)| is significant. The main disadvantages of Metropolis are its inefficiency as currently implemented and the "sign problem." Let us define the denominator of Eq. (2.21) by Φ : For the Metropolis algorithm, we take random steps in the fields time-slice by time-slice, following a sweeping procedure introduced by Sugiyama and Koonin [4]. For the Monte Carlo results to be valid, one requires that the points in the random walk be both distributed as the weight function and be statistically independent. The first requirement translates into starting the fields in a region of statistically significant weight; if the initial configuration has a small weight, a number of initial "thermalization" sweeps are usually needed to relax the fields to this region. The second requirement means that the walker must sweep through the fields many times to de-correlate the samples.
The number of thermalization sweeps and decorrelation sweeps increases with system size, and the choice of random walk procedure greatly affects the sampling efficiency. In the early stage of investigation, we allowed all fields σ α at one time-slice to change with equal probability within a certain step size. The acceptance is then determined by the ratio of the weight | exp(−S)| of the old and new configurations, and we adjusted the step size to give an average acceptance probability of 0.5. The calculations of sd-shell nuclei described below then needed approximately 2000 thermalization sweeps and up to 200 decorrelation sweeps between samples. In our later work, we used another random walk algorithm, in which a fixed number of fields at one time slice are randomly chosen for updating; those chosen are generated from the Gaussian distribution in Eq. (2.9b) while the others are kept at their previous values, and the acceptance is determined by the ratio of ζ in the new and old configurations. This alternative algorithm needs only some 200 thermalization sweeps and 10 decorrelation sweeps; i.e., it is 10 times more efficient than the previous method. In addition, boundaries where | exp(−S)| = 0 can confine the first random walk algorithm, while they do not affect the second.
The random walk can thermalize faster if it starts from a configuration where the weight W (σ) is large. Usually we start from the static mean fields, given by σ α,n =σ α , n = 1, N t ; (6.3a) One can show that for canonical and grand-canonical calculations, the static mean fieldσ α is a saddle point of the weight function, For the zero-temperature approach, we also choose the self-consistent field solutionσ α , although Eq. (6.4) is not rigorously satisfied. This is preferable to starting the fields at zero, which may be far from configurations of significant weight.
The mean field solution (6.3) can be found iteratively For the systems we have tested, Eq. (6.5) usually converges within 100 interactions starting fromσ σ = 0. For larger systems and at lower temperature, convergence is slow and sometimes unstable or barely stable. In that case taking for a starting configuration also leads to a shorter thermalization time than the choice of σ = 0.

B. Examples
We now describe several examples to demonstrate our methods for calculating nuclear properties. Two major considerations arise in the implementation: the choice of decomposition scheme and the choice of ensemble. Different decomposition schemes involve different dimensions of matrices and numbers of fields, which control the computational speed. Also, the rates of convergence as ∆β → 0 are different and determine the number of time slices to be used. Finally, and perhaps most importantly, this choice also affects the "sign problem" associated with the statistical stability of the calculation. Different decomposition schemes will be compared in Examples 1 and 2 below.
The choice of ensemble among zero-temperature, canonical, and grand canonical ensembles usually do not affect the issues noted above. Instead, it depends on the kind of properties to be calculated. The zero-temperature formalism with a good trial wavefunction is effective in projecting out the ground state and is suitable for calculating ground-state static observables. For calculating finite temperature properties, the canonical ensemble is physically most relevant but also is most difficult to implement. Examples 1, 2, and 3 below demonstrate the grand canonical, zero-temperature, and canonical ensembles, respectively.
Finally, we choose a particular nucleus, 20 Ne, to demonstrate the calculation of the dynamical responses for different operators and recover the strength function by the MaxEnt Method. The examples shown below were done with 3000 to 6000 samples generated on the passively parallel Touchstone Gamma and Delta computers at Caltech.

Example 1: Monopole pairing interaction in the sd-shell
We have described the considerable flexibility in writing the two-body interaction in quadratic form; e.g., density versus pairing decomposition and direct versus exchange decomposition. To illustrate this flexibility, we consider protons only in the sd-shell (N s = 12), and keep only the J = 0 terms in Eq. (3.5); the values of V J=0 and single particle energies are taken from the Wildenthal interaction [1]. We first recouple the Hamiltonian into a quadratic form in the density operators in Eq. (3.21); because all possible density operators are required, there are 144 fields for each time slice. We then use the pairing decomposition in Eq. (3.41); after diagonalization only 6 fields are required. Finally, we write the two-body interaction asĤ 2 = 1 2Ĥ 2 + 1 2Ĥ 2 and decompose the first half using densities and the second half using pairs; the total number of fields in this case is 150. It turns out that this system is 99% free of the sign problem (i.e., Real(Φ i ) > 0 99% of the time), independent of the decomposition.
All three calculations were performed in the grand-canonical ensemble using a Gaussian weight function around the static mean field, at an inverse temperature of β = 1 (here and henceforth, we measure all physical energies in MeV) and fixed chemical potential. The expectation value of the proton number, energy, and J 2 are given in Fig. 1 as a function of ∆β. As the number of time slices increases (i.e., ∆β → 0), all three decompositions converge to the exact answer, demonstrating their mutual equivalence in the continuum limit. Note, however, the different rates of convergence.

Example 2: 24 Mg with schematic forces
Next we consider sd-shell nuclei with a schematic Pairing + Multipole density interaction, where the multipole force is separable; it is the same interaction used in [7]: where This Hamiltonian was also decomposed in several ways. We first decomposed the pairing interaction as pair operators and the multipole force as density operators. In this way, the number of fields is kept to a minimum (only 21 per time slice). The pair-operator P Tz =0 mixes protons and neutrons and therefore one matrix representing the mixed neutron and proton wavefunction is needed. In addition, the pairing decomposition requires matrices whose dimension is 2N s , so that the matrices involved in this calculation are 48 × 48.
We calculated 24 Mg (4 valence protons and 4 valence neutrons) in the zero temperature formalism; i.e., using the evolution operator at large β to project out the ground state from a trial state ψ 0 . Sinceĥ is hermitian (hereρ has the propertyρ † K,M =ρ K,−M (−1) M ), in the static path ζ = ψ 0 | exp(−βĥ)|ψ 0 is always positive definite and Φ i = 1. However, for calculations with more than one time slice Φ becomes very small, so that we can only obtain sensible results in the SPA. These results turned out to be extremely good, relaxing to the right energy and angular momentum (Fig. 2). However, the success of the static path is case-specific and not well understood. For example, we have also tried the case of just multipole interactions among 4 protons in the sd-shell ( 20 Mg), and find that the SPA relaxes to an energy 2 MeV higher than the ground state.
In a second scheme, we transformed the pairing part of the interaction (6.7) into a quadratic form in the density operators. In this transformation, we used only density operators that conserve proton and neutron numbers (3.24-3.34), and chose all E K,T =1 elements in Eq. (3.27) to be zero so that only isoscalar density operators were present in the quadratic form.
In this case, the interaction is in a much more complicated form due to the Pandya transformation of the pairing interaction, and one needs 144 fields per time slice (as compared to 21 fields needed in the first decomposition). However, an advantage lies in the separation of the Slater-determinants for protons and neutrons since only density operators that separately conserve neutron and proton numbers are present. In the sd-shell, the dimension of matrix involved is only 12 × 12. For this particular interaction, an even more desirable property of the second scheme is that the eigenvalues λ K,α in Eq. (3.34) found by diagonalizing E K in Eq. (3.27) satisfy sign(λ K ) = (−1) K+1 . We prove in the next section that this property guarantees ζ(σ) to be positive definite for even-even nuclei if a suitable trial state is chosen. This allows calculations with any number of time slices that are free from any sign problem.
We performed the calculation using the zero-temperature formalism at different β and ∆β values, choosing first the Hartree Slater-determinant as the trial wavefunction. The SPA calculation and that with ∆β = 0.125 are shown in Fig. 2. We have also performed calculations at ∆β = 0.5 and ∆β = 0.25. These results are not shown but lie between the SPA and the ∆β = 0.125 results, and show a convergence to the result at ∆β = 0.125. At ∆β = 0.125, the energy relaxes to the right energy, whereas the SPA also relaxes, but to a much higher energy.
We repeated these last calculations with a different trial wavefunction: a Slater determinant of the orbitals (j, m) = ( 5 2 , ± 1 2 ), ( 5 2 , ± 3 2 ). Although a different relaxation curve is traced out by the results at ∆β = 0.125, it also converges to the same exact result (Fig. 2). The choice of the trial state is therefore important for determining the rate of relaxation of the zero-temperature approach, although the final result is independent of the trial state as expected. In this case, although the Hartree state is lower in energy than the maximal prolate state (compare H at β = 0), it contains some high angular momentum components (compare J 2 at β = 0), so that it relaxes more slowly and reaches the ground state at a larger value of β.

Example 3: Canonical calculations of 20 Ne and 24 Mg
Next, we demonstrate the canonical ensemble for the same interaction (6.7) using the pure density decomposition as described in Example 2. We calculate 20 Ne because it is small enough to allow for an exact diagonalization to give all the states ofĤ, since we are concerned with both ground-state and thermal properties. The canonical trace for this path-integral in also positive definite (see Section VII), allowing the calculations to be done using many time slices.
The results for calculations with ∆β = 0.25, 0.125, and 0.0625 are shown in Fig. 3. The convergence as a function of ∆β is also apparent. Note that particular attention should be given to the extrapolation at high temperature. However, it is not hard to increase the number of time slices to decrease the error at high temperature. For example, for β = 0.5, ∆β = 0.0625 amounts to only 8 time slices. Similar results on 24 Mg in the canonical ensemble are shown in Fig. 4. The relaxation to the ground state can be compared with the zero-temperature result in Fig. 2; however, now the results at small values of β are also physically significant.
In the calculation of 20 Ne, the activity expansion method is numerically stable. However, instabilities appear for sd-nuclei when the number of proton or neutron valence particles (or holes) is greater than 4. The instability is signalled by the deviation of N p and N n from integers. In those cases, the expansion in Eq. (4.27) or (4.35) gives the canonical trace as the small difference between very large numbers, so that mid-shell nuclei cannot be calculated directly by those equations. We have developed an alternative Fourier expansion technique to extract the canonical trace that gives satisfactory results [12].
The real coherent-state method for the canonical trace for 24 Mg gives the same results as the activity expansion. However, Φ i is not always unity due to the need to integrate over the coherent states in Eq. (4.40) (as will be explained in the next section). It changes from 0.70 to 0.23 for β changing from 0.5 to 1.0 at ∆β = 0.25.
We have also studied rotating nuclei in the canonical ensemble by adding a Lagrange multiplier to the Hamiltonian, Ĥ ′ =Ĥ − ωĴ z , whereĴ z is the z-component of the angular momentum and ω is the cranking frequency. The calculations at β = 1 for 20 Ne with ∆β = 0.125, 0.0625, and 0.03125 are shown in Fig. 5, where the convergence to the exact results can be seen. The moments of inertia fitted from the three sets of data are 5.10, 5.30, and 4.95, compared to 4.74, the value from the exact curve. By adding the term linear inĴ z , we break the time reversal symmetry ofĥ(σ), which is related to the sign properties of the weight function ζ(σ). Φ decreases from 1 to 0.55 when ω is increased from 0 to 1.5 at ∆β = 0.125 while it decreases from 1 to 0.52 at ∆β = 0.03125.

Example 4: Response and strength functions for 20 Ne
Finally, we demonstrate calculation of the imaginary-time response functions and the reconstruction of the strength functions. The calculation for 20 Ne is done by the activity expansion method in a pure density decomposition; the Hamiltonian is that of Eq. (6.7). The canonical ensemble is more suitable than the zero temperature formalism for measuring the dynamical response, because in the latter case many "boundary" time slices are needed to project out the ground state on both the left and the right, and extra time slices would have to be introduced in the middle to measure the response. In contrast, the cyclic property of the trace allows full use to be made of every time slice in the canonical ensemble. We choose β = 2.5 (the system has already reached the ground state at this low temperature, as can be seen from Fig. 3) and calculate the response functions at ∆β = 0.125 and ∆β = 0.0625 for several one-body operators: the isovector-and the isoscalar-quadrupole operatorsQ v =Q p −Q n ,Q s =Q p +Q n , and the isoscalar and isovector angular momentum operatorsĴ s =Ĵ p +Ĵ n ,Ĵ v =Ĵ p −Ĵ n . We choose this latter 1 + operator purely out of convenience;Ĵ s is just the total angular momentum, which we verified to produce a constant response equal to Ĵ 2 , which follows from the rotational invariance ofĤ.
The response functions are shown in a semi-log plot in Fig. 6(a,c,e). For these Hermitian operators, R(τ ) is symmetric about τ = β/2, so it is shown only for τ ≤ β/2. The slope of the plot is approximately the energy of the dominant strength. TheĴ v andQ v responses relax more rapidly than does theQ s response, indicating that the two isovector operators couple to states with higher excitation energy than does theQ s operator. Since 20 Ne has a J = 0, T = 0 ground state, the states excited by an operator will carry the J and T quantum numbers of the operator. The plots are therefore consistent with the existence of a low-lying 2 + state.
The MaxEnt reconstructions of the most probable strength function for the different one-body operators are shown in Fig. 6(b,d,e) Table 1 for different operators and different ∆β. The extrapolated moments (∆β → 0) and the exact results for the ground state transitions are also shown.
The single-particle pick-up and stripping response functions for different orbitals are given in a semi-log plot in Fig. 7(a,c,e) and Fig. 8. The strength functions extracted from these responses are related to the excitation spectrum in the neighboring nucleus with one additional particle or hole. The stripping and pick-up responses are the same for protons and neutrons as the ground state of 20 Ne is isoscalar, and the final states have the angular momentum of the single nucleon that is added or removed. We see from Fig. 7 that both the pick-up and stripping responses for the j = 5 2 orbital and the pick-up response for the j = 1 2 orbital converge to the exact results as ∆β becomes small; the MaxEnt reconstruction of the corresponding strength functions in Fig. 7(b,d,e) also show a convergence to the exact results. The extracted moments are listed in Table 1. However, the responses for the j = 3 2 orbital show an anomalous behavior: they are close to the exact result at τ = 0, and then, with a sudden change in slope, follow the responses of the j = 5 2 orbital. A possible reason is that angular momentum is not conserved sample-by-sample in the calculation, but rather only statistically. The J = 5 2 states for 19 Ne and 20 Ne nuclei are much lower in energy than the J = 3 2 states (because the j = 5 2 orbital is lower than the j = 3 2 orbital by 5.6 MeV), so that if a small admixture of the J = 5 2 states "leaks" into the intermediate states for the j = 3 2 response, it will dominate with an exponentially growing correlation function. (The j = 1 2 orbital is much closer to the j = 5 2 orbital in comparison, only 0.8 MeV higher.)

VII. THE SIGN PROBLEM AND FUTURE DEVELOPMENTS
We have seen in the examples above that the sign properties for the function G(σ)ζ(σ) are crucial to the stability of the calculation, and may frustrate attempts to apply the Monte Carlo path integral to any general two-body interaction. In general, if we choose an arbitrary two-body interaction and arbitrarily decompose it into pair or density operators, Φ vanishes rapidly as β increases or as the number of time slices N t increases at a fixed β. For example, with the Wildenthal interaction and for any calculation for β > 1 with more than three time slices, the noise due to the sign completely swamps the signal. This "sign-problem" is well-known in auxiliary-field path-integral calculations [20] and fermion quantum Monte Carlo more generally. If there is no explicit symmetry to enforce the positivity of ζ σ , Φ decays exponentially as a function of β, and the problem is more severe for smaller ∆β, so that it is very difficult to calculate low-temperature properties.
Only a handful of interacting fermion systems are known to give rise to positive definite path-integral: the onedimensional Hubbard Model, the half-filled Hubbard Model, and the attractive Hubbard Model at any dimension and filling [3]. We will show in the following that an important class of interactions for the nuclear shell model has a positive-definite path-integral representation for even-even nuclei. It arises from the symmetry between time-reversed single-particle orbitals and may serve as a starting point in understanding and controlling the 'sign problem'.
We first define the 'time-reversed' partner of each creation and annihilation operator to bẽ due to the spin-half statistics.
The class of Hamiltonians that give rise to a positive-definite path integral are of the form where χ α > 0, e α can be generally complex, and ρ α is a general density operator of the form Note the requirement of a negative coupling constant for the density operator with its time-reversed partner and the time-reversal invariant form of the remaining one-body part.
For application to the shell model, we refer to Eq. (3.17), so that in the density decomposition, The requirement of a negative coupling constant for the density operators leads to a sign rule for the λ Kα , namely As we will show below, the Hamiltonian (7.3) gives rise to a one-body hamiltonian in the path integral,ĥ σ , that is symmetric in time-reversed orbitals. Time-reversed pairs of single-particle orbitals thus acquire complex-conjugate phases in the propagation, guaranteeing a positive definite overlap function ζ.
After the HS transformation on Eq. (7.3), the linearized Hamiltonian iŝ so that ρ andρ couple to complex-conjugate fields. (If some part of the interaction does not satisfy (7.6), then there will be terms inĥ that are of the form i(σ α + iτ α )ρ α + i(σ α − iτ α )ρ α , and the above statement is not true.) We represent the single-particle wavefunction by a vector of the form jm jm , m > 0 , (7.8) with states with m > 0 in the first half of the vector and their time reversed orbitals in the second half. Due to Eq. (7.3) and the fact that time reversed operators are coupled to complex conjugate fields, the matrix M i representingĥ i has the structure and one can easily verify that the total evolution matrix, is of the same form. Here, A, B, P, Q are matrices of dimension N s /2. One can show that this matrix has pairs of complex conjugate eigenvalues ǫ, ǫ * with respective eigenvectors u v and −v * u * . In the case where ǫ is real, it is doubly degenerate since the two eigenvectors are distinct.
For the grand-canonical ensemble, the overlap function is given by If only particle-type (neutron-proton) conserving density operators are present in Eq. (7.3), each type of nucleon is represented by a separate Slater-determinant having the structure (7.10), and therefore ζ = ζ p × ζ n > 0, since ζ p > 0 and ζ n > 0.
In the zero-temperature formalism, if the trial wavefunction for an even number of particles is chosen to consist of time-reversed pairs of single-particle states, where a, b are matrices with dimension ( Ns 2 × Nv 2 ), then Ψ † 0 U Ψ 0 is a N v × N v matrix with the structure (7.10), and the overlap function again satisfies If only particle-type conserving operators are present, then time-reversed pairs of trial wavefunctions can be chosen for both protons and neutrons in an even-even nucleus, giving rise to ζ = ζ p × ζ n > 0. Note that while the overlap function is positive definite in the grand-canonical ensemble for any chemical potential (and thus any average number of particles), it is true only for an even (or even-even if involving a type-conserving decomposition) system in the zero-temperature formalism with a suitable trial wavefunction.
In the canonical ensemble for N particles, a fixed-number trace is involved and therefore (7.14) We have found empirically that ζ is positive definite for even-even systems, although we lack a rigorous proof. A special case of the Hamiltonian (7.2) exists in which the overlap function is positive definite also for odd-odd N = Z nuclei. The required condition is that only isoscalar-density operators are present in Eq. (7.2). This leads to a further symmetry that protons and neutrons couple to the same field in Eq. (7.7) and therefore the evolution matrices satisfy U p = U n ≡ U . In the zero-temperature formalism, if we choose the trial wavefunction for proton and neutrons to be time reversed partners of each other, so that and ζ = ζ p × ζ n > 0. On the other hand, for the canonical ensemble, one can prove that Eq. (7.13) is real. Given that ζ n = ζ p , ζ is a square and is therefore positive.
For a system with general form (7.3), if we perform the canonical trace by integration over real coherent states, then there is an extra operator, exp(X ph a † p a h ), multiplying the evolution operator. Time-reversed states couple to different real fields in the extra operator, and a sign problem arises, as seen in Example 3 above.
Cranking also causes the sign to depart from unity. In cranking, the Lagrange multiplier term ωJ z is added toĥ in Eq. (7.7), destroying the property that time-reversed operators are coupled to complex-conjugate numbers (becausẽ J z = −J z ). Notice, however, that cranking with an imaginary Lagrange multiplier iωJ z will preserve the time-reversal symmetry and will give rise to positive path-integral.
In summary, for a Hamiltonian of the form (7.3), the above proof guarantees the overlap function to be positive for any nucleus in the grand canonical ensemble, and for even-even nuclei in either the canonical ensemble or zerotemperature formalism (with suitable trial wavefunction). It also guarantees positivity for odd-odd N = Z systems when only isoscalar density operators are involved.
The Hamiltonian (6.7) satisfies (7.6) upon decomposition of the paring interaction using density operators and involves only isoscalar operators, so that ζ is positive for even-even and N = Z nuclei. For other systems, Φ decreases as a function of β. At ∆β = 0.0625, Φ = 0.4 at β = 1.5 for 24 Na and Φ = 0.2 at β = 2.0 for 23 Na. Thus, even for odd-A nuclei, the sign properties of (6.7) are still much better than that of a general interaction violating the criteria (7.6). For example, using the Wildenthal interaction, Φ drops to several percent at β = 1 for any nucleus.
Arbitrary two-body interactions do not satisfy the sign rule (7.6), and when the rule is violated, Φ rapidly decreases as a function of β. Note that monopole pairing plays an important role here. The pairing interaction can be written as It produces a constant shift in every χ α in Eq. (7.3), as may be seen also from the multipole decomposition of the pairing force: Therefore, if the pairing interaction is strong enough compared to the remaining interactions, the sign rule can be satisfied.

VIII. CONCLUSION
We have developed a general framework for carrying out auxiliary-field Monte Carlo calculations of the nuclear shell model. In this framework we evaluate ground state or thermal observables, using pairing or density fields or both.
Although the use of pairing fields naturally embodies important aspects of the residual interaction, these calculations are more difficult due to the larger matrix dimension needed and the extra effort to keep track of the sign of the overlap function as the wavefunction is propagated. Furthermore, for calculations with multiple time slices, the Monte Carlo method with a pairing decomposition suffers from severe sign problems. However, pairing fields are suitable for carrying out static path or two-time-slice calculations where the linearized Hamiltonian is hermitian, thereby enforcing the positivity of the overlap function. This can be easily verified by observing that for hermitianĥ,ĥ a , andĥ b , with real In these cases, there is no sign problem and also there is no need to keep track of the evolution of the sign.
For the density decomposition, we have found a class of interactions which give rise to a positive definite integrand upon the HS transformation. For these interactions, stable calculations can be carried out for many time slices to extrapolate to the exact results (∆β → 0). This class of interactions includes the phenomenological pairing-plusmultipole interaction used widely. We have carried out calculations with such interactions in the sd-shell, demonstrating the power of the method in calculating both static and dynamical properties in the ground state and at finite temperature; high spin nuclei were also studied by cranking. The calculations converge to the exact results (as found by direct diagonalization) with increasing number of time slices. Although the nuclear wavefunction is not found explicitly in these calculations, many nuclear properties can be obtained. For general shell model interactions, it appears that the sign or phase property of the integrand is the major factor determining successful application of the Monte Carlo sampling. Successful calculations are usually confined to high temperature studies. We have demonstrated the freedom one has in the decomposition scheme of the two-body interaction and found that the behavior of the sign can be different in the various schemes. The next crucial step is to explore whether we can manipulate these degrees of freedom to enable stable calculations of nuclei using general forces.

APPENDIX A: DERIVATION OF THE OVERLAP FORMULAE FOR PAIRING FIELDS
We consider operators of the form where eachĥ t is a quadratic operator Without loss of generality, we choose ∆, Λ to be anti-symmetric (∆ = −∆ T , etc.). We follow the development of Berezin [13], who considered the special caseÛ = exp(−iĥt) withĥ Hermitian; we take the general case.
We calculate the grand-canonical traceT (the sum is over all states of all particle number) by using the fermion coherent state (FCS) representation of unity Here ξ α are Grassman variables and the |ξ are fermion coherent states. Then We need the FCS representation ofÛ . In what immediately follows we show that if U is the matrix representation ofÛ , that is where then with and TrΘ(n)).
In this case the trace becomes a Gaussian integral over Grassman variables; the result is given beginning with equation (A29) below. However, to come to that point we must derive equations (A8-A10).
To this end, we employ the standard rules for operating on |ξ , ξ| with a † , a: Next, we derive expressions for a αÛ , a † αÛ , as linear combinations ofÛ a α ,Û a † α . Then, with the ansatz (A8) for ξ|Û |ξ as a Gaussian in the Grassman variables ξ, ξ * , we use (A11) to derive the elements B of the Gaussian given in (A9).
To do this, we introduce the operators b,b (which are not necessarily Hermitian conjugates) Then and we seek b,b as linear combinations of a † , a.
Solving the differential equation (A15), and so in general where the summation on γ is implicit. Upon inserting (A18) in (A12b), and using the ansatz (A8), one can straightforwardly derive the B's in terms of the U 's as given in (A9).
Although we do not show it in detail, we note that B is antisymmetric (B 11T = −B 11 , B 21T = −B 12 , etc.), which can be proved using and Now we must show the normalization C is of the form (A10). To do so, we find a differential equation for C. Lettinĝ we defineÛ soÛ n+1 =Û n (∆β). Taking the expectation value of (A22) between FCS's, invoking (A11) and equating the parts independent of ξ, ξ * , one obtains d dt ln C n (t) = −Tr Λ(n + 1)B 22 n (t) .
The integral (A5) is straightforward (see Berezin [13]; the magnitude iŝ The phase ofTrÛ , though critical, is more difficult to obtain (see Appendix B for details).
One can rewrite (A27) in a more compact form. The constant C contains the factor det(S 22 ) 1 2 , which can be written using Upon introducing (A28) into (A27), performing some algebra, and using relationships from (A20), one arrives at where, again, U is the matrix in (A6) representing the evolution operatorÛ .
As for the density case, one can introduce an activity expansion to project out an exact particle number Because the matrices 1 and U are of dimension 2N s × 2N s , one can write λ Ns/2 = det 1 0 0 λ 1 2

and (A30) becomeŝ
Tr λNÛ = det This can be expanded into a polynomial in λ, which then gives the canonical ensemble.
Finally, we give the expectation value of ψ t |Û |ψ t .First we note that the vacuum expectation value 0|Û |0 is the term in (A31) independent of λ, Any quasi-particle excitations can be represented as Hartree-Fock-Bogoliubov vacua for properly defined new quasiparticle operators, which corresponds to doing a similarity transformation on the matrices. If |ψ t is the vacuum to the quasi-particle annihilation operator β i ; i.e., The formulae derived in Appendix A for calculating the overlap ζ σ (β) = ψ L |Û (β, 0; σ)|ψ R in the zero-temperature formalism andTr[Û σ (β, 0)] in the thermal formalisms all involve the square root of a determinant, leaving the phase of ζ σ (β) undetermined by a factor of ±1. This ambiguity is irrelevant for the Monte Carlo random walk as we typically take |ζ| as the weight function, but the phase must be known unambiguously for calculation of observables, as important cancellations may result.
We determine the phase by following the evolution of ζ σ (τ ) and its derivatives with respect to τ as τ goes from 0 to β. For example, if ζ σ (τ ) were purely real, zero-crossing (with ζ ′ (τ ) = 0) would indicate a change in the phase by -1. The initial phase at τ = 0 is real and positive. Following the evolution is computationally expensive, but as most of the time is spent on the random walk, where the phase is irrelevant, the overall computational time is negligible.
In what follows we give the formulae for up to fourth derivatives for each of the different formalisms.

Grand Canonical Ensemble
Define then Using the abbreviation M = M k+1 , Θ = Θ k+1 , let The derivatives of ln(f ) can be expressed in terms of the matrices G and M , Then ζ σ (k∆β + t) is given by, Zero temperature formalism In this case Then and so on and all formulae are the same as in the grand canonical ensemble (B6-9), except that the matrix G is now different. In fact, in both cases G can be shown to be the matrix for the Green's function; i.e., Canonical Ensemble In Eq. (4.35), the undetermined sign involves only the vacuum expectation 0|Û (β, 0)|0 . Once it is determined, the sign forTr[Û (β, 0)] is known. We can use the equations for the zero temperature formalism to obtain the sign of 0|Û (β, 0)|0 .

APPENDIX C: MAXIMUM ENTROPY EXTRACTION OF THE STRENGTH FUNCTION
We use the Maximum entropy (MaxEnt) method to reconstruct the strength function from the response function. Here we give a brief description of the Classic MaxEnt Method; details can be found in the paper by Gull [19].
The MaxEnt method is a Bayesian approach for reconstruction of positive additive images f from noisy data. In our case the image is the strength function f (ω). The noisy data, D j = d j + η j , are the measurements of the response function R(τ ) at the discrete imaginary times, where d j = R(j∆β) and η j is the noise in the data. In the absence of any data, the most probable image is chosen to be a default model m. Skilling [18] proved that, in that case, the only consistent choice for the probability of an image f is determined up to a parameter α: pr(f ) = exp(αS(f, m))/Z(α, m) , where S is the entropy of image f relative to the default model. If the image is discretized to f i (i = 1, . . . , r), then and where the last step follows from a change of variable u i = √ f i . (Z S has no relation to the nuclear partition function.) In the presence of data, we gain some knowledge about the image. Assuming a Gaussian distribution of errors in the data D i , the probability for f is where G ij = η i η j is the correlation matrix of the errors in the data, and For a given choice of α, the most probable image can be found by maximizing αS − 1 2 χ 2 , giving rise to the term "Maximum Entropy" method. However, α is not predetermined. Rather, it is varied until χ 2 at the maximum of αS − 1 2 χ 2 is approximately equal to the total number of data D i . But this assignment of α is ad hoc and, according to Gull [19], usually leads to an underfitting of the data. In the classic MaxEnt, α is fixed by maximizing the probability of α given the data set D i and the default model m: where After α is fixed, we can find the most probable f by maximizing αS − 1 2 χ 2 , or we can find the average f by Monte Carlo sampling using the integrand of (C4) as weight function. Information about the uncertainty in f can also be obtained from Monte Carlo sampling. In the particular application at hand, we would like to know, for example, the uncertainty in the location of the peaks in f , or in the moments M n = f (ω)ω n dω/ f (ω)dω.
Let us now return to the problem of extracting the strength function. Since the imaginary time response function R(τ ) is given at discrete times, we can only allow a limited amount of parameters in the strength functionf , which we do by discretizingf (ω) to f i at ω i . To allow for some smoothness, we choose a Gaussian function centered around each ω i , rather than a delta function. Given the imaginary time response function, we bound the range of ω over which the f is significant by ω min and ω max , and choose the ω i to be evenly distributed between them. The number of f n , n ω , should not exceed the number of data D i we have, which is the total number of time slices, N t. We choose the width of each Gaussian, ∆ω, to be half of the spacing between the ω i 's. For a non-hermitian operatorÔ, the strength function f (ω) is related to f i by while for an hermitian operator, f (−ω) = e −βω f (ω) in the canonical ensemble, and we choose f (ω) = i f i e β/2(ω−ωi) e −1/2(ω−ωi) 2 /∆ω 2 + e −1/2(ω+ωi) 2 /∆ω 2 1 ∆ω √ 2π . (C10) The response R(τ ) generated by f (ω) is, for nonhermitian operators, while for hermitian operators it is We choose the default model m i , i = 1, . . . , n ω to be a constant fixed by R(0), which is related to the total strength. The data D j are the values of R(τ ) for τ = j∆β measured from the Monte Carlo sampling. The error correlation function can also be measured in these calculations.
To maximize the probability (C7), we have to know the dependence of Z S and Z Q on α. Some simplification can be had by calculating Z S in the saddle point approximation. There is a saddle point in S(u 2 , m) at u 2 i = m i with the second derivatives ∂S ∂ui∂uj | u 2 i =mi = −4δ ij , which leads to the approximate integral The condition ∂p ∂α = 0 then becomes and the average image f i is given by We do these integrals by Monte Carlo sampling of u i with exp(αS − 1 2 χ 2 ) as the weight function. The value α is determined by the self-consistent condition S α = − r 2α . When α is known, the average distribution f i = u 2 i and the uncertainty δf i = f 2 i − f i 2 is also given in the course of the Monte Carlo evaluation of (C15).