{'fraction_non_alphanumeric': 0.0424368486511987, 'fraction_numerical': 0.015182056010636017, 'mean_word_length': 4.191472781921407, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Perturbed Static Path Approximation at Finite Temperature: Observables and Strength Functions

We present an approximation scheme for calculating observables and strength functions of finite fermionic systems at finite temperature such as hot nuclei. The approach is formulated within the framework of the Hubbard-Stratonovich transformation and goes beyond the static path approximation and the RPA by taking into account small amplitude time-dependent fluctuations around each static value of the auxiliary fields. We show that this perturbed static path approach can be used systematically to obtain good approximations for observable expectation values and for low moments of the strength function. The approximation for the strength function itself, extracted by an analytic continuation from the imaginary-time response function, is not always reliable, and we discuss the origin of the discrepancies and possible improvements. Our results are tested in a solvable many-body model.


I. INTRODUCTION
Mean-field approximations [1,2], such as the finite temperature Hartree-Fock (HF), are standard for describing nuclei at finite temperature. Collective excitations of the system are obtained by linearizing the time-dependent HF equations around the static HF solution, leading to the finite temperature random phase approximation (RPA) [3]. However, this treatment is inadequate for situations where various nuclear configurations compete and have comparable free energies, which causes observables to fluctuate widely about their mean values. This is the case, for example, in the vicinity of a shape transition from a spherical to a deformed nucleus, where various discontinuities predicted in the mean-field approach are smoothed in the finite nuclear system by the presence of fluctuations. Large fluctuations in the nuclear shape also play an important role in giant dipole resonances (GDR) whose frequency is strongly coupled to the quadrupole deformation. In the adiabatic approximation, the observed GDR strength function is obtained by integrating the strength function that corresponds to each quadrupole configuration over all possible configurations weighted by their respective Boltzmann factor [4].
It has been shown [2,5] that a systematic description of fluctuations can be obtained in the framework of the auxiliary-field path integral (Hubbard-Stratonovich transformation [6]). In this framework the original many-body propagator (in imaginary time) is decomposed into a superposition of one-body propagators describing non-interacting fermions moving in a fluctuating external time-dependent potential (auxiliary fields). The expectation value of any observable is represented as a weighted average of its expectation values in the corresponding non-interacting systems. Applying the method of steepest descent one can obtain selfconsistent mean-field approximations at finite temperature. By taking into account small amplitude time-dependent fluctuations of the auxiliary fields around the mean-field solution in the Gaussian approximation one rederives the RPA [5,7].
The path integral representation also constitutes a starting point for exact numerical solution of the many-body problem where the auxiliary field integration is performed by Monte Carlo techniques that were recently developed for strongly correlated electron systems [8] and for the interacting nuclear shell model [9]. In the latter case a practical solution to the Monte Carlo sign problem [10], which is generic to all fermionic systems, enables the study of nuclear properties in medium and heavy mass nuclei using realistic effective interactions. The auxiliary field formulation also suggests new approximation schemes which are non-perturbative, the simplest of which is the static path approximation (SPA) [11] - [17]. Here, the path integral is approximated by summing over the time-independent fields only, weighted by the appropriate Boltzmann factor. This amounts to averaging over all possible static mean-field configurations rather than the self-consistent ones alone. The SPA has been used to calculate free energies and level densities in nuclei and was found to be superior to the mean-field approximation. In particular, it accounts for the enhancement of level density due to thermal fluctuations of the shape [13]. More recently the approximation has also been applied to the calculation of strength functions [16]. At high temperatures the SPA partition function approaches the exact result. However, as the temperature decreases the SPA becomes inaccurate since time-dependent fluctuations about the mean-field configurations can no longer be neglected. This is manifested especially in the strength function where even its first moment is significantly underestimated [16].
Recently, a method to improve the SPA has been proposed for the partition function [18][19][20][21]. In this approach contributions from time-dependent fields in the neighborhood of each static field are incorporated perturbatively to the second order in their amplitudes. When one considers such time-dependent fluctuations around the equilibrium configuration only, one obtains the RPA corrections to the partition function. However, in the proposed approach such small amplitude time-dependent fluctuations are taken around each static configuration of the auxiliary fields and the static integration is still fully retained. At even lower temperatures some of these time-dependent fluctuations can become unstable and the approximation breaks down. However, this happens only at temperatures below the phase transition, when the temperature drops below the largest imaginary RPA frequency. We remark that as a saddle point develops in the nuclear free energy surface below the transition temperature, it is still possible that the imaginary RPA frequencies are small enough in magnitude that the time-dependent fluctuations are all stable (although there is an unstable direction in the static free energy surface). This approximation scheme has been applied to free energy and level density calculations in simple models and shown to work well down to low temperatures.
It is interesting to investigate whether incorporating small time-dependent fluctuations provides a significant improvement over the SPA evaluation of quantities other than the partition function. This paper will explore the validity and applicability of this scheme, which we term the perturbed static-path approximation (PSPA), for the calculation of expectation values of observables at finite temperature as well as of strength functions. Previous work employed a formulation of the PSPA based on ordinary quantum-mechanical perturbation theory [19,20], which was specialized for free energy calculations and is not easily extended to other quantities. We therefore reformulate it in a general way, using many-body methods, and apply it to the calculation of observables and strength functions. By testing the PSPA in a solvable model we find that the PSPA results agree closely with the exact solution for observable expectation values and for low moments of the strength function. The strength function itself is approximated well at high temperatures and also at low temperatures for a certain regime of the model's parameters, although the low-temperature results generally do not improve much on the SPA. We discuss possible improvements. This paper is organized as follows. In Section II we review the auxiliary-field path integral formalism and in Section III we formulate the PSPA for the partition function. Section IV discusses the application of the PSPA to the calculation of expectation values of one-and two-body operators. Finally, in Section V we present a treatment of the strength function and its lowest two moments in this framework. We illustrate the PSPA and compare it to other approximations in a simple model whose exact solution is given in Appendix B.

II. AUXILIARY FIELD PATH INTEGRAL
In this Section we briefly review the auxiliary-field path integral representation of the imaginary-time evolution operator for an Hamiltonian with two-body interactions (Hubbard-Stratonovich transformation) [2,5]. We consider a system of interacting fermions, and assume for simplicity an Hamiltonian of the form where K and V are one-body operators a † i , a i are creation and annihilation operators for a set of single particle states i, and χ is a coupling constant. Our results can be easily generalized to a superposition of such separable interactions We remark that an arbitrary two-body interaction can be brought into the form (3) by diagonalizing u αβ , with the number q of separable interactions V 2 α equals at most to the square of the number of single particle states. In order to obtain a path integral representation of U = exp(−βH) we divide the imaginary-time interval [0, β) into N sub-intervals of length ǫ = β/N and write At each time slice n we use the identity valid for any operatorâ with a bounded spectrum and λ > 0 (we shall comment on the case λ < 0 below), to replace e 1 2 χV 2 ǫ by an integral over an auxiliary variable ξ n . This is the Hubbard-Stratonovich (HS) transformation [6], resulting in In the limit N → ∞, ξ n becomes a field ξ(τ = nǫ) = ξ n and we obtain the continuous version of the HS transformation where U ξ is the (imaginary time) propagator for a time-dependent one-body Hamiltonian and T denote time ordering. Eq.(8) provides a representation of the many-body evolution operator as an average over one-body evolutions U ξ which correspond to non-interacting particles moving in a fluctuating time-dependent field ξ(τ ), weighted by a Gaussian factor. Following [22] it is advantageous to describe the field ξ(τ ) in terms of its Fourier components, where σ −r = σ * r to keep ξ(τ ) real and ω r = 2πr/β are the Matsubara frequencies.
Here we assume the number of time slices N to be odd. Rewriting the functional integral over the field in terms of σ r , we have where The measure in (11) is given by where σ r = σ ′ r + iσ ′′ r . The one-body Hamiltonian H σ whose corresponding propagator is U σ separates into static and time-dependent parts It is useful to introduce the interaction picture representation U σ of the one-body propagator U σ with respect to the static part of the Hamiltonian h 0 where h 1 (τ ) = e τ h 0 h 1 e −τ h 0 . We then rewrite (11) as with the measure Eq. (16) describes the evolution operator for the two-body H as a Gaussian-weighted average of one-body evolutions e −βh 0 corresponding to a static field σ 0 , multiplied by a correction factor which represents the contribution from time-dependent fluctuations of the ξ-field about σ 0 . In fact, as we demonstrate below, any quantity of interest can be written as an integral over its static-field value times a correction factor. The objective of this paper is to approximate this factor for various quantities by evaluating the small-amplitude fluctuation contribution to the integral over σ r (r > 0) and to explore the validity of this approximation. For a superposition of separable interactions (3) one introduces auxiliary-field variables σ α r corresponding to each V α but the preceding development remains unchanged. We point out that for χ < 0 in (1) one should use (λ > 0) instead of (6).

III. PARTITION FUNCTION
We now consider the partition function Z ≡ Tre −βH . We work in the grand canonical ensemble and set the chemical potential µ = 0 to keep the notation simple. It is convenient to choose a σ 0 -dependent basis for the Fock space in which h 0 is diagonal: We can write Z using (11) as or in a form more convenient for our purpose (using (16)): Here is the partition function corresponding to the static part h 0 and is the correction factor to ζ 0 due to the time-dependent fluctuations of the field about σ 0 . We use the notation O 0 ≡ Tr(e −βh 0 O)/ζ 0 to denote the thermal average of an observable O with respect to h 0 . The effective static-field free energy F 0 is defined by (21) to be The representation (21) of Z is a starting point for various approximations. The meanfield approximation (MFA) is obtained when the contribution of the time-dependent paths to Z is neglected by setting h 1 = 0, implying ζ ′ 0 = 1, and the integration over the static fields σ 0 is performed in the method of steepest descent. This amounts to approximating the path integral by the contributions of the static pathsσ 0 that minimize the free energy F 0 in (24). It is easy to show [2,5] that this minimization condition is where f i = (1 + e βǫ i ) −1 are the Fermi occupation numbers, and that the solutionσ 0 of (25) is the Hartree mean field. One can improve on the MFA result by performing the integration over σ r in the expression (23) for ζ ′ 0 also by steepest descent. The saddle point is now given byσ 0 , σ r = 0 and one obtains the finite temperature RPA corrections to the partition function [5].
The static-path approximation (SPA) [11] - [17] is obtained by again setting h 1 = 0 but now the integration over σ 0 in (21) is performed exactly, thus approximating the path integral by the contributions from all static paths: The SPA is expected to become exact at high temperatures since one can use the one time slice approximation in (7) with an error of O(β 2 ) that vanishes as T → ∞. This method is advantageous to the MFA since it takes into account exactly large amplitude static fluctuations around the mean field. However, it neglects the time-dependent fluctuations which constitute the RPA corrections. This shortcoming of the SPA can be remedied if the exact integration over the static paths is supplemented by evaluation of ζ ′ 0 in the limit of small amplitude time-dependent fluctuations about each static value σ 0 . This scheme, the perturbed static-path approximation (PSPA), is the focus of our work. It was introduced in Refs. [18] - [21] for the partition function Z. Our approach is different and has the advantage that it can be generalized to the calculation of observables and response functions as shown in the following sections. We first illustrate our method by rederiving the corresponding expression for the partition function. We expand log U σ 0 in (23) log with V (τ ) = e τ h 0 V e −τ h 0 in the interaction picture. The time-ordered averages are calculated using the finite-temperature Wick's theorem [1,2] where the subscript c means that only the connected diagrams should be summed up in the diagrammatic representation of this expansion. The first term in (27) vanishes since V (τ ) c is τ -independent. For the second term Wick's theorem gives where the unperturbed temperature Green's function g 0 i is given by with frequencies ν k = (2π + 1)k/β. Using (28-29) we obtain for the double integral in (27) The infinite sum over k in (30) is calculated using the frequency summation technique [1,2]. The essence of this method lies in the observation that the points z = iν k are the poles of the function −β/(e βz + 1) with residue of one and the sum therefore equals the contour integral where C encircles the imaginary axis. This contour can be continuously transformed into a circle centered at the origin of an arbitrarily large radius, which is deformed at two places to include the poles at z = ǫ i and z = ǫ j − iω r . The residue theorem then gives with ∆ ij = ǫ i − ǫ j . Using (32) in (27) we have with The correction factor ζ ′ 0 in (23) is now given by a Gaussian integral and we obtain The second equality in (35) defines the frequencies Ω ν (β, σ 0 ) through [20] 1 where the prime in ij ′ and ij ′ restricts the product or sum to pairs (i, j) that satisfy i < j and ∆ ij = 0. Note that there are as many Ω ν in the numerator of (36) as there are ∆ ij in the denominator. The third equality in (35) uses the infinite product representation sinh x = x r>0 (1 + x 2 /π 2 r 2 ). Together with (21) and (22) we finally have This is a closed-form expression which corresponds to the limit N → ∞ (N is the number of imaginary-time slices, see (5)) since the infinite product in (35) has been performed exactly. Hence the PSPA result does not contain errors originating from a discretization of [0, β) into sub-intervals of a finite length, as is the case in a Monte Carlo evaluation of the path integral.
It can be shown [20] that for the value σ 0 =σ 0 (β) which minimizes F 0 (β; σ 0 ) (see Eqs. (24) and (25)), these roots ±Ω ν are the RPA frequencies, i.e. the frequencies of small amplitude oscillations around the equilibrium configuration. For an arbitrary static σ 0 , the frequencies ω = ±Ω ν solve the generalized finite temperature RPA equations obtained by replacing the equilibrium configuration by the arbitrary σ 0 (see Appendix A). We note that the Gaussian approximation in (23) leads to a convergent integral only if 1 + a r > 0 for all r, i.e. −Ω 2 ν < ω 2 r for all r and ν. When all RPA frequencies are real, these conditions are always met. In particular this is the case if σ 0 is a local minimum of the static free energy i.e. a stable mean-field configuration σ 0 (Thouless theorem [23]). At zero temperature ω r = 0 and an imaginary RPA frequency would lead to a breakdown of the Gaussian approximation. However, at finite temperature imaginary RPA frequencies do not necessarily lead to instability. If the largest modulus of all imaginary RPA frequencies is below ω 1 = 2πT , then the quadratic fluctuations in σ r are still stable. An instability in the σ 1 direction occurs when the magnitude of one of the imaginary RPA frequencies crosses 2πT . The RPA frequencies depend both on β and σ 0 . In particular, for any temperature β −1 there exists a static field σ 0 = σ ′ 0 (β) such that Ω 2 ν (β, σ ′ 0 ) = −ω 2 1 = −(2π/β) 2 for some ν, causing the correction factor ζ ′ 0 to diverge. Note that at T = 0 instability occurs as soon as some Ω ν becomes imaginary whereas at T > 0 an instability develops only for a sufficiently imaginary Ω ν . These instabilities can be ignored when they occur at fields for which the static free energy is large. In practical applications the PSPA breaks down only at temperatures that are significantly below the shape transition temperature, when the saddle point in the static free energy surface becomes unstable to small amplitude time-dependent fluctuations.
In the case (3) where the interaction is a sum of several separable interactions we have an auxiliary-field variable σ α r for each V α . Assuming χ α = χ (this is always possible by redefining V α ), expression (33) becomes The correction factor ζ ′ 0 is then generalized from (35) into where a r is the q-dimensional matrix defined in (39). As shown in Appendix A where Ω ν are again the RPA frequencies at finite temperature. We then obtain for the PSPA partition function an expression similar to (37), except that the integral over σ 0 is replaced by an integration over q static fields σ α 0 . When one or more of the χ α in (3) are negative and a representation of the type (18) is used in the HS transformation, the one-body Hamiltonian h 0 in (19) becomes non-hermitean, so that its eigenvalues ǫ i (σ 0 ) are in general complex and the associated one-body partition function can be negative. It was recently pointed out [24] that the static fields which correspond to these "repulsive" terms in the interaction do not represent large amplitude thermal fluctuations, and can simply be treated in a saddle point approximation. This amounts to keeping the exact integration over the static "attractive" fields but taking the mean-field solution for the static "repulsive" fields for every configuration of the "attractive" fields in the integrand.
An alternative way to approach "repulsive" interactions is to subtract from the Hamiltonian a term proportional toN 2 (N is the particle number operator) with a coefficient large enough so that the two-body part of the Hamiltonian becomes a negative-definite quadratic form (as in the case where all χ α in (3) are positive). However, by doing this we also increase the magnitude of "attractive" terms and it remains to be investigated whether the SPA and PSPA would still work well for such modified Hamiltonians. Note that the subtracted term is just a constant for a fixed number of particles.
To illustrate our results we apply the formalism to a simple many-body model (a variant of the Lipkin model [25]) based on a U(2) algebra which is defined and solved in Appendix B. At each σ 0 we have Thus the single-particle Hamiltonian corresponding to h 0 has two g-fold degenerate levels and the matrix where Ω 2 = 4ǭ 2 − (4gǫ 2 χ/ǭ) tanh(βǭ/2), so there are only two RPA frequencies ±Ω. The saddle-point mean-field equation (25) is and has a single spherical stable solution σ 0 = 0 for β < β c , which becomes unstable for β > β c and bifurcates into two deformed minima σ 0 = ±σ 0 (β). Here β c is the transition temperature The dimensionless parameter κ thus determines the mean-field behavior of the model. For κ < 1 there is no phase transition and σ 0 = 0 is the only solution at all temperatures.
In Fig. 1 we present the MFA, SPA and PSPA results for the free energy F (β) ≡ −β −1 log Z in our U(2) model using ǫ = 1 and g = 10 and compare them against the exact result. We consider three cases characterized by different values of κ. For κ = 0.5 there is no mean-field transition whereas for κ = 1.5 and κ = 3.0 the signature of a transition at temperatures β c = 1.61 and β c = 0.693, respectively, is seen in the mean-field curve. The SPA overestimates the exact result especially at low T whereas the PSPA is quite accurate at all temperatures. All three approximations converge to the exact result as T → ∞. We find that for β sufficiently small the divergence of ζ ′ 0 occurs only at a very large σ 0 = σ ′ 0 (β), where the other factors in the integrand of (37) are vanishingly small, and therefore does not constitute a practical problem. The breakdown of the PSPA occurs when β becomes large enough such that σ ′ 0 (β) ∼ 0. This happens only at temperatures far below the transition temperature T c = 1/β c .
In order to study the effect of the time-dependent fluctuations we compare in Fig. 2 the effective free energy F 0 (β; σ 0 ) in (24) as a function of the static field σ 0 with the SPA free energy (ζ ′ 0 = 1 in (24)) at different temperatures. We consider the case κ = 1.5. Above the mean-field transition temperature β c = 1.61, ζ ′ 0 ≈ 1 and the two approximations yield similar results. As the temperature is lowered, time-dependent fluctuations deepen the freeenergy minimum at σ 0 = 0. Below T c these fluctuations also lower the barrier between the two mean-field configurations σ 0 = ±σ 0 . Thus the PSPA result for the free energy improves significantly the SPA result especially below the transition temperature (see the middle panel of Fig. 1). The MFA result for the free energy F (β) is the most inaccurate especially near the transition where the free energy has comparable values over a broad range of configurations σ 0 .
We remark that in spite of its simplicity, the U(2) model has features that are generic to more realistic nuclear interactions (e.g. quadrupole interaction). In particular, the phase transition in its mean-field theory is analogous to the shape transitions from deformed to spherical shapes that occur in deformed nuclei [26]. We therefore expect the PSPA to improve significantly on the SPA also for more realistic nuclear shell model interactions. For interactions that include both pairing and multipole components (e.g. pairing plus quadrupole model), the SPA works better when a mixed pairing-density decomposition is used rather than just a pure density decomposition [9]. A mixed decomposition is thus preferable when the PSPA is applied to such interactions.

IV. THERMAL EXPECTATION VALUES OF OBSERVABLES
In this Section we consider an observable O and treat its expectation value at finite temperature O = Tr(e −βH O)/Tr(e −βH ) in the framework of the PSPA. Although we discuss observables of the form O = D and O = D 2 for a one-body operator D, any n-body O can be treated in a similar fashion.

A. One-Body Observables
In constructing an auxiliary-field path integral representation for D we can use Eqs. (11) and (20) to get Alternatively, to facilitate our calculations we can use (16) and (21) to write where The effective static-field free energy F 0 (β; σ 0 ), defined in (24), has been calculated in the previous Section (see Eqs. (22) and (35)). Our aim here is to calculate D 0 . In a diagrammatic representation of the perturbation series for U σ D 0 the connected diagrams can be easily shown to factor out at each order [1,2]: Thus we can write D 0 itself as a weighted average of U σ D c : U σ 0 is given in (33) and (34), whereas for U σ D c we expand Notice the absence of first-order terms in σ r since they vanish upon the integration in (52). Before continuing we introduce the following notation. Using (23) and (33) we define a matrix A(σ 0 ) by where the vector σ has N − 1 components (σ ′ 1 , σ ′ 2 , · · · , σ ′′ 1 , σ ′′ 2 , · · ·) and A is diagonal with elements (1 + a 1 , 1 + a 2 , · · · , 1 + a 1 , 1 + a 2 , · · ·) with a r given in (34). We also use (53) to define b(σ 0 ) and a matrix B(σ 0 ) by Using this notation we find Tr Recall that at each σ 0 we choose a different basis (19) for the Fock space in which we write The zeroth-order term is then In order to obtain an expression for Tr(A −1 B) in (56) we calculate the double integral in (53) using Wick's theorem and the frequency summation technique [1,2]. Notice that since A is diagonal it is sufficient to consider the case s = −r. The calculation is similar to the one performed in the previous Section and yields the result where In the cases ǫ j = ǫ i etc. it is understood that I ijk r ≡ lim ǫ j →ǫ i I ijk r . The sum over r = 0 in (59) which consists of N − 1 elements can be expressed in a closed form for N → ∞ by applying once more the frequency summation technique. Writing a −1 r = a −1 (iω r ) and I ijk r = I ijk (iω r ), we observe that the sum is carried out over the points z = iω r which are the poles of the function β/(e βz − 1) with residue of one. Therefore where the contour C encircles the imaginary axis. To evaluate the integral we transform C into an arbitrarily large circle, deformed to include the poles at z = ∆ ij , z = −∆ jk , and z = ±Ω ν for all ν, and obtain Here and the prime in lm ′ restricts the product to pairs (l, m) satisfying l < m and ∆ lm = 0.
In the case of q > 1 separable interactions in (3) the calculation is more complicated since the matrix A(σ 0 ) in (56) is block-diagonal with blocks of dimension 2q (see Eqs. (38)-(40)). However, the general form (59) including the summation over r still holds and the subsequent use of frequency summations to get a closed-form result can be carried out in this case as well.
We test the approximations (62) by applying them to the calculation of J z in our U(2) model ( J x = J y = 0 both exactly and in the above approximations). The matrix corresponding to D = J z is block-diagonal with g 2 × 2-blocks d ij = 1 2 cos 2φ sin 2φ sin 2φ − cos 2φ (64) in the notation of (43)-(44). Thus The more complicated expression for J z (P SP A) is similarly obtained using (43)-(45) and (62)-(64). The SPA and PSPA results are presented in Fig. 3 together with the meanfield calculation derived from a steepest-descent evaluation of (65). Above the transition temperature all approximations agree fairly well with the exact result (B8). For lower temperatures the MFA and SPA significantly overestimate the exact result whereas the PSPA works well down to the breakdown temperature.

B. Two-Body Observables
For simplicity we assume a two-body observable of the form O = D † D where D is a one-body operator. The calculation of the expectation value of such an observable, albeit more complicated, can be carried out along the same lines. As for D we have where Following similar steps we get an expression of the form (56) with the same A(σ 0 ) but different B(σ 0 ), b(σ 0 ). After further manipulations we obtain the final result ijkl in (69) are given by To test the approximations (68) we apply them to the calculation of J 2 x , J 2 y and J 2 z in our U(2) model. The matrices corresponding to D = J z , J y are block-diagonal with g 2 × 2-blocks is given in (64). We therefore have Expressions for D 2 (P SP A) can similarly be derived using (43)-(45) and (68-71). These results are presented in Fig. 4 together with mean-field calculations for different values of κ. The MFA and SPA exhibit large deviations from the exact results in most cases even above the transition temperature and when no transition occurs. The PSPA shows the best agreement and in the only case where its deviation is appreciable ( J 2 y at κ = 1.5), the other approximations give qualitatively wrong results. Note that the quantitative differences between the various approximations depend on the observable being calculated. In the notransition case κ = 0.5 all approximations agree for J 2 z (as they do for J z ) but not for J 2 x , J 2 y . The expectation value of the Hamiltonian H = 2ǫJ z − 2χJ 2 x (B1) itself is shown in Fig.  5 (note that the additivity of the average, O 1 + O 2 = O 1 + O 2 , is preserved in MFA, SPA and PSPA). The PSPA result is still superior but the agreement of the SPA with the exact result is also quite good, contrary to what it predicts for J z and J 2 x when taken separately.

A. Linear Response Theory
In this Section we are interested in the response of the system to an external perturbation D as reflected by its thermal strength function where | n and E n are the many-body eigenstates and corresponding energies. G(ω) is the Fourier transform of the real-time response function G(t − t ′ ) which originates in linearresponse theory [1,2] where the subscript H refers to the Heisenberg picture D H (t) = e iHt De −iHt . If an equilibrated system with a Hamiltonian H is perturbed by some external potential V (t) at t = t 0 , the change in the thermal expectation value of an observable D with respect to the equilibrium situation is where is the retarded real-time response function. The Fourier transform of the latter is related to that of δ D (t) by The strength function can also be obtained directly from G R (ω): This suggests that in order to approximate G(ω) we should approximate G(t) or G R (t) and Fourier-transform the result. This has been done in [16] where G(ω) has been obtained from a static-path calculation of G(t). However, it is very difficult to incorporate the contribution of time-dependent paths into the real-time formalism. In contrast, the imaginary-time framework discussed in the previous Sections is quite suitable for this task. Thus instead of G(t) we consider the imaginary-time response function [1,2] where D H (τ ) = e τ H De −τ H . G is periodic with period β (G(τ + β) = G(τ )), hence we are interested in its Fourier coefficients where ω n = 2πn/β are the Matsubara frequencies [1]. G R (ω) and G n are related by an analytic continuation since their Lehmann representations are determined by the same weight function ρ(ω). Defining it is easily verified that as η → 0 + . We can therefore calculate G R (ω) by obtaining an expression for G n with an explicit dependence on iω n , then perform the analytic continuation by formally replacing iω n → ω + iη. We can then calculate the strength function G(ω) from (79). In general, an analytic continuation of G n is not unique since it is based on extrapolating from a discrete set of points iω n to a continuum ω. The sum rule ∞ −∞ dωρ(ω) = D † D , however, selects a continuation that satisfies Γ(z) ∼ D † D /z as z → ∞ which is unique [1,2]. In this Section we present an auxiliary-field path integral treatment of G(iω n ) ≡ G n from which we extract several approximations for G(ω).

B. Strength Function in the PSPA
In the auxiliary-field path integral representation Eqs. (80) and (81) can be written in the form where The time-dependence in D is understood to be in the interaction picture with respect to h 0 , D(τ ) = e τ h 0 De −τ h 0 . Our aim is to calculate G 0 (iω n ). The connected diagrams factor out as in (51) [1,2]: where terms proportional to δ n,0 , such as are omitted since they do not contribute to the analytic continuation. Similarly to (56) we can rewrite G 0 (iω n ) in the form Tr where A(σ 0 ) is defined in (54) and b(σ 0 , iω n ), B(σ 0 , iω n ) are given implicitly by (86). The where the second sum is over the values σ 0 = σ λ 0 satisfying ∆ ij (σ λ 0 ) = ω. In the case of our U(2) model ∆ ij = 0, ±2 ǫ 2 + χ 2 σ 2 0 , hence for | ω |< 2ǫ the strength function G (SP A) (ω) = 0 and the transitions at this ω-range are not reflected.
It is interesting to compare the shell-model Monte Carlo (SMMC) methods with the PSPA. While in the SMMC the auxiliary-field path integral is evaluated exactly (except for statistical errors), the problem of extracting the strength function from the imaginarytime response function is quite difficult due to statistical noise. The strength function is calculated in the Monte Carlo using a maximal entropy reconstruction method [8,9], but this method work well only in some cases. In contrast, the PSPA strength function is extracted by exact analytic continuation (although only within the approximation). An additional advantage of the PSPA is that the infinite-discretization limit of the imaginary-time interval [0, β) is taken exactly, whereas in the SMMC it is necessary to extrapolate the finite time step to zero. The validity of the PSPA for strength functions is tested in Sections V.C and V.D.

C. Moments of the Strength Function
The moments of the strength function provide another measure of the quality of the approximations we develop. Rather than integrating over the expressions we have for G(ω), it is more convenient to obtain the moments directly in terms of H and D. From G(ω) being the Fourier transform of G(t) it follows that and differentiating the definition of G(t) (second Eq. in (74)) we get for the lowest moments Thus the zeroth moment (total strength) is simply a two-body expectation value, which was discussed in the previous Section and found to be well reproduced by the PSPA (contrary to the SPA) in our U(2) model (see Fig. 4). The first moment M 1 (M 1 /M 0 is average transition energy) is also given by the expectation value of a two-body operator. Taking D = J x in our U(2) model we can exploit the angular-momentum commutation relations to get M 1 = 2iǫ J x J y . The latter is calculated by a generalization of the results for D † D obtained in the previous Section to the case D † 1 D 2 : whereas d ij in (90) is the matrix corresponding to D 2 , d † ij should be taken to be that of D † 1 . The results for M 1 are shown in Fig. 6 for different values of the mean-field parameter κ in (47). The PSPA agrees well with the exact result. The SPA is good at high temperatures but worsens appreciably near the mean-field transition and below it, except for the no-transition case (κ = 0.5) where it remains a good approximation also at low temperatures. This comparison is interesting in view of the results of Ref. [16] where two versions of the static approximation for the strength function were tested. The first version (called the adiabatic approximation) is identical to our SPA. The second version (called the static-path approximation in Ref. [16]) consists of estimating the path integral representation of G(τ ) in (80) by step function paths ξ(τ ′ ) which have a discontinuity at τ ′ = τ , rather than by the constant ones (σ 0 = ρ 0 ) alone. We shall discuss this approximation further below. The PSPA includes the effect of discontinuous paths only to the second order in σ 0 − ρ 0 through the correction factor from small oscillations about the average value σ 0 τ /β + (ρ 0 − σ 0 )(1 − τ /β). It was found in Ref. [16] that although the inclusion of the discontinuous paths (102) provided an improvement over the SPA, a good agreement with the exact results was not achieved, contrary to the situation in the PSPA case. This suggests that the important contribution to the moments of the strength function beyond the SPA comes from small-amplitude oscillatory paths rather than from step function paths with a large discontinuity.

D. Strength Function in a Simple Model: Results and Discussion
We test the PSPA for G(ω) in our U(2) model with perturbing operators D = J x and D = J y . The results are presented in Fig. 7 and Fig. 8, respectively. We consider the same three cases studied previously characterized by κ = 0.5, 1.5, 3.0 at different temperatures, chosen to be at the mean-field transition in each case as well as above and below it (a transition does not occur for κ = 0.5). Note that the exact result (B9) consists of a sum of δ-functions in the limit η → 0 + and is therefore singular, as is the case for the MFA result. This is not the situation in the SPA and PSPA which involve an integration over the static field σ 0 . In order to facilitate a meaningful comparison we keep η small but finite both in the exact result, which becomes a sum of Lorentzians of width η, and in the approximations. Finally, we plot G(ω) only for positive ω since G(−ω) = e −βω G(ω) for an Hermitean D.
The shortcoming of the SPA expressions (97)-(98) is manifested clearly in Figs. 7 and 8. Since we use ǫ = 1 the SPA strength function cannot reflect transitions with ω < 2. Hence it vanishes in the range ω < 2 even if most of the strength is concentrated there, as is the case for κ = 0.5, 1.5 at β = 1.7, 3.0 (where the SPA result has a shifted peak near ω = 2 to the right of the exact peak). Furthermore, since the PSPA result consists of an additive correction to the SPA expression (see Eq.(89)), this shifted SPA peak leaves its trace in the PSPA strength function. For κ = 0.5, for instance, even though the main PSPA peak is in excellent agreement with the exact one, it is accompanies by a small additional (false) peak to its right, left over from the SPA, which becomes larger as the temperature decreases. In general, the PSPA works quite well for small κ but decreases in quality and becomes comparable with the SPA as κ increases or the temperature decreases. The MFA generates a sharp peak located near the middle of the broader SPA peak, consistently with its origin as a steepest-descent approximation of the SPA integral.
We mentioned above that in Ref. [16] G(ω) was calculated in a modified static-path approach which included constant paths with a discontinuity at τ in the path integral representation of G(τ ). It is interesting to note that the resulting strength function in [16] is accurate for large values of κ and deteriorates as κ decreases, in contrast with the PSPA result. This suggests that for small κ the major contribution to the strength function beyond static fields comes from small oscillations about them. However, for large κ these small oscillations are negligible and large imaginary-time discontinuities in the static fields become important. It is therefore desirable to have an approximation scheme which takes both contributions into account. In the following we discuss an approach to this problem in the imaginary-time framework and the difficulties it encounters.
We start from the auxiliary-field path integral representation of the imaginary-time response function However, rather than use the Fourier decomposition (10) of ξ(τ ′ ) over the entire interval [0, β) to obtain the form (84), (85) from which the SPA and PSPA are derived, we divide the intervals [0, τ ) and [τ, β) into N and M sub-intervals, respectively, and use separate decompositions in each. We have with the reality condition σ −r = σ * r , ρ −r = ρ * r and τ -dependent frequencies ω r = 2πr/τ , ν r = 2πr/(β − τ ). In terms of the new variables the representation (103) becomes The discontinuous static-path approximation (DSPA), originally introduced in [16] using the real-time framework, is now obtained by neglecting the contribution of the oscillations about the static paths (102) which results in the two-dimensional integral The partition function in this approximation in given by Note that Z (DSP A) acquires a τ -dependence.
In order to carry out the imaginary-time technique of Fourier-transforming G(τ ) to get G(iω n ) and extract the strength function G(ω) by an analytic continuation as was done above, it is necessary to obtain the functional dependence on τ in (106) analytically. However, the τ -dependence of the traces involved is non-trivial and had to be studied numerically in Ref. [16] even for the simple U(2) model; note that the analogous situation in the real-time framework (namely that the t-dependence of G (DSP A) (t) is not given analytically) does not pose a problem since a numerical Fourier transform produces G (DSP A) (ω) directly. Furthermore, unlike the SPA case where we had a static Hamiltonian h 0 and Wick's theorem could be used to calculate the traces, here we have a discontinuous Hamiltonian for which Wick's theorem is not applicable. In particular, without a generalization of Wick's theorem to this situation, it would be difficult to use our methods to calculate corrections due to small oscillations about the discontinuous paths.

VI. CONCLUSION
In this paper we present an approximation scheme, the PSPA, for the calculation of thermodynamic quantities and finite temperature response functions in finite fermionic systems. The approximation is derived in the framework of the auxiliary-field path integral. We use an imaginary-time formulation which facilitates the extension of this approximation to physical quantities beyond the free energy and the level density to which it was previously limited.
Testing the PSPA in a simple many-body model, we find that it improves on the SPA and is a good approximation for expectation values of observables as well as for low moments of strength functions. This indicates that the contribution of time-dependent fluctuations about the static fields (neglected in the SPA) is significant. The required computational work involved in the PSPA includes a q-dimensional numerical integration over the static fields σ α 0 and a diagonalization of a q × q matrix at each σ α 0 -point, where q is the number of separable interactions in the Hamiltonian. This approximation breaks down at low temperatures when the small-oscillation correction factor diverges for the dominant static fields, indicating that large time-dependent fluctuations become important. However, the breakdown occurs at temperatures well below the mean-field transition and does not affect the usefulness of the PSPA except at very low temperatures. For the strength function itself the PSPA results become less reliable when the contribution of static paths with large discontinuity (at the time where the response function is calculated) is important. Further improvement would require the inclusion of both discontinuous static paths and small time-dependent oscillations around them.
It would be interesting to test the PSPA methods for more realistic nuclear interactions, such as pairing plus multipole interactions. For these interactions, the SPA works better in a mixed pairing-density decomposition than in a pure density decomposition. Thus, it would be useful to extend the present PSPA techniques to such a mixed pairing-density decomposition in the HS representation.
This work was supported in part by the Department of Energy Grant DE-FG02-91ER40608.

APPENDIX: RPA FREQUENCIES AT FINITE TEMPERATURE
The finite temperature RPA equations for a general two-body interaction u ijkl are given by [3] − where ∆ ij = ǫ i − ǫ j and f i are the Fermi-Dirac occupation numbers f i = (1 + e βǫ i ) −1 . The solutions of (A1) are the RPA frequencies Ω ν and ξ ν are the associated RPA amplitudes. The single-particle energies ǫ i in (A1) correspond to the mean-field solutionσ 0 but in the following we replaceσ 0 by a general static field σ 0 . For a separable interaction as in (1) with V Hermitean, or more generally for an interaction as in (3) which is a sum of q such separable terms, u iljk = − α χ α v α ij v α lk , and Eq. (A1) can be rewritten as (A2) can be converted to a set of coupled equations for the q quantities by multiplying Eq. (A3) by v α ′ ji (f i − f j ) and summing over ij for each α ′ = 1, . . . , q. We get For (A4) to have a non-trivial solution (i.e. not all ζ ν α = 0), we require that the determinant of the coefficient matrix vanishes Regarding the l.h.s. of Eq. (A5) as a function of ω, where we have substituted ω for Ω ν , we notice that its roots are ±Ω ν , while its poles are ±∆ ij . It then follows that .

APPENDIX: THE MODEL
The formalism developed in this paper is illustrated and tested in a simple Fermionic system, a variant of a model introduced in [25] which is based on a U(2) algebra and is therefore solvable. This is a two-level system where each level is g-fold degenerate and may therefore contain between zero and 2g Fermions. The Hamiltonian is given in terms of quasi-angular momentum operators which has the form (1) with K = 2ǫJ z and V = 2J x . The quasi-angular momentum operators are given by a † 1i a 1i + a † 2i a 2i (B2) withN being the particle-number operator. The 2 2g states are arranged in U(2)-multiplets (n, j) where the quantum numbers are the number of particles n = 0, ..., 2g and the quasiangular momentum j = 0, ..., n/2 or 1/2, ..., n/2 (depending on whether n is even or odd). Each multiplet (n, j) contains 2j + 1 states labeled by | njm with m = −j, ..., j are the eigenvalues of J z . To find the number d n (j) of (n, j)-multiplets we first observe that the number of states with given n and m is N n (m) = g n/2 − m g n/2 + m .
Since every multiplet (n, j ′ ) with j ′ ≥ j contributes a single state with J z = j, we also have N n (j) = n/2 j ′ =j d n (j ′ ) .
(B4) Therefore d n (j) = N n (j) − N n (j + 1) = g n/2 − j g n/2 + j − g n/2 − j − 1 g n/2 + j + 1 , which checks to give n/2 j=0(1/2) (2j + 1)d n (j) = 2g n , the total number of states with n particles. The Hamiltonian matrix in this basis is block-diagonal with d n (j) identical blocks of dimension 2j + 1 for each pair (n, j), whose diagonalization gives the energies E jm and corresponding eigenstates. The exact partition function in the grand canonical ensemble is then Z = The strength function associated with an operator O is given by which reduces to a sum over δ-functions as η → 0 + .