text:
# Complex molecular gas structure in the Medusa merger

## 1 Introduction

The Medusa merger, NGC 4194, (Table 1; Fig. 1) belongs to a class of lower luminosity ($`L_{\mathrm{FIR}}=8.5\times 10^{10}`$ L at $`D`$=39 Mpc) mergers, an order of magnitude fainter than well known Ultra Luminous Galaxies (ULIRGs ($`L_{\mathrm{FIR}}\genfrac{}{}{0pt}{}{_>}{^{}}10^{12}`$ L)) such as Arp 220, but with an extended region of intense star formation (Prestwhich et al. 1994; Armus et al. 1990). As the name indicates, the morphology of the Medusa merger is as spectacular as that of the brighter objects. The “hair” or “tentacles” are a patchy tidal tail stretching out 60′′ (11.3 kpc) north of the main body, which appears twisted. There also is a sharp, curved feature (and a very faint tail) to the south.

NGC 4194 has been mapped in 12CO 1-0 and 2-1 with a single dish telescope by Casoli et al. (1992) and detected in the more highly excited 12CO 3-2 line by Devereux et al. (1994). It is one of the handful known mergers (e.g. Arp 220, IC 694, NGC 6240) with unusually high 12CO/13CO 1–0 intensity ratios, $`>`$20 (e.g. Aalto et al. 1991b; Casoli et al. 1992). In contrast, the majority of large galaxies exhibit $``$ in the range 5 — 15 (e.g. Aalto et al. 1995). One feature these HR (High Ratio) mergers have in common is a high FIR $`f(60\mu \mathrm{m})/f(100\mu \mathrm{m})`$ flux ratio $`\genfrac{}{}{0pt}{}{_>}{^{}}0.8`$ – which likely is an indication of a high dust temperature (e.g. Aalto et al. 1991a, 1991b). High resolution 12CO maps of several luminous ($`L_{\mathrm{IR}}\genfrac{}{}{0pt}{}{_>}{^{}}10^{11}`$ L) HR mergers reveal that nearly all the molecular gas is concentrated within a few hundred parsecs of the center (e.g. Scoville et al. 1989; Bryant & Scoville 1996; Downes & Solomon 1998).

The resulting large central surface densities require high pressures in hydrostatic equilibrium and the molecular ISM becomes highly turbulent. Large turbulent line widths in molecular clouds will result in reduced optical depths in the 12CO 1–0 line and the 12CO/13CO 1–0 line ratio will be elevated compared to the optically thick ratio (for a discussion see Aalto et al. 1995). So for the luminous HR mergers this could be one important mechanism behind the large observed 12CO/13CO line ratios (Aalto et al. 1995), aside from the additional effect from large kinetic temperatures. Will the lower luminosity HR galaxy NGC 4194 also follow this pattern of a compact nuclear starburst, or will an extended 12CO distribution require a different dominant mechanism behind the elevated $``$?

Most of the high resolution imaging effort of 12CO in mergers has been directed toward the ULIRGs and not much is known, statistically, of the more modest luminosity merging galaxies. Apart from addressing the question of the compactness of the 12CO emission in an HR merger, a 12CO image of NGC 4194 will improve the statistics of the 12CO distribution in moderate luminosity mergers.

The elevated values of $``$ in hot dust mergers have also been interpreted as a global underabundance in 13CO with respect to 12CO and not an effect of reduced optical depths in 12CO (e.g. Casoli et al. 1992). However, the moderate optical depth ($`\tau 1`$) in 12CO 1–0 towards centers of galaxies in general, and ULIRGs in particular, is now a reasonably well established fact (e.g. Aalto et al. 1997; Dahmen et al. 1998; Downes & Solomon 1998).

We have imaged NGC 4194 in the 12CO 1–0 line with OVRO to study the distribution and kinematics of the molecular gas in an HR merger galaxy of moderate luminosity. We have also obtained single dish spectra of the 1–0 transitions of 12CO, 13CO and HCN. The observations are presented in Sect. 2; the results in Sect. 3; the consequences of these findings are discussed in terms of merger history, central gas surface densities and final fate of the molecular gas in Sect. 4.

## 2 Observations

We have obtained high resolution maps of 12CO 1–0 using the Caltech six-element Owens Valley Radio Observatory (OVRO) millimeter array. Two tracks were taken in the low resolution configuration in September and October 1996 and two in the high resolution configuration in December 1996. The naturally weighted synthesized beam size is $`2.^{\prime \prime }5\times 2.^{\prime \prime }0`$ ($`472\times 378`$ pc for $`D`$=39 Mpc) and the beam position angle (BPA) is –66. The sensitivity of this map is 11 mJy beam-1 channel-1, corresponding to 0.19 K channel-1. For a beamsize of $`2`$$`.^{\prime \prime }`$3 a brightness temperature ($`T_\mathrm{B}`$) of 1 K corresponds to 57 mJy at $`\lambda `$=2.6 mm.

For the (naturally weighted) high resolution data only, the beamsize is $`2.^{\prime \prime }0\times 1.^{\prime \prime }7`$ and BPA is –71. For a resolution of $`1`$$`.^{\prime \prime }`$85 a brightness temperature of 1 K corresponds to 36 mJy.

System temperatures were 500–600 K. The quasar 1150+497 was used for phase calibration and Neptune and Uranus for absolute flux calibration. The digital correlator, centered at $`V_{\mathrm{LSR}}=2560`$ $`\mathrm{km}\mathrm{s}^1`$, provided a total velocity coverage of 1170 $`\mathrm{km}\mathrm{s}^1`$. Data were binned to 4 MHz resolution, corresponding to 10 $`\mathrm{km}\mathrm{s}^1`$, and to construct the map, the resolution was reduced to 20 $`\mathrm{km}\mathrm{s}^1`$.

Single dish spectra of the 1–0 transition of 12CO, 13CO and HCN were obtained toward the center of NGC 4194 with the OSO 20m telescope. The observations were carried out in October 1999, and the system temperatures were typically 400 K for 12CO and 13CO and 200–300 K for HCN. The pointing accuracy was checked on SiO masers and was found to be about 3′′. Beam sizes and efficiencies are listed in Table 2, together with the results.

## 3 Results

### 3.1 Morpholology and kinematics of the molecular distribution

Integrated intensity: The naturally weighted 12CO emission (Fig. 2a) is distributed on a scale of $`30^{\prime \prime }`$ (5.7 kpc). The bulk of the emission (67 %) is concentrated in an 10′′ (1.9 kpc) sized oval structure centered on the optically bright main body. Fainter, patchy emission curves out to the north from the brighter emission.

An overlay of the 12CO distribution on an R-band image (Fig. 3) reveals that the general 12CO morphology is quite different from that of the main optical body of the galaxy. The molecular gas extends from the center to the north-eastern part of the object, and the major axis of the central 12CO intensity distribution is almost perpendicular to that of the optical body. The position angle (PA) of the 12CO oval intensity distribution of the bright emission is $`49^{}`$ and the FWHM source size is $`4.^{\prime \prime }2\times 2`$$`.^{\prime \prime }`$5 The position angle of the optical emission is 146.

The curved northern 12CO emission reaches into the base of the tidal tail, and some faint emission is also detected near a bright patch (indicated with an arrow in Fig. 3) in the tail itself. In Fig. 3, the 12CO emission appears to be inside of an optical protrusion, or arm, in the east (marked with an arrow in Fig. 3). Armus et al. (1990) describe this feature as a “spur” or a “jet”. The blue band image in Arp’s Catalog of Peculiar Galaxies (1966) gives instead a strong indication of a foreground eastern dust lane blocking the background emission and creating the impression of a protruding feature separate from the main body. The molecular emission traces the space between the “spur” and the main body which is consistent with the interpretation of a dust lane.

The total 12CO flux in NGC 4194 detected by OVRO is 170 Jy $`\mathrm{km}\mathrm{s}^1`$ which, for a standard 12CO to H2 mass conversion factor, translates to $`2.3\times 10^9`$ M ($`0.92\times 10^4S\mathrm{\Delta }VD^2`$ M for $`X`$ = N(H2)/I(12CO) = $`2.3\times 10^{20}`$ $`\mathrm{cm}^2`$ (K $`\mathrm{km}\mathrm{s}^1`$)-1) of molecular gas. Within calibrational error bars, we recover all of the single dish flux. The peak flux is 27 Jy beam-1 $`\mathrm{km}\mathrm{s}^1`$ (16% of the total flux) at $`\alpha `$: $`12^\mathrm{h}11^\mathrm{m}41.22^\mathrm{s}`$ ; $`\delta `$: $`54^{}48^{}14.^{\prime \prime }4`$. which translates to a molecular mass of $`3.8\times 10^8`$ M within a radius of 240 pc. However, we suspect that the conversion factor overestimates the molecular mass in this system — in particular in the center (see Sect. 4.3.1.).

The velocity field: In Fig. 2b we present the 12CO velocity field which shows the complex dynamics of NGC 4194. Close to the center there seems to be organized rotation of the gas (even though the velocity contours are far from a regular spider diagram) with a position angle of 160 — close to the optical major axis. Thus, the major axis of the rotation of the gas follows the PA of the optical light distribution, rather than the distribution of the 12CO emission. A position-velocity (pV) cut at PA 160 (Fig. 4) indicates a solid-body rotation curve which starts to flatten (or fall) at a radius of 2′′ (380 pc). The Figure also shows that the 12CO peak is offset from the kinematic center by $`1^{\prime \prime }`$. The projected rotational velocity is 90 $`\mathrm{km}\mathrm{s}^1`$. Thus, for an inclination of 40 (from the optical isophotes) the deprojected rotational velocity is 140 $`\mathrm{km}\mathrm{s}^1`$ in the inner 16′′ (3 kpc). The systemic velocity at the kinematic center is 2510 $`\mathrm{km}\mathrm{s}^1`$.

The extended northern gas in general does not move according to the pattern established in the inner region. The eastern part of the emission is at higher velocities than the western emission and velocity discontinuities can be traced along the curved ridge at the transition between the velocity systems. The highest velocity in the system, 2700 $`\mathrm{km}\mathrm{s}^1`$, can be found on the eastern dust lane at $`\delta `$: $`54^{}48^{}22^{\prime \prime }`$. The velocity of the gas in the eastern dust lane then gradually drops along the lane, until it reaches the tidal tail where the velocity of the gas clumps is 2580 $`\mathrm{km}\mathrm{s}^1`$.

The Dispersion Map: Fig. 5 displays the dispersion map with a few selected spectra. The one dimensional dispersion is quite moderate, $`\sigma _{1\mathrm{d}}\genfrac{}{}{0pt}{}{_<}{^{}}10`$ $`\mathrm{km}\mathrm{s}^1`$ in large parts of the galaxy. Regions of high dispersion are found where the velocity field shows discontinuities or irregularities, and the sample spectra reveal that here the enhanced velocity dispersion is caused by multiple spectral features rather than genuinely broad lines. The velocity shifts between the narrow peaks is 150–200 $`\mathrm{km}\mathrm{s}^1`$, i.e. the typical rotational velocity of a spiral galaxy. There is also a continuous curved region of elevated dispersion (i.e. multiple spectra) in the center of the map. Even the width of the spectrum at the intensity peak is quite low (lower panel, right hand side) with $`\sigma _{1\mathrm{d}}=50`$ $`\mathrm{km}\mathrm{s}^1`$, which is narrow for the center of a merging system.

Channel Maps: Dividing the emission into three velocity bins, (2380 – 2485 $`\mathrm{km}\mathrm{s}^1`$), (2506 – 2570 $`\mathrm{km}\mathrm{s}^1`$), (2590 – 2720 $`\mathrm{km}\mathrm{s}^1`$), (Figs. 6a,b,c) shows the emission being distributed in long, arm-like features of quite different appearence in the three bins. The low and high velocity bins were chosen to include each one of the two features seen in the double peaked spectra (e.g. middle right panel in Fig. 5). The middle bin includes emission at intermediate velocities near the systemic value.

The low velocity map (Fig. 6a) shows a patchy distribution extending over 20′′ to the north-south. Closer to the center, brighter emission curves 5′′ in the east-west direction.

At mid-velocities (Fig. 6b), the distribution is dominated by a striking, elongated curved shape, approximately 15′′ (2.8 kpc) in length. This structure is the main cause of the discrepancy between the 12CO emission distribution and the R-band light distribution. It traces a central prominent dust lane. To show this, we have overlayed the mid-velocity 12CO emission on an archival HST image in Fig. 7. Another, fainter, 20′′ (3.8 kpc) long curved structure is found in the north. This feature may be a continuation of the central dust lane and the gap between them a depression in the 12CO brightness.

The high velocity emission bin (Fig. 6c) shows a bright peak in the center with a source size of 2$`.^{\prime \prime }5\times 1.^{\prime \prime }\text{8}`$. A prominent, narrow curved feature stretches out to the north-east, (20′′ to the north) and follows most of the eastern obscuring dust lane seen in Arp’s B-band image. We suggest the northern, curved feature in the mid-velocity bin is the continuation of this curve. As mentioned above, the central dust lane may connect to this feature as well. In Fig. 6d the features in the three velocity bins are marked with dashed lines.

#### 3.1.1 The High Resolution Map — a close-up on the center

For a closer look at the structure of emission in the inner 10′′, we selected the high resolution data only, and thus increased the resolution to $`1.^{\prime \prime }7`$ (but lost some sensitivity and $`uv`$-coverage). The integrated intensity map (Fig. 8a) of the inner region consists of 5 peaks, $`a`$, $`b`$, $`c`$, $`d`$ and $`e`$ (Table 3) and extended emission surrounding them. The total flux of the inner $`10^{\prime \prime }\times 10^{\prime \prime }`$ is 64 Jy km/s. The brightest peak, $`a`$ (17 Jy $`\mathrm{km}\mathrm{s}^1`$ beam-1), is at the same position as in the lower resolution map, i.e. close to the radio continuum peak, but shifted to the south by $`1.^{\prime \prime }5`$ (280 pc). Feature $`a`$ has a size of $`1.^{\prime \prime }8\times 1.^{\prime \prime }3`$ and a position angle of $`45^{}`$. The fitted FWHM sizes of features $`a`$$`e`$ are typically 200–400 pc, which is comparable to the sizes of Giant Molecular Associations (GMAs) found in other galaxies. Features $`b`$ and $`c`$ can be identified with knots in the H$`\alpha `$ distribution (Armus et al. (1990)). The peak intensity in the high resolution map is 200 mJy, which corresponds to 5.4 K.

Fig. 9 is an overlay of the integrated high resolution 12CO emission on an archival HST image (enhanced to show the dust lane, same as Fig. 5). The peaks $`a`$, $`b`$, $`d`$ and $`e`$ all lie on the central dust lane.

The velocity field (Fig. 8b) is slightly more twisted than in the lower resolution map and covers the same velocity range. Features $`d`$ and $`c`$ are causing twists in the velocity field because they are at too low velocities to fit a spider diagram. The maximum dispersion (100 $`\mathrm{km}\mathrm{s}^1`$) in Fig. 8c does not occur in $`a`$, but in an elongated (east-west) region located north of the 12CO peak, close to feature $`b`$. This is caused by the crossing of the dust lane, i.e. by a superposition of multiple spectral features.

Fig. 10 shows an overlay of the 12CO emission on the 1.49 GHz radio continuum structure (Condon et al. 1990). There is a rough agreement between the CO and the radio-continuum distribution which seems to break down at smaller scales (i.e. less than 1 kpc). The 12CO peak (feature $`a`$), for example, is located $`1^{\prime \prime }`$ south of the radio continuum peak. There is an arm-like structure curving from the center to the west which has a 12CO counterpart in feature $`c`$. The one 12CO feature entirely lacking in the continuum distribution is $`e`$ in the north. Comparisons between radio continuum and CO distributions in starburst galaxies often show a correlation between the two distributions (e.g. Yun et al. 1994, Bajaja et al. 1995). This effect is often attributed to cosmic ray heating of the molecular gas (e.g. Adler et al. 1991; Allen 1992). We will discuss this correlation further in Sect. 4.2.

### 3.2 Single dish results

Results are presented in Table 2 and 12CO, 13CO spectra in Fig. 12 (HCN was not detected, and its spectrum is not shown). The 12CO flux detected by OVRO is, within calibrational error bars, the same as measured by the OSO 20m telescope, which is approximately 60% more than the flux Devereux et al. (1994) find with the FCRAO 14m telescope.

We note that the integrated 12CO/13CO intensity ratio we find is lower, $`20`$, than previously reported by Casoli et al. (1992), $`55`$, although their error bar of $`\pm 22`$ includes our measurement. Our result means that NGC 4194 is not an extreme HR merger and 13CO 1–0 is quite detectable, even if the ratio is clearly above what is usually found for starburst galaxies. There is a hint of a line shape difference in the 12CO and 13CO lines with an extra narrow component at lower velocities in the 13CO spectrum. The signal to noise ratio in the 13CO spectrum is too low to warrant further speculation on this matter here, but a high resolution 13CO image may reveal interesting excitation changes across the galaxy, as has already been shown for the merger starburst system Arp 299 (Aalto et al. 1997; Casoli et al. 1999). The HCN emission is at least as faint as the 13CO emission with the 12CO/HCN intensity ratio $`\genfrac{}{}{0pt}{}{_>}{^{}}25`$. In comparison to other FIR luminous galaxies (e.g. Solomon et al. 1992) the HCN is relatively faint.

## 4 Discussion

### 4.1 Merger timescale and 12CO morphology

#### 4.1.1 The Medusa as a major disk-disk merger

Joseph & Wright (1985) suggested that NGC 4194 is a collision between two disk galaxies and put it fairly late in a merger sequence: after Arp 220 and before NGC 3310. They used the apparent diffuseness of the tidal tail and the coalesced optical main body as principal arguments for this placement in the evolutionary scheme.

However, despite the suggested advanced stage of the NGC 4194 merger, the molecular gas structure is complex and distributed on a scale of several kpc. In the late stages of a merger, simulations suggest that the bulk of the molecular gas would have sunk to the inner 200 pc of the merger, regardless of the previous morphology/orientation of the precursor galaxies (e.g. Barnes & Hernquist 1996). The corresponding size scale for NGC 4194 is an order of magnitude greater than that. The funneling of molecular gas toward the center of NGC 4194 has been less effective than for ULIRGs, such as Arp 220 or NGC 6240, which is interesting since both these mergers appear younger than NGC 4194. Apart from the ULIRGs being more luminous by factors of at least 10, they also have broader line widths, and deeper potentials. If NGC 4194 is the result of the merger of two disk galaxies of comparable mass, i.e. a ‘major merger’, either the merger is less advanced than previously thought, or the intrinsic and/or orbital properties of the NGC 4194 progenitors are, after all, important factors.

One can use the presence and separation of a double nucleus as an indication of merger age. In NGC 4194 there is no obvious sign of another nucleus, but there is a fairly bright extension to the south of the radio continuum emission, which is the location of the 12CO peak, $`a`$. If $`a`$ is another nucleus in the merging system, then the projected separation between them is $`1.^{\prime \prime }5`$, or 280 pc. However, the non-compact appearance (and lower lumniosity) of $`a`$ indicates instead that it is a bright star forming region rather than a second nucleus, in which case the upper limit of the nuclear separation becomes even more strict. It is perhaps also possible that merger orientation hides the second nucleus behind the first one, but this should be a rather unlikely scenario.

We find it difficult to reconcile the extended, unrelaxed appearance of the molecular gas distribution with a major merger scenario where the lack of a double nucleus indicates a very advanced merger. Together with the morphological appearance of the Medusa, this lead us to suggest another merger scenario:

#### 4.1.2 The Medusa as a young shell-galaxy

The southern sharp, curved feature seen in the optical image may well be a shell in its early stages of formation, thus raising questions on whether an elliptical/early-type galaxy was one of the progenitors. The total molecular mass of $`2\times 10^9`$ M is indeed less than what one would expect from the collision of two gas-rich disks.

Quinn (1984) found that shells can form through “phase-wrapping” of dynamically cold disk material in a rigid, spherical potential. The similarity between the overall Medusa morphology and simulations of a radial, planar collision between a disk galaxy and a more massive spherical, is striking (e.g. Quinn 1984, Fig. 1 timestep 6; Kojima & Noguchi (KN) 1997, Fig. 4a timestep 4.7, 6.2). The first shell is forming on the opposite side of the single, broad tidal tail creating an almost “cometary” structure. Although, we note that the implied molecular mass of NGC 4194 suggests that the disk galaxy involved was more massive than a dwarf — which is often assumed in the simulations.

KN also simulated the behaviour of the gaseous component during the encounter and found that the gas becomes scattered much like the disk stars. This is in sharp contrast to the models of Weil & Hernquist (WH)(1993) who instead conclude that the gas from the disk galaxy quickly sinks to the center of the merger. This also leads to different views on the occurence of starbursts in these encounters. In the KN approach, the scattering of gas clouds prevents star formation, while WH suggest that starbursts are likely to occur in the central gas concentration. In the WH model the shells are expected to be gas-poor because of the separation between gas and stars. The difference between the KN and WH models is largely caused by their ISM represenation: KN using sticky particles and WH a hydrodynamical model.

Now, if we adopt the notion that NGC 4194 is indeed a radial collision between a spiral and a more massive, gas-poor spherical (i.e. an unequal mass, or ‘minor’ merger), we can make the following observations: No molecular gas is detected in the southern shell in accordance with the WH model, but perhaps more sensitive observations can find gas here. Both HI and CO has been detected in the shells of Centaurus A (Schiminovich et al. 1994; Charmandaris, Combes & van der Hulst 2000). The molecular gas is not scattered all over the system, nor can all of it be found in the very inner region. The distribution is somewhere in between the two models, and the gas in the inner region seems to be going through a starburst on a fairly large scale (see Sect.4̇.2). The multiple spectral components may reflect gas on different orbits and in different stages of its dance around the center of the elliptical. The velocity shift between the two components is roughly the same as the typical rotational velocity for a spiral galaxy. Thus the existence of two spectral components does not necessarily imply that the molecular gas is coming from two different progenitor galaxies. For both the major and minor merger models, however, the presence of two spectral components suggests that the gas has some time to go before it reaches relaxation.

However, the elliptical-spiral model leaves us with (at least) two quesions: 1) where is the massive elliptical? 2) what is the significance of the (very) faint tail to the south (which can only be seen in deep images)? The first question can better be answered once there is a near infrared K-band image that can see through the dust lanes and dust enshrouded starburst currently hiding the center. The southern tail is very diffuse, and may be a remnant of an earlier interaction. Its significance can be better judged by sensitive HI observations.

We have made a search through galaxy morphology catalogues looking for Medusa-like mergers. We found three more objects with similar morphology (i.e. one prominent tidal tail and on the opposite side of the galaxy a sharp curved feature): ESO 240-G01, NGC 7135 and NGC 5018. So, the Medusa morphology is not unique, but fairly rare, which could be an effect of the rather short timescales involved in the radial encounters. The Medusa still stands out among these objects because of its bright FIR emission and hot dust. Perhaps the Medusa is in a short-lived stage of intense star formation that will consume a large fraction of the molecular gas, leaving the galaxy to look more like NGC 7135 once the burst is over.

### 4.2 The starburst

The starburst of NGC 4194 is extended on a scale of $``$ 2 kpc (e.g. Wynn-Williams & Becklin 1993; Armus et al. 1990) where also 67% of the total 12CO emission can be found. We suggest that overall solid body rotation of the inner gas (see Fig. 4) provided a dynamical environment favourable to a burst of star formation (Keel 1993). At least 90% of the FIR luminosity can be attributed to a burst of star formation (Prestwich et al. 1994; Rowan-Robinson & Crawford 1989) and the high $`f`$(60$`\mu `$m)/$`f`$(100$`\mu `$m) flux ratio (Table 1) indicates that the starburst is heating the dust clouds to high temperatures.

Suchkov et al. 1993 argue that the correlation between CO and radio continuum distributions found in galaxies is an effect of cosmic ray heating of the molecular gas. They suggest that changing values in the ratio of CO and radio continuum intensity are indicators of the starburst age, with young starbursts having a higher value than old bursts since the cosmic rays are more confined to the molecular layer if the burst is young. In this model, the (low energy) cosmic rays originate in supernove remnants (SNRs). Their electron component is responsible for the radio continuum emission at cm wavelengths, while their proton component heats the gas. This provides a natural explanation for the rough correlation between continuum and molecular line emission – both are tied to star formation activity. This relation is, however, also expected if the (largely non-thermal) continuum emission simply traces SNRs and UV radiation from massive stars plays the major role in the heating of the dust and gas.

It is of course also possible to find molecular gas without a strong radio continuum counterpart if the gas is quiescent and star formation has not (yet) set in. The molecular peak $`e`$ in Fig. 10 may be an example for such a region.

Bonatto et al. (1999) compare NGC 4194 with the starburst merger NGC 3256. Their population synthesis results suggest that both galaxies have very similar stellar populations. Starbursts occured both at intermediate (1–2 Gyr) and young ages. NGC 3256 also has a kpc scale starburst (Graham et al. 1984) and is an HR merger with relatively faint 13CO 1–0 emission (e.g. Aalto et al. 1991a) and hot dust.

### 4.3 The central gas surface densities and the molecular line ratios

The 12CO surface brightness of NGC 4194 is lower by an order of magnitude compared to other HR mergers such as Arp 220 and Mrk 231. The average gas surface density for the extended starburst region of NGC 4194 is 500 — 1000 M pc-2, i.e. more than an order of magnitude below the surface density found in ULIRGs for the same CO–H2 conversion factor (e.g. Bryant & Scoville 1999). Addressing the question of the reason for the high $``$ posed in the introduction, we therefore conclude that the turbulent, extreme high pressure scenario for compact nuclear starburst regions (Aalto et al. 1995) is not well suited for explaining the unusual molecular line ratio of NGC 4194. Instead, we adopt a similar explanation as for NGC 3256: high gas temperatures depopulating the lower levels coupled with a disrupted and fragmented ISM (Aalto et al. 1991). One way of testing this scenario is to measure the 12CO/13CO intensity ratio for higher rotational transitions, and to obtain high resolution images of 12CO and 13CO like those for the Arp 299 merger. We can also conclude that NGC 4194 is the exception to the general trend of increasing central concentration of the molecular cloud distribution with increasing $`f`$(60$`\mu `$m)/$`f`$(100$`\mu `$m) flux ratio (Bryant 1997). A study of the intrinsic properties of compact and kpc scale HR galaxies is underway (Aalto et al., in preparation).

The star formation rate (SFR) in NGC 4194 has been estimated to be $``$40 M yr-1 (Storchi-Bergmann et al. 1994), based on H$`\alpha `$ emission. This is not much lower than the estimate for the HCN-bright ULIRG Arp 220, $``$110 M yr-1 for Arp 220 (Smith et al. 1998). Assuming a similar conversion factor for both galaxies, Arp 220 has about an order of magnitude more molecular gas in its central region. This would imply a considerably higher star formation efficiency in NGC 4194. However, the SFR for Arp 220 was determined based on its FIR luminosity (see Scoville & Soifer 1991 for the method). Thus, the estimates are done at different wavelenghts and should be viewed with some caution. Indeed, if the SFR of NGC 4194 is determined using the FIR luminosity, the result is much lower (6 – 7 M yr-1), resulting in star formation efficiencies for NGC 4194 and Arp 220 that are comparable.

The ratio of LFIR/LCO ratio is indeed similar in both galaxies implying that both starbursts are equally efficient in producing IR photons per unit molecular mass (again assuming that CO traces mass in roughly the same way in both objects).

Even if the star formation efficiencies and the LFIR/LCO are similar, there may be significant differences in the star formation properties of moderate luminosity mergers and ULIRGs, since the ISM in the two types of galaxies has a different structure. The starburst of NGC 4194 takes place in a low pressure environment compared to the starbursts at the heart of ULIRGs. The lower gas pressure results in a lower mass fraction of dense clumps and may allow destructive UV photons to penetrate further into the molecular clouds, where they can heat the dust and dissociate the molecules. Then, the young, massive stars in NGC 4194 would have a larger impact on the surrounding gas, disrupting it more effciently. This may also account for the relatively high brightness in H$`\alpha `$ and the hot dust. NGC 4194 has near solar metallicity, so the apparently efficient heating of the dust cannot be explained by a gas-to-dust ratio that is unusually low (e.g. Bonatto et al. 1999; Storchi-Bergmann et al. 1994).

### 4.4 The Fate of the Molecular Gas

The elevated 12CO/13CO 1–0 molecular line ratio together with the IR properties indicate that the molecular gas is being strongly affected by the newborn stellar population. The estimated star formation rate of NGC 4194 is about 40 M yr-1 which, if it continues, will have consumed the available gas in the center within 40 million years. The unusual molecular cloud properties may well be an indication that the standard conversion factor is inapplicable, so the real gas consumption time scale may be lower. The cloud properties will likely be quite similar to those in the NGC 3256 model (Aalto et al. 1991a) where the kinetic gas temperatures are found to be close to 100 K and the clouds are quite small, a few pc in size. Such a molecular ISM will radiate more brightly per unit mass than the cool disk clouds for which the conversion factor is calibrated. Experience from the NGC 3256 model leads us to suggest that the real molecular mass in the starburst region may be a factor of $``$5 lower than the mass based on the ‘standard’ conversion factor.

If all of the gas is to be consumed, also the gas at higher galactic radii needs to be transported to the center. The morphology of the central dust lane with respect to the 12CO peak is very suggestive of gas transport to the center. A position-velocity cut (Fig. 12) along the central dust lane, from $`e`$ to $`a`$, shows that the gas on the dust lane may be moving continuously into the 12CO peak $`a`$. On both sides of $`a`$ along the dust lane the velocities are lower by 50 $`\mathrm{km}\mathrm{s}^1`$ than in $`a`$, which could imply that gas is flowing into $`a`$. It is of course also possible that $`a`$ does not lie on the dust lane, but further to the core of the galaxy, in which case this may simply be a projection effect. Then again, the fact that $`a`$ is offset from the kinematic center argues against this last possibility. There are interesting morphological features on the dust lane just east and west of $`a`$ where it seems to flare out into a triangular shape, and darken (brighten in the negative Figs.). From Fig. 12 it is also clear that at the position of $`b`$ the emission splits into a brighter part at higher velocities and fainter emission at lower velocities.

The final appearance of NGC 4194 may well be that of an elliptical with shells and a remaining minor axis dust lane.

## 5 Summary

We have obtained a high resolution map of the Medusa merger, NGC 4194, in the 12CO 1–0 line. The main conclusions we draw from this map and single dish 12CO, 13CO and HCN observations are as follows:

1. The molecular distribution in NGC 4194 is extended over almost 5 kpc. Thus, this moderate luminosity merger shows clear morphological differences to typical ULIRGs, where an FIR luminosity higher by an order of magnitude usually goes along with a very compact molecular distribution on a scale of $`\genfrac{}{}{0pt}{}{_<}{^{}}`$ 1kpc.
2. The structure revealed by our map is complex, consisting of a rotating central concentration, low level emission extending north almost into the base of the tidal tail and very prominent arcs and filaments that can be identified with two dust lanes visible in the optical.
3. Along the central dust lane, we can identify several clumps that may be Giant Molecular Associations (GMAs). Some of these GMAs are positionally coincident with radio continuum features. However the strongest molecular concentration visible in 12CO is offset from the radio continuum peak. This peak may represent a second pre-merger nucleus in the system or – more likely – a region of exceptionally strong gas concentration and star formation.
4. NGC 4194 is clearly detected in a single dish 13CO spectrum. Thus, it has a high but not extreme 12CO/13CO line intensity ratio of $`20`$. The extent of the molecular distribution and the relative central surface density of the gas indicate that the physical mechanism causing the high $``$ is different in NGC 4194 from ULIRGs with high values of $``$. We suggest that a scenario with small clouds having high kinetic temperatures may account for the high $``$ found in NGC 4194, but cannot exclude that diffuse gas may also play a role.
5. The non-detection of the high density tracer HCN indicates that the amount of dense gas in NGC 4194 is lower than in ULIRGs. However, the star formation efficiency may be nearly the same, with a gas consumption time of less than 40 million years. Alternatively, massive young stars may have a higher impact on the fragmented ISM we suggest for NGC 4194.
6. It is likely that gas is at present flowing along the dust lanes toward the nucleus. However, the extended molecular distribution is at odds with the classification of NGC 4194 as an advanced major merger. The two properties can only be reconciled if the gas flow to the center of the system was inefficient compared to ULIRGs. Alternatively, NGC 4194 can be interpreted as an emerging shell galaxy resulting from a merger beween an early type and a spiral galaxy. Comparing simulations of such encounters with the morphology of NGC 4194 lead us to favour this latter interpretation of the Medusa merger history.

To further our understanding of the merger history of NGC 4194, the system needs to be mapped in more molecular transitions as well as H I and H$`\alpha `$. It represents a class of little studied moderate luminosity mergers. More generally, knowledge about the parameters that govern the extent and efficiency of the starburst resulting from a merger and its connection to the progenitors can only be gained if the moderate luminosity systems are investigated along with the high luminosity ones. If galaxy formation at high $`z`$ happened through mergers, it is essential to understand the differences between different types of mergers in the local universe, where the systems can be resolved with present techniques.

###### Acknowledgements.

Many thanks to Peter Bryant for his help in initiating this study. We thank J. Mazzarella for generously sharing his R-band image of NGC 4194 with us. We are thankful for the useful suggestions by an anonymous referee. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The OVRO mm-array is supported in part by NSF grant AST 9314079, and by the K.T. and E.L. Norris Foundation.


format: text

text:
# Coulomb exchange and pairing contributions in nuclear Hartree–Fock–Bogoliubov calculations with the Gogny force

## 1 Introduction

The mean field approximation is the backbone of many-body calculations in Nuclear Physics, either as a zero order approach to the problem or as a basis for more elaborated theories like the Random-Phase approximation or the Generator Coordinate Method. The modern and fast computer facilities have made possible the use of effective density dependent interactions, as the Skyrme  or the Gogny force , in standard mean field calculations like the Hartree-Fock (HF) or the Hartree-Fock-Bogoliubov (HFB) approaches.

In spite of the relative simplicity of the mean field approach, some additional approximations have been commonly used when effective forces were used in the mean field calculations in order to make the problem more tractable. The complications in the calculations usually arise from the exchange terms, either from the Coulomb force or from some components of the nuclear force itself. Depending on the mean field approach (HF or HFB) and on the force (zero or finite range) different approximations have been used in the past. In the HF theory one has only one exchange term (the Fock term) which contributes to the HF field. Commonly, in this type of calculations, in order to retain the required simplicity of the model, the Fock term of the Coulomb force has been neglected, or treated in the Slater approximation. In ref.  the validity of different approaches for the Coulomb Fock contribution was tested in HF calculations finding the Slater approximation a good one. In the HFB approach, and for any two-body interaction, one has two exchange terms, the Fock term which contributes to the HF field and the particle-particle (pp) term which contributes to the pairing field. In the HFB approach one must distinguish between the Skyrme and the Gogny force. Since the Skyrme force is assumed to provide only the particle-hole (ph) part of the force, all terms of the pp type are neglected even those arising from the Coulomb force. With the Gogny force the situation is different because with its finite range it has been designed to provide both, the ph and the pp part of the nuclear force.

A test similar to the one performed for HF in ref.  has never been done for HFB calculations. Exact HFB calculations with the Gogny force have been performed for spherical nuclei , but the effects of the approximation of neglecting the exchange terms never have been analyzed. For the more complex triaxial HFB calculations it has been assumed that the findings of ref.  for the HF approach, concerning the exchange terms, will also apply to this case. This point, however, never has been checked.

HFB calculations with triaxial symmetry using density–dependent interactions like the Gogny force, used by different groups , are usually performed using the following approximations:

* First, the contributions to the pairing field stemming from the center of mass correction, Coulomb and spin–orbit terms are neglected. This is motivated by the fact that their contribution, in the limit of axial deformations, is assumed (but not proved) to be negligible small . Thus, under these conditions the pairing field arises only from the Brink–Boeker term.
* Second, since the exchange (Fock) Coulomb term contribution to the HF field requires a large CPU time, its contribution to the energy is not taken into account or it is included in the Slater approximation.

On the other hand in the HFB theory, pairing is an important degree of freedom and it is important to analyze how the properties related to this part of the field change when some approximations are assumed. In this paper, for the first time, we calculate in the triaxial basis the contributions mentioned above . We shall investigate the effect of including all the terms in the fields as compared with the calculations which are normally done for deformed nuclei. Mainly the D1S parameterization of the Gogny force  is used in the calculations, though some calculations with the D1  parameterization will be done. We shall study nuclei which have been analyzed by the approximate triaxial HFB, like some nuclei in the rare earth region , in the superdeformed A$`190`$ region  and in the actinide region . In the same way, we shall also analyze these contributions in spherical nuclei in order to understand how the results obtained depend on the basis used. A brief review of the method is given in Sect. 2. The effect of the different contributions in spherical nuclei is studied in Sect. 3.1 for the N$`=50`$ region and in Sect. 3.2 for the N$`=82`$ y N$`=126`$ regions. The rare earth region will be analyzed in the Sect. 4.1. The effect of the new terms in the nucleus 240Pu (actinide region) will be treated in Sect. 4.2. The behavior of a nucleus in the superdeformed region $`A190`$ will be studied in Sect. 5. Finally an study of energy surfaces with considerations of the exchange terms will be done in Sec. 6. A discussion and the conclusions will be presented at the end of the paper.

## 2 Theory

The HFB theory  unifies the self–consistent description of nuclear orbitals, as given by the Hartree–Fock (HF) approach and the Bardeen–Cooper–Schrieffer (BCS) pairing theory into a single variational theory. As trial wave function an independent quasiparticle state $`|\mathrm{\Phi }`$ is considered. This wave function is a linear combination of independent multi-particle states representing various possibilities of occupying single–particle states. The quasiparticle operators are defined by

$$\alpha _\mu =\underset{i}{}U_{i\mu }^{}c_i+V_{i\mu }^{}c_i^{},$$
(1)
with $`c_k^{},c_k`$ the particle creation and annihilation operators in the harmonic oscillator basis and $`U`$ and $`V`$ the Bogoliubov wave functions to be determined by the Ritz variational principle. The corresponding equations have been solved by the Conjugate Gradient Method . The Bogoliubov transformation (1) leads to wave functions $`|\mathrm{\Phi }`$ that are not eigenstates of the particle number operator. Furthermore, to generate a wave function $`|\mathrm{\Phi }_I`$ corresponding to a rotational state with angular momentum $`I`$ we shall use the cranking prescription. Therefore, in the minimization process one has to add Lagrange multipliers to satisfy the pertinent constraints. That means, the wave function $`|\mathrm{\Phi }_I`$ is determined by the condition

$$\delta \mathrm{\Phi }_I|\widehat{H}\omega \widehat{J}_x\lambda _N\widehat{N}\lambda _Z\widehat{Z}|\mathrm{\Phi }_I=0$$
(2)
with the Lagrange multipliers $`\omega ,\lambda _N`$ and $`\lambda _Z`$ determined by the usual constraints $`\widehat{J}_x=\sqrt{I(I+1)}`$, $`\widehat{N}=N`$ and $`\widehat{Z}=Z`$, in an obvious notation.

In the second quantization formalism, the nuclear Hamiltonian is given by

$$\widehat{H}=\underset{k_1k_2}{}t_{k_1k_2}c_{k_1}^{}c_{k_2}+\frac{1}{2}\underset{k_1k_2k_3k_4}{}v_{k_1k_2k_3k_4}c_{k_1}^{}c_{k_2}^{}c_{k_4}c_{k_3}.$$
(3)
where $`t`$ is the kinetic energy and $`v`$ the two-body interaction, $`\overline{v}`$ if antisymmetrized. The Wick theorem allows the evaluation of the expectation value of a two-body operator in a simple way:

$$c_{k_1}^{}c_{k_2}^{}c_{k_4}c_{k_3}=c_{k_1}^{}c_{k_3}c_{k_2}^{}c_{k_4}c_{k_2}^{}c_{k_3}c_{k_1}^{}c_{k_4}+c_{k_1}^{}c_{k_2}^{}c_{k_4}c_{k_3},$$
(4)
where we have introduced the shorthand notation $`\mathrm{\Phi }|\widehat{O}|\mathrm{\Phi }\widehat{O}`$ for the expectation value of any operator $`\widehat{O}`$. The first term on the r.h.s. of eq. (4) is called the Hartree term, the second one the Fock term and the third one the pairing term. Using the expressions

$$\rho _{k_1,k_2}=c_{k_2}^{}c_{k_1},\kappa _{k_1,k_2}=c_{k_2}c_{k_1},$$
(5)
for the density matrix $`\rho `$ and the pairing tensor $`\kappa `$, and

$$\mathrm{\Gamma }_{k_1k_3}=\underset{k_2k_4}{}\overline{v}_{k_1k_2k_3k_4}\rho _{k_4k_2},\mathrm{\Delta }_{k_1k_2}=\frac{1}{2}\underset{k_3k_4}{}\overline{v}_{k_1k_2k_3k_4}\kappa _{k_3k_4},$$
(6)
for the HF field $`\mathrm{\Gamma }`$ and the pairing field $`\mathrm{\Delta }`$, we can express the expectation value of the nuclear Hamiltonian by

$$E=\widehat{H}=Tr(t\rho )+\frac{1}{2}Tr\left(\mathrm{\Gamma }\rho \right)\frac{1}{2}Tr\left(\mathrm{\Delta }\kappa ^{}\right).$$
(7)
The second addend in this expression arises from the contributions of the Hartree and the Fock terms of eq. (4) while the third one stems from the pairing term of eq. (4). The two-body effective nuclear interaction $`v`$ used in our calculations, see Appendix A, is made up from the Gogny force, the Coulomb (C) interaction and the two–body correction of the kinetic energy (TK). The Gogny force itself has the following terms: Brink–Boeker (BB), spin–orbit (SO) and density–dependent (DD) contributions. Then, the total HFB energy can be written as

$$E=Tr(t\rho )+V_{BB}+V_{SO}+V_{DD}+V_{TK}+V_C$$
(8)
where we have written separately the different energy contributions. Each contribution, $`V_L`$, is calculated from the corresponding Hartree–Fock, $`\mathrm{\Gamma }_L`$, and pairing, $`\mathrm{\Delta }_L`$, fields (L $``$ BB, SO, DD, TK or C) and is given by

$$V_L=\frac{1}{2}Tr\left(\mathrm{\Gamma }_L\rho \mathrm{\Delta }_L\kappa ^{}\right)=V_L^H+V_L^F+V_L^P$$
(9)
where $`\mathrm{\Gamma }_L`$ and $`\mathrm{\Delta }_L`$ are given by eq. (6) but considering, instead of the full interaction $`v`$, only the part $`L`$ of the force. On the r.h.s. of this expression we have further separated $`V_L`$ into the three contributions obtained from each of the terms of eq. (4), namely, the Hartree-, the Fock- and the pairing term.

As mentioned above, in HFB calculations with the Gogny force and a triaxial basis several approximations have been done. In particular the contributions to the pairing field from the spin–orbit, two–body correction of the kinetic energy and Coulomb terms have been usually neglected in the self–consistent procedure. Furthermore, the exchange contribution of the Coulomb potential (Fock term) to the Hartree–Fock field has not been taken into account in most of the cases or treated in the Slater approximation. Thus, in the simplest HFB approximation, which we shall denote HFBs ($`s`$ for standard), the HFB energy used in the variational process is given by

$$E_s=Tr(t\rho )+V_{BB}+V_{DD}+\frac{1}{2}Tr\left((\mathrm{\Gamma }_{SO}+\mathrm{\Gamma }_{TK}+\mathrm{\Gamma }_C^H)\rho \right)$$
(10)
A first order perturbation theory correction to this approach, which we shall denote HFBs+, is to calculate the neglected terms just at the end of the iterative procedure and to add them to the total energy. We shall test, in different mass regions, how good this approach is with respect to the exact calculation (7), denoted HFBe ($`e`$ for exact), and with respect to other approaches that we shall introduce later on. Details on the calculation of the terms neglected in the HFBs approximation but included in the exact one are given in Appendix B.

The triaxial basis used in the calculations is spanned by harmonic oscillator states with quantum numbers $`\{n_x,n_y,n_z\}`$ which fulfill the condition, $`a_xn_x+a_yn_y+a_zn_zN_0`$. The coefficients $`a_i`$ are related to the axis of the nuclear matter distribution by $`a_x=(qp)^{1/3}`$, $`a_y=q^{1/3}p^{2/3}`$ and $`a_z=p^{1/3}q^{2/3}`$ .

## 3 Spherical nuclei

In this section we concentrate mainly on the effect of the Coulomb contribution in HFB calculations for spherical nuclei. We are interested to know the behavior of the pairing and total energies when different approximations for the Coulomb terms are used. At the end of sect. 3.1, a short discussion is devoted to the contributions to the total energy of the pairing field of both, the spin–orbit term and the two–body correction of the kinetic energy term.

### 3.1 N$`=50`$ region

In this subsection we shall study the medium heavy nuclei from the region $`N=50`$. We shall use the basis determined by $`q=1.0`$, $`p=1.0`$ and $`N_0=11.1`$. In the calculations the D1S parameterization of the Gogny force has been used. First we analyze in detail the Coulomb terms and ignore for the moment the SO and TK pairing exchange terms.

The Coulomb interaction contributes to the energy with three terms : The Hartree (H / h) or direct term, the Fock (F / f) or exchange term and the pairing (P / p) or Bogoliubov term. Each term can be considered either selfconsistently, i.e., taken into account in the variational process, or non-selfconsistently (perturbationally), i.e., ignored during the variation but its value being added to the total energy at the end of the variational process. In the first case we shall use capital letters (H, F or P) and in the second one small letters (h, f or p) to label the approximation used for the respective terms. For example the approximation $`Hfp`$ means that the Hartree Coulomb term has been taken into account selfconsistently and the Fock and Pairing Coulomb terms only perturbationally. To understand which contributions are more important and the convenience of considering in the variational process the pertinent terms we shall investigate the following approximations:

In the simplest calculation the three Coulomb terms mentioned above, i.e., even the Coulomb Hartree term, are neglected. In this approximation the energy, which we shall call reference energy, $`E_{ref}`$, given by

$$E_{ref}=Tr(t\rho )+V_{BB}+V_{DD}+\frac{1}{2}Tr\left((\mathrm{\Gamma }_{SO}+\mathrm{\Gamma }_{TK})\rho \right)$$
(11)
is considered in the minimization process. The next approximations will be

$`E_{hfp}`$ $`=`$ $`E_{ref}+[V_C]`$ (12)
$`E_{Hfp}`$ $`=`$ $`E_{ref}+V_C^H+[V_C^F+V_C^P]`$ (13)
$`E_{HFp}`$ $`=`$ $`E_{ref}+V_C^H+V_C^F+[V_C^P]`$ (14)
$`E_{HFP}`$ $`=`$ $`E_{ref}+V_C`$ (15)
We use squared brackets to indicate the terms considered perturbationally, the others, are used in the variational equations. Notice, therefore, that only the latter ones determine the wave function. To remember the meaning of each approximation one just needs to look at the subindex of the energy and to keep in mind the convention introduced above about capital and small letters. We shall also investigate the Slater approximation used by several groups. In this case the energy looks like

$$E_{H\stackrel{~}{F}p}=E_{ref}+V_C^H+\stackrel{~}{V}_C^F+[V_C^P]$$
(16)
In $`E_{H\stackrel{~}{F}p}`$ the tilde on $`F`$ means that the Fock term has been calculated in the Slater approximation but it has been considered in the variational process. In these calculations the best one is, obviously, $`HFP`$ where all Coulomb terms are taken into account selfconsistently.

In the upper part of Table 1 we present the binding energies of five spherical nuclei with $`N=50`$ calculated in the different approximations. The experimental binding energies are from ref. . As expected from a variational method, the more terms are treated selfconsistently the deeper gets the energy minimum. The largest energy decrease is provided by the inclusion of the Hartree term in the variational equations. The largest deviation of the $`HFP`$ energies from the experimental values is about 2 MeV, but already the selfconsistent treatment of the Hartree term provides a good approximation to it. We also find that the Slater approximation provides energies close to the $`Hfp`$ approximation. These effects can be more clearly seen in Fig. 1a where we plot the quantity $`\delta E_{app.}`$ which represents the percentage of error of a given approximation with respect to the calculation $`HFP`$ and is defined by

$$\delta E_{app.}=100\times \frac{E_{app.}E_{HFP}}{E_{HFP}}$$
(17)
where $`E_{app.}`$ represents any of the approximations that we have considered. Here, we see that the $`hfb`$ approximation, where none of the Coulomb terms is treated selfconsistently, deviates the most from the exact $`HFP`$ results and that the deviation increases with the proton number. The other approaches provide a good approximation to the complete $`HFP`$.

In the second section of Table 1 we show the proton pairing energies $`E_{app.}^{pp}`$ for the ground state of the same nuclei. The neutron pairing is zero since we have a major shell closure for these nuclei. These energies are defined by the expression $`E_{app.}^{pp}=\frac{1}{2}Tr\left(\mathrm{\Delta }\kappa ^{}\right)`$. One has to notice that from the three Coulomb terms only the pairing term gives a contribution to $`E_{app.}^{pp}`$, the Hartree and Fock terms only give an indirect contribution through the changes introduced in the wave function by their consideration in the variational equations. We have also included in the table $`E_{ref}^{pp}`$, which gives us the pairing energy without any Coulomb contribution. The difference between $`E_{ref}^{pp}`$ and $`E_{hfp}^{pp}`$ gives us an indication of the size of the pairing Coulomb energy, which is about 1 MeV for all nuclei. When the direct Coulomb term is included in the HFB calculation, $`Hfp`$, we observe a slight decrease ( ranging from $`0.150`$ MeV to $`0.670`$ MeV) in the pairing energies with respect to the calculation $`hfp`$. The inclusion of the Fock Coulomb term, $`HFp`$, produces the opposite effect, we get a slight increase (up to 1.0 MeV) of the pairing energy, again with respect to $`hfp`$. Finally, treating all the contributions of the Coulomb force in a self–consistent way, $`HFP`$, we observe an important increase (in absolute value, a decrease) in the values of the proton pairing energy of up to $`4`$ MeV for the nucleus 90Zr. It is interesting to notice that though, for this nucleus, the pairing Coulomb term itself is about 650 keV large, when considered selfconsistently it grows up to 4 MeV. Concerning the Slater approximation, $`H\stackrel{~}{F}p`$, we find again that it provides results very close to the $`Hfp`$ ones. We observe that the maximum effect of the three contributions occurs for $`Z=40`$ ( semi-magic nucleus). This can be easily seen in Fig. 1b where a definition similar to eq. (17) has been used for the pairing energy. Comparing a) and b), we find a clear difference: the percentage of error in the pairing energy is much bigger than in the binding energy. We obtain that even the $`Hfp`$ approximation gives us a good result for the total energy, while this approximation is not good for the pairing energy. As a matter of fact, an important increase ( $`1.5`$ MeV for all nuclei) is obtained for the pairing energy when we compare $`HFP`$ with all other approaches considered.

In the third section of Table 1 we present the lowest two-quasiproton energies calculated in different approximations. In this case we have also included the exact solution, HFBe. In line with the decrease of the pairing energies just discussed, we observe a decrease of the excitation energies as more terms are included in the variational procedure. This decrease is particularly strong for the pairing exchange Coulomb term. We may conclude from this discussion that neglecting the Fock and pairing Coulomb terms in the variational equations does not affect significantly the total energy but it does influence the pairing properties of the wave function. Then, if we are interested in properties of spherical nuclei for which pairing correlations play an important role, like excited states or transition energies, one must be careful with the Coulomb exchange terms. Notice that not only the two-quasiparticle excited states will be affected but also any calculation based on a HFB solution, like the Random Phase Approximation, will also experience a lowering in the excitation energies of the low-lying collective states. We can see in Table 1 that the most important reduction in the pairing energy is due to the introduction of the Coulomb pairing term in a self–consistent way, though the contribution of this term itself is not very important, about $`1`$ MeV for the considered nuclei.

From the results just discussed one could conclude that the Coulomb interaction produces an antipairing effect in nuclei. In order to understand this effect we shall now analyze in a closer detail the behavior of the most relevant Coulomb term concerning pairing correlations, namely the pairing Coulomb term. We have done an additional test for the nucleus 90Zr. First of all, an exact HFB calculation (including the three Coulomb terms and the spin–orbit and two–body kinetic mass correction contributions to the pairing field, i.e., HFBe) has been done for this nucleus. Then, we do a progressive switching off of the Coulomb pairing contribution by a parameter $`\eta `$. For example, $`\eta =1.0`$ means that we take into account the full Coulomb pairing selfconsistently, $`\eta =0.6`$ means that 60% of it is taken selfconsistently and 40% perturbationally and so on. In fig. 2a we display the proton–proton correlation energy versus the $`\eta `$ parameter for the nucleus 90Zr. We find a linear behavior, which means that the change in the pairing energy due to the switching off of the Coulomb pairing contribution is produced in a continuous way. The solution that we obtain for $`\eta =1.0`$ is very similar to the solution for $`\eta =0`$ except for a decrease in the proton–proton pairing energy.

In fig. 2b) we present the proton single–particle energies, $`\epsilon (nlj)`$, outside the Z$`=28`$ shell closure in the canonical basis, and in fig. 2c),d) the occupation $`\nu (nlj)`$ for each level $`(nlj)`$, which is given by

$$\nu (nlj)=\frac{1}{2j+1}\underset{m}{}\mathrm{\Phi }|c_{nljm}^{}c_{nljm}|\mathrm{\Phi }$$
(18)
In this figure we observe that the single–particle energies $`\epsilon `$ barely change with $`\eta `$. With respect to the occupation numbers, we observe that the Coulomb pairing contribution goes in the direction to fill up the levels $`0f_{5/2}`$, $`1p_{3/2}`$ and $`1p_{1/2}`$ and to empty the level $`0g_{9/2}`$. That means, the Coulomb pairing term tries to increase the occupancy of the levels below the Fermi level at the cost of decreasing the occupancy of the levels above.

In the forth section of Table 1 we show the two–proton separation energies $`S^{2p}(N,Z)=E(N,Z2)E(N,Z)`$. As we can see in the table we obtain similar results in all approximations and the agreement with the experiment is good in all cases. We can understand this result if we consider that the binding energies are well described in all approximations and that that $`S^{2p}`$ is not an observable related with the pairing properties of the nuclei, which are really sensitive to the different approximations.

Finally, to complete the study of the neglected terms we shall now investigate the contributions of the spin-orbit term and of the two-body correction of the kinetic energy to the pairing field. In the fifth section of Table 1 we show the binding energies $`E_{so,tk}`$ obtained by adding these terms at the end of the minimization of eq. (15), i.e.,

$$E_{so,tk}=E_{HFP}+[V_{SO}^P+V_{TK}^P]$$
(19)
in a similar notation as introduced above. These terms are repulsive and a few hundred keV large. If these two terms are now taken into account selfconsistently we obtain the exact solution, i.e., HFBe, of the full HFB equations, eq. (8). The resulting energies are on the average 100 keV smaller than if the terms are considered perturbationally. The last two entries of the table are the pairing energies $`E_{so,tk}^{pp}`$, defined as before. We observe that these terms provide an additional quenching of the pairing correlations. The behavior is similar to the one observed with the other terms, i.e., the largest quenching takes place for the semi-magical nucleus 90Zr. For this particular nucleus we have calculated the effect of taking into account selfconsistently only one term and the other one perturbationally, i.e., we have performed the following calculations

$$E_{SO,tk}=E_{HFP}+V_{SO}^P+[V_{TK}^P]$$
(20)
$$E_{so,TK}=E_{HFP}+V_{TK}^P+[V_{SO}^P].$$
(21)
We obtain $`E_{SO,tk}^{pp}=1.671`$ MeV, $`E_{SO,tk}=784.360`$ MeV and $`E_{so,TK}^{pp}=1.355`$ MeV, $`E_{so,TK}=784.354`$ MeV, here again capital letters mean selfconsistent and small letters perturbational approaches. As we can see each term provides about half of the total effect.

### 3.2 The N$`=82`$ and N$`=126`$ regions

In this section we focus on the chains of the N$`=82`$ and N$`=126`$ isotones, again with neutron shell closure. Then, only proton pairing correlations are going to be important in this region. Two approximations are used in this case: the approximation HFBs, see eq. (10), which neglects all the terms mentioned in the introduction, and the exact calculation HFBe, see eq. (8), which takes into account all the contributions from the different terms of the force to the fields. In the HFBs approximation, the missing terms are not added at the end of the variational procedure. In order to be able to compare the total binding energies in the HFBs approach with the exact ones one must add these terms perturbativaly, i.e, the HFBs+ approximation. In the N$`=82`$ region, the basis is determined by $`q=1.0`$, $`p=1.0`$ and $`N_0=11.1`$ and in this case the parameterization D1 is used. In the N$`=126`$ region, the basis $`q=1.0`$, $`p=1.0`$ and $`N_0=12.1`$ has been used and for the force parameterization we have chosen D1S.

In the first two columns of Table 2 we display the total binding energies. We find that the approach HFBs+ provides a very good approximation to the exact results. In the last three columns of the same table we show the proton–proton pairing energies obtained in the approximations HFBs, HFBs+ and the exact ones HFBe. Both approximations provide more proton pairing correlations that the exact calculation HFBe. We can see that, contrary to the total binding energy case, the perturbational calculation HFBs+ is not a good approximation for the pairing energy. We do not find any noticeable difference between both mass regions or the Gogny parameterizations.

From these results we may conclude that, for spherical nuclei, the approximation of ignoring some terms of the hamiltonian, HFBs, presents some differences with the exact calculation. The approximation of considering these terms at least in first order perturbation theory, HFBs+, works well for observables related to, or that can be obtained from, total binding energies. On the other hand, the wave function content may be different from the one obtained in the exact calculation.

## 4 Normal deformed nuclei

In this section we will study the effects neglecting some terms in the hamiltonian for deformed nuclei. In particular we shall investigate their behavior at high angular momentum as well as the impact of the exchange terms on the neutron system.

### 4.1 Rare earth nuclei

In this region, we study six nuclei: the three Erbium isotopes 164,166,168Er, and the three Ytterbium isotopes, 162,164,166Yb. As in the preceding subsection we shall investigate only the approximation HFBs (or HFBs+ when the last one does not make sense) to compare with the exact calculation HFBe and the experimental results when available. For these calculations we use the basis with q$`=1.3`$, p$`=1.0`$ and N$`{}_{0}{}^{}=11.1`$ and the D1S parameterization of the force.

We shall first investigate the ground state properties of these nuclei. In the upper part of Table 3 we show the binding energies. We find that the approximate calculation, HFBs+, provides results very close to the exact ones. In the middle of the table we display the pairing energies calculated in both approaches. We find that the exact calculation provides proton pairing energies much smaller, in absolute value, than the approximate calculation, HFBs. The global effect that we observe in the HFBe calculations is a decrease of the absolute value of the proton pairing correlations compared to the approximate calculation. For 164Yb and 166Yb these correlations even vanish. The behavior of the neutron pairing energy is different depending on the nucleus: for the last two nuclei $`|E_e^{pp,\nu }|`$ is increased by about $`1`$ MeV and for the other nuclei it is decreased by approximately the same amount as compared to $`|E_s^{pp,\nu }|`$. In general, we can say that the global effect of including the Coulomb exchange and pairing contributions in the self–consistent procedure is a strong reduction of the absolute value of the proton pairing correlations in all the cases. The circumstance that the inclusion of the neglected terms influences more the proton pairing than the neutron one is obviously due to the Coulomb terms that do not affect the neutron system. With respect to the quadrupole moments we observe that the relative variation between both calculations is relative small, approximately $`0.5\%`$ for all the nuclei except for 164Yb where Q increases in the exact calculation by about $`10\%`$.

The fact that the total binding energies are very close in both approaches is the reason why the non self–consistent and perturbational approximation has been considered in the past as a good approximation in spite of the fact that the wave functions are different.

The proton pairing energies of the nuclei 164Yb and 166Yb show an anomalous behavior: For these nuclei, the proton pairing correlations vanish completely in the HFBe approach. Taking into account the behavior of the Er isotopes in Table 3, it would be expected to obtain for the nuclei 164Yb and 166Yb approximately $`4`$ MeV for the proton pairing energy. In order to understand this quenching we have performed a calculation for 164Yb similar to the one done for the nucleus 90Zr switching off progressively the Coulomb pairing contribution. In Fig. 3a) the proton and neutron pairing energies are displayed versus the $`\eta `$ parameter. As it can be seen from this figure, increasing the $`\eta `$ value up to $`\eta =0.8`$ produces a linear increase of the proton pairing energy. However, for $`\eta =0.9`$ a sharp increase is obtained, going to the unpaired solution. By extrapolation of the small $`\eta `$ values to $`\eta =1.0`$ the expected proton pairing energy for this value would be around $`2.5`$ MeV. However, the selfconsistent solution at $`\eta =1.0`$ corresponds to a non-paired proton system, which seems to indicate that the mean field approximation breaks down collapsing to the non-correlated regime. In Fig. 3b) we show the proton single particle energies in the canonical basis versus the $`\eta `$ parameter. Here we can observe the effect of the sharp phase transition on the single particle energies. In Fig. 3c) we can see how the occupations of the levels above (below) the Fermi level get empty (full). This sharp phase transition will explain the anomalous behavior of these two nuclei. The phase transition itself has more to do with the mean field approach than with the exchange terms. This can be seen very clearly in Fig.4a) where we show potential energy curves as a function of the proton fluctuation for different $`\eta `$ values. That means we have solved the HFB equations with an additional constraint on the w.f., namely the fluctuation in the proton number. For 164Yb and $`\eta =0`$ we find a well defined minimum around the superfluid solution, $`(\mathrm{\Delta }N)^2_\pi =5`$, i.e., the mean field approach makes sense. For $`\eta =0.9`$, however, the energy surface is quite flat indicating that a theory able to include further correlations would be more appropriate. For 162Yb, where no pairing collapse takes place, we find well defined minima for all $`\eta `$ values, see fig. 4b).

In order to investigate the high spin properties of these nuclei we have solved the cranking HFB equations, eq. (2), in both approximations for the same nuclei. In Fig. 5, we display the pairing energies, $`E^{pp}`$, as a function of the angular momentum. We find large differences between the approximate and the exact results. The absolute values of the HFBs proton pairing energies are much larger than the exact ones, the difference being not only quantitative but also qualitative for the nuclei 164Yb and 166Yb. These nuclei remain super-fluid for the whole spin range in the approximate solution and normal fluid in the exact one. For the neutron system, the behavior is qualitatively similar in both approaches, the small differences observed depend on the nucleus considered. The strong quenching of the pairing energies obtained in the exact calculations will affect some observables, for example, the transition energies along the yrast band. In Table 4 we show the E$`{}_{\gamma }{}^{}(2^+0^+)`$ values obtained in the approximations HFBs, HFBs+ and in the exact calculation. As it can be observed, to consider the exchange Coulomb term and all the contributions to the pairing field leads to a decrease in the transition energies, with a relative change around $`20\%`$. Since the perturbational calculation up to first order (HFBs+) provides a good approximation to the total binding energies we expect that the transition energies calculated in this way will also be close to the exact ones. As we can see in the table this is the case. On the other hand, by comparing our results with the experimental ones, which are included in the same table, we observe that the HFBe method gives a worse agreement with the experiment than the HFBs+ one.

It is interesting to know what happens with the D1 parameterization of the Gogny force. We have done the calculation with this parameterization for the nucleus 164Er, we obtain the value $`E_\gamma =95.7`$ keV with HFBs and $`E_\gamma =79.0`$ keV with HFBe, with a relative variation between both calculations of $`21.1\%`$, the same relative variation as obtained with the D1S parameterization, as we can see in the Table 4. The two values of the transition energy Eγ obtained with the D1 are closer to the experiment, but again the approximate calculation HFBs gives the best agreement with the experiment.

The behavior of the gamma ray energies at high spin is displayed in Fig. 6 for the six analyzed nuclei, together with the experimental data. We find that the results of the approximation HFBs are rather different from the exact ones as one expects from the disparity of the pairing energies found in the two calculations. For the nuclei 164Er, 168Er and 164Yb we have also included the approximation HFBs+. As expected the agreement with the exact calculation is much better than the one obtained in the HFBs calculation. We have also checked the Slater approximation for deformed nuclei and at high angular momentum. That means, we have evaluated the Fock exchange Coulomb term in the Slater approximation and added it to the energy $`E_s`$ of eq. (10) before the variation. The results of these calculations, quadrupole moments, pairing energies and gamma ray energies along the Yrast band, practically coincide with the results of the plain HFBs calculation.

From the results discussed up to now, it seems that it is enough to consider the neglected terms in first order perturbation theory. This, however is not quite true since some observables can not be expressed as a function of the binding energies. For instance, the transition probabilities, magnetic moments, nuclear radii, etc. It is interesting to investigate the predictions for some of these observables in the exact and in the approximated calculation. In Fig. 7, we show the beta deformations as a function of the angular momentum for the six nuclei under study. The beta deformation parameter is mainly affected by the alignment effects and the pairing correlations. In general, larger pairing correlations imply smaller deformations and larger alignments smaller $`\beta `$-deformations (Coriolis anti-stretching effect). In the exact calculations, in general, we have smaller absolute values of the pairing correlations and larger alignment. In principle, one could expect a compensation and only in cases where one of the effects is much larger than the other one we should obtain a change in deformation corresponding to the largest effect. In the figure we observe that in the nuclei 164Er, 168Er and 162Yb such a compensation takes place. In 166Er the compensation takes place for small spins but not for large ones and in 164Yb and 166Yb there is no compensation at all. Looking in Fig. 6, we find that in 166Er, there is clearly a stronger alignment in the exact solution than in the approximated one, which causes an smaller deformation in the exact calculation at high spins. In 164Yb and 166Yb the large quenching observed in the pairing energy explains the distinct behavior of the exact solution from the approximate one.

In Fig. 8 we present the reduced transition probabilities B(E2) for the nuclei we are analyzing. Though these probabilities are related with the $`\beta `$-deformations we expect larger differences between the two calculations in the B(E2) case because this quantity has to do only with protons while to the $`\beta `$-deformations both the proton and the neutron systems contribute. We should remember that neglecting the Coulomb exchange terms affects only the proton system. We find this supposition to be right, specially at high spins. For the nucleus 164Yb there are differences even at zero angular momentum due to the change of deformation in the ground state (see Table 3). Though we find some differences between both calculation, in the comparison with the experimental results it does not matter which one we take.

In Fig. 9 we finally present the gyromagnetic factors versus the angular momentum. For this observable we expect only small differences between the two calculations for the neutron gyromagnetic factor, $`g_n`$, for the same reasons as before, and somewhat larger for the proton gyromagnetic factor $`g_p`$. We find this supposition to be right. In the HFBs approximation, the general behavior of these nuclei is, first, a more or less sharp neutron (i13/2) alignment take place at medium spin values and then a proton (h11/2) alignment take place at high spins. In the HFBe calculation the situation is qualitatively similar with the exception of the nucleus 166Er which shows a very delayed neutron alignment, contrary to the experimental situation .

### 4.2 Actinide region: The nucleus 240Pu.

In this subsection we shall now discuss an example of a heavy deformed nucleus, to see if there is a possible mass dependence of the terms neglected in the approximate calculations. As a representative of the actinide region we have chosen the nucleus 240Pu. In the calculations we have used the basis q$`=1.3`$, p$`=1.0`$ and N$`{}_{0}{}^{}=14.0`$ which is large enough to guarantee the convergence, and the D1S parameterization of the force. Let us first discuss the wave function content of the HFBs approximation and of the exact theory HFBe. In Fig. 10, panel (a), we show the pairing energies $`E^{pp}`$ in both theories as a function of the angular momentum. As before, we find a strong quenching for the values of the proton system in the HFBe as compared to the HFBs one. The values for the neutron system, on the contrary, are very similar in both approaches. In panel (b) we display the electric quadrupole moment. At small angular momentum both predictions are very close but at large angular momentum they differ by up to $`2\%`$ at the largest spin considered. The fact that the HFBe prediction is smaller can be understood looking at panel (c) where we show the gyromagnetic factors. As expected the $`g_n`$ are very similar in both approaches for all spin values. The $`g_p`$, in the HFBe prediction indicates a larger proton alignment than in the HFBs one, which will cause the observed anti-stretching effect in the HFBe approach.

With respect to observables related to energy differences, we can consider the approaches HFBs, HFBs+ and HFBe. In the same figure, we present the transition energies $`E_\gamma `$, panel d), first moment of inertia $`𝒥^1`$, panel (e), and second moment of inertia $`𝒥^2`$, panel (f), in the three approaches, as well as the experimental data. The agreement of the plain HFBs with the experiment for the three observables is spectacular. The exact solution HFBe, according to the less pairing-correlated wave function, see panel (a), provides too large moments of inertia and too small gamma-ray energies. For the perturbative approach HFBs+, we obtain values similar to the exact calculation HFBe, though for the more sensitive quantity $`𝒥^2`$ some differences are found at the smallest and the largest angular momenta.

It is very remarkable that the HFBs+ binding- and $`\gamma `$-ray energies are very close to the HFBe ones. This fact also happens in the other nuclei calculated so far. To gain some insight in the reasons for this agreement we have calculated separately all contributions to the total energy in both approaches. That means, with the wave function $`|\mathrm{\Phi }_s`$ determined by minimization of $`E_s`$, we have evaluated the Fock and pairing term of the Coulomb force, as well as the pairing terms stemming from the spin-orbit and the two-body kinetic energy terms. We have also separately calculated the mentioned terms with the wave function $`|\mathrm{\Phi }_e`$ determined by minimization of the exact energy $`E_e`$. In the upper part of Table 5 we show the energy $`E_s`$ of Eq. (10) evaluated with the wave function $`|\mathrm{\Phi }_s`$ (approach $`s`$) and with $`|\mathrm{\Phi }_e`$ (approach $`e`$), as well as the terms just mentioned (we use the notation of Eq. (9)) for the ground state of 240Pu. The last column, denoted $`E_{tot}`$, corresponds to the total energy, $`E_{s+}`$ in the $`s`$ approach and $`E_e`$ in the $`e`$ one. As expected the binding energy is lower in the exact than in the perturbative calculation. On the other hand, since the wave function $`|\mathrm{\Phi }_s`$ has been determined by minimizing $`E_s`$, it is also obvious that this quantity is lower in the $`s`$ approximation than in the $`e`$ one. Concerning the other terms, each of them separately provides a bigger energy lowering in the $`e`$ than in the $`s`$ appoach. We also observe that the relative largest difference between both approaches corresponds to the pairing Coulomb term. In the lower part of the same Table we show, as a function of the angular momentum, the energy differences $`\mathrm{\Delta }\xi =\xi (I)\xi (I2)`$, where $`\xi `$ stands for each energy entry of the upper part of the Table, i.e., $`\mathrm{\Delta }\xi `$ represents the contribution of each term to the gamma ray energy from the state $`I`$ to $`I2`$. These energy differences have been calculated again in the approaches $`s`$ and $`e`$. Let us first concentrate on the $`I=2\mathrm{}`$ row. The prediction for the gamma ray energy taking into account only the terms of the standard approach are almost identical independently of the fact that we evaluate them with $`|\mathrm{\Phi }_s`$ or with $`|\mathrm{\Phi }_e`$ and the same happens with the other terms. Obviously, if one adds all these contributions, the prediction for the total $`\gamma `$-ray energy with all terms is very similar in both approaches. It is also interesting to note that the contributions from the terms not considered in the standard approach are all negative and lead to a more compressed spectrum. We also observe that the largest contribution by far is the one from the Coulomb pairing term, though the other contributions are not at all negligible. At higher angular momentum we observe the same trend as for $`I=2\mathrm{}`$, i.e., the energy differences for a given $`I`$ are very similar in both approaches. It seems therefore that the reason, for the good agreement of the gamma ray energies in both approaches, is that the dependence of the exchange terms from the angular momentum is almost identical in the HFBs+ and in the HFBe approaches. For $`I=0\mathrm{}`$ the binding energies differ only by about 500 keV. The exchange terms are rather different in both approaches, however, when the $`\gamma `$-ray are calculated, the differences cancel and one gets similar values in both approaches.

## 5 Superdeformed nuclei in the A $`190`$ region

This region of superdeformation has already been studied using the HFBs method with the Gogny force D1S . In our calculations we use the basis q$`=1.5`$, p$`=1.0`$ and N$`{}_{0}{}^{}=12.5`$ and the D1S interaction. Taking into account the results obtained in the other regions, we do not expect a big influence of the Coulomb exchange terms because for the superdeformed states the proton pairing energies either vanish or are very small even with the HFBs method. A large change of the pairing energies due to the Coulomb exchange and the other terms cannot exist, therefore, in this states. As an example of superdeformation we have chosen the nucleus 190Hg. From this example we will learn how the other pairing contributions (spin–orbit and two–body correction) affect the neutron pairing energies and the other nuclear properties through the selfconsistency. As before, we shall first look for the wave function content of the HFBs and the HFBe approaches. In Fig. 11, panel (a), we can see the behavior of the pairing energies in both methods. As mentioned, the proton pairing energies along the band vanish in both calculations. The neutron pairing energies are close to each other in both calculations, though somewhat smaller values are obtained for the HFBe at low angular momentum. The charge quadrupole moment is depicted in panel (b). Only small differences between both approaches are found at large angular momentum. The gyromagnetic factors, panel (c), are again very close to each other. We shall now turn to a comparison between the HFBs, HFBs+ and HFBe predictions for the $`\gamma `$-ray energies, panel (d), first moments of inertia, panel (e) and second moments of inertia, panel (f). For these observables we do not find remarkable differences between the three approaches.

From this section we may conclude that the Coulomb exchange and pairing contribution barely modify the properties of the nucleus when the approximate calculation HFBs has vanishing proton pairing energies. Then, the approach HFBs is a good approximation to the completely self–consistent calculation, HFBe, even for the wave function.

## 6  Energy Surfaces

In this section we shall finally investigate the effect of the exchange terms as a function of a collective variable at zero angular momentum. We have used an axially symmetric HFB code that allows reflexion asymmetric shapes. We have performed two different calculation. In the first one we have neglected all exchange terms, i.e., the HFBs approach, but we have included, the Fock Coulomb term in the Slater approximation, we shall call this approach HFBs,Sla. The second approach includes all exchange terms with the exception of the contribution of the spin-orbit term to the pairing field,111We do not include this term in the calculations just because it is not relevant for the problem we are considering, besides the fact that it is very cumbersome to include in the axial HFB code. we shall call this approach HFBe,nso. Notice that since we do not add at the end of the calculation HFBs,Sla the pairing terms of the Coulomb and the two-body kinetic energy (both contributions are repulsive) it is possible to find HFBs,Sla solutions with an energy deeper than the HFBe,nso approach. That means we are more interested in the general behavior of the energy surface than in the absolute values. The most common case of energy dependence as a function of a collective variable appears in the calculation of fission barriers. In this calculations we use the basis q$`=1.5`$, p$`=1.0`$ and N$`{}_{0}{}^{}=15.1`$ and the D1S parameterization. In fig. 12, on panel a), we display the binding energy of 254No as a function of the constrained quadrupole moment in both approaches. We find a similar behavior in both calculations, though small differences are noticeable. For example, the height of the first barrier is slightly larger in the HFBe,nso approach than in the HFBs,Sla. The second one, however, has more or less the same height in both approaches. In panel b) we display the pairing energies along the fission path. The neutron pairing energy is practically the same in both calculations. The proton pairing energies in the HFBe,nso calculations show the Coulomb antipairing effect already discussed in other nuclei. They are in absolute value 5 to 10 MeV smaller than in the calculation without Coulomb pairing term.

We shall now investigate the effect of the exchange terms in cases where the reflexion symmetry is broken i.e. the wave functions do not have a good parity. In these calculations we use the same theoretical approaches as in the 254No ones, we use the basis q$`=1.3`$, p$`=1.0`$ and N$`{}_{0}{}^{}=13.1`$ and the D1S interaction. In panel c) of the same figure, we present the binding energy of the nucleus 226Th versus the constrained octupole moment. As with the quadrupole operator we do not find any remarkable difference between both calculations. The depth and shape of the octupole deformed minimum are almost identical in both approaches. In panel d), finally, we present the dipole moment versus the constrained octupole momentum. Again, small differences between both calculations are found but without further relevance.

## 7 Discussion

The main outcome of our investigation is that the exchange terms usually neglected affect mainly the nuclear pairing properties. More specifically, they affect mostly the proton nuclear pairing properties because of all missing terms the largest contribution to the pairing field arises from the Coulomb potential. Observables little or not dependent on the pairing correlations, are, therefore, not affected by these terms. On the contrary, observables like the moments of inertia or two-quasiparticle states energies, strongly dependent on pairing correlations, are very sensitive to these terms.

We have found that the results of calculations performed with the plain HFBs approximation differ significantly, for some observables, from those of the exact solution HFBe, in nuclei where the proton pairing energies are relatively large. On the contrary, the HFBs+ results agree, in general, very well with those of the exact solution. In cases of very small proton correlation energies all three approaches, HFBs, HFBs+ and HFBe practically coincide.

For those situations where proton correlations are relevant, the predictions of the HFBs approximation agree better with the experimental data than the exact ones or than the ones from the HFBs+ approximation. Paradoxically, the Gogny force was fitted  using the exact HFB approach to calculate properties of spherical nuclei. The pairing properties of the Gogny force were adjusted by fitting the odd-even mass difference of the $`Sn`$ isotones. The $`Sn`$ isotones, however, have a major shell closure in protons, Z$`=50`$, and in this case, practically, it does not matter whether one performs approximate or exact calculations, i.e., effectively the Coulomb exchange terms were not considered in this fit. Furthermore, the fact that including the exchange terms in the calculations worsens the good agreement of the theoretical results with the experimental data, reinforces the conjecture that the actual fits of the Gogny force should be used, in calculations with open proton shells, without the mentioned terms.

Of course the Gogny force is isospin invariant but the Coulomb force is not. The Hartree-Fock potential that one obtains after solving the HFB equations neither is isospin invariant. It is not clear, therefore, that a different parameterization (concerning the pairing part of the interaction) of the force would not be obtained if nuclei with open proton shells had been considered in the fit. The whole issue of adjusting the pairing properties of any force is itself rather fuzzy. Most of the fits have been done to reproduce the experimentally observed odd-even mass difference. This energy difference can, however, not be attributed only to pairing correlations . From the theoretical point of view the calculation of this energy difference presents also problems because binding energies of odd nuclei are not easy to calculate. Blocking effects and angular momentum projection should be performed in order to evaluate properly this energy, besides additional couplings. However, in order to save CPU time during the fit, approximations are used, in particular in order not to break the time reversal symmetry the equal filling approximation is used.

Taking into account the complication of the neglected exchange terms, see Appendix B, the fact that the HFBs approximation, commonly used in most HFB calculations, provides such good agreement with the experimental data is very fortunate.

## 8 Conclusions

For the first time we have performed an exact self–consistent HFB calculation in a triaxial basis with the Gogny force. The exact Coulomb exchange term and the contributions of the spin–orbit, two–body correction of the kinetic energy and Coulomb terms to the pairing field have been included in the calculations. An exhaustive analysis of each term has been performed for spherical nuclei in different mass regions. We find the Coulomb contribution to the pairing field to be the most relevant one. The other terms, though not negligible, give smaller contributions.

For deformed nuclei we have studied the commonly used HFBs approximation, the perturbative HFBs+, and the exact solution. We have been concerned with different mass regions at zero and high spins. As in the case of spherical nuclei, we observe a reduction of the proton pairing energies in the exact calculations as compared to the HFBs ones. Ground state energies are very similar in the perturbative and in the exact calculations. Energies of excited states are somewhat different in the HFBs approximation than in the exact one. We have also investigated superdeformed states. Since the proton pairing correlations of these nuclei are small we do not find any significant difference between the approximate and the exact solution. Lastly, we have analyzed the effect of the most relevant exchange terms in the calculation of energy surfaces. The energy as a function of the quadrupole (octupole) mass operator is practically identical in the approximate and in the exact calculations.

We also have found, that the Slater approximation for the Fock Coulomb term is a good one for all HFB calculations, in agreement with ref. for the HF case.

For all nuclei and situations analyzed the following general conclusions apply:

1. In cases where the proton pairing correlations do not play an important role, anyone of the approximations considered practically coincide with the exact calculation.
2. When proton pairing correlations are relevant, then:

* For energy-related observables the perturbative approach HFBs+ is a good approximation to the exact one.
    * For those observables whose evaluation requires the wave functions, one should perform the exact calculation.
    * If the HFB wave functions are thought as a basis for more elaborated theories, like RPA or GCM, then the exact calculation is required.
3. The way in which the fits of the Gogny force were performed favors neglecting the mentioned exchange terms in HFB calculations using these fitted parameter sets.

It remains to investigate the real importance of this terms by performing new fits of the Gogny force taking into account explicitly all terms discussed here.

We thank J.F. Berger, M. Girod and S. Peru for discussions. This work was supported in part by DGICyT, Spain under Project PB97-0023. One of us (M.A) acknowledges a grant from the Spanish Ministery of Education under Project PB94-0164.

Appendix A: The interaction

As an effective interaction we use the Gogny force

$`v_{12}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{2}{}}}e^{(\stackrel{}{r}_1\stackrel{}{r}_2)^2/\mu _i^2}(W_i+B_iP_\sigma H_iP_\tau M_iP_\sigma P_\tau )+`$ (22)
$`+`$ $`W_{LS}(\stackrel{}{\sigma }_1+\stackrel{}{\sigma }_2)\stackrel{}{k}\times \delta (\stackrel{}{r}_1\stackrel{}{r}_2)\stackrel{}{k}+`$
$`+`$ $`t_3(1+x_0P_\sigma )\delta (\stackrel{}{r}_1\stackrel{}{r}_2)\rho ^{1/3}\left({\displaystyle \frac{1}{2}}(\stackrel{}{r}_1+\stackrel{}{r}_2)\right),`$
and the Coulomb force

$$v_{12}^C=(1+2\tau _{1z})(1+2\tau _{2z})\frac{e^2}{|\stackrel{}{r}_1\stackrel{}{r}_2|}.$$
(23)
Additional contributions taken into account in the calculations arise from the one–body and two–body center of mass corrections

$$\widehat{T}=\underset{i}{}\frac{\stackrel{}{p}_{i}^{}{}_{}{}^{\mathrm{\hspace{0.17em}2}}}{2m}\left(1\frac{1}{A}\right)\frac{1}{Am}\underset{i>j}{}\stackrel{}{p}_i\stackrel{}{p}_j.$$
(24)
In the calculations we use the parametrizations D1S  and D1 .

Appendix B: Calculation of the Coulomb Hartree-Fock field and pairing tensor.

To compute the Hartree-Fock field and pairing tensor for the Coulomb interaction we have followed the idea used in . It amounts to use the identity

$$\frac{1}{|\stackrel{}{r}_1\stackrel{}{r}_2|}=\frac{2}{\sqrt{\pi }}_0^{\mathrm{}}\frac{d\mu }{\mu ^2}e^{\left(\stackrel{}{r}_1\stackrel{}{r}_2\right)^2/\mu ^2},$$

to compute the Coulomb matrix elements in terms of those of Gaussians with range $`\mu `$. The main advantage of this approach is that the matrix elements of the Gaussians are factorizable in terms of quantities defined for each Cartesian direction $`x`$, $`y`$ and $`z`$. Unfortunately, the price to pay is the $`\mu `$ integration. Introducing the HF field $`\mathrm{\Gamma }_{k_1k_2}^G(\mu )`$ and the pairing tensor $`\mathrm{\Delta }_{k_1k_2}^G(\mu )`$, see eq.(6), of a Gaussian of range $`\mu `$ we can express the corresponding quantities for the Coulomb potential as

$$\mathrm{\Gamma }_{k_1k_2}^C=\frac{2e^2}{\sqrt{\pi }}_0^{\mathrm{}}\frac{d\mu }{\mu ^2}\mathrm{\Gamma }_{k_1k_2}^G(\mu )$$

and

$$\mathrm{\Delta }_{k_1k_2}^C=\frac{2e^2}{\sqrt{\pi }}_0^{\mathrm{}}\frac{d\mu }{\mu ^2}\mathrm{\Delta }_{k_1k_2}^G(\mu )$$

To transform the integration interval to a finite one, one makes the change of variables $`u^2=\left(1+\alpha \mu ^2\right)^1`$, with $`\alpha `$ a parameter to be determined below. We finally obtain

$$\mathrm{\Gamma }_{k_1k_2}^C=2e^2\sqrt{\frac{\alpha }{\pi }}_0^1\frac{du}{\left(1u^2\right)^{3/2}}\mathrm{\Gamma }_{k_1k_2}^G\left(\frac{\sqrt{\frac{1}{\alpha }\left(1u^2\right)}}{u}\right)$$

and

$$\mathrm{\Delta }_{k_1k_2}^C=2e^2\sqrt{\frac{\alpha }{\pi }}_0^1\frac{du}{\left(1u^2\right)^{3/2}}\mathrm{\Delta }_{k_1k_2}^G\left(\frac{\sqrt{\frac{1}{\alpha }\left(1u^2\right)}}{u}\right).$$

The integrations are carried out numerically with the Gauss-Legendre method. The optimal choice of $`\alpha `$, turns out to be equal to the maximum of $`(\alpha _x,\alpha _y,\alpha _z)`$ where $`\alpha _i=\frac{1}{2b_i^2}`$ and $`b_i`$ are the oscillator lengths (see ref. ).

In the numerical computation of the Coulomb exchange and pairing terms we have used 20 points for the Gauss-Legendre integration. This number of points has proven to be enough for a precise determination of the exchange field and pairing tensor. The number of integration points used makes the evaluation of the Coulomb terms rather costly in terms of CPU time: it usually takes a factor 5.5 longer to compute those factors than to calculate the Fock and pairing parts of the Brink-Boecker interaction (in this case, we have two gaussians and we have to compute them for both protons and neutrons). The calculation of the exchange and pairing parts of the BB interaction is already the most costly part of each iteration as it usually takes between a factor 10 and 15 longer than the evaluation of the direct terms. In a standard workstation and for a configuration space with parameters $`p=1`$, $`q=1.3`$ and $`N_0=11.1`$ the evaluation of the direct field takes 1.9 seconds, the combined Fock and pairing terms of the BB part takes 19 seconds and the evaluation of the Fock and pairing parts of the Coulomb interaction takes 103 seconds. One has to keep in mind that the previous numbers are for each iteration of the minimization process. If we increase the size of the basis by using $`q=1.6`$ and $`N_0=14.0`$ the previous numbers climb up to 7, 92 and 540 seconds respectively. We clearly see the tremendous slow down introduced in the calculations stemming from the consideration of the Fock and pairing terms of the Coulomb interaction. This tremendous slow down can be even worse for the practitioners of the Skyrme interaction as their CPU time consumption per iteration has to be of the order of the one needed to calculate the direct field of the Gogny force.

Concerning the computation of the pairing field associated to the spin-orbit interaction we have followed the ideas already presented in . An explicit expression for all these terms can be found in . The numerical computation of this field takes almost the same time as the computation of the HF field and therefore, its inclusion is not significant in terms of CPU time.