text:
# Formalism and quality of a proper motion link with extragalactic objects for astrometric satellite missions

## 1 Introduction

Since 1998 the reference frame for astrometric measurements is represented by the International Celestial Reference Frame (ICRF, Ma et al. 1998), which is the realisation of the International Celestial Reference System (ICRS, Arias et al. 1995), defined by the International Earth Rotation Service (IERS). The ICRF consists of 212 VLBI observed quasars as defining extragalactic radio sources. These objects are assumed to have no proper motion and no overall global rotation – the quasars represent a quasi inertial system.

In the optical, the ICRS is represented today by the Hipparcos catalogue (ESA 1997). To be consistent, the Hipparcos catalogue had to be related to the ICRF. This means that it had to be assured that $`1^{\mathrm{st}}`$, the coordinate axes of the Hipparcos system are aligned with those of the ICRF and $`2^{\mathrm{nd}}`$, that the coordinate axes have no additional time dependent rotation component. Differences would show up as a global offset in positions and in a fictitious proper motion component. The transformation between the two systems is in general given by a rotation of the positional and proper motion system. This is the so-called *link algorithm*.

For the Hipparcos satellite it was not possible to link the catalogue directly with extragalactic objects, because the only extragalactic object in the Hipparcos catalogue is the quasar 3C 273 (HIP 60936). Consequently intermediate observations were needed (Kovalevsky et al. 1997): VLBI, MERLIN and VLA observations of radio emitting stars, proper motions determined with the HST fine guidance system, and photographic observations with respect to extragalactic objects made for example in Bonn (Geffert et al. 1997; Tucholke et al. 1997). The final link solutions were accurate to $`0.6\mathrm{mas}`$ for the positions and $`0.25\mathrm{mas}\mathrm{yr}^1`$ for the proper motion system (Kovalevsky et al. 1997).

The present study was initiated for the astrometric satellite mission AMEX111The missions DIVA (Röser et al. 1997) and FAME (Horner et al. 1998) are both cancelled, but the new satellite mission AMEX is a joint American-German merger of them., now proposed by the US Naval Observatory, and the ESA GAIA mission (Perryman et al. 2001), with a further step in astrometric accuracy, foreseen to be launched no later than 2012.

With the expected performance of AMEX some hundreds or even thousands of QSOs will be detectable, for GAIA even about 500 000 QSOs. If extragalactic objects are detectable it is possible to calibrate the proper motion system without further observations, so this would be a kind of self-calibration.

The aim of this paper is to determine the accuracy of the link of the proper motions of the astrometric satellite missions. We present Monte-Carlo simulations of observations of QSOs with the performance expected for the AMEX satellite and also for the GAIA satellite mission. These simulations allow us to estimate the accuracy of a proper motion link using extragalactic objects.

In Sect. 2 the link algorithm will be discussed. Some details about the implementation are described in Sect. 2.2. In Sect. 3 we present a search program for galaxies that could also be used for a proper motion link. In Sect. 4.1 we give a brief discussion of the underlying model to simulate the observations and the testing is described in Sect. 4.2. Finally, in Sect. 5 we present the results of our simulations.

## 2 The link algorithm

In this section we will give a complete mathematical description of the approximation of a small angle rotation and how it can be used to determine the components of the rotation vectors.

The link algorithm has been used in several publications like Frœschlé & Kovalevsky (1982), Arias et al. (1988) or Lenhardt et al. (1997). In Lindegren & Kovalevsky (1995) a theoretical study of the underlying exact model can be found. However, in most of these papers an approximation algorithm was used and the correction terms that will be described later (Eq. 13) were not taken into account, especially for the Hipparcos link program. We will show that with improved accuracy of astrometric data it is important to include these correction terms.

The basic idea of the algorithm is that in general the transformation between any two coordinate systems with the same origin is given by a rotation. If the coordinate axes of these systems are nearly coincident a first order approximation of the rotation can be used. It is required that the catalogues representing these two systems have a high degree of internal consistency. This is the case for the ICRF, for the Hipparcos catalogue, and it will be for future astrometric catalogues.

### 2.1 The first-order approximation

Let $`𝒙`$, $`𝒚`$ and $`𝒛`$ be the axes of a cartesian coordinate system: $`𝒙`$ points towards the origin of right ascension, $`𝒛`$ towards the celestial north pole and $`𝒚=𝒛\times 𝒙`$. Following Murray (1983), the coordinates of a unit vector $`𝒓`$ towards a celestial object in this coordinate system are given in general by

$$𝒓=\left(\begin{array}{c}\mathrm{cos}\delta \mathrm{cos}\alpha \\ \mathrm{cos}\delta \mathrm{sin}\alpha \\ \mathrm{sin}\delta \end{array}\right).$$
(1)
The two orthogonal vectors $`𝒑`$ and $`𝒒`$ which point towards positive right ascension and declination, respectively, are given by

$$𝒑=𝒛\times 𝒓=\left(\begin{array}{c}\mathrm{sin}\alpha \\ \mathrm{cos}\alpha \\ 0\end{array}\right)$$
(2)
and

$$𝒒=𝒓\times 𝒑=\left(\begin{array}{c}\mathrm{sin}\delta \mathrm{cos}\alpha \\ \mathrm{sin}\delta \mathrm{sin}\alpha \\ \mathrm{cos}\delta \end{array}\right)$$
(3)
where the brackets $``$ denotate vector normalisation.

Any small rotation about an arbitrary axis $`𝒆`$ with value $`|\epsilon |1`$ is given by

$$\mathrm{\Delta }𝒓=𝜺\times 𝒓$$
(4)
with the rotation vector $`𝜺=|\epsilon |𝒆`$. As shown in Murray (1983) the vector $`𝜺`$ can be split into a rotation $`\mathrm{\Delta }\alpha `$ about $`𝒛`$ and a rotation $`\mathrm{\Delta }\delta `$ about $`𝒑`$:

$$𝜺=\mathrm{\Delta }\alpha 𝒛\mathrm{\Delta }\delta 𝒑$$
(5)
$`\mathrm{\Delta }\alpha `$ and $`\mathrm{\Delta }\delta `$ are observable quantities: differences between right ascension and declination of two astrometric catalogues. It then follows

$$\mathrm{\Delta }𝒓=\mathrm{\Delta }\alpha \mathrm{cos}\delta 𝒑+\mathrm{\Delta }\delta 𝒒$$
(6)
and with the Eqs. (2) and (3), one easily finds

$`\mathrm{\Delta }\alpha \mathrm{cos}\delta `$ $`=`$ $`\mathrm{\Delta }𝒓𝒑`$ (7)
$`=`$ $`\mathrm{sin}\delta \mathrm{cos}\alpha \epsilon _1\mathrm{sin}\delta \mathrm{sin}\alpha \epsilon _2+\mathrm{cos}\delta \epsilon _3`$
and

$$\mathrm{\Delta }\delta =\mathrm{\Delta }𝒓𝒒=\mathrm{sin}\alpha \epsilon _1\mathrm{cos}\alpha \epsilon _2$$
(8)
where $`\epsilon _i`$ denotes the component of the vector $`𝜺`$. For the proper motion components one finds

$`\mathrm{\Delta }\mu _\alpha `$ $`=`$ $`\mathrm{cos}^2\delta \mathrm{cos}\alpha \mu _\delta \epsilon _1+\mathrm{tan}\delta \mathrm{sin}\alpha \mu _\alpha \epsilon _1`$ (9)
$`\mathrm{cos}^2\delta \mathrm{sin}\alpha \mu _\delta \epsilon _2\mathrm{tan}\delta \mathrm{cos}\alpha \mu _\alpha \epsilon _2`$
$`\mathrm{tan}\delta \mathrm{cos}\alpha \omega _1\mathrm{tan}\delta \mathrm{sin}\alpha \omega _2+\omega _3`$
and

$`\mathrm{\Delta }\mu _\delta `$ $`=`$ $`\mathrm{cos}\alpha \mu _\alpha \epsilon _1+\mathrm{sin}\alpha \mu _\alpha \epsilon _2`$ (10)
$`+\mathrm{sin}\alpha \omega _1\mathrm{cos}\alpha \omega _2`$
with $`\dot{\alpha }=\mu _\alpha `$, $`\dot{\delta }=\mu _\delta `$ and $`\dot{\epsilon _i}=\omega _i`$.

Finally the Eqs. (7) – (10) can be written in a matrix equation that can be solved for the components of the rotation vector $`𝜺`$ and the spin vector $`𝝎`$

$$\left(\begin{array}{c}\mathrm{\Delta }\alpha \mathrm{cos}\delta \\ \mathrm{\Delta }\delta \\ \mathrm{\Delta }\mu _\alpha \mathrm{cos}\delta \\ \mathrm{\Delta }\mu _\delta \end{array}\right)=\left(\begin{array}{cc}A_{11}& 0\\ A_{21}& A_{11}\end{array}\right)\left(\begin{array}{c}\epsilon _1\\ \epsilon _2\\ \epsilon _3\\ \omega _1\\ \omega _2\\ \omega _3\end{array}\right)$$
(11)
with the sub-matrices

$$A_{11}=\left(\begin{array}{ccc}\mathrm{sin}\delta \mathrm{cos}\alpha & \mathrm{sin}\delta \mathrm{sin}\alpha & \mathrm{cos}\delta \\ \mathrm{sin}\alpha & \mathrm{cos}\alpha & 0\end{array}\right)$$
(12)
and

$$A_{21}=\left(\begin{array}{ccc}\left(\begin{array}{cc}\frac{\mathrm{cos}\alpha }{\mathrm{cos}\delta }\mu _\delta +& \\ \mathrm{sin}\delta \mathrm{sin}\alpha \mu _\alpha & \end{array}\right)& \left(\begin{array}{cc}\frac{\mathrm{sin}\alpha }{\mathrm{cos}\delta }\mu _\delta & \\ \mathrm{sin}\delta \mathrm{cos}\alpha \mu _\alpha & \end{array}\right)& 0\\ \mathrm{cos}\alpha \mu _\alpha & \mathrm{sin}\alpha \mu _\alpha & 0\end{array}\right)$$
(13)
The Equation (13) is important: this term corrects for the orientation difference of the vectors $`𝒑`$ and $`𝒒`$ on slightly different places on the celestial sphere. This is important because these are the directions in which the proper motion components are measured.

For each celestial object one gets the 4 observational properties and the matrix in Eq. (11) which can be concatenated to a single matrix equation. For values $`\mathrm{\Delta }X`$ we use the convention (reference $``$ satellite) $`\mathrm{\Delta }X=X_{\mathrm{ref}}X_{\mathrm{sat}}`$ to calculate the matrix we use the reference catalogue.

### 2.2 Solving the equation

The matrix equation (11) can be solved by a weighted linear least-squares method. We use the singular value decomposition (SVD) because this algorithm is numerically the most appropriate (Press et al. 2002). For the SVD, vector, and matrix operations we use the C++ library newmat10A by R. Davies222Source code available for free from http://www.robertnz.net/nm10.htm. Details about the software implementation can be found in Metz (2003).

## 3 Galaxies for a proper motion link

In principle galaxies are, besides quasars, also good candidates for a proper motion link if they are sufficiently far away and well detectable for the instrument. Ideal candidates would be galaxies that have a star-like core region for which the position can be measured by the satellite instrument.

While – as mentioned above – only one extragalactic object can be found in the Hipparcos catalogue, the Tycho measurements should contain data for several galaxies. Since the Tycho-2 Catalogue (Høg et al. 2000) has a completeness of 90% at $`V=11.5\mathrm{mag}`$, one would expect some bright galaxies in this catalogue – however, not yet recognised as galaxies. On the other hand these objects may be important, since their positions from the Tycho-2 mission may be used as first epoch material. In this paragraph we have identified these galaxies and discuss the possibility of using them for the proper motion link of future satellite missions.

### 3.1 Search program

To identify galaxies in the Tycho-2 Catalogue we have compared the positions of all objects in the *Tycho-2 main catalogue* and the *Tycho-2 supplement 1* (Høg et al. 2002) with catalogues of galaxies and quasars available at the CDS FTP-Server (ftp://cdsarc.u-strasbg.fr/). The most important were the *Updated Zwicky Catalogue* (UZW), the *Catalogue of principal galaxies* (PGC), the *Third Reference Catalogue of Galaxies* (RC3) and the *Arcsecond positions of UGC galaxies* (AUGC) catalogue.

All objects found were further reviewed: we checked each single object with Digitized Sky Survey (DSS) images, whenever possible in different colors. In almost all cases miss-identifications due to stars in front of galaxies could be eliminated. All other objects were further checked with the NASA/IPAC Extragalactic Database (NED, http://nedwww.ipac.caltech.edu/) and Simbad (http://simbad.u-strasbg.fr/Simbad/).

Finally the galaxies were classified in two categories: “galaxies in the Tycho-2 Catalogue” and “*uncertain* galaxies in the Tycho-2 Catalogue”. Objects were classified as uncertain in the case of: galaxies in dense star fields, a comment in the database that a star may be in front of the galaxy, or a clearly visible positional offset in the DSS images. Galaxies found in the Tycho-2 Supplement 1 were not classified as certain/uncertain.

### 3.2 Galaxies in the Tycho-2 Catalogue

A total number of 181 galaxies were found: 116 galaxies in the Tycho-2 Catalogue, 35 uncertain galaxies in the Tycho-2 Catalogue and 30 galaxies in the Tycho-2 Supplement 1. The tables with all galaxies are available at CDS and the celestial distribution of the galaxies is shown in Fig. 1. Beside the galaxies, the quasar 3C 273 that was included in the Hipparcos catalogue was found as well.

All galaxies found have a PGC number but not all were found by the position cross correlation with the PGC due to the improper positional data in the PGC (and other galaxy catalogues as well). So there may be more galaxies in the Tycho-2 Catalogue which could not be identified by the applied method without more accurate galaxy catalogues.

However the galaxies found are hardly usable as first epoch material for a proper motion determination since they have mean position errors of $`\overline{\sigma _\alpha }=117\pm 22\mathrm{mas}`$ and $`\overline{\sigma _\delta }=127\pm 24\mathrm{mas}`$. With second epoch material e.g. from GAIA the resulting proper motion error would be of the order of $`\sigma _\mu 6\mathrm{mas}\mathrm{yr}^1`$. The large errors are mainly due to the fact that the galaxies are close to the limit of the Tycho-2 Catalogue with a mean magnitude of $`V_\mathrm{T}=12.5\mathrm{mag}`$. On the other side one would expect that a few months of observations with GAIA would produce results with a higher accuracy than it could be possible with help of the Tycho-2 Catalogue data.

### 3.3 Astrometric Observations

For three galaxies we made observations at our 1 m-Cassegrain telescope at the Observatorium Hoher List with the CCD camera HoLiCam (Sanner et al. 1998) in October 2002: PGC 7223, 11404, and 65866. For these galaxies we determined the internal position accuracy with a software package developed by Geffert et al. (1997), by fitting a star profile to the galaxy position. As reference we used the Tycho-2 Catalogue. For PGC 7223 and PGC 65866 the positions could be determined equally good as for stars in the fields with mean errors $`\begin{array}{c}<\\ \end{array}15\mathrm{mas}`$. Both galaxies have a very bright, star-like core region. For PGC 11404 we found a superimposed faint star but with a position clearly offset from the centre of the galaxy and from the Tycho-2 Catalogue position.

The results give a hint that in principle the positions of most of the galaxies found could be determined quite well even if a star profile were fitted, as we did.

## 4 Monte-Carlo Simulations: Models and tests

### 4.1 Models for the simulations

To analyse the possibility of linking the proper motion system of an astrometric satellite with QSOs we made Monte-Carlo simulations of observations of QSOs and for comparison also of observations of radio stars. The steps of the simulation were as follows:

1. Creation of a (synthetic) reference catalogue.
2. Creation of a satellite catalogue by applying an exact rotation model.
3. Applying an error model for the reference catalogue.
4. Applying an error model for the satellite catalogue.
5. Using the link algorithm to obtain the rotation and spin vectors.
6. Comparison with the vectors used in step B.

To create a reference catalogue we could either use a totally synthetic model that creates equally distributed apparent brightnesses and proper motions or make use of existing catalogues to create lists of objects with adopted apparent brightness and proper motion distributions. The first model is useful for testing purposes, the second one for simulations of observations. To apply a rotation to the synthetic reference catalogue we used the exact rotation equations (Lindegren & Kovalevsky 1995).

Errors could be simulated with the following models: no errors; normally distributed errors for testing; errors depending on the apparent brightness of a simulated object should be used to simulate the performance of an astrometric satellite (see Sect. 5); error distribution of existing catalogues could be used to simulate a realistic error distribution for example of radio observations. The adopted errors were applied to the catalogues by assuming each single error to be normally distributed around zero.333This does not imply that the total error distribution must be normal distributed!

We allowed different parameters for the simulations: total number of simulations $`N_{\mathrm{sim}}`$, number of simulated catalogues $`N_{\mathrm{cat}}`$, and number of simulated objects $`N_{\mathrm{obj}}`$. For each simulation we calculated the RMS of the difference between input and derived rotation and spin vector, $`\sigma (\epsilon _i)`$ and $`\sigma (\omega _i)`$, which are a measure for the mean expected errors of each vector component.

### 4.2 Testing the algorithm

We made detailed tests of the simulations and the algorithm. The accuracy of the link algorithm itself could be tested with simulations without errors. We varied the following parameters:

* number of objects $`N_{\mathrm{obj}}`$
* absolute value of the rotation and spin vectors $`|\epsilon |`$ and $`|\omega |`$
* the absolute value of the simulated proper motion of objects $`\overline{\mu }`$

Varying the number of objects, $`N_{\mathrm{obj}}`$, does not lead to a change in accuracy.

The effects of the correction terms of Eq. (13) were tested by varying $`|\epsilon |`$ and $`|\omega |`$ with fixed mean proper motions, which is shown in Fig. 2. The error of the position link $`\sigma (\epsilon _i)`$ increases with increasing $`|\epsilon |`$. The value $`\sigma (\omega _i)`$ increases similarly if $`|\epsilon |`$ and $`|\omega |`$ increase simultaneously. If in contrast the correction terms are *not* taken into account, the error $`\sigma (\omega _i)`$ is larger by a factor $`500`$. The error $`\sigma (\omega _i)`$ increases also if only the rotation vector $`|\epsilon |`$ is varied whereas the value $`\sigma (\epsilon _i)`$ is not affected by any variation of $`|\omega |`$.

The different behaviour for small values of $`|\epsilon |`$ and $`|\omega |`$ is a numerical effect. Important for the numerical accuracy of the link algorithm are differences between two values, for example $`\mathrm{\Delta }\alpha `$. Whereas for the proper motions one has to store small differences of small values, for positions one has to store small differences of large values. The smallest absolute value $`|\mathrm{\Delta }\alpha |`$ that could be stored in a double precision floating point number for a object with $`\alpha >15\stackrel{h}{.}3`$ is $`|\mathrm{\Delta }\alpha |=2^{352}410^7\mathrm{mas}`$.444For a IEC/IEEE 64 bit floating point number, e.g. a double in C++ on i386 architecture.

Also visible in the plot are little “bumps”. These appear because the RMS varies due to the Poisson-like distribution while a median would be more appropriate. However, in the plot we show the RMS because this it better comparable with the further results and shows the same basic behaviour as a median.

Further tests showed that the mean errors are independent of the value of the mean proper motion $`\overline{\mu }`$ if the correction terms are used. If the simulated objects have no proper motion the results do not differ whether the correction terms are taken into account or not.

We also made simulations with normally distributed errors for positions and proper motions of both catalogues. If we fix the mean errors of both catalogues and vary the number of link objects $`N_{\mathrm{obj}}`$ the resulting error could, as one would expect, be fitted by a function $`1/\sqrt{N_{\mathrm{obj}}}`$.

The relation between the resulting link mean error and the mean errors of the catalogue is shown in Fig. 3. We plotted the error $`\sigma (\omega _i)`$ depending on the size of the mean errors $`\mathrm{\Delta }`$ of the simulated reference (ref) and satellite (cat) catalogues. The behaviour is as follows: if we fix the mean error of one catalogue, for example $`\mathrm{\Delta }_{\mathrm{cat}}=0\mathrm{mas}\mathrm{yr}^1`$, and increase the other one, the resulting error for the proper motion link $`\sigma (\omega _i)`$ increases linearly. If we now fix $`\mathrm{\Delta }_{\mathrm{ref}}=4\mathrm{mas}\mathrm{yr}^1`$ and increase $`\mathrm{\Delta }_{\mathrm{cat}}`$, one sees that the link error is further on dominated by the much larger error $`\mathrm{\Delta }_{\mathrm{cat}}`$ until the value of $`\mathrm{\Delta }_{\mathrm{ref}}`$ gets in the order of $`\mathrm{\Delta }_{\mathrm{cat}}`$. Finally, both contribute equally to the resulting error.

If we assume normally distributed errors for both catalogues a $`\chi ^2`$-Test could be applied, but it fails if we use models for astrometric satellites (see next section), because the errors are no longer normally distributed, which is a requirement for the $`\chi ^2`$-Test (Brandt 1999).

## 5 Simulations for astrometric satellite missions:
AMEX like satellites and GAIA

In Hirte (1998) the relation between the apparent brightness of a star and the parallax accuracy for the DIVA satellite is presented, which then could be transformed in a position and proper motion accuracy. We fitted two exponential functions through the data points (Hirte 2002) to get a smooth relation and we use this data as estimation for AMEX like satellites (AMEX model). For GAIA we used data from Perryman et al. (2001) and fitted a single exponential function (GAIA model). We then used these models to simulate a proper motion link with different kinds of objects.

### 5.1 Simulations with quasars

To simulate observations of quasars we used the catalogue *Quasars and Active Galactic Nuclei (QAGN)* (Véron-Cetty & Véron 2001) to create the reference catalogue with a luminosity distribution of quasars as realistic as possible. We choose different limiting magnitudes to determine how the accuracy behaves with satellite performance. All proper motions and their errors were set to $`0\mathrm{mas}\mathrm{yr}^1`$ for the reference catalogue. This means that the resulting accuracy is only depending on the expected errors of the satellite observations.

The parameters for the AMEX like satellite simulation were a limiting magnitude of $`V_{\mathrm{max}}=15\mathrm{mag}`$ ($`N_{\mathrm{obj}}=88`$ QSOs) and we made $`N_{\mathrm{sim}}=\mathrm{10\hspace{0.17em}000}`$ simulations. In Fig. 4 the results of one of those simulations are shown. It is clearly visible that the values used to create the satellite catalogue (marked by the cross) are found on average but the results spread over a large range.

The results of all simulations performed are given in Table 1. The expected errors of the proper motion link for GAIA are, at the same limiting magnitude and the same number of objects, a factor of $`400`$ – 800 smaller than for the AMEX model, based on the expected GAIA performance (Perryman et al. 2001).

The mean error of the proper motion link for the AMEX model is $`\sigma (\omega _i)=0.61\mathrm{mas}\mathrm{yr}^1`$ at a limiting magnitude of $`V_{\mathrm{max}}=15\mathrm{mag}`$. If we include more QSOs by using fainter objects the mean error decreases only slightly. This is due to the fact that simultaneously the error of the proper motion determination increases for fainter objects. So the profit of using more objects is essentially cancelled out by the worsened performance. For GAIA the error decreases steadily.

To study the influence of the individual parameters we have also varied the number of QSOs and the length of the mission, both increased by a factor of 1.5. The results were as expected: if we increase the number of objects by a factor of 1.5 the error of the proper motion link decreases by a factor of $`\sqrt{1.5}`$, whereas a 1.5 times longer mission leads to a link error that is by a factor $`\sqrt{1.5^3}`$ smaller (Hirte 1998).

In Kovalevsky (2003) the effect of the curvature of the Solar system barycentric motion around the centre of our galaxy is described. This appears as an apparent proper motion of up to $`4\mu \mathrm{as}\mathrm{yr}^1`$. Our simulations show that this affects the proper motion link by an error that is a factor of $``$ 10% larger for all limiting magnitudes for the GAIA model if the data is not corrected for the aberration effect (which can be done before the link is applied). The apparent proper motion has no significant influence on the AMEX model with quasars.

It has to be taken into account that the QAGN is not complete for faint objects. If we extrapolate the number of quasars brighter than $`V_{\mathrm{max}}=16.5\mathrm{mag}`$ in the *SDSS Quasar Catalog I* (Schneider et al. 2002) we find that at this limiting magnitude about $`57\%\pm 16\%`$ of the quasars are in the QAGN whereas only $`5.4\%\pm 0.1\%`$ brighter than $`V_{\mathrm{max}}=19\mathrm{mag}`$ are in the QAGN. So for GAIA we underestimate the number of quasars in our simulations by a large factor.

### 5.2 Simulations with galaxies

As described in Sect. 3 we identified 116 galaxies in the Tycho-2 Catalogue. We used the position and magnitude data from the Tycho-2 Catalogue to simulate a proper motion link for AMEX like satellites with galaxies. All proper motions and their errors were set to $`0\mathrm{mas}\mathrm{yr}^1`$ for the reference catalogue, for the satellite errors we used the AMEX model. The resulting mean errors of the proper motion link were $`\sigma (\omega _i)=0.16\mathrm{mas}\mathrm{yr}^1`$.

This result is a crude estimation. Firstly we made the assumption that the astrometric errors only depend on the magnitude but they may be affected by the source not being point-like. Secondly we used Tycho-2 magnitudes as estimation of the brightness of the core region. On the other side it is likely that more galaxies may be detected (see Sect. 3.2).

### 5.3 Simulations with radio stars

To compare the results of the simulations with quasars and galaxies with other methods we made simulations with radio-emitting stars (Lestrade et al. 1999) which were used for the Hipparcos link program. Radio stars are suitable link objects, because their proper motions may be determined with high accuracy within short time scales. For the construction of the synthetic reference catalogue the astrometric data from Lestrade et al. (1999) were used, magnitudes were taken from the Hipparcos catalogue. We varied the number of objects $`N_{\mathrm{obj}}`$ and in a second run the VLBI errors were reduced by a factor of 2.

The results are given in Table 2. The errors are mainly caused by the VLBI errors. The accuracy of the proper motion link could be increased by using more link objects which could be described by a function $`1/\sqrt{N}`$ as shown in Sect. 4.2 for normal distributed errors.

## 6 Conclusions and outlook

The main goal of the proper motion link of future astrometric satellite missions is to establish a proper motion system with an accuracy that is better than the accuracy of the single proper motions of one star. Our simulations show that for astrometric missions like GAIA this may be achieved by the use of QSOs alone. For astrometric satellite missions with AMEX like performance the use of quasars alone is not sufficient and galaxies have to be included. However, an accurate proper motion link requires external data from VLBI observations of radio-emitting stars or HST observations of stars with respect to extragalactic objects.

###### Acknowledgements.

We thank the referee Jean Kovalevsky for his useful comments as well as Ulrich Bastian and Klaas S. de Boer for useful discussions.


format: text

text:
# Fourfold oscillations and anomalous magnetic irreversibility of magnetoresistance in the non-metallic regime of Pr1.85Ce0.15CuO4

## Abstract

Using magnetoresistance measurements as a function of applied magnetic field and its direction of application, we present sharp angular-dependent magnetoresistance oscillations for the electron-doped cuprates in their low-temperature non-metallic regime. The presence of irreversibility in the magnetoresistance measurements and the related strong anisotropy of the field dependence for different in-plane magnetic field orientations indicate that magnetic domains play an important role for the determination of electronic properties. These domains are likely related to the stripe phase reported previously in hole-doped cuprates.

In the low doping regime of high temperature superconductors (HTSC), resistivity crosses over most often from a high temperature metallic-like state to a low temperature non-metallic behaviorAndo et al. (1996). Interestingly, the crossover temperature observed by in-plane resistivity does not correspond to the Néel temperature ($`T_N`$), the onset of long range antiferromagnetic order. In fact, there is no hint of this transition in $`\rho _{xx}(T)`$ as metallic-like behavior persists sometimes below $`T_N`$ Lavrov et al. (1999); Ando et al. (2001). Ando and coworkers argued that this independence of $`\rho _{xx}(T)`$ from the underlying magnetic order is a sign of phase separation, probably with conducting stripes separated by antiferromagnetic domains. Such self-organized structures of the charge carriers was shown to lead to intriguing behaviors for several physical propertiesAndo et al. (1999, 2003); Lavrov et al. (2001); Sun et al. (2003a). An anisotropic twofold response to magnetic field applied along the copper-oxygen (CuO2) planes was reported for instance using in-plane magnetoresistance (MR) in the non-metallic regime of YBa2Cu3O6.33 (YBCO)Ando et al. (1999) and in lightly doped La1.99Sr0.01CuO4 (LSCO)Ando et al. (2003).

To avoid the possible contributions of orthorhombic distortions or phases (as in YBCO and LSCO), we have chosen the tetragonal electron-doped cupratesFontcuberta and Fabrega (1996) as good candidates to make a definitive test of anisotropic transport properties. We focus on the results obtained for magnetic fields always applied parallel to the CuO2 planes. In Figure 1, we present an example of raw in-plane and out-of-plane resistance data as a function of angle with respect to the in-plane crystal a-axis. These data were measured deep into the non-metallic regime of non-superconducting Pr1.85Ce0.15CuO4 (PCCO). We observe clear fourfold and symmetrical oscillations of the resistance with sharp maxima \[minima\] for $`R_{xx}`$ \[$`R_c`$\] for fields applied along the in-plane crystal axis (i.e. along the CuO bonds of the CuO2 planes) and broad minima \[maxima\] for field applied along the diagonals. This effect can be emphasized using the polar plot of Fig. 1(c).

In this Letter, we present this magnetoresistance anisotropy, for current along the CuO2 planes and along the c-axis (thickness) of non-superconducting electron-doped single crystals. We show for the first time that sharp fourfold oscillations persist over the whole non-metallic regime for both in-plane ($`\rho _{xx}`$) and out-of-plane ($`\rho _c`$) resistivities. More importantly, the field dependence of $`\rho _c`$ presents irreversibilities which cannot be explained by conventional band theory. We ascribe the new transport signatures to the presence of magnetic domains. Their origin is probably related to the presence of stripe domains (order) in the cuprates.

The single crystals of Re2-xCexCuO4 (Re = Pr, Nd, Sm and Eu) used for this study were grown by the directional flux growth method in alumina and high purity magnesia cruciblesPeng et al. (1991); Kaneko et al. (1999). The as-grown crystals are known to be non-superconducting even for x = 0.15, and only a reduction process allows them to become superconducting with a maximum transition temperature between 18 and 25K depending on the rare-earth ionFontcuberta and Fabrega (1996). In the present work, we study only as-grown non-superconducting crystals with x = 0.15 in order to extend as much as possible the temperature range for which a non-metallic regime is observed.

Silver epoxy contacts were applied directly onto the as-grown crystals (typical in-plane size : 2mm x 1 mm, and 30 $`\mu m`$ thickness along the c-axis) in two different configurations insuring uniform current density for the measurement of the in-plane ($`\rho _{xx}`$) and c-axis ($`\rho _c`$) resistivity as shown in the inserts of Figs. 1(a) and (b). The samples were then mounted on sapphire supports and attached onto a specially designed Physical Property Measurement System (PPMS) rotator chip, such that the applied magnetic field can be rotated over a full 360. The direction of the in-plane crystal axis (a-axis direction relative to the edge of the rectangular crystals) were determined using Laue x-ray diffraction and by Raman scattering. Several tests were performed to rule out contact configuration, thermometry and misorientation problems.

In Figure 2(a), we present the typical temperature dependence of in-plane and c-axis resistivity components for non-superconducting Pr1.85Ce0.15CuO4 single crystals. Both components show a non-metallic regime at low temperatures where the MR oscillations are observed. The in-plane resistivity $`\rho _{xx}`$ is metallic-like at room temperature while it shows an upturn starting at $`T_{min,xx}125K`$ which is very sensitive to the oxygen content. Below $`T_{min,xx}`$, $`\rho _{xx}`$ approaches a $`\mathrm{ln}(T)`$ behavior (not shown), similar to that reported for non-superconducting NCCO and PCCOFournier et al. (2000). This temperature behavior and the strong anisotropy illustrated in Fig. 2(b) have been interpreted as a signature of two-dimensional weak localization (2DWL) by disorderLee and Ramakrishnan (1985). For the c-axis resistivity $`\rho _c`$ , a non-metallic trend is observed at room temperature for several crystals, followed by a maximum (at $`T200K`$)com (a), then by a metallic-like regime. It is finally followed by a non-metallic regime, also approaching $`\mathrm{ln}(T)`$ , below $`T_{min,c}70K`$ .

The field dependence of the in-plane magnetoresistance oscillations at 10K are presented in Fig. 3(a). The MR oscillations evolve steadily from broad features at 2T to sharp maxima for 4T whenever the field is applied along the Cu-O bonds. For fields approaching 9T, new maxima develop for fields applied approximately along the diagonal directions. The relative proportion of both oscillations varies also steadily with temperature as evidenced by the same 45 features barely appearing at 5T and 2K in Fig. 3(b) (see arrows). The amplitude of these oscillations in $`\rho _{xx}`$ represents only a small fraction ($`\mathrm{\Delta }\rho _{osc}/\rho 0.05\%`$) of the total negative MR \[$`\mathrm{\Delta }\rho _{tot}/\rho 1\%`$ : see Fig. 2(b)\]. Thus, unlike YBCO and LSCO, $`\rho _{xx}`$ oscillations are very weak and have almost perfect fourfold symmetry.

In Fig. 3(b), we present the evolution of these $`\rho _{xx}`$ oscillations (at 5T) with temperature. A similar effect is obtained if one decreases the temperature or increases the applied field. We observe that the weak oscillations disappear quickly with increasing temperature, vanishing close to $`T_{min,xx}`$. This seems to imply that this anisotropic behavior can only be observed in the non-metallic regime of these materials (within our sensitivity).

In Fig. 4(a), we present the angular dependence of $`\rho _c`$ in several magnetic fields at 5K. We observe very sharp minima in the c-axis resistivity, developing particularly for intermediate fields (3 to 6T). These sharp cusps appear for magnetic fields applied along the Cu-O bonds. As the applied magnetic field is further increased, these anomalies are gradually replaced by oscillations of smaller amplitude. In Figs. 4(b) and (c), we illustrate the strong temperature dependence of these features as they seem once again to vanish as the sample reaches temperatures close to the crossover to the metallic-like state. As will be shown below, the magnitude of the oscillations observed for $`\rho _c`$ are comparable to the total positive MR, in sharp contrast with those observed with $`\rho _{xx}`$. Therefore, they are less sensitive to mis-orientation of the crystals and thus easily observed.

In Figure 5 (a), we present the c-axis magnetoresistance at T = 5K as a function of in-plane magnetic field applied along three different directions (0, 15 and 45). Contrary to the in-plane resistivity, the c-axis magnetoresistance is positive and presents an unusual field dependence. For the field applied along the a-axis (at 0), the resistivity remains remarkably flat at low fields until a threshold field is reached. At this point, the resistance varies sharply, as if the system was crossing a transition. Beyond this threshold, the magnetoresistance resumes a high field $`B^2`$ behavior. As soon as the field direction deviates from the a-axis \[for 15 and 45 in Fig. 5(a)\], the magnetoresistance presents a very sharp positive increase at low fields, quickly reaching the saturation $`B^2`$ regime. The large magnitude of the oscillations in $`\rho _c`$ at about 4 - 5T in Fig. 4 (a) can easily be explained by the strong variations in resistance at 5T for different field orientations underlined by the dashed line in Fig. 5(a). We should emphasize here that $`\rho _{xx}`$ displays its sharpest peaks in the same range of applied field.

Interestingly, upon decreasing the applied magnetic field, the resistance for $`\theta =0^{}`$ shows a similar transition-like pattern, but the resistance presents a clear sign of irreversibility observed in Fig. 5(b). The arising of magnetic hysteresis is a clear demonstration of the influence of magnetic history, i.e. the presence of magnetic domains. These domains (and their domain walls) are strongly affected by the direction and strength of the in-plane rotating field, and in turn they affect the conductivitycom (b). We can interpret the data in Fig. 5(a) as evidence of easy magnetization axis along the diagonal directions, while the Cu-O bond directions (a and ”b” axis) correspond to hard axis. Since the c-axis resistivity is a measure of interplane tunneling between CuO2 planesSuzuki et al. (1999), the alignment of magnetic domains at high fields is detrimental to interplane tunneling causing an increase of c-axis resistivity. This effect is present up to temperatures approaching $`T_{min,c}`$ \[see Fig. 5(c)\]. Inversely, in-plane resistivity $`\rho _{xx}`$ presents peaks for fields applied along the hard axis and flat minima for field applied along the diagonal easy axis.

The behavior of the magnetoresistance anisotropy for PCCO is significantly different from YBCOAndo et al. (1999) and LSCOAndo et al. (2003). For YBCOAndo et al. (1999), there is a clear sign change when magnetic field is applied parallel and transverse to the applied in-plane current. In LSCO, no sign change, but the presence of four asymmetrical lobes was attributed to twinning. In our case, we never observe a sign change in the MR with angle (it always remains negative) and the oscillations are not sensitive to the direction of the applied current as they remain about the same height for the field parallel  and transverse  to the current, except for the additional features observed at high magnetic fields in $`\rho _{xx}`$. The magnitude of the in-plane MR oscillations is much smaller for PCCO than for YBCO and LSCO. In both cases however, the oscillations disappear close the non-metal to metal crossover. We suspect that orthorhombic distortions could be a key player in promoting the differences in the magnitude and the anisotropy of these MR oscillations in YBCO, LSCO and PCCO.

MR oscillations of c-axis resistivity in cuprates were first reported in strongly overdoped Tl-based cupratesHussey et al. (1996) with $`T_c25K`$. In this case, the fourfold oscillations are large (0.33$`\%`$ of the total resistance at 40K and 10T) and correspond to a doping region of the phase diagram where both $`\rho _{xx}`$ and $`\rho _c`$ remain metallic-like over the whole temperature range. Because the mean-free path (MFP) is fairly large in this strongly overdoped regime, Dragulescu et al. argued that the oscillations in Tl-based cuprates could be due to angular-dependent magnetoresistance oscillationsDragulescu et al. (1998). In our case, the very small MFP for this dopingFournier et al. (1997); Gollnik and Naito (1998) precludes such interpretation.

Recent experiments indicate that the carriers injected into the cuprates through chemical substitution have a strong tendency to distribute non-uniformly in the copper-oxygen planes, segregating into phase-separated regions, clusters and stripesCarlson et al. (2002). Several experimental observations, including the MR oscillations obtained with YBCOAndo et al. (1999) and LSCOAndo et al. (1999), fit into this possible scenario. However, only a recent report by Sun et al. indicates the possible presence of stripes in the electron-doped cupratesSun et al. (2003b). Assuming their existence in the electron-doped cupratescom (c), we should expect the c-axis resistivity to decrease whenever stripes in adjacent CuO2 planes are directly on top of each other and aligned in the same direction : in this particular case, electrons can tunnel more easily between planes, thus decreasing resistivity. However, our data show that $`\rho _c`$ is higher when the field is high (for well aligned domains). We suggest here that a low density of stripes in adjacent planes makes it more difficult to have stripes on top of each other whenever domains are well aligned (high fields), in particular if domains on adjacent planes are not or weakly correlated (as in LSCOMatsuda et al. (2000)). At low fields, random orientation of stripe domains could lead to a better overlap from weakly correlated adjacent planes, and thus to a better conductivity. For in-plane MR \[see Fig. 1(a)\], the application of in-plane magnetic field could promote a partial displacement of domain walls (the stripes)Kivelson et al. (1998), enough to change the resistivity and improve the channeling of electrons along longer conducting ”rivers of charges”. Because the in-plane MR oscillations are so small, the changes in domain wall configurations are probably scarse. Our data would suggest that such modifications of the wall configuration are hard to develop only along the hard axis (the Cu-O bonds), leading to maxima in $`\rho _{xx}`$.

We should mention that this simple scenario ignores completely the possibility of interaction with the underlying rare earth magnetic order (different in PrCeCuO, NdCeCuO, SmCeCuO)Hundley et al. (1989). Our preliminary data on NdCeCuO and SmCeCuO showed no significant difference for the field and angular dependence of $`\rho _{ab}`$ and $`\rho _c`$ , thus ruling out a direct implication of rare earth magnetism. Moreover, it remains unclear how an electronic system could present at the same time signatures consistent with 2DWL by disorder and one dimensional features like stripes, unless the spin stripes induce only a partial charge density wave in the CuO2 planes. This aspect will need further exploration, possibly through the doping dependence of the observed oscillations.

In summary, we presented anomalous magnetoresistance oscillations for non-superconducting electron-doped cuprates in their non-metallic regime for magnetic field applied along the copper-oxygen planes. We showed that sharp fourfold oscillations persist over the whole non-metallic regime for both in-plane ($`\rho _{xx}`$) and out-of-plane ($`\rho _c`$) resistivities. The field dependence of $`\rho _c`$ presents irreversibilities which can be explained by the presence of magnetic domains. Their origin is probably related to the presence of stripe domains in the cuprates.

We thank A.-M. Tremblay, C. Bourbonnais, K. LeHur and M. Li for several discussions and S. Pelletier for technical assistance. We acknowledge the support of CIAR, CFI, NSERC and the Fondation of the Université de Sherbrooke.


text:
# Non-deterministic random bit generator based on electronics noise

## 1 Motivation

Today’s computers are Turing machines governed by deterministic laws. It has been noticed that such machines can not solve certain class of problems because of their inability to produce random numbers. For example, it is impossible to simulate a simple radioactive decay. A computer program that would be able to do this must be able to produce a sequence of times of decays which is by no means deterministic.

Consequently, a complete machine that is believed to be the ultimate universal computing machine is referred to as ”Turing machine with a random number generator” or ”probabilistic Turing machine” QCbook ; Wikipedia-pTuring .

Since computers play an increasingly important role in modern science and life, it is also increasingly important to develop ”good” random number generators.

The classical approach is to approximate truly random, non-deterministic generator of random numbers by a carefully chosen mathematical function which produces approximately random numbers. Such a function of course can be calculated on Turing machines. Without loss of generality we will only consider random bit generators (RBG). The basic technique is as follows. First, one chooses at will an integer number $`x_0`$ (so called seed) from a large set of numbers that are known to be valid seeds for a given function $`f`$. (Good coverage of the theory and practice of pseudo-random generators may be found in CRC ; NRinC ; Knuth ; RANLUX ; RANLUX2 ; Hellekalek ; rreview ; Miller ; Marsag90 .) Then one calculates:

$`x_1`$ $`=`$ $`f(x_0)`$ (1)
Iterating this as follows:

$`x_{i+i}`$ $`=`$ $`f(x_i),i=1,2,3\mathrm{}`$ (2)
results in a sequence of numbers:

$$x_1,x_2,x_3,\mathrm{}$$
(3)
These numbers converted to bits and put together side-by-side form a long sequence of pseudo-random bits. Sometimes, not all bits from the numbers $`x_i`$ are used to form the sequence but only a pre-selected part (for example the most significant half).

Variants of this basic technique exists but the underlying property of all pseudo-random generators is that they must accept a seed, a form of mathematical initial state that completely determines the sequence of bits produced thereafter. The idea behind all this is that to someone who does not know the seed and/or function $`f`$ (or does not care about them), the sequence of bits produced appears to be random.

Pseudo-random numbers can be produced efficiently and are used a lot in simulations of stochastic processes like passage of particles through matter Geant  or speeding up of calculation of exact problems such as multi-dimensional integration or primality testing primality .

While good generator functions which approximate true randomness quite well for most purposes RANLUX ; RANLUX2 ; Hellekalek ; Marsag90 , it is absolutely clear that the entropy of the whole pseudo-random sequence, no matter how long, is equal to the entropy of the seed only, which is usually not more than 16 to 64 bits. This is most clearly seen by noting that the optimal compression of a sequence of type (3) would result only in the seed $`x_0`$ and a decompression routine that equals $`f`$. Since $`f`$ is usually publicly known it does not represent a useful source of entropy, thus the only entropy left is that of the seed. And even this tiny entropy, namely the entropy of the seed, has to be provided by somebody or something that has nothing to do with the pseudo-random generator itself. Therefore the conclusion that pseudo-random generators do not generate entropy seems unavoidable.

For applications in cryptography this may be a killing property. For example provability of information-theoretic security of Quantum cryptography protocols such as BB84 BB84 , or Maurer’s SKAPD protocol Mau93  assumes existence of generators of truly random (non-deterministic) numbers. Even much less complicated applications such as PIN or TAN number generators may not be proven secure if the numbers generated are based on a deterministic procedure. Cases are known where seeding by a low-entropy source (such as clock) led to serious compromise of the subsequent cryptographic protocol. An example is attack to the Netscape’s 40 bit RC4-40 rc4  challenge data and encryption keys which could be revealed in a minute or so dobbs96 . The authors of that article suspect that 128 bit version RC4-128 would not be much harder to break either if seeding is done in a similar fashion. (This example also shows how a perfectly good cryptographic method may be ruined by a non-educated or malicious implementator.)

All these cases boil down to the fact that it is hard to get truly random data form today’s computing machines (usually PC’s). Some ad-hoc methods exist Davis ; rfc1750 ; urandom  to extract a handful of random bits per second, but the only way to get a lot of random bits of guaranteed quality in a short time is to add a non-deterministic random bit generator to the computer, thus effectively realizing the probabilistic Turing machine.

Problem, however, is to construct a good enough generator, whose output can be for all practical purposes considered as being truly unpredictable and random.

We describe, in this article, a theory and one practical realization of a reliable and fast random bit generator which should not be expensive to produce.

## 2 The method

Since it is obvious that only physical reality may provide true randomness, the non-deterministic RBG should be based on a repeated measurements of identical, independent physics processes and associated method for extraction of random bits from these measurements. An RBG should be characterized by the following:

1. its output can not be predicted regardless of the amount of knowledge about it (just as it would be impossible to predict outcome of flipping a fair coin by knowing exactly how it looks);
2. two identical generators may not be synchronized to produce the same sequence of bits.

The second requirement could well be regarded as the definition of a non-deterministic generator. This requirement may be reformulated so that a non-deterministic generator may not accept an initial state. This is opposite to pseudo-random generators which must accept an initial state.

Since we do not expect that measuring of independent physics processes, that is making independent identical experiments whose outcome is random, may form any pseudo-random pattern we are left with only two problems associated with non-deterministic random bit generators. These are:

1. statistical bias, defined as:

$`b=p(1)0.5`$ (4)
where $`p(1)`$ stands for probability of ones;
2. correlations among bits, among which the most important serial correlation coefficient is defined as Knuth :

$`a`$ $`=`$ $`{\displaystyle \frac{n_{i=0}^{n1}b_ib_{(i+1)\mathrm{mod}n}\left(_{i=0}^{n1}b_i\right)^2}{n\left(_{i=0}^{n1}b_i^2\right)\left(_{i=0}^{n1}b_i\right)^2}}`$ (5)
where $`b_1,i=0\mathrm{}n1`$ denotes $`i`$-th bit in the sequence.

Both are measures of imperfections that are inevitable in practical realizations of generators. Correlations appear if the measured physics processes (experiments) are not completely statistically independent of each other, whereas statistical bias is mainly associated with imperfections in the measuring equipment (eg. electronics). Serial correlation coefficient is a measure of the extent to which a bit (in a sequence) depends on a previous bit (Knuth ), and takes on a value between -1 and 1. We suppose that correlations between bits further apart (corresponding to experiments further apart in time) are smaller if not negligible. For truly random sequences of bits, of course, both statistical bias and correlations tend to zero, when the length of the sequence goes to infinity.

The method presented here solves both problems. It consists of counting events which are result of measurement of some random physical process. For example events may be a radioactive decay or taping of rain drops on a tin roof. These events appear at random but measured for a long time they have some mean period, $`\mu `$. The state of the counter is examined at regular time intervals with the period $`T\mu `$. If the counter state is found to be odd then the output from the generator is ”1”, otherwise it is ”0”. Unlike some other methods which require exponential distribution of time intervals between events vincent70  or sampling of an analog white noise source bucci01 , we will show that this method, in the limit of slow sampling ($`T/\mu \mathrm{}`$) leads to vanishing bias and correlations regardless of the distribution of the random process being measured. Therefore it may used to extract random bits from a large variety of processes.

## 3 Practical realization of a random bit generator

Useful non-deterministic random bit generator built in hardware should satisfy several requirements:

1. Bit sampling method must not rely strongly on any property of the measured physical process other than its randomness, at least in some easily achievable limit;
2. Generator should withstand reasonable tolerance in components and operating conditions (eg. supply voltage, temperature, EM noise) without the need for (re)calibration or compensation;
3. Possible malfunctions during the lifetime of the generator should be foreseen and checked for at each generator restart (for example like recommended in FIPS140-1 ) or even continuously;
4. Sequences of bits produced by the generator should pass, with a high probability, any known statistical randomness test.

In the generator described here physical processes in Zener diode serve as a source of randomness. The noise voltage at terminals of the diode is ”measured” by a special electronic circuit at regular time intervals. Each measurement results in a random bit.

The circuit presented in Fig. 1 follows the method of operation described in the previous section. It consists of the following five blocks:

1. A source of electric noise;
2. a DC decoupling capacitor $`C`$;
3. a circuit for digitizing the noise voltage consisting of a comparator controlled with a rough automatic zero-bias correction circuitry;
4. a counting circuit (JK flip-flop);
5. a sampling circuit which delivers a bit upon an external request.

It is a well known fact that a Zener diode operating in a reverse polarity and the current strength near the knee, can serve as a noise generator. For example, a 6.8 Volt commercial Zener diode can produce a noise fluctuation of an amplitude of 30 to 50 mV (peak to peak) with a mean frequency of zero crossings of the order of 10 MHz Somlo75 .

It is also well known that a Zener diode with a knee voltage of less than 6.2 V operates mainly in the Quantum Mechanical (tunneling) regime, while the diodes with the knee above that operate mainly in the micro plasma regime Somlo75 . Both regimes have ideal properties of unpredictability needed for a truly random noise voltage source.

The best temperature stability of the noise amplitude is also obtained in the Zener diodes with the knee voltage of approx. 6.2 V, because in this case the two regimes with opposite temperature coefficients are in equilibrium Somlo75 . The same condition is also optimal from the point of the long term stability.

The generator functions in the following way. Noisy voltage from the Zener diode is AC decoupled from the digitizing stage which consists of a comparator COMP controlled with a rough automatic zero-bias feedback control. The role of the comparator is to convert tiny analog noise to the digital signal suitable for further processing.

The capacitor $`C`$ in the series with the output resistance $`R`$ of the noise source and the effective input resistance $`R_i`$ of the comparator forms a major contribution to an unwanted ”memory”. The timely persistence $`\tau _e`$ of this memory is equal to:

$`\tau _e=C(R+R_i).`$ (6)
Voltage amplitudes of any two noise variations which happen within the period $`\tau _e`$ will be mutually correlated because of the electric charge in the capacitor $`C`$ which has no time to discharge through the resistance in the system. If physical events are not completely independent of each other there will be another persistence $`\tau _p`$ giving rise to the total memory:

$`\tau =\tau _e+\tau _p`$ (7)
Luckily, this ”memory” effect dies off exponentially with the time distance between the two variations, thus one can conclude that any two variations that are distant enough in time may be considered as statistically independent. Whenever this applies, our method is valid, as will be explained later. Nevertheless, the ”memory” of the circuit limits frequency bandwidth of the noise and sets an absolute upper limit to the bit extraction rate from the generator, which limit is independent of the latter bit extraction method.

The negative input of the comparator COMP is connected to a suitable DC reference voltage $`V_r`$. Between the two inputs of the comparator (positive and negative) a small DC ”offset” voltage can be induced by virtue of the resistor $`R_i`$ and the control current $`I_c`$ that flows through it. As the result, the positive input of the comparator ”sees” the sum of the offset voltage and the noise voltage. Whenever the sum exceeds the reference voltage $`V_r`$, the output of the comparator goes into the high logical state ”1”, whereas when the sum goes below the $`V_r`$, the output goes into the low logical state ”0”. This is illustrated in the Fig. 2

By setting the control current $`I_c`$ a little lower or higher, the output of the comparator COMP spends respectively a little more or less time in the logical state ”1”. At certain strength of the control current $`I_c`$, one can achieve that the comparator COMP spends approximately equal amount of time in either of the logical states, the point at which the comparator has the highest efficiency of detecting voltage variations. In other words, the duty cycle of the comparator COMP would then on average be close to 0.5. This is ensured by setting $`V_{1/2}=\frac{1}{2}V_{DD}`$. If, for any reason, the duty cycle ever gets changed, the feedback network closed by comparator COMP, operational amplifier OPA and low-pass filter LP (Fig. 1) will respond by changing the control current $`I_c`$ so that the duty cycle of 0.5 will be restored. This behavior effectively solves the requirement 2. Already here, at the output of the comparator COMP, there is an approximately equal chance that the state is at logical ”0” or at logical ”1”, i.e. the bias of the sampled output would be close to zero. However, due to technical reasons it is quite difficult to keep the bias below 1/1000 for long term and even this would be not possible without some fine-tuning. Tuning the bias consumes a lot of time (due to its statistical nature) and would therefore be unfavorable for a mass production of the generator.

Present method for extracting random bits eliminates the need to tune the bias to zero value, and actually allows to achieve bias as low as desired, without the need to do any modifications to the circuit.

Event counting is the crucial point of the method. Namely, the JK-type flip-flop (Fig. 1) performs a continued counting of the events, but keeps track only whether the count is even or odd. The result of the counting appears at the output Q of the JK flip-flop and represents a new random bit sequence with highly suppressed bias and serial correlation. To understand how this works it is important to keep in mind that the output Q of the JK flip-flop shall be sampled periodically in time. This is a condition for good operation of this generator. The sampling is done by the D-type flip-flop.

Without the loss of generality, let us suppose that at the time zero ($`t=0`$) the Zener diode Z exhibits a voltage breakdown. (Voltage breakdown is only one of the processes that may cause fluctuations of the voltage across the diode that are manifested as noise. But to simplify the language we will refer to ”voltage breakdown” as a synonym for a sudden positive jump of the voltage across the Zener diode.) This causes the output of the comparator COMP to exhibit a positive going transition (and later a negative one). Every such event triggers the JK flip-flop to reverse its state at the output Q. This is illustrated in the Fig. 2.b and 2.c.

Let us suppose that the voltage breakdowns happen at times $`t_0`$ = 0, $`t_1`$, $`t_2`$, $`t_3`$, $`t_4`$, etc. Tiny intervals of time between neighboring voltage breakdowns $`\mathrm{\Delta }t_k=t_kt_{k1}`$ where $`k`$ = 1, 2, 3… are distributed according to some statistical distribution. The mean value $`\mu `$ of duration of this tiny intervals is defined as:

$`\mu =\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}t_k`$ (8)
For example, in the pure Quantum Mechanical regime of the Zener diode, intervals between neighboring voltage breakdowns $`\mathrm{\Delta }t_k`$ follow the Exponential distribution. Mixing of the Quantum Mechanical model with the plasma noise, noise sources other than a Zener diode or other models including effective memory and/or filtering of the noise signal prior to digitization, may lead to more centered distributions such as Poisson or even Uniform-like distribution. The power of this method lies in the fact that the knowledge of the actual distribution of the time intervals between neighboring voltage breakdowns is irrelevant. The only important thing is that the shape of the distribution stays stable over a period of time substantially larger than the sampling period $`T`$.

The method presented in this work consists in sampling (i.e. reading out) of the output Q of the JK flip-flop periodically, at times $`T`$, $`2T`$, $`3T\mathrm{}`$ and generally at times $`i\times T`$ with $`i=1,2,3\mathrm{}`$ . The sampling is made by a D-type flip-flop (Fig. 1). Namely, when a positive-going edge of the sampling signal (that is transition from ”0” to ”1”) appears at the ”Request” input, state that is present at the input D gets sampled. Sampled value is memorized (frozen) and displayed at the output Q of the D-type flip-flop, and stays there until a new positive-going edge appears. According to the definition of the D-type flip-flop, the sampling is done almost instantaneously so its output stays frozen for almost the whole sampling period $`T`$ which makes possible for other devices (such as a computer) to read it safely.

Let us now suppose that the sampling period $`T`$ is chosen large enough so that many voltage breakdowns occur during each sampling period. In another words, we suppose:

$`T\mu .`$ (9)
This situation is illustrated in the Fig. 2 where several voltage breakdowns in the noise signal (FIG.2.a) happen during any sampling period of length $`T`$ (FIG 2.d).

Furthermore, let us denote that the JK flip-flop changes its state at the output Q: $`N_0`$ times from $`t=0`$ till $`t=T`$, $`N_1`$ times from $`t=T`$ till $`t=2T`$, and generally $`N_i`$ times from $`t=iT`$, till $`t=(i+1)T`$, where $`i=1,2,3,\mathrm{}`$ .

Now we can conclude that every sampling period (of length $`T`$) approximately equals the sum of many tiny intervals $`\mathrm{\Delta }t_k`$:

$`T`$ $``$ $`{\displaystyle \underset{k=1}{\overset{N_0}{}}}\mathrm{\Delta }t_k`$
$`T`$ $``$ $`{\displaystyle \underset{k=N_0+1}{\overset{N_0+N_1}{}}}\mathrm{\Delta }t_k`$
$`\mathrm{etc}.`$ (10)
The above equalities hold only approximately because the first and the last of the tiny intervals in a given sum may be only partially happening during the respective sampling period $`T`$. In such a case, a part of the first tiny interval may have actually happened during the previous period, while a part of the last tiny interval may have actually happened during the next period (sum). However, approximation can be made arbitrarily good by setting the large enough sampling period $`T`$.

The Central limit theorem of the Statistics states that the sum of a large number (say $`N`$) of independent random variables (for example $`\mathrm{\Delta }t_k`$), which follow some (any) distribution, approaches Normal distribution in the limit of large $`N`$. This means that the statistical variable $`x`$ defined as:

$`x={\displaystyle \underset{k=1}{\overset{N}{}}}\mathrm{\Delta }t_k,`$ (11)
follows approximately the Normal distribution for large $`N`$. However, in our case, defined by equations (10), the ”variable” $`x`$ is fixed (being just $`T`$), and the relevant statistical random variable becomes the number of summands, $`N`$ defined as:

$`{\displaystyle \underset{k=1}{\overset{N}{}}}\mathrm{\Delta }t_kx<{\displaystyle \underset{k=1}{\overset{N+1}{}}}\mathrm{\Delta }t_k.`$ (12)
In the special case mentioned above, when the tiny intervals $`\mathrm{\Delta }t_k`$ are distributed according to the Exponential distribution, random variable $`N`$ follows, by definition, the Poisson distribution:

$`P(N)={\displaystyle \frac{(T/\mu )^N}{N!}}exp(T/\mu ).`$ (13)
Quite generally, regardless of the distribution function of tiny intervals $`\mathrm{\Delta }t_k`$, it can be shown Stip  that the integer random variable $`N`$ is distributed according to the Binomial distribution $`B(N;2T/\mu ,p)`$, where $`p`$ quickly approaches 0.5 as $`T/\mu `$ tends to infinity. Given the general expression for the Binomial distribution:

$`B(r;n,p)=\left(\begin{array}{c}n\\ r\end{array}\right)p^r(1p)^{nr}.`$ (16)
distribution of the variable $`N`$ can be written as:

$`B(N;2T/\mu ,p)=\left(\begin{array}{c}2T/\mu \\ N\end{array}\right)p^N(1p)^{2T/\mu N}.`$ (19)
This is not in a contradiction with the previous result because the Poisson distribution (13) becomes equal to the Binomial distribution (19) in the limit $`T/\mu \mathrm{}`$.

Even without a rigorous proof it is easy to understand why $`N`$ follows symmetric Binomial distribution (19) in the limit of slow sampling or, equivalently, large $`N`$. Following eq. (12), for large $`N`$ we can write:

$`x{\displaystyle \underset{k=1}{\overset{N}{}}}\mathrm{\Delta }t_k.`$ (20)
Approximating the definition (8) with:

$`\mu {\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}\mathrm{\Delta }t_k`$ (21)
to eliminate the sum in (20) we may conclude:

$`N{\displaystyle \frac{1}{\mu }}x.`$ (22)
Since $`x`$ follows Normal p.d.f. so should $`N`$. But since $`N`$ is an integer number it actually follows a discrete version of the Normal p.d.f., i.e. symmetric Binomial distribution with $`p`$=0.5, Q.E.D. This is a very important result because it means that (19) holds regardless of the distribution of tiny intervals $`\mathrm{\Delta }t_k`$. This solves the requirement 1.

It can be shown Stip  that as a special consequence of the equation (19) the probability that the random variable $`N`$ is even becomes equal to the probability that it is odd, in the limit $`T/\mu \mathrm{}`$. But, when the JK flip-flop changes its state an even number of times its output Q returns to the initial state (say ”0”) whereas when it changes an odd number of times it returns the complementary state (say ”1”). Thus, in the limit of slow sampling, the probability of producing ”0” would be equal to the probability of producing ”1” at the output Sample (Fig. 1). This fact can be expressed like this:

$`T/\mu \mathrm{}b0.`$ (23)
Thus it has been shown that the generator and the corresponding method of sampling random phenomena solve the first problem stated in the section 2, i.e. that the bias should vanish in some limit.

The second problem, namely the statistical independence of the logical states (bits) at the output of the generator, is solved as follows. First we realize that the subsequent voltage breakdowns in the Zener diode are statistically independent because the breakdown itself is an unpredictable physical phenomenon. Secondly, the internal capacitive memory that appears in the circuit (which for example may be introduced by imperfections in electronics design, or be a part of physical process governing the alternate current noise source) ”dies off” exponentially with the half-life $`\tau `$ (Eq. (7)). Effect of the memory on correlations can be made as small as desired in the limit of slow sampling, that is by making T/$`\tau `$ large enough. Namely, in the limit $`T/\tau \mathrm{}`$ sums of the type (10) become mutually statistically independent. The statistical independence of two subsequent sums, in that limit, takes place because the two sums contain large portion of summands that are ”out of the range” of the memory effect of each other. If we consider two sums that are not subsequent, then the memory effects are even smaller. This means that a bit corresponding to a sum becomes statistically independent of neighboring and even more all other bits, in the limit of slow sampling. It can be concluded that the generator and the corresponding method of sampling random phenomena also solve the second problem stated in the section 2.

Functioning of the whole circuit in Fig. 1 can be best understood by looking at signals shown in the Fig. 2a-e. One can see that with each positive-going edge at the output of the comparator COMP (Fig. 2.b), the JK flip-flop changes state at its output Q (Fig. 2.c). The sampling signal (Fig. 2.d) isn’t, of course, in any synchronization with these signals. Whenever it exhibits a positive-going transition, the output Q of the JK flip-flop gets sampled and the state thus obtained is the output bit from the generator (Fig. 2.e).

It is important to understand that the output Q of the JK flip-flop always produces ”1” after ”0” and ”0” after ”1”, therefore its output is not random at all. However a duration of zeros and ones varies randomly (Fig. 2c) and periodic sampling of the said output by the D-type flip-flop provides a good quality random bits, in the limit of slow sampling ($`T\mu `$ and $`T\tau `$).

## 4 Simulations

In order to check the theory of operation described above, a series of simulations of the circuit in Fig. 1 were made. Simulations should provide an insight on how the type of distribution of tiny intervals between subsequent physical events and ratios $`T/\tau `$ and $`T/\mu `$ influence the quality of generated bits. The circuit was simulated by a specifically designed computer program which simulates a random digital signal obeying a specified distribution which would be present at the output of the comparator COMP, and assuming that this signal gets processed by ideal flip-flops. The approximation of considering the two flip-flops to be ideal is very good because their fall and rise times (of approximately 7 nanoseconds) are much shorter than the changes in the noise being processed.

Two distributions of tiny intervals were simulated: Exponential and Uniform. Exponential distribution is natural distribution for many processes. Uniform distribution assumes that the intervals between physical events may be between zero and 1 (in arbitrary time units, for example microseconds) and evenly distributed between the two limits. The two distributions are chosen because they are very different: for example Exponential d. is completely non-symmetric and peaked at one value (zero) whereas Uniform d. is symmetric and has no peak, etc. Most processes should obey distributions that are somewhere between those two extremes.

The effect of memory was simulated as a pile up: any two breakdowns closer (in time) than $`\tau `$ were considered as being just one breakdown happening at the end of the last breakdown. This in effect may chain up more than two events to be seen as one. This is a very faithful model for a noise source which we use, and a reasonable approximation for any other noise source because memory (by definition) always piles up events in a way that events which happen to close to one another get piled up and be detected as just one event.

The method for extracting bits does not prefer zeros nor ones, therefore we do not expect the bias to be significant even in the case when both ratios $`T/\tau `$ and $`T/\mu `$ are small. What happens however is that when sampling is being made too fast there is an excess of longer patterns of ones or zeros. This in turn means that the probability of single ones (pattern ”010”) a well as that of single zeros (pattern ”101”) will be smaller than expected (expected vales are: $`p(010)=p(101)=1/8`$).

For each choice of distribution, ratio $`T/\tau `$ and ratio $`\tau /\mu `$, the program calculates (by simulation) three output values: probability of ones $`p(1)`$, probability of single ones $`p(010)`$ and serial correlation coefficient $`a`$. Results of simulations for the two distributions and a range of ratios $`T/\tau `$ and $`\tau /\mu `$ are shown in the tables 3-6. Values shown are bias ($`b=p(1)0.5`$), second order bias ($`b_{010}=p(010)0.125`$), and serial correlation coefficient. Number of bits generated in each simulation is 6,250,000. According to that statistics, 1 sigma errors are: 0.00020 for $`b`$, 0.00007 for $`b_{010}`$ and 0.00040 for $`a`$. As a source of random numbers the program used semi-hardware ”/dev/urandom” Linux kernel random bit generator. The result were cross-checked using the standard implementation of ”secure” software random number generator rand() that is part of the Linux C compiler (we used gcc version 2.95). No discrepancies were found in the results within the error margins.

As expected, bias is always consistent with zero. Absolute values of the second order bias and the serial correlation coefficient quickly drop to zero (within the error margins) when the sampling period becomes large enough and memory effects small enough. The Exponential distribution seems slightly more favorable in this respect than the Uniform. From these results it can be concluded that in the worst case good bits are achieved by generators with the memory effect $`\tau /\mu <0.33`$ and the ratio $`T/\mu >9`$. This means that using a standard Zener diode with 3 MHz mean frequency of voltage turnovers it should be possible to build a generator that is capable of generating up to several hundred kilobits of good-quality random bits per second.

## 5 Prototype testing

We have constructed a prototype based on the block-diagram in Fig. 1 using a Zener diode, standard analog chips and CMOS logic. Bits were produced at 300 kbit/sec. We built four identical circuits and fed their outputs to a PC computer, thus obtaining a total of 1.2 Mbit/sec. Thanks to such a high speed we were able to easily produce long sequences of bits for the subsequent testing. The sequences of bits produced by the physical generator were tested by three batteries of statistical tests: Walker’s ENT ENT , Marsaglia’s DIEHARD DIEHARD  and NIST’s Statistical Test Suite sts-1.50 .

The ENT consists of battery of standard tests: entropy, Chi-square, mean value, Monte Carlo pi value test and the serial correlation test. These tests may be evaluated for short sequences thus making possible a fast check of the prototype in the development stage. The most important of them are the mean value and the serial correlation test because they best grasp problems in hardware. The mean value addresses the inequality between number of ones and zeros which is very hard to keep close to zero with imperfect hardware. Non zero serial correlation reflects existence of a short-term memory in the system, as explained in the section 3. Indeed, at the development stage of the prototype a slight consistent bias (positive, of the order of $`10^4`$) has been noticed on all four gadgets which was attributed to the fact that comparator had quite a different rise time from the fall time at its output. This was corrected by taking another comparator after which the bias became undetectable.

Once the generator was debugged and ENT tests passed well, we proceeded by generating long sequences (approx. 10<sup>8</sup> bits) needed by the DIEHARD, probably the stringiest battery of tests known today. Table 7 shows a typical result of testing of $`8\times 10^7`$ bits (10 megabyte) of random data obtained from the generator.

According to Marsaglia DIEHARD  sequence has failed a test if the p-value of that test is very close to 1 ($`p>0.999999`$), otherwise the sequence has passed the test. There are 15 different tests in the DIEHARD, but some are run more than once with different parameters (eg. binary rank test) or on independent parts of the sequence under test (eg. count the 1’s). Results in the table 7 indicate a random sequence that has passed all the tests.

It is interesting to note that efforts to standardize requirements for random bit generators for use with cryptographic software have been taken, notably the recommendations of NIST FIPS140-1 . A sequence of 10<sup>9</sup> bits divided in 100 equally long pieces (files) has been subjected to the NIST’s battery of randomness tests, version sts-1.50 sts-1.50 . We used the recommended test parameters. Results are summarized in the Table 8.

According to the NIST analysis program, the minimum pass rate for each statistical test with the exception of the random excursion (variant) test is approximately 0.9602 whereas the minimum pass rate for the random excursion (variant) test is approximately 0.9527. Taking this into account we conclude that all 16 tests were passed with excellent marks. Passing all three batteries of tests shows that this generator conforms also to the general requirement 4 stated in the section 3.

As one more check of the sampling method itself bits were generated from a cardiogram data of a healthy patient. A record of 1190 consequent heartbeats was made in course of an ergometric stress test. During the test patient is required to walk on a treadmill whose speed varies according to a standardized procedure. Test lasts approximately 12 minutes. The times at which main heartbeat peaks occur were considered as physical ”events”. Taking $`T/\mu 5`$ we obtained a sequence of 216 bits using our method. Results of ENT test shown in the table 9 indicate excellent randomness. This demonstrates the power of the method which is able to extract good random bits form a sequence of heartbeats which are nearly periodic but have only a slight random jitter.

The conclusion is that heartbeat timings may very well be used as a base for a non-deterministic random generator !

## 6 Conclusion

A method for extraction of non-deterministic random bits from random physics processes (such as nuclei decay or noise voltage variations) is presented. The method is characterized by the fact that it ensures a good quality random bits in an easily obtainable limit, that is the limit of slow sampling, regardless of the distribution of times between adjacent processes. A physical random bit generator which makes use of the method and Zener diode noise was built and successfully tested. The generator is shown to conform to general requirements for non-deterministic generators stated in section 3. A sample of approximately $`10^9`$ bits produced by the physical generator prototype is subjected to a series of well-known statistical tests showing no weaknesses.

## 7 Acknowledgments

The heartbeat recording was provided by courtesy of M. Martinis from Rudjer Bošković Institute and the Institute for cardiac deseases, Zagreb. These data were taken with approval of the patient.

Parts of the method and apparatus for producing random bits are subjected to a patent procedure.


text:
# Electromagnetic Mechanism for the Origin of Knots in Parsec Scale Jets

### Acknowledgments.

The authors are grateful to D.C. Gabuzda for helpful discussions. This work was done under the partial support of the Russian Foundation for Fundamental Research (grant number 02-02-16762). Support from DOE grant DE-FG02-00ER54600 is acknowledged.

## References

Attridge J.M., Roberts D.H., & Wardle J.F.C. 1999, ApJ, 518, L87

Beresnyak, A.R., Istomin, Ya.N., & Pariev, V.I. 2003, A&A, 403, 793

Gabuzda D.C. 2000, in Conway J.E., Polatidis A.G., Booth R.S., Pihlström Y.M., eds, EVN Symposium 2000, Proceedings of the 5th European VLBI Network Symposium. Onsala Space Observatory, Onsala, p.53

Gabuzda, D.C., & Murray, É. 2003, this volume

Harris D.E., et al. 2003, ApJ, 586, L41

Istomin Ya.N., Pariev V.I. 1994, MNRAS, 267, 629

Istomin Ya.N., Pariev V.I. 1996, MNRAS, 281, 1

Leppänen K.J., Zensus J.A., & Diamond P.J. 1995, AJ, 110, 2479

Lister M.L. 2001, ApJ, 562, 208

Marshall H.L. et al. 2001, ApJ, 549, L167

Pariev, V.I., Istomin, Ya.N., & Beresnyak, A.R. 2003, A&A, 403, 805


text:
# On the extension of the Bethe-Weizsäcker mass formula to light nuclei

## Abstract

Some general features of the Bethe-Weizsäcker mass formula recently extended to light nuclei have been explored. Though this formula improves fits to the properties of light nuclei and it does seem to work well in delineating the positions of all old and new magic numbers found in that region, yet it is not well tuned for predicting finer details. The mass predictions have also been found to be less accurate compared to those by the macroscopic-microscopic calculations. It is concluded that such semi-empirical mass formulae can not be a substitute for more fundamental mass formulae having its origin based upon the basic nucleon-nucleon effective interaction.

In a recent article  the Bethe-Weizsäcker (BW) mass formula has been extended to light nuclei and some new shell closures have been identified. This modified BW mass formula explains the shapes of the binding energy versus neutron number curves of most of the elements from $`Li`$ to $`Bi`$. For heavier nuclei this modified BW mass formula approaches the old BW mass formula. This modified BW formula suggests additional stabilities for few light nuclei at neutron numbers N=16 (Z=7,8), N=14 (Z=7-10), Z=14 (N=13-19), N=6 (Z=3-8) and loss of magicity for nuclei with neutron numbers N=8 (Z=4) and N=20 (Z=12-17), most of which were suggested earlier from the systematics. The magicity at neutron number N=6 (Z=3-8) suggested by the new formula is also supported by the experimental r.m.s. radii values that show a quenching at N=6 (Z=3-8). The new formula supports already known $`{}_{}{}^{32}Ne`$, $`{}_{}{}^{35}Na`$, $`{}_{}{}^{38}Mg`$ and $`{}_{}{}^{41}Al`$ as the last bound isotopes of neon, sodium, magnesium and aluminium, respectively. The existing mass formulae  support the existence of old magic numbers but fail to explain the general trend of the binding energy versus neutron number curves for several light nuclei near the drip lines and evidence of new magic number at neutron number N=16. Wherever there is a signature of shell closure in experimental data, the BW as well as the new modified BW formula both show marked deviations delineating clearly the positions of old magic numbers at 2, 8, 20, 28, 50, 82 and 126.

The binding energy B(A,Z) of any nucleus of mass number A and atomic number Z was obtained from a phenomenological search and was given by a more genaralized BW formula

$$B(A,Z)=a_vAa_sA^{2/3}a_cZ(Z1)/A^{1/3}a_{sym}(A2Z)^2/[(1+e^{A/k})A]+\delta _{new},$$
(1)
where the asymmetry term and the pairing term had been modified to obtain great improvements to the fits of the binding enegy versus neutron number curves of the light nuclei. Fitted value of the asymmetry term modifying constant $`k=17`$, while the modified pairing term $`\delta _{new}`$ can be expresssed as

$$\delta _{new}=(1e^{A/c})\delta wherec=30,$$
(2)
and $`\delta `$ being the old pairing term is given by

$`\delta =`$ $`a_pA^{1/2}forevenNevenZ,`$ (3)
$`=`$ $`a_pA^{1/2}foroddNoddZ,`$ (4)
$`=`$ $`0foroddA,`$ (5)
while the other constants

$$a_v=15.79MeV,a_s=18.34MeV,a_c=0.71MeV,a_{sym}=23.21MeVanda_p=12MeV,$$
(7)
remaining the same as the corresponding values of the constants for the old BW formula. Such modifications, however, do not alter significantly the results for heavier nuclei.

In order to investigate the features of this newly extended Bethe-Weizsäcker mass formula, the one and the two proton separation energies, the one and the two neutron separation energies and the atomic mass excesses have been calculated using the formula given by Eq.(1) with the $`\delta _{new}`$ quantity described by Eq.(2). The best fit values used for the new constants are $`k=17`$ and $`c=30`$ while the values of the older constants used for calculating the binding energies are given in Eq.(4). The one and two proton separation energies $`S_{1p}`$ and $`S_{2p}`$ defined as the energies required to remove one proton and two protons from a nucleus are given by

$$S_{1p}=B(A,Z)B(A1,Z1)andS_{2p}=B(A,Z)B(A2,Z2)$$
(8)
respectively, while the one neutron and two neutron separation energies $`S_{1n}`$ and $`S_{2n}`$ defined as the energies required to remove one neutron and two neutrons from a nucleus are given by

$$S_{1n}=B(A,Z)B(A1,Z)andS_{2n}=B(A,Z)B(A2,Z)$$
(9)
respectively. The reason that the atomic mass is considered rather than the nuclear mass is that historically, the former has been the actual experimentally measured quantity, whereas the latter is less accurate bacause its extraction requires a knowledge of binding energy of the Z atomic electrons. However, recent developments now allow the nuclear masses to be measured directly . For those applications where it is necessary to know the actual mass of the nucleus itself, its value (in $`MeV`$) can be found from the atomic mass excess or from nuclear binding energy by use of the relationship

$`M_{nucleus}(A,Z)=`$ $`Au+\mathrm{\Delta }M_{A,Z}Zm_e+a_{el}Z^{2.39}`$ (10)
$`=`$ $`Zm_p+(AZ)m_nB(A,Z),bydefinitionofnuclearbindingenergyB(A,Z)`$ (11)
where $`m_p`$, $`m_n`$ and $`m_e`$ are the masses of proton, neutron and electron respectively, the atomic mass unit u is 1/12 the mass of $`{}_{}{}^{12}C`$ atom, $`\mathrm{\Delta }M_{A,Z}`$ is the atomic mass excess of an atom of mass number A and atomic number Z and the electronic binding energy constant $`a_{el}=1.433\times 10^5`$ MeV. Hence from Eqs.(7) the atomic mass excess is given by

$`\mathrm{\Delta }M_{A,Z}=`$ $`Zm_p+(AZ)m_n+Zm_eAua_{el}Z^{2.39}B(A,Z)`$ (12)
$`=`$ $`Z(m_p+m_eu)+(AZ)(m_nu)a_{el}Z^{2.39}B(A,Z)`$ (13)
$`=`$ $`Z\mathrm{\Delta }m_H+(AZ)\mathrm{\Delta }m_na_{el}Z^{2.39}B(A,Z)`$ (14)
where $`\mathrm{\Delta }m_H=m_p+m_eu`$ = 7.289034 MeV and $`\mathrm{\Delta }m_n=m_nu`$ = 8.071431 MeV.

Though such modifications certainly cause impovements over earlier semi-empirical mass formulae, but it is still beset with many difficulties. The new formula suggests that the $`{}_{}{}^{31}F`$ nucleus is unstable while experimental evidence of its existence has been found. However, similar instability was predicted by other models also. The formula wrongly predicts $`{}_{}{}^{5}Li`$, $`{}_{}{}^{6}Be`$, $`{}_{}{}^{15}F`$ and $`{}_{}{}^{16}Ne`$ as stable , and puts the neutron drip line at $`{}_{}{}^{17}B`$, $`{}_{}{}^{20}C`$ and $`{}_{}{}^{26}O`$  whereas $`{}_{}{}^{19}B`$, $`{}_{}{}^{22}C`$ are in fact also stable, but $`{}_{}{}^{25,26}O`$ are not. While it may be used to provide gross features like predicting approximate positions of the drip lines, it must be realised that a semi-empirical formula such as the Bethe-Weizsäcker mass formula can not be well tuned for predicting finer details.

The calculated theoretical atomic mass excesses using the macroscopic-microscopic model  and the semi-empirical modified Bethe-Weizsäcker (BW) mass formula  denoted by $`\mathrm{\Delta }M_{MS}`$ and $`\mathrm{\Delta }M_{BW}`$ respectively, have been calculated using Eq.(8) while the atomic mass excesses denoted by $`\mathrm{\Delta }M_{Ex}`$ are the experimentally measured quantities. The chi-square per degree of freedom $`\chi ^2/F`$ can be used for a relative comparison between different theoretical models which is defined as

$$\chi ^2/F=(1/N)[(\mathrm{\Delta }M_{Th}\mathrm{\Delta }M_{Ex})/\mathrm{\Delta }M_{Ex}]^2$$
(15)
where $`\mathrm{\Delta }M_{Th}`$ is the theoretically calculated atomic mass excess and the summation extends to $`N`$ data points for which experimental atomic mass excesses are known. The $`\chi ^2/F`$ calculated for the semi-empirical Bethe-Weizsäcker (BW) mass formula using 1742 experimental atomic mass excesses has been found to be 12.936. The $`\chi ^2/F`$ for the macroscopic-microscopic model using 1727 experimental atomic mass excesses has been found to be 1.821. Only one point corresponding to the $`{}_{}{}^{12}C`$ has been excluded from the summation in both the calculations for $`\chi ^2/F`$ because the experimental atomic mass excess for the $`{}_{}{}^{12}C`$ atom is, by definition, equal zero and can not be accomodated inside the summation given by Eq.(9). It is evident that the macroscopic-microscopic model  which uses momentum and density dependent Seyler-Blanchard  type effective interaction provide better mass estimates than the semi-empirical modified Bethe-Weizsäcker (BW) mass formula.

To conclude, the extension of Bethe-Weizsäcker mass formula to light nuclei is beset with some difficulties. There are predictions about the stability of some light nuclei which are not in agreement with the experimental observations. Mass estimates also show clear superiority for more fundamental mass formulae having its origin based upon the basic nucleon-nucleon effective interaction. The need for improvements of mass formula having its origin based upon the basic nucleon-nucleon effective interaction is stressed rather than modifying semi-empirical mass formula, such as the Bethe-Weizsäcker, which can not be tuned to reproduce finer details. Ideally, a density dependent effective interaction  which provides unified descriptions of the nuclear scattering , cluster radioactivity  and nuclear matter  may be used for further improvements in the macroscopic-microscopic model.