# Calculation of exciton densities in SMMC

## Abstract

We develop a shell-model Monte Carlo (SMMC) method to calculate densities of states with varying exciton (particle-hole) number. We then apply this method to the doubly closed-shell nucleus <sup>40</sup>Ca in a full $`0s`$-$`1d`$-$`0f`$-$`1p`$ shell-model space and compare our results to those found using approximate analytic expressions for the partial densities. We find that the effective one-body level density is reduced by approximately 22% when a residual two-body interaction is included in the shell model calculation.

1. Introduction

Particle-hole, or exciton, level densities enter into the description of partial decay rates in nuclear preequilibrium emission. These level densities have been modeled using analytic expressions that describe nuclear excitations in terms of the number of particles, $`p`$, and holes, $`h`$, measured from the Fermi surface, with the exciton number $`N_e=(p+h)/2`$.

For a single species of particles, Williams  derived an expression for the partial density of states given by

$$\rho _{N_e}(E)=g\frac{\left(gEG\right)^{2N_e1}}{p!h!\left(2N_e1\right)!},$$
(1)   
where $`E`$ is the excitation energy measured above the ground-state configuration, $`g`$ is the single-particle density of states, and $`G/g`$ plays the role of an effective Pauli energy with $`G=(p^2+h^2)/4+(ph)/4+h/2`$. There exist more complicated expressions that distinguish between protons (with single-particle density $`g_p`$) and neutrons ($`g_n`$) in a given nucleus, but we will not quote them here. In the most naive picture, a uniform spacing of single-particle states, $`d`$, is assumed, in which case the single-particle level density is $`g=1/d`$, measured in units of MeV<sup>-1</sup>.

Equation 1 and its neutron/proton counterpart suffer from several deficiencies including the assumption of an unlimited number of single-particle states, an inexact treatment of the Pauli principle, and the assumption of a uniform single-particle level spacing. Extensions to the basic model that ameliorate some of these effects have also been pursued. As examples, we mention Bogilia et al.  in which the energy dependence of the single-particle level spacing was included in a general way; using the equidistant single-particle picture, Kalbach  and Zhang and Yang  considered Pauli principle corrections to the state densities; and De and Hua  considered the effects of pairing in addition to the Pauli blocking on the state densities. These effects were combined and extended to non-uniform level spacings by, for example, Harangozo et al. .

A further difficulty with the simple formula is that the residual two-body interaction, which is present beyond the nuclear mean field, and which includes important contributions beyond $`J=0`$ pairing, is not incorporated. While some progress has been made to approximate the effects of the residual interaction  on the partial level densities, no interacting shell-model calculations have been performed in large model spaces that would indicate the effect of the two-body interaction on the single-body density parameter, $`g`$, nor have there been any partial density-of-state calculations in the interacting shell model.

In this article, we describe calculations that study the effects of the residual two-body interaction on the partial level densities, and present results for partial level densities in <sup>40</sup>Ca. Our approach is to study a related quantity, $`Y_{N_e}(\beta )`$, which is the ratio of the particle-hole partition functions, $`Z_{N_e}(\beta )`$, to the full partition function, $`Z_A(\beta )`$, as a function of the inverse temperature, $`\beta `$, (measured in MeV<sup>-1</sup>) in the system. We perform our calculations in a full $`0s`$-$`1d`$-$`0f`$-$`1p`$ model space using shell-model Monte Carlo (SMMC) techniques , and an interaction that describes reasonably well the low-lying spectral properties of nuclei in the $`sd`$-$`fp`$ region. We will compare our results with those obtained from eq. (1) and its proton/neutron counterparts. Finally, we will show partial densities of states for several exciton numbers. In section 2, we give an overview of our calculational method. We present results in section 3, and conclude with a brief summary in section 4.

2. Calculation of Excitons in SMMC

Investigations into both ground-state and thermal properties of nuclei have been described using the SMMC technique . This method offers an alternative description of nuclear structure properties in the shell-model context that is complementary to direct diagonalization. SMMC is designed to give thermal or ground-state expectation values for various one- and two-body observables. Indeed, for larger nuclei, SMMC may be the only way to obtain information on the thermal properties of the system from a shell-model perspective. In this method, we make use of the imaginary time many-body propagator $`\widehat{U}=\mathrm{exp}(\beta \widehat{H})`$ to calculate the expectation values. For example, the excitation energy of a nucleus is $`E(\beta )=\widehat{H}(\beta )\widehat{H}(\mathrm{})`$, where $`\widehat{H}(\mathrm{})`$ is the ground-state energy. In order to find the excitation energy of a nucleus with particle number $`A`$, we must then calculate

$$\widehat{H}=\frac{\text{Tr}\widehat{P}_A\widehat{U}\widehat{H}}{\text{Tr}\widehat{P}_A\widehat{U}},\frac{\text{Tr}_A\widehat{U}\widehat{H}}{\text{Tr}_A\widehat{U}},$$
(2)   
where $`\widehat{P}_A=\delta (\widehat{N}A)`$ projects the trace over all many-body states in the system to those states that have the desired particle number.

Two-body terms in $`\widehat{H}`$ are linearized through the Hubbard-Stratonovich transformation, which introduces auxiliary fields over which one must integrate to obtain physical answers. Since $`\widehat{H}`$ contains many two-body terms that do not commute, one must discretize $`\beta =N_t\mathrm{\Delta }\beta `$. The method can be summarized as

$`Z_A`$ $`=`$ $`\text{Tr}_A\widehat{U}=\text{Tr}_A\mathrm{exp}(\beta \widehat{H})\text{Tr}_A\left[\mathrm{exp}(\mathrm{\Delta }\beta \widehat{H})\right]^{N_t}`$ (3)   
$``$ $`{\displaystyle 𝒟[\sigma ]G(\sigma )\text{Tr}_A\underset{n=1}{\overset{N_t}{}}\mathrm{exp}\left[\mathrm{\Delta }\beta \widehat{h}(\sigma _n)\right]},`$ (4)   
where $`\sigma _n`$ are the auxiliary fields at a given imaginary time-step $`\mathrm{\Delta }\beta `$ (there is one $`\sigma `$-field for each two-body matrix-element in $`\widehat{H}`$ when the two-body terms are recast in quadratic form), $`𝒟[\sigma ]`$ is the measure of the integrand, $`G(\sigma )`$ is a Gaussian in $`\sigma `$, and $`\widehat{h}`$ is a one-body Hamiltonian. Thus, the shell-model problem is transformed from the diagonalization of a large matrix to one of large dimensional quadrature. Dimensions of the integral can reach up to 10<sup>5</sup> for systems of interest, and it is thus natural to use Metropolis random walk methods to sample the space. Such integration is most efficiently performed on massively parallel computers. Further details are discussed in Koonin et al. .

In order to obtain density-of-state information, we calculate in SMMC the expectation of the energy and integrate the thermodynamic relationship

$$E(\beta )=\frac{d\mathrm{ln}Z_A(\beta )}{d\beta }$$
(5)   
to obtain

$$\mathrm{ln}Z_A(\beta )=_0^\beta 𝑑\beta ^{}E(\beta ^{})\mathrm{ln}Z_A(0),$$
(6)   
where $`Z_A(0)=\text{Tr}_A\mathrm{𝟏}`$ is the total number of $`A`$-particle states in the system. $`Z_A(\beta )`$ and $`\rho (E)`$ are related by the inverse Laplace transform

$$Z_A(\beta )=_{\mathrm{}}^{\mathrm{}}𝑑E\mathrm{exp}(\beta E)\rho (E),$$
(7)   
which can be solved in a saddle point approximation to yield

$`\rho (E)`$ $`=`$ $`{\displaystyle \frac{\mathrm{exp}(S)}{\sqrt{2\pi \beta ^2C}}}`$ (8)   
$`S`$ $`=`$ $`\beta E+\mathrm{ln}Z(\beta );\beta ^2C={\displaystyle \frac{dE}{d\beta }}.`$ (9)   
In this expression $`C`$ is the heat capacity of the system.

For our discussion, we study <sup>40</sup>Ca. The non-interacting ground state is a filled $`sd`$-shell with no particles in the $`fp`$-shell. Our excitons are then enumerated with respect to the filled $`sd`$-shell. We may excite both protons ($`\pi `$) and neutrons ($`\nu `$) so that $`N_e=(p_\pi +h_\pi +p_\nu +h_\nu )/2=N_{fp}`$, where $`N_{fp}`$ ($`N_{sd}`$) gives the number of particles in the $`fp`$\- ($`sd`$-) shell. For example, $`N_e=2`$ includes the following particle-hole excitations: ($`0p_\pi 0h_\pi `$, $`2p_\nu 2h_\nu `$) ($`1p_\pi 1h_\pi `$, $`1p_\nu 1h_\nu `$) ($`2p_\pi 2h_\pi `$, $`0p_\nu 0h_\nu `$). Furthermore, $`0N_e24`$ since, at most, 24 particles can be excited from the $`sd`$-shell into the $`fp`$-shell.

The ratio of the partition function for $`N_e`$ excitons to the full partition function for the $`A`$-particle system may be found by introducing a second number projection operator,

$$\widehat{P}_{N_e}=\delta (N_{sd}\widehat{N}_{sd})\delta (N_{fp}\widehat{N}_{fp}),$$
(10)   
provided that $`A=N_{sd}+N_{fp}`$. In reality we perform this projection for both protons and neutrons simultaneously, but this only complicates notation and will not be discussed here. We calculate the ratio of partition functions as

$$Y_{N_e}(\beta )=\frac{Z_{N_e}(\beta )}{Z_A(\beta )}=\frac{\text{Tr}\widehat{P}_{N_e}\widehat{U}}{\text{Tr}\widehat{P}_A\widehat{U}}.$$
(11)   
Therefore, $`_{N_e}Y_{N_e}(\beta )=1`$ which we use as a convenient numerical check. We also extract the energy of the particle-hole excitations, $`E_{N_e}`$, as

$$E_{N_e}(\beta )=\frac{d\mathrm{ln}Y_{N_e}(\beta )}{d\beta }+E(\beta ).$$
(12)   
We may now employ eq. (9) for the partial density of states:

$`\rho _{N_e}(E)`$ $`=`$ $`{\displaystyle \frac{\mathrm{exp}(S_{N_e})}{\sqrt{2\pi \beta ^2C_{N_e}}}}`$ (13)   
$`S_{N_e}`$ $`=`$ $`\beta _{N_e}E_{N_e}+\mathrm{ln}Z_{N_e}(\beta );\beta _{N_e}^2C_{N_e}={\displaystyle \frac{dE_{N_e}}{d\beta }}.`$ (14)   
Here $`\beta _{N_e}=\beta _{N_e}(E_{N_e})`$ is determined by inverting the relation $`E_{N_e}=E_{N_e}(\beta _{N_e})`$, and $`C_{N_e}`$ is the heat capacity for the particular exciton number.

3. Results

We now turn to a description of our <sup>40</sup>Ca calculation in the $`0s`$-$`1d`$-$`0f`$-$`1p`$ shell model space. Our starting point for an appropriate interaction is taken from ref. . In order to obtain a microscopic effective interaction, we begin with a free nucleon-nucleon interaction which is appropriate for a description of low-energy nuclear structure. The choice made in ref.  was to work with the charge-dependent version of the Bonn nucleon-nucleon potential model as found in ref. . Standard perturbation techniques, as discussed in ref. , were employed to obtain an effective interaction in the full $`sd`$-$`fp`$ model space. Finally, the interaction was modified in the monopole terms using techniques developed by Zuker and co-workers .

SMMC calculations for realistic interactions typically have a Monte Carlo sign problem which can be overcome by an extrapolation technique discussed in ref. , and successfully applied to the $`sd`$-$`fp`$ region in . This extrapolation technique was also applied to thermal properties of nuclei , but the statistical error inherent in the energy upon extrapolation prevents a full description of the density of states unless one has good justification to spend the computational resources to reduce the statistical error. It was recently demonstrated that a good reproduction of the experimental density of states could be obtained for nuclei in the $`0f1p`$-$`0g_{9/2}`$ shell  using a pairing-plus-quadrupole interaction that was free from the sign problem. In this work, we fit our realistic two-body interaction discussed above to a pairing-plus-multipole interaction given by

$`\widehat{H}_2=g_0\pi \widehat{P}_{00}^{}\widehat{P}_{00}+4\pi {\displaystyle \underset{\nu \mu }{}}\chi _\nu :{\displaystyle \underset{\mu }{}}()^\mu \widehat{Q}_{\nu \mu }\widehat{Q}_{\nu \mu }:,`$ (15)   
where $`::`$ denotes normal ordering and $`\widehat{P}_{\lambda \mu }^{},\widehat{Q}_{\nu \mu }`$ are pair and quadrupole operators given by

$`\widehat{P}_{\lambda \mu }^{}`$ $`=`$ $`{\displaystyle \underset{ab}{}}()^\mathrm{}_b(j_a𝒴_\nu j_b)[\widehat{a}_{j_a}^{}\times \widehat{a}_{j_b}^{}]_{\lambda \mu }`$ (16)   
$`\widehat{Q}_{\nu \mu }`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\nu +1}}}{\displaystyle \underset{ac}{}}(j_ar^\nu 𝒴_\nu j_c)[\widehat{a}_{j_a}^{}\times \widehat{\stackrel{~}{a}}_{j_c}]_{\nu \mu }.`$ (17)   
In eq. (16) $`an\mathrm{}j`$ denotes a single-particle orbit and $`\widehat{\stackrel{~}{a}}_{jm}=()^{j+m}\widehat{a}_{jm}`$. We fit $`g_0`$ and the $`\mu =2,4,6`$ multipoles to the realistic interaction from . A least-squares fit gives an interaction which indeed has a good Monte Carlo sign. After some minor adjustments to the pairing strength in order to obtain a better gap between ground states and first excited states in several light nuclei, we use the following parameter set: $`g_0=0.63`$ MeV, $`\chi _2=0.047`$ MeV fm<sup>-4</sup>, $`\chi _4=0.001`$ MeV fm<sup>-8</sup>, and $`\chi _6=0.17\times 10^3`$ MeV fm<sup>-12</sup>. (Large enhancements of the $`ar^\nu b`$ matrix elements as $`\nu `$ increases is the reason for the decrease in $`\chi _\nu `$ values, although contributions to two-body matrix elements arising from the higher multipoles is significant.) Our single-particle energies are 0.0, 5.36, 0.64, 8.21, 14.21, 10.14, and 12.07 MeV for the $`0d_{5/2}`$, $`0d_{3/2}`$, $`1s_{1/2}`$, $`0f_{7/2}`$, $`0f_{5/2}`$, $`1p_{3/2}`$, and $`1p_{1/2}`$ orbitals, respectively. We do not correct for center-of-mass motion in these calculations, although such a contamination to the $`Y_{N_e}`$ should be fairly small in this system. Furthermore, we do not include odd multipoles which upon fitting were found to give coefficients that cause Monte Carlo sign problems. Thus our negative parity states are probably less well described by this choice of Hamiltonian.

We use the results of the non-interacting case to demonstrate the validity of our technique for finding the partial partition functions. In order to show this, we calculate by enumeration the total number of many-body states for each $`N_e`$. This can most easily be done by using eq. (6) to find $`Z_A(\beta =0)`$. We also find $`Y_{N_e}(0)`$ by an extrapolation from small, but finite, $`\beta `$. We show in fig. 1 our results for the number of states as a function of the exciton number. The SMMC results are compared to an exact counting of the number of states of a given exciton number. The agreement is excellent. The total number of calculated SMMC states is 3.834$`\times `$10<sup>16</sup> as compared to the exact value of 5.095$`\times `$10<sup>16</sup>.

We show in fig. 2 a comparison of the non-interacting (left) and interacting (right) calculation. The $`N_e=0`$ calculation gives some indication of the thermal freeze-out of the ground state. The non-interacting calculation requires fairly large $`\beta `$ to fully reach the ground state, since the first excited state is only 0.64 MeV above the ground state. We pursued these calculations to $`\beta =4.0`$ MeV<sup>-1</sup>, for which $`H`$ is 0.139 MeV from the ground-state value. Since it takes more thermal energy to overcome the pairing interaction and to excite them from the ground-state configuration, the interacting curves representing different $`N_e`$ excitations are compressed in $`\beta `$ relative to the non-interacting curves. (The excitation energy of the first excited state is approximately 3.5 MeV.) However, the same features remain. Clearly, as thermal energy is decreased, and $`\beta `$ becomes larger (temperature decreases), it is more difficult to produce large particle-hole excitations in the system. The converse is also true at higher temperatures, where it is difficult to obtain only $`N_e=2`$, for example. In the interacting case, the low temperature tail of the $`N_e=2`$ exciton tends to spread further in $`\beta `$ than does the $`N_e=1`$ tail. As we shall see, this has direct consequences on the partial densities of states. Furthermore, since $`_{N_e}Y_{N_e}(\beta )=1`$, we can interpret $`Y_{N_e}`$ as a measure of likelihood to find $`N_e`$ excitons at a given temperature. Since the excitation energy is a monotonic function of the temperature, one expects the density of states to be dominated by excitons of a particular type in a given energy range. As we shall see, this is indeed the case.

We compare our results to those obtained from eq. (1) and its proton/neutron equivalent by adjusting $`g`$ to obtain a fit to the calculated $`Y_{N_e}`$ curves. This is demonstrated for the $`N_e=4`$ curve in fig. 3. We fit each of our SMMC curves to eq. (1) for an effective one-body density parameter $`g_{\mathrm{eff}}=g_p+g_n`$ , and the effective level density parameter is $`a=(\pi ^2/6)g_{\mathrm{eff}}`$. A uniform Fermi gas yields $`aA/15`$ (=2.67 for A=40) MeV<sup>-1</sup>, a harmonic oscillator potential yields $`aA/10`$ MeV<sup>-1</sup> (=4.0), and the empirical value is $`A/8`$ MeV<sup>-1</sup>. We obtain $`a=5.43`$ ($`g_{\mathrm{eff}}=3.3`$) MeV<sup>-1</sup> for the non-interacting case, and $`a=4.27`$ ($`g_{\mathrm{eff}}=2.6`$) MeV<sup>-1</sup> in the interacting case, rather independent (within 0.02 MeV<sup>-1</sup>) of $`N_e2`$. For $`N_e=1`$ the comparison cannot be made as the Blann-Williams formula breaks down. Thus, $`g_{\mathrm{eff}}`$ is reduced by 22% in the presence of an interaction.

The decomposition of the $`N_e=4`$ case into the various proton-neutron components is also shown in fig. 3. The ($`0p_\pi 0h_\pi `$, $`4p_\nu 4h_\nu `$) plus its neutron counterpart carries very little weight here, while the ($`2p_\pi 2h_\pi `$, $`2p_\nu 2h_\nu `$) component of $`Y_4`$ carries the most weight. As expected, in all $`N_e`$ cases the largest component of $`Y_{N_e}`$ is the one in which the number of excited neutrons equals the number of excited protons.

Finally we present our result for the calculation of $`\rho _{N_e}`$ derived by using $`Z_{N_e}(\beta )`$ in eq. (6) and $`E_{N_e}`$ from eq. (12). The natural log of $`\rho _{N_e}`$ is shown in fig. 4 for the $`N_e`$=1–6 excitations. We also include in the figure the total state density as a function of the excitation energy in the system. The saddle point approximation breaks down in regions where there are very few states, which makes it difficult to describe well-separated states in the low-lying spectrum ($`E^{}<3`$ MeV) for the full density or for the individual exciton densities. We also propagated our statistical error bars through the calculation of $`\rho _{N_e}`$, but, as can be seen, they are very small except in the case of the $`N_e=1`$ excitons. Note that the majority of states, for example at $`E^{}=25`$ MeV, are $`N_e=4`$ states, while at $`E^{}=35`$ MeV the $`N_e=5`$ states contribute most. This localization in excitation energy of excitons was reflected in our earlier discussion of the behavior of $`Y_{N_e}`$. Interestingly, the $`N_e=1`$ density of states begins in energy slightly above 2p2h state density. Recall that experimentally the first excited state of <sup>40</sup>Ca is a 0<sup>+</sup> 2p2h state (at 3.3 MeV), and that the first negative parity state (a 3<sup>-</sup>) occurs at a slightly higher energy of 3.7 MeV. Our Hamiltonian fairly closely gives the correct relative starting energies for these two exciton configurations, although, due to the breakdown of the saddle-point approximation for low state densities, we cannot precisely determine the excitation energy of the first excited 0<sup>+</sup> level. We also note an interesting pairing effect that shows up in the partial densities. Note that the $`N_e=1`$ state density starts about 0.6 MeV above the $`N_e=2`$ case. The $`N_e=3`$ state density begins approximately 2 MeV above the $`N_e=4`$ case. This is a manifestation of pairing in the system. It takes more excitation energy to produce an odd particle-hole excitation than it does to produce an even particle-hole excitation since energy must be expended to break pairs for odd excitations. This effect was already apparent from the discussion of the low-temperature behavior of the $`Y_{N_e}`$, as indicated above.

4. Conclusion

We have discussed in this article how one may obtain information on particle-hole excitations using SMMC methods for calculations of a nuclear system. Our technique uses the ratio of the particle-hole partition function to the full partition function of the system. The method incorporates exact Pauli blocking, non-equidistant single-particle energies, and gives the exact partial densities for a given nuclear effective interaction, within statistical errors. It also has a well-defined energy scale. One drawback of the present calculation is that the space size is limited to two major oscillator shells, although this can be rather easily overcome. The projection operator introduced in this may be applied in any Monte Carlo technique where ratios of partition functions are needed. Our results also indicate that the effective $`g`$-parameter used in eq. (1) should be reduced by approximately 22% to account for the inclusion of the two-body interaction which acts to correlate the nucleus beyond the simple mean-field or pairing prescription.

The method we have described here could be further advanced in two ways. One may increase the model space used, thus allowing for a broader range of energies and excitation modes to be explored. The method is also applicable to $`mp`$-$`nh`$ excitations if we extend our studies to open-shell nuclei such as, e.g., <sup>42</sup>Ca. This may be pursued in future work.

###### Acknowledgements.

This work was supported in part through grant DE-FG02-96ER40963 from the U.S. Department of Energy. Oak Ridge National Laboratory (ORNL) is managed by Lockheed Martin Energy Research Corp. for the U.S. Department of Energy under contract number DE-AC05-96OR22464. We also acknowledge support from the U.S. National Science Foundation under Grants No. PHY-9722428, PHY94-12818 and PHY94-20470.