text:
# GMRT Observations of NGC 3079

## 1 Introduction

NGC 3079 is unusual amongst spiral galaxies in that it displays two well-defined radio lobes extending $``$ 30′′ (2.4 kpc, assuming a distance of 16.5 Mpc as listed in Table 1) from the nucleus (Duric et al. 1983), which are clearly out of the plane of this edge-on (84) system. Within the eastern radio lobe is a “bubble” of H$`\alpha `$ emission extending to 550 pc and displaying outflowing velocities up to 1000 km s-1, as shown by Hubble Space Telescope (HST) observations (Cecil et al. 2001 and references therein). Recent Chandra spacecraft images have also revealed X-ray emitting gas that is spatially correlated with the emission line gas (Cecil et al. 2002).

The high IR luminosity (Table 1) and bubble/cone morphology of the emission line gas has led to models of wind-like outflow related to a starburst. However, the galaxy also harbours an active galactic nucleus (AGN) and jet and about half of the IR luminosity could be due to recycled AGN (rather than hot stellar) photons (Iyomoto et al. 2001). The core/jet nature of the nuclear region is confirmed by VLBI detections (Irwin & Seaquist 1988; Trotter et al. 1998) as well as the recent measurement of an increase in separation between the compact components implying a jet outflow velocity of 0.13 $`\pm `$ 0.03 c (Kondratko et al. 2002). H2O maser emission in the disk indicates that a core mass of 106 M is present within a radius of 1 pc. Recent BeppoSax observations (Iyomoto et al. 2001) show that X-ray emission from the core is highly obscured and that after correction for absorption, the luminosity is 1042 to 1043 ergs s-1. In the inner parsec, this is 1% to 10% of the Eddington luminosity. NGC 3079 is therefore as active as the Seyfert galaxies and quasars. While it is still unclear to what extent the radio core powers the various components of the kpc-scale outflow, the radio lobes themselves appear to be a nearby analogue of the distant, powerful, extra-galactic radio sources, but within the dense, high angular momentum ISM of a nearby spiral.

Given the strength of the radio continuum emission in NGC 3079 and its peculiar nature, this galaxy represented a good target for low frequency observations using the Giant Metrewave Radio Telescope (GMRT) during its commissioning phase. We were particularly interested in searching for broad scale emission, especially as might be related to the nuclear outflow or to some of the other neutral hydrogen (HI) extensions and supershells seen in this galaxy (Irwin & Seaquist 1990). Multi-frequency observations were also important as a means to test for consistency and to enable the determination of minimum energy parameters such as magnetic field strength, cosmic ray energy density, and particle lifetimes. These are presented in Sec. 3. The basic data for NGC 3079 are presented in Table 1.

## 2 Observations and Data Reduction

The GMRT consists of 30 45-m antennas in an approximate ‘Y’ shape similar to the Very Large Array but each antenna cannot be moved from a fixed position. Twelve antennas are randomly placed within a central 1 km by 1 km square (the “Central Square”) and the remainder form the irregularly shaped Y (6 on each arm) over a total extent of about 25 km. Further details about the array can be found at the GMRT website at http://www.gmrt.ncra.tifr.res.in. We present details of the observations in Table 2. The observations were made in the standard fashion, with each source observation interspersed with observations of the phase calibrator. Data are collected in spectral-line mode at the GMRT over 128 channels in 2 Stokes (RR and LL) and one IF. The data are first written in a GMRT-based “lta” format and then converted into FITS format to facilitate processing in the standard fashion using the Astronomical Image Processing System (AIPS) of the National Radio Astronomy Observatory (NRAO).

Because no on-line flagging occurred in real time during GMRT observations, the data were first “pre-edited” using the program, GMRED, written by one of us (JAI) to run within AIPS. The purpose of GMRED is to search through each channel, each stokes and each IF separately and edit out bad data prior to normal editing and calibration. It combines clipping (flagging by a maximum and minimum) as well as median filtering to eliminate gross outliers. The maximum and minimum, the time interval for determination of the median, and the deviation allowed from the median before flagging can all be set by the user after inspection of the raw UV data. This allowed for removal of points such as might be obtained when an antenna is receiving data but is not yet pointed exactly on source, points that might be poor at the beginning of scans, and spikes from interference or electronics. GMRED is freely available with instructions for installation from http://www.astro.queensu.ca/$`\stackrel{~}{}`$irwin.

We then carried out a calibration in time in the standard fashion but using a single good channel of the phase calibrator rather than the inner 75% of the band. This was possible since the calibrator sources were quite strong at each frequency and avoided potential variations with frequency.

Since the data were acquired in spectral-line mode, we also carried out a bandpass calibration using the phase calibrator. The data from a single good channel were then Fourier Transformed and cleaned, first using a wide field to identify all sources and then of smaller multiple fields so that the outlying sources could be cleaned with a beam for the respective positions. All channels were then self-calibrated several times (and again edited, if required) using the single channel map. Cleaned maps were then made of all channels (i.e. a cube was created) and plots of the source flux density and map rms noise as a function of channel number were made. Loss of signal at the ends of the channels and occasional interference spikes or hardware problems required further editing of a number of channels at this stage. The final edited bandwidths and band centres are given in Table 2.

We then checked for positional offsets and variations in source flux density across the band. We found that the positions of NGC 3079 (central peak) and other point sources within the field agreed in position from one edge of the band to the other to within a single pixel (about 1/4 of a beam). The standard deviation in the flux density of NGC 3079 across the band (measuring all channels) in comparison to the flux density is given in Table 2. These errors can be compared to the bootstrapped flux density error on the phase calibrator as well as the estimated absolute error on the flux density calibrators (the latter, a few percent) to arrive at an estimate of the final calibration error listed in Table 2.

The UV data were then averaged over all good channels to make a single-channel UV data set and various maps were made and cleaned from these data, using a variety of UV weighting functions and tapers. Several maps experienced small improvements from further minor editing and self-calibration at this stage. We then corrected for the primary beam at the various frequencies for all maps. However, after these corrections, the increase in source flux density was much less than the absolute calibration error and the point-by-point change in flux density per beam was less than the rms map noise for all maps. Thus we present the final maps without this correction.

We then checked and made minor corrections to the map registration by comparing the positions of the core of NGC 3079 and point sources around it (via gaussian fits) with those of the NRAO VLA Sky Survey (NVSS) and Faint Images of the Radio Sky at Twenty-centimeters (FIRST) maps (see Figs. 4b and 4c, respectively). We estimate the final positional accuracy to be within 1 arcsec. We also compared our 326 MHz flux densities with those measured from the Westerbork Northern Sky Survey (WENSS, Rengelink et al. 1997) map (Fig. 4a). After smoothing to the same resolution, the WENSS map has a peak and integrated flux density that are 10% lower than our GMRT flux densities. Given the 10% accuracy of our GMRT 326 MHz flux density and the 5% accuracy of the WENSS flux density, these values are in agreement.

## 3 Results

### 3.1 Total Intensity Maps

The total-intensity maps, shown over a variety of spatial scales, can be seen in Figs. 1 (326 MHz), 2 (615 MHz), and 3 (1281 MHz). Their map properties are listed in Table 2. In these figures, we show positive contours to low intensity levels (1.5$`\sigma `$) but caution in the interpretation of the lowest levels and summarize what we consider to be real features in Sect. 3.1.4. Total flux densities from these maps are given in Table 3 and are plotted on the spectrum shown in Fig. 5. The new GMRT flux densities are consistent with the known spectrum at these frequencies. (The one discordant low point at 365 MHz is from the Texas Survey using data taken between the years 1974 and 1983; see Douglas et al. 1996.) The new data are clearly and firmly in the non-thermal part of the spectrum.

#### 3.1.1 326 MHz

Our low-resolution GMRT image (Fig. 1a) shows strong emission from NGC 3079, but no emission from either of its two companion galaxies, namely NGC 3073 and MCG-9-17-9, the upper limit to their flux densities being 9 mJy and 5 mJy, respectively.

At 326 MHz (Fig. 1) an envelope of emission is seen around the galaxy which is more extensive than previously observed (cf. Duric & Seaquist 1988, Duric et al. 1983, de Bruyn 1977, and Fig. 4) with the possible exception of the low-resolution image of Condon (1987). The emission envelope follows the disk of the galaxy and is likely due to distributed sources in the disk similar to other star-forming spirals.

To determine the vertical scale height of this thick disk, we created a model galaxy with the geometry of NGC 3079, whose projected intensity at any point is proportional to the line of sight distance through the galaxy at that point. The model galaxy was then convolved with a beam equivalent to that of Fig. 1a. We then matched the extent of the model major axis with that of NGC 3079 at the 2$`\sigma `$ level and varied the vertical height until a match was achieved in the perpendicular direction (the western extension was omitted, see below). The best fit was found for a vertical height of 60′′ (4.8 kpc). Since the intrinsic brightness distribution should fall off with height from mid-plane, the vertical scale height will be larger than 4.8 kpc.

The 326 MHz emission also shows much substructure with improving resolution. The most obvious feature is an extension towards the west which ends in a knot at RA = 10 01 38.1, DEC = 55 41 14. A feature is also seen at a similar position in the WENSS image (Fig. 4a) at the same frequency though with somewhat different structure, given the different uv coverage and spatial resolution. This feature resolves into two distinct extensions (Fig. 1c) which roughly flank the western nuclear radio lobe. These two extensions have counterparts in HI, best seen in the moment map of Irwin & Seaquist (1987).

A large loop of emission is visible at the north end of the galaxy disk extending towards the east of it. The north edge of this loop makes the major axis on the north side appear to bend around to the east. The central minimum is at RA = 10 01 58.3, DEC = 55 43 01 and the loop diameter is 47′′ (3.8 kpc) parallel to the major axis and extends to 90′′ (7.2 kpc) from the major axis. This loop is also coincident with an HI (open-topped) loop seen at a velocity of 1011 km s-1 in Irwin & Seaquist (1991, their Fig. 1).

Another interesting result is the curvature of the emission towards the east away from the optical disk at the southern tip of the major axis. This behaviour again mimics that seen in the neutral hydrogen distribution (Irwin & Seaquist 1991). The NVSS image (Fig. 4b) also hints at an extension in this direction. Since both the radio continuum and HI emission are disconnected from the optical disk, the curvature is not tidal but is likely due to ram pressure affecting only the ISM components. This indicates that there may be more extensive intergalactic gas in this region through which NGC 3079 has some local motion. Such motion would be (in projection) towards the west or south-west or could be rotation (tumbling) with the southern major axis moving towards the west. At the present time, the only extended X-ray emission observed appears to be associated with the nuclear outflow (Pietsch et al. 1998, Cecil et al. 2002). Thus deeper X-ray observations are warranted.

#### 3.1.2 615 MHz

The 615 MHz data (Fig. 2) reveal the largest radio continuum loop which has yet been associated with this galaxy. The loop extends to the east, has a central minimum located at RA = 10 02 07.4, DEC = 55 42 37, has a diameter of 64′′ (5.1 kpc) parallel to the major axis and extends to 2.2 (11 kpc) from the major axis. This is much farther than the minimum vertical scale height measured from the 326 MHz data and confirms that a very extensive radio continuum halo exists around NGC 3079. The loop is not centered over the eastern nuclear radio lobe, but rather offset from it to the north. However, there is much substructure in the feature, especially on its southern edge which is closer to the position of the radio lobe. It is possible that several features are blended along the line of sight in this edge-on system. Alternatively, the feature may be related to the nuclear outflow but source precession, motion of the galaxy with respect to an IGM, or flows into previously emitted gas may have affected its position over time. This feature is also hinted at in the WENSS image (Fig. 4a).

The same two western extensions on the west side of the disk (as noted from the 326 MHz data) are seen prominently at 615 MHz as well. These (Fig. 2b) reach $``$ 2 (9.6 kpc) from the plane.

The distinctive cross shape and markedly stronger emission of the nuclear region in this galaxy is especially visible in Fig. 2c.

#### 3.1.3 1281 MHz

At the highest frequency (Fig. 3), we again see emission related to the northern loop at RA $``$ 10 01 55, DEC $``$ 55 43 18, a feature related to the giant eastern loop at RA $``$ 10 02 03, DEC $``$ 55 42 00, and the two western extensions. At these higher resolutions, we also detect details of the nuclear radio lobes. Fig. 3b, for example, can be compared to the 1420 MHz FIRST image of Fig. 4c which is quite comparable in structure, except that our 1281 MHz image has a slightly lower resolution and shows some additional emission away from the plane. (Note that the FIRST image contains a noticeable stripe.) Fig. 3c reveals the lobes’ substructure in exquisite detail, especially the loop-like nature of the eastern lobe (cf. Duric & Seaquist 1988). The H$`\alpha `$ emission occurs roughly in the region of the gap between the eastern loop and the nucleus (Cecil et al. 2001).

#### 3.1.4 Summary of Radio Continuum Structure

Each data set set has its own limitations, via uv coverage, spatial resolution, and possible residual errors which could not be eliminated. We will therefore consider a feature to be real only if it is observed at more than one frequency (note that all data were taken on different dates, Table 2) or if it is observed in another independent observation using a different telescope. Thus we summarize the structures that we believe to be real as follows:

1. The western extension (Fig. 1a) which resolves into two distinct extensions (Fig. 1c, Fig. 2a, 2b, Fig. 3a, 3b, Fig. 4a, and is also seen in HI, Irwin & Seaquist 1987). These features extend to $``$ 10 kpc, roughly flank the western nuclear radio lobe and may be associated with it.

2. The northern loop (Fig. 1b, c, Fig. 3a and HI, Irwin & Seaquist 1991) of 3.8 kpc diameter and extending 7.2 kpc from the major axis. This loop appears to be associated with the disk.

3. The giant eastern radio loop (Fig. 2a, 2b, Fig. 3a, Fig. 4a) of 5.1 kpc diameter and extending 11 kpc from the major axis. It is offset from the eastern nuclear radio lobe but could still be related to larger scale and/or earlier nuclear outflow.

4. The offset, curved emission on the south tip of the major axis (Fig. 1a, 1b, 1c, Fig. 4b and HI, Irwin & Seaquist 1991), suggesting that there may be motion through broader scale IGM gas.

### 3.2 Spectral Index Maps

To make the spectral index maps, we took Fig. 1c, Fig. 2b, and Fig. 3a, which are all at similar spatial resolutions, smoothed them to a common resolution, and blanked all points below a level of 5 sigma. Spectral index maps ($`S_\nu \nu ^\alpha `$) were then formed and are shown in Fig. 6. Fig. 6a shows the 326-615 MHz spectral index and Fig. 6b shows the spectral index map between 615 and 1281 MHz. To estimate the error in the spectral index maps, at each freuency we extracted a map of noise from a region in which there was no emission. We then added this noise map to the map of the emission region of NGC 3079. This resulted in 3 maps of emission with added noise typical of that data set. New spectral index maps were then formed from the noise-added images. We then formed a map of the difference between the noise-added spectral index maps and the original spectral index maps. The rms of these difference maps is typical of the error on any individual point in Figs. 6a and 6b, i.e. $``$ 0.1.

In Table 3, we list values of the spectral index, as determined from the total fluxes listed in Table 3 (the Global spectral index) and as determined from the mean of the spectral index maps over equivalent regions (the Mean spectral index). Note that the mean spectral index, which applies to high S/N regions, shows very little change with changing frequency and the global spectral index is consistent with these values between 615 and 1281 MHz ($`\alpha `$ = $``$0.92 to $``$1.02). However, the 326$``$615 MHz global value is flatter at $``$0.75. This suggests that there may be more faint extended emission at 326 MHz which remained undetected. There may indeed be some trend, however, for stronger changes in spectral index in the high latitude gas than in the disk. For example, extrapolating the 615 MHz emission in the broad extended halo region using $`\alpha `$ = $``$0.98 to 326 MHz should have revealed as extended a halo in Fig. 1c above the noise. A trend towards a flatter spectral index at lower frequencies can also be observed in the global spectrum (Fig. 5) when other independent data points are considered. Flattening of the spectral index at lower frequencies in the non-thermal domain has been observed before in other galaxies. Israel & Mahoney (1990) attribute such flattening to thermal absorption but this is unlikely to play a major role in our data since the mean spectral index, which applies mainly to the galaxy disk where most of the thermal gas should occur, shows only very minor flattening. Energy losses at the higher frequencies seems more likely (Hummel 1991).

In this regard, it is interesting to examine the distribution of spectral index as shown in Fig. 6a and Fig. 6b for differences. We note that there are differences between the maps that are greater than a typical error on any point. For example, a map of the difference between Fig. 6a and 6b (not shown) has an rms of 0.49 whereas a typical error bar is about 0.2. However, closer examination shows that the larger differences occur near the edges of the maps (where the S/N is lower) or in regions in which there are strong intensity gradients (where 1′′ errors in registration can make larger changes in the measured value of $`\alpha `$). Therefore, in the region shown in Fig. 6c, we see no convincing evidence for changes in the spectral index with frequency above the error bars. Therefore, we averaged Fig. 6a and 6b to create a map of the average spectral index, shown in Fig. 6c.

Fig. 6c has a mean of $`\alpha `$ = $``$0.97 and an rms of 0.18 with a typical error less than 0.1. There are hints that the spectral index may flatten at the bases of extensions (e.g. on the west side). The region centered on the nucleus clearly has a flatter spectral index than in other regions of the disk, e.g. a mean of $``$0.73 in the central $``$ 30′′ diameter region, in agreement with previous observations of this region (Duric & Seaquist 1983). The flatter values near the nucleus are consistent with the fact that this is a region in which cosmic ray acceleration is occurring and particles are younger.

### 3.3 Computation of Minimum Energy Parameters

From the mean spectral index map (Fig. 6c) and one total intensity map (the 615 MHz map that was used as to create the images of Fig. 6), we can compute the following minimum energy parameters (see Pacholczyk 1970, Duric 1991): the cosmic ray energy density, $`u_{CR}`$, the cosmic ray electron diffusion length, $`L_D`$, magnetic field strength, $`B`$, and the particle lifetime, $`t`$. These quantities are computed for each pixel and represent means along a line of sight. A line of sight distance, and therefore a geometry, is required for these calculations. We also require assumptions regarding the lower and upper frequency cutoffs of the spectrum, $`\nu _1`$ and $`\nu _2`$, and the ratio of the heavy particle energy to the electron energy, $`k`$. (See Irwin et al. 1999 for a previous example.) Values which have been used are specified via several models which are listed in Table 4. In Fig. 7, we show the results for Model 1 which uses an inclined disk, $`k=\mathrm{\hspace{0.17em}40}`$, $`\nu _1=\mathrm{\hspace{0.17em}10}^7`$Hz and $`\nu _2=\mathrm{\hspace{0.17em}10}^{11}`$Hz. Table 4 lists the means calculated over these maps as well as those of the other models.

The results for Models 1 to 4 (Table 4) show that changes in geometry (inclined disk or oblate spheroid) and changes in the upper frequency cutoff make relatively small differences in the results. The main differences are due to the assumption of heavy particle to electron energy (cf. Model 1, Model 2) which introduce variations in the minimum energy parameters of order a factor of $``$ 1.5. Thus, we expect Model 1 (Fig. 7) to be a good representation of the data to within a factor of $``$ 2. Note that this choice affects only the absolute scale of the maps of Fig. 7 and not the point to point variation.

The fact that the measured flux densities for NGC 3079 fall within the non-thermal regime of the spectrum (Fig. 5) indicates that the thermal fraction should be very small. However, point to point variation in thermal fraction could change the appearance of the minimum energy parameter maps. Condon (1992) has provided an estimate (thought to be correct to within a factor of $``$ 2) of the thermal fraction in spiral galaxies, as measured globally, viz. $`S/S_T\mathrm{\hspace{0.17em}1}+\mathrm{\hspace{0.17em}10}\left(\nu /GHz\right)^{(0.1|\alpha _{NT}|)}`$, where $`\alpha _{NT}`$ is the non-thermal spectral index, $`S`$ is the total flux density and $`S_T`$ is the flux density due to the thermal component. Taking $`\alpha _{NT}`$ to be $``$0.97 (from the mean spectral index map of Fig. 6c), we find that only 3.6%, 6%, and 11% of the total emission could be due to thermal components at 326, 615, and 1281 MHz, respectively. To see what such a contribution might do to the calculation of the minimum energy parameters, we subtracted the above fractions from the total intensity images at the 3 frequencies, recomputed the spectral index map and recomputed the minimum energy maps. These results are also listed in Table 4 (Model 5) and show that there is a negligible difference in the minimum energy parameters as a result of possible thermal contributions.

In Fig. 7, we also show a map of power, $`P=U/t`$, where $`U`$ is the cosmic ray energy density, $`u_{CR}`$ (Fig. 7a), integrated along a line of sight and $`t`$ is the particle lifetime (Fig. 7d). The result is identically, $`P=(1+k)L`$, where $`L`$ is the observed luminosity at a point. Thus the map of $`P`$ closely resembles the map of total flux density but does not match it exactly because the computation of $`L`$ requires an integration over frequency which is dependent on spectral index, and the spectral index is different at different points in the map. The map of $`P`$ represents the rate at which cosmic rays must be accelerated in order to maintain equilibrium. Integrated over the map, the total power is $`P_T=\mathrm{\hspace{0.17em}9.2}\times \mathrm{\hspace{0.17em}10}^{41}`$ ergs s-1. The total energy in cosmic rays is $`U_T=\mathrm{\hspace{0.17em}1.1}\times \mathrm{\hspace{0.17em}10}^{57}`$ ergs. These values are minima since not all of the galaxy’s emission is represented in the regions shown in Fig. 7. In comparison, the values for the Galaxy are $`P_T\mathrm{\hspace{0.17em}2}\times \mathrm{\hspace{0.17em}10}^{40}`$ ergs s-1 and $`U_T\mathrm{\hspace{0.17em}10}^{55}`$ ergs (Condon 1992). Thus the cosmic ray energy generation rate and total energy are two orders of magnitude higher in NGC 3079 than in the Galaxy. The spectral power of Table 3 suggests a supernova rate of 0.27 yr-1 (see Condon 1992) if all emission has this origin. This is clearly an upper limit since some of the emission is due to the AGN and associated outflow.

Although the total cosmic ray energy and power are higher than the Galaxy, they are not atypical of spiral galaxies (Duric 1991). Our average magnetic field strength of $``$ 10 $`\mu `$Gauss, for example, agrees with the 12 $`\mu `$Gauss found for NGC 6946 (Ehle & Beck 1993). In addition, the distribution of values is as expected. The magnetic field strength, energy density, and power are all higher near the nucleus where there is much activity and correspondingly, the lifetime of the particles and the diffusion length are shorter. We note, in addition, that the peak power of $`8\times \mathrm{\hspace{0.17em}10}^{38}`$ ergs s-1 is considerably smaller than the power of the X-ray core in this galaxy ($`10^{4243}`$ ergs s-1, Sec. 1).

## 4 Discussion

NGC 3079 is a unique nearby laboratory for studying the presence of radio jets in a spiral galaxy and the effect that these might have on the global characteristics of the galaxy. Although little else has been seen in the way of such well-defined radio lobes in a spiral, recent observations of Abell 428 have also revealed a classical double-lobed radio source associated with a disk galaxy (Ledlow et al. 1998, Ledlow et al. 2001), suggesting that NGC 3079 is no longer unique. Our new GMRT observations have revealed not only the details of the nuclear radio structure (e.g. Fig. 3c) at high (3′′) resolution, but broad scale radio continuum emission (e.g. Fig. 1a) at low (55′′) resolution and much substructure.

It appears that much of the extended emission around NGC 3079 is related to the disk. This is seen via the envelope of emission whose morphology follows the disk (Fig. 1a) and also from at least one feature which forms a kpc-scale loop to the east of the northern tip of the galaxy (Fig. 1b) which has an HI counterpart. Thus, the radio continuum features are similar to those found in NGC 5775 in which extraplanar features and supershells are detected in every component of the ISM (Lee et al. 2001) or NGC 2613 in which a large open-topped radio continuum arc with some associated HI spans much of the galaxy (Chaves & Irwin 2001). Neither of the latter two galaxies are known to harbour an AGN.

The nuclear radio emission and radio lobes have been well-delineated by our observations. In addition, large radio extensions and loops are seen which may be associated with it. This includes the two western extensions and possibly the giant eastern loop best seen at 615 MHz which is somewhat offset from a direct extension of the eastern nuclear lobe. Fig. 8 shows the western extensions and eastern loop in relation to the high resolution nuclear radio lobes. Indeed the placement of these features and their similar sizes are suggestive of a nuclear origin. If the two large features on either side of the galaxy are indeed associated with the nuclear outflow then they indicate regions of previous outflow and suggest that the nuclear activity in NGC 3079 is episodic. If we extrapolate a particle lifetime of $`10^8`$ yrs (Fig. 7d) and use the outer distance of the giant eastern loop (11 kpc, Sect. 3.1.2) then the required average bulk speed of the particles from the disk to this distance is $``$ 100 km s-1.

The presence of offset radio continuum and HI emission at the southern tip of the major axis argues for some motion with respect to an IGM and this may also be affecting the location of the eastern loop. However, other effects may also be important. For example, the fact that the VLBI jet direction is offset from the direction of the nuclear radio lobes argues either that the source outflow direction is changing with time or that the nuclear jets are bending through the galaxy’s ISM and being focussed along the minor axis until the outflow moves from an ISM-dominated region to an IGM dominated region. We speculate that the emission around NGC 3079 may be similar to what has been observed at low frequencies around M 87 in which the galaxy is embedded in a complex IGM in its immediate vicinity which is due to its own nuclear outflow (Owen et al. 2000).

## 5 Conclusions

We have observed the double-lobed spiral galaxy, NGC 3079 using the GMRT at three frequencies and detected structure over resolutions from 3′′ to 56′′. At high resolution, the structure of the nuclear radio emission is shown in exquisite detail (Fig. 3c). In addition, we have detected a number of new intermediate and broad scale features, namely:

1. Two western extensions above and below the western nuclear radio lobe which reach $``$ 10 kpc from the galaxy’s disk. These features have associated HI emission.

2. A giant eastern radio loop, best seen in Fig. 2a or Fig. 8, which is slightly offset from a line extending through the nucleus and the eastern radio lobe. The giant loop extends to 11 kpc from the disk.

3. A northern loop, best seen in Fig. 1b. There is a large HI arc at this position. It is 3.8 kpc in diameter and extends 7.2 kpc from the major axis. This loop appears to have originated in the disk.

4. A broad halo of emission with a vertical scale height of at least 4.8 kpc, which follows the disk of the galaxy.

5. Emission at the southern tip of the major axis that is offset from the optical major axis suggesting that there may be motion through broader scale IGM gas.

No emission has been detected from either companion galaxy.

We suggest that Features 1 and 2 above may be associated with earlier nuclear outflow in which case the average outflow velocity would be $``$ 100 km s-1. In addition, the IGM that is implied in 5, above, could be a local IGM, formed via outflow from NGC 3079 itself.

We have also computed spectral index maps over a restricted spatial region and have computed maps of minimum energy parameters for this galaxy, including the cosmic ray energy density, the diffusion length, the magnetic field strength, the particle lifetime (each of these representing averages along a line of sight) and finally the power. The total energy content is several orders of magnitude larger than that of our Milky Way but the parameters, in general, are typical of spiral galaxies.

## Acknowledgments

JAI wishes to thank the scientists and staff of NCRA, Pune, for graciously allowing her to work at this institute during a sabbatical year. The GMRT is a national facility operated by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. This research has made use of the NASA/IPAC extragalactic database (NED) which is operated by the Jet Propulsion Laboratory, Caltech, under contract with the National Aeronautics and Space Administration.


format: text

text:
# Optimal State Discrimination Using Particle Statistics

## Abstract

We present an application of particle statistics to the problem of optimal ambiguous discrimination of quantum states. The states to be discriminated are encoded in the internal degrees of freedom of identical particles, and we use the bunching and antibunching of the external degrees of freedom to discriminate between various internal states. We show that we can achieve the optimal single-shot discrimination probability using only the effects of particle statistics. We discuss interesting applications of our method to detecting entanglement and purifying mixed states. Our scheme can easily be implemented with the current technology.

Recently, there has been an emerging interest in the use of particle statistics (both bosonic and fermionic) for quantum information processing . In fact, it was shown that a useful task such as entanglement concentration could be accomplished, even if non-optimally, using only the effects of quantum statistics, without the need for any other interactions . The above investigations differ significantly from some previous suggestions, where either anyonic statistics  or the effects of electronic statistics in conjunction with other interactions  were used for quantum information processing. Schemes using only particle statistics  would be very useful for tasks implemented with identical particles that interact very weakly or not at all with each other, such as photons or neutrons. This weak interaction can be beneficial for information processing as it may reduce the unwanted coupling to the environment. Such schemes are also extremely general in the sense of being independent of the actual particle species. It is however not known whether such schemes can accomplish quantum information processing efficiently. Here we present a particular quantum information processing task involving two qubits and show that it can be performed optimally using only quantum statistics. Moreover, we point out how the task of discriminating quantum states can be applied to detecting entanglement and purifying mixed states. We also discuss how to generalize these tasks to $`N`$ qubits and argue that quantum statistics could be used to perform even this generalized task optimally. While the two qubits and other small $`N`$ qubits versions of our protocol can be tested with photons, electrons, neutrons or atoms, the large number of qubits versions could have interesting implementations in optical lattices .

One of the striking aspects of quantum mechanics is that it is not possible to perfectly discriminate between two states unless they are orthogonal. Suppose someone prepares two qubits encoded in the internal degrees of freedom of two identical particles — say, in the spin of two electrons or the polarization of two photons — in one of the following two possible states:

* spins aligned (parallel) and pointing in an arbitrary direction;
* spins anti-aligned and pointing in an arbitrary direction.

We will assume that the two states are equally likely for the sake of simplicity, but all the presented results will be valid for any *a priori* distribution. Note also that, unless otherwise stated, whenever we mention spin we will actually refer to any two-dimensional internal degree of freedom, be it for fermions or for bosons. As we do not have any knowledge of the direction of alignment in both of the above cases, the overall states are mixed. They are described respectively by the following density operators:

$$\rho _2=\frac{1}{4\pi }𝑑\mathrm{\Omega }|\mathrm{\Omega }\mathrm{\Omega }||\mathrm{\Omega }\mathrm{\Omega }|,$$
(1)
and

$$\sigma _2=\frac{1}{4\pi }𝑑\mathrm{\Omega }|\mathrm{\Omega }\mathrm{\Omega }||\mathrm{\Omega }^{}\mathrm{\Omega }^{}|.$$
(2)
Here, the subscripts indicate that we are considering two particles. The ket $`|\mathrm{\Omega }`$ represent the spin-up state along the axis defined by the angle $`\mathrm{\Omega }`$, while $`|\mathrm{\Omega }^{}`$ is the orthogonal state spin-down. Thus, $`\rho _2`$ represents an equal mixture of spins aligned along an arbitrary axis in space, while $`\sigma _2`$ is the equal mixture of anti-aligned spins of two spin-$`\frac{1}{2}`$ particles. It is impossible to discriminate between these states perfectly because they are not orthogonal.

Optimal results are known for the discrimination of any two given quantum states $`\eta `$ and $`\eta ^{}`$ . The maximal probability of ambiguously discriminating between two *a priori* equally likely quantum states in a single-shot measurement is given by the Helstrom formula:

$$P_H(\eta ,\eta ^{})=\frac{1}{2}+\frac{1}{4}Tr|\eta \eta ^{}|.$$
(3)
We now present a procedure for discriminating between $`\rho _2`$ and $`\sigma _2`$, both for fermions and bosons, based only on the effects of particle statistics. To use these effects we interfere particles at a beam splitter. Here we use *beam splitter* in a generic sense, referring not only to the common optical element (partially silvered mirror) used with photons, but also to any device presenting an analog behavior for other kinds of particles, as already suggested for electrons . In order to distinguish the states we rely on path measurements that discriminate between bunching and antibunching, and that are performed on our particles after letting them pass simultaneously through a $`50/50`$ beam splitter. Note that in such balanced beam splitters two indistinguishable particles will always bunch if they are bosons , and always antibunch if they are fermions  (see also ).

For fermions, our guess in the case of the antibunching result is that the input state was $`\rho _2`$, while in the case of bunching is that it was $`\sigma _2`$. The probability of success of our procedure is then:

$$P_{BS}(\rho _2,\sigma _2)=\frac{3}{4}.$$
(4)
This probability can easily be calculated by noticing that the only case for which our guess could be incorrect is when we have the antibunching result. Then, we conclude that the input state was $`\rho _2`$, while it could have actually been $`\sigma _2`$. On the other hand, in the case of the bunching result we know for sure that the input state was $`\sigma _2`$, since — according to the Pauli exclusion principle — two particles with aligned spins cannot end up in the same output arm of the beam splitter. When our input state is $`\sigma _2`$, the antibunching happens with probability $`1/2`$, giving in total a probability of incorrect inference of $`1/4`$. In the case of bosons, our protocol is exactly the opposite of the fermionic one, but yields precisely the same efficiency: this time the antibunching results stand for $`\sigma _2`$, whereas the bunching ones stand for $`\rho _2`$, but the probability of success coincides with Eq. (4).

Interestingly, the Helstrom formula gives the same result for the maximal probability of discriminating between these two states:

$$P_H(\rho _2,\sigma _2)=\frac{3}{4}.$$
(5)
We can thus conclude that our procedure is optimal for both fermions and bosons.

Helstrom’s probability can still be achieved in other cases using the effects of quantum statistics. We now introduce a case of special interest that can be applied to other tasks such as entanglement detection and state purification, as will be shown later in the paper. Suppose that we have to discriminate between the following two states:

* spins aligned (parallel) and pointing in an arbitrary direction (the same as in Eq. (1));
* each spin in the maximally mixed state.

The latter state is represented by the following operator:

$`\tau _2={\displaystyle \frac{1}{4}}I^2={\displaystyle \frac{1}{4}}\left(|00|+|11|\right)\left(|00|+|11|\right),`$ (6)
where $`|0`$ and $`|1`$ are any two orthogonal spin states.

The strategy now is exactly the same as before. It relies on the fact that if the state is $`\rho _2`$ the particles can, due to their indistinguishability, give only one result (antibunching in the case of fermions or bunching in the case of bosons), while if the input state is $`\tau _2`$ both results are possible. This time, the probability of success is:

$$P_{BS}(\rho _2,\tau _2)=\frac{5}{8},$$
(7)
and this coincides with the Helstrom result

$$P_H(\rho _2,\tau _2)=\frac{5}{8},$$
(8)
so our procedure is optimal in this case too. The calculations leading to the above results are analogous to the ones in the previous discrimination case.

One of the most interesting applications of our approach can be found in the detection of entanglement. We illustrate this in the case of pure states of two particles. Suppose, for example, that we are to discriminate any maximally entangled state of two qubits from any disentangled (product) state. As before, the qubits are supported by the internal degrees of freedom of two identical particles (labelled, say, $`A`$ and $`B`$). Suppose in addition that we are given two identical copies of the state. In order to detect entanglement, we take the same particle (either $`A`$ or $`B`$) from each pair and interfere them at a beam splitter, as shown in Fig. 1. The crux of the argument is that if the state is entangled, then the reduced states of these particles will be maximally mixed as in the state $`\tau _2`$. On the other hand, if the state is separable, then the two interfering particles are in the state $`\rho _2`$. This is the same as in our discrimination procedure above, and so there is a probability of $`5/8`$ to detect entanglement. This example can be generalized to other entangled pure states, and, more interestingly, to more particles. We note that there is a close analogy between this method and our entanglement concentration scheme in . In particular, if the two states were less than maximally entangled, then by detecting entanglement we would actually also amplify it (see  for more details).

Another interesting application is in mixed state purification, as in . Suppose that we start with two qubits, each in some mixed state. We would like to make the state of these qubits purer (in the sense of having lower linear entropy), but also to preserve their original direction in the Bloch sphere. The optimal way of doing so (as proven in ) is to project the joint state onto the symmetric subspace, in which case the resulting mixed state is purer and yet preserves the original direction. If the projection is unsuccessful, the qubits are thrown away. This is exactly the same as our probabilistic discrimination with a beam splitter.

We would now like to investigate the generalization of the above results to $`N`$ particles. For this, we use a generalized (N-port) balanced beam splitter, as shown in Fig. 2, which acts only on the spatial degrees of freedom of the input particles. This action is given by a unitary matrix $`U_N`$, with elements:

$$u_{mn}=\frac{1}{\sqrt{N}}e^{i\frac{2\pi }{N}(m1)(n1)}.$$
(9)
(Note that there exist alternative descriptions of balanced multiport beam splitters .) The square of the norm of each element in the matrix represents the probability that the particle in the $`m`$-th input arm of the beam splitter ends up in the $`n`$-th output arm. Since all these elements have norm $`\frac{1}{\sqrt{N}}`$, we have a representation of a balanced N-port beam splitter.

The aim is now to discriminate between the N-particle generalizations of $`\rho _2`$ and $`\tau _2`$. Those states are given by:

$$\rho _N=\frac{1}{4\pi }𝑑\mathrm{\Omega }\left(|\mathrm{\Omega }\mathrm{\Omega }|\right)^N,$$
(10)
and

$$\tau _N=\frac{1}{2^N}\left(|00|+|11|\right)^N.$$
(11)
To calculate the Helstrom probability, one has to diagonalize the matrix $`\rho _N\tau _N`$. This turns out to be straightforward once we notice that $`\rho _N`$ can also be represented as an equal mixture of all possible symmetric states within the basis $`\{|S_i:i=1,\mathrm{},N+1=|S|\}`$ of $`N`$ qubits. Then, we have:

$$\rho _N=\frac{1}{|S|}\underset{i=1}{\overset{|S|}{}}|S_iS_i|.$$
(12)
It is now easy to calculate the Helstrom formula by expanding $`\rho _N`$ in a basis consisting of the union of a basis of the symmetric sub-space and a basis of its orthogonal complement. The result is:

$$P_H(\rho _N,\tau _N)=1\frac{(N+1)}{2^{(N+1)}}.$$
(13)
Alternatively, we can calculate the average probability of success to distinguish states $`\rho _N`$ and $`\tau _N`$ using the following expression:

$$P_H(\rho _N,\tau _N)=\frac{1}{2}1+\frac{1}{2}p.$$
(14)
Here, the $`1/2`$ factors refer to the fact that the two states are prepared with equal probability. The term $`1`$ comes from the fact that the state $`\rho _N`$, supported on the symmetric subspace (of dimension $`d_S=N+1`$), is always identified reliably as such, and $`p`$ is the probability of identifying the other state, $`\tau _N`$. Since $`\tau _N`$ is maximally mixed, it is uniformly distributed over the whole space of $`N`$ qubits (of dimension $`d=2^N`$). In this case $`p=\frac{1d_S}{d}`$, which after substitution in Eq. (14) gives Eq. (13) right away.

The problem to apply our discrimination scheme to $`N`$ particles is that it becomes exponentially hard to calculate $`P_{BS}(\rho _N,\tau _N)`$ as $`N`$ increases. Moreover, it is not clear which inference strategy should be followed. For fermions, the natural generalization seems to be to associate the antibunching results with $`\rho _N`$ and the others with $`\tau _N`$. For bosons, on the other hand, a more subtle strategy may be needed. This is because, loosely speaking, for bosons there is no clear analogue of the Pauli exclusion principle. Furthermore, even without having the complete calculations for $`N>2`$, we would like to emphasize the remarkable fact that the Helstrom probability $`P_H(\rho _N,\tau _N)`$ is equal to the probability of success of a fermionic beam splitter strategy described above, if calculated under the assumption that the particles are classical (i.e., always distinguishable by some arbitrary label), but obey a constraint equivalent to the Pauli exclusion principle (not allowing more than two particles in the same internal state to share the same output arm of the beam splitter). The overall probability is then calculated by summing up the probabilities of all the possible outcomes rather than the amplitudes, as it would be done in the quantum case. However, in the case of three fermions ($`N=3`$), we have performed the full quantum calculations (i.e., taking properly into account the effects of statistics) for a three-port balanced beam splitter and obtained:

$$P_{BS}(\rho _3,\tau _3)=\frac{3}{4},$$
(15)
which is equal to $`P_H(\rho _3,\tau _3)`$. We believe this result of obtaining the optimal discrimination probability using the effects of particle statistics (in multiports) can be generalized to an arbitrary $`N`$, both for fermions and bosons, and we continue research in this direction. For now this remains a conjecture. In an optical lattice, with one particle in each lattice site, a multiport beam splitter could probably be simulated by dissolving N potential wells and then creating a new set of N wells . Of course, if the particles interact, then the effective beam splitter will be modified, and here we only point out the plausibility of creating multiport beamsplitters (or multiparticle interference) in an optical lattice.

In this paper we have shown that it is possible to perform an optimal quantum information processing task using only the effects of particle statistics. In particular, we have presented a strategy for discriminating between two non-orthogonal states of two qubits (encoded in the internal degrees of freedom of identical particles) using beam splitters. We have considered two discrimination scenarios and in each of them our strategy differs (symmetrically) between fermions and bosons, but offers the same efficiency. We also pointed out how our discrimination scheme can be applied to detecting entanglement and purifying mixed states. In addition, we have calculated the Helstrom probability for $`N`$ qubits in one of our discrimination scenarios. We have shown by explicit calculation that this probability can be achieved in a fermionic three-port beam splitter strategy and that it is the same as the fermionic strategy for general $`N`$ if the fermions are considered as classical particles that obey the Pauli exclusion principle as the only additional constraint. An advantage of our method is that it can also be easily implemented with the current technology.

Our work suggests a number of interesting research directions. One problem is to prove the optimality of the beam splitter strategy in the case of $`N`$ qubits and its application to multiparty entanglement detection. This, we hope, will answer the question of weather the symmetry between fermions and bosons in our strategy will be preserved for a generalized beam splitter. It may also lead to a simple and physically intuitive selection principle governing bosonic behavior. Another possible direction is to classify all the pairs of states that can be optimally discriminated with our scheme. Finally, our results suggest that it would be worth to further explore the role of particle statistics in quantum information tasks, in particular in efficient quantum computation.

After the conclusion of this paper, an interesting related work was pointed out to us . The authors would like to thank G. Castagnoli on many useful discussions on the link between computational complexity and particle statistics. His ideas have stimulated much of the work in this paper. We also would like to thank F. A. Bovino and P. Varisco for useful discussions on the role of particle statistics in quantum information processing. S.B. acknowledges support from the NSF under Grant Number EIA-00860368. Y.O. acknowledges support from Fundação para a Ciência e a Tecnologia from Portugal and wishes to thank the Institute for Quantum Information at Caltech for their hospitality. N.P. and V.V. thank Elsag S.p.A. for financial support and Centro de Física de Plasmas at IST in Lisbon for their hospitality. V.V. acknowledges support from Hewlett-Packard company, EPSRC and the European Union projects EQUIP and TOPQIP. This work has been partly supported by the QUIPROCONE collaboration grant number 044.


text:
# Systematic variation of the residual of PSR B1937+21, the fingerprint of its companion star?

## 1 Introduction

PSR B1937+21, the first millisecond pulsar and also the fastest rotating pulsar up to date, was discovered by Backer $`\&`$ Kulkarni et al. (1982). Which has been measured as a singular pulsar. The systemic long-term variation in the residual of time of arrivals (TOAs) has been measured by different authors, Kaspi et al. (1994), Lommen $`\&`$ Backer (2001). The variation may be caused by planet around the pulsar. However, such possibility was constrained stringently by Thorsett $`\&`$ Phillips (1992).

This paper provides an alternative model to interpret the systematic variation in TOAs, the orbital precession induced time delay of a binary pulsar system. Three numerical solutions are obtained through fitting the residual of TOAs, all of which demands small companion star mass, $`10^2M_{}`$, and in turn small projected semi-major axis, $`x10^4`$ or $`10^3`$. Therefore, the corresponding radial velocity of the pulsar to the line of sight (LOS) is so small that it cannot be resolved from the timing noise of TOAs. Which might explain why the companion star has not been observed from timing measurement.

This paper shows that although the companion star of PSR B1937+21 cannot be measured directly from the velocity curve, however, the orbital parameters could still be extracted through the long-time variation of TOAs.

By the orbital precession model, when the fitted orbital parameters are put into the standard timing model, the behavior of PSR B1937+21 as binary pulsar should be similar to PSR J2051-0827 (Doroshenko et al. 2001) and PSR B1957+21 (Arzoumanian et al. 1994), which have similar companion masses and orbital period. Moreover, when the additional time delay is also added to the standard timing model, the systemic variation of TOAs of PSR B1937+21 should be eliminated, so are PSR J2051-0827 and PSR B1957+21. PSR B1937+21 might be another millisecond pulsar whose companion star is evaporated by the radiation of the pulsar as PSR B1957+21.

## 2 The orbital precession model

Barker $`\&`$ O’Connell (1975) gives the first gravitational two-body equation with spins, in which the orbital angular momentum, vector, $`𝐋`$, is expressed as precesses around a vector combined by $`𝐒_1`$ and $`𝐒_2`$, and the magnitude is small 2 Post Newtonian (PN), typically $`1^{}10^4`$yr-1, which is insignificant in pulsar timing measurement.

On the other hand, the precession of orbital angular momentum vector, $`𝐋`$, around the total angular momentum vector, $`𝐉`$, in the special case that ignores the spin angular momentum of one star in a binary pulsar system, has been studied by many authors (Apostolatos et al. 1994, Kidder 1995, Wex $`\&`$ Kopeikin 1999). In this case the precession velocity of the orbit plane is significant, 1.5 PN, typically $`1^{}`$yr-1 for binary pulsars with orbital period around 10 hours. Which is measurable in pulsar timing.

Therefore, the two kinds of expressions on the precession velocity of $`𝐋`$ seems contrary each other. To investigate their relationship, the first and most important question to be answered is: which vector should $`𝐋`$ precess around ? The two kinds of expressions both agree that $`𝐒_1`$ and $`𝐒_2`$, as well as the vector combined by them precess rapidly (1.5PN) around the total angular momentum, $`𝐉`$. Since $`𝐉`$ is at rest relative to LOS (after counting out the proper motion), so the small precession velocity of $`𝐋`$ around the vector combined by $`𝐒_1`$ and $`𝐒_2`$ doesn’t mean that precession velocity of $`𝐋`$ around $`𝐉`$ is also small, and in turn relative velocity LOS is small.

The situation is analogy to the following case. If the binary system we observe is replaced by a solar system, one cannot say that the velocity of planet A is ignorable because its velocity with respect to planet B is very small. Instead only when the velocity of planet A is very small relative to the Baryon center of the solar system which is at rest to the observer is very small, can one conclude that this velocity is ignorable.

Obviously the role of $`𝐉`$ in a binary system is equivalent that baryon center in the solar system. And only the velocity relative to $`𝐉`$ make sense in observation. Therefore, the precession of $`𝐋`$ should be expressed as around $`𝐉`$ to compare with observation.

In the one spin case the precession velocity of $`𝐋`$ must be around $`𝐉`$ and the magnitude must be significant 1.5PN (Apostolatos et al. 1994, Kidder 1995, Wex $`\&`$ Kopeikin 1999). Gong develop the one spin to the general two spins case, and obtains that still the precession velocity of $`𝐋`$ must be around $`𝐉`$ and must be significant 1.5PN (Gong1 2003).

The difference is that in the two spins case, the magnitude of the precession velocity of $`𝛀_0`$ varies rapidly due to the variation of the sum the spin angular momenta of the two stars, $`𝐒`$, which can lead to significant secular variabilities in binary pulsars. Whereas, the one spin case predicts a constant magnitude of $`𝛀_0`$, and thus it cannot explain the significant secular variabilities measured in binary pulsars.

Since the spin angular momenta of the two stars in binary pulsar, $`S_1`$ and $`S_2`$, are much smaller than that of the orbital angular momentum, $`L`$, therefore, the sum of the two spins, $`𝐒=𝐒_1+𝐒_2`$, is also much smaller than $`L`$, or, $`S/L1`$ must be satisfied for a general binary pulsar system. Moreover since $`𝐉=𝐋+𝐒`$ ($`𝐉`$, $`𝐋`$ and $`𝐒`$ forms a triangle), the constraint, $`S/L1`$, means the misalignment angle between $`𝐉`$ and $`𝐋`$ must be very small, $`\lambda _{LJ}1`$. The conservation of the total angular momentum of a binary pulsar system, $`\dot{𝐉}=0`$, leads to (Barker $`\&`$ O’Connell 1975),

$$𝛀_0\times 𝐋=𝛀_1\times 𝐒_1𝛀_2\times 𝐒_2,$$
(1)
the magnitude of the left hand side of Eq(1) is $`\mathrm{\Omega }_0L\mathrm{sin}\lambda _{LJ}`$. From which one can see that for a given torque caused by $`S_1`$ and $`S_2`$ at the right hand side of Eq(1), there are two choices for $`\mathrm{\Omega }_0`$. The first one is that $`\mathrm{\Omega }_0`$ is small in the case that $`\lambda _{LJ}1`$, the second one is that $`\mathrm{\Omega }_0`$ is large in the case that $`\lambda _{LJ}1`$. As analyzed above, not matter from the point of view of observation, or from the constraint that must be satisfied by a binary pulsar ($`S/L1`$), $`\mathrm{\Omega }_0`$ have to follow the second choice. Which means $`𝐋`$ must precess around $`𝐉`$ with a very small angle of precession cone, $`\lambda _{LJ}1`$, and the magnitude of $`𝛀_0`$ must be significant.

By Eq(1), the velocity of the orbital precession, $`𝛀_0`$, which is along $`𝐉`$ and including two spins can be obtained (Gong1 2003).

$$\mathrm{\Omega }_0=\rho \mathrm{\Omega }_2\mathrm{sin}\lambda _{LS}+\rho (\mathrm{\Omega }_1\mathrm{\Omega }_2)\frac{S_1^{}}{S}\mathrm{sin}\lambda _{LS_1}^{},$$
(2)
where $`\rho =\frac{1}{\mathrm{sin}\lambda _{JS}}`$, and $`S_1^{}=S_1\mathrm{cos}\eta _{SS1}`$ denotes the component of $`𝐒_1`$ in the plane determined by $`𝐒`$ and $`𝐉`$. $`\lambda _{LS_1}^{}`$ is the angle between $`S_1^{}`$ and $`𝐋`$. Note that $`L\mathrm{sin}\lambda _{LJ}\rho S`$ is used in Eq(2), since $`S/L1`$. The right-hand side of Eq(2) can as well be written by replacing subscribes 1 with 2 and 2 with 1. The precession velocity of the two stars, $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$, are given by Barker $`\&`$ O’Connell’s equation.

Actually Eq(2) can be regarded as transforming the direction of $`𝛀_0`$ of the Barker $`\&`$ O’Connell’s equation (two spins) to $`𝐉`$; and it can also be regarded as add one more spin into the former one spin case given by Apostolatos et al. (1994), Kidder (1995), Wex $`\&`$ Kopeikin (1999).

The scenario of the motion of a binary system has been discussed by authors, Smarr $`\&`$ Blandford (1976), Hamilton $`\&`$ Sarazin (1982). In which $`𝐋`$, $`𝐒_1`$ and $`𝐒_2`$ all precess around $`𝐉`$ rapidly (1.5PN). And the relative velocities of $`𝐋`$ to $`𝐒_1`$ and $`𝐒_2`$ are very small. From which a slow precession of $`𝐋`$ around $`𝐉`$ is impossible because of the triangle constraint, $`𝐉=𝐋+𝐒`$. The scenario can solve all the puzzle concerning the motion of the three vectors, $`𝐋`$, $`𝐒_1`$ and $`𝐒_2`$.

Since $`𝐒_1`$ and $`𝐒_2`$ precess with different velocities, $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ respectively ($`m_1m_2`$), then $`𝐒`$ varies in both magnitude and direction ($`𝐒_1`$, $`𝐒_2`$ and $`𝐒`$ form a triangle), then from the triangle of $`𝐒`$, $`𝐋`$ and $`𝐉`$, in react to the variation of $`𝐒`$, $`𝐋`$ must vary in direction ($`|𝐋|=`$const), which means the variation of $`\lambda _{LJ}`$ ($`𝐉`$ is invariable).

The change of $`\lambda _{LJ}`$ means that the orbital plane tilts back and forth. In turn, both $`\lambda _{LS}`$ and $`\lambda _{JS}`$ vary with time. Therefore, from Eq(2), the derivative of the rate of orbital precession can be given by,

$$\dot{\mathrm{\Omega }}_0=\rho \dot{\mathrm{\Omega }}+\dot{\rho }\mathrm{\Omega },$$
(3)
where $`\mathrm{\Omega }=\mathrm{\Omega }_0/\rho `$ and $`\dot{\mathrm{\Omega }}`$ is given by

$$\dot{\mathrm{\Omega }}=\mathrm{\Omega }_2\mathrm{\Omega }_{12}X_3X_4\mathrm{\Omega }_{12}X_1(\mathrm{\Omega }_{01}X_2+\mathrm{\Omega }_{12}X_3)+\mathrm{\Omega }_{12}X_1X_5,$$
(4)
where $`\mathrm{\Omega }_{12}=\mathrm{\Omega }_1\mathrm{\Omega }_2`$, $`\mathrm{\Omega }_{01}=\mathrm{\Omega }_1\mathrm{\Omega }_0`$, $`X_1=\frac{S_1^{}}{S}\mathrm{sin}\lambda _{LS_1}^{}`$, $`X_2=\mathrm{tan}\eta _{ss1}`$, $`X_3=\frac{S_{V1}S_{V2}}{S^2}\frac{\mathrm{sin}\eta _{s1s2}}{\alpha \mathrm{sin}\lambda _{JS}}`$, $`X_4=\frac{\mathrm{cos}^2\lambda _{LS}}{\mathrm{sin}\lambda _{LS}}`$, $`X_5=\mathrm{cot}\lambda _{LS_1}^{}\dot{\lambda }_{LS_1}^{}`$ and $`\dot{\rho }=\frac{X_3X_4\mathrm{\Omega }_{12}}{\mathrm{sin}^2\lambda _{JS}}`$. In which $`\alpha =\mathrm{sin}\lambda _{JS}+\frac{\mathrm{cos}^2\lambda _{LS}}{\mathrm{sin}\lambda _{LS}}`$, $`S_{V1}=S_1\mathrm{sin}\lambda _{JS1}`$ and $`S_{V2}=S_2\mathrm{sin}\lambda _{JS2}`$ represent components of $`𝐒_\mathrm{𝟏}`$ and $`𝐒_\mathrm{𝟐}`$ that are vertical to $`𝐉`$.

Note that $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ are unchanged when ignoring the orbital decay (2.5PN), and $`\lambda _{LS_\alpha }`$ are unchanged also, since they decay much slower than the orbital decay (Apostolatos et al. 1994). Through Eq(3), $`\ddot{\mathrm{\Omega }}_0`$ can be easily obtained (Gong1 2003) .

The S-L coupling in a binary pulsar system results the following effects. (1) the apsidal motion of the orbital plane which can be absorbed by $`\dot{\omega }`$ (Smarr $`\&`$ Blandford 1976)

$$\dot{\omega }^{obs}=\dot{\omega }^{GR}+\mathrm{\Omega }_0\mathrm{cos}\lambda _{LJ}\dot{\omega }^{GR}+\mathrm{\Omega }_0.$$
(5)
Notice that $`\mathrm{\Omega }_0`$ is a function of time due to the variation of $`S`$, as shown in Eq(2), whereas, $`\dot{\omega }^{GR}`$, the GR prediction of the advance of periastron is a constant.

(2) the precession of the orbital plane which can be absorbed by $`\dot{x}`$,

$$\dot{x}=x\mathrm{\Omega }_0\mathrm{sin}\lambda _{LJ}\mathrm{sin}\eta _0\mathrm{cot}i.$$
(6)
$`|\dot{x}|`$ of Eq(6) is usually much larger than $`|\dot{x}^{GR}|`$, which is caused by the gravitational wave induced orbital decay.

(3) the nutation of the orbital plane which can be absorbed by $`\dot{P}_b`$. The variation in the precession velocity of the orbit results in a variation of orbital frequency ($`\nu _b=2\pi /P_b`$), $`\nu _b^{}\nu _b=\dot{\mathrm{\Omega }}_0\mathrm{\Delta }t`$ . Then $`\dot{\nu }_b=\dot{\mathrm{\Omega }}_0`$, which lead to the variation of $`P_b`$ (Gong1 2003),

$$\dot{P}_b=\frac{\dot{\mathrm{\Omega }}_0P_b^2}{2\pi }.$$
(7)
Notice that the coupling of the spin induced quadrupole moment with the orbit (Q-L coupling) can also cause apsidal motion which can explain the secular variabilities on $`\dot{\omega }`$ and $`\dot{x}`$ measured in binary pulsars. However, such effect cannot explain $`\dot{P}_b`$, and the second and third order of derivatives of parameters like, $`\ddot{x}`$, $`\ddot{P}_b`$ and $`P_b^{(3)}`$. Since the Q-L coupling is similar to that of the S-L coupling in the one spin case, which corresponds to a static precession or apsidal motion of the orbital plane. In other words, the velocity, $`\mathrm{\Omega }_0`$ is a constant in these two cases. Whereas, in the two spins S-L coupling, $`\mathrm{\Omega }_0`$ is a function of time, as given by Eq(2), which can not only explain the parameters that can be explained by the Q-L coupling or the one spin S-L coupling, but also explain parameters that cannot be explained by them.

## 3 Application to PSR B1937+21

The essential transformation relating solar system barycentric time $`t_b`$ to pulsar proper time $`T`$ is given by the equation (Manchester and Taylor 1977)

$$t_bt_0=T+\mathrm{\Delta }_R+\mathrm{\Delta }_E+\mathrm{\Delta }_S+\mathrm{\Delta }_A,$$
(8)
where $`\mathrm{\Delta }_R`$ is the ”Roemer time delay”, which corresponds to the propagation time across the binary orbit; $`\mathrm{\Delta }_E`$ and $`\mathrm{\Delta }_S`$ are the orbital Einstein and Shapiro delays; and $`\mathrm{\Delta }_A`$ is a time delay related with aberration caused by rotation of the pulsar. The dominant term concerning the orbital effects is the Roemer delay, $`\mathrm{\Delta }_R`$, given by (Taylor $`\&`$ Weisberg 1989)

$$\mathrm{\Delta }_R=xF(\omega +u),$$
(9)
where

$$F(\omega +u)=\mathrm{sin}\omega [\mathrm{cos}ue(1+\delta _r)]+[1e^2(1\delta _\theta )^2]^{1/2}\mathrm{cos}\omega \mathrm{sin}u,$$
(10)
where $`u`$ is the eccentric anomaly, and $`\omega `$ is the angular distance periastron from the node. In calculation, the small quantities $`\delta _r`$ and $`\delta _\theta `$ due to aberration are ignored.

When the contribution of orbital precession to the time of arrival is considered, the dominant term that absorbs the additional time delay is also the Roemer delay. Which can deviate from the standard one, $`\mathrm{\Delta }_R`$, and therefore, leads to an additional time delay,

$$\mathrm{\Delta }=xF(\omega ^{}+u^{})xF(\omega +u)=2x\mathrm{sin}(\mathrm{\Omega }_0t/2)\mathrm{cos}\theta x\mathrm{\Omega }_0t\mathrm{cos}\theta .$$
(11)
where $`\theta =(\dot{\omega }+\mathrm{\Omega }_0)t/2+\omega _i`$. The additional time delay, $`\mathrm{\Delta }`$, given by Eq(11), is a function of time due to the variation of $`\mathrm{\Omega }_0`$. Which can contribute the TOAs via Eq(8). And such dynamic effect might be interpreted as effect of propagation, like $`\dot{DM}`$ in some binary pulsars. The effect of $`\dot{\mathrm{\Omega }}_0`$, $`\ddot{\mathrm{\Omega }}_0`$ can be absorbed by such parameters as $`\dot{P}_b`$, $`\ddot{P}_b`$, $`\ddot{x}`$, $`x^{(3)}`$, $`\ddot{DM}`$, $`DM^{(3)}`$, etc (Gong2, Gong3, 2003). The orbital precession model predicts that $`\frac{P_b^{(3)}}{\ddot{P}_b}\frac{\ddot{P}_b}{\dot{P}_b}\mathrm{\Omega }_2`$ for a same binary pulsar, so is the derivatives of $`x`$ and $`DM`$. Evidences of such correlation can be found in binary pulsars, such as PSR J2051-0827 (Doroshenko et al. 2001) and PSR B1957+21 (Arzoumanian et al. 1994), PSR B1534+12 (Konacki et al. 2003). Further more the orbital precession model predicts that binary pulsars with long orbital period, i.e., for days, correspond to much smaller $`\mathrm{\Omega }_0`$ and $`\mathrm{\Omega }_2`$, and in turn much smaller derivatives of $`x`$, $`P_b`$ and $`DM`$ by Eq(6), Eq(7) and Eq(11). Which make the corresponding derivatives of such binary pulsar more difficult to observe than that of binary pulsars with small $`P_b`$, i.e., a few hours. Such correlation can also be find in binary pulsars

As shown in Fig 1, the amplitude of the measured time delay is about $`10^6`$s. By simple estimation, if $`x\mathrm{\Omega }_0\mathrm{cos}\theta 10^{13}`$ in Eq(11), then the orbital precession induced time delay can be $`\mathrm{\Delta }10^6`$s also. There are two possibilities to satisfy this relation: (1) small $`\mathrm{\Omega }_0`$ and $`\dot{\omega }`$ like PSR B1855+09 (large $`P_b`$), and large $`x`$, however in such case, it will take much longer time to finish a variation like that of Fig 1. And (2) large $`\mathrm{\Omega }_0`$ and $`\dot{\omega }`$, i.e., $`10^{}`$yr-1, and small $`x`$, obviously case (2) can not only satisfy the amplitude but also the periodicity of variation. Therefore, the orbital precession model implies that PSR B1937+21 should have small $`x`$, and small $`P_b`$ (which corresponds to large $`\mathrm{\Omega }_0`$ and $`\dot{\omega }`$).

The vectors $`𝐒_1`$, $`𝐒_2`$ and $`𝐒`$ are studied in the coordinate system of the total angular momentum, in which the z-axis directs to $`𝐉`$, and the x- and y-axes are in the invariance plane. $`𝐒`$ can be represented by $`S_P`$ and $`S_V`$, the components parallel and vertical to the z-axis, respectively:

$$S=(S_V+S_P)^{1/2}.$$
(12)
$`S_P`$ and $`S_V`$ can be expressed (recall $`S_{V1}`$, $`S_{V2}`$ and $`S_V`$ form a triangle) as

$$S_P=S_1\mathrm{cos}\lambda _{JS1}+S_2\mathrm{cos}\lambda _{JS2},$$

$$S_V=(S_{V1}^2+S_{V2}^22S_{V1}S_{V2}\mathrm{cos}\eta _{S1S2})^{1/2},$$
(13)
where $`\eta _{S1S2}`$ is the misalignment angle between $`S_{V1}`$ and $`S_{V2}`$, which can be written as

$$\eta _{S1S2}=(\mathrm{\Omega }_1\mathrm{\Omega }_2)t+\varphi _i.$$
(14)
Therefore, by the variation of $`S`$ as function of time (in the case of one spin, $`S=`$const), one can obtain $`\mathrm{\Omega }_0`$ as function of time through Eq(2). Thus the measured $`\mathrm{\Delta }^{obs}`$ of Fig 1 can be fitted step by step via the orbital precession induced $`\mathrm{\Delta }`$ given by Eq(11),

$$\mathrm{\Delta }^{obs}(t_{k+1})\mathrm{\Delta }^{obs}(t_k)=\mathrm{\Delta }(t_{k+1})\mathrm{\Delta }(t_k).$$
(15)
Note that both $`\mathrm{\Omega }_0`$ and $`\mathrm{cos}\theta `$ are functions of time in Eq(11). Which are responsible for variation of residual in Fig 1.

Three numerical solutions can be obtained by the Monte Carlo solution of Eq(15). The solutions are shown in Table 1,3,5, respectively. The three solutions shows that $`P_b`$ is from $`0.0708`$d to $`0.10`$d, and $`x`$ from $`3.7310^4`$s to $`4.2510^3`$s, and $`m_2`$ from $`0.0126M_{}`$ to $`0.0373M_{}`$. The possibility of other solutions cannot be excluded.

The standard DD timing model uses 12 parameter concerning the orbit to fit the delays of Eq(8) (Taylor $`\&`$ Weisberg 1989). The total number of binary parameters of the orbital precession model as shown in Table 1 is also 12. There are some parameters in common, such as, $`e`$, $`P_b`$, $`\omega _0`$ and $`T_0`$. The difference is that the latter can cause an additional time delay, as shown in Eq(11), through the parameters of the second row of Table 1 (including $`m_1`$ and $`m_2`$ of the first row). Whereas, the former doesn’t have this delay. In other words, the latter only adds one more delay term to the right hand side of Eq(8), and all other terms (effects) are the same as the former models. The relation of the orbital precession model with the DDGR model can be given by $`\dot{P}_b^{obs}=\dot{P}_b^{GR}+\dot{P}_b^{OP}`$, and by Eq(5).

By the results of Table 1,3,5 and the orbital precession model, one can easily obtain the following Post-Keplerian parameters. Notice that Table 2,4,6 correspond to the solution 1,2,3, respectively.

## 4 discussion

The orbital precession model predicts that PSR B1937+21 is a binary pulsar with very small orbital period, $`P_b=1.7`$hr. Then there is a question, why it has been measured as a singular pulsar? Actually it can be answered by the fitted orbital parameters of Table 1. The radial velocity (to line of sight) of the pulsar is given (Shapiro $`\&`$ Teukolsky 1983),

$$𝐯_1𝐧=\frac{2\pi a_1\mathrm{sin}i}{P_b(1e^2)^{1/2}}[\mathrm{cos}(\omega +u)+e\mathrm{cos}\omega ],$$
(16)
where $`a_1=\frac{m_2a}{m_1+m_2}`$. Then the observed pulsar period deviates from the intrinsic period by a small amount,

$$P^{obs}=P^{int}(1+\frac{𝐯_1𝐧}{c})P^{int}(1+\frac{2\pi x}{P_b}).$$
(17)
Therefore, the contribution of the orbital precession effect to $`P^{obs}`$ is about $`\delta P^{OP}2\pi xP/P_b1.610^{20}`$s (corresponds to solution 1). Whereas, the last decimal of pulsar period, $`P`$, is $`4(2)10^{18}`$s (Kaspi et al. 1994), which is approximately two order of magnitude larger than $`\delta P^{OP}`$. In other words, the modulation of the radial velocity of the pulsar (induced by binary motion) to the TOAs is ignorable. Moreover the small orbital inclination angle ($`i10^2`$) given by solution 1 and 2, means that the orbital plane of the binary may be face on, then it is also impossible to observe the eclipse of this binary pulsar. These factor may make it difficult to observe the effects of the companion star.

By the orbital precession model, a binary pulsar with very small orbital period, i.e., a few hours, can produce a long-term time delay i.e., years, as given by Eq(11). Therefore, the companion star of PSR B1937+21, although unmeasurable in the orbital period time scale, due to small $`x`$, might imprint in the long-term measurement of TOAs.

The results of the three numerical solutions can be tested by TEMPO. Put the orbital parameters, such as $`P_b`$, $`x`$, $`e`$ of Table 1 (3,5), as well as Post-Kepler parameters $`\dot{P}_b`$, $`\dot{\omega }`$ of Table 2 (4,6) into the corresponding terms in DD model. Then behavior of PSR B1937+21 should be very similar to PSR B1957+20 and PSR J2051-0827, since they have similar companion mass and close orbital period. Moreover, when the orbital precession induced time delay of Eq(11) is added into Eq(8), then the systematic variation of residuals of TOAs should be eliminated. Since different initial time corresponds to different phase angles, so the initial phases, like $`\omega _i`$, $`\varphi _i`$ may need to be adjusted when fitted by TEMPO.

Assuming PSR J2051-0827 and PSR B1957+20 are fitted as a singular pulsar, then they should have similar the residual of TOAs (after counting out DM variation) as that of PSR B1937+21. And binary pulsars such as, PSR J0437-4715, with much larger $`P_b`$, should take much longer time to finish the similar variation in residuals of TOAs like PSR B1937+21, in the case that PSR J0437-4715 is fitted as a singular pulsar.

The evolutionary links between millisecond pulsars and their binary progenitor systems involve a number of exotic astrophysical phenomena. In their late stages of evolution, neutron stars in low-mass X-ray and pulsar binaries may evaporate their companions through the strength of their radiation (Alpar et al. 1982, Bhattacharya $`\&`$ van den Heuvel 1992, King et al 2003). PSR B1957+20 has provided strong evidences to that scenario, the presumed binary pulsar PSR B1937+21 might provide another one.

By the spin angular momentum of the pulsar given in Table 1,3,5, one can easily obtain its moment of inertia, $`I_1=1.5I_{45}`$, $`2.6I_{45}`$, $`0.44I_{45}`$ for solution 1,2,3, respectively. (since pulsar period is $`P=1.5578`$ms), which are consistent with the prediction of neutron star equation of state. It is interesting that PSR B1937+21 having the smallest pulsar period, $`P`$, might also has the smallest $`x`$ of all millisecond pulsars.

I would like to thank Prof. T.Lu for constructive suggestions. I also thank Prof.Z.R.Wang, X.D.Li, Z.Y.Li, Z.D.Dai, D.M.Wei for helpful discussion on the manuscript.


text:
# Conductance distribution in nanometer-sized semiconductor devices due to dopant statistics

## I Introduction

In semiconductor physics, the influence of doping is generally accounted for by a homogeneous shift of the Fermi-level, caused by the introduction of free carriers in the semiconductor. However, when the dimensions of a device are small compared to the average distance between individual dopants, the discrete nature of doping must be taken into account. Each individual ionized dopant introduces a Coulomb potential well in the semiconductor, locally distorting the potential landscape. When the number of dopants in the volume of semiconductor that determines the transport characteristics of a device gets small, these random potential fluctuations cause atypical behavior of semiconductor devices Calvet et al. (2002). This is commonly viewed as one of the fundamental limits in the ongoing size-reduction of CMOS-technology Meindl et al. (2001); Hoeneisen and Mead (1972). Mapping the positions of individual dopants Voyles et al. (2002) and the potential fluctuations they induce Modesti et al. (2003); Rau et al. (1999) has been performed experimentally with various techniques. Furthermore, the influence of statistical fluctuations due to random dopants on device behavior has been subject of simulations Sano and Tomizawa (2001).

In this work, we experimentally investigate the effects of the discreteness of doping on the transport properties of small diodes by comparing many identically prepared devices. We find that statistical fluctuations caused by randomly positioned individual dopant atoms do not average out for very small devices. In contrast, fluctuations dominate the electrical transport properties of the smallest devices and cause large differences in the conductance of nominally equal devices. Furthermore, we demonstrate that at low temperature the Coulomb well of a single dopant gives rise to a resonant tunneling channel.

## II Measurements

In order to allow for measuring the transport characteristics of many identically prepared diodes, we use self-assembly methods to fabricate epitaxial CoSi2-diodes. The tip of a scanning tunneling microscope (STM) is used to characterize and access the devices individually.

All experiments are performed in an ultra-high vacuum (UHV) system with a base pressure of $`510^{11}`$ mbar. Self assembled CoSi2-islands are grown on Si-substrates (resistivity $`0.015\mathrm{\Omega }`$cm, doping concentration around $`210^{18}`$ cm-3) by evaporating a sub-monolayer of cobalt onto a clean $`7\times 7`$-reconstructed Si(111) surface, followed by an anneal at 800 C for about 5 min. The resulting hexagon-shaped epitaxial CoSi2-islands have heights in the range of 2–10 nm and diameters in the range of 15–80 nm. The inset of Fig. 2 displays an STM-image of a typical island. Each island is regarded as a nanometer sized epitaxial metal-semiconductor diode. The inter-island distances are roughly ten times larger than the island diameters. To minimize the effect of surface related transport channels, the $`7\times 7`$ surface reconstruction surrounding the islands is destroyed by exposing it to atomic hydrogen for 10 min at a substrate temperature of 400 C.

Current-voltage ($`IV`$) measurements are performed by positioning the STM-tip over an island. After switching off the feedback loop, the tip is lowered by 15 Å, which is found to be sufficient to make contact to the island. Then the current is measured while ramping the voltage, yielding well-reproducible $`IV`$-curves that reflect the properties of the metal-semiconductor contact. The sample preparation and measurements techniques are described in more detail elsewhere Smit et al. (2002a).

### II.1 Room temperature measurements

Some typical room-temperature measurements are displayed in the inset of Fig.1, showing weakly rectifying $`IV`$-curves, compatible with the $`p`$-type substrate-doping. The overall shape of the curves is very similar. To investigate the dependence of the diode’s conductances on their size, in Fig. 1 the zero-bias conductance per unit area is plotted as a function of island area for more than 40 different islands, grown on two similar $`p`$-type (boron-doped) samples. One would expect that devices of equal size yield the same value for this quantity. Indeed, for larger devices the conductance per unit area falls within a narrow range. However, towards smaller areas there is no definite behavior. Instead, the scatter in the measured values increases rapidly and nominally equal devices yield very different results.

To study the increasing scatter in more detail, the standard deviation of the conductance per unit area is plotted in Fig. 2 as a function of diode area, calculated from the data of Fig. 1. We clearly observe that the standard deviation increases with decreasing island size.

Measurement inaccuracies are not the source of this random scatter, as proved by the low noise level (typically less than 10 pA around zero bias) and good reproducibility of measurements on the same device. Moreover, this scatter is absent in similar measurements on low-doped samples (not shown). Schottky barrier inhomogeneities Tung (1992), which have been observed at non-epitaxial interfaces Talin et al. (1994); Meyer and von Känel (1997); Jäger et al. (2003), do not play a role here, as they are due to variations in the atomic arrangement at the metal-semiconductor interface itself. The metal-semiconductor interfaces in our devices are perfectly epitaxial and mono-crystalline, which is supported by cross-sectional transmission electron micrographs of similar structures Bennett et al. (2001). This also rules out interface defects or grain boundaries as the origin of the fluctuations. In addition, BEEM measurements of CoSi2-films on undoped Si(111) reveal a perfectly homogeneous Schottky barrier height Meyer and von Känel (1997). Only islands showing an atomically flat and defect-free surface in the STM-images were included in the analysis.

We will show that the increased scatter in the measured data is caused by the presence of randomly distributed dopant atoms in the Schottky barrier. An ionized dopant atom locally distorts the barrier, giving rise to a local barrier reduction and thus a high conductance spot (see the inset of Fig. 3). Then, the observed spread in Fig. 2 is directly related to the spread in the number of dopant atoms $`N`$ in the barrier of a device. From Poisson statistics, the relative spread is given by

$$\frac{\mathrm{sd}(N)}{N}=\frac{1}{\sqrt{N}},$$

where $`\mathrm{sd}(N)`$ is the standard deviation of $`N`$. Because $`N`$ is proportional to the island area $`A`$, this gives that the standard deviation of the number of dopants per unit area increases when the device area decreases. Indeed, Fig. 2 shows that a function of the form $`C/\sqrt{A}`$ describes our observations appropriately.

### II.2 Low temperature measurements

To further investigate the influence of individual dopant atoms on the transport properties of small diodes, we performed similar experiments in a low-temperature UHV-STM operating at 4.5 K. In the upper part of Fig. 3, several $`IV`$-curves of two different islands are plotted. Each represents a full and independent measurement cycle (stop scanning, make contact, ramp $`V`$, release contact, resume scanning), demonstrating their reproducibility. The current is plotted on a logarithmic scale, to make the features of the curves visible over several orders of magnitude in current. Note that the noise-level is below 1 pA.

From the measured $`IV`$-curves, the quantity $`(V/I)(\mathrm{d}I/\mathrm{d}V)`$ (the so-called normalized differential conductance) has been calculated. It was calculated from the average of several $`IV`$-curves per island and plotted for the same two islands in the lower part of Fig. 3. This is a well-known approach in the field of scanning tunnelling spectroscopy Stroscio et al. (1986) and is very useful as it reduces the overall exponential behavior and enhances bias-dependent features. Indeed, the weak features in the $`IV`$-curves turn into clear peaks and the peaks show a one-to-one correspondence to the features in the $`IV`$-curves. As can be seen from Fig. 3, features appear at one bias polarity, only. This is true for all acquired curves on the same sample. For substrates with the opposite doping type, the peaks occur at the opposite bias polarity (not shown).

The peaks can be explained by resonant tunneling through a discrete energy level of a dopant atom, occurring when the Fermi-level at either side of the barrier lines up with an energy level of the dopant’s potential well. The resulting resonant channel produces a feature in the measurement curve. Resonances are expected at bias voltages from zero upto the barrier height (roughly 0.5–0.7 V, depending on the dopant type). This is consistent with the observations. The actual bias voltage at which a resonance occurs is predominantly determined by the distance of the dopant atom to the interface. The number of peaks in a typical spectrum (1 to 4) corresponds to the expected number of dopants in the devices, as we will show later. This confirms our hypothesis that individual dopants influence the conduction path.

## III Analysis

Motivated by the foregoing observations, we present a model that links the conductance fluctuations observed at room temperature to (random) dopant positions. Moreover, the parameters in the model are directly related to device parameters.

Because the dopants in the substrate are randomly distributed, the number of dopants $`N`$ in the barrier of a given island is Poisson-distributed with parameter $`\lambda `$. This means that the probability $`P_\lambda (k)`$ that $`N`$ equals $`k`$ for a certain device is given by

$$P_\lambda (k)=\frac{\lambda ^k}{k!}e^\lambda ,$$
(1)
where $`\lambda =N`$ is the mean value of the number of dopants in the barrier of this particular device. The parameter can be expressed as $`\lambda =AtN_\mathrm{d}`$, where $`A`$ is the area of the island, $`t`$ the effective thickness of the barrier and $`N_\mathrm{d}`$ the average doping concentration in the barrier. As shown in Fig. 2, the scatter in the number of dopants in the barrier per unit area satisfies $`\mathrm{sd}(N)/A1/\sqrt{A}`$ and therefore is large for small islands.

We will now show that this effect also creates scatter in conductance measurements. We assume that dopants located in a certain region close to the metal-semiconductor interface induce a local barrier lowering that gives rise to a low-resistance transport channel 111Note that here the term ‘channel’ refers to a weak spot in a barrier and is unrelated to the conduction channels occurring in the field of quantum transport.. The conductance of a transport channel induced by a dopant atom depends on the dopant’s distance to the interface. Since the distance is a random variable, the conductance $`G_1`$ of a single channel is not the same for all channels, but has a certain probability distribution. The distribution is given by a probability density function $`f_1(g)`$, meaning that $`f_1(g)\mathrm{d}g`$ is the probability for a certain channel to have a conductance $`G_1`$ between $`g`$ and $`g+\mathrm{d}g`$. Both the distribution of the position of a dopant in the barrier and the dependence of the conductance of a channel on that position are contained in $`f_1(g)`$.

Assuming that the values of the conductance of the individual channels are independent and characterized by the same distribution given by $`f_1(g)`$, the total conductance of a device is given by the sum of the contributions of the individual channels. Here, we neglect the background conduction and assume that the conductance is dominated by these channels. If there are $`k`$ channels, the total conductance is $`G_k=_{i=1}^kG_1^{(i)}`$, where the conductance of the individual channels is denoted by $`G_1^{(i)}`$. The density function $`f_k(g)`$ of $`G_k`$ can be calculated explicitly by taking the $`k`$-fold convolution of $`f_1(g)`$ Grimmett and Stirzaker (1992)

$$f_k(g)=\underset{k\mathrm{times}}{\underset{}{f_1(g)\mathrm{}f_1(g)}}.$$

Finally, taking into account the Poisson distribution of the number of channels (Eq. 1), the density function $`f(g,\lambda )`$ of the total conductance $`G`$ of a device can be computed as

$$f(g,\lambda )=\underset{k=0}{\overset{\mathrm{}}{}}P_\lambda (k)f_k(g).$$
(2)
In other words, for given device parameters $`N_\mathrm{d}`$, $`A`$ and $`t`$, $`f(g,\lambda )\mathrm{d}g`$ is the probability that the total conductance $`G`$ due to dopant atoms in the diode is between $`g`$ and $`g+\mathrm{d}g`$.

To illustrate its behavior, Fig. 4 shows a map of $`f(g,\lambda )`$ (from Eq. 2) as a function of $`g/\lambda `$ (proportional to the conductance per unit area) and $`\lambda `$ (which is proportional to $`A`$). For a fixed value of the device area (that is a fixed value of $`\lambda `$ and a vertical line in the plot) the color scale gives the probability density to find a device with a particular conductance per unit area $`g/\lambda `$. For large values of $`\lambda `$, all density is concentrated around $`g/\lambda =\frac{1}{2}`$, while for smaller $`\lambda `$ it is spread over an increasingly wide range of values. In this figure $`f_1(g)`$ was chosen as a uniform distribution. However, as we will show next, the general properties of $`f(g,\lambda )`$ are not strongly dependent on the particular choice of $`f_1(g)`$.

Without making any assumptions on the choice of $`f_1(g)`$, we can compute the moments of $`f(g,\lambda )`$ in terms of those of $`f_1(g)`$. First, the mean value of $`G`$ satisfies

$$\begin{array}{cc}\hfill G=& gf(g,\lambda )dg=\underset{k=0}{\overset{\mathrm{}}{}}P_\lambda (k)gf_k(g)dg=\hfill \\ & \\ \hfill =& \underset{k=0}{\overset{\mathrm{}}{}}P_\lambda (k)kG_1=\lambda G_1.\hfill \end{array}$$

In fact, to allow for easy comparison with the data, we consider the total conductance *per unit area* $`\sigma `$. By making the substitution $`\sigma =G/A`$. It follows directly that the average value of $`\sigma `$ satisfies

$$\sigma =tN_\mathrm{d}G_1.$$
(3)
This is a very intuitive result, since $`\sigma `$ equals the average number of dopants in the barrier per unit area multiplied by the average conductance per channel.

To obtain the standard deviation we first observe that

$$\begin{array}{cc}\hfill G_k^2=& (G_1^{(1)}+\mathrm{}+G_1^{(k)})^2=\hfill \\ & \\ \hfill =& kG_1^2+k(k1)G_1^2,\hfill \end{array}$$

when the $`G_1^{(i)}`$ are independently identically distributed. Using this fact, we find that

$$\begin{array}{cc}\hfill G^2=& g^2f(g,\lambda )dg=\underset{k=0}{\overset{\mathrm{}}{}}P_\lambda (k)g^2f_k(g)dg=\hfill \\ & \\ \hfill =& \underset{k=0}{\overset{\mathrm{}}{}}P_\lambda (k)(kG_1^2+k(k1)G_1^2)=\hfill \\ & \\ \hfill =& \lambda G_1^2+\lambda ^2G_1^2.\hfill \end{array}$$

Finally, this yields

$$\mathrm{sd}(G)=\sqrt{G^2G^2}=\sqrt{\lambda G_1^2}.$$

By making the substitution $`\sigma =G/A`$ once more, we find that the standard deviation $`\mathrm{sd}(\sigma )`$ (which can be interpreted as the spread in $`\sigma `$) is given by

$$\mathrm{sd}(\sigma )=\sqrt{\frac{tN_\mathrm{d}}{A}G_1^2}.$$
(4)
The most important observation from this equation is that $`\mathrm{sd}(\sigma )`$ is proportional to $`1/\sqrt{A}`$, showing that the area-dependence of $`\mathrm{sd}(N)/A`$ leads to a similar behavior of spread in $`\sigma `$. This also justifies the choice of the fit-function in Fig. 2.

## IV Application to the data

Our simple model captures the general features of the data and yields reasonable values for the parameters. To demonstrate this, a least square fit of the function $`C/\sqrt{A}`$ (cf. Eq. (4); $`C`$ is a fit-parameter) to the standard deviations has been performed. The result is the dashed line in Fig. 2. The fit gives a good description of the data, showing that the spread in the data is consistent with the prediction of the model.

From the fit, we make some estimates for the physical quantities in the model. The value of the fit-parameter $`C=1.8`$ nA V-1 nm-1 should be equal to $`\sqrt{tN_\mathrm{d}G_1^2}`$ (according to Eq. (4)). By looking at the large area values in Fig. 1, we find that $`tN_\mathrm{d}G_10.04`$ nA V-1 nm-2 (Eq. (3)). Combining these numbers and assuming that $`G_1^2G_1^2`$ we find values for parameters in the model. First, $`1/tN_\mathrm{d}1/(2500\mathrm{nm}^2)`$, which corresponds to an average of one dopant per 2500 nm2 device area (indicated in Fig. 2). This number is consistent with e.g. an average doping concentration at the interface 222The doping concentration close to the silicon surface can be more than an order of magnitude lower than in the bulk due to out-diffusion during the sample preparation, as confirmed by simulations. around $`N_\mathrm{d}=10^{17}`$ cm-3 and an effective barrier thickness of $`t=2.5`$ nm. Note that because $`t`$ is the thickness of the barrier region where dopants influence the barrier height, it is thinner than the total Schottky barrier thickness (a few tens of nanometers in this case). Second, we find $`G_1100`$ nA/V for the average conductance per channel. In order to achieve this value, it is necessary to have e.g. a small patch with a local barrier height of approximately $`0.15`$ eV and an area 333This estimate was made using the thermionic emission model \[S. M. Sze, *Physics of Semiconductor Devices* (Wiley, New York, 1981), 2nd ed.\]. of approximately 10 nm2. These numbers are consistent with the actual sample parameters.

On average, the conductance per unit area increases for decreasing island area, as can be seen from Fig. 1. This can be explained by the scaling mechanism described in Ref. Smit et al., 2002b, which predicts that the barrier thickness decreases with device size for devices that are smaller than a few times the Debye length. We note that the increasing value of $`\sigma `$ is not the cause of the increment of $`\mathrm{sd}(\sigma )`$, since similar measurements on low-doped samples do not exhibit an increased scatter (not shown).

As demonstrated by our measurements, the discreteness of doping is easily observed in highly doped samples. In lower doped samples, where the barrier is thicker, the effect is expected to be much weaker. The local distortion of the potential landscape due to the presence of a dopant is roughly as large as its effective Bohr-radius, which equals about 3 nm in silicon. When the barrier thickness is much larger than the dopant’s potential well, the effective barrier lowering will be negligible. Hence the induced local barrier lowering is the most pronounced in thin barriers.

## V Conclusions

In conclusion, we have shown that individual dopant atoms dominate the transport characteristics of epitaxial nanometer sized metal semiconductor diodes. Room temperature data show increasingly large device-to-device conductance fluctuations towards smaller device sizes. Measurements at 4.5 K reveal pronounced structure in the $`IV`$-curves. A statistical analysis based on the assumption of randomly positioned individual dopant atoms leads to a good description of the experimental data.

###### Acknowledgements.

We thank S. G. Lemay for making the low-temperature measurements possible. One of us, S.R., acknowledges the Royal Netherlands Academy of Arts and Sciences for financial support. This work is part of the research program of the Stichting voor Fundamenteel Onderzoek der Materie, which is financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek.