text:
# Dark matter and visible baryons in M33

## 1 Introduction

M33 is a low luminosity spiral galaxy in the Local Group. Its large angular extent and well determined distance make it ideal for a detailed study of the radial distribution of visible and dark matter. The interplay between the dark and visible matter as well as the nature and the radial distribution of dark matter are still open crucial issues with strong implications for galaxy formation and evolution theories. M33 has three advantages over its brighter neighbor M31: (i) being of lower luminosity, it is a dark matter dominated spiral in which it is easier to disentangle the dark and luminous mass components (Persic, Salucci $`\&`$ Stel 1996, Corbelli $`\&`$ Salucci 2000 hereafter CS); (ii) it has no prominent black hole or bulge component in the center (Gebhardt et al. 2001) and therefore it can be used to test if the innermost dark matter density profile is consistent with cosmological models which predict a central density cusp; (iii) being at lower declination, high sensitivity 21-cm observations with the Arecibo radio telescope made possible a velocity field map out to large galactocentric radii. CS have used this velocity field map, together with aperture synthesis data of Newton (1980) and with the tilted ring model results of Corbelli $`\&`$ Schneider (1997), to derive the rotation curve of M33 out to 13 disk scalelengths. This unveiled large scale properties of the mass distribution in M33, such as the lack of correlation between the dark and visible matter, the dark halo mass and the possible radial slope of the dark matter density. However some ambiguities are left, especially in the inner regions, where the low spatial resolution of 21-cm observations and the possible role of a stellar nucleus or of the molecular gas component have not allowed to establish unambiguously the steepness of the dark matter density profile.

Several authors, e.g. Bolatto et al. (2002), have shown the advantages of using CO observations for determining the central densities of dark matter halos. A full map of the CO emission in M33 would not only provide higher spatial resolution data on the velocity field but it would also give a first estimate of the molecular mass distribution in this galaxy. The very blue color of M33 and the presence of prominent HII regions indicate vigorous star formation activity in the nuclear region and along the spiral arms: here the molecular gas component might represent a substantial fraction of the total disk mass. CO observations of the innermost 500 pc of M33 (Wilson $`\&`$ Scoville 1989), have confirmed in fact that here the molecular gas has a higher surface density than the atomic gas. Recently Heyer, Corbelli $`\&`$ Schneider (2003) have fully mapped the CO emission in M33 using the FCRAO 14-m telescope. A detailed analysis of the molecular gas distribution and emission features will be presented elsewhere. In this paper we derive only the average radial distribution of the molecular gas mass (Section 2) and analyze the kinematics of the CO gas (Section 3). The implications for the dark matter halo models are shown in Section 4. An accurate knowledge of the velocity field and of the mass surface densities over a large disk area are useful for understanding the processes which trigger the star formation in this galaxy, such as disk instabilities, spiral arms, cooling and heating mechanisms. A possible link between the radial extent of the observed star formation region and the gravitationally unstable region in the disk is discussed in Section 5.

## 2 The molecular data

The nearby spiral galaxy M33 has been fully mapped in the CO J=1-0 transition with a velocity resolution of 0.81 km s-1 and a spatial resolution of 45 arcsec using the FCRAO 14-m telescope (Heyer, Corbelli $`\&`$ Schneider 2003). A 22.5 arcsec spatial sampling is obtained over an area of about 900 arcmin2 centered on RA 01h33m50s.89 and DEC 3039’36”.7. For a distance to M33 of 0.84 Mpc (Freedman et al. 2001), which will be used throughout this paper, this means that the total disk area observed is about 44 kpc2 at a spatial resolution of 165 pc. Here we shall use the data at the original spatial and spectral resolution to derive the kinematical information on the molecular gas. We analyze in detail each spectrum by assigning a spectral window $`\mathrm{\Delta }V=V^fV^i`$ to any detectable emission feature. There are a few cases where more than one spectral window is needed but secondary signals are always quite low. If no emission is clearly present in a given spectrum, a spectral window is set based on the diffuse emission detectable when averaging nearby positions. We compute the integrated brightness temperature $`F`$ in K km s-1 by integrating over $`\mathrm{\Delta }V`$ the telescope temperature observed at each position in the sky, and dividing this number by the main beam efficiency of the telescope (which is 0.45 at 115 GHz). In order to derive the average radial distribution of the molecular gas we need a concentric ring model. 45 concentric rings of width 150 pc cover all the M33 surface area observed. They are tilted as to best fit the neutral hydrogen distribution (Corbelli $`\&`$ Schneider 1997). Inside this area the position angle of the best fitting ring model changes only slightly, and therefore we consider it constant and equal to 22o. The inclination of the rings increases from 50o to 58o moving radially outwards, as described in details in Section 3 of CS. We use the concentric ring model for assigning a face on radial distance to each observed position, and for deriving the inclination corrected average flux in each ring. Following Wilson’s results (Wilson 1995) we use the standard conversion factor between the CO luminosity and the line of sight column density of $`H_2`$ ($`S_{H_2}=2.8\times 10^{20}\times F`$ cm-2 = 4.5$`\times F`$ M pc-2).

The resulting average surface density of molecular hydrogen, perpendicular to the plane of the galaxy (i.e. corrected for the disk inclination respect to our line of sight), is shown in Figure 1. The estimated total mass of molecular gas is $`2\times 10^8`$ M. This is a factor 6 higher than measured by Wilson $`\&`$ Scoville (1989) in the smaller nuclear region. The molecular hydrogen mass is about 10$`\%`$ the neutral hydrogen mass ($`2.2\times 10^9`$ M for D=0.84 Mpc). The surface density distribution does not peak at the center but shows a peak at a few hundred parsecs from the nucleus, as noted previously by Wilson $`\&`$ Scoville (1989). The most prominent CO emission features are located in the innermost 1 kpc region and along the spiral arms. This can be seen already in the 1.3 $`\times `$ 7 kpc2 preliminary emission map (Corbelli 2000) oriented along the M33 major axis. An exponential is fitted to the azimuthally averaged surface density data from 0 to 6 kpc. We leave out from the fit the last 6 points shown in Figure 1 which seem to indicate a plateau at the level of 0.6 M pc-2. In reality we don’t have enough sensitivity and coverage in our survey to determine the slope of the $`H_2`$ distribution so far out. The best fit to S$`_{H_2}`$ is shown by the straight line in Figure 1 and is given by

$$S_{H_2}=6.3\times \text{exp}(\frac{R}{2.5\text{ kpc}})\text{M}_{}\text{pc}^2$$
(1)
Therefore the CO radial scalelength in M33 is larger than the B and K band scalelengths (1.9 and 1.3 kpc respectively, Regan $`\&`$ Vogel 1994). In Figure 2 we show the HI as well as the total (HI + H2) gas distribution in this galaxy. The total gas surface density can be described by a gaussian profile with half width at half maximum of about 6 kpc. Note that the bump of the spiral arm feature around 2 kpc is visible in the azimuthal averages of the molecular and atomic gas shown in Figure 2.

## 3 The CO rotation curve

As Wilson $`\&`$ Scoville (1989) have already pointed out from a comparison of interferometric and single dish data, small clouds and diffuse emission in M33 contain a non negligible fraction of the total molecular mass. Half of the molecular gas mass that we estimate in M33 resides in fact in regions with integrated flux less than 3$`F_\sigma `$ 1.33 K km s-1. $`F_\sigma `$ is the integrated noise defined in terms of the spectral window $`\mathrm{\Delta }V`$ as:

$$F_\sigma =\sigma _{rms}\times \sqrt{\mathrm{\Delta }V\times 0.81}$$
(2)
where $`\sigma _{rms}`$ is the noise in the spectrum at the original resolution of 0.81 km s-1. We fit single or double gaussians to all signals with integrated flux $`F>F_{5\sigma }`$ with a success rate of 75$`\%`$. There is no significant shift between the gaussian and the flux weighted mean velocities (difference is $`\stackrel{<}{}2`$ km s-1 and is never above 5 km s-1). In deriving kinematical informations we shall consider only spectra with $`F>5F_\sigma `$, well fitted by gaussian profiles, from positions at angles $`\alpha <45^{}`$ from the major axis. These conditions restrict the analysis to a smaller radial range than observed. Luckily the prominent CO features come from locations close to the major axis and therefore we do have many points to average for building up a rotation curve. Corrections for the inclination of the disk and for the angle $`\alpha `$ between the spectrum position and the major axis are applied according to the tilted ring model described in the previous Section. Radial bin averages of the flux weighted mean velocities are plotted in Figure 3, also for the northern and southern side of the galaxy separately. Error bars in each bin show the dispersion of the velocities about the mean.

In the innermost 1 kpc the south and north rotation curves coincide reasonably well. Beyond this radius there are some differences both in the dispersions and in the mean values of $`V`$. Velocities in the southern half of M33 have higher dispersion than in the northern half. Also the southern side shows more clearly wiggles which can be interpreted as signs of density wave perturbations associated with the spiral arms. The detection of such features in the northern side is more doubtful. Residual velocities, obtained by subtracting the average rotational velocity $`V(R)`$ from the rotational velocity measured in each spectrum, $`V_{obs}(R)`$, are shown in Figure 4. In the left hand panel of Figure 4 the map of residual velocities is superimposed to the used tilted ring model, after rotating the galaxy by 22o clockwise in the plane of the sky. Cross and square symbols are for positive and negative residuals respectively, and their size is proportional to the amplitude of the residual velocities. Since the approaching side of the galaxy is the northern side ($`y>0`$), a positive residual there means that the local velocity is less than the average rotational velocity. For each quadrant of the galaxy, the right hand panels of Figure 4 show the residual velocities as function of the galactocentric radius. It hard to judge whether there are signs of non circular motions and spiral density waves in the velocity residuals map. Part of the difficulty is due to the fact the bulk of the CO emission comes from localized regions and is not smoothly distributed throughout the disk. A better attempt can be made in the future with the use of a numerical simulation for any given model assumption.

Comparison of the rotation curve in Figure 3 with the 21-cm data of Newton (1980) enlightens a discrepancy in the northern inner region, where at $`R\stackrel{<}{}0.5`$ kpc the HI velocities are lower by about 10 km s-1 than the CO velocities. This might be due to the lack of atomic gas in the innermost region, which together with the larger 21-cm beam, limits the accuracy of the atomic line measure. Finite angular resolution is not affecting significantly the velocities measured in our CO survey. To estimate any possible “beam smearing” effect we assume that the true velocity, for a given mass model, is that resulting from rotation curve and perform a numerical simulation to derive the observed velocities (van den Bosch $`\&`$ Swaters 2001). The observed velocities at each position on the plane of the sky are the flux weighted means of the rotational velocities along the line of sight, corrected for the assumed ring model and convolved with the shape of the telescope beam. We then compare the binned data with equally binned values of the velocities resulting from the simulation (we exclude points which form an angle $`\alpha >45^{}`$ with the major axis). Differences between the observed and the true velocities, for mass models which give a good fit to the data, are always smaller than 0.6 km s-1 except in the innermost 0.5 kpc where the beam smearing effect can lower the observed velocity as much as 2 km s-1.

## 4 The mass components and constrains on cosmological halo models

Rotational velocities, which we use to derive the visible and dark matter distributions in M33, are shown in Figure 5. The innermost point at $`R150`$ pc is from H$`\alpha `$ measurements of Rubin $`\&`$ Ford (1985). Our CO data extending over the interval $`0.2R<3.4`$ kpc have been averaged into the subsequent 11 bins. For $`3.4R<5.5`$ kpc data are taken from 21-cm observations of Newton (1980) at the Cambridge Half-Mile telescope while for larger radii we use Arecibo 21-cm data (CS and references therein). The plotted error bars in each bin are twice the dispersion of the velocities about the mean.

We consider the dynamical contribution of two visible mass components: stars, distributed in a disk and in the nucleus, and gas, in molecular and in neutral atomic form. For the stellar disk we consider that the mass follows the light distribution in the K-band, well fitted by an exponential with scalelength $`R_d1.3\pm 0.2`$ kpc (Regan $`\&`$ Vogel 1994). After subtracting the exponential disk from the surface brightness, a central emission excess remains at $`R\stackrel{<}{}0.5`$ kpc (Bothun 1992, Regan $`\&`$ Vogel 1994) which has been attributed to spiral arms, to an extended semistellar nucleus or to a diffuse halo or bulge. We shall refer to this excess as the “nucleus” of M33 despite the fact that in the literature the word “nucleus” is often used for the light excess observed only at $`R\stackrel{<}{}10`$ pc (Kormendy $`\&`$ McClure 1993). Nuclear rotational velocities have been observed in H$`\alpha `$ along the major axis between 20 and 200 pc from the center (Rubin $`\&`$ Ford 1985), and seem roughly constant at about $`27`$ km s-1. This flatness of the rotation curve is quite puzzling and needs further investigation once uncertainties on the nuclear rotational velocities will be established together with a knowledge of the velocity field along different directions. In this paper we model the rotation curve at $`R\stackrel{>}{}150`$ pc only, considering two different possibilities for the nuclear component $`V_n`$ at these radii:

$`(1)`$ $`V_nR^{0.5}/(R+a)`$. This rotational velocity approximates the dynamical contribution of a spheroidal mass distribution whose surface density can be fitted by a de Vaucouleurs $`R^{1/4}`$ law (Hernquist 1990). This scaling law has been used by Regan $`\&`$ Vogel (1994) to fit the J-band nuclear surface brightness. The compatibility of the fit with photometric data in other bands (e.g. Bothun 1992) is somewhat uncertain, as is the value of the scalelength $`a`$. Suggested values of $`a`$ lie in the range $`0.1a1`$ kpc, which we will consider in this paper.

$`(2)`$ $`V_nR^{0.115}`$, as if core collapse has happened or is happening in the compact stellar nucleus (Cohn 1980, Kormendy $`\&`$ McClure 1993). The corresponding power-law mass density distribution develops up to a limiting outer radius $`R_c`$, which should also be given. For $`R>R_c`$, $`V_nR^{0.5}`$. We consider as possible $`R_c`$ values the radius of each data point of the rotation curve inside the optical region i.e. $`R_c<8`$ kpc. We consider also the case $`R_c150`$ pc, which is appropriate if most of the nuclear mass resides at very small radii and the nucleus can be treated as a point source throughout the extent of the present rotation curve.

Hereafter, we will refer to these models as nuclear model (1) and (2) respectively. It has been established by CS that the visible baryons in M33 are only a small fraction of the total galaxy mass, and that a dark matter halo is in place. We will use the M33 rotation curve presented here to test in details the consistency of the required halo density profile with theoretical models which predict a well defined dark matter distribution around galaxies. Namely, in addition to the widely used non-singular isothermal dark halo model we shall consider a spherical halo with a dark matter density profile as originally derived by Navarro Frenk and White (1996, 1997, hereafter NFW) for galaxies forming in a Cold Dark Matter scenario (hereafter CDM). We consider also the Burkert and the Moore dark matter density profile (Burkert 1995, Moore et al. 1998). All models have two free parameters to be determined from the rotation curve fit. The density distribution for a non-singular isothermal dark halo can be written as:

$$\rho (R)=\frac{\rho _{iso}}{1+(\frac{R}{R_{iso}})^2}$$
(3)
The density profile proposed by Burkert (1995) has also a flat density core:

$$\rho (R)=\frac{\rho _B}{(1+\frac{R}{R_B})\left(1+(\frac{R}{R_B})^2\right)}$$
(4)
and a strong correlation between the two free parameters $`\rho _B`$ and $`R_B`$ has been found by fitting the rotation curve of low mass disk galaxies.

The NFW density profile is:

$$\rho (R)=\frac{\rho _{NFW}}{\frac{R}{R_{NFW}}\left(1+\frac{R}{R_{NFW}}\right)^2}$$
(5)
Numerical simulations of galaxy formation in CDM models find a correlation between $`\rho _{NFW}`$ and $`R_{NFW}`$ which depends on the cosmological model (Avila-Reese et al. 2001, Eke Navarro and Steinmetz 2001). Often this correlation is expressed using the concentration parameter $`C_\mathrm{\Delta }R_\mathrm{\Delta }/R_{NFW}`$ and $`M_\mathrm{\Delta }`$ or $`V_\mathrm{\Delta }`$. $`R_\mathrm{\Delta }`$ is the radius of a sphere containing a mean density $`\mathrm{\Delta }`$ times the critical density, $`V_\mathrm{\Delta }`$ and $`M_\mathrm{\Delta }`$ are the characteristic velocity and mass inside $`R_\mathrm{\Delta }`$. We use a Hubble constant $`H_0=65`$ km s-1 Mpc-1 and show results for $`\mathrm{\Delta }=100`$, as in CDM halo models of Eke Navarro and Steinmetz (2001), unless otherwise specified.

A similar profile but with a different central slope has been suggested instead by numerical simulations of Moore et al. (1998) for galaxies forming in CDM scenarios:

$$\rho (R)=\frac{\rho _M}{(\frac{R}{R_M})^{1.5}\left(1+(\frac{R}{R_M})^{1.5}\right)}$$
(6)
For each dark halo model and each nuclear model we shall determine five free parameters applying the least-square method to the observed rotation curve. These free parameters are: the dark halo core radius $`R_{iso,NFW,M,B}`$, the dark matter density $`\rho _{opt}`$ computed at $`R=R_{opt}3.2R_d4.1`$ kpc, the ratio between the stellar disk mass $`M_d`$ and the blue luminosity of the galaxy $`L`$ ($`L5.7\times 10^9`$ L), the nuclear mass $`M_n`$, and the nuclear scale length $`a`$ or the cut-off radius $`R_c`$. The 1-$`\sigma `$, 2-$`\sigma `$ and 3-$`\sigma `$ ranges are determined using the reduced $`\chi ^2`$ and assuming Gaussian statistics. They are computed by determining in the free parameters space the projection ranges, along each axis, of the hypersurfaces corresponding to the 68.3$`\%`$, 95.4$`\%`$ and 99.7$`\%`$ confidence levels.

Both the non-singular isothermal dark halo model and the Burkert dark matter density profile give good fits to the rotation curve since the minimum reduced $`\chi ^2`$ is 0.7 for the nuclear model (1) and 1.2 for the nuclear model (2). The best fitting values of the free parameters are: $`a=0.4`$ kpc, $`M_n=3\times 10^8`$ M, $`M_d/L=1`$M/L, $`R_{iso}=7`$ kpc or $`R_B=12`$ kpc, $`\rho _{opt}=5\times 10^{25}`$ g cm-3 for the nuclear model (1), and $`R_c=0.7`$ kpc, $`M_n=10^8`$ M, $`M_d/L=1`$ M/L, $`R_{iso}=8`$ kpc or $`R_B=13`$ kpc, $`\rho _{opt}=5\times 10^{25}`$ g cm-3 for the nuclear model (2). The 3-$`\sigma `$ ranges for the free parameters of both the isothermal and the Burkert model are:

$$0.8\frac{\text{M}_d/\text{L}}{\text{M}_{}/\text{L}_{}}1.1;5\times 10^7\frac{\text{M}_n}{\text{M}_{}}8\times 10^8;4\times 10^{25}\frac{\rho _{opt}}{\text{g cm}^3}7\times 10^{25};$$

$$0.4\frac{R_c}{\text{kpc}}<8;7\frac{R_B}{\text{kpc}}19\text{or}4\frac{R_{iso}}{\text{kpc}}13.$$

The resulting 3-$`\sigma `$ range for the nuclear scale length $`a`$ in the nuclear model (1) covers the whole range we have considered i.e. $`0.1a1`$ kpc. For the nuclear model (2) instead we find a lower limit to $`R_c`$. This means that the nucleus cannot be considered a point source for any possible dark halo models we are examining. This ‘dynamical’ conclusion agrees with photometric data (Bothun 1992, Minniti, Olszewski $`\&`$ Rieke 1993, Regan $`\&`$ Vogel 1994) showing an excess of light inward of 0.5 kpc after subtracting the exponential disk. Unfortunately estimates of the nuclear luminosity $`L_n`$ are still quite uncertain: in the B-band $`7\times 10^7L_{n,B}3\times 10^8`$ L⊙,B (Bothun 1992, Regan $`\&`$ Vogel 1994), and we cannot limit the nuclear stellar mass range further based on possible unrealistic values of the mass to light ratio (Bell $`\&`$ de Jong 2001). An extended nucleus is not required instead if the density of the dark matter scales as $`R^{1.3\pm 0.1}`$ throughout the halo; these type of solutions has been investigated by CS and provides also good fits to the rotation data presented in this paper. Both the Burkert and the isothermal model require a very extended core region of nearly constant dark matter density. Figure 5$`(a)`$ shows the best fit using the Burkert profile with the nuclear model (1). The 1-$`\sigma `$ and 3-$`\sigma `$ probability contours for $`a=0.4`$ kpc in the $`\rho _{opt}R_B`$ plane are displayed in Figure 6$`(a)`$. The correlation found by Burkert between $`\rho _B`$ and $`R_B`$ implies a relation between $`\rho _{opt}`$ and $`R_B`$ shown with crosses in Figure 6$`(a)`$. We can see that the M33 acceptable parameters lie close to, but outside, the Burkert correlation line. So either the found scaling relation between $`\rho _{opt}`$ and $`R_B`$ should be slightly modified to include M33, or it changes for values of $`R_B`$ larger than those considered by Burkert. A Burkert model which fits M33 is indistinguishable from an isothermal model since both require a large constant density core region and, even at the outermost fitted radius, the dark matter density does not decline yet as an $`R^2`$ or $`R^3`$ power law. The dark matter density profile suggested by Moore fails to fit the observed rotation curve in M33.

For the NFW dark matter density profile for $`\mathrm{\Delta }=100`$ or 200 the minimum $`\chi ^2`$ is 0.7 for the nuclear model (1) and 1.0 for the nuclear model (2). The best fitting values of the free parameters for the nuclear model (1) are: $`a=0.4`$ kpc, $`M_n=2\times 10^8`$ M, $`M_d/L=0.8`$ M/L, $`R_{NFW}=140`$ kpc, $`\rho _{opt}=6\times 10^{25}`$ g cm-3. For the nuclear model (2) $`R_c=0.7`$ kpc, $`M_n=8\times 10^7`$ M, $`M_d/L=0.8`$ M/L, $`R_{NFW}=300`$ kpc, $`\rho _{opt}=6\times 10^{25}`$ g cm-3. The 3-$`\sigma `$ ranges for the free parameters are:

$$0.5\frac{\text{M}_d/\text{L}}{\text{M}_{}/\text{L}_{}}0.9;3\times 10^7\frac{\text{M}_n}{\text{M}_{}}4\times 10^8;5\times 10^{25}\frac{\rho _{opt}}{\text{g cm}^3}8\times 10^{25};$$

$$0.3\frac{R_c}{\text{kpc}}<8;\frac{R_{NFW}}{\text{kpc}}>15$$

The corresponding 3-$`\sigma `$ upper limit for the concentrations is $`C_{100}<12.5`$ for nuclear model (1) and $`C_{100}<11`$ for nuclear model (2), while the corresponding best fits to the rotation curve require $`C_{100}4`$ and $`C_{100}2.5`$ respectively. For all the dark halos we are examining the nuclear model (2) provides a flatter rotation curve than model (1) inward of the center-most data point fitted, in agreement with the preliminary measures of Rubin $`\&`$ Ford (1985). The central steepness of the NFW dark matter density profile lowers the nuclear stellar mass required to fit the rotational data. Figure 6$`(b)`$ shows the 1-$`\sigma `$ and 3-$`\sigma `$ areas in the $`M_{100}C_{100}`$ plane. The dashed and dashed-dotted lines in the same figure are the fits to numerical simulation results for structures formation in the cosmological models $`W_8`$ and $`S_{0.9}`$ respectively of Eke, Navarro $`\&`$ Steinmetz (2001). Both model assume matter density parameters $`\mathrm{\Omega }_0=0.3`$ and $`\mathrm{\Lambda }_0=0.7`$, a power spectrum in the form given by Bardeen et al. (1986) with $`\sigma _8`$=0.9 and shape parameter $`\mathrm{\Gamma }=0.2`$. $`S_{0.9}`$ corresponds to this ‘standard’ $`\mathrm{\Lambda }`$CDM spectrum, while $`W_8`$ mimics a Warm Dark Matter power spectrum because in addition its power is reduced on scales smaller than that of the characteristic free streaming mass $`M_f=8\times 10^{11}`$ M. It can be easily seen in Figure $`6(b)`$ that the possible $`C_\mathrm{\Delta }`$ and $`M_\mathrm{\Delta }`$ values for M33 are consistent with a Wark Dark Matter cosmology but lie below the relation found in simulated standard CDM halos. These $`M_\mathrm{\Delta }C_\mathrm{\Delta }`$ relations have however an associated scatter. Even a small scatter, $`\sigma `$(lg $`C_\mathrm{\Delta }`$)=0.1 as suggested for galaxies which did not undergo major mergers (Wechsler et al. 2002), makes the M33 dark halo still compatible with standard $`\mathrm{\Lambda }`$CDM galaxy formation models. We display in Figure 5$`(b)`$ the fit to the M33 rotation curve using a NFW dark halo profile with $`C_{100}=8`$, compatible with a standard $`\mathrm{\Lambda }`$CDM and a Wark Dark Matter cosmology. The nuclear model (2) is used. The reduced $`\chi ^2`$ and the values of the free parameters found from the best fit are: $`\chi ^2=1.5`$, $`R_c=0.7`$ kpc, $`M_n=7.5\times 10^7`$ M, $`M_d/L=0.7`$ M/L, $`R_{NFW}=35`$ kpc, $`\rho _{opt}=7\times 10^{25}`$ g cm-3.

## 5 Surface densities of visible matter and the global stability of the disk

In the previous Section we have determined the dark matter and stellar contribution to the M33 rotation curve. Despite the fact that the radial scaling law of the dark matter density required for fitting the rotation curve is not unique, the total baryonic and dark mass up to the last measurable point (R=17 kpc) are determined with small uncertainties. The total stellar mass in the M33 disk is estimated to be 3-6$`\times `$ 109 M. The mass of the dark halo at the outermost observable radius is about a factor 10 higher than the stellar disk mass. The neutral hydrogen plus the molecular gas mass is 2.4$`\times 10^9`$ M. We don’t have a firm estimate of the nuclear stellar mass $`M_n`$, but all three possible dark matter models discussed in this paper require an extended distribution. Its mass lies in the range $`5\times 10^7M_n8\times 10^8`$ M if the dark halo has a large size core of constant density, while $`3\times 10^7M_n4\times 10^8`$ M if the dark halo has a central cusp as in NFW model. We can estimate the total surface density of visible baryons in M33 neglecting at the moment the small contribution from the ionized gas. Figure 7 shows the total baryonic surface density computed by adding the helium to the molecular and the neutral atomic gas (a 1.33 factor correction), and adding to this the nuclear and stellar disk surface density. We use the mass model displayed in Figure 5$`(b)`$, for which we show also the dark matter surface density within 0.5 kpc from the plane of the galaxy. We extrapolate the fits of the stellar and molecular gas surface density out to 16 kpc but it is likely they experience a truncation at smaller radii. For the HI surface density only we know that it declines sharply beyond 15-16 kpc due to ionization effects by an extragalactic background (Corbelli $`\&`$ Salpeter 1993). The gas distribution is likely to extend further out but the dominant phase becomes ionized atomic hydrogen. Assuming a gaseous vertical extent of 0.5 kpc above and below the galactic plane, we can say that the surface density of visible baryons falls below the dark matter surface density contained in the M33 disk for $`R\stackrel{>}{}10`$ kpc. A similar conclusion applies also for the Burkert best fit model.

Using our accurate measurement of the rotational velocity gradient and of the surface density of the baryonic and non baryonic components, we can check if the onset of massive star formation in M33 is strongly related to a global instability of the disk. One of the most popular tools used to assess the gravitational stability of gaseous disks is the Toomre stability criterion:

$$f\mathrm{\Sigma }_g<\mathrm{\Sigma }_T\frac{\kappa c_g}{\pi G}$$
(7)
where $`\mathrm{\Sigma }_g`$ is the gas surface mass density, $`f`$ is the disk thickness factor, $`\kappa `$ is the epicyclic frequency (Toomre 1964), G is the gravitational constant, $`c_g`$ is the gas velocity dispersion. We shall use both $`f=1`$ as for a thin disk and $`f=0.7`$ to mimic instead a disk with a finite thickness (Romeo 1992). The criterion described by equation (7) is not at all related to the star formation efficiency but it correctly predicts the occurrence of a star formation threshold radius which, in many spiral galaxies, has been determined from the drop of the H$`\alpha `$ surface brightness (e.g. Martin $`\&`$ Kennicutt 2001). In M33 Kennicutt (1989) observed such a drop at $`R_{}6.4`$ kpc. However, since the early work of Kennicutt it is well known that the simple Toomre stability criterion fails to account for the active star forming region in the disk of low mass spirals such as M33 (Elmegreen 1993, van Zee et al. 1997, Wong $`\&`$ Blitz 2002). An alternative stability criterion includes the shear rate of the disk and it may be more appropriate for irregular galaxies (Elmegreen 1993, Hunter Elmegreen $`\&`$ Baker 1998). It is based on the consideration that making new stars is possible only when there is enough time for the instability to grow, before shear stretches the perturbation apart. This corresponds to the following condition:

$$f\mathrm{\Sigma }_g<\mathrm{\Sigma }_A\frac{2.5Ac_g}{\pi G}$$
(8)
where the Oort constant $`A=0.5Rd\mathrm{\Omega }/dR`$ is the local shear rate. In Figure 8 we display $`\mathrm{\Sigma }_T`$ and $`\mathrm{\Sigma }_A`$, using $`c_g=6`$ km s-1. This value is at the lowest end of the cold gas velocity dispersion range estimated in M33 (6-9 km s-1, Warner at al. 1973, Thilker 2000). The thin and thick continuous lines in Figure 8 show the mass surface densities for $`f=1`$ and for $`f=0.7`$ respectively. The arrow indicates $`R_{}`$. We can see in Figure 8$`(a)`$ that at $`RR_{}`$ the azimuthally averaged surface gas density for a thin disk falls below the threshold for a gravitational instability to develop according to the shear rate criterion. We shall investigate now if considering additional sources of gravity in the disk, a modified version of the classical Toomre criterion can still give correct predictions for the star formation threshold radius.

Since the masses of stars and gas in the M33 disk are comparable, we are encouraged to use a stability parameter for a two components fluid. Following Wang $`\&`$ Silk (1994) we define $`\mathrm{\Sigma }_{}\mathrm{\Sigma }_s\times c_g/c_s`$, where $`\mathrm{\Sigma }_s`$ and $`c_s`$ are the surface density and velocity dispersion of the stars in the disk. Using the epicyclic frequency or the shear rate of the disk, the stability condition reads:

$$f(\mathrm{\Sigma }_g+\mathrm{\Sigma }_{})<\mathrm{\Sigma }_T,\mathrm{\Sigma }_A$$
(9)
Figure 8$`(b)`$ shows $`f(\mathrm{\Sigma }_g+\mathrm{\Sigma }_{})`$, using $`M_d/L=0.7`$ M/L for deriving $`\mathrm{\Sigma }_{}`$ (as for the mass model shown in Fig. 5$`(b)`$), and $`c_g/c_s=0.3`$. For $`c_g=6`$ km s-1 this last ratio implies $`c_s=20`$ km s-1, a value observed in the central regions of M33 (Kormendy and McClure 1993, Gebhardt et al. 2001). We can see that with a value of $`M_d/L`$ slightly higher than 0.7, the Toomre criterion for a thin disk predicts correctly a gravitationally unstable disk for $`R<R_{}`$. If finite thickness effects are taken into account or if $`c_g`$ is higher than 6 km s-1 at $`RR_{}`$ we need the shear rate criterion.

Following Hunter, Elmegreen $`\&`$ Baker (1998) we consider also the possible influence of dark matter gravity. If $`\mathrm{\Sigma }_{dm}`$ is the dark matter surface density within the vertical extent of the gas, we can rewrite equation (9) as

$$f(\mathrm{\Sigma }_g+\mathrm{\Sigma }_{dm}+\mathrm{\Sigma }_{})<\mathrm{\Sigma }_T,\mathrm{\Sigma }_A$$
(10)
In Figure 8$`(c)`$ we show $`f(\mathrm{\Sigma }_g+\mathrm{\Sigma }_{dm}+\mathrm{\Sigma }_{})`$ assuming a gaseous extent of 0.5 kpc above and below the galactic plane and a NFW dark matter profile as used for the rotation curve fit displayed in Figure 5$`(b)`$. In this case finite thickness effects, or equivalently a velocity dispersion as high as 9 km s-1, make the disk gravitationally unstable, according to the Toomre criterion, exactly for $`R<R_{}`$. This solution is not unique. If the gas is in hydrostatic equilibrium for example, the gas vertical extension and the velocity dispersion are not independent variables. An approximation to the scale height of the gas for a galactic disk of stars and gas in hydrostatic equilibrium in a dark matter halo is given by

$$h_g=\frac{c_g}{\pi G}\left(\frac{\mathrm{\Sigma }_g+\mathrm{\Sigma }_{dm}}{c_g}+\frac{\mathrm{\Sigma }_s}{c_s}\right)^1$$
(11)
Since $`\mathrm{\Sigma }_{dm}`$ is the dark matter surface density between $`+h_g`$ and $`h_g`$ the above formula should be solved numerically for determining $`h_g`$ and $`\mathrm{\Sigma }_{dm}`$. In Figure 8$`(d)`$ we plot $`f(\mathrm{\Sigma }_g+\mathrm{\Sigma }_{dm}+\mathrm{\Sigma }_{})`$, for $`\mathrm{\Sigma }_{dm}`$ between $`\pm h_g`$ for the same dark matter density used in Figure 8$`(c)`$. We find that $`h_g`$ varies between 20 pc and 1 kpc through the radial extent of the disk analyzed in this paper. This model reduces the radii of the unstable star forming regions respect to the constant scale height model displayed in Figure 8$`(c)`$, making both $`\mathrm{\Sigma }_A`$ and $`\mathrm{\Sigma }_T`$ suitable indicators of the growth of instabilities in the disk. Which one it is better to use for defining a threshold radius in M33 is therefore still unclear since this depends on the radial variations of the involved parameters such as sound speed, thickness etc. which are not well determined yet. It is important to notice however that since the densities of visible baryons and dark matter are well constrained at intermediate radii, the intersections between $`\mathrm{\Sigma }_A`$ or $`\mathrm{\Sigma }_T`$ and $`f(\mathrm{\Sigma }_g+\mathrm{\Sigma }_{dm}+\mathrm{\Sigma }_{})`$ do not depend on which mass model we use, as soon as it is compatible with the rotational velocity data. The Burkert best fitting mass model, for example, gives the same solutions displayed in Figure 8 for a NFW dark matter halo.

## 6 Summary

A full map of the nearby galaxy M33 in the CO J=1-0 transition has given additional data for improving the accuracy of the rotation curve of this galaxy and the knowledge of the visible baryons distribution. The molecular gas peaks in the central regions with a surface density of about 7 M pc-2 and declines exponentially with a radial scalelength of about 2.5 kpc down to about 0.6 M pc-2 at $`R6`$ kpc. Our sensitivity was not sufficient to outline the molecular gas radial distribution further out and in particular to detect any possible sudden drop of the azimuthally averaged molecular gas surface density. The molecular gas is not smoothly distributed in the inner disk since there are remarkable mass concentrations along the spiral arms. Some wiggles of about 10 km s-1 are visible in the rotation curves of the two separate halves of the galaxy, which are reminiscent of density wave perturbations. Details of atomic and ionized gas distribution in M33 can be found in Hoopes $`\&`$ Walterbos (2000) and in Thilker, Braun $`\&`$ Walterbos (2001). The rotation curve derived using H$`\alpha `$ data for the innermost point, the CO data and the HI 21-cm line data, has been used to constrain the visible baryonic mass and dark matter halo models. The total gas mass (HI+H2+He) is $`3.2\times 10^9`$ M. This is of the same order of the stellar disk mass, estimated to be between 3$`\times 10^9`$ and 6$`\times 10^9`$ M. The total visible baryonic mass of M33 can therefore be as high as $`10^{10}`$ M. The dark halo mass out to the last measured point ($`17`$ kpc from the center) is $`5\times 10^{10}`$ M.

Both the non-singular isothermal halo and the Burkert halo model, give good fits to the M33 rotation curve when the flat density core is extended. The pure $`R^2`$ or $`R^3`$ outer scaling law applies only beyond the outermost observed radius. Dark matter density profiles with an inner $`R^{1.5}`$ cusp (Moore et al. 1998) are too steep to be consistent with the rotation of the central regions of M33 and are ruled out. Halo models with milder central cusps, such as the NFW radial profile, fit well all the observed rotational velocities. The mass of the virialized NFW dark halo for this galaxy is $`5\times 10^{11}`$ M and the halo size comparable to the distance between M33 and its bright companion M31 ($`180`$ kpc). The baryon fraction in M33 is then $`0.02`$. The required concentrations $`C_\mathrm{\Delta }`$ are below the scaling laws found by numerical simulations of structures formation in a standard Cold Dark Matter scenario ($`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Lambda }_0=0.7`$, $`\sigma _8`$=0.9) and favor cosmological models which predicts less concentrated dark halos (Eke, Navarro $`\&`$ Steinmetz 2001, Zentner $`\&`$ Bullock 2002, van den Bosch, Mo $`\&`$ Yang 2003). We have shown however that if an associated scatter in the numerical simulation results is considered (Wechsler et al. 2002), the M33 dark halo is still compatible with galaxy formation models in a standard $`\mathrm{\Lambda }`$CDM cosmology. The finding that dark matter halos are consistent with constant density core models or with central cusps less steep than $`R^{1.5}`$, is known as the cusp-core degeneracy (van den Bosch $`\&`$ Swaters 2001). It is remarkable that this degeneracy holds also in M33 which is a galaxy with small uncertainties about its distance, its gaseous and stellar content, and with circular velocities sampled over a wide range of radii: from galactocentric distances of a few parsecs out to 13 disk scalelengths. For this galaxy it seems that the degeneracy holds mostly because of a nuclear stellar component, which regulates the dynamics of the innermost 0.5 kpc region. This extended ‘stellar nucleus’, with mass between $`3\times 10^7`$ and $`8\times 10^8`$ M, is not only dynamically required by the rotational data but it is also evident from photometric data (Bothun 1992, Minniti, Olszewski $`\&`$ Rieke 1993, Regan $`\&`$ Vogel 1994). We have considered a core collapse model and a de Vaucouleurs $`R^{1/4}`$ law for the nuclear stellar density but additional data on the kinematics and light distribution of the central few hundreds parsec region are needed to constrain better the mass distribution there.

Having a very accurate rotation curve and a good estimate of the gas, stellar and dark matter surface density in the M33 disk, we have addressed the question on the possible coincidence of the active star formation region in this galaxy with a gravitationally unstable region. Azimuthal averages of the observed gas surface densities are above the critical densities in the star forming region of M33 only if one considers the minimum acceptable value of the gas velocity dispersion and a stability criterion based on the shear rate of the disk. However, the simple Toomre condition predicts a star formation threshold radius in agreement with the observed drop in the H$`\alpha `$ surface brightness if the additional compression due to the stellar disk or to the dark matter surface density is considered.

I acknowledge Steve Schneider and Mark Heyer for having given me the opportunity to use the CO data prior to publication. I am also very thankful to the anonymous referee who made several excellent suggestions for improving the quality of the work presented in the original manuscript.


format: text

text:
# Radio Observations of Two Isolated Neutron Stars, RX J0720.4–3125 and RX J0806.4–4132

## 1 Introduction

A new class of neutron stars has emerged with the discovery of the so-called Isolated Neutron Stars (INS) in the X-ray band with ROSAT. Their identification as neutron stars is secure as they have extremely high X-ray to optical flux ratios. They are isolated in the sense that they have no binary companions which rules out accretion from the companion as a mechanism for generating the X-rays. They are all extremely close, probably only a few hundred parsecs (at most) from Earth. There is no known radio emission (pulsed or unpulsed) associated with these objects although deep searches have not yet been reported in the literature. A review of these objects can be found in Treves et al. (2000).

There are 3 possibilities as to the nature of these neutron stars. First, they might be old stars accreting from the interstellar medium. These still rotate, albeit slowly, but do not produce ‘conventional’ radio or X-ray emission. The second possibility is that they are ‘standard’ (old) neutron stars, with magnetic fields of order 1012 G which pulsate in X-rays and pulsate in the radio but may or may not be beamed in our direction. Finally, they may be ‘magnetars’, which are young neutron stars with extremely high magnetic fields (greater than 1014 G). These have X-ray pulsations but no radio pulsations. The mechanism for producing X-ray pulsations is completely different to conventional pulsars; they are powered by magnetic rather than rotational energy.

When the INS were first discovered a few years ago, it seemed unlikely that they were standard radio pulsars \[Heyl & Kulkarni 1998\]. No radio pulsars at that time had periods greater than 4 seconds and magnetic fields more than 1013 G. However, in 1999, PSR J2144$``$3933 was discovered to have a period of 8.5 s (Young, Manchester & Johnston 1999) and Camilo et al. (2000) discovered PSR J1814–1744 with a long period and a high magnetic field similar (in parameter space at least) to the so-called ‘magnetars’. PSR J2144$``$3933 has a very low dispersion measure, placing it about 200 pc distant, a similar distance to the known INS. It therefore now seems more likely that the INS are normal pulsars rather than some new exotic species.

One way of determining the nature of these objects is detection of periods and period derivatives. The combination of period and period derivative naturally yields the magnetic field strength, which is the key discriminator in the scenarios outlined above. Four of the 7 known INS have X-ray pulsations with rotation periods ranging from 5.2 to 22 seconds. Of interest here are RX J0720.4–3125 and RX J0806.4–4132 which have spin periods of 8.4 and 11.4 seconds respectively \[Haberl & Zavlin 2002, Kaplan et al. 2002, Zane et al. 2002\]. Timing analysis has so far not yielded a period derivative (and hence magnetic field strength) for these pulsars. Zane et al. (2002) claim a period derivative for RX J0720.4–3125 but their phase connected solution has been called into question by Kaplan et al. (2002). However, the upper limit on the period derivative seems to rule out the magnetar model. For RX J0806.4–4132 the upper limit on the period derivative does not allow discrimination between the models.

The question then is - are there radio pulses and/or continuum radio emission (e.g. from a wind nebula) from these objects? If radio pulsations were detected then it becomes relatively easy to measure a period derivative and hence determine the nature of these objects. If detected, they would be among the closest radio pulsars, and likely the lowest luminosity pulsars known. This would have implications for the birth rate of pulsars and the low-luminosity cut-off. It is a non-trivial task, however, to detect these long period pulsars using conventional search techniques as the red-noise in the Fourier transform of the time series seriously reduces the sensitivity. Knowledge of the period a prori helps the search process.

## 2 Observations and Data Reduction

Observations of RX J0720.4–3125 and RX J0806.4–4132 were made at the Australia Telescope Compact Array (ATCA) on 2002 December 4 and 5 using the 6A array configuration. In this configuration the shortest baseline is 627 m and the longest 5940 m. In each observation, two continuum frequencies were observed simultaneously with 32 4-MHz channels across each; all four Stokes parameters were recorded. The centre frequencies were 1384 and 1704 MHz. The pointing centres used were RA (J2000) 07h 20m 24.96, Dec (J2000) –31° 25′ 19.6″ and RA (J2000) 08h 06m 23.47, Dec (J2000) –41° 22′ 02.3″ respectively, both $``$30″ offset from the X-ray location of the INS. The total integration time was 10 hr per source. The ATCA is also capable of splitting each correlator cycle into bins corresponding to different phases of a pulsar’s period, and, in this case, the pulse period divided into 32 time bins. The flux density scale of the observations was determined by observations of PKS B1934–638, while polarization and antenna gains and phases were calibrated using 3 min observations every 50 min of PKS B0614$``$349 for RX J0720.4–3125 and PKS B0823$``$500 for RX J0806.4–4132. As a system test, observations were made of the known 8.5 s radio pulsar PSR J2144$``$3933 using an identical setup to that described above. The source was observed for 30 min.

Data were edited and calibrated using the MIRIAD package. At each frequency the field of interest was imaged using conventional techniques, after first collapsing the phase bins.

The data were also folded to produce a pulse profile. The dispersion smearing was assumed to be negligable compared to the duration of one phase bin. The data were folded according to the ephemeris given in Young et al. (1999) for PSR J2144$``$3933, assuming a barycentric period of 8.391115 s for RX J0720.4–3125 as given by Kaplan et al. (2002), and assuming a barycentric period of 11.3714 s for RX J0806.4–4132 as given in Haberl & Zavlin (2002). Note that in the case of RX J0720.4–3125, if the star has a period derivative in excess of $`3.8\times 10^{13}`$, this would be sufficient to broaden the pulse profile by more than one bin over the course of the 10 hr observation.

## 3 Results and Discussion

### 3.1 PSR J2144$``$3933

The pulse profile of PSR J2144$``$3933 at 1.4 GHz is shown in Fig. 1. The pulsar is known to have a short duty cycle and hence its pulse fits entirely within one phase bin in this instance. The polarization of the pulsar is low, less than 10% of the total intensity. This is consistent with observations at lower frequencies \[Manchester, Han & Qiao 1998\]. The flux density at 1.4 GHz is $``$0.8 mJy; the flux density at 436 MHz is 4 mJy (Young et al. 1999) which yields a (fairly typical) spectral index of –1.4. This observation shows the viability of detecting these long period pulsars using an interferometer.

### 3.2 RX J0720.4–3125

Continuum images of the region around RX J0720.4–3125 are shown in Figs. 2 and 3. The position of RX J0720.4–3125 as given by the Chandra observations of Kaplan et al. (2002) is shown by the + symbol. The extended source seen in Fig. 3 does not appear to have an optical counterpart in the Deep Sky Survey nor in the deeper images of Kulkarni & van Kerwijk (1998). It is likely that it is a background extragalactic source. No radio emission is detected from RX J0720.4–3125 with a 3-$`\sigma `$ upper limit of 0.2 mJy. A more stringent limit can be placed on the flux density if we assume that the pulse is contained within one phase bin (as is the case for PSR J2144$``$3933 above). The upper limit is then reduced by a factor $`\sqrt{32}`$ to 0.05 mJy.

If RX J0720.4–3125 were a conventional pulsar (i.e. not a magnetar), limits on the efficiency of the coupling of the relativistic wind with the interstellar medium can be determined. Following Gaensler et al. (2000), if a putative wind nebulae was a point source, the efficiency parameter, $`ϵ`$ is given by

$$ϵ=\frac{d^2S}{2.1\times 10^2\dot{E}_{31}}$$
(1)
where $`d`$ is the distance to the source in kpc, $`S`$ is the flux density upper limit in mJy and $`\dot{E}`$ is the spin-down energy in units of $`10^{31}`$ erg s-1. For RX J0720.4–3125, therefore, with $`d=0.2`$ kpc and $`\dot{E}<2\times 10^{31}`$ erg s-1, $`ϵ`$ would be $`<5\times 10^5`$, a value not atypical of other middle-aged pulsars. Therefore we cannot conclude from the lack of a radio nebula that this is not a neutron star powered by spin-down.

Unfortunately, the lack of radio emission from RX J0720.4–3125 does little to help our understanding of the nature of this source. The beaming fraction of these long period pulsars is expected to be less than 10% on theoretical grounds \[Tauris & Manchester 1998\] and the radio luminosity is low, as seen in PSR J2144$``$3933. The upper limit to pulsed emission from this source is a factor 10 better than for PSR J2144$``$3933 and, assuming a similar distance of 180 pc makes it unlikely that it is beaming in our direction but is below the sensitivity limit. It therefore remains likely that this is a conventional pulsar in which the radio beams do not intersect Earth \[Brazier & Johnston 1999\]. The question then remains as to whether it is an old pulsar, with an age similar to that of PSR J2144$``$3933 or a much younger pulsar. The former implies that the age derived from the blackbody fits to the X-ray spectrum was incorrect whereas the latter implies that it must have been born with a very long period, unlike the bulk of the known pulsars.

### 3.3 RX J0806.4–4132

A continuum image of a 2′ region around RX J0806.4–4132 is shown in Fig. 4. The position of RX J0806.4–4132 as given by the XMM-Newton observations of Haberl & Zavlin (2002) is shown. The area covered by the + symbol is 3 times larger than the error box given by Haberl & Zavlin (2002). Only 1 source is seen in the field. It has a flux density of 0.4 mJy and is offset from the position of RX J0806.4–4132 by 34″. A search for pulsations at the known period failed to reveal any significant detection. Assuming the pulsed flux is contained within one phase bin implies a lower limit on the pulsed flux of 0.05 mJy.

The pulsed fraction in X-rays is very low for RX J0806.4–4132, only 6%. Haberl & Zavlin (2002) speculate that perhaps the neutron star has its magnetic axis and rotational axis nearly aligned. If this is the case, then radio emission might also be expected. The lack of a detection in the radio either implies that this proposition is not correct, that the radio pulsar is an order of magnitude less luminous than PSR J2144$``$3933, already one of the lowest luminosity pulsars known or perhaps that it has crossed the death-line and is incapable of producing radio emission.

### 3.4 Speculation on the nature of the INS

Table 1 lists the 7 known radio pulsars and the X-ray pulsar Geminga within 400 pc of Earth (excluding the millisecond pulsars) ranked in descending order of $`\dot{E}`$ along with the two INS observed here. Of these, the four with the highest $`\dot{E}`$ have been detected in the X-ray band, essentially due to the sensitivity of the ROSAT all-sky survey. Nearby, spin-powered pulsars, must have $`\dot{E}`$ greater than $``$5$`\times 10^{31}`$ erg s-1 in order to be detected. Only PSR B1133+16 of the non-detected pulsars has an $`\dot{E}`$ slightly in excess of this value. PSR J2144$``$3933 has an $`\dot{E}`$ which is three orders of magnitude less than this and therefore not only is the lack of X-ray emission from this pulsar not surprising, but it also seems unlikely that the INS share the same characteristics.

It therefore seems most likely that the INS will turn out to have $`\dot{E}`$ in excess of 5$`\times 10^{31}`$ erg s-1, consistent with the inferred $`\dot{E}`$ of RX J1856.5–3754 \[van Kerkwijk & Kulkarni 2001, Kaplan, van Kerkwijk & Anderson 2002\] and marginally consistent with the upper limit derived by Kaplan et al. (2002) for RX J0720.4–3125. For RX J0806.4–4132, the combination of an $`\dot{E}`$ of 5$`\times 10^{31}`$ erg s-1 and the spin period of 11.4 s would imply a high value of the period derivative ($`2\times 10^{12}`$) and a magnetic field in the magnetar range ($`2\times 10^{14}`$ G) and even more extreme parameters for the 22.7 s INS RX J0420.0–5022. These could then be magnetars where the energy source derives from the magnetic field rather than the spin-down. In summary, the INS appear to be one of the manifestations of neutron stars born with high magnetic fields. Their number suggests they must form a substantial fraction of the population of neutron stars in the Galaxy but that their radio emission may be inhibited.

## Acknowledgements

The Australia Telescope is funded by the Commonwealth of Australia for operation as a National Facility managed by the CSIRO. I thank Warwick Wilson for his efforts beyond the call of duty in programming the correlator to ensure maximum efficiency for these long period pulsars.


text:
# Reply to “Comment on ‘Tunnel Window’s Imprint on Dipolar Field Distributions’ ”

Stamp and Tupitsyn (ST) have recently posted a Comment ST  on a paper of ours ourPRL . They claim ST  that, in Ref. ourPRL , “Short-time relaxation was caused by a simple noise field, acting uniformly over an ‘energy window’ of width $`2\delta h_{hf}`$ (to simulate the fluctuating nuclear spin bias)”. We see no basis for this statement.

Whereas previous work of Prokof’ev and Stamp PSPRL ; PS1  deals with nuclear spin bath mechanisms that give rise to a tunnel energy window (TEW) for single-molecule magnets, such as Fe8, the work we report in Ref. ourPRL  is independent of such mechanisms. The only definition of the TEW, $`\delta h_{hf}`$, we make use of in Ref. ourPRL  is the one we give therein in the abstract: “we allow spin flips only if the corresponding energy change is less than some $`2\delta h_{hf}`$”. This is further explained below Eq. (2) in Ref. ourPRL , where the system under consideration is defined as an Ising model of spins (that clearly stand for the electronic cluster spins) interacting through dipolar fields with the flipping rule just mentioned. Thus, the definition of $`\delta h_{hf}`$ stands on its own without need for reference to the nuclear bath.

In addition, we mention hyperfine interactions only in passing, in an introductory remark, but the phrase “typical hyperfine coupling”, which PT ascribe to us, does not appear in Ref. ourPRL . We do refer once to $`\delta h_{hf}`$ as “a typical hyperfine energy”, but make no use whatsoever of this non-definition.

We can similarly see no basis for the statement towards the end of the Comment by ST ST : “… but here using the kinetic equation of Ref. , rather than just a phenomenological noise.”, since the Monte Carlo simulations we perform in Ref. ourPRL  happen to be exp  a stochastic implementation of “the kinetic equation”, i.e., Eq. (4) of Ref.PSPRL .

Finally, at the risk of being repetitious, we cannot understand the closing remark in the preceeding Comment, “…so we predict a steep fall of $`\xi _0`$, which would be very hard to understand in the approach of ref. ”. Indeed, there is nothing to understand about this in the approach of Ref. ourPRL , because the physics of nuclear baths which is behind TEW’s is not at all within the scope of Ref. ourPRL .

What we do show in Ref.ourPRL  is that a hole of half-width $`W`$ develops in the dipole field distribution of an Ising model of spins interacting through dipolar fields with the flipping rule just mentioned, and that $`W0.75\delta h_{hf}`$, and we know of no previously published work establishing a relation between these two quantities.

Juan J. Alonso1 and Julio F. Fernández2
1Física Aplicada I, Universidad de Málaga, 29071-Málaga, Spain. 2ICMA, CSIC and Universidad de Zaragoza, 50009-Zaragoza, Spain


text:
# Semigroup techniques for the efficient classical simulation of optical quantum information

## 1 Introduction

Information processing using the rules of quantum mechanics may allow tasks that cannot be performed using classical laws . The efficient factorization algorithm of Shor  and secure quantum cryptography  are two examples. Of the many possible realizations of quantum information processes, optical realizations have the advantange of negligible decoherence: light does not interact with itself, and thus a quantum state of light can be protected from becoming entangled with the environment. Several proposed optical schemes  offer significant potential for quantum information processing.

In order to prove theorems regarding the possibilities and limitations of optical quantum computation, one must construct a framework for describing all types of physical processes (unitary transformations, projective measurements, interaction with a reservoir, etc.) that can be used by an experimentalist to perform quantum information processing. Most frameworks currently employed (e.g., ) are restricted to describing only unitary transformations. However, such transformations are a subset of all possible physical processes. Non-unitary transformations such as dissipation, noise, and measurement must also be described within a complete framework. The new results of Knill et al  show that photon counting measurements allow for operations that are “difficult” with unitary transformations alone; thus, non-unitary processes may be a powerful resource in quantum information processing and must be considered in any framework that attempts to address the capabilities of quantum computation with optics.

In this paper, we show that unitary transformations, measurements and any other physical process can be described in the unified formalism of completely positive (CP) maps. Also, a broad class of these maps which includes linear optics and squeezing transformations, noise processes, amplifiers, and measurements with feedforward that are typical to quantum optics experiments can be described within the framework of a *Gaussian semigroup*. This framework allows us to place limitations on the potential power of certain quantum information processing tasks.

One important goal is to identify classes of processes that can be efficiently simulated on a classical computer; such processes cannot possibly be used to provide any form of “quantum speedup”. The Gottesman-Knill (GK) theorem  for qubits and the CV classical simulatability theorems of Bartlett et al  provide valuable tools for assessing the classical complexity of a quantum optical process. It is shown here that semigroup techniques provide a powerful formalism with which one can address issues of classical simulatability. In particular, a classical simulatability result is presented for a general class of quantum optical operations, and thus a no-go theorem for quantum computation with optics is proven using semigroup techniques.

## 2 Semigroup Description of Gaussian operations

Consider an optical quantum information process involving $`n`$ coupled electromagnetic field modes, with each mode described as a quantum harmonic oscillator. The two observables for the (complex) amplitudes of a single field mode serve as canonical operators for this oscillator. A system of $`n`$ coupled oscillators, then, carries an irreducible representation of the Heisenberg-Weyl algebra hw($`n`$), spanned by the $`2n`$ canonical operators $`\{q_i,p_i,i=1,\mathrm{},n\}`$ along with the identity operator $`I`$. These operators satisfy the commutation relations $`[q_i,p_j]=i\mathrm{}\delta _{ij}I`$. We express the $`2n`$ canonical operators in the form of a phase space vector $`\mathrm{z}`$ with components $`z_i=q_i`$ and $`z_{n+i}=p_i`$ for $`i=1,\mathrm{},n`$. These operators satisfy $`[z_i,z_j]=i\mathrm{}\mathrm{\Sigma }_{ij}`$, with $`\mathrm{\Sigma }`$ the skew-symmetric $`2n\times 2n`$ matrix

$$\mathrm{\Sigma }=\left(\begin{array}{cc}0& I_n\\ I_n& 0\end{array}\right)$$
(1)
and $`I_n`$ the $`n\times n`$ identity matrix. For a state $`\rho `$ represented as a density matrix, the *means* of the canonical operators is a vector defined as the expectation values $`𝝃=\mathrm{z}\mathbf{}_𝝆`$, and the *covariance matrix* is defined as

$$\mathrm{\Gamma }=(\mathrm{z}\mathbf{}𝝃\mathbf{)}\mathbf{(}\mathrm{z}\mathbf{}𝝃\mathbf{)}^{\mathbf{}}\mathbf{}_𝝆\mathbf{}\mathrm{i}\mathrm{\Sigma }\mathbf{.}$$
(2)
A *Gaussian state* (a state whose Wigner function is Gaussian and thus possesses a quasiclassical description) is completely characterized by its means and covariance matrix . Coherent states, squeezed states, and position- and momentum-eigenstates are all examples of Gaussian states.

We define $`𝒞_n`$ to be the group of linear transformations of the canonical operators $`\{z_i\}`$ ; this group corresponds to the infinite-dimensional (oscillator) representation of the “Clifford group” employed by Gottesman . For a system of $`n`$ oscillators, it is the unitary representation of the group ISp($`2n,`$) (the inhomogeneous linear symplectic group in $`2n`$ phase space coordinates)  which is the semi-direct product of phase-space translations (the Heisenberg-Weyl group HW($`n`$)) plus one- and two-mode squeezing (the linear symplectic group Sp($`2n,`$)). Phase space displacements are generated by Hamiltonians that are linear in the canonical operators; a displacement operator $`X(𝜶)`$ HW($`n`$) is defined by a real $`2n`$-vector $`𝜶`$. A symplectic transformation $`M(A)`$ Sp($`2n,`$), with $`A`$ a real matrix satisfying $`A^{}\mathrm{\Sigma }A=\mathrm{\Sigma }`$, is generated by a Hamiltonian that is a homogeneous quadratic polynomial in the canonical operators. A general element $`C𝒞_n`$ can be expressed as a product $`C(𝜶,A)=X(𝜶)M(A)`$, and transforms the canonical operators as

$$C(𝜶,A):\mathrm{z}\mathbf{}\mathrm{z}^{\mathbf{}}\mathbf{=}\mathrm{z}\mathrm{A}\mathbf{+}𝜶\mathbf{.}$$
(3)
The group $`𝒞_n`$ consists of unitary transformations that map Gaussian states to Gaussian states; however, unitary transformations do not describe all physical processes. In the following, we include other (non-unitary) CP maps that correspond to processes such as dissipation or measurement. We define the *Gaussian semigroup*, denoted $`𝒦_n`$, to be the set of Gaussian CP maps  on $`n`$ modes: a Gaussian CP map takes any Gaussian state to a Gaussian state. Because Gaussian CP maps are closed under composition but are not necessarily invertible, they form a semigroup. A general element $`T𝒦_n`$ is defined by its action on the canonical operators as

$$T(𝜶,A,G):\mathrm{z}\mathbf{}\mathrm{z}^{\mathbf{}}\mathbf{=}\mathrm{z}\mathrm{A}\mathbf{+}𝜶\mathbf{+}𝜼\mathbf{,}$$
(4)
where $`𝜶`$ is a real $`2n`$-vector, $`A`$ and $`G`$ are $`2n\times 2n`$ real matrices, and $`A`$ is no longer required to be symplectic. Eq. (4) includes the transformations (3) plus additive noise processes  described by quantum stochastic noise operators (the vector $`𝜼`$) with expectation values equal to zero and covariance matrix

$$𝜼𝜼^{}_{\rho _R}i\mathrm{\Sigma }=GiA^{}\mathrm{\Sigma }A.$$
(5)
Here, $`\rho _R`$ is a Gaussian ‘reservoir’ state which, in order to define a CP map, must be chosen such that the noise operators satisfy the quantum uncertainty relations. This condition is satisfied if the noise operators define a positive definite density matrix, which leads to the condition

$$G+i\mathrm{\Sigma }iA^{}\mathrm{\Sigma }A0.$$
(6)
The group $`𝒞_n`$ is recovered for $`G=0`$.

The action of the Gaussian semigroup on the means and covariance matrix is straightforward and given by

$$T(\alpha ,A,G):\{\begin{array}{cc}𝝃𝝃^{}=𝝃A+𝜶\hfill & \\ \mathrm{\Gamma }\mathrm{\Gamma }^{}=A^{}\mathrm{\Gamma }A+G.\hfill & \end{array}$$
(7)
Because the means and covariance matrix completely define a Gaussian state, the resulting action of the Gaussian semigroup on Gaussian states can be easily calculated via this action.

The Gaussian semigroup $`𝒦_n`$ represents a broad framework to describe several important types of processes in a quantum optical circuit. The group $`𝒞_n𝒦_n`$ comprises the unitary transformations describing phase–space displacements and squeezing (both one– and two–mode). Introduction of noise to the circuit (e.g., via linear amplification) is also in $`𝒦_n`$. Furthermore, the Gaussian semigroup describes certain measurements in the quantum circuit. These include measurements where the outcome is discarded (thus evolving the system to a mixed state) or retained (where the system follows a specific quantum trajectory defined by the measurement record ). Finally, the Gaussian semigroup includes Gaussian CP maps conditioned on the outcome of such measurements. For details and examples of all of these types of Gaussian semigroup transformations, see .

## 3 Classical Simulation of Gaussian Semigroup Processes

Using the framework of the Gaussian semigroup, it is straightforward to prove the classical simulatability result of Bartlett and Sanders .

Theorem: Any quantum information process that initiates in a Gaussian state and that performs only Gaussian semigroup maps can be *efficiently* simulated using a classical computer.

Proof: Recall that any Gaussian state is completely characterized by its means and covariance matrix. For any quantum information process that initiates in a Gaussian state and involves only Gaussian semigroup maps, one can follow the evolution of the means and the covariance matrix rather than the quantum state itself. For a system of $`n`$ coupled oscillators, there are $`2n`$ independent means and $`2n^2+n`$ elements in the (symmetric) covariance matrix; thus, following the evolution of these values requires resources that are polynomial in the number of coupled systems. QED

Because most current experimental techniques in quantum optics are describable by Gaussian semigroup maps, this theorem places a powerful constraint on the capability of achieving quantum computational speedups (tasks that are not efficient on any classical machine) using quantum optics.

## 4 Conclusions

Semigroup techniques provide a powerful tool for constructing and assessing new quantum information protocols using quantum optics. These techniques have been used to show that algorithms or circuits consisting of only Gaussian semigroup maps can be efficiently simulated on a classical computer, and thus do not provide the ability to perform quantum information processing tasks efficiently that cannot be performed efficiently on a classical machine. Eisert et al  use related techniques to show that local Gaussian semigroup transformations are insufficient for distilling entanglement: an important process for quantum communication and distributed quantum computing. Most current quantum optics experiments consist only of Gaussian semigroup transformations; thus, the challenge is to exploit this semigroup to prove new theorems, limitations and possibilities for quantum information processing using optics.

This project has been supported by Macquarie University and the Australian Research Council. The author thanks B. C. Sanders for helpful discussions.

## References