text:
# THE NATURE OF THE COMPACT SUPERNOVA REMNANTS IN STARBURST GALAXIES

## 1 INTRODUCTION

Radio observations of starburst galaxies have revealed a population of compact supernova remnants that are more luminous than Cas A, the most luminous supernova remnant in our Galaxy. This population has been best studied in M82 (Kronberg, Biermann, & Schwab 1985), and similar populations have been found in NGC 253 (Antonucci & Ulvestad 1988), NGC 3448 (Noreau & Kronberg 1987), and Arp 220 (Smith et al. 1998). Although there is general agreement that these nonthermal radio sources are the result of supernovae, there has been no clear sense of their evolutionary status. One possibility is that they involve interaction with the winds from the progenitor massive star. The radio supernova phenomenon can be interpreted in terms of supernova interaction with the free wind of the progenitor (Chevalier 1982; Weiler et al. 1986). In this case, the radio properties do not depend on the interstellar environment. Alternatively, the remnants are interacting with interstellar material, in which case the starburst environment is likely to be important.

Recent observations of the compact remnants in M82 have yielded clues about their nature. Kronberg et al. (2000) found that 18 of the sources (75% of the sample) showed no discernible flux evolution over a 12 year period, which suggests a flux variability $`<0.1`$% yr-1 (see also Ulvestad & Antonucci 1994). The evolution timescale is then $`>10^3|n|`$ yr if the luminosity evolves as $`Lt^n`$. An estimate of $`n`$ can be obtained from the radio $`\mathrm{\Sigma }D`$ (surface brightness – diameter) relation, assuming that it represents an evolutionary trend; for $`\mathrm{\Sigma }D^{3.5}`$ (Huang et al. 1994), we have $`|n|=1.5`$. Muxlow et al. (1994) resolved most the remnants in M82, finding diameters of $`0.64`$ pc, and suggested that they are supernova remnants rather than radio supernovae. In some cases, shell structure was observable, as expected for a supernova remnant. The combination of size and timescale yields a velocity $`10^3\mathrm{km}\mathrm{s}^1`$, which implies that these are evolved remnants. This forms the basis of our treatment of the compact remnants in M82 (§ 2) and in other galaxies (§ 3). We also briefly discuss supernova-driven superwinds in § 3. We take a distance to M82 of 3.63 Mpc (Freedman et al. 1994).

## 2 THE COMPACT REMNANTS IN M82

Of the 24 radio sources in M82 followed by Kronberg et al. (2000), two, 41.5+59.7 and 41.9+57.5, show especially rapid evolution. The first one was observed in 1981, but not in 1982. The rapid evolution was consistent with the radio emission from a Type Ib or Ic supernova (Kronberg et al. 2000), and we agree that this is a plausible interpretation. The emission is from a shock wave propagating in the free wind from the progenitor star. The other source, which has been the brightest of the compact sources, had an evolution that is compatible with that of a bright radio supernova like SN 1986J, but recent radio imaging shows that it has a bipolar structure (McDonald et al. 2001). We do not include 41.9+57.5 in our discussion.

The duration of the initial radio supernova phase is limited by the extent of the freely expanding progenitor wind. In the case of a starburst region, the wind expansion is likely to be limited by the high pressure of the interstellar medium. Based on energy input from supernovae, Chevalier & Clegg (1985) estimated a pressure in the central region of M82: $`p/k10^7`$ cm-3 K, where $`k`$ is Boltzmann’s constant. The ram pressure of a wind equals the external pressure at a radius $`r_w=0.2\dot{M}_4^{1/2}v_{w1}^{1/2}p_7^{1/2}`$ pc, where $`\dot{M}_4`$ is the mass loss rate in units of $`10^4M_{}\mathrm{yr}^1`$, $`v_{w1}`$ is the wind velocity in units of $`10\mathrm{km}\mathrm{s}^1`$, and $`p_7`$ is the external pressure in units of $`10^7`$ cm-3 K. Except for 41.9+57.5, the radio remnants have $`r0.5`$ pc (Muxlow et al. 1994), which is too large for interaction with a red supergiant wind. A Wolf-Rayet progenitor wind could expand to a larger radius because of the larger wind velocity ($`v_{w1}100`$), but the density is low, leading to a faint source. The majority of the compact sources in M82 cannot be interpreted as standard radio supernovae.

The interstellar medium in the center of M82 is believed to have a hot component that drives a galactic superwind (Chevalier & Clegg 1985; McCarthy, van Breugel, & Heckman 1987; Strickland & Stevens 2000). In the reference wind model of Chevalier & Clegg (1985) with a power of one $`10^{51}`$ erg supernova every $`35`$ years and a mass loss rate of $`1M_{}\mathrm{yr}^1`$, the central temperature is $`T_c5\times 10^8`$ K and central H density is $`n_H0.05`$ cm-3. The models of Strickland & Stevens (2000) are cooler, and appear to have $`n_H0.1`$ cm-3 in the central region. In Chandra observations, Griffiths et al. (2000) found an energetic X-ray component in the central kpc of M82 that is more luminous than predicted by the thermal emission in the wind models, but much of the emission could be from the inverse Compton mechanism. The radius at which a supernova has swept up its own mass is $`R_s=8.9(M/10M_{})^{1/3}(n_H/0.1\mathrm{cm}^3)^{1/3}`$ pc, where $`M`$ is the ejecta mass, and the radius at which a $`10^{51}`$ erg supernova comes into pressure equilibrium with its surroundings is $`R_p=16p_7^{1/3}`$ pc. Under these conditions, the observed radio remnants, with diameters $`4`$ pc (Muxlow et al. 1994), would be primarily in free expansion, with characteristic velocities $`5,000\mathrm{km}\mathrm{s}^1`$.

In view of the large evolutionary times for the compact remnants, we propose that most of the compact remnants are interacting with a dense, possibly molecular, interstellar medium. The analogy in our Galaxy would be remnants evolving in the interclump medium of molecular clouds (Chevalier 1999), which has a H density $`n_H10`$ cm-3. Galactic remnants like IC 443 and W44, which have high radio surface brightnesses and radii $`10`$ pc, can be interpreted as recently having entered the radiative phase of evolution. The remnants in M82 can also be interpreted as recently having entered the radiative phase if the surrounding density is $`n_H10^3`$ cm-3. This is a plausible density for the interclump medium of molecular clouds in M82 if the high pressure in the central region is taken into account. In our Galaxy, the pressure of the interclump medium of clouds is maintained at $`10^5`$ cm-3 K by nonthermal support, presumably related to magnetic fields. In M82, a similar situation is likely. Shen & Lo (1995) find a CO line width of $`7\mathrm{km}\mathrm{s}^1`$ and up in molecular clouds in M82. When this line width is combined with a density $`n_H10^3`$ cm-3, the result is a nonthermal pressure $`p/k10^7`$ cm-3 K, as is inferred in the central region of M82. In addition, there is direct evidence for densities $`10^3`$ H atoms cm-3 in the interclump medium of molecular clouds in starburst regions (Mao et al. 2000; Solomon 2001). The location of the radio remnants in M82 overlaps the regions of strong CO emission (e.g., Fig. 4 of Lord et al. 1996) and the gas column density to the remnants is high (Mattila & Meikle 2001), which is consistent with interaction with molecular gas.

A possible complication is that the analysis of the dense gas emission in terms of photodissociation regions implies that the gas has a filling factor $`0.050.1`$ and is in many clouds with radii $`1`$ pc (Lord et al. 1996; Mao et al. 2000). If this were the case, the observed radio emission is presumably from shock fronts in the dense gas, but some of the expansion is in a lower density medium. However, an LVG (large velocity gradient) analysis of the CO data yields cloud radii of $`150`$ pc (Mao et al. 2000) and Shen & Lo (1995) identify $`60`$ CO “clouds” with half-power diameters of $`10100`$ pc. In any case, the dense gas is likely to have complex structure. The radio images of the remnant do show considerable structure in the nonthermal emission (Muxlow et al. 1994).

An upper limit on the filling factor $`f`$ of $`10^3`$ cm-3 gas can be obtained by assuming that all of the central gas is at this density. Assuming a disk of diameter 600 pc, a total height of 100 pc, and an interclump density of $`n_H=10^3`$ cm-3, the gas mass is $`9.8\times 10^8fM_{}`$. The total gas mass within a radius of 400 pc, including H2, HI, and HII, is $`1.3\times 10^8M_{}`$ (Götz et al. 1990), yielding $`f0.13`$. The filling factor appears to be low, but we expect the supernovae to come from massive stars that are born in the molecular gas, so that the supernovae are correlated with the presence of molecular gas.

Although a late red supergiant phase is unlikely to significantly affect the surrounding medium, as discussed above, photoionization and a wind during the main sequence phase could have an effect on the surroundings. Draine & Woods (1991) considered these effects and found that for a surrounding density of $`10^3`$ cm-3, stars with an initial mass of $`25M_{}`$ or less have supernovae that enter the radiative blast wave phase in that medium. This should apply to most of the massive star supernovae in M82.

Wheeler, Mazurek, & Sivaramakrishnan (1980) find that a remnant in a medium with $`n_H=10^3`$ cm-3 becomes radiative at a radius $`R_c=0.78`$ pc for a $`10^{51}`$ erg explosion. For the same parameters, Draine & Woods (1991) find $`R_c=1.37`$ pc; the difference can be attributed to the fact that their cooling function is a factor $`2`$ smaller than that used by Wheeler et al. (1980). Except for the sources 41.9+57.5 and 44.0+59.6, which may be an active galactic nucleus (Seaquist, Frayer, & Frail 1997), Muxlow et al. (1994) find that the diameters of the compact remnants in M82 are in the range 1.1 to 4 pc. We suggest that the remnants in M82 have recently entered the radiative phase of evolution. The evidence in our Galaxy is that remnants which have recently become radiative (e.g., W28, W44) have higher radio surface brightnesses for a given diameter than those that are non-radiative with a similar diameter (e.g., SN 1006).

To estimate the observed X-ray fluxes, we use an analytical model for the evolution of the supernova remnants, as in Draine & Woods (1991), for different values of the ambient density, $`n_H`$. From the shock temperature and density we calculate the cooling rate and luminosity of the remnant as function of time. Using a simple bremsstrahlung form for the spectrum and including the absorption by the interstellar medium, we estimate the resulting X-ray luminosity in the 0.1 – 2 keV and 2 – 10 keV bands. In the upper panels of Fig. 1 we show the unabsorbed luminosities in the two bands, while the lower panels show the observed luminosities for a column density of $`4.3\times 10^{22}\mathrm{cm}^2`$ (solid lines), and for a lower column density of $`2\times 10^{22}\mathrm{cm}^2`$ (dashed lines). Mattila & Meikle (2001) find that the average H column density to the remnants, including atomic and molecular H, is $`4.3\times 10^{22}`$ cm-2 ($`\sigma 1.7\times 10^{22}`$ cm-2). We terminate the calculations at the point when the internal pressure of the remnant is equal to the external pressure of the interstellar medium, here taken to be $`p/k=10^7\mathrm{K}\mathrm{cm}^3`$. This mainly affects the cases with $`n_H10^2\mathrm{cm}^3`$. The peaks in the unabsorbed X-ray luminosity, especially apparent in the light curves of the 0.1 – 2 keV band, occur at the time of the transition from the adiabatic to the radiative phase.

It is seen that in the 0.1 – 2 keV band the observed luminosities are $`3\times 10^{37}\mathrm{erg}\mathrm{s}^1`$ for $`n_H=10^3\mathrm{cm}^3`$, and considerably below this for lower densities. For the low value of the column density the luminosity approaches $`1\times 10^{38}\mathrm{erg}\mathrm{s}^1`$. The observed luminosities in the 2 – 10 keV band are basically unaffected by the interstellar medium absorption. For $`n_H=10^3\mathrm{cm}^3`$, the peak luminosity is $`3\times 10^{38}\mathrm{erg}\mathrm{s}^1`$, scaling roughly with the ambient density.

A Chandra image of M82 (Griffiths et al. 2000) shows many compact X-ray sources, but their positions are $`>1^{\prime \prime }`$ from the positions of the compact radio remnants, except in one case. Most of the sources are likely to be X-ray binaries, or supernova remnants that are not strong radio emitters. The luminosity limit of the image is $`10^{37}\mathrm{erg}\mathrm{s}^1`$. While especially the 2 – 10 keV luminosities would be observable with Chandra and in possible conflict with the small number of X-ray sources coinciding with the compact radio sources, we note that the peak of the X-ray flux occurs considerably before the radiative transition when we expect an enhanced radio flux. The radio and X-ray bright remnants may therefore represent two different populations, with the former more numerous because of their greater ages. An uncertainty in the predicted X-ray emission is that it is strongest at a small radius when the effects of winds from the progenitor may be a factor in the surrounding density distribution.

There are possible problems with our picture. VLBI observations by Pedlar et al. (1999) show that the remnant 43.3+59.2 is expanding at $`9850\pm 1500\mathrm{km}\mathrm{s}^1`$, as opposed to the $`500\mathrm{km}\mathrm{s}^1`$ that we expect. This remnant shows a relatively constant flux (Kronberg et al. 2000), so it belongs to the main class of remnant. If further VLBI measurements confirm the initial result, and find similar rapid expansion in other remnants, our model will have to be abandoned. Also, Muxlow et al. (1994) find the number – diameter relation for the remnants to be consistent with free expansion; the number of remnants with size less than diameter $`D`$ increases linearly with $`D`$. However, the relation may be due to decelerating remnants evolving in an interstellar medium with density variations (Berkhuijsen 1987). In our model, the remnants would be brightest when they first become radiative, which is a function of density.

There are arguments, in addition to the expected rapid evolution, against free expansion for the remnants at a velocity $`10^4\mathrm{km}\mathrm{s}^1`$. The high radio luminosities of the remnants suggests substantial thermalization and deceleration of the ejecta. Huang et al. (1994) argued that the surface brightness – diameter relation that applies to the M82, LMC, and Galactic remnants could be approximately understood by assuming that some constant fraction of the supernova energy goes into relativistic particles and magnetic fields. The small diameter remnants are not observable in the Galaxy and LMC because they have not interacted with a sufficiently dense medium to thermalize the ejecta energy.

There is also the rate of formation of the radio remnants. If the remnants expand at $`10^4\mathrm{km}\mathrm{s}^1`$ and there are 28 with diameters $`<3`$ pc (Huang et al. 1994), the formation rate is 1/(5 yr) (see also Muxlow et al. 1994). In our model, the remnants have mean expansion velocities about an order of magnitude smaller and a correspondingly lower formation rate. No new remnants have been found above the flux level observed by Kronberg et al. (1985) in the past 20 years, so the observations are more consistent with the slow expansion model.

An individual radiative remnant is expected to have a peak total luminosity $`10^{40}\mathrm{erg}\mathrm{s}^1`$, most of which is expected to be at infrared wavelengths (Wheeler et al. 1980; Draine & Woods 1991). The problem with observing this emission is that the starburst nucleus of M82 has an infrared luminosity $`10^{44}\mathrm{erg}\mathrm{s}^1`$. For the high absorption toward the compact remnants, the observable emission is likely due to dust continuum emission and infrared fine structure lines. The \[OI\] 63$`\mu `$m line dominates the postshock cooling below 5,000 K and should be prominent (Hollenbach & McKee 1989). For a newly cooling remnant in a medium with $`n_H=10^3`$ cm-3 and shock velocity $`v_{sh}=426\mathrm{km}\mathrm{s}^1`$ and $`R=1.37`$ pc (Draine & Woods 1991), the expected \[OI\] 63 $`\mu `$m luminosity is $`1\times 10^{37}\mathrm{erg}\mathrm{s}^1`$ (Chevalier 1999). The total observed \[OI\] 63 $`\mu `$m luminosity from the starburst region of M82 is $`1\times 10^{41}\mathrm{erg}\mathrm{s}^1`$, most of which is thought to arise from warm neutral gas photodissociated by radiation from OB stars (Lord et al. 1996). The expected line emission is faint, but might be observable because the line is broad and is expected in localized regions. However, we estimate that the line is somewhat below the expected detection limit of the AIRES echelle spectrograph on SOFIA.

Greenhouse et al. (1997) have detected 6 sources of \[Fe II\] 1.2 and 1.6 $`\mu `$m line emission that they identify as shock emission from supernova remnants. The emission does not coincide with that from the radio remnants. We attribute this partially to the high extinction to the radio supernova remnants. Mattila & Meikle (2001) find an average $`A_V24`$ to the radio remnants, while Greenhouse et al. (1997) find $`A_V=69`$ to the \[Fe II\] remnants. Taking $`A_V=7.5`$ to the \[Fe II\] remnants and the extinction curve of Fitzpatrick (1999), the additional extinction would reduce the observed flux of the radio remnants by a factor of $`13`$ compared to the \[Fe II\] remnants. In addition, both theoretical shock models (Hollenbach & McKee 1989) and observations of a Galactic remnant (Oliva et al. 1999) indicate that the luminosity ratio $`L([\mathrm{FeII}]1.65\mu \mathrm{m})/L([\mathrm{OI}]63\mu \mathrm{m})<1`$, leading to a low expected \[Fe II\] 1.65 $`\mu `$m line luminosity. The same considerations indicate that the \[Fe II\] remnants should be strong \[O I\] 63 $`\mu `$m line sources, easily detectable with SOFIA. The uncertainties in this argument are that the expected shock velocities in the radio remnants are higher than have been modeled or observed, and the shock velocities in the \[Fe II\] remnants are unknown.

Three of the \[Fe II\] remnants are resolved with diameters of $`40`$ pc and the remainder are unresolved at a diameter of 22 pc (Greenhouse et al. 1997), an order of magnitude larger than the radio remnants. Expansion in a $`10^3`$ cm-3 medium would lead to substantial energy loss and a lower \[Fe II\] luminosity, yet the observed high luminosity requires interaction with relatively dense gas. We suspect that the remnants are primarily interacting with a low density gas, and the \[Fe II\] emission is from radiative shock waves in clouds. The positions of the remnants extend to a larger distance from the center of M82 than do the radio remnants (Greenhouse et al. 1997). The \[Fe II\] remnants should be weak nonthermal radio sources; based on the surface brightness – diameter relation approximately followed by radio remnants (Huang et al. 1994), we estimate radio fluxes of $`0.010.04`$ mJy. As noted by Greenhouse et al. (1997), the high radio background in M82 makes faint sources difficult to detect.

## 3 OTHER GALAXIES AND SUPERWINDS

Ulvestad & Antonucci (1997) have studied compact radio sources in NGC 253 with a resolution of 1 pc. They find that both HII regions and supernova remnants are present (based on their spectral indices), that some of the sources are resolved, and that there is little evidence for flux variation over a timespan of 8 yr. The upper limits on the rate of flux decrease of the strongest sources are $`12`$% yr-1. These properties indicate that the remnants are quite similar to those in M82.

The sources found by Smith et al. (1998) in the NW nucleus of the Arp 220 merging system have somewhat different properties: they are unresolved with a $`1.1\times 2.9`$ pc beam and they are more luminous than the M82 sources. Their luminosities are comparable to that of the luminous radio supernova SN 1986J and Smith et al. (1998) identify the objects with this class of radio source. This requires that essentially all of the supernovae in the nuclear region of Arp 220 produce very luminous radio supernovae. Events like SN 1986J are thought to involve interaction with an especially dense wind and are rare among supernova explosions. However, the interstellar conditions in the nuclear region of Arp 220 are probably more extreme than those in M82 and may account for the difference in the radio source properties. The projected area covered by the radio sources is more compact: $`75\times 150`$ pc in Arp 220 as compared to $`100\times 600`$ pc in M82. The infrared luminosity in Arp 220 is $`10^{12}L_{}`$ as compared to $`4\times 10^{10}L_{}`$ in M82, which suggests that the supernova rate is correspondingly higher in Arp 220. In the spherically symmetric galactic wind theory of Chevalier & Clegg (1985), the central pressure is $`p_c\dot{E}^{1/2}\dot{M}_w^{1/2}R^2`$, where $`\dot{E}`$ is the total power input from supernova, $`\dot{M}_w`$ is the mass loss rate in the galactic wind, and $`R`$ is the radius of the region of power input. The pressure in the starburst region of Arp 220 is plausibly $`10`$ greater than that in M82 and the density in the interclump region of molecular clouds may have a correspondingly higher density. Solomon (2001) in fact finds that the mean density in the starburst region of Arp 220 is $`10^4`$ cm-3. For $`n_H=10^4`$ cm-3, the cooling radius for a supernova remnant is $`0.5`$ pc when its age is $`100`$ yr (Draine & Woods 1991). We thus predict that the remnants are only somewhat smaller than the current size limits and that they evolve more rapidly than the remnants in M82.

In our model for the compact remnants, they are radiative and thus do not efficiently provide mechanical energy to heat gas in the starburst region. In addition, our estimate of ages of the objects is larger than in models where they expand at a high constant velocity (e.g., Muxlow et al. 1994) and the rate of supernova explosions is correspondingly lower. We hypothesize that there is a population of massive star supernovae that are able to escape from their parent molecular cloud and explode in a low density, hot medium ($`n_H1`$ cm-3). These supernovae are responsible for driving the galactic winds from M82 and other starburst galaxies (Chevalier & Clegg 1985; Heckman, Armus, & Miley 1990). Unfortunately, the remnants of these supernovae are difficult to observe; Fig. 1 shows that remnants occurring in the low density component should have luminosities $`<10^{36}\mathrm{erg}\mathrm{s}^1`$ and are not detectable with Chandra. However, the initial supernovae may be directly observable (e.g., Bregman, Temi, & Rank 2000). Supernovae in the low density medium may have smaller line-sight column densities that the ones in molecular clouds.

RAC is grateful to R. O’Connell and J. Ulvestad for useful conversations, and we thank an anonymous referee for a very helpful report. Support for this work was provided in part by NASA grant NAG5-8232.


format: text

text:
# Dark Matter in the center of galaxies and galaxy clusters: ruling out the CDM scenario?

## 1. Introduction

Cold Dark Matter (CDM) models provide a solid framework capable to explain most of the properties of the universe at large scales. In the last few years a very large effort has been spent in the investigation of the detailed structure of the virialized haloes within the CDM scenario, in which dark particles are assumed to be cold and collisionless. Using N-body simulations it was showed that virialized haloes are well described by an universal density profile, the Navarro, Frenk, $`\&`$ White model (1997; hereafter NFW). This profile diverges at the center showing a cuspy core: $`\rho 1/r`$. Recent high-resolution N-body simulations show that as the numerical resolution is increased, the resulting dark density profile goes as $`\rho 1/r^{1.5}`$ (Moore et al. 1999). On galactic scales the theory seem to be in conflict with the observations, since the rotation curves of LSB and dwarf galaxies rise too slowly with respect to the predictions of the CDM model. Indeed, the rotation curves of these galaxies call for a finite central density (soft core), in conflict with the cuspy cores predicted by the model. Note that LSB and dwarf galaxies are systems strongly dark matter dominated, so their rotation curves are good tracers of the underlying dark haloes gravitational potential and therefore good candidates to explore the innermost shape of the dark matter distribution. Recently, however, some authors have challenged the existence of soft cores in centers of galaxies measuring H$`\alpha `$ rotation curves of late-type dwarf and LSB galaxies. Basically they find two results: HI rotation curves for most of the galaxies are affected by beam-smearing and a good agreement with concentration values predicted by the NFW model (van den Bosch and Swaters 2000).

## 2. Soft core observational evidence

We have carried out H$`\alpha `$ rotation curves at the TNG of dwarf and LSB galaxies. The spatial resolution in the central regions is a factor 10 better than with respect the HI data published in literature for the same galaxies. In Fig.1 we reproduce the circular velocities for the galaxies UGC4325, UGC11861 and F571-8 (empty circles are our H$`\alpha `$ observations and filled circles are HI data). From a comparison between optical and radio data it is interesting to note that HI observations for these galaxies are not affected (or are very little) by beam smearing (see Marchesini et al. at this conference, for details). To discriminate between the NFW and King models as representative profile of the primeval dark halo, we built a fiducial galaxy. The dotted lines in the top panels represent a NFW profile while the dotted lines in the bottom panels are due to a King profile. Accounting for the disc formation (dashed lines) we computed next the adiabatic contraction of the DM component due to the cooling of the baryons into the virialized haloes. The final rotation curve is then obtained (solid lines) and compared with the H$`\alpha `$ data. The final curve is best fitted for all galaxies assuming that the dark matter in the primeval halo is distribuited as a King profile. On galaxy cluster scale, there are at least two clusters with evidence of soft cores: CL0024+1654 as inferred by strong lensing technique (Tyson et al. 1998) and A1705 as revealed by the Chandra data (Ettori et al. at this conference). The central densities of the galaxies analyzed here or taken from the literature and of CL0024+1654 seem to confirm the evidence that the central density is independent on the halo mass (Firmani et al. 2000). Future observations will increase the sample and put these observational findings on a more robust basis.

## 3. Implications on the DM nature

The question is: how can we justify the existence of soft cores in galaxies and in galaxy clusters? Warm dark matter (WDM) has been proposed in order to solve the conflict with CDM models. However, even if interesting, this scenario is unable to produce soft cores at galactic scales (see e.g. Avila-Reese et al. 2001). A possible solution is to assume a weakly self-interacting cold dark matter (SIDM) (Spergel $`\&`$ Steinhardt 2000). It is well known that the NFW profile is far from thermal equilibrium, and shows, in fact, a inner positive temperture gradient. Collisions between dark particles would modify the velocity distribution to a Maxwellian distribution generating a constant velocity dispersion. We have developed different approaches to explore the SIDM regime: thermodynamic (Firmani et al. 2000) and dynamic models building a cosmological Boltzmann code (Firmani, D’Onghia, $`\&`$ Chincarini 2000, astro-ph/0010497). Here, we propose a new numerical procedure applied to N-body simulations (D’Onghia, Firmani, $`\&`$ Chincarini 2001, in preparation). Our technique is capable to simulate a single halo of any size and mass in a cosmological framework coupling semi-numerical technique and N-body simulations. The cosmological initial conditions are fixed by the mass aggregation history of the halo using semi-analytical models and only the dynamical evolution of the halo is followed by the N-body code. To achieve this, the HYDRA code has been ad hoc modified and our version can describe simulations both of collisional and collisionless dark haloes. In our simulations we assume a cross section inversely proportional to the particle velocity dispersion: $`\sigma /m_xv_{100}=10^{24}`$ cm2/Gev with $`v_{100}`$ the halo velocity dispersion in units of 100 km/s. The top panels of Fig. 2 show the dark density profiles of haloes of $`M=10^{11}M_{}`$ (left panel) and $`10^{15}M_{}`$ (right panel) as obtained in a SIDM model. Each halo have been modeled assuming the cross section value above (in the figure the corresponding density profiles have empty circles) and a cross section three times larger (filled circles). In the same plot solid lines are the NFW profiles. Within a Hubble time a modest self-interaction cross section value creates soft cores in haloes of any size and a three times stronger cross section value produces larger soft cores. Bottom panels of Fig. 2 show the corresponding halo radial dispersion velocity profiles (symbols have the same meanings as above). This is the case in which any trend towards the core collapse is avoided by the competition between a mass accretion determined by the halo merging history and a thermalization process by collisions.

## 4. Conclusions

H$`\alpha `$ rotation curves carried out at TNG of two late-type dwarf galaxies and one LSB galaxy show evidence of soft cores as the mass distribution of CL0024+1654 on galaxy cluster scales. Integrating the information from galaxies to clusters of galaxies the halo central density seems to be independent of the halo mass. WDM in unable to produce visible soft cores, but SIDM in a very weak self-interacting regime ($`\sigma /m_xv_{100}=10^{24}`$ cm2/Gev) produces soft cores in agreement with the available observations on different scales.

## 5. Acknowledgments

E.D. and D.M. are grateful to Richard Bower and Oleg Gnedin for stimulating discussions about this topic during the congress.

## 6. References

Avila-Reese V.,Colin P.,Valenzuela O.,D’Onghia E.,Firmani C., 2001,(astro-ph/0010525)
Firmani C.,D’Onghia E.,Avila-Reese V.,Chincarini G., Hernández X.,2000,MNRAS,315,L29
Moore B., Quinn T., Governato F., Stadel J., Lake G., 1999, MNRAS, 310, 1147
Navarro J., Frenk C.S., White S.D.M., 1997, ApJ, 490, 493
Spergel D., Steinhardt P., 2000, PhRvL, 84, 3760
Tyson J.A., Kochanski G.P., Dell’Antonio I.P., 1998, ApJ, 498, L107
van den Bosch F.C., Swaters R.A., 2001, MNRAS, 325, 1017


text:
# Relativistic static fluid spheres with a linear equation of state

## I Introduction

Einstein’s equations for static and spherically symmetric perfect fluids have been studied by many authors \[1-3\]. The unknown functions are more than the equations, therefore, either an equation of state should be prescribed or some simplifying ansatz assumed. Schwarzschild found the first solution  in 1916. It was a degenerate case ($`n=0`$) of the linear equation of state

$$\rho =np+\rho _0$$
(1)
where $`\rho `$ and $`p`$ are the density and the pressure of the fluid and $`n`$, $`\rho _0`$ are constant parameters. A year later, two other solutions satisfying Eq. (1) were found, describing static cosmological models . Tolman  used eight simplifying ansatze for the metric components to find new solutions. Some limiting cases of these satisfy Eq. (1) too, leading to a simple, but singular solution for arbitrary $`n`$ and $`\rho _0=0`$ . Klein rediscovered this solution first for the case of incoherent radiation $`n=3`$  and then for general $`n`$  by studying systematically the equation of state (1) with vanishing $`\rho _0`$ (the so-called $`\gamma `$-law). This Klein-Tolman (KT) solution was rediscovered later three times more \[10-12\]. Meanwhile, two new solutions with $`\rho _00`$ appeared in 1968, developing further the idea of Schwarzschild’s incompressible sphere. Whittaker argued that the gravitational mass density $`\rho +3p`$ should be constant , while Buchdahl and Land demanded that the speed of sound should be the largest possible, namely, equal to the speed of light .

The main stream of papers, however, utilized simplifying assumptions, which often lead to complicated equations of state, if any. Already Tolman used the ansatz

$$e^\lambda =a+br^2$$
(2)
where $`a`$ and $`b`$ are constants, $`0a1`$, and $`e^\lambda `$ is the $`g_{rr}`$ component of the metric. A thorough study of this ansatz was performed by Kuchowicz  who found seven different cases in total, many of them having rather involved expressions for the pressure and the density. He didn’t study whether his solutions satisfy Eq. (1).

What happens when we impose on the field equations both Eqs. (1) and (2)? An overdetermined system of differential equations is obtained which may be devoid of any non-trivial solutions. The real answer, derived in this paper, is somewhat surprising. All known explicit solutions with a linear equation of state satisfy also Eq. (2), except for the Whittaker solution. We find an additional solution with $`n=3`$ instead, which appears to be new.

It is much easier to impose both requirements directly on the field equations, than to study case by case the results of Ref. . In this way we derive in a unified and simple manner almost all known solutions with linear equation of state, including a new hybrid between the Einstein static universe (ESU) and the KT solution. This programme is executed in the following section. The last section contains some discussion.

## II Field equations and their solutions

The metric of a static spherically symmetric spacetime reads

$$ds^2=e^\nu dt^2e^\lambda dr^2r^2d\mathrm{\Omega }^2$$
(3)
where $`d\mathrm{\Omega }^2`$ is the metric on the two-sphere and $`\nu `$, $`\lambda `$ depend on $`r`$. The system of equations is

$$\rho =\frac{\lambda ^{}}{r}e^\lambda +\frac{1}{r^2}\left(1e^\lambda \right)$$
(4)
$$p=\frac{\nu ^{}}{r}e^\lambda \frac{1}{r^2}\left(1e^\lambda \right)$$
(5)
$$p^{}+\frac{1}{2}\left(\rho +p\right)\nu ^{}=0$$
(6)
where the prime means a $`r`$-derivative and one of the field equations has been replaced by the contracted Bianchi identity (6). Let us impose Eq. (2). Then Eq. (4) gives a simple expression for the density

$$\rho =3b+\frac{1a}{r^2}$$
(7)
which is regular at the centre if $`a=1`$. When Eq. (1) holds, Eq. (6) may be integrated. Two different cases arise, $`n=1`$ and $`n1`$.

Let us discuss first the case $`n=1`$. Then Eq. (6) yields

$$p=\frac{1}{2}\rho _0\nu +h$$
(8)
with $`h`$ being an integration constant. Plugging Eq. (8) into Eq. (1) and the resulting expression for $`\rho `$ into Eq. (7) gives

$$\frac{1}{2}\rho _0\nu =h\rho _03b+\frac{1a}{r^2}$$
(9)
There are two possibilities. If $`\rho _0=0`$ then from Eq. (8) $`p=h`$, from Eq. (1) $`\rho =p`$ and from Eq. (9) $`a=1`$, $`b=h/3`$. Therefore

$$e^\lambda =1+\frac{h}{3}r^2$$
(10)
while Eq. (5) gives $`\nu `$ after a simple integration: $`\nu =\lambda `$. This is the de Sitter solution.

If $`\rho _00`$ then the insertion of Eqs. (2,8,9) into Eq. (5) yields the relation

$$\rho _0\left(\rho _0+2b\right)r^4=4\left(1a\right)\left(a+br^2\right)$$
(11)
Obviously, $`a=1`$ and $`b=\rho _0/2`$. From Eq. (7) $`\rho =3\rho _0/2`$ and from Eq. (1) $`p=\rho _0/2`$. Eq. (8) gives $`\nu =1+2h/\rho _0=const`$, while Eq. (2) becomes

$$e^\lambda =1\frac{\rho _0}{2}r^2$$
(12)
The pressure and the density satisfy, in fact, the equation of state $`\rho +3p=0`$. Hence, we obtain ESU, although initially $`\rho =p+\rho _0`$.

Let us continue with the generic case $`n1`$. Now the integration of Eq. (6) brings a different answer

$$p=p_1e^{\frac{n+1}{2}\nu }p_0$$
(13)
where $`p_1`$ is a constant of integration, while $`p_0=\rho _0/\left(n+1\right)`$. Proceeding like we have done with Eq. (8) the analog of Eq. (9) is obtained

$$np_1e^{\frac{n+1}{2}\nu }np_0+\rho _0=3b+\frac{1a}{r^2}$$
(14)
Again, there are two possibilities. If $`n=0`$ we have $`a=1`$, $`b=\rho _0/3`$ and from Eq. (1) $`\rho =\rho _0=const`$. The metric component $`\nu `$ is obtained not from Eq. (14) but from Eq. (5), the last equation to be satisfied. For a general $`n`$ it becomes a linear equation for $`y=e^{\frac{n+1}{2}\nu }:`$

$$\frac{2}{n+1}r\left(a+br^2\right)y^{}=\left[1a\left(p_0+b\right)r^2\right]y+p_1r^2$$
(15)
In our case it reduces to

$$2\left(1\frac{\rho _0}{3}r^2\right)y^{}=\frac{2}{3}r\rho _0y+p_1r$$
(16)
and is easily integrated

$$e^\nu =\left(\frac{3p_1}{2\rho _0}+Ce^{\lambda /2}\right)^2$$
(17)
$$e^\lambda =1\frac{\rho _0}{3}r^2$$
(18)
which, after some change of notation, coincides with solution III of Tolman. This is the Schwarzschild interior solution.

If $`n0`$ both Eqs. (14) and (15) determine $`\nu `$ and we must study their compatibility. Eq. (14) may be written as

$$\frac{1}{y}=A+\frac{B}{r^2}$$
(19)
$$A=\frac{1}{n\left(n+1\right)p_1}\left[\rho _0+3\left(n+1\right)b\right]$$
(20)
$$B=\frac{1a}{np_1}$$
(21)
When these relations are put into Eq. (15), a polynomial of second degree in $`r^2`$ arises, which must vanish identically. Three compatibility conditions result

$$\left(1a\right)\left[\frac{4n+\left(n+1\right)^2}{\left(n+1\right)^2}a1\right]=0$$
(22)
$$A\left[\rho _0+\left(n+3\right)b\right]=0$$
(23)
$$B\left[\left(n+1\right)\rho _0+4nb+\left(n+1\right)\left(n+3\right)b\right]\left(n+1\right)^2\left(1a\right)A=0$$
(24)
Eq. (22) has two solutions for $`a`$. Let $`a=1`$. Then $`B=0`$, $`\nu `$ is a constant and $`\rho =3b`$. Eq. (24) is also satisfied, while Eq. (23) provides two subcases. If $`A=0`$ then from Eq. (19) $`y=\mathrm{}`$ which is unacceptable. If $`\rho _0=\left(n+3\right)b`$ then Eq. (20) shows that $`A0`$ when $`b0`$. Eq. (1) gives $`p=b`$, hence the equation of state is $`\rho +3p=0`$ for any $`n`$. We again obtain the ESU with constant $`\nu `$ and

$$e^\lambda =1+br^2$$
(25)
This is the same as Eq. (12), since both can be written as $`e^\lambda =1+pr^2`$.

The last case which remains is

$$a=\frac{\left(n+1\right)^2}{4n+\left(n+1\right)^2}$$
(26)
where $`a1`$ since $`n0`$. Eq. (23), like before, offers two choices.

a) $`A=0`$ ($`\rho _0=3\left(n+1\right)b`$). Then Eq. (24) simplifies to $`\left(n1\right)b=0`$. There are two subcases:

a1) $`n=1`$, $`b0`$. Then $`\rho _0=6b`$, Eq. (26) gives $`a=1/2`$, Eq. (21) gives $`B=1/2p_1`$ and Eqs. (2,19,7,1) give respectively

$$e^\lambda =\frac{1}{2}\left(1\frac{\rho _0}{3}r^2\right)$$
(27)
$$e^\nu =2p_1r^2$$
(28)
$$\rho =\frac{1}{2}\left(\rho _0+\frac{1}{r^2}\right)$$
(29)
$$\rho =p+\rho _0$$
(30)
This is precisely the solution of Buchdahl and Land  after $`e^\nu `$ is rescaled.

a2) $`b=0`$. Then $`\rho _0=0`$ and we easily find

$$e^\lambda =\frac{\left(n+1\right)^2}{4n+\left(n+1\right)^2}$$
(31)
$$e^{\frac{n+1}{2}\nu }=\frac{4}{\left[4n+\left(n+1\right)^2\right]p_1r^2}$$
(32)
$$\rho =np=\frac{4n}{\left[4n+\left(n+1\right)^2\right]r^2}$$
(33)
This is the KT solution in curvature coordinates \[3,8-12\]. Eq. (31) shows that the solution exists when $`4n+\left(n+1\right)^2>0`$. This condition is fulfilled for $`n`$ outside the interval $`\left(5.83=32\sqrt{2},3+2\sqrt{2}=0.17\right)`$. When $`n<5.83`$ or $`0.17<n<0`$ the solution has $`p`$ and $`\rho `$ of different signs. The physically realistic range is $`1n\mathrm{}`$, which includes the important cases of stiff fluid ($`n=1`$) and incoherent radiation ($`n=3`$). The KT solution is singular at $`r=0`$ but it can serve as an interior solution for $`r>r_0>0`$ in multi-layered perfect fluid models. It also represents a focal point for a regular solution, which has been studied numerically .

b) $`\rho _0=\left(n+3\right)b`$. Since $`a1`$, Eq. (24) reduces to $`\left(n+3\right)b=0`$. The two subcases are

b1) $`b=0`$. Then $`A=0`$ and we go back to case a).

b2) $`n=3`$. Then $`\rho _0=0`$, $`a=1/2`$, $`\rho +3p=0`$ and

$$e^\lambda =\frac{1}{2}+br^2$$
(34)
$$p_1e^\nu =b\frac{1}{2r^2}$$
(35)
$$\rho =3b+\frac{3}{2r^2}$$
(36)
The constant $`p_1`$ may be set to one by a time rescaling. This appears to be a new solution - a hybrid between ESU with its constant pressure and density and the KT solution. Its pressure, density and metric component $`g_{tt}`$ are linear combinations of those two solutions, although Eqs. (4-6) form a highly non-linear system. When $`b=0`$ we obtain the KT solution for $`n=3`$ which, in fact, does not exist. When $`b0`$, $`e^\lambda `$ and $`e^\nu `$ are negative for small $`r`$. However, if we choose $`r_0^2=1/2b`$ and $`b>0`$, then the solution exists for $`r>r_0`$, has positive pressure, negative density and represents a regular companion of ESU.

Whittaker has solved the general case $`n=3`$ when $`\rho _0`$ may be different from zero and the ansatz (2) is not imposed . However, the above solution has been missed. The reason is the following. Eqs. (1,4,5) lead to the relation

$$\left(e^{\lambda +\nu }\right)^{}=kre^{2\left(\lambda +\nu \right)}$$
(37)
where $`k`$ is a constant. The integration of this equation yields

$$e^{\lambda +\nu }=\frac{2}{ckr^2}$$
(38)
$`c`$ being a constant. Whittaker assumes that $`k=c\alpha `$, where $`\alpha `$ is yet another constant. This choice excludes the above solution which has $`c=0`$, $`k=2p_10`$ and in principle may be regained as a limiting case. However, when one sets $`\rho _0=0`$ in the results of Ref.  only ESU is obtained. It should be noticed that the singular term in Eq. (35) is different from the universal Schwarcschild term $`m/r`$ which was correctly discarded from the interior solution.

## III Discussion

In this paper we have studied the gravitation of static spherically symmetric perfect fluid solutions with a linear equation of state. We imposed the additional requirement (2), making the system of equations overdetermined. Amazingly, it gives six non-trivial solutions: the Schwarzschild interior solution for an incompressible sphere, the de Sitter cosmological solution, the Einstein static universe, the stiff fluid solution of Buchdahl and Land, the Klein-Tolman singular solution and a new solution with $`\rho +3p=0`$. Among the known solutions with linear equation of state only the Whittaker solution does not satisfy Eq. (2). The reason is that it satisfies Eq. (38) and if Eq (2) is required too, $`\nu `$ has the only option to remain a constant, which is too restrictive. These solutions have been found originally starting from a variety of physical considerations and following quite different mathematical approaches. Here they are derived in a unified way as a result of straightforward calculations and a new solution is added. Roughly speaking, it is a linear combination between ESU and the KT solution, which is unusual for a non-linear system of differential equations. It closes the list of solutions satisfying both Eqs. (1) and (2). The equation of state $`\rho +3p=const`$ is unrealistic because the speed of sound is negative. However, models which satisfy it have interesting properties, are often integrable and occupy important places in different classification schemes. This concerns one whole class of algebraically special perfect fluids of Petrov type II or D  (Eq.(29.20)), a rigidly rotating solution of type D , the well-known Wahlquist solution \[21-22\] which has the Whittaker solution as a static limit and the Kerr metric as a vacuum subcase, and the recently found Wahlquist-Newman solution  which reduces to a number of well-known solutions in appropriate limits. An interesting phenomenon is the effective change of the equation of state by which ESU is obtained first from $`\rho =p+\rho _0`$ and then from $`\rho =np+\rho _0`$. When $`\rho _0=0`$ a one-parameter sequence of KT solutions exists which serve as focal points for regular numerical solutions. When $`\rho _00`$ only two discrete cases satisfy Eq. (2) for $`n=0`$ and $`n=1`$, the third known one ($`n=3`$) drops out of this scheme. Presumably when one replaces the ansatz (2) by a different one, other integrable cases will appear, although the system is overdetermined.


text:
# Improved tensor–product expansions for the two–particle density matrix

## Abstract

We present a new density–matrix functional within the recently introduced framework for tensor–product expansions of the two–particle density matrix. It performs well both for the homogeneous electron gas as well as atoms. For the homogeneous electron gas, it performs significantly better than all previous density-matrix functionals, becoming very accurate for high densities and outperforming Hartree-Fock at metallic valence electron densities. For isolated atoms and ions, it is on a par with generalized gradient approximations to density-functional theory. We also present analytic results for the correlation energy in the low density limit of the free electron gas for a broad class of such functionals. PACS 71.10.-w 71.15.Mb, submitted to Physical Review Letters

The long standing quest for an accurate description of electronic correlation with effective single–particle theories has renewed interest in density–matrix functionals in recent years. The success of density functional approaches, based on Hohenberg–Kohn–Sham theory , is well established . Unfortunately, the known approximations to the requisite universal energy functional of the density are not sufficiently accurate to predict the rates of chemical and physical processes at room temperatures. This theory can easily be generalized to show that there is a universal correlation energy functional of the single–particle density matrix as well . A density matrix functional requires an expression of the two–body matrix in terms of the one-body density matrix. An ansatz for the two–body matrix containing one free parameter was proposed by Müller , and the case where this parameter is set to 1/2 has been discussed by Buijse . In addition, Buijse attempted to perform self-consistent calculations but optimized only the lowest two natural orbitals and occupation numbers and only for the H2 molecule. Independently, Goedecker and Umrigar  proposed a functional based on a very similar, but self-interaction corrected, two–body matrix and derived the Euler–Lagrange formalism necessary for the numerical minimization of the functional. Optimizing typically 50 natural orbitals and occupation numbers, necessary to achieve convergence, for a variety of atoms and ions it was found that the Goedecker-Umrigar (GU) functional yielded energies that were comparable or better than those from the generalized gradient approximation (GGA) in density functional theory and densities that were better than those from GGA. Following these first encouraging numerical results several groups of authors have recently proposed and studied the accuracy of various new approximations to the density-matrix functional. However, numerical studies  demonstrate that density–matrix functionals, to date, give poor correlation energies for the homogeneous electron gas.

It has been shown analytically that the per particle correlation energy of the homogeneous electron gas in the low density limit approaches a non-zero constant for these initial functionals. Building on our tensor–product formalism, below we extend this proof and demonstrate this shortcoming for a broad class of functionals. As regions of low densities contribute relatively little to the total energies of solids and molecules, we focus our attention on improving the results for the high and intermediate density regimes. We present a new density–matrix functional which is more than an order of magnitude more accurate in the high–density limit for the homogeneous electron gas than those proposed previously and also represents a significant improvement in the intermediate regime. In atoms, the new functional performs similarly to the generalized gradient approximation (GGA) to density functional theory but not as well as the GU functional. Because selfinteraction corrections complicate numerical calculations, we do not enforce them on our new functional in the present work.

Notation — For a system of $`N`$ electrons, the exact energy functional of the two–body density matrix $`\gamma (𝐱_1𝐱_2,𝐱_1^{}𝐱_2^{})\mathrm{\Psi }|\widehat{\psi }^{}(𝐱_1^{})\widehat{\psi }^{}(𝐱_2^{})\widehat{\psi }(𝐱_2)\widehat{\psi }(𝐱_1)|\mathrm{\Psi }/2`$, is

$`E`$ $`=`$ $`{\displaystyle 𝑑𝐱_1\left(\left[\frac{1}{2}_{𝐫_1}^2+U(𝐫_1)\right]n(𝐱_1,𝐱_1^{})\right)}|_{𝐱_1^{}=𝐱_1}`$ (2)
$`+{\displaystyle 𝑑𝐱_1𝑑𝐱_2\frac{\gamma (𝐱_1𝐱_2,𝐱_1𝐱_2)}{|𝐫_1𝐫_2|}}.`$
Here and throughout, we work in atomic units, and the $`𝐱_i`$ are compound coordinates representing both position $`𝐫_i`$ and spin $`s_i`$, so that integration over $`𝐱_i`$ represents integration over space and summation over spin channels. In the above expression, the external potential and one–body density matrix are $`U(𝐫)`$ and $`n(𝐱_1,𝐱_1^{})\mathrm{\Psi }|\widehat{\psi }^{}(𝐱_1^{})\widehat{\psi }(𝐱_1)|\mathrm{\Psi }`$, respectively. Finally, we note that simple integration of the definitions reveals the sum–rule which connects the one– and two–body matrices,

$$\left(\frac{N1}{2}\right)n(𝐱_1,𝐱_1^{})=𝑑𝐱_2\gamma (𝐱_1𝐱_2,𝐱_1^{}𝐱_2).$$
(3)
As in , we expand $`\gamma `$ into Hartree and exchange–correlation contributions,

$$\widehat{\gamma }=\widehat{\gamma }_\text{H}+\widehat{\gamma }_{\text{XC}}.$$
(4)
By separation of variables, it is always possible to expand the unknown exchange–correlation part as a (potentially infinite) sum of tensor products,

$$\widehat{\gamma }_{\text{XC}}=\underset{i}{}\widehat{u}_i\widehat{v}_i.$$
(5)
In this context we view the two–body density matrix as a two–body operator (a function of four variables) and $`\widehat{u}_i`$ and $`\widehat{v}_i`$ as one–body operators (functions of two variables). Note that throughout we employ the twisted tensor product (denoted as “Type III” product in ),

$`[uv](𝐱_1𝐱_1^{},𝐱_2𝐱_2^{})u(𝐱_1,𝐱_2^{})v(𝐱_1^{},𝐱_2),`$

to combine the one–body operators.

Remarkably, given such an expansion truncated to finite order, the combination of the symmetry constraints on $`\gamma `$ with the sum–rule is sometimes sufficient to determine the full four–variable function $`\gamma `$ directly in terms of the two–variable density matrix $`n`$. Under these conditions, one can then evaluate the energy as an explicit functional of the density matrix using (2). Minimizing the density–matrix functional over all physical density matrices then yields the ground–state of the system under this approximation.

New Functional — Considerable freedom remains in the terms for the expansion (5). Previously, we explored two choices, “Corrected Hartree–Fock” (CHF), and the “Corrected Hartree” (CH) functional,

$`\widehat{\gamma }_{\text{XC}}^{\text{CHF}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\widehat{n}\widehat{n}){\displaystyle \frac{1}{2}}(\sqrt{\widehat{n}(1\widehat{n})}\sqrt{\widehat{n}(1\widehat{n})}),`$
$`\widehat{\gamma }_{\text{XC}}^{\text{CH}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\sqrt{\widehat{n}}\sqrt{\widehat{n}}),`$
respectively. The latter is identical to the one proposed by Buijse  and a special case of the one proposed by Müller . The GU functional is the CH functional with a selfinteraction correction. In the case of the homogeneous electron gas at high densities the GU and CH functionals become identical since the selfinteraction correction vanishes but at low densities they could differ if the natural orbitals are localized. One way to view these functionals is in terms of the coefficient of the Fock term $`\widehat{n}\widehat{n}`$, which is $`1/2`$ in the CHF functional and is zero in the CH functional. Under truncation of the expansion to one additional term, the value of this coefficient, through the sum–rule, then uniquely determines the form of the remaining tensor product.

Homogeneous electron gas — Figure 1 shows the correlation energies of the resulting functionals for the spin–unpolarized homogeneous electron gas. Corrected Hartree–Fock leads to no improvement in the correlation energy over Hartree–Fock, and approaches zero for high electron density. The CH functional is a little better, and represents an improvement by a factor of two over Hartree–Fock at high densities, where it is over-correlated. Both functionals are over-correlated at low densities, where they perform badly, with the per particle correlation energy approaching a non–zero constant. This property of the CH functional has been demonstrated analytically in .

Having seen the behaviour of the above two functionals in the high density regime, it is natural to look for a form with a Fock coefficient intermediate to those above. A coefficient of $`1/4`$ for the Fock term results in the following functional,

$`\widehat{\gamma }_{\text{XC}}={\displaystyle \frac{1}{4}}\left(\widehat{n}\widehat{n}+\sqrt{\widehat{n}(2\widehat{n})}\sqrt{\widehat{n}(2\widehat{n})}\right).`$ (6)
Figure 1 shows that this new functional results in a dramatic improvement at high densities (low $`r_s`$) and significant improvement at metallic valence densities ($`1<r_s<6`$). At low density, the behaviour of all three functionals becomes the same, resulting in unrealistically large correlation energies per particle. Note that because this happens at low density, the contribution of these regions to the total energy is likely to be small in an inhomogeneous system. We now explore the reasons for this behaviour to ascertain the feasibility of making improvements in this regime as well.

Low density limit — To understand the approach of the three functionals shown in Fig. 1 to one another in the low density limit, we recall that Cioslowski and Pernal have determined the predictions of the CH functional analytically for $`r_s`$   $`>`$  $``$    $`5.769`$. Their analysis applies, provided that $`\widehat{\gamma }_{\text{XC}}`$ is represented as a single tensor product of a specific homogeneous form. As an illustration of the anomalous correlation behaviour which these functionals exhibit in the low density limit, Figure 2 compares the known many–body momentum distribution with numerical results for our new functional and the analytic distribution from  derived for the CH functional for $`r_s=10`$. Both density–matrix functionals have very low occupation numbers, and lose the characteristic signature of the sharp drop near the Fermi wave vector evident in the exact curve.

We now show that, as a consequence of the sum–rule, all of the above functionals reduce to this particular form in the limit of high $`r_s`$. Fourier transforming (2), the total energy in the homogeneous electron gas becomes

$$E=2\frac{Vd𝐤}{(2\pi )^3}\frac{𝐤^2}{2}f_𝐤\frac{1}{V}\frac{V^2d𝐤d𝐤^{}}{(2\pi )^6}g_{\mathrm{𝐤𝐤}^{}}\frac{4\pi }{|𝐤^{}𝐤|^2},$$
(7)
where $`g_{\mathrm{𝐤𝐤}^{}}`$ is the momentum space representation of $`\widehat{\gamma }_{\text{XC}}`$. In the case that $`\widehat{\gamma }_{\text{XC}}`$ is expanded as a sum of tensor products of single particle operators, its momentum space representation will be a sum of tensor products of the momentum space representation of those operators. Because our one–body operators are explicit functions of the density matrix, which for the homogeneous electron gas is diagonal in the momentum representation, their momentum space representations are also diagonal. Thus, their tensor products become simple products, and we have

$$g_{\mathrm{𝐤𝐤}^{}}=\mathrm{\Omega }_1[n(𝐤),n(𝐤^{})]+\mathrm{\Omega }_2[n(𝐤),n(𝐤^{})]+\mathrm{},$$
(8)
where the $`\mathrm{\Omega }`$ are the same functions of the density matrix which appear in the tensor products expressing $`\widehat{\gamma }_{\text{XC}}`$ but with the tensor–product operator replaced with simple multiplication. For example, for the CH functional $`g_{\mathrm{𝐤𝐤}^{}}=\frac{1}{2}\sqrt{n(𝐤)n(𝐤^{})}`$.

We now show that, in the low–density limit, the analytic solution presented in  for the case when $`g_{\mathrm{𝐤𝐤}^{}}`$ consists of a single $`\mathrm{\Omega }`$ function of the homogeneous form

$$\mathrm{\Omega }[\lambda n(𝐤),\lambda n(𝐤)]=\lambda ^\beta \mathrm{\Omega }[n(𝐤),n(𝐤^{})],$$
(9)
applies more generally to cases with multiple tensor products. The analysis of  begins with the scaling ansatz

$`n(𝐤)=\rho ^\mathrm{\Delta }\eta (\rho ^\nu 𝐤),`$

where $`\rho `$ is the total electron density per spin channel, $`\nu `$ is a constant, $`\eta `$ is a universal function that does not depend on the density explicitly, and $`\mathrm{\Delta }=13\nu `$ due to the constraint of the total particle density. Note that as $`\rho 0`$, $`n(𝐤)0`$ uniformly (provided that $`\mathrm{\Delta }>0`$). This behaviour, evident in Figure 2, forms the crux of our argument; we verify it self–consistently below.

In the above limit, all of the $`\mathrm{\Omega }_i`$ in the above functionals approach homogeneous functions with their own exponents $`\beta _i`$. Substituting $`n(𝐤)`$ into (8), we then find that a single homogeneous function dominates $`g_{\mathrm{𝐤𝐤}^{}}`$, the one with the smallest $`\beta _i`$ ($`\beta _{\mathrm{min}}`$). Therefore, all the conditions for the analysis of  apply, in particular that

$`\nu `$ $`=`$ $`{\displaystyle \frac{1\beta _{\mathrm{min}}}{3\beta _{\mathrm{min}}2}}.`$
One can show that not only the above functionals, but any functional satisfying the sum–rule and dominated by a single homogeneous term in the low–density limit, will exhibit $`\beta _{\mathrm{min}}=1`$ in that limit. If indeed a single tensor product term dominates $`\widehat{\gamma }`$, then the symmetry arguments presented in  imply that $`\widehat{\gamma }_{\text{XC}}`$ must take the form

$$\widehat{\gamma }_{\text{XC}}=\widehat{u}[\widehat{n}]\widehat{u}[\widehat{n}].$$
(10)
The sum–rule (3) then implies

$$𝑑𝐱_2^{}\widehat{\gamma }_{\text{XC}}=\widehat{n},$$
(11)
which, using (10), can be written as an operator equation

$$(u[\widehat{n}])^2=\widehat{n},$$
(12)
which implies that $`u[]`$ must be a homogeneous function with exponent $`1/2`$ and therefore that $`\widehat{\gamma }_{\text{XC}}`$ is dominated by a term of homogeneous exponent $`\beta _{\mathrm{min}}=1`$. Therefore $`\nu =0`$ and $`\mathrm{\Delta }=1`$; the fact that $`\mathrm{\Delta }>0`$ confirms the self–consistency and the validity of the above analysis.

Ref.  analyzes the case $`\beta =1`$ in detail, finding for the correlation energy per particle the result

$`ϵ_c`$ $`{\displaystyle \frac{1}{8}},`$
as $`\rho 0`$. Figure 1 verifies that, in the low density limit, indeed the correlation energy of all three functionals approach one another and, in particular, that they all approach a non–zero constant correlation energy of $`1/8`$ Hartree per particle.

Atoms — The electrons in solids exhibit both free–electron like and localized atomic–orbital behaviour. Thus, another important limiting case for assessing new functionals is their behavior in atoms. Table I presents energies obtained from the new functional for light atoms and ions. For comparison, we include energies obtained from the CH and GU functionals as well as from the local density approximation (LDA), and the Perdew–Burke–Ernzerhof (PBE)  generalized gradient approximation.

The table reveals that, on the whole, the new functional performs for atoms quite comparably to both the CH functional and to PBE. In particular, for the two–electron series, the new functional is somewhat worse, but for the four–electron series, it is somewhat better. All density–matrix functionals are generally much better than LDA. For the neutral atoms, CH and the present functional are not quite as good as PBE, but, for the ions both are better or in some cases similar to PBE. Note that the CH functional always overestimates the correlation energy, which is not the case for the new functional. In the GU functional this overcorrelation is corrected by the selfinteraction correction and this functional yields the best energies for atoms and ions. Thus, we feel that the density–matrix functionals provide an adequate description of atom–like systems.

Conclusions — We have introduced a new density–matrix functional, which represents a major improvement over existing functionals of the density matrix. It is very accurate in the high density regime of the homogeneous electron gas, significantly improves the correlation energies at typical valence densities, and it is comparable to the generalized gradient approximation in atoms. Although it is not as accurate as the GU functional for atoms, we expect it to perform better than the GU functional for solids. We have also shown that, in our tensor–product expansion of the two–body density matrix, little further improvement can be expected at the present truncation of two terms. Our test systems span the range of environments encountered in solids and molecules, so we conclude that this new functional is a good candidate to be used in electronic structure calculations of condensed matter.

We thank Cyrus Umrigar for interesting discussions and helpful comments on the manuscript.


text:
# Long-Term RXTE Monitoring of the Anomalous X-ray Pulsar 1E 1048.1-5937

## Introduction

The nature of anomalous X-ray pulsars (AXPs) has been a mystery since the discovery of the first example some 20 years ago. Although it is clear that AXPs are young neutron stars, it is not clear why they are observable. In particular, they show no evidence for possessing a binary companion, making conventional accretion problematic. Furthermore, given their spin periods and period derivatives, their rate of loss of rotational kinetic energy is orders of magnitude too small for these sources to be rotation-powered. One important clue is that two AXPs (and one AXP candidate) are clearly associated with supernova remnants. Although only five AXPs are known, their origin is likely to be of great importance to our understanding of the fate of massive stars and the basic properties of the young neutron star population. For an excellent recent review of these objects, see mer99 .

Currently there are two models to explain AXPs. One model proposes that AXPs are young, isolated, highly magnetized neutron stars or “magnetars” dt92a ; td96a . High magnetic fields ($`10^{14}10^{15}`$G) are inferred from their spin-down under the assumption of magnetic dipole braking, as well as by association with the soft gamma repeaters (SGRs) which show AXP-like pulsations in quiescence kds+98 ; ksh+99 , and are thought to have high magnetic fields for independent reasons td96a . The second model of AXP emission is that they are powered by accretion from a fall-back disk of material remaining from the supernova explosion chn00 .

One way to distinguish between these classes of models may be through timing observations. In the magnetar model, relatively smooth spin-down should be expected, punctuated by occasional abrupt spin-up or spin-down events or “glitches,” as well as low-level, long-time-scale deviations from simple spin-down, or “timing noise.” Both phenomena are well known among young radio pulsars (e.g. lyn96 ), although their physical origins in magnetars may be different given the much larger inferred magnetic field. However, according to the magnetar model, no extended spin-up should be seen. On the other hand, accretion power is usually associated with much noisier timing behavior, which can be correlated with spectral, luminosity, and pulse morphology changes. In addition, some accreting binary systems undergo extended ($``$years) episodes of spin-up, although these generally seem to alternate with long intervals of spin-down bcc+97 .

1E 1048.1$``$5937 is a 6.4 s AXP in the Carina region scs86 . It exhibits no evidence for any binary companion, as no Doppler shifts of the pulse period are seen mis98 , and no optical counterpart to a limiting magnitude of $`m_V20`$ has been detected mcb92 . The pulsar’s spectrum, like those of other AXPs, is well described with a two component model consisting of a soft black body with a power-law tail opmi98 . Occasional monitoring observations over more than 20 years show that the pulsar is spinning down, though significant deviations from a simple spin-down model have been noted opmi98 ; pkdn00 ; bss+00 . The paucity of data thus far has made it impossible to unambiguously identify the origin of the deviations.

Here we report on our monthly Rossi X-Ray Timing Explorer (RXTE) monitoring of 1E 1048.1$``$5937 in which we have attempted long-term phase-coherent timing like that achieved for other AXPs kcs99 . The results described here are reported in more detail in kgc+01 .

## Observations and Results

The observations we report on were made with RXTE’s Proportional Counter Array (PCA; jsg+96 ). Observations of 3–6 ks in length of 1E 1048.1–5937 were made on a monthly basis during 1996 November – 1997 December and 1999 January – 2000 August. In addition, we used archival observations from 1996; these generally had longer integration times than the other data sets. To minimize use of telescope time, our monitoring data consist of brief (usually 3 ks) snapshots of the pulsar. These snapshots suffice to measure pulse arrival times for a phase-coherent timing analysis to good precision. However, for any one epoch, the measured period has typical uncertainty $`3`$ ms, quite large by normal timing standards. Thus, our snapshot method of measuring pulse arrival times can determine spin parameters with extremely high precision only when phase coherence can be maintained. For details regarding this timing procedure, see kgc+01 . The snapshot observations are always, however, useful for monitoring the source pulsed flux and pulse morphology (see below).

### Timing

We maintained unambiguous phase coherence for 1E 1048.1$``$5937 in our monthly observations from 1999 January 23 through 1999 November 15. We required a fourth-order polynomial to characterize the 17 pulse arrival times obtained in this span. These results alone clearly imply that the rotational behavior of 1E 1048.1$``$5937 is quite different from that of AXPs 1E 2259+586 and RXS J170849.0$``$400910. Those AXPs exhibit much more stable rotation on comparable and even longer time scales, that is, terms of higher order than $`\dot{\nu }`$ are very small or negligible for those pulsars on time scales of over a year kcs99 . The span 1999 January through November represents the longest over which we can phase-connect timing data from 1E 1048.1$``$5937. Investigating archival RXTE data going back to 1997 for 1E 1048.1$``$5937, we find timing results that are similar those obtained in our recent monitoring program, namely we are able to maintain phase coherence only over few-month intervals.

We can compare our pulse ephemerides with measurements of pulse frequency made over the past 20 yr in order to look for long-term trends. Figure 1 shows the spin history of 1E 1048.1$``$5937 with previously measured spin frequencies plotted as points with their corresponding 1$`\sigma `$ error bars. Data were taken from a variety of sources opmi98 ; pkdn00 ; bss+00 . Our RXTE timing results are plotted as lines representing separate, short, phase-connected segments. The dotted line represents an extrapolation of the $`\nu `$ and $`\dot{\nu }`$ from the 1999 coherent fit. The lower plot shows the same data set with the linear term subtracted off. This ephasizes deviations from the simple linear trend.

### Pulsed Flux and Pulse Morphology

In accreting systems in which the neutron star is undergoing spin-up, changes in torque should be correlated with changes in X-ray flux. AXPs are spinning down, however. Chatterjee et al. (2000) suggest that AXPs might be spinning down in the propeller regime due to accretion from a fall-back disk. In that case, although the physics of the propeller regime is not well understood, it is still likely that $`L_x`$ should be correlated with torque kgc+01 .

Given the large field-of-view of the PCA and that the bright, nearby, but unrelated source $`\eta `$ Carinae exhibited large flux changes over the course of our observations, direct flux measurements of 1E 1048.1$``$5937 could not be made with our RXTE data. Instead, we determined the pulsed component of the flux by using off-pulse emission as a background estimator. This renders our analysis insensitive to changes in the fluxes of other sources in the field-of-view.

The results are shown in Figure 2. We find no large pulsed flux variations. The $`\chi ^2`$ strictly speaking does suggest some low-level variability; longer individual observations are clearly necessary to verify this is the case. However, as we discuss below, the pulsed flux is certainly much more stable than previous analyses have suggested opmi98 .

We have also used the RXTE data to search for pulse profile changes, as many accretion-powered pulsars exhibit significant changes in their average pulse profiles. Such changes can be correlated with the accretion state, and hence accretion torque and timing behavior bcc+97 . Furthermore, X-ray pulse profiles from the SGRs 1806$``$20 and 1900+14 have shown differences at different epochs depending on time since outburst kds+98 ; ksh+99 . However, we find no significant changes in the pulse profile morphology in any of the RXTE observations of 1E 1048.1$``$5937.

## Discussion and Conclusions

Long-term RXTE monitoring of the AXP 1E 1048.1$``$5937 has shown it to be a much less stable rotator than other AXPs, yet its pulse profile and pulsed flux are stable. Previously, Oosterbroek et al. (1998) compiled flux data from a variety of different X-ray instruments that observed 1E 1048.1$``$5937. That compilation suggested that the pulsar shows variability by over a factor of $``$5 on time scales of a few years. The reality of those flux changes is not supported by our results. One caveat is that we measure pulsed flux, while they report flux, so the results could be reconciled if the pulsed fraction is variable.

In the context of the magnetar model, we note that the timing behavior of 1E 1048.1$``$5937 is somewhat similar to that observed for the soft gamma repeaters SGR 1806$``$20 and 1900+14 mrl99 ; wkf+00 ; wkp+00 . However, as the stable flux time series (Fig. 2) for the AXP shows, it has not undergone any outbursts. This can perhaps be understood in terms of persistent seismic activity and small-scale crustal fractures tb98  or low amplitude toroidal modes resulting in angular momentum loss following crustal twisting fractures dun00 .

1E 1048.1$``$5937 is unusual among AXPs for reasons other than just its timing behavior. In particular, it shows the highest ratio of blackbody to total flux (once energy band is accounted for), and the largest pulsed fraction. In addition, it has the lowest photon index for the power-law tail in its spectrum of any AXP, which makes it the closest to those measured in the X-ray band for SGRs 1806$``$20 and 1900+14. Further, the thermal component of 1E 1048.1$``$5937’s spectrum has the highest temperature (0.64 keV) of any AXP. This temperature is comparable to that seen for SGR 1900+14 post-burst, 0.62 keV wkp+00 . It therefore could be the case that 1E 1048.1$``$5937 is a transition object between the populations of AXPs and SGRs, and the AXP most likely to one day undergo an outburst.

In the context of accretion models, perhaps the best source with which to compare 1E 1048.1$``$5937 is 4U 1626–67, a 7.7 s accreting pulsar in a 42-min binary with a low-mass companion. Although 1E 1048.1$``$5937 is noisy by AXP timing standards, its noise is comparable in strength to that of 4U 1626$``$67 cbg+97 . Still, we regard the case for 1E 1048.1$``$5937 as an accreting binary, even with a very low-mass companion, as weak, given the other evidence against this hypothesis, namely, its much softer spectrum than other accreting binaries, the absence of pulsed flux or pulse morphology changes correlated with the timing behavior, and the spin-down over some 20 yr. It is more difficult to dismiss the possibility that 1E 1048.1$``$5937 is accreting from a supernova fall-back disk, since there is not yet a consensus on the properties such a disk would have or on the expected timing and variability properties of the pulsar. However, one expectation is that such a disk would be a significant emitter in the optical and infrared, Future optical/IR observations following a more precise localization using the Chandra X-ray Observatory could test the fallback disk model.