text:
# Birth rates of different types of neutron star and possible evolutions of these objects

## 1 Introduction

Existence of single neutron stars with different physical properties is very well known. But the amount of data about different types of neutron star except radio pulsars is small. Actually, the change in the beaming factor together with the value of and the change in the angle between the magnetic field and the rotation axes are not well known. It is also not clear how the relation between the characteristic time ($`\tau `$) and the real age changes in time. Therefore, the birth rate, initial periods and values of the real magnetic field are not known well. So, it is difficult to understand what the evolutionary tracks of pulsars on the P-Ṗ diagram must be and where these tracks end. As there are not many available observational data for the other types of neutron star and since they show some exotic phenomena such as $`\gamma `$-ray bursts, it is difficult to understand their nature. Below, we analyse birth rates of different types of neutron star, their locations on the P-Ṗ diagram and possible evolutions of these objects.

## 2 Analysis of the data about the birth rates of DITNSs and supernova rate around the Sun

According to Yakovlev et al. (2002) and Kaminker et al. (2002) neutron stars cool down to T=4.6$`\times `$105 K (or kT=40 eV) in t$``$106 yr. According to Haberl (2003, and references therein), there are 7 X-ray dim isolated thermal neutron stars (DITNSs) located within 120 pc around the Sun and only one DITNS, namely RX J1836.2+5925, is located at a distance of 400 pc. All of these 8 radio-quiet objects have T$`>`$4.6$`\times `$105 K. On the other hand, there are only 5 radio pulsars with characteristic times smaller than 107 yr in the cylindrical volume around the Sun with a radius of 400 pc and height 2$``$z$``$=400 pc, where $``$z$``$ is the distance from the Galactic plane. Among these 5 pulsars, Vela has $`\tau `$$``$104 yr and the other 4 pulsars have 3$`\times `$106$`<`$$`\tau `$$`<`$5$`\times `$106 yr. The real ages of 1-2 of these 4 radio pulsars may be smaller than 106 yr. We have included Geminga pulsar (B0633+1748) in Table 1 to compare it with different types of neutron star, since this pulsar has properties in between the properties of DITNSs and of young radio pulsars.

In order to make a reliable estimation of the number density of the sources in each considered volume, it is necessary to consider how these objects move in the Galaxy. Practically, for pulsars with such ages the space velocity must not decrease in the gravitational field of the Galaxy. In order to demonstrate this in a simple and reliable way, we can use the observational data about the scale height and the average peculiar velocity of old stars near the Sun.

The average values of $``$z$``$ and $``$Vz$``$ for population II stars which belong to the old disk, the intermediate and the halo are, respectively, 400 pc and 16 km/s, 700 pc and 25 km/s, and 2000 pc and 75 km/s (Allen 1991). Ages of these stars are close to 1010 years. They are in dynamical equilibrium and their ages are large enough to make many oscillations in the gravitational field of the Galaxy.

Let us consider, in the first approximation, their oscillations in a homogeneous field in the direction perpendicular to the Galactic plane. In this case

$$V_z=A\sqrt{z}cm/sec$$
(1)
and the period of oscillation is

$$T=\frac{8V_z}{A^2}sec$$
(2)
where A is a constant which is related to the gravitational field intensity near the Sun. If we use the average values of $``$z$``$ and $``$Vz$``$ for different subgroups of population II stars as written above, then we find values of A very close to each other showing the reliability of the approximation. Therefore, we can adopt a value of A=5.2$`\times `$10-5 cm1/2s-1, which is the average value for the 3 subclasses of population II objects, for all objects located up to about 2000 pc away from the Galactic plane and we can also use this value for pulsars. If we put the average value of $``$Vz$``$=170-200 km/s for pulsars in eqn.(1), then we find T=(1.2-1.7)$`\times `$109 yr. As we see the age values of young radio pulsars are very small compared to the average period of pulsar oscillation. So, we can adopt that pulsars with ages 106-107 yr move practically with constant velocity which they gain at birth.

Since the progenitors of neutron stars have scale height of about 60 pc, most of the neutron stars are born close to the Galactic plane. Therefore, they can go 200 pc away from the plane in about 106 yr. As the Sun is located very close to the Galactic plane, considering the value of z of the Sun in the analysis does not change the results and so can be neglected.

There are 23 supernova remnants (SNRs) within about 3.2 kpc around the Sun with surface brightness $`>`$10-21 Wm-2Hz-1sr-1 and the ages of these SNRs are not greater than (3-4)$`\times `$104 yr (Green 2001; Ankay et al. 2003). So, the number of neutron stars which were born in the last 106 yr within the region of 400 pc around the Sun must be about 11. Only 2 or 3 of these 11 neutron stars can go more than 200 pc away from the Galactic plane. Therefore, the considered region may contain about 8-9 neutron stars. But the birth rate of radio pulsars is up to 3-4 times smaller than the supernova explosion rate in our galaxy (Ankay et al. 2003). So, we can expect only 3 radio pulsars with an age up to 106 yr located in this region and this is in accordance with the data mentioned above.

According to Haberl (2003) (see also the references therein), surface temperatures of DITNSs are 40-95 eV and their ages must not be more than 106 yr (Kaminker et al. 2002; Yakovlev et al. 2002) as they are in the cooling stage. Guseinov et al. (2003a) include most of these sources (for which there exist considerably more information) in the list of isolated radio quiet neutron stars which were observed in X-ray band. We have included in Table 1 only a few very important data taken from the tables given in Guseinov et al. (2003a) together with some new data. In some cases, the distance values of DITNSs are about 2-3 times larger than the ones given in Haberl (2003) and in some other papers (see also the references in Guseinov et al. 2003a and Haberl 2003). Is it acceptable to adopt such smaller values of distance as given by Haberl (2003) for these objects?

The statistics about the SNRs and the radio pulsars in the considered region are very poor, but the number density of DITNSs for the case of small distances (Haberl 2003) is very large so that there exists a contradiction. If the data about the distances and ages of DITNSs were reliable, then their birth rate would turn out to be about 4 times more than the supernova explosion rate.

Since the measured space velocities of all types of neutron star are very large compared to the space velocities of O and B-type stars, there is no doubt about the origin of neutron star formation which is due to the collapse of the progenitor star together with supernova explosion. For example, RX J0720.4-3125 at a distance of 200 pc (which may actually be $``$300 pc) has a tangential velocity of about 100 km/s calculated from the proper motion $`\mu `$=97$`\pm `$12 mas/yr (Motch et al. 2003). PSR J0538+2817 has proper motion $`\mu `$=67 mas/yr and has a transverse velocity in the interval 255-645 km/s at a distance of 1.2 kpc (Kramer et al. 2003). Therefore, the birth rate of neutron stars can not exceed the rate of supernova explosion. How can we explain the contradiction between the neutron star birth rate and the supernova explosion rate? Note that Haberl (2003) does not assume all the DITNSs to be only cooling neutron stars but also discusses other possibilities, for example the accretion from interstellar gas. But as a rule, today the cooling origin of the X-ray radiation is considered and this is reliable (Haberl 2003). In order to solve the problem about the difference in the birth rates it is necessary to adopt either larger distance values or longer lifetimes for DITNSs.

## 3 The distances of DITNSs

Since the theory of cooling of neutron stars has been well developed and all the DITNSs have T$`>`$4.6$`\times `$105 K, we can say that these neutron stars have ages less than 106 yr (Pavlov et al. 2002a,b; Kaplan et al. 2003a,b; Pavlov & Zavlin 2003; Yakovlev et al. 2003). As the ages of these objects are small, their number density must be small. Therefore, it is better to adopt larger values of distance for these neutron stars.

Before beginning to discuss how to adopt distance values, it is necessary to analyse the data of DITNSs. These data vary a lot from one observation to another. The most recent data about DITNSs are given in Haberl (2003) and Guseinov et al. (2003a).

For all DITNSs there exist data about their temperatures and it is known that they have approximately blackbody radiation (Pavlov et al. 2002b; Haberl 2003). In this work, we have used the temperature values given in Pavlov et al. (2002b) and Haberl (2003) which are more in accordance with the other data that they are more reliable.

The cooling curves of neutron stars, which are obtained from the fit of the data of different pulsars on the surface temperature versus age diagram, are given in various articles (see for example Yakovlev et al. 2002). Since the general form of the cooling curves given by different authors is approximately the same in all the works, we will use the data given in Yakovlev et al. (2002).

DITNS RX J1605.3+3249 (Motch et al. 1999) is located in a region of sky (l=53o, b=48o) where it is easier to determine the temperature with small uncertainty. There are 2 different distance estimations for this source, 0.1 kpc (Haberl 2003) and 0.3 kpc (Kaplan et al. 2003b). We may adopt reliable distance values for any source which has a spectrum close to the blackbody using luminosity values. The temperatures of DITNS RX J1605.3+3249 and of radio pulsars J1057-5226 and J0659+1414 which have similar soft X-ray spectrums are 0.092, 0.070 and $``$0.092 keV, respectively (Yakovlev et al. 2003). Therefore, the luminosity of RX J1605.3+3249 must not be smaller than the luminosity of these radio pulsars (Becker & Aschenbach 2002; Guseinov et al. 2003b), i.e. it must not be less than about 1033 erg/s. The luminosity of RX J1605.3+3249 at 0.1 kpc is 1.1$`\times `$1031 erg/s (Haberl 2003), so that, its luminosity at 0.3 kpc must be about 1032 erg/s. On the other hand, Geminga pulsar (B0633+17), which has temperature of about 0.045 keV, has Lx=1.05$`\times `$1031 erg/s in 0.1-2.4 keV band (Becker & Trumper 1997). This also suggests to adopt a larger Lx value for RX J1605.3+3249. So, we can adopt a distance value of 0.3 kpc (Kaplan et al. 2003b) or even a larger value which is more reliable.

Following the same path, we can adopt distance values for other DITNSs; the distances of J0420.0-5022, J0720.4-3125, J0806.4-4123, J1308.8+2127, and J214303.7+065419 must be about 3-4 times larger than the distance values given in Haberl (2003). All the new distance values and other reliable data are represented in Table 1.

## 4 Where must DITNSs be located on the P-Ṗ diagram?

There are 6 single radio pulsars (J1952+3252, J0659+1414, J0117+5914, J1057-5226, J0358+5413 and J0538+2817) with characteristic times in the range 105-6$`\times `$105 yr from which X-ray radiation has been observed (Becker & Aschenbach 2002; Possenti et al. 2002). These pulsars are located up to 2 kpc from the Sun. The luminosities of these radio pulsars are in the range 1031-5$`\times `$1033 erg/s (Guseinov et al. 2003b). There are also 3 radio pulsars (J0826+2637, J0953+0755 and J1932+1059) with characteristic times 3$`\times `$106$`<`$$`\tau `$$`<`$2$`\times `$107 yr (it is necessary to remember that value of $`\tau `$ may considerably exceed the real age value). The distances of these 3 radio pulsars are less than 400 pc and their luminosities are $`<`$1030 erg/s (Possenti et al. 2002).

As seen from Table 2, the X-ray luminosities of all the 8 DITNSs and Geminga are in the interval 1030-1.7$`\times `$1032 erg/s. Among these sources, Geminga has the lowest luminosity (Teff value of Geminga is also very small, see Table 1) and the $`\tau `$ value of Geminga is 3.5$`\times `$105 yr. Taking these data into consideration, the ages of DITNSs can be adopted as 105-106 yr in accordance with the age values calculated from the cooling models. Moreover, these sources are nearby objects and there is not any pulsar wind nebula around them nor any SNR to which they are connected; this also shows that the ages of these objects must be greater than 105 yr. On the other hand, pulsar wind nebula is present around the neutron stars with the rate of rotational energy loss Ė$`>`$5$`\times `$1035 erg/s and with Lx(2-10 keV)$`>`$5$`\times `$1032 erg/s (Guseinov et al. 2003b). Naturally, it may be possible to observe pulsars with smaller values of Ė and Lx (2-10 keV) which are located closer to the Sun. Taking these facts into account, we can assume that Ė values of DITNSs must be less than about 3$`\times `$1035 erg/s (constant Ė=3$`\times `$1035 erg/s line is shown on the P-Ṗ diagram, see Figure 1).

It is well known that the difference between $`\tau `$ and the real age can be significant for very young pulsars (Lyne & Graham-Smith 1998). On the other hand, for single-born old pulsars $`\tau `$ must be approximately equal to the real age if the evolution takes place under the condition B=constant. But none of the DITNSs is connected to a SNR, so that, DITNSs must be located in the belt between $`\tau `$=3$`\times `$104 yr and $`\tau `$=106 yr lines on the P-Ṗ diagram if the condition B=constant is satisfied.

For 5 of the 9 DITNSs (including Geminga) represented in Table 1, the spin periods (P) have been measured. Four of these objects have spin periods greater than 8 s and Geminga pulsar has P=0.237 s. As known, the period values of anomalous X-ray pulsars (AXPs) and soft gamma repeaters (SGRs) are greater than 5 s (see for example Mereghetti 2001; Guseinov et al. 2003a and the references therein). The small value of Ṗ=5.4$`\times `$10-13 s/s belongs to AXP 1E2259+586 which has P=6.98 s. Therefore, all single neutron stars with P$`>`$10 s must be related to AXPs and SGRs. The P=10 s line is displayed in Figure 1 to show the two separate regions in which DITNSs can be located. The DITNSs with $`\tau `$$`<`$106 yr may have values of P$`>`$10 s (see Table 1) and this is possible if the values of Ṗ are very large. Therefore, these objects may be the evolutionary continuations of SGRs and AXPs. Naturally, their birth rates should be in agreement with each other for this assumption to be true.

The locations of pulsars Geminga, RX J0720.4-3125 and RX J1308.8+2127 on the P-Ṗ diagram are shown in Figure 1. From the position of Geminga pulsar it is seen that this pulsar evolves similar to the radio pulsars with B=1012-1013 G. The position of pulsar RX J0720.4-3125 is not within our chosen interval of $`\tau `$=104-106 yr, but as this pulsar has a high value of kTeff (see Table 1) its real age (according to the cooling models) must be smaller than its $`\tau `$ value. So, there may be magnetic field decay or some other reason for this pulsar (i.e. n$`>`$3, where n is the braking index). Pulsar RX J1308.8+2127 is located in the SGR/AXP region on the P-Ṗ diagram so that this pulsar seems to have a relation with the SGR/AXP class of neutron stars. Pulsars RX J0806.4-4123 and RX J0420.0-5022 with ages $`<`$106 yr have large values of P. Although their Ṗ values are not known, they must be located on the upper part of the P-Ṗ diagram and their positions must not be lower than the position of J0720.4-3125, if the P values are correct.

## 5 Birth rates of different types of neutron star near the Sun

In section 2, we have mentioned that the birth rate for all the types of neutron star, in other words the supernova rate, must be about 11 in 106 yr in the region up to 400 pc from the Sun. In the same region, the birth rate of radio pulsars can be about 3-4 in 106 yr.

In Figure 1, we have plotted all the 9 radio pulsars with $`\tau `$$``$4$`\times `$104 yr which are connected to SNRs and located at distances up to 3.5 kpc (Guseinov et al. 2003c). There are also 2 other radio pulsars, J1048-5832 and J1837-0604, in this region with such values of $`\tau `$ but without any connection with SNRs (Guseinov et al. 2003b). From these data it follows that the birth rate of radio pulsars is 3.6 in 106 yr in the region up to 400 pc. It is necessary to take into consideration that for such young radio pulsars the beaming factor is close to 1 and the influence of the luminosity function on the estimations is small. Note that the searches of pulsars near the Sun in the central regions of SNR shells and the searches of pulsar wind nebulae are considerably better than the searches of pulsars under the surveys. Also note that in some cases pulsars have been found after observing point X-ray sources in the central parts of SNRs. Therefore, we can adopt that the birth rate in the region up to 400 pc from the Sun is about 3-4 radio pulsars in 106 yr.

It is necessary to take into account that in the region with distance up to 3.5 kpc from the Sun, there are also 6 dim radio quiet neutron stars (DRQNSs) which are connected to SNRs (Table 1). The locations of 2 of these objects, 1E1207.4-5209 and RXJ0002+6246, are shown in Figure 1. These objects have considerably large values of $`\tau `$ compared to the ages of the SNRs in which they are located, so that, they have different evolutionary tracks compared to other radio pulsars which are connected to SNRs. If we assume that all DRQNSs are also radio pulsars with 1012$`<`$B$`<`$1014 G, then they must have low radio luminosities and/or the direction of their radio radiation does not pass through the line of sight. On the other hand, all of them may have large P and $`\tau `$ values and significantly different evolutionary tracks compared to other pulsars with the same magnetic field, because there exist some other important differences between DRQNSs and most of the radio pulsars. Practically, all the radio pulsars which are connected to SNRs and which have similar ages as DRQNSs have pulsar wind nebula (Guseinov et al. 2003b). None of the DRQNSs has such property. Therefore, these DRQNSs have Ė$`<`$3$`\times `$1035 erg/s. The positions of 1E1207.4-5209 and RX J0002+6246 on the P-Ṗ diagram require magnetic field decay or some additional ideas for the pulsar models (in Fig.1 the location of RX J0002+6246 has been found using the condition $`\tau `$ = age of the SNR). Therefore, the number of SNRs which contain ordinary and other types of pulsars with similar properties for d$``$3.5 kpc is 17. By this approach, we give the upper limit for pulsar birth rate in the region up to 400 pc which is 5.5 in 106 yr.

In section 3, we have adopted up to 3-4 times larger distance values for DITNSs and this gives us the possibility to estimate the birth rate of this type of neutron star as 9 in 106 yr in the same region with a radius of 400 pc. These sources have P$`>`$0.1 s, Ė$`<`$3$`\times `$1035 erg/s and, according to the cooling theories, ages between 3$`\times `$104-106 yr that they must be in later stages of the evolution of single neutron stars with initial magnetic field B$`>`$1012 G (see the locations of these sources in Figure 1). But how can we explain such a large birth rate for these sources which is comparable with the supernova rate? First, note that the statistical data are poor. Second, the actual ages of the SNRs may be not up to 4$`\times `$104 yr but up to 3$`\times `$104 yr (see Table 1). On the other hand, the rate of supernova must be a little more if we take into account the SNRs which have low surface brightness during their evolution. In this case, the birth rate of DITNSs must roughly be equal to the birth rate of radio pulsars and DRQNSs together.

As seen from Figure 1, there exist 5 single radio pulsars which have been detected in X-ray band and have 105$`<`$$`\tau `$$`<`$6$`\times `$105 yr and 5$`\times `$1011$`<`$B$`<`$1.1$`\times `$1012 G. They are located in the region up to 2 kpc from the Sun and 2 of them, J0659+1414 and J0538+2817, are most probably connected to S type SNRs (Kramer et al. 2003; Guseinov et al. 2003b and the references therein). If we also consider that 2 such pulsars may go far away from the Galactic plane that they can be missed in the surveys, the total number of the pulsars in the considered volume turns out to be 7. Therefore, the birth rate of this type of pulsar in the region with d$``$400 pc is not more than 0.6 in 106 yr.

From the estimations done in this section, we see that approximately 60% of the SNRs with surface brightness $`>`$10-21 Wm-2Hz-1sr-1 are connected to normal pulsars and DRQNSs. The rate of birth for DITNSs is also approximately equal to 60% of the rate of supernova explosion. Therefore, the neutron stars with ages approximately $`<`$5$`\times `$105 yr, which show themselves as radio pulsar or DRQNS in SNRs, may mainly transform to DITNSs.

The numbers of radio pulsars with effective values of magnetic field B$``$1013 G and B$``$3$`\times `$1012 G which have $`\tau `$$`<`$106 yr and d$``$3.5 kpc are 5 and 32, respectively (ATNF pulsar catalogue 2003; Guseinov et al. 2002). From these data, the birth rates of radio pulsars with B$``$1013 G and B$``$3$`\times `$1012 G located up to 400 pc from the Sun must be considerably more than 0.06 and 0.4, respectively, in 106 yr, because the $`\tau `$ values may be several times larger than the real ages.

In the region up to 8 kpc from the Sun, there are 4 AXPs and one SGR and the ages of these objects must not be larger than 5$`\times `$104 yr (Mereghetti 2001; Guseinov et al. 2003a). Therefore, the birth rate of AXPs and SGRs in the considered cylindrical volume with the radius of 400 pc in 106 yr is not less than 0.15. This is about 60 times smaller than the supernova rate in the same volume and the time interval, but it approximately coincides with the birth rate of radio pulsars with effective value of B$``$1013 G. As seen from Figure 1, only 3 pulsars, J1740-3015, J1918+1444 and J1913+0446 (in the order of increasing value of P), are located in the volume with a radius of 3.5 kpc and with $`\tau `$$``$105 yr.

## 6 Discussion and Conclusions

It is known that, different types of neutron star which have different properties are born as a result of supernova explosion: radio pulsars, DRQNSs, DITNSs, AXPs and SGRs. There exist PWN and often SNR shell around very young single radio and X-ray pulsars which have Ė$`>`$5$`\times `$1035 erg/s (Guseinov et al 2003b). There may exist only the shell around pulsars with Ė$`<`$3$`\times `$1035 erg/s, DRQNSs and AXPs the ages of which are less than 105 yr. No shell nor PWN has been found around some of the very young pulsars (e.g. J1702-4310 and J1048-5832) and SGRs. Absence of the shell or PWN around DITNSs must be considered normal because they have considerably large ages.

In section 4, we have adopted $``$3-4 times larger distance values, compared to the values usually adopted (Haberl 2003), for most of the DITNSs and this gives the possibility to decrease their birth rate down to the sum of the birth rates of radio pulsars and DRQNSs. Birth rate of these 3 types of neutron star together is approximately equal to the rate of supernova. The combined birth rate of these 3 types of neutron star may be more than the supernova rate, because some of the DITNSs may be formed as a result of the evolution of DRQNSs and radio pulsars. Birth rate of AXPs and SGRs, which belong to the same class of objects, is about 60 times smaller than the supernova rate and it is about same as the birth rate of radio pulsars with effective values of B$``$1013 G.

As seen from Table 1, the period values of 4 DITNSs are very large, though their ages are smaller than 106 yr. From this situation, there arises a possibility of a relation between some of the DITNSs and AXPs/SGRs. Naturally, we must take into account the position of each DITNS on the P-Ṗ diagram to show the relation between some of these objects and AXPs/SGRs.

Existence of radio pulsars with n$`<`$3 and with real ages smaller than $`\tau `$ (for young pulsars) show that the condition B=constant is not satisfactory in all cases and this is well seen in Figure 1. Most of the pulsars with 105$`<`$$`\tau `$$`<`$107 yr are not in the belt B=1012-1013 G where most pulsars are born in; often, there occurs magnetic field decay. But the evolutions of AXPs/SGRs and some of the DITNSs according to the field decay approach (the magnetar model) lead to bimodality in the number of neutron stars versus the magnetic field distribution. On the other hand, the time scale of the magnetic field decay must be very short. This shows that the large effective B values of these objects and the shape of their evolutionary tracks must be related mainly to the masses of and the density distributions in the neutron stars and also to the activity of the neutron star. This is also necessary to understand the different positions of radio pulsars which are connected to SNRs and of DRQNSs on the P-Ṗ diagram despite the fact that they have similar ages. If the evolution under the condition B=constant were true, then they would be located along $`\tau `$=constant belt, but not along the constant magnetic field belt.

We think that there is a possibility to get rid of these difficulties and to understand the large X-ray luminosities and also the bursts of AXPs/SGRs. It is necessary to assume the birth of neutron stars with masses about half of the maximum mass values found from the given equations of state and rotational moment. In principle, it must be easy to identify such smaller mass neutron stars as they are far away from hydrodynamical equilibrium. They must have an elipsoidal shape due to rotation and possibly they do not rotate as a rigid body. In this case, the young pulsar may especially demonstrate itself when the angle between the magnetic field and the rotation axes is close to 90o. In such a case, a considerably larger effective value of magnetic field can be produced as compared to the real magnetic field.

| Table 1 - The data of DRQNSs and DITNSs. Ages of the | | | | | | | |
| --- | --- | --- | --- | --- | --- | --- | --- |
| SNRs connected to these objects are given in the 4th | | | | | | | |
| column. References are shown at the end of the table. | | | | | | | |
| Some data adopted in this work are given without reference. | | | | | | | |
| Names | P | $`\tau `$ | t | d | kTeff | L$`_\text{x}`$1032 | L$`_\text{x}`$ |
|  | s | kyr | kyr | kpc | keV | erg/s |  |
| 1E1207.4-5209 | 0.424 | 340- | 7-20 | 1.8 | 0.11 | 10d$`{}_{}{}^{2}{}_{2}{}^{}`$ | $``$0.1 |
| G296.5+10.0 |  | 480 |  | 2 |  | (0.5-6) |  |
| S |  |  |  |  | 0.25 |  |  |
|  |  |  |  | 2.1 | (0.5-6) |  |  |
|  |  |  |  |  |  |  |  |
| 1E0820-4247 |  |  | 3-4 | 2 | 0.15 | 12 |  |
| RXJ0822-4300 |  |  |  |  |  | (0.1-2.4) |  |
| Puppis A |  |  |  |  | 0.44 ? |  |  |
| G260.4-3.4 |  |  |  |  |  |  |  |
| S |  |  |  |  |  |  |  |
| CXOJ2323+5848 |  |  | 0.32 | 3.2 |  |  |  |
| Cas A |  |  |  |  |  |  |  |
| G111.7-2.1 |  |  |  |  |  |  |  |
| RXJ0002+6246 | 0.2418 |  | 10-20 | 3.5 | 0.10 | 2 | 0.0002- |
| $`\gamma `$-ray source |  |  |  |  |  | (0.5-2) | 0.0006 |
| G117.7+0.6 |  |  |  |  |  |  |  |
| RXJ0007.0+7302 |  |  | 10-24 | 1.4 |  | 0.15 |  |
| G119.5+10.2 |  |  | \[5,8-10\] |  |  | (0.1-2.4) |  |
| CTA 1 |  |  |  |  |  |  |  |
| RXJ2020.2+4026 |  |  | 6-10 | 1.5 |  | 9 | 0.0009- |
| $`\gamma `$-ray source |  |  | \[4-6\] |  |  |  | 0.00009 |
| G78.2+2.1 |  |  |  |  |  |  |  |
| $`\gamma `$ Cyg |  |  |  |  |  |  |  |
| CXOJ0852-4615 |  |  | 0.7-2 | 1 | 0.40 | 2.3 |  |
| G266.2-1.2 |  |  |  |  |  | (0.4-6) |  |
|  |  |  |  |  |  | d=1kpc |  |
|  |  |  |  |  |  |  |  |
| Geminga | 0.237 | 350 |  | 0.16 | 0.048 | 0.01 | 0.00003 |
| B0633+1748 |  |  |  |  |  | (0.6-5) |  |
|  |  |  |  |  |  |  |  |

| Table 1 (continued) | | | | | | | |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Names | P | $`\tau `$ | t | d | kTeff | L$`_\text{x}`$1032 | L$`_\text{x}`$ |
|  | s | kyr | kyr | kpc | keV | erg/s |  |
| RXJ1836.2+5925 |  |  |  | 0.43 | 0.043 | 0.054 |  |
|  |  |  |  |  |  |  |  |
| RXJ1856.5-3754 |  |  |  | 0.13 | 0.062 | 0.17 |  |
|  |  |  |  |  |  |  |  |
| RXJ0720.4-3125 | 8.39 | 3000 |  | 0.3 | 0.085 | 1.7 | 55 |
|  |  |  |  |  |  | 0.9d$`{}_{}{}^{2}{}_{0.3}{}^{}`$ |  |
|  |  |  |  |  |  |  |  |
| RXJ0420.0-5022 | 22.7 ? |  |  | 0.4 | 0.057 | 0.43 |  |
|  |  |  |  |  |  |  |  |
| RX J0806.4-4123 | 11.37 ? |  |  | 0.4 | 0.095 | 0.91 |  |
|  |  |  |  |  |  |  |  |
| RX J1308.8+2127 | 10.3 | 7-14 |  | 0.4 | 0.091 | 0.82 | $``$0.2-0.6 |
| RBS 1223 |  |  |  | 0.67 |  |  |  |
|  |  |  |  |  |  |  |  |
| 1RXSJ214303.7 |  |  |  | 0.3- | 0.09 | 1-1.7 |  |
| +065419 |  |  |  | 0.4 |  |  |  |
| RBS 1774 |  |  |  |  |  |  |  |
| RXJ1605.3+3249 |  |  |  | 0.3- | 0.092 | 1-1.7 |  |
|  |  |  |  | 0.4 |  |  |  |
|  |  |  |  | 0.3 |  |  |  |
|  |  |  |  |  |  |  |  |
|  Yakovlev et al. 2002;  Guseinov et al. 2003c;  Haberl 2003; | | | | | | | |
|  Bignami et al. 2003;  De Luca et al. 2002;  Vasisht et al. 1997; | | | | | | | |
|  Pavlov et al. 2002a;  Petre et al. 1996;  Brazier & Johnston 1999; | | | | | | | |
|  Pavlov et al. 2002b;  Hailey & Craig 1995;  Brazier et al. 1996; | | | | | | | |
|  Kargaltsev et al. 2002;  Zane et al. 2002;  Kaplan et al. 2003a; | | | | | | | |
|  Hambaryan et al. 2002;  Haberl & Zavlin 2002;  Haberl et al. 2003; | | | | | | | |
|  Kaplan et al. 2002;  Kaplan et al. 2003b;  McLaughlin et al. 1999; | | | | | | | |
|  Haberl et al. 1999;  Halpern & Wang 1997;  Pavlov et al. 2001; | | | | | | | |
|  Guseinov et al. 2003d. | | | | | | | |

Figure Caption
Figure 1: Period versus period derivative diagram for different types of pulsar. The ’+’ signs denote the radio pulsars with d$``$3.5 kpc which are connected to SNRs. The ’X’ signs show the positions of the radio pulsars with d$``$3.5 kpc and 105$`<`$$`\tau `$$`<`$2$`\times `$107 yr which have been detected in X-rays. The locations of 3 radio pulsars which have d$``$3.5 kpc and $`\tau `$$`<`$105 yr are shown with ’circles’ to make a comparison between the birth rates (see text). DITNSs are represented with ’stars’ and DRQNSs are displayed with ’empty squares’. The ’filled squares’ show the positions of all AXPs/SGR in the Galaxy. Names of DITNSs, DRQNSs, 2 of the AXPs, and some of the radio pulsars are written. Constant lines of B = 1011-15 G, $`\tau `$ = 103-9 yr, and Ė = 1029, 1032, 1035, 3$`\times `$1035 and 1038 erg/s are shown. P=10 s line is also included (see text).


format: text

text:
# Molecular gas in the galaxy M 83

## 1 Introduction

In disk galaxies the mass of the molecular (H2) gas is usually similar to, or a bit less than, that of the atomic (HI) gas (Young & Scoville, 1991; Casoli et al., 1998). However, their spatial distributions are completely different, both in the radial and in the vertical directions. The molecular gas is concentrated to a relatively thin disk, which often has a radial extent similar to that of the optical disk and a vertical scale height which is smaller than that of the young, massive stars. The radius of the HI disk is usually much larger than that of the optical disk (up to 5 times the Holmberg radius), and the atomic gas is also less concentrated to the disk in the vertical direction (Knapp, 1987).

If the mass surface density in a region exceeds a certain critical density, massive star formation is triggered (Kennicutt, 1989). So, naturally, the study of interstellar gas, and especially the molecular gas phase, plays a key role in understanding the structure and evolution of galaxies. Unfortunately, direct studies (in emission) of the cold ($``$10 K) H2 gas in molecular clouds are very difficult, since the energy of the first excited rotational level is of the order of 500 K (corresponding to a wavelength of 28 $`\mu `$m), and therefore much larger than the average kinetic temperatures in the clouds. Carbon monoxide (CO) is the second most common molecule with a relative abundance with respect to H2 of about 10-4. The heavy, asymmetrical CO molecule has rotational levels which are easily excited by collisions with H2, and it radiates at mm wavelengths where observations can be carried out using radio telescopes. However, the conversion factor which is used to derive the H2 column density from the velocity-integrated CO intensity, $`X_{\mathrm{CO}}`$, is not well-known, and this is a source of major uncertainty (Young & Scoville, 1991; Combes, 1991).

CO emission at mm wavelengths has been detected in many galaxies ranging in distance from the nearby irregular galaxies The Large and The Small Magellanic Clouds to quasars at a redshift of more than $`z`$$``$ 5, and in essentially all galaxy types. A large survey of 300 galaxies (of which 236 were detected) of different types was done by Young et al. (1995) using the FCRAO 14 m telescope. The analysis of the data is still ongoing, but some interesting results have been published, e.g., the conversion factor ($`X_{\mathrm{CO}}`$) appears relatively constant, at least within each galaxy type (Young & Rownd, 2001), and the global star formation efficiency decreases with increasing galaxy size (Young, 1999). In this survey a number of galaxies were also mapped and radial distributions of the molecular gas can be examined. Recently, in a more detailed study, 28 nearby galaxies were mapped with the 45 m Nobeyama telescope (Kuno et al., 2000; Nishiyama et al., 2001). Among other things, this showed that barred and non-barred galaxies differ in both the distribution and the kinematics of the molecular gas. Another large survey is the BIMA SONG project (Regan et al., 2001; Helfer et al., 2003), in which 44 nearby ($`v_{\mathrm{}}<2000`$ km s-1) galaxies have been observed with the BIMA array and the NRAO 12m telescope. Early results show that many galaxies have multiple peaks at the nucleus, that the CO distribution usually cannot be fitted with a single exponential disk model, and that the ratio of the scale length in the CO disk and the scale length in the stellar disk is on average 1 (but the scatter around this average is large).

Detailed morphological studies of selected galaxies are more rare, partly because there are so few nearby galaxies for which the limited angular resolution of radio telescopes gives an adequate spatial resolution. Nakai et al. (1994) outlined a number of important observational conditions that should be fulfilled in order to be able to extract useful information on the properties of the molecular gas in different regions of a galaxy: (1) a spatial resolution better than one kpc, which is the typical width of molecular spiral arms and the spacing between them, (2) high sensitivity in order to detect the weak inter-arm emission, (3) mapping of a large field in order to get statistically meaningful results, and (4) filled aperture telescope data, since it is important to detect the diffuse emission in the inter-arm regions. To this list we add a fifth condition, (5) the observation should be carried out in (at least) two lines per molecular species (and preferably also partially observed in a few isotopomers) in order to gain insight into the density, temperature, and chemistry of the molecular gas.

Single-dish measurements are therefore essential to our determination of the molecular gas distribution, its relation to star formation activity, its location with respect to dynamical resonances, and the total reservoir of gas capable of forming future generations of stars. The most detailed single-dish maps of molecular gas in spiral galaxies are the mappings of M 51, using the IRAM 30m telescope (Garcia-Burillo et al., 1993) and the Nobeyama 45m telescope (Nakai et al., 1994), and two projects currently underway to map our neighbor M 31 (Loinard et al., 1999; Nieten et al., 2000). In M 51 the molecular spiral arms were resolved, and streaming motions were detected. Also, the CO emission was found to peak on the inside of the spiral arms, giving support to a scenario where molecular gas agglomerates and forms new stars, which in their turn produce H$`\alpha `$ emission and other star formation tracers, close to the spiral arms. Support was also found for the proposal that gravitational instability in the disk is the main formation mechanism of GMAs (Galactic Molecular Associations of mass $``$107 M). Other large spiral galaxies that have been completely surveyed in CO with a single-dish telescope are IC 342 (Crosthwaite et al., 2001), M 101 (Kenney et al., 1991), NGC 253 (Houghton et al., 1997), and NGC 6946 (Tacconi & Young, 1989). Some large scale CO mapping projects, such as the OVRO maps of M 51 (Aalto et al., 1999) and M 83 (Rand et al., 1999), have been done using mm-wave interferometers. Although this gives a superior angular resolution and important unique information, a significant part of the emission is resolved out (especially the diffuse component in the inter-arm regions).

We have mapped the CO($`J`$=1–0) and CO($`J`$=2–1) emission in the barred, grand-design, spiral galaxy M 83 located at a distance of 4.5 Mpc. The observations were done with the 15 m Swedish-ESO Submillimetre Telescope (SEST) providing beam widths of 45″ (1–0) and 23″ (2–1), which corresponds to 980 pc and 500 pc, respectively, at the adopted distance. M 83 is among the brightest galaxies in terms of CO emission, and it is also one of the most nearby barred spirals. Its low inclination ($``$24) and proximity make it a perfect candidate to investigate the correlation between various tracers of stellar activity.

In this paper we present the CO($`J`$=1–0) and CO($`J`$=2–1) observational data. We show the distribution of the molecular gas, as inferred from the CO data, and compare it to the distributions of optical light, obtained in different filters, and the HI gas. We examine the CO ($`J=21`$)/($`J=10`$) line ratio and its variation over the disk. In two forth-coming papers we will present and discuss the kinematics of the molecular gas, and discuss the relation between the molecular gas and star formation.

## 2 The galaxy M83

M 83 is a relatively nearby, barred, grand-design, spiral galaxy (see Table 1). It is viewed almost face-on and the bar is aligned with the major axis. Color images show clumpy, well-defined spiral arms, evidently rich in young blue stars. It is also fairly symmetrical and has no nearby massive optical companions and no immediate evidence of interaction or outflows.

Due to its favorable properties it is one of the most well-observed galaxies: HI (Rogstad et al., 1974; Huchtmeier & Bohnenstengel, 1981; Tilanus & Allen, 1993), radio continuum (Ondrechen, 1985; Neininger et al., 1993), CO (see Table 2), IR continuum (Adamson et al., 1987; Fitt et al., 1992; Elmegreen et al., 1998), optical continuum (Jensen et al., 1981), Balmer lines (Comte, 1981; de Vaucouleurs et al., 1983; Tilanus & Allen, 1993), and in X-ray emission (Immler et al., 1999; Soria & Wu, 2002).

The distance to M 83 is a long-standing issue. Sandage & Tammann (1974) used a relation between the luminosity class and the angular size of the three largest Hii regions to deduce a distance of 8.9 Mpc, while de Vaucouleurs (1979) used a number of methods to derive a distance of 3.7 Mpc. Recently, Cepheids have been observed in M 83 by Thim et al. (2003) using the VLT. Twelve Cepheid candidates were observed and the distance was estimated to be 4.5$`\pm `$0.3 Mpc. This distance is also supported by studies that seem to indicate that M 83 interacted with NGC 5253 some 1–2 $`\times `$ 109 years ago (Kobulnicky & Skillman, 1995; Calzetti et al., 1999). This small, metal-poor dwarf galaxy about 2° SE of M 83 has roughly the same systemic velocity, and shows a highly disturbed HI velocity field and intense star formation. There has been a number of distance measurements to NGC 5253 (Saha et al., 1995; Parodi et al., 2000; Gibson et al., 2000), and they all favor a distance of about 4 Mpc.

## 3 Observations and data reduction

### 3.1 Observations

The CO($`J=10`$ and $`J=21`$) observations of M 83 were done using the 15 m SEST111The Swedish-ESO Submillimetre Telescope is operated jointly by ESO and the Swedish National Facility for Radio Astronomy, Chalmers University of Technology during two epochs, 1989-1994 and 1997-2001 \[a description of the telescope is given by Booth et al. (1989)\]. Between 1988 and 1994 we observed CO($`J`$=1–0) in the inner disk of M 83. A segment of these observations were presented in Wiklind et al. (1990). In 1998 and 1999 we completed this map by observing the outer disk. During the latter observations it was possible to do simultaneous observations in two frequency bands. We used this option to observe the CO($`J`$=2–1) line in parallel with the CO($`J`$=1–0) line. The map spacing (11″) and integration times were set by the requirements on the CO($`J`$=1–0) data. In 2000 and 2001 our primary task was to fill in the central area of the CO($`J`$=2–1) map. For these observations we chose a grid spacing of 7″. We used the other receiver to observe the 13CO($`J`$=1–0) line, but these data will be presented elsewhere.

The receivers were centered on the CO($`J`$=1–0) and CO($`J`$=2–1) lines (115.271 and 230.538 GHz, respectively), where the full half-power beam widths (HPBWs) are 45$`\mathrm{}`$ and 23$`\mathrm{}`$, respectively. Observations performed during 1988-1994 were made with a cooled Schottky-diode mixer receiver, while the subsequent observations were done with SIS receivers. All receivers operated in the single-sideband mode, and the average system temperatures, corrected for the atmospheric and antenna ohmic losses, were 660 K and 390 K for the CO($`J`$=1–0) observations in the first and second epoch, respectively, and 340 K for the CO($`J`$=2–1) observations. The backends were acousto-optical spectrometers, with channel separations of 0.69 MHz (1.8 and 0.9 km s-1 at the two rest frequencies) and total bandwidths varying between 0.5 and 1 GHz.

We used the dual beam switching mode with a throw of $``$12$`\mathrm{}`$ in azimuth. In dual beam switching the source is placed first in the signal beam and then in the reference beam. The two spectra produced are then subtracted to generate a spectrum with a minimum of baseline variation. Great care was taken to ensure the quality of the data. No observations were done below the elevation of 30° or above 84° (M 83 passes within 1° from zenith at SEST), and pointing was regularly checked and updated using the SiO($`v`$=1, $`J=21`$) maser line of the nearby AGB-star W Hya ($``$3 from M 83). Total pointing errors were typically less than 3″, and data with pointing offsets suspected to be larger than 4″ were disregarded in the subsequent data reduction (when possible, map points with such data were reobserved). The intensity scale was calibrated using the conventional chopper wheel method, and the internal calibration errors in the corrected antenna temperatures, $`T_\mathrm{A}^{}`$, is within $`\pm `$10% according to the SEST manual. In this paper we use the main beam brightness temperature scale $`T_{\mathrm{mb}}`$, which is defined by $`T_{\mathrm{mb}}T_\mathrm{A}^{}/\eta _{\mathrm{mb}}`$. According to the SEST manual the main beam efficiencies \[$`\eta _{\mathrm{mb}}`$(115 GHz)=0.7, $`\eta _{\mathrm{mb}}`$(230 GHz)=0.5\] have been constant throughout the period of our observations.

The observations were centered on the coordinates $`13^\mathrm{h}36^\mathrm{m}59\stackrel{s}{.}4`$, $`29\mathrm{°}52\mathrm{}05\mathrm{}`$ (J2000), which is the optical center given by de Vaucouleurs et al. (1976). For CO($`J`$=1–0) we have 3634 spectra in 1900 positions with 11″ spacing, and the coverage is complete out to a radius of 4′20″. The grid extends parallel to the equatorial coordinate system, and the dense spacing (1/4 of the HPBW) was chosen in order to facilitate the use of deconvolution techniques. The integration time per position was chosen to render an rms around 70 mK ($`T_{\mathrm{mb}}`$-scale) at the original velocity resolution (1.8 km s-1). During normal weather conditions we had a typical on-source integration time of 60 seconds. Among the final 1900 spectra we have $``$ 3$`\sigma `$ detections in 1201 positions, i.e. a detection rate of 63%. The average rms noise in our data set is 74 mK, which should be compared to the peak temperature which usually lies around 0.3 K ($`T_{\mathrm{mb}}`$-scale) in the spiral arms in the outskirts of the map.

In CO($`J`$=2–1) we have two maps with different spacings covering two, partly overlapping, regions. The inner 5′$`\times `$3′ is covered with a 7″ grid, while the rest of the optical disk (out to a radius of 4$`\stackrel{}{.}`$5) is covered with an 11″ grid. In total, we have spectra in 2574 positions, out of which we have $``$ 3$`\sigma `$ detections in 1898 positions (73% detection rate). The spectra have an average rms of about 90 mK at the original velocity resolution (0.9 km s-1). The observational parameters are summarized in Table 3.

### 3.2 Data reduction

The data reduction was carried out with class222Continuum and Line Analysis Single-dish Software from the Observatoire de Grenoble and Institut de Radio-Astronomie Millimétrique , drp333Data Reduction Package developed at Onsala Space Observatory, and xs444xs: Spectral Line Reduction Package for Astronomy developed at Onsala Space Observatory. During the CO($`J`$=1–0) observations, the number of channels and the frequency resolution varied between the spectrometers. We therefore resampled all data to 1.8 km s-1 (0.69 MHz) and by selecting the central 700 channels this gave us a (fully sufficient) total bandwidth of 1260 km s-1.

After flagging bad channels, and replacing them with a value obtained from interpolation between adjacent channels, all spectra at the same position were averaged using a weight determined by the noise level. The baselines were stable, and only linear polynomial fits were subtracted. From the spectra we created a fits data cube for each transition, where the first two axes represent spatial coordinates, and the third axis represents frequency (or velocity).

We adopted a relatively dense grid spacing, compared to the beamwidth, in order to be able to increase the spatial resolution of the data using a mem-deconvolution routine (maximum entropy method); the Statistical Image Analysis (SIA) routine (Rydbeck, 2000). This procedure is particularly effective on a data set as densely spaced as ours (spacing in the range 1/4–1/2 of the HPBW), due to the high redundancy. The routine starts by deconvolving a velocity-integrated intensity map, where the velocity range is chosen to be large enough to encompass all the emission. Once this is done, the procedure splits the initial velocity range into three equally wide regions, calculates their velocity-integrated intensities, and then deconvolves these three maps, using the deconvolved map from the previous step as the input guess. The deconvolution continues in this hierarchical manner, successively narrowing the velocity interval, until a velocity range of 5 km s-1 is reached. To deconvolve beyond this point is not meaningful since the S/N-ratios in the (narrow velocity range) maps drop below a critical value. For each of the two CO transitions we created 81 (34) maps per transition (each covering about 5 km s-1) which we assembled into a data cube per transition. By convolving the mem-deconvolved data cube by a Gaussian with a FWHM equal to the observational HPBW and comparing with the raw spectra, we verified that the mem-results are reliable, and that the raw data set is homogeneous, i.e., free from variations in gain and pointing errors. The angular resolutions of the mem-cubes are estimated to be $``$22″ and $``$14″ for the CO($`J`$=1–0) and CO($`J`$=2–1) data, respectively, and the velocity resolution is 5 km s-1.

We have also constructed convolved data sets by convolving the raw data cubes with a HPBW beam size of 20$`\mathrm{}`$ and 15$`\mathrm{}`$ for the CO($`J`$=1–0) and CO($`J`$=2–1) data, respectively (the spacing in these cubes are 11″). Even if this degrades the resolution by $``$10–20%, the average rms noise in these cubes decreases from 74 mK to 25 mK and from 90 mK to 31 mK, respectively. The dramatically decreased noise level is due to the fact that the spacing is small with respect to the HPBW of the convolution kernels. Since the peak intensity in the spectra are relatively unaffected by this process, the average S/N-ratio in the spectra increases by a factor of $``$3. As a result of this, the detection rate increases to 92% and 89% for the CO($`J`$=1–0) and CO($`J`$=2–1) data, respectively (cf Sect. 3.1). These cubes were used to cross-check the results we got from the mem-deconvolved data. Some results were also derived directly from the convolved cubes, especially in cases where the velocity resolution was more important than the spatial resolution. Furthermore, we convolved the raw CO($`J`$=2–1) cube with a beam size of 44″ and regridded it to 11″ spacing in order to compare with the convolved CO($`J`$=1–0) cube at the same resolution ($``$49″).

In order to avoid problems with baseline irregularities when calculating the zero moment (over velocity) of the intensities we used a sliding window technique where at each position integration was performed inside a velocity window whose center and width depend on the position within the galaxy. We checked carefully that the window always covered all the emission.

### 3.3 Correction for the error beam

A closer inspection of the CO($`J`$=2–1) spectra reveals, apart from the expected main component seen in the CO($`J`$=1–0) spectra, a broad (velocity width 400 km s-1) and faint (peak intensity around 0.03 K) component. Despite the low peak intensity of this feature, its integrated intensity is quite appreciable due to its large width. Typically, it makes up to 20%–35% of the total integrated intensity. This feature was discovered in the mem-deconvolved data, which are noise-free, but it can be seen also in spectra produced by convolving the raw CO($`J`$=2–1) data cube with a 45″ beam. The appearance of this broad component changed with position in the disk, and we found that its shape could, in any given position, be reproduced very well by convolving the data set with a relatively large Gaussian beam. This indicates that the SEST picks up emission from a much larger area than the HPBW of 23”, i.e., at 230 GHz the SEST suffers from a non-negligible error beam.

We have corrected the CO($`J`$=2–1) data sets for this error using a method described in Appendix A. We have found no trace of an error beam contribution in the CO($`J`$=1–0) data (not unexpectedly, considering the higher main beam efficiency), and consequently we have not applied any correction to these data.

## 4 The data

### 4.1 Spectra

In general, the spectra have a smooth and symmetric shape, but in some areas they show deviations. The spectra close to the nucleus are asymmetric, show extended wings, and a few of them have multiple peaks. The main reason for this is a large velocity gradient over spatial scales small compared to the telescope beam. Further out in the disk there are a few (interarm) regions which show multiple-peak spectra, such as the “Gould-belt structure” located about 3$`\mathrm{}`$ SE of the galaxy center (Comerón, 2001). However, in the vast majority of the observed positions the spectral feature can be well described by a single, Gaussian profile.

#### 4.1.1 The CO($`J`$=1–0) spectra

In the arms we find that the peak temperatures are typically 0.2–0.4 K, and in the inter-arm regions they are about 0.1 K. The spectra are widest in the nuclear region with (Gaussian) widths of 80–100 km s-1. Further out the linewidths are of the order of 7–14 km s-1 in the arms, with no systematic difference between arm and inter-arm regions.

The global CO($`J`$=1–0) spectrum (left panel in Fig. 1) has a velocity-integrated intensity of 8.18$`\pm `$0.04 K km s-1 (excluding the calibration error) and a noise level of 1.8 mK (all errors in this paper are 3$`\sigma `$ unless otherwise noted). The receding part of the galaxy (SW) has more CO($`J`$=1–0) emission than the approaching part (NE): the velocity interval 510–610 km s-1 is 18% more luminous than the 410–510 km s-1 interval. An excess of higher-velocity emission has also been seen in the atomic gas. Huchtmeier & Bohnenstengel (1981) found that the higher-velocity part of the HI emission is 16% more luminous than the lower-velocity part.

The peak intensity spectrum (Fig. 1, right panel) has been obtained by selecting, for every channel, the highest intensity in the map. Therefore, the noise level in this spectrum reflects the noise in the overall worst case. From studies of the channel maps we conclude that the peaks in this spectrum arise from the nuclear region (487, 514, and 534 km s-1), the western bar end (580 km s-1), the eastern bar end (460 km s-1), and the NE outer spiral arm (437 km s-1).

By adding the CO($`J`$=1–0) spectra (from the raw data cube) in concentric, 20″ wide annuli (compensated for inclination and position angle and centered on the IR nucleus, see Table 1), we produced the spectra in Fig. 3. They show how the average CO($`J`$=1–0) velocity-integrated intensity changes with galactocentric radius (which is indicated, in arc-seconds, in the upper right corner of each panel). The coverage in the map is complete to a radius of 229″.

#### 4.1.2 The CO($`J`$=2–1) spectra

The global and peak intensity CO($`J`$=2–1) spectra show characteristics very similar to those of the CO($`J`$=1–0) data, Fig. 2 \[since the CO($`J`$=2–1) spectra were obtained with two different spacings, we regridded these data to a common 11″ spacing before producing the spectra in Fig. 2\]. The main difference is that the peak intensity spectrum is brighter than the corresponding CO($`J`$=1–0) spectrum. This is most likely an effect of the emitting objects being smaller than the CO($`J`$=1–0) beam.

In Table 4 we summarize the characteristics of the global spectra. The global line intensity ratio $`R_{21}T_{\mathrm{mb}}(21)𝑑v/T_{\mathrm{mb}}(10)𝑑v`$ is 0.77$`\pm `$0.01, based on the errors from Table 4. Note that, including errors in the calibration and the main beam and error beam efficiencies, the result is 0.77$`\pm `$0.16. However, it should be noted that the relative calibration of the CO($`J`$=1–0) and CO($`J`$=2–1) data may be better than the absolute calibration and hence the error in their ratio can be smaller than indicated by the absolute error.

### 4.2 The maps

#### 4.2.1 CO($`J`$=1–0)

Figure 4 shows maps of the CO($`J`$=1–0) velocity-integrated intensity (hereafter $`I_{10}`$) and the peak intensity (hereafter $`T_{10}`$) for both the convolved data set (left column) and the mem-deconvolved data set (right column). The CO($`J`$=1–0) emission is concentrated to the nucleus, the central part of the bar, the bar ends, and the spiral arms. The velocity-integrated intensity at the nucleus and the eastern and western bar ends are: 73, 24, and 30 K km s-1, respectively, which is very close to the corresponding values in Crosthwaite et al. (2002): 73, 22, and 28 K km s-1. In the mem-deconvolved data set these features are brighter due to the increased spatial resolution: at the eastern bar end we measure 49 K km s-1, at the western bar end 52 K km s-1, and at the nucleus the maximum $`I_{10}`$ is 235 K km s-1 ($`T_{10}`$=2.8 K). By fitting an ellipse to the mem-deconvolved $`I_{10}`$-distribution we find that the CO peak coincides with the IR nucleus within the absolute position accuracy of our maps (3″). The axes ratio of this ellipse is 60″:29″ with a position angle of 36°. In the mem-deconvolved maps the spiral arms are resolved, and they can be followed for almost 360°. The emission in the spiral arms breaks up into regularly spaced maxima, which we interpret to be giant molecular associations (GMAs), a phenomenon previously seen in e.g. M 51 (Nakai et al., 1994). The outer arm, to the SE of the nucleus, seems somewhat disturbed in the sense that in the mem-deconvolved $`T_{10}`$ map it is possible to discern two parallel arms. Further toward the west, this arm appear to almost merge with the inner arm, which also appears to be the case on the opposite side of the galaxy (NE of the nucleus).

Rand et al. (1999) (hereafter RLH) mapped the CO($`J`$=1–0) emission in the inner eastern arm using the OVRO Millimeter Array (spatial resolution 6$`\stackrel{}{.}`$5 $`\times `$ 3$`\stackrel{}{.}`$5). In Fig. 5 we show their data overlayed on our deconvolved $`I_{10}`$ map. The components close to the nucleus (components 1–5 using their nomenclature), are not prominent in our data. They are distributed along the leading edge of the bar. The reason is probably that a shock in this region creates substructures (in temperature, in density, and/or in velocity) in the CO distribution which elsewhere is smooth and therefore undetectable by an interferometer (RLH found that only 2–5% of the single-dish emission was recovered in their interferometer observations). Our eastern bar end is located between components 6 and 7, and the spiral arm segment which continues southward is located between components 9 and 10.

#### 4.2.2 CO($`J`$=2–1)

The CO($`J`$=2–1) maps show structures very similar to those in the CO($`J`$=1–0) maps, however with a higher angular resolution, Fig. 6. The deconvolved $`I_{21}`$ map shows the “double nucleus” previously seen in the CO($`J=32`$) and CO($`J=43`$) lines by Petitpas & Wilson (1998). This “double nucleus” is apparent also in our raw CO($`J`$=2–1) data, but it is not possible to see in the convolved $`I_{21}`$ map, since this map is made using data which was both convolved (with a 15″ beam) and regridded (to 11″ spacing). Since the gray-scale and contour levels in Figs. 4 and 6 are the same, we can directly infer that the CO($`J`$=2–1) emission is less concentrated to the arms than the CO($`J`$=1–0) emission.

## 5 Results and discussion

### 5.1 Molecular mass in M 83

The H2 column density in external galaxies is commonly derived from the velocity-integrated CO($`J`$=1–0) intensity ($`I_{\mathrm{CO}}=T_{\mathrm{mb}}𝑑v`$) using the expression

$$N(\mathrm{H}_2)=X_{\mathrm{CO}}\times I_{\mathrm{CO}},$$
(1)
where $`X_{\mathrm{CO}}`$ is a conversion factor. The actual value of this factor, and its applicability in different environments have been widely discussed, see e.g. Wilson (1995); Arimoto et al. (1996); Combes (2000); Young & Rownd (2001). In this paper we have adopted the value $`X_{\mathrm{CO}}`$= 2.3$`\times 10^{20}`$ (K km s-1)-1 cm-2 (Strong et al., 1988), which is supported by the analysis done in RLH. In order to convert the measured intensities into mass surface density we used the expression

$$\mathrm{\Sigma }(\mathrm{H}_2)=3.7\times I_{\mathrm{CO}}\times \mathrm{cos}i=3.3\times I_{\mathrm{CO}}\mathrm{M}_{}\mathrm{pc}^2,$$
(2)
where i is the inclination of the galaxy, and $`I_{\mathrm{CO}}`$ is given in K km s-1 (This expression can be obtained by integrating Eq. 1). Using this conversion factor to convert the integrated intensity of the global CO($`J`$=1–0) spectrum (8.18 K km s-1, Tab. 4) we find that the average (over the entire optical disk, and including the bar and the nucleus) surface mass density is 27.2 M pc-2, and that the total molecular gas mass within the surveyed area is 3.3$`\times 10^9`$ M.

### 5.2 The radial distribution

#### 5.2.1 Molecular gas

We have, using the sliding-window technique described above, calculated the velocity-integrated intensity, centroid velocity, peak intensity, and center and width of a Gaussian emission profile fitted to the spectrum for all spectra in all data cubes. The left panel of Fig. 7 shows the CO($`J`$=1–0) velocity-integrated intensities (obtained from the convolved data set) as a function of the (inclination-corrected) distance from the galaxy center. The scatter is dominated by the arm-interarm contrast. The data in the right panel of Fig. 7 were obtained from the integrated intensities of spectra averaged within concentric annuli (see Fig. 3 and Table 5) by converting the CO($`J`$=1–0) integrated intensity into mass surface density as described above. In this procedure we used the raw data set, since we needed neither the higher spatial resolution of the mem-deconvolved spectra nor the lower noise of the convolved spectra. The error bars reflect the uncertainty in the azimuthally-averaged radial distribution. To these data we fitted an H2 mass surface density distribution which is a combination of a Gaussian component and an exponential function (for the disk),

$`\mathrm{\Sigma }(\mathrm{H}_2)`$ $`=`$ $`216\mathrm{exp}\left[\left({\displaystyle \frac{r}{649\mathrm{pc}}}\right)^2\right]+`$ (3)
$`120\mathrm{exp}\left[{\displaystyle \frac{r}{2265\mathrm{pc}}}\right]\mathrm{M}_{}\mathrm{pc}^2`$
The Gaussian shape of the inner part of the distribution depends mainly on the convolution of the true distribution with the beam (the scale length 649 pc in Eq. 3 corresponds to a FWHM of 50″). The right panel of Fig. 7 shows also the mass surface density of the HI gas. These data were obtained from Tilanus & Allen (1993) (hereafter TA). The column densities in TA were obtained with the VLA, and the authors noted that their estimated total HI mass within the Holmberg radius ($`0.87\times 10^9\mathrm{M}_{}`$) was lower than the corresponding masses estimated by other authors using single dish observations, e.g., both Rogstad et al. (1974) and Huchtmeier & Bohnenstengel (1981) obtained about $`1.6\times 10^9\mathrm{M}_{}`$ (all masses recalculated to the distance 4.5 Mpc). The most probable explanation for the deficiency was, according to TA, the absence of short-spacing information. In order to compensate for this we have added the “missing” mass surface density (2.4 M pc-2) to the HI data of TA. The resulting HI mass surface densities are still low and typically H2 gas constitute about 70% of the total mass of neutral hydrogen (H2 \+ HI) in the disk inside 6.8 kpc.

Based upon the H2 mass surface density distribution (Eq. 3) we estimate a total molecular gas mass within a radius of 7.3 kpc (the extent of the map) of $`3.5\times 10^9\mathrm{M}_{}`$, and extrapolated to the Holmberg radius (7$`\stackrel{}{.}`$3 or 9.6 kpc) the mass becomes $`3.9\times 10^9\mathrm{M}_{}`$ (which agrees well with the result of Crosthwaite et al. (2002), $`3.7\times 10^9\mathrm{M}_{}`$ (cf. also to the result $`3.3\times 10^9\mathrm{M}_{}`$ that was derived in Sect. 5.1 using the velocity-integrated intensity in the global spectrum). This is more than twice the total HI mass within the same galactocentric radius. It should, however, be noted that 80% of the HI mass in M 83 is found outside the Holmberg radius (Huchtmeier & Bohnenstengel, 1981). As a comparison, the stellar mass of M 83 is $``$ 37$`\times 10^9\mathrm{M}_{}`$, given V$`{}_{}{}^{0}{}_{T}{}^{}`$=7.37 (de Vaucouleurs et al., 1992) and assuming a mass-to-light ratio of 2.

In the Gaussian/exponential disk decomposition the major part of the molecular gas mass within the Holmberg radius lies in the disk, $`3.6\times 10^9\mathrm{M}_{}`$. Lord (1987) made CO($`J`$=1–0) observations in 21 positions in the central parts of M 83, and he estimated the total molecular gas mass within 115″ to be 1.2$`\times 10^9\mathrm{M}_{}`$ (we have compensated for the difference in $`X_{\mathrm{CO}}`$ and distance). Our H2 distribution gives 1.4$`\times 10^9\mathrm{M}_{}`$.

#### 5.2.2 Comparison with data at other wavelengths

In Fig. 8 we compare the CO radial distributions with corresponding photometry data in the B, K, and H$`\alpha `$ bands. The B and H$`\alpha `$ data come from Talbot et al. (1979), and the K data from Adamson et al. (1987). The scale lengths of the various tracers compare well over most of the radial range, but the comparatively low spatial resolution affects the CO curves in the central region and at the bar ends. The latter are responsible for the peaks at $``$120″. We have fitted a single exponential to the data in the radial range 140″–240″, which corresponds to the region outside the bar but inside the area where our CO survey is complete. The scale lengths are: 135$`\pm `$14″, 126$`\pm `$16″, 112$`\pm `$20″, 100$`\pm `$14″ and 132$`\pm `$35″ for CO($`J`$=1–0), CO($`J`$=2–1), K, B, and H$`\alpha `$, respectively (the errors of the latter three do not include the measurement errors). Thus, the disk scale lengths are, within the errors, the same for these star formation tracers.

The CO intensities decrease exponentially over the entire surveyed area, and we do not observe the sharp cut-off at 4$`\stackrel{}{.}`$5 seen by Crosthwaite et al. (2002). In fact, we find CO emission with $`I_{10}`$$``$ 1 K km s-1 in selected spectra as far from the nucleus as 5$`\stackrel{}{.}`$5.

### 5.3 Brightness distributions

We find that in the mem-deconvolved map the average CO($`J`$=1–0) intensity in the arms is about 2.5 times higher than in the interarm regions. Assuming a direct proportionality between CO emission and mass surface density, this agrees well with the lower limit of the arm–interarm mass surface density ratio, 2.3, found by RLH using purely kinematical arguments.

In Fig. 9 we show the mem-deconvolved CO($`J`$=1–0) velocity-integrated intensity map superposed on a true-color image. This image is a composite of three images in the filters B, V, and R, rendered as blue, green, and red, respectively. These images have been obtained at the ESO Danish 1.54 m telescope by Sören Larsen and the data are described in Larsen & Richtler (1999). On the large scale the CO($`J`$=1–0) emission follows the dust lanes, both in and between spiral arms, but at some locations the CO($`J`$=1–0) emission is displaced toward the regions where star formation takes place (Wiklind et al., 1990). Most notably, in the inner spiral arm on the eastern side, the CO($`J`$=1–0) ridge is located downstream of the dust lane with a separation of about 10–15″ (also shown in RLH), as are the HI and H$`\beta `$ emissions as shown by TA. We note that in general the OB associations straddle the rims of the CO($`J`$=1–0) emission concentrations, presumably due to extinction in combination with the photodissociation by the young stellar population. The outer and inner spiral arms appear to connect at several locations, and these connection regions often show an increased number of OB associations. The “Gould belt structure”, at (1$`\stackrel{}{.}`$75,–2$`\stackrel{}{.}`$25) in the coordinate system of Fig. 9, is such a region (Comerón, 2001).

Figure 10 shows the CO($`J`$=2–1) velocity-integrated intensity as contours on an RGB image, where the colors represent star formation (H$`\alpha `$; red), infrared emission (I band; green), and dust lanes (filters V–I; blue). In the disk, the observed CO($`J`$=2–1) brightness distribution strengthens the conclusions drawn from the CO($`J`$=1–0) emission. In the central region the higher angular resolution at the CO($`J`$=2–1) frequency reveals additional details. There are two central nuclear components relatively symmetrically placed on opposite sides of the IR center. At these locations the dust lanes, on the leading edges of the bar, attach to the nuclear ring seen in J–K by Elmegreen et al. (1998). Their color difference map (V–I) shows an area of high dust extinction at the location of the NE component, whereas at the SW component there appears to be relatively small amounts of dust. It appears as if the NE CO component lies above the disk (or the bulge) and the SW CO component lies below or behind the same. This was also concluded by Sofue & Wakamatsu (1994) based on similar data. Along the bar, the CO($`J`$=2–1) emission follows the dust lanes on the leading edges. In the disk the CO emission traces the H$`\alpha `$ emission very well, although there are regions with a clear anti-correlation, i.e., substantial CO emission but none, or very weak, H$`\alpha `$ emission, such as the peak at about half an arcminute west of the nucleus.

Figure 11 shows the CO($`J`$=1–0) velocity-integrated intensity on a gray-scale image of the HI column density (data kindly provided by R.P.J. Tilanus and R.J. Allen). The angular resolution of the HI map is slightly better than that of the CO map. The CO emission (which is assumed to be linearly proportional to the H2 mass surface density intensity) and the HI column density follow each other closely, with one outstanding exception: the nucleus, where very little HI is present. At some places in the disk the H2 and HI maxima are displaced from each other, but there is no apparent systematic trend that H2 is outside HI or vice verse. In general, the separations are of the order of 3-9$`\mathrm{}`$, which is small compared to the resolution in these maps (23$`\mathrm{}`$).

Local minima in the HI and H2 gas surface densities also correlate well. At (50$`\mathrm{}`$,180$`\mathrm{}`$) in the coordinate system of Fig. 11, marked minima can be seen in both H2 and HI. The total gas surface density in this area is about 3–6 M pc-2, of which only $``$1 M pc-2 is in the form of molecular gas. Such low mass surface densities should not give rise to any massive star formation, and indeed, the H$`\alpha `$ image shows very little emission here.

Figure 12 compares the H2 (as probed by the CO($`J`$=2–1) emission) and the HI gas in more detail, the resolution being $``$ 10$`\mathrm{}`$. In the main disk the CO($`J`$=2–1) emission and HI emission follow each other closely and several of the concentrations in the CO($`J`$=2–1) emission are also seen in the HI distribution. Again the exception is the nuclear region. The mass surface density of HI at the two nuclear components differ markedly: the SW component has a peak mass surface density of 6 M pc-2, while at the center of the NE component HI is actually seen in absorption against a nuclear continuum source, which again speaks in favor of the idea that the SW and NE components lie on different sides of the disk (or a bulge). The mass surface density of HI at the center is very much lower than that of H2, where the highest value in the $`I_{10}`$ map corresponds to $``$750 M pc-2. In the $`I_{10}`$ map it is not possible to distinguish the separate nuclei due to the relatively low spatial resolution and the coarse grid. However, using the average 2–1/1–0 line ratio (see Sect. 5.4), the mass surface density is $``$1100 M pc-2 for each of these components. As the physical properties of the molecular gas in the center may differ from that in the disk, application of the standard conversion factor may lead to an erroneous result. Indeed, it is likely that at least part of the emission in the center comes from a diffuse and non-virialized component. This would lead to an overestimate of the molecular gas mass (Papadopoulos & Allen, 2000).

The explanation for the nice correlation may be that the HI is mostly a dissociation product, as TA discuss. A smooth distribution of ‘primordial’ HI gas would not be detected in the interferometer observations (in which roughly 50% of the HI gas was detected) while HI produced as an effect of dissociation of molecular gas in the neighborhood of star formation would be localized, and therefore detected by the interferometer.

### 5.4 The CO ($`J`$=2–1)/($`J`$=1–0) line intensity ratio

We have created a map of the CO($`J`$=2–1)/($`J`$=1–0) line intensity ratio, $`R_{21}`$, at 49″ resolution. The mean value over this map is 0.72$`\pm `$0.19, which is within the errors the value derived from the velocity-integrated intensities in Table 4. Crosthwaite et al. (2002) estimate a ratio of 1.1$`\pm `$0.2, which we suspect is higher than ours because they have not corrected their CO($`J`$=2–1) data for the error beam contribution (considering the main beam efficiency of the NRAO 12 m telescope at this frequency, $``$0.56, the error beam is likely to be comparable to that of the SEST). In the nuclear region and along the bar the ratio is 0.83$`\pm `$0.04. The $`R_{21}`$ map shows also that on the arms, the line intensity ratio is lower than in the interarm regions. This is illustrated in Figure 13, where $`R_{21}`$ is plotted as a function of radius in a $`45\mathrm{°}`$ wide sector of the galaxy (position angle 90°–135°). In order to quantify and compare the ratio on the arms to that in the interarm regions, we have selected two sectors where the arms are well separated. The position angle of these sectors are: 90°–135° (SE) and 270°–315° (NW). In these sectors we selected spectra belonging to the inner, outer and inter-arm regions, and averaged them separately. Finally, “arm”-spectra were created by taking the average of the inner and outer arm spectra in the two respective sectors. From these eight spectra, four ratios were calculated. Assuming our error beam is correct (Sect. A), the line ratio is significantly lower on the arms: 0.59$`\pm `$0.04 (SE) and 0.69$`\pm `$0.04 (NW) compared to the ratio in the interarm region 0.88$`\pm `$0.17 (SE) and 0.98$`\pm `$0.21 (NW) (the errors are 3$`\sigma `$ and are estimated from the noise level in the spectra). We believe that this is an effect of a (partly) different cloud component in the interarm regions, where the heating of the gas is more efficient than in the dense arm cloud complexes. As a note: We have also calculated the $`R_{21}`$ value for the mem-deconvolved CO($`J`$=1–0) and CO($`J`$=2–1) maps at the common resolution 30″. The arm value is roughly the same (0.6$`\pm `$0.1), but the interarm value is higher, 1.3$`\pm `$0.3.

By averaging $`R_{21}`$ in concentric annuli, we find a trend in the line intensity ratio in the sense that it decreases with galactocentric distance, Fig. 14 (this is also indicated in Fig. 8). In particular, the lower envelope of the $`R_{21}`$-values is clearly decreasing. A similar difference between the $`R_{21}`$ line ratio in the disk and towards the nucleus has previously been shown to exist in our Galaxy (Sawada et al., 2001). The wiggles in Fig. 14 can be associated with features in the $`R_{21}`$ map. The nucleus appears to be surrounded by a ring of excited gas, which is reflected as a bump in the radial distribution. The bar ends and inner arms are reflected as a flat region at $``$120–150$`\mathrm{}`$ , while the low $`R_{21}`$-values in the outer arms suppresses the average $`R_{21}`$-ratio in the range 200–300″.

## 6 Conclusions

We have mapped, in over-sampling mode, the CO($`J`$=1–0) and CO($`J`$=2–1) brightness distributions over the entire optical disk of M83. Significant corrections for an error beam have been applied to the CO($`J`$=2–1) data. A mem-method has been used to increase the angular resolution of the data.

We find that:

i) the CO emission is strongly peaked toward the nucleus, which splits up into two components in the CO($`J`$=2–1) data. Also the bar ends are prominent. The CO emission follows the leading edges of the bar.

ii) Molecular gas spiral arms are clearly identified, and they trace, in most cases, the dust lanes. There are frequent “bridges” between spiral arms, and in some areas the arms are clearly disturbed.

iii) An average arm–interarm brightness contrast of about 2.5 is found for the CO($`J`$=1–0) line.

iv) The CO ($`J`$=2–1)/($`J`$=1–0) line ratio, which is about 0.77 on average, differ between the arm and interarm regions. It is significantly higher in the latter.

v) Regularly spaced molecular gas concentrations of mass $``$107 M lie along the arms.

vi) The estimated total molecular gas mass is 3.9$`\times `$109 M, and within 7.3 kpc the H2 mass dominates over that of HI by a factor of more than two. The estimated gas/stellar mass ratio is $``$0.1 in the optical disk.

vii) The CO and HI emissions are very well correlated in the optical disk, and frequently concentrations of the two gas components coincide.

viii) The CO radial brightness distribution in the disk follows that of other starformation tracers as H$`\alpha `$ emission, and continuum light in the B and K filters. The estimated scale length, of a fitted exponential, is about 120$`\mathrm{}`$, corresponding to 2.6 kpc at the adopted distance.

###### Acknowledgements.

We are very grateful toward the SEST staff for their support during observations and Swedish Natural Science Research Council for travel expenses support. We also wish to thank Remo Tilanus and Ron Allen for letting us use their HI data and Sören Larsen for the use of the optical images and the referee for valuable input.

## Appendix A Error beam

The beam pattern of a radio telescope is affected by its surface accuracy. This means that the power received in the error beam may very well be comparable to the power received in the main beam. For sources that are compact relative to the extent of the beam, it is relatively easy to compensate for this underestimate of the true flux: The measured antenna temperature ($`T_\mathrm{A}^{}`$) is simply divided by the main beam efficiency, $`\eta _{\mathrm{mb}}`$, to attain the main beam brightness temperature ($`T_{\mathrm{mb}}`$). The main beam efficiency is usually well known and can be found in the user manual of the different telescopes. For extended sources, the situation is more complex, since the antenna may receive power not only from the main beam. If the size of the error beam is comparable to the extent of the source, a simple multiplication of $`(\eta _{\mathrm{mb}})^1`$ will overestimate the flux at the position where the measurement is done. If the antenna pattern is well known, it is possible to calculate how the beam couples to the source and thus compensate for this effect. Unfortunately, the antenna pattern for SEST has never been measured properly, so this was not an option for us. Instead we devised a procedure to remove the impact of the error beam, using our own data set to estimate the strength and size of the error beam pattern.

Nearly all CO($`J`$=2–1) spectra did show two components: one narrow line, that agrees with the width and center velocity of the CO($`J`$=1–0) spectra at the same position and another, broader, component that is picked up by the error beam. By studying these broad wings in a number of spectra in the deconvolved CO($`J`$=2–1) data set and by comparing with spectra in data cubes convolved with different FWHPs, we estimated that the size of the error beam at 230 GHz is of the order of 200″. This is in good agreement with the error beam size of 3′, which was assumed by Johansson et al. (1998).

In order to try to estimate the impact of the error beam we produced a data set where we convolved the mem-deconvolved data set with a 200″ Gaussian profile. In each position an automated process, using the sliding window technique, calculated the velocity-integrated intensity in the region of the extended wings \[i.e. outside the velocity range where we have seen emission in the CO($`J`$=1–0) data\] in both the deconvolved cube and in the 200″ convolved cube. The ratio between these intensities corresponds to the fraction of the available emission which is picked up by the error beam. We will refer to this as the error beam efficiency. In principle, the error beam efficiency must lie in the interval 0.1–0.5, since the main beam efficiency is 0.5 and the moon beam efficiency is 0.9. Some values of our error beam efficiency estimates are higher than 0.5. This happens in regions where the emission in the main beam falls outside the expected velocity range (e.g., in regions with exceptionally strong streaming motions or wide emission profiles), leading to an over-estimation of the flux in the wings. An attempt to make the window wider had a negative impact on the results since we were left with fewer channels to calculate the velocity-integrated emission in the wings. Neglecting values above 0.5, we found that the error beam efficiency was on average 0.27$`\pm `$0.11.

Figure 15 shows the original spectrum, the 200″ convolved spectrum (scaled with the error beam efficiency), and the resulting (cleaned) spectrum after subtraction of the two previous spectra, in four typical positions. The convolved spectra reproduces the shape of the wings very well, and the shape of the cleaned spectra resembles the corresponding CO($`J`$=1–0) spectra (not shown).

We have corrected the convolved and the mem-deconvolved data sets for this error by subtracting the error beam contribution for each individual spectrum. Based on the previous result, we took a conservative approach and fixed the error beam contribution to 0.27. In principle, we could have let this ratio vary from spectrum to spectrum by fitting the 200″-convolved spectrum to the extended wings, but this would have introduced other, not easily controlled, errors.


text:
# Testing the locality of transport in self-gravitating accretion discs

## 1 Introduction

One of the basic unknowns in accretion disc theory is the physical mechanism ultimately responsible for angular momentum transport and energy dissipation in the disc. The usual way to overcome this difficulty is to assume that transport is dominated by some kind of viscous process, and to give an ad hoc prescription for the $`r\varphi `$ component of the stress tensor $`T`$ (usually vertically integrated, under the assumption that the disc is geometrically thin). The most widely used prescription is the so-called $`\alpha `$-prescription (Shakura and Sunyaev, 1973), according to which:

$$T_{r\varphi }=\left|\frac{\text{d}\mathrm{ln}\mathrm{\Omega }}{\text{d}\mathrm{ln}r}\right|\alpha \mathrm{\Sigma }c_s^2,$$
(1)
where $`\mathrm{\Omega }`$ is the angular velocity of the disc, $`c_s`$ is the thermal speed, $`\mathrm{\Sigma }`$ is the surface density, and $`\alpha `$ a free parameter, expected to be smaller than unity.

It has recently been recognized that accretion discs threaded by a weak magnetic field are subject to MHD instabilities (see Balbus and Hawley (1998) and references therein), that can induce turbulence in the disc, thereby being able to transport angular momentum outwards and to promote the accretion process. However, in many astrophysically interesting cases, such as the outer regions of protostellar discs, the ionization level is expected to be very low, reducing significantly the effects of magnetic fields in determining the dynamics of the disc. A possible alternative source of transport in cold discs is provided by gravitational instabilities.

The axisymmetric stability of a thin disc with respect to gravitational disturbances is determined by the condition (Toomre, 1964):

$$Q=\frac{c_s\kappa }{\pi G\mathrm{\Sigma }}>\overline{Q}1,$$
(2)
where $`\kappa `$ is the epicyclic frequency. The external regions of many observed systems are likely to be gravitationally unstable. If we consider, for example, the case of AGN accretion discs, it can be shown (Lodato and Bertin, 2003) that $`Q`$ falls below unity already at a distance of $`10^2`$ pc from the central black hole. In the case of protostellar discs, models of the outburst in FU Orionis systems (a class of young stellar objects that experience a phase of enhanced accretion) show that the disc is marginally stable already at a distance of $`1`$ AU from the central young star (Bell and Lin, 1994).

It has been suggested (Paczyński, 1978; Bertin and Lodato, 1999) that, in a cold enough disc, where condition (2) is violated, gravitational instabilities would provide a source of energy dissipation, therefore heating up the disc and leading to a self-regulated state, with the disc close to marginal stability. Gravitationally induced spiral structures are also able to transport angular momentum. In this paper we examine the issue of characterizing the angular momentum transport and energy dissipation in self-gravitating discs. In particular, we would like to test the extent to which the transport due to gravitational instabilities can be described as a local, viscous process, amenable to a description analogous to Eq. (1), as proposed by Lin and Pringle (1987).

Balbus and Papaloizou (1999) have argued that the long-range nature of the gravitational field would preclude a local description of transport, as implicitly assumed in a viscous scenario. In particular, they have shown how the energy flux contains extra terms, related to wave transport, which are “anomalous” with respect to a viscous flux. Gammie (2001) has performed local, shearing sheet, zero-thickness simulations of self-gravitating accretion discs, and has concluded that a local description is adequate in such “razor-thin” discs, extrapolating his results to discs for which $`H/r<0.1`$. However, Gammie’s simulations are not appropriate to test global effects, since locality is set up from the beginning, and they are only valid for infinitesimally thin discs. Rice et al. (2003a, b) have already shown using global, 3D simulations, how global effects can be important in the dynamics of self-gravitating discs, with respect to the related issue of fragmentation. Here we follow up the work of Rice et al., in order to quantify under which conditions a local model for transport in self-gravitating discs is justified.

## 2 Numerical setup

The simulations presented here were performed using smoothed particle hydrodynamics (SPH), a Lagrangian hydrodynamics code (see Monaghan, 1992). Our code uses a tree structure to calculate gravitational forces and the nearest neighbours of particles. The particles are advanced with individual time-steps (Bate et al., 1995), resulting in an enormous saving in computational times when a large range of dynamical timescales are involved.

We consider a system comprising a central star, modelled as a point mass $`M_{}`$, surrounded by a disc with mass $`M_{disc}`$. We have performed several simulation using different values of the parameter $`q=M_{disc}/M_{}`$. The initial surface density profile was taken to be a power-law $`\mathrm{\Sigma }r^1`$. The disc is initially gravitationally stable and the minimum value of $`Q`$ is $`Q_{min}=2`$. Our calculations are essentially scale-free. In dimensionless units, the disc extends from $`r_{in}=0.25`$ to $`r_{out}=25`$. We have used $`N=250000`$ particles.

We use an adiabatic equation of state, with adiabatic index $`\gamma =5/3`$. We explicitly solve the energy equation for our system, including heating from shocks and from artificial viscosity and introducing a cooling term, simply parameterized as:

$$\left(\frac{\text{d}u_i}{\text{d}t}\right)_{cool}=\frac{u_i}{t_{cool}},$$
(3)
where $`u_i`$ is the internal energy of a representative particle, and, as in Gammie (2001) and in Rice et al. (2003b), the cooling time-scale $`t_{cool}=\beta \mathrm{\Omega }^1`$, where $`\beta `$ is a free parameter, independent of radius. Gammie (2001) and Rice et al. (2003b) have shown that if cooling is too fast, then the disc will fragment as a result of the gravitational instabilities. In our simulations we have used $`\beta =7.5`$, in which case none of our discs was found to fragment.

## 3 Results

We have run our simulations for a few thermal timescales, in order to reach thermal equilibrium. In all cases the disc initially cools down until the value of $`Q`$ becomes small enough for gravitational instabilities to become effective and to provide a source of heating to the disc. At later stages the disc settles into a self-regulated state with $`Q1`$ over a large portion of the disc.

The $`r\varphi `$ component of the stress tensor (integrated in the $`z`$-direction) associated with self-gravity is given by (Lynden-Bell and Kalnajs, 1972):

$$T_{r\varphi }^{grav}=\text{d}z\frac{g_rg_\varphi }{4\pi G},$$
(4)
where $`𝐠`$ is the gravitational field of the disc. The hydrodynamic (Reynolds) contribution to the stress tensor is given by:

$$T_{r\varphi }^{Reyn}=\mathrm{\Sigma }\delta v_r\delta v_\varphi ,$$
(5)
where $`\delta 𝐯=𝐯𝐮`$ is the velocity fluctuation, $`𝐯`$ is the fluid velocity and $`𝐮`$ is the mean fluid velocity.

After averaging the stress tensor azimuthally and radially, over a small region $`\mathrm{\Delta }r=0.1r`$, we have computed the corresponding value of $`\alpha `$ (see Eq. (1)):

$$\alpha (r)=\left|\frac{\text{d}\mathrm{ln}\mathrm{\Omega }}{\text{d}\mathrm{ln}r}\right|^1\frac{<T_{r\varphi }^{grav}>+<T_{r\varphi }^{Reyn}>}{\mathrm{\Sigma }c_s^2}.$$
(6)
The resulting radial profiles of $`\alpha `$ are shown in Fig. 2 for the three cases $`q=0.05`$, $`q=0.1`$, and $`q=0.25`$. The upper panels show separately the hydrodynamic and gravitational contributions to $`\alpha `$, while the bottom panels show the sum of the two. The results are also averaged in time, over the last 500 timesteps of the simulations.

The simple parameterization of the cooling time adopted here allows us to have a direct insight on the problem of characterizing transport properties and dissipation in the disc, by using the following result (Pringle, 1981), which follows directly from the energy balance equation: if the disc is in thermal equilibrium, and if heating is dominated by a purely viscous process, described by Eq. (1), then the viscosity coefficient $`\alpha `$ and the cooling time-scale $`t_{cool}`$ satisfy the following relation:

$$\alpha =\left|\frac{\text{d}\mathrm{ln}\mathrm{\Omega }}{\text{d}\mathrm{ln}r}\right|^2\frac{1}{\gamma (\gamma 1)}\frac{1}{t_{cool}\mathrm{\Omega }}$$
(7)
In our simulations $`t_{cool}=7.5\mathrm{\Omega }^1`$, so from Eq. (7) we would expect that in thermal equilibrium, if energy dissipation can be described locally, $`\alpha `$ should be approximately constant with radius (this is indicated with the dotted line in the bottom panels of Fig. 2). We can interpret the bottom panels of Fig. 2 in the following way: the solid line gives a measure of the actual torque induced by gravitational instabilities in our simulations, while the dotted line is the torque that a viscous process should exert in order to dissipate the amount of energy which is actually dissipated in the simulations. Therefore, when the computed value of $`\alpha `$ is larger than the expected one, gravitational instabilities dissipate less energy than a viscous process, whereas in the opposite case, they dissipate more energy than a viscous process. Our results indicate that, up to a disc mass $`M_{disc}=0.25M_{}`$, gravitationally induced transport is fairly well described within a local framework. This result could also be anticipated, based on Fig. 1, that shows how the disc dynamics is dominated by high-$`m`$ modes, that dissipate on a short length-scale.

As a separate test for the locality of transport, we have also computed $`\alpha _{part}(r,d)`$, which is defined as the gravitational part of $`\alpha (r)`$, where the gravitational field $`𝐠`$ is computed taking into account only those particles which are inside a spherical radius $`d`$ from the radial point $`r`$ at which we are computing the stress. This quantity gives us a measure of the size of the region that mostly contributes to the stress at a given point.

Fig. 3 shows the results for $`r=15`$ for the three simulations. The upper panel shows clearly that the more massive the disc is, the more far regions contribute to the stress. However, it should be noted that, since $`Q1`$, a more massive disc is also hotter, and therefore thicker. The bottom panel of Fig. 3 shows $`\alpha _{part}`$ as a function of $`d/H`$, where $`H`$ is the disc thickness. In all cases, more than 80% of the contribution to the gravitational torque come from a region with size $`\mathrm{\Delta }10H`$. Transport is going to be local if $`\mathrm{\Delta }r`$. We can therefore conclude that, as long as $`M_{disc}<0.25M_{}`$, a disc with $`H/r0.1`$ will be reasonably described in terms of a local model. This is basically the global analogous of the results obtained by Gammie (2001) based on a local, 2D model.

Note that detailed vertical structure models of FU Orionis objects (Bell and Lin, 1994) lead to disc thickness $`H/r0.1`$, so that in this case we could expect global effects to be important in determining the emission properties of the disc, as suggested by Lodato and Bertin (2001).

## 4 Conclusions

We have performed global, three-dimensional SPH simulations of self-gravitating accretion discs. Our results show how the heating provided by gravitational instabilities can balance the cooling (that we have parameterized in a simplified way), leading to a self-regulated state where $`Q1`$. With a cooling rate $`t_{cool}=7.5\mathrm{\Omega }^1`$, none of our simulations was found to fragment. We have characterized the transport properties of the disc, by computing the stress tensor associated with gravitational instabilities and comparing the corresponding viscous dissipation rate to the actual heating rate in the disc. We have found that, if the disc is less massive than $`0.25M_{}`$, gravitational disturbances are dominated by high-$`m`$ modes, so that the transport properties of the disc can be described reasonably well in terms of a local, viscous formalism. In particular, we have found that the gravitational torque comes from a region of size $`10H`$, where $`H`$ is the disc thickness. We can therefore expect that discs with an aspect ratio $`H/r<0.1`$ can be treated locally.

The simulations presented here were performed at the UK Astrophysical Fluids Facility (UKAFF).