text:
# Box- and peanut-shaped bulges Based on observations collected at ESO/La Silla (61.A-0143), DSAZ/Calar Alto, and TIRGO/Gornergrat.

## 1 Introduction

Our statistical study of $``$ 1350 edge-on galaxies derived from the RC3 (*Third Reference Catalogue of Bright Galaxies*, de Vaucouleurs et al. rc3 (1991)) (Lütticke et al. 2000a , hereafter Paper I) has revealed that 45 % of all bulges (S0-Sd galaxies) are box- and peanut-shaped (b/p). This high frequency must be considered in a theory of bulge formation and especially in scenarios explaining the formation process of b/p bulges.

Currently discussed theories of bulge formation are primordial or secular evolution scenarios. While monolithic collapse (Eggen et al. ELS (1962) (ELS)), clumpy collapse (Kauffmann et al. kau (1993), Baugh et al. bau (1996)), and inside-out formation (van den Bosch bosch (1998), Kepner kep (1999)) are proposed for primordial bulge formation, merging processes (Wyse et al. wys (1997) and references therein) and dynamical evolution triggered by bars (e.g. Combes et al. com90 (1990)) are the basis of secular evolution scenarios. However, combinations of both scenarios are also possible (Combes com99 (2000)).

In the last 20 years several working hypothesis about the origin of b/p bulges have been developed. External cylindrically symmetric torques can act on an initially non-rotating triaxial spheroidal stable system so that it becomes box-shaped (May et al. may (1985)). However, such torques are “problematical to relate to some astrophysical counterparts” (Combes et al. com90 (1990)). Mergers of two disk galaxies are proposed by Binney & Petrou (bin (1985)) and Rowley (row (1988)) as a possible formation process of b/p bulges. However, such a scenario requires very special conditions such as the precise alignment of the spin and orbital angular momenta of the two galaxies (Bureau bur (1998)). A more common process to form b/p bulges could be the accretion of satellite galaxies, although an oblique impact angle of the satellite is needed (Binney & Petrou bin (1985), Whitmore & Bell whi (1988)) and the mass of the accreted satellite is constrained. If it is too large, the stellar disk would be disrupted (Barnes bar (1992), Hernquist her (1993)). The origin of b/p bulges by bars was first proposed by Combes & Sanders (com81 (1981)). This theory was strongly supported ten years later by the studies of Combes et al. (com90 (1990)), Raha et al. (rah (1991)), and Pfenniger & Friedli (pfe (1991)) dealing with N-body simulation of disk galaxies.

The connection between b/p bulges and the presence of a bar is difficult to prove by surface photometry, because b/p bulges are observable only in almost edge-on galaxies. However, there is photometric evidence in a few edge-on galaxies from cuts parallel to the major axis (de Carvalho & da Costa deca (1987), Dettmar & Ferrara df (1996)) and in two intermediately inclined galaxies (Bettoni & Galletta bet (1994), Quillen et al. qui (1997)) pointing to bars. Additionally, kinematic observations of around 20 b/p bulges with characteristic bar signatures in velocity fields and in orbit analyses (e.g. Kuijken & Merrifield kui (1995)) are reported. In addition, our statistics (Paper I) have shown that the overall frequency of barred galaxies ($``$ 55 %, Paper I, Knapen et al. kna2000 (2000)) is able to explain the high fraction of b/p bulges.

On the other side observational evidence for accreted material in b/p bulges is found by Bureau & Freeman (bf (1999)) for a few b/p bulges. Indirect evidences are substructures and asymmetries in some special b/p bulges (Dettmar det89 (1989), Dettmar & Barteldrees 1990a , Lütticke et al. 2000b , hereafter Paper III). Combes et al. (com90 (1990)) and Shaw (sha93 (1993)) conclude from their investigation of b/p bulges, listed by Jarvis (jar (1986)), de Souza & dos Anjos (sa (1987)), and Shaw (sha93 (1993)), that they show less evidence for significant merging and are likely not the result of recent accretion events. However, b/p bulges can indirectly originate from interaction and finally accretion by triggering the formation of bars (Noguchi nog (1987), Gerin et al. ger (1990), Walker et al. wal (1996)), which then actually build up the b/p structures (Fisher et al. fis (1994), Mihos et al. mih (1995)). Therefore the search for bars is one of the best ways to discriminate between the proposed evolutionary scenarios of b/p bulges.

The essential advantage of investigations in the NIR is the small extinction by dust at these wavelengths (Rieke & Lebowsky rie (1985), Knapen et al. kna91 (1991)). Therefore the influence of the dust lane in the plane of edge-on galaxies is largely reduced by NIR observations, especially in the $`K`$-band. In this way the structure of a galaxy near the plane in radial as well as in vertical direction can be analysed and thereby bars detected (Dettmar & Ferrara df (1996)). In addition it is even possible to determine parameters of bars, e.g. their thickness and projected length. Therefore it is possible to investigate the connection between bars and shapes of bulges.

## 2 Observation and data reduction

Observations were carried out at the 1.5m TIRGO/Gornergrat in 1994/95 (run identification: TI) with NICMOS3 ARNICA (0.92$`\mathrm{}`$/pixel), the 1.2m telescope at Calar Alto in June 1996 (CA4) and June 1997 (CA2) with MAGIC (Rockwell NICMOS3 array) using the high resolution mode (1.2$`\mathrm{}`$/pixel), the 2.2m telescope at Calar Alto during two runs in January 1997 (CA3) and March 1998 (CA1) also equipped with MAGIC (hr mode: 0.64$`\mathrm{}`$/pixel), and at ESO/La Silla with the 2.2m telescope in May 1998 (E1) with IRAC-2b (NICMOS3 array) using the LC objective (0.49$`\mathrm{}`$/pixel). To the resulting sample of 56 galaxies four galaxies kindly provided by Emsellem (E2: ESO/La Silla 2.2m telescope, April 1998, IRAC-2b, LC objective) are added. The sample was selected from the list of edge-on galaxies with bulge classification (Paper I) according to the allocated observing time, matching the field size, and covering all bulge types. Filters and typical on-source integration times for each observing run are listed in Table 1 and the final sample is given in Table 2. At ESO/La Silla and TIRGO the $`K^{}`$\- resp. $`K`$-band filters were used. Due to large thermal emission in $`K`$ for MAGIC we observed at Calar Alto mainly in $`H`$. Using different filters we can show that the bulge shape in the NIR is independent from the observed wavelengths (Paper I).

Data reduction is performed with standard techniques for NIR images: dark subtraction, sky background subtraction, bad pixel correction, and flatfielding (CA1-4: sky flats, E1+2: dome flats) (see e.g. MAGIC Observer’s Guide 111http://www.caha.es/CAHA/Instruments or ESO-MIDAS User Guide, Volume B 222http://www.eso.org/projects/esomidas). The sky background frames are determined by offsets, after each object frame, integrating sky images with the same time as used for the object images. Combining the individually reduced object frames leads to the final image of the galaxy. For the photometric calibration of the images infrared standard stars (Elias et al. eli (1982); Persson et al. per (1998)) are used giving a photometric error of $`\mathrm{\Delta }0.3`$ mag. The data obtained at TIRGO are calibrated using published integrated aperture data.

## 3 Analysis of the data

Following Wakamatsu & Hamabe (wak (1984)) “bumps”, which are symmetric to the minor axis, in cuts along or parallel to the major axis of edge-on disk galaxies can be associated with bars. This is indicated by the fact that the bump becomes more and more insignificant at larger distances from the plane (D’Onofrio et al. dono (1999)). Using this technique bars in edge-on galaxies NGC 4762 (Wakamatsu & Hamabe wak (1984)), NGC 1381 (de Carvalho & da Costa deca (1987)), NGC 5170 (Dettmar & Barteldrees 1990b ), NGC 4302 (Dettmar & Ferrara df (1996)), NGC 3250B, NGC 5047, NGC 6771, IC 4767 (Lütticke lue99 (1999)), and at a marginal level in NGC 128 (D’Onofrio et al. dono (1999)) have been found. However, for edge-on galaxies these radial cuts are obviously hampered by dust near the plane. Therefore previous studies using optical observations are limited to early type galaxies with low dust content or the resulting evidence for a bar remain uncertain (NGC 5170). The influence of dust is significantly reduced in the NIR enabling the identification of bars also in late type spirals (Dettmar & Ferrara df (1996)).

Therefore we have searched for bumps, which are signatures of a bar, in radial cuts of the galaxies in our NIR-sample. Additionally, the shape of the bumps in the cuts along and parallel to the major axis can be divided in profiles with a flat part and in profiles decreasing from the center of the galaxy throughout. For these purposes we have determined the surface brightness at the radial position of the top of the bump ($`R_{\mathrm{tb}}`$) and at the end of the central bulge ($`R_{\mathrm{CBU}}`$). We call the feature in the cuts “strong bar signature”, if the gradient

$$_{\mathrm{bump}}=\left|\frac{\mu (R_{\mathrm{tb}})\mu (R_{\mathrm{CBU}})}{R_{\mathrm{tb}}R_{\mathrm{CBU}}}\right|$$

is smaller than 0.015 mag arcsec-1 (Fig. 1 and 5). Otherwise the bump is called “weak bar signature” (Fig. 2). The error of this gradient is $`\pm `$0.003 mag arcsec-1. An example for a galaxy without a bar signature is given in Figure 3.

For galaxies showing a bar signature several characteristic quantitative parameters are derived from the cuts parallel to the major axis (for illustration see Fig. 1):

* BAL: The projected bar length is marked in the cut along the major axis by an increasing light distribution towards the center compared to the radial exponential light distribution (Pohlen et al. poh (2000)) of the disk.
* BAT: The bar thickness is defined by the vertical extent of cuts showing the bar signature.
* CBU: We introduced the length of the central bulge as the length of the structure being brighter than the bar signature.
* BUL: The bulge length is marked by the increasing light distribution over the exponential disk well above the bar.
* BPL: The length of the b/p structure is distance between the maxima of the b/p distortion. This parameter is obviously only definable for peanut bulges.

However, taking into account seeing effects the measured values are only upper limits.

For 17 galaxies of our sample it is not possible to search for bar signatures. Sometimes the heavy dust extinction near the plane, even in the NIR, prevents a final analysis. In other galaxies extra features (e.g. stars or possible spiral arms), asymmetrical structures, or low signal-to-noise ratios would lead to a too uncertain verification of a bar. Therefore the resulting sample is reduced to 43 galaxies.

## 4 Results

The derived bar signatures and parameters for our sample are presented in Table 2. Furthermore, we have computed characteristic ratios of the parameters (Table 2).

In Paper I we have visually classified bulge shapes in three types of b/p bulges (1: peanut-shaped, 2: box-shaped, 3: nearly box-shaped) and one elliptical type (4). Regarding the ratio of the projected bar length to the bulge length (BAL/BUL) we measure a clear trend with the bulge type (Fig. 4): The more prominent the b/p type, the larger the ratio. The mean value for type 1 bulges is BAL/BUL $`=1.9\pm 0.3`$, type 2: $`1.6\pm 0.2`$, type 3: $`1.4\pm 0.2`$, and type 4: $`1.1\pm 0.1`$. However, we are dealing with small number statistics, $`\mathrm{N}=21`$, and the scatter is relatively high. This could be explained by the different morphological types contributing to these statistics, because early type galaxies have smaller ratios than late types (Fig. 4). Additionally, according to Martin (mar (1995)), Elmegreen & Elmegreen (2elm95 (1995)), and Elmegreen et al. (elm (1996)) the relative bar length ($`R_{\mathrm{bar}}/R_{25}`$) has also a high scatter pointing to a large range of bar lengths.

Two other characteristic ratios, CBU/BPL and BAL/BPL, are only related to peanut bulges. They can be compared with predictions of N-body simulations and thereby test the evolution scenario of b/p bulges proposed in these simulations (see next section). The mean value of BAL/BPL is $`2.7\pm 0.3`$ (Table 2). A surprising result for CBU/BPL is the mean value of 1.0 with the small scatter of 0.1 (Table 2 and Fig. 1).

Five out of the six galaxies without a bar signature (Table 3), but with a b/p bulge are thick boxy bulges (TBBs) (Dettmar & Lütticke dl (1999)). This group of galaxies shows remarkable morphological features such as large asymmetries and substructures and the box-shape is more prominent in the outer part of the bulge, possibly pointing to a different formation scenario. Therefore this group of b/p bulges will be studied separately in a forthcoming paper (Paper III). The remaining galaxy (NGC 3762) without bar signature and b/p bulge is classified as a weak type 3 bulge on an optical image (Paper I). Excluding the five TBBs from Table 3 there is an almost one-to-one correspondence between b/p bulges and bar signatures. The fraction of b/p bulges with such a signature is 95 % whereas over 80 % of the non-b/p bulges show no bar signature. Additionally, it should be mentioned that the three galaxies (NGC 4026, NGC 4762, and UGC 9858) with a bar signature and without a b/p bulge have the smallest BAL/BUL ratios in our sample (Table 2 and Fig. 4). Table 3 reveals a significant correlation on one side between strong bar signatures and prominent b/p bulges (type 1 + 2) and on the other side between weak bar signatures and less pronounced b/p bulges of type 3. This strong dependence is shown on a 1.0 % significance-level by a $`\chi ^2`$-test. We do not find a correlation between the Hubble type and the weak or strong bar signature.

For NGC 1175 (Fig. 5) and NGC 3203 we measure a strong bar signature. Therefore the classification as unbarred systems within the RC3 can be replaced. Since the bar classification for edge-on galaxies is difficult and lead to misclassifications in galaxy catalogues based on photographic plates, our method of bar detection points at a solution for this problem.

## 5 N-body simulation

### 5.1 Model

The origin of b/p bulges based on the evolution of bars is shown in several N-body simulations (e.g. Combes et al. com90 (1990), Pfenniger & Friedli pfe (1991)). To check the consistency of simulations with the results of our observations, we have investigated a representative simulation of a barred galaxy. For this purpose a particle position file of the final stage of a N-body simulation of a barred galaxy in three dimensions was kindly made available by Daniel Pfenniger. The simulation uses a particle-mesh code with 200.000 particles described in Pfenniger & Friedli (pfe (1991)). They choose axisymmetric initial conditions which are likely to develop a bar instability on their own and an initial mass distribution corresponding to a superposition of two axisymmetric Miyamoto-Nagai disks. The provided binary file represents the final time step of their simulation at T = 5000 ($``$ 5 Gyrs) in which a bar is formed. The length unit is chosen to be one kpc (Pfenniger & Friedli pfe (1991)). For the transformation of the three-dimensional mass distribution of the particles to a two-dimensional surface brightness distribution we assume a constant M/L ratio and take into account the inclination ($`i`$) of the galaxy and the aspect angle ($`\varphi `$) of the bar. An edge-on bar is denoted with $`\varphi =90^{}`$, while $`\varphi =0^{}`$ refers to an end-on bar.

### 5.2 Classification of the bulges

In the investigated N-body simulation the peanut shape of bulges can be detected down to a galaxy inclination of $``$$`75^{}`$, if the bar is edge-on. If the galaxy is exactly edge-on, peanut bulges can be seen down to an aspect angle of the bar of $``$$`77^{}`$. In general the simulation shows that the degree of boxiness of the bulge, reflected by the bulge type, is weakened with decreasing inclination of the galaxy (Fig. 6, Table 4). At inclinations lower than $``$ 75 it is not possible to distinguish between structures perpendicular to the galactic plane and the bar itself. This value was also previously found by Shaw et al. (sdb (1990)) as limit for the observation of b/p bulges. Regarding galaxy inclinations above this limit the b/p morphology is detected in the simulation whenever $`\varphi `$ is larger than $``$ 20. In this way $``$ 77 % of projection angles of bars result in b/p bulges (Table 4). Combes et al. (com90 (1990)) find in their simulations a similar percentage (85 %) for “favourable” projection angles to see b/p bulges in galaxies with an inclination between 80 and 90.

Since the frequencies of the simulated bulge types (Table 4, Col. 4) are derived only from a barred galaxy, these frequencies must be modified to compare them with our statistics of bulge types (Paper I) derived from all edge-on galaxies, with bar and without. Taking this into account the percentages for the bulge types of the simulated galaxy have to be multiplied with the observed fraction of barred galaxies. Adding afterwards the fraction of unbarred galaxies to the percentage of bulge type 4 — N-body simulations (Pfenniger & Friedli pfe (1991)) suggest that unbarred galaxies have such a bulge type — leads to the modified frequencies of the bulge types (Table 4, Col. 5 and 6). The errors of these frequencies depend on the uncertainties of the limits in the aspect angle of the bar ($`\varphi `$) (up to $`\pm 4^{}`$) and have values of $``$ 2.5 % for type 1 and 4 bulges and $``$ 5 % for type 2 and 3 bulges. The observed frequencies of the bulge types (Paper I; type 1: 4.1 %, type 2: 15.7 %, type 3: 25.2 %, and type 4: 55.0 %) are very similar (inside the errors) to the modified frequencies computed with a fraction of 55 % barred galaxies. This percentage is revealed by studies from optical surveys (e.g. RC3) (Ho et al. ho (1997), Knapen et al. kna2000 (2000), Paper I). However, the higher percentage of 70 % derived from NIR studies (Mulchaey & Regan mul (1997), Knapen et al. kna2000 (2000), Eskridge et al. esk (2000)) seems to be more reliable, since the NIR-light traces the dominant, old stellar population and is relatively unaffected by dust. Using this fraction for computing the modified frequencies, there is a difference of $``$ 10 % between observed and modified frequencies of galaxies with b/p bulges (type 1-3). There are less b/p bulges observed as computed from the investigated simulation coupled with NIR observations.

The existence of galaxies with non-b/p bulges, but with bars (Tab. 2) are consistent with the simulation where almost end-on orientated bars lead to elliptical bulges.

### 5.3 Quantitative parameters

Along the bar major axis of the simulated galaxy, measured in the face-on case ($`i=0^{}`$), the decrease of the surface density is exponential (Pfenniger & Friedli pfe (1991)) and the length of the bar can be determined to $``$ 20 kpc. This length exactly fits inside the radius of corotation which occurs at 10 kpc (Pfenniger & Friedli pfe (1991)). The axis ratio of the bar is $`b/a=0.7\pm 0.1`$ (Fig. 4, top right) and is inside the range of observed values (0.2 - 1) (Martin mar (1995)).

Regarding the galaxy edge-on ($`i=90^{}`$) with an edge-on orientated bar ($`\varphi =90^{}`$) is the best way to derive the quantitave bar parameters (introduced in section 3) of the simulated galaxy, because bar and b/p structures are most prominent at this geometry. However, the cuts along and parallel to the major axis do not reveal the typical bar signature (Fig. 7) which is visible in the NIR observations. The three components bulge, bar, and central part of the bulge cannot be distinguished as it is possible in the observations (cf. Fig. 1, 2, and 5). Due to the lack of a bar signature the bar thickness as well as the central bulge length (CBU) cannot be determined. The ratio of BAL/BUL $`=1.25\pm 0.2`$, taking into account BAL $`=20`$ kpc and BUL $`=16`$ kpc (Fig. 7), is not consistent with the observed ratios for peanut bulges ($`1.9\pm 0.3`$) (Fig. 4), while BAL/BPL $`=2.5\pm 0.2`$ (BPL $``$ 8 kpc; Fig. 7) is consistent with observed ratios in peanut bulges (BAL/BPL $`=2.7\pm 0.3`$). While in observations the positions of the maxima of the b/p distortion coincide with the limits of the central bulge, in the N-body simulation the positions of the maxima are not related to any other features in the profiles.

### 5.4 Relevance of the simulated galaxy for observations

The differences between the galaxy derived from the N-body simulation and the observed galaxies likely result from the initial conditions of the simulated model. Pfenniger & Friedli (pfe (1991)) are just using two merged disks with the same rotation properties instead of an independent bulge component in their simulation. The inner disk represents the “bulge” and is chosen with the main purpose “to have sufficient freedom to shape the initial rotation curve”. Therefore the initial resemblance with a real bulge is only approximate. Using such a bulge component in the simulation could also be the reason for the lack of the bar signature when the galaxy is seen edge-on, since the inner disk could cover the signature. Another explanations for the fact that no bar signature is visible could lay in the limited spatial resolution or in the lack of gas in the simulation, which influences the bar structure. It is also possible that the exponential decrease of the bar surface brightness and the smooth end of the bar prevent the detection of a bar signature.

Additionally, the investigated simulated galaxy is only a single, although representative, example of a barred galaxy. Variations of the initial conditions such as the bulge-to-disk ratio, Q parameter, gas fraction and distribution, and halo structure end up in different bar strenghts, shapes, thickness, lenghts, and density profiles in N-body simulations (e.g. Friedli & Benz friben (1993)). Therefore conclusions drawn from comparisons between our statistics of observed bulge types (Paper I) and results from the simulated galaxy, which is investigated in this study, are limited. Future studies have to analyse recently published simulations of barred galaxies (e.g. Pfenniger pfe2000 (1999)) to get more relevant comparisons with observations.

## 6 Discussion

### 6.1 Bars and b/p bulges

The strong correlation of the bar signature with b/p bulges in our study shows that b/p bulges are well explained by dynamical processes in bar potentials. The ratio of the bar length to the length of the b/p distortion (BPL) ($`2.5\pm 0.2`$) in the simulated barred galaxy and the observed ratio of BAL/BPL ($`2.7\pm 0.3`$) in galaxies with peanut bulges is consistent with such a scenario.

The correlation of the bar signature with b/p bulges is previously claimed, using the same detection method, by the investigations of individual galaxies, namely NGC 1381 (de Carvalho & da Costa deca (1987)), NGC 5170 (Dettmar & Barteldrees 1990b ), and NGC 4302 (Dettmar & Ferrara df (1996)). Many further evidences for the connection between bars and b/p bulges are reported in the literature: Bettoni & Galletta (bet (1994)) and Quillen et al. (qui (1997)) find them in their investigation of the intermediately inclined galaxies NGC 4442 resp. NGC 7582. Furthermore, N-body simulations (e.g. Combes et al. com90 (1990), Raha et al. rah (1991), Pfenniger & Friedli pfe (1991)) show that the evolution of a b/p bulge is coupled with the development of a bar.

Independent evidence is obtained by observing the gas kinematics in N II and H$`\alpha `$. The characteristic “figure-of-eight” rotation curve which is a strong signature for the presence of a bar and is detected in many b/p bulges: e.g. NGC 5746, NGC 5965 (Kuijken & Merrifield kui (1995)), NGC 2683 (Merrifield mer96 (1996)), and additional ten further galaxies by Bureau & Freeman (bf (1999)) and five by Merrifield & Kuijken (mer (1999)). Control samples of non-b/p bulges show no kinematical signature of a bar (Bureau & Freeman bf (1999), Merrifield & Kuijken mer (1999)). In two b/p bulges (ESO 597- 36 and NGC 3390) without kinematic bar signature Bureau & Freeman (bf (1999)) find indications for accretion events. However, since the positioning of the slit was difficult (Bureau & Freeman bf (1999)), dust in the galactic plane could affected their observations. Therefore bars can not be excluded in these two galaxies. A merger-induced bar formation with a resulting b/p bulge (Mihos et al. mih (1995)) could be a possible scenario for these galaxies. The four detected non-b/p bulges: NGC 5907 (Miller & Rubin mil (1995)), IC 5096, ESO 240- 11 (Bureau & Freeman bf (1999)), and NGC 681 (Merrifield & Kuijken mer (1999)) with a bar signature in the gas kinematics can be explained by our analysis of the N-body simulation with bars which are nearly end-on (see next subsection).

Signatures of the gas kinematics in barred potentials are also detected in retrograde orbits of NGC 128 (Emsellem & Arnault ems (1997)), the galaxy with one of the most prominent b/p bulges (D’Onofrio et al. dono (1999)). Direct kinematic evidence for bar streaming motions in two galaxies with b/p bulges (NGC 3079 and NGC 4388) is presented by Veilleux et al. (vei (1999)). Furthermore, in H I the presence of a bar is also verified in one galaxy with a peanut bulge analysing the position velocity diagram (IC 2531: Bureau & Freeman bur97 (1997)). The gas kinematics in the central region obtained in CO ($`\mathrm{J}=21`$ line) in the case of NGC 3628 (Reuter et al. reu (1991)), NGC 891 (Garcia-Burillo & Guelin gar (1995)), and NGC 4013 (Garcia-Burillo et al. gar99 (1999)) again point to the presence of a bar.

The conclusion (often given in the literature, e.g. Pfenniger & Friedli pfe (1991)), that b/p bulges are thick bars, cannot be supported. From the observational point of view bars and bulges, if b/p-shaped or not, are obviously individual components with different lengths (e.g. Fig. 1). Stellar orbits in a barred potential build up b/p structures and the bar itself. Therefore b/p bulges are the composition of a spheriodal bulge and a b/p structure representing the thick part of the bar.

The non-detection of bars (Table 3) in all observed TBBs ($``$ 5 % of all b/p bulges belong to this group, Paper III) is an indication for the hypothesis that the boxiness of these bulges is not related to the the presence of a bar. This idea gets more convincing support from the fact that the extent of the box shape is too large to result from an usual bar potential (Paper III). However, the lack of a bar signature does not necessarily mean that there is no bar, since weak bars, which have exponential surface brightness profiles (Elmegreen & Elmegreen 2elm89 (1989)), could remain undetecable by the here presented search for a bar signature.

In this way also the missing bar signature in the simulated barred galaxy is explainable, because the surface density distribution along the simulated bar is exponential and additionally, the bar axis ratio of this galaxy points to a weak bar ($`b/a>0.6`$, Martinet & Friedli marfri (1997)).

### 6.2 Explanations for various bulge shapes

The connection between b/p bulges and bars can be understood as interplay between resonant bending fed by vertical diffusion of orbits and instabilities similar to the fire-hose one (Bureau bur (1998), Friedli fri (1999)). While there is evidence for this connection in the most galaxies with b/p bulges, the question for the different shape of the b/p bulges remains to be answered.

A correlation between a strong signature of a bar, marked by constant surface brightness in radial direction (cf. section 3), and prominent b/p bulges is found in our study. Thereby the strong bar signature can be associated with a flat bar, which is known to be strong (Elmegreen & Elmegreen 2elm89 (1989)). Therefore the correlation can be interpreted by the assumption that strong barred galaxies are responsible for prominent b/p bulges, while galaxies with a weak bar have bulges with a lower degree of boxiness. N-body simulations by Friedli (fri (1999)) support this explanation for different bulge shapes. In this way it is also possible to explain partly the observed correlation of BAL/BUL with the bulge type (cf. section 4) by larger bars in galaxies with prominent b/p bulges, because flat bars tend to be larger (Elmegreen & Elmegreen 2elm95 (1995), Elmegreen et al. elm (1996)). Since this correlation holds for early as well as for late types, it cannot be explained by the fact that the same bar will produce more pronounced b/p bulges, the smaller the pre-existing spheroidal bulge is. However, having this in mind it is not surprising that the weak bar of the here investigated simulated galaxy can produce a peanut bulge. Furthermore, one could expect that peanut bulges more often occur in late type spirals, but the observed fraction of peanut bulges in galaxies in the range of S0-Scd is rather constant and equal to zero in spirals later than Sc (Paper I). This result could be interpreted by a dependence of the bar strength on the Hubble type. Such a correlation of bars is found by Elmegreen & Elmegreen (2elm (1985)) and Elmegreen et al. (elm (1996)). However, in contrast to them Seigar & James (sei (1998)) do not find that flat (strong) bars are more common in early Hubble types and exponential (weak) bars in later types. This statement is in agreement with our results for the profiles of bar signature concerning the Hubble type, but it is uncertain, if their results concerning bar morphology detected in face-on galaxies can be transformed to bar signatures in edge-on galaxies.

Additionally, for the explanation of bulge shapes it has to be taken into account that gas mass and a central black hole effect the shape of bulges. A significant fraction of gas mass (Berentzen et al. ber98 (1998)) or a supermassive central black hole can weaken the b/p structure (Friedli fri (1999)).

Explaining bulge shapes by aspect angles of bars was introduced by N-body simulations (Combes et al. com90 (1990), Pfenniger & Friedli pfe (1991)). Kinematical bar diagnostics in edge-on spiral galaxies using simulations of families of periodic orbits and hydrodynamical simulations are also consistent with observations. They confirm the dependence of the boxiness on the aspect angle of the bar (Bureau & Athanassoula ba (1999), Athanassoula & Bureau ab (1999)). An exact determination of aspect angles is not possible, because the envelope of the signature of the outer bar region in the position-velocity diagram is too faint (Athanassoula & Bureau ab (1999)).

The strong correlation of BAL/BUL with the bulge type and the agreement of the observed ratio BAL/BUL only of the most prominent b/p bulges, $`1.9\pm 0.3`$, with that ratio obtained for face-on galaxies, $`2.3\pm 0.7`$ (Athanassoula & Martinet ath (1980), Baumgart & Peterson baum (1986), Martin mar (1995)), shows that the differences in BAL/BUL can be interpreted by different aspect angles of bars. Therefore our analysis points also to a dependence of the shape of a bulge on the aspect angle of the bar. In this way and consistent with our investigation of the simulated barred galaxy peanut bulges result from bars seen nearly edge-on (large projected bar length, BAL), boxy bulges from intermediately inclined bars, and elliptical ones from end-on bars (small BAL) or they have no bar at all. Thereby bars detected in non-b/p bulges can be explained as well as differences between the fraction of barred galaxies (70 % in the NIR) and galaxies with b/p bulges (45 %), although also other parameters might play a role (weak bar, large gas mass, or dominant bulge component). The dependence of the bulge shape on the aspect angle of the bar is also supported in our study by the modified frequencies of the bulge types derived in the simulation. These frequencies are consistent with the observed frequencies of the bulge types, if one keeps in mind that only a single simulated barred galaxy is compared with all morphological types of disk galaxies. Larger bulge components in the simulation would weaken the b/p structures and thereby reduce the modified frequency of b/p bulges. Stronger bars would have the contrary effect on the modified frequency. Additionally, many galaxies of our investigated sample with peanut bulge and strong bar signature show features at the ends of the bumps, which can be associated with spiral arms extending outward from the ends of a bar (Baggett & MacKenty bag (1996)). Therefore it can be suggested that bars of these galaxies are indeed nearly edge-on.

There exists only a few rough estimates of aspect angles of bars in edge-on galaxies from observations. Wakamatsu & Hamabe (wak (1984)) adopt for NGC 4762 (bulge type 3 in the optical, type 4 in the NIR; Paper I) an aspect angle of the bar of $`\varphi =33^{}\pm 10^{}`$ (although a lower limit of 21 is also possible from their argumentation) and Garcia-Burillo & Guelin (gar (1995)) for NGC 891 (bulge type 2) $`\varphi 45^{}`$. These results are in agreement with a dependence of the bulge shape on the bar aspect angle.

### 6.3 A new correlation

The coincidence of the positions of the maxima of the b/p distortion with the limits of the central bulge, which is found here for the first time, cannot finally be explained ($`\mathrm{CBU}/\mathrm{BPL}=1.0\pm 0.1`$). Selection effects due to different aspect angles of the bar, which would influence the length of the b/p distortion (BPL), can be neglected, since we have shown that peanut bulges, for which BPL is defined only, can be associated with nearly edge-on bars. The coincidence points also to a possible connection between the bar length and the length of the central bulge (CBU), since BPL is connected with BAL.

Therefore possible explanations could be related to secondary bars (e.g. Friedli fri (1999), Greusard et al. gre (2000)). This idea is supported by the result from N-body simulations (Friedli & Martinet frimar (1993)) that the corotation of the inner bar roughly coincides with the location of the horizontal inner Lindblad resonance (ILR). At this position also the maximum of the elevation of the peanut structure occurs, because the vertical ILR is very close to the horizontal ILR and the conjunction of these two resonances build up the peanut structure (Combes et al. com90 (1990)). However, the observed ratio of the bar lengths of the primary to the secondary, which ends near corotation (Greusard et al. gre (2000)), varies between 4.0 and 13.4 (Friedli et al. fri96 (1996)), respectively 2.7 and 7.1 (Greusard et al. gre (2000)), while the ratio of bar length to length of the central bulge is much smaller and varies only between 2.4 and 3.1.

The central bulge could also be associated with a primordial bulge, which influences the bar structure and the position of the horizontal ILR. Since at this location also the peanut structures is most prominent, the extent of the b/p distortion would depend on the mass and length of the primordial bulge. Such a hypothesis is supported by N-body simulations in which the size of the evolving bar is directly related to the initial size of the bulge (Pfenniger pfe2000 (1999)). In this way the nearly constant ratio of bar length or BPL to CBU would be a natural result.

The correlation between CBU and BPL and especially the CBU/BPL ratio of 1 could also be explained by nuclear rings. They correspond to the horizontal ILR of bars (Buta & Crocker but (1993), Buta & Combes butcom (1996)) and therefore occur at the same radial position as the maximum of the elevation of the peanut structure. A nuclear ring could be responsible for the increase of the surface brightness, which we observe in radial cuts near the plane (Fig. 1, 2, 5), or the ring is only the border for the central bulge. However, we have not seen any structure inside the central bulge pointing to rings, but it is likely that our resolution is too low to detect them.

### 6.4 Bar thickness

Since we know now with high probability the fact that bars in peanut bulges are edge-on orientated, it is for the first time possible to determine the thickness of bars by observations. From the mean value of BAL/BAT for the galaxies with peanut bulges a ratio of bar length to thickness of $`14\pm 4`$ (taking the error in BAL and BAT into account) is derived. This value is at the upper limit of the values estimated in the literature (ranging from 4 to 15, Sellwood & Wilkinson sel (1993)). It must be mentioned that bars are more extended in the NIR than in the optical (Friedli et al. fri96 (1996), Chapelon et al. cha (1999)). However, from a sample of nine barred galaxies (Wozniak et al. woz (1995), Friedli et al. fri96 (1996)) investigated in the $`I`$\- and $`K`$-band we have computed that this effect is only of the order of $``$ 5%. Therefore our rather large value for the ratio is not an effect of the observed wavelengths.

## 7 Conclusion

We have analysed a sample of edge-on disk galaxies in the NIR and images derived from a N-body simulation of a barred galaxy using cuts parallel to the major axis and our classification of bulge shapes introduced in Paper I. Our main conclusions are:

* Cuts along and parallel to the major axis of edge-on galaxies show a strong correlation between bars and b/p bulges. This correlation is for the first time verified in a larger sample of galaxies (N $`=60`$).
* The investigation of the bar signatures in the NIR reveals a correlation between prominent b/p bulges and strong bar signatures pointing to a dependence of the boxiness of a galaxy on the bar strength.
* The ratio of the projected bar length (BAL) to the length of the bulge (BUL) correlates with the bulge types, i.e. with the bulge shape. This correlation is the first evidence from surface photometry for an explanation of the different degree of boxiness by the aspect angle of a bar. The maximum of this ratio for galaxies with peanut bulges is consistent with ratios observed in face-on galaxies.
* The ratio of BAL to the length of the b/p distortion (BPL) has a mean value for the observed peanut bulges of $`2.7\pm 0.3`$. This value is consistent with the investigated N-body simulation where a ratio of $`2.5\pm 0.2`$ can be derived for an edge-on galaxy with an edge-on bar.
* The frequencies of bulge types in the simulated galaxy, seen under different aspect angles of a bar and under high inclinations, are consistent with the observed frequencies corrected for the total fraction of barred galaxies. These results support the interpretation that bars viewed under different angles in edge-on galaxies play a major role for the bulge shape. However, the viewing angle is not the only parameter, since the size of the pre-existing spheroidal bulge and various bar shapes and types also influence the bulge shape.
* Interpreting galaxies with peanut bulges as galaxies with a bar seen edge-on the ratio of bar length to thickness, $`\mathrm{BAL}/\mathrm{BAT}=14\pm 4`$, can be directly measured for the first time.
* From the observational point of view bars and b/p bulges are different components revealed by their different lengths. The inner part of barred disk galaxies can be described by three components: the spheroidal bulge (including the central bulge), the bar (=thin bar), and the b/p structure (=thick bar).
* We have found a new size relation of the positions of the maxima of the b/p distortion with the extent of the central bulge ($`\mathrm{CBU}/\mathrm{BPL}=1\pm 0.1`$). This size relation is most likely related to the in-plane ILR of the bar.

###### Acknowledgements.

Part of this work was supported by the *Deutsche Forschungsgemeinschaft, DFG*. Furthermore the authors want to thank Daniel Pfenniger who provided the particle position file, Andreas Schröer for helping with IDL programming, and the referee F. Combes for constructive comments.


format: text

text:
# Many-body perturbation calculation of spherical nuclei with a separable monopole interaction: I. Finite Nuclei

## I Introduction

The central problem of nuclear structure theory is the solution of the many-body Schrödinger equation (MBSE). For Hamiltionans of interest in the nuclear case, an analytic solution is impossible, and one is compelled to use some approximation, either in the numerical solution of the equation or the specification of the Hamiltonian, or both.

Approaching the problem with the aim of using as realistic a representation of the potential as possible usually means fitting a combination of a meson exchange and phenomenological interaction to low-energy nucleon-nucleon scattering data and properties of few-body systems. To get good agreement with experiment both two- and three-body forces seem to be necessary. Recent examples of such potentials include the Bonn , the Argonne two-body  with Urbanna 3-body, Nijmegen  and Moscow  potentials, the last of which also incorporates quark degrees of freedom. These forces share the property of having a hard repulsive core which is a natural consequence of meson-exchange. It is this hard core which presents the difficulty in solving the MBSE. For instance, Hartree-Fock (HF) mean-field calculations with such interactions result in unbound nuclei. Treating corrections beyond the HF approximation order-by-order in perturbation theory is also unsuccessful since the interactions used are non-perturbative. One has to solve the full MBSE numerically in as exact a way as possible using techniques such as Variational Monte-Carlo, Green’s Function Monte Carlo, the coupled-cluster method, and the Fermion Hyper-netted chain model. Using effective interactions derived from realistic potentials, no-core shell-model calculations have been made in light nuclei and heavier nuclei close to closed shells have been treated.

The computational difficulty of performing numerically exact solutions of the MBSE has limited the techniques to light nuclei, for instance $`A=8`$ results have been published recently using the Argonne $`v18`$ and Urbanna IX potentials in the GFMC framework. In this work, it is seen that although the lightest nuclei are reproduced very well, the quantitative comparison of theory to data gets worse as $`A`$ increases. This may be due to the necessarily phenomenological nature of the three–body potential, a problem which may be overcome with re-fitting. On the other hand, it is not obvious that higher-body forces will not prove necessary or that the concept of a bare interaction between nucleons is valid for small distances.

Attempts were made in the late sixties primarily by the Kerman group at MIT to parameterize the NN interaction in such a way that it is weak in the sense of being perturbative. Such a weak interaction allows one to perform Hartree-Fock calculations to obtain a reasonable approximation to the full wavefunction and then to calculate corrections in perturbation theory. While this technique seems very attractive, the result obtained were only moderately successful at reproducing experimental data, a fact which was presumed to be due to inadequacies in the potentials used. The efficacy of developing a suitable interaction when similar, though more complicated, techniques were available for realistic interactions has been questioned and no better interaction was developed. Separable parameterizations, particularly the quadrupole-quadrupole interaction, have retained currency, as residual interactions . Even when the interactions are too strong for regular perturbation theory, separable interactions requiring solution of Brückner Hartree-Fock equations have proved fruitful because of their simplicity.

On the other hand, interactions have been developed which are not intended for use in the full MBSE, but rather to give good results with a Hartree-Fock calculation alone. Good quantitative success came with the zero-range density-dependent force of Ehlers and Moszkowski and Skyrme’s interaction, used in HF calculations by Vautherin and Brink  and subsequently by many others, and also Gogny’s finite-range interaction. Skyrme’s interaction has been particularly successful, in part due to its simple form, that of a delta function, which leads to easy calculation, even of the exchange part of the force. This computational simplicity has allowed extensive study of the properties of nuclei to be made with the Skyrme interaction across the entire range of nuclei in the periodic table . Related somewhat to the Skyrme-Hartree-Fock model is the Relativistic Mean Field (RMF) approach, which also gives single-particle motion in a mean field, but as a solution to the Dirac equation as opposed to the Schrödinger equation. The RMF approach has some nice features such as the natural occurrence of the spin-orbit splitting without recourse to an assumed spin-orbit interaction.

These mean-field models are inherently single-particle in nature. The forces used are not intended to be used in the MBSE, nor are explicit corrections beyond the mean-field part of the framework, although Skyrme’s interaction can be considered as a phenomenological G-matrix equivalent, and in that sense includes a subset of possible correlation effects in the mean-field. Although appropriate for use in mean-field calculations, Skyrme’s interaction would actually diverge in perturbation theory because of the zero range. This compels one to use a different interaction to obtain correlation beyond the mean field than was used to create it. Extra residual forces are used, such as pairing, or shell model interactions to allow for more general wavefunctions and obtain more accurate reproduction of the physics, or certain approximations are used such as RPA, which can describe certain observables, particularly those of giant resonance states and other forms of collective motion, but not others. It is thought that correlations should be particularly important in nuclei at the limits of stability, where the nucleons nearest the Fermi level couple strongly to the continuum.

We revisit the idea that it is possible to parameterize a nuclear interaction in such a way that it is weak enough with which to perform perturbation theory, thereby allowing correlated wavefunctions and observables to be calculated across the entire range of nuclei. Using the separable ansatz of previous “weak” interactions we have developed a density dependent interaction which we hope will provide some insight into the correlation structure of nuclear wavefunctions while retaining the quantitative power of contemporary effective interactions used in Hartree-Fock. In contrast to previous work, the interaction is designed to be an effective interaction with parameters fitted to the properties of finite nuclei within the calculation framework for which it is intended.

The paper is organized as follows. In Section 2 we briefly review basics of many-body perturbation theory. The separable interaction, used in the present work is given and discussed in Section 3. Results of the calculation for doubly magic nuclei are summarized in Section 4. Derivation of the HF energy and potential is outlined in Appendix A.

## II Single particle models and many-body perturbation theory

For standard perturbation theory to be successful, the Hamiltonian must be separated into two parts, one which is solvable ($`H_0`$), and another which is “small” ($`H_1`$). The Hamiltonian for a many-fermion problem in which the particles interact via one- and two-body interactions may be written schematically as

$$H=\underset{ij}{}i|U|ja_i^{}a_j+\frac{1}{4}\underset{ijkl}{}ij|V|kla_i^{}a_j^{}a_la_k.$$
(1)
Typically, the one-body part of the Hamiltonian is just the kinetic energy. This is a unsuitable choice for $`H_0`$ in the case of nuclei since the eigenstates of this operator, namely plane wave states, are not close enough to the exact eigenstates of the full Hamiltonian and $`H_1`$ is thus not small. By adding and subtracting a one-body term, $`U^{}`$, of ones choice, the $`H_0`$ part of the Hamiltonian may be solvable to give a wavefunction close to the solution of the full problem, thereby making the corrections from the $`H_1`$ part small enough for perturbation theory to succeed. Schematically the Hamiltonian is now split up into two terms;

$$H=\underset{H_0}{\underset{}{\underset{ij}{}i|U+U^{}|ja_i^{}a_j}}+\underset{H_1}{\underset{}{\frac{1}{4}\underset{ijkl}{}ij|V|kla_i^{}a_j^{}a_la_k\underset{ij}{}i|U^{}|ja_i^{}a_j}}.$$
(2)
One practical choice of $`U^{}`$ is a simple, analytically-solvable external potential, typically – in the nuclear case – of harmonic oscillator or Woods-Saxon form. In such a case, all of the physics of the nuclear interaction is in the residual part, $`H_1`$. Unfortunately, for sensible choices of two-body interaction, the residual part of the Hamiltonian is non-perturbative in the basis obtained from the solution of the one-body part of the Hamiltonian, and one has to perform more exact and exacting interaction shell-model calculations to arrive at a meaningful result. Traditionally, this has involved the diagonalisation of large matrices, or more recently, auxiliary field Monte-Carlo calculations in the SMMC.

Another popular choice for $`U^{}`$ is that of the Hartree-Fock mean-field. The HF mean-field potential is usually derived from the full Hamiltonian by a variational principle. Viewed in this way, it is the one-body potential whose occupied eigenstates form the lowest energy Slater-Determinant many-body wavefunction possible for the full many-body Hamiltonian. Using such a one-body potential, some physics of the two-body interaction is included in the single-particle problem. With a judicious choice of interaction, one ought to be able to produce very good results, since the approximation of the nucleus as a system of non-interacting particles in a mean-field is known to be a good one. Using the HF potential for $`U^{}`$ has the added attraction that it leads to vanishing first order corrections to the energy in perturbation theory. In this case, using Wick’s theorem, one can re-write (1) as

$$H=\underset{H_0}{\underset{}{E_0+\underset{i<ϵ_F}{}\epsilon _i:a_i^{}a_i:}}+\underset{H_1}{\underset{}{\frac{1}{4}\underset{ijkl}{}ij|V|kl:a_i^{}a_j^{}a_la_k:}}$$
(3)
where $`H_0`$ is the HF mean-field Hamiltonian and $`H_1`$ is the perturbing Hamiltonian. Note that $`H_1`$ is just the full two-body interaction matrix elements with a time-ordered product of creation and annihilation operators. The perturbation series for the energy is ordered by the number of matrix elements of the potential. Usually a diagrammatic representation is used. One can write down the number of diagrams for any particular order of perturbation theory. We have evaluated these for the vacuum amplitude up to order seven. The results are presented in Table I. One sees that, for the method of direct evaluation of diagrams by order to be effective, the series must be sufficiently converged by fourth, or perhaps fifth, order.

In this work, the second and third order diagrams for the vacuum amplitude, as shown in Figs. 1 and 2, are calculated. Their algebraic form is given here as

$`E_2`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{ab<ϵ_F}{}}{\displaystyle \underset{rs>ϵ_F}{}}{\displaystyle \frac{|ab|\stackrel{~}{V}|rs|^2}{ϵ_a+ϵ_bϵ_rϵ_s}}`$ (4)
$`E_3^{(pp)}`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \underset{ab<ϵ_F}{}}{\displaystyle \underset{rs>ϵ_F}{}}{\displaystyle \underset{tu>ϵ_F}{}}{\displaystyle \frac{ab|\stackrel{~}{V}|rsrs|\stackrel{~}{V}|tutu|\stackrel{~}{V}|ab}{(ϵ_a+ϵ_bϵ_rϵ_s)(ϵ_a+ϵ_bϵ_tϵ_u)}}`$ (5)
$`E_3^{(hh)}`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \underset{ab<ϵ_F}{}}{\displaystyle \underset{cd<ϵ_F}{}}{\displaystyle \underset{rs>ϵ_F}{}}{\displaystyle \frac{ab|\stackrel{~}{V}|rscd|\stackrel{~}{V}|abrs|\stackrel{~}{V}|cd}{(ϵ_a+ϵ_bϵ_rϵ_s)(ϵ_c+ϵ_dϵ_tϵ_u)}}`$ (6)
$`E_3^{(ph)}`$ $`=`$ $`{\displaystyle \underset{abc<ϵ_F}{}}{\displaystyle \underset{rst>ϵ_F}{}}{\displaystyle \frac{ab|\stackrel{~}{V}|rscr|\stackrel{~}{V}|atst|\stackrel{~}{V}|cb}{(ϵ_a+ϵ_bϵ_rϵ_s)(ϵ_b+ϵ_cϵ_sϵ_t)}}`$ (7)
in which the tildes over the potential indicate the matrix element is antisymmetrized. The state vectors label HF single-particle states, whose energies are given by the subscripted $`ϵ`$.

It is important to note that our interaction is not intended to fit scattering data, having, as it does, density dependence. On the level of the perturbation theory it is necessary to treat the density functions as just the spatial form of the interaction, rather than a representation of a many-body force. This is to be considered a part of the present model. To do otherwise would be to surrender the simplifications our weak, separable potential affords.

In the present work, the Hartree-Fock problem is solved in a basis of spherical harmonic oscillator states. This yields, along with the hole states, a large number of particle states, with which to directly evaluate the sums of the perturbation series. A sufficient number of states is used so that the particle states are oscillatory over the size of the nucleus and that both the HF solution and the perturbation corrections are reasonably converged.

## III Interaction

We have developed an interaction written in the form of a sum of separable terms, which is to say it is in the form $`V(r_1,r_2)g(r_1)g(r_2)`$. The functions $`g`$ carry no angular momentum ($`l=0`$), and the force is dubbed a monopole-monopole interaction. For future applications, it is in intended to include higher multipole forces, with $`l=1,2,\mathrm{}`$, within our framework as these will presumably be necessary for calculation of excited states and deformed nuclei. Although higher multipole forces will contribute to spherical nuclei from the exchange term in Hartree Fock order and via correlations in perturbation theory, they are not included in the present calculation since it seems unwise to attempt to fit the parameters of such forces to spherical nuclei alone.

In coordinate space, the monopole interaction is written as

$`V(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ $`=`$ $`W_af_a\rho ^{\beta _a}(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_2)(1+a_a(\tau _1^+\tau _2^{}+\tau _1^{}\tau _2^+)+b_at_{1z}t_{2z})`$ (8)
$`+`$ $`W_rf_r\rho ^{\beta _r}(\stackrel{}{r}_1)\rho ^{\beta _r}(\stackrel{}{r}_2)(1+a_r(\tau _1^+\tau _2^{}+\tau _1^{}\tau _2^+)+b_rt_{1z}t_{2z})`$ (9)
$`+`$ $`k_1^2\rho (\stackrel{}{r}_1)_2^2\rho (\stackrel{}{r}_2),`$ (10)
where the function $`f_\xi `$ is defined as

$$f_\xi =\left[\mathrm{d}^3\stackrel{}{r}\rho ^{\alpha _\xi }(\stackrel{}{r})\right]^1,$$
(11)
for subscripts $`\xi =a`$ and $`\xi =r`$. Throughout this work, the three terms in (10) are referred to, in the order they appear in the above expression, as the attractive, repulsive and derivative terms.

In addition, the spin-orbit force is taken to be

$$V_{\mathrm{s}\mathrm{o}}(r)=c\frac{1}{r}\frac{\rho }{r}\stackrel{}{l}\stackrel{}{s},$$
(12)
which is similar to that used in the modified delta interaction .

The parameters $`W_a`$, $`\alpha _a`$, $`\beta _a`$, $`a_a`$, $`b_a`$, $`W_r`$, $`\alpha _r`$, $`\beta _r`$, $`a_r`$, $`b_r`$, $`k`$ and $`c`$ are to be fitted to experimental data.

One notices that the two-body interaction consists of a sum of terms, each of which is separable in form and that the expressions for the attractive and repulsive terms in (10) differ only by the values of their parameters

The energy, $`E_{\mathrm{pot}}`$, due to the interaction (10) in the Hartree-Fock approximation is derived in Appendix A (A4), (A10), (A41) and is presented here;

$`E_{HF}`$ $`=`$ $`T+E_{\mathrm{coul}}+E_{\mathrm{pot}}=T+E_{\mathrm{coul}}+{\displaystyle \underset{\xi =a,r}{}}\{{\displaystyle \frac{1}{2}}W_\xi f_\xi N_\xi ^2{\displaystyle \frac{1}{2}}W_\xi f_\xi M_\xi +{\displaystyle \frac{1}{2}}W_\xi b_\xi f_\xi (\mathrm{\Delta }N_\xi )^2`$ (14)
$`{\displaystyle \frac{1}{2}}W_\xi f_\xi [b_\xi M_\xi +a_\xi M_\xi ^{(\tau \overline{\tau })}]\}+{\displaystyle \frac{1}{2}}kN_d^2+cN_w`$
where $`T`$ is the kinetic energy, $`E_{\mathrm{coul}}`$ is the direct Coulomb energy plus exchange in the Slater approximation. The following quantities have been defined;

$`N_\xi `$ $`=`$ $`{\displaystyle \mathrm{d}^3\stackrel{}{r}\rho ^{\beta _\xi +1}(\stackrel{}{r})}`$ (15)
$`M_\xi `$ $`=`$ $`{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\left[\rho _p(\stackrel{}{r}_1,\stackrel{}{r}_2)\rho ^{\beta _\xi }(\stackrel{}{r}_1)\rho ^{\beta _\xi }(\stackrel{}{r}_2)\rho _p(\stackrel{}{r}_1,\stackrel{}{r}_2)+\rho _n(\stackrel{}{r}_1,\stackrel{}{r}_2)\rho ^{\beta _\xi }(\stackrel{}{r}_1)\rho ^{\beta _\xi }(\stackrel{}{r}_2)\rho _n(\stackrel{}{r}_1,\stackrel{}{r}_2)\right]}`$ (16)
$`\mathrm{\Delta }N_\xi `$ $`=`$ $`{\displaystyle \mathrm{d}^3\stackrel{}{r}\rho ^{\beta _\xi }(\stackrel{}{r})\left[\rho _p(\stackrel{}{r})\rho _n(\stackrel{}{r})\right]}`$ (17)
$`M_\xi ^{(\tau \overline{\tau })}`$ $`=`$ $`{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\left[\rho _p(\stackrel{}{r}_1,\stackrel{}{r}_2)\rho ^{\beta _\xi }(\stackrel{}{r}_1)\rho ^{\beta _\xi }(\stackrel{}{r}_2)\rho _n(\stackrel{}{r}_1,\stackrel{}{r}_2)+\rho _n(\stackrel{}{r}_1,\stackrel{}{r}_2)\rho ^{\beta _\xi }(\stackrel{}{r}_1)\rho ^{\beta _\xi }(\stackrel{}{r}_2)\rho _p(\stackrel{}{r}_1,\stackrel{}{r}_2)\right]}`$ (18)
$`N_d`$ $`=`$ $`{\displaystyle \mathrm{d}^3\stackrel{}{r}\rho (\stackrel{}{r})^2\rho (\stackrel{}{r})}`$ (19)
$`N_w`$ $`=`$ $`{\displaystyle \mathrm{d}^3\stackrel{}{r}\frac{1}{r}\frac{\rho }{r}\rho _w(\stackrel{}{r})}`$ (20)
and the following densities are used:

$`\rho (\stackrel{}{r})`$ $`=`$ $`\rho _p(\stackrel{}{r})+\rho _n(\stackrel{}{r})={\displaystyle \underset{i<ϵ_Fp}{}}\phi _i^{}(\stackrel{}{r})\phi _i(\stackrel{}{r})+{\displaystyle \underset{i<ϵ_Fn}{}}\phi _i^{}(\stackrel{}{r})\phi _i(\stackrel{}{r})`$ (21)
$`\rho (\stackrel{}{r}_1,\stackrel{}{r}_2)`$ $`=`$ $`\rho _p(\stackrel{}{r}_1,\stackrel{}{r}_2)+\rho _n(\stackrel{}{r}_1,\stackrel{}{r}_2)={\displaystyle \underset{i<ϵ_Fp}{}}\phi _i^{}(\stackrel{}{r}_1)\phi _i(\stackrel{}{r}_2)+{\displaystyle \underset{i<ϵ_Fn}{}}\phi _i^{}(\stackrel{}{r}_1)\phi _i(\stackrel{}{r}_2)`$ (22)
$`\rho _w(\stackrel{}{r})`$ $`=`$ $`{\displaystyle \underset{i<ϵ_F}{}}{\displaystyle \frac{1}{2}}\left(j_i(j_i+1)l_i(l_i+1)3/4\right)\phi _i^{}(\stackrel{}{r})\phi _i(\stackrel{}{r}).`$ (23)
The variation of the total energy is carried out in Appendix A (see A28, A31, A38, A45) . The resulting local Hartree-Fock potential is

$`U_{L,\tau }(x)`$ $`=`$ $`{\displaystyle \underset{\xi =a,r}{}}\{W_\xi f_\xi [N_\xi (\beta _\xi +1)+b_\xi \mathrm{\Delta }_\tau N_\xi ]\rho ^{\beta _\xi }(x)`$ (24)
$``$ $`W_\xi (\alpha _\xi /2)f_\xi ^2\left[N_\xi ^2+b_\xi (\mathrm{\Delta }N_\xi )^2(1+b_\xi )M_\xi a_\xi M_\xi ^{(\tau \overline{\tau })}\right]\rho ^{\alpha _\xi 1}(x)`$ (25)
$``$ $`W_\xi f_\xi \beta _\xi \left[(1+b_\xi )G_\xi (x)+a_\xi G_\xi ^{(\tau \overline{\tau })}(x)\right]\rho ^{\beta _x1}(x)`$ (26)
$`+`$ $`\left[W_\xi b_\xi \beta _\xi f_\xi \mathrm{\Delta }N_\xi \right]\rho ^{\beta _x1}(x)\delta \rho (x)\}`$ (27)
$`+`$ $`2kN_d^2\rho (x)`$ (28)
which differs for protons ($`\tau =p`$) and neutrons ($`\tau =n`$) through the function

$$\mathrm{\Delta }_\tau N_\xi =\{\begin{array}{cc}\mathrm{\Delta }N_\xi ,& \tau =p\hfill \\ \mathrm{\Delta }N_\xi ,& \tau =n.\hfill \end{array}$$
(29)
The other newly-introduced functions in (28) are

$`G_\xi (x)`$ $`=`$ $`G_\xi ^{(pp)}(x)+G_\xi ^{(nn)}(x)={\displaystyle \mathrm{d}^3\stackrel{}{r}\left[\rho _p(\stackrel{}{r},\stackrel{}{x})\rho ^{\beta _\xi }(r)\rho _p(\stackrel{}{x},\stackrel{}{r})+\rho _n(\stackrel{}{r},\stackrel{}{x})\rho ^{\beta _\xi }(\stackrel{}{r})\rho _n(\stackrel{}{x},\stackrel{}{r})\right]}`$ (30)
$`G_\xi ^{(\tau \overline{\tau })}(x)`$ $`=`$ $`G_\xi ^{(pn)}(x)+G_\xi ^{(np)}(x)={\displaystyle \mathrm{d}^3\stackrel{}{r}\left[\rho _p(\stackrel{}{r},\stackrel{}{x})\rho ^{\beta _\xi }(r)\rho _n(\stackrel{}{x},\stackrel{}{r})+\rho _n(\stackrel{}{r},\stackrel{}{x})\rho ^{\beta _\xi }(\stackrel{}{r})\rho _p(\stackrel{}{x},\stackrel{}{r})\right]}.`$ (31)
In addition, the non-local component to the mean field is (see A39)

$$U_{NL,\tau }(\stackrel{}{x},\stackrel{}{x}^{})=\underset{\xi =a,r}{}W_xf_x\rho ^{\beta _\xi }(\stackrel{}{x})\rho ^{\beta _\xi }(\stackrel{}{x}^{})\left\{(1+b_\xi )\rho _\tau (\stackrel{}{x},\stackrel{}{x}^{})+a_\xi \rho _{\overline{\tau }}(\stackrel{}{x},\stackrel{}{x}^{})\right\}.$$
(32)
and there is a state-dependent potential from the spin-orbit interaction of the form

$$U_{\mathrm{so}}(\stackrel{}{x})\phi _b(x)=c\left(w_b\frac{1}{x}\frac{\rho }{x}\frac{1}{x}\frac{\rho _w}{x}\frac{1}{x^2}\rho _w(x)\right)\phi _b(x)$$
(33)
where $`w_b=1/2(j_b(j_b+1)l_b(l_b+1)3/4)`$ is the spin-orbit weight factor.

Note that the one-body spin-orbit term could be taken as either a one-body force, or as a one-body potential deriving from a two-body force. Since the latter approach would render the perturbation calculation problematic due to the absence of a suitable of form of the two-body force, we choose the former approach. Hence, since the force is density-dependent, we have also included the rearrangement contribution to the HF potential. Only the non-rearrangement term actually gives rise to the spin-orbit splittings, but the rearrangement terms, coming as they do from a variational principle, result in a lowering of the HF energy. Combining the potentials (28), (32) and (33) gives us the HF equation

$$U_{L,\tau }(\stackrel{}{x})\phi _b(\stackrel{}{x})+\mathrm{d}^3\stackrel{}{x}^{}U_{NL,\tau }(\stackrel{}{x},\stackrel{}{x}^{})\phi _b(\stackrel{}{x}^{})+U_{\mathrm{so}}(\stackrel{}{x})\phi _b(\stackrel{}{x})=\epsilon _b\phi _b(\stackrel{}{x})$$
(34)
In this potential, as well as in the expression for the total energy (14), the exchange contribution from the derivative term is omitted. While it would, in principle, be desirable to include this term, the calculational complexity involved in doing so has forced the omission in the present case. However, for the main attractive and repulsive terms, the exchange part is much smaller than the direct in all nuclei, and the direct derivative term gives a rather small contribution to the mean-field and the binding energy in comparison to the other direct terms, so it is not considered an unwarranted approximation to neglect the effects of this term.

It might be objected that the form of the interaction is too unrealistic. For instance, since it is separable it cannot satisfy Galilean invariance. Furthermore, the unusual form of the isospin operators gives rise to a different effective interaction for protons and neutrons. As for the second point, since this is a density-dependent interaction, the effect of the difference between neutrons and protons comes from the densities as well as the isospin operators. This being the case, it seems reasonable to allow for the “stretching” of the isospin operator as we have done, to allow greater freedom in fitting the interaction to data. As for the first issue, the unusual form has been motivated by the desire to keep the force “weak” and is justified by the quality of the results produced. A possible method for restoring Galilean invariance within our framework is discussed in Section IV.

The choice of omitting a spin-spin ($`\sigma \sigma `$) type force yet having an isospin-isospin type force is motivated by the nuclei under study. All the closed-shell nuclei are spin-saturated and would contribute only through the exchange term in the HF order. For this separable interaction, the space-exchange terms are rather small and a spin-spin force would add little to the results. In addition, even if the effects in closed-shell nuclei are important, it does not seem reasonable to fit this term to closed-shell nuclei alone. It remains an open question whether such a force will prove necessary or useful in open-shell nuclei.

It is interesting to compare the leading terms in the HF mean field to that of other models. The first line of equation (28) gives us this as

$$U(x)\underset{\xi =a,r}{}C_\xi (f_\xi N_\xi )\rho ^{\beta _\xi }(x)$$
(35)
where $`C_\xi `$ is a combination of constants. The product $`f_\xi N_\xi `$ is

$$f_\xi N_\xi =\frac{\mathrm{d}^3\stackrel{}{r}\rho ^{\beta _\xi +1}(\stackrel{}{r})}{\mathrm{d}^3\stackrel{}{r}\rho ^{\alpha _\xi }(\stackrel{}{r})}.$$
(36)
If $`\beta _\xi +1=\alpha _\xi `$ then the product $`f_\xi N_\xi `$ is constant and the leading mean field terms go like

$$U(x)𝒞_a\rho ^{\beta _a}(\stackrel{}{x})+𝒞_r\rho ^{\beta _r}(\stackrel{}{x})$$
(37)
which, for the special case $`\beta _a=1`$, are the same as the terms in the Skyrme and Gogny mean-field proportional to the parameters $`t_0`$ and $`t_3`$, which give the bulk of the binding energy and saturation properties. In this work, we do not strictly keep $`\beta _\xi +1=\alpha _\xi `$, thus allowing for some A-dependence of the coefficients in the mean-field potential. It has been found, however, that one can not go too far away from the equality and still obtain reasonable results.

That one can get similar results in a mean-field calculation from two very different interactions is reflected in the different constitution of the resulting perturbative part of the Hamiltonian, $`H_1`$.

## IV Doubly (Semi-) Magic Nuclei

In order to find the best set of parameters for the interaction (10), calculations have been made of 14 doubly closed-shell nuclei across the periodic table. They are <sup>16</sup>O, <sup>34</sup>Si, <sup>40,48</sup>Ca, <sup>48,56,68,78</sup>Ni, <sup>90</sup>Zr, <sup>100,114,132</sup>Sn, <sup>146</sup>Gd and <sup>208</sup>Pb. The nuclei represent a selection of doubly-closed (sub-)shell nuclei both close to and far from stability. There is limited experimental information about <sup>48</sup>Ni  and <sup>100</sup>Sn . <sup>78</sup>Ni has yet to be discovered.

The ability to reproduce the properties of such exotic nuclei will be important for applications of our technique and discrepancies will help direct refinements.

A Hartree-Fock code assuming spherical symmetry and representing wavefunctions in a basis of spherical harmonic oscillator states was used to calculate uncorrelated wavefunctions. Perturbation corrections to the binding energy were directly evaluated using the results of the HF calculation. The results presented here were obtained in a basis of 12 expansion coefficients per single-particle wavefunction and iterated until the HF energy had converged to within 1keV. The parameters of the force were fitted to binding energies to second order and charge radii, charge density distributions, single-particle energies and spin-orbit splittings to HF order of the nuclei listed above where experimental data were available, and are presented in Table II.

The results of the the calculated energies, in HF order and in each order of perturbation theory are presented in Table III. The differences between the HF energy $`E_{\mathrm{cal}}=(\text{14})E_{\mathrm{HF}}`$ and the experimental ground state energy $`E_{\mathrm{exp}}`$ and between HF plus second order perturbation correction (4) $`E_{\mathrm{cal}}=E_{\mathrm{HF}}+E_2`$ and experiment are shown in Figure 3. The experimental energies are taken from the mass table of Audi and Wapstra  with two exceptions. An estimate of the mass of the recently-discovered nuclei <sup>48</sup>Ni  and the measured mass of <sup>100</sup>Sn . The energy for <sup>78</sup>Ni was taken from  in which extrapolated values are given, which are thought to be in error by less than $`0.2\%`$.

One sees from Figure 3 that most of the nuclei fit the binding energy to within $``$2%. The most obvious exception is <sup>16</sup>O which is quite under-bound. This may be due to the omission of a center-of-mass correction which would undoubtedly go a long way to close the discrepancy in the energy. It was not calculated in this case since a rigorous microscopic correction would destroy the mean-field which provides the essential basis for the perturbation calculation. A phenomenological correction could have been calculated, but perhaps the most suitable method in our framework would be to include with our multipole forces an isoscalar dipole force which could be fitted to restore the translational invariance of the many-body Hamiltonian, and evaluated exactly in perturbation theory.

A general trend can be seen in which lighter nuclei are somewhat over-bound and the heaviest are under-bound. It is the exceptions which conspire to stop the fitting algorithm from doing better, but the somewhat systematic nature of this discrepancy suggests that a better mass or isospin dependence may improve matters. It is unclear as yet the extent to which multipole correlations or a spin-spin force would improve the fit to spherical nuclei. That question awaits the study of deformed nuclei and excited states.

In Table IV a comparison is made of the quality of the fit to the binding energy to properties of the same nuclei calculated with a selection of Skyrme parameterizations. The parameterizations used are SIII, SkP, SLy4 and SkI4. In this comparison, it is seen that the energies from the different Skyrme parameterizations are of a similar quality, all reproducing the binding energies of closed-shell nuclei very well, with only a few binding energies being reproduced no better than 1% – including <sup>48</sup>Ni whose experimental value is in any case not well known. It is clear that the results from the separable force are somewhat worse. Particularly problematic is <sup>16</sup>O, whose large under-binding was mentioned above, and also <sup>48</sup>Ni which is, as with the Skyrme parameterizations, over-bound, although more so with the separable interaction.

Results for one-body properties are also presented. Comparison of the charge density results to experiment  and to the selection of Skyrme interactions is made in Table VI. One-body observables are generally reproduced better in the HF calculation alone than the binding energies. The comparison of the calculated charge radii with experiment is generally more favourable than the energy data. The radius of oxygen is too large by about 5% which is consistent with its under-binding. It can be seen that the agreement with experiment is of the same level as the Skyrme interactions. Perturbative corrections to the one-body observables, such as the densities, and hence radii, will be calculated in future work.

Figures 5-9 show the electron scattering form-factors of a selection of the nuclei compared to experiment . The proton density was corrected for the finite proton size by folding with a Gaussian to give the charge density, from which the radii and form-factors were calculated. The form-factors agree with experiment quite well, which is expected given the generally correct radii.

Some typical single-particle energies for light and heavy nuclei are shown in Figures 10 and 11 for <sup>40</sup>Ca and <sup>208</sup>Pb respectively. The single-particle energies of a density-dependent Hartree-Fock calculation do not directly correspond to an experimental observable, so caution should be used in comparing values. It can be seen that the level spacings and shell closures are better reproduced in <sup>208</sup>Pb, which is true of heavy nuclei in general. The comparatively poorer results in light nuclei seems to be common to mean-field approaches . In the case of <sup>40</sup>Ca, the gap at the fermi level may be widened by the inclusion of multipole forces which will link the occupied levels with the $`f_{7/2}`$ states, which is otherwise “inert”.

Table V shows spin-orbit splittings for some cases where the experimental values are known. The “experimental” data presented represents that used in previous work for fitting effective interactions to data . Clearly the splittings are all systematically small. This could be remedied by an increase in the spin-orbit coefficient, $`c`$. In a previous work , a value 10% higher than ours was used for the same spin-orbit interaction, and hence the spin-orbit splittings were more realistic. The lower value used in our work is the result of a compromise between the reproduction of the spin-orbit splittings and the total binding energies. This is a further indication that a more suitable spin-orbit potential needs to be sought.

The perturbation correction to the energy are seen to be rather small in all nuclei considered. This is consistent with our goal that the mean field solution should be close to the exact solution of the MBSE. The size of the second order correlation is roughly constant across the periodic table. It is characterized by a dimensionless strength parameter, $`\kappa `$, defined as

$$\kappa =\frac{1}{4}\underset{abϵ_F}{}\underset{rs>ϵ_F}{}\frac{ab|\stackrel{~}{V}|rsrs|\stackrel{~}{V}|ab}{\left(ϵ_a+ϵ_bϵ_rϵ_s\right)^2}.$$
(38)
It is related to the “wound integral” and is proportional to the number of 2p2h states excited due to the second order perturbation.

Figure 4 shows the correlation structure from the second order correction in the $`N=Z`$ nucleus <sup>40</sup>Ca, and the nuclei <sup>48</sup>Ca and <sup>208</sup>Pb. The contribution to the second-order energy is defined as a function of one of the particle states $`r`$;

$$E_2(|r)=\frac{1}{4}\underset{abϵ_F}{}\underset{s>ϵ_F}{}\frac{\left|ab|\stackrel{~}{V}|rs\right|^2}{ϵ_a+ϵ_bϵ_rϵ_s}.$$
(39)
The plot shows the contribution to the total second order energy correction as a function of the single particle energy $`ϵ_r`$, in 5 Mev wide bins. In all three cases particles are dominantly excited to low-lying states above the Fermi level. This results in a ground state with occupation probabilities similar to those which result from pairing forces. It is also a further indication that perturbation theory makes sense for our interaction since it does not predict excitation of particles in the ground state to extremely high energies.

Since the calculations were made only using a monopole force, the correlation structure is not expected to be complete. Only corrections involving simultaneous $`l=0`$ scattering of two particles is included. An indication of this is seen in the difference between the results for $`N=Z`$ and $`NZ`$ nuclei. The second order correction in <sup>48</sup>Ca is much larger than that in <sup>40</sup>Ca, due to the possibility of an $`f_{7/2}`$ neutron exciting to the $`f_{7/2}`$ proton state while another proton excites to a neutron state. This extra excitation is the labelled peak in Fig. 4. As well as having large wavefunction overlaps, the energy denominator in this case is much smaller than in any other possible monopole excitation which must excite any particle across major shells to keep all angular and isospin quantum numbers the same. When general excitations are permitted by higher multipole forces $`l=1,2\mathrm{}`$, this difference between correlations in $`N=Z`$ and $`NZ`$ nuclei will be smoothed out. For this reason, too, the correlation energies should not be considered too quantitatively at this stage, but rather as an indication of the perturbative properties of the interaction.

Comparing the form of the interaction to that of Skyrme suggests other possible sources of improvement to the model. One such may come from a better parameterization of the spin-orbit interaction. A two-body form which fits the philosophy of the separable effective interaction has not been found, but may be necessary to give the correct contribution to the binding energy. In any case, the simple form used in the present work which depends on a radial derivative will need modification if it is to be applied to deformed nuclei. It may also prove fruitful to explore a more general term dependent upon the derivatives of the density than the single term with parameter $`k`$, such as is found in the Skyrme interaction with two terms proportional to $`t_1`$ and $`t_2`$, which often carry further exchange parameters $`x_1`$ and $`x_2`$.

## V Conclusion

We have presented a new effective nuclear interaction which is designed for use in calculations which go beyond the mean-field. The technique of using perturbation theory to build correlations on top of the Hartree-Fock result is applicable to our interaction and results in small corrections to the single-particle behavior. A monopole-monopole force alone gives reasonable results for the ground state properties of spherical doubly-magic nuclei. It is expected that the addition of multipole forces will improve these results, particularly through the completion of the correlation structure. Such multipole forces will also presumably be important in giving the correct shapes of deformed nuclei, which are the subject of a forthcoming study and in the correct description of excited states.

## VI Acknowledgements

The authors would like to acknowledge useful discussions with D. Vautherin, P.-G. Reinhard, D. M. Brink and D. J. Dean. This research was sponsored by the Division of Nuclear Physics, U.S. Dept. of Energy under contract DE-AC05-00OR 22725 managed by UT–Battelle, LLC, and by research the UK EPSRC, and by US DOE grant no. DE-FG02-94ER40834.

## A HF Energy and Potential

The interaction is given in equation (10). Its expectation value is the contribution it makes to the total energy and is

$$E_{\mathrm{pot}}=\frac{1}{2}\underset{ij<ϵ_F}{}ij|V\left\{|ij|ji\right\}.$$
(A1)
If we consider just the attractive term - that is, the term whose parameters have the subscript “a” - then we will obtain the contribution from the repulsive term by simply substituting the subscript “a” for “r”.

$`E_\mathrm{a}`$ $`=`$ $`{\displaystyle \frac{1}{2}}W_af_a{\displaystyle \underset{ij<ϵ_F}{}}ij|\rho ^{\beta _a}(r_1)\rho ^{\beta _a}(r_2)\left\{|ij|ji\right\}`$ (A2)
$`+`$ $`{\displaystyle \frac{1}{2}}W_af_aa_a{\displaystyle \underset{ij<ϵ_F}{}}ij|\rho ^{\beta _a}(r_1)\rho ^{\beta _a}(r_2)(\tau _1^+\tau _2^{}+\tau _1^{}\tau _2^+)\left\{|ij|ji\right\}`$ (A3)
$`+`$ $`{\displaystyle \frac{1}{2}}W_af_ab_a{\displaystyle \underset{ij<ϵ_F}{}}ij|\rho ^{\beta _a}(r_1)\rho ^{\beta _a}(r_2)t_{1z}t_{2z}\left\{|ij|ji\right\}.`$ (A4)
Taking the first line, the matrix element is represented in space (and spin and isospin) coordinates;

$`{\displaystyle \frac{1}{2}}W_af_a{\displaystyle \underset{ij<ϵ_F}{}}{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\varphi _i^{}(\stackrel{}{r}_1)\varphi _j^{}(\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_2)\varphi _i(\stackrel{}{r}_1)\varphi _j(\stackrel{}{r}_2)}`$ (A6)
$`{\displaystyle \frac{1}{2}}W_af_a{\displaystyle \underset{ij<ϵ_F}{}}{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\varphi _i^{}(\stackrel{}{r}_1)\varphi _j^{}(\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_2)\varphi _j(\stackrel{}{r}_1)\varphi _i(\stackrel{}{r}_2)}`$
$`=`$ $`{\displaystyle \frac{1}{2}}W_af_a{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\rho ^{\beta _a+1}(\stackrel{}{r}_1)\mathrm{d}^3\stackrel{}{r}_2\rho ^{\beta _a+1}(\stackrel{}{r}_2)}`$ (A9)
$`{\displaystyle \frac{1}{2}}W_af_a{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\rho _p(\stackrel{}{r}_1,\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_2)\rho _p(\stackrel{}{r}_2,\stackrel{}{r}_1)}`$
$`{\displaystyle \frac{1}{2}}W_af_a{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\rho _n(\stackrel{}{r}_1,\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_2)\rho _n(\stackrel{}{r}_2,\stackrel{}{r}_1)}`$
$`=`$ $`{\displaystyle \frac{1}{2}}W_af_aN_a^2{\displaystyle \frac{1}{2}}W_af_aM_a`$ (A10)
where quantities defined in section III are used. Note that integrals include sums over spinors and isospinors where appropriate, and the coordinates include spin and isospin coordinates where appropriate. Where densities are used, the summing over isospin states has already been done and where densities do not carry isospin labels, the isoscalar density is assumed. See equations (21) and (22) for definitions.

The second line of (A4) contains “isospin-flipping” operators whose action is to turn an isospin state where particle one is a proton and particle two a neutron, $`|pn`$, into $`|np`$ and vice-versa. The direct contribution, in which the labels in the bra and the ket are in the same order is zero since all proton states are orthogonal to all neutron states. The exchange term is similar to that in (A10) but with different isospin combination of the density matrices;

$`{\displaystyle \frac{1}{2}}W_af_aa_a{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\rho _p(\stackrel{}{r}_1,\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_2)\rho _n(\stackrel{}{r}_2,\stackrel{}{r}_1)}`$ (A12)
$`{\displaystyle \frac{1}{2}}W_af_aa_a{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\rho _n(\stackrel{}{r}_1,\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_2)\rho _p(\stackrel{}{r}_2,\stackrel{}{r}_1)}`$
$`=`$ $`{\displaystyle \frac{1}{2}}W_af_aa_aM_a^{(\tau \overline{\tau })}.`$ (A13)
The third line of A4 contains isospin-projection operators which have a value $`+1`$ when $`i`$ and $`j`$ are like particles and $`1`$ when they are unlike. In the direct term this gives an energy of

$`{\displaystyle \frac{1}{2}}W_af_ab_a{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\rho _p(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_1)\mathrm{d}^3\stackrel{}{r}_2\rho _p(\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_2)}`$ (A17)
$`+{\displaystyle \frac{1}{2}}W_af_ab_a{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\rho _n(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_1)\mathrm{d}^3\stackrel{}{r}_2\rho _n(\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_2)}`$
$`{\displaystyle \frac{1}{2}}W_af_ab_a{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\rho _p(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_1)\mathrm{d}^3\stackrel{}{r}_2\rho _n(\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_2)}`$
$`{\displaystyle \frac{1}{2}}W_af_ab_a{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\rho _n(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_1)\mathrm{d}^3\stackrel{}{r}_2\rho _p(\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_2)}`$
$`=`$ $`{\displaystyle \frac{1}{2}}W_af_ab_a\left[{\displaystyle \mathrm{d}^3\stackrel{}{r}\rho _p(\stackrel{}{r})\rho ^{\beta _a}(\stackrel{}{r})}\right]^2+{\displaystyle \frac{1}{2}}W_af_ab_a\left[{\displaystyle \mathrm{d}^3\stackrel{}{r}\rho _n(\stackrel{}{r})\rho ^{\beta _a}(\stackrel{}{r})}\right]^2`$ (A19)
$`W_af_ab_a{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\rho _n(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_1)\mathrm{d}^3\stackrel{}{r}_2\rho _p(\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_2)}`$
$`=`$ $`{\displaystyle \frac{1}{2}}W_af_ab_a\left[{\displaystyle \mathrm{d}^3\stackrel{}{r}(\rho _p(\stackrel{}{r})\rho _n(\stackrel{}{r}))\rho ^{\beta _a}(\stackrel{}{r})}\right]^2={\displaystyle \frac{1}{2}}W_af_ab_a(\mathrm{\Delta }N_a)^2.`$ (A20)
The exchange term gives a contribution only when $`i`$ and $`j`$ have the same isospin quantum number, which is just like the case for having no isospin operator there, so the contribution is like that in (A10)

$$\frac{1}{2}W_af_ab_aM_a.$$
(A21)
The HF mean-field is obtained by varying the total energy with respect to the single particle states. This gives, for the case of the attractive term, without the explicit isospin dependence:

$`{\displaystyle \frac{\delta }{\delta \varphi _b(\stackrel{}{x})}}\left({\displaystyle \frac{1}{2}}W_af_aN_a^2{\displaystyle \frac{1}{2}}W_af_aM_a\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}W_a{\displaystyle \frac{\delta f_a}{\delta \varphi _b^{}(\stackrel{}{x})}}N_a^2+W_af_aN_a{\displaystyle \frac{\delta N_a}{\delta \varphi _b^{}(\stackrel{}{x})}}`$ (A22)
$``$ $`{\displaystyle \frac{1}{2}}W_a{\displaystyle \frac{\delta f_a}{\delta \varphi _b^{}(\stackrel{}{x})}}M_a{\displaystyle \frac{1}{2}}W_af_a{\displaystyle \frac{\delta M_a}{\delta \varphi _b^{}(\stackrel{}{x})}}.`$ (A23)
The variation of the function $`f_a`$ is given by

$`{\displaystyle \frac{\delta f_a}{\delta \varphi _b^{}(\stackrel{}{x})}}`$ $`=`$ $`{\displaystyle \frac{\delta }{\delta \varphi _b^{}(\stackrel{}{x})}}\left[{\displaystyle \rho ^{\alpha _a}(\stackrel{}{r})\mathrm{d}^3\stackrel{}{r}}\right]^1`$ (A24)
$`=`$ $`f_a^2{\displaystyle \mathrm{d}^3\stackrel{}{r}\frac{\rho ^{\delta \alpha _a}(\stackrel{}{r})}{\delta \varphi _b^{}(\stackrel{}{x})}}`$ (A25)
$`=`$ $`f_a^2\alpha _a{\displaystyle \mathrm{d}^3\stackrel{}{r}\rho ^{\alpha _a1}(\stackrel{}{r})\frac{\delta \rho (\stackrel{}{r})}{\delta \varphi _b^{}(\stackrel{}{x})}}`$ (A26)
$`=`$ $`f_a^2\alpha _a\rho ^{\alpha _a1}(\stackrel{}{x})\varphi _b(\stackrel{}{x})`$ (A27)
so that the contributions of the two terms in (A23) involving the variation of $`f_a`$ give a contribution to the HF mean–field of

$$W_a(\alpha _a/2)f_a^2(N_a^2M_a)\rho ^{\alpha _a1}(x)$$
(A28)
The function $`N_a`$ is similar in form to $`f_a`$ and the functional variation proceeds in a similar manner;

$`{\displaystyle \frac{\delta N_a}{\delta \varphi _b^{}(\stackrel{}{x})}}`$ $`=`$ $`{\displaystyle \frac{\delta }{\delta \varphi _b^{}(\stackrel{}{x})}}{\displaystyle \mathrm{d}^3\stackrel{}{r}\rho ^{\beta _a+1}(\stackrel{}{r})}`$ (A29)
$`=`$ $`(\beta _a1)\rho ^{\beta _a}(\stackrel{}{x})\varphi _b(\stackrel{}{x})`$ (A30)
and the contribution from the second term in (A23) to the mean-field is

$$W_af_aN_a(\beta _a1)\rho ^{\beta _a}(\stackrel{}{x}).$$
(A31)
Finally, the functional variation of the exchange matrix element $`M_a`$ is

$`{\displaystyle \frac{\delta M_a}{\delta \varphi _b^{}(\stackrel{}{x})}}`$ $`=`$ $`{\displaystyle \frac{\delta }{\delta \varphi _b^{}(\stackrel{}{x})}}\left({\displaystyle \underset{ij<ϵ_F}{}}{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\varphi _i^{}(\stackrel{}{r}_1)\varphi _j^{}(\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_2)\varphi _j(\stackrel{}{r}_1)\varphi _i(\stackrel{}{r}_2)}\right)`$ (A32)
$`=`$ $`2{\displaystyle \underset{ij<ϵ_F}{}}{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\frac{\delta \varphi _i^{}}{\delta \varphi _b^{}(\stackrel{}{x})}\varphi _j(\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{r}_1)\rho ^{\beta _a}(\stackrel{}{r}_2)\varphi _j(\stackrel{}{r}_1)\varphi _i(\stackrel{}{r}_2)}`$ (A33)
$`+`$ $`2{\displaystyle \underset{ij<ϵ_F}{}}{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\mathrm{d}^3\stackrel{}{r}_2\varphi _i^{}(\stackrel{}{r}_1)\varphi _j^{}(\stackrel{}{r}_2)\frac{\delta \rho ^{\beta _a}(\stackrel{}{r}_1)}{\delta \varphi _b^{}(\stackrel{}{x})}\rho ^{\beta _a}(\stackrel{}{r}_2)\varphi _j(\stackrel{}{r}_1)\varphi _i(\stackrel{}{r}_2)}`$ (A34)
$`=`$ $`2{\displaystyle \underset{j<ϵ_F}{}}{\displaystyle \mathrm{d}^3\stackrel{}{r}_2\varphi _j^{}(\stackrel{}{r}_2)\rho ^{\beta _a}(\stackrel{}{x})\rho ^{\beta _a}(\stackrel{}{r}_2)\varphi _j(\stackrel{}{x})\varphi _b(\stackrel{}{r}_2)}`$ (A35)
$`+`$ $`2{\displaystyle \underset{ij<ϵ_F}{}}{\displaystyle \mathrm{d}^3\stackrel{}{r}_2\varphi _i^{}(\stackrel{}{x})\varphi _j(\stackrel{}{r}_2)\beta _a\rho ^{\beta _a1}(\stackrel{}{x})\varphi _b(\stackrel{}{x})\rho ^{\beta _a}(\stackrel{}{r}_2)\varphi _j(\stackrel{}{x})\varphi _i(\stackrel{}{r}_2)}`$ (A36)
where use is made of the symmetry of the integral to combine the four terms into two. The last term in (A36) gives rise to a local term in the mean-field of

$`W_af_a\beta _a{\displaystyle \underset{i<ϵ_F}{}}\left({\displaystyle \mathrm{d}^3\stackrel{}{r}\rho (\stackrel{}{r},\stackrel{}{x})\rho ^{\beta _a}(\stackrel{}{r})\rho (\stackrel{}{x},\stackrel{}{r})}\right)\rho ^{\beta _a1}(\stackrel{}{x})`$ (A37)
$`=`$ $`W_af_a\beta _aG_a(\stackrel{}{x})\rho ^{\beta _a1}(\stackrel{}{x})`$ (A38)
where $`G_a(\stackrel{}{x})`$ has been defined as in (31).

The other term in (A36) gives rise to a truly non-local Fock term in the mean-field:

$$U(\stackrel{}{x},\stackrel{}{x}^{})\varphi _b(\stackrel{}{x})=W_af_a\underset{i<ϵ_F}{}\rho ^{\beta _a}(\stackrel{}{x})\varphi _i(\stackrel{}{x})\left[\mathrm{d}^3\stackrel{}{r}\varphi _i^{}(\stackrel{}{r})\rho ^{\beta _a}(\stackrel{}{r})\varphi _b(\stackrel{}{r})\right].$$
(A39)
This completes the non-isospin-dependent part of the attractive force, and so also the repulsive by change of subscript. The isospin-dependent terms are obtained in an analogous way, except that when the variation applies to the density of a single nucleon species, so the contribution to the mean-field applies only to that species.

For the final term in equation (10), the so-called derivative term, only the direct part of the energy is at present considered. It is

$`E_{\mathrm{deriv}.}`$ $`=`$ $`{\displaystyle \frac{1}{2}}k{\displaystyle \underset{ij<ϵ_F}{}}ij|_1^2\rho (\stackrel{}{r}_1)_2^2(\stackrel{}{r}_2)|ij`$ (A40)
$`=`$ $`{\displaystyle \frac{1}{2}}k{\displaystyle \mathrm{d}^3\stackrel{}{r}_1\rho (\stackrel{}{r}_1)_1^2\rho (\stackrel{}{r}_1)\mathrm{d}^3\stackrel{}{r}_2\rho (\stackrel{}{r}_2)_2^2\rho (\stackrel{}{r}_2)}={\displaystyle \frac{1}{2}}kN_d^2.`$ (A41)
The functional variation proceeds as

$`{\displaystyle \frac{\delta E_{\mathrm{deriv}.}}{\delta \varphi _b^{}(\stackrel{}{x})}}`$ $`=`$ $`kN_d{\displaystyle \frac{\delta }{\delta \varphi _b^{}(\stackrel{}{x})}}{\displaystyle \mathrm{d}^3\stackrel{}{r}\rho (\stackrel{}{r})^2\rho (\stackrel{}{r})}`$ (A42)
$`=`$ $`kN_d{\displaystyle \mathrm{d}^3\stackrel{}{r}\left\{\frac{\delta \rho (\stackrel{}{r})}{\delta \varphi _b^{}(\stackrel{}{x})}\right\}^2\rho (\stackrel{}{r})}+kN_d{\displaystyle \mathrm{d}^3\stackrel{}{r}\rho (\stackrel{}{r})\left\{\frac{\delta }{\delta \varphi _b^{}(\stackrel{}{x})}^2\rho (\stackrel{}{r})\right\}}.`$ (A43)
The first term gives a contribution to the mean field of

$$kN_d^2\rho (\stackrel{}{x}).$$
(A44)
By integrating the second term by parts twice, one in fact gets exactly the same contribution to the mean-field again, so that the total contribution to the mean-field from the direct term of the derivative interaction is

$$2kN_d^2\rho (\stackrel{}{r})$$
(A45)
The exchange part of this term is not calculated.


text:
# 1 The ”classical” Tip of the Red Giant Method

## 1 The ”classical” Tip of the Red Giant Method

Figure 1 shows the (I, V-I) and (I, B-I) CM diagrams of Antlia. As far as the data reduction is concerned, the interested reader is referred to Piersimoni et al. (1999). Data plotted in these CMDs present two key features: a well-defined Red Giant Branch, and a sample of bright blue stars, belonging to a young stellar component. By adopting the new calibration of the Tip of the Red Giant Branch (TRGB) provided by Salaris & Cassisi (1998) we estimated that the Antlia distance modulus is (m-M)0=25.98 $`\pm `$0.10 mag, equivalent to D=1.51$`\pm `$0.07 Mpc.

The main aim of this investigation is to present a new calibration of the TRGB method based on NIR magnitudes which can allow us to supply reliable distance determinations for both metal-poor and metal-rich stellar populations. In fact empirical evidence suggest not only that the the mean metallicity of LG galaxies range from $`[Fe/H]=2.2\pm 0.1`$ (Ursa Minor) to $`[Fe/H]=0.8\pm 0.1`$ (NGC 205), but also that the average metallicity spread is of the order of 0.5 dex (Mateo 1998). Therefore, the TRGB method cannot be used to supply a homogeneous distance scale in the LG, since both theory and observations support the evidence that I magnitude of the TRGB stars is constant within 0.1 mag, for stellar populations with metal contents \[Fe/H\] ranging from -2.2 to -0.7 (Lee, Freedman, & Madore 1993; Salaris & Cassisi 1997,1998).

In order to overcome this drawback and to improve the TRGB method in such a way that could be safely used to estimate the distances of stellar systems which include both metal-poor and metal-rich stars, we investigated on theoretical grounds the properties of TRGB stars in the NIR bands. We focused our attention to NIR passbands because they present a negligible dependence on interstellar reddening and also a marginal dependence on metallicity.

## 2 Discussion

The behavior of TRGB magnitude in NIR bands (I,J,H,K) was investigated by adopting the theoretical predictions by Salaris & Cassisi (1997, 1998) which were transformed into the observational plane using the bolometric corrections and the color-temperature relations provided by Wood & Bessel (1994).

Panel (a) of Figure 2 shows the well known behavior of TRGB magnitudes in the I band, i.e. the I magnitude of the tip is not a reliable standard candle for metallicities greater that \[Fe/H\]=-0.7. Panel (b) of the same figure shows the J band magnitude of TRGB stars. Data plotted in this panel clearly show that the TRGB magnitude in the J band cannot be adopted as a standard candle, since it shows a steady decrease when moving from metal-poor to metal rich structures. The same outcome applies for H and K magnitudes (see panels c and d). A closer inspection, however, reveals that the behavior on this two filters is quite similar and in particular that they present the same slope. This evidence prompted us to estimate the ratio MH/MK as a function of metallicity.

Interestingly enough, we found that the ratio MH/MK shows a maximum variance of approximately 0.01 mag over the entire metallicity range (see bottom panel) we have taken into account ($`2.5[M/H]0.2`$). This result suggests that HK ratio can be adopted as a standard candle, since it presents a negligible dependence on metallicity.

## 3 A possible improvement of the TRGB method

On the basis of the result of the previous section, we can assume that MH/MK $``$ const. = a. Then, by accounting for the evidence that the extinction in the H and K bands are quite similar (for EB-V = 0.1, we have AH-AK$``$0.02 mag), we obtain

mH \- MH $``$ mK \- MK.

At this point, we can simply derive

MK $``$ (mH \- mK) / ( a - 1) = ( H - K ) / ( a - 1 ).

This relation suggests that dereddened measurements of the TRGB H and K magnitudes, allow us to estimate the distance modulus. Once this ”new” method is supported by empirical evidence it could be adopted to estimate distances of LG galaxies whose stellar populations are older few Gyrs.

Let us briefly mention some drawbacks of the distance indicators currently adopted:

(1) The TRGB cannot be safely applied to all LG galaxies, mainly because of the already quoted ”metallicity problem”.

(2) Due to the fact that the Horizontal Branch (HB) is intrinsically several magnitude fainter than TRGB, distances based on $`M_V`$(HB) can be derived only for nearby stellar systems.

(3) Distances based on ”Cepheid PL relations” can be derived only for young stellar systems.

On the other hand, this new approach can be adopted in LG stellar systems characterized by a wide range of ages and metallicities. In fact, modern NIR detectors aboard on VLT (ISAAC) and on NTT (SOFI) allow us to reach H and K magnitudes roughly equal to 19 mag. By assuming a S/N=20, and a seeing of 1.0 arcsec, a limit magnitude K=19 can be reached within time exposures ranging from 650 sec (ISAAC) to 3000 sec (SOFI). It turns out that the TRGB can be measured in galaxies with distance moduli up to 24-25. This means that this approach can be applied for approximately 70% of southern LG galaxies.

The ”a” parameter of the previous equation has to be tested and calibrated on NIR data of stellar systems in the LG. Moreover, to test the intrinsic accuracy and reliability of this method we need to collect new NIR data for galactic globular clusters and LG dwarf galaxies whose distance is already well known from other distance indicators.